parameter ichip iflat

if ichip>0
  string outfield '{field}%i1' ichip
else
  string outfield {field}
end_if

! Set some options
opt wa=0 ma=highbad fi=fi ps=ps fw=fw ip=ip
ifake=0

! Error traps
if nfilt>4
  printf 'cant handle nfilt>4 right now'
  pause
end_if

! Read in the frames and sum up all frames to be used for finding, as
!  well as making sum frames for each color
! Load up {frames} and {exp} strings
! Create links for individual frame PSFs 
i=1
unmask
string frames ' '
string exp 'exp='
string gain 'gain='
string rn 'rn='
string col 'col='
do ifilt=1,nfilt
  string filt 'filt%i1' ifilt
  string var 'nframe%i1' ifilt
  do iframe=1,{var}
    string var2 'name%i1%i1' ifilt iframe
    if ichip>0
      string name '{{var2}}%i1' ichip
    else
      string name '{{var2}}'
    end_if
    string psf '{field}%i2.2' i
    if ip==-1
      $ln -fs {name}var.lib psf/{psf}.lib
    else
      $ln -fs {name}.lib psf/{psf}.lib
    end_if
    if i>13
      printf 'ERROR: multistar currently compile limited to 13 frames'
      pause
    end_if
    if i>20
      printf |
          'ERROR: phot.pro and register currently compile limited to 20 frames'
      pause
    end_if
    rd $i '{name}' {pictype}

!  Load up approximate frame offsets if appropriate coordinate files exist
    if i==1
      string ref {name}
    else
      err continue
      register 1 $i load ref=pos/{ref}.coo dat=pos/{name}.coo trans
    end_if 

!  If this is WFPC2 data, correct for pixel areas, and flatten if we
!    are using unflattened data
    if wfpc==1
     if iflat==0
      $echo {{filt}} | awk '\{print toupper($0)\}' > /tmp/tmp.lis
      open tmp /tmp/tmp.lis
      string ufilt {tmp}
      close tmp

      string flat 'iflat_{ufilt}_%i1' ichip
      rd 109 wfpc2cal/flats/{flat}
      mul $i 109
     end_if
     string file 'pixarea_%i1' ichip
     rd 109 wfpc2cal/distort/{file}
     mul $i 109

    end_if

    ! get proper nomalization for each frame, which is exposure time and
    !   possibly extinction correction (if extinction coef has been defined).
    !   load normalization factor into scaleII variables
    string scvar 'scale%i2.2' i
    string extvar 'ext%i2.2' i
    string k 'k%i1' ifilt
    err continue
    hjd $i silent
    err goto noair
    {extvar}=10^(0.4*{k}*{i:am})
    {scvar}={?i:exptime}*{extvar}
    air0=0
    goto haveair
    noair:
    {extvar}=1
    {scvar}={?i:exptime}
    air0=-1
    haveair:

    ! Create find buffer 101
    if (i==1)
      copy 101 1
      sub 101 101
      get mask=/usr/tmp/blank
      masktoim 102 nr=nr[i] sr=sr[i] nc=nc[i] sc=sc[i]
    end_if
   !get bad pixel mask 
    get mask=msk/{name}.msk
    masktoim 110 nr=nr[i] sr=sr[i] nc=nc[i] sc=sc[i]
    unmask
    clip maskonly 110 min=0.5

    string find 'find%i1' ifilt
    ! Add image into find frame if requested to do so
    if {find}==1
      add 101 $i
      ! Use mask for instrument plus any one defined for this specific frame
      err continue
      mul 102 110
    end_if

    ! Make individual frames for each filter to be used for aperture phot
    string var 'tot%i1' ifilt
    if iframe==1
      copy $105+ifilt $i
      {var}={scvar}/exp0
    else
      add $105+ifilt $i
      {var}={var}+{scvar}/exp0
    end_if

    ! Set bad pixels to lowbad-100, but after creating 105+ifilt frames since
    !  aperture photometry is done on these and will bomb with lots of bad
    !  pixels
    clip $i vmask=lowbad-100

    string frames '{frames} %i2.2' i

    exp1={scvar}
    if i>1
      string exp '{exp},' 
      string col '{col},'
      string gain '{gain},'
      string rn '{rn},'
    end_if
    if exp1/exp0<10
        string exp '{exp}%f5.3' exp1/exp0
    else_if exp1/exp0<100
        string exp '{exp}%f6.3' exp1/exp0
    else_if exp1/exp0<1000
        string exp '{exp}%f7.3' exp1/exp0
    else
        string exp '{exp}%f8.3' exp1/exp0
    end_if
    string col '{col}%i1' ifilt
    if (ichip>0) 
      call wfpcgain i ichip
      if gain>13
        rn=9
      else
        rn=6
      end_if
    end_if
    if gain>9.99
      string gain '{gain}%f5.2' gain
    else
      string gain '{gain}%f4.2' gain
    end_if
    if rn>9.99
      string rn '{rn}%i2' rn
    else
      string rn '{rn}%i1' rn
    end_if

    i=i+1
  end_do
