;+
; NAME:
;	
; PURPOSE:
;	
; CATEGORY:
;	
; CALLING SEQUENCE:
;	MAP_GEO, HEADBYTES, IMAGE, HEIGHT, PIMG $
;	       , NEW_CAM = NEW_CAM, VRBG = VRBG, FULLSCALE = FULLSCALE $
;	       , COLOR = COLOR, ZOOM = ZOOM, ANNOTATE = ANNOTATE $
;              , GRID = GRID, ROT = ROT, POLE = POLE, SAT = SAT
; INPUTS:
;	HEADBYTES: The header bytes from the image to be transformed.
;       IMAGE:  The image to be transformed.  The procedure is to
;		read the images from a .ray or .bgs file.
;	HEIGHT:	The assumed altitude of emission (km).
; OPTIONAL INPUTS:
;	None.
; KEYWORD PARAMETERS:
;       /NEW_CAM: Set this keyword if MAP_GEO has already been
;               called in this session of IDL, but with a different
;               camera/filter combination.  This keyword forces
;               execution of code that would otherwise be bypassed.
;       VRBG:   Gives the zenith brightness (R) of the assumed van
;               Rhijn background emission layer to be subtracted
;		from the images.
;       FULLSCALE: A 2-element array giving the minimum and maximum
;               brightnesses (R) to be included in the color table.
;               If not specified, each image is autoscaled to the
;		0.1 and 99.9 percentiles for display.
;       COLOR:  The element of the color table to use in annotating
;		the images.
;       /ANNOTATE: If set, the image is annotated with date/time/
;		camera/filter/exposure time information, and color 
;		scale limits.
;       GRID:   Takes STRING value.  If it contains [cC], the image
;               is decorated with coastlines.  If it contains [gG],
;		the image is decorated with a coordinate grid.
;       /ROT:   If set, the image is rotated to put Eureka in the
;               correct orientation w.r.t. the PACE geomagnetic
;               coordinate pole at 'height'.  Eureka is always
;		located at the center of the image, unless the POLE
;		keyword is set; see below.
;       /POLE:  If set, this keyword forces the PACE geomagnetic
;               pole at the specified HEIGHT to be the center of
;               the map plot.  This keyword presupposes that ROT
;               is set, and in fact forces it to be set if the user
;		has not already done so.
;       SAT:    Contains one or more names (in standard 2-line
;               element form) of satellites whose footprints at
;               the specified emission altitude are to be overplotted
;               onto the mapped image.  The footprint across the
;               whole image is diplayed, and the satellite positions
;               at the start and end of the image exposure time
;		are highlighted.
;	/ZOOM:	This keyword has been disabled.
; OUTPUTS:
;       PIMG:   A 383x383 BYTE array containing the mapped image,
;		along with annotation and decoration if selected.
; OPTIONAL OUTPUTS:
;	None.
; COMMON BLOCKS:
;       COORDS: COORDHT duplicates HEIGHT, and flags whether the
;                       present call of MAP_GEO is the initial one;
;               PLAT    is used only to pass information on plot
;			resolution back to the calling program;
;		PLONG	isn't used at all in MAP_GEO;
;               GLAT    is a rectangular grid of latitudes to which
;                       IMAGE is initially mapped to provide input
;			for MAP_IMAGE;
;		GLONG	is the corresponding grid of longitudes;
;		ASIZE	is the array descriptor returned by SIZE(GLAT);
;               AZ      is the array of azimuths for the coordinates
;                       (GLONG,GLAT) seen at altitude HEIGHT from
;			the Polar Camera location;
;		ZA	is the corresponding array of zenith angles;
;		COL	is the corresponding array of image columns;
;		ROW	is the corresponding array of image rows;
;               VRLAYER is the synthesized van Rhijn layer (of unit
;                       zenith brightness) for the given camera at
;			the assumed HEIGHT;
;		POCALAT	is the latitude of the Polar Camera;
;		POCALONG is the corresponding longitude;
;		POCAHT	is the altitude (km) of the Polar Camera;
;		HOR_SIZE is the width in pixels of the mapped image;
;		VERTSIZE is the corresponding height;
;               EUR_MLT is the PACE Magnetic Local Time at Eureka
;			at the time of the current IMAGE;
;               PACEPOLECOORDS is a catch-all array, the contents
;			of which are as follows:
;               0:      geographic latitude of the PACE pole at
;			'height';
;               1:      corresponding geographic longitude;
;		2:	'height';
;		3:	azimuth of Eureka from the PACE pole;
;               4:      CW angle by which to rotate image so Eureka
;                       is at correct MLT w.r.t. PACE pole; e.g.,
;                       if MLT is 1.0 h, Eureka is at azimuth 180
;                       - 1*15 = 165 from PACE pole in rotated plot;
;               5:      X distance (in DEVICE units) across a circle
;                       centered at the Polar Camera and spanning
;                       60 deg zenith angle at all azimuths, at
;			'height';
;               6:      corresponding Y distance;
;               7-10:   longitudes at 60 deg zenith angle and four
;                       cardinal directions from Eureka at 'height';
;		11-14:  corresponding latitudes;
;
;               UNROTLIMIT is an 8-element array of latitudes and
;                       longitudes of the bounds of the coordinate
;                       system created by MAP_SET if the image is
;                       plotted with North at the top; it is in
;                       the format required for the LIMIT keyword
;                       to MAP_SET, and is used to determine the
;                       appropriate map limits when the image is
;			rotated (see /ROT above);
;	DATE:	YEAR	is the (UT) year of the current image;
;		MONTH	is the corresponding month;
;		DAY	is the corresponding day;
; SIDE EFFECTS:
;       A pixmap is created to hold the mapped image and annotation.
;       The contents of this pixmap are read back using TVRD to
;       generate the image returned to the calling program.  An
;       information file called map_geo.parms is created in ~cnsr3
;	to assist in debugging the code.
; RESTRICTIONS:
;	None known.  Give me time.
; PROCEDURE:
;	Later.
; EXAMPLE:
;	See VIDEO.PRO and POCAFLIC.PRO for a couple of different uses.
; SEE ALSO:
;	MAP_PGM.PRO
; MODIFICATION HISTORY:
; 	Written by:	D P Steele, ISR, Feb.-Mar. 1995.
;-

