;|-----------------------------------------------------------------|
;|                                                                 |
;|    Subroutine HDFREAD                                           |
;|                                                                 |
;|    reads DUCT and IDM data from UTD HDF files.                  |
;|                                                                 |
;|-----------------------------------------------------------------|

pro hdfread

@var.inc

filename=' '
name=' '
 readf,1,name

 name=strmid(name,1,7) 

 filename='/usr35/dynexp/hdf/duct/DE2HDF'+strtrim(number,2)+'/D'+strtrim(name,2)+'.HDF'
 print,filename

 HDF_DFSD_GETINFO, filename, dims = ddimsizes
 print,'data dimensions are ',ddimsizes
 ddata  = fltarr(ddimsizes(0),ddimsizes(1))
 HDF_DFSD_GETDATA, filename, ddata
 

 filename='/usr35/dynexp/hdf/drift/DE2HDF'+strtrim(number,2)+'/I'+strtrim(name,2)+'.HDF'
 print,filename

 HDF_DFSD_GETINFO, filename, dims = idimsizes
 print,'data dimensions are ',idimsizes
 idata  = fltarr(idimsizes(0),idimsizes(1))
 HDF_DFSD_GETDATA, filename, idata

 end;---------------------------------------------------------------------------------------

;|-----------------------------------------------------------------|
;|                                                                 |
;|    Subroutine FINDWATS                                          |
;|                                                                 |
;|    calculates which WATS file is needed.                        |
;|                                                                 |
;|-----------------------------------------------------------------|


pro findwats

@var.inc

; Calculate name for wats ascii files

year=1981
date=ddata(1,0)-81000.
if date gt 1000 then begin
	date=date-1000.
	year=1982
	endif
if date gt 1000 then begin
	date=date-1000.
	year=1983
	endif
date=fix(date)

end;-----------------------------------------------------------------------



;|-----------------------------------------------------------------|
;|                                                                 |
;|    Subroutine ASCIIREAD                                         |
;|                                                                 |
;|    reads DUCT and IDM data from UTD HDF files.                  |
;|                                                                 |
;|-----------------------------------------------------------------|

pro asciiread

@var.inc


tempdat=fltarr(11)
wdata=fltarr(22,10000)

j=0
while not eof(2) do begin

	readf, 2, tempdat, format ="(i6,i9,i4,i3,e12.5,f7.1,f7.1,f8.1,i6,i5,i7)"
	wdata(0:10,j)=tempdat(*)
	readf, 2, tempdat, format ="(f9.1,f8.1,i6,f7.1,f6.1,f7.1,f6.2,f6.2,f5.2,f6.1,f7.1)"
	wdata(11:21,j)=tempdat(*)
	j=j+1
endwhile
wcolumns=j

end;----------------------------------------------------------------------------------



;|-----------------------------------------------------------------|
;|                                                                 |
;|    Subroutine DUCT                                              |
;|                                                                 |
;|    extracts basic parameters from Ddata.                        |
;|                                                                 |
;|-----------------------------------------------------------------|

pro duct

@var.inc


time=fltarr(512*ddimsizes(1))
timehr=lonarr(512*ddimsizes(1))
timemin=lonarr(512*ddimsizes(1))
timesec=lonarr(512*ddimsizes(1))
N=fltarr(512*ddimsizes(1))
LogN=fltarr(512*ddimsizes(1))
DUTS=strarr(512*ddimsizes(1))


dut=fltarr(ddimsizes(1))
duthr=lonarr(ddimsizes(1))
dutmin=lonarr(ddimsizes(1))
dutsec=lonarr(ddimsizes(1))
dutx=strarr(ddimsizes(1))
;DUT=fltarr(ddimsizes(1))
DILAT=fltarr(ddimsizes(1))
DMLT=fltarr(ddimsizes(1))
DALT=fltarr(ddimsizes(1))
DGLAT=fltarr(ddimsizes(1))
DGLON=fltarr(ddimsizes(1))

