parameter ichip iflat ifake

if wfpc==1
  subclip=0.75
  if ifake==0
    cte=1
    string cte 'cte'
  else
    cte=0
    string cte ' '
  end_if
!subclip=-1
!cte=0
else
  subclip=-1
  string cte ' '
  cte=0
end_if
if mosaic==-1
  string xtype 'cubic'
else
  string xtype 'trans'
end_if

if ichip>0
  if wfpc==1|acs==1
    string outfield '{field}%i1' ichip
  else
    string outfield '{field}_%i1' ichip
  end_if
else
  string outfield {field}
end_if

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

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

! Buffers for find frames and default work buffers
work1=100 work2=99 work3=96
findbuf=109 findnorm=110
subbuf=97 subnorm=98

! 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 phot 'phot='
string gain 'gain='
string rn 'rn='
string col 'col='
string sky ' '
tot=0
do ifilt=1,nfilt
  string filt 'filt%i1' ifilt
  string var 'nframe%i1' ifilt
  string sumbuf 'sumbuf%i1' ifilt
  if {var}>1
    {sumbuf}=100+ifilt
  else
    {sumbuf}=i
  end_if
  do iframe=1,{var}
    if i>20&ip>=-9
      printf 'ERROR: multistar currently compile limited to 20 frames'
      pause
    end_if
    if i>20&ip>=-9
      printf |
          'ERROR: phot.pro and register currently compile limited to 20 frames'
      pause
    end_if

    ! get the appropriate file name
    string var2 'name%i1%i2.2' ifilt iframe
    if ichip>0
      if wfpc==1|acs==1
        string name '{{var2}}%i1' ichip
        if ifake==0
          string file '{{var2}}'
          string nim 'nim=%i1' ichip
        else
          string file {name}
        end_if
      else
        string name '{{var2}}_%i1' ichip
        string file {name}
        string nim ' '
      end_if
    else
      string name '{{var2}}'
    end_if

    ! Make links to the appropriately named PSF files
    string psf '{field}%i2.2' i
    if ip==-1
      $ln -fs {name}var.lib psf/{psf}.lib
    else_if ip>-10
      $ln -fs {name}.lib psf/{psf}.lib
    end_if

    ! Read in the data
    if ifake>0
      rd $i 'tmp/{file}' {pictype} {nim}
    else_if archival==1|acs==1
      rd $i '{file}_m0f' {pictype} {nim}
!    else_if acs==1
!      call acsread i {file} ichip
    else
      rd $i '{name}' {pictype}
    end_if

    !  Load up approximate frame offsets if appropriate coordinate files exist
    if archival==1|acs==1
     string var2 'dx%i1%i2.2' ifilt iframe
     dx1={var2}
     string var2 'dy%i1%i2.2' ifilt iframe
     dy1={var2}
     if wfpc==1
       call wfpcoff ichip dx1 dy1
     else
       dx=-dx1 dy=-dy1
     end_if
     if ishift==1
       dx=0 dy=0
     end_if
     fits $i float=x(01,1) -dx
     fits $i float=x(01,2) -dy
     do ix=2,10
       string var 'x(%i2.2,1)' ix
       fits $i float={var} 0.
       string var 'x(%i2.2,2)' ix
       fits $i float={var} 0.
     end_do
    else
     if i==1
      string ref {name}
     end_if
     err continue
     register 1 $i load dat=pos/{ref}.coo ref=pos/{name}.coo {xtype}
     dx=0 dy=0
     err continue
     dx={i:x(01,1)} 
     err continue
     dy={i:x(01,2)}
    end_if
    fits $i int=dxpos int[dx]
    fits $i int=dypos int[dy]