PRO map_geo,headbytes,image,height,pimg $
	   ,new_cam=new_cam,vRbg=vRbg,fullscale=fullscale $
	   ,color=color,zoom=zoom,annotate=annotate,grid=grid $
	   ,rot=rot,pole=pole,sat=sat

COMMON coords,coordht,plat,plong,glat,glong,asize,az,za,col,row $
	     ,vRlayer $
	     ,pocalat,pocalong,pocaht,hor_size,vertsize,eur_mlt $
	     ,PACEpolecoords,unrotlimit
COMMON date,year,month,day

; Save the current graphics device and set the device to the Z
; buffer.

olddev=!D.NAME
TVLCT,rx,gx,bx,/GET	; save the X-windows RGB color table
oldord=!ORDER
!ORDER=0

SET_PLOT,'Z'
DEVICE,Z_BUFFERING=0
loadct,0,/silent
IF N_ELEMENTS(color) EQ 0 THEN color=!D.N_COLORS-1

; If this is the first invocation, read in the coordinate arrays

IF NOT KEYWORD_SET(coordht) THEN coordht=0
IF (height NE coordht) THEN BEGIN
    PRINT,'Generating geographic coordinates for '+STRTRIM(height,2)+' km'
    coordht=FIX(ROUND(height))
    maxza=70.
    header=gethd(headbytes)
    loc=pocaloc(header.misc.dt.year,header.misc.dt.month,header.misc.dt.day)
    pocalat=loc(0)
    pocalong=loc(1)-360.
    pocaht=loc(2)

    lowlat=72.9		; These limits encompass the range of latitudes
    highlat=88.2	; and longitudes covered by images from both
    lowlong=234.1-360.	; cameras when mapped to 300 km altitude.  They
    highlong=320.1-360.	; should be adequate for all purposes.

    ; These empirically determined values are close to the actual
    ; values set by the routine.
    IF KEYWORD_SET(pole) THEN mapres=4.2 ELSE mapres=3.5

    rb2ll,pocalong,pocalat,mapres/(relc(pocalat)+height),0.,Nlong,Nlat
    dellat=Nlat-pocalat
    rb2ll,pocalong,pocalat,mapres/(relc(pocalat)+height),90.,Elong,Elat
    IF Elong GT 180 THEN Elong=Elong-360.
    dellong=Elong-pocalong

    neast=FIX((pocalong-lowlong)/dellong)
    nwest=FIX((highlong-pocalong)/dellong)
    glong=(FINDGEN(neast+1+nwest)-neast)*dellong+pocalong

    nsouth=FIX((pocalat-lowlat)/dellat)
    nnorth=FIX((highlat-pocalat)/dellat)
    glat=(FINDGEN(nsouth+1+nnorth)-nsouth)*dellat+pocalat

    tmp=glat
    glat=REPLICATE(1.,N_ELEMENTS(glong))#tmp
    glong=glong#REPLICATE(1.,1,N_ELEMENTS(tmp))
    tmp=0B

    asize=SIZE(glat)

