images (version 2.8)
You can also look at:
#! /usr/bin/ksh
###########################################################################
#
# SYNTAX: images
#
# BRIEFLY: Generate images and exposure maps
#
# DESCRIPTION: This routine creates exposure maps and sky images for all
# DESCRIPTION: of the filtered event files.
# DESCRIPTION: Images are generated in three energy bands: lo < 2keV,
# DESCRIPTION: hi > 2keV, and all.
# DESCRIPTION: The actual PI channels used for the energy cut are
# DESCRIPTION: %gis_picut_128 , %gis_picut_256 , or %gis_picut_1024
# DESCRIPTION: for GIS spectra with 128, 256, or 1024 PI channels, respectively.
# DESCRIPTION: For SIS spectra the PI channel is %sis_picut .
# DESCRIPTION: The exposure maps take into account the region filtering
# DESCRIPTION: of GIS files to remove the noisy outer regions, but they
# DESCRIPTION: do not correct for GIS deadtime. Note that periods of saturated
# DESCRIPTION: telemetry are removed from the SIS filtered event files,
# DESCRIPTION: so no correction is needed for this effect.
# DESCRIPTION: The exposure maps do not account for vignetting in either
# DESCRIPTION: instrument.
# DESCRIPTION: <P>
# DESCRIPTION: The individual images for each filtered event file are summed
# DESCRIPTION: so that there are two files of each type - one combining
# DESCRIPTION: GIS2 and GIS3 and one combining SIS0 and SIS1.
# DESCRIPTION: <P>
# DESCRIPTION: Finally <TT>ximage</TT> is used to produce an
# DESCRIPTION: exposure-corrected color image from the "all" energy
# DESCRIPTION: summed images.
# DESCRIPTION: <TT>ximage</TT> produces a PPM format image which is converted to
# DESCRIPTION: GIF format for easier viewing.
#
# VERSION: 2.8
#
# HISTORY: 0.0 -> 1.0 9/19/96
# HISTORY: Now deletes individual exposure maps.
# HISTORY:
# HISTORY: 1.0 -> 1.1 9/27/96
# HISTORY: Now makes all energy counts map with xselect along with the hi and
# HISTORY: low energy images. This is because ascaexpo does not put in all
# HISTORY: the keywords that the extractor does. Also ascaexpo images have
# HISTORY: floating point pixels instead of integer pixels (in case this
# HISTORY: matters).
# HISTORY:
# HISTORY: 1.1 -> 1.2 10/17/96
# HISTORY: Don't generate individual exposure maps for BRIGHT2 event files
# HISTORY: any more. This is because they are not used in the summed maps
# HISTORY: and probably won't be kept anyway. Also they are redundant with
# HISTORY: the corresponding BRIGHT mode files. Before this the BRIGHT2
# HISTORY: exposure maps were not being deleted and weren't included in
# HISTORY: the distribution or trend products.
# HISTORY:
# HISTORY: 1.2 -> 1.3 12/04/96
# HISTORY: Corrected TOTCNTS to TOTCTS for the total counts keyword
# HISTORY: in the image files. This bug was creating a new keyword
# HISTORY: instead of updating the value of the existing one.
# HISTORY: Also started updating instrument keyword to
# HISTORY: GIS2,GIS3 or SIS0,SIS1 as appropriate.
# HISTORY:
# HISTORY: 1.3 -> 1.4 2/4/97
# HISTORY: GIF format image from ximage is now exposure-corrected.
# HISTORY: 2/5/97 Changed from floatadd to equals STOOLS.
# HISTORY: Added some quotes to add_fits_keyword calls.
# HISTORY:
# HISTORY: 1.4 -> 1.5 3/18/97
# HISTORY: Now ignores event files which have ontime less than a tolerance.
# HISTORY: This is to avoid large regions of very short exposure time which
# HISTORY: can screw up exposure-corrected images and source detection.
# HISTORY: Note: This is only done for SIS event files.
# HISTORY:
# HISTORY: 1.5 -> 1.6 3/26/97
# HISTORY: Undid the last change. The same thing has been achieved in
# HISTORY: the sources routine by clipping individual low-exposure
# HISTORY: pixels when smoothing and exposure-correcting.
# HISTORY:
# HISTORY: 1.6 -> 1.7 5/5/97
# HISTORY: Made compatible with GIS region files that depend on the
# HISTORY: spread discriminator.
# HISTORY:
# HISTORY: 1.7 -> 1.8 5/15/97
# HISTORY: Corrected read/exposure error in ximage script.
# HISTORY: Should instead be read/exp_map.
# HISTORY:
# HISTORY: 1.8 -> 1.9 7/11/97
# HISTORY: Added a check for exposure maps <60 seconds. Ximage can't handle
# HISTORY: such sort exposures, so no GIF is produced for them.
# HISTORY:
# HISTORY: 1.9 -> 2.0 8/5/97
# HISTORY: Now uses a region filter to exclude the (1,1) pixel where ascalin
# HISTORY: dumps photons outside the image coordinates. This was done because
# HISTORY: it could screw up ximage if the corner pixel had the highest
# HISTORY: value in the image. Note that no info is lost by doing this,
# HISTORY: since the event files still retain these events.
# HISTORY:
# HISTORY: 2.0 -> 2.1 12/30/97
# HISTORY: Now explicitly sets the equinox in ximage to 2000.
# HISTORY: Previously the equinox was not specified and it was likely
# HISTORY: defaulting to 1950 (although it could have been anything).
# HISTORY:
# HISTORY: 2.1 -> 2.2 1/20/98
# HISTORY: Now deletes the "all" energy individual images. Before these were
# HISTORY: saved for rev1 compatibility, although they were assigned
# HISTORY: a trend class of "trash" and thus deleted anyway.
# HISTORY: This saves some work in post-processing.
# HISTORY:
# HISTORY: 2.2 -> 2.3 4/13/98
# HISTORY: Now skip making exposure maps and images for FAST mode files.
# HISTORY: Previously a bug elsewhere kept FASt mode filtered event files
# HISTORY: from being produced. Also explicitly skip MOC mode files even though
# HISTORY: the current software does not allow filtered MPC mode "event"
# HISTORY: files.
# HISTORY:
# HISTORY: 2.3 -> 2.4 1998-09-17
# HISTORY: Modified to use xanadu initialization. Now interlace GIFs
# HISTORY: using ppmtogif -interlace instead of using giftool
# HISTORY:
# HISTORY: 2.4 -> 2.5 1999-03-24
# HISTORY: Took out $UTIL/generate_filename call for rate file since
# HISTORY: the rate file is not used by ascaexpo
# HISTORY:
# HISTORY: 2.5 -> 2.6 1999-04-07
# HISTORY: Changed the ximage script to override the default title
# HISTORY: in the GIF plots. This is because ximage by default does
# HISTORY: not correctly handle INSTRUME="inst1,inst2".
# HISTORY:
# HISTORY: 2.6 -> 2.7 1999-11-30
# HISTORY: Fixed a bug where $inst was given as the argument to
# HISTORY: $UTIL/generate_filename instead of $inst_type. This caused
# HISTORY: problems when there were GIS images with more than one dimension.
# HISTORY:
# HISTORY: 2.7 -> 2.8 2000-02-02
# HISTORY: Changed "read/exp_map" to "read/expo" in ximage command file
# HISTORY: to comply with the syntax used in XANADU 5.0. Also dropped the
# HISTORY: explicit title. The title had been added in images 2.6. However
# HISTORY: The problem which prompted that change has been fixed and
# HISTORY: using an explicit lower title activates another bug which makes
# HISTORY: ximage not read the exposure correctly.
#
############################################################################
#DEBUG=1
$UTIL/milestone "Generating images and exposure maps"
####################################################
# Make a region file which removes the corner point
####################################################
rm -f corner_point.reg
echo "-POINT(1,1)" > corner_point.reg
########################################
# Loop through all filtered event files
########################################
list=$($UTIL/generate_filename event any any any any)
list=$(ls $list 2>/dev/null)
for event in ${list} ; do
#############
# File names
#############
inst=$( $UTIL/parse_filename inst $event)
index=$( $UTIL/parse_filename index $event)
mode=$( $UTIL/parse_filename mode $event)
bitrate=$($UTIL/parse_filename bitrate $event)
#################################################################
# Skip BRIGHT2 mode files since their time is included in BRIGHT
# mode (converted) files
#################################################################
if [ "$mode" = "12" ]; then
continue
fi
#######################################
# Skip non imaging FAST and MPC modes
#######################################
if [ "$mode" = "03" -o "$mode" = "71" ]; then
continue
fi
teldef=$($UTIL/generate_filename teldef $inst $index $mode $bitrate)
expo=$( $UTIL/generate_filename expo $inst $index $mode $bitrate)
rm -f $expo
###################
# Debugging output
###################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: inst=$inst
echo ${0##*/}: index=$index
echo ${0##*/}: mode=$mode
echo ${0##*/}: bitrate=$bitrate
echo ${0##*/}: event=$event
echo ${0##*/}: teldef=$teldef
echo ${0##*/}: rate=$rate
echo ${0##*/}: expo=$expo
echo ${0##*/}: sky=$sky
fi
##############################################################
##############################################################
##
## Exposure maps
##
##############################################################
##############################################################
$UTIL/log_entry "Generating exposure map $expo"
##########################################
# Extraction region map for exposure maps
##########################################
case "$inst" in
g[23])
rawxbins=$($UTIL/read_fits_keyword ${event}[1] RAWXBINS)
spread=$( $UTIL/read_fits_keyword ${event}[0] S_DSCR)
extractregion=$($UTIL/read_parfile $PARFILE \
${inst}region${rawxbins}$spread)
dir=$($UTIL/read_parfile $PARFILE regiondir)
$UTIL/fetch_file $dir $extractregion
regionmap=${extractregion}_regionmap.tmp
if [ ! -a "$regionmap" ]; then
####################################
# Haven't already generated the map
####################################
$STOOL/regionmap $regionmap $extractregion $rawxbins \
>stdout.log 2>stderr.log
$UTIL/stool_test regionmap $? $0 2 stdout.log stderr.log
fi
imath="mul"
;;
s[01])
regionmap="none"
imath="none"
;;
*) $UTIL/exception $0 1 "Unknown instrument $inst"
esac
###################
# Debugging output
###################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: imath=$imath
echo ${0##*/}: regionmap=$regionmap
fi
####################
# Fetch teldef file
####################
caldir=$($UTIL/read_parfile ./$PARFILE caldir)
$UTIL/fetch_file $caldir $teldef
###############
# Run ascaexpo
###############
$UTIL/log_entry "Running ascaexpo for $event"
$UTIL/setup_parfile $FTOOL/ascaexpo.par \
datafile=$event \
calfile=$teldef \
attitude=default \
defattpath=./ \
ratefile=none \
satfcol="SAT" \
instfile="$regionmap" \
imath="$imath" \
expofile=$expo \
imapfile=none \
skyfile=none \
attstep=5 \
rebin=-1 \
binoff=0.0
$FTOOL/ascaexpo >stdout.log 2>stderr.log
$UTIL/ftool_test ascaexpo $? $0 2 stdout.log stderr.log
##############################################################
##############################################################
##
## Low, High, and All Energy Images
##
##############################################################
##############################################################
$UTIL/log_entry "Generating low and high energy images for $event"
command=command.tmp
###########################################
# Determine PI bin for high/low energy cut
###########################################
case "$inst" in
g*)
pibins=$($UTIL/read_fits_keyword ${event}[0] PHA_BINS)
picut=$($UTIL/read_parfile $PARFILE gis_picut_$pibins);;
s*)
pibins=2048
picut=$($UTIL/read_parfile $PARFILE sis_picut);;
esac
losky=$($UTIL/generate_filename sky $inst ${index} $mode $bitrate lo )
hisky=$($UTIL/generate_filename sky $inst ${index} $mode $bitrate hi )
sky=$($UTIL/generate_filename sky $inst ${index} $mode $bitrate all)
rm -f $losky
rm -f $hisky
rm -f $sky
#############################
# Setup xselect command file
#############################
rm -f $command
echo "${event%.*}" >>$command
echo " " >>$command
echo "%ED%off" >>$command
echo "set mission ASCA" >>$command
echo "set datadir ." >>$command
echo "read events $event" >>$command
echo "set image sky" >>$command
echo "filter region corner_point.reg" >>$command
echo "extract image" >>$command
echo "save image outfile=${sky}" >>$command
echo "filter pha_cut 0 $picut" >>$command
echo "extract image" >>$command
echo "save image outfile=${losky}" >>$command
echo "clear pha_cut" >>$command
echo "filter pha_cut $picut $pibins" >>$command
echo "extract image" >>$command
echo "save image outfile=${hisky}" >>$command
echo "exit save_session=no" >>$command
##############
# Run xselect
##############
rm -f ${losky}
rm -f ${hisky}
$FTOOL/xselect @${command} </dev/null > stdout.log 2> stderr.log
$UTIL/xselect_test $? $0 stdout.log stderr.log
rm -f $command
done
rm -f *_regionmap.tmp
rm -f corner_point.reg
########################################
########################################
##
## Combine Images
##
########################################
########################################
datalist=datalist.tmp
listfile=list.tmp
ppm=pgplot.ppm
export DISPLAY=$($UTIL/read_parfile $PARFILE display)
ppm2gif=$($UTIL/read_parfile $PARFILE ppm2gif)
datalist=datalist.tmp
listfile=listfile.tmp
for inst_type in sis gis; do
$UTIL/log_entry "Summing $inst_type images"
#########################################################
# Select BRIGHT mode for the SIS and PH mode for the GIS
#########################################################
case "$inst_type" in
sis) mode=02;;
gis) mode=70;;
esac
###########################################################
# List the files and dimensions like this: NAXIS1 filename
###########################################################
list=$($UTIL/generate_filename expo ${inst_type%is}? any $mode any)
list=$(ls $list 2>/dev/null)
rm -f dimenlist.tmp
touch dimenlist.tmp
for file in $list ; do
naxis=$($UTIL/read_fits_keyword ${file}[0] NAXIS1)
inst=$( $UTIL/parse_filename inst $file)
index=$( $UTIL/parse_filename index $file)
bitrate=$($UTIL/parse_filename bitrate $file)
echo "$naxis ${inst}|${index}|${bitrate}" >>dimenlist.tmp
done
#####################################
# Loop through each unique dimension
#####################################
for dimen in $(awk '{print $1}'<dimenlist.tmp | sort -n | uniq); do
rm -f $datalist
touch $datalist
awk '$1 == word {print $2}' word=$dimen dimenlist.tmp \
>$datalist 2>/dev/null
###################
# Debugging output
###################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: dimen=$dimen
fi
############################################
# Loop over the types of files to be summed
############################################
for type in expo all lo hi; do
case "$type" in
expo)
summed=$($UTIL/generate_filename \
totexpo $inst_type $dimen $mode any);;
*)
summed=$($UTIL/generate_filename \
totsky $inst_type $dimen $mode any $type);;
esac
###################
# Debugging output
###################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: type=$type
echo ${0##*/}: summed=$summed
fi
###########################
# Assemble a list of files
# and sum the exposure
###########################
totexposure=0.0
totcounts=0
rm -f $listfile
for thing in $(cat $datalist); do
#####################
# Determine filename
#####################
inst=$( echo $thing | awk -F'|' '{print $1 }')
index=$( echo $thing | awk -F'|' '{print $2 }')
bitrate=$(echo $thing | awk -F'|' '{print $3 }')
case "$type" in
expo)
file=$($UTIL/generate_filename \
expo $inst $index $mode $bitrate);;
*)
file=$($UTIL/generate_filename \
sky $inst $index $mode $bitrate $type);;
esac
###################
# Debugging output
###################
if [ -n "$DEBUG" ]; then
echo ${0##*/}: file=$file
fi
########################
# Echo filename to list
########################
if [ -a "$file" ]; then
echo "$file" >>$listfile
exposure=$($UTIL/read_fits_keyword ${file}[0] EXPOSURE)
totexposure=$($STOOL/equals $exposure +$totexposure)
if [ "$type" != "expo" ]; then
counts=$($UTIL/read_fits_keyword ${file}[0] TOTCTS)
let totcounts=$totcounts+$counts
fi
else
$UTIL/exception $0 2 "Missing $file for image sum"
fi
done # loop over files
#################
# Add the images
#################
rm -f $summed
$UTIL/merge_images $summed $listfile
#################################################
# Modify INSTRUME, EXPOSURE, and TOTCTS keywords
#################################################
case $inst_type in
s*)
$UTIL/add_fits_keyword ${summed}[0] INSTRUME "SIS0,SIS1" ;;
g*)
$UTIL/add_fits_keyword ${summed}[0] INSTRUME "GIS2,GIS3" ;;
esac
$UTIL/add_fits_keyword ${summed}[0] EXPOSURE "$totexposure"
if [ "$type" != "expo" ]; then
$UTIL/add_fits_keyword ${summed}[0] TOTCTS "$totcounts"
fi
#############################################
# Delete individual images and exposure maps
#############################################
for file in $(cat $listfile); do
rm -f $file
done
rm -f $listfile
done #loop over types
rm -f $datalist
############
# GIF image
############
totexpo=$($UTIL/generate_filename totexpo $inst_type $dimen \
$mode any)
totsky=$( $UTIL/generate_filename totsky $inst_type $dimen \
$mode any all)
index=${dimen#dimen}
gif=$($UTIL/generate_filename gif $inst_type $dimen $mode any)
rm -f $gif
####################################################################
# Check if exposure is large enough.
# Exposures less than 60 seconds aren't handled correctly by ximage
####################################################################
exposure=$($UTIL/read_fits_keyword $totexpo[0] EXPOSURE)
if [ $($STOOL/floatcmp $exposure 60.0) = "lt" ]; then
$UTIL/log_entry "Skipping $gif since EXPOSURE = $exposure"
continue
fi
#####################################################
# extract keyword values to use in image title
#####################################################
object=$( $UTIL/read_fits_keyword $totsky+0 OBJECT )
mission=$($UTIL/read_fits_keyword $totsky+0 TELESCOP)
inst=$( $UTIL/read_fits_keyword $totsky+0 INSTRUME)
date=$( $UTIL/read_fits_keyword $totsky+0 DATE-OBS)
date=$($UTIL/iso_date_to_words $date)
exposure=$($STOOL/equals "int($exposure*10)/10")
title1="$object"
title2="$mission $inst $date Exposure: $exposure s"
$UTIL/log_entry "Running XIMAGE to create $gif"
rm -f $command
echo "read/fits $totsky" >>$command
echo "read/expo $totexpo" >>$command
echo "smooth" >>$command
echo "cpd /ppm" >>$command
echo "cey 2000" >>$command # equinox 2000
# echo "title \"$title1\" " >>$command
# echo "title/lower \"$title2\"">>$command
echo "disp/correct/noframe" >>$command
echo "grid" >>$command
echo "scale" >>$command
echo "exit" >>$command
$XANBIN/bin/ximage @$command </dev/null >stdout.log 2>stderr.log
$UTIL/stool_test ximage $? $0 2 stdout.log stderr.log
########################
# Convert to GIF format
########################
$ppm2gif -interlace < $ppm > $gif 2>> stderr.log
$UTIL/exit_test $? $0 2 "While converting a PPM file to GIF format"
rm -f $ppm
rm -f $command
done # loop over dimen
rm -f dimenlist.tmp
done # loop over inst
exit 0