!  If this is WFPC2 data, flatten if we are using unflattened data
    if wfpc==1
      call wfpcgain i ichip
      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 $work1 wfpc2cal/flats/{flat}
        mul $i $work1
      end_if
    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
    string hjd 'hjd%i1' ifilt
    string hjd0 'hjd0%i1' ifilt
    err continue
    hjd $i 
    err goto noair
    a={k}*{i:am}
    a=a+{hjd}*({i:hjd}-2444000-{hjd0})
    {extvar}=10^(0.4*a)
    {scvar}={?i:exptime}*{extvar}
    air0=0
    goto haveair
    noair:
    {extvar}=1
    {scvar}={?i:exptime}
    air0=-1
    haveair:

    ! set a whole bunch of string variables for reduction
    string frames '{frames} %i2.2' i
    exp1={scvar}
    if i>1
      string exp '{exp},' 
      string col '{col},'
      string gain '{gain},'
      string rn '{rn},'
      string phot '{phot},'
    end_if
    if exp1/exp0<0.1
      string exp '{exp}%f7.5' exp1/exp0
    else_if exp1/exp0<1
      string exp '{exp}%f6.4' exp1/exp0
    else_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 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
    string photvar 'phot%i1%i2.2' ifilt iframe
    string phot '{phot}%i1' {photvar}

    ! Create find buffer findbuf consisting of all frames in filters on 
    ! which find
    ! is specified. Also create summed buffers 100+ifilt for each filter, to
    ! be used for initial aperture photometry guesses. All of these frames
    ! will be normalized to exposure time exp0

    unmask
    err continue
    ! instrument bad pixels
    get mask=/usr/tmp/blank
    masktoim $work1 nr=nr[i] sr=sr[i] nc=nc[i] sc=sc[i]
    ! get bad pixel mask for this frame
    err continue
    get mask=msk/{name}.msk
    ! Add in bad pixels, e.g. zapped CRs from PICCRS
    clip $i min=lowbad maskonly
    masktoim $work2 nr=nr[i] sr=sr[i] nc=nc[i] sc=sc[i]
    ! Combination of all sources of bad pixels will be in buffer 99
    mul $work2 $work1
  
    ! Dont include bad pixels in summed frames
    copy $work1 $i
    mul $work1 $work2

    ! Normalize to correct gain
    mul $work1 c=gain/gain0

    ! Final normalization will be by exposure time
    mul $work2 c={scvar}/exp0
 
    string find 'find%i1%i2.2' ifilt iframe
    ! Add image into find frame if requested to do so
    if (i==1)
      copy $findbuf 1
      sub $findbuf $findbuf
      copy $findnorm $findbuf
    end_if
    if {find}>=1
      ! Dont include bad pixels
      dx=0 dy=0
      err continue
      dx={i:dxpos} 
      err continue
      dy={i:dypos}
      add $findbuf $work1 dr=int[dy] dc=int[dx]
      add $findnorm $work2 dr=int[dy] dc=int[dx]
      tot=tot+{scvar}/exp0
    end_if

    ! Correct for pixel areas now. Dont do this before finding to insure that
    !  the sky stays flat...
    if wfpc==1
      string file 'pixarea_%i1' ichip
      rd $work3 wfpc2cal/distort/{file}
      mul $i $work3
      call wfpcgain i ichip
    end_if
    if acs==1
      string file 'wfc%i1_pam' acscam
      rd $work3 acscal/distort/{file}.fits.gz
      mul $i $work3
    end_if
    if abs[mosaic]==1
      if mosaic==-1
        string file 'pixarea_%i1' {i:IMAGEID}
      else
        string file 'pixarea_%i1' ichip
      end_if
      rd $work3 distort/{file}
      mul $i $work3
    end_if

    ! Make individual frames for each filter to be used for aperture phot
    string var 'tot%i1' ifilt
    if {sumbuf}~=i
      if iframe==1
        copy ${sumbuf} $i
        sub ${sumbuf} ${sumbuf}
        copy ${sumbuf}+nfilt ${sumbuf}
      end_if
      add ${sumbuf} $work1 dr=int[dy] dc=int[dx]
      add ${sumbuf}+nfilt $work2 dr=int[dy] dc=int[dx]
    end_if

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

    i=i+1
  end_do
end_do


! Normalize find frame - dont include pixels which appear in less than 1/4
!  total exposure time
!clip $findnorm min=tot/4 vmin=0
div $findbuf $findnorm
dispose $work1 $work2 

! 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 sumbuf 'sumbuf%i1' ifilt
  string var 'nframe%i1' ifilt
  if {var}>1
    div ${sumbuf} ${sumbuf}+nfilt
    dispose ${sumbuf}+nfilt
  end_if
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 $findbuf 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 $100+ifilt red/{file} maxtry=0
!end_do

