1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
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"
|
