;
; IDL routine plots the magnetic field components and the total field 
; strength stored in CDF files.
;
; 1/12/98 by SHC, Raytheon STX
;
;- functions and procedures -------------------------------------------
;
function fun_epoch,year,doy,tod
yr=year
if n_params() eq 0 then begin
   print,'function,epoch0:',year,doy,tod
   return,-1
endif
;ep0 is epoch at start of 20th cent.
;tod - time of doy in seconds
;doy - day of year starting at 1
;year - ie 1974 is 74
if(yr gt 1000) then yr=yr-1900
ep0 = 5.99582304d13
return,(yr*365d0 + fix((yr-1)/4)+ doy-1)*8.64d7 +tod*1000d0 + ep0 ;in millisec
end
;
;------------------------------------------------------------------------
; to transfer epoch0 time in msec to yr doy hr min sec
;
pro yrdhms,epoch0,yr,doy,hr,min,sec 		; in milliseconds
ep0 = 5.99582304d13 		;  1900 jan 1 00:00:00 UT
ep1=ep0+3.15576d12		; year of 2000
if(epoch0 lt ep0) or (epoch0 gt ep1) then begin
	print,'func-yrdhms: out of range'
	return
endif
depoch=0.001d0*(epoch0-ep0)
yr=long(depoch/31536000.d0)	; start at the year of 1900
doy=long(depoch/86400.d0)-yr*365-fix((yr-1)/4)+1
if(doy le 0) then begin
	yr=yr-1
	doy=365+doy
endif
days=yr*365+fix((yr-1)/4)+doy-1
hr=long(depoch/3600.d0-24.d0*days)
min=long(depoch/60.d0-1440.d0*days-60.d0*hr)
sec=depoch-86400.d0*days-hr*3600.d0-min*60.d0
return
end
;------------------------------------------------------------------------
pro ical,yr,doy,month,dom,idoy = idoy ;set idoy=1 to compute day of year
if( n_params() eq 0 ) then begin
        print,'positional param: yr,doy,month,dom
        print,'keyword: idoy = 0-doy to month and dom, 1-month dom to doy'
        return
endif
if( n_elements(idoy) eq 0 ) then idoy = 0
days=[31,28,31,30,31,30,31,31,30,31,30,31]
yr0=yr
if yr0 lt 100 then yr0=1900+yr0
if( ((yr0-1900) mod 4) eq 0)then begin
        days(1) = 29 ;leap year
        ;print,'year:',yr0,' is a leap year'
endif else begin
        days(1) = 28
endelse
if(idoy eq 1)then begin
        doy =dom
        for i= 0,month-2 do doy = doy + days(i)
        return
endif else begin
        dom = doy
        for month = 1,12 do begin
                if( dom le days(month-1)) then goto, jump
                dom = dom - days(month-1)
        endfor
endelse
print,' error in date conversion'
jump:return
end
;
;------------------------------------------------------------------------
;
function TTICKS, offset,tvalue,doy0,yr0,day,year
tvalue_buf=tvalue+offset ;in minutes
hour=fix(tvalue_buf/60.)
minute=fix(tvalue_buf-60*hour)
sec=fix((tvalue_buf-fix(tvalue_buf))*60.)
day=hour/24
hour=hour-day*24
day=day+doy0
year=yr0
Check_Again:
if(day gt 365) then begin
    if((year-fix(year/4)*4) eq 0) then begin
	day=day-366
	year=year+1
    endif else begin
	day=day-365
	year=year+1
    endelse
    goto, Check_Again
endif

if(sec ne 0) then $
	return, string(hour, minute, sec, format="(i2.2,':',i2.2,':',i2.2)") $
else return, string(hour, minute, format="(i2.2,':',i2.2)")
end
;
;-------------------------------------------------------------------------
;
pro get_position,luntra,time,epoch0,act_time,start,dt,xgsm,ygsm,zgsm,mlat,mlt
re=6371.;
epoch=0.d0
temp0=0.
temp1=0.
temp2=0.
temp3=0.
temp4=0.
temp5=0.
temp6=0.
epoch0=time*60000.d0+start
Read_again:
if(not EOF(luntra)) then begin
;       readu,21,epoch,r,mlat,mlt,invlat,xgsm,ygsm,zgsm
        readu,luntra,epoch,temp0,temp1,temp2,temp3,temp4,temp5,temp6
endif else begin
;	print,'end of the file 21'
	xgsm=999.0
	ygsm=999.0
	zgsm=999.0
	mlat=999.0
	mlt=999.0
	return
endelse
if(epoch ge epoch0-dt) and (epoch le epoch0+dt) then begin
	xgsm=temp4/re
	ygsm=temp5/re
	zgsm=temp6/re
	mlat=temp1
	mlt=temp2
endif else begin
if(epoch gt epoch0+dt) then begin
	xgsm=999.0
	ygsm=999.0
	zgsm=999.0
	mlat=999.0
	mlt=999.0
