1 files changed, 109 insertions, 0 deletions
diff --git a/funtest/nwcs b/funtest/nwcs
new file mode 100755
index 0000000..cc21ba7
--- /dev/null
+++ b/funtest/nwcs
@@ -0,0 +1,109 @@
+#!/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"
+