; The following is a kludge to make MAP_GEO look like MAP_PGM for
; determining latitude resolution of the mapped image.

    plat=FLTARR(asize(1)/2+1,asize(2))

; Center the map on what location? (Part 1)
    IF KEYWORD_SET(pole) THEN rot=1B

; Determine geographic coordinates of PACE pole at given altitude.
    IF KEYWORD_SET(rot) THEN BEGIN
	pacepole,height,gltp,glnp
        PACEpolecoords=[glnp,gltp,height]
	ENDIF ELSE PACEpolecoords=[280.73,80.57,300.]

; Center the map on what location? (Part 2)
    IF KEYWORD_SET(pole) THEN BEGIN
        ctrlong=PACEpolecoords(0)
        ctrlat=PACEpolecoords(1)
        ENDIF ELSE BEGIN
        ctrlong=pocalong
        ctrlat=pocalat
        ENDELSE
    ctrlong=ctrlong-(ctrlong GT 180.)*360.

; Define appropriate limits for map.  The following code sets up
; the 8-element LIMIT vector used by the MAP_SET command.  See IDL
; Reference Guide, p. 1-149, for details.

    lat1=MAX(glat,MIN=lat3)
    lon0=MIN(glong,MAX=lon2)
    unrotlimit=[ctrlat,lon0,lat1,ctrlong,ctrlat,lon2,lat3,ctrlong]
    tmp=REFORM(unrotlimit,2,4)
    unrotlons=REFORM(tmp(1,*))
    unrotlats=REFORM(tmp(0,*))
    ll2rb,ctrlong,ctrlat,unrotlons,unrotlats,dist,azi

    pdist=dist*COS((azi-SHIFT(FINDGEN(4)*90.,1))*!DTOR)
    pdist=REPLICATE(mean(pdist)<(6./!RADEG),N_ELEMENTS(pdist))
    rb2ll,ctrlong,ctrlat,pdist,SHIFT(FINDGEN(4)*90.,1),edgelon,edgelat
    edgelon=edgelon-(edgelon GT 180)*360.
    unrotlimit=REFORM([TRANSPOSE(edgelat),TRANSPOSE(edgelon)],8)

    hdist=(relc(ctrlat)+height)*pdist
    perpdist=hdist*COS((azi-SHIFT(FINDGEN(4)*90.,1))*!DTOR)

    xpixels=ROUND((perpdist(0)+perpdist(2))/mapres)
    ypixels=ROUND((perpdist(1)+perpdist(3))/mapres)
    xpixmax=383		; choose this size for compatibility with 
    ypixmax=383		; images mapped to PACE geomagnetic coords.
    grow=(FLOAT(xpixmax)/FLOAT(xpixels))<(FLOAT(ypixmax)/FLOAT(ypixels))
    PRINT,mapres/grow $
	 ,FORMAT='("Actual map resolution is ",F3.1," km at Polar Camera.")'

