;+
;PROCEDURE:	plot_fa_crossing
;
;PURPOSE:
;
;	Plots magnetic footprint of spacecraft across earth. (FLAT,FLNG)
;	Shows the night/day terminator, auroral ovals, apogee and perigee
;	footprints, etc.  (See examples below.)
;
;       Keywords are grouped into categories: TIMESPAN, VIEW, DISPLAY,
;       AURORAL ZONE, MISCELLANEOUS.
;
; TIMESPAN KEYWORDS:
;
;	ORBIT	The orbit to plot. If unset, will plot over interval
;		containing present time and show craft position right now
;		(unless TMIN and TMAX are set).  ORBIT cannot be greater than
;		the last orbit listed in the predicted orbit almanac file.
;  TMIN, TMAX	Time points to be included on the chart.  Use these if
;               you want to display timespan in distant future.
;               Should not be more than a couple hours apart.  Format
;               of TMIN and TMAX must be the type of string accepted
;               by time_double() or a double float in seconds since
;               1970.  ORBIT must not be set if these keywords are to
;               be used.  These times will be labeled on the map as t1
;               and t2.  (Good for showing AOS and LOS taken from
;               contact schedules.)
;	XMARK	Set this to a time (string or double float) to have that point
;               labeled on the map as a big X.  If none of ORBIT,
;               TMIN, TMAX are set then this time will be used as the
;               reference time for which to create the plot.  (Good
;               for showing conjunctions.)
;
; VIEW KEYWORDS:
;
;   KIRUNA, POKER, WALLOPS, MCMURDO, CANBERRA, SANTIAGO, BERKELEY:
;		Set any of these for an overhead view of the station.
;   VIEWPOINT   This allows arbitrary views of the Earth. Should be a
;               3-element array designating lattitude, longitude, and
;               rotation in degrees.
;	WHOLE	If set, will not confine plot to polar regions.
;	SOUTH	Set to view earth from directly over South Pole.
;     MAGPOLE   View from above the magnetic pole with magnetic local
;               noon at the top of the plot.  Gridlines are still
;               geographic.  Poles are not symmetric because of the
;               eccentric dipole.  (When this keyword is not set
;               the display defaults to above the geographic poles,
;               with geographic local noon at the top.)
;        ZOOM   Enlarges the map. 1 is normal, 2 is twice as big.
;               See WINSIZE keyword for zooming in postscript mode.
;
; DISPLAY KEYWORDS:
;
;     WINSIZE	The width of the plot window in pixels. Default width is
;               640.  Height is scaled automatically so that lattitude
;               lines will be circular when output on a Tek printer:
;               H = W * 1.031.  In PostScript mode (POST keyword set),
;               WINSIZE acts like a magnification factor; WINSIZE=10
;               is normal.
;	  PC	Set this to nonzero if you are using a windowing system
;		that does not provide backing store to retain hidden
;		windows.  This keyword also sets the VECTOR_FONTS keyword.
;	GREY	Set this keyword to output the plot in greyscale.
;		(Run @startup to restore colortable.)
;       FILL    Fills oceans and continents with solid color to make a
;               pretty plot.  Setting this keyword when printing to
;               a color printer is discouraged.
;  VECTOR_FONTS Disable switch to device (hardware) font from
;               (default) Hershey vector-drawn fonts.  Device fonts
;               may be prettier, but vector fonts are device
;               independent.
;	POST    Set this to a filename to direct the ouput to an 8-bit color
;               postscript file instead of the graphics window.  When
;               viewing the postscript file, be sure your viewer is
;               set to 8 bits, not 24, and that it is displaying
;               "perfect colors".  (Do not add ".ps" to the name.)
;        GIF    Captures the image on the output graphics window to a
;               GIF file.  (Add ".gif" to the name yourself.)
;
; AURORAL ZONE KEYWORDS:
;
;    ACTIVITY	Set the Activity Index of aurora. ( 0-6, default 3)
;     SSCZONE   Show the auroral zone used by the Satellite Situation
;               Center instead of the one from Holzworth and Meng.
;               (THIS IS NOT PERFECTED YET.) See auroral_zone_ssc.pro
;               and notes below.
;
; MISCELLANEOUS KEYWORDS:
;
;      POLAR    Plots the magnetic footprint of the POLAR spacecraft
;               and information about its closest approach to that of
;               FAST.  Good for finding conjunctions or displaying
;               known ones. (For a description of how POLAR orbit data
;               is obtained, see the documentation for
;               get_po_orbit.pro.)
;   DRAG_PROP   If set, the orbit propagator will include the effects
;               of atmospheric drag.  The orbit track may be
;               inaccurate by a few degrees for propagations as short
;               as a few weeks.  This keyword minimizes this error.
;ALMANAC_INFO   Prints a line under the plot title telling the orbit
;               file used, the last epoch in that orbit file, and
;               whether or not drag was included in the orbit propagation.
;
;   EXAMPLES:   Show what FAST is doing now, and make it pretty:
;
;               IDL>  plot_fa_crossing, /fill
;
;               Show what FAST is doing for Christmas, make output
;               greyscale, and copy output to a GIF file:
;
;               IDL>  plot_fa_crossing, xmark='98-12-25/00:00:00', $
;                     gif='~/images/northpole.gif', /grey
;
;               Show the polar crossing of orbit 3000 from a
;               bird's-eye view over Canberra, and send the output to
;               an 8-bit color postscript file:
;
;               IDL>  plot_fa_crossing, /can, post='~/images/orb3000', $
;                     orbit=3000
;
;NOTES ON THE
;AURORAL ZONE:	The auroral zone is drawn using a procedure by Jim
;               Clemens (suggested by Holzworth & Meng, and corrected
;               by Mike Temerin).  The source for this procedure is
;               "Mathematical Representation of the Auroral Oval",
;               R.H. Holzworth and C.-I. Meng, Geophysical Research
;               Letters, Vol.2, No.9, Sept 1975.
;
;               The mathematical representation of the auroral ovals
;               is found by fitting Feldstein statistical ovals in
;               corrected geomagnetic coordinates to a 7-parameter
;               Fourier series:
;	
;		theta =  A1  +  A2 cos(phi + A3)
;		      +  A4 cos(2phi + 2A5)
;		      +  A6 cos(3phi + 3A7)
;
;		theta = corrected geomagnetic co-lattitude
;		phi = 2(pi)(MLT)/(24hrs)
;
;               The best fit constants are found for each value of
;               Q=0..6, where Q is the activity index which describes
;               how quiet (0) or active (6) the Feldstein auroral
;               ovals are.  The characteristic radius A1 of the ovals
;               increases monotonically with Q. (Source: Holzworth and
;               Meng.)
;
;               The function azonloc.pro generates the southern
;               auroral oval by reflection of the northern oval
;               through the magnetic equator.
;
;               The formula in Holzworth and Meng gives the ovals in
;               corrected geomagnetic coordinates.  The procedure
;               transform_mag_geo.pro converts them to geograhic
;               coordinates using the eccentric dipole model.
;
;CREATED BY:    J.Rauchleiba            96-12-20
;-
pro plot_fa_crossing, $
	ORBIT=orbit, $
	TMIN=time1, $
	TMAX=time2, $
        XMARK=xmarkin, $
        ZOOM=zoom, $
        DRAG_PROP=drag_prop, $
        ALMANAC_INFO=almanac_info, $
	WHOLE=whole, $
	VIEWPOINT=pov, $
	KIRUNA=kir, $
	POKER=pok, $
	WALLOPS=wal, $
	MCMURDO=mcm, $
	CANBERRA=can, $
	SANTIAGO=san, $
	BERKELEY=ber, $
        SOUTH=south, $
        MAGPOLE=magpole, $
        SSCZONE=ssc, $
	ACTIVITY=act, $
	WINSIZE=win, $
	PC=pc, $
        GREY=grey, $
        FILL=fill, $
        VECTOR_FONTS=vdf, $
        POST=psfile, $
        GIF=giffile, $
                      POLAR=polar, $
                      USER_OVERPLOT = uop

