;+
; NAME:
;	IMF_PLOT
; PURPOSE:
;	Plot interplanetary magnetic field data from CDF files of IMP 8
;	vector magnetometer data.
; CATEGORY:
;	
; CALLING SEQUENCE:
;	IMF_PLOT, YEAR, MONTH, DAY, START, FINISH [, NAME] $
;	       [, /POLAR] [, /NONORM] [, /PS] $
;	       [, OUT = output_filename ] [, _EXTRA = e]
; INPUTS:
;	YEAR,
;	MONTH,
;	DAY:	The date of the observations to be displayed.
;	START,
;	FINISH: The first and last UT times (in hours) to be displayed.
; OPTIONAL INPUTS:
;	NAME:	The name of the CDF file from which to read the data.
; KEYWORD PARAMETERS:
;       /POLAR: Set this keyword if you want a plot of the magnetic
;               field magnitude, longitude, and latitude.  The
;               default is to plot 3 cartesian components and
;		magnitude.
;       /NONORM: Set this keyword if you want to plot the field
;               components or other quantities on fixed-range
;               scales.  The default is to adjust the axis limits
;		to show the maximum detail.
;	/PS:	Set this keyword if you want to produce a PostScript file 
;		ofthe plot, and output the file to queue SYS$PS. 
;	OUT:	Specifies the filename of a text file to which the plotted
;		data are to be written for later use.
;       _EXTRA: Any keywords applicable to the PLOT procedure may
;               be supplied as keywords to IMF_PLOT.  They will be
;		supplied to the PLOT commands in the procedure.
; OUTPUTS:
;	None.
; OPTIONAL OUTPUTS:
;	None.
; COMMON BLOCKS:
;	None.
; SIDE EFFECTS:
;       A CDF file is read.  A plot is generated on the specified
;	graphics device.
; RESTRICTIONS:
;	No more than 24 hours of data may be plotted at one time.
;	The data are read and plotted in GSM coordinates.
; PROCEDURE:
;	Straightforward.
; EXAMPLE:
;	
; SEE ALSO:
;	
; MODIFICATION HISTORY:
; 	Written by:	Your name here, Date.
; Modified 94/06/21 by D P Steele to read data for a given day from all 
; possible files, if the user doesn't specify a single file to read.
;-
PRO imf_plot,yr,month,day,start,finish,name $
            ,polar=polar,nonorm=nonorm,ps=ps,out=output_filename,_EXTRA=e

IF N_PARAMS() EQ 0 THEN BEGIN
    doc_library,'imf_plot'
    RETURN
    ENDIF

; First make sure the year is in the correct format

year=yr
IF yr LT 100 THEN IF yr LT 50 THEN year=year+2000 ELSE year=year+1900

; Make sure the time limits make sense

IF finish LT start then swap,start,finish

; Calculate the Julian date for the desired day, for comparison with the
; start dates of the available files 

jdobs=julday(month,day,year)
doy=jdobs-julday(1,1,year)+1

IF N_PARAMS() EQ 6 THEN name=[name] ELSE BEGIN

; If the user hasn't specified a file to read, check all the available files 
; to see which ones might have data for the desired day.

    cdf=FINDFILE('I8_K0_MAG_199%%%%%_V0%.CDF',COUNT=ncdf)
    jddiff=INTARR(ncdf)
    FOR i=0,ncdf-1 DO BEGIN
        p=STRPOS(cdf(i),'199')
        cdfy=FIX(STRMID(cdf(i),p,4))
        cdfm=FIX(STRMID(cdf(i),p+4,2))
        cdfd=FIX(STRMID(cdf(i),p+6,2))
        jddiff(i)=jdobs-julday(cdfm,cdfd,cdfy)
        ENDFOR

; Restrict the search to files starting on or before the desired day.

    ok=WHERE(jddiff GE 0,ncdf)
    IF ncdf LT 1 THEN BEGIN
        MESSAGE,'No appropriate data files found; quitting',/INFORMATIONAL
        RETURN
        ENDIF

    jddiff=jddiff(ok)
    cdf=cdf(ok)
    n=MIN(jddiff,ix)

; Create an array of the file name(s) that might have useful data.

    IF n GT 0 THEN name=[cdf(ix)] ELSE name=cdf(ix-1:ix)
    PRINT,'Checking the following file(s) for data:'
    FOR k=0,N_ELEMENTS(name)-1 DO PRINT,k+1,'  ',name(k)
    ENDELSE

; Now check all the files selected.

FOR j=0,N_ELEMENTS(name)-1 DO BEGIN
    PRINT,name(j)
    id=cdf_open(name(j))
    cdf_info=cdf_inquire(id)
    len=cdf_info.maxrec
    t=LONARR(3,len)
    b=FLTARR(3,len)
    gsmpos=b

