PRO modlat,image,h,sat,device=device,shrink=shrink
;
; This procedure takes the name of an image file, an assumed emission
; height (km), and a satellite name (strict SpaceTrack format), and
; plots the specified image in geographic coordinates, with the
; trajectory of the specified satellite overplotted on the projected
; image.  The satellite name is optional.  If keyword device is
; set, the plot goes to the specified device, otherwise the default
; device is 'X'.  If keyword shrink is set, the image is shrunk by
; the factor shrink, otherwise the full 256x256 image is used.
;
; Note that variable names that correspond to large arrays having
; outlived their usefulness are reassigned to short strings to save
; memory.
;
t0=SYSTIME(1)
; set 'decimation' factor to speed up processing
IF KEYWORD_SET(shrink) THEN ff=FLOAT(shrink) ELSE ff=1.0
; edge length of 'decimated' image
l=LONG(292L/ff+0.5)
; read in image and subtract dark frame
PRINT,'Read image and subtract dark frame...'
rdkimg,image,hb,i
head=gethd(hb)
cam=head.misc.cam
exp=head.misc.exp_scale*head.exp.exposure/1000.
rdmeandk,cam,exp,dk1
i=(i-dk1)>0
ibord=INTARR(292,292)
ibord(18:273,18:273)=i
dk1='empty'
; 'decimate' the image
IF ff EQ FIX(ff) THEN BEGIN
    bi=REBIN(ibord,l,l)
ENDIF ELSE BEGIN
    bi=CONGRID(ibord,l,l,/INTERP,/MINUS_ONE)
ENDELSE
i='empty'

; identify useful min and max values for byte-scaling
prad=150/ff
dddd=SHIFT(dist(l),l/2,l/2)
wprad=WHERE(dddd LE prad)
base=WHERE(dddd GT prad)
lmaxi=pcentile(bi(wprad),0.999)

; if viewing sky image, take dark 'skirt' as min value
bh=HISTOGRAM(bi(base),MIN=0,BIN=1)
maxbh=MAX(bh,lmaxbh)                ; bottom = most prob. skirt value
bi=BYTSCL(bi,MIN=lmaxbh,MAX=lmaxi)

; plot the 'decimated' image for comparison
loadct,13
whitetop
IF (NOT KEYWORD_SET(device)) OR (STRUPCASE(device) EQ 'X') THEN BEGIN
    WINDOW,2,XSIZE=l,YSIZE=l
    TVSCL,ROTATE(bi,2)
ENDIF
PRINT,SYSTIME(1)-t0,' seconds elapsed...'

; create arrays giving the column and row numbers of all pixels
PRINT,'Make azimuth and zenith angle maps...'
w=REFORM(TEMPORARY(INDGEN(l^2)),l,l)
one2two,w,bi,col,row
w='empty'

; cr2az,...

; parameters of the star fit for camera 0
;c0=130.575/ff
;r0=126.297/ff
;r=203.723/ff
;k=0.626410
;ca=0.229637
;c0=130.302/ff
;r0=126.375/ff
c0=148.302/ff
r0=144.375/ff
r=199.895/ff
k=0.638308
ca=0.227637
; apply the inverse transformation, from col/row to az/za
dc=col-REPLICATE(c0,l,l)
col='empty'
dr=row-REPLICATE(r0,l,l)
row='empty'
rho=SQRT(dc^2+dr^2)<r
phi=ATAN(dc,-dr)
dc='empty'
dr='empty'
xp=COS(phi)
yp=SIN(phi)
phi='empty'
cc=COS(ca)
sc=SIN(ca)
x=xp*cc+yp*sc
y=-xp*sc+yp*cc
z=(ASIN(rho/r)/k)<(89.5*!DTOR)

; zero out pixels of image lying outside field of view, just before
; variable 'rho' disappears
bi(WHERE(rho GT 0.70*r))=0
rho='empty'
az=!RADEG*ATAN(-y,x)
x='empty'
y='empty'
PRINT,SYSTIME(1)-t0,' seconds elapsed...'

; define Rabbit Lake coordinates
PRINT,'Make latitude and longitude maps...'
ree=6378.150D0
rep=6356.7505D0
RLlat=58.233
RLlon=256.33-360.0
sR=SIN(RLlat*!DTOR)
cR=COS(RLlat*!DTOR)
; calculate Earth radius at Rabbit Lake
RLrad=SQRT((ree*ree*rep*rep)/(ree*ree*sR*sR+rep*rep*cR*cR))
; angular distance between R.L. and observed point, user-specified
; altitude
gcdist=z-ASIN((RLrad*SIN(z))/(RLrad+h))
z='empty'
; transform image coordinates from (az,za) to (longitude,latitude)
rb2ll,RLlon,RLlat,gcdist,az,imglon,imglat
az='empty'
gcdist='empty'
; make imglat and imglon single-precision to save space
imglat=FLOAT(imglat)
imglon=FLOAT(imglon)-REPLICATE(360.,l,l)
; restrict plotted coordinates to points inside lat/lon grid
clat=58.25
clon=-103.5
dlatdh=2.25D0/110.
dlondh=5.5D0/110.
wedge=FLOAT(clon-dlondh*h)
eedge=FLOAT(clon+dlondh*h)
nedge=FLOAT(clat+dlatdh*h)
sedge=FLOAT(clat-dlatdh*h)
PRINT,wedge,eedge,sedge,nedge $
     ,FORMAT='("Plot limits: ",F7.2," - ",F7.2," longitude, ",F6.2," - ",F6.2," latitude")'
; set up device
IF KEYWORD_SET(device) AND (STRUPCASE(device) EQ 'PS') THEN BEGIN
    psopen,'tmp.ps',/landscape,/color
    DEVICE,BITS_PER_PIXEL=8
    loadct,13
