;+
; NAME:
;	RDSSCOUT
; PURPOSE:
;       To read and plot (!) segments of satellite orbits that pass
;       within the Polar Camera field of view at Eureka.  The input
;       data files are generated by the 'Query' function of the
;       NODIS Satellite Situation Center at GSFC, using the 'Trace'
;	option and a 2000 km distance window centered on Eureka.
; CATEGORY:
;	?
; CALLING SEQUENCE:
;	RDSSCOUT, YYMMDD, SAT, START, FINISH [, MAXLAT = maxlat] $
;              [, /PS] [, FILE=file]
; INPUTS:
;       YYMMDD: Date of orbit segments to be plotted; a 6-digit
;               number with the first 2 digts giving the year, the
;               next two giving the month number, and the last 2
;		giving the day of the month.
; OPTIONAL INPUTS:
;       SAT:    Specifies the satellite(s) for which orbits are to
;               be plotted.  If omitted, all satellites represented
;               in the file are plotted, each on its own dial plot.
;               Strings must match the satellite names in the file
;               exactly, ignoring case.  A satellite name of '*'
;               will also get plots for all available satellites.
;               If more than one name is specified, the parameter
;               must be presented as a STRING array, i.e.,
;		['first','second'].
;	START,
;       FINISH: If set, specify the UT hours bounding the period
;               in which orbits are to be plotted, i.e., START=20
;               and END=24 will produce plots of all orbits taking
;		place between 20:00:00 and 24:00:00 UT.
; KEYWORD PARAMETERS:
;       MAXLAT: Specifies the maximum PACE geomagnetic latitude
;               that must be attained by the satellite ground track
;               in order to be plotted.  The default is 85 deg,
;               which (usually) puts the satellite F-region footprint
;               within the Polar Camera field of view from Eureka.
;       /PS:    If set, plots will be generated to a B/W PostScript
;               file, which will be queued for printing after it
;               is complete.  Otherwise, plots appear in a large
;		X-window.
;       FILE:   May be used to specify a particular data file from
;               which the orbit data are to be read.  By default,
;               the routine looks for a data file in
;               /aragorn/user2/steele/ssc with a name of the format
;               YYDD1DD2.*, where YY gives the year of the first
;               entry in the file, and DD1 and DD2 are the days of
;               the year for the first and last entries.  If the
;               file spans a year boundary, DD2 is greater than
;               365.  This parameter would be needed for NOAA12
;               data from winter 1993-1994, which is collected in
;               a file names NOAA9394.EUR.
; OUTPUTS:
;	None.
; OPTIONAL OUTPUTS:
;	None.
; COMMON BLOCKS:
;	None.
; SIDE EFFECTS:
;       Unless the FILE keyword is set, /aragorn/user2/steele/ssc
;       is searched for all files beginning with digits.  A selected
;       file is read, if uncompressed, or copied to
;       /jasper/cnsr3_data1/tmp_file.gz, if compressed, and the
;       data are read from the uncompressed version of the latter
;       file, which is deleted after the data are successfully
;       read.  A plot is created in a PostScript file in ~steele/ssc
;       if /PS is included, or in Window 0.  The CONVERT.PRO routine
;       is used to generate PACE geomagnetic coordinates for the
;       specified geographic coordinates of the field-line-traced
;	satellite positions.
; RESTRICTIONS:
;	The orbit segment file for the desired date must have been
;	created.
; PROCEDURE:
;	Later.
; EXAMPLE:
;	
; SEE ALSO:
;	CONVERT.PRO
; MODIFICATION HISTORY:
; 	Written by:	D P Steele, ISR, 950207-9.
;-

PRO rdsscout,date,sat,start,finish,maxlat=maxlat,ps=ps,file=file

IF N_PARAMS() LT 1 THEN BEGIN
    doc_library,'rdsscout'
    RETURN
    ENDIF

; Set MAXLAT to default value if necessary.
IF N_ELEMENTS(maxlat) EQ 0 THEN maxlat=85.

t0=SYSTIME(1)

; Decode the date parameter

year=date/10000L
month=(date-10000*year)/100
year=year+(year GE 50)*1900+(year LT 50)*2000
day=date MOD 100
datestr=ymd2date(year,month,day,format='y$-n$-d$')
WHILE STRPOS(datestr,' ') NE -1 DO STRPUT,datestr,'0',STRPOS(datestr,' ')

