! Procedure to make a PSF from image in buffer ibuf, with name fname

parameter ibuf string=frame inter

cubic=0

! Can do all tasks or just partial tasks as specified by inter:
!  inter=0  do everything
!  inter=1  do non-interactive preprocessing, i.e. zap frames, find stars,
!           aperture photometry, group stars
!  inter=2  do interactive marking of stars only, previous step must have been
!           done
!  inter=3  do post-interactive processing, i.e. make the actual PSFs, previous
!           step must have been done. Still displays images for inspection
!  inter=4  ?show PSF subtractions?


! Photometry parameters
opt fw=fw ip=1 ps=psfrad ng=1 fi=fi
opt ma=maxbad wa=0

again:

! Find stars
if inter==0|inter==2|inter==4
  daosky $ibuf 
  tv $ibuf l=max[10.,10.*skysig] z=sky-3*skysig
  !mn $ibuf pix=9 silent
  !tv $ibuf l=2.*mean
end_if
fits $ibuf float=ronoise rn
fits $ibuf float=gain gain

! Setup file names
$mkdir psf
daofiles coo=psf/{frame}.coo mag=psf/{frame}.ap 
daofiles psf=psf/{frame} pro=psf/{frame} 
daofiles grp=psf/{frame}

! Do the finding, ignoring masked pixels
if inter<2
  ! Use mask for instrument plus any one defined for this specific frame
  unmask
  err continue
  get mask=/usr/tmp/blank
  masktoim 101 nr=nr[ibuf] sr=sr[ibuf] nc=nc[ibuf] sc=sc[ibuf]
  get mask=msk/{frame}.msk
  masktoim 102 nr=nr[ibuf] sr=sr[ibuf] nc=nc[ibuf] sc=sc[ibuf]
  mul 101 102
  unmask
  clip 101 min=0.5 maskonly
  clip $ibuf maskonly max=maxbad
  clip $ibuf maskonly min=lowbad
  clip $ibuf vmask=lowbad-10

  sky $ibuf mask
  thresh=3*skysig
  copy 101 $ibuf
  do iter=1,5
    zap 101 sig=5 size=3
  end_do
  if inter==0
    tv $ibuf old noerase
  end_if
  if inter==0|inter==1
    find 102 101 thresh=thresh lowbad=lowbad relerr=relerr
    find 102 101 thresh=thresh*relerr lowbad=lowbad
    ! Find bright stars
    daofiles coo=psf/{frame}bright.coo
    thresh=10*skysig
    find 102 101 thresh=thresh*relerr lowbad=lowbad nwant=100,200
    daofiles coo=psf/{frame}.coo
  end_if
  unmask
end_if

err goto nozap
open zapfile msk/{frame}.zap
close zapfile
lzap $ibuf list=msk/{frame}.zap surgical sig=0
nozap:

! Display current star list
if inter==0|inter==2
  get coo=psf/{frame}.coo sr=sr[ibuf]-1 sc=sc[ibuf]-1
!  mark box=1 exit
end_if

! Do first pass with grouping by separation only, then second pass using
!  the predetermined PSF for better grouping by critical ratio. Actually,
!  grouping by separation is more conservative, so in most cases, just use
!  single pass!
npass=1
do ipass=1,npass

 buf=ibuf
 ! Various subtracted iterations will be placed into buffers ibuf+1...
 sbuf=ibuf+1

 ! Aperture photometry with relatively smaller aperture to be used for
 !  grouping (in pass 2)
 if inter<2
   photom $ibuf rad=r1 skyrad=skyrad1,skyrad2 gain=gain rn=rn
   daofiles coo=psf/{frame}bright mag=psf/{frame}bright grp=psf/{frame}bright
   photom $ibuf rad=r1 skyrad=skyrad1,skyrad2 gain=gain rn=rn

   ! Get candidate PSF stars by grouping of master list
   if ipass==1
     opt ig=1
     critval=critsep
     fudge=1.
    critval=2*ps
     regroup:
     group crit=critval*fudge
     daofiles file=psf/{frame}bright.grp file2=psf/{frame}psf.grp
     select size=1,1
