#!/bin/tcsh -f

if ($#argv != 2 ) then

 cat <<EOF
 
Usage DOA.csh MECS2.evt MECS3.evt

This program determines the average optical axis position 
of the MECS2 and MECS3 combined in both X,Y and DETX,DETY.  It
gives the distance of the mean in arcmin.

EOF
 exit
endif

set evt1 = $1
set evt2 = $2

set X1 = `fdump ${evt1} STDOUT - - prdata=no page=no | grep 'OPTIC9  ' | awk '{print$3}'`
set Y1 = `fdump ${evt1} STDOUT - - prdata=no  page=no | grep 'OPTIC10 ' | awk '{print$3}'`
set X2 = `fdump ${evt2} STDOUT - - prdata=no page=no  | grep 'OPTIC9  ' | awk '{print$3}'`
set Y2 = `fdump ${evt2} STDOUT - - prdata=no page=no  | grep 'OPTIC10 ' | awk '{print$3}'`
set MX = `echo "scale=10;(${X1}+${X2})/2" | bc -l`
set MY = `echo "scale=10;(${Y1}+${Y2})/2" | bc -l`

# Since the Detector coordinates of the MECS2 is mirrored L-R
# A shift left on the sky is a shift right on the Detector
set DX1 = `echo "scale=10;(${X1}-${MX})" | bc -l`
set DY1 = `echo "scale=10;(${MY}-${Y1})" | bc -l`
set DX2 = `echo "scale=10;(${MX}-${X2})" | bc -l`

# Since the Detector coordinates of the MECS3 is mirrored U-D
# A shift left on the sky is a shift right on the Detector
set DY2 = `echo "scale=10;(${Y2}-${MY})" | bc -l`
set DETX1 = `echo "scale=10;(256+(${DX1}))" | bc -l`
set DETY1 = `echo "scale=10;(256+(${DY1}))" | bc -l`
set DETX2 = `echo "scale=10;(256+(${DX2}))" | bc -l`
set DETY2 = `echo "scale=10;(256+(${DY2}))" | bc -l`


set DEC = `fdump ${evt1} STDOUT - - prdata=no  page=no | grep 'DEC_PNT ' | awk '{printf"%10.10f\n",$3}'`
set DPP = `fdump ${evt1} STDOUT - - prdata=no  page=no | grep 'TCDLT10 ' | awk '{print$3}'`
set SPP = `echo "scale=10;${DPP}*3600" | bc -l`
set RASPP = `echo "scale=10;${SPP}/c(${DEC}*pi/180)" | bc -l`
#set RASPP = ${SPP}
set PDIST = `echo "scale=10;sqrt((${DX1}*${RASPP})^2+(${DY1}*${SPP})^2)" | bc -l`

# NOW DETERMINE THE RA AND DEC OF THE NEW CENTER
set RA  = `fdump ${evt1} STDOUT - - prdata=no page=no | grep 'TCRVL9  ' | awk '{printf"%10.10f\n",$3}'`
set DEC = `fdump ${evt1} STDOUT - - prdata=no page=no | grep 'TCRVL10 ' | awk '{printf"%10.10f\n",$3}'`
set X0  = `fdump ${evt1} STDOUT - - prdata=no page=no | grep 'TCRPX9  ' | awk '{printf"%10.10f\n",$3}'`
set Y0  = `fdump ${evt1} STDOUT - - prdata=no page=no | grep 'TCRPX10 ' | awk '{printf"%10.10f\n",$3}'`
set deltaX = `echo "scale=10;${X0} - ${MX}" | bc -l`
set deltaY = `echo "scale=10;${Y0} - ${MY}" | bc -l`
set deltad = `fdump ${evt1} STDOUT - - prdata=no page=no | grep 'TCDLT9  ' | awk '{printf"%10.10f\n",$3}'`
set deltar = `fdump ${evt1} STDOUT - - prdata=no page=no | grep 'TCDLT10 ' | awk '{printf"%10.10f\n",$3}'`
set deltar = `echo "scale=10;${deltar}/c(${DEC}*pi/180)" | bc -l`
set RA0 = `echo "scale=10;($RA + (${deltar})*(${deltaX}))/15" | bc -l`
set HRS = `echo "scale=0;${RA0}/1" | bc -l`
set MIN = `echo "scale=0;((${RA0}-${HRS})*60)/1" | bc -l`
set SEC = `echo "scale=4;((((${RA0}-${HRS})*60)-${MIN})*60)/1" | bc -l`
set DEC0 = `echo "scale=10;$DEC + (${deltad})*(${deltaY})" | bc -l`
set DEG = `echo "scale=0;${DEC0}/1" | bc -l`
if (${DEG} < 0 ) then
 set DMN = `echo "scale=0;((${DEC0}*(-1)+(${DEG}))*60)/1" | bc -l`
 set DSC = `echo "scale=4;((((${DEC0}*(-1)+(${DEG}))*60)-${DMN})*60)/1" | bc -l`
else
 set DMN = `echo "scale=0;((${DEC0}-(${DEG}))*60)/1" | bc -l`
 set DSC = `echo "scale=4;((((${DEC0}-${DEG})*60)-${DMN})*60)/1" | bc -l`
endif
set col = ':'
set RAOUT = `echo "scale=10;${RA0}*15" | bc -l`
cat <<EOF

 The mean center in sky coordinates:
 ${MX} ${MY}

 For ${evt1} this corresponds to:
 ${DETX1} ${DETY1}
 
 For ${evt2} this corresponds to:
 ${DETX2} ${DETY2}
 
 The mean is:
 ${PDIST} arcsecs
 from the optical axes of each event file.
 
 The RA and DEC of this point is:
 ${RAOUT},${DEC0}
 ${HRS}${col}${MIN}${col}${SEC}  ${DEG}${col}${DMN}${col}${DSC}

EOF
 
 
