#! /bin/tcsh -f

# t_veron 4.0

# Creates a smoothed GIS sky image with veron2001
# objects shown.

# Paul O'Neill

# Date in final form: 17 December 2004
 

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

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


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

set obsid = $1           # observation sequence number
set SEQ = $2             # root directory for current sequence
set SCRIPTS = $3         # directory containing script files
set VERON = $4           # 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


########################################################################################
# Check for necessary and existing files 
########################################################################################


#### check for files ####

foreach file ( $WORK/${obsid}_gis2_sky_cent.img $WORK/${obsid}_gis3_sky_cent.img $WORK/${obsid}_gis2.evt $WORK/${obsid}_gis3.evt $VERON )
    if ( ! -e $file ) then
	echo "gnt_fatal (t_veron) $file not found"
	exit 1
    endif
end


#### get name of attitude file ####

set attfiles = `ls $AUX/fa*`

if ( $#attfiles == "0" ) then
    echo "gnt_fatal (t_veron) attitude file not found"
    exit 1
else
    set attfile = $attfiles[1]
endif


#### remove existing files ####

if ( -e $WORK/giscentre.temp ) rm $WORK/giscentre.temp
if ( -e $WORK/${obsid}_veronobj.fits ) rm $WORK/${obsid}_veronobj.fits
if ( -e $WORK/${obsid}_veronobj_fov.fits ) rm $WORK/${obsid}_veronobj_fov.fits
if ( -e $WORK/${obsid}_veronobj_fov.temp ) rm $WORK/${obsid}_veronobj_fov.temp
if ( -e $WORK/${obsid}_gis_veron.img ) rm $WORK/${obsid}_gis_veron.img
if ( -e $WORK/${obsid}_gis_veron.ppm ) rm $WORK/${obsid}_gis_veron.ppm
if ( -e $WORK/${obsid}_gis_veron.gif ) rm $WORK/${obsid}_gis_veron.gif
if ( -e $WORK/${obsid}_gis_veron_large.gif ) rm $WORK/${obsid}_gis_veron_large.jpeg
if ( -e $WORK/${obsid}_gis_veron_thumb.gif ) rm $WORK/${obsid}_gis_veron_thumb.jpeg


########################################################################################
# Check is object was detected
########################################################################################

cd $WORK

if ( -e sds.temp ) then
    set s = `head -1 sds.temp`
    if ( $#s == "0" ) then
	echo "gnt_infrm (t_veron) the sds.temp file is empty"
	set detect = 0
    else
	set detect = `echo $s | awk '{print ($1 >= 5.0)}'`  # = 1 if $s >= 5.0
    endif
else
    set detect = 0
endif

if ( $detect == 1 ) then
    echo "gnt_infrm (t_veron) object detected with significance $s"
else
    echo "gnt_infrm (t_veron) object was not detected or no sds.temp file"
endif


########################################################################################
# Find centre of GIS2 and GIS3 fields 
########################################################################################

echo "gnt_infrm (t_veron)"
echo "gnt_infrm (t_veron) finding centre of combined GIS2 and GIS3 field-of-view"


#### centre of GIS2 field ####

set gis2X = "128.5"
set gis2Y = "128.5"

xy2sky infile=${obsid}_gis2.evt+1 xpix=$gis2X ypix=$gis2Y

set gis2ra = `pget xy2sky xsky`
set gis2dec = `pget xy2sky ysky`

echo "gnt_infrm (t_veron)"
echo "gnt_infrm (t_veron) centre of GIS2 FOV in (sky pixels): $gis2X $gis2Y"
echo "gnt_infrm (t_veron) centre of GIS2 FOV in (RA DEC degrees): $gis2ra $gis2dec"


#### centre of GIS3 field ###

set gis3X = "128.5"
set gis3Y = "128.5"

xy2sky infile=${obsid}_gis3.evt+1 xpix=$gis3X ypix=$gis3Y

set gis3ra = `pget xy2sky xsky`
set gis3dec = `pget xy2sky ysky`

echo "gnt_infrm (t_veron)"
echo "gnt_infrm (t_veron) centre of GIS3 FOV in (sky pixels): $gis3X $gis3Y"
echo "gnt_infrm (t_veron) centre of GIS3 FOV in (RA DEC degrees): $gis3ra $gis3dec"


#### get mean of GIS2 and GIS3 ####

set gisra = `echo "scale=7;( ${gis2ra} + ${gis3ra} ) / 2.0" | bc`
set gisdec = `echo "scale=7;( ${gis2dec} + ${gis3dec} ) / 2.0" | bc`

echo "gnt_infrm (t_veron)"
echo "gnt_infrm (t_veron) mean GIS FOV centre (RA DEC degrees): $gisra $gisdec"


#### write to file ####

touch giscentre.temp
echo "$gisra $gisdec" >> giscentre.temp
echo "$gis2X $gis2Y $gis2ra $gis2dec" >> giscentre.temp
echo "$gis3X $gis3Y $gis3ra $gis3dec" >> giscentre.temp


########################################################################################
# Get RA and DEC of target object in field
########################################################################################

echo "gnt_infrm (t_veron)"
echo "gnt_infrm (t_veron) determing location of object"


#### get nominal target coordinates from attitude file ####

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

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

set cenra = $nomra      # default values for centroid position
set cendec = $nomdec

echo "gnt_infrm (t_veron) nominal RA and DEC (degrees): $nomra $nomdec"


#### look for the object centroid position ####

set havecent = "no"      # do not have a centroid position yet

if ( ($detect == 1) && (-e $WORK/centroid.temp) ) then     # only use centroid if object was detected

    set centroid = `cat $WORK/centroid.temp`

    if ( $#centroid != "20" ) then   # there should be 20 numbers in centroid.temp
	echo "gnt_infrm (t_veron)"
	echo "gnt_infrm (t_veron) centroid file corrupt - showing nominal position"
    else
	set cenra = $centroid[3]
	set cendec = $centroid[4]
	echo "gnt_infrm (t_veron)"
	echo "gnt_infrm (t_veron) centroid RA and DEC (degrees): $cenra $cendec"
	set havecent = "yes"
    endif

else

    echo "gnt_infrm (t_veron)"
    echo "gnt_infrm (t_veron) centroid not available - showing nominal position"

endif


########################################################################################
# Find VERON2001 objects in field 
########################################################################################

cd $WORK

#### setup a VERON2001 FITS file in the WORK directory ####

cp $VERON $WORK/${obsid}_veronobj.fits

set searchrad = 22

fkeypar fitsfile=${obsid}_veronobj.fits+1 keyword=NAXIS2
set totnumobj = `pget fkeypar value`

echo "gnt_infrm (t_veron)"
echo "gnt_infrm (t_veron) finding VERON2001 objects within $searchrad arcmin GIS field-of-view radius"
echo "gnt_infrm (t_veron) number of objects in VERON2001 catalog: $totnumobj"


#### calculate offsets from FOV centre and centroid ####

# uses the Haversine formula to calculate the offset in arcmin 
# between each veron2001 object and the FOV and centroid

# FOV

set expr = "(10800.0 / #pi) * 2.0 * arcsin(min(1.0,sqrt((sin((DECRAD-(#deg * ${gisdec}))/2.0d0)**2) + ( cos((#deg * ${gisdec})) * cos(DECRAD) * (sin((RARAD-(#deg * ${gisra}))/2.0d0)**2)))))"

echo "gnt_infrm (t_veron)"
echo "gnt_infrm (t_veron) using Haversine formula"

fcalc infile=${obsid}_veronobj.fits outfile=${obsid}_veronobj.fits clname=OFFSETF  expr="${expr}" clobber=yes

# centroid

set expr = "(10800.0 / #pi) * 2.0 * arcsin(min(1.0,sqrt((sin((DECRAD-(#deg * ${cendec}))/2.0d0)**2) + ( cos((#deg * ${cendec})) * cos(DECRAD) * (sin((RARAD-(#deg * ${cenra}))/2.0d0)**2)))))"

fcalc infile=${obsid}_veronobj.fits outfile=${obsid}_veronobj.fits clname=OFFSETC  expr="${expr}" clobber=yes


#### check for veron objects in FOV ####

fselect infile=${obsid}_veronobj.fits outfile=${obsid}_veronobj_fov.fits expr="OFFSETF.LE.${searchrad}" clobber=yes

fkeypar fitsfile=${obsid}_veronobj_fov.fits+1 keyword=NAXIS2
set fovnumobj = `pget fkeypar value`

echo "gnt_infrm (t_veron)"
echo "gnt_infrm (t_veron) number of veron2001 objects in FOV: $fovnumobj"


#### sort according to centroid offset and dump ####

if ( $fovnumobj != "0" ) then
    fsort infile=${obsid}_veronobj_fov.fits+1 columns=OFFSETC method=heap ascend=yes
    fdump infile=${obsid}_veronobj_fov.fits+1 outfile=${obsid}_veronobj_fov.temp prhead=no prdata=yes showcol=no showunit=no showrow=yes pagewidth=200 columns="NAME,RA,DEC,REDSHIFT,TYPE,OFFSETF,OFFSETC" clobber=yes rows="-"
endif


########################################################################################
# Create combined smoothed GIS image
########################################################################################

# use exposure corrected image if available, otherwise combine the GIS2 and GIS3 
# sky images that are created for every sequence

if ( -e $WORK/${obsid}_gis_sky_sm.img ) then

    echo "gnt_infrm (t_veron) using existing smoothed GIS image"

    cp $WORK/${obsid}_gis_sky_sm.img $WORK/${obsid}_gis_veron.img

else

    echo "gnt_infrm (t_veron) creating combined smoothed GIS image"

    fparimg 0.0 $WORK/${obsid}_gis2_sky_cent.img 1,1
    fparimg 0.0 $WORK/${obsid}_gis3_sky_cent.img 1,1

    farith infil1=$WORK/${obsid}_gis2_sky_cent.img infil2=$WORK/${obsid}_gis3_sky_cent.img \
	outfil=$WORK/${obsid}_gis_veron.tmp ops=ADD clobber=yes

    fgauss infile=$WORK/${obsid}_gis_veron.tmp outfile=$WORK/${obsid}_gis_veron.img \
	sigma=2 datatype=e clobber=yes

endif


########################################################################################
# Make sky image of GIS for web page
########################################################################################

echo "gnt_infrm (t_veron) Making web page sky image with veron2001 objects"


#### start ximage script ####

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

#### show the object centroid/nominal position ####

if ( $havecent == "yes" ) then
    echo "label/ra=${cenra}/dec=${cendec}"'/symbol=25/symcolor=3/symcsize=1.0/symlwidth=3.0/just=left/text="C"'"/color=3/csize=1.5" >> gisveron.xco
else
    echo "label/ra=${cenra}/dec=${cendec}"'/symbol=25/symcolor=3/symcsize=1.0/symlwidth=3.0/just=left/text="N"'"/color=3/csize=1.5" >> gisveron.xco
endif


#### show each object in field-of-view ####

@ objnum = 1
while ( $objnum <= $fovnumobj )

    ftabpar fitsfile=${obsid}_veronobj_fov.fits+1 column=NAME row=$objnum
    set objname = "`pget ftabpar value`"

    ftabpar fitsfile=${obsid}_veronobj_fov.fits+1 column=RADEG row=$objnum
    set objra = "`pget ftabpar value`"

    ftabpar fitsfile=${obsid}_veronobj_fov.fits+1 column=DECDEG row=$objnum
    set objdec = "`pget ftabpar value`"

    if ( $objnum == "1" ) then
	echo "label/ra=${objra}/dec=${objdec}"'/symbol=18/symcolor=4/symcsize=2.0/symlwidth=2.0/just=center/text="'"${objname}"'"'"/color=4/csize=1.5/font=roman" >> gisveron.xco
    else
    echo "label/ra=${objra}/dec=${objdec}"'/symbol=6/symcolor=4/symcsize=1.5/symlwidth=1.5/just=center/text="'"${objname}"'"'"/color=4/csize=0.9/font=roman" >> gisveron.xco
    endif

    @ objnum ++

end


#### finish script and run ####

echo "quit" >> gisveron.xco
ximage @gisveron.xco


#### convert ppm to gif ####

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


#### create images for web page ####

if ( -e ${obsid}_gis_veron.gif ) then
    convert -shave 175x90 gif:${obsid}_gis_veron.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_veron_large.jpeg
    convert -quality 80 gif:${obsid}_sky_3.gif jpeg:${obsid}_gis_veron_thumb.jpeg
    rm ${obsid}_sky_1.gif
    rm ${obsid}_sky_3.gif
else
    echo "gnt_fatal (t_vernon) ${obsid}_gis_veron.gif does not exist"
    exit 1
endif


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

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