; Read all of the time and GSM magnetic field information from the file.

    FOR i=0,len-1 DO BEGIN
        cdf_varget,id,'Time_PB5',time,rec_start=i
        t(*,i)=time
	IF time(1) EQ doy THEN BEGIN
            cdf_varget,id,'B_GSM_c',imf,rec_start=i
            b(*,i)=imf
            cdf_varget,id,'SC_pos_sm',pos,rec_start=i
            gsmpos(*,i)=pos
	    ENDIF
	IF ((i+1) MOD 100) EQ 0 THEN PRINT,i+1,FORMAT='(I5,$)'
	IF ((i+1) MOD 1600) EQ 0 THEN PRINT,FORMAT='()'
        ENDFOR

; Close the CDF file before going any further.

    cdf_close,id
    PRINT,FORMAT='()'
    PRINT,len,' records read',FORMAT='(I0,A)'

; Now restrict attention to the desired day.

    ok=WHERE(t(1,*) EQ doy,nok)
    IF nok EQ 0 $ 
	THEN MESSAGE,'No data found for day '+STRTRIM(doy,2),/INFORMATIONAL $
        ELSE BEGIN
        t=t(2,ok)

; Reformat B to a 1-d array so array concatenation will do what we want.

        b=REFORM(b(*,ok),3*N_ELEMENTS(ok))
	gsmpos=REFORM(gsmpos(*,ok),3*N_ELEMENTS(ok))
        PRINT,nok,' records from day '+STRTRIM(doy,2),FORMAT='(I0,A)'

; Now incorporate the times and magnetic field data found into the arrays 
; giving all the data for the day.

        IF N_ELEMENTS(allt) EQ 0 THEN BEGIN
            allt=REFORM(t)
            allb=b
	    allpos=gsmpos
            ENDIF ELSE BEGIN
            allt=[allt,REFORM(t)]
            allb=[allb,b]
	    allpos=[allpos,gsmpos]
            ENDELSE

        ENDELSE

    ENDFOR

IF N_ELEMENTS(allt) EQ 0 THEN BEGIN
    MESSAGE,'No data available for day '+STRTRIM(doy,2),/INFORMATIONAL
    RETURN
    ENDIF ELSE BEGIN

; Now put ALLB back in the original format.

    allb=REFORM(allb,3,N_ELEMENTS(allb)/3)
    allpos=REFORM(allpos,3,N_ELEMENTS(allpos)/3)

; Sort the data by time.

    s=SORT(allt)
    t=allt(s)
    b=allb(*,s)
    pos=allpos(*,s)/6378.15		; convert to Earth radii

; Convert the time into UT seconds from midnight

    uts=t/1000.D0

; Select the data to be displayed.

    futs=3600.*start
    luts=3600.*finish
    wdisp=WHERE((futs LE uts) AND (uts LE luts),nwdisp)
    IF nwdisp GT 0 THEN BEGIN
        t=t(wdisp)
        uts=t/1000.D0
        b=b(*,wdisp)
        pos=pos(*,wdisp)
        ENDIF ELSE MESSAGE,'No data for specified interval!'

; If a polar plot is desired, generate the needed angles.

    bmag=TRANSPOSE(SQRT(TOTAL(b^2,1)))
    IF KEYWORD_SET(polar) THEN BEGIN
        blong=(!RADEG*ATAN(b(1,*),b(0,*))+360.) MOD 360.
        blat=90.-!RADEG*ACOS(b(2,*)/bmag)
        ENDIF

; Write the data to a file if desired.

    IF N_ELEMENTS(output_filename) NE 0 THEN BEGIN
        IF KEYWORD_SET(polar) THEN bout=[bmag,blat,blong] ELSE bout=b
        HELP,uts,bout
        STOP
        OPENW,outunit,output_filename,/GET_LUN
        FOR iout=0L,N_ELEMENTS(uts)-1 DO $
            PRINTF,outunit,uts(iout),bout(*,iout) $
                  ,FORMAT='(F7.1,3(2X,F5.1))'
        FREE_LUN,outunit
        ENDIF

; Set up for plotting.

    !Y.MARGIN=[7,2]
    !P.FONT=-1
    nhrs=finish-start
    CASE 1 OF
                   nhrs LE 3: interv=0.5
 (3 LT nhrs) AND (nhrs LE 6): interv=1.0