icount=long(-1)
incr=1000./64. 
l_step=36             ; ddimsizes(1)..L_step is the Legend Step
Orbit='Orbit #: '+strcompress(fix(ddata(4,0)),/remove_all)
 
max3=-1.e3
min3=+1.e3

;print,'data dimensions are ',ddimsizes

for j=0,ddimsizes(1)-1 do begin
  for i=167,678 do begin
        icount=icount + 1
	time(icount)= ddata(2,j)+ddata(165,j)+(i-167)*incr
        timehr(icount)=fix(time(icount)/3600000.)
        timemin(icount)=(time(icount)/1000.-timehr(icount)*3600)/60
        timesec(icount)=time(icount)/1000.-(timehr(icount)*3600+timemin(icount)*60)
	DUTS(icount)=strcompress(timehr(icount),/remove_all)+':'+strcompress(timemin(icount),/remove_all)+':'+strcompress(timesec(icount),/remove_all)
        N(icount)= ddata(i,j)
	if (N(icount) gt max3) then max3=N(icount)
	if (N(icount) lt min3) then min3=N(icount)
  endfor
endfor 

for j=0,ddimsizes(1)-1 do begin
        dut(j)=ddata(2,j)
        duthr(j)=fix(dut(j)/3600000.)
        dutmin(j)=(dut(j)/1000.-duthr(j)*3600)/60
        dutsec(j)=dut(j)/1000.-(duthr(j)*3600+dutmin(j)*60)
	dutx(j)=strcompress(duthr(j),/remove_all)+':'+strcompress(dutmin(j),/remove_all)+':'+strcompress(dutsec(j),/remove_all)
        
      

endfor




for j=0,ddimsizes(1)-1 do begin
      DUT(j)=ddata(2,j)    ; UT
      DALT(j)=ddata(9,j)   ; Altitude
      DGLAT(j)=ddata(10,j) ; Latitude
      DGLON(j)=ddata(11,j) ; Longitude
      DMLT(j)=ddata(12,j)  ; MLT
      DILAT(j)=ddata(13,j) ; Invariant Latitude

endfor

    

; Get rid of these and see if it still works the same way...

DUTS=DUTS(0:icount-1)    ;String  (12:22:34)
N=N(0:icount-1)
;time=time(0:icount-1)   ; UT in Msec
;time2=time(81578:122367)

;  Keep the rest please 

LogN = alog10(N)
index=where(time gt 21100000 and  time lt 21500000 )
;index2=where(i gt 35 and i lt 85)
;      plot,time(index),N(index),yrange=[2400000,2800000],xrange=[21100000,21500000]
;      xyouts,21500000,2780000,Orbit, charsize=1.5,alignment=1

end;---------------------------------------------------------------------------------------


;|-----------------------------------------------------------------|
;|                                                                 |
;|    Subroutine IDM                                               |
;|                                                                 |
;|    extracts basic parameters from Idata.                        |
;|                                                                 |
;|-----------------------------------------------------------------|

pro idm

@var.inc


;Get the no. of icolumns in the HDF file

icolumns=idimsizes(1)
orbitnumber=fix(idata(4,0))
OrbitNum='IDM Orbit NO.'+ strcompress(fix(idata(4,0)))
orbitdate='Date :'+strcompress(idata(1,0))
date=round(idata(1,0))

catch,errorstat
if errorstat ne 0 then begin
	print,'!!!!!'
        print,'Error is ',errorstat
endif


;
;First extract Horinzontal and vertical SDS DATA which is stored in rows 46-109 but reduce 1 since IDL starts arrays with row number 0. 
;That is, we will extract 64 rows between 45 and 108.We'll get all icolumns.
;