ENDIF ELSE WINDOW,0
; create map title
PRINT,SYSTIME(1)-t0,' seconds elapsed...'
; create title for map
date=INTARR(3)
FOR j=0,2 DO date(j)=head.misc.dt.(j)
dst=ymd2date(date(2),date(1),date(0),format='d$ N$ Y$')
time=INTARR(3)
FOR j=0,2 DO time(j)=head.misc.tm.(j)
ut=time(2)+60.*(time(1)+60.*time(0))
tst=strsec(ut)
PRINT,'Enter format string for camera wavelength'
wvst=''
READ,wvst
mtitle=dst+' '+tst+' '+wvst
; plot map
PRINT,'MAP_SET...'
MAP_SET,RLlat,RLlon,/ORTHOGRAPHIC,LIMIT=[sedge,wedge,nedge,eedge] $
       ,TITLE=mtitle,/NOBORDER,XMARGIN=[3,3],YMARGIN=[3,3] $
       ,CHARSIZE=0.75,/LABEL,LATALIGN=0.0,LATDEL=0.5 $
       ,LATLAB=-108,LONLAB=56,LONALIGN=1.0,LONDEL=1.0
; interpolate Polar Camera image to a rectangular grid in latitude
; and longitude
PRINT,SYSTIME(1)-t0,' seconds elapsed...'
PRINT,'TRIANGULATE...'
TRIANGULATE,imglon(wprad),imglat(wprad),tr,b
PRINT,SYSTIME(1)-t0,' seconds elapsed...'
PRINT,'TRIGRID...'
lonres=0.05*(h/110.)
latres=0.025*(h/110.)
gridi=FIX(TRIGRID(imglon(wprad),imglat(wprad),bi(wprad),tr $
		 ,[lonres,latres],[wedge,sedge,eedge,nedge]))
wprad='empty'
imglon='empty'
imglat='empty'
; display interpolated image
PRINT,SYSTIME(1)-t0,' seconds elapsed...'
PRINT,'MAP_IMAGE...'
mapi=FIX(MAP_IMAGE(gridi,sx,sy,xz,yz,LATMIN=sedge,LATMAX=nedge $
                  ,LONMIN=wedge,LONMAX=eedge,MISSING=0,COMPRESS=2))
!ORDER=0
PRINT,SYSTIME(1)-t0,' seconds elapsed...'
PRINT,'Plot the warped image...'
TVSCL,mapi,sx,sy
; overplot the latitude/longitude grid
MAP_GRID,/LABEL,LATALIGN=0.0,LATDEL=0.5,LATLAB=-108 $
	,LONLAB=56,LONALIGN=1.0,LONDEL=1.0,CHARSIZE=0.75
; show Rabbit Lake position
prov
PLOTS,rllon,rllat,PSYM=1
XYOUTS,rllon-0.4,rllat,'Rabbit Lake',ALIGNMENT=1.
PRINT,SYSTIME(1)-t0,' seconds elapsed...'
; overplot the horizon at Rabbit Lake
PRINT,'Overplot zenith angle circles...'
nza=4
hang=!DTOR*15.*(FINDGEN(nza)+1)
gcrad=hang-ASIN((RLrad*SIN(hang))/(RLrad+h))
naz=37
gaz=FINDGEN(naz)*360./(naz-1)
hlon=FLTARR(naz,nza)
hlat=hlon
FOR j=0,naz-1 DO BEGIN
    rb2ll,RLlon,RLlat,gcrad,gaz(j),alon,alat
    hlon(j,*)=alon
    hlat(j,*)=alat
ENDFOR
FOR j=0,nza-1 DO BEGIN
    PLOTS,hlon(*,j)-360.,hlat(*,j)
    XYOUTS,hlon(0,j),hlat(0,j)+0.05,ALIGNMENT=0.5,CHARSIZE=0.75 $
          ,STRTRIM(FIX(!RADEG*hang(j)+.001),2)+'!9%!x'
ENDFOR
t1=SYSTIME(1)
PRINT,t1-t0,' seconds elapsed...'
; now get the footprint trajectory of the desired satellite and
; overplot it
IF KEYWORD_SET(sat) THEN BEGIN
    PRINT,'Overplot satellite trajectory...'
    t2=SYSTIME(1)
    nut=FIX((ut+30.)/60.)*60.
    sechms,nut-300.,h,m,s,sh,sm,ss
    stime=(h*100+m)*100L+s
    sechms,nut+300.,h,m,s,sh,sm,ss
    etime=(h*100+m)*100L+s
    indate=STRING(date(2),date(1),date(0),FORMAT='(I2.2,I2.2,I2.2)')
    plotorb,sat,indate,stime,etime,0.25,'x' $
	   ,/noplot,ptime,plat,plon
    USERSYM,[1,0,-1,0],[0,1,0,-1],/FILL
    PLOTS,plon,plat,PSYM=-8,NOCLIP=0
    atmin=WHERE(LONG(ptime) MOD 100 EQ 0)
    XYOUTS,plon(atmin),plat(atmin)+0.1,STRTRIM(LONG(ptime(atmin)),2) $
	  ,ALIGNMENT=0.5,CHARSIZE=0.75
    itime=strsec(ut+30.)
    XYOUTS,RLlon,sedge-(nedge-sedge)/22.86,itime,ALIGNMENT=0.5
    PRINT,SYSTIME(1)-t2+t1-t0,' seconds elapsed...'
ENDIF
; end plot
IF !D.NAME EQ 'PS' THEN BEGIN
    psclose
    SPAWN,'rmps3 tmp.ps' 
ENDIF
END