end_do

! Normalize summed frames in each filter to get aperture photometry normalized
!  to default output exposure time for best first guess
do ifilt=1,nfilt
  string var 'tot%i1' ifilt
  div $105+ifilt c={var}
end_do

! Read in prexisting summed frames and find frame if they exist
if ichip>0
  string file '{field}%i1find' ichip
else
  string file '{field}find' ichip
end_if
err continue
rd 101 red/{file} maxtry=0
do ifilt=1,nfilt
  string filt 'filt%i1' ifilt
  if ichip>0
    string file '{field}{{filt}}%i1' ichip
  else
    string file '{field}{filt}' ichip
  end_if
  err continue
  rd $105+ifilt red/{file} maxtry=0
end_do

unmask
clip maskonly 102 min=0.5
! Mask out edges
box 1 nr=10 sr=1 nc=nc[101] sc=sc[101]
clip 101 maskonly max=-1e10 box=1
box 1 nr=10 sr=nr[101]-9 nc=nc[101] sc=sc[101]
clip 101 maskonly max=-1e10 box=1
box 1 nc=10 sc=1 nr=nr[101] sr=sr[101]
clip 101 maskonly max=-1e10 box=1
box 1 nc=10 sc=nc[101]-9 nr=nr[101] sr=sr[101]
clip 101 maskonly max=-1e10 box=1

! Get sky value
daosky 101 mask

! zap for CRs (on find frame only)
zap 101 size=3 sig=3
zap 101 size=3 sig=3

! Loop over multiple photometry passes
npass=2
do ipass=1,npass

 if ipass==1
  string pass 'a'
 else
  string pass 'b'
 end_if
 daofiles coo=./{outfield}{pass}
 if ipass==1
  ! First pass for geometric transformations only, consider just bright stars
  fudge=1
  refind:
  find 102 101 thresh=10*skysig*fudge lowbad=lowbad nfound=nfound
  if nfound>1000
    fudge=fudge*1.1
    type fudge
    goto refind
  else_if nfound<100
    fudge=fudge/1.1
    type fudge
    goto refind
  end_if
  string dpos 'mindpos=1 header'
 else
  find 102 101 thresh=3*skysig lowbad=lowbad relerr=relerr
  find 102 101 thresh=3*skysig*relerr lowbad=lowbad
  string dpos 'dposmax=1 header'
