Index: trunk/Ohana/src/relphot/Makefile =================================================================== --- trunk/Ohana/src/relphot/Makefile (revision 35104) +++ trunk/Ohana/src/relphot/Makefile (revision 35416) @@ -40,8 +40,10 @@ $(SRC)/load_catalogs.$(ARCH).o \ $(SRC)/load_images.$(ARCH).o \ +$(SRC)/synthetic_mags.$(ARCH).o \ $(SRC)/plot_scatter.$(ARCH).o \ $(SRC)/plotstuff.$(ARCH).o \ $(SRC)/reload_catalogs.$(ARCH).o \ $(SRC)/relphot.$(ARCH).o \ +$(SRC)/relphot_images.$(ARCH).o \ $(SRC)/relphot_objects.$(ARCH).o \ $(SRC)/select_images.$(ARCH).o \ @@ -69,4 +71,5 @@ $(SRC)/args.$(ARCH).o \ $(SRC)/help.$(ARCH).o \ +$(SRC)/synthetic_mags.$(ARCH).o \ $(SRC)/plotstuff.$(ARCH).o \ $(SRC)/liststats.$(ARCH).o \ Index: trunk/Ohana/src/relphot/include/relphot.h =================================================================== --- trunk/Ohana/src/relphot/include/relphot.h (revision 35104) +++ trunk/Ohana/src/relphot/include/relphot.h (revision 35416) @@ -4,4 +4,13 @@ # include # include + +# define MARKTIME(MSG,...) { \ + gettimeofday (&stopTimer, (void *) NULL); \ + float dtime = DTIME (stopTimer, startTimer); \ + fprintf (stderr, MSG, __VA_ARGS__); } + +# define INITTIME \ + struct timeval startTimer, stopTimer; \ + gettimeofday (&startTimer, (void *) NULL); /* # define GRID_V1 */ @@ -193,4 +202,5 @@ int UpdateAverages; int ApplyOffsets; +int SyntheticPhotometry; char *PhotcodeList; @@ -353,4 +363,5 @@ void write_coords PROTO((Header *header, Coords *coords)); int relphot_objects (int hostID, char *hostpath); +int relphot_images (void); void relphot_usage (void); @@ -390,2 +401,5 @@ void SetMrelInfoInit (SetMrelInfo *results, int allocLists); + +int init_synthetic_mags (); +int add_synthetic_mags (AverageTiny *average, SecFilt *secfilt, MeasureTiny *measure, off_t *Nmeasure, off_t *Nm); Index: trunk/Ohana/src/relphot/src/MosaicOps.c =================================================================== --- trunk/Ohana/src/relphot/src/MosaicOps.c (revision 35104) +++ trunk/Ohana/src/relphot/src/MosaicOps.c (revision 35416) @@ -152,10 +152,4 @@ } -# define MARKTIME(MSG,...) { \ - float dtime; \ - gettimeofday (&stopTimer, (void *) NULL); \ - dtime = DTIME (stopTimer, startTimer); \ - fprintf (stderr, MSG, __VA_ARGS__); } - /* find mosaic frames (unique time periods) (NOTE : we do NOT require matching photcodes...) this function will also identify the images NOT in the subset which belong to a selected mosaic @@ -167,9 +161,7 @@ char *pname; - struct timeval startTimer, stopTimer; - if (!MOSAIC_ZEROPT) return; - gettimeofday (&startTimer, (void *) NULL); + INITTIME; /* a 'mosaic' in relphot is (unlike relastro) a virtual concept: there is no @@ -264,5 +256,5 @@ } } - fprintf (stderr, "%d total images, %d overlap skyregion, %d do not overlap, but match overlapping mosaics\n", (int) Nimage, (int) Nskip, (int) Nmark); + fprintf (stderr, "%d total images, %d overlap skyregion & match selection criteria, %d do not overlap, but match overlapping mosaics\n", (int) Nimage, (int) Nskip, (int) Nmark); MARKTIME("find unselected images matching selected mosaics: %f sec\n", dtime); Index: trunk/Ohana/src/relphot/src/args.c =================================================================== --- trunk/Ohana/src/relphot/src/args.c (revision 35104) +++ trunk/Ohana/src/relphot/src/args.c (revision 35416) @@ -251,4 +251,8 @@ remove_argument (N, &argc, argv); MOSAIC_ZEROPT = TRUE; + if (!strcasecmp (MOSAICNAME, "none")) { + fprintf (stderr, "mosaic astrometry selected by MOSAICNAME not defined\n"); + exit (2); + } } @@ -315,4 +319,11 @@ remove_argument (N, &argc, argv); DophotSelect = TRUE; + } + + SyntheticPhotometry = FALSE; + if ((N = get_argument (argc, argv, "-synthphot"))) { + remove_argument (N, &argc, argv); + SyntheticPhotometry = TRUE; + init_synthetic_mags(); } @@ -516,4 +527,11 @@ } + SyntheticPhotometry = FALSE; + if ((N = get_argument (argc, argv, "-synthphot"))) { + remove_argument (N, &argc, argv); + SyntheticPhotometry = FALSE; + init_synthetic_mags(); + } + ImagSelect = FALSE; if ((N = get_argument (argc, argv, "-instmag"))) { Index: trunk/Ohana/src/relphot/src/bcatalog.c =================================================================== --- trunk/Ohana/src/relphot/src/bcatalog.c (revision 35104) +++ trunk/Ohana/src/relphot/src/bcatalog.c (revision 35416) @@ -179,4 +179,5 @@ REALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE); } + myAssert (Nmeasure < NMEASURE, "realloc failure"); } @@ -204,4 +205,17 @@ continue; } + + // for w-band photometry (& other cases?) convert gri photometry + // to a synthetic w-band magnitude + add_synthetic_mags (&subcatalog[0].averageT[Naverage], + &subcatalog[0].secfilt[Nsecfilt*Naverage], + &subcatalog[0].measureT[Nmeasure], + &Nmeasure, + &Nm); + if (Nmeasure == NMEASURE) { + NMEASURE += 1000; + REALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE); + } + myAssert (Nmeasure < NMEASURE, "realloc failure"); subcatalog[0].averageT[Naverage].Nmeasure = Nm; Index: trunk/Ohana/src/relphot/src/global_stats.c =================================================================== --- trunk/Ohana/src/relphot/src/global_stats.c (revision 35104) +++ trunk/Ohana/src/relphot/src/global_stats.c (revision 35416) @@ -1,9 +1,3 @@ # include "relphot.h" - -# define MARKTIME(MSG,...) { \ - float dtime; \ - gettimeofday (&stopTimer, (void *) NULL); \ - dtime = DTIME (stopTimer, startTimer); \ - fprintf (stderr, MSG, __VA_ARGS__); } void global_stats (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) { @@ -11,10 +5,8 @@ StatType stN, stX, stS, imN, imX, imM, imD, msM, msX, msN, msD; - // struct timeval startTimer, stopTimer; + // INITTIME; fprintf (stderr, "\n"); fprintf (stderr, "STATS median mean sigma min max Nmeas\n"); - - // gettimeofday (&startTimer, (void *) NULL); int Ns; Index: trunk/Ohana/src/relphot/src/load_catalogs.c =================================================================== --- trunk/Ohana/src/relphot/src/load_catalogs.c (revision 35104) +++ trunk/Ohana/src/relphot/src/load_catalogs.c (revision 35416) @@ -148,12 +148,14 @@ char tmpline[1024]; - if (VERBOSE) { snprintf (tmpline, 1024, "%s -v", command); strcpy (command, tmpline); } - if (VERBOSE2) { snprintf (tmpline, 1024, "%s -vv", command); strcpy (command, tmpline); } - if (RESET) { snprintf (tmpline, 1024, "%s -reset", command); strcpy (command, tmpline); } - if (RESET_ZEROPTS) { snprintf (tmpline, 1024, "%s -reset-zpts", command); strcpy (command, tmpline); } - if (!KEEP_UBERCAL) { snprintf (tmpline, 1024, "%s -reset-ubercal", command); strcpy (command, tmpline); } - if (DophotSelect) { snprintf (tmpline, 1024, "%s -dophot %d", command, DophotValue); strcpy (command, tmpline); } - if (ImagSelect) { snprintf (tmpline, 1024, "%s -instmag %f %f", command, ImagMin, ImagMax); strcpy (command, tmpline); } - if (MaxDensityUse) { snprintf (tmpline, 1024, "%s -max-density %f", command, MaxDensityValue); strcpy (command, tmpline); } + if (VERBOSE) { snprintf (tmpline, 1024, "%s -v", command); strcpy (command, tmpline); } + if (VERBOSE2) { snprintf (tmpline, 1024, "%s -vv", command); strcpy (command, tmpline); } + if (RESET) { snprintf (tmpline, 1024, "%s -reset", command); strcpy (command, tmpline); } + if (RESET_ZEROPTS) { snprintf (tmpline, 1024, "%s -reset-zpts", command); strcpy (command, tmpline); } + if (!KEEP_UBERCAL) { snprintf (tmpline, 1024, "%s -reset-ubercal", command); strcpy (command, tmpline); } + if (DophotSelect) { snprintf (tmpline, 1024, "%s -dophot %d", command, DophotValue); strcpy (command, tmpline); } + if (ImagSelect) { snprintf (tmpline, 1024, "%s -instmag %f %f", command, ImagMin, ImagMax); strcpy (command, tmpline); } + if (MaxDensityUse) { snprintf (tmpline, 1024, "%s -max-density %f", command, MaxDensityValue); strcpy (command, tmpline); } + if (SyntheticPhotometry) { snprintf (tmpline, 1024, "%s -synthphot", command); strcpy (command, tmpline); } + if (TimeSelect) { char *tstart = ohana_sec_to_date (TSTART); Index: trunk/Ohana/src/relphot/src/load_images.c =================================================================== --- trunk/Ohana/src/relphot/src/load_images.c (revision 35104) +++ trunk/Ohana/src/relphot/src/load_images.c (revision 35416) @@ -1,9 +1,3 @@ # include "relphot.h" - -# define MARKTIME(MSG,...) { \ - float dtime; \ - gettimeofday (&stop, (void *) NULL); \ - dtime = DTIME (stop, start); \ - fprintf (stderr, MSG, __VA_ARGS__); } // This function generates a subset of the images based on selections. Input db has already @@ -14,5 +8,4 @@ off_t Nimage, Nsubset, Nchar; off_t *LineNumber; - struct timeval start, stop; char *inSubset; @@ -20,5 +13,5 @@ SkyList *skylist = NULL; - gettimeofday (&start, (void *) NULL); + INITTIME; // load the current sky table (layout of all SkyRegions) @@ -93,5 +86,5 @@ } -# if (1) +# if (0) FILE *ftest = fopen ("skydump.dat", "w"); int i; Index: trunk/Ohana/src/relphot/src/relphot.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot.c (revision 35104) +++ trunk/Ohana/src/relphot/src/relphot.c (revision 35416) @@ -1,19 +1,5 @@ # include "relphot.h" -# define MARKTIME(MSG,...) { \ - float dtime; \ - gettimeofday (&stop, (void *) NULL); \ - dtime = DTIME (stop, start); \ - fprintf (stderr, MSG, __VA_ARGS__); } - int main (int argc, char **argv) { - - int i, status, Ncatalog; - Catalog *catalog = NULL; - FITS_DB db; - struct timeval start, stop; - SkyList *skylist = NULL; - - gettimeofday (&start, (void *) NULL); /* get configuration info, args */ @@ -36,242 +22,5 @@ } - /* register database handle with shutdown procedure */ - set_db (&db); - db.mode = dvo_catalog_catmode (CATMODE); - db.format = dvo_catalog_catformat (CATFORMAT); - - /* lock and load the image db table */ - status = dvo_image_lock (&db, ImageCat, 60.0, (UPDATE ? LCK_XCLD : LCK_SOFT)); - if (!status && UPDATE) { - fprintf (stderr, "error\n"); - Shutdown ("ERROR: failure to lock image catalog %s", db.filename); - } - // if the file is missing, db.dbstate will have a value of either: - // LCK_EMPTY (if UPDATE) or LCK_MISSING (if !UPDATE) - if ((db.dbstate == LCK_EMPTY) || (db.dbstate == LCK_MISSING)) { - - // 25.0 is the nominal zero point for measurements in the database - // for all measurements except PHOT_REF: - // measure.M contains: -2.5*log(DN) + 2.5*log(exptime) + 25.0 - // for measurements using PHOT_REF: - // measure.M contains: -2.5*log(DN) + 2.5*log(exptime) + True Zero Point (but exptime and ZP are known) - dvo_image_create (&db, 25.0); - - // XXX why are we running relphot on an empty database? - } else { - if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename); - } - MARKTIME("-- load image data: %f sec\n", dtime); - - /* load regions and images based on specified sky patch and/or catalog */ - skylist = load_images (&db, UserCatalog, &UserPatch); - MARKTIME("-- load images: %f sec\n", dtime); - - /* unlock, if we can (else, unlocked below) */ - if (!UPDATE) dvo_image_unlock (&db); - - // load the flat correction table (if defined) - char flatcorrfile[256]; - sprintf (flatcorrfile, "%s/flatcorr.fits", CATDIR); - FlatCorrectionTable *flatcorr = FlatCorrectionLoad (flatcorrfile, VERBOSE); - - if (NLOOP > 0) { - /* load catalog data from region files (hostID is 0 since we are not a client */ - catalog = load_catalogs (skylist, &Ncatalog, 0, NULL); - MARKTIME("-- load catalog data: %f sec\n", dtime); - - /* add in a loop over the catalogs calling dvo_catalog_chipcoords */ - - /* match measurements with images, mosaics */ - initImageBins (catalog, Ncatalog, TRUE); - MARKTIME("-- make image bins: %f sec\n", dtime); - - initMosaicBins (catalog, Ncatalog, TRUE); - initGridBins (catalog, Ncatalog); - initMrel (catalog, Ncatalog); - - findImages (catalog, Ncatalog, TRUE); - MARKTIME("-- set up image indexes: %f sec\n", dtime); - - findMosaics (catalog, Ncatalog, TRUE); /* also sets Grid values */ - MARKTIME("-- set up mosaic indexes: %f sec\n", dtime); - - SAVEPLOT = FALSE; - - setExclusions (catalog, Ncatalog, TRUE); - - global_stats (catalog, Ncatalog, flatcorr); - - if (PLOTSTUFF) { - plot_star_coords (catalog, Ncatalog); - // plot_mosaic_fields (catalog); - } - - // if we are measuring the flat-field correction grid, we need to perform a number of iterations first: - if (USE_GRID) { - int star_toofew; - - // until we finish the grid analysis, do not reject stars out-of-hand based on ID_STAR_FEW - // XXX this is kind of poor: need to have a better distinctions about STAR_BAD in setMrel vs getMrel - star_toofew = STAR_TOOFEW; - STAR_TOOFEW = 0; - for (i = 0; i < NGRID; i++) { - STAR_BAD = ID_STAR_POOR; - setMrel (catalog, Ncatalog, flatcorr); - STAR_BAD = ID_STAR_POOR | ID_STAR_FEW; - setMgrid (catalog, flatcorr); - } - STAR_BAD = ID_STAR_POOR | ID_STAR_FEW; - STAR_TOOFEW = star_toofew; - } - - /* determine fit values */ - for (i = 0; i < NLOOP; i++) { - if (PLOTSTUFF) { - plot_scatter (catalog, Ncatalog, flatcorr); - } - setMrel (catalog, Ncatalog, flatcorr); // threaded - if (PLOTSTUFF) { - plot_scatter (catalog, Ncatalog, flatcorr); - } - setMcal (catalog, FALSE, flatcorr); - setMmos (catalog, FALSE, flatcorr); - setMgrid (catalog, flatcorr); - MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime); - - if (PLOTSTUFF) { - plot_scatter (catalog, Ncatalog, flatcorr); - plot_grid (catalog, flatcorr); - plot_mosaics (); - plot_images (); - plot_stars (catalog, Ncatalog); - plot_chisq (catalog, Ncatalog); - } - // if (i < NLOOP - 1) rationalize_mosaics (catalog, Ncatalog); - // if (i % 6 == 1) rationalize_images (); - - // NOTE : in the past, I was not iterating enough before cleaning. make sure we do - // at least 8 loops first -- that should get the systematic errors down to the ~1% - // level, even in cases where we have an even split between photometric data and - // data with 1 mag of extinction. - - if ((i > 8) && (i % 8 == 2)) clean_measures (catalog, Ncatalog, FALSE, flatcorr); - if ((i > 8) && (i % 8 == 3)) clean_stars (catalog, Ncatalog); - if ((i > 8) && (i % 8 == 5)) clean_mosaics (); - if ((i > 8) && (i % 8 == 5)) clean_images (); - - // if ((i == 1) || (i == 5) || (i == 9) || (i == 13)) clean_measures (catalog, Ncatalog, FALSE); - // if ((i == 2) || (i == 6) || (i == 10) || (i == 14)) clean_stars (catalog, Ncatalog); - // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_mosaics (); - // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_images (); - if (i % 3 == 2) global_stats (catalog, Ncatalog, flatcorr); - MARKTIME("-- finished loop %d: %f sec\n", i, dtime); - } - - if (PLOTSTUFF) { - plot_scatter (catalog, Ncatalog, flatcorr); - plot_grid (catalog, flatcorr); - plot_mosaics (); - plot_images (); - plot_stars (catalog, Ncatalog); - plot_chisq (catalog, Ncatalog); - } - - if (USE_GRID) dump_grid (); - - /* set Mcal & Mmos for bad images */ - setMcal (catalog, TRUE, flatcorr); - setMmos (catalog, TRUE, flatcorr); - MARKTIME("-- finalize Mcal values: %f sec\n", dtime); - - setMcalFinal (); // copy per-mosaic calibrations to the images - - if (SAVE_IMAGE_UPDATES) { - - FITS_DB dbX; - off_t *LineNumber; - Image *image, *subset; - off_t Nimage, Nsubset, i, Nx; - char filename[1024]; - - gfits_db_init (&dbX); - dbX.lockstate = LCK_XCLD; - dbX.timeout = 60.0; - dbX.mode = db.mode; - dbX.format = db.format; - - snprintf (filename, 1024, "%s.bck", ImageCat); - if (!gfits_db_lock (&dbX, filename)) { - fprintf (stderr, "can't lock backup image image catalog\n"); - return (FALSE); - } - - // apply the changes from the image subset to the full table: - - // convert database table to internal structure (binary to Image) - // 'image' points to the same memory as db->ftable->buffer - image = gfits_table_get_Image (&db.ftable, &Nimage, &db.swapped); - if (!image) { - fprintf (stderr, "ERROR: failed to read images\n"); - exit (2); - } - gfits_scan (db.ftable.header, "NAXIS1", OFF_T_FMT, 1, &Nx); - - subset = getimages (&Nsubset, &LineNumber); - for (i = 0; i < Nsubset; i++) { - if (LineNumber[i] == -1) continue; - memcpy (&image[LineNumber[i]], &subset[i], Nx); - } - - // copy Images.dat data (all or only the subset / vtable elements?) I think I need to dump - // the entire Image table, but with the updates in place I think this says: re-work the - // ftable/vtable usage in this program to be more sensible... - gfits_copy_header (&db.header, &dbX.header); - gfits_copy_matrix (&db.matrix, &dbX.matrix); - gfits_copy_header (&db.theader, &dbX.theader); - gfits_copy_ftable (&db.ftable, &dbX.ftable); - - dbX.ftable.header = &dbX.theader; - dbX.virtual = FALSE; - - // save Images.dat using the copied structure - if (UPDATE_CATFORMAT) { - // ensure the db format is updated - dbX.format = dvo_catalog_catformat (UPDATE_CATFORMAT); - char photcodeFile[1024]; - sprintf (photcodeFile, "%s/Photcodes.dat", CATDIR); - SavePhotcodesFITS (photcodeFile); - } - dvo_image_save (&dbX, VERBOSE); - dvo_image_unlock (&dbX); - MARKTIME("-- save Image.dat.bck: %f sec\n", dtime); - } - - /* at this point, we have correct cal coeffs in the image/mosaic structures */ - for (i = 0; i < Ncatalog; i++) { - free_tiny_values (&catalog[i]); - dvo_catalog_free (&catalog[i]); - } - freeImageBins (Ncatalog, TRUE); - freeMosaicBins (Ncatalog, TRUE); - freeGridBins (Ncatalog); - - } // (NLOOP > 0) : this loop determines the offsets per chip - - reload_images (&db); - - // only change the real database files if -update is requested - if (!UPDATE) exit (0); - - /* load catalog data from region files, update Mrel include all data */ - reload_catalogs (skylist, flatcorr, 0, NULL); - MARKTIME("-- updated all catalogs: %f sec\n", dtime); - - if (UPDATE_CATFORMAT) { - // ensure the db format is updated - db.format = dvo_catalog_catformat (UPDATE_CATFORMAT); - } - dvo_image_update (&db, VERBOSE); - dvo_image_unlock (&db); + relphot_images (); exit (0); Index: trunk/Ohana/src/relphot/src/relphot_images.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_images.c (revision 35416) +++ trunk/Ohana/src/relphot/src/relphot_images.c (revision 35416) @@ -0,0 +1,253 @@ +# include "relphot.h" + +// 25.0 is the nominal zero point for measurements in the database +// for all measurements except PHOT_REF: +// measure.M contains: -2.5*log(DN) + 2.5*log(exptime) + 25.0 +// for measurements using PHOT_REF: +// measure.M contains: -2.5*log(DN) + 2.5*log(exptime) + True Zero Point (but exptime and ZP are known) + +int relphot_images () { + + int i, status, Ncatalog; + Catalog *catalog = NULL; + FITS_DB db; + SkyList *skylist = NULL; + + INITTIME; + + /* register database handle with shutdown procedure */ + set_db (&db); + db.mode = dvo_catalog_catmode (CATMODE); + db.format = dvo_catalog_catformat (CATFORMAT); + + /* lock and load the image db table */ + status = dvo_image_lock (&db, ImageCat, 60.0, (UPDATE ? LCK_XCLD : LCK_SOFT)); + if (!status && UPDATE) { + fprintf (stderr, "error\n"); + Shutdown ("ERROR: failure to lock image catalog %s", db.filename); + } + + // if the file is missing, db.dbstate will have a value of either: + // LCK_EMPTY (if UPDATE) or LCK_MISSING (if !UPDATE) + if ((db.dbstate == LCK_EMPTY) || (db.dbstate == LCK_MISSING)) { + Shutdown ("ERROR: database %s contains no image data", CATDIR); + } + + // read data from Image.dat file + if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename); + MARKTIME("-- load image data: %f sec\n", dtime); + + /* load regions and images based on specified sky patch and/or catalog */ + skylist = load_images (&db, UserCatalog, &UserPatch); + MARKTIME("-- load images: %f sec\n", dtime); + + /* unlock, if we can (else, unlocked below) */ + if (!UPDATE) dvo_image_unlock (&db); + + // load the flat correction table (if defined) + char flatcorrfile[256]; + sprintf (flatcorrfile, "%s/flatcorr.fits", CATDIR); + FlatCorrectionTable *flatcorr = FlatCorrectionLoad (flatcorrfile, VERBOSE); + + if (NLOOP > 0) { + /* load catalog data from region files (hostID is 0 since we are not a client */ + catalog = load_catalogs (skylist, &Ncatalog, 0, NULL); + MARKTIME("-- load catalog data: %f sec\n", dtime); + + /* add in a loop over the catalogs calling dvo_catalog_chipcoords */ + + /* match measurements with images, mosaics */ + initImageBins (catalog, Ncatalog, TRUE); + MARKTIME("-- make image bins: %f sec\n", dtime); + + initMosaicBins (catalog, Ncatalog, TRUE); + initGridBins (catalog, Ncatalog); + initMrel (catalog, Ncatalog); + + findImages (catalog, Ncatalog, TRUE); + MARKTIME("-- set up image indexes: %f sec\n", dtime); + + findMosaics (catalog, Ncatalog, TRUE); /* also sets Grid values */ + MARKTIME("-- set up mosaic indexes: %f sec\n", dtime); + + SAVEPLOT = FALSE; + + setExclusions (catalog, Ncatalog, TRUE); + + global_stats (catalog, Ncatalog, flatcorr); + + if (PLOTSTUFF) { + plot_star_coords (catalog, Ncatalog); + // plot_mosaic_fields (catalog); + } + + // if we are measuring the flat-field correction grid, we need to perform a number of iterations first: + if (USE_GRID) { + int star_toofew; + + // until we finish the grid analysis, do not reject stars out-of-hand based on ID_STAR_FEW + // XXX this is kind of poor: need to have a better distinctions about STAR_BAD in setMrel vs getMrel + star_toofew = STAR_TOOFEW; + STAR_TOOFEW = 0; + for (i = 0; i < NGRID; i++) { + STAR_BAD = ID_STAR_POOR; + setMrel (catalog, Ncatalog, flatcorr); + STAR_BAD = ID_STAR_POOR | ID_STAR_FEW; + setMgrid (catalog, flatcorr); + } + STAR_BAD = ID_STAR_POOR | ID_STAR_FEW; + STAR_TOOFEW = star_toofew; + } + + /* determine fit values */ + for (i = 0; i < NLOOP; i++) { + if (PLOTSTUFF) { + plot_scatter (catalog, Ncatalog, flatcorr); + } + setMrel (catalog, Ncatalog, flatcorr); // threaded + if (PLOTSTUFF) { + plot_scatter (catalog, Ncatalog, flatcorr); + } + setMcal (catalog, FALSE, flatcorr); + setMmos (catalog, FALSE, flatcorr); + setMgrid (catalog, flatcorr); + MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime); + + if (PLOTSTUFF) { + plot_scatter (catalog, Ncatalog, flatcorr); + plot_grid (catalog, flatcorr); + plot_mosaics (); + plot_images (); + plot_stars (catalog, Ncatalog); + plot_chisq (catalog, Ncatalog); + } + // if (i < NLOOP - 1) rationalize_mosaics (catalog, Ncatalog); + // if (i % 6 == 1) rationalize_images (); + + // NOTE : in the past, I was not iterating enough before cleaning. make sure we do + // at least 8 loops first -- that should get the systematic errors down to the ~1% + // level, even in cases where we have an even split between photometric data and + // data with 1 mag of extinction. + + if ((i > 8) && (i % 8 == 2)) clean_measures (catalog, Ncatalog, FALSE, flatcorr); + if ((i > 8) && (i % 8 == 3)) clean_stars (catalog, Ncatalog); + if ((i > 8) && (i % 8 == 5)) clean_mosaics (); + if ((i > 8) && (i % 8 == 5)) clean_images (); + + // if ((i == 1) || (i == 5) || (i == 9) || (i == 13)) clean_measures (catalog, Ncatalog, FALSE); + // if ((i == 2) || (i == 6) || (i == 10) || (i == 14)) clean_stars (catalog, Ncatalog); + // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_mosaics (); + // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_images (); + if (i % 3 == 2) global_stats (catalog, Ncatalog, flatcorr); + MARKTIME("-- finished loop %d: %f sec\n", i, dtime); + } + + if (PLOTSTUFF) { + plot_scatter (catalog, Ncatalog, flatcorr); + plot_grid (catalog, flatcorr); + plot_mosaics (); + plot_images (); + plot_stars (catalog, Ncatalog); + plot_chisq (catalog, Ncatalog); + } + + if (USE_GRID) dump_grid (); + + /* set Mcal & Mmos for bad images */ + setMcal (catalog, TRUE, flatcorr); + setMmos (catalog, TRUE, flatcorr); + MARKTIME("-- finalize Mcal values: %f sec\n", dtime); + + setMcalFinal (); // copy per-mosaic calibrations to the images + + if (SAVE_IMAGE_UPDATES) { + + FITS_DB dbX; + off_t *LineNumber; + Image *image, *subset; + off_t Nimage, Nsubset, i, Nx; + char filename[1024]; + + gfits_db_init (&dbX); + dbX.lockstate = LCK_XCLD; + dbX.timeout = 60.0; + dbX.mode = db.mode; + dbX.format = db.format; + + snprintf (filename, 1024, "%s.bck", ImageCat); + if (!gfits_db_lock (&dbX, filename)) { + fprintf (stderr, "can't lock backup image image catalog\n"); + return (FALSE); + } + + // apply the changes from the image subset to the full table: + + // convert database table to internal structure (binary to Image) + // 'image' points to the same memory as db->ftable->buffer + image = gfits_table_get_Image (&db.ftable, &Nimage, &db.swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } + gfits_scan (db.ftable.header, "NAXIS1", OFF_T_FMT, 1, &Nx); + + subset = getimages (&Nsubset, &LineNumber); + for (i = 0; i < Nsubset; i++) { + if (LineNumber[i] == -1) continue; + memcpy (&image[LineNumber[i]], &subset[i], Nx); + } + + // copy Images.dat data (all or only the subset / vtable elements?) I think I need to dump + // the entire Image table, but with the updates in place I think this says: re-work the + // ftable/vtable usage in this program to be more sensible... + gfits_copy_header (&db.header, &dbX.header); + gfits_copy_matrix (&db.matrix, &dbX.matrix); + gfits_copy_header (&db.theader, &dbX.theader); + gfits_copy_ftable (&db.ftable, &dbX.ftable); + + dbX.ftable.header = &dbX.theader; + dbX.virtual = FALSE; + + // save Images.dat using the copied structure + if (UPDATE_CATFORMAT) { + // ensure the db format is updated + dbX.format = dvo_catalog_catformat (UPDATE_CATFORMAT); + char photcodeFile[1024]; + sprintf (photcodeFile, "%s/Photcodes.dat", CATDIR); + SavePhotcodesFITS (photcodeFile); + } + dvo_image_save (&dbX, VERBOSE); + dvo_image_unlock (&dbX); + MARKTIME("-- save Image.dat.bck: %f sec\n", dtime); + } + + /* at this point, we have correct cal coeffs in the image/mosaic structures */ + for (i = 0; i < Ncatalog; i++) { + free_tiny_values (&catalog[i]); + dvo_catalog_free (&catalog[i]); + } + freeImageBins (Ncatalog, TRUE); + freeMosaicBins (Ncatalog, TRUE); + freeGridBins (Ncatalog); + + } // (NLOOP > 0) : this loop determines the offsets per chip + + reload_images (&db); + + // only change the real database files if -update is requested + if (!UPDATE) exit (0); + + /* load catalog data from region files, update Mrel include all data */ + reload_catalogs (skylist, flatcorr, 0, NULL); + MARKTIME("-- updated all catalogs: %f sec\n", dtime); + + if (UPDATE_CATFORMAT) { + // ensure the db format is updated + db.format = dvo_catalog_catformat (UPDATE_CATFORMAT); + } + dvo_image_update (&db, VERBOSE); + dvo_image_unlock (&db); + + exit (0); +} + Index: trunk/Ohana/src/relphot/src/select_images.c =================================================================== --- trunk/Ohana/src/relphot/src/select_images.c (revision 35104) +++ trunk/Ohana/src/relphot/src/select_images.c (revision 35416) @@ -17,10 +17,4 @@ off_t getRegionStartByRA (double R, double *Rref, off_t Nregions); -# define MARKTIME(MSG,...) { \ - float dtime; \ - gettimeofday (&stop, (void *) NULL); \ - dtime = DTIME (stop, start); \ - fprintf (stderr, MSG, __VA_ARGS__); } - Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, char *inSubset, off_t **LineNumber, off_t *Nimage, SkyRegion *region) { @@ -31,5 +25,4 @@ double Ri[5], Di[5], Xi[5], Yi[5]; Coords tcoords; - struct timeval start, stop; if (skylist[0].Nregions < 1) { @@ -40,5 +33,5 @@ } - gettimeofday (&start, (void *) NULL); + INITTIME; // the comparison is made in the catalog local projection. below we set crval1,2 @@ -82,5 +75,5 @@ } -# if (1) +# if (0) // XXX quick test of the dec bands: FILE *fout = fopen ("dec.bands.dat", "w"); @@ -132,5 +125,5 @@ // this adds 1.3 sec for 3M images if (!FindMosaicForImage (timage, Ntimage, i)) { - fprintf (stderr, "cannot find mosaic for "OFF_T_FMT"\n", i); + if (VERBOSE2) fprintf (stderr, "cannot find mosaic for "OFF_T_FMT"\n", i); continue; } Index: trunk/Ohana/src/relphot/src/synthetic_mags.c =================================================================== --- trunk/Ohana/src/relphot/src/synthetic_mags.c (revision 35416) +++ trunk/Ohana/src/relphot/src/synthetic_mags.c (revision 35416) @@ -0,0 +1,71 @@ +# include "relphot.h" + +static int NSr = -1; +static int NSi = -1; +static int W_SYNTH = -1; + +int init_synthetic_mags () { + int Pr = GetPhotcodeCodebyName ("r"); + int Pi = GetPhotcodeCodebyName ("i"); + + NSr = GetPhotcodeNsec (Pr); + NSi = GetPhotcodeNsec (Pi); + + W_SYNTH = GetPhotcodeCodebyName ("SYNTH.w"); + + return TRUE; +} + +int add_synthetic_mags (AverageTiny *average, SecFilt *secfilt, MeasureTiny *measure, off_t *Nmeasure, off_t *Nm) { + + if (!SyntheticPhotometry) return TRUE; + + myAssert (NSr >= 0, "failed to call init_synthetic_mags"); + myAssert (NSi >= 0, "failed to call init_synthetic_mags"); + myAssert (W_SYNTH >= 0, "failed to call init_synthetic_mags"); + + // XXX for now, hard-wire the photometric relationship + + // we want to create a synthetic w-band measurement based on the observed mean r & i photometry + + // from Tonry + // (w-r)_obs = 0.042 + 0.166 (r-i)_obs - 0.398 (r-i)_obs^2, (r-i)_obs < 0.5 + // (w-r)_obs = 0.268 - 0.435 (r-i)_obs - 0.078 (r-i)_obs^2, (r-i)_obs > 0.5 + // thus, for (r-i < 0.5): + // (r - i < 0.5) : w = r + 0.042 + 0.166*(r-i) - 0.398(r-i)^2 + // (r - i > 0.5) : w = r + 0.268 - 0.435 (r-i) - 0.078(r-i)^2 + + float Mr = secfilt[NSr].M; + float Mi = secfilt[NSi].M; + + if (!isfinite(Mr)) return FALSE; + if (!isfinite(Mi)) return FALSE; + + // apply a color correction + // XXX this is very GPC1 specific and hard-wired -- be very afraid! + // select the filter; default to fixed photcode and mag limit otherwise + + float Mri = Mr - Mi; + + // saturate at some valid range limits + Mri = MAX (MIN(Mri, 2.0), -0.2); + + float Mw = NAN; + if (Mri < 0.5) { + Mw = Mr + 0.042 + 0.166*(Mri) - 0.398*SQ(Mri); + } else { + Mw = Mr + 0.268 - 0.435*(Mri) - 0.078*SQ(Mri); + } + + dvo_measureT_init (measure); + + measure->M = Mw; + measure->dM = 0.001; + measure->photcode = W_SYNTH; + measure->dbFlags |= ID_MEAS_SYNTH_MAG; + + (*Nmeasure) ++; + (*Nm) ++; + + return TRUE; +}