sdsstart=45.0
sdsrows=64.0
sdsvel=EXTRAC(idata,sdsstart,0,sdsrows,icolumns)
isdsvel=fix(sdsvel)
iflag=round(100.0*abs(sdsvel-isdsvel))
tsdshvel=fltarr(sdsrows/2,icolumns)
tsdsvvel=fltarr(sdsrows/2,icolumns)
sdstimest=fltarr(icolumns)
sdstimeinc=fltarr(sdsrows)
sdstime=fltarr(sdsrows,icolumns)
tsdshtime=fltarr(sdsrows/2,icolumns)
tsdsvtime=fltarr(sdsrows/2,icolumns)
tsdsdate=fltarr(icolumns)
tsdsut=fltarr(icolumns)
tsdsilat=fltarr(icolumns)
tsdsmlt=fltarr(icolumns)
tsdsalt=fltarr(icolumns)
tsdsglat=fltarr(icolumns)
tsdsglon=fltarr(icolumns)
tsdshline_fit=fltarr(2,icolumns)
tsdsvline_fit=fltarr(2,icolumns)
tsdshvel_fit=fltarr(sdsrows/2,icolumns)
tsdsvvel_fit=fltarr(sdsrows/2,icolumns)
jcnt=0
max2=-1.e3
min2=+1.e3
tmax=1.
tmin=1.e10
;
;Big loop
;
For j=0,icolumns-1 do begin



	Flagsize=Size(where(iflag(*,j) eq 60))
	sdstimest(j)=idata(2,j)
	hcnt=0
	vcnt=0
;
;skip low invariant lats (above), need 32 60s
;

	If (idata(13,j) gt 45) and Flagsize(1) eq 32  then begin

		For i=0,sdsrows-1 do begin

			sdstimeinc(i)=i*4*31.25+75
			sdstime(i,j)=sdstimest(j)+sdstimeinc(i)

			If iflag(i,j) eq 60 then begin    ;HORIZONTAL SDS

				If i eq 0 then begin
					tsdsdate(jcnt)=idata(1,j)
				        tsdsut(jcnt)  =idata(2,j)
				        tsdsalt(jcnt) =idata(9,j)
				        tsdsglat(jcnt)=idata(10,j)
				        tsdsglon(jcnt)=idata(11,j)
				        tsdsmlt(jcnt) =idata(12,j)
				        tsdsilat(jcnt)=idata(13,j)
				Endif
				tsdshvel(hcnt,jcnt)=fix(sdsvel(i,j))


				tsdshtime(hcnt,jcnt)=sdstime(i,j)
      ;tsdshtime(hcnt,jcnt)=sdstime(i,j)/1000.
				hcnt=hcnt+1

			Endif else if iflag(i,j) eq 61 then begin ;vertical SDS
				tsdsvvel(vcnt,jcnt)=fix(sdsvel(i,j))
      ;tsdsvtime(vcnt,jcnt)=sdstime(i,j)/1000.
				tsdsvtime(vcnt,jcnt)=sdstime(i,j)
				vcnt=vcnt+1
				endif

				if (sdsvel(i,j) gt max2) then max2=sdsvel(i,j)
				if (sdsvel(i,j) lt min2) then min2=sdsvel(i,j)
                                if (sdstime(i,j) gt tmax) then tmax=sdstime(i,j)
                                if (sdstime(i,j) lt tmin) then tmin=sdstime(i,j)

			
			
 
		Endfor

		jcnt=jcnt+1

	Endif ;end of the big if 

Endfor ;end of the big loop
 
;tmin=sdstime(0,0)
;tmax=sdstime(sdsrows-1,icolumns-1)

sdsjcnt=jcnt

if jcnt eq 0 then return

;clean up the zeros