thm_init
loadct2,43
cols=get_colors()
window,10,xsize=500,ysize=500

; FIDDLE WITH THE COLORS

;n_colors = !d.n_colors		; Number of colors already loaded
n_colors = !d.table_size		; Number of colors already loaded
black = 0			; Assuming loadct2, 39
white = n_colors - 1
blue = fix(.33*float(n_colors))
if keyword_set(grey) then begin		; Make Greyscale colortable
    g1 = indgen(256)
    g2 = indgen(256)
    g3 = indgen(256)
    g1 = congrid(g1, n_colors)	; Sample n_colors of new table
    g2 = congrid(g2, n_colors)	; so size of new table equal to
    g3 = congrid(g3, n_colors)	; size of old table.
    tvlct, g1, g2, g3
    col_seas = white		; Seas
    col_cont = fix(n_colors*4/5) ; Continents
    col_fast = fix(n_colors/4)	; Spacecraft
    col_term = fix(n_colors/2)	; Terminator line
    col_oval = fix(n_colors/3)	; Auroral Ovals
    col_tags = black		; Various Labels
    erase, white                ; PS device non-supportive
endif else begin
    col_seas = blue
    col_cont = fix(.81*float(n_colors))
    col_fast = 6
    if keyword_set(fill) then col_term=5 else col_term=blue
    col_oval = 4
    if keyword_set(fill) then col_tags=white else col_tags=black
    ;; PS device does not support erase command.
    ;; Must use polyfill if /FILL in this case.
    if !d.name EQ 'PS' AND keyword_set(fill) then begin
        polyfill, [0,1,1,0,0], [0,0,1,1,0], /normal, col=black
    endif else begin 
        if keyword_set(fill) then erase, black else erase, white
    endelse