!     open psf psf/{frame}psf.grp
!     stat npsf=count[psf]
!     close psf
!     if (npsf<maxpsf)
!       fudge=fudge/1.1
!       goto regroup
!     end_if
   else
     opt ig=0
     critval=100.
     daofiles psf=psf/{frame}.psf
     daofiles grp=psf/{frame}psf.grp
     group crit=critval fi=fi ps=ps psf
   end_if
   daofiles file=psf/{frame}psf.grp file2=psf/tmp.grp
   sort index=4 norenum
   $mv psf/tmp.grp psf/{frame}psf.grp 
   daofiles coo=psf/{frame} mag=psf/{frame} grp=psf/{frame}

 ! Before making final PSF, redo aperture photometry with big aperture 
 !  radius to get better PSF normaliziations between frames
   if (ipass==npass)
     photom $ibuf rad=r0 skyrad=skyrad1,skyrad2 gain=gain rn=rn
   end_if
 end_if

 if inter==1
   goto done
 end_if

 if inter<3
   ask 'Enter maximum number of PSFs (# to display individually): ' maxpsf
   !Show the maxpsf brightest semi-isolated stars, then let the user manually
   !  select (if inter==2)
   get coo=psf/{frame}psf.grp sr=sr[ibuf]-1 sc=sc[ibuf]-1
   open psf psf/{frame}psf.grp
   stat ngrp=count[psf]
   ngrp=ngrp-3
   read psf
   read psf
   read psf
   maxpsf=min[ngrp,maxpsf]
   npsf=0
   do i=1,maxpsf
     read psf
     m=@psf.4
     if (i==1)
       bright=m
     else
       m=@psf.4
       if (m-bright>3) 
         goto next
       end_if
     end_if
     if (m-bright<1)
       color=2
     else_if (m-bright<2)
       color=1
     else
       color=3
     end_if
     box 1 cr=@psf.3 cc=@psf.2 nr=ps*2 nc=ps*2
     tvbox box=1 color=color
     string rvar 'r%i2.2' i
     {rvar}=@psf.3
     string cvar 'c%i2.2' i
     {cvar}=@psf.2
     npsf=npsf+1
     next: 
   end_do
   close psf

   ! User selection of stars
   printf 'To select star, mark PSF stars with C key in TV window. '
   printf 'Enter E to go to next star. '
   string new 'new'
   j=1
   do i=1,npsf
     string rvar 'r%i2.2' i
     string cvar 'c%i2.2' i
     box 1 cr={rvar} cc={cvar} n=ps*4
     create 101 box=1
     srow=sr[101] scol=sc[101] erow=sr[101]+nr[101]-1 ecol=sc[101]+nc[101]-1
     srow0=sr[ibuf] scol0=sc[ibuf]
     erow0=sr[ibuf]+nr[ibuf]-1 ecol0=sc[ibuf]+nc[ibuf]-1
     if (srow<srow0)
       srow=srow0
     else_if (erow>erow0)
       erow=erow0
     end_if
     if (scol<scol0)
       scol=scol0
     else_if (ecol>ecol0)
       ecol=ecol0
     end_if
     box 1 sr=srow sc=scol nr=erow-srow+1 nc=ecol-scol+1
     abx $ibuf 1 silent high=high
   
     ! move image subsection to another buffer to avoid showing stars
     ! that arent on the image 
     fits $ibuf int=cnpix1 1
     fits $ibuf int=cnpix2 1
     copy 199 $ibuf  box=1
     tv 199 z=sky-3*skysig l=max[20.*skysig,(high-sky)/100]
     get coo=psf/{frame}.coo sr=sr[ibuf]-1 sc=sc[ibuf]-1
     mark box=1 exit 
     err continue
     get coo=psf/psf sr=0 sc=0
     mark circ=ps nobox surgical sig=0 {new} id=j |
          zapfile=msk/{frame}.zap append
     err goto nextstar
     save coo=psf/psf sr=0 sc=0
     string new ' '
     j=j+1
     nextstar:
   end_do
   tv $ibuf l=max[10.,10.*skysig] z=sky-3*skysig
   printf 'Add more stars, E to continue...'
   mark circ=ps id=j surgical sig=0 zapfile=msk/{frame}.zap append
