sources (version 2.4)
You can also look at:
#! /usr/bin/ksh
###########################################################################
#
# SYNTAX: sources
#
# BRIEFLY: Find sources in the observed fields of view
#
# DESCRIPTION: This routine finds X-ray sources from the summed images
# DESCRIPTION: and exposure maps.
# DESCRIPTION: The positions and extraction
# DESCRIPTION: radii for the sources are determined in a four step
# DESCRIPTION: process:
# DESCRIPTION: <OL>
# DESCRIPTION: <LI> Smooth the summed sky and exposure maps
# DESCRIPTION: with the <TT>vsmooth</TT> STOOL. This does a variable size
# DESCRIPTION: boxcar smooth using the minimum square box
# DESCRIPTION: containing at least 50 events. An identical filter
# DESCRIPTION: is applied to the counts image and to the exposure map.
# DESCRIPTION: The results are divided to produce a smoothed,
# DESCRIPTION: exposure-corrected image. Pixels with 0 exposure are
# DESCRIPTION: assigned a pixel value of 0. Note that
# DESCRIPTION: each combination of instrument (SIS or GIS) and
# DESCRIPTION: image resolution (256 or possibly 64 for the GIS), has three
# DESCRIPTION: images corresponding to "lo" (<2keV), "hi" (>2keV) and "all"
# DESCRIPTION: energy bands.
# DESCRIPTION: <P>
# DESCRIPTION: <LI> Detect sources in the smoothed images.
# DESCRIPTION: This is done using the <TT>ascasource</TT> STOOL
# DESCRIPTION: which reports the
# DESCRIPTION: center pixel of each source along with two radii.
# DESCRIPTION: The first is the minimum radius which will include the
# DESCRIPTION: entire source, and the second is the maximum
# DESCRIPTION: radius which will avoid confusion with other sources and
# DESCRIPTION: the background. Sources are detected with the following
# DESCRIPTION: algorithm:
# DESCRIPTION: <OL>
# DESCRIPTION: <LI> Locate the 20 brightest local peaks in the image.
# DESCRIPTION: A local peak is a pixel which is as bright or brighter
# DESCRIPTION: than the eight pixels surrounding it.
# DESCRIPTION: <LI> Examine the peaks to see if they are distinct
# DESCRIPTION: and significant in the following way.
# DESCRIPTION: Choose the brightest peak. For each of the
# DESCRIPTION: other peaks, compare
# DESCRIPTION: the intensity of the lower peak
# DESCRIPTION: to the minimum intensity along the
# DESCRIPTION: line segment connecting the peaks.
# DESCRIPTION: Delete the lower peak
# DESCRIPTION: if it is less than 2.5 times brighter
# DESCRIPTION: than the trough separating it from the brighter peak.
# DESCRIPTION: Repeat this for the second brightest surviving peak,
# DESCRIPTION: etc.
# DESCRIPTION: <LI> Determine the maximum extraction radius.
# DESCRIPTION: For each source calculate the distance to
# DESCRIPTION: the lowest point along a straight line to each other
# DESCRIPTION: source. Set the maximum radius to the nearest
# DESCRIPTION: such trough.
# DESCRIPTION: <LI> Determine the radial distribution of azimuthally
# DESCRIPTION: averaged intensity around each source.
# DESCRIPTION: Pixels with zero intensity are not included in the
# DESCRIPTION: average. The lowest intensity in the radial
# DESCRIPTION: distrubution is taken as the background flux.
# DESCRIPTION: <LI> Delete all sources whose peak intensity is less than
# DESCRIPTION: four times the background intensity.
# DESCRIPTION: <LI> Determine the "suggested" extraction radius.
# DESCRIPTION: This is the minimum radius at which the
# DESCRIPTION: azimuthally averaged intensity is more than 1/5
# DESCRIPTION: of the peak intensity and more than twice the
# DESCRIPTION: background intensity. The suggested radius is
# DESCRIPTION: never greater than the maximum radius.
# DESCRIPTION: </OL>
# DESCRIPTION: <P>
# DESCRIPTION: <LI> Determine the extraction radius for each source.
# DESCRIPTION: The default extraction radius is
# DESCRIPTION: 6 arc minutes for GIS and 4 arc minutes for SIS.
# DESCRIPTION: If the suggested radius is less than the default radius
# DESCRIPTION: (true for point sources), use the default radius.
# DESCRIPTION: Otherwise, use the suggested radius. In no case is the
# DESCRIPTION: extraction radius larger than the maximum radius.
# DESCRIPTION: <P>
# DESCRIPTION:
# DESCRIPTION: <LI> Reconcile the sources detected in the three energy bands
# DESCRIPTION: and eliminate duplicates.
# DESCRIPTION: This is done by giving a list of all source positions and
# DESCRIPTION: radii to the <TT>sift_sources</TT> STOOL.
# DESCRIPTION: The sources are listed in order of decreasing
# DESCRIPTION: intensity with all of the "all" energy band sources
# DESCRIPTION: listed first, followed by the "hi" energy sources, and the
# DESCRIPTION: "lo" energy sources. <TT>sift_sources</TT> does the
# DESCRIPTION: following:
# DESCRIPTION: <OL>
# DESCRIPTION: <LI> Take the first source in the list.
# DESCRIPTION: <LI> Examine all of the other sources in the list.
# DESCRIPTION: <LI> If the second source is within the extraction radius
# DESCRIPTION: of the first source, eliminate it.
# DESCRIPTION: <LI> If the extraction radii of the two sources overlap,
# DESCRIPTION: reduce the radius of the second source so that
# DESCRIPTION: they just touch.
# DESCRIPTION: <LI> Repeat these steps for each of the remaining
# DESCRIPTION: sources.
# DESCRIPTION: </OL>
# DESCRIPTION: <P>
# DESCRIPTION: </OL>
# DESCRIPTION: <P>
# DESCRIPTION: After the sources are detected, their positions, radii,
# DESCRIPTION: and total counts are stored in FITS format source catalogs.
# DESCRIPTION: <P>
# DESCRIPTION: Finally, this routine creates <TT>SAOimage</TT>/<TT>xselect</TT>
# DESCRIPTION: style region filters for each extracted source.
# DESCRIPTION: The sources are detected in sky coordinates, but the region
# DESCRIPTION: filters must be in detector coordinates in order for
# DESCRIPTION: spectra to be extracted properly. The conversion from sky to
# DESCRIPTION: detector coordinates is done in the <TT>sky2det</TT>
# DESCRIPTION: utility routine. <TT>sky2det</TT>
# DESCRIPTION: calculates the mean detector coordinate
# DESCRIPTION: position of the events within a given radius of the peak
# DESCRIPTION: in sky cordinates. First a radius of
# DESCRIPTION: 2 pixels is tried. If there are no events within that radius
# DESCRIPTION: (e.g. a dim extended source or one not in
# DESCRIPTION: the field of view of one detector), then
# DESCRIPTION: the conversion is retried using the extraction radius.
# DESCRIPTION: If there are no events within the extraction radius, this
# DESCRIPTION: routine gives an error and does not produce a region filter.
# DESCRIPTION: Note that there are separate region filters for individual
# DESCRIPTION: instruments, since the conversion from sky to detector
# DESCRIPTION: coordinates may not be the same for the two SISs or GISs.
# DESCRIPTION: The routine also creates "source 0" region files which
# DESCRIPTION: <em>exclude</em> all of the detected sources. Note that
# DESCRIPTION: the exclusion radius is never smaller than the standard
# DESCRIPTION: (4 or 6 arc minute) radius.
#
# VERSION: 2.4
#
# HISTORY: 0.0 -> 1.0 10/4/96
# HISTORY: Added code to create a FITS BINTABLE list of sources.
# HISTORY:
# HISTORY: 1.0 -> 1.1 12/4/96
# HISTORY: Now copies all relevant keywords from the totsky image headers
# HISTORY: to the smoothed headers, and adds the SEQNUM keyword.
# HISTORY: This is so that the colorpic tool can read some of this info
# HISTORY: to label the plot it creates.
# HISTORY: Modified to accomodate the new version of colorpic, which
# HISTORY: creates a postscript plot.
# HISTORY: Also stopped deleting the smoothed images.
# HISTORY: Modified the source catalog format by deleting the spectral
# HISTORY: fitting columns and adding the sequence column.
# HISTORY:
# HISTORY: 1.1 -> 1.2 12/19/96
# HISTORY: Now generates region files for each source in this routine.
# HISTORY:
# HISTORY: 1.2 -> 1.3 1/24/97
# HISTORY: Supplied a missing "let" when calculating the default extraction
# HISTORY: radius for gis data. This resulted in gis sources with the
# HISTORY: default radius being dropped. Also deleted a few unused reads from
# HISTORY: the parfile.
# HISTORY: 1/28/97 Supplied code for the background region file in the case
# HISTORY: where there are no foreground sources. Also added comments to the
# HISTORY: top of each region file.
# HISTORY: 2/4/97 Modified the source catalog format.
# HISTORY: 2/7/97 Began adding the same keywords to the sourcecat as to the
# HISTORY: last smoothed image.
# HISTORY:
# HISTORY: 1.3 -> 1.4 3/26/97
# HISTORY: Began clipping smoothed image pixels with low but non-zero
# HISTORY: exposures. This was to eliminate the spurious one-pixel
# HISTORY: sources which can crop up around the edges of GIS images
# HISTORY: of very bright sources.
# HISTORY: Note that this requires STOOLS.3.2 or higher. Also, ascasource
# HISTORY: in STOOLS.3.2 corrects a pixel numbering problem in earlier versions.
# HISTORY: Pixels are now numbered from "1" instead of "0" in the source
# HISTORY: catalog.
# HISTORY:
# HISTORY: 1.4 -> 1.5 3/31/97
# HISTORY: Now includes extended/compact flag in sourcecat. Note this requires
# HISTORY: sift-sources in STOOLS 3.3 or higher.
# HISTORY: Changed grouping when calculating defaultr.
# HISTORY: Now creates region files in detector coordinates.
# HISTORY: Source-free region files with no sources are now a BOX region
# HISTORY: instead of a CIRCLE.
# HISTORY:
# HISTORY: 1.5 -> 1.6 5/14/97
# HISTORY: Now increments the sources counter when there are no counts
# HISTORY: within the extraction radius for sky to det coordinate conversion.
# HISTORY: Now, when either low or high energy images contain no counts,
# HISTORY: issues an exception and makes a black and white colorpic from
# HISTORY: the all-energy image. This can happen when the GIS low energy
# HISTORY: discriminator (keyword LE_DS) is set to something other than zero.
# HISTORY:
# HISTORY: 1.6 -> 1.7 6/3/97
# HISTORY: Accomodated shortened "no photons" message from $UTIL/sky2det.
# HISTORY:
# HISTORY: 1.7 -> 1.8 6/11/97
# HISTORY: Now skips source detection for smoothed images with no counts.
# HISTORY: Also fixed a bug that was hiding the exit code from ascasource.
# HISTORY:
# HISTORY: 1.8 -> 1.9 7/14/97
# HISTORY: Fixed a problem with source-zero region files when there are
# HISTORY: no detected sources. Previously the size of the extraction box
# HISTORY: was only half as large as it was supposed to be due to a
# HISTORY: misinterpretation of the region syntax.
# HISTORY:
# HISTORY: 1.9 -> 1.10 8/6/97
# HISTORY: Now deletes empty region files which can occur when a given source
# HISTORY: is outside the field of view of an instrument. Before these were
# HISTORY: left empty, and so would have bad side effects when they were
# HISTORY: used for extraction.
# HISTORY:
# HISTORY: 1.10 -> 2.0 9/18/97
# HISTORY: Now uses the sky2pix STOOL to calculate ra and dec. This tool takes
# HISTORY: into account corrections for the projection from spherical
# HISTORY: coordinates to a plane. Before this a linear mapping was assumed.
# HISTORY: The errors seem to be on the order of .001 degrees.
# HISTORY:
# HISTORY: 2.0 -> 2.1 9/25/97
# HISTORY: Made it more silent when there are no source-zero region files.
# HISTORY:
# HISTORY: 2.1 -> 2.2 10/9/97
# HISTORY: Now specifies the primary extension of the totsky when extracting
# HISTORY: keywords to be used in smoothed image and source catalog.
# HISTORY: If there is only one screened event file for a given instrument
# HISTORY: resolution, etc., then merge_images will retain the GTI extension
# HISTORY: from the original file. If the primary extension isn't specified
# HISTORY: explicitly, the GTI extensions will be dumped too, and things
# HISTORY: like the TFORM keywords will screw the FITS files.
# HISTORY:
# HISTORY: 2.2 -> 2.3 2/6/98
# HISTORY: Added "1" "repeat count" to the TFORM specifications for the
# HISTORY: source catalog file. Apparently leaving it out is valid FITS,
# HISTORY: but can give trouble to IDL FITS readers.
# HISTORY:
# HISTORY: 2.3 -> 2.4 1999-10-14
# HISTORY: Now delete sources with extraction radius < %minsrcradius
#
###########################################################################
#DEBUG=1
columns=columns.tmp
data=data.tmp
header=header.tmp
$UTIL/milestone "Detecting sources in summed images"
sourcelist=source.tmp
expoclipfrac=$($UTIL/read_parfile $PARFILE expoclipfrac)
minsrcradius=$($UTIL/read_parfile $PARFILE minsrcradius)
for inst in gis sis; do
###################################
# Imaging mode for each instrument
###################################
case "$inst" in
sis) mode=02;;
gis) mode=70;;
*) $UTIL/exception $0 1 "Unknown instrument $inst"
esac
###################
# Debugging output
###################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: inst=$inst
echo ${0##*/}: mode=$mode
fi
#################################
# Loop over all image dimensions
#################################
list=$($UTIL/generate_filename totsky $inst any $mode any all)
list=$(ls $list 2>/dev/null)
for imagefile in $list; do
dimen=$($UTIL/read_fits_keyword ${imagefile}[0] NAXIS1)
###################
# Debugging output
###################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: imagefile=$imagefile
echo ${0##*/}: dimen=$dimen
fi
#########################
# Loop over all energies
#########################
rm -f $sourcelist
touch $sourcelist
for nrg in all hi lo; do
totsky=$($UTIL/generate_filename totsky \
$inst $dimen $mode any $nrg)
totexpo=$($UTIL/generate_filename totexpo $inst $dimen $mode any)
smooth=$($UTIL/generate_filename smooth \
$inst $dimen $mode any $nrg)
########################################################
# Read exposure keyword and set exposure clipping level
########################################################
exposure=$($UTIL/read_fits_keyword $totexpo[0] EXPOSURE)
clip_expo=$($STOOL/equals "$exposure * $expoclipfrac")
###################
# Debugging output
###################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: nrg=$nrg
echo ${0##*/}: totsky=$totsky
echo ${0##*/}: totexpo=$totexpo
echo ${0##*/}: smooth=$smooth
echo ${0##*/}: exposure=$exposure
echo ${0##*/}: clip_expo=$clip_expo
fi
#############################
# Create a smoothed flux map
#############################
$UTIL/log_entry "Smoothing $totsky with $totexpo"
$UTIL/log_entry "Clipping exposures below $clip_expo seconds"
$STOOL/vsmooth $totsky $totexpo $smooth $clip_expo \
>stdout.log 2>stderr.log
$UTIL/stool_test vsmooth $? $0 2 stdout.log stderr.log
#################################################
# Copy image header keywords into smoothed image
#################################################
$UTIL/setup_parfile $FTOOL/fdump.par infile=$totsky[0] \
outfile=out.tmp \
columns=" " \
rows=" " \
prhead=yes \
prdata=no \
showcol=no \
showunit=no \
showrow=no \
showscale=no
rm -f out.tmp
rm -f header.tmp
$FTOOL/fdump >stdout.log 2>stderr.log
$UTIL/ftool_test fdump $? $0 2 stdout.log stderr.log
grep -v SIMPLE out.tmp | \
grep -v NAXIS | \
grep -v BITPIX | \
grep -v COMMENT | \
grep -v HISTORY | \
grep -v CONTENT | \
grep -v FILIN | \
grep -v CREATOR | \
grep -v END >$header
rm -f out.tmp
sequence=$($UTIL/read_parfile $JOBPAR sequence)
echo "SEQNUM = $sequence" >>$header
$UTIL/setup_parfile $FTOOL/fmodhead.par infile=${smooth}[0] \
tmpfil=$header
$FTOOL/fmodhead >stdout.log 2>stderr.log
$UTIL/ftool_test fmodhead $? $0 2 stdout.log stderr.log
###################################
# Detect sources in smoothed image
###################################
counts=$($UTIL/read_fits_keyword $smooth[0] TOTCTS)
if [ $counts -eq 0 ]; then
################################################
# Image is empty, so don't run source detection
################################################
$UTIL/log_entry \
"$smooth has no counts, so skipping source detection"
else
$UTIL/log_entry "Detecting sources in $smooth"
$STOOL/ascasource $smooth >stdout.log 2>stderr.log
code=$?
cat stdout.log >>$sourcelist
$UTIL/stool_test ascasource $code $0 2 stdout.log stderr.log
fi
done
##########################################################
# Default radius based on instrument and image resolution
##########################################################
case "$inst" in
sis)
let defaultr=38*$dimen/320
let compactr=10*$dimen/320
;;
gis)
let defaultr=24*$dimen/256
let compactr=10*$dimen/256
;;
esac
$UTIL/log_entry "Determining extraction radii"
##############################################
# Determine extraction radius for each source
##############################################
rm -f z.tmp
touch z.tmp
nsources=$(wc -l <$sourcelist)
typeset -i source=1
while [ $source -le $nsources ]; do
i0=$( awk ' NR == line {print $1}' line=$source $sourcelist)
j0=$( awk ' NR == line {print $2}' line=$source $sourcelist)
maxr=$(awk ' NR == line {print $4}' line=$source $sourcelist)
minr=$(awk ' NR == line {print $5}' line=$source $sourcelist)
###################
# Determine radius
###################
radius=$defaultr
if [ $minr -gt $defaultr ]; then
radius=$minr
fi
if [ $radius -gt $maxr ]; then
radius=$maxr
fi
##################
# Is it extended?
##################
if [ $minr -gt $compactr ]; then
extended=T
else
extended=F
fi
#################################
# Record the position and radius
#################################
echo "$i0 $j0 $radius $extended" >> z.tmp
let source=$source+1
done
rm -f $sourcelist
#################################################
# Sift through the sources to eliminate overlaps
#################################################
$UTIL/log_entry "Eliminating redundant sources"
$STOOL/sift_sources <z.tmp >stdout.log 2>stderr.log
code=$?
nawk '$3 >= min' min=$minsrcradius stdout.log >$sourcelist
$UTIL/stool_test sift_sources $code $0 2 stdout.log stderr.log
$UTIL/log_entry "Sources with radius >= $minsrcradius"
$UTIL/file_to_log $sourcelist
rm -f z.tmp
######################
# Make a nice picture
######################
losky=$($UTIL/generate_filename smooth $inst $dimen $mode any lo)
allsky=$($UTIL/generate_filename smooth $inst $dimen $mode any all)
hisky=$($UTIL/generate_filename smooth $inst $dimen $mode any hi)
sourcepic=$($UTIL/generate_filename sourcepic $inst $dimen $mode any)
################################################
# Check if there are counts in all three images
################################################
lo_counts=$($UTIL/read_fits_keyword $losky[0] TOTCTS)
hi_counts=$($UTIL/read_fits_keyword $hisky[0] TOTCTS)
rm -f $sourcepic
if [ $lo_counts -eq 0 -o $hi_counts -eq 0 ]; then
########################################################
# If the low or high energy image has no counts
# just make a black and white image from the all energy
# image
########################################################
$UTIL/exception $0 1 "Either $losky or $hisky has no counts"
awk '{print $1" "$2" "$3}' $sourcelist | \
$STOOL/colorpic $sourcepic $allsky $allsky $allsky \
>stdout.log 2>stderr.log
$UTIL/stool_test colorpic $? $0 2 stdout.log stderr.log
else
#####################
# Make a color image
#####################
awk '{print $1" "$2" "$3}' $sourcelist | \
$STOOL/colorpic $sourcepic $losky $allsky $hisky \
>stdout.log 2>stderr.log
$UTIL/stool_test colorpic $? $0 2 stdout.log stderr.log
fi
#########################################
# Make a FITS format list of all sources
#########################################
sourcecat=$($UTIL/generate_filename sourcecat $inst $dimen $mode any)
$UTIL/log_entry "Creating fits catalog of sources: $sourcecat"
rm -f $sourcecat
rm -f $data
rm -f $columns
#######################################
# Column info for source catalog files
#######################################
rm -f $columns
echo "sequence 1J" >> $columns
echo "index 1I" >> $columns
echo "xpixel 1I pixels" >> $columns
echo "ypixel 1I pixels" >> $columns
echo "rpixel 1I pixels" >> $columns
echo "RA 1E degrees" >> $columns
echo "Dec 1E degrees" >> $columns
echo "rdeg 1E degrees" >> $columns
echo "hardcnts 1J counts" >> $columns
echo "softcnts 1J counts" >> $columns
echo "extended 1L " >> $columns
########################
# Read coordinate scale
########################
totsky=$($UTIL/generate_filename totsky $inst $dimen $mode any all)
totsky_lo=$( $UTIL/generate_filename totsky $inst $dimen $mode any lo)
totsky_hi=$( $UTIL/generate_filename totsky $inst $dimen $mode any hi)
cdelt1=$($UTIL/read_fits_keyword ${totsky}[0] CDELT1)
################################
# List sources to an ascii file
################################
nsources=$(wc -l <$sourcelist)
touch $data
typeset -i source=1
while [ $source -le $nsources ]; do
x=$( awk ' NR == line {print $1}' line=$source $sourcelist)
y=$( awk ' NR == line {print $2}' line=$source $sourcelist)
r=$( awk ' NR == line {print $3}' line=$source $sourcelist)
extended=$(awk ' NR == line {print $4}' line=$source $sourcelist)
echo "$x $y" | $STOOL/pix2sky $totsky | read ra dec
rdeg=$($STOOL/equals "abs(${r}*${cdelt1})")
softcnts=$($STOOL/counts_in_circle $totsky_lo $x $y $r)
hardcnts=$($STOOL/counts_in_circle $totsky_hi $x $y $r)
echo "$sequence $source $x $y $r $ra $dec $rdeg" \
"$softcnts $hardcnts $extended" >> $data
let source=$source+1
done
##################################################
# Clear some unwanted things from the header file
##################################################
rm -f z.tmp
grep -v RPIX $header | \
grep -v RVAL | \
grep -v DELT | \
grep -v DUNIT | \
grep -v CROTA | \
grep -v CTYPE | \
grep -v CHANMIN | \
grep -v CHANMAX | \
grep -v CHANTYPE | \
grep -v TOTCTS | \
grep -v USER | \
grep -v HDUCLAS | \
grep -v EXPOSURE >z.tmp
mv z.tmp $header
#################################
# Create the source catalog file
#################################
$UTIL/setup_parfile $FTOOL/fcreate.par cdfile=$columns \
datafile=$data \
headfile=$header \
outfile=$sourcecat \
tbltype=binary
$FTOOL/fcreate >stdout.log 2>stderr.log
$UTIL/ftool_test fcreate $? $0 2 stdout.log stderr.log
rm -f $data
rm -f $columns
rm -f $sourcelist
rm -f header.tmp
done #loop over image dimension
done #loop over instrument
####################################################################
# Generate region files for all instruments, dimensions and sources
####################################################################
$UTIL/log_entry "Generating region files"
regiondir=$($UTIL/read_parfile $PARFILE regiondir)
########################
# Loop over instruments
########################
for inst in s0 s1 g2 g3; do
case "$inst" in
s?)
sourceinst=sis
display=SIS${inst#s}
;;
g?)
sourceinst=gis
display=GIS${inst#g}
;;
esac
############################
# Loop over source catalogs
############################
list=$($UTIL/generate_filename sourcecat $sourceinst any any any)
list=$(ls $list 2>/dev/null)
for sourcecat in $list; do
dimen=$($UTIL/parse_filename index $sourcecat)
mode=$( $UTIL/parse_filename mode $sourcecat)
nsources=$($UTIL/read_fits_keyword ${sourcecat}[1] NAXIS2)
#######################
# Image binning factor
#######################
case "$inst" in
s?)
bin=4
let defaultr=38*$dimen/320
;;
g?)
bin=1
let defaultr=24*$dimen/256
;;
esac
###################
# Debugging output
###################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: sourcecat=$sourcecat
echo ${0##*/}: inst=$inst
echo ${0##*/}: dimen=$dimen
echo ${0##*/}: mode=$mode
echo ${0##*/}: bin=$bin
echo ${0##*/}: defaultr=$defaultr
echo ${0##*/}: nsources=$nsources
fi
################################
# Header for source-free region
################################
freeregion=$($UTIL/generate_filename sourceregion $inst $dimen \
$mode none 0)
rm -f $freeregion
echo "# This is the $display, source 0 region file" >$freeregion
echo "# at $dimen pixel resolution." >>$freeregion
echo "# It gives the region NOT occupied by any of the" >>$freeregion
echo "# detected sources. " >>$freeregion
echo "# This region file is in DETECTOR coordinates" >>$freeregion
case "$inst" in
s?)
echo "# Data extracted from this region may be" >>$freeregion
echo "# used cautiously for background subtraction " >>$freeregion
echo "# of SIS data for which vignetting is small." >>$freeregion
echo "#" >>$freeregion
;;
g?)
echo "# Vignetting makes it unwise to use data" >>$freeregion
echo "# extracted from this region for background" >>$freeregion
echo "# subtraction." >>$freeregion
;;
esac
if [ $nsources -eq 0 ]; then
#######################################################
# No foreground sources, so background is entire field
#######################################################
echo "# There are no foreground sources, so this" >>$freeregion
echo "# region occupies the entire field of view" >>$freeregion
echo "#" >>$freeregion
center=$($STOOL/equals "$dimen/2")
echo " BOX($center,$center,$dimen,$dimen)" >>$freeregion
fi
####################
# Loop over sources
####################
source=1
while [ $source -le $nsources ]; do
#####################
# Source region file
#####################
region=$($UTIL/generate_filename sourceregion $inst $dimen \
$mode none $source)
rm -f $region
echo "# This is the $display, source $source region file" >$region
echo "# at $dimen pixel resolution" >>$region
echo "# This region file is in DETECTOR coordinates" >>$region
echo "#" >>$region
x=$($UTIL/read_bintable $sourcecat xpixel $source)
y=$($UTIL/read_bintable $sourcecat ypixel $source)
r=$($UTIL/read_bintable $sourcecat rpixel $source)
##################################################
# Convert sky coordinates to detector coordinates
##################################################
let x=$x*$bin
let y=$y*$bin
#############################
# First try 2.0 pixel radius
#############################
result=$($UTIL/sky2det $x.0 $y.0 $inst $mode 2.0)
if [ "$result" = "No photons" ]; then
############################
# Try the extraction radius
############################
$UTIL/log_entry "No photons in 2.0 pixel radius"
result=$($UTIL/sky2det $x.0 $y.0 $inst $mode $r.0)
if [ "$result" = "No photons" ]; then
##############################################
# No photons in full radius, somethings wrong
##############################################
$UTIL/log_entry \
"No photons for inst $inst, dimen $dimen, source $source"
rm -f $region
let source=$source+1
continue
fi
fi
######################################################
# Extract detector coordinates and record region file
######################################################
detx=$(echo "$result" | awk '{print $1}')
dety=$(echo "$result" | awk '{print $2}')
detx=$($STOOL/equals "int(($detx-1)/$bin + 1.5)")
dety=$($STOOL/equals "int(($dety-1)/$bin + 1.5)")
echo " CIRCLE($detx,$dety,$r)" >>$region
##########################
# Source-free region file
##########################
if [ $r -lt $defaultr ]; then
r=$defaultr
fi
echo "-CIRCLE($detx,$dety,$r)" >>$freeregion
let source=$source+1
done # loop over sources
done # loop over source catalogs
done # loop over instruments
########################################################################
# Remove any empty source 0 region files
# these files can be empty is all detected sources are out of the field
# of view for a given detector.
########################################################################
list=$($UTIL/generate_filename sourceregion any any any none 0)
list=$(ls $list 2>>/dev/null)
for region in $list; do
count=$(grep -v "^#" $region | wc -l)
if [ $count -eq 0 ]; then
$UTIL/log_entry "Removing empty region file $region"
rm -f $region
fi
done
rm -f include.tmp
rm -f exclude.tmp
exit 0