endelse

; FONT SETUP

oldfont = !p.font
if NOT keyword_set(vdf) and NOT keyword_set(pc) then !p.font = 0

; VIEWPOINT AND PLOT LIMITS

case 1 of
	keyword_set(kir): 	pov = [67.883,21.067,0 ]
	keyword_set(pok): 	pov = [64.8,-147.85,0 ]
	keyword_set(wal):	pov = [37.9,-75.4,0]
	keyword_set(mcm):	pov = [-77.85,166.62,0]
	keyword_set(can):	pov = [-35.35,149.167,0]
	keyword_set(san):	pov = [-33.433,-70.667,0]
        keyword_set(ber):	pov = [37.9,-122.1,0]
        else:
endcase
if keyword_set(whole) then latlim=[-90, 90] else latlim = [40, -40]
if keyword_set(south) then begin
    lat_ref = -90.
endif else begin
    if keyword_set(pov) then lat_ref = pov(0) else lat_ref = 90.
endelse
;;ns_sign = (lat_ref/abs(lat_ref) + 1)/2  ; 1:0 ::  N:S

; GET THE MAGNETIC PRINT DATA OF FAST

;get_fa_orbit, tmin, tmax, /all, orbit_file=orbit_file, drag_prop=drag_prop,/no_sto, struc=fa_orbit, status=status
fa_orbitrange,orbit
trange=timerange()
orbits=fa_time_to_orbit(trange)
if orbits[0] eq orbits[1] then orbit_string = string(orbits[0]) else orbit_string=string(orbit[0])+'-'+string(orbit[1])

fa_k0_load,'orb'
;if status NE 0 then message, 'Error returned by get_fa_orbit.'
get_data, 'flat', data=tmp1	
get_data, 'flng', data=tmp2
clock = tmp1.x
flat = tmp1.y
flng = tmp2.y
tmin = trange(0)
tmax = trange(1)
;clock = fa_orbit.TIME
;flat = fa_orbit.FLAT
;flng = fa_orbit.FLNG
Mtheta = flat
Mphi = flng

; Set time_N to when craft nearest N/S pole if orbit keyword set
; or median of tmin and tmax if they are specified
; or right now if no time specification.
; middle_time will be set to either of the latter two if orbit not set.
; This variable used by both the terminator and auroral oval generators.

if keyword_set(orbit) then begin
    dummy = min(abs(lat_ref - Flat), sub_N) 
    time_N = clock(sub_N) 
endif else time_N=middle_time


; GET LAT AND LNG OF SUN 

center_time = time_N
if (!VERSION.RELEASE LE '5.4') then begin
    hr_min = (str_sep((str_sep(time_string(center_time),'/'))(1),':'))(0:1)
endif else begin
    hr_min = (strsplit((strsplit(time_string(center_time),'/',/EXTRACT))(1),':',/EXTRACT))(0:1)
endelse
gmt_hrs = float(hr_min(0)) + float(hr_min(1))/60.0
noon_long  = (!pi - gmt_hrs*(!pi/12.)) MOD (2*!pi)
t0 = time_double('96-12-21/0:00') ;
ang = (center_time - t0)/(365.25d*24d*3600d)*2d*!pi
tilt = .410152
noon_lat = -(tilt * cos(ang))

; HANDLE SPECIAL POLE VIEWS