; Determine amount of rotation needed to put Eureka at correct 'MLT'
; w.r.t. PACE pole at specified altitude

    IF KEYWORD_SET(rot) THEN BEGIN
	ll2rb,glnp,gltp,pocalong,pocalat,Edist,Eazi	; get azimuth of
	Eazi=Eazi(0)					; Eureka from 
	year=header.misc.dt.year			; PACE pole
	month=header.misc.dt.month
	day=header.misc.dt.day
        hour=header.misc.tm.hour
	minute=header.misc.tm.minute
	second=header.misc.tm.second
	mconvert,year,month,day,hour,minute,second $
		,pocalat,pocalong,height,1,mlat,mlong,eur_mlt
        rotangle=(15.*(12.-eur_mlt)-Eazi+720.) MOD 360.
	ENDIF ELSE BEGIN
	ll2rb,280.73,80.57,pocalong,pocalat,Edist,Eazi
	Eazi=Eazi(0)
	rotangle=0
	eur_mlt=-99.0	; dummy value
	ENDELSE
    PACEpolecoords=[PACEpolecoords,Eazi,rotangle]

; Rotate map limits to be consistent with rotangle.

    bazi=(SHIFT(FINDGEN(4)*90.,1)-rotangle+360) MOD 360.
    rb2ll,ctrlong,ctrlat,pdist,bazi,rotlons,rotlats
    rotlons=rotlons-(rotlons GT 180)*360.
    rotlimit=REFORM([TRANSPOSE(rotlats),TRANSPOSE(rotlons)],8)

; Adjust window dimensions until aspect ratio is close enough to 1.0

; Which window dimension needs adjustment?

    modix=(xpixels GT ypixels)
    winsize=[ROUND(grow*xpixels),ROUND(grow*ypixels)]
    fasp=1.0

; Define the reference points for distance comparison.

    dist=za2dist(60.,height,/degrees,radius=relc(ctrlat)) $
        /(relc(ctrlat)+height)
    rb2ll,ctrlong,ctrlat,dist,FINDGEN(4)*90.-rotangle,blong,blat
    blong=blong-(blong GT 180)*360.

    REPEAT BEGIN	; until map aspect ratio = 1.00
	winsize(modix)=ROUND(fasp*winsize(modix))

;       WINDOW,/FREE,XSIZE=winsize(0),YSIZE=winsize(1) $
;             ,PIXMAP=(N_PARAMS() EQ 4)
;       dwin=!D.WINDOW

        DEVICE,SET_RESOLUTION=winsize

        MAP_SET,ctrlat,ctrlong,rotangle,/STEREOGRAPHIC,/NOBORDER $
               ,LIMIT=rotlimit,XMARGIN=[0,0],YMARGIN=[0,0],/WHOLE_MAP

; Check the aspect ratio of the resulting plot: will a circular
; field of view look circular?

        devco=CONVERT_COORD(blong,blat,/DATA,/TO_DEVICE)
        xdiam=devco(0,1)-devco(0,3)
        ydiam=devco(1,0)-devco(1,2)
	CASE modix OF
	 0: fasp=ydiam/xdiam
	 1: fasp=xdiam/ydiam
	    ENDCASE
;       WDELETE,dwin

        DEVICE,/CLOSE			; close the Z-buffer

	ENDREP UNTIL (0.995 LT fasp) AND (fasp LT 1.005)

    hor_size=winsize(0)
    vertsize=winsize(1)
    PACEpolecoords=[PACEpolecoords,xdiam,ydiam,blong,blat]
    PRINT,"Window is "+STRTRIM(hor_size,2)+" by "+STRTRIM(vertsize,2)
    PRINT,fasp $
	 ,FORMAT='("Map aspect ratio is ",F5.3," (hor/vert)")'
    chsz=1.