; Ensure a partial argument list is interpreted correctly

IF N_PARAMS() EQ 1 THEN sat='*'
IF N_PARAMS() EQ 3 THEN BEGIN	; SAT omitted, START and FINISH specified
    finish=start		; shift parameters to correct places
    start=sat			; and set SAT to wild-card value.
    sat='*'
    ENDIF
IF N_ELEMENTS(start) EQ 0 THEN start=0. ELSE start=start(0)
IF N_ELEMENTS(finish) EQ 0 THEN finish=24. ELSE finish=finish(0)

; No file specified, so check out all available file names

IF N_ELEMENTS(file) EQ 0 THEN BEGIN
    sscdir='/aragorn/user2/steele/ssc'
    allfiles=rfindfile(sscdir+'/*',count=nfiles)
    IF nfiles EQ 0 THEN BEGIN
        MESSAGE,'No files found!',/INFORMATIONAL
        RETURN
        ENDIF ELSE BEGIN

	; Exclude .PS files

        allfiles=allfiles(WHERE(STRPOS(STRUPCASE(allfiles),'PS') EQ -1))
        nfiles=N_ELEMENTS(allfiles)

        fjstart=LONARR(nfiles)
        fjend=fjstart
        jday=julday(month,day,year)
	ig=0
        FOR i=0,nfiles-1 DO BEGIN
	    dd=rstrpos(allfiles(i),'/')
	    filepart=STRMID(allfiles(i),dd+1,STRLEN(allfiles(i))-dd-1)
	    IF isnumber(STRMID(filepart,0,2)) NE 0 THEN BEGIN

		; File name begins with at least 2 digits

	        yr=FIX(STRMID(filepart,0,2))
	        yr=yr+(yr GE 50)*1900+(yr LT 50)*2000
	        doy=FIX(STRMID(filepart,2,3))
	        fjstart(ig)=julday(1,1,yr)+doy-1
	        doy=FIX(STRMID(filepart,5,3))
	        fjend(ig)=julday(1,1,yr)+doy-1
		ig=ig+1
		ENDIF	; (isnumber(STRMID(filepart,0,2)) NE 0)
	    ENDFOR	; (i=0,nfiles-1)

	IF ig LT nfiles THEN BEGIN
	    fjstart=fjstart(0:ig-1)
	    fjend=fjend(0:ig-1)
	    ENDIF

        whin=WHERE((fjstart LE jday) AND (jday LE fjend),nwhin)
	IF nwhin EQ 0 THEN BEGIN
	    MESSAGE,'Can''t find a file for day '+STRTRIM(jday,2) $
		   ,/INFORMATIONAL
	    RETURN
	    ENDIF ELSE BEGIN
	    file=allfiles(whin(0))
	    PRINT,file
	    PRINT,STRTRIM(date,2)
	    ENDELSE

        IF STRMID(file,STRLEN(file)-1,1) EQ 'z' THEN BEGIN	; 'gzip'ped file!
	    zflag=1B
	    SPAWN,'/usr/bsd43/bin/cp '+file+' tmp_file.gz'
	    SPAWN,'/usr/local/bin/gunzip -f tmp_file.gz'
	    file='tmp_file'
	    ENDIF ELSE zflag=0B
        ENDELSE		; (nfiles EQ 0)
    ENDIF		; (N_ELEMENTS(file) EQ 0)

inf=qgetfile(file)	; read in all lines of the text file
nl=N_ELEMENTS(inf)
name=STRARR(nl)
time=name
dstr=name
ut=LONARR(nl)
slat=FLTARR(nl)
slon=slat
srad=slat
tlat=slat
tlon=slat
trad=slat

; Now process all data lines that begin with the right date string

j=0
FOR i=0,nl-1 DO BEGIN
    CASE 1 OF