sdshvel=tsdshvel(*,0:jcnt-1)
sdsvvel=tsdsvvel(*,0:jcnt-1)
sdshtime=tsdshtime(*,0:jcnt-1)
sdsvtime=tsdsvtime(*,0:jcnt-1)
sdsdate=tsdsdate(0:jcnt-1)
sdsut=tsdsut(0:jcnt-1)
sdsmlt=tsdsmlt(0:jcnt-1)
sdsilat=tsdsilat(0:jcnt-1)
sdsalt=tsdsalt(0:jcnt-1)
sdsilat=tsdsilat(0:jcnt-1)
sdsmlt=tsdsmlt(0:jcnt-1)
sdsalt=tsdsalt(0:jcnt-1)
sdsglat=tsdsglat(0:jcnt-1)
sdsglon=tsdsglon(0:jcnt-1)



end;--------------------------------------------------------------------------------------


;|-----------------------------------------------------------------|
;|                                                                 |
;|    Subroutine WATS                                              |
;|                                                                 |
;|    extracts basic parameters from Wdata.                        |
;|                                                                 |
;|-----------------------------------------------------------------|


pro wats


@var.inc


; Now you can use Neu_Data and Columns in this code

neu_orbitnum=intarr(wcolumns)
nflag = intarr(wcolumns)
neu_vel=fltarr(wcolumns)
sc_vel=fltarr(wcolumns)
neu_flag_row = 2
neu_vel_row = 11
;sc_vel_row = 7

neu_orbitnum = fix(wdata(13,0:wcolumns-1))
neu_vel = wdata(neu_vel_row,0:wcolumns-1)
;sc_vel = wdata(sc_vel_row,0:wcolumns-1)
nflag = fix(wdata(neu_flag_row,0:wcolumns-1))/100
Orbitnumber=fix(idata(4,0))

; This next section does the same thing (keep it for future ref)

;neu_vel=EXTRAC(wdata, neu_vel_row, 0, 1, wcolumns)
;nflag=fix(EXTRAC(wdata, neu_flag_row, 0, 1, wcolumns)) /100
;neu_orbitnum=fix(EXTRAC(wdata, 13, 0, 1, wcolumns))

; Get Horizontal and Vertical Components of Neutral Velocity

tneu_hvel = FLTARR(wcolumns)
tneu_vvel = FLTARR(wcolumns)
tneu_schvel = FLTARR(wcolumns)
tneu_scvvel = FLTARR(wcolumns)
tneu_hut = FLTARR(wcolumns)
tneu_vut = FLTARR(wcolumns)
tneu_hglat = FLTARR(wcolumns)
tneu_hilat = FLTARR(wcolumns)
tneu_vglat = FLTARR(wcolumns)
tneu_vilat = FLTARR(wcolumns)

max1=-1.e3
min1=+1.e3
;tmax=1.
;tmin=1.e10

hcnt =1
vcnt =1

;a=""
for j=0,(wcolumns-1) do begin

;neutral_vel(j)=wdata(7,j)
;nflag(j)=fix((wdata(2,j))/100)



 if (neu_orbitnum(j) eq orbitnumber) then begin

   ; Set boundaries
        if neu_vel(j) ne 9999. then begin
		if neu_vel(j) gt max1 then max1=neu_vel(j)
		if neu_vel(j) lt min1 then min1=neu_vel(j)
		;if wdata(1,j) gt tmax then tmax=wdata(1,j)
		;if wdata(1,j) lt tmin then tmin=wdata(1,j)

     ;Horizontal neutral velocity

		if ((nflag(j) eq 3) or (nflag(j) eq 4)) and (neu_vel(j) gt -9000) then begin
        		tneu_hvel(hcnt)=neu_vel(j)
        		tneu_schvel(hcnt)=sc_vel(j)
        		tneu_hglat(hcnt)=wdata(15,j) 
        		tneu_hilat(hcnt)= wdata(20,j)
        		tneu_hut(hcnt)=wdata(1,j)    ;/3600000.
        		if abs(tneu_hglat(hcnt)) lt abs(tneu_hglat(hcnt-1)) then begin 
        			tneu_hvel(hcnt) =- tneu_hvel(hcnt) 
        			endif
        		;print,hcnt
      			hcnt=hcnt +1
      			endif 
     
     ;Vertical neutral velocity

     		if (nflag(j) eq 5) or (nflag(j) eq 6) then begin
        		tneu_vvel(vcnt)=neu_vel(j)
        		;tneu_scvvel(vcnt)=sc_vel(j)
       			tneu_vglat(vcnt)=wdata(15,j) 
		        tneu_vilat(vcnt)= wdata(20,j)
		        tneu_vut(vcnt)=wdata(1,j)   ;/3600000.
           		if abs(tneu_vglat(vcnt)) lt abs(tneu_vglat(vcnt-1)) then begin 
               			tneu_vvel(vcnt) =- tneu_vvel(vcnt) 
           			endif
        		;print,vcnt
        		vcnt=vcnt +1
     			endif 
		endif


 endif

