;reads dat files, creates alfvenicity graph
;includes a day interval input over readheliosdatdays-cj.pro

;initialize variable types
filein = '' & filebase = '' & scid=1 & fileyear=1 & ndaybeg=1 & ndayend=1 & choice=-1

device,decomposed=0
tek_color
makevb = 0  ; 1 to make ViSBARD files
close,1
read, scid,prompt = 'Enter Helios number (1 or 2)'
read, fileyear, prompt = 'Enter year 19XX'
read, ndaybeg,prompt = 'Enter first day'
read,ndayend,prompt = 'Enter last day'
;print, '1 - Original values'
;print, '2 - Absolute values'
;read, choice, prompt = 'Enter choice'
choice = 1
;for ii = 76, 76 do begin
    ;fileyear = strcompress(string(ii),/remove_all)

    i=0L
    npts=0L
    for kk = ndaybeg, ndayend do begin
       counter=1
       ;in if statement, if kk eq to first day tested
       if kk eq ndaybeg then counter=0 else counter=1

       if (kk lt 10) then fileday='00'+strcompress(string(kk),/remove_all)
       if (kk ge 10 and kk lt 100) then fileday='0'+strcompress(string(kk),/remove_all)
       if (kk ge 100) then fileday = strcompress(string(kk),/remove_all)
       filebase='H'+strcompress(string(scid),/remove_all)+strcompress(string(fileyear),/remove_all)+'_'+fileday
       ;filein = '\Documents and Settings\Maria Luckyanova\My Documents\NASA\HeliosData\'+filebase+'.dat'
       filein = '/Users/hancj/Documents/IDL/Data/Helios' +strcompress(string(scid),/remove_all)+'_Dat_full/'+filebase+'.dat'
       ;fileout= '/Data/ViSBARD/Helios'+scid+'_hires/'+filebase+'vb.dat'
       on_ioerror, endloops
       openr,1,filein

       m=0L
       dummy = ' '
       while not eof(1) do begin
         readf,1,dummy
         i=i+1
         m=m+1
       endwhile

       nptsm = m-1
       npts=npts+nptsm
       close, 1

       print,npts, nptsm



       if (counter eq 0) then begin
         stringin = ' '
         datain = strarr(npts)
         data = fltarr(34,npts)

       openr,1,filein
       readf,1,stringin
       vars=strsplit(stringin,' ',/extract)

       readf,1,datain

         for j = 0, nptsm-1 do begin
          test=strsplit(datain(j),' ',/extract)
          data(*,j) = float(test)
         endfor

       endif

       if (counter eq 1) then begin
         stringin = ' '
         datain = strarr(nptsm)
         data2 = fltarr(34,nptsm)
         ;data = data2

              openr,1,filein
       readf,1,stringin
       vars=strsplit(stringin,' ',/extract)

       readf,1,datain
       endif



       if (counter eq 1) then begin
         for j = 0, nptsm-1 do begin
              test=strsplit(datain(j),' ',/extract)
              data2(*,j) = float(test)
         endfor
       endif

       if counter eq 1 then data = [[data],[data2]]

       close,1

       ; Assign data values to the variable names.
       for k = 0, 33 do begin
          qq=execute(vars(k)+'=data('+string(k)+',*)')

       endfor

       vr = vp1r
       vt = vp1t
       vn = vp1n
       den = np1
       temp = tp1
       r = rh

       print, year(0,0), day(0,0), npts

       endloops: continue
    endfor ; on kk
;endfor ; on ii

nside = 5


int = lindgen(npts)

indat=fltarr(14,i)

;openr,1,ifile
;readf,1,indat
;close,1

br=fill(br)
bt=fill(bt)
bn=fill(bn)
vr=fill(vr)
vt=fill(vt)
vn=fill(vn)


jj=1L

 jj=where(indat eq 999.9 or indat eq 9999. or indat eq 99.99999 or indat eq 9999999.,count)
 ;if count ne 0 then indat(jj) = 0

decday = day + dechr/24

br=transpose(br)
vr=transpose(vr)

bn=transpose(bn)
vn=transpose(vn)

bt=transpose(bt)
vt=transpose(vt)


dbr = br-smth(br,nside)
dvr = vr-smth(vr,nside)

dbt = bt-smth(bt,nside)
dvt = vt-smth(vt,nside)

dbn = bn-smth(bn,nside)
dvn = vn-smth(vn,nside)

dbr = dbr*21.8/sqrt(fill(den(int)))
dbt = dbt*21.8/sqrt(fill(den(int)))
dbn = dbn*21.8/sqrt(fill(den(int)))
;dvx = dvx*21.8/sqrt(fill(den1(int)))
topr = dbr*dvr
bottomr = dbr^2+dvr^2
sigcr = 2*smth(topr,nside)/(smth(bottomr,nside)+0.00001)

if(choice eq 1) then begin
;	window, /free
;	plot,transpose(decday), sigcr, yrange = [-2,2], title = 'Alfvenicity for Year 19' + $
;		strcompress(string(fileyear), /remove_all) + ' Day ' + $
;		strcompress(string(ndaybeg), /remove_all) + ' to Day ' + $
;		strcompress(string(ndayend), /remove_all)

	top = dbr*dvr + dbt*dvt + dbn*dvn
	bottom = dbr^2 + dvr^2 + dbn^2 + dvn^2 + dbt^2 + dvt^2
	sigc = 2*smth(top,nside)/(smth(bottom,nside)+0.00001)
;	oplot,transpose(decday),sigc,color=2
endif else begin
;	window, /free
;	plot,transpose(decday), abs(sigcr), yrange = [-2,2], title = 'Alfvenicity for Helios ' + $
;		strcompress(string(scid), /remove_all) + ' Year 19' + $
;		strcompress(string(fileyear), /remove_all) + ' Day ' + $
;		strcompress(string(ndaybeg), /remove_all) + ' to Day ' + $
;		strcompress(string(ndayend), /remove_all)

	top = dbr*dvr + dbt*dvt + dbn*dvn
	bottom = dbr^2 + dvr^2 + dbn^2 + dvn^2 + dbt^2 + dvt^2
	sigc = 2*smth(top,nside)/(smth(bottom,nside)+0.00001)
;	oplot,transpose(decday),abs(sigc),color=2
endelse

decdayt=transpose(decday)


;plots,[decdayt[0],decdayt[npts-1]],[0,0], linestyle=2

if (makevb eq 1) then begin
   phi = HGIlong*!pi/180.0
   theta = clat*!pi/180.0

   bx=br*cos(phi)*cos(theta)-bt*sin(phi)-bn*sin(theta)*cos(phi)
   by=br*sin(phi)*cos(theta)+bt*cos(phi)-bn*sin(phi)*sin(theta)
   bz=br*sin(theta)+bn*cos(theta)

   vx=vr*cos(phi)*cos(theta)-vt*sin(phi)-vn*sin(theta)*cos(phi)
   vy=vr*sin(phi)*cos(theta)+vt*cos(phi)-vn*sin(phi)*sin(theta)
   vz=vr*sin(theta)+vn*cos(theta)

   x=r*cos(theta)*cos(phi)
   y=r*cos(theta)*sin(phi)
   z=r*sin(theta)

   outdat = [year,day,dechr,x,y,z,bx,by,bz,den,temp,vx,vy,vz]

   openw,2,fileout
   printf,2,format='(14(f12.4,1x))',outdat

   print, 'Wrote ', fileout

   close,2
endif

print, 'Done'
end