#! /bin/tcsh 
##### NOTE: USER MAY HAVE TO EDIT THIS LINE TO FIND LOCAL tcsh
#
#
# 1.0 AMR 27/05/03 - conversion of BGrebinimage2SKY to 4" - for exposure maps
# 1.1 AMR 27/01/05 - improvement in transform, removed references to local 
#     disks etc.
# 1.2 AMR 08/03/06 - Slight changes in (pn) coordinates for new 2006 BGWG 
#     (JAC/AMR) files + name change to BGrebinimage2SKY_2006
# 1.3 AMR 14/03/06 - Includes routine to catch and repair NAN values
# 1.4 JBB 07/02/07 - Accepts any image size, new scaling parameter
# 1.4 JBB 03/10/07 - Adapt NAN correction to Solaris output
# 1.5 AMR 04/10/07 - Change nawk to awk & correct no-NaN error
set version=1.5

#nice +19

cat <<EOF 

--------------- $0 VERSION $version `date +%x`  ---------

EOF

if ( $#argv != 6 ) then
 cat <<EOF 

Script to convert a detector (DETX/DETY) map into a sky (X/Y) image.
The input DETX/DETY map is rebinned and reprojected onto the sky to
the same spatial scale and sky position of a template image. The user
provides this template image (constructed using the evselect or
xmmselect SAS tasks), containing the attitude information, and an
event file (ideally, the one used to create the template image file,
though any may do [just used for general header purposes]). A rebinned
map is produced with the same spatial resolution, and at the same sky
position as the template image. The radius of the interpolation circle
considered for the re-griding can be given. Note that the routine
requires SAS and ftools. Also the SAS environment variables SAS_CCF
and SAS_CCFPATH need to be correctly set. This script calls the
interpolation routine BGrebinimage2SKY_interp which must be present in
the current working directory.

Specific notes: 

- Because of the size of the images, the task is CPU-intensive, and
may take a long time (5 minutes for a 250x250 input map, 20 minutes
for a 500x500 input map).

- Scaling of the map intensity values to account for changes in bin size is
performed optionally.

- The user may have to edit the first line of this script to access
their local tcsh installation.

- Some trigonometric calls within awk may not work on some unix
versions/platforms. The user may wish to substitute awk with nawk,
a newer version of the awk pattern scanning and processing language.

-----------------------------------------------------------------------------

Use: $0 image imtemplate evfile rinterp outimage scaling

 image      - Input map in DETX/DETY coordinates
 imtemplate - Template image file containing attitude info (high-res SKY image)
 evfile     - Event file, ideally the one used to create imtemplate
 rinterp    - Interpolation radius [arcmin] 
 outimage   - Output image (rebinned map in X/Y sky coordinates)
 scaling    - 0 means no scaling, 1 means rescaling according to pixel size

e.g. $0 EXPim.fits myimage.fits evfile.fits 0.15 outexp.fits 0

EOF
 exit
endif

set im=$1
set templ=$2
set ev=$3
set rinterp=$4
set out=$5
set scaling=$6

set radian=57.29577951

if (! -e $im) then 
  echo 'Error: Input image does not exist... exiting...'
  exit
endif
if (! -e $templ ) then 
 echo "Error: Template file does not exist... exiting..."
 exit
endif
if (! -e $ev ) then 
 echo "Error: Events file does not exist... exiting..."
 exit
endif
if (! -e $SAS_CCF ) then 
 echo "Error: SAS_CCF environment variable not set... exiting..."
 exit
endif

echo 'Creating events list from image...'
fim2lst $im temptbl2 copyall=Y clobber=Y
fdump $im tempdmp2 X - prdata=no showcol=no showunit=no showrow=no clobber=yes
set inx=`grep NAXIS1 tempdmp2 | head -1 | awk '{print $3}'`
set iny=`grep NAXIS2 tempdmp2 | head -1 | awk '{print $3}'`
echo 'Image is '$inx' by' $iny' pixels...'

set ctype=`grep CTYPE1 tempdmp2 | head -1 | awk '{print $3}'`
if ($ctype == "'DETX_ANG'") then
  set detbin=`grep CDELT1 tempdmp2 | head -1 | awk '{print int(72001*$3)}'`
  set detpixx=`grep CRPIX1 tempdmp2 | head -1 | awk '{print $3}'`
  set dxmin=`echo $detpixx $detbin | awk '{print int((1-$1)*$2-0.49)}'`
  set detpixy=`grep CRPIX2 tempdmp2 | head -1 | awk '{print $3}'`
  set dymin=`echo $detpixy $detbin | awk '{print int((1-$1)*$2-0.49)}'`
else if ($ctype == "'DETX") then
  set detbin=`grep CDELT1 tempdmp2 | head -1 | awk '{print int($3+0.49)}'`
  set dxmin=`grep CRVAL1 tempdmp2 | head -1 | awk '{print int($3-0.49)}'`
  set dymin=`grep CRVAL2 tempdmp2 | head -1 | awk '{print int($3-0.49)}'`
else
  echo 'Unrecognised file type ...exiting...' $ctype
  exit
endif

set dxmax=`echo $dxmin $detbin $inx | awk '{print int($1+$2*($3-1)-0.99)}'`
set dymax=`echo $dymin $detbin $iny | awk '{print int($1+$2*($3-1)-0.99)}'`
set dxmin2=`echo $dxmin $detbin | awk '{print $1-$2}'`
set dymin2=`echo $dymin $detbin | awk '{print $1-$2}'`
fcalc temptbl2 temptbl2 DETX "${dxmin2} + (X*${detbin})" clobber=yes
fcalc temptbl2 temptbl2 DETY "${dymin2} + (Y*${detbin})" clobber=yes
fdump temptbl2 templst2 columns="DETX,DETY,VALUE" rows=- prhead=no prdata=yes showcol=no showunit=no showrow=no clobber=yes

if (-e tempevl2) rm -f tempevl2
if (-e tempfcr2) rm -f tempfcr2
touch tempfcr2
echo "DETX I" >> tempfcr2
echo "DETY I" >> tempfcr2
echo "VAL D" >> tempfcr2
echo "" >> tempfcr2
fcreate tempfcr2 templst2 tempevl2 extname="EVENTS" clobber=yes
cphead $ev tempevl2 
fcalc tempevl2 tempevl2 TIME "1.0" clobber=yes

if (-e tempfmd2) rm -f tempfmd2
touch tempfmd2
echo "TLMIN1  =               ${dxmin} / Minimum for DETX" >> tempfmd2
echo "TLMAX1  =               ${dxmax} / Maximum for DETX" >> tempfmd2
echo "TCRPX1  =                    0 / Center for DETX" >> tempfmd2
echo "TCDLT1  = 1.38888888888889E-05 / [degrees] Pixel size for DETX" >> tempfmd2
echo "TLMIN2  =               ${dymin} / Minimum for DETY" >> tempfmd2
echo "TLMAX2  =               ${dymax} / Maximum for DETY" >> tempfmd2
echo "TCRPX2  =                    0 / Center for DETY" >> tempfmd2
echo "TCDLT2  = 1.38888888888889E-05 / [degrees] Pixel size for DETY" >> tempfmd2
echo "" >> tempfmd2
fmodhead tempevl2+1 tempfmd2

if (! -e $templ ) then 
 echo "Error: Template file does not exist... exiting..."
 exit
else 
 echo "Getting attitude information from template file..."
 if (-e tempdmp2) rm -f tempdmp2
 fdump $templ+0 tempdmp2 1 1-1
 set ra_pnt=`grep RA_PNT tempdmp2 | head -1 | awk '{print $3}'`
 set dec_pnt=`grep DEC_PNT tempdmp2 | head -1 | awk '{print $3}'`
 set pa_pnt=`grep PA_PNT tempdmp2 | head -1 | awk '{print $3}'`
 set ranom=`grep REFXCRVL tempdmp2 | head -1 | awk '{print $2}'`
 set decnom=`grep REFYCRVL tempdmp2 | head -1 | awk '{print $2}'`
endif

echo "Calculating sky coordinates... calling XMM-SAS attcalc..."
attcalc -w 0 -V 0 eventset=tempevl2 attitudelabel=fixed fixedra=$ra_pnt fixeddec=$dec_pnt fixedposangle=$pa_pnt withatthkset=N refpointlabel=user nominalra=$ranom nominaldec=$decnom

fdump $templ tempdmp3 X - prdata=no showcol=no showunit=no showrow=no clobber=yes
set tnx=`grep NAXIS1 tempdmp3 | head -1 | awk '{print $3}'`
set tny=`grep NAXIS2 tempdmp3 | head -1 | awk '{print $3}'`
echo 'Image is '$tnx' by '$tny' pixels...'
set tpx_o=`grep CDELT1 tempdmp3 | head -1 | awk '{print $3}'`
set tpy_o=`grep CDELT2 tempdmp3 | head -1 | awk '{print $3}'`
set tpx=`echo 72000 $tpx_o | awk '{print $1*sqrt($2*$2)}'| awk '{print int($1+0.5)}'`
set tpy=`echo 72000 $tpy_o | awk '{print $1*sqrt($2*$2)}'| awk '{print int($1+0.5)}'`
echo 'Image pixels are '$tpx' by '$tpy' sky pixels...'

set cpx=`grep CRPIX1 tempdmp3 | head -1 | awk '{print $3}'`
set cpy=`grep CRPIX2 tempdmp3 | head -1 | awk '{print $3}'`
#echo 'Image ref. pixel is at '$cpx','$cpy'...'
set cvx=`grep CRVAL1 tempdmp3 | head -1 | awk '{print $3}'`
set cvy=`grep CRVAL2 tempdmp3 | head -1 | awk '{print $3}'`
#echo 'Image ref. coords are '$cvx','$cvy'...'

set trx=`grep REFXCRPX tempdmp3 | head -1 | awk '{print $2}'`
set try=`grep REFYCRPX tempdmp3 | head -1 | awk '{print $2}'`
#echo 'Image event file sky coordinates centred at X,Y='$trx','$try'...'
#echo 'Image event file central sky coordinates: '$ranom','$decnom
set tdx=`grep REFXCDLT tempdmp3 | head -1 | awk '{print $2}'`
set tdy=`grep REFYCDLT tempdmp3 | head -1 | awk '{print $2}'`

set ty0=`echo $cvy $tpy_o $cpy 0.5 | awk '{print $1 - ($2*($3-$4))}'`
# cos does not work with awk on Solaris
set tx0=`echo $cvx $tpx_o $cpx 0.5 $ty0 $radian | awk '{print $1 - ($2*($3-$4)/cos($5/$6))}'`
#echo 'Image pixel (1,1) is at '$tx0','$ty0'... '

set px0=`echo $trx $ranom $tx0 $ty0 $radian $tdx | awk '{print $1-( ($2-$3)*(cos($4/$5))/$6 )}'`
set py0=`echo $try $decnom $ty0 $tdy | awk '{print $1-( ($2-$3)/$4 )}'`
#echo $px0 $py0
set px0=`echo $px0 | awk '{print int($1+0.5)}'`
set py0=`echo $py0 | awk '{print int($1+0.5)}'`

if (-e tempint2) rm -f tempint2
echo 'Calling interpolation and rebinning routine BGrebinimage2SKY_interp...'
./BGrebinimage2SKY_interp tempevl2 $px0 $py0 $tnx $tny $tpx $tpy $rinterp tempint2

if (-e tempfcr3) rm -f tempfcr3
touch tempfcr3
echo "X J" >> tempfcr3
echo "Y J" >> tempfcr3
echo "VALUE D" >> tempfcr3
echo "" >> tempfcr3

echo 'Checking for NAN dropouts... If there are many, your interpolated image may be invalid...'
if (-e tempnan1) rm -f tempnan1
grep NAN tempint2 > tempnan1
grep NaN tempint2 >> tempnan1   # for Solaris

set nx=`tail -1 tempint2 | awk '{print $1}'`
set ny=`tail -1 tempint2 | awk '{print $2}'`

if (-e tempint3) rm -f tempint3
touch tempint3

set nan=`wc -l tempnan1 | awk '{print $1}'`
echo "No of NANs found : "$nan
if ($nan > 0) then
# On Solaris NaN is not separated from y by a blank. Remove it first
 set nan2=`grep NaN tempnan1 | wc -l`
 if ($nan2 > 0) then
  sed '1,$ s/NaN//' tempnan1 > tempnan2
  mv tempnan2 tempnan1
 endif
 echo " ... attempting to repair..."
 echo $nx $ny $nan
 set n = 0
 set oldrow = 0
 set rowm1 = 0 
 while ($n != $nan)
  @ n = $n + 1
  set x = `head -$n tempnan1 | tail -1 | awk '{print $1}'`
  set y = `head -$n tempnan1 | tail -1 | awk '{print $2}'`
  set row=`echo $x $y $ny | awk '{print (($1 - 1) * $3) + $2}'`
  @ rowm1 = $row - 1
  set rowt = `echo $rowm1 $oldrow | awk '{print $1 - $2}'`
  head -$rowm1 tempint2 | tail -$rowt >> tempint3 
  @ oldrow = $row  

  set nc = 0
  set sumval = 0
  if ($x > 1) then 
    set new = `echo $row $ny | awk '{print $1 - $2}'`
    set val = `head -$new tempint2 | tail -1 | awk '{print $3}'`
    set sumval = `echo $sumval $val | awk '{print $1 + $2}'`  
    @ nc = $nc + 1
  endif
  if ($x < $nx) then 
    set new = `echo $row $ny | awk '{print $1 + $2}'`
    set val = `head -$new tempint2 | tail -1 | awk '{print $3}'`
    set sumval = `echo $sumval $val | awk '{print $1 + $2}'`  
    @ nc = $nc + 1
  endif
  if ($y > 1) then 
    set new = `echo $row | awk '{print $1 - 1}'`
    set val = `head -$new tempint2 | tail -1 | awk '{print $3}'`
    set sumval = `echo $sumval $val | awk '{print $1 + $2}'`  
    @ nc = $nc + 1
  endif
  if ($y < $ny) then 
    set new = `echo $row | awk '{print $1 + 1}'`
    set val = `head -$new tempint2 | tail -1 | awk '{print $3}'`
    set sumval = `echo $sumval $val | awk '{print $1 + $2}'`  
    @ nc = $nc + 1
  endif
  set av=`echo $sumval $nc | awk '{print $1/$2}'`
  echo $n $x $y $av
  echo $x $y $av >> tempint3
  if ($n == $nan) then 
    set row=`echo $x $y $ny | awk '{print (($1 - 1) * $3) + $2}'`
    set left=`echo $nx $ny $row | awk '{print ($1 * $2) - $3}'`
    tail -$left tempint2 >> tempint3  
  endif
 end
else
 if (-e tempint3) rm -f tempint3
 mv tempint2 tempint3
endif

echo 'Creating image from interpolated and rebinned data...'
flst2im tempint3 tempimg2 xrange=1,$tnx yrange=1,$tny rows='-' clobber=yes

if ($scaling > 0) then
  echo 'Counts map... rescaling on pixel size...'
  set scale=`echo $detbin $tpx $tpy | awk '{print $2*$3/$1/$1}'`
  fcarith tempimg2 $scale tempimg3 MUL clobber=yes
else
  echo 'Exposure map... no scaling of image...'
  mv tempimg2 tempimg3
endif

echo 'Copying template image header...'
cphead $templ tempimg3 scale=yes

if (-e $out) then 
  echo 'Output file exists... overwriting...'
  rm -f $out
endif
mv -f tempimg3 $out

rm -f templst2 tempfmd2 tempevl2 tempfcr2 tempint2 tempint3 tempfcr3 tempdmp3 tempimg2 tempdmp2 temptbl2 tempimg3 tempnan1

echo $0' completed.'
exit