;procedure to plot given row and column coordinates onto an image.
;called by starid.pro to plot the positions of stars onto 256X256 images.

;inputs: starcols,starrows - row and column positions of points to be plotted
;	 file2 - file name of image to be overplotted onto
; Modified 94/06/03 by D P Steele: now reads mean camera transformation
; parameters from standard location using rdmnparm.pro
; Modified May 8, 1994 by T A Oliynyk: can now be used on unix of dos machines.

PRO stplot,starcols,starrows,file2

@isitdos
!ORDER = 1
file3=''

;read in image
rdkimg,file2,head,img1
kh=gethd(head)
cam=kh.misc.cam

;find time of exposure of image
getdatetime,kh,date,time
yr = date(0)+1900 & mo=date(1) & da=date(2) & hr=time(0) & min=time(1) & sec=time(2) 

;rabbit lake latitude and longitude
PRINT,'Choose Polar Camera site for this image'
PRINT,'1. Rabbit Lake'
PRINT,'2. Eureka'
READ,site
IF site EQ 1 THEN BEGIN
    lat=58.233
    lon=256.330 
    ENDIF ELSE BEGIN
    lat=80.053
    lon=273.584
    ENDELSE
cam=kh.misc.(0) & filt=kh.misc.(1) & exp=kh.exp.(1)

;find positions of all visible bright stars at time of exposure
stephem,yr,mo,da,hr,min,sec,lat,lon,n,desig,az,za,mag
PRINT,STRTRIM(n,2),' stars read from catalog'

;read in mean dark frame to subtract - not really necessary
rdmeandk,cam,60,dk60
img=(img1-dk60)>0

;rebin image to 512X512 for easier viewing when points are plotted 
bimg=rebin(img,512,512,/sample)
pc99=percentile(bimg,0.99)
bimgmin=MIN(bimg)
wsat=WHERE(bimg GT 1.5*pc99)
bimg(wsat)=bimgmin

;plot only stars with magnitude less than or equal to 4.00
bright=where(mag le 4.00)
numst = N_ELEMENTS(bright)
mag = mag(bright)
desig=desig(bright)
az = az(bright)
za =za(bright)
PRINT,STRTRIM(numst,2),' stars available to match'

;to get correct stretch and warp paramaters, read in data from .meanparms files.  
rdmnparm,cam,c0mean,r0mean,rmean,kmean,camean,site=site
camaz = -camean*!RADEG
IF site EQ 2 THEN camaz=0.
warp = kmean
stretch = rmean
ccol = c0mean
crow=r0mean
PRINT,ccol,crow,stretch,warp,camaz*!RADEG

;find column and row positions of the stars read in using stephem.pro 
;from their azimuth and zenith angles 
polrec3d,REPLICATE(1.,numst),za,-az,x,y,z,/degrees
rot_3d,3,x,y,z,camaz,xp,yp,zp,/degrees
rho=stretch*SIN(warp*za*!DTOR)
phi=+ATAN(yp,xp)
cl=2.*(ccol+rho*sin(phi))
rw=2.*(crow-rho*COS(phi))
rw = REPLICATE(512.,numst)-rw
;window,0

;It is not useful to plot stars outside of the fov.  Only take
;stars within a certain radius of the center
notedge = where(SQRT((starcols-ccol)^2. + (starrows-crow)^2.) LT 135)
ct=starcols & rt=starrows
starcols = starcols(notedge)
starrows = starrows(notedge)
PRINT,STRTRIM(N_ELEMENTS(notedge),2),' candidate stars found'
dd=shift(dist(512,512),256,256)
x=-x


;dfile is the file to be written to to save all relevant data.  STRMID(file2,30,1)
;is the image type identification letter (`a' or `g' here)
dfile=STRING(root,cam,filt,STRMID(file2,STRLEN(file2)-12,9),hr,min,sec $
     ,FORMAT='(A,"/starpos/stpos.c",I1,".f",I1,A,3I2.2)')
OPENW,unit,dfile,/GET_LUN
j=-1
k=0
FOR i=0,N_ELEMENTS(starcols)-1 DO BEGIN

  IF (i NE 0) AND (k LE 1) THEN BEGIN
      oldcol=2.*starcols(i-1)
      oldrow=2.*starrows(i-1)
      PRINT,oldcol,oldrow,' patch'
      piece=extrac(bimg,oldcol-10,oldrow-10,21,21)
      TV,BYTSCL(piece,MIN=bimgmin,MAX=pc99,TOP=!D.N_COLORS-1) $
	,oldcol-10,511-oldrow-10
      ENDIF ELSE tvscl,bimg<pc99

  ; choose colour for symbols plotted, depending on background
  meanbg=mean(extrac(bimg,2.*starcols(i),2.*starrows(i),21,21))
  IF meanbg GT (0.6*bimgmin+0.4*pc99) THEN pcolor=0 ELSE pcolor=!D.N_COLORS-1

  ;plot position of bright point to be identified, as a diamond shape
  PRINT,STRTRIM(2.*starcols(i),2),',',STRTRIM(2.*starrows(i:i),2)
  PLOT,2.*starcols(i:i),2.*(-starrows(i:i)+256) $
      ,CLIP=[0,0,511,511],/DEVICE,XMARGIN=[0,0],YMARGIN=[0,0],/NOERASE $
      ,SYMSIZE=2,POSITION=[0,0,511,511],PSYM=4,COLOR=pcolor $
      ,XRANGE=[0,511],YRANGE=[0,511],XSTYLE=5,YSTYLE=5

  j=-1
  k=0
  ans=''
  ; plot positions of known stars until the user has found one in the 
  ; same position as the bright point
  ps=1
  REPEAT BEGIN
      j=j+1
      IF SQRT((2.*starcols(i)-cl(j))^2 $
	      +(2.*(-starrows(i)+256)-rw(j))^2) LT 40 THEN BEGIN
	  k=k+1
          PRINT,desig(j),mag(j)
          PLOT,cl(j:j),rw(j:j),CLIP=[0,0,511,511],/DEVICE $
              ,XMARGIN=[0,0],YMARGIN=[0,0],/NOERASE,SYMSIZE=2 $
              ,POSITION=[0,0,511,511],PSYM=ps,COLOR=pcolor $
              ,XRANGE=[0,511],YRANGE=[0,511],XSTYLE=5,YSTYLE=5
          ans = '' & READ,ans
          ENDIF
      IF ((ans EQ '') AND (j EQ numst-1)) THEN BEGIN	; No matches; on to the next candidate
	  IF (k EQ 0) THEN BEGIN
	      PRINT,'No stars match; go on to next star'
	      ans='p'
	      ENDIF ELSE BEGIN
	      j=-1	;back to start of list (-1 because j=j+1 at start of loop)
	      ps=ps+1
	      IF ps EQ 3 THEN ps=4
	      ENDELSE
          ENDIF

      ENDREP UNTIL (ans EQ 'y') OR (ans EQ 'p')
      ; stop plotting positions when the user agrees to the position 
      ; or finds no match and calls the spot a planet
  
  IF ans EQ 'y' THEN BEGIN
     print,desig(j),mag(j)
     printf,unit,strmid(desig(j),0,4),az(j)*!DTOR,za(j)*!DTOR,starcols(i),starrows(i),$
     FORMAT='(I5,2f12.6,2I4)'  ; write to stpos file the position of the star

     print,desig(j:j),az(j),za(j)
  ENDIF
ENDFOR
close,unit & FREE_LUN,unit

end