endif else goto, Read_again	
endelse
act_time=epoch
return
end
;-------------------------------------------------------------------------
;
; searching input files
;
pro findf,path,filename,k,timef=timef,timel=timel
    list=findfile(path)
    nf=-1
    k=-1
    if list(0) ne '' then begin
      if n_elements(timef) ne 0 then begin
	if(n_elements(timel) eq 0) then  begin
	  filename=strarr(1)
          for i=0,n_elements(list)-1 do begin
            parts=str_sep(list(i),'/')
            fname=str_sep(parts(n_elements(parts)-1),'_')
            time=fname(3)
            if(time eq timef) then begin
              filename=list(i)
              ;print,filename
	      k=0
              return
            endif
          endfor
	endif else begin
          reads,format='(i4,i3)',timef,yr_start,doy_start
          reads,format='(i4,i3)',timel,yr_end,doy_end
	; if((yr_end-yr_start) gt 1 ) then begin
	;	print,'The time interval is too long...returned to main'
	;	return
	; endif
	  if(((yr_start/4)*4-yr_start) eq 0) then leapyr=1 $
	  else leapyr=0
	  filename=strarr((yr_end-yr_start)*365+doy_end-doy_start+leapyr)
          for i=0,n_elements(list)-1 do begin
            parts=str_sep(list(i),'/')
            fname=str_sep(parts(n_elements(parts)-1),'_')
            time=fname(3)
            if(time ge timef) and (time le timel) and $
			(nf lt n_elements(filename)) then begin
              ;print,list(i)
	      nf=nf+1
              filename(nf)=list(i)
	    endif else begin
	      if(time gt timel) or (nf ge n_elements(filename)) then begin
	        k=nf
                return
	      endif
            endelse
          endfor
	endelse
      endif
    endif else begin
	print,' no files found in the directory:',path
	;filename=strarr(1)
	;read,'Type in the name including the full path:',filename
    endelse
end
;--------------------------------------------------------------------------
;
; see what's available
;
pro whatf,path
    list=findfile(path)
    if list(0) ne '' then begin
        time=strarr(n_elements(list))
        for i=0,n_elements(list)-1 do begin
          parts=str_sep(list(i),'/')
          fname=str_sep(parts(n_elements(parts)-1),'_')
          time(i)=fname(3)
        endfor
        ;print,format='(6(" ",i7.7))',time
    endif
end
;
;-------------------------------------------------------------------------
;-------------------------------------------------------------------------
; Main program
;-------------------------------------------------------------------------
;-------------------------------------------------------------------------
; Multiple panel plotting routines for IDL 
;
iframe=0
idevice=1
N_dimen=1500000 	; the dimension of the buffers containing data points
tframeleng=0.	; time frame length in minutes
yyddd=74100	; Hawkeye year of launch
infile=strarr(1)
file21=strarr(1)
ps_out=strarr(1)
ymax=fltarr(6,5)
ymin=fltarr(6,5)
;
;*********************************Input file request*****************
;
close,/all
ascii_binary=3
orbital_info=1  ; =0: no orbital information; =1 have...
yr_start=0
mon_start=0
day_start=0
doy_start=0
hr_start=0
min_start=0
sec_start=0
yr_end=0
mon_end=0
day_end=0
doy_end=0
hr_end=0
min_end=0
sec_end=0
;path_hk='/ncf/nssdc/pub/www/hawkeye/cdf/data/hk_h0*v02.cdf'
;openr,lun,'/ncf/nssdc/pub/www/hawkeye/idl/tplot_hk_B_cdf.feed',/get_lun
;openw,lun_write,'/ncf/nssdc/pub/www/hawkeye/scratch_images/status.asc',/get_lun
path_hk='/usr/local/httpd/htdocs/nssdc/hawkeye/cdf/data/hk_h0*v02.cdf'
openr,lun,'/usr/local/httpd/htdocs/nssdc/hawkeye/idl/tplot_hk_B_cdf.feed',/get_lun
openw,lun_write,'/usr/local/httpd/htdocs/nssdc/hawkeye/pub_scratch_images/status.asc',/get_lun
;
time_format=0
if(time_format eq 0) then begin
    readf,lun,yr_start,mon_start,day_start,hr_start,min_start,sec_start
    ical,yr_start,doy_start,mon_start,day_start,idoy=1
    readf,lun,yr_end,mon_end,day_end,hr_end,min_end,sec_end
    ical,yr_end,doy_end,mon_end,day_end,idoy=1
endif else begin
    readf,lun,yr_start,doy_start,hr_start,min_start,sec_start
    ical,yr_start,doy_start,mon_start,day_start
    readf,lun,yr_end,doy_end,hr_end,min_end,sec_end
    ical,yr_end,doy_end,mon_end,day_end
endelse