if keyword_set(magpole) then begin ; MAG pole view
    if keyword_set(south) then begin
        pov = [-75.3, 118.6, 0.0] ; override
        rotcoeff = -1
        rot_add = 0.0
    endif else begin
        pov = [82.7, -92.0, 0.0]  ; override
        rotcoeff = +1
        rot_add = !pi
    endelse
    mnoon_long = noon_long -pov(1)*!dtor ; Sm pole is 92 deg W of 0.
    rotate = (mnoon_long)*rotcoeff + rot_add
    rerange, rotate, /deg
    pov(2) = rotate
endif else if NOT keyword_set(pov) then begin ; GEO pole view
    if keyword_set(south) then begin
        pov = [-90., 0., 0.]
        rotcoeff = -1
        rot_add = 0.0
    endif else begin
        pov = [90., 0., 0.]
        rotcoeff = +1
        rot_add = !pi
    endelse
    rotate = noon_long*rotcoeff + rot_add
    rerange, rotate, /deg
    pov(2) = rotate
endif

; MAKE THE MAP (This sets up the plotting axes.)

if keyword_set(orbit) $
  then maptitle = 'ORBIT ' + strtrim(string(orbit), 2)+'  '+time_string(tmin) $
else maptitle = time_string(tmin) + '  ' + time_string(tmax)

maptitle = '!3' + maptitle
if keyword_set(zoom) then scale = 81.25e6/float(zoom)/(win/640.)
map_set, /ortho, pov(0), pov(1), pov(2), /noerase, /noborder, scale=scale, $
  horizon=keyword_set(fill), e_horizon={fill:keyword_set(fill), color:blue}

xyouts, .5, .975, /norm, align=.5, color=col_tags, charsize=1.3, maptitle
map_continents, fill=fill, color=col_cont
map_grid, glinestyle=0, color=0, label=2, latlab=-20, latalign=0.5, $
  lats = [-90,-75,-60,-45,-30,-15,0,15,30,45,60,75,90], $
  lonlab=.75*lat_ref, lonalign=0.5, $
  lons = [-180,-135,-90,-45,0,45,90,135,135,180]

; If ORBIT PROPAGATION INFO IS DESIRED

if keyword_set(almanac_info) then begin
    if (!VERSION.RELEASE LE '5.4') then begin
        ofpieces = str_sep(orbit_file, '/')
    endif else begin
        ofpieces = strsplit(orbit_file, '/', /EXTRACT)
    endelse
    ofname = ofpieces(n_elements(ofpieces) - 1)
    if keyword_set(drag_prop) then drag_inclu='with drag.' $
    else drag_inclu='without drag.'
    last_epoch = time_string(find_last_epoch(orbit_file))
    propinfo = 'Propagation: '+ofname+' ('+last_epoch+') '+drag_inclu
    xyouts, .29, .955, propinfo, /norm, color=col_tags, charsize=.90
endif

; SHOW SOME GROUND STATIONS AND OTHER POINTS OF INTEREST

citylng = [-147.85,21.067,-70.667,-75.4,166.62,149.167,-122.2,-92.0,118.6]
citylat = [64.8,67.883,-33.433,37.9,-77.85,-35.35,37.9,82.7,-75.3]
citytag = ['!3 Poker',' Kiruna',' Santiago',' Wallops',' McMurdo',$
           ' Canberra',' Berkeley',' !6N!3',' !6S!3']
oplot, citylng, citylat, psym=2, color=black
xyouts, citylng, citylat, citytag, color=black

;;xyouts, -147.85, 64.8, '!3*Poker', color=black
;;xyouts, 21.067, 67.883, '*Kiruna', color=black
;;xyouts, -70.667, -33.433, '*Santiago', color=black
;;xyouts, -75.4, 37.9, '*Wallops', color=black
;;xyouts, 166.62, -77.85, '*McMurdo', color=black
;;xyouts, 149.167, -35.35, '*Canberra', color=black
;;xyouts, -122.2, 37.9, '*Berkeley', color=black
;;xyouts, -92.0, 82.7, '*!6N!3', color=black  ; Magnetic North
;;xyouts, 118.6, -75.3, '*!6S!3', color=black ; Magnetic South

; GET POSITION OF FAST IN GSE TO FIND WHEN IN SHADOW

;;get_data, 'fa_pos', data=tmp	
get_data,'r',data=tmp	
;pos_arr = fa_orbit.fa_pos/6372.1 ; Position in GEI, km --> Re
store_data, 'fa_pos_re', data={x:clock, y:tmp.y/6372.1}
cotrans, 'fa_pos_re', 'fa_pos_gse', /GEI2GSE
get_data, 'fa_pos_gse', data=gse_stc
store_data, 'fa_pos_re', /delete
store_data, 'fa_pos_gse', /delete