! Loop over multiple photometry passes
if ip>-10
  npass=3
else
  npass=1
end_if
do ipass=1,npass

 if ipass==1
  string pass 'a'
 else_if ipass==2
  string pass 'b'
  string oldpass 'a'
 else_if ipass==3
  string pass 'c'
  string oldpass 'b'
 else_if ipass==4
  string pass 'd'
  string oldpass 'c'
 end_if
 unmask
 ! Mask out edges
 box 1 nr=10 sr=1 nc=nc[findbuf] sc=sc[findbuf]
 clip $findbuf maskonly max=-1e10 box=1
 box 1 nr=10 sr=nr[findbuf]-9 nc=nc[findbuf] sc=sc[findbuf]
 clip $findbuf maskonly max=-1e10 box=1
 box 1 nc=10 sc=1 nr=nr[findbuf] sr=sr[findbuf]
 clip $findbuf maskonly max=-1e10 box=1
 box 1 nc=10 sc=nc[findbuf]-9 nr=nr[findbuf] sr=sr[findbuf]
 clip $findbuf maskonly max=-1e10 box=1

 ! Mask out undefined regions
 clip $findnorm min=0.00001 maskonly

 if ipass==1
  ! Get sky value
  daosky $findbuf mask
  if (skysig<0)
    skysig=sqrt[sky/gain+rn^2/gain^2]
  end_if
  if ifake>0
    open thresh {outfield}/{outfield}.thresh
    read thresh
    skysig=@thresh.1  hmin=@thresh.2
    close thresh
  end_if
  
  ! zap for CRs (on find frame only)
  zap $findbuf size=3 sig=5
  zap $findbuf size=3 sig=5

  daofiles coo=./{outfield}{pass}
  ! First pass for geometric transformations only, consider just bright stars
  if ifake==0
    if ip>-10
     find $work1 $findbuf thresh={thresh1} |
      lowbad=lowbad maxbad=1e10 nwant=100,1000 hmin=hmin
    else
     find $work1 $findbuf thresh={thresh2} lowbad=lowbad maxbad=1e10 
     hmin={thresh2}
    end_if
    string dpos 'mindpos=1 header'
    printf '%2f12.4' skysig hmin >{outfield}.thresh
  else
    find $work1 $findbuf thresh=hmin |
      lowbad=lowbad maxbad=1e10 
    string dpos 'dposmax=1 header'
  end_if

!  if nofind==1
!    tv $findbuf 50.
!    mark new
!    save coo=./{outfield}{pass} sr=0 sc=0
!  end_if

  string fitphot '{phot}'
 else
  zap $subbuf size=3 sig=5
  zap $subbuf size=3 sig=5
  daofiles coo=./{outfield}sub{pass}
  if subclip>0
   ! First find new stars on subtracted frame
   copy $work1 $subbuf
   copy $work2 $findbuf
   sub $work1 c=sky+{thresh2}
   sub $work2 c=sky+{thresh2}
   clip $work1 vmin=0
   clip $work2 vmin=0
   div $work1 $work2
   clip $work1 min=subclip vmin=0
   smooth $work1 fw=fw
   clip $work1 min=0.01 vmin=0 max=0.01 vmax=1
   smooth $work1 fw=1
   copy $work2 $findbuf
   sub $work2 c=sky
   mul $work1 $work2
   add $work1 c=sky
  else
   copy $work1 $subbuf
  end_if
 
  find $work2 $work1 thresh={thresh2} lowbad=lowbad relerr=relerr
  find $work2 $work1 thresh={thresh2}*relerr lowbad=lowbad

  daofiles file=./{outfield}{filt1}{oldpass}.nst file2=./{outfield}sub{pass}.coo
  daofiles file3=./{outfield}{pass}.coo
  append append

  error goto nomaster
  if master==1
    string master 'master%i1' ichip
    error goto nomaster
    open master ../{master}.coo
    close master
    $'rm' master.coo
    $head -3 {outfield}{pass}.coo >./master.coo
    $cat ../{master}.coo >>./master.coo
    $mv master.coo {outfield}{pass}.coo
  end_if

  nomaster:
  daofiles coo=./{outfield}{pass} 
  string dpos 'dposmax=1 header'
  string fitphot ' '
