Changeset 23924 for trunk/extsrc/gpcsw/gpcsrc/fits/burntool

Timestamp:

Apr 18, 2009, 11:11:15 AM (17 years ago)

Author:

eugene

Message:

import updates from camera group (fits table outputs, fixes to persistent trails); modify import to avoid libpsf (fortran)

Location:

trunk/extsrc/gpcsw/gpcsrc/fits/burntool

Files:

• Makefile (modified) (3 diffs)
• burnparams.h (modified) (3 diffs)
• burntool.c (modified) (23 diffs)
• burntool.h (modified) (6 diffs)
• man/burntool.1 (modified) (4 diffs)
• persist_fits.c (added)
• persist_fits.h (added)
• persistfix.c (modified) (1 diff)
• persistio.c (modified) (1 diff)
• psfstats.c (modified) (3 diffs)
• trailfit.c (modified) (1 diff)

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

@@ -5 +5 @@
 CCDEFS+=$(VERSIONDEFS)
 CCLINK += -lm
-
 # $(EXECNAME): $(OBJS) libfh.a libpsf.a libpscoords.a
 $(EXECNAME): $(OBJS) libfh.a libpscoords.a
@@ -13 +12 @@
 # Dependencies by Make.Common $Revision: 2.17 $
 
-$(OBJ)/basement.o: basement.c
-
 $(OBJ)/burnfix.o: burnfix.c burntool.h burnparams.h
 
 $(OBJ)/burntool.o: burntool.c fh/fh.h fhreg/general.h \
   fhreg/macros.h fhreg/gpc_detector.h \
-  burntool.h burnparams.h
+  burntool.h burnparams.h persist_fits.h
 
 $(OBJ)/burnutils.o: burnutils.c burntool.h burnparams.h
 
 $(OBJ)/cheapfits.o: cheapfits.c
+
+$(OBJ)/persist_fits.o: persist_fits.c fh/fh.h burntool.h \
+  burnparams.h persist_fits.h
 
 $(OBJ)/persistfix.o: persistfix.c burntool.h burnparams.h
@@ -34 +34 @@
 $(OBJ)/psfstats.o: psfstats.c burntool.h psf/psf.h
 
+$(OBJ)/sort.o: sort.c
+
 $(OBJ)/stardetect.o: stardetect.c burntool.h burnparams.h
 

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

@@ -1 +1 @@
 /* A global buffer into which routines can accumulate detections */
+#ifndef _INCLUDED_burnparams_
+#define _INCLUDED_burnparams_
 
 /* '#define EXTERN' in just one file, burntool.c, to declare variables */
@@ -28 +30 @@
 EXTERN int MAX_READ_NOISE;      /* Maximum believable read noise (ADU) */
 EXTERN double MIN_EADU;         /* Minimum believable e/ADU */
+EXTERN int SAT4SURE;            /* Ignore pixels above for noise estimate */
+EXTERN double MIN_BLAST_PASS;   /* Allow blasted cells if they have BIG satfrac */
+EXTERN double MAX_BLAST_PASS;   /* But not if it's all wiped out! */
 
 EXTERN int BURN_THRESH  ;       /* Threshold for onset of burning */
@@ -48 +53 @@
 
 EXTERN double NEGLIGIBLE_TRAIL; /* Don't sweat less than this * sigma */
+EXTERN int EXPIRE_TRAIL_TIME;   /* Expire trails after this interval */
 
-
+#endif /* _INCLUDED_burnparams_ */

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

@@ -20 +20 @@
 #define EXTERN  /* Define EXTERN to declare variables in params.h */
 #include "burnparams.h"
+#include "persist_fits.h"
 
 int
@@ -28 +29 @@
    HeaderUnit ehu;
    const char* ifilename = "-";
-   const char* ofilename = NULL;
    char extname[FH_MAX_STRLEN+1], otaposn[FH_MAX_STRLEN+1];
    int cellmask[MAXCELL];
@@ -34 +34 @@
    int nextend, cellxy, cell, cellcode;
    int ext, update, restore, psfsize, psfavg;
-   int outfd;
    IMTYPE *buf;
-   const char *burnfile=NULL,  *persistfile=NULL;
+   const char *burnfile=NULL,  *persistfile=NULL, *persistfitsfile=NULL;
+   const char *deltablefitsfile=NULL;
    const char *psffile=NULL, *psfstatfile=NULL;
    CELL OTA[MAXCELL];   /* Cell structure for entire OTA */
 