; Convert the geographic coordinate arrays to look directions from
; Eureka

    PRINT,'Converting coordinates to look directions'
    geo2az,glat,glong,REPLICATE(height,asize(1),asize(2)) $
	  ,pocalat,pocalong,pocaht,az,za

; If van Rhijn correction is desired, generate the van Rhijn
; background image to be subtracted from all images

    IF vRbg GT 0 THEN BEGIN
	vRlayer=vanRhijn(za,height)
	PRINT,vRbg $
	     ,FORMAT='("van Rhijn layer brightness ",F4.0," R")'
	ENDIF ELSE vRlayer=INTARR(256,256)

; Initialize the file to hold the mapping parameters for later
; inspection.

    OPENW,parm,'map_geo.parms',/GET_LUN
    PRINTF,parm,'Day HH:MM:SS RotAngl Xsz Ysz Xst Yst Xdi    ' $
	       +'Ydi    Xct    Yct     Max'
    FREE_LUN,parm

    ENDIF

; Extract information from the image header

header=gethd(headbytes)
year=header.misc.dt.year
month=header.misc.dt.month
day=header.misc.dt.day
hour=header.misc.tm.hour
minute=header.misc.tm.minute
second=header.misc.tm.second
PCsite=PoCasite(year,month,day)
camera=header.misc.cam
filter=header.misc.fw
exptime=header.misc.exp_scale*header.exp.exposure/1000

chsz=1.
IF (N_ELEMENTS(col) EQ 0) OR KEYWORD_SET(new_cam) THEN BEGIN

; Transform the look directions to image coordinates (column,row)

    PRINT,'Converting look directions to column/row'
    az2cr,year,month,day,camera,filter,az,za,col,row $
	 ,/small,/deg,site=PCsite

; Avoid impossible or useless image coordinates

    crej=WHERE((col LE 0) OR (255 LE col),ncrej)
    rrej=WHERE((row LE 0) OR (255 LE row),nrrej)
    IF ncrej GT 0 THEN BEGIN
	col(crej)=0
	row(crej)=0
	ENDIF
    IF nrrej GT 0 THEN BEGIN
	col(rrej)=0
	row(rrej)=0
	ENDIF
    ENDIF

; Sample the mapped image from the observed image

IF N_PARAMS() LT 4 THEN PRINT,'Sampling image'
map=FIX(ROUND((image(col,row)-vRbg*vRlayer(col,row))>0))
gt0=WHERE(map,ngt0)
IF KEYWORD_SET(fullscale)  $
    THEN lims=fullscale $
    ELSE lims=pcentile(map(gt0),[0.001,0.999])

; Ensure pixel above Eureka is at full scale brightness

PCaz=[0.]
PCza=PCaz
az2cr,year,month,day,camera,filter,PCaz,PCza,PCcol,PCrow $
     ,/small,/deg
map(ROUND(PCcol(0)),ROUND(PCrow(0)))=lims(1)

; Determine amount of rotation needed to put Eureka at correct 'MLT'
; w.r.t. PACE pole at specified altitude

IF KEYWORD_SET(rot) THEN BEGIN
    year=header.misc.dt.year
    month=header.misc.dt.month
    day=header.misc.dt.day
    hour=header.misc.tm.hour
    minute=header.misc.tm.minute
    second=header.misc.tm.second
    mconvert,year,month,day,hour,minute,second $
            ,pocalat,pocalong,height,1,mlat,mlong,eur_mlt
    Eazi=PACEpolecoords(3)
    rotangle=(15.*(12.-eur_mlt)-Eazi+720.) MOD 360.
    ENDIF ELSE BEGIN
    rotangle=0
    eur_mlt=-99.0
    ENDELSE
PACEpolecoords(4)=rotangle

; Center the map on what location? (Part 2)
IF KEYWORD_SET(pole) THEN BEGIN
    ctrlong=PACEpolecoords(0)
    ctrlat=PACEpolecoords(1)
    ENDIF ELSE BEGIN
    ctrlong=pocalong
    ctrlat=pocalat
    ENDELSE