!  string dpos 'header'
  dispose $work2
 end_if
 dispose $work1 

 ! 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}
  string sumbuf 'sumbuf%i1' ifilt
  photom ${sumbuf} rad=r1 skyrad=skyrad1,skyrad2 gain=gain rn=rn {3sig}
  shortap 
  ! Separation here doesn't matter, it is done again during multistar, here
  !  we just need it to create group files
  opt ig=1
  oldgroup crit=1
 end_do

 if ip<-9
   goto apdone
 end_if
 
 ! 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
 if ipass<npass
   opt ps=ps/2
 end_if
 if cte==1&ipass==npass
  daofiles psf=psf/{field} pro=./{outfield}{pass}nocte
  multistar {frames} |
   {sky} |
   {col} |
   {exp} |
   {fitphot} |
   {gain} |
   {rn} |
   {dpos} drift=5 exp0=exp0 gain0=gain0 air0=air0 write |
   skyrad=skyrad1,skyrad2 {3sig} crit=crit
  do ifilt=1,nfilt
    string filt 'filt%i1' ifilt
    string file {outfield}{{filt}}{pass}
    $mv {file}.nst {file}nocte.nst
  end_do
 end_if
 daofiles psf=psf/{field} pro=./{outfield}{pass}
 multistar {frames} |
   {sky} |
   {col} |
   {exp} |
   {fitphot} |
   {gain} |
   {rn} |
   {dpos} drift=5 exp0=exp0 gain0=gain0 air0=air0 write |
   skyrad=skyrad1,skyrad2 {3sig} {cte} crit=crit 
 opt ps=ps

 ! Get the geometric transformation if first pass
 if (ipass==1)
   if ifake==0
     register {frames} load dpos=./{outfield}{pass}.dpos silent |
      niter=10 maxrms=0.1 cubic 
   else
     open shift {outfield}/{outfield}{pass}.shift
   end_if
   $'rm' ./{outfield}{pass}.shift
   open fitdata ./{outfield}{pass}.frm
   read fitdata
   i=1
   string exp 'exp='
   do ifilt=1,nfilt
     string var 'nframe%i1' ifilt
     do iframe=1,{var}
       if (ifake>0)
         read shift
         do j=1,10
           string xvar 'X(%i2.2,1)' j
           string yvar 'X(%i2.2,2)' j
           read shift
           fits $i float={xvar} @shift.1
           fits $i float={yvar} @shift.2
         end_do
       end_if
       printf '%i3' i >>./{outfield}{pass}.shift
       do j=1,10
         string xvar 'X(%i2.2,1)' j
         string yvar 'X(%i2.2,2)' j
         printf '{i:{xvar}} {i:{yvar}}' >>./{outfield}{pass}.shift
       end_do

       string scvar 'scale%i2.2' i
       read fitdata
       photscal=@fitdata.2
       exp1={scvar}
       if i>1
         string exp '{exp},' 
       end_if
       if exp1/exp0<0.1
         string exp '{exp}%f7.5' exp1/exp0
       else_if exp1/exp0<1
         string exp '{exp}%f6.4' exp1/exp0
       else_if exp1/exp0<10
         string exp '{exp}%f5.3' exp1/exp0*photscal
       else_if exp1/exp0<100
         string exp '{exp}%f6.3' exp1/exp0*photscal
       else_if exp1/exp0<1000
         string exp '{exp}%f7.3' exp1/exp0*photscal
       else
         string exp '{exp}%f8.3' exp1/exp0*photscal
       end_if
       i=i+1
     end_do
   end_do
   if ifake>0
     close shift
   end_if
 end_if

 ! Sort the output files by magnitude, weed out bad objects and 
 !    subtract stars from each frame, and create new summed frames for
 !    finding stars
 i=1
 tot=0
 string sky 'sky='
 do ifilt=1,nfilt
  string filt 'filt%i1' ifilt
  if cte==1&ipass==npass
    string file {outfield}{{filt}}{pass}nocte
    daofiles file=./{file}.nst file2=./{file}2.nst
    sort index=4 norenum
    weed in={file}2.nst out={file}.nst
    $'rm' {file}2.nst
  end_if
  string file {outfield}{{filt}}{pass}
  daofiles file=./{file}.nst file2=./{file}2.nst
  sort index=4 norenum
  weed in={file}2.nst out={file}.nst
  $'rm' {file}2.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}
    if (wfpc==1) 
      call wfpcgain i ichip
    end_if
    if i>1
      string sky '{sky},'
    end_if
    string sky '{sky}%i2' i+20

    string pro '{outfield}{pass}%i2.2' i
    daofiles pro=./{pro}
    if wfpc==1
      sub* $i+20 square {cte} gain=gain
    else
      sub* $i+20
    end_if

    string var2 'name%i1%i2.2' ifilt iframe
    if ichip>0
      if wfpc==1|acs==1
        string name '{{var2}}%i1' ichip
      else
        string name '{{var2}}_%i1' ichip
      end_if
    else
      string name '{{var2}}'
    end_if
    if ipass==npass
      wd $i+20 tmp/{name}s
    end_if

    string find 'find%i1%i2.2' ifilt iframe
    ! Add image into find frame if requested to do so
    if (i==1)
      copy $subbuf 1
      sub $subbuf $subbuf
      copy $subnorm $subbuf
    end_if
    if {find}>=ipass+1
      unmask
      err continue
      ! instrument bad pixels
      get mask=/usr/tmp/blank
      masktoim $work1 nr=nr[i] sr=sr[i] nc=nc[i] sc=sc[i]
      ! get bad pixel mask for this frame
      err continue
      get mask=msk/{name}.msk
      ! Add in bad pixels, e.g. zapped CRs from PICCRS
      clip $i min=lowbad maskonly
      masktoim $work2 nr=nr[i] sr=sr[i] nc=nc[i] sc=sc[i]
      ! Combination of all sources of bad pixels will be in buffer 99
      mul $work2 $work1
  
      ! Dont include bad pixels in summed frames
      copy $work1 $i+20
      mul $work1 $work2

      ! Normalize to correct gain
      mul $work1 c=gain/gain0

      ! Final normalization will be by exposure time
      mul $work2 c={scvar}/exp0

      ! Dont include bad pixels
      dx=0 dy=0
      err continue
      dx={i:dxpos}
      err continue
      dy={i:dypos}
      add $subbuf $work1 dr=int[dy] dc=int[dx]
      add $subnorm $work2 dr=int[dy] dc=int[dx]
    end_if


    i=i+1
  end_do
 end_do

 div $subbuf $subnorm