-   if (argc == 1)  {
-     syntax(argv[0]);
-     exit(EXIT_SUCCESS);
-   }
-
    if(argc > 1 && strncmp(argv[1], "help", 4) == 0) {
-      syntax(argv[0]);
-      exit(EXIT_SUCCESS);
-   }
-
-   if(argc > 1 && strncmp(argv[1], "--help", 6) == 0) {
-      syntax(argv[0]);
-      exit(EXIT_SUCCESS);
-   }
-
-   if(argc > 1 && strncmp(argv[1], "-h", 6) == 0) {
       syntax(argv[0]);
       exit(EXIT_SUCCESS);
@@ -71 +56 @@
       fprintf(stderr, "\rerror: `%s' is not a multi-extension FITS\n",
               ifilename);
+#ifndef JT2DHACK
       exit(EXIT_FAILURE);
+#endif
    }
 
@@ -91 +78 @@
    MAX_READ_NOISE = 20;         /* Maximum believable read noise (ADU) */
    MIN_EADU = 0.3;              /* Minimum believable e/ADU */
+   SAT4SURE  = 60000;           /* Ignore pixels above for noise estimate */
+   MIN_BLAST_PASS = 0.1;        /* Allow blasted cells if they have BIG satfrac */
+   MAX_BLAST_PASS = 0.9;        /* But not if it's all wiped out! */
 
    BURN_THRESH  = 30000;        /* Threshold for onset of burning */
@@ -108 +98 @@
 
    NEGLIGIBLE_TRAIL = 0.5;      /* Don't sweat less than this * sigma */
-
+   EXPIRE_TRAIL_TIME = 2000;    /* Expire a persist after this [sec] */
 
 /* Parse the args */
@@ -155 +145 @@
          burnfile = argv[i] + 9;
 
-/* Output file for burn streaks */
+/* Output text file for burn streaks */
       } else if(strncmp(argv[i], "out=", 4) == 0) {     /* out=fname */
          burnfile = argv[i] + 4;
 
-# if (0)
-/* XXX disable for now: is not yet working */
-/* Alternate Output file for deburned image */
-      } else if(strncmp(argv[i], "outimage=", 9) == 0) {        /* outimage=fname */
-        ofilename = argv[i] + 9;
-# endif
-
-/* Input file for previous burn persistence streaks */
-      } else if(strncmp(argv[i], "trailin=", 8) == 0) { /* trailin=fname */
+/* Input text file for previous burn persistence streaks */
+      } else if(strncmp(argv[i], "infits=", 8) == 0) {  /* in=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 */
+         persistfitsfile = argv[i] + 12;
 
 /* Input file for previous burn persistence streaks */
       } else if(strncmp(argv[i], "in=", 3) == 0) {      /* in=fname */
          persistfile = argv[i] + 3;
+
+/* 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;
 
 /* Output file for PSF gallery */
@@ -244 +235 @@
       } else if(strncmp(argv[i], "thrpsf=", 7) == 0) {  /* thrpsf=thresh */
          if(sscanf(argv[i]+7, "%d", &PSF_THRESH) != 1) {
+            fprintf(stderr, "\rerror: cannot get value from `%s'\n", argv[i]);
+            exit(EXIT_FAILURE);
+         }
+
+      } else if(strncmp(argv[i], "expire=", 7) == 0) {  /* expire=nsec */
+         if(sscanf(argv[i]+7, "%d", &EXPIRE_TRAIL_TIME) != 1) {
             fprintf(stderr, "\rerror: cannot get value from `%s'\n", argv[i]);
             exit(EXIT_FAILURE);
@@ -295 +292 @@
          }
 
+      } else if(strncmp(argv[i], "deltables=", 8) == 0) { /* trailin=fname */
+         deltablefitsfile = argv[i] + 10;        
+
       } else if(strncmp(argv[i], "help", 4) == 0) {     /* help output */
          syntax(argv[0]);
@@ -306 +306 @@
    }
 