ctrlong=ctrlong-(ctrlong GT 180.)*360.

; Rotate map limits to be consistent with rotangle.

tmp=REFORM(unrotlimit,2,4)
unrotlons=REFORM(tmp(1,*))
unrotlats=REFORM(tmp(0,*))
ll2rb,ctrlong,ctrlat,unrotlons,unrotlats,dist,azi
dist=REPLICATE(mean(dist)<(6./!RADEG),N_ELEMENTS(dist))
bazi=(azi-rotangle+360.) MOD 360.
rb2ll,ctrlong,ctrlat,dist,bazi,rotlons,rotlats
rotlons=rotlons-(rotlons GT 180)*360.
rotlimit=REFORM([TRANSPOSE(rotlats),TRANSPOSE(rotlons)],8)

; Transform image into appropriate map system and display, with
; coordinate grid and continental outlines.  If desired, rotate
; image and map to give the correct MLT orientation.

;WINDOW,/FREE,XSIZE=hor_size,YSIZE=vertsize,PIXMAP=(N_PARAMS() EQ 4)
;dwin=!D.WINDOW

DEVICE,SET_RESOLUTION=[hor_size,vertsize]

lat1=MAX(glat,MIN=lat3)
lon0=MIN(glong,MAX=lon2)
MAP_SET,ctrlat,ctrlong,rotangle,/STEREOGRAPHIC,/NOBORDER $
       ,LIMIT=rotlimit,XMARGIN=[0,0],YMARGIN=[0,0],/WHOLE_MAP

; Is the aspect ratio and image scaling preserved?

dist=za2dist(60.,height,/degrees,radius=relc(ctrlat)) $
    /(relc(ctrlat)+height)
rb2ll,ctrlong,ctrlat,dist,FINDGEN(4)*90.-rotangle,blong,blat
blong=blong-(blong GT 180)*360.
devco=CONVERT_COORD(blong,blat,/DATA,/TO_DEVICE)
xdiam=devco(0,1)-devco(0,3)
ydiam=devco(1,0)-devco(1,2)
PACEpolecoords(5)=xdiam		; keep current values of these
PACEpolecoords(6)=ydiam		; parameters for access by
PACEpolecoords(7:10)=blong	; calling program
PACEpolecoords(11:14)=blat

; Now cobble up a pseudo-plat for use by the calling routine if it
; wants to.

dist=ACOS(SIN(blat(0)*!DTOR)*SIN(blat(2)*!DTOR) $
	 +COS(blat(0)*!DTOR)*COS(blat(2)*!DTOR) $
	 *COS((blong(0)-blong(2))*!DTOR))*!RADEG
del=dist/ydiam
yd2=FLOOR(ydiam/2.)
plat(asize(1)/2,0)=(FINDGEN(1,asize(2))-asize(2)/2)*del+ctrlat

; Now map the image into the coordinate system.

mmap=MAP_IMAGE(map,sx,sy,xs,ys $
	      ,LATMIN=lat3,LATMAX=lat1,LONMIN=lon0,LONMAX=lon2,COMPRESS=1)
mmapmax=MAX(mmap)
mmap(0,0)=lims(1)

; Crop by hand - this shouldn't be necessary!
IF (sx+xs GT hor_size) THEN mmap=mmap(0:hor_size-sx-1,*)
IF (sy+ys GT vertsize) THEN mmap=mmap(*,0:vertsize-sy-1)

TVSCL,mmap<lims(1),sx,sy
PLOTS,ctrlong,ctrlat,PSYM=1,COLOR=0.7*!D.N_COLORS
IF ctrlong EQ pocalong $
    THEN PLOTS,PACEpolecoords(0),PACEpolecoords(1),PSYM=3,COLOR=!D.N_COLORS-1 $
    ELSE PLOTS,pocalong,pocalat,PSYM=3,COLOR=!D.N_COLORS-1

; Locate Resolute Bay and Qaanaaq as well

PLOTS,-94.88,74.78,PSYM=7,COLOR=!D.N_COLORS-1	; Resolute
PLOTS,-69.20,77.58,PSYM=7,COLOR=!D.N_COLORS-1   ; Qaanaaq