! Get PSF residuals in first pass, full subtraction and aperture photometry
!  on last pass
 if (wfpc==1&ipass==1)|ipass==npass
   ! 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 (wfpc==1)
         call wfpcgain i ichip
       end_if

       string pro '{outfield}{pass}%i2.2' i
       daofiles pro=./{pro}
!      zap $i+20 sig=10 size=7 gain=gain rn=rn noisemod

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

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

         string file '{outfield}_%i2.2appsf' i
         daofiles file=./{file}.nst
         subphot $i+20 $i rad=r1 gain=gain rn=rn psfcor cte |
          skyrad=skyrad1,skyrad2 oscale=exp1/exp0*gain0/gain {3sig} {cte}
       else_if psfres==1
         rd 201 psf/{psf}.lib nim=1 silent
         nsamp={201:nlib}
         subphot $i+20 $i rad=r1 gain=gain rn=rn |
           skyrad=skyrad1,skyrad2 oscale=exp1/exp0*gain0/gain {3sig} {cte} |
           res=202 nsamp=nsamp errmax=0.05
         fits 201 int=naxis3 25
         wd 201 psf/{psf}restmp.lib headonly
         do ilib=1,25
           rd 201 psf/{psf}.lib nim=ilib silent
           add 201 202 dc=64 dr=64
           wd 201 psf/res.lib nohead notail full
           $cat psf/res.lib >> psf/{psf}restmp.lib
         end_do

         do ilib=1,25
           rd 201 psf/{psf}.lib nim=ilib+25 silent
           wd 201 psf/res.lib nohead notail full
           $cat psf/res.lib >> psf/{psf}restmp.lib
         end_do
         $mv psf/{psf}restmp.lib psf/{psf}res.lib
         ! set things up to use the new libraries in next iteration
         $ln -fs {psf}res.lib psf/{psf}.lib

