#! /bin/tcsh -f

# t_cent 4.0

# Finds the centroids of the SIS and GIS images and writes
# a region file.
# Uses code from do_cent (v3.2) script by K.Nandra/D.Turcan

# Adaptation by Paul O'Neill

# Calls the scripts: t_zipc t_detimg do_inst_sky_cent.xco do_gis2_sky_bin2.xco
#                    do_gis3_sky_bin2.xco t_veron do_cent_extr_spec.xco
#                    do_gis2_sky_init.xco do_cent_extr.xco

#### NOTES REGARDING THE PERMANENT DIRECTORY FILES ####

# The centroiding will ALWAYS be performed. However, during the centroiding process,
# the $PERM directory will be checked for any manual files as follows:
#
#             ${inst}_centcus.txt   -   Sky pixel coordinates to use as nominal centroid position,
#                                       instead of the X-Y coordinates that correspond to the RA and DEC
#                                       in the attitude file.
#
#             gis2_srccus.reg      -   Contains the custom region to use when extracting a GIS2 source
#             gis2_bgdcus.reg          and background spectrum. This spectrum is used to determine if the 
#                                      object was detected. In sky pixels.
#
#             ${obsid}_${inst}_manual_${reg}.reg  -  Custom regions to use for all subsequent source and background
#                                                    extractions. In detector coordinates.
#
# These manual files are NEVER over-written or deleted by the Tartarus scripts 


# The centroids in detector coordinates are written to the files detcent_${inst}.temp. These files are
# read by ascaarf (t_spec) to determine the position of the object. Therefore, the centroids given in
# the custom region files (_manual_) must be the same as in detcent_${inst}.temp. This is expected to ALWAYS be
# the case because custom regions are only ever used to remove other sources in the field. The centroiding process below
# is the best way to determine the source centroid. t_cent checks the custom region files to make sure that
# they each contain a line having the coordinates as found by the centroiding. t_cent exits with a fatal error if no
# such line is found.

# The script t_veron creates a smoothed GIS sky image with an overlay of the locations of
# any Veron objects in the field. The script t_Detimg creates detector images with or
# without an overlay of the extraction regions.

# Date in final form: 17 December 2004


########################################################################################
# Welcome 
########################################################################################

echo ""
echo "gnt_infrm (t_cent) ++++++++++++++++++++++++++++++++++++++++++++++++++++++"
echo "gnt_infrm (t_cent)  t_cent 4.0  STARTED "`date`
echo "gnt_infrm (t_cent) ++++++++++++++++++++++++++++++++++++++++++++++++++++++"
echo ""


########################################################################################
# Get parameters
########################################################################################

set obsid = $1           # observation sequence number
set SEQ = $2             # root directory for current sequence
set PERM = $3            # permanent directory for this sequence
set SCRIPTS = $4         # directory containing tartarus scripts        
set XCO = $5             # directory with xselect/Xspec scripts
set REGS = $6            # directory with region files
set VERON = $7           # path and name of text file of VERON2001 catalog


########################################################################################
# Set directory variables
########################################################################################

# NOTE: aux, screened, unscreened, and work directories MUST be sub-directories
# of the sequence directory (SEQ). This structure is hard wired.

set AUX = $SEQ/aux            # contains house keeping data
set SCR = $SEQ/screened       # contains screened data
set UNSCR = $SEQ/unscreened   # contains unscreened data
set WORK = $SEQ/work          # a work area


########################################################################################
# Set XIMAGE colour table
########################################################################################

# the images are created using this colour table

set coltab = "bb"


########################################################################################
# Check and prepare directories
########################################################################################


#### check permanent directory ####

echo "gnt_infrm (t_cent) checking permanent directory"

if ( -e $PERM ) then    # if a PERM directory exists (it should) then unzip files
    $SCRIPTS/t_zipc $PERM
else
    mkdir $PERM     
endif


#### check work directory ####

echo "gnt_infrm (t_cent) checking work directory"

if ( -e $WORK ) then    # if a WORK directory exists (it should) then unzip files
    $SCRIPTS/t_zipc $WORK
else
    echo "gnt_fatal (t_cent) work directory $WORK does not exist"
    exit 1
endif


########################################################################################
# Make sure event and attitude files exist
########################################################################################


#### look for event files #####

foreach inst ( sis0 sis1 gis2 gis3 )
    set file = ${obsid}_$inst.evt
    if  ( ! -e $WORK/$file ) then
	echo "gnt_fatal (t_cent) event file $file does not exist"
	exit 1
    endif
end


#### look for attitude file ####

set attfiles = `ls -1 $AUX/fa*`