; Array shade holds indices of data that are in shadow
; Craft in shade if x < 0 and y^2 + z^2 < Re^2 (1)

shade = where( gse_stc.y(*,0) LT 0 AND $
	(gse_stc.y(*,1)*gse_stc.y(*,1) +  gse_stc.y(*,2)*gse_stc.y(*,2)) LE 1 )
n_datapts = n_elements(gse_stc.x)
n_shadepts = n_elements(shade)
if shade(0) NE -1 then begin		; If there are  points in shade
	; The index to the shade array of the last shade point before sunrise:
	shade_break = (where(shade NE indgen(n_shadepts)+shade(0))) (0) - 1
	; data_break is index of last data point in shade before sunrise.
	; If the data pts ref'd by shade not all adjacent, 
	; set data_break to the index of the last point in the continuous segment, 
	; otherwise just set it to last point ref'd by shade.
	; (If there was no break in the shade then shade_break equals -2)
	if shade_break NE -2 then data_break = shade(shade_break) $
	else data_break = shade(n_shadepts - 1)
	; Assume, if shaded (shade(0) ne -1), data_break is either the last point
	; in continuous shade before sunlight (shade_break EQ -2), or is the last
	; point in the first interval of shade where the shade is broken up into
	; two intervals (shade_break NE -2).  data_break is never the last point
	; in the data array unless daybreak happens just before timespan.  In this
	; case, data_break = n_datapts - 1 AND shade_break EQ -2.
endif
	
; PLOT THE PATH OF THE CRAFT IN N AND S POLAR REGIONS (UNLESS /WHOLE).

oplot, Mphi, Mtheta, min_value=latlim(0), color=col_fast, linestyle=0, thick=3

; If we didn't already plot all the necessary points with the last plot command

if not keyword_set(whole) then $
oplot, Mphi, Mtheta, max_value=latlim(1), color=col_fast, linestyle=0, thick=3

;if defined(uop) then begin ; user-defined overplot routine
if keyword_set(uop) then begin ; user-defined overplot routine
    call_procedure, uop
endif


; OVERPLOT A BLACK LINE WHERE CRAFT IS IN SHADOW

if shade(0) NE -1 then begin		; if there are some points in shade
	oplot, Mphi(shade(0):data_break), Mtheta(shade(0):data_break), $
	min_value=latlim(0), color=black, linestyle=0, thick=2
	if not keyword_set(whole) then begin	; if we didn't aready plot all
		oplot, Mphi(shade(0):data_break), Mtheta(shade(0):data_break), $
		max_value=latlim(1), color=black, linestyle=0, thick=2
	endif

	; We must plot the remaining shady points if there was a break in the shade

	if data_break LT (n_datapts - 1) AND shade_break NE -2 then begin
		last_shade = shade(n_elements(shade) - 1)
		oplot, Mphi(shade(shade_break+1):last_shade), $
		Mtheta(shade(shade_break+1):last_shade), $
		min_value=latlim(0), color=black, linestyle=0, thick=2
		if not keyword_set(whole) then begin
			oplot, Mphi(shade(shade_break+1):last_shade), $
			Mtheta(shade(shade_break+1):last_shade), $
			max_value=latlim(1),color=black,linestyle=0,thick=2
		endif
	endif	
endif


if 0 then begin			; this is no longer needed

; Show print of craft right now if "now" is between tmin, tmax.
; Write text to plot window telling UT, ORBIT, FLAT, FLNG right now.

if (plot_craft) then begin
	dt = min( abs(clock - now_time), pos_ind)	; pos_ind=nowdataindex
	oplot, [Flng(pos_ind)], [Flat(pos_ind)], psym=6, syms=3, col=col_fast
        ;;get_data, 'ORBIT', data=orbit_nums
	get_data