(STRMID(inf(i),0,9) EQ datestr): $
	BEGIN
	    s=inf(i)
	    nw=nwrds(s)
	    wix=0
	    dstr(j)=getwrd(s,wix) & wix=wix+1
	    time(j)=getwrd(s,wix) & wix=wix+1
	    ut(j)=LONG(secstr(time(j)))
	    name(j)=getwrd(s,wix)

	    IF isnumber(name(j)) NE 0 $		; sat. name field blank!
		THEN IF j GT 0 $		; not first line for date
		    THEN name(j)=name(j-1) $	; use last name read
		    ELSE name(j)=getwrd(inf(i-1),wix) $	; read line above
		ELSE wix=wix+1			; update word index

	    slat(j)=FLOAT(getwrd(s,wix)) & wix=wix+1
	    slon(j)=FLOAT(getwrd(s,wix)) & wix=wix+1
	    srad(j)=FLOAT(getwrd(s,wix)) & wix=wix+1
	    tlat(j)=FLOAT(getwrd(s,wix)) & wix=wix+1
	    tlon(j)=FLOAT(getwrd(s,wix)) & wix=wix+1
	    trad(j)=FLOAT(getwrd(s,wix))

	    j=j+1
	    END
  ELSE:
	ENDCASE
    ENDFOR

IF zflag THEN SPAWN,'/usr/bsd43/bin/rm '+file

t1=SYSTIME(1)
PRINT,'Reading ',nl,' lines of data took ',t1-t0,' sec' $
     ,FORMAT='(A,I0,A,F6.2,A)'

; If satellite unspecified, choose all satellites represented in file

IF sat(0) EQ '*' THEN sat=name(uniq(name,SORT(name(0:j-1))))
nsat=N_ELEMENTS(sat)

; Set up plot

IF KEYWORD_SET(ps) THEN BEGIN
    IF nsat GT 1 THEN satstr='all' ELSE satstr=sat(0)
    psfile=file+'.'+STRTRIM(date,2)+'.'+satstr+'.ps'
    psopen,psfile
    DEVICE,XSIZE=6,YSIZE=9,XOFFSET=1,YOFFSET=1,/INCHES
    psres=300
    mlr=22.5
    ENDIF ELSE BEGIN
    WINDOW,0,XSIZE=640,YSIZE=960
    mlr=17.5
    ENDELSE

!P.NOCLIP=0
!P.NOERASE=1
usersymbol,'star'	; define plot symbol for PSYM=8

FOR isat=0,nsat-1 DO BEGIN	; loop through all satellites selected

    wsat=WHERE(name EQ STRUPCASE(sat(isat)),nl)

    IF nl GT 0 THEN BEGIN	; identify relevant data
        PRINT,STRUPCASE(sat(isat))
        sdstr=dstr(wsat)
        stime=time(wsat)
        sut=ut(wsat)
        sname=name(wsat)
        sslat=slat(wsat)
        sslon=slon(wsat)
        ssrad=srad(wsat)
        stlat=tlat(wsat)
        stlon=tlon(wsat)
        strad=trad(wsat)

        tmlat=FLTARR(nl)
        tmlt=tmlat
        FOR i=0,nl-1 DO BEGIN	; convert GEO coordinates to PACE GM
	    tj=sut(i)
            hour=tj/3600L
            minute=(tj/60) MOD 60
            convert,year,month,day,hour,minute,0,stlat(i),stlon(i),0.0 $
	           ,mlat,mlt
            tmlat(i)=mlat
            tmlt(i)=mlt
	    ENDFOR

        ENDIF ELSE BEGIN
        MESSAGE,'No entries for '+STRUPCASE(sat)+' on '+STRTRIM(date,2) $
	       ,/INFORMATIONAL
        RETURN
        ENDELSE


; Determine plot limits in DEVICE coordinates

    IF KEYWORD_SET(ps) $
	THEN posn=[scl((isat MOD 2)*900,res=psres) $
	          ,scl(((5-isat)/2)*900,res=psres) $
	          ,scl(((isat MOD 2)+1)*900-1,res=psres) $
	          ,scl((((5-isat)/2)+1)*900-1,res=psres)] $
        ELSE posn=[(isat MOD 2)*320 $
	          ,((5-isat)/2)*320 $
	          ,((isat MOD 2)+1)*320-1 $
	          ,(((5-isat)/2)+1)*320-1]
    !P.CLIP=posn

; Establish DATA coordinates

    PLOT,[0,20],[-180,180]*!DTOR,/POLAR,/NODATA,/DEVICE $
        ,POSITION=posn,XRANGE=[-20,20],XSTYLE=5 $
        ,YRANGE=[-20,20],YSTYLE=5