!         copy $i+40 $i
!         daofiles psf=psf/{psf}restmp.lib
!         sub* $i+40 square
!         daofiles psf=psf/{psf}.lib
       end_if

       i=i+1
     end_do
   end_do
!   pause
 end_if


 apdone:
! Done with this pass
end_do

if ip<-9
  do ifilt=1,nfilt
   string filt 'filt%i1' ifilt
   string file {outfield}{{filt}}{pass}
   $mv {file}.nst {outfield}{{filt}}.nst
  end_do
  goto apdone2
end_if

! Clean up files
$'rm' {outfield}*.ap {outfield}*.grp
$gzip -f {outfield}?.coo {outfield}sub?.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
   $mv {file}nocte.nst {outfield}{{filt}}nocte.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
printf '%i6' cte >>./average.inp
$(cd src; make gaverage)
$bin/gaverage < average.inp > /dev/null

!Create a .nst for big aperture photometry
do ifilt=1,nfilt
  string filt 'filt%i1' ifilt
  $head -3 {outfield}{{filt}}ap.nst >./{outfield}{{filt}}ap0.nst
  $awk 'NR>3 \{printf("%6d%9.2f%9.2f%9.3f%9.3f%9.3f%9.1f%9.2f%9.3f\n",$1,$2,$3,$4-$8>99?99.999:$4-$8,$5,$6,$7,1.,$9)\}' {outfield}{{filt}}ap.nst >>./{outfield}{{filt}}ap0.nst
  if cte==1
    $head -3 {outfield}{{filt}}ap.nst >./{outfield}{{filt}}ap0nocte.nst
    $awk 'NR>3 \{printf("%6d%9.2f%9.2f%9.3f%9.3f%9.3f%9.1f%9.2f%9.3f\n",$1,$2,$3,$4-$8>99?99.999:$4-$8,$5,$6,$7,1.,$9)\}' {outfield}{{filt}}apnocte.nst >>./{outfield}{{filt}}ap0nocte.nst
  end_if
end_do

apdone2:

!combine
string file '{outfield}{filt1}'
string file2 '{outfield}{filt2}'
combine cmb=./{file}.nst dat=./{file2}.nst dist=0 id
if ip>-9
  combine cmb=./{file}ap.nst dat=./{file2}ap.nst dist=0 id
  combine cmb=./{file}appsf.nst dat=./{file2}appsf.nst dist=0 id
  combine cmb=./{file}ap0.nst dat=./{file2}ap0.nst dist=0 id
  if cte==1
    combine cmb=./{file}nocte.nst dat=./{file2}nocte.nst dist=0 id
    combine cmb=./{file}apnocte.nst dat=./{file2}apnocte.nst dist=0 id
    combine cmb=./{file}appsfnocte.nst dat=./{file2}appsfnocte.nst dist=0 id
    combine cmb=./{file}ap0nocte.nst dat=./{file2}ap0nocte.nst dist=0 id
  end_if
end_if
if (nfilt>=3)
  string file2 '{outfield}{filt3}'
  combine cmb=./{file} dat=./{file2}.nst dist=0 id
  if ip>-9
    combine cmb=./{file}ap dat=./{file2}ap.nst dist=0 id
    combine cmb=./{file}appsf dat=./{file2}appsf.nst dist=0 id
    combine cmb=./{file}ap0 dat=./{file2}ap0.nst dist=0 id
    if cte==1
      combine cmb=./{file}nocte dat=./{file2}nocte.nst dist=0 id
      combine cmb=./{file}apnocte dat=./{file2}apnocte.nst dist=0 id
      combine cmb=./{file}appsfnocte dat=./{file2}appsfnocte.nst dist=0 id
      combine cmb=./{file}ap0nocte dat=./{file2}ap0nocte.nst dist=0 id
    end_if
  end_if