;        orbit_nums = fa_orbit.ORBIT
        orbit_nums = fa_orbit.ORBIT
	current_orb = orbit_nums(pos_ind)
        if (!VERSION.RELEASE LE '5.4') then begin
	  xyouts, .8, .95, /norm, color=col_tags, $
            'CURRENT FAST INFO' + $
            '!C------------' + $
            '!CUT: ' + (str_sep(time_string(now_time), '/'))(1) + $
            '!CORBIT: ' + strtrim(string(current_orb), 2) + $
            '!CFLAT: ' + string(format='(F6.1)', Flat(pos_ind)) + $
            '!CFLNG: ' + string(format='(F6.1)', Flng(pos_ind))
        endif else begin
	  xyouts, .8, .95, /norm, color=col_tags, $
            'CURRENT FAST INFO' + $
            '!C------------' + $
            '!CUT: ' + (strsplit(time_string(now_time), '/', /EXTRACT))(1) + $
            '!CORBIT: ' + strtrim(string(current_orb), 2) + $
            '!CFLAT: ' + string(format='(F6.1)', Flat(pos_ind)) + $
            '!CFLNG: ' + string(format='(F6.1)', Flng(pos_ind))
        endelse
endif

endif

;; ADD TIME TICKS, LABELS TO PATH

label_foot_ticks, time=clock, latitude=Flat, longitude=Flng, interval=300, $
  latlim=45, color=col_tags

; LABEL TIME1 AND TIME2 IF SUPPLIED EXPLICITLY BY USER

if keyword_set(time1s) AND keyword_set(time2s) AND NOT keyword_set(orbit) $
then begin
	dt1 = min( abs(clock - time1s), time1_ind )
	dt2 = min( abs(clock - time2s), time2_ind )
	xyouts, Flng(time1_ind), Flat(time1_ind), '.t1', $
		align=0.0, color=col_tags
	xyouts, Flng(time2_ind), Flat(time2_ind), '.t2', $
		align=0.0, color=col_tags
endif

; PUT AN X MARK AT XMARK 

if keyword_set(xmark) then begin
    dtx = min( abs(clock - xmark), x_ind )
    oplot, [Flng(x_ind)], [Flat(x_ind)], psym=7, symsize=3.0, color=col_tags
endif

; PLOT THE TERMINATOR

terminator, time=time_N, tlat, tlng	; tlat, tlng as yet undefined
; The use of max_value and min_value below is a patch to remove ugly line
if pov(0) EQ 90 then min_term = 10 else min_term = -90
if pov(0) EQ -90 then max_term = -10 else max_term = 90
oplot, tlng, tlat, thick=3, color=col_term, min_val=min_term, max_val=max_term

; APOGEE, PERIGEE

;;get_data, 'ALT', data=alt
;;get_data, 'LAT', data=Glat
get_data,'alt',data=tmp & alt=tmp.y
get_data,'ilat',data=tmp & ilat=tmp.y
;alt = fa_orbit.ALT
;Glat = fa_orbit.LAT
max_alt = max(alt, app_sub)
min_lat = min(alt, per_sub)
xyouts, Flng(app_sub), Flat(app_sub), 'A', /data, color=col_tags
xyouts, Flng(per_sub), Flat(per_sub), 'P', /data, color=col_tags
;apogee_lat = Glat(app_sub)	; Geographic lattitude of apogee
apogee_lat = ilat(app_sub)	; Invariant lattitude of apogee
perigee_lat = ilat(per_sub)	; Invariant lattitude of apogee

; AURORAL OVALS

IF n_elements(act) NE 1 then act = 3		; Default activity index

; Get Corrected Geomagnetic Lattitudes for all MLT

mlt = findgen(241)/10.          ; [.1, .2, ..., 23.9, 24.0]
if keyword_set(ssc) then begin
    nelat = auroral_zone_ssc(mlt, poleward=nplat)
    selat = auroral_zone_ssc(mlt, poleward=splat, /south)
endif else begin
    nelat = auroral_zone(mlt,act,/lat) ; Corr. Geomag. Lat for all MLT
    nplat = auroral_zone(mlt,act,/lat,/pole) ; Northern poleward
    selat = auroral_zone(mlt,act,/lat,/south) ; Southern equatorward
    splat = auroral_zone(mlt,act,/lat,/south,/pole) ; Southern poleward
endelse

if (!VERSION.RELEASE LE '5.4') then begin
    date_time = str_sep(time_string(time_N), '/')	; time_N is reference time
    year = fix(strmid(date_time(0), 0, 4))		; 4-digit integer
    N_hms = fix(str_sep(date_time(1),':'))		; [hh,mm,ss]
    UT_hrs = N_hms(0) + N_hms(1)/60.		; Hours into UT day