; Draw and label polar dial plot in MLAT/MLT

    rt=28.
    tt=315.*!DTOR
    XYOUTS,rt*COS(tt),rt*SIN(tt),'!17'+sdstr(0) $
          ,ALIGNMENT=1,CHARTHICK=2
    tt=225.*!DTOR
    XYOUTS,rt*COS(tt),rt*SIN(tt),'!17'+sname(0) $
          ,ALIGNMENT=0,CHARTHICK=2
    r2=SQRT(2.)
    len=[20,20*r2,20,20*r2,20,20*r2,20,20*r2]
    FOR i=0,7 DO BEGIN
        ang=45*i-90
        radii,0,len(i),ang
        XYOUTS,mlr*COS(ang*!DTOR),mlr*SIN(ang*!DTOR) $
              ,STRTRIM(FIX(3*i),2),ALIGNMENT=1
        ENDFOR
    FOR i=1,5 DO BEGIN
        arcs,5*i,0,360,NOCLIP=0
        ang=(45.-15./i)*!DTOR
        XYOUTS,5*i*COS(ang),5*i*SIN(ang),STRTRIM(90-5*i,2) $
              ,ALIGNMENT=0
        ENDFOR

; Locate the beginnings of passes, by the time gaps between data lines

    dut=(sut-SHIFT(sut,1))>0
    startpass=[0,WHERE(dut GT 3600),nl]

    ips=0
    FOR i=0,N_ELEMENTS(startpass)-2 DO BEGIN
        IF (start*3600. LE sut(startpass(i))) AND $
           (sut(startpass(i)) LT finish*3600.) THEN BEGIN
	    j1=startpass(i)
	    j2=startpass(i+1)-1
	    IF j2 GT j1 THEN BEGIN	; ensure no degenerate orbit segments

; Calculate PACE coordinates of Eureka at mid-time of orbit segment

	        utmid=(sut(j1)+sut(j2))/2
	        hour=utmid/3600L
	        minute=(utmid/60L) MOD 60
	        second=utmid MOD 60
	        convert,year,month,day,hour,minute,second,80.05,273.58,0.6 $
	               ,emlat,emlt
                r=90.-[tmlat(j1),tmlat(j2),emlat]
                theta=([tmlt(j1),tmlt(j2),emlt]*15.-90.)*!DTOR

; How close does the satellite get to Eureka?

	        x=r*COS(theta)
	        y=r*SIN(theta)
	        a=[x(0),y(0)]	; Cartesian coordinates of first point on footprint
	        b=[x(1),y(1)]	;     "          "      "  last    "
	        e=[x(2),y(2)]	;     "          "      "  Eureka

	        xn=b(1)-a(1)	; X-component of normal vector
	        yn=a(0)-b(0)	; Y-component of normal vector
	        normvec=[xn,yn]/SQRT(xn*xn+yn*yn)	; normalized ( |N| = 1 )
	        c=-TOTAL(a*normvec)

; Points on orbit track in MLAT/MLT satisfy NORMVEC . P + C = 0

	        q=TOTAL(normvec*e)+c
	        closest=e-q*normvec

; Vector from Eureka to closest point on orbit is parallel to NORMVEC

	        clmlat=90.-SQRT(TOTAL(closest*closest))
	        clmlt=((ATAN(closest(1),closest(0))*!RADEG+90.)/15.+24.) MOD 24
	        PRINT,strsec(utmid),clmlat,clmlt,FORMAT='(A,2(2X,F5.2),$)'

	        IF clmlat GT maxlat THEN BEGIN	; orbit gets close enough to
	    					; Eureka to be worth plotting
                    PRINT,' X'

; Set up a 1-min resolution time grid during the orbit, and
; interpolate MLAT/MLT positions on the grid

	            ptime=[stime(j1),stime(j2)]
                    tabut=secstr(stime(j1:j2))	; UT of all available points
		    npts=(tabut(N_ELEMENTS(tabut)-1)-tabut(0))/60.+1
		    orbut=FINDGEN(npts)*60+tabut(0)	; regular grid of UT values
		    tabr=90.-tmlat(j1:j2)		; 2-D polar coords for all
		    tabtheta=(tmlt(j1:j2)*15.-90.)*!DTOR	; available points
		    tabx=tabr*COS(tabtheta)		; convert polar coords 
		    taby=tabr*SIN(tabtheta)		; -> Cartesian coords

