PRO tolatlon,image,h,sat,device=device
;
; 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'.
;
; 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
ff=2
; edge length of 'decimated' image
l=256L/ff
; 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
dk1='empty'
; 'decimate' the image
bi=REBIN(i,l,l)
i='empty'
; plot the 'decimated' image for comparison
loadct,13
IF (NOT KEYWORD_SET(device)) OR (STRUPCASE(device) EQ 'X') THEN BEGIN
    WINDOW,2
    TVSCL,REBIN(ROTATE(bi,2),512,512)
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'
; 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
; 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.72*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 psopen,'/tmp/tmp.ps' $
    ELSE WINDOW,0
; create map title
mtitle=image
PRINT,SYSTIME(1)-t0,' seconds elapsed...'
; 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,imglat,tr,b
PRINT,SYSTIME(1)-t0,' seconds elapsed...'
PRINT,'TRIGRID...'
gridi=FIX(TRIGRID(imglon,imglat,bi,tr,[0.05,0.025] $
		 ,[wedge,sedge,eedge,nedge]))
; 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
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)
    date=INTARR(3)
    FOR j=0,2 DO date(j)=head.misc.dt.(j)
    indate=((date(2) MOD 100)*100+date(1))*100L+date(0)
    time=INTARR(3)
    FOR j=0,2 DO time(j)=head.misc.tm.(j)
    ut=time(2)+60.*(time(1)+60.*time(0))
    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
    plotorb,sat,indate,stime,etime,0.25,'x' $
	   ,/noplot,ptime,plat,plon
    USERSYM,[1,0,-1,0],[0,1,0,-1],/FILL
    stop
    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/tmp.ps' 
ENDIF
END