#!/bin/sh
# set -x

if [ x"$3" = x ]; then
  echo "usage: $0 ievents iwcs oimage"
  echo "where:"
  echo "  ievents       # input events file"
  echo "  iwcs          # input image file containing new WCS header"
  echo "  oimage        # output image file"
  echo " "
  echo "The x, y position of each event in the input file is transformed to"
  echo "the a new x, y based on the specified WCS and the result is written"
  echo "to the specified output image."
  exit 1
fi

IEVENTS="$1"
IWCS="$2"
OIMAGE="$3"

WEDIT=wcs.edit
WCALC=nwcs.calc
TEVENTS=tevents.fits
TIMAGE=timage.fits

echo "retrieving WCS information from $IWCS ..."
funhead "$IWCS" | egrep "^EPOCH|^RADECSYS|^EQUINOX|^CTYPE1|^CTYPE2|^CRVAL1|^CRVAL2|^CRPIX1|^CRPIX2|^CDELT1|^CDELT2|^CROTA1|^CROTA2|^CUNIT1|^CUNIT2|^CD1_1|^CD1_2|^CD2_1|^CD2_2|^PV1_0|^PV1_1|^PV1_2|^PV1_3|^PV1_4|^PV1_5|^PV1_6|^PV1_7|^PV1_8|^PV1_9|^PV2_0|^PV2_1|^PV2_2|^PV2_3|^PV2_4|^PV2_5|^PV2_6|^PV2_7|^PV2_8|^PV2_9|^LTM1_1|^LTM1_2|^LTM2_1|^LTM2_2|^LTV1|^LTV2|^DTM1_1|^DTM1_2|^DTM2_1|^DTM2_2|^DTV1|^DTV2" > $WEDIT

echo "converting events from $IEVENTS using new WCS ..."
cat <<END > $WCALC
local
double tli2p (double di, double tlmin, double binsiz, int type);
char *wcsfile=NULL;
char *s;
int i;
int init=0;
int got=0;
int offscl;
double tlmin1, tlmin2;
double dval1, dval2;
double nx, ny;
void *wcs1;
void *wcs2;
Fun fun2;
FITSCard card;
end
before
  if( ARGC < 1 ){
    gerror(stderr, "requires -a wcsfile\n");
    exit(1);  
  }
  wcsfile = ARGV(0);
  if( !(fun2 = FunOpen(wcsfile, "r", NULL)) ){
    gerror(stderr, "can't FunOpen wcs file (or find extension): %s\n",
	   wcsfile);
  }
  FunInfoGet(fun2, FUN_WCS, &wcs2, 0);
  FunInfoGet(fun,  FUN_WCS, &wcs1, FUN_MIN1, &tlmin1, FUN_MIN2, &tlmin2, 0);
  xdim=FunParamGeti(fun2, "NAXIS1", 0, -1, &got);
  if( !got ){
    gerror(stderr, "NAXIS1 missing from WCS image\n");
    exit(1);
  }
  ydim=FunParamGeti(fun2, "NAXIS2", 0, -1, &got);
  if( !got ){
    gerror(stderr, "NAXIS2 missing from WCS image\n");
    exit(1);
  }
  pix2wcs(wcs2, (double)0, (double)0, &dval1, &dval2);
  wcs2pix(wcs1, dval1, dval2, &nx, &ny, &offscl);
  nx = tli2p(nx, tlmin1, 0, 'D');
  ny = tli2p(ny, tlmin2, 0, 'D');
  FunParamPuti(ofun, "NWCSX", 1, (int)nx, "lower corner x of wcs image", 1);
  FunParamPuti(ofun, "NWCSY", 1, (int)ny, "lower corner y of wcs image", 1);
  pix2wcs(wcs2, (double)xdim-1, (double)ydim-1, &dval1, &dval2);
  wcs2pix(wcs1, dval1, dval2, &nx, &ny, &offscl);
  nx = (double)tli2p(nx, tlmin1, 0, 'D');
  ny = (double)tli2p(ny, tlmin2, 0, 'D');
  FunParamPuti(ofun, "NWCSX", 2, (int)nx, "upper corner x of wcs image", 1);
  FunParamPuti(ofun, "NWCSY", 2, (int)ny, "upper corner y of wcs image", 1);
end
pix2wcs(wcs1, (double)cur->x, (double)cur->y, &dval1, &dval2);
wcs2pix(wcs2, dval1, dval2, &nx, &ny, &offscl);
/* cur->x = (int)(nx+0.5); */
/* cur->y = (int)(ny+0.5); */
cur->x = nx;
cur->y = ny;
END
funcalc -n -f "$WCALC" -a "$IWCS" "$IEVENTS" "$TEVENTS" > foo.c
funcalc -f "$WCALC" -a "$IWCS" "$IEVENTS" "$TEVENTS"

XDIM=`funhead $IWCS | egrep NAXIS1 | awk '{print $3}'`
YDIM=`funhead $IWCS | egrep NAXIS2 | awk '{print $3}'`
XCEN=`expr $XDIM / 2`
YCEN=`expr $YDIM / 2`
echo "found wcs image dimensions $XCEN@$XDIM, $YCEN@$YDIM ..."

echo "converting new event file to image file $OIMAGE ..."
funimage "$TEVENTS[$XDIM@$XCEN,$YDIM@$YCEN]" "$TIMAGE"

echo "adding new WCS to output image $OIMAGE ..."
funhead "$TIMAGE" "$OIMAGE" "$WEDIT"

echo "determine lower and upper corners of $IWCS image in event coords:"
funhead "$OIMAGE" | egrep "corner . of wcs image"

echo "cleaning up ..."
rm -f "$WEDIT" "$WCALC" "$TEVENTS" "$TIMAGE"