if(yr_start lt 100) then $
timef=string(format='("19",i2.2,i3.3)',yr_start,doy_start) $
else timef=string(format='(i4.4,i3.3)',yr_start,doy_start)
;
if(yr_end lt 100) then $
timel=string(format='("19",i2.2,i3.3)',yr_end,doy_end) $
else timel=string(format='(i4.4,i3.3)',yr_end,doy_end)
;print,timef,' ',timel
;
;Request1:
;    read,' Input file? ',infile
;    openr,lun,infile(0),/get_lun
;    read,' File type?(1 - ascii/Card; 2 - binary/Card; 3 - binary/Spher) ',$
;	ascii_binary
;
;if(ascii_binary eq 3) then $
;    read,' Show orbital information? (0 - no; 1- yes)', orbital_info $
;else orbital_info=0
;
;if(orbital_info eq 1) then begin
;    read,' Input file for orbital info ? ',file21
;    openr,luntra,file21(0),/get_lun
;    reads,infile(0),format='(6x,i5,i2,1x,i5,i2)',yyddd0,hr0,yyddd1,hr1
;    file21(0)=string(format='("hktra_",i5.5,i2.2,"_",i5.5,i2.2,".bin")',$
;	yyddd0,hr0,yyddd1,hr1)
;    openr,21,file21(0)
;    print,' The orbital information will be read from ',file21(0)
;endif
;
;
Request1:
    ;print,'Input file for Hawkeye magnetic field:'
    findf,path_hk,infile,index_mgf,timef=timef,timel=timel
    ;nfiles_mgf=n_elements(infile)
    nfiles_mgf=index_mgf+1
    ;print,'number of files:',nfiles_mgf
    printf,lun_write,index_mgf ; if index_mgf < 0 : no data found
    close,lun_write
    free_lun,lun_write
;
; exit if no data found
;
    if(index_mgf lt 0) then exit
;
;    lun_cdf_mgf=cdf_open(infile(0))
;    cdf_info=cdf_inquire(lun_cdf_mgf)
;    cdf_max_rec_mgf=cdf_info.maxrec
;    cdf_control,lun_cdf_mgf,variable='Epoch',get_var_info=x
;    cdf_max_rec_mgf=x.maxrec
    ;print,'Maximum number of records: ',cdf_max_rec_mgf
;    print,'Pad value: ',x.padvalue
    ;print,'Number of rVariables: ',cdf_info.nvars
    ;print,'Number of zVariables: ',cdf_info.nzvars
;
;   time resolution of the data in minute
;
    dt_mgf=1.9/60.0;
    dt_pos=3*dt_mgf;
;
;    cdf_varget1,lun_cdf_mgf,'Epoch',t0,rec_start=5
;    cdf_varget1,lun_cdf_mgf,'Epoch',t1,rec_start=6
;    dt_mgf=(t1-t0)/60000.d0


;read,' Plot in GSM : 0 or SM : 1 coordinate? ',index_mgf
index_mgf=0
readf,lun,index_mgf
;read,lun,' Time frame length (minutes)? ',tframeleng
;read,' Number of major-tick intervals (0 for automatic selection)? ',xticks
xticks=0 ; to be determined later in the program
;if(xticks gt 0) then read,' Number of minor-tick intervals? ',xminor
xminor=1
yon=1
;read,' Define your own y scales? (yes 0, or no 1)',yon
if(yon eq 0) then begin
        temp1=0.0 ;float
        temp2=0.0 ;float
        temp3=0 ; integer
        temp4=0 ; integer
	yi=5
	read,'Bx min, max, # y.ticks, # y.minor? ',temp1,temp2,temp3,temp4
	ymin(5,1)=temp1
	ymax(5,1)=temp2
	!y.ticks=temp3
	!y.minor=temp4
	read,'By min, max? ',temp1,temp2
	ymin(5,2)=temp1
	ymax(5,2)=temp2
	read,'Bz min, max? ',temp1,temp2
	ymin(5,3)=temp1
	ymax(5,3)=temp2
	read,'Bt min, max? ',temp1,temp2
	ymin(5,4)=temp1
	ymax(5,4)=temp2
endif else yi=2
;
;every_npts=1
;read,' Draw every how many points?',every_npts
;if (every_npts eq 0) or (every_npts gt N_dimen) then every_npts=1
;
thick=1
color=191
;
;read,' Curve thickness (>0)? ',thick
;read,' Output device (0 - postcript file ; 1 - X window ; 2 - gif file (z-buffer))? ',idevice
readf,lun,idevice
readf,lun,b_range
if(b_range gt 0) then begin
    highcut=b_range
    lowcut=-highcut
endif else begin
    highcut=60000.0
    lowcut=-60000.0
endelse
ps_out=strarr(1)
gif_out=strarr(1)
readf,lun,gif_out
;
if (idevice eq 0) then begin 
;*** for poscript ***************************************************
	;print,' '
;	read,' Output file name ? ',ps_out
	ps_out(0)='idl_out.ps'
;        print,'The graphic output file will be ',ps_out(0)
	set_plot,'PS'
	device,filename=ps_out(0),xoffset=2,xsize=17,yoffset=2.5,ysize=22,$
		/HELVETICA,/BOLD,FONT_SIZE=12
        charsize=2.25
        xticklen=0.1
        !p.charthick=4.
        !p.charsize=charsize
        !x.thick=4.
        !y.thick=4.
