#!/bin/csh -f
#
# Script to automatically align CCD images with IRAF and imstar.
# Requires WCSTOOLS for image coregistration (i.e. alignment). 
# DS9 is also useful but not required.
#
# GSS 2006.05.19-21, split from reduce.csh 2006.09.23
#
# to do: cosmic ray randling, better selection of nstatpix and minrad

set GAIN=1.
set RDNOISE=0.
# settings for IMSTAR for star detection
set minrad=2
set nstatpix=200
# number of iterations for IRAF to use when calculating shifts
set niter=10
# pixel rejection algorithm for combined images
set reject="none"
set pclip=-0.6

if !( -e login.cl ) then
 echo "login.cl not found. Run mkiraf and edit login.cl appropriately."
 exit
endif

if !( -X cl ) then
 echo "IRAF not found in path, aborting."
 exit
endif

if !( -e obj.lst ) then
  ls *.{f,F}{i,I}{t,T}{,s,S} | egrep -v "tmp|bias|dark|flat|Zero|Dark|Flat|BPM|ref|test|average|median|sum" > obj.lst
endif

\rm -f A_*.fits aligned.lst >& /dev/null

if !( -X imstar ) then
 echo "Cannot find IMSTAR in path, so cannot do realignment. Aborting."
 exit
endif

if !( -e ref.fits ) then
  echo -n "ref.fits not found, analyzing images; "
cl >& /dev/null << EOF 

imutil
imstat @obj.lst fields=image,kurtosis format=no > imstat.lst
logout

EOF
  set ref=`cat imstat.lst | sort -g -k 2 | awk '{ print $1}' | head -1`
  echo "using ${ref} as reference image."
  cp ${ref} ref.fits
endif

echo "Calculating shifts for coregistration (alignment) (may take some time)."
imstar ref.fits minrad=$minrad nstatpix=$nstatpix | egrep -v IMAGE | sort -n -k 4 -r | awk '{ if ( NF==7 ) { a=5; b=6 } else { a=1; b=2 }; print $a, $b }' > coords.lst
\rm -f shifts.lst to_shift.lst no_stars.lst
set x=`head -1 coords.lst | awk '{ print $1 }'`
set y=`head -1 coords.lst | awk '{ print $2 }'`
echo "    (Reference coordinates: ${x} ${y}.)"
foreach f (`cat obj.lst`)
  echo -n "     $f "
  set shift=`imstar $f minrad=$minrad nstatpix=$nstatpix |& egrep -v "IMAGE|ListStars" | awk -v X=${x} -v Y=${y} '{ if ( NF==7 ) { a=5; b=6 } else { a=1; b=2 }; print ((X-$a)^2)+((Y-$b)^2), X-$a, Y-$b}' | sort -g | awk '{ print $2, $3 }'| head -1`
  if ( $#shift ) then 
    echo $shift | tee -a shifts.lst
    echo $f >>! to_shift.lst
  else
    echo "< no stars found, check image >"
    echo $f >>! no_stars.lst
  endif
end
cat to_shift.lst | sed -e 's/^/A_/g' > aligned.lst

echo -n "Realigning images with IRAF ..."
# blank line before IRAF commands required (cl bug?)
cl >& /dev/null << EOF

imalign ("@to_shift.lst", referenc="ref.fits", coords="coords.lst",\
   output="@aligned.lst", shifts="shifts.lst", boxsize=11,\
   bigbox=31, negative=no, background=INDEF, lower=INDEF, upper=INDEF,\
   niterate=${niter}, tolerance=0, maxshift=INDEF, shiftimages=yes,\
   interp_type="spline3", boundary_typ="nearest", constant=0., trimimages=no,\
   verbose=no)

# Removing physical-image shift fields from header (LTV1,LTV2)
hedit ("@aligned.lst", "LTV1,LTV2", "", add=no, addonly=no, delete=yes,\
 verify=no, show=no, update=yes)

imcomb ("@aligned.lst",\
   "median", headers="", bpmasks="", rejmasks="", nrejmasks="", expmasks="",\
 sigmas="",logfile="/dev/null",combine="median",reject="${reject}",project=no,\
   outtype="real", outlimits="", offsets="none", masktype="none", maskvalue=0.,\
   blank=0., scale="none", zero="none", weight="none", statsec="", expname="",\
   lthreshold=INDEF, hthreshold=INDEF, nlow=1, nhigh=1, nkeep=1, mclip=yes,\
   lsigma=3., hsigma=3.,\
   rdnoise="${RDNOISE}",gain="${GAIN}",\
   snoise="0.", sigscale=0.1, pclip=${pclip}, grow=0.)

imcomb ("@aligned.lst",\
   "average", headers="", bpmasks="", rejmasks="", nrejmasks="", expmasks="",\
 sigmas="",logfile="/dev/null",combine="average",reject="${reject}",project=no,\
   outtype="real", outlimits="", offsets="none", masktype="none", maskvalue=0.,\
   blank=0., scale="none", zero="none", weight="none", statsec="", expname="",\
   lthreshold=INDEF, hthreshold=INDEF, nlow=1, nhigh=1, nkeep=1, mclip=yes,\
   lsigma=3., hsigma=3.,\
   rdnoise="${RDNOISE}",gain="${GAIN}",\
   snoise="0.", sigscale=0.1, pclip=${pclip}, grow=0.)

imcomb ("@aligned.lst",\
   "sum",     headers="", bpmasks="", rejmasks="", nrejmasks="", expmasks="",\
  sigmas="",logfile="/dev/null",combine="sum",   reject="${reject}",project=no,\
   outtype="real", outlimits="", offsets="none", masktype="none", maskvalue=0.,\
   blank=0., scale="none", zero="none", weight="none", statsec="", expname="",\
   lthreshold=INDEF, hthreshold=INDEF, nlow=1, nhigh=1, nkeep=1, mclip=yes,\
   lsigma=3., hsigma=3.,\
   rdnoise="${RDNOISE}",gain="${GAIN}",\
   snoise="0.", sigscale=0.1, pclip=${pclip}, grow=0.)

logout
EOF

echo " done."

if !( -d unaligned ) mkdir unaligned
foreach f (`cat obj.lst`) 
 \mv -f  $f unaligned
end

echo -n "Generating phot.lst ..."
imstar average.fits minrad=$minrad nstatpix=$nstatpix | egrep -v IMAGE | awk '{ if ( NF==7 ) { a=5; b=6 } else { a=1; b=2 }; print $a, $b }' > phot.lst
cat phot.lst | awk '{ print "circle", $1, $2, "15" }' > ds9.reg
echo " done."

cat << EOF

----------------------------------------------------------------------------

Reduction complete.  DS9 will now be started (if available) and  display the 
average image.  Stars  automatically detected  for photometry will be marked
with circles using the region file "ds9.reg". 

The list of stars detected for photometry is in the file "phot.lst". Compare
this to the image and  add or delete stars as required. Also, make sure that
your main target star is the first on the list (the stars are sorted more or 
less by brightness).

Once everything is OK run phot.csh.

EOF
if ( -X ds9 ) then
 ds9 -zscale -linear average.fits -view image yes -region ds9.reg &
else
 echo "(DS9 not found.)"
endif
echo -n ""