!  string dpos 'header'
 end_if

 ! Aperture photometry and grouping
 do ifilt=1,nfilt
  string filt 'filt%i1' ifilt
  string file {outfield}{{filt}}{pass}
  daofiles mag=./{file} pro=./{file} grp=./{file}
  photom $105+ifilt rad=r1 skyrad=skyrad1,skyrad2 gain=gain rn=rn {3sig}
  shortap 
  oldgroup crit=crit
 end_do
 
 ! Set up group file names for multistar
 do ifilt=1,nfilt
  string filt 'filt%i1' ifilt
  string file {outfield}{{filt}}{pass}
  if ifilt==1
    daofiles grp=./{file}.grp 
  else_if ifilt==2
    daofiles file=./{file}.grp 
  else_if ifilt==3
    daofiles file2=./{file}.grp 
  else_if ifilt==4
    daofiles file3=./{file}.grp 
  end_if
 end_do

 ! Do the profile photometry
 daofiles psf=psf/{field} pro=./{outfield}{pass}
 opt ig=4
 multistar {frames} |
   {col} |
   {exp} |
   {gain} |
   {rn} |
   {dpos} drift=5 exp0=exp0 gain0=gain0 air0=air0 skyrad=skyrad1,skyrad2 write

 ! Get the geometric transformation if first pass
 if (ipass==1)
   register {frames} load dpos=./{outfield}{pass}.dpos silent |
      niter=10 maxrms=0.1 cubic
 end_if

 ! Sort the output files by magnitude, and subtract stars from each frame
 i=1
 do ifilt=1,nfilt
  string filt 'filt%i1' ifilt
  string file {outfield}{{filt}}{pass}
  daofiles file=./{file}.nst file2=./{file}2.nst
  sort index=4 norenum
  $mv {file}2.nst {file}.nst
  daofiles pro=./{file}.nst 

  ! Subtract the stars
  string var 'nframe%i1' ifilt
  do iframe=1,{var}
    string psf '{field}%i2.2' i
    daofiles psf=psf/{psf}.lib
    copy $i+20 $i
    string scvar 'scale%i2.2' i
    exp1={scvar}
!    sub* $i+20 shift scale=exp1/exp0

    string pro '{outfield}{pass}%i2.2' i
    daofiles pro=./{pro}
    sub* $i+20 

    string var2 'name%i1%i1' ifilt iframe
    if ichip>0
      string name '{{var2}}%i1' ichip
    else
      string name '{{var2}}'
    end_if
    if ipass==npass
      wd $i+20 tmp/{name}s
    end_if
    i=i+1
  end_do
 end_do

! Done with this pass
end_do

! Now do aperture photometry on subtracted frame
i=1
do ifilt=1,nfilt
  string filt 'filt%i1' ifilt
  string file {outfield}{{filt}}{pass}
  daofiles pro=./{file}.nst 

  string var 'nframe%i1' ifilt
  do iframe=1,{var}
    string scvar 'scale%i2.2' i
    exp1={scvar}
    if (ichip>0) 
      call wfpcgain i ichip
      if gain>13
        rn=9
      else
        rn=6
      end_if
    end_if

    string pro '{outfield}{pass}%i2.2' i
    daofiles pro=./{pro}

    string psf '{field}%i2.2' i
    daofiles psf=psf/{psf}.lib
    string file '{outfield}_%i2.2ap' i
    daofiles file=./{file}.nst
!    subphot $i+20 $i rad=r1 scale=exp1/exp0*gain0/gain shift gain=gain rn=rn |
!       skyrad=skyrad1,skyrad2  {3sig}
    subphot $i+20 $i rad=r1 gain=gain rn=rn |
       skyrad=skyrad1,skyrad2 oscale=exp1/exp0*gain0/gain {3sig}

    string file '{outfield}_%i2.2ap0' i
    daofiles file=./{file}.nst
!    subphot $i+20 $i rad=r0 scale=exp1/exp0*gain0/gain shift gain=gain rn=rn |
!       skyrad=skyrad1,skyrad2  {3sig}
    subphot $i+20 $i rad=r0 gain=gain rn=rn |
       skyrad=skyrad1,skyrad2 oscale=exp1/exp0*gain0/gain {3sig}
    i=i+1
  end_do
end_do

! Clean up files
$'rm' {outfield}*.ap {outfield}*.grp
$gzip -f {outfield}a.coo {outfield}b.coo

