parameter i1 i2 nointer 

! Gets shifts between images
! Sorts files into stacks for CR rejections by pointing only
! Sorts files into stacks for combination by pointing and exptime
open input ./fields.lis
stat n=count[input]
string field {input}
call {field}setup

if acs==1
  n=2
end_if
if wfpc==1
  n=4
end_if
if mosaic==1
  n=8
end_if

if i1==0
  i1=1
  i2=n
else_if i2==0
  i2=i1
end_if

$mkdir oldcomb
$mv comb_*.out oldcomb/

! Do two passes. First pass will just use header shifts and will then
!   call PICCRS for CR rejection. Second pass will use the CR-rejected
!   frames and get shifts by cross-correlation

$'rm' getshift.dat
do ipass=1,2
 if ipass==1
   i1=1 i2=1
 else
   i1=1 i2=4
 end_if
 dxtot=0 dytot=0 wtot=0
 do ifield=i1,i2
   ! read in PC frames
   call look ipass-1 ifield
   tmax=0
   ! Use longest exptime in first filter as reference
   do i=1,nframe1
     t={i:exptime}
     if (t>tmax) 
       iref=i
       tmax=t
     end_if
   end_do
   i=0
   do ifilt=1,nfilt
     string nf 'nframe%i1' ifilt
     string com 'crmask 1 3 0'
     do j=1,{nf}
       i=i+1
       if (i==iref)
         dx=0 dy=0 dxhead=0 dyhead=0 w=1
       else_if ipass==1
         call getshift iref i 0 ifield
       else
         call getshift iref i 1 ifield
	 printf '%2i3 %4f8.3 %e10.3' i ifield dx dy dxhead dyhead w >>./getshift.dat
       end_if
       if ipass==1
         call getfilt i
         string dx 'dxh%i1%i2.2' ifilt j
         string dy 'dyh%i1%i2.2' ifilt j
         {dx}=dx {dy}=dy
         string var 'name%i1%i2.2' ifilt j
         string com '{com} {{var}}_c0f 0. 0. 0. %2f8.2' dx dy
         $cp archive/{{var}}_c0f.fits ./
         $cp archive/{{var}}_c1f.fits ./
         $mask {{var}}_c0f.fits {{var}}_c1f.fits
         $getsky {{var}}_c0f 1
       else
         if ifield==1
           w=w*2
         end_if
         string dx 'dx%i1%i2.2' ifilt j
         string dy 'dy%i1%i2.2' ifilt j
         string wtot 'wtot%i1%i2.2' ifilt j
         if ifield==1
           {dx}=dx*w {dy}=dy*w {wtot}=w
         else
           {dx}={dx}+dx*w {dy}={dy}+dy*w {wtot}={wtot}+w
           if ifield==i2
             {dx}={dx}/{wtot} {dy}={dy}/{wtot}
           end_if
         end_if 
       end_if 
     end_do
     if ipass==1
       ${com}
     end_if
   end_do
 end_do
end_do

$'rm' cross.dat getshift.err
do ifield=1,1
  if ifield>=i1&ifield<=i2
   if wfpc==1|mosaic==1|acs==1
    call {field}setup ifield
   else
    call {field}setup
   end_if

   ! Load up headers (to get filter right below!)
   call look 

   i=0
   do ifilt=1,nfilt 

     do isort=1,2
       if isort==1
       ! For CR reject, sort by pointing only
         tmin=0 tmax=1e10 
         dmax=0.3
       else_if isort==2
         ! For combining , sort by pointing and exposure time
         if wfpc==1
           tmin=0.5 tmax=2.
         end_if
         dmax=0.5
       end_if
       npoint=0 
       string nf 'nframe%i1' ifilt
       do j=1,{nf}
         i=i+1
         string dx 'dx%i1%i2.2' ifilt j
         string dy 'dy%i1%i2.2' ifilt j
         dx={dx} dy={dy}
         string dx 'dxh%i1%i2.2' ifilt j
         string dy 'dyh%i1%i2.2' ifilt j
         dxhead={dx} dyhead={dy}
         if (i==1)
           dx0=dx dy0=dy dxhead0=dxhead dyhead0=dyhead
         end_if
         texp={i:exptime}
         call getfilt i
         string var 'cfilt%i1' ifilt
         string {var} {filt}
!         string filt {i:filtnam1}
         string var 'name%i1%i2.2' ifilt j
         if isort==1
           printf '%4f8.2 {filt} {i:exptime} {{var}} %i' dx-dx0 dy-dy0 |
              dxhead-dxhead0 dyhead-dyhead0 ifield >>./cross.dat
           string rout 'crmask_{filt}'
         else
           string rout 'comb_{filt}'
         end_if
         jpoint=0 new=0
         if npoint>0
           do ipoint=1,npoint
             string xvar 'pdx%i' ipoint
             string yvar 'pdy%i' ipoint
             string tvar 'pexp%i' ipoint
             if (abs[dx-{xvar}]<dmax&abs[dy-{yvar}]<dmax&texp/{tvar}>tmin&texp/{tvar}<tmax)
               jpoint=ipoint
             end_if
           end_do
         end_if
         if jpoint==0
           new=1
           npoint=npoint+1
           string var 'pdx%i' npoint
           {var}=dx
           string var 'pdy%i' npoint
           {var}=dy
           string var 'pexp%i' npoint
           {var}=texp
           jpoint=npoint
         end_if
         string pdx 'pdx%i' jpoint
         string pdy 'pdy%i' jpoint
         string out '{rout}%i' jpoint
         string var 'name%i1%i2.2' ifilt j
         if isort==2
           if new==1
             printf '{{var}}_c0f %2f8.2' dx-dx0 dy-dy0 >./{out}.out
           else
             printf '{{var}}_c0f %2f8.2' dx-dx0 dy-dy0 >>./{out}.out
           end_if
         else
           if new==1
             printf '{{var}}_c0f %2f8.2' dx-{pdx} dy-{pdy} >./{out}.out
           else
             printf '{{var}}_c0f %2f8.2' dx-{pdx} dy-{pdy} >>./{out}.out
           end_if
           if j==1
             printf '{{var}}_c0f %2f8.2' dx-dx0 dy-dy0 >./{rout}.all
           else
             printf '{{var}}_c0f %2f8.2' dx-dx0 dy-dy0 >>./{rout}.all
           end_if
         end_if
       end_do
       if isort==1
         i=i-{nf}
         ! Find out if any of the pointings only have one image, in which case we will
         !   want to CR-reject them along with another pointing
         do jpoint=1,npoint
           string var 'n%i' jpoint
           string out '{rout}%i' jpoint
           open inpoint ./{out}.out
           stat {var}=count[inpoint]
           close inpoint
         end_do
       end_if 
     end_do
   end_do

  end_if
end_do

$mkdir newcomb
$cp comb_*.out newcomb/
$diff oldcomb newcomb >./{field}.diffcomb

end