endif else begin
    date_time = strsplit(time_string(time_N), '/', /EXTRACT)	; time_N is reference time
    year = fix(strmid(date_time(0), 0, 4))		; 4-digit integer
    N_hms = fix(strsplit(date_time(1),':', /EXTRACT))		; [hh,mm,ss]
    UT_hrs = N_hms(0) + N_hms(1)/60.		; Hours into UT day
endelse

; Convert MLT to MLNG using UT and get MLNG of Sun

sun_lng = [!pi - UT_hrs*(!pi/12.)] ; LNG of Sun
sun_lat = [0.]
mag_to_geo, sun_lat, sun_lng, /mag ; MLNG of Sun
azon_mlng = mlt*(!pi/12.) + sun_lng(0) - !pi ; MLNG of oval
nelon = azon_mlng
nplon = nelon
selon = nelon
splon = nelon

; Convert auroral ovals from MAG to GEO

transform_mag_geo, nelat, nelon, tnelat, tnelon, year=year
transform_mag_geo, nplat, nplon, tnplat, tnplon, year=year
transform_mag_geo, selat, selon, tselat, tselon, year=year
transform_mag_geo, splat, splon, tsplat, tsplon, year=year

rerange, tnelon, tnelat, /deg
rerange, tnplon, tnplat, /deg
rerange, tselon, tselat, /deg
rerange, tsplon, tsplat, /deg

; Plot the ovals

oplot, tnelon, tnelat, color=col_oval, thick=2
oplot, tnplon, tnplat, color=col_oval, thick=2
oplot, tselon, tselat, color=col_oval, thick=2
oplot, tsplon, tsplat, color=col_oval, thick=2

; Put a dot on the ovals at MLT=12, 24

dot_lng = [tnelon(240), tnplon(240), tselon(240), tsplon(240), $
           tnelon(120), tnplon(120), tselon(120), tsplon(120)]
dot_lat = [tnelat(240), tnplat(240), tselat(240), tsplat(240), $
           tnelat(120), tnplat(120), tselat(120), tsplat(120)]
oplot, psym=4, symsize=.5, dot_lng, dot_lat, color=black

; COLOR-CODED LEGEND, LABELS, ETC.

leg_loc = convert_coord(.02, .95, /norm, /to_device)
line_xpos = make_array(8, /int, value=leg_loc(0))
line_ypos = intarr(8)
if !d.name EQ 'PS' then linewidth=500 else linewidth=12
for p=0,7 do line_ypos(p) = leg_loc(1) - p*linewidth
colors = [col_tags, col_tags, $
          col_fast, col_fast, $
          col_tags, col_tags, $
          col_term, col_oval]