; Interpolate in Cartesian coordinates at 1-min intervals, and save
; interpolated values for later reference

		    orbx=interpol(tabx,tabut,orbut)
		    orby=interpol(taby,tabut,orbut)
		    IF ips EQ 0 THEN BEGIN
			allx=orbx
			ally=orby
			alltime=ptime
			ENDIF ELSE BEGIN
			allx=[allx,orbx]
			ally=[ally,orby]
			alltime=[alltime,ptime]
			ENDELSE

		    orbmlat=90.-SQRT(orbx*orbx+orby*orby)
		    orbmlt=((ATAN(orby,orbx)*!RADEG+90.)/15.+24.) MOD 24

; Plot (interpolated) satellite positions at 1-min intervals

		    ips=(ips+1) MOD 9	; recycle if too many passes
		    IF ips EQ 3 THEN ips=ips+1
                    OPLOT,90.-orbmlat,(orbmlt*15.-90.)*!DTOR $
		         ,/POLAR,PSYM=-ips,SYMSIZE=0.5

; Plot line joining Eureka to closest point on orbit

                    OPLOT,90.-[emlat,clmlat],([emlt,clmlt]*15.-90.)*!DTOR $
		         ,/POLAR
                    devco=CONVERT_COORD([TRANSPOSE(r*COS(theta)) $
				        ,TRANSPOSE(r*SIN(theta))] $
                                       ,/DATA,/TO_DEVICE)
	            algn=[1,0]
	            IF devco(0,0) GT devco(0,1) THEN algn=[0,1]
	            offdir=0.5-algn
	            FOR k=0,1 DO $
                        XYOUTS,devco(0,k)+offdir(k)*!D.X_CH_SIZE $
		              ,devco(1,k)-0.2*!D.Y_CH_SIZE $
                              ,STRMID(ptime(k),0,5) $
			      ,/DEVICE,ALIGNMENT=algn(k) $
			      ,CHARSIZE=0.5
	            PLOTS,devco(0,2),devco(1,2),PSYM=6 $
			 ,/DEVICE,SYMSIZE=0.5
                    ENDIF ELSE BEGIN
		    PRINT,FORMAT='()'
		    IF clmlat LT 69.0 THEN BEGIN
		        SAVE,/ALL,FILENAME='rdsscout.debug.save'
		        RETURN		; quit here
			ENDIF		; (clmlat LT 69.0)
		    ENDELSE		; (clmlat GT maxlat)

	        ENDIF	; (j2 GT j1)

            ENDIF	; (start*3600. LE sut(i)) ...
        ENDFOR		; (i=0,N_ELEMENTS(startpass)-2)

; Make a table of orbit times in the least populated quadrant of
; the dial plot

    IF ips GT 0 THEN BEGIN	; at least one orbit was plotted!
        quad=INTARR(4)	; number of orbit points in each quadrant
        quad(0)=TOTAL((allx GE 0) AND (ally GE 0))
        quad(1)=TOTAL((allx LT 0) AND (ally GE 0))
        quad(2)=TOTAL((allx LT 0) AND (ally LT 0))
        quad(3)=TOTAL((allx GE 0) AND (ally LT 0))
        qmin=MIN(quad,qi)	; choose least populated quadrant
        IF (qi EQ 0) OR (qi EQ 3) $
	    THEN xt=(posn(0)+posn(2)+1)/2 $
	    ELSE xt=posn(0)
        IF (qi EQ 2) OR (qi EQ 3) $
	    THEN yt=(posn(1)+posn(3)+1)/2 $
	    ELSE yt=posn(3)
        xt=xt+2*!D.X_CH_SIZE	; initialize text positions
        yt=yt-3*!D.Y_CH_SIZE
        psa=[1,2,4,5,6,7,8]

; Show symbol and start, finish times for each orbit plotted

        FOR i=0,N_ELEMENTS(alltime)/2-1 DO BEGIN
	    PLOTS,xt+0.5*!D.X_CH_SIZE,yt+0.5*!D.X_CH_SIZE $
	         ,PSYM=-psa(i),/DEVICE
            XYOUTS,xt+2*!D.X_CH_SIZE,yt $
	          ,STRMID(alltime(2*i),0,5)+' - ' $
	          +STRMID(alltime(2*i+1),0,5) $
	          ,/DEVICE
            yt=yt-1.25*!D.Y_CH_SIZE
	    ENDFOR		; (i=0,N_ELEMENTS(alltime)/2-1)
        ENDIF

    ENDFOR		; (isat=0,nsat-1)

IF KEYWORD_SET(ps) THEN BEGIN
    psclose
    SPAWN,'lpr -l -r -Pps '+psfile
    ENDIF

RETURN
END