PRO plotstar,name,threshmag,maxza,sobel=sobel,grid=grid
;
; Procedure to read in a named image, get the time and date from
; the header, load the bright star catalog, calculate their positions
; at the given date and time from Rabbit Lake, SK, select stars
; brighter than the specified threshold magnitude and within the
; specified maximum zenith angle (if desired), transform those
; positions to image locations, plot the image, and overplot the
; stars.  If the sobel keyword is set, work with SOBEL(img).  If
; the grid keyword is set, overlay an azimuth/zenith angle grid 
; on the image.
;
; First get the requested image
;
rdkimg,name,hb,img
header=gethd(hb)
;
; Get the required dark frame and subtract it
;
cam=header.misc.cam
exp=header.misc.exp_scale*header.exp.exposure/1000
rdmeandk,cam,exp,dk
img=(img-dk)>0
IF KEYWORD_SET(sobel) THEN img=SOBEL(img)
;
; Then get the date and time of the image
;
getdatetime,header,date,time
;
; Then get the topocentric star coordinates for Rabbit Lake
;
RLlat=58.233
RLlon=256.330
stephem,date(0)+1900,date(1),date(2),time(0),time(1),time(2) $
       ,RLlat,RLlon,n,stid,az,za,mag
;
; Reject stars fainter than the specified threshold magnitude
;
IF NOT KEYWORD_SET(maxza) THEN maxza=71.
keep=WHERE((mag LT threshmag) AND (za LE maxza))
n=N_ELEMENTS(keep)
stid=stid(keep)
az=az(keep)
za=za(keep)
mag=mag(keep)
; big=WHERE(STRMID(stid,11,3) EQ 'UMA')
; n=N_ELEMENTS(big)
; stid=stid(big)
; az=az(big)
; za=za(big)
; mag=mag(big)
;
; Define the parameters of the transformation
; (present values derived from camera 0, filter 1 only)
;
cam_azimuth=-13.530
warp=0.615301
stretch=207.322
center_col=130.847
center_row=126.224
;
; Transform star coordinates to image locations
;
rdum=REPLICATE(1.,n)
polrec3d,rdum,za,-az,x,y,z,/degrees
rot_3d,3,x,y,z,cam_azimuth,xp,yp,zp,/degrees
rho=stretch*SIN(warp*za*!DTOR)
phi=+ATAN(yp,xp)
col=center_col+rho*SIN(phi)
row=center_row-rho*COS(phi)
stop
;
; Define the approximate maximum useful DN to display
;
ih=HISTOGRAM(img,MIN=0,BIN=1)
cih=ih & FOR i=1,N_ELEMENTS(cih)-1 DO cih(i)=cih(i)+cih(i-1)
ceiling=MAX(WHERE(cih LT 0.999*N_ELEMENTS(img)))
; IF KEYWORD_SET(sobel) THEN ceiling=20*ceiling ELSE ceiling=5*ceiling
;
; Rebin the image and adjust the col and row arrays to suit
;
zoom=2
img=REBIN(img,zoom*256,zoom*256,/SAMPLE)
col=zoom*col
row=zoom*row
;
; Plot the image with the first row at the top, and change the row
; array to suit
;
!ORDER=1
row=REPLICATE(zoom*256,n)-row
;
; Now plot the image and overplot the stars
;
TVSCL,img<ceiling
IF KEYWORD_SET(sobel) THEN psym=7 ELSE psym=4
PLOT,col,row,/DEVICE,/NOERASE,POSITION=[0,0,zoom*256-1,zoom*256-1] $
    ,XRANGE=[0,zoom*256-1],XSTYLE=5,YRANGE=[0,zoom*256-1],YSTYLE=5 $
    ,PSYM=psym
;
; Generate the grid and display it, if requested.
;
IF KEYWORD_SET(grid) THEN BEGIN
    ncirc=7
    circint=10
    nrays=12
    devradii=FLTARR(ncirc)
    FOR i=1,ncirc DO BEGIN
	devradii(i-1)=zoom*stretch*SIN(warp*i*circint*!DTOR)
	tvcirc,zoom*center_col,zoom*center_row,devradii(i-1)
    ENDFOR
    FOR i=1,nrays DO radii,devradii(0),devradii(ncirc-1) $
		          ,90+cam_azimuth-30*i,zoom*center_col $
			  ,zoom*center_row,/device
ENDIF
;
; We're done!
;
RETURN
END