endfor

;
;clearing up the zeros at the end
;sometimes the neutral data is not available so the if statement is necessary
;

if vcnt gt 1 then begin
  neu_hvel=tneu_hvel(1:hcnt-1)
  neu_schvel=tneu_schvel(1:hcnt-1)
  neu_hut=tneu_hut(1:hcnt-1)
  neu_hilat=tneu_hilat(1:hcnt-1)
  neu_hglat=tneu_hglat(1:hcnt-1)
endif else if vcnt eq 1 then begin 
  neu_hvel= tneu_hvel
  neu_schvel= tneu_hvel
  neu_hut = tneu_hut
  neu_hglat = tneu_hglat
  neu_hilat = tneu_hilat
  ;hcnt = wcolumns
endif 

if vcnt gt 1 then begin
  neu_vvel=tneu_vvel(1:vcnt-1)
  neu_scvvel=tneu_scvvel(1:vcnt-1)
  neu_vut=tneu_vut(1:vcnt-1)
  neu_vilat=tneu_vilat(1:vcnt-1)
  neu_vglat=tneu_vglat(1:vcnt-1)
endif else if vcnt eq 1 then begin
  neu_hvel= tneu_hvel
  neu_scvvel=tneu_scvvel
  neu_hut = tneu_hut
  neu_vilat=tneu_vilat
  neu_vglat=tneu_vglat
  ;hcnt = wcolumns
endif

end;--------------------------------------------------------------------------------------


;|-----------------------------------------------------------------|
;|                                                                 |
;|    Subroutine FILEPLOT                                          |
;|                                                                 |
;|    plots three graphs in a systematic fashion.                  |
;|                                                                 |
;|-----------------------------------------------------------------|

pro fileplot;

@var.inc

;set_plot,'PS'
;device,/inches,yoffset=10.5,ysize=7.0,xsize=10.0,/landscape,/color,filename='/usr33/users/klenzing/images/'+strtrim(name,2)+'.ps'

;!p.thick = 2
;!x.thick = 2
;!y.thick = 2

print,'Plotting file'

!P.Multi=[0,1,3]
!Y.Charsize=1.2
;!X.Charsize=1.2

;Plotting Neutral velocities

	!P.position=[0.1,0.7,0.9,0.9]
	plot,neu_hut,neu_hvel, xstyle=1,ystyle=1,symsize=.5,$
	title='Neutral vel,Ion vel and Ion Density plots for orbit '+strtrim(Orbitnumber,2),$
	ytitle ='Neutral Vel(m/sec)',xtitle='UT',Yrange=[min1,max1],xrange=[tmin,tmax],$
	xcharsize=0.01,psym=1
	;xyouts,21200000,80,OrbitNum, charsize=1.1,alignment=1
	;xyouts,21800000,400,'NVvel: ', charsize=1.0,alignment=1
	   ;OVERPLOT
	oplot,neu_vut,neu_vvel,psym=2,symsize=.5 