+   /* If we're told to remove the table from a FITS file, then
+    * do nothing else. */
+   if(deltablefitsfile) {
+     fprintf(stderr, 
+             "\rRemoving burn tables from %s, all other options ignored.\n",
+             ifilename);
+     if(persist_fits_remove_tables(ihu, deltablefitsfile) != FH_SUCCESS) {
+       exit(EXIT_FAILURE);
+     }
+     else {
+       exit(EXIT_SUCCESS);
+     }
+   }
+   
+   /* If there is no other persistence info supplied, try getting
+    * it from the input FITS file. */
    if(restore && persistfile == NULL) {
-      fprintf(stderr, "\rerror: must specify an input file for restore\n");
-      exit(EXIT_FAILURE);
+     if(persistfitsfile == NULL) persistfitsfile = ifilename;
    }
 
@@ -322 +337 @@
       if(persist_read(OTA, persistfile)) exit(EXIT_FAILURE);
    }
+   else if(persistfitsfile != NULL) {
+     if(persist_fits_read(OTA, persistfitsfile) != FH_SUCCESS) exit(EXIT_FAILURE);
+   }
 
 /* Which OTA is this??? */
@@ -333 +351 @@
    }
 
-   if (ofilename) {
-     outfd = creat (ofilename, S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | S_IROTH | S_IWOTH);
-     if (outfd == -1) {
-         fprintf(stderr, "\rerror: Cannot open %s for output\n", ofilename);
-         exit(EXIT_FAILURE);
-      }
-
-     if (fh_write(ihu, outfd)) {
-         fprintf(stderr, "\rerror: Trouble writing PHU to %s\n", ofilename);
-         exit(EXIT_FAILURE);
-      }
+   /* Warn that what you're doing might not be a good idea. */
+   {
+     fh_bool burn_applied;
+     
+     if(fh_get_bool(ihu, PHU_NAME_BURN_APPLIED, 
+                    &burn_applied) != FH_SUCCESS) {
+        if(VERBOSE > 0) {
+           fprintf(stderr, 
+                   "warning: Unable to determine whether burn correction "
+                   "already applied - unable to find %s in primary header\n",
+                   PHU_NAME_BURN_APPLIED);
+        }
+        burn_applied = FH_FALSE;
+     }
+     if(restore && (burn_applied == FH_FALSE)) {
+       fprintf(stderr, 
+               "warning: Restoring old burns, but header indicates no burns previously corrected.\n");             
+     }
+     else if (update && (burn_applied == FH_TRUE)) {
+       fprintf(stderr, 
+               "warning: Applying burn correction, but header indicates burns previously corrected.\n");             
+       }
    }
 
 /* Look at all the MEF's extensions */
+#ifndef JT2DHACK
    for (ext = 1; ext <= nextend; ext++) {
+#else
+   for (ext = 1; ext <= MAX(1,nextend); ext++) {
+#endif
       int naxis, naxis1, naxis2, naxis3;
       int prescan1, ovrscan1, ovrscan2, pontime;
@@ -353 +386 @@
       char xtension[FH_MAX_STRLEN + 1];
 
+#ifdef JT2DHACK
+      if(nextend == 0) {
+         ehu = ihu;
+         sprintf(extname, "%s", "xy00");
+      } else {
+#endif
       if (!(ehu = fh_ehu(ihu, ext)) || 
           fh_get_str(ehu, "EXTNAME", extname, sizeof(extname)) != FH_SUCCESS) {
@@ -377 +416 @@
          fprintf(stderr,
                  "warning: Skipping non-image extension %s\n", extname);
+#ifdef JT2DHACK
+      }
+#endif
 /*
  * %%% When reading from a pipe, does the data still need
@@ -448 +490 @@
       if (fh_get_BZERO(ehu, &bzero_d) == FH_SUCCESS) BZERO = bzero_d;
 
+      if(VERBOSE & VERB_NORM) {
+         printf("nx=%d ny=%d prex=%d postx=%d posty=%d BZERO=%.1f\n", 
+                naxis1, naxis2, prescan1, ovrscan1, ovrscan2, BZERO);
+      }
+
       buf = (IMTYPE*)malloc(naxis1*naxis2*naxis3*sizeof(short));
       if (fh_read_padded_image(ehu, fh_file_desc(ehu), buf,
@@ -471 +518 @@
 /* Get bias and sky levels */
          err = cell_stats(nx, ny, naxis1, naxis2, imbuf, OTA+cell);
+         if(err) {
+            if(VERBOSE > 0) {
+               fprintf(stderr, "logonly: error getting bias/sky/rms for cell %d, skipping...\n", cell);
+            }
+            continue;
+         }
+
+/* Does this cell get a pass because it's heavily blasted? */
+/* Turn it in to a relaxation factor for the "read noise" */
+         i = OTA[cell].satfrac > MIN_BLAST_PASS &&
+             OTA[cell].satfrac < MAX_BLAST_PASS &&
+             OTA[cell].sky < TRAIL_THRESH ? 30 : 1;
 
 /* Does this cell look kosher? */
          if(OTA[cell].rms*OTA[cell].rms > 
-            OTA[cell].sky/MIN_EADU + MAX_READ_NOISE*MAX_READ_NOISE) {
+            OTA[cell].sky/MIN_EADU + i*i*MAX_READ_NOISE*MAX_READ_NOISE) {
             if(VERBOSE > 0) {
                fprintf(stderr, "logonly: cell %d is unreasonably noisy, skipping...\n", cell);
+               fprintf(stderr, "logonly: bias = %d sky = %d rms = %d satfrac = %.2f\n", 
+                       OTA[cell].bias, OTA[cell].sky, OTA[cell].rms, OTA[cell].satfrac);
             }
             continue;
@@ -506 +567 @@
 
          } else {
-
 /* Restore the old burns */
             burn_restore(naxis1-ovrscan1, naxis2-ovrscan2, naxis1, 
@@ -514 +574 @@
 /* Write the corrected data back to the FITS. */
          if(update) {
-           if (ofilename) {
-             fh_ehu(ehu, 0);    /* Seek back to the start of data */
-
-             if (fh_write(ehu, outfd)) {
-               fprintf(stderr, "\rerror: Trouble writing EXT %d to %s\n", ext, ofilename);
-               exit(EXIT_FAILURE);
-             }
-
-             if (fh_write_padded_image(ehu, outfd, buf,
-                                       naxis1*naxis2*naxis3*sizeof(short),
-                                       FH_TYPESIZE_16) != FH_SUCCESS) {
+            fh_ehu(ehu, 0);     /* Seek back to the start of data */
+            if (fh_write_padded_image(ehu, fh_file_desc(ehu), buf,
+                                      naxis1*naxis2*naxis3*sizeof(short),
+                                      FH_TYPESIZE_16) != FH_SUCCESS) {
                fprintf(stderr, "\rerror: failed to re-write image data for extension `%s'.\n",
                        extname);
                free(buf);
                exit(EXIT_FAILURE);
-             }
-           } else {
-             fh_ehu(ehu, 0);    /* Seek back to the start of data */
-             if (fh_write_padded_image(ehu, fh_file_desc(ehu), buf,
-                                       naxis1*naxis2*naxis3*sizeof(short),
-                                       FH_TYPESIZE_16) != FH_SUCCESS) {
-               fprintf(stderr, "\rerror: failed to re-write image data for extension `%s'.\n",
-                       extname);
-               free(buf);
-               exit(EXIT_FAILURE);
-             }
-           }
-         }
-      }
-      
+            }
+         }
+      }
+
       free(buf);
    }
 
-   if (ofilename) {
-     if (close(outfd)) {
-       fprintf(stderr, "\rerror: trouble writing data to output file %s\n", ofilename);
-       exit(EXIT_FAILURE);
-     }
-   }
-      
 /* Dump out the postage stamp file */
    if(psffile != NULL) {
@@ -564 +599 @@
    }
 
+/* Write burn info to FITS file. */
+   if(update) persist_fits_write(OTA, ihu);
+
+   
+   if(restore) {    
+     /* 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_destroy(ihu);
 
 /* Write the persistence data for the next image */
    if(burnfile != NULL) persist_write(OTA, burnfile);
+
 
    exit(EXIT_SUCCESS);
@@ -613 +660 @@
                       CELL *cell)
 {
-   int i, j, k, n;
+   int i, j, k, n, nsat;
 
 /* Get bias stats */
@@ -624 +671 @@
 
 /* Get sky stats */
-   for(k=n=0; k<nx*ny; k+=((617*nx)/1000)) {
+   for(k=n=nsat=0; k<nx*ny; k+=((617*nx)/1000)) {
       i = k % nx;
       j = k / nx;
-      median_buf[n++] = data[i+NX*j];
-   }
+      if(data[i+NX*j] > SAT4SURE) {
+         nsat++;
+      } else if(data[i+NX*j] != NODATA) {
+         median_buf[n++] = data[i+NX*j];
+      }
+   }
+
+   if(n < 20) {         /* Better have hit at least 20! */
+      cell->sky = cell->rms = 0.0;
+      return(-1);
+   }
+   cell->satfrac = ((double)nsat) / n;
+
 /* First pass at sky and quartile */
    cell->sky = int_median(n, median_buf);
@@ -656 +714 @@
    printf(" restore={t|f}  Restore the input MEF by adding input fits?\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");
+   printf("                   in table extensions).  If both this and the 'in' input\n");
+   printf("                   file option are specified, then 'in' takes precedence.\n");
    printf(" out=fname      Output file for burn streaks\n");
-   printf(" outimage=imname Output file for MEF image\n");
+
 //   printf(" trailin=fname  Input file for previous burn persistence streaks\n");
+//   printf(" trailinfits=fname  Input FITS file for previous burn persistence streaks\n");
 //   printf(" trailout=fname Output file for burn streaks\n");
+   printf(" deltables=fname Copy mef_file to new FITS file 'fname' with burn streak\n");
+   printf("                   tables removed. NOTE: if specified, all other options\n");
+   printf("                   will be ignored!\n");
    printf(" psf=fname      Output file for PSF FITS stamp gallery\n");
    printf(" psfstat=fname  Output file for PSF statistics listing\n");
@@ -676 +741 @@
    printf(" rmask=X        Diameter growth factor of burned spots\n");
    printf(" bmask=X        Box size growth of burn/star boxes\n");
+   printf(" expire=N       Retire a blasted burn after N seconds\n");
    printf(" quiet={t|f}    Quiet?\n");
    printf(" verbose=N      Set verbosity bits:\n");

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

@@ +1 @@
+#ifndef _INCLUDED_burntool_
+#define _INCLUDED_burntool_
+
 /* Burn correction routines:
    -------------------------  
@@ -63 +66 @@
 
 #define MAXCELL 64              /* Max cells in an OTA */
-#define MAXSIZE 700             /* Maximum vertical cell size */
+#define MAXSIZE 2048            /* Maximum vertical cell size */
 
 #define STAR_RADIUS  4          /* Radius over which a star ctr must be max */
@@ -87 +90 @@
 #define BURN_PWR  1             /* Power law */
 #define BURN_EXP  2             /* Exponential */
-#define PSF_STAR  3             /* Unfitted: good psf star */
+#define BURN_BLASTED 3          /* Blasted top to bottom: flag only for IPP */
+#define PSF_STAR  9             /* Unfitted: good psf star */
 
 /* Fit error codes */
@@ -93 +97 @@
 #define FIT_TOP_ERROR  2        /* Saturation extends to top: no points */
 #define FIT_SLOPE_ERROR  3      /* Unreasonable fit */
+#define FIT_EXPIRED 9           /* Don't carry any more as persistent */
 
 /* Fit parameters */
@@ -152 +157 @@
       int rms;                  /* RMS in the sky */
       int time;                 /* PON time of this cell */
+      double satfrac;           /* Fraction of SAT4SURE saturated pixels */
       int nburn;                /* Number of trails left by sat stars */
       OBJBOX *burn;             /* Stars which we think left a trail */
@@ -212 +218 @@
 int write_2ddata(int nx, int ny, int *ntot, IMTYPE *data, int fd);
 int write_3dend(int *ntot, int fd);
+
+#endif /* _INCLUDED_burntool_ */

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

@@ -53 +53 @@
         power laws; downward burns as exponentials.
         
+        Burntool also identifies really blasted areas which are saturated from
+        top to bottom.  These cannot be fitted, but are carried along for
+        EXPIRE_TRAIL_TIME seconds (default 2000) in the persistence list.  The
+        idea is to give IPP a clue that something bad happened, and of course
+        the persistence info tells when it happened.  This expiration time
+        for carrying the region can be set with "expire=N".
+
         Burntool finishes by rewriting the modified MEF and all of the fits in
         an output file ("out=file_name").  If desired, burntool can be run
@@ -112 +119 @@
            npsf=3          Number of contributing PSF stars
            fwhm=5.38       PSF FWHM [pix] (0.0 if fit fails)
-           m2=11.08        Second moment of distribution [pix^2] 
-           qp=5.154        Plus quadrupole qxx-qyy [pix^2]       
-           qc=0.560        Cross quadrupole 2*qxy [pix^2]        
-           qt=-5.201       Tangential quadrupole [pix^2]          
+           fmed=5.76       FWHM averaged over psfavg (0.0 if fit fails)
+           m2=11.08        Second moment of distribution [pix^2] (-99.99 if fails)
+           qp=5.154        Plus quadrupole qxx-qyy [pix^2] (-99.99 if fails)
+           qc=0.560        Cross quadrupole 2*qxy [pix^2] (-99.99 if fails)
+           qt=-5.201       Tangential quadrupole [pix^2] (-99.99 if fails)
                                 q+ * cos(2*phi) + qx * sin(2*phi)
+           qpm=5.154       Plus quadrupole averaged over psfavg
+           qcm=0.560       Cross quadrupole averaged over psfavg
+           qtm=-5.201      Tangential averaged over psfavg
 
         Helpful utilities include:
@@ -148 +159 @@
                 Input file for previous burn persistence streaks
 
+        infits=fname       
+                Input FITS file for previous burn persistence streaks
+                (burn streaks stored in table extensions). If both
+                this and the 'in' input file option are specified, then
+                'in' takes precedence.
+
         out=fname      
                 Output file for burn streaks
+
+        deltables=fname
+                Copy mef_file to new FITS file 'fname' with burn streak 
+                tables removed. NOTE: if specified, all other options 
+                will be ignored!
 
         psf=fname      
@@ -199 +221 @@
         bmask=X        
                 Box size growth of burn/star (BMASK_GROW)
+
+        expire=N
+                Retire a blasted burn after N seconds
 
         quiet={t|f}    

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

@@ -27 +27 @@
 /* Fix up all the persistence streaks */
    for(k=0; k<cell->npersist; k++) {
+
+/* Is this just a blasted area being carried for IPP? */
+      if((cell->persist)[k].func == BURN_BLASTED) {
+         if(cell->time - (cell->persist)[k].time > EXPIRE_TRAIL_TIME) {
+            (cell->persist)[k].fiterr = FIT_EXPIRED;
+         }
+         continue;
+      }
 
 /* Fit the trail */

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

@@ -233 +233 @@
    fprintf(fp,  "# Cell: %d  sky= %d   rms= %d   bias= %d\n", 
            cell[0].cell, cell[0].sky, cell[0].rms, cell[0].bias);
-   fprintf(fp, "#Cell time    cx  cy  max   sx  sy   ex  ey  y0m y0p y1m y1p x0m x0p x1m x1p     slope    zero\n");
+   fprintf(fp, "#Cell time    cx  cy  max   y0   sx  sy   ex  ey  y0m y0p y1m y1p x0m x0p x1m x1p F up    slope nfit sxf exf\n");
 
    for(j=0; j<MAXCELL; j++) {

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

@@ -1 @@
-/* psfstats.c - calculate some statistics about a PSF */
+/* dummy version of psfstats.c - disable use of the psf */
 #include <stdio.h>
 #include <stdlib.h>
@@ -12 +12 @@
 #include "psf/psf.h"
 
-//#define PSFSTATS              /* ignore psf stats so we can skip fortan */
-//#define TEST          /* Dump all the PSF info? */
-
-static unsigned short int *psfbuf=NULL;
-static int npsfbuf=0;
-#define PHONYSKY 1000
-
 /****************************************************************/
 /* psf_stats(): Tell us about this PSF star */
@@ -24 +17 @@
                      double *fwhm, double *q)
 {
-
-#ifdef PSFSTATS
-
-   int i, j, r, r2, err, i0, j0;
-   double d, qxx, qxy, qyy, mx, my, sum;
-   PSF_PARAM psfout;      /* Results of fit */
-   PSF_EXTRA psfextra;    /* Extra results from fit */
-   double pi=4*atan(1.0);
-
-   if(nx*ny > npsfbuf) {
-      if(psfbuf != NULL) free(psfbuf);
-      psfbuf = (ushort *)calloc(nx*ny, sizeof(ushort));
-      npsfbuf = nx*ny;
-   }
-   for(i=0; i<nx*ny; i++) {
-      j = ((int)data[i]) - bias + PHONYSKY;
-      if(j >= 0 && j < 65535) psfbuf[i] = j;
-      else psfbuf[i] = 0;
-   }
-
-/* 2D PSF fit */
-   err = psf(nx, ny, psfbuf, PSF_2DIM, NULL, &psfout, &psfextra);
-
-#ifdef TEST
-   printf("err = %d\n", err);
-#endif
-
-   if(!err) {
-#ifdef TEST
-      printf("err = %d\n", err);
-      printf("bias = %d\n", bias);
-      printf("ix, iy = %d %d\n", psfout.ix, psfout.iy);
-      printf("x0, y0 = %.2f %.2f\n", psfout.x0, psfout.y0);
-      printf("fwhm = %.2f\n", psfout.fwhm);
-      printf("peak, bckgnd = %.1f %.1f\n", psfout.peak, psfout.bkgnd);
-      printf("flux, S/N = %.1f %.2f\n", psfout.flux, psfout.sn);
-      printf("weight = %.2f\n", psfout.weight);
-      printf("FWmaj, FWmin = %.2f %.2f\n", psfextra.majfw, psfextra.minfw);
-      printf("theta = %.1f\n", psfextra.thfw*180/pi);
-      printf("wpeak, wbkgnd = %.1f %.1f\n", psfextra.wpeak, psfextra.wbkgnd);
-      printf("dflux, dbckgnd = %.1f %.1f\n", psfextra.dflux, psfextra.dbkgnd);
-      printf("rmsbkgnd = %.1f\n", psfextra.rmsbkgnd);
-#endif
-      fwhm[0] = psfextra.majfw;
-      fwhm[1] = psfextra.minfw;
-      fwhm[2] = psfextra.thfw;
-
-/* Moments */
-      i0 = MIN(psfout.ix, nx-psfout.ix);
-      j0 = MIN(psfout.iy, ny-psfout.iy);
-   } else {
-      fwhm[0] = fwhm[1] = fwhm[2] = 0.0;
-      i0 = nx/2;
-      j0 = ny/2;
-   }
-
-   r = MIN(i0, j0) - 1;
-   sum = mx = my = qxx = qxy = qyy = 0.0;
-   for(j=j0-r; j<j0+r; j++) {
-      for(i=i0-r; i<i0+r; i++) {
-         r2 = (i-i0)*(i-i0) + (j-j0)*(j-j0);
-         if(r2 > r*r) continue;
-         d = (double)data[i+j*nx] - bias;
-         sum += d;
-         mx += (i+0.5) * d;
-         my += (j+0.5) * d;
-         qxx += (i+0.5) * (i+0.5) * d;
-         qxy += (i+0.5) * (j+0.5) * d;
-         qyy += (j+0.5) * (j+0.5) * d;
-      }
-   }
-   if(sum > 0) {
-      mx /= sum;
-      my /= sum;
-      qxx = qxx / sum - mx*mx;
-      qxy = qxy / sum - mx*my;
-      qyy = qyy / sum - my*my;
-   } else {
-      q[0] = q[1] = q[2] = 0.0;
-   }
-
-   q[0] = (qxx + qyy) / 2;
-   q[1] = qxx - qyy;
-   q[2] = 2 * qxy;
-
-/* phi is the angular position of the OTA around the center [rad] */
-//   qtang = (qxx - qyy) * cos(2*phi) + 2 * qxy * sin(2*phi);
-
-#ifdef TEST
-   printf("Sum= %.1f mx= %.2f my= %.2f qxx= %.4f qxy= %.4f qyy= %.4f\n", 
-          sum, mx, my, qxx, qxy, qyy);
-#endif
-
-# else
-   printf("psfstats ignored in this build\n");
-   return (1);
-#endif /* PSFSTATS */
-   return(0);
+   return(1);
 }

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

@@ -148 +148 @@
    }
 
+/* Burn extends all the way to the top */
    if(up && y1 >= ny-1) {
       box->slope = -1.0;
-      box->nfit = xe - xs + 1;
       box->fiterr = FIT_TOP_ERROR;
+/* No fit here because the burn is to top, but the persist may be fitable */
+      if(box->sy > FIT_EDGE) {
+         box->nfit = xe - xs + 1;
+/* However, if burn extends to bottom, persist will not be able to fit */
+      } else {
+         box->nfit = 1;         /* Just enough to keep it alive */
+         box->func = BURN_BLASTED;
+      }
       return(0);
    }