Index: branches/eam_branches/20090522/extsrc/gpcsw/Makefile.in =================================================================== --- branches/eam_branches/20090522/extsrc/gpcsw/Makefile.in (revision 24238) +++ branches/eam_branches/20090522/extsrc/gpcsw/Makefile.in (revision 24557) @@ -2,5 +2,5 @@ default: all -all: burntool +all: burntool tiltystreak install: all @@ -24,4 +24,12 @@ svn export https://svn.ifa.hawaii.edu/gpc/repos/sw/trunk/fits/libfhreg gpcsrc/fits/libfhreg.new svn export https://svn.ifa.hawaii.edu/gpc/repos/sw/trunk/fits/burntool gpcsrc/fits/burntool.new + svn export https://svn.ifa.hawaii.edu/gpc/repos/sw/trunk/fits/tiltystreak gpcsrc/fits/tiltystreak.new +# ensure the target directory exists + mkdir -p gpcsrc/analysis/libpscoords + mkdir -p gpcsrc/analysis/libpsf + mkdir -p gpcsrc/fits/libfh + mkdir -p gpcsrc/fits/libfhreg + mkdir -p gpcsrc/fits/burntool + mkdir -p gpcsrc/fits/tiltystreak # replace new files with those from the imported directories cp -rf gpcsrc/analysis/libpscoords.new/* gpcsrc/analysis/libpscoords @@ -30,4 +38,5 @@ cp -rf gpcsrc/fits/libfhreg.new/* gpcsrc/fits/libfhreg cp -rf gpcsrc/fits/burntool.new/* gpcsrc/fits/burntool + cp -rf gpcsrc/fits/tiltystreak.new/* gpcsrc/fits/tiltystreak # remove the temp directory rm -rf gpcsrc/analysis/libpscoords.new @@ -36,4 +45,5 @@ rm -rf gpcsrc/fits/libfhreg.new rm -rf gpcsrc/fits/burntool.new + rm -rf gpcsrc/fits/tiltystreak.new # modify burntool to avoid libpsf (which requires fortran) cp -f burntool.nopsf/Makefile gpcsrc/fits/burntool @@ -48,4 +58,5 @@ ln -sf ../Make.Common gpcsrc/fits/libfhreg/Make.Common ln -sf ../Make.Common gpcsrc/fits/burntool/Make.Common + ln -sf ../Make.Common gpcsrc/fits/tiltystreak/Make.Common burntool: @@ -54,2 +65,7 @@ cd gpcsrc && make -C analysis/libpscoords install cd gpcsrc && make -C fits/burntool install + +tiltystreak: + cd gpcsrc && make -C fits/libfh install + cd gpcsrc && make -C fits/libfhreg install + cd gpcsrc && make -C fits/tiltystreak install Index: branches/eam_branches/20090522/extsrc/gpcsw/burntool.nopsf/psfstats.c =================================================================== --- branches/eam_branches/20090522/extsrc/gpcsw/burntool.nopsf/psfstats.c (revision 24238) +++ branches/eam_branches/20090522/extsrc/gpcsw/burntool.nopsf/psfstats.c (revision 24557) @@ -10,10 +10,9 @@ #include "burntool.h" #include "math.h" -#include "psf/psf.h" /****************************************************************/ /* psf_stats(): Tell us about this PSF star */ -STATIC int psf_stats(int nx, int ny, IMTYPE *data, int bias, - double *fwhm, double *q) +STATIC int psf_stats(int nx, int ny, IMTYPE *data, int bias, + double *fwhm, double *q) { return(1); Index: branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/analysis/libpscoords/man/libpscoords.1 =================================================================== --- branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/analysis/libpscoords/man/libpscoords.1 (revision 24557) +++ branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/analysis/libpscoords/man/libpscoords.1 (revision 24557) @@ -0,0 +1,136 @@ +.nf +NAME: + libpscoords - convert all the various Pan-STARRS coordinate systems + +API: + pscoords [args] < infile > outfile + +DESCRIPTION: + +There are five coordinate systems used by Pan-STARRS: + + * Sky coordinates of conventional RA, Dec, and PA; + * Tangent plane coordinates (meaning the angles at which rays enter the + optical system, relative to the optical system itself); + * Focal plane coordinates (meaning physical position of an image on the + focal plane); + * OTA coordinates (meaning the pixel position of an image on the silicon + of a given OTA); and + * Cell coordinates (meaning the pixel position of an image on a given + cell). + +libpscoords has routines to convert between these various coordinate +systems. + +Normal Application +------------------ + +To go from a position (a,d) relative to (a0,d0,pa) pointing to cell coords: +--------------------------------------------------------------------------- + psc_tproject(a, d, a0, d0, pa, density, pi/2-alt, vertical, &x, &y); + psc_psoptics2(x/sec2rad, y/sec2rad, dx, dy, dpa, pscale, d2, d3, &xfp, &yfp); + psc_fp_to_pixel(xfp, yfp, &ota_xid, &ota_yid, &xota, &yota); + psc_pixel_to_cell(xota, yota, &cell_xid, &cell_yid, &xcell, &ycell); + +(skipping the d2 distortion term one can also use) + psc_psoptics(x/sec2rad, y/sec2rad, dx, dy, dpa, pscale, d3, &xfp, &yfp); + +To go from cell coords to a position (a,d) relative to (a0,d0,pa) pointing: +--------------------------------------------------------------------------- + psc_cell_to_pixel(cell_xid,cell_yid, xcell,ycell, &xota,&yota); + psc_pixel_to_fp(ota_xid, ota_yid, xota, yota, &xfp, &yfp); + psc_invoptics2(xfp, yfp, dx, dy, dpa, pscale, d2, d3, &x, &y); + psc_tplonglat(x*sec2rad, y*sec2rad, a0, d0, pa, density, pi/2-alt, vertical, &a, &d); + +(skipping the d2 distortion term one can also use) + psc_invoptics(xfp, yfp, dx, dy, dpa, pscale, d3, &x, &y); + +See the prototypes for suggested parameter values. Note that for +Pan-STARRS the "vertical" argument is the ROT angle (in rad). + + +Prototypes +---------- + +/* Enable application of chip offsets? */ +int psc_do_chipoff(int doit); /* Apply chip offsets? (0/1, default=0) */ + +/* Convert a focal plane position to OTA otax,otay, Cell cellx,celly, Pixel */ +int psc_fp_to_pixel(double xfp_um, double yfp_um, /* FP position (um) */ + int *ota_xid, int *ota_yid, /* OTA ID (0:7) */ + double *xota_pix, double *yota_pix); /* OTA position (pix) */ + +/* Convert a pixel position in OTA ota_xid,ota_yid to FP */ +int psc_pixel_to_fp(int ota_xid, int ota_yid, /* OTA ID (0:7) */ + double xota_pix, double yota_pix, /* OTA position (pix) */ + double *xfp_um, double *yfp_um); /* FP position (um) */ + +/* Convert an OTA pixel position to cell ID and cell coords */ +int psc_pixel_to_cell(double xota_pix, double yota_pix, /* OTA position (pix) */ + int *cell_xid, int *cell_yid, /* Cell ID (0:7) or -1 */ + double *xcell_pix, double *ycell_pix);/* Cell position (pix) */ + +/* Convert a cell position and ID to OTA coords */ +int psc_cell_to_pixel(int cell_xid, int cell_yid, /* Cell ID (0:7) */ + double xcell_pix, double ycell_pix, /* Cell position (pix) */ + double *xota_pix, double *yota_pix); /* OTA position (pix) */ + +/* Convert tangent plane (x,y) at (a0,d0,pa) to RA,Dec (a,d) [radians] */ +/* y is rotated CCW by pa from N; NOTE: x is west when PA = 0 (pos parity) */ +/* Set density=0 to apply no correction for differential refraction, otherwise P/T in STP */ +int psc_tplonglat(double x_rad, double y_rad, double a0_rad, double d0_rad, double pa_rad, + double density/*0.71 for HK*/, double zenith_rad/*pi/2-ALT*/, + double vertical_rad/*ROT*/, double *a_rad, double *d_rad); + +/* Project RA,Dec (a,d) to the tangent plane (x,y) at (a0,d0,pa) [radians] */ +/* y is rotated CCW by pa from N; NOTE: x is west when PA = 0 (pos parity) */ +/* Set density=0 to apply no correction for differential refraction, otherwise P/T in STP */ +int psc_tproject(double a_rad, double d_rad, double a0_rad, double d0_rad, double pa_rad, + double density/*0.71 for HK*/, double zenith_rad/*pi/2-ALT*/, + double vertical_rad/*ROT*/, double *x_rad, double *y_rad); + +/* Evaluate an inverse optical model: microns -> arcsec notionally */ +/* Offset, scale with *small* distortion, rotate */ +int psc_invoptics(double xfp_um, double yfp_um, double dx_sec, double dy_sec, double dpa_rad, + double pscale/*38.86um/sec*/, double d3/*1.49e-10sec^-2*/, + double *x_sec, double *y_sec); + +/* Evaluate an optical model: arcsec -> microns notionally */ +/* Offset, scale with distortion, rotate */ +int psc_psoptics(double x_sec, double y_sec, double dx_sec, double dy_sec, double dpa_rad, + double pscale/*38.86um/sec*/, double d3/*1.49e-10sec^-2*/, + double *xfp_um, double *yfp_um); + +/* Evaluate an inverse optical model: microns -> arcsec notionally */ +/* Offset, scale with *small* distortion, rotate */ +int psc_invoptics2(double xfp_um, double yfp_um, double dx_sec, double dy_sec, double dpa_rad, + double pscale/*38.793um/sec*/, + double d2/*9.8e-7sec^-2*/, double d3/*3.1e-11sec^-2*/, + double *x_sec, double *y_sec); + +/* Evaluate an optical model: arcsec -> microns notionally */ +/* Offset, scale with distortion, rotate */ +int psc_psoptics2(double x_sec, double y_sec, double dx_sec, double dy_sec, double dpa_rad, + double pscale/*38.793um/sec*/, + double d2/*9.8e-7sec^-2*/, double d3/*3.1e-11sec^-2*/, + double *xfp_um, double *yfp_um); + +/* Enable application of chip offsets (um, mrad)? */ +int psc_load_otaoff(const char *fname); + + + + + + +BUGS: + 090516: Conceivably a tiny bug in Sky<->TP projections very near + the pole, equally conceivably a bug in libm. + +SEE ALSO: + pscoords + +AUTHORS: + John Tonry + +VERSION: Index: branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/analysis/libpscoords/pscoords.c =================================================================== --- branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/analysis/libpscoords/pscoords.c (revision 24238) +++ branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/analysis/libpscoords/pscoords.c (revision 24557) @@ -24,9 +24,47 @@ rcsid[] = "$Id: pscoords.c,v 1.2 2008/09/06 03:26:48 jt Exp jt $"; +#define MRAD2 (0.001) /* Convert from mrad to rad, used for PSC_OFFROT_T */ + /* Overall focal plane offset and rotation (um, rad) */ -const PSC_OFFROT_T psc_psc_fpoff={0.,0.,0.}; +PSC_OFFROT_T psc_psc_fpoff={0.,0.,0.}; /* Offsets and rotations of each OTA (um, millirad) */ -const PSC_OFFROT_T psc_otaoff[PSC_NX*PSC_NY]={ /* Offsets and rotations derived o4885 090228 */ +PSC_OFFROT_T psc_otaoff[PSC_NX*PSC_NY]={ /* Offsets and rotations derived o4885 090516 */ + { 0.0, 0.0, 0.00}, { -6.8, -32.4, -1.30}, /* OTA00 , OTA10 */ + { 15.7,-130.7, 0.07}, { -9.4, -43.5, -1.34}, /* OTA20 , OTA30 */ + { 17.8, 71.7, -0.48}, { 7.2, 41.3, -0.82}, /* OTA40 , OTA50 */ + { 22.1, -48.3, -1.74}, { 0.0, 0.0, 0.00}, /* OTA60 , OTA70 */ + { -23.0, -12.9, 0.68}, { -52.0, -52.2, 0.88}, /* OTA01 , OTA11 */ + { 28.8, -10.9, -0.73}, { -22.9, 18.1, -0.89}, /* OTA21 , OTA31 */ + { -13.8, 45.9, 0.90}, { -12.2, 49.7, -0.32}, /* OTA41 , OTA51 */ + { -31.9, -9.8, -1.60}, { 1.9, 8.0, 0.53}, /* OTA61 , OTA71 */ + { -27.2, 42.6, 2.73}, { 6.5, -52.6, 0.48}, /* OTA02 , OTA12 */ + { 10.1, -37.9, -1.53}, { 3.1,-139.0, 1.75}, /* OTA22 , OTA32 */ + { -21.4, 25.6, -1.12}, { -13.8, 0.9, 0.02}, /* OTA42 , OTA52 */ + { -1.8, -16.1, -2.06}, { 76.0, 31.7, 1.70}, /* OTA62 , OTA72 */ + { -15.7, 24.1, 0.14}, {-104.2, -76.8, 1.97}, /* OTA03 , OTA13 */ + { 75.4, 31.8, 1.04}, { -18.7, 14.6, 0.51}, /* OTA23 , OTA33 */ + { -30.6, 53.4, -0.78}, { 12.6, 24.4, 0.43}, /* OTA43 , OTA53 */ + { 8.0, -34.9, 1.44}, { 6.4, -73.2, -0.93}, /* OTA63 , OTA73 */ + { -31.2, 11.2, 0.23}, { -10.3, -2.1, 0.14}, /* OTA04 , OTA14 */ + { 33.2, 207.9, -0.40}, { -75.7,-105.9, 0.93}, /* OTA24 , OTA34 */ + { 8.1, -15.9, -0.03}, { -16.1, -13.3, -0.06}, /* OTA44 , OTA54 */ + { 4.1, 47.6, 1.26}, { 71.6, -6.0, 0.42}, /* OTA64 , OTA74 */ + { -23.4, 46.7, 2.65}, { 31.1, 5.5, -0.64}, /* OTA05 , OTA15 */ + { -43.2, -6.5, 1.68}, { -26.0, -53.2, 0.00}, /* OTA25 , OTA35 */ + { -43.0, 9.1, -0.69}, { 30.1, 69.7, 0.62}, /* OTA45 , OTA55 */ + { -48.5, 61.6, -0.51}, { 65.5, -30.1, 0.69}, /* OTA65 , OTA75 */ + { 2.5, 57.2, 0.17}, { 0.6, 25.4, -0.91}, /* OTA06 , OTA16 */ + { -70.8, -17.4, 1.36}, { 13.6, -24.2, 0.53}, /* OTA26 , OTA36 */ + { -6.8, 28.1, -1.35}, { 20.9, 42.0, 0.19}, /* OTA46 , OTA56 */ + { 6.2, 13.2, -1.50}, { 36.5, 48.8, -0.36}, /* OTA66 , OTA76 */ + { 0.0, 0.0, 0.00}, { 20.4, 1.7, 0.09}, /* OTA07 , OTA17 */ + { 83.8, 11.8, -0.56}, { -26.5, -21.2, 0.56}, /* OTA27 , OTA37 */ + { -70.5, 7.2, -1.05}, { -41.2, -3.6, -1.89}, /* OTA47 , OTA57 */ + { -35.0, 49.5, -0.40}, { 0.0, 0.0, 0.00}}; /* OTA67 , OTA77 */ + +/* Offsets and rotations of each OTA (um, millirad) */ +#ifdef RUN_THREE_V1 +PSC_OFFROT_T psc_otaoff[PSC_NX*PSC_NY]={ /* Offsets and rotations derived o4885 090228 */ { 0.0, 0.0, 0.00}, { 20.2, -59.5, -1.23}, /* OTA00 , OTA10 */ { 28.1,-138.1, 0.18}, { 2.0, -43.6, -1.22}, /* OTA20 , OTA30 */ @@ -61,4 +99,5 @@ { -74.8, 0.3, -0.92}, { -44.6, -1.6, -1.78}, /* OTA47 , OTA57 */ { -50.3, 72.1, -0.47}, { 0.0, 0.0, 0.00}}; /* OTA67 , OTA77 */ +#endif #ifdef RUN_TWO @@ -136,5 +175,5 @@ #endif -static int DOOFFSET; /* Do these offsets or no? */ +static int DOOFFSET=1; /* Do these offsets or no? */ #if 0 @@ -143,5 +182,29 @@ #endif -/* Enable application of chip offsets? */ +/* Enable application of chip offsets (um, mrad)? */ +int psc_load_otaoff(const char *fname) +{ + int i; + char rbuf[1024]; + FILE *fp; + if( (fp=fopen(fname, "r")) == NULL) { + fprintf(stderr, "error: Error opening OTA offset file '%s'\n", fname); + return(-1); + } + for(i=0; i microns notionally */ +/* Offset, scale with quadratic and cubic distortion, rotate */ +int psc_psoptics2(double x, double y, double dx, double dy, double dpa, + double pscale, double d2, double d3, double *xfp, double *yfp) +{ + double w2=((x-dx)*(x-dx)+(y-dy)*(y-dy)), w; + w = sqrt(w2); + x = pscale * (1.0 + d2*w + d3*w2) * (x-dx); + y = pscale * (1.0 + d2*w + d3*w2) * (y-dy); +/* Apply any extra TP-FP rotation */ + *xfp = x * cos(dpa) + y * sin(dpa); + *yfp = -x * sin(dpa) + y * cos(dpa); + return(0); +} + +/* Evaluate an inverse optical model: microns -> arcsec notionally */ +/* Offset, scale with *small* quadratic and cubic distortion, rotate */ +int psc_invoptics2(double xfp, double yfp, double dx, double dy, double dpa, + double pscale, double d2, double d3, double *x, double *y) +{ + double xp, yp, w; +/* Undo any extra TP-FP rotation */ + xp = (xfp * cos(dpa) - yfp * sin(dpa)) / pscale; + yp = (xfp * sin(dpa) + yfp * cos(dpa)) / pscale; + w = sqrt(xp*xp+yp*yp); + *x = xp / (1.0 + d2*w + d3*(xp*xp+yp*yp)) + dx; + *y = yp / (1.0 + d2*w + d3*(xp*xp+yp*yp)) + dy; + return(0); +} + + +/* Evaluate an optical model: arcsec -> microns notionally */ /* Offset, scale with distortion, rotate */ int psc_psoptics(double x, double y, double dx, double dy, double dpa, @@ -318,5 +412,4 @@ return(0); } - /* Project RA,Dec (a,d) to the tangent plane (x,y) at (a0,d0,pa) [radians] */ Index: branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/analysis/libpscoords/pscoords.h =================================================================== --- branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/analysis/libpscoords/pscoords.h (revision 24238) +++ branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/analysis/libpscoords/pscoords.h (revision 24557) @@ -2,21 +2,28 @@ #define NINT(x) (x<0?(int)((x)-0.5):(int)((x)+0.5)) +#define OTA_CTE (5.0e-6) /* Coefficient of thermal expansion of Si+Mo */ +#define CFFP_CTE (0.5e-6) /* Coefficient of thermal expansion of CFFP */ +#define DELTA_T (-100.0) /* Operating temp minus metrology temp */ + +#define OTA_SCALE (1.0+OTA_CTE*DELTA_T) /* Room temp um to operating um of Si+Mo */ +#define CFFP_SCALE ((1.0+CFFP_CTE*DELTA_T)*25400) /* Room temp inches to um on CFFP */ + /* Silicon layout of OTAs (um) */ -#define PSC_PIXEL 10.0 /* Pixel size of an OTA */ -#define PSC_HCELL 5900.0 /* Horizontal cell size */ -#define PSC_VCELL 5980.0 /* Vertical cell size */ -#define PSC_HSTREET 120.0 /* Horizontal street between cells (HS) */ -#define PSC_VSTREET 180.0 /* Vertical street between cells (VS) */ -#define PSC_LBORDER 442.0 /* Left border (to left-most cell) (LB) */ -#define PSC_RBORDER 418.0 /* Right border (to right-most vertical street) (RB) */ -#define PSC_TBORDER 354.5 /* Top border (to top-most cell) (TB) */ -#define PSC_BBORDER 975.5 /* Bottom border (to lowest horizontal street) (BB) */ -#define PSC_HDIE 49500.0 /* Horizontal size of silicon die */ -#define PSC_VDIE 50130.0 /* Vertical size of silicon die */ +#define PSC_PIXEL (10.0*OTA_SCALE) /* Pixel size of an OTA */ +#define PSC_HCELL (5900.0*OTA_SCALE) /* Horizontal cell size */ +#define PSC_VCELL (5980.0*OTA_SCALE) /* Vertical cell size */ +#define PSC_HSTREET (120.0*OTA_SCALE) /* Horizontal street between cells (HS) */ +#define PSC_VSTREET (180.0*OTA_SCALE) /* Vertical street between cells (VS) */ +#define PSC_LBORDER (442.0*OTA_SCALE) /* Left border (to left-most cell) (LB) */ +#define PSC_RBORDER (418.0*OTA_SCALE) /* Right border (to right-most vertical street) (RB) */ +#define PSC_TBORDER (354.5*OTA_SCALE) /* Top border (to top-most cell) (TB) */ +#define PSC_BBORDER (975.5*OTA_SCALE) /* Bottom border (to lowest horizontal street) (BB) */ +#define PSC_HDIE (49500.0*OTA_SCALE) /* Horizontal size of silicon die */ +#define PSC_VDIE (50130.0*OTA_SCALE) /* Vertical size of silicon die */ /* Mechanical layout of OTAs (um) */ -#define PSC_HMECH (1.957*25400) /* Horizontal mech spacing of OTA placement */ -#define PSC_VMECH (2.025*25400) /* Vertical mech spacing of OTA placement */ -#define PSC_VMOFF (PSC_VDIE-(1.568+0.375)*25400) /* Vert die offset between sides */ +#define PSC_HMECH (1.957*CFFP_SCALE) /* Horizontal mech spacing of OTA placement */ +#define PSC_VMECH (2.025*CFFP_SCALE) /* Vertical mech spacing of OTA placement */ +#define PSC_VMOFF (PSC_VDIE-(1.568+0.375)*CFFP_SCALE) /* Vert die offset between sides */ #define PSC_REFRACT_CONST 60.0 /* Standard refraction ("/tanz) at STP */ @@ -123,11 +130,22 @@ /* Overall focal plane offset and rotation (um, rad) */ -const PSC_OFFROT_T psc_fpoff; +PSC_OFFROT_T psc_fpoff; /* Offsets and rotations of each OTA (um, millirad) */ -const PSC_OFFROT_T psc_otaoff[PSC_NX*PSC_NY]; - +PSC_OFFROT_T psc_otaoff[PSC_NX*PSC_NY]; + + +/* This is from Run 3, 090516 */ +#define PS_scale 38.7932 /* Default PS plate scale [um/arcsec] */ +#define PS_d2 9.78e-7 /* Default quadratic PS distortion [arcsec^-1] */ +#define PS_d3 3.16e-11 /* Default cubic PS distortion [arcsec^-2] */ + + +#ifdef RUN_THREE_V1 /* This is from Run 3, 090228 */ #define PS_scale 38.860 /* Default PS plate scale [um/arcsec] */ -#define PS_d3 1.49e-10 /* Default PS distortion [arcsec^-2] */ +#define PS_d3 1.49e-10 /* Default cubic PS distortion [arcsec^-2] */ +#endif + + #define PS_airdens 0.71 /* Default PS1 qir density (Haleakala) */ @@ -139,7 +157,11 @@ * --------------------------------------------------------------------------- * psc_tproject(a, d, a0, d0, pa, density, pi/2-alt, vertical, &x, &y); - * psc_psoptics(x/sec2rad, y/sec2rad, dx, dy, dpa, pscale, d3, &xfp, &yfp); + * psc_psoptics2(x/sec2rad, y/sec2rad, dx, dy, dpa, pscale, d2, d3, &xfp, &yfp); * psc_fp_to_pixel(xfp, yfp, &ota_xid, &ota_yid, &xota, &yota); * psc_pixel_to_cell(xota, yota, &cell_xid, &cell_yid, &xcell, &ycell); + * + * can also use + * psc_psoptics(x/sec2rad, y/sec2rad, dx, dy, dpa, pscale, d3, &xfp, &yfp); + * * * To go from cell coords to a position (a,d) relative to (a0,d0,pa) pointing: @@ -147,6 +169,9 @@ * psc_cell_to_pixel(cell_xid,cell_yid, xcell,ycell, &xota,&yota); * psc_pixel_to_fp(ota_xid, ota_yid, xota, yota, &xfp, &yfp); + * psc_invoptics2(xfp, yfp, dx, dy, dpa, pscale, d2, d3, &x, &y); + * psc_tplonglat(x*sec2rad, y*sec2rad, a0, d0, pa, density, pi/2-alt, vertical, &a, &d); + * + * can also use * psc_invoptics(xfp, yfp, dx, dy, dpa, pscale, d3, &x, &y); - * psc_tplonglat(x*sec2rad, y*sec2rad, a0, d0, pa, density, pi/2-alt, vertical, &a, &d); */ @@ -155,5 +180,5 @@ /**************/ /* Enable application of chip offsets? */ -int psc_do_chipoff(int doit); /* Apply chip offsets? (0/1, default=0) */ +int psc_do_chipoff(int doit); /* Apply chip offsets? (0/1, default=1) */ /* Convert a focal plane position to OTA otax,otay, Cell cellx,celly, Pixel */ @@ -202,2 +227,19 @@ double pscale/*38.86um/sec*/, double d3/*1.49e-10sec^-2*/, double *xfp_um, double *yfp_um); + +/* Evaluate an inverse optical model: microns -> arcsec notionally */ +/* Offset, scale with *small* distortion, rotate */ +int psc_invoptics2(double xfp_um, double yfp_um, double dx_sec, double dy_sec, double dpa_rad, + double pscale/*38.793um/sec*/, + double d2/*9.8e-7sec^-2*/, double d3/*3.1e-11sec^-2*/, + double *x_sec, double *y_sec); + +/* Evaluate an optical model: arcsec -> microns notionally */ +/* Offset, scale with distortion, rotate */ +int psc_psoptics2(double x_sec, double y_sec, double dx_sec, double dy_sec, double dpa_rad, + double pscale/*38.793um/sec*/, + double d2/*9.8e-7sec^-2*/, double d3/*3.1e-11sec^-2*/, + double *xfp_um, double *yfp_um); + +/* Enable application of chip offsets (um, mrad)? */ +int psc_load_otaoff(const char *fname); Index: branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/burntool/burnfix.c =================================================================== --- branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/burntool/burnfix.c (revision 24238) +++ branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/burntool/burnfix.c (revision 24557) @@ -105,4 +105,19 @@ if(!cell->persist[k].fiterr) sub_fit(nx, ny, NX, buf, &(cell->persist[k]), -1); + } + return(0); +} + + +/****************************************************************/ +/* burn_apply(): Subtract the trail fits from the image */ +STATIC int burn_apply(int nx, int ny, int NX, IMTYPE *buf, CELL *cell) +{ + int k; + +/* Restore all the burns */ + for(k=0; knpersist; k++) { + if(!cell->persist[k].fiterr) + sub_fit(nx, ny, NX, buf, &(cell->persist[k]), +1); } return(0); Index: branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.c =================================================================== --- branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.c (revision 24238) +++ branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.c (revision 24557) @@ -33,5 +33,5 @@ int otanum; int nextend, cellxy, cell, cellcode; - int ext, update, restore, psfsize, psfavg; + int ext, update, restore, apply, tableonly, psfsize, psfavg; IMTYPE *buf; const char *burnfile=NULL, *persistfile=NULL, *persistfitsfile=NULL; @@ -102,6 +102,8 @@ /* Parse the args */ cellxy = -1; - update = 1; - restore = 0; + update = 1; /* Calc fits, apply fits, write img and table */ + restore = 0; /* Restore previous fit only, write img */ + apply = 0; /* Apply previous fit only, write img */ + tableonly = 0; /* Calc fits (apply fits), write table only */ psfsize = 32; psfavg = 0; @@ -137,7 +139,15 @@ update = argv[i][7] == 'y' || argv[i][7] == '1' || argv[i][7] == 't'; -/* Modify the input MEF by subtracting fits? */ +/* Modify the input MEF by adding back fits? */ } else if(strncmp(argv[i], "restore=", 8) == 0) { /* restore={t|f} */ restore = argv[i][8] == 'y' || argv[i][8] == '1' || argv[i][8] == 't'; + +/* Modify the input MEF by subtracting previously calculated fits? */ + } else if(strncmp(argv[i], "apply=", 6) == 0) { /* apply={t|f} */ + apply = argv[i][6] == 'y' || argv[i][6] == '1' || argv[i][6] == 't'; + +/* Calculate and write tables only? */ + } else if(strncmp(argv[i], "tableonly=", 10) == 0) {/* tableonly={t|f} */ + tableonly = argv[i][10] == 'y' || argv[i][10] == '1' || argv[i][10] == 't'; /* Output file for burn streaks */ @@ -150,9 +160,9 @@ /* Input text file for previous burn persistence streaks */ - } else if(strncmp(argv[i], "infits=", 8) == 0) { /* in=fname */ + } else if(strncmp(argv[i], "trailin=", 8) == 0) {/* trailin=fname */ persistfile = argv[i] + 8; /* Same thing, but information is stored in tables in a FITS file. */ - } else if(strncmp(argv[i], "trailinfits=", 8) == 0) { /* infits=fname */ + } else if(strncmp(argv[i], "trailinfits=", 12) == 0) { /* trailinfits=fname */ persistfitsfile = argv[i] + 12; @@ -162,6 +172,6 @@ /* Same thing, but information is stored in tables in a FITS file. */ - } else if(strncmp(argv[i], "trailinfits=", 8) == 0) { /* trailin=fname */ - persistfitsfile = argv[i] + 12; + } else if(strncmp(argv[i], "infits=", 7) == 0) { /* infits=fname */ + persistfitsfile = argv[i] + 7; /* Output file for PSF gallery */ @@ -322,6 +332,6 @@ /* If there is no other persistence info supplied, try getting * it from the input FITS file. */ - if(restore && persistfile == NULL) { - if(persistfitsfile == NULL) persistfitsfile = ifilename; + if((restore || apply) && persistfile == NULL) { + if(persistfitsfile == NULL) persistfitsfile = ifilename; } @@ -334,9 +344,9 @@ /* Read the persistence data for this OTA */ - if(persistfile != NULL) { + if(persistfile != NULL) { /* Text data file */ if(persist_read(OTA, persistfile)) exit(EXIT_FAILURE); - } - else if(persistfitsfile != NULL) { - if(persist_fits_read(OTA, persistfitsfile) != FH_SUCCESS) exit(EXIT_FAILURE); + } else if(persistfitsfile != NULL) { /* FITS table */ + if(persist_fits_read(OTA, persistfitsfile) != FH_SUCCESS) + exit(EXIT_FAILURE); } @@ -369,5 +379,5 @@ "warning: Restoring old burns, but header indicates no burns previously corrected.\n"); } - else if (update && (burn_applied == FH_TRUE)) { + else if ((update||apply) && (burn_applied == FH_TRUE)) { fprintf(stderr, "warning: Applying burn correction, but header indicates burns previously corrected.\n"); @@ -491,5 +501,5 @@ if(VERBOSE & VERB_NORM) { - printf("nx=%d ny=%d prex=%d postx=%d posty=%d BZERO=%.1f\n", + printf("nx=%d ny=%d prex=%d postx=%d posty=%d BZERO=%d\n", naxis1, naxis2, prescan1, ovrscan1, ovrscan2, BZERO); } @@ -542,7 +552,17 @@ } - - if(!restore) { - + if(apply) { +/* Use the table-driven fits instead of calculating new ones */ + burn_apply(naxis1-ovrscan1, naxis2-ovrscan2, naxis1, + buf, OTA+cell); +/* Tell us about it? */ + if(VERBOSE & VERB_NORM) burn_blab(OTA+cell); + + } else if(restore) { +/* Restore the old burns */ + burn_restore(naxis1-ovrscan1, naxis2-ovrscan2, naxis1, + buf, OTA+cell); + + } else { /* Fix up the burns */ burn_fix(naxis1-ovrscan1, naxis2-ovrscan2, naxis1, naxis2, buf, @@ -559,19 +579,15 @@ if(VERBOSE & VERB_NORM) burn_blab(OTA+cell); -/* Fix up the streaks */ - persist_fix(naxis1-ovrscan1, naxis2-ovrscan2, naxis1, buf, - OTA+cell); - +/* Fix up the streaks (don't bother if table only) */ + if(!tableonly) { + persist_fix(naxis1-ovrscan1, naxis2-ovrscan2, naxis1, buf, + OTA+cell); /* Tell us about it? */ - if(VERBOSE & VERB_NORM) persist_blab(OTA+cell); - - } else { -/* Restore the old burns */ - burn_restore(naxis1-ovrscan1, naxis2-ovrscan2, naxis1, - buf, OTA+cell); + if(VERBOSE & VERB_NORM) persist_blab(OTA+cell); + } } /* Write the corrected data back to the FITS. */ - if(update) { + if(!tableonly && update) { fh_ehu(ehu, 0); /* Seek back to the start of data */ if (fh_write_padded_image(ehu, fh_file_desc(ehu), buf, @@ -590,10 +606,10 @@ /* Dump out the postage stamp file */ - if(psffile != NULL) { + if(psffile != NULL && !restore && !apply) { psf_write(psfsize, psfsize, OTA, otanum, psffile); } /* Dump out the PSF stats */ - if(psfstatfile != NULL) { + if(psfstatfile != NULL && !restore && !apply) { psf_write_stats(psfsize, psfsize, OTA, otanum, psfstatfile, psfavg); } @@ -601,11 +617,15 @@ /* Write burn info to FITS file. */ if(update) persist_fits_write(OTA, ihu); - - if(restore) { + if(restore || tableonly) { /* Indicate in the header that the burns are not applied. */ - fh_set_bool(ihu, FH_AUTO, PHU_NAME_BURN_APPLIED, - FH_FALSE, PHU_COMMENT_BURN_APPLIED); - fh_rewrite(ihu); + fh_set_bool(ihu, FH_AUTO, PHU_NAME_BURN_APPLIED, + FH_FALSE, PHU_COMMENT_BURN_APPLIED); + fh_rewrite(ihu); + } else if(apply) { + /* Indicate in the header that the burns have been applied. */ + fh_set_bool(ihu, FH_AUTO, PHU_NAME_BURN_APPLIED, + FH_TRUE, PHU_COMMENT_BURN_APPLIED); + fh_rewrite(ihu); } @@ -711,6 +731,8 @@ printf(" cell=N Work on just one cell? Cell count [0:63] mode.\n"); printf(" mask=0101... 64 digits to work on cells 0:63.\n"); - printf(" update={t|f} Modify the input MEF by subtracting fits?\n"); + printf(" update={t|f} Modify the input MEF writing table and subtracting fits?\n"); printf(" restore={t|f} Restore the input MEF by adding input fits?\n"); + printf(" apply={t|f} Modify the input MEF by subtracting previously calculated fits?\n"); + printf(" tableonly={t|f} Calculate fits but do *not* modify the input MEF images, only write tables\n"); printf(" in=fname Input file for previous burn persistence streaks\n"); printf(" infits=fname Input FITS file for previous burn persistence streaks (stored\n"); Index: branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.h =================================================================== --- branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.h (revision 24238) +++ branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.h (revision 24557) @@ -173,4 +173,5 @@ MTYPE *mask, OBJBOX *box); STATIC int burn_restore(int nx, int ny, int NX, IMTYPE *buf, CELL *cell); +STATIC int burn_apply(int nx, int ny, int NX, IMTYPE *buf, CELL *cell); STATIC int persist_read(CELL *cell, const char *infile); STATIC int persist_write(CELL *cell, const char *outfile); Index: branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/burntool/man/burntool.1 =================================================================== --- branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/burntool/man/burntool.1 (revision 24238) +++ branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/burntool/man/burntool.1 (revision 24557) @@ -66,7 +66,54 @@ inhibited from subtracting the fits by "update=f". + Two other modes can be invoked using "tableonly=t" and "apply=t". + The first mode is just like "update=t", except that it does not write + the modified images, only the FITS table and/or tabulated fit data. + The second mode is like "restore=t", except that it subtracts the + tabulated fits instead of adding them back in. Note that this requires + slightly different specifications for the "in=" file: + + update or tableonly: in=previous_file_for_persistence + restore or apply: in=this_file_for_burn_and_persistence + + Thus a sequence of images could be processed (leaving off usual + arguments such as mask=, etc) to immediately fit and subtract: + + burntool o4991g0001o24.fits \ + out=o4991g0001o24.burn \ + psfstat=o4991g0001o24.stat psf=o4991g0001o24.psf + + burntool o4991g0002o24.fits \ + in=o4991g0001o24.burn out=o4991g0002o24.burn \ + psfstat=o4991g0002o24.stat psf=o4991g0002o24.psf + + burntool o4991g0003o24.fits \ + in=o4991g0002o24.burn out=o4991g0003o24.burn \ + psfstat=o4991g0003o24.stat psf=o4991g0003o24.psf + + or as a two step process of generating tables: + + burntool o4991g0001o24.fits tableonly=t \ + out=o4991g0001o24.burn \ + psfstat=o4991g0001o24.stat psf=o4991g0001o24.psf + + burntool o4991g0002o24.fits tableonly=t \ + in=o4991g0001o24.burn out=o4991g0002o24.burn \ + psfstat=o4991g0002o24.stat psf=o4991g0002o24.psf + + burntool o4991g0003o24.fits tableonly=t \ + in=o4991g0002o24.burn out=o4991g0003o24.burn \ + psfstat=o4991g0003o24.stat psf=o4991g0003o24.psf + + and then later on carrying out the subtractions of the fits: + + burntool o4991g0001o24.fits apply=t + burntool o4991g0002o24.fits apply=t + burntool o4991g0003o24.fits apply=t + The identification of significant (but unsaturated) stars is an important component of burntool's mask generation, and if requested burntool can assemble a gallery of postage stamps of suitable stars. + This does not work with "restore=t" or "apply=t" of course, since + these work only on previous tables and do not examine the images. Any star whose maximum exceeds PSF_THRESH above sky but is not deemed @@ -125,7 +172,16 @@ qt=-5.201 Tangential quadrupole [pix^2] (-99.99 if fails) q+ * cos(2*phi) + qx * sin(2*phi) + q3c=-1.193 Cosine trefoil measure [pix^2] (-99.99 if fails) + q3s=-0.038 Sine trefoil measure [pix^2] (-99.99 if fails) qpm=5.154 Plus quadrupole averaged over psfavg qcm=0.560 Cross quadrupole averaged over psfavg qtm=-5.201 Tangential averaged over psfavg + q3cm=-1.230 Cosine trefoil averaged over psfavg + q3sm=-0.079 Sine trefoil averaged over psfavg + + The net, polar coordinate quadrupole can be assembled from + + Qmagnitude = sqrt(qpm^2+qcm^2) + Qtheta = 0.5 * atan2(qcm, qpm) Helpful utilities include: @@ -151,8 +207,18 @@ update={t|f} - Modify the input MEF by subtracting fits (default t)? + Calculate fits and modify the input MEF images by subtracting + the trail fits and writing a FITS table (default t). If + false the MEF is not altered, although PSF can be calculated. restore={t|f} - Restore the input MEF by adding input fits (default f)? + Restore the input MEF by adding input fits (default f) + + apply={t|f} + Use previously calculated fits to modify the input MEF image + by subtraction (default f) + + tableonly={t|f} + Calculate fits but do *not* modify the input MEF images, + only write the fit data as a FITS and/or text table (default f) in=fname Index: branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/burntool/psfstamp.c =================================================================== --- branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/burntool/psfstamp.c (revision 24238) +++ branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/burntool/psfstamp.c (revision 24557) @@ -136,5 +136,5 @@ int i, k, l, nstar, fdout, otacx, otacy, xid, yid, ntot=0, sumax; int cellcount, ota_xid, ota_yid; - double scale, phi, fwhm[3], q[3], qt, xfp, yfp, pi=4*atan(1.0); + double scale, phi, fwhm[3], q[5], qt, xfp, yfp, pi=4*atan(1.0); IMTYPE *median_image; CELL *cell; @@ -203,6 +203,7 @@ qt = -q[1] * cos(2*phi) - q[2] * sin(2*phi); - printf("N= %d PSFmaj= %.2f min= %.2f theta= %.1f m2= %.2f q+= %.3f qx= %.3f qt= %.3f\n", - nstar, fwhm[0], fwhm[1], fwhm[2]*180/pi, q[0], q[1], q[2], qt); + printf("N= %d PSFmaj= %.2f min= %.2f theta= %.1f m2= %.2f q+= %.3f qx= %.3f qt= %.3f q3c= %.3f q3s= %.3f\n", + nstar, fwhm[0], fwhm[1], fwhm[2]*180/pi, q[0], q[1], q[2], qt, + q[3], q[4]); if(CONCAT_FITS) { @@ -258,9 +259,12 @@ double m2[MAXPSFMEDIAN], qp[MAXPSFMEDIAN], qc[MAXPSFMEDIAN]; double qt[MAXPSFMEDIAN], fwavg[MAXPSFMEDIAN]; + double q3c[MAXPSFMEDIAN], q3s[MAXPSFMEDIAN]; double qpavg[MAXPSFMEDIAN], qcavg[MAXPSFMEDIAN], qtavg[MAXPSFMEDIAN]; + double q3cavg[MAXPSFMEDIAN], q3savg[MAXPSFMEDIAN]; int nstar[MAXCELL], nfw[MAXCELL]; double fwmed[MAXCELL], m2med[MAXCELL]; double qpmed[MAXCELL], qcmed[MAXCELL], qtmed[MAXCELL]; - double qpmacro, qcmacro, qtmacro, fwmacro; + double q3cmed[MAXCELL], q3smed[MAXCELL]; + double qpmacro, qcmacro, qtmacro, fwmacro, q3cmacro, q3smacro; FILE *fp; CELL *cell; @@ -304,4 +308,6 @@ qp[nstar[k]] = q[1]; qc[nstar[k]] = q[2]; + q3c[nstar[k]] = q[3]; + q3s[nstar[k]] = q[4]; /* Get the position in the focal plane and therefore the qt statistic */ psc_cell_to_pixel(cellx, celly, 0.5*PSC_HCELL/PSC_PIXEL, @@ -325,6 +331,9 @@ qcmed[k] = double_median(nstar[k], qc); qtmed[k] = double_median(nstar[k], qt); + q3cmed[k] = double_median(nstar[k], q3c); + q3smed[k] = double_median(nstar[k], q3s); } else { m2med[k] = qpmed[k] = qcmed[k] = qtmed[k] = -99.99; + q3cmed[k] = q3smed[k] = -99.99; } /* Toss these results into the macrocell median hopper */ @@ -337,4 +346,6 @@ qcavg[nqavg] = qc[l]; qtavg[nqavg] = qt[l]; + q3cavg[nqavg] = q3c[l]; + q3savg[nqavg] = q3s[l]; nqavg++; } @@ -346,5 +357,7 @@ qcmacro = double_median(nqavg, qcavg); qtmacro = double_median(nqavg, qtavg); - if(nqavg == 0) qpmacro = qcmacro = qtmacro = -99.99; + q3cmacro = double_median(nqavg, q3cavg); + q3smacro = double_median(nqavg, q3savg); + if(nqavg == 0) qpmacro = qcmacro = qtmacro = q3cmacro = q3smacro = -99.99; fwmacro = double_median(nfwave, fwavg); @@ -357,9 +370,9 @@ cell = ota + k; - fprintf(fp, "ext=xy%1d%1d bias=%d sky=%d rmssky=%d npsf=%d fwhm=%.2f fwmed=%.2f m2=%.2f qp=%.3f qc=%.3f qt=%.3f qpm=%.3f qcm=%.3f qtm=%.3f\n", + fprintf(fp, "ext=xy%1d%1d bias=%d sky=%d rmssky=%d npsf=%d fwhm=%.2f fwmed=%.2f m2=%.2f qp=%.3f qc=%.3f qt=%.3f q3c=%.3f q3s=%.3f qpm=%.3f qcm=%.3f qtm=%.3f q3cm=%.3f q3sm=%.3f\n", cellx, celly, cell->bias, cell->sky, cell->rms, nstar[k], fwmed[k], fwmacro, m2med[k], - qpmed[k], qcmed[k], qtmed[k], - qpmacro, qcmacro, qtmacro); + qpmed[k], qcmed[k], qtmed[k], q3cmed[k], q3smed[k], + qpmacro, qcmacro, qtmacro, q3cmacro, q3smacro); } } Index: branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/burntool/psfstats.c =================================================================== --- branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/burntool/psfstats.c (revision 24238) +++ branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/burntool/psfstats.c (revision 24557) @@ -10,10 +10,9 @@ #include "burntool.h" #include "math.h" -#include "psf/psf.h" /****************************************************************/ /* psf_stats(): Tell us about this PSF star */ -STATIC int psf_stats(int nx, int ny, IMTYPE *data, int bias, - double *fwhm, double *q) +STATIC int psf_stats(int nx, int ny, IMTYPE *data, int bias, + double *fwhm, double *q) { return(1); Index: branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/libfh/fh.c =================================================================== --- branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/libfh/fh.c (revision 24238) +++ branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/libfh/fh.c (revision 24557) @@ -200,5 +200,5 @@ /* Functions that do file operations WITH RETRY: */ static int write_file(int fd, const void* buf, int len); -static int read_file(int fd, const void* buf, int len); +static int read_file(int fd, void* buf, int len); static int seek_file(int fd, off_t off); /* Performs lseek(SEEK_SET) */ static int fh_lock_file(HeaderUnit hu, int fd, fh_mode mode); @@ -234,4 +234,5 @@ * this information in debug mode. */ +const char* fh_rcs_version(void); /* proto to make gcc happy */ const char* fh_rcs_version(void) { return rcs_id; } @@ -274,5 +275,5 @@ static void -null_log(const char* s) { if (s) ; } +null_log(const char* s) { if (s) return; } /* * It's also possible to pass a NULL to fh_log_xxx() and rely only on the @@ -309,5 +310,5 @@ HeaderUnit -fh_create() +fh_create(void) { HeaderUnitStruct* rtn; @@ -1941,5 +1942,6 @@ HeaderUnitStruct* list = FH_HU(hu); char writebuf[FH_BLOCK_SIZE]; - int i = 0, writelen; + int i = 0; + unsigned int writelen; int blocks_left; fh_result result; @@ -2446,5 +2448,5 @@ else /* byte-swapping case requires a temporary buffer */ { - unsigned char outbuf[32768]; + unsigned char* outbuf[32768]; while (bytes_left) @@ -2811,5 +2813,32 @@ while (len) { - rtn = write(fd, (char*)buf, len); + rtn = write(fd, (const char*)buf, len); + + switch (rtn) + { + case -1: if (errno == EINTR || errno == EAGAIN) + continue; /* retry */ + return -1; /* permanent failure */ + case 0: return count; + default: + { + len -= rtn; + count += rtn; + buf = (const char*)buf + rtn; + } + } + } + return count; +} + +static int +read_file(int fd, void* buf, int len) +{ + int rtn; + int count = 0; + + while (len) + { + rtn = read(fd, (char*)buf, len); switch (rtn) @@ -2831,31 +2860,4 @@ static int -read_file(int fd, const void* buf, int len) -{ - int rtn; - int count = 0; - - while (len) - { - rtn = read(fd, (char*)buf, len); - - switch (rtn) - { - case -1: if (errno == EINTR || errno == EAGAIN) - continue; /* retry */ - return -1; /* permanent failure */ - case 0: return count; - default: - { - len -= rtn; - count += rtn; - buf = (char*)buf + rtn; - } - } - } - return count; -} - -static int seek_file(int fd, off_t off) { @@ -2969,5 +2971,5 @@ fh_compare(const void* a, const void* b) { - double diff = (*((FitsCard**)a))->idx - (*((FitsCard**)b))->idx; + double diff = (*((FitsCard* const*)a))->idx - (*((FitsCard* const*)b))->idx; if (diff < 0) return -1; Index: branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/libfh/fh_table.c =================================================================== --- branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/libfh/fh_table.c (revision 24238) +++ branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/libfh/fh_table.c (revision 24557) @@ -180,8 +180,10 @@ /* Write the format specifiers for each column. */ - for(i=0; i < table->num_cols; i++) - { + for (i=0; i < table->num_cols; i++) + { + int n; + idx = 5.0 + i/10.; - int n = i + 1; /* FITS numbering scheme starts from 1. */ + n = i + 1; /* FITS numbering scheme starts from 1. */ snprintf(keyword, FH_NAME_SIZE + 1, "TFORM%d", n); if((table->cols[i].format == FH_TABLE_FORMAT_FLOAT) || @@ -214,8 +216,10 @@ * nice to have from a self-documenting perspective, isn't strictly * necessary to interpret the actual table values. */ - for(i=0; i < table->num_cols; i++) - { + for (i=0; i < table->num_cols; i++) + { + int n; + idx = 7.0 + i/10.; - int n = i + 1; /* FITS numbering scheme starts from 1. */ + n = i + 1; /* FITS numbering scheme starts from 1. */ snprintf(keyword, FH_NAME_SIZE + 1, "TTYPE%d", n); fh_set_str(hu, idx, keyword, table->cols[i].name, Index: branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/libfhreg/gpc_telescope.h =================================================================== --- branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/libfhreg/gpc_telescope.h (revision 24238) +++ branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/libfhreg/gpc_telescope.h (revision 24557) @@ -89,4 +89,5 @@ MK_PFL( 523.3, COMALT ,6, "Commanded telescope pointing altitude (degrees)" ) MK_PFL( 523.4, COMROT ,6, "Commanded telescope rotator angle (degrees)" ) +MK_PFL( 523.5, MOONANG ,6, "Angular distance to moon (degrees)" ) MK_CMT( 530.00,CMTTELO1, "NOTE: Telescope RA DEC or ALT AZ already include") MK_CMT( 530.10,CMTTELO2, " the following offsets. Do not re-apply!" ) Index: branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/tiltystreak/Make.Common =================================================================== --- branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/tiltystreak/Make.Common (revision 24557) +++ branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/tiltystreak/Make.Common (revision 24557) @@ -0,0 +1,1 @@ +link ../Make.Common Index: branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/tiltystreak/Makefile =================================================================== --- branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/tiltystreak/Makefile (revision 24557) +++ branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/tiltystreak/Makefile (revision 24557) @@ -0,0 +1,16 @@ +include ../Make.Common + +VERSION=0.00 +CCWARN+=$(WERROR) +CCDEFS+=$(VERSIONDEFS) +$(EXECNAME): $(OBJS) libfh.a + +include ../Make.Common + +# Dependencies by Make.Common $Revision: 2.17 $ + +$(OBJ)/sort.o: sort.c + +$(OBJ)/tiltystreak.o: tiltystreak.c fh/fh.h \ + fhreg/general.h fhreg/macros.h \ + fhreg/gpc_detector.h Index: branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/tiltystreak/linearfit.c =================================================================== --- branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/tiltystreak/linearfit.c (revision 24557) +++ branches/eam_branches/20090522/extsrc/gpcsw/gpcsrc/fits/tiltystreak/linearfit.c (revision 24557) @@ -0,0 +1,129 @@ +/* Routine to fit a linear function */ +/* 060304 v1.2 add wlinearfit() */ +/* 060130 v1.1 replaced static allocations with calloc's */ +/* 040208 v1.0 John Tonry (derived from fortran) */ +#include +#include +#include + +#define ABS(x) ((x)<0?(-(x)):(x)) + +static int linsolve(int n, double *y, double *a, double *x); + +int linearfit(int npt, double *y, int nx, double *x, double *param) +{ +/* Y = NPT data being fitted as a linear function of X[*][NX] + * PARAM = returned parameters of fit: + * Y[*] = PARAM[0]*X[*][0], + * + PARAM[1]*X[*][1], + * ... + * + PARAM[NX-1]*X[*][NX-1] + * + */ + int i, j, k; + double *v, *a; + v = (double *)calloc(nx, sizeof(double)); + a = (double *)calloc(nx*nx, sizeof(double)); + + if(npt < nx) { + fprintf(stderr, "linearfit: Too few points provided: %d %d\n", npt, nx); + return(-2); + } + +#if 0 /* Show us data to solve? */ + for(i=0; i rat) { + k = i; + rat = ABS(a[j+k*n]); + } + } + rowstatus[k] = 1; + row[j] = k; + if(rat == 0.0) { + fprintf(stderr, "WHOA: singular matrix, iter %d, row %d\n", j, k); + free(rowstatus); + free(row); + return(1); + } +/* Subtract away matrix below diagonal */ + for(i=0; i=0; j--) { + x[j] = y[row[j]]; + for(k=j+1; k +#include +#include + +#define MAXSTACK 256 +//#define NSTOP 15 +#define NSTOP 8 + +/* Quicksort program 1980 (John Tonry) */ +int qsort_dbl(int n, double *x) +{ + double key, kl, kr, km, temp; + int l, r, m, lstack[MAXSTACK], rstack[MAXSTACK], sp; + register int i, j, k; + int mgtl, lgtr, rgtm; + + sp = 0; + lstack[sp] = 0; + rstack[sp] = n-1; + + while(n > NSTOP && sp >= 0) { +/* Sort a subrecord off the stack */ + l = lstack[sp]; + r = rstack[sp]; + sp--; + m = (l + r) / 2; +/* Set KEY = median of X(L), X(M), X(R) */ + kl = x[l]; + km = x[m]; + kr = x[r]; + mgtl = km > kl; + rgtm = kr > km; + lgtr = kl > kr; +#ifndef STD_KEY + if(mgtl ^ rgtm) { + if(mgtl ^ lgtr) key = kr; + else key = kl; + } else { + key = km; + } +#else +/* Curiously enough, this non-median seems to work as well or better */ + if(mgtl ^ rgtm) { + key = km; + } else { + if(mgtl ^ lgtr) key = kl; + else key = kr; + } +#endif + i = l; + j = r; + while(1) { +/* Find a big record on the left */ + while(x[i] < key) i++; + +/* Find a small record on the right */ + while(x[j] > key) j--; + + if(i >= j) break; +/* Exchange records */ + temp = x[i]; + x[i] = x[j]; + x[j] = temp; + i++; + j--; + } + +/* Subfile is partitioned into two halves, left .le. right */ +/* Push the two halves on the stack */ + if(j-l+1 > NSTOP) { + sp = sp + 1; + lstack[sp] = l; + rstack[sp] = j; + } + if(r-j > NSTOP) { + sp = sp + 1; + lstack[sp] = j+1; + rstack[sp] = r; + } + if(sp >= MAXSTACK) { +// fprintf(stderr,"QSORT8: Fatal error from stack overflow\n"); +// fprintf(stderr,"Fall back on sort by insertion\n"); + break; + } + } + +/* Sorting routine that sorts the N elements of single precision */ +/* array X by straight insertion between previously sorted numbers */ + for(j=n-2; j>=0; j--) { + for(i=j+1, k=j; i NSTOP && sp >= 0) { +/* Sort a subrecord off the stack */ + l = lstack[sp]; + r = rstack[sp]; + sp--; + m = (l + r) / 2; +/* Set KEY = median of X(L), X(M), X(R) */ + kl = x[l]; + km = x[m]; + kr = x[r]; + mgtl = km > kl; + rgtm = kr > km; + lgtr = kl > kr; +#ifndef STD_KEY + if(mgtl ^ rgtm) { + if(mgtl ^ lgtr) key = kr; + else key = kl; + } else { + key = km; + } +#else +/* Curiously enough, this non-median seems to work as well or better */ + if(mgtl ^ rgtm) { + key = km; + } else { + if(mgtl ^ lgtr) key = kl; + else key = kr; + } +#endif + i = l; + j = r; + while(1) { +/* Find a big record on the left */ + while(x[i] < key) i++; + +/* Find a small record on the right */ + while(x[j] > key) j--; + + if(i >= j) break; +/* Exchange records */ + temp = x[i]; + x[i] = x[j]; + x[j] = temp; + i++; + j--; + } + +/* Subfile is partitioned into two halves, left .le. right */ +/* Push the two halves on the stack */ + if(j-l+1 > NSTOP) { + sp = sp + 1; + lstack[sp] = l; + rstack[sp] = j; + } + if(r-j > NSTOP) { + sp = sp + 1; + lstack[sp] = j+1; + rstack[sp] = r; + } + if(sp >= MAXSTACK) { +// fprintf(stderr,"QSORT8: Fatal error from stack overflow\n"); +// fprintf(stderr,"Fall back on sort by insertion\n"); + break; + } + } + +/* Sorting routine that sorts the N elements of single precision */ +/* array X by straight insertion between previously sorted numbers */ + for(j=n-2; j>=0; j--) { + for(i=j+1, k=j; i +#include +#include + +#include "math.h" +#include "fh/fh.h" +#include "fhreg/general.h" +#include "fhreg/gpc_detector.h" + +#define MAX(a,b) (((a) > (b)) ? (a) : (b)) +#define MIN(a,b) (((a) < (b)) ? (a) : (b)) +#define NINT(x) (x<0?(int)((x)-0.5):(int)((x)+0.5)) +#define ABS(a) (((a) > 0) ? (a) : -(a)) + +#define NODATA 0 /* Marker for *No Data* */ +#define NCELL 8 /* Number of cells in a row */ +#define SAT4SURE 45000 /* Blasted, we don't care about pixels above */ +#define NPOLY 3 /* N terms in of flattening polynomial */ + +//#define DEBUG /* Turn on debug output? */ + +/* Info for a cell */ +typedef struct cell_info { + HeaderUnit hu; /* Header unit for this cell */ + int cell; /* Cell number */ + int mefxt; /* Where is it in the MEF? */ + int bias; /* Bias level */ + int sky; /* Sky level */ + int rms; /* RMS in the sky */ + int nsat; /* Saturated pixel count */ + double ampl; /* Streak amplitude */ + short *data; /* Image data */ + int *mask; /* Mask bits */ +} CELL; + +/* Prototypes */ +static int cell_stats(int sx, int nx, int ny, int NX, int NY, int bzero, + CELL *cell); +static int int_median(int n, int *key); +int qsort_dbl(int n, double *x); +int qsort_int(int n, int *x); +static int streak_ratio(int sx, int nx, int ny, int NX, CELL otarow[]); +static int linear_streak_fit(int sx, int nx, int ny, int NX, CELL otarow[]); +static int bilinear_streak_fit(int sx, int nx, int ny, int NX, CELL otarow[]); +static int parabola_streak_fit(int sx, int nx, int ny, int NX, CELL otarow[]); +int linearfit(int npt, double *y, int nx, double *x, double *param); +static void syntax(const char *prog); + +static int nmedian_buf=0; +static int *median_buf=NULL; +static double *Abuf=NULL, *Bbuf=NULL, *Xbuf=NULL; + +int main(int argc, const char* argv[]) +{ + int i, j, row, col, ext; + HeaderUnit ihu = fh_create(); + HeaderUnit ehu = NULL; + const char* ifilename; + char extname[FH_MAX_STRLEN+1], xtension[FH_MAX_STRLEN+1], cellsought[5]; + int nextend; + CELL otarow[NCELL]; + int naxis, naxis1=0, naxis2=0; + int n1, n2; + int prescan1, ovrscan1, ovrscan2; + double bzero_d; + int bzero = 32768, firstime = 1, algo=2; + int cellmask[NCELL*NCELL]; + + if(argc < 2) { + syntax(argv[0]); + fprintf(stderr, "error: syntax: tiltystreak Uncompressed_GPC1_MEF.fits\n"); + exit(EXIT_FAILURE); + } + ifilename = argv[1]; + +/* Defaults */ + for(i=0; i nextend) { + fprintf(stderr, + "error: Could not find extension `%s'\n", cellsought); + exit(EXIT_FAILURE); + } + +/* Positioned at the desired cell now */ + if (fh_get_PRESCAN1(ehu, &prescan1) != FH_SUCCESS || + fh_get_OVRSCAN1(ehu, &ovrscan1) != FH_SUCCESS || + fh_get_OVRSCAN2(ehu, &ovrscan2) != FH_SUCCESS || + fh_get_NAXIS(ehu, &naxis) != FH_SUCCESS || + fh_get_NAXIS1(ehu, &n1) != FH_SUCCESS || + fh_get_NAXIS2(ehu, &n2) != FH_SUCCESS) { + fprintf(stderr, "error: Cannot get NAXIS* in extension %d row %d col %d\n", ext, row, col); + exit(EXIT_FAILURE); + } +/* Cell isn't an image or is 32bpp */ + if(naxis != 2) { + otarow[col].mefxt = -1; + continue; + } + +/* Make some space now that we finally know the size */ + if(firstime) { + naxis1 = n1; + naxis2 = n2; + firstime = 0; + for(i=0; i 0) i++; + if(i == 0) continue; /* Nope, entirely wiped out */ + +/* Process this row of cells */ + +/* What's the stats in each cell, generate mask */ + for(col=0; colbias = 0; + } else { + for(j=1; jdata[sx+nx+(NX-nx-sx)/2+j*NX] + bzero; + cell->bias = int_median(ny-2, median_buf); + } + +/* Get sky stats */ + for(k=n=0; kdata[sx+i+NX*j]+bzero < SAT4SURE && + cell->data[sx+i+NX*j]+bzero != NODATA) { + median_buf[n++] = cell->data[sx+i+NX*j] + bzero; + } + } + + if(n < 20) { /* Better have hit at least 20! */ + cell->sky = cell->rms = 0.0; + return(-1); + } + +/* First pass at sky and quartile */ + cell->sky = int_median(n, median_buf); + cell->rms = 1.33*(cell->sky - median_buf[n/4]); + +/* Clip at 3 sigma */ + for(j=0; median_buf[j] < cell->sky - Q_CLIP*cell->rms && j cell->sky + Q_CLIP*cell->rms && k>n/2; k--); + cell->sky = median_buf[(j+k)/2]; + cell->rms = 1.33 * (cell->sky - median_buf[(3*j+k)/4]); + cell->sky -= cell->bias; + +/* Create mask */ + thresh = cell->bias + cell->sky + MASK_CLIP*cell->rms; + cell->nsat = 0; + for(j=0; jmask[i+NX*j] = 1; + if(cell->data[i+NX*j]+bzero > thresh) cell->mask[i+NX*j] = 0; + if(cell->data[i+NX*j]+bzero > SAT4SURE) (cell->nsat)++; + } + } + + return(0); +} + +static int int_median(int n, int *key) +{ + if(n == 0) return(0); + qsort_int(n, key); + return((key[n/2]+key[(n-1)/2]) / 2); +} + +/* Determine the relative streakiness of all cells in this row */ +/* + * Compare each cell's pixels in the difference which are greater + * than MINSTREAK with the average using pixels surviving the + * intersection of all the masks + */ +#define MINSTREAK 50 +#define EDGE 3 +static int streak_ratio(int sx, int nx, int ny, int NX, CELL otarow[]) +{ + int i, j, k, nalive=0, ncontrib; + double diff; + + for(k=0; k 0) { + nalive++; + } else { + otarow[k].ampl = 0.0; + } + } + +/* Loop over all the pixels */ + ncontrib = 0; + for(j=1+EDGE; j MINSTREAK)) { + for(k=0; k 0) { + k = int_median(n, medimg) + Abuf[j-1] * (i-nx*sx); + sy += k; + s1 += i - (sx+nx); + } + } + Abuf[j] = 0.0; + if(s1 != 0) Abuf[j] = sy / s1; + for(k=0, n=1; k 0) { + k = int_median(n, medimg); + if(i < sx+nx/3) { + s1 += k + Abuf[j-1]; + n1++; + } else { + s2 += k + Bbuf[j-1]; + n2++; + } + } + } + Abuf[j] = Bbuf[j] = 0.0; + if(n1 > 0) Abuf[j] = s1/n1; + if(n2 > 0) Bbuf[j] = s2/n2; + for(k=0, n=1; k 0) { + k = int_median(n, medimg); + sy += k; + s1 += i - (sx+nx); + } + } + if(s1 != 0) A = sy / s1; + +/* Fit a full line z = B*(i-nx-sx)+C to the left third */ + s0 = s1 = s2 = sm = B = C = 0.0; + npt = 0; + for(i=sx+EDGE; i 0) { + k = int_median(n, medimg); + s1 += k; + s0 += (i-sx-nx); + s2 += (i-sx-nx)*(i-sx-nx); + sm += (i-sx-nx)*k; + npt++; + } + } + if(npt > 0) { + B = (npt*sm-s0*s1) / (npt*s2-s0*s0); + C = (s2*s1 - s0*sm) / (npt*s2-s0*s0); + } + + if(A != B) { + i0 = C / (A-B) + sx+nx; + } else { + i0 = 0; + } + +#ifdef DEBUG +// if(otarow[3].mefxt == 4) printf("%4d %4d %9.3f %9.3f %9.3f\n", +// j, i0, A, B, C); +#endif + +/* Correct all the cells */ + for(k=0; k