endif else begin 
    if(idevice eq 1) then begin
;
;*** for X-window ***************************************************
	set_plot,'X'
	device,cursor_standard=34    ;cross hair 34, hand 60,  
        charsize=1.5
        xticklen=0.1
        !p.background=255
        !p.color=0
        !p.charthick=1
        !p.charsize=charsize
        !x.thick=2.
        !y.thick=2.
;
;*** for z-buf ****
;
endif else begin
        set_plot,'Z'
        charsize=2
        xticklen=0.1
        !p.background=255
        !p.color=0
        !p.charthick=1
        !p.charsize=charsize
        !x.thick=2.
        !y.thick=2.

endelse
endelse
;*********************************************************************
;
;epoch1900=5.99582304d10 	; in seconds
DegtoRad=3.141592653589793/180.
start=0.d0
start0=0.d0
;
; Setting starting time of the plots. ****
;
toffset=hr_start*60.0+min_start
temp=toffset*60.
start0=fun_epoch(yr_start,doy_start,temp)
temp=hr_end*3600.+min_end*60.
end0=fun_epoch(yr_end,doy_end,temp)
;print,format='(" start0=",f16.1,"      end0=",f16.1)',start0,end0
;
tframeleng=float((end0-start0))/60000.
;
years=yr_start
days=doy_start
epoch0=0.
dt=57000.d0*0.5  ; 5.7 s is the time res. of orbital information 
;
npts_aframe=1200 ; maximum number of points in a frame
every_npts_mgf=fix(tframeleng/dt_mgf/npts_aframe+0.5)
if(every_npts_mgf le 0) then every_npts_mgf=1
;
; the dimension of the buffers containing data points
;
N_dimen=fix(npts_aframe*1.1)
dt_mgf=every_npts_mgf*dt_mgf; time resolution of field showing on the plot
every_npts_pos=3*every_npts_mgf
dt_pos=dt_pos*every_npts_pos
;
time_mgf=dblarr(N_dimen)
bgsm=fltarr(3,N_dimen)
bt=fltarr(N_dimen)
;
N_dimen_pos=N_dimen/3+1
time_pos=dblarr(N_dimen_pos)
posx=fltarr(N_dimen_pos)
posy=fltarr(N_dimen_pos)
posz=fltarr(N_dimen_pos)
posr=fltarr(N_dimen_pos)
;
; Plotting routine --------------------------------------------------
;
npanels=4 	;number of panels.
!p.multi=[0,1,npanels]
;!p.charsize=2
;!p.charthick=2.
;
; To define the characteristics of the time axis
;
if(tframeleng le 0) then tframeleng=30.0 ;in minutes
;
if(xticks le 0) then begin
    if(tframeleng ge 2880.0) then !x.ticks=fix(tframeleng/360.0) $
    else if(tframeleng ge 60.) then !x.ticks=fix(tframeleng/60.+0.1) $
    else !x.ticks=fix(tframeleng/5)
    if(!x.ticks eq 0) then !x.ticks=4
  checkxticks:
    if(!x.ticks lt 4) then begin
	!x.ticks=!x.ticks*2
	goto,checkxticks
    endif else begin
    if(!x.ticks gt 14) then begin
	!x.ticks=!x.ticks/2
	goto,checkxticks
    endif 
    endelse
    if(tframeleng ge 60.) then !x.minor=fix(tframeleng/!x.ticks/5) $
    else !x.minor=fix(tframeleng/!x.ticks)
    if(!x.minor eq 0) then !x.minor=1
  checkxminor:
    if(!x.minor gt 6) then begin
	!x.minor=!x.minor/2
	goto, checkxminor
    endif
    xticks=!x.ticks
    xminor=!x.minor
endif else begin
    !x.ticks=xticks
    !x.minor=xminor
endelse
!x.style=9
;
; To define the default settings for y axises
;
;yi=2
!y.style=1
!y.minor=1
if (yi eq 0) then yi=1
if(yi eq 1) then begin
	ymax(yi,1)=60.
	ymax(yi,2)=60.
	ymax(yi,3)=60.
	ymax(yi,4)=120.
	ymin(yi,1)=-60.
	ymin(yi,2)=-60.
	ymin(yi,3)=-60.
	ymin(yi,4)=0.
endif else begin
if(yi eq 2) then begin 
	ymax(yi,1)=40.
	ymax(yi,2)=40.
	ymax(yi,3)=40.
	ymax(yi,4)=80.
	ymin(yi,1)=-40.
	ymin(yi,2)=-40.
	ymin(yi,3)=-40.
	ymin(yi,4)=0.
endif else begin
if(yi eq 3) then begin
	ymax(yi,1)=80.
	ymax(yi,2)=80.
	ymax(yi,3)=80.
	ymax(yi,4)=160.
	ymin(yi,1)=-80.
	ymin(yi,2)=-80.
	ymin(yi,3)=-80.
	ymin(yi,4)=0.