!   do i=1,99
!     string rvar 'r%i2.2' i
!     string cvar 'c%i2.2' i
!     if i<10
!       string idvar 's%i1' i
!     else
!       string idvar 's%i2' i
!     end_if
!     {rvar}=0
!     {cvar}=0
!     {idvar}=0
!   end_do
!   npsf=0
!   mark circ=ps nobox new surgical sig=0
   ! Save list of selected PSF stars
   save coo=psf/psf sr=0 sc=0

   !Get original list back
   open injunk psf/psf.coo
 else
   err goto useauto
   open injunk psf/{frame}psf.coo
 end_if

 stat npsf=count[injunk]
 npsf=npsf-3
 read injunk
 read injunk
 read injunk
 do i=1,npsf
     read injunk
     string rvar 'r%i2.2' i
     string cvar 'c%i2.2' i
     {rvar}=@injunk.3 
     {cvar}=@injunk.2
 end_do
 close injunk

 useauto:
 get coo=psf/{frame} sr=0 sc=0

 open injunk psf/{frame}.coo
 string line '{injunk}'
 printf '{line}' >psf/{frame}psf.coo 
 string line '{injunk}'
 printf '{line}' >>psf/{frame}psf.coo 
 string line '{injunk}'
 printf '{line}' >>psf/{frame}psf.coo 
 close injunk
 do i=1,npsf
   string rvar 'r%i2.2' i
   string cvar 'c%i2.2' i
   if i<10
     string idvar 's%i1' i
   else
     string idvar 's%i2' i
   end_if
   {idvar}=id[{rvar},{cvar}]
   if i==1
     string psfs '{idvar}'
   else
     string psfs |
'{psfs},{idvar}'
   end_if
   printf '%i6%2f9.2    0.000    0.000    0.000       0.     0.00    1.000' |
      {idvar} {cvar} {rvar} >>psf/{frame}psf.coo
 end_do

 if inter==2
   goto done
 end_if

 if inter==4
   goto showpsf
 end_if

 ! Now make the PSF
 ! We return here to successively iterate the PSF with neighbors subtracted
 old:

 ! Make individual PSF files for each star for field-dependent PSF
 if (sbuf-ibuf)==4
  nnpsf=0
  do i=1,npsf
   if i<10
     string idvar 's%i1' i
   else
     string idvar 's%i2' i
   end_if
   string cvar 'c%i2.2' i
   string rvar 'r%i2.2' i
   err goto nextpsf
   daofiles psf=psf/{frame}.psf
   psf $buf stars={idvar}
   nnpsf=nnpsf+1
   !  Turn it into a library and normalize it
   daolib n=4
   rd 101 psf/{frame}.lib
   tot={101:psftot}
   if tot<1
     nnpsf=nnpsf-1
     goto nextpsf
   end_if
   mul 101 c=1/tot
   fits 101 float=psftot 1.
   fits 101 float=row {rvar}-nr[buf]/2
   fits 101 float=col {cvar}-nc[buf]/2
  row0=nr[101]/2+4 col0=nc[101]/2+4
  copy 102 101
  axes 102
  shift 102 dr=row0-axr dc=col0-axc
  sub 101 101
  add 101 102
   string pframe '{frame}_%i2.2' nnpsf
   wd 101 psf/{pframe}.lib full
   nextpsf:
  end_do
  npsf=nnpsf
 end_if

 ! Make single PSF (after individuals so we have the correct PSG file)
 daofiles psf=psf/{frame}.psf
 psf $buf stars={psfs}

 !  Turn it into a library and normalize it
 daolib n=4
 rd 101 psf/{frame}.lib
 tot={101:psftot}
 mul 101 c=1/tot
 fits 101 float=psftot 1.
 wd 101 psf/{frame}.lib full

 !  Show results
 get coo=psf/{frame} sr=sr[ibuf]-1 sc=sc[ibuf]-1
 printf ' PSF stars are marked with biggest boxes, neighbors with '
 printf '    somewhat smaller boxes '
 if inter==0
   mark circ=ps star={psfs} color=2
   mark box=4 exit
 end_if

 !subtract neighbors, or all stars after final pass (to inspect for PSF
 ! variation) 
 opt ip=0 ig=4
 daofiles grp=psf/{frame}.psg
 multistar $ibuf exp=1 fi=fi skyrad=skyrad1,skyrad2
 opt ig=0
 copy $sbuf $ibuf
 daofiles pro=psf/{frame}.nst
 if (sbuf-ibuf)<=4
   sub* $sbuf exclude={psfs} subonce
   ! Redo aperture photometry for better normalizations
   daofiles coo=psf/{frame} mag=psf/{frame}
   photom $sbuf rad=r0 skyrad=skyrad1,skyrad2 gain=gain rn=rn
   ! Display current iteration
   if inter==0
     tv $sbuf old noerase
     mark circ=ps star={psfs} color=2
     mark box=4 exit
   end_if
 end_if
 opt ip=1

 !  Iterate again?
 !ask 'Enter a 0 to end, else redo subtraction: ' cont
 if ((sbuf-ibuf)<4)
   buf=sbuf
   sbuf=sbuf+1
   goto old
 end_if