! Average aperture mags
printf '{outfield}' >./average.inp
printf '%i6' nfilt >>./average.inp
string exp ' '
ntot=0
i=1
do ifilt=1,nfilt
   string filt 'filt%i1' ifilt
   string file {outfield}{{filt}}{pass}
   $mv {file}.nst {outfield}{{filt}}.nst
   string var 'filt%i1' ifilt
   printf '{{var}}' >>./average.inp
   string var 'istart%i1' ifilt
   {var}=ntot+1
   string var 'nframe%i1' ifilt
   ntot=ntot+{var}
   do iframe=1,{var}
     string scvar 'scale%i2.2' i
     exp1={scvar}
     string exp '{exp} %f7.3' exp1/exp0
     i=i+1
   end_do
   string var 'iend%i1' ifilt
   {var}=ntot
end_do

one=1
if nfilt==2
   printf '%i6 %i6 %i6 %i6 %i6' ntot istart1 istart2 iend1 iend2 |
       >>./average.inp
else_if nfilt==3
   printf '%i6 %i6 %i6 %i6 %i6 %i6 %i6' |
           ntot istart1 istart2 istart3 iend1 iend2 iend3 >>./average.inp
else_if nfilt==4
   printf '%i6 %i6 %i6 %i6 %i6 %i6 %i6 %i6 %i6' |
           ntot istart1 istart2 istart3 istart4 iend1 iend2 iend3 iend4 |
           >>./average.inp
end_if

printf '{exp}' >>./average.inp
printf '%i6' ifake >>./average.inp
$gaverage < average.inp > /dev/null


!combine
string file '{outfield}{filt1}'
string file2 '{outfield}{filt2}'
combine cmb=./{file}ap.nst dat=./{file2}ap.nst dist=0 id
combine cmb=./{file}.nst dat=./{file2}.nst dist=0 id
if (nfilt>=3)
  string file2 '{outfield}{filt3}'
  combine cmb=./{file}ap dat=./{file2}ap.nst dist=0 id
  combine cmb=./{file} dat=./{file2}.nst dist=0 id
end_if
if (nfilt>=4)
  string file2 '{outfield}{filt4}'
  combine cmb=./{file}ap dat=./{file2}ap.nst dist=0 id
  combine cmb=./{file} dat=./{file2}.nst dist=0 id
end_if

! Weed out bad observations based on CHI
weed in=./{file}.cmb out=./tmp.cmb
weed in=./{file}ap.cmb out=./tmpap.cmb

if (nfilt==2)
    string out '{outfield}ap'
    magaver norm id name={filt1},{filt2} in=./tmpap out=./{out}.mag |
     ncolor=2 meas=1,1 meas=2,2 short
    string out '{outfield}'
    magaver norm id name={filt1},{filt2} in=./tmp out=./{out}.mag ncolor=2 |
     meas=1,1 meas=2,2 short
else_if (nfilt==3)
    string out '{outfield}ap' 
    magaver norm id name={filt1},{filt2},{filt3} |
     in=./tmpap out=./{out}.mag ncolor=3 |
     meas=1,1 meas=2,2 meas=3,3 short
    string out '{outfield}'
    magaver norm id name={filt1},{filt2},{filt3} |
     in=./tmp out=./{out}.mag ncolor=3 |
     meas=1,1 meas=2,2 meas=3,3 short
else_if (nfilt==4)
    string out '{outfield}ap' 
    magaver norm id name={filt1},{filt2},{filt3},{filt4} |
     in=./tmpap out=./{out}.mag ncolor=4 |
     meas=1,1 meas=2,2 meas=3,3 meas=4,4 short
    string out '{outfield}'
    magaver norm id name={filt1},{filt2},{filt3},{filt4} |
     in=./tmp out=./{out}.mag ncolor=4 |
     meas=1,1 meas=2,2 meas=3,3 meas=4,4 short
end_if

$'rm' {outfield}*.cmb.old {outfield}*. tmp.cmb tmpap.cmb
$gzip -f {outfield}*.nst {file}.cmb {file}ap.cmb
$mkdir {outfield}
$mv {outfield}* {outfield}/

end