;+
; NAME:
;	Medit
; PURPOSE:
;	Edit Mathematica program fragment created by StPlotFix.PRO to
;       correct syntax errors. Optionally, tack on the rest of the
;       Mathematica code needed to create a full program.
; CATEGORY:
;	
; CALLING SEQUENCE:
;	Medit, file, test=test, append=append
; INPUTS:
;	file:	Fully-qualified path name of the Mathematica program
;		fragment file created by StPlotFix.PRO.
; OPTIONAL INPUTS:
;	None.
; KEYWORD PARAMETERS:
;	/test:	Write the new file to file 'Medit.m' to allow comparison
;		with the original file.
;       /append: Append the rest of the Mathematica code necessary to
;               have a full Mathematica program. Ths option is used when
;               the procedure is invoked by StPlotFix.
; OUTPUTS:
;	None.
; OPTIONAL OUTPUTS:
;	None.
; COMMON BLOCKS:
;	None.
; SIDE EFFECTS:
;	The specified file is opened and read. After being edited
;	internally a new file of the same name is written over the
;	original file, unless the /test keyword is set.
; RESTRICTIONS:
;	The named file must exist and have the expected structure.
; PROCEDURE:
;	The file is read as a byte string. The erroneous byte patterns
;	are located using StrPos, and the correct byte patterns are
;	inserted. Optionally, the remaining Mathematica code from
;       'starfit.m' is read and appended to the byte array before it is
;       written out.
; EXAMPLE:
;	Medit,'c:\data\starpos\g0601060.m',/test
; SEE ALSO:
;	StarID.PRO, StPlotFix.PRO
; MODIFICATION HISTORY:
; 	Written by:	DPSteele, Jan. 30/97.
;-

PRO Medit, file, test=test, append=append

@isitdos

IF Keyword_Set(test) THEN Stop

; Read in the file.

OpenR, lun, file, /Get_LUN
fs=FStat(lun)
nb=Fs.SIZE
bytes=BytArr(nb)
ReadU,lun,bytes
Free_LUN,lun

; Define the erroneous and correct byte patterns.

patt1=[123B, 13B, 10B, 123B]
len1=4
patt2=[125B, 44B, 13B, 10B, 125B]
len2=5

rpatt1=[123B, 32B, 123B]
rpatt2=[125B, 32B, 125B]

; Locate and replace the offending bytes.

pos1=StrPos(String(bytes),String(patt1))
IF pos1 LT 0 THEN BEGIN
    Message,'patt1 not found; aborting',/INFORMATIONAL
    Return
    ENDIF ELSE bytes=[bytes(0:pos1-1),rpatt1,bytes(pos1+len1:*)]

pos2=StrPos(String(bytes),String(patt2))
IF pos2 LT 0 THEN BEGIN
    Message,'patt2 not found; aborting',/INFORMATIONAL
    Return
    ENDIF ELSE bytes=[bytes(0:pos2-1),rpatt2,bytes(pos2+len2:*)]

; Load and append the rest of the Mathematica code from 'starfit.m'

IF Keyword_Set(append) THEN BEGIN
    mfile=starroot+dd+'starfit.m'
    mfile=rFindFile(mfile,count=nmf)
    IF nmf EQ 0 THEN BEGIN
        Message,mfile(0)+' not found',/INFORMATIONAL
        Return
        ENDIF ELSE mfile=mfile(0)
    OpenR,mlun,mfile,/Get_LUN
    mfs=FStat(mlun)
    mnb=mfs.SIZE
    mbytes=BytArr(mnb)
    ReadU,mlun,mbytes
    Free_LUN,mlun

    bytes=[bytes,mbytes]
    ENDIF

; Write out the corrected file.

IF Keyword_Set(test) THEN file=starroot+dd+'Medit.m'
OpenW,lun,file,/Get_LUN
WriteU,lun,bytes
Free_LUN,lun

Return
END

; ----------------------------------------------------------------------

PRO stplotfix,starcols,starrows,amp,file2,wkimg

;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
;        NOTE: These values are actual positions for the given image and thus
;              depend on the image's actual dimensions: 128x128 or 256x256.
;         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.



@isitdos

;read in image
rdkihd,file2,head
kh=gethd(head)
cam=kh.misc.cam
imcols=kh.exp.cols/kh.exp.pbin
imrows=kh.exp.rows/kh.exp.sbin

; Normalize positions of observed stars to a standard 256x256 image
normfac=256/imcols
starcols=normfac*starcols
starrows=normfac*starrows

;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
    dksize=SIZE(dk60)
    imsize=SIZE(img)
    IF dksize(1) NE imsize(1) THEN dk60=binf(dk60,256/imsize(1))
    img=(img1-dk60)>0
    ENDELSE

;rebin image to 512X512 for easier viewing when points are plotted
bimg=rebin(img,512,512,/sample)
pc99=pcentile(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)
; The following two lines map the stars onto a 512x512 image.
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*pc99/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'

PRINT,'Type Y to accept, P to go on to next star in image, Q to quit,'
PRINT,'or any other key to try the next candidate star in the catalogue.'

window,2,XSIZE=512,YSIZE=512

;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)
IF dos $
    THEN dfile=STRING(starroot,dd,STRMID(file2,STRLEN(file2)-12,9) $
                     ,FORMAT='(A,A1,A,"stp")') $
    ELSE 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
PRINTF,unit,'stardat = {'
j=-1
k=0
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 & ans=STRLOWCASE(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') OR (ans EQ 'q')
      ; stop plotting positions when the user agrees to the position
      ; or finds no match and calls the spot a planet

  CASE ans OF
  'y': BEGIN
           print,desig(j),mag(j)
        ; write to stpos file the position of the star
           printf,unit,strmid(desig(j),0,4),az(j)*!DTOR,za(j)*!DTOR $
                      ,starcols(i),starrows(i) $
                      ,FORMAT='("{",I5,2(",",f12.6),2(",",I4),"},")'
           print,desig(j:j),az(j),za(j)
           END
  'q': GOTO, DONE
 ELSE:
      ENDCASE
ENDFOR

DONE:
PRINTF,unit,"}"
FREE_LUN,unit

; Now edit and augment the file so it is ready for input to Mathematica.

Medit,dfile,/append

Return
END