end_do

! Make field dependent PSF basis functions
string vpsfs ' '
if cubic==1
  npar=10
  string cubic 'cubic' 
else
  npar=6
  string cubic ' '
end_if
do i=1,npsf
   string pframe '{frame}_%i2.2' i
   rd $sbuf+i psf/{pframe}.lib full
   if i<=npar
     copy $100+i $sbuf+i
     j=100+i
   else
     j=sbuf+i
   end_if
   string vpsfs '{vpsfs} %i3' j
   $'rm' psf/{pframe}.lib
end_do
psfinterp {vpsfs} header {cubic}
unfit 101 card=row
unfit 101 card=col
fits 101 int=naxis 3
fits 101 int=naxis3 npar
fits 101 float=x0 nc[buf]/2
fits 101 float=y0 nr[buf]/2
wd 101 psf/{frame}var.lib noauto notail full
do i=2,npar
   wd $100+i psf/junk.lib noauto nohead notail full
   $cat psf/junk.lib >>psf/{frame}var.lib
end_do

! Determine the seeing
get coo=psf/{frame}psf.coo sr=0 sc=0
starplot 1 phot load scale=scale noplot silent load gauss
printf '%f10.2' fwavg >psf/{frame}.seeing

! Now fit and subtract using field-dependent PSF to compare
opt ip=-1 ig=4
$cp psf/{frame}.psg psf/{frame}var.psg
daofiles psf=psf/{frame}var.lib grp=psf/{frame}var.psg
multistar $ibuf exp=1 fi=fi skyrad=skyrad1,skyrad2

! Come here to display residuals
showpsf:
get coo=psf/{frame}.coo sr=0 sc=0
sbuf=sbuf+1
copy $sbuf $ibuf
opt ip=0
daofiles pro=psf/{frame}.nst psf=psf/{frame}.lib
sub* $sbuf subonce
if inter==4
tv $sbuf old noerase
end_if

sbuf=sbuf+1
copy $sbuf $ibuf
daofiles pro=psf/{frame}var.nst psf=psf/{frame}var.lib
opt ip=-1
sub* $sbuf subonce
if inter==4
tv $sbuf old noerase
mark circ=ps star={psfs} color=2
mark box=4 exit
end_if

opt ip=1

daofiles none

!$'rm' psf/{frame}.ap psf/{frame}.dpos psf/{frame}.grp psf/{frame}.nst 
!$'rm' psf/{frame}.psf
!$'rm' psf/{frame}psf.grp 
!$gzip -f psf/{frame}.coo

done:
end