if ( $#attfiles == "0" ) then
    echo "gnt_fatal (t_cent) attitude file not found    # this file is required
    exit 1
endif

set attfile = $attfiles[1]      # use first fa* file in $AUX as the attitude file


########################################################################################
# Get some info from the attitude file 
########################################################################################


#### find target coordinates ####

echo "gnt_infrm (t_cent) getting RA and DEC from attitude file"

set ra = `fkeyprint ${attfile}+1 RA | grep "RA      =" | awk '{print $3}'`
set dec = `fkeyprint ${attfile}+1 DEC | grep "DEC     =" | awk '{print $3}'`


#### Get object name ####

if ( ! -e $WORK/object.temp ) then
    echo "gnt_infrm (t_cent) Getting object name from attitude file"
    set object = `fkeyprint ${attfile}+1 OBJECT | grep "OBJECT  =" | sed s/\\[/"\("/g | sed s/\\]/"\)"/g`
    set object = `echo "$object[3] $object[4]" | sed s/"'"//g | sed s/"\/"//g | sed s/" "//g`
    echo "$object" >! $WORK/object.temp
endif

set object = `head -1 $WORK/object.temp`


#### Get otime from attitude file ####

echo "gnt_infrm (t_cent) Getting OTIME from attitude file"
set otime = `fkeyprint ${attfile}+1 MTIME0 | grep "MTIME0  =" | awk '{print int($3)}'`
echo $otime >! $WORK/otime.temp


#### show some info ####

echo "gnt_infrm (t_cent) Object name                       : $object"
echo "gnt_infrm (t_cent) Object RA                         : $ra"
echo "gnt_infrm (t_cent) Object DEC                        : $dec"
echo "gnt_infrm (t_cent) Observation start in mission time : $otime"


########################################################################################
# Delete unwanted files
########################################################################################

cd $WORK

set files = `ls -1 $WORK/*.tmp`
if ( $#files != "0" ) rm -f $WORK/*.tmp    # if there are .tmp files to delete

set ra_dec = ""
set cent = ""


########################################################################################
# Create some detector images
########################################################################################

# runs the script t_detimg
# detector images .img were produced by t_scrn. Now produce some gif and jpeg
# without the extract regions (which don't yet exist). t_detimg
# is run again at the end of t_cent, if the objects is detected, to
# produce images showing the regions.

$SCRIPTS/t_detimg $obsid $SEQ


########################################################################################
# Extract sky images for each instrument
########################################################################################

# These sky images will have the nominal object position shown with a circle. The images
# may be used for a visual inspection of the sources in the field-of-view of each instrument.
# The images will be used to show the region used for centroiding.


#### start loop through each instrument ####

foreach inst ( sis0 sis1 gis2 gis3 )

    echo "gnt_infrm (t_cent) extracting sky image for $inst"


#### find nominal centroid for the inst - this will be shown on the image ####

    set coord = `sky2xy $WORK/${obsid}_${inst}.evt+1 xsky=$ra ysky=$dec | grep pixel | sed s/","//g`
   
    if ( $inst == sis0 || $inst == sis1 ) then
	set rad = 33
	set x = `echo $coord[4] | awk '{print int($1/4)}'`   # nominal RA in rebinned sky pixels
	set y = `echo $coord[5] | awk '{print int($1/4)}'`   # nominal DEC in rebinned sky pixels
    else    # GIS
	set rad = 14
	set x = `echo $coord[4] | awk '{print int($1)}'`     # nominal RA in unbinned sky pixels 
	set y = `echo $coord[5] | awk '{print int($1)}'`     # nominal DEC in unbinned sky pixels
    endif


#### show the custom centroid circle if being used to find centroid ####

    if ( -e $PERM/${inst}_centcus.txt ) then
	echo "gnt_infrm (t_cent) **** SEQUENCE $obsid - using custom centroiding circle file $PERM/${inst}_centcus.txt ****"
	echo "gnt_warng (t_cent) **** SEQUENCE $obsid - using custom centroiding circle file $PERM/${inst}_centcus.txt ****"
	set centcus = `cat $PERM/${inst}_centcus.txt`
	set x = $centcus[1]
	set y = $centcus[2]
    endif


#### extract the sky image using XSELECT ####

    echo "gnt_infrm (t_cent) creating sky image for $inst"

    xselect @$XCO/do_inst_sky_cent.xco $obsid $inst



#### Replace bad events in pixel 1,1 with a zero value ####

    fparimg 0.0 $WORK/${obsid}_${inst}_sky_cent.img 1,1 


#### create an image using XIMAGE ####

    set linewidth = 1

    echo "gnt_infrm (t_cent) creating ximage script for $inst"

    echo "read/${obsid}_${inst}_sky_cent.img" >! $WORK/${inst}centsky.xco
    echo "cey 2000" >> $WORK/${inst}centsky.xco
    echo "cpd ${obsid}_${inst}_sky_cent.ppm/ppm" >> $WORK/${inst}centsky.xco
    echo "cct/set $coltab" >> $WORK/${inst}centsky.xco
    echo "levels/num=128" >> $WORK/${inst}centsky.xco
    echo 'title " " ' >> $WORK/${inst}centsky.xco
    echo 'title/lower " " ' >> $WORK/${inst}centsky.xco

# SIS0 and SIS1
	
    if ( ($inst == "sis0") || ($inst == "sis1") ) then

        set sisx = `echo $x | awk '{print $1*4.0}'`
	set sisy = `echo $y | awk '{print $1*4.0}'`
	set sisrad = `echo $rad | awk '{print $1*4.0}'`

	echo "disp/sqrt nobox v1=0.1 v2=0.9 v3=0.1 v4=0.9" >> $WORK/${inst}centsky.xco
	echo "draw/circle $sisx $sisy $sisrad color=3 lstyle=1 lwidth=$linewidth" >> $WORK/${inst}centsky.xco
	echo "quit" >> $WORK/${inst}centsky.xco
	    
	echo "gnt_infrm (t_cent) running ximage $inst"

	ximage @$WORK/${inst}centsky        # run XIMAGE

	echo "gnt_infrm (t_cent) converting image $inst"

	if ( -e $WORK/${obsid}_${inst}_sky_cent.ppm ) then
	    convert ppm:${obsid}_${inst}_sky_cent.ppm gif:${obsid}_${inst}_sky_cent.gif
	    rm ${obsid}_${inst}_sky_cent.ppm
	endif

	if ( ! -e $WORK/${obsid}_${inst}_sky_cent.gif ) then
	    echo "gnt_fatal (t_cent) centering $inst sky image not produced"
	    exit 1
	else
	    convert -shave 235x150 gif:${obsid}_${inst}_sky_cent.gif gif:${obsid}_sky_1.gif
	    convert -resize 500x500 gif:${obsid}_sky_1.gif gif:${obsid}_sky_2.gif
	    convert -resize 75x75 gif:${obsid}_sky_1.gif gif:${obsid}_sky_3.gif
	    convert -quality 80 gif:${obsid}_sky_2.gif jpeg:${obsid}_${inst}_sky_cent_large.jpeg
	    convert -quality 80 gif:${obsid}_sky_3.gif jpeg:${obsid}_${inst}_sky_cent_thumb.jpeg
	    rm ${obsid}_sky_1.gif
	    rm ${obsid}_sky_2.gif
	    rm ${obsid}_sky_3.gif
	endif

    endif


# GIS2 and GIS3

    if ( ($inst == "gis2") || ($inst == "gis3") ) then

	echo "disp/sqrt nobox v1=0 v2=1.0 v3=0.0 v4=1.0" >> $WORK/${inst}centsky.xco
	echo "draw/circle $x $y $rad color=3 lstyle=1 lwidth=$linewidth" >> $WORK/${inst}centsky.xco
	echo "quit" >> $WORK/${inst}centsky.xco
	    
	echo "gnt_infrm (t_cent) running ximage $inst"

	ximage @$WORK/${inst}centsky

	echo "gnt_infrm (t_cent) converting image $inst"

	if ( -e $WORK/${obsid}_${inst}_sky_cent.ppm ) then
	    convert ppm:${obsid}_${inst}_sky_cent.ppm gif:${obsid}_${inst}_sky_cent.gif
	    rm ${obsid}_${inst}_sky_cent.ppm
	endif

	if ( ! -e $WORK/${obsid}_${inst}_sky_cent.gif ) then
	    echo "gnt_fatal (t_cent) centering $inst sky image not produced"
	    exit 1
	else
	    convert -shave 175x90 gif:${obsid}_${inst}_sky_cent.gif gif:${obsid}_sky_1.gif
	    convert -resize 200x200 gif:${obsid}_sky_1.gif gif:${obsid}_sky_3.gif
	    convert -quality 80 gif:${obsid}_sky_1.gif jpeg:${obsid}_${inst}_sky_cent_large.jpeg
	    convert -quality 80 gif:${obsid}_sky_3.gif jpeg:${obsid}_${inst}_sky_cent_thumb.jpeg
	    rm ${obsid}_sky_1.gif
	    rm ${obsid}_sky_3.gif
	endif

    endif

end


########################################################################################
# Create combined gis2 gis3 sky image with XY binning of 2
########################################################################################

# This sky images will have the nominal object position shown with a circle. The image
# may be used for a visual inspection of the sources in the field-of-view of each instrument.


#### extract the sky images using XSELECT and combine images ####

echo "gnt_infrm (t_cent) creating rebinned combined gis2 gis3 sky image"

echo "gnt_infrm (t_cent) extracting gis2 rebinned sky image"
xselect @$XCO/do_gis2_sky_bin2.xco $obsid

echo "gnt_infrm (t_cent) extracting gis3 rebinned sky image"
xselect @$XCO/do_gis3_sky_bin2.xco $obsid

echo "gnt_infrm (t_cent) adding gis2 and gis3 rebinned sky images"
farith infil1=$WORK/${obsid}_gis2_sky_bin2.img infil2=$WORK/${obsid}_gis3_sky_bin2.img \
    outfil=$WORK/${obsid}_gis_sky_bin2.img ops=ADD clobber=yes

fparimg 0.0 $WORK/${obsid}_gis_sky_bin2.img 1,1    # replace this bad pixel with zero


#### find nominal centroid for GIS this will be shown on the image ####

set coord = `sky2xy $WORK/${obsid}_gis2.evt+1 xsky=$ra ysky=$dec | grep pixel | sed s/","//g`
set rad = 14
set x = `echo $coord[4] | awk '{print int($1)}'`     # nominal RA in unbinned sky pixels 
set y = `echo $coord[5] | awk '{print int($1)}'`     # nominal DEC in unbinned sky pixels
	

#### show the custom centroid circle if being used to find centroid ####

if ( -e $PERM/${inst}_centcus.txt ) then
    echo "gnt_infrm (t_cent) **** SEQUENCE $obsid - using custom centroiding circle file $PERM/${inst}_centcus.txt ****"
    echo "gnt_warng (t_cent) **** SEQUENCE $obsid - using custom centroiding circle file $PERM/${inst}_centcus.txt ****"
    set centcus = `cat $PERM/${inst}_centcus.txt`
    set x = $centcus[1]
    set y = $centcus[2]
endif


#### create an image using XIMAGE ####

set linewidth = 1

echo "gnt_infrm (t_cent) creating ximage script for gis_sky_bin2"

echo "read/${obsid}_gis_sky_bin2.img" >! $WORK/gis_sky_bin2.xco
echo "cey 2000" >> $WORK/gis_sky_bin2.xco
echo "cpd ${obsid}_gis_sky_bin2.ppm/ppm" >> $WORK/gis_sky_bin2.xco
echo "cct/set $coltab" >> $WORK/gis_sky_bin2.xco
echo "levels/num=128" >> $WORK/gis_sky_bin2.xco
echo 'title " " ' >> $WORK/gis_sky_bin2.xco
echo 'title/lower " " ' >> $WORK/gis_sky_bin2.xco
echo "disp/sqrt nobox v1=0 v2=1.0 v3=0.0 v4=1.0" >> $WORK/gis_sky_bin2.xco
echo "draw/circle $x $y $rad color=3 lstyle=1 lwidth=$linewidth" >> $WORK/gis_sky_bin2.xco
echo "quit" >> $WORK/gis_sky_bin2.xco
	    
echo "gnt_infrm (t_cent) running ximage for gis_sky_bin2"

ximage @$WORK/gis_sky_bin2

echo "gnt_infrm (t_cent) converting image for gis_sky_bin2"

if ( -e $WORK/${obsid}_gis_sky_bin2.ppm ) then
    convert ppm:${obsid}_gis_sky_bin2.ppm gif:${obsid}_gis_sky_bin2.gif
    rm ${obsid}_gis_sky_bin2.ppm
endif

if ( ! -e $WORK/${obsid}_gis_sky_bin2.gif ) then
    echo "gnt_fatal (t_cent) gis_sky_bin2 sky image not produced"
    exit 1
else
    convert -shave 175x90 gif:${obsid}_gis_sky_bin2.gif gif:${obsid}_sky_1.gif
    convert -resize 200x200 gif:${obsid}_sky_1.gif gif:${obsid}_sky_3.gif
    convert -quality 80 gif:${obsid}_sky_1.gif jpeg:${obsid}_gis_sky_bin2_large.jpeg
    convert -quality 80 gif:${obsid}_sky_3.gif jpeg:${obsid}_gis_sky_bin2_thumb.jpeg
    rm ${obsid}_sky_1.gif
    rm ${obsid}_sky_3.gif
endif


########################################################################################
# Create Veron image - first call
########################################################################################

# create an image before checking for a detection

echo "gnt_infrm (t_cent) creating veron image"

$SCRIPTS/t_veron $obsid $SEQ $SCRIPTS $VERON


########################################################################################
# Extract source and background spectra from GIS2
########################################################################################

echo "gnt_infrm (t_cent) extracting GIS2 source and background spectra"

#### create nominal extraction regions ####

set coord = `sky2xy $WORK/${obsid}_gis2.evt+1 xsky=$ra ysky=$dec | grep pixel | sed s/","//g`

set x = `echo $coord[4] | awk '{print int($1)}'`
set y = `echo $coord[5] | awk '{print int($1)}'`

echo "gnt_infrm (t_cent) RA in sky pixels  (X) = $x"
echo "gnt_infrm (t_cent) DEC in sky pixels (Y) = $y"

set rad = 14
echo " circle($x,$y,$rad)" >! $WORK/gis2_src_reg.tmp      # write source and bg region file
echo " circle($x,$y,55.0)" >! $WORK/gis2_bgd_reg.tmp      # in sky pixels
echo "-circle($x,$y,35.0)" >> $WORK/gis2_bgd_reg.tmp


#### check for custom regions and use if present ####

if ( -e $PERM/gis2_srccus.reg ) then
    echo "gnt_infrm (t_cent) **** SEQUENCE $obsid - using custom source region for GIS2 detection: $PERM/gis2_srccus.reg ****"
    echo "gnt_warng (t_cent) **** SEQUENCE $obsid - using custom source region for GIS2 detection: $PERM/gis2_srccus.reg ****"
    cp $PERM/gis2_srccus.reg $WORK/gis2_src_reg.tmp
endif

if ( -e $PERM/gis2_bgdcus.reg ) then
    echo "gnt_infrm (t_cent) **** SEQUENCE $obsid - using custom background region for GIS2 detection: $PERM/gis2_bgdcus.reg ****"
    echo "gnt_warng (t_cent) **** SEQUENCE $obsid - using custom background region for GIS2 detection: $PERM/gis2_bgdcus.reg ****"
    cp $PERM/gis2_bgdcus.reg $WORK/gis2_bgd_reg.tmp
endif


#### copy extraction regions to permanent directory ####

cp $WORK/gis2_src_reg.tmp $PERM
cp $WORK/gis2_bgd_reg.tmp $PERM


#### extract source and background spectra ####

xselect @$XCO/do_cent_extr_spec.xco $obsid   # .xco script uses energy channel 60-900 = 0.5-10 keV
                                             # extracts a source and bg spectrum using the rg files


########################################################################################
# Create unsmoothed GIS2 sky image - show the extraction region
########################################################################################


#### create unsmoothed image ####

echo "gnt_infrm (t_cent) extracting initial GIS2 sky image"
echo "-circle(1,1,10)" >! gis2_excl_corner.reg   # don't want this contamination to affect colour scale
xselect @$XCO/do_gis2_sky_init.xco $obsid 


#### start the XIMAGE script ####

echo "gnt_infrm (t_cent) creating ximage script"

set linewidth = 1

echo "read/${obsid}_gis2_sky_init.img" >! $WORK/gis2initsky.xco
echo "cey 2000" >> $WORK/gis2initsky.xco
echo "cpd ${obsid}_gis2_sky_init.ppm/ppm" >> $WORK/gis2initsky.xco
echo "cct/set $coltab" >> $WORK/gis2initsky.xco
echo "levels/num=128" >> $WORK/gis2initsky.xco
echo 'title " " ' >> $WORK/gis2initsky.xco
echo 'title/lower " " ' >> $WORK/gis2initsky.xco
echo "disp/sqrt nobox v1=0 v2=1.0 v3=0.0 v4=1.0" >> $WORK/gis2initsky.xco


#### draw the extraction regions - get coordinates from the region files ####

foreach regfile ( $WORK/gis2_src_reg.tmp $WORK/gis2_bgd_reg.tmp )
    
    set line = 1
    set xreg = `sed s/"("/" "/g $regfile | sed s/")"/" "/g | sed s/","/" "/g | sed -n ${line}p`

    while ( "$#xreg" != "0" )
	
        if ( "$xreg[1]" == "box" ) then
	    if ( "$#xreg" < "6" ) set xreg = ($xreg 0.0)
	    set xreg = ( `echo "scale=3;$xreg[2]" | bc` `echo "scale=3;$xreg[3]" | bc` \
			 `echo "scale=3;$xreg[4]" | bc` `echo "scale=3;$xreg[5]" | bc` \
			 `echo "scale=3;$xreg[6]" | bc` )
	    echo "draw/box $xreg[1] $xreg[2] $xreg[3] $xreg[4] $xreg[5] color=2 lstyle=1 lwidth=$linewidth" >> $WORK/gis2initsky.xco
	endif


	if ( "$xreg[1]" == "-box" ) then
	    if ( "$#xreg" < "6" ) set xreg = ($xreg 0.0)
	    set xreg = ( `echo "scale=3;$xreg[2]" | bc` `echo "scale=3;$xreg[3]" | bc` \
			 `echo "scale=3;$xreg[4]" | bc` `echo "scale=3;$xreg[5]" | bc` \
			 `echo "scale=3;$xreg[6]" | bc` )
	    echo "draw/box $xreg[1] $xreg[2] $xreg[3] $xreg[4] $xreg[5] color=7 lstyle=2 lwidth=$linewidth" >> $WORK/gis2initsky.xco
	endif


	if ( "$xreg[1]" == "circle" ) then
	    set xreg = ( `echo "scale=3;$xreg[2]" | bc` `echo "scale=3;$xreg[3]" | bc` \
			 `echo "scale=3;$xreg[4]" | bc` )
	    if ( $regfile == "$WORK/gis2_bgd_reg.tmp" ) then  # outside of bgd annulus
		echo "draw/circle $xreg[1] $xreg[2] $xreg[3] color=2 lstyle=1 lwidth=$linewidth" >> $WORK/gis2initsky.xco
	    else
		echo "draw/circle $xreg[1] $xreg[2] $xreg[3] color=3 lstyle=1 lwidth=$linewidth" >> $WORK/gis2initsky.xco
	    endif
	endif


	if ( "$xreg[1]" == "-circle" ) then
	    set xreg = ( `echo "scale=3;$xreg[2]" | bc` `echo "scale=3;$xreg[3]" | bc` \
		         `echo "scale=3;$xreg[4]" | bc` )
	    echo "draw/circle $xreg[1] $xreg[2] $xreg[3] color=7 lstyle=2 lwidth=$linewidth" >> $WORK/gis2initsky.xco
	endif


	if ( "$xreg[1]" == "-annulus" ) then
	    set xreg = ( `echo "scale=3;$xreg[2]" | bc` `echo "scale=3;$xreg[3]" | bc` \
		         `echo "scale=3;$xreg[4]" | bc` )
	    echo "draw/circle $xreg[1] $xreg[2] $xreg[3] color=7 lstyle=2 lwidth=$linewidth" >> $WORK/gis2initsky.xco
	endif


	@ line ++
	set xreg = `sed s/"("/" "/g $regfile | sed s/")"/" "/g | sed s/","/" "/g | sed -n ${line}p`

    end

end


#### finish the XIMAGE script ####

echo "quit" >> $WORK/gis2initsky.xco


#### run XIMAGE script ####

echo "gnt_infrm (t_cent) running ximage"

ximage @$WORK/gis2initsky


#### convert ppm image to jpeg for web page ####

echo "gnt_infrm (t_cent) converting image"

if ( -e $WORK/${obsid}_gis2_sky_init.ppm ) then
    convert ppm:${obsid}_gis2_sky_init.ppm gif:${obsid}_gis2_sky_init.gif
    rm ${obsid}_gis2_sky_init.ppm
endif

if ( ! -e $WORK/${obsid}_gis2_sky_init.gif ) then
    echo "gnt_fatal (t_cent) initial GIS2 sky image not produced"
    exit 1
else
    convert -shave 175x90 gif:${obsid}_gis2_sky_init.gif gif:${obsid}_sky_1.gif
    convert -resize 200x200 gif:${obsid}_sky_1.gif gif:${obsid}_sky_3.gif
    convert -quality 80 gif:${obsid}_sky_1.gif jpeg:${obsid}_gis2_sky_init_large.jpeg
    convert -quality 80 gif:${obsid}_sky_3.gif jpeg:${obsid}_gis2_sky_init_thumb.jpeg
    rm ${obsid}_sky_1.gif
    rm ${obsid}_sky_3.gif
endif


########################################################################################
# Create smoothed GIS2 sky image - show the extraction region
########################################################################################


#### smooth image ####

echo "gnt_infrm (t_cent) smoothing initial sky image"

fgauss infile=$WORK/${obsid}_gis2_sky_init.img outfile=$WORK/${obsid}_gis2_sky_init_sm.img \
    sigma=1.0 datatype=e clobber=yes


#### start XIMAGE script ####

echo "gnt_infrm (t_cent) creating ximage script"

echo "read/${obsid}_gis2_sky_init_sm.img" >! $WORK/gis2initsky_sm.xco
echo "cey 2000" >> $WORK/gis2initsky_sm.xco
echo "cpd ${obsid}_gis2_sky_init_sm.ppm/ppm" >> $WORK/gis2initsky_sm.xco
echo "cct/set $coltab" >> $WORK/gis2initsky_sm.xco
echo "levels/num=128" >> $WORK/gis2initsky_sm.xco
echo 'title " " ' >> $WORK/gis2initsky_sm.xco
echo 'title/lower " " ' >> $WORK/gis2initsky_sm.xco
echo "disp/sqrt nobox v1=0 v2=1.0 v3=0.0 v4=1.0" >> $WORK/gis2initsky_sm.xco


#### draw the extraction regions - get coordinates from the region files ####

foreach regfile ( $WORK/gis2_src_reg.tmp $WORK/gis2_bgd_reg.tmp )
    
    set line = 1
    set xreg = `sed s/"("/" "/g $regfile | sed s/")"/" "/g | sed s/","/" "/g | sed -n ${line}p`

    while ( "$#xreg" != "0" )
	
	if ( "$xreg[1]" == "box" ) then
	    if ( "$#xreg" < "6" ) set xreg = ($xreg 0.0)
	    set xreg = ( `echo "scale=3;$xreg[2]" | bc` `echo "scale=3;$xreg[3]" | bc` \
			    `echo "scale=3;$xreg[4]" | bc` `echo "scale=3;$xreg[5]" | bc` \
			    `echo "scale=3;$xreg[6]" | bc` )
	    echo "draw/box $xreg[1] $xreg[2] $xreg[3] $xreg[4] $xreg[5] color=2 lstyle=1 lwidth=$linewidth" >> $WORK/gis2initsky_sm.xco
	endif


	if ( "$xreg[1]" == "-box" ) then
	    if ( "$#xreg" < "6" ) set xreg = ($xreg 0.0)
	    set xreg = ( `echo "scale=3;$xreg[2]" | bc` `echo "scale=3;$xreg[3]" | bc` \
			    `echo "scale=3;$xreg[4]" | bc` `echo "scale=3;$xreg[5]" | bc` \
			    `echo "scale=3;$xreg[6]" | bc` )
	    echo "draw/box $xreg[1] $xreg[2] $xreg[3] $xreg[4] $xreg[5] color=7 lstyle=2 lwidth=$linewidth" >> $WORK/gis2initsky_sm.xco
	endif


	if ( "$xreg[1]" == "circle" ) then
	    set xreg = ( `echo "scale=3;$xreg[2]" | bc` `echo "scale=3;$xreg[3]" | bc` \
			    `echo "scale=3;$xreg[4]" | bc` )
	    if ( $regfile == "$WORK/gis2_bgd_reg.tmp" ) then  # outside of bgd annulus
		echo "draw/circle $xreg[1] $xreg[2] $xreg[3] color=2 lstyle=1 lwidth=$linewidth" >> $WORK/gis2initsky_sm.xco
	    else
		echo "draw/circle $xreg[1] $xreg[2] $xreg[3] color=3 lstyle=1 lwidth=$linewidth" >> $WORK/gis2initsky_sm.xco
	    endif
	endif


	if ( "$xreg[1]" == "-circle" ) then
	    set xreg = ( `echo "scale=3;$xreg[2]" | bc` `echo "scale=3;$xreg[3]" | bc` \
			    `echo "scale=3;$xreg[4]" | bc` )
	    echo "draw/circle $xreg[1] $xreg[2] $xreg[3] color=7 lstyle=2 lwidth=$linewidth" >> $WORK/gis2initsky_sm.xco
	endif


	if ( "$xreg[1]" == "-annulus" ) then
	    set xreg = ( `echo "scale=3;$xreg[2]" | bc` `echo "scale=3;$xreg[3]" | bc` \
			    `echo "scale=3;$xreg[4]" | bc` )
	    echo "draw/circle $xreg[1] $xreg[2] $xreg[3] color=7 lstyle=2 lwidth=$linewidth" >> $WORK/gis2initsky_sm.xco
	endif


	@ line ++
	set xreg = `sed s/"("/" "/g $regfile | sed s/")"/" "/g | sed s/","/" "/g | sed -n ${line}p`


    end

end


#### finish the XIMAGE script ####

echo "quit" >> $WORK/gis2initsky_sm.xco


#### run XIMAGE script ####

echo "gnt_infrm (t_cent) running ximage"

ximage @$WORK/gis2initsky_sm


#### convert ppm image to jpeg for web page ####

echo "gnt_infrm (t_cent) converting image"

if ( -e $WORK/${obsid}_gis2_sky_init_sm.ppm ) then
    convert ppm:${obsid}_gis2_sky_init_sm.ppm gif:${obsid}_gis2_sky_init_sm.gif
    rm ${obsid}_gis2_sky_init_sm.ppm
endif

if ( ! -e $WORK/${obsid}_gis2_sky_init_sm.gif ) then
    echo "gnt_fatal (t_cent) smoothed initial GIS2 sky image not produced"
    exit 1
else
    convert -shave 175x90 gif:${obsid}_gis2_sky_init_sm.gif gif:${obsid}_sky_1.gif
    convert -resize 200x200 gif:${obsid}_sky_1.gif gif:${obsid}_sky_3.gif
    convert -quality 80 gif:${obsid}_sky_1.gif jpeg:${obsid}_gis2_sky_init_sm_large.jpeg
    convert -quality 80 gif:${obsid}_sky_3.gif jpeg:${obsid}_gis2_sky_init_sm_thumb.jpeg
    rm ${obsid}_sky_1.gif
    rm ${obsid}_sky_3.gif
endif


########################################################################################
# Copy the gis2 detection regions
########################################################################################

# these files will be included in the final products

cp $WORK/gis2_src_reg.tmp $WORK/gis2_det_src.reg
cp $WORK/gis2_bgd_reg.tmp $WORK/gis2_det_bgd.reg


########################################################################################
# Check for a detection in GIS2 
########################################################################################

#### determine number of counts in each spectrum ####

    set tot_src = "undef"
    set tot_bgd = "undef"

    set tot_src = `fkeyprint $WORK/gis2_spec_src.tmp+1 TOTCTS | grep "TOTCTS  =" | awk '{print $3}'`
    set tot_bgd = `fkeyprint $WORK/gis2_spec_bgd.tmp+1 TOTCTS | grep "TOTCTS  =" | awk '{print $3}'`

    echo "gnt_infrm (t_cent) total number of counts in source region     : $tot_src"
    echo "gnt_infrm (t_cent) total number of counts in background region : $tot_bgd"


#### determine the number of image pixels used for each spectrum ####

    fimgstat $WORK/gis2_spec_src.tmp threshlo=-0.000001 threshup=999999
    set num_src = `pget fimgstat num | awk '{print int($1)}'`    # number of source pixels

    fimgstat $WORK/gis2_spec_bgd.tmp threshlo=-0.000001 threshup=999999
    set num_bgd = `pget fimgstat num | awk '{print int($1)}'`    # number of background pixels

    set F = `echo "scale=3;$num_src / $num_bgd" | bc`

    echo "gnt_infrm (t_cent) total number of pixels in source region     : $num_src"
    echo "gnt_infrm (t_cent) total number of pixels in background region : $num_bgd"
    echo "gnt_infrm (t_cent) scaling factor                              : $F"


#### calculate the significance of detection number (SDS) for the source ####

    set Ts = $tot_src                              # counts in source region
    set Tb = `echo "scale=3; $tot_bgd * $F" | bc`  # normalised bgd counts

    echo "gnt_infrm (t_cent) estimated number of background counts in source region : $Tb"

    set Tbu = `echo $tot_bgd $F | awk '{print ((($1*$2*$2)+($1*$2))^0.5)}'` 
						# estimate of the uncertainty in the estimation of bgd in source region
                                                # this takes into account the fact that the bgd in the source region
                                                # is independent of the bgd in the background region, except
                                                # that they have the same population mean, having taken into 
                                                # account the scaling factor F
                                                                             
    echo "gnt_infrm (t_cent) estimated uncertainty of background counts in source region : $Tbu"

    set src = `echo $Ts $Tb | awk '{print ($1-$2)}'`              # estimated source counts in source region
    set sc = `echo $Ts $Tb $Tbu | awk '{print (($1-$2)/($3))}'`   # signal-to-bg noise
    set s = `echo $sc | awk '{print int($1*100+0.5)/100}'`        # S/N to 2 decimal places
    set S = `echo $s | awk '{print int($1*100)}'`                 # S/N times 100 (an integer)

    echo "gnt_infrm (t_cent) number of source counts in source region: $src"
    echo "gnt_infrm (t_cent) signal-to-background-noise of source (sc): $sc"
    echo "gnt_infrm (t_cent) signal-to-background-noise of source (s): $s"
    echo "gnt_infrm (t_cent) signal-to-background-noise of source (S): $S"

    echo $s >! $WORK/sds.temp           # create new file


#### calculate 3-sigma upper limit if object not detected ####

# the upper limit is defined as the best estimate of the source counts
# plus 3 times the propagated errors in the total count (Ts) and estimate of the
# mean of the background (F * tot_bgd)

    if ( $S < 500 ) then
	set diff = `echo $src | awk '{print ($1 < 0.0)}'`  # = 1 if $src < 0 
	if ( $diff == 1 ) set src = 0    # src < 0
	set ulim = `echo $src $Ts $tot_bgd $F | awk '{print ($1+(3.0*($2+($4^2)*$3)^0.5))}'`
	echo $ulim >! ulim.temp
    endif


#### report upper limit if object was not detected in GIS2 and exit ####

    if ( $S < 500 ) then
	echo "gnt_warng (t_cent) **** SEQUENCE $obsid - source was not detected with any confidence because S < 5  (S=$s) ****"
	echo "gnt_infrm (t_cent) source upper limit is $ulim counts"
	echo "gnt_fatal (t_cent) exiting because source was not detected"
	exit 1
    else
	echo "gnt_infrm (t_cent) source was detected with confidence because S = $s"
	echo "gnt_infrm (t_cent) the scripts will find centroids from smoothed sky images."
	echo "gnt_infrm (t_cent) the source coordinates will be X and Y of the centroid."
    endif


########################################################################################
# Get centroid for each instrument
########################################################################################

# centroiding is only perform if the object was detected


    foreach inst ( sis0 sis1 gis2 gis3 )

	echo "gnt_infrm (t_cent) finding centroid for $inst"


#### find initial centroid for the inst ####

	set coord = `sky2xy $WORK/${obsid}_${inst}.evt+1 xsky=$ra ysky=$dec | grep pixel | sed s/","//g`
   
	if ( $inst == sis0 || $inst == sis1 ) then
	    set rad = 33
	    set x = `echo $coord[4] | awk '{print int($1/4)}'`   # nominal RA in rebinned sky pixels
	    set y = `echo $coord[5] | awk '{print int($1/4)}'`   # nominal DEC in rebinned sky pixels
	else    # GIS
	    set rad = 14
	    set x = `echo $coord[4] | awk '{print int($1)}'`     # nominal RA in unbinned sky pixels 
	    set y = `echo $coord[5] | awk '{print int($1)}'`     # nominal DEC in unbinned sky pixels
	endif


	if ( -e $PERM/${inst}_centcus.txt ) then
	    echo "gnt_infrm (t_cent) **** SEQUENCE $obsid - using custom centroiding circle file $PERM/${inst}_centcus.txt ****"
	    echo "gnt_warng (t_cent) **** SEQUENCE $obsid - using custom centroiding circle file $PERM/${inst}_centcus.txt ****"
	    set centcus = `cat $PERM/${inst}_centcus.txt`
	    set x = $centcus[1]
	    set y = $centcus[2]
	endif

	echo "gnt_infrm (t_cent) centroiding region is circle($x,$y,$rad)"
	echo " circle($x,$y,$rad)" > ! $WORK/init_reg.tmp


#### Create and smooth sky image - make raw image, even for SIS ####

	echo "gnt_infrm (t_cent) creating sky image for $inst"

	xselect @$XCO/do_cent_extr.xco $obsid $SCRIPTS $inst

	echo "gnt_infrm (t_cent) smoothing sky image for $inst"

	fgauss infile=$WORK/${inst}_sky.tmp outfile=$WORK/${inst}_sky_sm.tmp \
	    sigma=5 datatype=e clobber=yes


#### Find centroid of smoothed sky image ####
 
	echo "gnt_infrm (t_cent) finding maximum pixel (centroid) for $inst"

	fimgstat infile=$WORK/${inst}_sky_sm.tmp threshlo=1.e-3 threshup=10000

	set x_sky = `pget fimgstat xmax | awk '{print int($1)}'`   # X sky coord of max pixel in smoothed image
	set y_sky = `pget fimgstat ymax | awk '{print int($1)}'`   # Y sky coord of max pixel in smoothed image


#### Now convert to detector co-ordinates ####

	echo "gnt_infrm (t_cent) initial sky posns are $x_sky, $y_sky"
	echo "gnt_infrm (t_cent) converting to det co-ordinates"

	set subevt = $WORK/${obsid}_${inst}_subevt.tmp 
	echo $subevt

	if ( $inst == "sis0" || $inst == "sis1" ) then

	    set x_sel = `echo $x_sky | awk '{print int($1*4)}'`
	    set y_sel = `echo $y_sky | awk '{print int($1*4)}'`
	    set selsize = 20

	else if ( $inst == "gis2" || $inst == "gis3" ) then

	    set x_sel = `echo $x_sky | awk '{print int($1)}'`
	    set y_sel = `echo $y_sky | awk '{print int($1)}'`
	    set selsize = 5

	else    # NO region

	    echo "gnt_fatal (t_cent) unknown instrument $inst - exiting"
	    exit 1
   
	endif 


	fselect infile=$WORK/${obsid}_${inst}.evt+1 outfile=$subevt expr="circle($x_sel,$y_sel,$selsize,x,y)"

	fstatistic infile=$subevt+1 colname=DETX rows=-
	set x_tmp = `pget fstatistic mean`

	fstatistic infile=$subevt+1 colname=DETY rows=-
	set y_tmp = `pget fstatistic mean`


	if ( $inst == sis0 || $inst == sis1 ) then    # SIS region
      
	    set x_det = `echo $x_tmp | awk '{print int($1/4)}'`
	    set y_det = `echo $y_tmp | awk '{print int($1/4)}'`
	    set csize = 45
 
	else if ( $inst == gis2 || $inst == gis3 ) then  # GIS region

	    set x_det = `echo $x_tmp | awk '{print int($1)}'`
	    set y_det = `echo $y_tmp | awk '{print int($1)}'`
	    set csize = 27
	
	else   # NO region

	    echo "gnt_fatal (t_cent) unknown instrument $inst - exiting"
	    exit 1

	endif

	echo "gnt_infrm (t_cent) $inst centroid DETX = $x_det"
	echo "gnt_infrm (t_cent) $inst centroid DETY = $y_det"


#### Construct the source region and detector centroid files ####

	if ( -e $PERM/${obsid}_${inst}_manual_src.reg ) then     # a custom region file is available

	    echo "gnt_infrm (t_cent) **** SEQUENCE $obsid - using custom region file $PERM/${obsid}_${inst}_manual_src.reg ****"
	    echo "gnt_warng (t_cent) **** SEQUENCE $obsid - using custom region file $PERM/${obsid}_${inst}_manual_src.reg ****"
	    cp $PERM/${obsid}_${inst}_manual_src.reg $WORK/${obsid}_${inst}_src.reg


 # confirm that custom region and centroid are the same

	    set centline = ( `grep "circle(${x_det},${y_det}" $WORK/${obsid}_${inst}_src.reg` )
	    if ( $#centline == 0 ) then
		echo "gnt_fatal (t_cent) The detector centroids just determined do not match the src region file"
		echo "gnt_fatal (t_cent) The custom region files might need to be changed - EXITING"
		exit 1
	    endif

	else

	    echo "gnt_infrm (t_cent) writing source region for $inst"
	    echo " circle($x_det,$y_det,$csize)" >! $WORK/${obsid}_${inst}_src.reg   # create new source region file
	    cat $REGS/${inst}_excl.reg >> $WORK/${obsid}_${inst}_src.reg

	endif

	echo "gnt_infrm (t_cent) writing detector centroid file for $inst"    
	echo "$x_det $y_det" >! $WORK/detcent_${inst}.temp        # *** this file is used by ascaarf ***


#### Construct the background region file ####

	if ( -e $PERM/${obsid}_${inst}_manual_bgd.reg ) then        # a custom region file is available

	    echo "gnt_infrm (t_cent) **** SEQUENCE $obsid - using custom region file $PERM/${obsid}_${inst}_manual_bgd.reg ****"
	    echo "gnt_warng (t_cent) **** SEQUENCE $obsid - using custom region file $PERM/${obsid}_${inst}_manual_bgd.reg ****"
	    cp $PERM/${obsid}_${inst}_manual_bgd.reg $WORK/${obsid}_${inst}_bgd.reg


# confirm that custom region and centroid are the same

	    if ( $inst == sis0 || $inst == sis1 ) then     
		set centline = ( `grep "-circle(${x_det},${y_det}" $WORK/${obsid}_${inst}_bgd.reg` )
	    else if ( $inst == gis2 || $inst == gis3 ) then 
		set centline = ( `grep "-circle(${x_det},${y_det}" $WORK/${obsid}_${inst}_bgd.reg` )
	    endif
	    if ( $#centline == 0 ) then   # if don't find the same coordinates
		echo "gnt_fatal (t_cent) The detector centroids just determined do not match the bgd region file"
		echo "gnt_fatal (t_cent) The custom region files might need to be changed - EXITING"
		exit 1
	    endif


# create new region file

	else    

	    if ( $inst == sis0 || $inst == sis1 ) then   # create new SIS region file   

		set csize = 47
		echo "gnt_infrm (t_cent) writing SIS background region for $inst"
		echo $x_det,$y_det
		cat $REGS/${inst}.reg >! $WORK/${obsid}_${inst}_bgd.reg
		echo "-circle($x_det,$y_det,$csize)" >> $WORK/${obsid}_${inst}_bgd.reg

	    else if ( $inst == gis2 || $inst == gis3 ) then   # create new GIS region file

		echo "gnt_infrm (t_cent) writing GIS background region for $inst"
		echo $x_det,$y_det
		echo " circle($x_det,$y_det,55.0)" >! $WORK/${obsid}_${inst}_bgd.reg
		echo "-circle($x_det,$y_det,35.0)" >> $WORK/${obsid}_${inst}_bgd.reg

	    endif

	    cat $REGS/${inst}_excl.reg >> $WORK/${obsid}_${inst}_bgd.reg  # exclusion region

	endif


#### Store sky image centroids and overall centroid ####
	
	if ( $S > 500 ) then
	    set ra_dec1 = `xy2sky ${inst}_sky_sm.tmp $x_sky $y_sky | grep SKY`
	else
	    set ra_dec1 = `echo 0 0 0 $ra $dec`
	endif

	set ra_dec = `echo $ra_dec $ra_dec1[4] $ra_dec1[5] | sed s/,//g`
	set cent = `echo $cent $x_sky $y_sky`


#### end the foreach inst loop ####

    end

    set ra_dec = `echo $ra_dec | awk '{print $1+0,$2+0,$3+0,$4+0,$5+0,$6+0,$7+0,$8+0}'`


########################################################################################
# Compare the SIS and GIS centroids if SOURCE WAS FOUND
########################################################################################


#### compare centroids ####

    echo "gnt_infrm (t_cent) comparing SIS and GIS centroids respectively"
	
    set gis23_ra = `echo $ra_dec[5] $ra_dec[7] | awk '{print int(($1-$2)*3600)}' | sed s/-//g`
    set gis23_dec = `echo $ra_dec[6] $ra_dec[8] | awk '{print int(($1-$2)*3600)}' | sed s/-//g`
    
    if ( $gis23_ra > 60 || $gis23_dec > 60 ) then
	echo "gnt_warng (t_cent) **** SEQUENCE $obsid - centroids for GIS intruments disagree more than 1 arcmin ****"
    endif

    set sis01_ra = `echo $ra_dec[1] $ra_dec[3] | awk '{print int(($1-$2)*3600)}' | sed s/-//g`
    set sis01_dec = `echo $ra_dec[2] $ra_dec[4] | awk '{print int(($1-$2)*3600)}' | sed s/-//g`
   
    if ( $sis01_ra > 60 || $sis01_dec > 60 ) then
	echo "gnt_warng (t_cent) **** SEQUENCE $obsid - centroids for SIS intruments disagree more than 1 arcmin ****"
    endif


#### Combine respective sky images and smooth ####

    echo "gnt_infrm (t_cent) combining respective sky images"
    farith infil1=gis2_sky.tmp infil2=gis3_sky.tmp outfil=gis_sky.tmp ops=ADD 
    fgauss infile=$WORK/gis_sky.tmp outfile=$WORK/gis_sky_sm.tmp \
	sigma=2 datatype=e clobber=yes


#### Find centroid of GIS and compare to both SIS0 and SIS1 ####

    echo "gnt_infrm (t_cent) comparing combined GIS centroid to SIS0 and SIS1"
	
    fimgstat infile=$WORK/gis_sky_sm.tmp threshlo=1.e-3 threshup=10000

    set gis_x = `pget fimgstat xmax`
    set gis_y = `pget fimgstat ymax`

    set gis_ra_dec = `xy2sky gis_sky_sm.tmp $gis_x $gis_y | grep SKY`

    set gis_sis0_ra = `echo $gis_ra_dec[4] $ra_dec[1] | sed s/,//g | awk '{print int(($1-$2)*3600)}' | sed s/-//g`

    set gis_sis0_dec = `echo $gis_ra_dec[5] $ra_dec[2] | sed s/,//g | awk '{print int(($1-$2)*3600)}' | sed s/-//g`

    set gis_sis1_ra = `echo $gis_ra_dec[4] $ra_dec[3] | sed s/,//g | awk '{print int(($1-$2)*3600)}' | sed s/-//g`

    set gis_sis1_dec = `echo $gis_ra_dec[5] $ra_dec[4] | sed s/,//g | awk '{print int(($1-$2)*3600)}' | sed s/-//g`
	
    if ( $gis_sis0_ra > 60 || $gis_sis0_dec > 60 ) then
	echo "gnt_warng (t_cent) **** SEQUENCE $obsid - SIS0 and GIS centroids disagree more than 1 arcmin ****"
    endif

    if ( $gis_sis1_ra > 60 || $gis_sis1_dec > 60 ) then
	echo "gnt_warng (t_cent) **** SEQUENCE $obsid - SIS1 and GIS centroids disagree more than 1 arcmin ****"
    endif


########################################################################################
# Write all centroids to file
########################################################################################

    echo $ra $dec >! $WORK/centroid.temp    # RA and DEC from attitude file      

    if ( $S > 500 ) then     # RA and DEC in degrees from maximum in combined GIS sky image
	echo $gis_ra_dec[4] $gis_ra_dec[5] | sed s/,//g | awk '{print $1+0,$2+0}' >> $WORK/centroid.temp   
    else
	echo $ra $dec >> $WORK/centroid.temp   # RA and DEC from attitude file 
        echo "gnt_infrm (t_cent) no centroid comparison."
    endif

    echo $ra_dec >> $WORK/centroid.temp     # RA and DEC in degrees of centroid  for all 4 instruments
    echo $cent >> $WORK/centroid.temp       # X and Y in sky pixels of centroid for all 4 instruments


########################################################################################
# Make copy of region files 
########################################################################################

# the custom regions in the $PERM directory are never over-written.

#### copy to permanent area ####

    echo "gnt_infrm (t_cent) copying .reg files to permanent area" 

    cp -f $WORK/centroid.temp $PERM
    cp -f $WORK/sds.temp $PERM
    
    foreach inst ( sis0 sis1 gis2 gis3 )
	foreach reg ( src bgd )
	    cp $WORK/${obsid}_${inst}_${reg}.reg $PERM
 	end
	cp $WORK/detcent_${inst}.temp $PERM
    end





########################################################################################
# Create Veron image - second call
########################################################################################

# create a new Veron image - show the centroid

echo "gnt_infrm (t_cent) creating veron image"

$SCRIPTS/t_veron $obsid $SEQ $SCRIPTS $VERON


########################################################################################
# Create some detector images with regions shown
########################################################################################

# runs the script t_detimg
# detector images .img were produced by t_scrn. Now produce some gif and jpeg
# with the extraction regions shown

$SCRIPTS/t_detimg $obsid $SEQ


########################################################################################
# Check if all files were created
########################################################################################

foreach file ( centroid sds object otime )
    if ( ! -e $WORK/$file.temp ) then
	echo "gnt_fatal (t_cent) file $file.temp was not created."
	exit 1
    endif
end


foreach inst ( sis0 sis1 gis2 gis3 )
    foreach reg ( src bgd )
	set file = ${obsid}_${inst}_${reg}.reg
	if ( ! -e $WORK/$file ) then
	    echo "gnt_fatal (t_cent) File $file was not created."
	    exit 1
	endif
    end
end


########################################################################################
# All done 
########################################################################################

echo ""
echo "gnt_infrm (t_cent) ++++++++++++++++++++++++++++++++++++++++++++++++++++++"
echo "gnt_infrm (t_cent)  t_cent 4.0  FINISHED "`date`
echo "gnt_infrm (t_cent) ++++++++++++++++++++++++++++++++++++++++++++++++++++++"
echo ""