Index: trunk/extsrc/gpcsw/gpcsrc/fits/burntool/Makefile =================================================================== --- trunk/extsrc/gpcsw/gpcsrc/fits/burntool/Makefile (revision 23564) +++ trunk/extsrc/gpcsw/gpcsrc/fits/burntool/Makefile (revision 23924) @@ -5,5 +5,4 @@ CCDEFS+=$(VERSIONDEFS) CCLINK += -lm - # $(EXECNAME): $(OBJS) libfh.a libpsf.a libpscoords.a $(EXECNAME): $(OBJS) libfh.a libpscoords.a @@ -13,15 +12,16 @@ # 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,4 +34,6 @@ $(OBJ)/psfstats.o: psfstats.c burntool.h psf/psf.h +$(OBJ)/sort.o: sort.c + $(OBJ)/stardetect.o: stardetect.c burntool.h burnparams.h Index: trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burnparams.h =================================================================== --- trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burnparams.h (revision 23564) +++ trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burnparams.h (revision 23924) @@ -1,3 +1,5 @@ /* 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,4 +30,7 @@ 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,4 +53,5 @@ EXTERN double NEGLIGIBLE_TRAIL; /* Don't sweat less than this * sigma */ +EXTERN int EXPIRE_TRAIL_TIME; /* Expire trails after this interval */ - +#endif /* _INCLUDED_burnparams_ */ Index: trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.c =================================================================== --- trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.c (revision 23564) +++ trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.c (revision 23924) @@ -20,4 +20,5 @@ #define EXTERN /* Define EXTERN to declare variables in params.h */ #include "burnparams.h" +#include "persist_fits.h" int @@ -28,5 +29,4 @@ HeaderUnit ehu; const char* ifilename = "-"; - const char* ofilename = NULL; char extname[FH_MAX_STRLEN+1], otaposn[FH_MAX_STRLEN+1]; int cellmask[MAXCELL]; @@ -34,26 +34,11 @@ 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,5 +56,7 @@ fprintf(stderr, "\rerror: `%s' is not a multi-extension FITS\n", ifilename); +#ifndef JT2DHACK exit(EXIT_FAILURE); +#endif } @@ -91,4 +78,7 @@ 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,5 +98,5 @@ 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,22 +145,23 @@ 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,4 +235,10 @@ } 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,4 +292,7 @@ } + } 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,7 +306,22 @@ } + /* 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,4 +337,7 @@ 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,19 +351,34 @@ } - 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,4 +386,10 @@ 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,4 +416,7 @@ fprintf(stderr, "warning: Skipping non-image extension %s\n", extname); +#ifdef JT2DHACK + } +#endif /* * %%% When reading from a pipe, does the data still need @@ -448,4 +490,9 @@ 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,10 +518,24 @@ /* 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,5 +567,4 @@ } else { - /* Restore the old burns */ burn_restore(naxis1-ovrscan1, naxis2-ovrscan2, naxis1, @@ -514,44 +574,19 @@ /* 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,8 +599,20 @@ } +/* 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,5 +660,5 @@ CELL *cell) { - int i, j, k, n; + int i, j, k, n, nsat; /* Get bias stats */ @@ -624,9 +671,20 @@ /* Get sky stats */ - for(k=n=0; k 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,8 +714,15 @@ 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,4 +741,5 @@ 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"); Index: trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.h =================================================================== --- trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.h (revision 23564) +++ trunk/extsrc/gpcsw/gpcsrc/fits/burntool/burntool.h (revision 23924) @@ -1,2 +1,5 @@ +#ifndef _INCLUDED_burntool_ +#define _INCLUDED_burntool_ + /* Burn correction routines: ------------------------- @@ -63,5 +66,5 @@ #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,5 +90,6 @@ #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,4 +97,5 @@ #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,4 +157,5 @@ 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,2 +218,4 @@ int write_2ddata(int nx, int ny, int *ntot, IMTYPE *data, int fd); int write_3dend(int *ntot, int fd); + +#endif /* _INCLUDED_burntool_ */ Index: trunk/extsrc/gpcsw/gpcsrc/fits/burntool/man/burntool.1 =================================================================== --- trunk/extsrc/gpcsw/gpcsrc/fits/burntool/man/burntool.1 (revision 23564) +++ trunk/extsrc/gpcsw/gpcsrc/fits/burntool/man/burntool.1 (revision 23924) @@ -53,4 +53,11 @@ 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,9 +119,13 @@ 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,6 +159,17 @@ 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,4 +221,7 @@ bmask=X Box size growth of burn/star (BMASK_GROW) + + expire=N + Retire a blasted burn after N seconds quiet={t|f} Index: trunk/extsrc/gpcsw/gpcsrc/fits/burntool/persist_fits.c =================================================================== --- trunk/extsrc/gpcsw/gpcsrc/fits/burntool/persist_fits.c (revision 23924) +++ trunk/extsrc/gpcsw/gpcsrc/fits/burntool/persist_fits.c (revision 23924) @@ -0,0 +1,1224 @@ +/* -*- c-file-style: "Ellemtel" -*- + * + * persist_fits.c - read and write persistence info to FITS tables. + * + */ + +#include +#include +#include +#include +#include +#include +#include + +#include "fh/fh.h" + +#include "burntool.h" +#include "burnparams.h" +#include "persist_fits.h" + +/* Default names for burn info table extensions. We'll use keyword information + * when reading to determine the names of the extensions, but we'll always + * write them with these names. */ +#define DEFAULT_EXTNAME_AREA_TABLE "burntool_areas" +#define DEFAULT_EXTNAME_FIT_TABLE "burntool_fits" + +/* Index of all the columns of the area table. */ +typedef enum +{ + AREA_TABLE_COL_CELL, + AREA_TABLE_COL_TIME, + AREA_TABLE_COL_CX, + AREA_TABLE_COL_CY, + AREA_TABLE_COL_MAX, + AREA_TABLE_COL_Y0, + AREA_TABLE_COL_SX, + AREA_TABLE_COL_SY, + AREA_TABLE_COL_EX, + AREA_TABLE_COL_EY, + AREA_TABLE_COL_Y0M, + AREA_TABLE_COL_Y0P, + AREA_TABLE_COL_Y1M, + AREA_TABLE_COL_Y1P, + AREA_TABLE_COL_X0M, + AREA_TABLE_COL_X0P, + AREA_TABLE_COL_X1M, + AREA_TABLE_COL_X1P, + AREA_TABLE_COL_FUNC, + AREA_TABLE_COL_UP, + AREA_TABLE_COL_SLOPE, + AREA_TABLE_COL_NFIT, + AREA_TABLE_COL_SXFIT, + AREA_TABLE_COL_EXFIT, + + /* Add new columns above this line. */ + AREA_TABLE_NUM_COLS +} AREA_TABLE_COL_T; + +/* All of the columns of the area table. */ +static fhTableCol +area_table_cols[] = +{ +/* Name Description Units Format Width Dec places */ + { "cell", "Cell that correction area is in", "", FH_TABLE_FORMAT_INT, 2, 0 }, + { "time", "PON time when area created", "secs", FH_TABLE_FORMAT_INT, 10, 0 }, + { "cx", "Center x (position of max)", "pixels", FH_TABLE_FORMAT_INT, 3, 0 }, + { "cy", "Center y (position of max)", "pixels", FH_TABLE_FORMAT_INT, 3, 0 }, + { "max", "Max data value (above sky)", "ADU", FH_TABLE_FORMAT_INT, 5, 0 }, + { "y0", "y origin for fit (stamp sy for stars)", "pixels", FH_TABLE_FORMAT_INT, 3, 0 }, + { "sx", "Left corner of area", "pixels", FH_TABLE_FORMAT_INT, 3, 0 }, + { "sy", "Bottom corner of area", "pixels", FH_TABLE_FORMAT_INT, 3, 0 }, + { "ex", "Right corner of area", "pixels", FH_TABLE_FORMAT_INT, 3, 0 }, + { "ey", "Top corner of area", "pixels", FH_TABLE_FORMAT_INT, 3, 0 }, + { "y0m", "Min y value at sx", "pixels", FH_TABLE_FORMAT_INT, 5, 0 }, + { "y0p", "Max y value at sx", "pixels", FH_TABLE_FORMAT_INT, 5, 0 }, + { "y1m", "Min y value at ex", "pixels", FH_TABLE_FORMAT_INT, 5, 0 }, + { "y1p", "Max y value at ex", "pixels", FH_TABLE_FORMAT_INT, 5, 0 }, + { "x0m", "Min x value at sy", "pixels", FH_TABLE_FORMAT_INT, 5, 0 }, + { "x0p", "Max x value at sy", "pixels", FH_TABLE_FORMAT_INT, 5, 0 }, + { "x1m", "Min x value at ey", "pixels", FH_TABLE_FORMAT_INT, 5, 0 }, + { "x1p", "Max x value at ey", "pixels", FH_TABLE_FORMAT_INT, 5, 0 }, + { "func", "What are we going to do about it?", "", FH_TABLE_FORMAT_INT, 3, 0 }, + { "up", "Trails up or down", "1=up 0=down", FH_TABLE_FORMAT_INT,1, 0 }, + { "slope", "Slope of fit", "", FH_TABLE_FORMAT_DOUBLE, 9, 6 }, + { "nfit", "Number of columns corrected", "", FH_TABLE_FORMAT_INT, 3, 0 }, + { "sxfit", "Starting column for fit", "pixels", FH_TABLE_FORMAT_INT, 3, 0 }, + { "exfit", "Ending column for fit", "pixels", FH_TABLE_FORMAT_INT, 3, 0 }, +}; + +/* Index of all of the columns in the fit table. */ +typedef enum +{ + FIT_TABLE_COL_CELL, + FIT_TABLE_COL_CX, + FIT_TABLE_COL_CY, + FIT_TABLE_COL_XFIT, + FIT_TABLE_COL_YFIT, + FIT_TABLE_COL_ZERO, + + /* Add new columns above this line. */ + FIT_TABLE_NUM_COLS +} FIT_TABLE_COL_T; + +/* All of the columns in the fit table. */ +static fhTableCol +fit_table_cols[] = +{ +/* Name Description Units Format Width Dec places */ + { "cell", "Cell that correction area is in", "", FH_TABLE_FORMAT_INT, 2, 0 }, + { "cx", "Center x (position of max)", "pixels", FH_TABLE_FORMAT_INT, 3, 0 }, + { "cy", "Center y (position of max)", "pixels", FH_TABLE_FORMAT_INT, 3, 0 }, + { "xfit", "x of each value of the start of correction", "pixels", FH_TABLE_FORMAT_INT,3, 0 }, + { "yfit", "x of each value of the start of correction", "pixels", FH_TABLE_FORMAT_INT,3, 0 }, + { "zero", "Zero of this fit", "", FH_TABLE_FORMAT_DOUBLE, 8, 4 }, +}; + +/* Structure describing area table. Stuff like the number of + * rows, the overall size, etc. are calculated at runtime. */ +static fhTable +area_table = +{ + .extname = DEFAULT_EXTNAME_AREA_TABLE, + .num_cols = AREA_TABLE_NUM_COLS, + .num_rows = 0, /* Don't know yet - need to rip through the data to figure this out. */ + .cols = area_table_cols, +}; + +/* This is the fit table, in all its glory. */ +static fhTable +fit_table = +{ + .extname = DEFAULT_EXTNAME_FIT_TABLE, + .num_cols = FIT_TABLE_NUM_COLS, + .num_rows = 0, /* Don't know yet - need to rip through the data to figure this out. */ + .cols = fit_table_cols, +}; + +/* + * Finds existing fit and area tables in a FITS file. + * + * Given the primary header of a FITS file, this function + * tries to find the area and fit table extensions. It does this + * by first looking at the header to see if the extension names + * of the tables are listed, and if so, looks for those by name. + * If that information is not available for some reason, the default + * names are used instead. + * + * Inputs: + * phu - primary header of a FITS file. + * + * Outputs: + * area_hu - header of area table stored here if found, NULL otherwise. + * area_name - if non-NULL, this is a buffer where the name of the + * area table extension will be copied. + * fit_hu - header of fit table stored here if found, NULL otherwise. + * fit_name - if non-NULL, this is a buffer where the name of the + * fit table extension will be copied. + * + * This function returns nothing. + */ +static void +find_existing_tables(HeaderUnit phu, + HeaderUnit * area_hu, char * area_name, + HeaderUnit * fit_hu, char * fit_name) +{ + char area_extname[FH_MAX_STRLEN+1]; + char fit_extname[FH_MAX_STRLEN+1]; + + /* Does the FITS file provide the names of the burn tables? + * If so, use those, otherwise revert to the defaults. */ + if(fh_get_str(phu, PHU_NAME_BURN_AREA, area_extname, sizeof(area_extname)) != FH_SUCCESS) + { + snprintf(area_extname, FH_MAX_STRLEN+1, DEFAULT_EXTNAME_AREA_TABLE); + } + if(fh_get_str(phu, PHU_NAME_BURN_FIT, fit_extname, sizeof(fit_extname)) != FH_SUCCESS) + { + snprintf(fit_extname, FH_MAX_STRLEN+1, DEFAULT_EXTNAME_FIT_TABLE); + } + + *area_hu = fh_ehu_by_extname(phu, area_extname); + if(area_name) + { + if(*area_hu) strcpy(area_name, area_extname); + else area_extname[0] = '\0'; + } + *fit_hu = fh_ehu_by_extname(phu, fit_extname); + if(fit_name) + { + if(*fit_hu) strcpy(fit_name, fit_extname); + else fit_name[0] = '\0'; + } +} + +/* + * Determines how many rows should be in the area and fit tables. + * + * Given all of the burn information for the whole OTA, this function + * figures out how many rows are needed in the area and fit tables to + * accomodate all of that info. + * + * Inputs: + * cell - array of 64 CELL structures, one for each cell in the OTA. + * + * Outputs: + * num_areas_found - the number of burn areas for the OTA stored here. + * num_fits_found - the number of fits across all of the burn + * areas for the OTA stored here. + * + * This function returns FH_SUCCESS if the row was successfully + * written. + */ +static void +calculate_rows(CELL *cell, int * num_areas_found, int * num_fits_found) +{ + int j, k; + int num_areas = 0; + int num_fits = 0; + + for(j=0; jcell)) != FH_SUCCESS) || + (fh_table_write_value(table, data, row, AREA_TABLE_COL_TIME, &(area->time)) != FH_SUCCESS) || + (fh_table_write_value(table, data, row, AREA_TABLE_COL_CX, &(area->cx)) != FH_SUCCESS) || + (fh_table_write_value(table, data, row, AREA_TABLE_COL_CY , &(area->cy)) != FH_SUCCESS) || + (fh_table_write_value(table, data, row, AREA_TABLE_COL_MAX, &(area->max)) != FH_SUCCESS) || + (fh_table_write_value(table, data, row, AREA_TABLE_COL_Y0, &(area->y0)) != FH_SUCCESS) || + (fh_table_write_value(table, data, row, AREA_TABLE_COL_SX, &(area->sx)) != FH_SUCCESS) || + (fh_table_write_value(table, data, row, AREA_TABLE_COL_SY, &(area->sy)) != FH_SUCCESS) || + (fh_table_write_value(table, data, row, AREA_TABLE_COL_EX, &(area->ex)) != FH_SUCCESS) || + (fh_table_write_value(table, data, row, AREA_TABLE_COL_EY, &(area->ey)) != FH_SUCCESS) || + (fh_table_write_value(table, data, row, AREA_TABLE_COL_Y0M, &(area->y0m)) != FH_SUCCESS) || + (fh_table_write_value(table, data, row, AREA_TABLE_COL_Y0P, &(area->y0p)) != FH_SUCCESS) || + (fh_table_write_value(table, data, row, AREA_TABLE_COL_Y1M, &(area->y1m)) != FH_SUCCESS) || + (fh_table_write_value(table, data, row, AREA_TABLE_COL_Y1P, &(area->y1p)) != FH_SUCCESS) || + (fh_table_write_value(table, data, row, AREA_TABLE_COL_X0M, &(area->x0m)) != FH_SUCCESS) || + (fh_table_write_value(table, data, row, AREA_TABLE_COL_X0P, &(area->x0p)) != FH_SUCCESS) || + (fh_table_write_value(table, data, row, AREA_TABLE_COL_X1M, &(area->x1m)) != FH_SUCCESS) || + (fh_table_write_value(table, data, row, AREA_TABLE_COL_X1P, &(area->x1p)) != FH_SUCCESS) || + (fh_table_write_value(table, data, row, AREA_TABLE_COL_FUNC, &(area->func)) != FH_SUCCESS) || + (fh_table_write_value(table, data, row, AREA_TABLE_COL_UP, &(area->up)) != FH_SUCCESS) || + (fh_table_write_value(table, data, row, AREA_TABLE_COL_SLOPE, &(area->slope)) != FH_SUCCESS) || + (fh_table_write_value(table, data, row, AREA_TABLE_COL_NFIT, &(area->nfit)) != FH_SUCCESS) || + (fh_table_write_value(table, data, row, AREA_TABLE_COL_SXFIT, &(area->sxfit)) != FH_SUCCESS) || + (fh_table_write_value(table, data, row, AREA_TABLE_COL_EXFIT, &(area->exfit)) != FH_SUCCESS)) + { + fprintf(stderr, "error: Error writing data to row %d of area table.\n", row); + return FH_BAD_VALUE; + } + + return FH_SUCCESS; +} + +/* + * Writes area data to a FITS file. + * + * This function iterates through all of the burn areas in all cells + * and populates the body of the area table in the provided buffer. + * + * Inputs: + * hu - header of fit table extension. + * data - points to buffer for table body where data will be written. + * table - describes structure of the table. + * cell - array of 64 CELL structures for all cells of an OTA. Burn + * data from these structures will be written to the table. + * + * Outputs: + * None. + * + * This function returns FH_SUCCESS if the table was successfully + * written. + */ +static fh_result +write_area_table(HeaderUnit hu, void * data, fhTable * table, CELL * cell) +{ + int j, k; + int row=0; + fh_result result; + + for(j=0; j MAXCELL)) + { + fprintf(stderr, "error: illegal cell %d in area table row %d\n", + cell_num, row); + boxbuf[row].cell = -1; + boxbuf[row].nfit = 0; + return FH_BAD_VALUE; + } + + boxbuf[row].cell = cell_num; + + if((fh_table_read_value(table, data, row, AREA_TABLE_COL_TIME, &(boxbuf[row].time)) != FH_SUCCESS) || + (fh_table_read_value(table, data, row, AREA_TABLE_COL_CX, &(boxbuf[row].cx)) != FH_SUCCESS) || + (fh_table_read_value(table, data, row, AREA_TABLE_COL_CY, &(boxbuf[row].cy)) != FH_SUCCESS) || + (fh_table_read_value(table, data, row, AREA_TABLE_COL_MAX, &(boxbuf[row].max)) != FH_SUCCESS) || + (fh_table_read_value(table, data, row, AREA_TABLE_COL_Y0, &(boxbuf[row].y0)) != FH_SUCCESS) || + (fh_table_read_value(table, data, row, AREA_TABLE_COL_SX, &(boxbuf[row].sx)) != FH_SUCCESS) || + (fh_table_read_value(table, data, row, AREA_TABLE_COL_SY, &(boxbuf[row].sy)) != FH_SUCCESS) || + (fh_table_read_value(table, data, row, AREA_TABLE_COL_EX, &(boxbuf[row].ex)) != FH_SUCCESS) || + (fh_table_read_value(table, data, row, AREA_TABLE_COL_EY, &(boxbuf[row].ey)) != FH_SUCCESS) || + (fh_table_read_value(table, data, row, AREA_TABLE_COL_Y0M, &(boxbuf[row].y0m)) != FH_SUCCESS) || + (fh_table_read_value(table, data, row, AREA_TABLE_COL_Y0P, &(boxbuf[row].y0p)) != FH_SUCCESS) || + (fh_table_read_value(table, data, row, AREA_TABLE_COL_Y1M, &(boxbuf[row].y1m)) != FH_SUCCESS) || + (fh_table_read_value(table, data, row, AREA_TABLE_COL_Y1P, &(boxbuf[row].y1p)) != FH_SUCCESS) || + (fh_table_read_value(table, data, row, AREA_TABLE_COL_X0M, &(boxbuf[row].x0m)) != FH_SUCCESS) || + (fh_table_read_value(table, data, row, AREA_TABLE_COL_X0P, &(boxbuf[row].x0p)) != FH_SUCCESS) || + (fh_table_read_value(table, data, row, AREA_TABLE_COL_X1M, &(boxbuf[row].x1m)) != FH_SUCCESS) || + (fh_table_read_value(table, data, row, AREA_TABLE_COL_X1P, &(boxbuf[row].x1p)) != FH_SUCCESS) || + (fh_table_read_value(table, data, row, AREA_TABLE_COL_FUNC, &(boxbuf[row].func)) != FH_SUCCESS) || + (fh_table_read_value(table, data, row, AREA_TABLE_COL_UP, &(boxbuf[row].up)) != FH_SUCCESS) || + (fh_table_read_value(table, data, row, AREA_TABLE_COL_SLOPE, &(boxbuf[row].slope)) != FH_SUCCESS) || + (fh_table_read_value(table, data, row, AREA_TABLE_COL_NFIT, &(boxbuf[row].nfit)) != FH_SUCCESS) || + (fh_table_read_value(table, data, row, AREA_TABLE_COL_SXFIT, &(boxbuf[row].sxfit)) != FH_SUCCESS) || + (fh_table_read_value(table, data, row, AREA_TABLE_COL_EXFIT, &(boxbuf[row].exfit)) != FH_SUCCESS)) + { + fprintf(stderr, + "error: Error reading values from row %d of burn area table\n", + row); + boxbuf[row].nfit = 0; + return FH_BAD_VALUE; + } + + /* Make space for fits (to be read later). */ + if(boxbuf[row].nfit > 0) + { + /* This shouldn't happen, but be nice and don't leak. */ + if(boxbuf[row].zero) free(boxbuf[row].zero); + if(boxbuf[row].xfit) free(boxbuf[row].xfit); + if(boxbuf[row].yfit) free(boxbuf[row].yfit); + + boxbuf[row].zero = (double *)calloc(boxbuf[row].nfit, sizeof(double)); + boxbuf[row].xfit = (int *)calloc(boxbuf[row].nfit, sizeof(int)); + boxbuf[row].yfit = (int *)calloc(boxbuf[row].nfit, sizeof(int)); + boxbuf[row].fiterr = 0; + } + + return FH_SUCCESS; +} + +/* + * Reads the area table from a FITS file. + * + * This function fills out boxbuf[] with all of the areas in the + * given table. + * + * Inputs: + * hu - header of the table extension for the area table. + * table - describes structure of the table. + * + * Outputs: + * None. + * + * This function returns FH_SUCCESS if the tables were successfully + * appended to the FITS file. + */ +static fh_result +read_area_table(HeaderUnit hu, fhTable * table) +{ + fh_result result = FH_INVALID; + int i; + int row_chars, num_cols, num_rows; + void * data; + + if((fh_get_int(hu, "NAXIS1", &row_chars) != FH_SUCCESS) || + (fh_get_int(hu, "NAXIS2", &num_rows) != FH_SUCCESS) || + (fh_get_int(hu, "TFIELDS", &num_cols) != FH_SUCCESS)) + { + fprintf(stderr, "error: Unable to find required keywords for area table dimensions\n"); + return FH_NOT_FOUND; + } + + /* Sanity check. */ + if(num_cols != table->num_cols) + { + fprintf(stderr, "error: %d-column area table found, expected %d cols.\n", + num_cols, table->num_cols); + return FH_BAD_VALUE; + } + + if(num_rows > MAXBURN) + { + fprintf(stderr, + "error: too many boxes in area table. Max is %d, got %d\n", + MAXBURN, num_rows); + return FH_BAD_VALUE; + } + + table->num_rows = num_rows; + table->table_size = num_rows * row_chars; + + + if(!(data = malloc(table->table_size))) + { + fprintf(stderr, + "error: Unable to allocate %d bytes for area table.\n", + table->table_size); + return FH_NO_MEMORY; + } + + if(fh_read_image(hu, fh_file_desc(hu), data, + table->table_size, FH_TYPESIZE_8) != FH_SUCCESS) + { + fprintf(stderr, "error: Unable to map area table body for reading.\n"); + free(data); + return FH_INVALID; + } + + for(i = 0; i < num_rows; i++) + { + if((result = read_area(table, data, i)) != FH_SUCCESS) break; + } + + free(data); + return result; +} + +/* + * Reads one row of the fit table from a FITS file. + * + * Given information about a fit table and a row, this function + * reads a single row of data from a buffer containing the table + * body. + * + * Inputs: + * hu - header of the table extension for the fit table. + * table - describes structure of the table. + * data - points to buffer of data from fit table. + * row - the row to be read. + * + * Outputs: + * cell - cell number stored here. + * cx - center x from table row stored here. + * cy - center y from table row stored here. + * xfit - x of each value of correction start stored here. + * yfit - y of each value of correction start stored here. + * zero - zero of this fit stored here. + * + * This function returns FH_SUCCESS if the table row was read + * successfully. + * + * Note: this function must not be called until after read_area_table() + * has read out the burn areas across all cells and prepared + * boxbuf[]. This function assumes that each area's fit + * information has been filled out in boxbuf[], including + * allocating space for the fits. + */ +static fh_result +read_fit(HeaderUnit hu, fhTable * table, void * data, int row, + int * cell, int * cx, int * cy, int * xfit, int * yfit, double * zero) +{ + if(fh_table_read_value(table, data, row, FIT_TABLE_COL_CELL, cell) != FH_SUCCESS) + { + fprintf(stderr, + "error: Unable to get cell number from row %d of fit table\n", + row); + return FH_BAD_VALUE; + } + + if((*cell < 0) || (*cell > MAXCELL)) + { + fprintf(stderr, "error: illegal cell %d in fit table row %d\n", + *cell, row); + + return FH_BAD_VALUE; + } + + if((fh_table_read_value(table, data, row, FIT_TABLE_COL_CX, cx) != FH_SUCCESS) || + (fh_table_read_value(table, data, row, FIT_TABLE_COL_CY, cy) != FH_SUCCESS) || + (fh_table_read_value(table, data, row, FIT_TABLE_COL_XFIT, xfit) != FH_SUCCESS) || + (fh_table_read_value(table, data, row, FIT_TABLE_COL_YFIT, yfit) != FH_SUCCESS) || + (fh_table_read_value(table, data, row, FIT_TABLE_COL_ZERO, zero) != FH_SUCCESS)) + { + fprintf(stderr, + "error: Error reading values from row %d of burn area table\n", + row); + return FH_BAD_VALUE; + } + + return FH_SUCCESS; +} + +/* + * Reads the fit table from a FITS file. + * + * This function copies the data from a given fit table to boxbuf[]. + * Once complete, the contents of boxbuf[] should be equivalent to + * if persist_read() had read from a text file. + * + * Inputs: + * hu - header of the table extension for the fit table. + * table - describes structure of the table. + * num_areas - the number of areas that were previously read from + * the area table and written to boxbuf. + * + * Outputs: + * None. + * + * This function returns FH_SUCCESS if the table was read + * successfully. + * + * Note: this function must not be called until after read_area_table() + * has read out the burn areas across all cells and prepared + * boxbuf[]. This function assumes that each area's fit + * information has been filled out in boxbuf[], including + * allocating space for the fits. + */ +static fh_result +read_fit_table(HeaderUnit hu, fhTable * table, int num_areas) +{ + fh_result result; + int area, i; + int row_chars, num_cols, num_rows; + void * data; + int fit_table_row; + + int cell, cx, cy, xfit, yfit; + double zero; + + if((fh_get_int(hu, "NAXIS1", &row_chars) != FH_SUCCESS) || + (fh_get_int(hu, "NAXIS2", &num_rows) != FH_SUCCESS) || + (fh_get_int(hu, "TFIELDS", &num_cols) != FH_SUCCESS)) + { + fprintf(stderr, "error: Unable to find required keywords for fit table dimensions\n"); + return FH_NOT_FOUND; + } + + /* Sanity check. */ + if(num_cols != table->num_cols) + { + fprintf(stderr, "error: %d-column fit table found, expected %d cols.\n", + num_cols, table->num_cols); + return FH_BAD_VALUE; + } + + table->num_rows = num_rows; + table->table_size = num_rows * row_chars; + + /* Make space for the body of the table and copy it. + * + * %%% It would be preferable to just memory map the body of the + * table in the file itself, but I don't seem to be able to get + * that to work for reads. */ + if(!(data = malloc(table->table_size))) + { + fprintf(stderr, + "error: Unable to allocate %d bytes for area table.\n", + table->table_size); + return FH_NO_MEMORY; + } + if(fh_read_image(hu, fh_file_desc(hu), data, table->table_size, FH_TYPESIZE_8) != FH_SUCCESS) + { + fprintf(stderr, "error: Unable to read fit table body.\n"); + free(data); + return FH_INVALID; + } + + /* For all areas, grab all of the fits for that area. This should + * end up consuming the whole table. */ + fit_table_row = 0; + for(area = 0; area < num_areas; area++) + { + for(i = 0; i < boxbuf[area].nfit; i++) + { + if((result = read_fit(hu, table, data, fit_table_row, + &cell, &cx, &cy, &xfit, &yfit, &zero)) != FH_SUCCESS) + { + free(data); + return result; + } + + /* This shouldn't happen, but just in case... */ + if((cell != boxbuf[area].cell) || + (cx != boxbuf[area].cx) || + (cy != boxbuf[area].cy)) + { + fprintf(stderr, + "error: Fit in table row %d does not match current area " + "(table cell=%d cx=%d cy=%d vs. area cell=%d cx=%d cy=%d)\n", + fit_table_row, + cell, cx, cy, boxbuf[area].cell, boxbuf[area].cx, boxbuf[area].cy); + free(data); + return FH_BAD_VALUE; + } + + boxbuf[area].xfit[i] = xfit; + boxbuf[area].yfit[i] = yfit; + boxbuf[area].zero[i] = zero; + + /* Move on to the next row of the table. */ + fit_table_row++; + } + } + + free(data); + return FH_SUCCESS; +} + +/* + * Retrieves burn information from a FITS file. + * + * This function reads the burn tables from an existing FITS file + * and builds up burntool's internal table of areas and fits from + * the information within. This function should be equivalent to + * persist_read(), except that it operates on a FITS file instead + * of text. + * + * Inputs: + * cell - array of 64 CELL structures that describe all of the + * burn info for the whole OTA. Burn info will be written + * here. + * filename - filename of FITS file to read from. + * + * Outputs: + * None. + * + * This function returns FH_SUCCESS if the tables were successfully + * appended to the FITS file. + * + * Note: this function assumes that the passed file is a MEF and that + * any required sanity checks on the file have already been performed. + */ +fh_result +persist_fits_read(CELL *cell, const char * filename) +{ + HeaderUnit phu; + HeaderUnit area_hu, fit_hu; + fh_result result; + + if(!(phu = fh_create())) + { + fprintf(stderr, + "error: unable to create header structure for FITS file \"%s\"\n", + filename); + return FH_NO_MEMORY; + } + + /* Try to open the FITS file. */ + if ((result = fh_file(phu, filename, FH_FILE_RDONLY)) != FH_SUCCESS) + { + fprintf(stderr, + "error: unable to open FITS file \"%s\"\n", + filename); + fh_destroy(phu); + return result; + } + + + /* Don't even try if this FITS file doesn't have any burn info in it. */ + find_existing_tables(phu, &area_hu, NULL, &fit_hu, NULL); + if(!area_hu || !fit_hu) + { + fprintf(stderr, + "error: Unable to find persistence info in FITS file.\n"); + fh_destroy(phu); + return FH_NOT_FOUND; + } + + fh_table_init(&area_table); + fh_table_init(&fit_table); + + fh_ehu_by_extname(phu, DEFAULT_EXTNAME_AREA_TABLE); + if((result = read_area_table(area_hu, &area_table)) != FH_SUCCESS) + { + free(area_table.strbuf); + free(fit_table.strbuf); + fh_destroy(phu); + return result; + } + + fh_ehu_by_extname(phu, DEFAULT_EXTNAME_FIT_TABLE); + if((result = read_fit_table(fit_hu, &fit_table, area_table.num_rows)) != FH_SUCCESS) + { + free(area_table.strbuf); + free(fit_table.strbuf); + fh_destroy(phu); + return result; + } + + /* Now copy over the boxes that were read to our array of cells. + This section is ripped off from persist_read(). */ + { + int i, k; + int nbox; + +/* Initialize the counts */ + for(k=0; k 0) { + cell[k].persist = (OBJBOX *)calloc(i, sizeof(OBJBOX)); + cell[k].npersist = 0; + } + } +/* Copy the results to the cells */ + for(i=0; inpersist; 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 */ Index: trunk/extsrc/gpcsw/gpcsrc/fits/burntool/persistio.c =================================================================== --- trunk/extsrc/gpcsw/gpcsrc/fits/burntool/persistio.c (revision 23564) +++ trunk/extsrc/gpcsw/gpcsrc/fits/burntool/persistio.c (revision 23924) @@ -233,5 +233,5 @@ 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 #include @@ -12,11 +12,4 @@ #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,101 +17,4 @@ 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= 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 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); } Index: trunk/extsrc/gpcsw/gpcsrc/fits/burntool/trailfit.c =================================================================== --- trunk/extsrc/gpcsw/gpcsrc/fits/burntool/trailfit.c (revision 23564) +++ trunk/extsrc/gpcsw/gpcsrc/fits/burntool/trailfit.c (revision 23924) @@ -148,8 +148,16 @@ } +/* 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); }