end_if
if (nfilt>=4)
  string file2 '{outfield}{filt4}'
  combine cmb=./{file} dat=./{file2}.nst dist=0 id
  if ip>-9
    combine cmb=./{file}ap dat=./{file2}ap.nst dist=0 id
    combine cmb=./{file}appsf dat=./{file2}appsf.nst dist=0 id
    combine cmb=./{file}ap0 dat=./{file2}ap0.nst dist=0 id
    if cte==1
      combine cmb=./{file}nocte dat=./{file2}nocte.nst dist=0 id
      combine cmb=./{file}apnocte dat=./{file2}apnocte.nst dist=0 id
      combine cmb=./{file}appsfnocte dat=./{file2}appsfnocte.nst dist=0 id
      combine cmb=./{file}ap0nocte dat=./{file2}ap0nocte.nst dist=0 id
    end_if
  end_if
end_if

! Weed out bad observations based on CHI
if ip>-9
  weed in=./{file}.cmb out=./tmp.cmb
  ! Aperture photometry files dont have CHI!
  !weed in=./{file}ap.cmb out=./tmpap.cmb
  $cp {file}ap.cmb tmpap.cmb
  $cp {file}ap0.cmb tmpap0.cmb
  $cp {file}appsf.cmb tmpappsf.cmb
  if cte==1
    weed in=./{file}nocte.cmb out=./tmpnocte.cmb
    $cp {file}apnocte.cmb tmpapnocte.cmb
    $cp {file}ap0nocte.cmb tmpap0nocte.cmb
    $cp {file}appsfnocte.cmb tmpappsfnocte.cmb
  end_if
else
  $cp {file}.cmb tmp.cmb
end_if

if (nfilt==2)
    string out '{outfield}'
    magaver norm id name={filt1},{filt2} in=./tmp out=./{out}.mag ncolor=2 |
     meas=1,1 meas=2,2 short outgain=gain0
    if ip>-9
      string out '{outfield}ap'
      magaver norm id name={filt1},{filt2} in=./tmpap out=./{out}.mag |
       ncolor=2 meas=1,1 meas=2,2 short outgain=gain0
      string out '{outfield}appsf'
      magaver norm id name={filt1},{filt2} in=./tmpappsf out=./{out}.mag |
       ncolor=2 meas=1,1 meas=2,2 short outgain=gain0
      string out '{outfield}ap0'
      magaver norm id name={filt1},{filt2} in=./tmpap0 out=./{out}.mag |
       ncolor=2 meas=1,1 meas=2,2 short outgain=gain0
      if cte==1
        string out '{outfield}nocte'
        magaver norm id name={filt1},{filt2} in=./tmpnocte out=./{out}.mag |
         ncolor=2 meas=1,1 meas=2,2 short outgain=gain0
        string out '{outfield}apnocte'
        magaver norm id name={filt1},{filt2} in=./tmpapnocte out=./{out}.mag |
         ncolor=2 meas=1,1 meas=2,2 short outgain=gain0
        string out '{outfield}appsfnocte'
        magaver norm id name={filt1},{filt2} in=./tmpappsfnocte out=./{out}.mag |
         ncolor=2 meas=1,1 meas=2,2 short outgain=gain0
        string out '{outfield}ap0nocte'
        magaver norm id name={filt1},{filt2} in=./tmpap0nocte out=./{out}.mag |
         ncolor=2 meas=1,1 meas=2,2 short outgain=gain0
      end_if
    end_if
else_if (nfilt==3)
    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 outgain=gain0
    if ip>-9
      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 outgain=gain0
      string out '{outfield}appsf' 
      magaver norm id name={filt1},{filt2},{filt3} |
       in=./tmpappsf out=./{out}.mag ncolor=3 |
       meas=1,1 meas=2,2 meas=3,3 short outgain=gain0
      string out '{outfield}ap0' 
      magaver norm id name={filt1},{filt2},{filt3} |
       in=./tmpap0 out=./{out}.mag ncolor=3 |
       meas=1,1 meas=2,2 meas=3,3 short outgain=gain0
      if cte==1
        string out '{outfield}nocte'
        magaver norm id name={filt1},{filt2},{filt3} |
         in=./tmpnocte out=./{out}.mag ncolor=3 |
         meas=1,1 meas=2,2 meas=3,3 short outgain=gain0
        string out '{outfield}apnocte' 
        magaver norm id name={filt1},{filt2},{filt3} |
         in=./tmpapnocte out=./{out}.mag ncolor=3 |
         meas=1,1 meas=2,2 meas=3,3 short outgain=gain0
        string out '{outfield}appsfnocte' 
        magaver norm id name={filt1},{filt2},{filt3} |
         in=./tmpappsfnocte out=./{out}.mag ncolor=3 |
         meas=1,1 meas=2,2 meas=3,3 short outgain=gain0
        string out '{outfield}ap0nocte' 
        magaver norm id name={filt1},{filt2},{filt3} |
         in=./tmpap0nocte out=./{out}.mag ncolor=3 |
         meas=1,1 meas=2,2 meas=3,3 short outgain=gain0
      end_if
    end_if
