;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 June 8, 1994 by T A Oliynyk: can now be used on Unix or DOS machines.
; Modified June 9, 1994 by T A Oliynyk: an additional parameter, wkimg, was 
;	added. If this procedure is passed an image file, it will use this
;	image to plot the stars on.
;         Also an amplification factor, amp, was added to scale the limits
;       of the image plot.


PRO mstplot,starcols,starrows,amp,file2,wkimg

@isitdos

;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. Calgary'
PRINT,'2. Rabbit Lake'
PRINT,'3. Eureka'
READ,site
IF site EQ 2 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'

;if an image file is passed to the procedure, then this image will be used. if no
;image is passed, then the image named by the variable file2 will be used.
IF N_PARAMS() EQ 5 THEN BEGIN
	img=wkimg
ENDIF ELSE BEGIN
	rdmeandk,cam,60,dk60
	img=(img1-dk60)>0
ENDELSE

;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)
PRINT,bimgmin,pc99
;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 parameters, read in data from meanparms.?? files.  
rdmnparm,yr,mo,da,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

;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

; Make a mock-up star map for the location and time of the actual image.

!ORDER = 0
stmap=BYTARR(256,256)
cell=BYTARR(5,5,5)
cell(0,0,0)=[[1,0,1,0,1],[0,1,1,1,0],[1,1,1,1,1],[0,1,1,1,0],[1,0,1,0,1]]
cell(0,0,1)=[[0,0,1,0,0],[0,0,1,0,0],[1,1,1,1,1],[0,0,1,0,0],[0,0,1,0,0]]
cell(0,0,2)=[[0,0,0,0,0],[0,1,1,1,0],[0,1,1,1,0],[0,1,1,1,0],[0,0,0,0,0]]
cell(0,0,3)=[[0,0,0,0,0],[0,1,0,1,0],[0,0,1,0,0],[0,1,0,1,0],[0,0,0,0,0]]
cell(2,2,4)=1
;Repeat for all layers of cell
FOR i=-1,3 DO BEGIN
    magi=WHERE(FLOOR(mag) EQ i,nmagi)
    IF nmagi GT 0 THEN BEGIN
        magn=BYTARR(256,256)
	magn(cl(magi)/2.,rw(magi)/2.)=1B
	magn=DILATE(magn,cell(*,*,i+1))
	stmap=stmap+magn
	ENDIF
    ENDFOR
TVSCL,stmap<1,2
!ORDER = 1

; Superimpose calculated star positions on the image.
bimg=bimg+REBIN(ROTATE(stmap,7),512,512,/sample)*amp*percentile(bimg,0.9)/2.
loadct,3,/silent

;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'

window,2

;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,dd,dd,cam,filt,STRMID(file2,STRLEN(file2)-12,9),hr,min,sec $
     ,FORMAT='(A,A1,"starpos",A1,"stpos.c",I1,".f",I1,A,3I2.2)')
OPENW,unit,dfile,/GET_LUN
j=-1
k=0
PRINT,'Y for a match, <RET> for next candidate, P if no matches'
FOR i=0,N_ELEMENTS(starcols)-1 DO BEGIN

  IF (i NE 0) AND (k LT 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=amp*pc99,TOP=!D.N_COLORS-1) $
	,oldcol-10,511-oldrow-10
      ENDIF ELSE tvscl,bimg<amp*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*amp*pc99) THEN pcolor=0 ELSE pcolor=!D.N_COLORS-1

  ;plot position of bright point to be identified, as a diamond shape
  PRINT,STRTRIM(ROUND(2.*starcols(i)),2),',',STRTRIM(ROUND(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