;Plotting Ion  velocities

	!P.position=[0.1,0.5,0.9,0.7]
	Plot,sdsvtime,sdsvvel,xstyle=1,ystyle = 1,psym=1,symsize=.5,ytitle='Sds Vvelocity ',$
	  xtitle='U.T.',yrange=[min2,max2],xrange=[tmin,tmax]
	;xyouts,21200000,50,IOrbitNUM, charsize=1.1,alignment=1
	;xyouts,21800000,150,'Vvel: ', charsize=1.0,alignment=1
	     ;OVERPLOT
	
	oplot,sdshtime,sdshvel,psym=2,symsize=.5

;Plotting Ion Densities (N)
 
	!P.position=[0.1,0.3,0.9,0.5]
	plot,time,N,psym=1,xcharsize=0.01,xrange=[tmin,tmax]
	 ;       xyouts,21200000,2700000,Orbit, charsize=1.1,alignment=1
	        ;xyouts,21300000,2750000,Orbit, charsize=1.5,alignment=1
	     ;   xyouts,time(43844),2.35e6,'UT', charsize=.75
	;	xyouts,time(43844),2.3e6,'ILAT', charsize=.75
	;	xyouts,time(43844),2.25e6, 'MLT', charsize=.75
	;	xyouts,time(43844),2.2e6, 'ALT', charsize=.75
	;	xyouts,time(43844),2.15e6,'GLAT', charsize=.75
	;	xyouts,time(43844),2.1e6,'GLON', charsize=.75
	
	;for j=1.,2.,.25 do begin
	;
	;xyouts,time(j*l_step*512.),2.35e6,dutx(j*l_step), charsize=.75,alignment=0.25
	;xyouts,time(j*l_step*512.),2.3e6,DILAT(j*l_step), charsize=.75,alignment=0.5
	;xyouts,time(j*l_step*512.),2.25e6,DMLT(j*l_step), charsize=.75,alignment=0.5
	;xyouts,time(j*l_step*512.),2.2e6,DALT(j*l_step), charsize=.75,alignment=0.5
	;xyouts,time(j*l_step*512.),2.15e6,DGLAT(j*l_step), charsize=.75,alignment=0.5
	;xyouts,time(j*l_step*512.),2.1e6,DGLON(j*l_step), charsize=.75,alignment=0.5
	;  
	;endfor

	
;device,/close
;set_plot,'x'
;!p.thick=0
;!x.thick=0
;!y.thick=0

end;------------------------------------------------------------------------------------

;|-----------------------------------------------------------------|
;|                                                                 |
;|   Begin Main Program                                            |
;|                                                                 |
;|-----------------------------------------------------------------|

@var.inc
maxrank=long(2)
ddimsizes=lonarr(maxrank)
idimsizes=lonarr(maxrank)

cyclefile=' '
number=' '

Print,'Enter Directory Number(01 - 14)'
read, number
cyclefile=strtrim(number,2)+'.txt'
openr,1,cyclefile
loop:error_status=0
;---------> Begins Main Loop
while not eof(1) do begin

	hdfread   ;-------> Call Subroutine HDFREAD
	findwats  ;-------> Call Subroutine FINDWATS

	catch, error_status
        if error_status eq -250 then goto,loop

	filename='/usr33/users/klenzing/data/wats/'+strtrim(year,2)+strtrim(date,2)+'_de2_wats_2s_v01.asc'
	print,filename
	openr, 2, filename

	asciiread ;-------> Call Subroutine ASCIIREAD

	CLOSE,2
	FREE_lun,2

	duct      ;-------> Call Subroutine DUCT
	idm       ;-------> Call Subroutine IDM

	if jcnt eq 0 then goto,loop

	wats      ;-------> Call Subroutine WATS

@psl.cmd
	fileplot  ;-------> Call Subroutine FILEPLOT
@xplt.cmd
endwhile
;---------> Ends Main Loop

close,1
free_lun,1




end;----------------------------------------------------------