endif else begin
if(yi eq 4) then begin
	ymax(yi,1)=120.
	ymax(yi,2)=120.
	ymax(yi,3)=120.
	ymax(yi,4)=240.
	ymin(yi,1)=-120.
	ymin(yi,2)=-120.
	ymin(yi,3)=-120.
	ymin(yi,4)=0.
endif
endelse
endelse
endelse

;
; To determine the x positions of panels.
;
xleft=0.15
xright=0.98
;
; To determine the y positions of panels.
;
ybottom=0.15  	; In normal coordinate [0,1].
ytop=0.95	; in normal coordinate.
ygap=0.05	; ratio of the height.
;
ypos=fltarr(npanels)
height=(ytop-ybottom)/npanels
ypos(0)=ybottom
for i=1,npanels-1 do ypos(i)=ypos(i-1)+height
ygap=ygap*height	; gaps in normal coordinate.
;
; Initialize
;
;
time_label=strarr(30)
xvalue=fltarr(30)
pos_label=fltarr(2)
zeros=fltarr(2)
nsixty=fltarr(2)
vertical=fltarr(2)
adays=fltarr(30)
at=dblarr(2)
at_ver=fltarr(2)
pos_ut=fltarr(2)
nxlabel=4
gsm_buf=fltarr(nxlabel)
gsm=fltarr(nxlabel,30)
xlabel_head=strarr(nxlabel+1)
;
zeros(0)=0.
zeros(1)=0.
nsixty(0)=-100.
nsixty(1)=-100.
vertical(0)=ytop-ygap
vertical(1)=ybottom
if(index_mgf eq 0) then begin
    gsm_sm='(GSM)'
    xlabel_head(0)='Xgsm(Re)'
    xlabel_head(1)='Ygsm'
    xlabel_head(2)='Zgsm'
    xlabel_head(3)='R'
    ;xlabel_head(4)='MLT (hrs)'
endif else begin
    gsm_sm='(SM)'
    xlabel_head(0)='Xsm(Re)'
    xlabel_head(1)='Ysm'
    xlabel_head(2)='Zsm'
    xlabel_head(3)='R'
    ;xlabel_head(4)='MLT (hrs)'
endelse
xlabel_head(nxlabel)='UT(h:m)'
;
;
k=0
i=0
it1=0L
it2=0L
iskip=0
bx_buf=0.0
by_buf=0.0
bz_buf=0.0
bt_buf=0.0
colat_buf=0.0
long_buf=0.0
t1=0.d0
t2=0.d0
act_time=0.d0
time_buf=0.d0
time_buf2=0.d0
;
;if(idevice eq 1) then begin
	;print,' '
	;print,'*** Click in the left margin to go to the next page'
	;print,'*** Click in the right margin to quite'
	;print,'*** Click in the panels to select time intervals'
	;print,' '
