$'rm' err.dat rms.dat
! Loop over each set of fake stars
do i=1,23

  ! Make files giving lists of all of the errors of the fake stars
  string file 'phoenix_%i3.3' i
  rd 1 {file}/{file} wfpc
  histlist 1 buf=i+100
  rotate $i+100 left
  box 1 nc=1 sc=1 nr=nr[i+100] sr=sr[i+100]
  print $i+100 box=1 >{file}/{file}.dat
  $awk 'NR>9 \{print {1:mag1}+21.823-25+0.1*rand(),sqrt($2*$2)+rand() \}' |
      {file}/{file}.dat >>./err.dat

 ! compute the rms error values 
 rmsx=0.1
 rmsy=0.1
 
 if i==1
   copy 2 1
   sub 2 2
   add 2 c=1 
   copy 3 2
   add 2 2 dc=1  ! x
   copy 4 2
   mul 4 4       ! x^2
   add 3 3 dr=1  ! y
   copy 5 3
   mul 5 5       ! y^2
 end_if
 copy 11 4
 mul 11 1
 copy 12 5
 mul 12 1
 copy 13 2
 copy 14 3
 mul 13 1
 mul 14 1

 xbar=201
 ybar=201
 rmsx=50
 rmsy=50

 sigclip=3 

 do iter=1,5
 
  nrow=sigclip*rmsy ncol=sigclip*rmsx crow=ybar ccol=xbar
  srow=max[1,crow-nrow/2]
  scol=max[1,ccol-ncol/2]
  erow=min[400,crow+nrow/2]
  ecol=min[400,ccol+ncol/2]
  nrow=erow-srow+1
  ncol=ecol-scol+1
  box 1 sr=srow nr=nrow sc=scol nc=ncol
 
  abx 11 1 total=sumx2i
  abx 12 1 total=sumy2i
  abx 13 1 total=sumxi
  abx 14 1 total=sumyi
  abx 1 1 total=sumi
 
  xbar=sumxi/sumi
  rmsx=(sumx2i-2*xbar*sumxi+xbar^2*sumi)/sumi
  rmsx=sqrt[rmsx]
  ybar=sumyi/sumi
  rmsy=(sumy2i-2*ybar*sumyi+ybar^2*sumi)/sumi
  rmsy=sqrt[rmsy]
 
 end_do
 
 mag1={?1:mag1}
 mag2={?1:mag2}
 nmatch={?1:nmatch}
 ntot={?1:ntot}
 
 printf '%i5 %5f12.4 ' i mag1-3.192 mag2-4.056 .01*rmsx .01*rmsy nmatch/ntot
 printf '%i5 %5f12.4 ' i mag1-3.192 mag2-4.056 .01*rmsx .01*rmsy nmatch/ntot |
      >>./rms.dat
 
end_do
 
END