ds = string("260B)              ;" degree symbol
legend = ['LEGEND', $
          '----------', $
;          '[] FAST Position Now', $
          '+ FAST Path', $
          'A Apogee (' + string(format='(F5.1)', apogee_lat) + ds + ' ilat)', $
          'P Perigee(' + string(format='(F5.1)', perigee_lat) + ds + ' ilat)', $
          '- Terminator', $
          '- Auroral Oval']
xyouts, line_xpos, line_ypos, /device, legend, color=colors

; Tell when craft enters and exits shade. Must have shade(0) NE -1
; True transition from shade to sunlight if data_break NE n_datapts-1
; True transition from sunlight to shade if shade(0) NE 0.

if shade(0) NE -1 then begin
    if (!VERSION.RELEASE LE '5.4') then begin
      if data_break NE (n_datapts - 1) $
        then ecl_ex = 'Eclipse Ext: '+(str_sep(time_string(gse_stc.x(data_break)), '/'))(1) $
      else ecl_ex = ''
      if shade(0) NE 0 $
        then ecl_en = 'Eclipse Ent: '+(str_sep(time_string(gse_stc.x(shade(0))), '/'))(1) $
      else ecl_en = ''
    endif else begin
      if data_break NE (n_datapts - 1) $
        then ecl_ex = 'Eclipse Ext: '+(strsplit(time_string(gse_stc.x(data_break)), '/', /EXTRACT))(1) $
      else ecl_ex = ''
      if shade(0) NE 0 $
        then ecl_en = 'Eclipse Ent: '+(strsplit(time_string(gse_stc.x(shade(0))), '/', /EXTRACT))(1) $
      else ecl_en = ''
    endelse
    ecl_label = [ecl_ex, ecl_en]
    leg_loc = convert_coord(.98, .05, /norm, /to_device)
    line_xpos = make_array(2, /int, value=leg_loc(0))
    line_ypos = intarr(2)
    if !d.name EQ 'PS' then linewidth=500 else linewidth=12
    for p=0,1 do line_ypos(p) = leg_loc(1) - p*linewidth
    xyouts, line_xpos, line_ypos, /device, ecl_label, color=col_tags, align=1.0
endif

; Write TMIN and TMAX, or XMARK

if keyword_set(time1) AND keyword_set(time2) then begin
    if (!VERSION.RELEASE LE '5.4') then begin
        t12_label = ['t1 = '+(str_sep(time_string(time1s),'/'))(1), $
                     't2 = '+(str_sep(time_string(time2s),'/'))(1)  ]
    endif else begin
        t12_label = ['t1 = '+(strsplit(time_string(time1s),'/',/EXTRACT))(1), $
                     't2 = '+(strsplit(time_string(time2s),'/',/EXTRACT))(1)  ]
    endelse
    leg_loc = convert_coord(.02, .05, /norm, /to_device)
    line_xpos = make_array(2, /int, value=leg_loc(0))
    line_ypos = intarr(2)
    if !d.name EQ 'PS' then linewidth=500 else linewidth=12
    for p=0,1 do line_ypos(p) = leg_loc(1) - p*linewidth
    xyouts, line_xpos, line_ypos, /device, t12_label, color=col_tags
endif
if keyword_set(xmark) then xyouts, .02, .05, /normal, col=col_tags, $
  'X = '+time_string(xmark)

; Add POLAR spacecraft footprint if requested

if keyword_set(polar) then begin
    fast_flat = double(flat*!dtor)
    fast_flng = double(flng*!dtor)
    fast_times = clock
    plot_po_overlay, tag_color=col_tags, $
      time=clock, FLAT=polar_flat, FLNG=polar_flng
    ;;get_data, 'FLAT', data=po_flat_stc
    ;;get_data, 'FLNG', data=po_flng_stc
    polar_flat = polar_flat*!dtor
    polar_flng = polar_flng*!dtor
    
    ;; Find closest approach of footprints
    
    fx = cos(fast_flat)*cos(fast_flng)
    px = cos(polar_flat)*cos(polar_flng)
    fy = cos(fast_flat)*sin(fast_flng)
    py = cos(polar_flat)*sin(polar_flng)
    fz = sin(fast_flat)
    pz = sin(polar_flat)
    linesep = min(sqrt((fx-px)^2 + (fy-py)^2 + (fz-pz)^2), cji)
    ;;angsep = asin(linesep/2d)*!radeg
    ;;print, 'Geocent Ang Sep: ', time_string(fast_times(cji)), angsep, ds
    
    ;; Label Conjunction
    
    if (!VERSION.RELEASE LE '5.4') then begin
        cnj_time = (str_sep(time_string(fast_times(cji)), '/'))(1)
    endif else begin
        cnj_time = (strsplit(time_string(fast_times(cji)), '/', /EXTRACT))(1)
    endelse

    d_lat = string((fast_flat(cji) - polar_flat(cji))*!radeg, format='(F5.1)')
    d_lng = string((fast_flng(cji) - polar_flng(cji))*!radeg, format='(F5.1)')
    cnj_loc = convert_coord(.02, .17, /normal, /to_device)
    cnj_xpos = make_array(4, /int, value=cnj_loc(0))
    cnj_ypos = intarr(4)
    for p=0,3 do cnj_ypos(p) = cnj_loc(1) - p*linewidth
    cnj_stg = ['CONJUNCTION:', cnj_time, '!4D!3Lat '+d_lat, '!4D!3Lng '+d_lng]
    xyouts, cnj_xpos, cnj_ypos, /device, cnj_stg, color=col_tags
    ;;oplot, [fast_flng(cji),polar_flng(cji)]*!radeg, $
    ;;  [fast_flat(cji),polar_flat(cji)]*!radeg,color=col_oval
endif

; Close postscript file if opened

if keyword_set(psfile) then pclose

; Write to GIF file if desired

if keyword_set(giffile) AND (!d.name EQ 'X' OR !d.name EQ 'Z') then begin
    tvlct, /get, roy, gee, biv
    write_gif, giffile, tvrd(), roy, gee, biv
endif

; Reset to original font system and plot window

!p.font = oldfont
if keyword_set(orig_plot_pos) then !p.position=orig_plot_pos

; Message

print, 'Plot_fa_crossing graphical interface: PFCGUI.PRO'

end
