!procedure to analyze stellar photometry completeness runs

unalias rd
unalias wd

do iap=1,2

  ! Note we must do profile first to get proper weeding of spurious detections!
  if iap==1
    string ap ' '
  else
    string ap 'ap'
  end_if

  ! {field}.fake gives all of the input fake star counts 
  open input ./{field}.fake
  stat n=count[input]
  string field2 {field}{ap}
  $'rm' {field2}.complete
  do jfake=1,n
    read input
    ifake=@input.2
    string fake '%i3.3' ifake
    printf '%i5' nfilt >./bincomplete.inp
    do ifilt=1,nfilt
      string var 'input.%i1' ifilt+2
      printf '%f12.5' @{var} >>./bincomplete.inp
    end_do
    string file '{field}_%i3.3{ap}' ifake
    printf '{file}' >>./bincomplete.inp
    do ichip=1,4
       cnts=@input.3
       ! Make fake.nst file with counts of fake stars
       call makefake ichip cnts
   
  !  Combine output photometry list with input fake star list 
       string file '{field}%i1{ap}' ichip 
       string dir '{field}%i1_{fake}' ichip
       string dir2 '{field}%i1' ichip
  
  !  Use the first color (F555W) to determine matches. The fake star
  !    mags we will use will be the PSF magnitudes
       ifilt=1
       string filt 'filt%i1' ifilt
       string file2 '{field}fake{{filt}}%i1_%i3.3' ichip ifake
       $mv fake.nst {file2}.nst
       combine cmb={dir}/{file}.mag dat=./{file2}.nst dist=0.5,1 |
           silent out={dir}/{file}.cmb
  
  !  Combine list with original photometry list from frame with no fake stars 
       string file2 '{field}%i1{ap}' ichip 
       string var 'nstar%i1' ichip
       combine cmb={dir}/{file} dat={dir2}/{file2}.mag dist=0.5,1,2 |
            silent imag=ifilt nstar={var}
  
  !  For aperture photometry, need to correct fake star magnitudes for
  !     aperture correction effects
       if iap==2
         string trn '{field}%i1' ichip
         open trn {trn}/{trn}.trn
         read trn
         read trn
         z1=@trn.6
         read trn
         z2=@trn.6
         close trn
         open trn {trn}/{trn}ap.trn
         read trn
         read trn
         a1=@trn.6
         read trn
         a2=@trn.6
         close trn
       else
         a1=0 z1=0 a2=0 z2=0
       end_if 

       printf '%i10 %2f10.3' {var} z1-a1 z2-a2 >>./bincomplete.inp
       printf '{dir}/{file}' >>./bincomplete.inp
    end_do

    string file '{field}fake_%i3.3{ap}' ifake
    $ (cd src; make bincomplete)
    $ bin/bincomplete < bincomplete.inp > {file}.complete
    m1=-2.5*log10[@input.3]+25
    m2=-2.5*log10[@input.4]+25
    string outmag '%2f12.5' m1 m2
    $ echo "{outmag}" >>./{field2}.complete
    $ cat {file}.complete >>./{field2}.complete
    string file '{field}_%i3.3{ap}' ifake
    string dir '{field}_{fake}' ichip
    $ mkdir {dir}
    $ mv {field}*_{fake}*.* {dir}/
    do ichip=1,4
       string file '{field}%i1{ap}' ichip ifake
       string dir '{field}%i1_{fake}' ichip
       $'rm' {dir}/{file}.cmb 
       $'rm' {dir}/{file}.cmb.old {file}.
       string file '{field}fake{filt1}%i1_%i3.3' ichip ifake
       $'rm' {file}.*
       string file '{field}%i1{ap}' ichip ifake
       $'rm' {file}.
       string temp 'temp%i1_%i3.3{ap}' ichip ifake
       $'rm' {temp}.mag 
    end_do
    nextfake:
  end_do

  ! Now take individual error frames and interpolate in between them and
  !  make "error histogram" frames
  iii=0
  do i=1,n-1
    istart=i
    iend=50
    string file '{field}_%i3.3' i
    rd $istart {file}/{file}{ap} wfpc
    string file '{field}_%i3.3' i+1
    rd $iend {file}/{file}{ap} wfpc
  
    m1start={?istart:mag1}
    m1end={?iend:mag1}
    dm1=(m1end-m1start)/5.
    m2start={?istart:mag2}
    m2end={?iend:mag2}
    dm2=(m2end-m2start)/5.
    mag=m1start
    mag2=m2start
    copy 99 $iend
    sub 99 $istart
    buf=istart
  !  tv $istart l=50
    do j=2,5
      buf=buf+1
      mag=mag+dm1
      mag2=mag2+dm2
      copy $buf $istart
      copy 98 99
      mul 98 c=(j-1)/5.
      add $buf 98
      fits $buf float=m1 mag
      fits $buf float=m2 mag2
      abx $buf total=tot
      fits $buf int=nmatch int[tot]
  !    tv $buf l=50
    end_do
  !  tv $iend l=50
    if i==n-2
      buf=buf+1
      copy $buf 50
      no=6
    else
      no=5
    end_if
    buf=istart-1
    mag=m1start-dm1
    mag2=m2start-dm2
    do j=1,no
      mag=mag+dm1
      mag2=mag2+dm2
      buf=buf+1
      histlist $buf buf=98 ntot={?buf:ntot}
  ! Final output image will by ntotal*2*(n-1)*5+1
      if i==1&j==1
  !   first time
        copy 97 98
        fits 97 int=naxis 3
        fits 97 int=naxis3 (n-1)*5+1
        fits 97 int=bitpix 16
        $'rm' tempcomplete.fits
        $touch tempcomplete.fits
      end_if
  
  !   all times
      iii=iii+1
      string card 'mag1_%i3.3' iii
      fits 97 float={card} mag
      string card 'mag2_%i3.3' iii
      fits 97 float={card} mag2
  
      if i==n-1&j==no
  !  last time write out header
        wd 97 ./{field2}complete headonly
        fits 97 float=bscale 0.01
        fits 97 float=bzero 0.0
      end_if
  
      wd 98 ./junk.fits nohead notail noauto 
      $cat ./junk.fits >>./tempcomplete.fits
      
    end_do

  end_do

  $cat tempcomplete.fits >>./{field2}complete.fits
  $'rm' tempcomplete.fits


! Now make single 3D output file with the raw error data, binned to 0.02 mag
!   errors

  $'rm' temp.fits
  $touch temp.fits
  do i=1,n
    string file '{field}_%i3.3' i
    rd $i {file}/{file}{ap} wfpc
    bin $i bin=2
    if i==1
      copy 99 $i
    end_if
    mag1={i:mag1}
    mag2={i:mag2}
    string card 'mag1_%i3.3' i
    fits 99 float={card} mag1
    string card 'mag2_%i3.3' i
    fits 99 float={card} mag2
    wd $i ./junk.fits nohead notail noauto
    $cat junk.fits >> temp.fits
  end_do
  fits 99 int=crval1 -100
  fits 99 int=crval2 -100
  fits 99 int=naxis 3
  fits 99 int=naxis3 n
  fits 99 float=cdelt1 .02
  fits 99 float=cdelt2 .02
  wd 99 ./{field2}_all{ap}  headonly
  $cat temp.fits >> ./{field2}_all.fits
  $'rm' temp.fits junk.fits

end_do

END