(6 LT nhrs) AND (nhrs LE 12): interv=2.0
                (12 LT nhrs): interv=4.0
                        ELSE:
        ENDCASE

    IF interv LE 1.0 THEN xminor=6 ELSE xminor=8
    labuts=makex(interv*CEIL(start/interv) $
                ,interv*FLOOR(finish/interv) $
                ,interv)*3600.
    labgsmx=STRTRIM(STRING(interpol(REFORM(pos(0,*)),uts,labuts),FORMAT='(F6.2)'),2)
    labgsmy=STRTRIM(STRING(interpol(REFORM(pos(1,*)),uts,labuts),FORMAT='(F6.2)'),2)
    labgsmz=STRTRIM(STRING(interpol(REFORM(pos(2,*)),uts,labuts),FORMAT='(F6.2)'),2)
    labhrmn=STRMID(strsec(labuts),0,5)
    xticks=N_ELEMENTS(labuts)-1
    xlab=REPLICATE(' ',xticks+1)

    IF KEYWORD_SET(polar) THEN BEGIN
        ytit=['!17B!x Longitude (deg)','!17B!x Latitude (deg)']
        ylim=[[0.,360.],[-90.,90.]]
        yticks=[4,4]
        ytickv=[[0,90,180,270,360],[-90,-45,0,45,90]]
        ytickname=[['','90','180','270',''],['','-45','0','45','']]
        yminor=[3,3]
        unit=' deg'
        ENDIF ELSE BEGIN
        ytit=['B!ix!n (nt)','B!iy!n (nT)','B!iz!n (nT)']
        ylim=REBIN([-20.,20.],2,3)
        yticks=[4,4,4]
        ytickv=REBIN([-20,-10,0,10,20],5,3)
        ytickname=[['','-10','0','10',''],['','-10','0','10',''] $
	          ,['','-10','0','10','']]
        yminor=[5,5,5]
        unit=' nT'
        ENDELSE
    np=N_ELEMENTS(ytit)
    ptitle='IMP-8 IMF Data for '+ymd2date(year,month,day,format="d$ n$ y$")
    ptitle=ptitle+' (Day '+STRTRIM(doy,2)+')'
    !X.TICKLEN=0.08

; Open a PostScript file, if necessary.

    IF KEYWORD_SET(ps) THEN BEGIN
        psname=STRING(year MOD 100,month,day $
                     ,FORMAT='(3I2.2,"_IMF.PS")')
        IF !VERSION.OS EQ 'vms' THEN open_ps $
                                ELSE psopen,psname,/land
        ENDIF
	
; Plot the magnetic field, starting with the X, Y, Z components.

    multiplot,[0,1,np+1,0,0],/init
    FOR i=0,np-1 DO BEGIN
        multiplot

    ; If polar plot desired, plot first longitude and then latitude.

        IF KEYWORD_SET(polar) THEN BEGIN
	    CASE i of
	     0: data=blong
	     1: data=blat
          ELSE:
	        ENDCASE
	    limform='(F6.1,A)'
	    ENDIF ELSE BEGIN
	    data=b(i,*)
	    limform='(F6.2,A)'
	    ENDELSE

    ; If fixed axis scales desired, set them up, else check the data.

        PRINT,MIN(data,MAX=datamax),datamax
        IF KEYWORD_SET(nonorm) THEN BEGIN
            minb=ylim(0,i)
            maxb=ylim(1,i)
            ENDIF ELSE BEGIN
	    maxb=MAX(data,MIN=minb)
	    ENDELSE
	maxp=CEIL(maxb)
	minp=FLOOR(minb)
        IF (NOT (KEYWORD_SET(polar) OR KEYWORD_SET(nonorm))) $
           AND (i EQ 2) THEN BEGIN
            maxp=maxp>0		; ensure Bz plot includes 0
            minp=minp<0		; to show Bz magnitude clearly
            ENDIF
	keep=WHERE((minp LE ytickv(*,i)) AND $
		   (ytickv(*,i) LE maxp),nkeep)
        PRINT,minb,maxb

        IF (KEYWORD_SET(polar) OR KEYWORD_SET(nonorm)) $
            THEN PLOT,uts/3600.,data $
	             ,XRANGE=[start,finish],/XSTYLE,XTITLE='' $
                     ,XTICKS=xticks,XTICKV=labuts/3600. $
                     ,XTICKNAME=xlab,XMINOR=xminor $
	             ,YRANGE=[minp,maxp],/YSTYLE,YTITLE=ytit(i) $
                     ,YTICKS=yticks(i),YTICKV=ytickv(keep,i) $
                     ,YTICKNAME=ytickname(keep,i),YMINOR=yminor(i) $
                     ,TITLE=ptitle,_EXTRA=e $
            ELSE PLOT,uts/3600.,data $
	             ,XRANGE=[start,finish],/XSTYLE,XTITLE='' $
                     ,XTICKS=xticks,XTICKV=labuts/3600. $
                     ,XTICKNAME=xlab,XMINOR=xminor $
	             ,YRANGE=[minp,maxp],/YSTYLE,YTITLE=ytit(i) $
                     ,TITLE=ptitle,_EXTRA=e
        IF i EQ 0 THEN ptitle=''

