#!/bin/csh
# a third attempt at making a Supernova search analysis program
#
# Input the file_id of the target .fts, .coo, and .cand file
# 

# The script does the following:
#   0. copies images to scratch files and repairs chip defects
#   1. determines the relative offset from fiducial stars  
#   2. determines the brightness scaling factor from the fiducial stars
#   3. shifts the target  to align the stars with the reference
#   4. smooths reference image if necessary
#  	5. scales the reference image  so that the star brightnesses are the same
#   6. subtracts the reference from the target
#   7. makes a search for stars
#   8. filters for a. too weak , b. fwhm too little, c. same coords as star

# An output file  is created
#      1. with the root name of the target with '.cand'
#      2. the header of the target.coo file
#      3. keywords describing the subtraction
#      4. photometric points of the residual stars 
# 12/23/92 - localize the names so two process can run in the same directory
# 1/29/93 - use 'cosmic' program to sort out small fwhm stars
# 4/1/93 - REF_DIR changed to use value in ait.config
# 5/5/93 - defect file now depends on the target
# 5/7/93 - put correct path for defect file 
# 5/13/93- increase the radius to 4 in original findmatch
# 7/13/93- explicitly put the arguments for cosmic program raise fwhm 2->2.25
# 7/27/93 - separate the locating from the subtraction
# 3/4/94 - replace sed with getkey
# 3/29/94 - candidates with skysigma > 100 rejected.
# 4/12/94 - skip candidates in the bright regions of the reference
# 8/19/94 - REF_DIR --> sn/ref
# 1/30/97 default for defect and galbox

if ($#argv < 1) then
	echo usage: sn_sub target 
	exit 1
endif

set target=$argv[1]
set REF_DIR=`picksym data_dir`/sn/ref
set ref=$REF_DIR/`picksym file=$target.fts object`
set DEFECT=$REF_DIR/`picksym file=$ref.ref default=nodefect.pht defect `

#preliminary definitions
set REF=ref$$			# local image scratch files
set DIFF=diff$$

set TMP=tmp$$
set TMP2=tmp$$2
set TMP3=tmp$$3
set TMP4=tmp$$4
set TARGET=target$$
set BETTER=better$$
set BOTH=both$$

set APER_SIZE=8				# photometry aperture size

#NOTE the smoothing has been disabled
set SMOOTH_THRESH=15.0		# minimum blur to smooth reference
#set SMOOTH_THRESH=4.0		# minimum blur to smooth reference
set MAX_MAG=18.5 			# faintest fiducial star
set MAX_FIDUC=15			# maximum number of fiducial stars
set MAX_CAND=10				# maximum number of candidates

#copy the original data into scratch files
cp $target.fts $TARGET.fts
cp $ref.fts $REF.fts
chmod +w $TARGET.fts $REF.fts

set dr=`picksym file=$target.cand dr`
set dc=`picksym file=$target.cand dc`
set magdiff=`picksym file=$target.cand magdiff`

#remove the sky term from reference image
#  the target subtraction can be commented out if desired
sub $TARGET const=`picksym file=$target.coo sky`
sub $REF const=`picksym file=$ref.ref sky`

#derive image shift and move the image file to match reference
imshift $TARGET dr=$dr dc=$dc 

#smooth the reference image if needed 
set tfw=`picksym file=$target.cand mfwhm default=1` 
set rfw=`picksym file=$ref.ref mfwhm default=20` 
set hyp=`val "(($tfw^2)-($rfw^2))"`
if (` val "$hyp > $SMOOTH_THRESH^2" `) then
	set fw=`val "sqrt($hyp)"`
	echo smoothing reference by fw=$fw
	smooth $REF fw=$fw
	if ($status == 0) then
		putsym file=$target.cand refsmooth=$fw
	endif
endif

#derive scale factor (scaled by 1000) and scale reference to data
set FACTOR=1000
set scale=`val "10^($magdiff/2.5)" format=%.2f`

mul $REF const=$scale\*$FACTOR div=$FACTOR

#subtract the reference image from the target
sub $TARGET $REF  outfile=$DIFF
if ($status) then
	exit 1
endif

#locate any stars within the full frame but shifted box
box 10 sr=10 nr=495 sc=10 nc=495 > /dev/null
boxshift 10 dr=$dr dc=$dc edge
sn_finder $DIFF box=10

#Filter  the .coo file for cosmic rays
#using the assumption of larger fwhm 
cosmic $DIFF infile=$DIFF.coo inverse sky=0 fwhm=2.25 > $BETTER.coo

#If there are any candidates remaining
#  1. do photometry on the candidates
#  2. remove any that are near stars in the reference or Chip Defects
#  3. sort in order of brightness
#  4. strip off the non-photometric ones
#  5. add the new candidates  to the output file
#  6. take those with sky sigma < 100.
#  6. take only the MAX_CAND brightest 

#calculate the reference photometry
subheader $ref.ref | phot $ref aper=$APER_SIZE > $REF.pht

if (-z $BETTER.coo == 0) then
	cat $BETTER.coo | phot $DIFF aper=$APER_SIZE > $BETTER.pht
	cat $REF.pht $DEFECT > $BOTH.pht
	findmatch $BETTER.pht $BOTH.pht noshift radius=4 blind=1.4 inverse  | \
		sort -n +5  |			\
		awk '($6 < 20) { print $0}'		| \
		awk '($5 < 100) { print $0}'		| \
		head -$MAX_CAND > $TMP2
	if (-z $TMP2 == 0) then
		#now check the photometry of the same locations in the reference image
		#and select only those targets which are brighter than
		#2.5 magnitudes below the reference
		# note a fixed (zero) sky value is used
		cat $TMP2 | phot $ref aper=$APER_SIZE sky=0 | awk '{print $6}' > $TMP3
		stderr BUG HERE WRONG SKY VALUE
		exit 1
		paste $TMP2 $TMP3  | \
			awk '( $6 < $8+2.5) {print $0}' | cut -f1 > $TMP4
		#we would like to paste on the fwhm so we do the following nonsense
		findmatch $DIFF.coo $TMP4 noshift list | \
		awk '{print $5}' | \
		paste $TMP4 - >> $target.cand
	endif
endif


# install the reference galaxy box and  the new box corrected for clipping 
set galbox=`picksym file=$ref.ref default="0 0 0 0" galbox` 
let box11=$galbox
boxshift 11 other=10
putsym file=$target.cand galbox=\'"$galbox"\'
putsym file=$target.cand newbox=\'"`let box11`"\'

#for diagnosis calculate the rejection ratio
cat $REF.pht | phot $DIFF aper=$APER_SIZE >$TMP
set frej=`findmatch $TMP $REF.pht noshift mag`
#echo rejection $frej
putsym file=$target.cand rejmag=`echo $frej | getkey magdiff`

#finally put a label for the candidates
putsym file=$target.cand comment="id row   col     sky  skysig mag   magerr    fwhm"

#store  some things for external inspection and clean up!
cp $target.cand diff.cand
mv $TARGET.fts target.fts
mv $REF.fts ref.fts
mv $DIFF.fts diff.fts
rm $TARGET* $REF* $DIFF* $BETTER* $BOTH* $TMP*

exit 0