else_if (nfilt==4)
    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 outgain=gain0
    if ip>-9
      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 outgain=gain0
      string out '{outfield}appsf' 
      magaver norm id name={filt1},{filt2},{filt3},{filt4} |
       in=./tmpappsf out=./{out}.mag ncolor=4 |
       meas=1,1 meas=2,2 meas=3,3 meas=4,4 short outgain=gain0
      string out '{outfield}ap0' 
      magaver norm id name={filt1},{filt2},{filt3},{filt4} |
       in=./tmpap0 out=./{out}.mag ncolor=4 |
       meas=1,1 meas=2,2 meas=3,3 meas=4,4 short outgain=gain0
      if cte==1
        string out '{outfield}nocte'
        magaver norm id name={filt1},{filt2},{filt3},{filt4} |
         in=./tmpnocte out=./{out}.mag ncolor=4 |
         meas=1,1 meas=2,2 meas=3,3 meas=4,4 short outgain=gain0
        string out '{outfield}apnocte' 
        magaver norm id name={filt1},{filt2},{filt3},{filt4} |
         in=./tmpapnocte out=./{out}.mag ncolor=4 |
         meas=1,1 meas=2,2 meas=3,3 meas=4,4 short outgain=gain0
        string out '{outfield}appsfnocte' 
        magaver norm id name={filt1},{filt2},{filt3},{filt4} |
         in=./tmpappsfnocte out=./{out}.mag ncolor=4 |
         meas=1,1 meas=2,2 meas=3,3 meas=4,4 short outgain=gain0
        string out '{outfield}ap0nocte' 
        magaver norm id name={filt1},{filt2},{filt3},{filt4} |
         in=./tmpap0nocte out=./{out}.mag ncolor=4 |
         meas=1,1 meas=2,2 meas=3,3 meas=4,4 short outgain=gain0
      end_if
    end_if
end_if

! Clean up files!
$'rm' {outfield}*.cmb.old {outfield}*. 
$'rm' tmp.cmb tmpap.cmb tmpap0.cmb tmpappsf.cmb 
$gzip -f {outfield}*.nst {file}.cmb {file}ap.cmb {file}ap0.cmb {file}appsf.cmb
if cte==1
  $'rm' tmpnocte.cmb tmpapnocte.cmb tmpap0nocte.cmb tmpappsfnocte.cmb 
  $gzip -f {file}nocte.cmb {file}apnocte.cmb {file}ap0nocte.cmb {file}appsfnocte.cmb
end_if
$gzip -f {outfield}*.dpos
if ifake>0
  string outdir '{outfield}_%i3.3' ifake
  do ifilt=1,nfilt
    string var 'nframe%i1' ifilt
    do iframe=1,{var}
      string var2 'name%i1%i2.2' ifilt iframe
      if ichip>0
        string name '{{var2}}%i1' ichip
      else
        string name '{{var2}}'
      end_if
      $'rm' tmp/{name}.*
    end_do
  end_do
else
  string outdir '{outfield}'
end_if
$mkdir {outdir}
if ifake==0
  $mv {outfield}*.nst.gz {outdir}/
  $mv {outfield}*.coo.gz {outdir}/
  $mv {outfield}*.cmb.gz {outdir}/
else
  $'rm' {outfield}*.nst.gz 
  $'rm' {outfield}*.coo.gz 
  $'rm' {outfield}*.cmb.gz 
end_if
$mv {outfield}*.mag {outdir}/
$mv {outfield}*.dpos.gz {outdir}/
$mv {outfield}*.frm {outdir}/
$mv {outfield}*.thresh {outdir}/
$mv {outfield}*.shift {outdir}/

end