; Write mapping parameters to a file for later inspection

ctrdevco=CONVERT_COORD(ctrlong,ctrlat,/DATA,/TO_DEVICE)
outs=STRING(day,hour,minute,second,rotangle,xs,ys,sx,sy,xdiam,ydiam $
	   ,ctrdevco(0),ctrdevco(1),mmapmax $
	   ,FORMAT='(2I3,2(":",I2.2),F8.2,4I4,4F7.2,I5)')
OPENU,parm,'map_geo.parms',/GET_LUN,/APPEND
PRINTF,parm,outs
FREE_LUN,parm

;MAP_CONTINENTS,COLOR=0.7*!D.N_COLORS

; Decorate rotated image with coastlines and coordinate grids

IF N_ELEMENTS(grid) GT 0 THEN BEGIN

; Use CIA World Map database (sampled at 10%) to plot continents,
; islands, and lakes, as well as political boundaries, onto map.
    map_dir='/jasper/cnsr3_data1/CIA_World_Map/maps/sol_2/'
;   t0=SYSTIME(1)
    IF STRPOS(STRUPCASE(grid),'C') NE -1 THEN BEGIN
        ndx=map_lndx(map_dir+'cil.ndx')
        map_plot,map_dir+'cil.dat',mmndx=ndx $	; plot continents, ...
                ,color=0.7*!D.N_COLORS
        ENDIF

    IF STRPOS(STRUPCASE(grid),'G') NE -1 $
        THEN mmap_grid,latdel=2,londel=10,/nice,color=color $
                      ,/label,latlab=-70.,lonlab=83.8,latalign=1 $
		      ,lonalign=0.5
    ENDIF

; Annotate image with labels

IF KEYWORD_SET(annotate) THEN BEGIN

    XYOUTS,sx+0.5*xs,sy,STRING(height,FORMAT='(F4.0," km")') $
          ,ALIGNMENT=0.0,/DEVICE,CHARSIZE=chsz,COLOR=color

    IF KEYWORD_SET(fullscale) THEN BEGIN
        loval=fullscale(0)
        hival=fullscale(1)
	ENDIF ELSE BEGIN
	loval=lims(0)
	hival=lims(1)
	ENDELSE
    labelimg,headbytes,hival,loval,sx,sy,xs,ys,chsz $
            ,color=color,/camfilt,/date,/sun,rot=rot $
	    ,rang=rotangle
    ENDIF

; Overplot a satellite footprint?