;endif
;
npoints_mgf=0L
npoints_pos=0L
;
; read magnetic field cdf
;
for i=0,nfiles_mgf-1 do begin
	;print,'file :',i,infile(i)
	lun_cdf_mgf=cdf_open(infile(i))
	cdf_info=cdf_inquire(lun_cdf_mgf)
	cdf_control,lun_cdf_mgf,variable='Epoch',get_var_info=x
	cdf_max_rec_mgf=x.maxrec
	;
	cdf_search_epoch,lun_cdf_mgf,start0,cdf_rec_num_mgf,time_buf_mgf,$
                cdf_max_rec=cdf_max_rec_mgf
	cdf_search_epoch,lun_cdf_mgf,end0,end_rec_num_mgf,act_epoch_mgf,$
                cdf_max_rec=cdf_max_rec_mgf
	rec_count_mgf=end_rec_num_mgf-cdf_rec_num_mgf+1L
	rec_count_mgf=long(rec_count_mgf/every_npts_mgf)+1
	rec_count_pos=rec_count_mgf/3+1

	;print,'reading data...'
        ;print,'start rec:',cdf_rec_num_mgf
	;print,'number of rec:',rec_count_mgf,rec_count_pos
	;print,'every npoints:',every_npts_mgf

        cdf_varget,lun_cdf_mgf,'Epoch',time_mgf_buf,count=1,$
            rec_start=cdf_rec_num_mgf,rec_interval=every_npts_mgf,$
            rec_count=rec_count_mgf
        cdf_varget,lun_cdf_mgf,'B_total',bt_buf,count=1,$
            rec_start=cdf_rec_num_mgf,rec_interval=every_npts_mgf,$
	    rec_count=rec_count_mgf
        cdf_varget,lun_cdf_mgf,'Epoch',time_pos_buf,count=1,$
            rec_start=cdf_rec_num_mgf,rec_interval=every_npts_pos,$
            rec_count=rec_count_pos
        cdf_varget,lun_cdf_mgf,'r',r_buf,count=1,offset=0,$
            rec_start=cdf_rec_num_mgf,rec_interval=every_npts_pos,$
            rec_count=rec_count_pos
        if(index_mgf eq 0) then begin
            cdf_varget,lun_cdf_mgf,'B_gsm',bgsm_buf,count=3,offset=0,$
                rec_start=cdf_rec_num_mgf,rec_interval=every_npts_mgf,$
                rec_count=rec_count_mgf
            cdf_varget,lun_cdf_mgf,'Pos_gsm',pgsm_buf,count=3,offset=0,$
                rec_start=cdf_rec_num_mgf,rec_interval=every_npts_pos,$
                rec_count=rec_count_pos
        endif else begin
            cdf_varget,lun_cdf_mgf,'B_sm',bgsm_buf,count=3,offset=0,$
                rec_start=cdf_rec_num_mgf,rec_interval=every_npts_mgf,$
                rec_count=rec_count_mgf
            cdf_varget,lun_cdf_mgf,'Pos_sm',pgsm_buf,count=3,offset=0,$
                rec_start=cdf_rec_num_mgf,rec_interval=every_npts_pos,$
                rec_count=rec_count_pos
        endelse

	jjj=0
        for jj=0,rec_count_mgf-1 do begin
		; convert to minutes
                time_mgf(npoints_mgf)=(time_mgf_buf(jj)-start0)*0.001/60.
		;print,N_dimen,npoints_mgf,time_mgf(npoints_mgf)
		bgsm(0,npoints_mgf)=bgsm_buf(0,jj)
		bgsm(1,npoints_mgf)=bgsm_buf(1,jj)
		bgsm(2,npoints_mgf)=bgsm_buf(2,jj)
		bt(npoints_mgf)=bt_buf(jj)
		if((jj mod 3) eq 0) then begin
                	time_pos(npoints_pos)=$
				(time_pos_buf(jjj)-start0)*0.001/60.
			posx(npoints_pos)=pgsm_buf(0,jjj)/6378.0
			posy(npoints_pos)=pgsm_buf(1,jjj)/6378.0
			posz(npoints_pos)=pgsm_buf(2,jjj)/6378.0
			posr(npoints_pos)=r_buf(jjj)/6378.0
			jjj=jjj+1
        		npoints_pos=npoints_pos+1
		endif
		npoints_mgf=npoints_mgf+1

        endfor ; loop for points in each files
 	cdf_close,lun_cdf_mgf
endfor ; loop for files

;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;
; Main Plotting Loop
;
;ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
;
while (time_buf lt end0) do begin
 	if(idevice eq 1)  then window, window_no,xsize=620,ysize=600.
;	if(cdf_rec_num_mgf ge (cdf_max_rec_mgf-1L)) then goto, TheEnd
;
	t1=start
	t2=t1+tframeleng
	t2_buf=t2*60000.d0+start0

        ;if(time_buf_mgf gt t2_buf) and (npoints_mgf gt 1) then $
        ;        npoints_mgf=npoints_mgf-1L

        ;print,'dt_mgf=',dt_mgf,'  dt_pos=',dt_pos

	it1=0L
	it2=npoints_mgf-1L
	if(it1 ge it2) then begin
		iskip=iskip+1
		goto, J102
	endif

	if(iskip gt 0) then begin
	    ;print,'*** Skipped ',iskip,$
	;	'  frames containing no or only one data point ***'
	    iskip=0
	endif

	i=i+1
	;print,' Frame ',i,'   No. of data pts: ',npoints_mgf
	at(0)=t1
	at(1)=t2
;
;	calculate the position of time ticks
;
	dtframeleng=tframeleng/!x.ticks
	xvalue(0)=0.0
	xvalue(!x.ticks)=tframeleng
	for kk=1,!x.ticks-1 do begin
	    xvalue(kk)=kk*dtframeleng
	endfor
;
	axis_scale,bgsm(1,it1:it2),npoints=npoints_mgf,yrange,accu=1,/linear,$
                upcutoff=highcut,lowcutoff=lowcut
	;nsixty(0)=ymin(yi,2)
	nsixty(0)=yrange(0)
	nsixty(1)=nsixty(0)
	plot,at,nsixty,$
		;xrange=[t1,t2],xticklen=0.03,xtick_get=xvalue,$
		xrange=[t1,t2],xticklen=0.03,$
		xtickname=replicate(' ',30),xticks=xticks,xminor=xminor,$
		ytitle='By',yticklen=0.008,$;yrange=[ymin(yi,2),ymax(yi,2)],$
		yrange=yrange,ystyle=1,$
		position=[xleft,ypos(2),xright,ypos(3)-ygap]