; Plot a dashed line at 0 nT and print the minimum and maximum component 
; values for the day.

        IF (!Y.CRANGE(0) LE 0) AND (!Y.CRANGE(1) GE 0) THEN $
            PLOTS,!X.CRANGE,[0,0],LINESTYLE=1
        IF NOT KEYWORD_SET(nonorm) THEN BEGIN
            XYOUTS,0.97*!X.CRANGE(0)+0.03*!X.CRANGE(1) $
                  ,0.93*!Y.CRANGE(0)+0.07*!Y.CRANGE(1) $
                  ,STRING(minb,unit,FORMAT=limform) $
                  ,ALIGNMENT=0.0
            XYOUTS,0.97*!X.CRANGE(0)+0.03*!X.CRANGE(1) $
                  ,0.16*!Y.CRANGE(0)+0.84*!Y.CRANGE(1) $
                  ,STRING(maxb,unit,FORMAT=limform) $
                  ,ALIGNMENT=0.0
            ENDIF
        ENDFOR

; Now plot the field magnitude.

    multiplot
    IF KEYWORD_SET(nonorm) THEN BEGIN
        minb=0.
        maxb=20.
        ENDIF ELSE maxb=MAX(bmag,MIN=minb)
    PLOT,uts/3600.,bmag $
        ,XRANGE=[start,finish],/XSTYLE,XTITLE='',XTICKS=xticks $
        ,XTICKV=labuts/3600.,XTICKNAME=xlab,XMINOR=xminor $
        ,YRANGE=[0.,maxb],/YSTYLE,YTITLE='|!17B!x| (nT)' $
        ,_EXTRA=e
    IF NOT KEYWORD_SET(nonorm) THEN BEGIN
        XYOUTS,0.97*!X.CRANGE(0)+0.03*!X.CRANGE(1) $
              ,0.93*!Y.CRANGE(0)+0.07*!Y.CRANGE(1) $
              ,STRING(minb,FORMAT='(F6.2," nT")') $
              ,ALIGNMENT=0.0
        XYOUTS,0.97*!X.CRANGE(0)+0.03*!X.CRANGE(1) $
              ,0.16*!Y.CRANGE(0)+0.84*!Y.CRANGE(1) $
              ,STRING(maxb,FORMAT='(F5.2," nT")') $
              ,ALIGNMENT=0.0
        ENDIF

; Now annotate time axis with times and s/c positions.

    ll=CONVERT_COORD([[start,0.,0.],[finish,5.,0.]],/DATA,/TO_DEVICE)
    hrperch=(finish-start)/((ll(0,1)-ll(0,0))/!D.X_CH_SIZE)
    ntperch=5./((ll(1,1)-ll(1,0))/!D.Y_CH_SIZE)
    XYOUTS,start-4*hrperch,0.-1.5*ntperch,'UT:',ALIGNMENT=1,/DATA
    XYOUTS,labuts/3600.,-1.5*ntperch,labhrmn,ALIGNMENT=0.5,/DATA
    XYOUTS,start-4*hrperch,0.-3.0*ntperch,'GSM X:',ALIGNMENT=1,/DATA
    XYOUTS,labuts/3600.,-3.0*ntperch,labgsmx,ALIGNMENT=0.5,/DATA
    XYOUTS,start-4*hrperch,0.-4.5*ntperch,'GSM Y:',ALIGNMENT=1,/DATA
    XYOUTS,labuts/3600.,-4.5*ntperch,labgsmy,ALIGNMENT=0.5,/DATA
    XYOUTS,start-4*hrperch,0.-6.0*ntperch,'GSM Z:',ALIGNMENT=1,/DATA
    XYOUTS,labuts/3600.,-6.0*ntperch,labgsmz,ALIGNMENT=0.5,/DATA

; If a PostScript plot was generated, close the file, remove the spurious 
; characters in the first line, and send it to the B/W printer.

    IF KEYWORD_SET(ps) THEN BEGIN
        IF !VERSION.OS EQ 'vms' THEN close_ps ELSE psclose
        psf=qgetfile(psname)
        psf(0)=STRMID(psf(0),0,2)
        putfile,psname,psf
        IF !VERSION.OS EQ 'vms' $
            THEN SPAWN,'PRINT '+psname+'/QUEUE=QMS2/DELETE'
        ENDIF

    ENDELSE

RETURN
END