IF KEYWORD_SET(sat) THEN BEGIN

    ; First see whether the specified satellite is in the neighbourhood
    ; within 10 minutes of the current time

    satdate=((year MOD 100)*100L+month)*100L+day
    openut=(hour*60L+minute)*60L-600
    closeut=openut+1260
    sechms,openut,openh,openm,opens
    opentime=(openh*100L+openm)*100L+opens
    sechms,closeut,closeh,closem,closes
    closetime=(closeh*100L+closem)*100L+closes
    plotorb,sat,satdate,opentime,closetime,1,otime,olat,olong,oht $
	   ,type=4,/noplot,fpalt=height

    ; If any point on the footprint is within the map, then show
    ; the footprint for 5 minutes either side of closest approach
 
    inmap=WHERE((MIN(rotlons) LE olong) AND (olong LE MAX(rotlons)) AND $
 		(MIN(rotlats) LE olat) AND (olat LE MAX(rotlats)),ninmap)
 
    IF ninmap GE 1 THEN BEGIN

        ctime=otime(inmap(ninmap/2))
	centut=sechhmmss(ctime)
	openut=centut-300
	closeut=openut+660
        sechms,openut,openh,openm,opens
	opentime=(openh*100L+openm)*100L+opens
	sechms,closeut,closeh,closem,closes
	closetime=(closeh*100L+closem)*100L+closes
	plotorb,sat,satdate,opentime,closetime,0.25,otime,olat,olong,oht $
	       ,type=4,/noplot,fpalt=height

	; Plot the footprint at 1-minute intervals and show time labels


        inmap=WHERE((MIN(rotlons) LE olong) AND (olong LE MAX(rotlons)) AND $
		    (MIN(rotlats) LE olat) AND (olat LE MAX(rotlats)),ninmap)
        otime=otime(inmap)
        olong=olong(inmap)
	olat=olat(inmap)
	OPLOT,olong,olat,COLOR=color	; plot all 15-s segments in map

	; Now select points at whole UT minutes

	wholemin=WHERE((ROUND(otime) MOD 100) EQ 0,nwholemin)
	IF nwholemin GT 0 THEN BEGIN
	    otime=otime(wholemin)
	    olat=olat(wholemin)
	    olong=olong(wholemin)
	    OPLOT,olong,olat,PSYM=4,COLOR=color	; plot whole minutes
	    outs=sechhmmss(otime)
	    dvpos=CONVERT_COORD(olong,olat,/DATA,/TO_DEVICE)
	    satnum=FIX(STRMID(sat,STRLEN(sat)-1,1))
	    algn=(2-((satnum MOD 3) EQ 0))/2.	; => 0.5 for 3 (F8), 6 (F11)
	    XYOUTS,dvpos(0,*) $
	          ,dvpos(1,*)+0.5*!D.Y_CH_SIZE $
	          ,STRMID(strsec(outs,/hour),0,5) $
	          ,/DEVICE,ALIGNMENT=algn,COLOR=color

	; Put the satellite name somewhere near the footprint

	    algn=1.0-algn
            IF STRPOS(sat,'DMSP') NE -1 THEN BEGIN
	        series_no=FIX(STRMID(sat,STRLEN(sat)-1,1))
	        satname=STRING(series_no+5,FORMAT='("DMSP F",I0)')
	        ENDIF ELSE satname=sat
	    xpos=dvpos(0,0)
	    up=(dvpos(1,1) GT dvpos(1,0))	; true if track has + slope
	    lend=xpos-algn*STRLEN(satname)*!D.X_CH_SIZE
	    IF lend LT !D.X_CH_SIZE THEN xpos=xpos+!D.X_CH_SIZE-lend
            XYOUTS,xpos $
	          ,dvpos(1,0)-1.5*!D.Y_CH_SIZE $
	          ,satname $
	          ,/DEVICE,ALIGNMENT=algn,COLOR=color

	    ; Now zero in on the footprint during this image

	    openut=(hour*60L+minute)*60L+second
	    closeut=openut+exptime
	    sechms,openut,openh,openm,opens
	    opentime=(openh*100L+openm)*100L+opens
	    sechms,closeut,closeh,closem,closes
	    closetime=(closeh*100L+closem)*100L+closes
	    plotorb,sat,satdate,opentime,closetime,0.1,otime,olat,olong,oht $
	           ,type=4,/noplot,fpalt=height

            inmap=WHERE((MIN(rotlons) LE olong) AND (olong LE MAX(rotlons)) AND $
		        (MIN(rotlats) LE olat) AND (olat LE MAX(rotlats)),ninmap)
	    IF ninmap GT 0 THEN BEGIN
                olong=olong(inmap)
	        olat=olat(inmap)
                OPLOT,olong,olat,THICK=3,COLOR=color
	        ENDIF

	    ENDIF	; (nwholemin GT 0)

	ENDIF	; (MIN(dist,loc) LT za2dist...)

    ENDIF	; (KEYWORD_SET(sat))

EMPTY

; If the mapped image is to be returned, get it from the display
; subsystem

IF (N_PARAMS() EQ 4) THEN BEGIN
    pimg=FIX(TVRD()) 
    pimg=extrac(TEMPORARY(pimg),0,0,383,383)
    ENDIF ELSE k=GET_KBRD(1)
;WDELETE,dwin

DEVICE,/CLOSE

; Restore the graphics device.

SET_PLOT,olddev
TVLCT,rx,gx,bx
!ORDER=oldord

RETURN
END