;
;	if(idevice eq 0) or (idevice eq 1) then begin
;	    t1_shade=0.d0
;	    t2_shade=0.d0
;	    shade=fltarr(2,4)
;	    openr,25,'hk_tint74184_shade.asc'
;	    while(not EOF(25)) do begin
;		readf,25,t1_shade,t2_shade,temp1
;		t1_shade=(t1_shade-start0)/60000.
;		t2_shade=(t2_shade-start0)/60000.
;		at_ver=convert_coord(t1_shade,0.,/to_normal)
;		shade(0,0)=at_ver(0)
;		shade(1,0)=ybottom
;		shade(0,1)=at_ver(0)
;		shade(1,1)=ytop-ygap
;		at_ver=convert_coord(t2_shade,0.,/to_normal)
;		shade(0,2)=at_ver(0)
;		shade(1,2)=ytop-ygap
;		shade(0,3)=at_ver(0)
;		shade(1,3)=ybottom
;		polyfill,shade(0:1,0:3),col=200,/normal
;	    endwhile
;	    close,25
;	endif

	oplot_gaps,time_mgf(it1:it2),bgsm(1,it1:it2),thick=thick,$
                highcut=60000.,lowcut=-60000.,dx=dt_mgf
	oplot,at,zeros
;	axis,t1,ymin(yi,2),xaxis=0,xtickname=replicate(' ',30),/data
;
	axis_scale,bgsm(2,it1:it2),npoints=npoints_mgf,yrange,accu=1,/linear,$
                upcutoff=highcut,lowcutoff=lowcut
	plot_gaps,time_mgf(it1:it2),bgsm(2,it1:it2),$
                highcut=60000.,lowcut=-60000.,dx=dt_mgf,xticks=xticks,$
		xminor=xminor,$
		xrange=[t1,t2],xticklen=0.03,xtickname=replicate(' ',30),$
		ytitle='Bz',yticklen=0.008,$;yrange=[ymin(yi,3),ymax(yi,3)],$
		yrange=yrange,ystyle=1,$
		position=[xleft,ypos(1),xright,ypos(2)-ygap],thick=thick
	oplot,at,zeros
	;
	;The bottom panel.
	;
        kk_buf=0
	itemp2=0
        gsm_buf(0)=999.
        gsm_buf(1)=999.
        gsm_buf(2)=999.
        gsm_buf(3)=999.
	for jj=0,!x.ticks do begin
		time_label(jj)=TTICKS(toffset,xvalue(jj),days,years,$
			itemp1,itemp2)
		adays(jj)=itemp1
                if(orbital_info eq 1) then begin
                    for kk=kk_buf,npoints_pos-1 do begin
                        if(time_pos(kk) ge (xvalue(jj)-dt_pos*0.5)) and $
                          (time_pos(kk) lt (xvalue(jj)+dt_pos*0.5)) then begin
			    temp1=posx(kk)
			    temp2=posy(kk)
			    temp3=posz(kk)
			    temp4=posr(kk)
			    ;print,jj,kk,temp1,temp2,temp3,temp4
			    kk_buf=kk+1
			    goto, label1
                        endif else begin
                            if(time_pos(kk) ge (xvalue(jj)+dt_pos)) then begin
				temp1=999.
				temp2=999.
				temp3=999.
				temp4=999.
				kk_buf=kk+1
				goto, label1
                            endif
                        endelse
                    endfor
            label1: if(jj eq 0) and (temp1 gt 990.) then begin
                        temp1=gsm_buf(0)
                        temp2=gsm_buf(1)
                        temp3=gsm_buf(2)
                        temp4=gsm_buf(3)
                    endif else begin
                        if(jj eq !x.ticks) then begin
                            gsm_buf(0)=temp1
                            gsm_buf(1)=temp2
                            gsm_buf(2)=temp3
                            gsm_buf(3)=temp4
                        endif
                    endelse
		    gsm(0,jj)=temp1
		    gsm(1,jj)=temp2
		    gsm(2,jj)=temp3
		    gsm(3,jj)=temp4
                endif
	endfor
	;
	axis_scale,bt(it1:it2),npoints=npoints_mgf,yrange,accu=1,/linear,$
		upcutoff=highcut,lowcutoff=-0.5
	plot_gaps,time_mgf(it1:it2),bt(it1:it2),$
		xticks=xticks,xminor=xminor,xticklen=0.05,$
                highcut=60000.,lowcut=-0.1,dx=dt_mgf,xtitle=' ',$
		xrange=[t1,t2],xtick_get=xvalue,xtickname=time_label,$
		ytitle='|B|',yticklen=0.008,$;yrange=[ymin(yi,4),ymax(yi,4)],$
		yrange=yrange,ystyle=1,$
		position=[xleft,ypos(0),xright,ypos(1)-ygap],thick=thick
	;
	; X - Labels below the primary time lablels  
	;
	xlabel_shift=0.53
