Changeset 24391

Timestamp:

Jun 12, 2009, 3:53:07 PM (17 years ago)

Author:

eugene

Message:

updates to burntool-and-related from camera svn

Location:

trunk/extsrc/gpcsw/gpcsrc

Files:

• analysis/libpscoords/man (added)
• analysis/libpscoords/man/libpscoords.1 (added)
• analysis/libpscoords/pscoords.c (modified) (8 diffs)
• analysis/libpscoords/pscoords.h (modified) (6 diffs)
• fits/burntool/Makefile (modified) (1 diff)
• fits/burntool/burnfix.c (modified) (1 diff)
• fits/burntool/burntool.c (modified) (14 diffs)
• fits/burntool/burntool.h (modified) (1 diff)
• fits/burntool/man/burntool.1 (modified) (3 diffs)
• fits/burntool/psfstamp.c (modified) (8 diffs)
• fits/libfh/fh.c (modified) (9 diffs)
• fits/libfh/fh_table.c (modified) (2 diffs)
• fits/libfhreg/gpc_telescope.h (modified) (1 diff)

trunk/extsrc/gpcsw/gpcsrc/analysis/libpscoords/pscoords.c

@@ -24 +24 @@
 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 +99 @@
    { -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 +175 @@
 #endif
 
-static int DOOFFSET;    /* Do these offsets or no? */
+static int DOOFFSET=1;  /* Do these offsets or no? */
 
 #if 0
@@ -143 +182 @@
 #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<PSC_NX*PSC_NY; i++) {
+      if(fgets(rbuf, 1024, fp) == NULL || 
+         sscanf(rbuf, "%lf %lf %lf", &psc_otaoff[i].dx, 
+                &psc_otaoff[i].dy, &psc_otaoff[i].rot) != 3) {
+         fprintf(stderr, "error: Error reading OTA offset file '%s' line %d\n",
+                 fname, i);
+         fclose(fp);
+         return(-1);
+      }
+   }
+   fclose(fp);
+   return(0);
+}
+
+/* Load up a new offset table? */
 int psc_do_chipoff(int doit)
 {
@@ -196 +259 @@
 /* Position wrt OTA center after correction for offset and rotation */
    if(DOOFFSET) {
-      x =  (xp-psc_otaoff[idx].dx)*cos(psc_otaoff[idx].rot*0.001) + 
-         (yp-psc_otaoff[idx].dy)*sin(psc_otaoff[idx].rot*0.001);
-      y = -(xp-psc_otaoff[idx].dx)*sin(psc_otaoff[idx].rot*0.001) + 
-         (yp-psc_otaoff[idx].dy)*cos(psc_otaoff[idx].rot*0.001);
+      x =  (xp-psc_otaoff[idx].dx)*cos(psc_otaoff[idx].rot*MRAD2) + 
+         (yp-psc_otaoff[idx].dy)*sin(psc_otaoff[idx].rot*MRAD2);
+      y = -(xp-psc_otaoff[idx].dx)*sin(psc_otaoff[idx].rot*MRAD2) + 
+         (yp-psc_otaoff[idx].dy)*cos(psc_otaoff[idx].rot*MRAD2);
    } else {
       x = xp;
@@ -228 +291 @@
 /* Position wrt nominal OTA center after its offset and rotation */
    if(DOOFFSET) {
-      xp = psc_otaoff[idx].dx + x*cos(psc_otaoff[idx].rot*0.001) - y*sin(psc_otaoff[idx].rot*0.001);
-      yp = psc_otaoff[idx].dy + x*sin(psc_otaoff[idx].rot*0.001) + y*cos(psc_otaoff[idx].rot*0.001);
+      xp = psc_otaoff[idx].dx + x*cos(psc_otaoff[idx].rot*MRAD2) - y*sin(psc_otaoff[idx].rot*MRAD2);
+      yp = psc_otaoff[idx].dy + x*sin(psc_otaoff[idx].rot*MRAD2) + y*cos(psc_otaoff[idx].rot*MRAD2);
    } else {
       xp = x;
@@ -292 +355 @@
 
 /* Evaluate an optical model: arcsec -> 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 +412 @@
    return(0);
 }
-
 
 /* Project RA,Dec (a,d) to the tangent plane (x,y) at (a0,d0,pa) [radians] */

trunk/extsrc/gpcsw/gpcsrc/analysis/libpscoords/pscoords.h

@@ -2 +2 @@
 #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 +130 @@
 
 /* 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 +157 @@
  * ---------------------------------------------------------------------------
  *   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 +169 @@
  *   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 +180 @@
 /**************/
 /* 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 +227 @@
                  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);

trunk/extsrc/gpcsw/gpcsrc/fits/burntool/Makefile

@@ -32 +32 @@
   pscoords/pscoords.h
 
-$(OBJ)/psfstats.o: psfstats.c burntool.h
+$(OBJ)/psfstats.o: psfstats.c burntool.h psf/psf.h
 
 $(OBJ)/sort.o: sort.c

trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burnfix.c

@@ -105 +105 @@
       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; k<cell->npersist; k++) {
+      if(!cell->persist[k].fiterr) 
+         sub_fit(nx, ny, NX, buf, &(cell->persist[k]), +1);
    }
    return(0);

trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.c

@@ -33 +33 @@
    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 +102 @@
 /* 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 +139 @@
          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 +160 @@
 
 /* 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 +172 @@
 
 /* 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 +332 @@
    /* 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 +344 @@
 
 /* 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 +379 @@
                "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 +501 @@
 
       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 +552 @@
          }
 
-
-         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 +579 @@
             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 +606 @@
 
 /* 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 +617 @@
 /* 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 +731 @@
    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");

trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.h

@@ -173 +173 @@
                      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);

trunk/extsrc/gpcsw/gpcsrc/fits/burntool/man/burntool.1

@@ -66 +66 @@
         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 +172 @@
            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 +207 @@
 
         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       

trunk/extsrc/gpcsw/gpcsrc/fits/burntool/psfstamp.c

@@ -136 +136 @@
    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 +203 @@
    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 +259 @@
    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 +308 @@
                      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 +331 @@
                   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 +346 @@
                      qcavg[nqavg] = qc[l];
                      qtavg[nqavg] = qt[l];
+                     q3cavg[nqavg] = q3c[l];
+                     q3savg[nqavg] = q3s[l];
                      nqavg++;
                   }
@@ -346 +357 @@
          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 +370 @@
                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);
             }
          }

trunk/extsrc/gpcsw/gpcsrc/fits/libfh/fh.c

@@ -200 +200 @@
 /* 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 +234 @@
  * 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 +275 @@
 
 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 +310 @@
 
 HeaderUnit
-fh_create()
+fh_create(void)
 {
    HeaderUnitStruct* rtn; 
@@ -1941 +1942 @@
    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 +2448 @@
    else /* byte-swapping case requires a temporary buffer */
    {
-      unsigned char outbuf[32768];
+      unsigned char* outbuf[32768];
 
       while (bytes_left)
@@ -2811 +2813 @@
    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 +2860 @@
 
 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 +2971 @@
 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;

trunk/extsrc/gpcsw/gpcsrc/fits/libfh/fh_table.c

@@ -180 +180 @@
 
    /* 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 +216 @@
     * 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, 

trunk/extsrc/gpcsw/gpcsrc/fits/libfhreg/gpc_telescope.h

@@ -89 +89 @@
 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!"  )