;	xlabel_shift=0.65
	xlabel_gap=0.18
	for kk=0,nxlabel-1 do begin 
	    factor=xlabel_shift-kk*xlabel_gap
	    for jj=!x.ticks,0,-1 do begin
		;pos_label=convert_coord(xvalue(jj),ymin(yi,4),/to_device)
		pos_label=convert_coord(xvalue(jj),yrange(0),/to_device)
		pos_label(1)=factor*pos_label(1)
		if (abs(gsm(kk,jj)) lt 990.) then begin
		    xlabel=string(gsm(kk,jj),format='(f5.2," ")')
		    xyouts,pos_label(0),pos_label(1),xlabel,$
				alignment=0.5,/device,charsize=charsize*0.5
		endif
	    endfor
	    if(kk eq 0) then begin
		pos_ut(0)=pos_label(0)*0.05
		pos_ut(1)=pos_label(1)*(xlabel_shift+1.5*xlabel_gap)/factor
	    endif
	    xyouts,pos_ut(0),pos_label(1),xlabel_head(kk),/device,$
			charsize=charsize*0.5
	endfor
	xyouts,pos_ut(0),pos_ut(1),xlabel_head(nxlabel),/device,$
			charsize=charsize*0.5
	;
	;The top panel.
	;
	if(adays(!x.ticks) ne adays(0)) then begin
	   title=string(format=$
	   '("Hawkeye Magnetometer ",A5 ," Day ",i3,"-",i3,", ",i4)',$
		gsm_sm,adays(0),adays(!x.ticks),itemp2)
	endif else begin 
	   title=string(format=$
	   '("Hawkeye Magnetometer ",A5," Day ",i4,", ",i4)',$
		gsm_sm,adays(0),itemp2)
	endelse

	axis_scale,bgsm(0,it1:it2),npoints=npoints_mgf,yrange,accu=1,/linear,$
		upcutoff=highcut,lowcutoff=lowcut
	plot_gaps,time_mgf(it1:it2),bgsm(0,it1:it2),$
                highcut=60000.,lowcut=-60000.,dx=dt_mgf,title=title,$
		xstyle=1,xrange=[t1,t2],xticklen=0.03,xticks=xticks,$
		xminor=xminor,xtickname=replicate(' ',30),$
		ytitle='Bx (nT)',yticklen=0.008,$;yrange=[ymin(yi,1),ymax(yi,1)],$
		yrange=yrange,ystyle=1,$
		position=[xleft,ypos(3),xright,ytop-ygap],thick=thick
	oplot,at,zeros
;
;	select and output intervals
;
	if (idevice eq 1) then begin
J101:		cursor,xcursor,ycursor,/data,/down
		if(xcursor gt t2) then goto,TheEnd
		while(xcursor ge t1) and (xcursor le t2) do begin
			at_ver=convert_coord(xcursor,0.,/to_normal)
			at_ver(1)=at_ver(0)
			plots,at_ver,vertical,/normal
			start_int=start0+xcursor*60000.
			;print,format=$
			    ;'("starting time = ",f9.3,"   epoch= ",f16.1)',$
			    xcursor,start_int
			cursor,xcursor,ycursor,/data,/down
			if(xcursor lt t1) or (xcursor gt t2) then goto, J101
			at_ver=convert_coord(xcursor,0.,/to_normal)
			at_ver(1)=at_ver(0)
			plots,at_ver,vertical,/normal
			end_int=start0+xcursor*60000.
			;print,format=$
			    ;'("ending time   = ",f9.3,"   epoch= ",f16.1)',$
			    xcursor,end_int
			duration=(end_int-start_int)*0.001
			;print,'duration      = ',duration,' s'
			;print,'--------------------------------------------'
			if(duration lt 0.) then begin
				buf=start_int
				start_int=end_int
				end_int=buf
			endif
			kpoints=0L
			iframe=iframe+1
			if(iframe gt 1) then $
				;printf,24,1.e34,1.e34,1.e34,1.e34,1.e34
			ll=0L
			for ll=it1,it2 do begin
			    buf=time(ll)*60000.+start0
			    if(buf ge start_int) and $
					(buf le end_int) then begin 
				kpoints=kpoints+1L
				;printf,24,$
				    ;format='(f16.1,1x,f8.3,f8.3,f8.3,f8.3)',$
					buf,bx(ll),by(ll),bz(ll),bt(ll)
			    endif else begin
				if(buf gt end_int) then goto, TheNext
			    endelse
			endfor
TheNext:
			    ;printf,23,format='(f16.1,1x,f16.1,1x,i8)',$
					start_int,end_int,kpoints
			cursor,xcursor,ycursor,/data,/down
		endwhile
	endif
J102:
	start=t2
	if((start*60000.+start0) ge end0) or (time_buf ge end0) $
		then goto, TheEnd

	for jj=0,nxlabel-1 do begin
	    for kk=0,!x.ticks do begin
		gsm(jj,kk)=999.0
	    endfor
	endfor
        window_no=window_no+1
        window_no=window_no-window_no/2*2
endwhile

TheEnd:

if(idevice eq 0) then begin
    device,/close
    ;print,'The graphic output file will be ',ps_out(0)
    set_plot,'X'
endif else begin
    bytemap=tvrd()
    ;print,'creating GIF output: idl.gif...'
    write_gif,gif_out(0),bytemap
    ;print,'Done.'
endelse

close,/all
end