Index: trunk/Ohana/src/relphot/Makefile =================================================================== --- trunk/Ohana/src/relphot/Makefile (revision 41601) +++ trunk/Ohana/src/relphot/Makefile (revision 41647) @@ -30,4 +30,6 @@ $(SRC)/ConfigInit.$(ARCH).o \ $(SRC)/GridOps.$(ARCH).o \ +$(SRC)/GridIO.$(ARCH).o \ +$(SRC)/ResetOps.$(ARCH).o \ $(SRC)/ImageOps.$(ARCH).o \ $(SRC)/ImageSubset.$(ARCH).o \ @@ -39,8 +41,4 @@ $(SRC)/Shutdown.$(ARCH).o \ $(SRC)/StarOps.$(ARCH).o \ -$(SRC)/WarpImageMaps.$(ARCH).o \ -$(SRC)/StackImageMaps.$(ARCH).o \ -$(SRC)/RepairWarpMeasures.$(ARCH).o \ -$(SRC)/myIndex.$(ARCH).o \ $(SRC)/psps_ids.$(ARCH).o \ $(SRC)/args.$(ARCH).o \ @@ -93,4 +91,6 @@ $(SRC)/ConfigInit.$(ARCH).o \ $(SRC)/GridOps.$(ARCH).o \ +$(SRC)/GridIO.$(ARCH).o \ +$(SRC)/ResetOps.$(ARCH).o \ $(SRC)/ImageOps.$(ARCH).o \ $(SRC)/ImageSubset.$(ARCH).o \ @@ -102,8 +102,4 @@ $(SRC)/Shutdown.$(ARCH).o \ $(SRC)/StarOps.$(ARCH).o \ -$(SRC)/WarpImageMaps.$(ARCH).o \ -$(SRC)/StackImageMaps.$(ARCH).o \ -$(SRC)/RepairWarpMeasures.$(ARCH).o \ -$(SRC)/myIndex.$(ARCH).o \ $(SRC)/psps_ids.$(ARCH).o \ $(SRC)/bcatalog.$(ARCH).o \ Index: trunk/Ohana/src/relphot/include/relphot.h =================================================================== --- trunk/Ohana/src/relphot/include/relphot.h (revision 41601) +++ trunk/Ohana/src/relphot/include/relphot.h (revision 41647) @@ -5,11 +5,14 @@ # include -/* # define GRID_V1 */ -# define GRID_V2 -# define NO_IMAGE -100 - # define ID_SECF_STACK_PRIMARY 0x00004000 # define NBOOTSTRAP 100 + +// MEAS_BAD is used in ImageOps.c, MosaicOps.c, TGroupOps.c to skip specific measurements. +// ID_MEAS_AREA is raised for measurements outside the user-specified region of the chips +// ID_MEAS_NOCAL is raised for measurements outside the valid time range or for inactive photcodes +// Note: these bits are only raised on the temporary MeasureTiny elements and are not saved +// to the database files. +# define MEAS_BAD (ID_MEAS_AREA | ID_MEAS_NOCAL) # ifndef MAX_INT @@ -53,9 +56,9 @@ typedef enum { - MODE_NONE = 0, - MODE_LOAD = 1, - MODE_UPDATE = 2, + MODE_NONE = 0, + MODE_LOAD = 1, + MODE_UPDATE = 2, MODE_UPDATE_OBJECTS = 3, - MODE_SYNTH_PHOT = 4, + MODE_SYNTH_PHOT = 4, } ModeType; @@ -70,4 +73,10 @@ PS1_w = 6, } PS1filters; + +typedef enum { + GRID_MEAN = 0, + GRID_STDEV = 1, + GRID_NPTS = 2, +} GridOutputMode; typedef struct { @@ -99,5 +108,5 @@ char skipCal; // if TRUE, this mosaic is incomplete and should not be calibrated char inTGroup; - Coords coords; + Coords coords; // this is only used to set the mosaic center for assign_images used by region hosts } Mosaic; @@ -140,4 +149,16 @@ int nCode; } TGTimes; + +typedef struct { + unsigned short photcode; + float **Mgrid; // grid of average corrections + float **dMgrid; // grid of stdev of corrections + int **nMgrid; // grid of number of stars for corrections + int Nx; // bin = int(Xchip * (Nx / NxChip)) + int Ny; // bin = int(Ychip * (Ny / NyChip)) + float dX; // bin = int(Xchip * dX), dX = Nx / NxChip + float dY; // bin = int(Ychip * dY), dY = Ny / NyChip + // NxChip, NyChip = 4900,4900 for now +} GridCorrectionType; typedef enum { @@ -366,4 +387,5 @@ int SKY_DEPTH; /** XXX EAM : depth of catalog tables, fix usage */ char *SYNTH_ZERO_POINTS; +char *GRID_MEANFILE; // globals for parallel region operations @@ -389,5 +411,7 @@ double IMAGE_OFFSET; double NIGHT_SCATTER; -double NIGHT_OFFSET; +double NIGHT_CHISQ; +double MOSAIC_SCATTER; +double MOSAIC_CHISQ; double STAR_SCATTER; double STAR_CHISQ; @@ -413,4 +437,10 @@ int MOSAIC_ZEROPT; int TGROUP_ZEROPT; +int GRID_ZEROPT; + +int GRID_BIN_HSC; +int GRID_BIN_CFH; +int GRID_BIN_GPC1; +int GRID_BIN_GPC2; int FREEZE_IMAGES; @@ -424,9 +454,10 @@ int MANUAL_ITERATION; int NLOOP; -int NGRID; int RESET; int RESET_ZEROPTS; +int RESET_FLATCORR; int REPAIR_WARPS; int PRESERVE_PS1; +int REQUIRE_PSFFIT; int UPDATE; int SAVE_IMAGE_UPDATES; @@ -436,13 +467,10 @@ char MOSAICNAME[256]; char STATMODE[32]; -int STAR_BAD; -int MEAS_BAD; int STAR_TOOFEW; int GRID_TOOFEW; int IMAGE_TOOFEW; double IMAGE_GOOD_FRACTION; -int IMAGE_BAD; int CALIBRATE_STACKS_AND_WARPS; -int USE_GRID; + int KEEP_UBERCAL; char *OUTROOT; @@ -457,8 +485,4 @@ char *PhotcodeList; -int RELPHOT_GRID_X; -int RELPHOT_GRID_Y; -int RELPHOT_GRID_BINNING; - int *photseclist; int Nphotcodes; @@ -490,4 +514,7 @@ char *UserCatalog; +enum {GRID_ZPT_MODE_NONE, GRID_ZPT_MODE_ALL, }; +int GRID_ZPT_MODE; + enum {TGROUP_ZPT_MODE_NONE, TGROUP_ZPT_MODE_GOOD_NIGHT, TGROUP_ZPT_MODE_ALL, }; int TGROUP_ZPT_MODE; @@ -504,12 +531,4 @@ RelphotStages STAGES; - -# ifdef GRID_V1 -int setGridMeasure (off_t meas, int cat, double X, double Y); -# endif - -# ifdef GRID_V2 -int setGridMeasure (off_t meas, int cat, double X, double Y, int ccdnum); -# endif /*** relphot prototypes ***/ @@ -523,5 +542,5 @@ int bcatalog PROTO((Catalog *subcatalog, Catalog *catalog, int Ncat)); void clean_images PROTO((void)); -void clean_measures PROTO((Catalog *catalog, int Ncatalog, int final, FlatCorrectionTable *flatcorr)); +void clean_measures PROTO((Catalog *catalog, int Ncatalog, int final)); void clean_mosaics PROTO((void)); void clean_tgroups PROTO((void)); @@ -549,5 +568,5 @@ int SetSignals (void); -void freeGridBins PROTO((int Ncatalog)); +void freeGridBins PROTO((void)); void freeImageBins PROTO((int Ncatalog, int doImageList)); void freeMosaicBins PROTO((int Ncatalog, int doMosaicList)); @@ -560,6 +579,7 @@ float getMcal PROTO((off_t meas, int cat, dvoMagClassType class)); -float getMflat PROTO((off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog)); -float getMgrid PROTO((off_t meas, int cat)); +float getMflat PROTO((off_t meas, int cat, Catalog *catalog)); +float getMgridTiny PROTO((MeasureTiny *measure)); +float getMgrid PROTO((Measure *measure)); float getMmos PROTO((off_t meas, int cat)); float getMgrp PROTO((off_t meas, int cat, float airmass, float *dZpt)); @@ -570,12 +590,18 @@ Image *getimages PROTO((off_t *N, off_t **LineNumber)); ImageSubset *getimages_subset PROTO((off_t *N)); -void global_stats PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr, int nloop)); -void initGrid PROTO((int dX, int dY)); -void initGridBins PROTO((Catalog *catalog, int Ncatalog)); +void global_stats PROTO((Catalog *catalog, int Ncatalog, int nloop)); +void initGrid PROTO((void)); +void setMflatFromGrid PROTO((Catalog *catalog)); +void initGridBins PROTO((void)); +GridCorrectionType *getGridCorrByCode PROTO((int code)); +GridCorrectionType *getGridCorrNext PROTO((int *Nlast)); +int GridCorrectionSave PROTO((void)); +int GridCorrectionSaveFile PROTO((char *filename, GridOutputMode mode)); +void GridCorrectionLoad PROTO((char *filename)); void initImageBins PROTO((Catalog *catalog, int Ncatalog, int doImageList)); void initImagesSubset PROTO((ImageSubset *input, off_t *line_number, off_t N)); void initImages PROTO((Image *input, off_t *LineNumber, off_t N)); void initMosaicBins PROTO((Catalog *catalog, int Ncatalog, int doMosaicList)); -void initMosaicGrid PROTO((Image *image, off_t Nimage)); +void initMosaicMcal PROTO((Image *image, off_t Nimage)); void initMosaics PROTO((Image *subset, off_t Nsubset, Image *image, char *inSubset, off_t Nimage)); void initTGroups PROTO((Image *subset, off_t Nsubset)); @@ -621,10 +647,9 @@ void plot_chisq PROTO((Catalog *catalog, int Ncatalog)); void plot_defaults PROTO((Graphdata *graphdata)); -void plot_grid PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr)); void plot_images PROTO((void)); void plot_list PROTO((Graphdata *graphdata, double *xlist, double *ylist, int N, char *label, char *format, ...) OHANA_FORMAT(printf, 6, 7) ); void plot_mosaic_fields PROTO((Catalog *catalog)); void plot_mosaics PROTO((void)); -void plot_scatter PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr)); +void plot_scatter PROTO((Catalog *catalog, int Ncatalog)); void plot_star_coords PROTO((Catalog *catalog, int Ncatalog)); void plot_stars PROTO((Catalog *catalog, int Ncatalog)); @@ -635,27 +660,29 @@ int open_graph PROTO((int N)); -void reload_catalogs PROTO((SkyList *skylist, FlatCorrectionTable *flatcorr, int hostID, char *hostpath)); +void reload_catalogs PROTO((SkyList *skylist, int hostID, char *hostpath)); int reload_catalogs_parallel PROTO((SkyList *sky)); int reload_images PROTO((FITS_DB *db)); int setExclusions PROTO((Catalog *catalog, int Ncatalog, int verbose)); -void setMgrid PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr)); +void setMgrid PROTO((Catalog *catalog, int Ncatalog)); +void resetMgrid PROTO((void)); + void setMcalFromMosaics PROTO((void)); void setMcalFromTGroups PROTO((void)); -int setMcalOutput PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr)); - -void setMcal PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr)); -int setMmos PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr)); -int setMgrp PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr)); - -int setMrel PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr)); -void setMrelFinal PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr, int simpleAverage)); -int setMrelOutput PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr)); - -void setMcalTest PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr)); +int setMcalOutput PROTO((Catalog *catalog, int Ncatalog)); + +void setMcal PROTO((Catalog *catalog)); +int setMmos PROTO((Catalog *catalog)); +int setMgrp PROTO((Catalog *catalog)); + +int setMrel PROTO((Catalog *catalog, int Ncatalog)); +void setMrelFinal PROTO((Catalog *catalog, int simpleAverage)); +int setMrelOutput PROTO((Catalog *catalog, int Ncatalog)); + +void setMcalTest PROTO((Catalog *catalog)); int setMave PROTO((Catalog *catalog, int Ncatalog)); void set_ZP PROTO((double ZERO)); int setrefcode PROTO((Image *image, off_t Nimage)); -void skip_measurements PROTO((Catalog *catalog, int pass, FlatCorrectionTable *flatcorr)); +void skip_measurements PROTO((Catalog *catalog, int pass)); void sortA PROTO((double *X, int N)); void sortB PROTO((double *X, double *Y, int N)); @@ -663,9 +690,9 @@ void sortD PROTO((double *X, double *Y, double *Z, int N)); StatType statsImageM PROTO((Catalog *catalog)); -StatType statsImageN PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr)); +StatType statsImageN PROTO((Catalog *catalog)); StatType statsImageX PROTO((Catalog *catalog)); StatType statsImagedM PROTO((Catalog *catalog)); StatType statsMosaicM PROTO((Catalog *catalog)); -StatType statsMosaicN PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr)); +StatType statsMosaicN PROTO((Catalog *catalog)); StatType statsMosaicX PROTO((Catalog *catalog)); StatType statsMosaicdM PROTO((Catalog *catalog)); @@ -674,5 +701,5 @@ StatType statsTGroupdM PROTO((Catalog *catalog)); -StatType statsStarN PROTO((Catalog *catalog, int Ncatalog, int Nsec, int seccode, FlatCorrectionTable *flatcorr)); +StatType statsStarN PROTO((Catalog *catalog, int Ncatalog, int Nsec, int seccode)); StatType statsStarS PROTO((Catalog *catalog, int Ncatalog, int Nsec)); StatType statsStarX PROTO((Catalog *catalog, int Ncatalog, int Nsec)); @@ -720,8 +747,10 @@ int print_measure_set_alt (Average *average, SecFilt *secfilt, Measure *measure); -int setMrel_catalog_alt (Catalog *catalog, int Nc, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt); -int setMrelAverageExposure (Catalog *catalog, int cat, off_t ave, int Nsecfilt, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results); -int setMrelAverageStack (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr); -int setMrelAverageForcedWarp (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr, SetMrelInfo *results); + +// in setMrelCatalog.c: +int setMrelCatalog (Catalog *catalog, int Nc, int isSetMrelFinal, SetMrelInfo *results, int Nsecfilt); +int setMrelAverageExposure (Catalog *catalog, int cat, off_t ave, int Nsecfilt, int isSetMrelFinal, SetMrelInfo *results); +int setMrelAverageStack (Catalog *catalog, int cat, off_t ave, int Nsecfilt); +int setMrelAverageForcedWarp (Catalog *catalog, int cat, off_t ave, int Nsecfilt, SetMrelInfo *results); int setGlobalObjStats (Average *average, Measure *measure); @@ -735,5 +764,5 @@ StatDataSet *StatDataSetAlloc (int Nsecfilt, int Nmax); -int init_synthetic_mags (); +int init_synthetic_mags (void); int add_synthetic_mags (AverageTiny *average, SecFilt *secfilt, MeasureTiny *measure, off_t *Nmeasure, off_t *Nm); @@ -796,5 +825,5 @@ int SynthZeroPointsLoad (char *filename); -SynthZeroPoints *SynthZeroPointsGet (); +SynthZeroPoints *SynthZeroPointsGet (void); int relphot_synthphot (SkyList *sky, int hostID, char *hostpath); @@ -815,18 +844,5 @@ void put_astrom_table (AstromOffsetTable *myTable); -void free_astrom_table (); - -int RepairWarpMeasures (Catalog *catalog); - -int FindWarpGroups (void); -int GetWarpSeq (Image *image, int obstime, unsigned short photcode, double Rave, double Dave, float X, float Y); -void FreeWarpGroups (void); - -myIndexType *myIndexInit (); -int myIndexFree (myIndexType *myIndex); -int myIndexUpdateLimits (myIndexType *myIndex, int value); -int myIndexSetRange (myIndexType *myIndex); -int myIndexSetEntry (myIndexType *myIndex, int value, int entry); -int myIndexGetEntry (myIndexType *myIndex, int value); +void free_astrom_table (void); uint64_t CreatePSPSObjectID(double ra, double dec); @@ -834,11 +850,5 @@ uint64_t CreatePSPSDetectionID(double tobs, int ccdid, int detID); -int MakeStackIndex (void); -void FreeStackGroups (void); -int GetStackSeq (Image *image, double Rstk, double Dstk, unsigned short photcode, float X, float Y); void sort_by_ra (double *R, double *D, int *I, int *S, int N); - -void RepairWarpMeasuresOpenLogfile (); -void RepairWarpMeasuresCloseLogfile (); void dump_tgroups (Catalog *catalog, int Npass); @@ -849,5 +859,5 @@ int GetZptIteration (void); -void SetZeroPointModes (void); +void SetZeroPointModes (Catalog *catalog, int Ncatalog); int UseStandardOLS (ZptFitModeType mode); int GetActivePhotcodeIndex (int photcode); @@ -857,2 +867,10 @@ int MagResidSave(char *filename, Catalog *catalog); + +void ResetAverageAndMeasure (Catalog *catalog); +void ResetAverageObjects (Catalog *catalog); +void ResetImages (Image *subset, off_t Nsubset); +void ResetMeasureZeroPoints (MeasureTiny *measure, off_t Nmeasure, off_t Ncat); +void ResetAverageActivePhotcodes (SecFilt *secfilt); + + Index: trunk/Ohana/src/relphot/src/BrightCatalog.c =================================================================== --- trunk/Ohana/src/relphot/src/BrightCatalog.c (revision 41601) +++ trunk/Ohana/src/relphot/src/BrightCatalog.c (revision 41647) @@ -12,8 +12,4 @@ myAssert (!strcmp(type, #TYPE), "wrong column type"); -// XXX double check the header extname values for each table -// XXX make sure we free things as we can -// XXX make sure we close files as we can -// XXX handle free and close on error return as well BrightCatalog *BrightCatalogLoad(char *filename) { @@ -622,12 +618,5 @@ int BrightCatalogSplitFree (CatalogSplitter *catalogs) { - // XXX don't free the catalogs : - // for (i = 0; i < NCATALOG; i++) { - // free (catalogs->catalog[i].averageT); - // free (catalogs->catalog[i].measureT); - // free (catalogs->catalog[i].secfilt); - // } - // free (catalogs->catalog); - + // don't free the catalogs : free elsewhere free (catalogs->catIDs); free (catalogs->NAVERAGE); @@ -652,5 +641,4 @@ catIDmax = MAX(catIDmax, bcatalog->average[i].catID); } - // XXX validate the measure ID range here int maxIDold = catalogs->maxID; @@ -746,5 +734,4 @@ // assign the measures to the corresponding catalog - // XXX what about averef and related links? Do I need them? for (i = 0; i < bcatalog->Nmeasure; i++) { int ID = bcatalog->measure[i].catID; Index: trunk/Ohana/src/relphot/src/ConfigInit.c =================================================================== --- trunk/Ohana/src/relphot/src/ConfigInit.c (revision 41601) +++ trunk/Ohana/src/relphot/src/ConfigInit.c (revision 41647) @@ -26,11 +26,16 @@ GetConfig (config, "SIGMA_LIM", "%lf", 0, &SIGMA_LIM); + DefConfig ("STAR_SCATTER", "%lf", 0.1, STAR_SCATTER); + DefConfig ("STAR_CHISQ", "%lf", 10.0, STAR_CHISQ); - GetConfig (config, "STAR_SCATTER", "%lf", 0, &STAR_SCATTER); + DefConfig ("NIGHT_SCATTER", "%lf", 0.1, NIGHT_SCATTER); + DefConfig ("NIGHT_CHISQ", "%lf", 10.0, NIGHT_CHISQ); + + DefConfig ("MOSAIC_SCATTER", "%lf", 0.1, MOSAIC_SCATTER); + DefConfig ("MOSAIC_CHISQ", "%lf", 10.0, MOSAIC_CHISQ); GetConfig (config, "IMAGE_SCATTER", "%lf", 0, &IMAGE_SCATTER); GetConfig (config, "IMAGE_OFFSET", "%lf", 0, &IMAGE_OFFSET); - GetConfig (config, "STAR_CHISQ", "%lf", 0, &STAR_CHISQ); GetConfig (config, "GRID_TOOFEW", "%d", 0, &GRID_TOOFEW); GetConfig (config, "STAR_TOOFEW", "%d", 0, &STAR_TOOFEW); @@ -54,6 +59,4 @@ DefConfig ("RELPHOT_IMFIT_SYS_SIGMA_LIM", "%lf", 0.01, IMFIT_SYS_SIGMA_LIM); - DefConfig ("NIGHT_SCATTER", "%lf", 0.05, NIGHT_SCATTER); - DefConfig ("NIGHT_OFFSET", "%lf", 0.50, NIGHT_OFFSET); DefConfig ("SKY_DEPTH", "%d", 2, SKY_DEPTH); @@ -63,8 +66,4 @@ GetConfig (config, "ZERO_PT", "%lf", 0, &ZERO_POINT); - - GetConfig (config, "RELPHOT_GRID_X", "%d", 0, &RELPHOT_GRID_X); - GetConfig (config, "RELPHOT_GRID_Y", "%d", 0, &RELPHOT_GRID_Y); - GetConfig (config, "RELPHOT_GRID_BINNING", "%d", 0, &RELPHOT_GRID_BINNING); GetConfig (config, "CAMERA_CONFIG", "%s", 0, CameraConfig); GetConfig (config, "MOSAICNAME", "%s", 0, MOSAICNAME); Index: trunk/Ohana/src/relphot/src/GridIO.c =================================================================== --- trunk/Ohana/src/relphot/src/GridIO.c (revision 41647) +++ trunk/Ohana/src/relphot/src/GridIO.c (revision 41647) @@ -0,0 +1,181 @@ +# include "relphot.h" + +int GridCorrectionSave () { + + if (!GRID_ZEROPT) return TRUE; + + char filename[1024], uniquer[12]; + int TIME = time(NULL); + int PID = getpid(); + + snprintf_nowarn (uniquer, 12, "%05d.%05d", PID, TIME % 100000); + + snprintf (filename, 1024, "%s/gridcorr.mean.%s.fits", CATDIR, uniquer); + if (!GridCorrectionSaveFile (filename, GRID_MEAN)) return FALSE; + GRID_MEANFILE = strcreate (filename); // save in the global so we can pass to e.g., reload_catalogs + + snprintf (filename, 1024, "%s/gridcorr.stdev.%s.fits", CATDIR, uniquer); + if (!GridCorrectionSaveFile (filename, GRID_STDEV)) return FALSE; + + snprintf (filename, 1024, "%s/gridcorr.npts.%s.fits", CATDIR, uniquer); + if (!GridCorrectionSaveFile (filename, GRID_NPTS)) return FALSE; + + return TRUE; +} + +int GridCorrectionSaveFile (char *filename, GridOutputMode mode) { + + /* open file for input */ + FILE *f = fopen (filename, "w"); + if (!f) { + gprint (GP_ERR, "can't open file for write : %s\n", filename); + return FALSE; + } + + // generate a PHU + + Header phu_header; + Matrix phu_matrix; + + gfits_init_header (&phu_header); + gfits_create_header (&phu_header); + gfits_modify_alt (&phu_header, "EXTEND", "%t", 1, TRUE); + gfits_create_matrix (&phu_header, &phu_matrix); + gfits_fwrite_header (f, &phu_header); + gfits_fwrite_matrix (f, &phu_matrix); + gfits_free_header (&phu_header); + gfits_free_matrix (&phu_matrix); + + int Nlast = -1; + + while (TRUE) { + // getGridCorreNext() takes a photcode and returns the *next* valid grid correction pointer (start at -1) + GridCorrectionType *GridCorr = getGridCorrNext (&Nlast); + if (!GridCorr) break; + + int filter = GetPhotcodeEquivCodebyCode (GridCorr->photcode); + char *filtname = GetPhotcodeNamebyCode (filter); // reference, do not free + char *photname = GetPhotcodeNamebyCode (GridCorr->photcode); // reference, do not free + + Header header; + Matrix matrix; + + gfits_init_header (&header); + header.bitpix = (mode == GRID_NPTS) ? 32 : -32; + header.Naxes = 2; + header.Naxis[0] = GridCorr->Nx; + header.Naxis[1] = GridCorr->Ny; + + /* create the appropriate header and matrix */ + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + + // create header & matrix + if (!gfits_modify (&header, "PHOTCODE", "%d", 1, GridCorr->photcode)) myAbort ("failed to set PHOTCODE"); + if (!gfits_modify (&header, "FILTER", "%s", 1, filtname)) myAbort ("failed to set FILTER"); + if (!gfits_modify (&header, "SEASON", "%d", 1, 0)) myAbort ("failed to set SEASON"); // XXX hard-wired for now + if (!gfits_modify (&header, "X_CHIP", "%d", 1, -1)) myAbort ("failed to set X_CHIP"); + if (!gfits_modify (&header, "Y_CHIP", "%d", 1, -1)) myAbort ("failed to set Y_CHIP"); + if (!gfits_modify (&header, "DX_CHIP", "%f", 1, GridCorr->dX)) myAbort ("failed to set DX_CHIP"); + if (!gfits_modify (&header, "DY_CHIP", "%f", 1, GridCorr->dY)) myAbort ("failed to set DY_CHIP"); + if (!gfits_modify (&header, "NX_CHIP", "%d", 1, GridCorr->Nx)) myAbort ("failed to set NX_CHIP"); + if (!gfits_modify (&header, "NY_CHIP", "%d", 1, GridCorr->Ny)) myAbort ("failed to set NY_CHIP"); + + int Nx = GridCorr->Nx; + + float *fvalue = (float *) matrix.buffer; + int *ivalue = (int *) matrix.buffer; + + for (int iy = 0; iy < GridCorr->Ny; iy++) { + for (int ix = 0; ix < GridCorr->Nx; ix++) { + switch (mode) { + case GRID_MEAN: fvalue[ix + iy*Nx] = GridCorr-> Mgrid[ix][iy]; break; + case GRID_STDEV: fvalue[ix + iy*Nx] = GridCorr->dMgrid[ix][iy]; break; + case GRID_NPTS: ivalue[ix + iy*Nx] = GridCorr->nMgrid[ix][iy]; break; + } + } + } + + char extname[85]; + switch (mode) { + // XXX If we add seasons, add here: + case GRID_MEAN: snprintf (extname, 85, "%s.%s", photname, "MEAN"); break; + case GRID_STDEV: snprintf (extname, 85, "%s.%s", photname, "STDEV"); break; + case GRID_NPTS: snprintf (extname, 85, "%s.%s", photname, "NPTS"); break; + } + gfits_modify (&header, "EXTNAME", "%s", 1, extname); + + gfits_primary_to_extended (&header, "IMAGE", "Image Extension"); + gfits_fwrite_header (f, &header); + gfits_fwrite_matrix (f, &matrix); + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + } + fclose (f); + + return TRUE; +} + +void GridCorrectionLoad (char *filename) { + + int photcode; + + Header phu; + Header header; + Matrix matrix; + + if (!GRID_ZEROPT) return; + + // allocate the basic containing structures: + initGridBins (); + + /* open file for input */ + FILE *f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + gprint (GP_ERR, "can't open file for read : %s\n", filename); + return; + } + + // read the PHU header + if (!gfits_fread_header (f, &phu)) return; + + // move to the first extension + int Nbytes = gfits_data_size (&phu); + fseeko (f, Nbytes, SEEK_CUR); + + gfits_free_header (&phu); + + while (gfits_fread_header (f, &header)) { + if (!gfits_fread_matrix (f, &matrix, &header)) myAbort ("failed to read matrix"); + + if (!gfits_scan (&header, "PHOTCODE", "%d", 1, &photcode)) myAbort ("failed to get PHOTCODE"); + GridCorrectionType *GridCorr = getGridCorrByCode (photcode); + GridCorr->photcode = photcode; + + if (!gfits_scan (&header, "DX_CHIP", "%f", 1, &GridCorr->dX)) myAbort ("failed to get DX_CHIP"); + if (!gfits_scan (&header, "DY_CHIP", "%f", 1, &GridCorr->dY)) myAbort ("failed to get DY_CHIP"); + if (!gfits_scan (&header, "NX_CHIP", "%d", 1, &GridCorr->Nx)) myAbort ("failed to get NX_CHIP"); + if (!gfits_scan (&header, "NY_CHIP", "%d", 1, &GridCorr->Ny)) myAbort ("failed to get NY_CHIP"); + + GridCorr->dMgrid = NULL; + GridCorr->nMgrid = NULL; + ALLOCATE (GridCorr->Mgrid, float *, GridCorr->Nx); + for (int ix = 0; ix < GridCorr->Nx; ix++) { + ALLOCATE (GridCorr->Mgrid[ix], float, GridCorr->Ny); + } + + int Nx = GridCorr->Nx; + + float *fvalue = (float *) matrix.buffer; + + for (int iy = 0; iy < GridCorr->Ny; iy++) { + for (int ix = 0; ix < GridCorr->Nx; ix++) { + GridCorr-> Mgrid[ix][iy] = fvalue[ix + iy*Nx]; + } + } + gfits_free_header (&header); + gfits_free_matrix (&matrix); + } + gfits_free_header (&header); +} Index: trunk/Ohana/src/relphot/src/GridOps.c =================================================================== --- trunk/Ohana/src/relphot/src/GridOps.c (revision 41601) +++ trunk/Ohana/src/relphot/src/GridOps.c (revision 41647) @@ -1,898 +1,350 @@ # include "relphot.h" -enum { - GRID_FITTED, - GRID_FROZEN, - GRID_REFERENCE, -}; - -static int Ngrid; // number of grid elements (gridX * gridY) -static float *gridM; // magnitude offset for this grid cell -static float *gridS; // stdev of the magnitude offset for this grid cell -static int *gridN; // number of stars used to measure the magnitude offset for this grid cell -static int *gridV; // data mode for this cell: fitted, frozen, reference -static int gridX; // number of grid elements in X direction -static int gridY; // number of grid elements in Y direction - -static int **bin; // link from catalog, measure to grid element -static int **Xmeas; // grid x-coordinate for a measurement -static int **Ymeas; // grid y-coordinate for a measurement - -static int **clist; // link from measurement on a cell to catalog containing measurement -static off_t **mlist; // link from measurement on a cell to measurement in a catalog -static off_t *Nlist; // number of measurements for each grid cell -static off_t *NLIST; // allocated number of measurements for each grid cell - -static struct { - int Nchip; - int Mx, My; /* mosaic size in chips */ - int Nx, Ny; /* chip size in pixels */ - int *Fx, *Fy; /* chip flip */ - int *Ox, *Oy; /* chip offset */ - int *valid; - char **ccdname; -} camera; - -static int *ccdnum; -static char *config; - -void initGrid (int dX, int dY) { - OHANA_UNUSED_PARAM(dX); - OHANA_UNUSED_PARAM(dY); - - int i, N, ccdnum_max, refX, refY, refBin; - char *p, field[64], line[256]; - int *Fx, *Fy; /* chip flip */ - int *Ox, *Oy; /* chip offset */ - char **ccdname; - - /* load camera config file */ - config = LoadConfigFile (CameraConfig); - if (config == (char *) NULL) { - fprintf (stderr, "ERROR: can't find camera config file %s\n", CameraConfig); - exit (1); - } - - /* load basic mosaic parameters */ - ScanConfig (config, "NCCD", "%d", 1, &camera.Nchip); - ScanConfig (config, "MOSAIC_X", "%d", 1, &camera.Mx); - ScanConfig (config, "MOSAIC_Y", "%d", 1, &camera.My); - ScanConfig (config, "NAXIS1", "%d", 1, &camera.Nx); - ScanConfig (config, "NAXIS2", "%d", 1, &camera.Ny); - - ScanConfig (config, "REFCELL.X", "%d", 1, &refX); - ScanConfig (config, "REFCELL.Y", "%d", 1, &refY); - - // temporary storage - ALLOCATE (Fx, int, camera.Nchip); - ALLOCATE (Fy, int, camera.Nchip); - ALLOCATE (Ox, int, camera.Nchip); - ALLOCATE (Oy, int, camera.Nchip); - ALLOCATE (ccdname, char *, camera.Nchip); - ALLOCATE (ccdnum, int, camera.Nchip); - - /* load per-chip parameters */ - for (i = 0; i < camera.Nchip; i++) { - ALLOCATE (ccdname[i], char, 256); - sprintf (field, "CCD.%d", i); - ScanConfig (config, field, "%s", 1, line); - // XXX get error status! - - sscanf (line, "%s %d %d %d %d", ccdname[i], &Ox[i], &Oy[i], &Fx[i], &Fy[i]); +/* + We define a 'grid correction', essentially the flat-field correction, as a correction per + photcode subdivided into an NxN array. The dimensions of the chips corresponding to a + photcode will need to be added to the photcode table. This means a schema update, + which I detest. For now (2021.05.16), I will hard-wire the GPC1 / GPC2 chip size and + worry about HSC & Megacam in the future. + + We have a collection of photcodes, with photcodeID limited by design to 64k (unsigned short). + + Thus we can generate an array of pointers to the grid correction structures and access them + by photcode. + + + */ + +static GridCorrectionType **GridCorr = NULL; +static int NGridCorr = 0; + +# if (0) +/* PS1 / PS2 values */ +# define NX_CHIP 4900 +# define NY_CHIP 4900 +# define NX_BIN 16 +# define NY_BIN 16 +# endif + +/* test values */ +# define NX_CHIP_DEFAULT 1000 +# define NY_CHIP_DEFAULT 1000 +# define NX_BIN_DEFAULT 2 +# define NY_BIN_DEFAULT 2 + +void initGridBins (void) { + + // allocate the full possible range of GridCorrectionType pointers, Nphotcode + // loop over all images to find actual existing photcodes + // generate initial grid values for each existing photcode + + if (!GRID_ZEROPT) return; // skip if we are ignoring the grid correction + + if (GridCorr) return; + + PhotCodeData *photcodes = GetPhotcodeTable(); + if (!photcodes) return; + + // we have photcodes->Ncodes actually loaded. loop over them and allocate + // array only as large as the max photcodes->code[i].code + + int maxCode = 0; + for (int i = 0; i < photcodes->Ncode; i++) { + maxCode = MAX(maxCode, photcodes->code[i].code); + myAssert (maxCode < 0x10000, "oops"); + } + + NGridCorr = maxCode + 1; + ALLOCATE (GridCorr, GridCorrectionType *, NGridCorr); + for (int i = 0; i < NGridCorr; i++) { + GridCorr[i] = NULL; + } + + fprintf (stderr, "Generating grid corrections for %d photcodes\n", NGridCorr); + return; +} + +void freeGridBins() { + + if (!GridCorr) return; + + for (int code = 0; code < NGridCorr; code++) { + if (!GridCorr[code]) continue; + + for (int ix = 0; ix < GridCorr[code]->Nx; ix++) { + FREE (GridCorr[code]-> Mgrid[ix]); + if (GridCorr[code]->dMgrid) FREE (GridCorr[code]->dMgrid[ix]); + if (GridCorr[code]->nMgrid) FREE (GridCorr[code]->nMgrid[ix]); + } + + FREE (GridCorr[code]-> Mgrid); + FREE (GridCorr[code]->dMgrid); + FREE (GridCorr[code]->nMgrid); + + FREE(GridCorr[code]); + } + + FREE (GridCorr); + GridCorr = NULL; + + return; +} + +GridCorrectionType *getGridCorrNext (int *Nlast) { + + if (GridCorr == NULL) return NULL; + + if (*Nlast >= NGridCorr) return NULL; + + if (*Nlast < 0) *Nlast = -1; + + GridCorrectionType *result = NULL; + for (int i = *Nlast + 1; i < NGridCorr; i++) { + if (GridCorr[i] == NULL) continue; + result = GridCorr[i]; + *Nlast = i; + break; + } + return result; +} + +GridCorrectionType *getGridCorrByCode (int code) { + + if (GridCorr == NULL) return NULL; + + if (code >= NGridCorr) return NULL; + if (code < 0) return NULL; + + ALLOCATE(GridCorr[code], GridCorrectionType, 1) + GridCorrectionType *result = GridCorr[code]; + return result; +} + +/* for GPC1, we have 60 chips, 5 filters, 16x16 grid cells = 2MB of memory for this stuff + if we go to 64x64 grid cells (~75 pixels), then it is still only 29MB */ +void initGrid (void) { + + if (!GRID_ZEROPT) return; + + off_t Nimages = 0; + Image *images = getimages (&Nimages, NULL); + + int NGridReal = 0; + + for (int i = 0; i < Nimages; i++) { + int code = images[i].photcode; + myAssert (code >= 0, "oops"); + myAssert (code < NGridCorr, "oops"); - p = ccdname[i]; - while (!isdigit(*p) && *p) p++; - if (*p == 0) continue; - ccdnum[i] = atoi (p); - } - - /* we now have the parameters loaded into a minimal length list; reshuffle to the a list of length 0 - MAX(ccdnum) */ - ccdnum_max = 0; - for (i = 0; i < camera.Nchip; i++) { - ccdnum_max = MAX(ccdnum_max, ccdnum[i]); - } - ccdnum_max ++; - - if (ccdnum_max < camera.Nchip) { - fprintf (stderr, "problem with camera config: duplicate ccd IDs\n"); - exit (1); - } - - if (ccdnum_max > 0x1000) { - fprintf (stderr, "problem with camera config: absurd max ccd ID number %d\n", ccdnum_max); - exit (1); - } - - ALLOCATE (camera.Fx, int, ccdnum_max); - ALLOCATE (camera.Fy, int, ccdnum_max); - ALLOCATE (camera.Ox, int, ccdnum_max); - ALLOCATE (camera.Oy, int, ccdnum_max); - ALLOCATE (camera.valid, int, ccdnum_max); - ALLOCATE (camera.ccdname, char *, ccdnum_max); - - for (i = 0; i < ccdnum_max; i++) { - camera.valid[i] = FALSE; - } - - for (i = 0; i < camera.Nchip; i++) { - N = ccdnum[i]; - camera.Fx[N] = Fx[i]; - camera.Fy[N] = Fy[i]; - camera.Ox[N] = Ox[i]; - camera.Oy[N] = Oy[i]; - camera.ccdname[N] = ccdname[i]; - camera.valid[N] = TRUE; - } - - /* define mosaic 2d correction grid: - * GRID_X is the number of grid pixels per chip in the X direction */ - gridX = RELPHOT_GRID_X * camera.Mx; - gridY = RELPHOT_GRID_Y * camera.My; - Ngrid = gridX * gridY; - - ALLOCATE (gridM, float, Ngrid); - ALLOCATE (gridS, float, Ngrid); - ALLOCATE (gridN, int, Ngrid); - ALLOCATE (gridV, int, Ngrid); - - // the grid bins may have one of three possible states: fitted, frozen, reference - // set the initial values to indicate that the bins are frozen - for (i = 0; i < Ngrid; i++) { - gridM[i] = 0.0; - gridS[i] = 0.0; - gridN[i] = 0; - gridV[i] = GRID_FROZEN; - } - - // refBin is the index of the grid cell which is kept at 0.0 (all others are relative to this) - refBin = refX + refY*gridX; - gridV[refBin] = GRID_REFERENCE; -} - -void initGridBins (Catalog *catalog, int Ncatalog) { - - int i, j; - - if (!USE_GRID) return; - - /* define cat,meas -> grid pointers */ - ALLOCATE (bin, int *, Ncatalog); - ALLOCATE (Xmeas, int *, Ncatalog); - ALLOCATE (Ymeas, int *, Ncatalog); - for (i = 0; i < Ncatalog; i++) { - ALLOCATE (bin[i], int, MAX (catalog[i].Nmeasure, 1)); - ALLOCATE (Xmeas[i], int, MAX (catalog[i].Nmeasure, 1)); - ALLOCATE (Ymeas[i], int, MAX (catalog[i].Nmeasure, 1)); - for (j = 0; j < catalog[i].Nmeasure; j++) bin[i][j] = -1; - } - - /* define grid -> cat,meas pointers */ - ALLOCATE (Nlist, off_t, Ngrid); - ALLOCATE (NLIST, off_t, Ngrid); - ALLOCATE (clist, int *, Ngrid); - ALLOCATE (mlist, off_t *, Ngrid); - - for (i = 0; i < Ngrid; i++) { - Nlist[i] = 0; - NLIST[i] = 100; - ALLOCATE (clist[i], int, NLIST[i]); - ALLOCATE (mlist[i], off_t, NLIST[i]); - } -} - -void freeGridBins (int Ncatalog) { - - int i; - - if (!USE_GRID) return; - - /* define cat,meas -> grid pointers */ - for (i = 0; i < Ncatalog; i++) { - free (bin[i]); - free (Xmeas[i]); - free (Ymeas[i]); - } - free (bin); - free (Xmeas); - free (Ymeas); - - /* define grid -> cat,meas pointers */ - for (i = 0; i < Ngrid; i++) { - free (clist[i]); - free (mlist[i]); - } - free (Nlist); - free (NLIST); - free (clist); - free (mlist); -} - -int setGridMeasure (off_t meas, int cat, double X, double Y, int ccdnum) { - - int ix, iy, Cx, Cy, i; - double x, y; - - /* X, Y are chip coords on chip ccdnum */ - - /* normalize X & Y */ - x = X; - if (camera.Fx[ccdnum]) x = camera.Nx - X; - y = Y; - if (camera.Fy[ccdnum]) y = camera.Ny - Y; - - /* grid coords on the chip */ - Cx = MIN (MAX ((x / camera.Nx) * RELPHOT_GRID_X, 0), RELPHOT_GRID_X - 1); - Cy = MIN (MAX ((y / camera.Ny) * RELPHOT_GRID_Y, 0), RELPHOT_GRID_Y - 1); - - /* coordinates in the grid */ - ix = Cx + camera.Ox[ccdnum]*RELPHOT_GRID_X; - iy = Cy + camera.Oy[ccdnum]*RELPHOT_GRID_Y; - - i = ix + iy*gridX; - - bin[cat][meas] = i; - Xmeas[cat][meas] = x + camera.Ox[ccdnum]*camera.Nx; - Ymeas[cat][meas] = y + camera.Oy[ccdnum]*camera.Ny; - clist[i][Nlist[i]] = cat; - mlist[i][Nlist[i]] = meas; - - // for the moment, add up the total grid count - gridN[i]++; - - Nlist[i] ++; - if (Nlist[i] == NLIST[i]) { - NLIST[i] += 100; - REALLOCATE (clist[i], int, NLIST[i]); - REALLOCATE (mlist[i], off_t, NLIST[i]); - } - return (TRUE); - - fprintf (stderr, "error: star out of grid\n"); - exit (1); -} - -int showGridCount() { - - int ix, iy, i; - - for (iy = 0; iy < gridY; iy++) { - for (ix = 0; ix < gridX; ix++) { - i = ix + iy*gridX; - fprintf (stderr, "%3d ", gridN[i]); - } - fprintf (stderr, "\n"); + // valid photcodes values (code) are in range 1 <= code < 0x10000 + // photcode == 0 are e.g., PHU (mosaic) images, and should be ignored here + if (!code) continue; + + int NX_CHIP = NX_CHIP_DEFAULT; + int NY_CHIP = NY_CHIP_DEFAULT; + int NX_BIN = NX_BIN_DEFAULT; + int NY_BIN = NY_BIN_DEFAULT; + if (isGPC1chip(code)) { NX_CHIP = 4900; NY_CHIP = 4900; NX_BIN = NY_BIN = GRID_BIN_GPC1; } + if (isGPC2chip(code)) { NX_CHIP = 4900; NY_CHIP = 4900; NX_BIN = NY_BIN = GRID_BIN_GPC2; } + if (isHSCchip(code)) { NX_CHIP = 2100; NY_CHIP = 4200; NX_BIN = NY_BIN = GRID_BIN_HSC; } + if (isCFHchip(code)) { NX_CHIP = 2100; NY_CHIP = 4200; NX_BIN = NY_BIN = GRID_BIN_CFH; } + + if (GridCorr[code]) continue; // already created this one + + NGridReal ++; + + ALLOCATE(GridCorr[code], GridCorrectionType, 1); + GridCorr[code]->photcode = code; + GridCorr[code]->Nx = NX_BIN; + GridCorr[code]->Ny = NY_BIN; + GridCorr[code]->dX = NX_BIN / (float) NX_CHIP; + GridCorr[code]->dY = NY_BIN / (float) NY_CHIP; + + // we are normally accessing this array randomly, so there is no advantage to + // doing Mgrid[y][x] vs Mgrid[x][y] + ALLOCATE (GridCorr[code]-> Mgrid, float *, NX_BIN); + ALLOCATE (GridCorr[code]->dMgrid, float *, NX_BIN); + ALLOCATE (GridCorr[code]->nMgrid, int *, NX_BIN); + for (int ix = 0; ix < NX_BIN; ix++) { + ALLOCATE (GridCorr[code]-> Mgrid[ix], float, NY_BIN); + ALLOCATE (GridCorr[code]->dMgrid[ix], float, NY_BIN); + ALLOCATE (GridCorr[code]->nMgrid[ix], int, NY_BIN); + } + } + resetMgrid(); // start with values of 0 + fprintf (stderr, "Init grid corrections for %d photcodes\n", NGridReal); +} + +// reset the values in the arrays to 0 +void resetMgrid () { + + if (!GRID_ZEROPT) return; + if (!GridCorr) return; + + for (int code = 0; code < NGridCorr; code++) { + if (!GridCorr[code]) continue; + + for (int ix = 0; ix < GridCorr[code]->Nx; ix++) { + for (int iy = 0; iy < GridCorr[code]->Ny; iy++) { + GridCorr[code]-> Mgrid[ix][iy] = 0.0; + GridCorr[code]->dMgrid[ix][iy] = 0.0; + GridCorr[code]->nMgrid[ix][iy] = 0; + } + } + } +} + +void setMgrid (Catalog *catalog, int Ncatalog) { + + // check if we are actually doing this step + if (!GRID_ZEROPT) return; + if (GRID_ZPT_MODE == GRID_ZPT_MODE_NONE) return; + + resetMgrid(); // start with values of 0 + + // loop over all measurements, accumulate Sum (in Mgrid), Sum2 (in dMgrid), and Npts + + int Nsecfilt = GetPhotcodeNsecfilt (); + + for (int nc = 0; nc < Ncatalog; nc++) { + for (int na = 0; na < catalog[nc].Naverage; na++) { + + int nm = catalog[nc].averageT[na].measureOffset; + for (int k = 0; k < catalog[nc].averageT[na].Nmeasure; k++, nm++) { + + // skip measurements marked by AREA or TIME + if (catalog[nc].measureT[nm].dbFlags & MEAS_BAD) continue; + + float Mcal = getMcal (nm, nc, MAG_CLASS_PSF); + if (isnan(Mcal)) continue; + + float Mgrp = getMgrp (nm, nc, catalog[nc].measureT[nm].airmass, NULL); + if (isnan(Mgrp)) continue; + + float Mmos = getMmos (nm, nc); + if (isnan(Mmos)) continue; + + float Mflat = getMflat (nm, nc, catalog); + if (isnan(Mflat)) continue; + + // Mrel* is the average magnitude for this star. For PS1 stacks, we have too much + // PSF variability. We need to calibrate the PSF magnitudes separately from the + // Aperture-like magnitues. (We have an option to use the kron magnitudes or the + // other apertures here). I basically need to do this analysis separately for each + // magnitude type + + float MrelPSF = getMrel (catalog, nm, nc, MAG_CLASS_PSF, MAG_SRC_CHP); + if (isnan(MrelPSF)) continue; + + float MsysPSF = PhotSysTiny (&catalog[nc].measureT[nm], &catalog[nc].averageT[na], &catalog[nc].secfilt[na*Nsecfilt], MAG_CLASS_PSF); + if (isnan(MsysPSF)) continue; + + float Moff = Mcal + Mgrp + Mmos + Mflat; + + // Msys = Mrel + Moff + Mgrid + // thus Mgrid = Msys - Mrel - Moff + + int code = catalog[nc].measureT[nm].photcode; + if (code <= 0) continue; + if (code >= NGridCorr) continue; // does not match one of our image, skip + + GridCorrectionType *grid = GridCorr[code]; + if (!grid) continue; // does not match one of our images, skip + + // edge effects could cause some positions to be slightly out of range + // probably should trap extreme outliers + int ix = MIN(MAX(0, (int)(catalog[nc].measureT[nm].Xccd * grid->dX)), grid->Nx - 1); + int iy = MIN(MAX(0, (int)(catalog[nc].measureT[nm].Yccd * grid->dY)), grid->Ny - 1); + + float dM = MsysPSF - MrelPSF - Moff; + + // XXX by the time we get here, we should have already mostly fixed up the zero + // points. reject measurements which are way off + if (fabs(dM) > 0.5) continue; + + grid-> Mgrid[ix][iy] += dM; + grid->dMgrid[ix][iy] += dM*dM; + grid->nMgrid[ix][iy] ++; + } } - return (TRUE); -} - -float getMgrid (off_t meas, int cat) { - - int i; - float value; - - if (!USE_GRID) return (0); - i = bin[cat][meas]; - if (i == -1) return (NAN); - - // during the grid annealing process, we skip over grid cells until they have enough - // valid stars to be fitted - if (gridV[i] == GRID_FROZEN) { - return (NAN); - } - - value = gridM[i]; - return (value); -} - -/* direct (non-iterative) solution for Mgrid values for all grid bins */ -void setMgridDirect (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) { - - int **gotstar, **gridmeas; - int i, j, k, Ngood, Nbad, Nmos, Ncal, Nrel, Nsys, Ngrp; - double **A, **B, *Mjx, *Wjx; - float Msys, Mcal, Mmos, Mgrp, Merr, Wsys; - double Mj, Wj; - - if (!USE_GRID) return; - - ALLOCATE (A, double *, Ngrid); - ALLOCATE (B, double *, Ngrid); - for (i = 0; i < Ngrid; i++) { - ALLOCATE (A[i], double, Ngrid); - ALLOCATE (B[i], double, 1); - memset (A[i], 0, Ngrid*sizeof(double)); - memset (B[i], 0, sizeof(double)); - } - - Ngood = Nbad = Ncal = Nmos = Nrel = Nsys = Ngrp = 0; - - ALLOCATE (gotstar, int *, Ncatalog); - for (i = 0; i < Ncatalog; i++) { - ALLOCATE (gotstar[i], int, catalog[i].Naverage); - } - - // set up gridmeas table : grid index for each measurement - ALLOCATE (gridmeas, int *, Ncatalog); - for (i = 0; i < Ncatalog; i++) { - ALLOCATE (gridmeas[i], int, catalog[i].Nmeasure); - for (j = 0; j < catalog[i].Nmeasure; j++) { - gridmeas[i][j] = -1; - } - } - for (i = 0; i < Ngrid; i++) { - for (j = 0; j < Nlist[i]; j++) { - int m, c; - m = mlist[i][j]; - c = clist[i][j]; - gridmeas[c][m] = i; - } - } - - // as we loop over the grid cells, we need to accumulate the values of Wjx for each star - ALLOCATE (Wjx, double, Ngrid); - ALLOCATE (Mjx, double, Ngrid); - - // accumulate the elements of the matrix equation. We have an equation of the form: Ax = B - // where x is the vector of grid cell values G_x (x = 0 - Ngrid), A is an Ngrid x Ngrid matrix, - // and B is an Ngrid vector. For a cell A(x,y), we need the following elements from each - // star which touches grid cell (x): - // - // Mj : sum over all measurements of Msys / dMsys^2 - // Wj : sum over all measurements of 1.0 / dMsys^2 - // Mjx : sum over all measurements which touch cell (x) of Msys / dMsys^2 - // Wjx : sum over all measurements which touch cell (x) of 1.0 / dMsys^2 - - // Mjx and Wjx can be calculated by summing over all measurements which touch the cell - // Mj requires looping over stars which touch (x) - - // this is tricky because we need to know both the measurements which touch a cell - // and the stars for which any measurement touches a cell. we need to accumulate - // sums for each star which touches as cell on both bases. - - int Nsecfilt = GetPhotcodeNsecfilt (); - int thisCode = photcodes[0][0].code; - int Nsec = GetPhotcodeNsec(thisCode); - - for (i = 0; i < Ngrid; i++) { - - for (j = 0; j < Ncatalog; j++) { - memset (gotstar[j], 0, catalog[j].Naverage*sizeof(int)); - } - - // we are looping over the stars, but doing so by looping over the set of measurements: - // every star which touches this grid cell has a measurement in Nlist[i] - for (j = 0; j < Nlist[i]; j++) { - - int mx, c, n, m0, Npts; - - mx = mlist[i][j]; - c = clist[i][j]; - n = catalog[c].measureT[mx].averef; - - // if we have already visited this star, skip the stuff below - if (gotstar[c][n]) continue; - gotstar[c][n] = TRUE; - - // skip stars marked as BAD - if (catalog[c].secfilt[n*Nsecfilt+Nsec].flags & STAR_BAD) { - Nrel ++; - continue; - } - - m0 = catalog[c].average[n].measureOffset; - - // we accumuate an entry for each cell - memset (Wjx, 0, Ngrid*sizeof(double)); - memset (Mjx, 0, Ngrid*sizeof(double)); - Npts = Mj = Wj = 0.0; - - // if we have not yet visited this star, accumulate the Mj, Wj entries - for (k = 0; k < catalog[c].average[n].Nmeasure; k++) { - - int m, Ng; - - m = m0 + k; - - // skip measurements marked as BAD - if (catalog[c].measureT[m].dbFlags & MEAS_BAD) { - Nbad ++; + } + + // now calculate Mgrid, dMgrid, nMgrid from Sum, Sum, Npt + for (int code = 0; code < NGridCorr; code++) { + if (!GridCorr[code]) continue; + + for (int ix = 0; ix < GridCorr[code]->Nx; ix++) { + for (int iy = 0; iy < GridCorr[code]->Ny; iy++) { + + // cells without sufficient coverage stay at 0.0 + if (GridCorr[code]->nMgrid[ix][iy] < 5) { + GridCorr[code]-> Mgrid[ix][iy] = 0.0; + GridCorr[code]->dMgrid[ix][iy] = 0.0; continue; } - // skip images marked as BAD - Mcal = getMcal (m, c, MAG_CLASS_PSF); - if (isnan(Mcal)) { - Ncal ++; - continue; - } - - // skip mosaics marked as BAD - Mmos = getMmos (m, c); - if (isnan(Mmos)) { - Nmos ++; - continue; - } - - // skip mosaics marked as BAD - Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass, NULL); - if (isnan(Mgrp)) { - Ngrp ++; - continue; - } - - // select the color- and airmass-corrected observed magnitude for this star - // XXX need to be able to turn off the color-correction until initial average mags are found - Msys = PhotCatTiny (&catalog[c].measureT[m], MAG_CLASS_PSF); - if (isnan(Msys)) { - Nsys++; - continue; - } - - // mag-error for this measurement - Merr = MAX (catalog[c].measureT[m].dM, MIN_ERROR); - - // disable Wsys for now - Wsys = TRUE ? 1.0 : 1.0 / SQ(Merr); - - Ng = gridmeas[c][m]; - if (Ng == -1) continue; // skip measurements which do not touch any cell - - Mj += Msys * Wsys; // we are only including measurements touching this cell - Wj += Wsys; // we are only including measurements touching this cell - Npts ++; - Ngood ++; - - Mjx[Ng] += Msys * Wsys; // we are only including measurements touching cell (x) - Wjx[Ng] += Wsys; // we are only including measurements touching cell (x) - } - - // some stars will not have any valid measurements, skip these - if (Npts == 0) continue; - - B[i][0] += Mj*Wjx[i]/Wj - Mjx[i]; - A[i][i] -= Wjx[i]; - for (k = 0; k < Ngrid; k++) { - A[i][k] += Wjx[i]*Wjx[k]/Wj; - // fprintf (stderr, "%3.0f ", Wjx[k]); + float Mgrid = GridCorr[code]-> Mgrid[ix][iy] / GridCorr[code]->nMgrid[ix][iy]; // average Mgrid + float Mgrid2 = GridCorr[code]-> dMgrid[ix][iy] / GridCorr[code]->nMgrid[ix][iy]; // average Mgrid^2 + + float r = GridCorr[code]->nMgrid[ix][iy] / (float) (GridCorr[code]->nMgrid[ix][iy] - 1.0); // pop -> sample stdev + + GridCorr[code]-> Mgrid[ix][iy] = Mgrid; + GridCorr[code]->dMgrid[ix][iy] = sqrt(r*(Mgrid2 - Mgrid*Mgrid)); // sample stdev + // fprintf (stderr, "grid code %d, %d x %d : %f +/- %f : %d\n", code, ix, iy, GridCorr[code]-> Mgrid[ix][iy], GridCorr[code]->dMgrid[ix][iy], GridCorr[code]->nMgrid[ix][iy]); } } } - - if (1) { - - FILE *f; - Header theader; - Matrix matrix; - - /* we are writing to this file */ - f = fopen ("matrix.fits", "w"); - if (f == (FILE *) NULL) { - fprintf (stderr, "cannot open matrix.fits for output\n"); - return; - } - - /* save grid mag values */ - gfits_init_header (&theader); - theader.Naxes = 2; - theader.Naxis[0] = Ngrid; - theader.Naxis[1] = Ngrid; - theader.bitpix = -32; - gfits_create_Theader (&theader, "IMAGE"); - gfits_modify (&theader, "EXTNAME", "%s", 1, "MATRIX"); - gfits_create_matrix (&theader, &matrix); - for (i = 0; i < Ngrid; i++) { - for (j = 0; j < Ngrid; j++) { - gfits_set_matrix_value (&matrix, i, j, (double) A[i][j]); - } - } - gfits_fwrite_header (f, &theader); - gfits_fwrite_matrix (f, &matrix); - gfits_free_matrix (&matrix); - - gfits_modify (&theader, "EXTNAME", "%s", 1, "TRPOSE"); - gfits_create_matrix (&theader, &matrix); - for (i = 0; i < Ngrid; i++) { - for (j = 0; j < Ngrid; j++) { - gfits_set_matrix_value (&matrix, i, j, (double) A[j][i]); - } - } - gfits_fwrite_header (f, &theader); - gfits_fwrite_matrix (f, &matrix); - gfits_free_matrix (&matrix); - fclose (f); - } - - dgaussjordan (A, B, Ngrid, 1); - - fprintf (stderr, "grid cells fitted (Ngood: %d, Nbad: %d, Nmos: %d, Ncal: %d, Nrel: %d, Nsys: %d)\n", Ngood, Nbad, Nmos, Ncal, Nrel, Nsys); - - for (i = 0; i < Ngrid; i++) { - gridM[i] = B[i][0]; - gridS[i] = sqrt(A[i][i]); - gridN[i] = Ngood; - } - - free (Wjx); - free (Mjx); - - for (i = 0; i < Ngrid; i++) { - free (A[i]); - } - free (A); - free (B); - - for (i = 0; i < Ncatalog; i++) { - free (gotstar[i]); - free (gridmeas[i]); - } - free (gotstar); - free (gridmeas); -} - -/* determine Mgrid values for all grid bins */ -void setMgrid (Catalog *catalog, FlatCorrectionTable *flatcorr) { - - int i, j, m, c, n, N, Nmax, Nbad, Nmos, Ngrp, Ncal, Nrel, Nsys, Nfit; - double *list, *dlist; - float Msys, Mrel, Mcal, Mmos, Mgrp; - - StatType stats; - liststats_setmode (&stats, "INNER_WTMEAN"); + return; +} + +float getMgrid (Measure *measure) { + + if (!GRID_ZEROPT) return 0.0; + if (GRID_ZPT_MODE == GRID_ZPT_MODE_NONE) return 0.0; + + int code = measure->photcode; + if (code <= 0) return 0.0; + if (code >= NGridCorr) return 0.0; // does not match one of our image, skip + + GridCorrectionType *grid = GridCorr[code]; + if (!grid) return 0.0; // does not match one of our images, skip - if (!USE_GRID) return; - - int Nsecfilt = GetPhotcodeNsecfilt (); - - Nmax = Nlist[0]; - for (i = 0; i < Ngrid; i++) { - Nmax = MAX (Nmax, Nlist[i]); - } - ALLOCATE (list, double, Nmax); - ALLOCATE (dlist, double, Nmax); - - Nbad = Ncal = Nmos = Ngrp = Nrel = Nsys = Nfit = 0; - - for (i = 0; i < Ngrid; i++) { - - N = 0; - for (j = 0; j < Nlist[i]; j++) { - - m = mlist[i][j]; - c = clist[i][j]; - - if (catalog[c].measureT[m].dbFlags & MEAS_BAD) { - Nbad ++; - continue; - } - Mcal = getMcal (m, c, MAG_CLASS_PSF); - if (isnan(Mcal)) { - Ncal ++; - continue; - } - Mmos = getMmos (m, c); - if (isnan(Mmos)) { - Nmos ++; - continue; - } - Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass, NULL); - if (isnan(Mgrp)) { - Ngrp ++; - continue; - } - Mrel = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); - if (isnan(Mrel)) { - Nrel ++; - continue; - } - - n = catalog[c].measureT[m].averef; - Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF); - if (isnan(Msys)) { - Nsys++; - continue; - } - list[N] = Msys - Mrel - Mcal - Mmos - Mgrp; - dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); - N++; - } - - // the reference Cell is forced to have a value of 0.0, and is never changed to GRID_FITTED - if (gridV[i] == GRID_REFERENCE) { - gridM[i] = 0.0; - gridS[i] = 0.0; - gridN[i] = N; - continue; - } - - // until we have enough valid measurements on this grid cell, skip it - if (N < GRID_TOOFEW) { - gridV[i] = GRID_FROZEN; - continue; - } - - liststats (list, dlist, NULL, N, &stats); - gridM[i] = stats.mean; - gridS[i] = stats.sigma; - gridN[i] = N; - gridV[i] = GRID_FITTED; - Nfit++; - } - - fprintf (stderr, "%d of %d grid cells fitted (+ reference cell) (Nbad: %d, Nmos: %d, Ncal: %d, Nrel: %d, Nsys: %d)\n", Nfit, Ngrid, Nbad, Nmos, Ncal, Nrel, Nsys); - - free (list); - free (dlist); -} - -void plot_grid (Catalog *catalog, FlatCorrectionTable *flatcorr) { - - int i, j, m, c, n, N, Narea; - float Msys, Mrel, Mcal, Mmos, Mgrp; - double *xlist, *Mlist, *dlist, *ylist; - Graphdata graphdata; - - if (!USE_GRID) return; - - int Nsecfilt = GetPhotcodeNsecfilt (); - - N = 0; - for (i = 0; i < Ngrid; i++) - N += Nlist[i]; - - ALLOCATE (xlist, double, N); - ALLOCATE (ylist, double, N); - ALLOCATE (Mlist, double, N); - ALLOCATE (dlist, double, N); - - Narea = 0; - N = 0; - for (i = 0; i < Ngrid; i++) { - for (j = 0; j < Nlist[i]; j++) { - - m = mlist[i][j]; - c = clist[i][j]; - - if (catalog[c].measureT[m].dbFlags & MEAS_BAD) { - Narea ++; - continue; - } - Mcal = getMcal (m, c, MAG_CLASS_PSF); - if (isnan(Mcal)) continue; - Mmos = getMmos (m, c); - if (isnan(Mmos)) continue; - Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass, NULL); - if (isnan(Mgrp)) continue; - Mrel = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); - if (isnan(Mrel)) continue; - - n = catalog[c].measureT[m].averef; - Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF); - - xlist[N] = Xmeas[c][m]; - ylist[N] = Ymeas[c][m]; - Mlist[N] = Msys - Mrel - Mcal - Mmos - Mgrp; - dlist[N] = Msys - Mrel - Mcal - Mmos - Mgrp - gridM[i]; - N++; - } - } - - fprintf (stderr, "skipped %d meas for area\n", Narea); - - plot_defaults (&graphdata); - graphdata.ymin = PlotdMmin; - graphdata.ymax = PlotdMmax; - plot_list (&graphdata, xlist, Mlist, N, "X vs dM raw", "%s.XdM.png", OUTROOT); - plot_list (&graphdata, xlist, dlist, N, "X vs dM corrected", "%s.XdMf.png", OUTROOT); - plot_list (&graphdata, ylist, dlist, N, "Y vs dM corrected", "%s.YdMf.png", OUTROOT); - - plot_defaults (&graphdata); - plot_list (&graphdata, xlist, ylist, N, "X vs Y", "%s.XY.png", OUTROOT); - - free (ylist); - free (xlist); - free (Mlist); - free (dlist); - -} - -void dump_grid () { - - off_t i, Nimage; - int j, Nbytes, Nformat; - FILE *f; - Header header, theader; - Matrix matrix; - Mosaic *refmosaic; - char *filename; - char formatline[32], key[32], value[64]; - - Nbytes = strlen (OUTROOT) + 6; - ALLOCATE (filename, char, Nbytes); - snprintf (filename, Nbytes, "%s.fits", OUTROOT); - - /* select reference mosaic image */ - // off_t *imlist = SelectRefMosaic (&refmosaic, &Nimage); return value ignored - SelectRefMosaic (&refmosaic, &Nimage); - - /* we are writing to this file */ - f = fopen (filename, "w"); - if (f == (FILE *) NULL) { - fprintf (stderr, "cannot open %s for output\n", filename); - free (filename); - return; - } - - /* create empty phu */ - gfits_init_header (&header); - header.extend = TRUE; - gfits_create_header (&header); - gfits_create_matrix (&header, &matrix); - gfits_modify (&header, "NEXTEND", OFF_T_FMT, 1, Nimage + 3); - gfits_modify (&header, "FILTER", "%s", 1, photcodes[0][0].name); // XXXX note that this expects a single photcode, enforced in initialize.d - gfits_modify_alt (&header, "COMMENT", "%S", 1, "Mosaic Photometry Grid Analysis"); - - // we need to add lines to the PHU to identify the camera and format; these are used by the ipp config system - // Note that config must have been loaded (and not freed) above. - ScanConfig (config, "NFORMAT", "%d", 1, &Nformat); - for (i = 1; i <= Nformat; i++) { - ScanConfig (config, "FORMAT", "%s", i, formatline); - sscanf (formatline, "%s %s", key, value); - gfits_modify (&header, key, "%s", 1, value); - } - - gfits_fwrite_header (f, &header); - gfits_fwrite_matrix (f, &matrix); - gfits_free_matrix (&matrix); - - /* save grid mag values */ - gfits_init_header (&theader); - theader.Naxes = 2; - theader.Naxis[0] = gridX; - theader.Naxis[1] = gridY; - theader.bitpix = -32; - gfits_create_Theader (&theader, "IMAGE"); - gfits_modify (&theader, "FILTER", "%s", 1, photcodes[0][0].name); - gfits_modify (&theader, "EXTNAME", "%s", 1, "MAG_OFFSET"); - gfits_create_matrix (&theader, &matrix); - for (i = 0; i < gridX; i++) { - for (j = 0; j < gridY; j++) { - gfits_set_matrix_value (&matrix, i, j, (double) gridM[i + j*gridX]); - } - } - write_coords (&theader, &refmosaic[0].coords); - gfits_fwrite_header (f, &theader); - gfits_fwrite_matrix (f, &matrix); - gfits_free_matrix (&matrix); - - /* save grid Nmeas values */ - gfits_modify (&theader, "EXTNAME", "%s", 1, "NMEAS"); - gfits_modify (&theader, "FILTER", "%s", 1, photcodes[0][0].name); - gfits_create_matrix (&theader, &matrix); - for (i = 0; i < gridX; i++) { - for (j = 0; j < gridY; j++) { - gfits_set_matrix_value (&matrix, i, j, (double) gridN[i + j*gridX]); - } - } - write_coords (&theader, &refmosaic[0].coords); - gfits_fwrite_header (f, &theader); - gfits_fwrite_matrix (f, &matrix); - gfits_free_matrix (&matrix); - - /* save grid sigma values */ - gfits_modify (&theader, "EXTNAME", "%s", 1, "SIGMA"); - gfits_modify (&theader, "FILTER", "%s", 1, photcodes[0][0].name); - gfits_create_matrix (&theader, &matrix); - for (i = 0; i < gridX; i++) { - for (j = 0; j < gridY; j++) { - gfits_set_matrix_value (&matrix, i, j, (double) gridS[i + j*gridX]); - } - } - write_coords (&theader, &refmosaic[0].coords); - gfits_fwrite_header (f, &theader); - gfits_fwrite_matrix (f, &matrix); - gfits_free_matrix (&matrix); - - /* calculate value for each CCD pixel, write out CCD images */ - /* grid pixels are tied to detector pixels, but are flipped to match focal plane */ - for (i = 0; i < camera.Nchip; i++) { - int N, ix, iy, x, y, X, Y, bin; - - N = ccdnum[i]; - - gfits_modify (&theader, "EXTNAME", "%s", 1, camera.ccdname[N]); - gfits_modify (&theader, "FILTER", "%s", 1, photcodes[0][0].name); - gfits_modify (&theader, "NX", "%d", 1, camera.Nx); - gfits_modify (&theader, "NY", "%d", 1, camera.Ny); - - theader.Naxis[0] = RELPHOT_GRID_X; - theader.Naxis[1] = RELPHOT_GRID_Y; - gfits_modify (&theader, "NAXIS1", "%d", 1, RELPHOT_GRID_X); - gfits_modify (&theader, "NAXIS2", "%d", 1, RELPHOT_GRID_Y); - gfits_create_matrix (&theader, &matrix); - - for (Y = 0; Y < RELPHOT_GRID_Y; Y++) { - for (X = 0; X < RELPHOT_GRID_X; X++) { - - /* normalize X & Y */ - x = X; - if (camera.Fx[N]) x = RELPHOT_GRID_X - X - 1; - y = Y; - if (camera.Fy[N]) y = RELPHOT_GRID_Y - Y - 1; - - /* coordinates in the grid */ - ix = x + camera.Ox[N]*RELPHOT_GRID_X; - iy = y + camera.Oy[N]*RELPHOT_GRID_Y; - - bin = ix + iy*gridX; - gfits_set_matrix_value (&matrix, X, Y, (double) gridM[bin]); - } - } - gfits_fwrite_header (f, &theader); - gfits_fwrite_matrix (f, &matrix); - gfits_free_matrix (&matrix); - } - - free (filename); -} - -void InterpolateGrid (float *buffer, int Nx, int Ny, Coords *ccd, Coords *gcoords) { - - int i, j; - double x, y, r, d, X, Y, dx, dy; - double V00, V01, V10, V11; - double wV00, wV01, wV10, wV11; - double dV00, dV01, dV10, dV11; - double v1, v2, value; - int ix, iy, N; - - for (i = 0; i < Nx; i++) { - for (j = 0; j < Ny; j++) { - x = i * RELPHOT_GRID_BINNING / 2; - y = j * RELPHOT_GRID_BINNING / 2; - XY_to_RD (&r, &d, x, y, ccd); - RD_to_XY (&X, &Y, r, d, gcoords); - - X = X / RELPHOT_GRID_BINNING; - Y = Y / RELPHOT_GRID_BINNING; - - ix = (int) X; - dx = X - ix; - iy = (int) Y; - dy = Y - iy; - - if (ix < 0) continue; - if (iy < 0) continue; - if (ix >= gridX) continue; - if (iy >= gridY) continue; - - N = ix + iy*gridX; - V00 = gridM[N]; - V10 = gridM[N + 1]; - V01 = gridM[N + gridX]; - V11 = gridM[N + gridX + 1]; - - dV00 = gridS[N]; - dV10 = gridS[N + 1]; - dV01 = gridS[N + gridX]; - dV11 = gridS[N + gridX + 1]; - - wV00 = (dV00 == 0) ? 0.0 : 1 / SQ(dV00); - wV01 = (dV01 == 0) ? 0.0 : 1 / SQ(dV01); - wV10 = (dV10 == 0) ? 0.0 : 1 / SQ(dV10); - wV11 = (dV11 == 0) ? 0.0 : 1 / SQ(dV11); - - v1 = wV00*V00*(1 + dx*dy - dx - dy) + - wV10*V10*(dx - dx*dy) + - wV01*V01*(dy - dx*dy) + - wV11*V11*(dx*dy); - - v2 = wV00*(1 + dx*dy - dx - dy) + - wV10*(dx - dx*dy) + - wV01*(dy - dx*dy) + - wV11*(dx*dy); - - value = v1 / v2; - buffer[j*Nx + i] = value; - } - } -} + // edge effects could cause some positions to be slightly out of range + // probably should trap extreme outliers + int ix = MIN(MAX(0, (int)(measure->Xccd * grid->dX)), grid->Nx - 1); + int iy = MIN(MAX(0, (int)(measure->Yccd * grid->dY)), grid->Ny - 1); + + float Mgrid = grid-> Mgrid[ix][iy]; + return Mgrid; +} + +float getMgridTiny (MeasureTiny *measure) { + + if (!GRID_ZEROPT) return 0.0; + if (GRID_ZPT_MODE == GRID_ZPT_MODE_NONE) return 0.0; + + int code = measure->photcode; + if (code <= 0) return 0.0; + if (code >= NGridCorr) return 0.0; // does not match one of our image, skip + + GridCorrectionType *grid = GridCorr[code]; + if (!grid) return 0.0; // does not match one of our images, skip + + // edge effects could cause some positions to be slightly out of range + // probably should trap extreme outliers + int ix = MIN(MAX(0, (int)(measure->Xccd * grid->dX)), grid->Nx - 1); + int iy = MIN(MAX(0, (int)(measure->Yccd * grid->dY)), grid->Ny - 1); + + float Mgrid = grid-> Mgrid[ix][iy]; + return Mgrid; +} + +// for historical reasons, Mflat and Mgrid have opposite signs +void setMflatFromGrid (Catalog *catalog) { + if (!GRID_ZEROPT) return; + for (off_t j = 0; j < catalog->Nmeasure; j++) { + float Mgrid = getMgrid (&catalog->measure[j]); + catalog->measure[j].Mflat += Mgrid; + } +} + Index: trunk/Ohana/src/relphot/src/ImageOps.c =================================================================== --- trunk/Ohana/src/relphot/src/ImageOps.c (revision 41601) +++ trunk/Ohana/src/relphot/src/ImageOps.c (revision 41647) @@ -319,54 +319,7 @@ } -int findCCD (off_t idx, off_t meas, int cat, MeasureTiny *measure) { - - int ccdnum; - double X, Y; - char *pname, *filter, *p, base[256]; - - /* identify the ccd on the basis of the photcode name */ - pname = GetPhotcodeNamebyCode (image[idx].photcode); - - // skip measurements which do not match one of the requested photcodes ( - // (do we not already exclude in bcatalog -- maybe needed for reload_objects? - int Ns = GetActivePhotcodeIndex (measure[0].photcode); - if (Ns < 0) return FALSE; - - filter = photcodes[Ns][0].name; - sprintf (base, "%s.%s.", MOSAICNAME, filter); - if (strncmp (pname, base, strlen (base))) return (FALSE); - p = pname + strlen(base); - - /* p is pointing at the DETECTOR ID part of the photcode: CAMERA.FILTER.DETECTOR - for now, we will let this be of the form SSSSnn where SSS is an arbitrary string - and nn is an integer -> ccdnum */ - - /// XXX we have an inconsistency here wrt dvo.layout and dvo.photcode - - while (!isdigit(*p) && *p) p++; - if (*p == 0) return (FALSE); - ccdnum = atoi (p); - - /* ccdnum is an integer, but not necessarily a sequence number. when the camera layout is - constructed, there will be null values for undefined ccdnums */ - - // old code to add this measurement to the grid cell for this chip - // ave = measureT[0].averef; - // ra = catalog[cat].averageT[ave].R - measureT[0].dR / 3600.0; - // dec = catalog[cat].averageT[ave].D - measureT[0].dD / 3600.0; - // RD_to_XY (&X, &Y, ra, dec, &image[i].coords); - - // XXX we can now use these values (but need to be careful about old formats) - X = measure[0].Xccd; - Y = measure[0].Yccd; - setGridMeasure (meas, cat, X, Y, ccdnum); - - return (TRUE); -} - void matchImage (Catalog *catalog, off_t meas, int cat, int doImageList) { off_t idx, ID; - int status; MeasureTiny *measure; @@ -381,12 +334,4 @@ catalog[cat].measureT[meas].myDet = TRUE; - if (USE_GRID) { - status = findCCD (idx, meas, cat, measure); - if (!status) { - if (VERBOSE2) fprintf (stderr, "failed to determine CCD for "OFF_T_FMT", %d\n", meas, cat); - return; - } - } - // index for (catalog, measure) -> image MeasureToImage[cat][meas] = idx; @@ -440,5 +385,7 @@ if (i == -1) return NAN; - if (image[i].flags & IMAGE_BAD) return NAN; + // allow context dependence?: in relphot_images, we want to exclude detections from poor images + // but in the final pass, we want to allow even bad detections + if (image[i].flags & (ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW)) return NAN; switch (class) { @@ -464,5 +411,4 @@ if (i == -1) return (1000); - if (image[i].flags & IMAGE_BAD) return (1000); distance = image[i].ubercalDist; // was dummy3 in structure return (distance); @@ -552,24 +498,4 @@ } -float getMflat (off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog) { - OHANA_UNUSED_PARAM(flatcorr); - - off_t i = MeasureToImage[cat][meas]; - if (i == -1) return (NAN); - - float offset = catalog[cat].measureT[meas].Mflat; - if (!isfinite(offset)) { - offset = 0.0; - } - - // to do this, I need to pass in the catalog and flatcorr pointers - // int flat_id = image[i].photom_map_id; - // if (flat_id) { - // offset = FlatCorrectionOffset (flatcorr, flat_id, catalog[cat].measureT[meas].Xccd, catalog[cat].measureT[meas].Yccd); - // } - - return (offset); -} - Coords *getCoords (off_t meas, int cat) { @@ -582,5 +508,5 @@ /* determine Mcal values for all images */ -void setMcal (Catalog *catalog, FlatCorrectionTable *flatcorr) { +void setMcal (Catalog *catalog) { off_t i, j, m, c, n; @@ -690,5 +616,5 @@ continue; } - float Mgrid = getMgrid (m, c); + float Mgrid = getMgridTiny (&catalog[c].measureT[m]); if (isnan(Mgrid)) { Ngrid++; @@ -711,7 +637,8 @@ // apply that offset as well here for this image (in other words, each detection is // being compared to the model, excluding the zero point, Mcal. The model includes - // the flat-correction. NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat) - - float Mflat = getMflat (m, c, flatcorr, catalog); + // the flat-correction. NOTE the sign of Mflat (Image.Mcal = Measure.Mcal + Mflat) + // this was inconsistent w.r.t. PhotRel pre-r41606 + + float Mflat = getMflat (m, c, catalog); n = catalog[c].measureT[m].averef; @@ -722,5 +649,5 @@ } - float Moff = Mmos + Mgrp + Mgrid - Mflat; + float Moff = Mmos + Mgrp + Mgrid + Mflat; PhotCode *code = GetPhotcodebyCode (catalog[c].measureT[m].photcode); @@ -763,6 +690,9 @@ if (mark) { image[i].flags |= ID_IMAGE_PHOTOM_FEW; + image[i].McalPSF = 0.0; + image[i].McalAPER = 0.0; + image[i].McalChiSq = NAN; image[i].dMcal = NAN; - image[i].McalChiSq = NAN; + image[i].dMagSys = NAN; image[i].nFitPhotom = 0; Nfew ++; @@ -835,5 +765,5 @@ /* determine McalTEST values for all images -- this is not used to set measure.Mcal, but only to test */ -void setMcalTest (Catalog *catalog, FlatCorrectionTable *flatcorr) { +void setMcalTest (Catalog *catalog) { off_t i, j, m, c, n; @@ -889,5 +819,5 @@ continue; } - float Mgrid = getMgrid (m, c); + float Mgrid = getMgridTiny (&catalog[c].measureT[m]); if (isnan(Mgrid)) { Ngrid++; @@ -905,7 +835,7 @@ // apply that offset as well here for this image (in other words, each detection is // being compared to the model, excluding the zero point, Mcal. The model includes - // the flat-correction. NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat) - - float Mflat = getMflat (m, c, flatcorr, catalog); + // the flat-correction. NOTE the sign of Mflat (Image.Mcal = Measure.Mcal + Mflat) + + float Mflat = getMflat (m, c, catalog); // get the PSF magnitude for thie measurement, with airmass slope applied @@ -917,5 +847,5 @@ } - float Moff = Mmos + Mgrp + Mgrid - Mflat; + float Moff = Mmos + Mgrp + Mgrid + Mflat; psfStars.alldata-> yVector[Nref] = MsysPSF - MrelPSF - Moff; @@ -968,5 +898,5 @@ // measure stats for Mcal and dMcal for (i = N = 0; i < Nimage; i++) { - if (image[i].flags & IMAGE_BAD) continue; + if (image[i].flags & (ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW)) continue; if (FREEZE_IMAGES && isMosaicChip(image[i].photcode)) continue; @@ -996,6 +926,15 @@ mark = FALSE; - image[i].flags &= ~ID_IMAGE_PHOTOM_POOR; - mark = (image[i].dMcal > MaxScatter) || (fabs(image[i].McalPSF - MedOffset) > MaxOffset); + + if (image[i].flags & ID_IMAGE_PHOTOM_FEW) { + mark = TRUE; + } + if (image[i].dMcal > MaxScatter) { + mark = TRUE; + } + if (fabs(image[i].McalPSF - MedOffset) > MaxOffset) { + mark = TRUE; + } + if (mark) { Nmark ++; @@ -1146,5 +1085,5 @@ } -StatType statsImageN (Catalog *catalog, FlatCorrectionTable *flatcorr) { +StatType statsImageN (Catalog *catalog) { off_t i, j, m, c, n, N; @@ -1163,5 +1102,5 @@ n = 0; for (i = 0; i < Nimage; i++) { - if (image[i].flags & IMAGE_BAD) continue; + if (image[i].flags & (ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW)) continue; if (FREEZE_IMAGES && isMosaicChip(image[i].photcode)) continue; @@ -1179,5 +1118,5 @@ Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass, NULL); if (isnan(Mgrp)) continue; - Mgrid = getMgrid (m, c); + Mgrid = getMgridTiny (&catalog[c].measureT[m]); if (isnan(Mgrid)) continue; N++; @@ -1212,5 +1151,5 @@ for (i = 0; i < Nimage; i++) { - if (image[i].flags & IMAGE_BAD) continue; + if (image[i].flags & (ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW)) continue; if (FREEZE_IMAGES && isMosaicChip(image[i].photcode)) continue; @@ -1245,5 +1184,5 @@ for (i = 0; i < Nimage; i++) { - if (image[i].flags & IMAGE_BAD) continue; + if (image[i].flags & (ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW)) continue; if (FREEZE_IMAGES && isMosaicChip(image[i].photcode)) continue; @@ -1278,5 +1217,5 @@ for (i = 0; i < Nimage; i++) { - if (image[i].flags & IMAGE_BAD) continue; + if (image[i].flags & (ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW)) continue; if (FREEZE_IMAGES && isMosaicChip(image[i].photcode)) continue; Index: trunk/Ohana/src/relphot/src/MagResidSave.c =================================================================== --- trunk/Ohana/src/relphot/src/MagResidSave.c (revision 41601) +++ trunk/Ohana/src/relphot/src/MagResidSave.c (revision 41647) @@ -192,4 +192,6 @@ CHECK_STATUS (!status, "ERROR: problem closing meanmags file %s\n", filename); + FitDataSetFree (&psfStars); + return TRUE; } @@ -213,5 +215,4 @@ off_t c = ImageToCatalog[j]; - // XXX drop this? if (catalog[c].measureT[m].dbFlags & MEAS_BAD) continue; @@ -225,5 +226,5 @@ if (isnan(Mgrp)) continue; - float Mgrid = getMgrid (m, c); + float Mgrid = getMgridTiny (&catalog[c].measureT[m]); if (isnan(Mgrid)) continue; @@ -231,7 +232,7 @@ // apply that offset as well here for this image (in other words, each detection is // being compared to the model, excluding the zero point, Mcal. The model includes - // the flat-correction. NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat) - - // float Mflat = getMflat (m, c, flatcorr, catalog); + // the flat-correction. NOTE the sign of Mflat (Image.Mcal = Measure.Mcal + Mflat) + + float Mflat = getMflat (m, c, catalog); // MrelPSF is the average magnitude for this star. @@ -246,8 +247,5 @@ if (isnan(MsysPSF)) continue; - // XXX skip Mflat for now - // float Moff = Mcal + Mmos + Mgrp + Mgrid - Mflat; - - float Moff = Mcal + Mmos + Mgrp + Mgrid; + float Moff = Mcal + Mmos + Mgrp + Mgrid + Mflat; // Mrel is the true average apparent magnitude of this star Index: trunk/Ohana/src/relphot/src/MosaicOps.c =================================================================== --- trunk/Ohana/src/relphot/src/MosaicOps.c (revision 41601) +++ trunk/Ohana/src/relphot/src/MosaicOps.c (revision 41647) @@ -1,5 +1,5 @@ # include "relphot.h" off_t findMosaic (unsigned int *startTimes, off_t Nmosaic, unsigned int start); -int save_test_mosaic_measures (FILE *fout, int Nmos, Catalog *catalog, FlatCorrectionTable *flatcorr); +int save_test_mosaic_measures (FILE *fout, int Nmos, Catalog *catalog); // see discussion in ImagesOps.c re: IDX_T @@ -50,104 +50,4 @@ # undef COMPARE -} - -/* find mosaic frames (unique time periods & photcode name matches mosaic) */ -void initMosaics_old (Image *image, off_t Nimage) { - - off_t i, j, status, found, NMOSAIC, *MosaicN_IMAGE; - unsigned int start, stop; - char *pname; - - if (!MOSAIC_ZEROPT) return; - - Nmosaic = 0; - NMOSAIC = 10; - ALLOCATE (mosaic, Mosaic, NMOSAIC); - - ALLOCATE (MosaicToImage, off_t *, NMOSAIC); - ALLOCATE (MosaicN_Image, off_t, NMOSAIC); - ALLOCATE (MosaicN_IMAGE, off_t, NMOSAIC); - - ALLOCATE (ImageToMosaic, off_t, Nimage); // mosaic to which image belongs - - /* a 'mosaic' in relphot is (unlike relastro) a virtual concept: there is no - * entry in the image table that represents this mosaic. Instead, it is an - * internal construct that defines a group of related images - */ - - /* generate list of unique mosaics */ - for (i = 0; i < Nimage; i++) { - ImageToMosaic[i] = -1; - - /* select valid mosaic images by photcode */ - pname = GetPhotcodeNamebyCode (image[i].photcode); - status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME)); - if (status) continue; - - /* set image time range */ - start = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1); - stop = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1); - - /* find existing mosaic with this time range */ - found = FALSE; - for (j = 0; !found && (j < Nmosaic); j++) { - if (stop < mosaic[j].start) continue; - if (start > mosaic[j].stop) continue; - found = TRUE; - - // add reference from image to mosaic - ImageToMosaic[i] = j; - - /* add image to mosaic image list */ - MosaicToImage[j][MosaicN_Image[j]] = i; - MosaicN_Image[j] ++; - if (MosaicN_Image[j] == MosaicN_IMAGE[j]) { - MosaicN_IMAGE[j] += 10; - REALLOCATE (MosaicToImage[j], off_t, MosaicN_IMAGE[j]); - } - - } - if (found) continue; - - /* a new mosaic, define ranges -- preserve the original values incase this image is not used */ - mosaic[Nmosaic].start = start; - mosaic[Nmosaic].stop = stop; - mosaic[Nmosaic].McalPSF = 0.0; // note : at the end, mosaic.Mcal is added back to the input images - mosaic[Nmosaic].McalAPER = 0.0; // note : mosaic stores only offsets relative to the original image values - mosaic[Nmosaic].dMcal = 0.0; // note : at the end, mosaic.Mcal is added back to the input images - mosaic[Nmosaic].dMsys = 0.0; - mosaic[Nmosaic].McalChiSq = 0.0;// NAN or 0.0? - mosaic[Nmosaic].flags = image[i].flags; - mosaic[Nmosaic].secz = image[i].secz; - mosaic[Nmosaic].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode); - - // XXX do we need to do something about flag consistency across a mosaic? - - /* add image to mosaic image list */ - MosaicN_IMAGE[Nmosaic] = 10; - MosaicN_Image[Nmosaic] = 1; - ALLOCATE (MosaicToImage[Nmosaic], off_t, MosaicN_IMAGE[Nmosaic]); - MosaicToImage[Nmosaic][0] = i; - - // add reference from image to mosaic - ImageToMosaic[i] = Nmosaic; - - Nmosaic ++; - if (Nmosaic == NMOSAIC) { - NMOSAIC += 10; - REALLOCATE (mosaic, Mosaic, NMOSAIC); - REALLOCATE (MosaicToImage, off_t *, NMOSAIC); - REALLOCATE (MosaicN_Image, off_t, NMOSAIC); - REALLOCATE (MosaicN_IMAGE, off_t, NMOSAIC); - } - } - - // free this or not? - free (MosaicN_IMAGE); - - initMosaicGrid (image, Nimage); - - fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nimage, (int) Nmosaic); - return; } @@ -318,5 +218,5 @@ free (startTimesMosaic); - initMosaicGrid (subset, Nsubset); + initMosaicMcal (subset, Nsubset); fprintf (stderr, "matched %d images to %d mosaics, %d simple chips not matched to mosaics\n", (int) (Nsubset - Nsimple), (int) Nmosaic, (int) Nsimple); @@ -460,5 +360,5 @@ if (mergeMcal) { - initMosaicGrid (image, Nimage); + initMosaicMcal (image, Nimage); } @@ -511,4 +411,7 @@ } +// this function sets mosaic->coords to the median of the individual chips. This +// coordinate frame is used by the parallel region analysis to assign exposures (mosaics) +// to specific machines by a single center (rather than individual chips) void setMosaicCenters (Image *image, off_t Nimage) { OHANA_UNUSED_PARAM(Nimage); @@ -599,95 +502,24 @@ mosaic[i].coords.cdelt1 = 1.0 / 3600.0; mosaic[i].coords.cdelt2 = 1.0 / 3600.0; - - mosaic[i].McalPSF = 0.0; - mosaic[i].McalAPER = 0.0; - mosaic[i].dMcal = 0.0; - mosaic[i].dMsys = 0.0; - mosaic[i].McalChiSq = 0.0; } return; } -void initMosaicGrid (Image *image, off_t Nimage) { +void initMosaicMcal (Image *image, off_t Nimage) { OHANA_UNUSED_PARAM(Nimage); - /* find max dR, dD range for all mosaics */ - /* define mosaic.coords to cover dR, dD */ - /* send results to initGridBins */ - - off_t i, j, m, NX, NY; - int dXmax, dYmax; - double dS, dX, dY; - double R, D, Rmin, Rmax, Dmin, Dmax; double McalPSF, McalAPER, dMcal, McalChiSq; fprintf (stderr, "*** moving Mcal from image.Mcal to mosaic.Mcal ***\n"); - dXmax = dYmax = 0.0; - for (i = 0; i < Nmosaic; i++) { - Dmin = Rmin = 360.0; - Dmax = Rmax = -360.0; - dS = 0.0; + for (off_t i = 0; i < Nmosaic; i++) { McalPSF = McalAPER = dMcal = McalChiSq = 0; - for (j = 0; j < MosaicN_Image[i]; j++) { - m = MosaicToImage[i][j]; - NX = image[m].NX; - NY = image[m].NY; - dS += hypot (image[m].coords.cdelt1*image[m].coords.pc1_1, image[m].coords.cdelt1*image[m].coords.pc2_1); - - OhanaProjection proj = GetProjection (image[m].coords.ctype); - if ((proj == PROJ_WRP) && !image[m].coords.mosaic) { - XY_to_LM (&R, &D, 0.0, 0.0, &image[m].coords); - } else { - XY_to_RD (&R, &D, 0.0, 0.0, &image[m].coords); - } - Rmin = MIN (Rmin, R); - Rmax = MAX (Rmax, R); - Dmin = MIN (Dmin, D); - Dmax = MAX (Dmax, D); - - if ((proj == PROJ_WRP) && !image[m].coords.mosaic) { - XY_to_LM (&R, &D, (double) NX, 0.0, &image[m].coords); - } else { - XY_to_RD (&R, &D, (double) NX, 0.0, &image[m].coords); - } - Rmin = MIN (Rmin, R); - Rmax = MAX (Rmax, R); - Dmin = MIN (Dmin, D); - Dmax = MAX (Dmax, D); - - if ((proj == PROJ_WRP) && !image[m].coords.mosaic) { - XY_to_LM (&R, &D, (double) NX, (double) NY, &image[m].coords); - } else { - XY_to_RD (&R, &D, (double) NX, (double) NY, &image[m].coords); - } - Rmin = MIN (Rmin, R); - Rmax = MAX (Rmax, R); - Dmin = MIN (Dmin, D); - Dmax = MAX (Dmax, D); - - if ((proj == PROJ_WRP) && !image[m].coords.mosaic) { - XY_to_LM (&R, &D, 0.0, (double) NY, &image[m].coords); - } else { - XY_to_RD (&R, &D, 0.0, (double) NY, &image[m].coords); - } - Rmin = MIN (Rmin, R); - Rmax = MAX (Rmax, R); - Dmin = MIN (Dmin, D); - Dmax = MAX (Dmax, D); - - // XXX : this probably does not handle mosaics at RA = 0,360 well - - /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */ - - // XXX: how does this work with UBERCAL? We want to keep the Mcal values supplied by ubercal, but - // solve for a single offset for each exposure (Mosaic.Mcal). - // we also want to keep the flat-field terms for each exposure (regardless of ubercal or not) - // if it helps, note that ubercal uses a single zp per exposure, so the mean of those values is the same as the value + for (off_t j = 0; j < MosaicN_Image[i]; j++) { + off_t m = MosaicToImage[i][j]; if (!isfinite(image[m].McalPSF)) { - image[m].McalPSF = 0.0; - image[m].McalAPER = 0.0; - image[m].dMcal = 0.0; + image[m].McalPSF = 0.0; + image[m].McalAPER = 0.0; + image[m].dMcal = 0.0; image[m].McalChiSq = 0.0; fprintf (stderr, "warning: resetting NAN value for Mcal %s\n", image[m].name); @@ -703,15 +535,5 @@ image[m].dMcal = NAN; image[m].McalChiSq = NAN; - - } - dS /= MosaicN_Image[i]; - - InitCoords (&mosaic[i].coords, "DEC--TAN"); - mosaic[i].coords.crval1 = Rmin; - mosaic[i].coords.crval2 = Dmin; - mosaic[i].coords.cdelt1 = dS; - mosaic[i].coords.cdelt2 = dS; - - RD_to_XY (&dX, &dY, Rmax, Dmax, &mosaic[i].coords); + } mosaic[i].McalPSF = McalPSF / MosaicN_Image[i]; @@ -720,9 +542,4 @@ mosaic[i].McalChiSq = McalChiSq / MosaicN_Image[i]; } - if (!USE_GRID) return; - - dXmax = MAX (dXmax, dX); - dYmax = MAX (dYmax, dY); - initGrid (dXmax, dYmax); return; } @@ -963,5 +780,7 @@ if (i == -1) return (0.0); - if (mosaic[i].flags & IMAGE_BAD) return (NAN); + // if the mosaic cannot be calibrated (too few measurements, skip it) + // XXX: if the mosaic is bad (few) but the chips are OK, I can used them. + // XXX: if (mosaic[i].flags & ID_IMAGE_PHOTOM_FEW) return (NAN); value = mosaic[i].McalPSF; return (value); @@ -1001,11 +820,10 @@ Catalog *catalog; Image *image; - FlatCorrectionTable *flatcorr; SetMmosInfo info; } ThreadInfo; -int setMmos_mosaic (Mosaic *mosaic, off_t Nmos, Image *image, Catalog *catalog, SetMmosInfo *info, FlatCorrectionTable *flatcorr); +int setMmos_mosaic (Mosaic *mosaic, off_t Nmos, Image *image, Catalog *catalog, SetMmosInfo *info); void *setMmos_worker (void *data); -int setMmos_threaded (Catalog *catalog, FlatCorrectionTable *flatcorr); +int setMmos_threaded (Catalog *catalog); void SetMmosInfoInit (SetMmosInfo *info, off_t Nmax, int allocLists) { @@ -1073,5 +891,5 @@ } -int setMmos (Catalog *catalog, FlatCorrectionTable *flatcorr) { +int setMmos (Catalog *catalog) { off_t i, N, Nmax; @@ -1086,5 +904,5 @@ // so do not run setMmos in threaded mode if PLOTSTUFF is set if (NTHREADS && !PLOTSTUFF) { - int status = setMmos_threaded (catalog, flatcorr); + int status = setMmos_threaded (catalog); return status; } @@ -1106,5 +924,5 @@ for (i = 0; i < Nmosaic; i++) { - setMmos_mosaic (&mosaic[i], i, image, catalog, &info, flatcorr); + setMmos_mosaic (&mosaic[i], i, image, catalog, &info); } SetMmosInfoFree (&info); @@ -1123,5 +941,5 @@ // 'mosaic' is a pointer to the current mosaic of interest (Nmos) -int setMmos_mosaic (Mosaic *myMosaic, off_t Nmos, Image *image, Catalog *catalog, SetMmosInfo *info, FlatCorrectionTable *flatcorr) { +int setMmos_mosaic (Mosaic *myMosaic, off_t Nmos, Image *image, Catalog *catalog, SetMmosInfo *info) { off_t j, NimageReal; @@ -1204,5 +1022,5 @@ snprintf (filename, 64, "test.%05d.%02d.dat", (int) Nmos, npass_output); fout = fopen (filename, "w"); - save_test_mosaic_measures (fout, Nmos, catalog, flatcorr); + save_test_mosaic_measures (fout, Nmos, catalog); fclose (fout); } @@ -1231,5 +1049,5 @@ continue; } - float Mgrid = getMgrid (m, c); + float Mgrid = getMgridTiny (&catalog[c].measureT[m]); if (isnan(Mgrid)) { info->Ngrid ++; @@ -1253,7 +1071,8 @@ // apply that offset as well here for this image (in other words, each detection is // being compared to the model, excluding the zero point, Mcal. The model includes - // the flat-correction. NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat) - - float Mflat = getMflat (m, c, flatcorr, catalog); + // the flat-correction. NOTE the sign of Mflat (Image.Mcal = Measure.Mcal + Mflat). + // this was inconsistent pre r41606 + + float Mflat = getMflat (m, c, catalog); off_t n = catalog[c].measureT[m].averef; @@ -1278,5 +1097,5 @@ assert (Nbright >= 0); - float Moff = Mcal + Mgrp + Mgrid - Mflat; + float Moff = Mcal + Mgrp + Mgrid + Mflat; psfStars->alldata-> yVector[N] = MsysPSF - MrelPSF - Moff; @@ -1370,5 +1189,5 @@ } -int save_test_mosaic_measures (FILE *fout, int Nmos, Catalog *catalog, FlatCorrectionTable *flatcorr) { +int save_test_mosaic_measures (FILE *fout, int Nmos, Catalog *catalog) { int Nsecfilt = GetPhotcodeNsecfilt (); @@ -1382,12 +1201,12 @@ float Mcal = getMcal (m, c, MAG_CLASS_PSF); // float Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass, NULL); - float Mgrid = getMgrid (m, c); + float Mgrid = getMgridTiny (&catalog[c].measureT[m]); float MrelPSF = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); - float Mflat = getMflat (m, c, flatcorr, catalog); + float Mflat = getMflat (m, c, catalog); off_t n = catalog[c].measureT[m].averef; float MsysPSF = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF); - float delta = MsysPSF - MrelPSF - Mcal - Mgrid + Mflat; + float delta = MsysPSF - MrelPSF - Mcal - Mgrid - Mflat; int isBad = (catalog[c].measureT[m].dbFlags & MEAS_BAD); @@ -1398,5 +1217,5 @@ } -int setMmos_threaded (Catalog *catalog, FlatCorrectionTable *flatcorr) { +int setMmos_threaded (Catalog *catalog) { int i; @@ -1435,5 +1254,4 @@ threadinfo[i].catalog = catalog; threadinfo[i].image = image; - threadinfo[i].flatcorr = flatcorr; // we do NOT allocate the arrays here, we only supply basic info (Nmax, Nloop) used in @@ -1481,4 +1299,5 @@ summary.Nskip); free (threadinfo); + // XXX SetMmosInfoFree (&summary); npass_output ++; @@ -1503,8 +1322,7 @@ Catalog *catalog = threadinfo->catalog; - FlatCorrectionTable *flatcorr = threadinfo->flatcorr; Image *image = threadinfo->image; - setMmos_mosaic (&mosaic[i], i, image, catalog, &results, flatcorr); + setMmos_mosaic (&mosaic[i], i, image, catalog, &results); SetMmosInfoAccum (&threadinfo->info, &results); } @@ -1512,205 +1330,4 @@ SetMmosInfoFree (&results); return NULL; -} - -// When we rationalize the images/mosaics, we are driving the negative cloud images back -// to 0.0. At the same time, we make a guess to the effective impact on all other images, -// driven by the coupling of common stars. -// XXX this function is no longer used because we force significantly negative clouds to 0.0 -int rationalize_mosaics (Catalog *catalog, int Ncatalog) { - - double *imageOffset, **starOffset; - int **starNcount, *seclist; - int **Slist, *NSlist, *NSLIST; - int i, j, k, m, nMos, Ns, found; - - off_t Nimage; - Image *image; - - // set a test value for now - float CLOUD_TOLERANCE = 0.025; - - if (!MOSAIC_ZEROPT) return (FALSE); - if (FREEZE_MOSAICS) return (FALSE); - - image = getimages (&Nimage, NULL); - - ALLOCATE (imageOffset, double, Nmosaic); - ALLOCATE ( Slist, int *, Nmosaic); // array of calibrated star indexes on this mosaic - ALLOCATE (NSlist, int, Nmosaic); // number of stars on mosaic - ALLOCATE (NSLIST, int, Nmosaic); // number of Slist entries allocated - memset (Slist, 0, Nmosaic*sizeof(int *)); - - // find the images / mosaics with negative clouds and save their offset - for (i = 0; i < Nmosaic; i++) { - - NSlist[i] = 0; - NSLIST[i] = 100; - ALLOCATE (Slist[i], int, NSLIST[i]); - - imageOffset[i] = 0.0; - - if (VERBOSE2 && (fabs(mosaic[i].McalPSF) < CLOUD_TOLERANCE)) { - fprintf (stderr, "cloud-free: %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].McalPSF); - } - if (VERBOSE2 && (mosaic[i].McalPSF < -CLOUD_TOLERANCE)) { - imageOffset[i] = -mosaic[i].McalPSF; - // NOTE the negative sign: down below, we are going to add in the negative of McalPSF - // to this image, and the propagated mean values for other images - fprintf (stderr, "anti-clouds: %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].McalPSF); - } - if (VERBOSE2 && (mosaic[i].McalPSF > CLOUD_TOLERANCE)) { - fprintf (stderr, "cloudy : %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].McalPSF); - } - } - - int Nsecfilt = GetPhotcodeNsecfilt (); - ALLOCATE (seclist, int, Nphotcodes); - for (Ns = 0; Ns < Nphotcodes; Ns ++) { - int thisCode = photcodes[Ns][0].code; - seclist[Ns] = GetPhotcodeNsec(thisCode); - } - - // allocate an array for star offsets - int Nstars = 0; - int NSTARS = 1000; - ALLOCATE (starOffset, double *, NSTARS); - ALLOCATE (starNcount, int *, NSTARS); - memset (starOffset, 0, NSTARS*sizeof(double *)); - memset (starNcount, 0, NSTARS*sizeof(int *)); - - // find the mean offset for each star - for (i = 0; i < Ncatalog; i++) { - for (j = 0; j < catalog[i].Naverage; j++) { - ALLOCATE (starOffset[Nstars], double, Nphotcodes); - ALLOCATE (starNcount[Nstars], int, Nphotcodes); - memset (starOffset[Nstars], 0, Nphotcodes*sizeof(double)); - memset (starNcount[Nstars], 0, Nphotcodes*sizeof(int)); - - m = catalog[i].averageT[j].measureOffset; - - // determine the mosaic for each measurement - for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) { - - // skip unused measurements - if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; - - // skip unused measurements - - // XXX replace with: - // Ns = GetActivePhotcodeIndex (photcode); - // if (Ns < 0) continue; - - int Nsec; - int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode); - found = FALSE; - for (Ns = 0; !found && (Ns < Nphotcodes); Ns ++) { - if (ecode == photcodes[Ns][0].code) { - found = TRUE; - break; - } - } - if (!found) continue; - Nsec = seclist[Ns]; - - // bad stars for this secfilt - if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; - - // skip REF measurements (not tied to an image) - if (getImageEntry (m, i) < 0) continue; - - // find the source of this measurement (skip unassigned measurements) - nMos = MeasureToMosaic[i][m]; - if (nMos == -1) continue; - - // XXX how can this not be true? - if (mosaic[nMos].photcode != ecode) { - fprintf (stderr, "*"); - } - - assert (Ns >= 0); - assert (Ns < Nphotcodes); - - // accumulate the offsets from the negative cloud images (others have 0.0 value) - starOffset[Nstars][Ns] += imageOffset[nMos]; - starNcount[Nstars][Ns] ++; - - // record the mosaic->star reference - Slist[nMos][NSlist[nMos]] = Nstars; - NSlist[nMos] ++; - if (NSlist[nMos] == NSLIST[nMos]) { - NSLIST[nMos] += 100; - REALLOCATE (Slist[nMos], int, NSLIST[nMos]); - } - } - Nstars ++; - if (Nstars == NSTARS) { - NSTARS += 1000; - REALLOCATE (starOffset, double *, NSTARS); - REALLOCATE (starNcount, int *, NSTARS); - memset (&starOffset[NSTARS-1000], 0, 1000*sizeof(double *)); - memset (&starNcount[NSTARS-1000], 0, 1000*sizeof(int *)); - } - } - } - - // find the mean offset of the images without negative clouds - for (i = 0; i < Nmosaic; i++) { - - found = FALSE; - for (Ns = 0; !found && (Ns < Nphotcodes); Ns ++) { - if (mosaic[i].photcode == photcodes[Ns][0].code) { - found = TRUE; - break; - } - } - if (!found) { - fprintf (stderr, "invalid photcode for mosaic?\n"); - abort(); - } - - // a negative cloud image (cloud: McalPSF > 0; anti-clouds: McalPSF < 0; imageOffset = -McalPSF) - if (imageOffset[i] > 0.0) continue; - - // we need to actually have cross-references to count - if (NSlist[i] < 2) continue; - - float dM = 0.0; - for (j = 0; j < NSlist[i]; j++) { - k = Slist[i][j]; - if (starNcount[k][Ns] > 1) { - dM += (starOffset[k][Ns] / starNcount[k][Ns]); - } - } - imageOffset[i] = dM / NSlist[i]; - // fprintf (stderr, "adjust image: %s : (%f %d) : %f\n", image[MosaicToImage[i][0]].name, dM, NSlist[i], imageOffset[i]); - } - - // for (i = 0; i < Nmosaic; i++) { - // fprintf (stderr, "correction: %s : %f\n", image[MosaicToImage[i][0]].name, imageOffset[i]); - // } - - // apply all offset values to the mosaics - // find the images / mosaics with negative clouds and save their offset - for (i = 0; i < Nmosaic; i++) { - mosaic[i].McalPSF += imageOffset[i]; - } - - for (i = 0; i < Nstars; i++) { - free (starOffset[i]); - free (starNcount[i]); - } - free (starOffset); - free (starNcount); - - free (seclist); - free (NSlist); - free (NSLIST); - for (i = 0; i < Nmosaic; i++){ - free (Slist[i]); - } - free (Slist); - free (imageOffset); - return (TRUE); } @@ -1731,5 +1348,5 @@ n = 0; for (i = 0; i < Nmosaic; i++) { - if (mosaic[i].flags & IMAGE_BAD) continue; + if (mosaic[i].flags & (ID_IMAGE_MOSAIC_POOR | ID_IMAGE_PHOTOM_FEW)) continue; if (mosaic[i].skipCal) continue; if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue; @@ -1764,5 +1381,5 @@ n = 0; for (i = 0; i < Nmosaic; i++) { - if (mosaic[i].flags & IMAGE_BAD) continue; + if (mosaic[i].flags & (ID_IMAGE_MOSAIC_POOR | ID_IMAGE_PHOTOM_FEW)) continue; if (mosaic[i].skipCal) continue; if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue; @@ -1781,5 +1398,5 @@ } -StatType statsMosaicN (Catalog *catalog, FlatCorrectionTable *flatcorr) { +StatType statsMosaicN (Catalog *catalog) { off_t i, j, m, c, n, N; @@ -1797,5 +1414,5 @@ n = 0; for (i = 0; i < Nmosaic; i++) { - if (mosaic[i].flags & IMAGE_BAD) continue; + if (mosaic[i].flags & (ID_IMAGE_MOSAIC_POOR | ID_IMAGE_PHOTOM_FEW)) continue; if (mosaic[i].skipCal) continue; if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue; @@ -1811,5 +1428,5 @@ Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass, NULL); if (isnan(Mgrp)) continue; - Mgrid = getMgrid (m, c); + Mgrid = getMgridTiny (&catalog[c].measureT[m]); if (isnan(Mgrid)) continue; Mrel = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); @@ -1847,5 +1464,5 @@ n = 0; for (i = 0; i < Nmosaic; i++) { - if (mosaic[i].flags & IMAGE_BAD) continue; + if (mosaic[i].flags & (ID_IMAGE_MOSAIC_POOR | ID_IMAGE_PHOTOM_FEW)) continue; if (mosaic[i].skipCal) continue; if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue; @@ -1866,6 +1483,8 @@ /* mark mosaic if: abs(Mcal - ) too large, dMcal too large */ -static float MinMaxChiSq = NAN; -static float MinMaxScatter = NAN; +static float MinChiSqLim = NAN; +static float MaxChiSqLim = NAN; +static float MinScatterLim = NAN; +static float MaxScatterLim = NAN; void clean_mosaics () { @@ -1877,4 +1496,7 @@ if (VERBOSE) fprintf (stderr, "marking poor mosaics\n"); + + if (isnan (MaxChiSqLim)) MaxChiSqLim = MOSAIC_CHISQ; + if (isnan (MaxScatterLim)) MaxScatterLim = MOSAIC_SCATTER; ALLOCATE (mlist, double, Nmosaic); @@ -1899,13 +1521,13 @@ liststats (mlist, NULL, NULL, N, &stats); float ChiSqUpper90 = stats.Upper90; - if (isnan (MinMaxChiSq)) MinMaxChiSq = 2.0*stats.median; - float MaxChiSq = MAX (MinMaxChiSq, ChiSqUpper90); + if (isnan (MinChiSqLim)) MinChiSqLim = 2.0*stats.median; // chi-square cut cannot fall below this value (even if this is > MaxChiSqLim) + float ChiSqLimit = MAX(MinChiSqLim, MIN(MaxChiSqLim, ChiSqUpper90)); // chi-square cut should be between MinChiSqLim and MaxChiSqLim liststats (slist, NULL, NULL, N, &stats); float ScatterUpper90 = stats.Upper90; - if (isnan (MinMaxScatter)) MinMaxScatter = 2.0*stats.median; - float MaxScatter = MAX (MinMaxScatter, ScatterUpper90); - - fprintf (stderr, "MOSAIC: Max ChiSq: %f, Max Scatter: %f | ChiSquare Upper 90: %f, Scatter Upper 90: %f\n", MaxChiSq, MaxScatter, ChiSqUpper90, ScatterUpper90); + if (isnan (MinScatterLim)) MinScatterLim = 2.0*stats.median; // scatter cut cannto fall below this value (even if this is > MaxScatterLim) + float ScatterLimit = MAX(MinScatterLim, MIN(MaxScatterLim, ScatterUpper90)); + + fprintf (stderr, "MOSAIC: ChiSqLimit: %f, ScatterLimit: %f | ChiSquare Upper 90: %f, Scatter Upper 90: %f\n", ChiSqLimit, ScatterLimit, ChiSqUpper90, ScatterUpper90); int Ntotal = 0, Npoor = 0, Nmark = 0, Nscatter = 0, Nfew = 0, Nchisq = 0; @@ -1930,9 +1552,9 @@ Nfew ++; } - if (mosaic[i].stdev > MaxScatter) { + if (mosaic[i].stdev > ScatterLimit) { mark = TRUE; Nscatter ++; } - if (mosaic[i].McalChiSq > MaxChiSq) { + if (mosaic[i].McalChiSq > ChiSqLimit) { mark = TRUE; Nchisq ++; @@ -1979,5 +1601,5 @@ c = MosaicToCatalog[i][j]; - if (catalog[c].measureT[m].dbFlags & (ID_MEAS_AREA | ID_MEAS_NOCAL)) continue; + if (catalog[c].measureT[m].dbFlags & MEAS_BAD) continue; // ave = catalog[c].measureT[m].averef; Index: trunk/Ohana/src/relphot/src/RepairWarpMeasures.c =================================================================== --- trunk/Ohana/src/relphot/src/RepairWarpMeasures.c (revision 41601) +++ (revision ) @@ -1,181 +1,0 @@ -# include "relphot.h" - -FILE *logfile = NULL; - -void RepairWarpMeasuresOpenLogfile () { - - char name[DVO_MAX_PATH]; - snprintf (name, DVO_MAX_PATH, "%s/repair.warp.log", HOSTDIR); - logfile = fopen (name, "w"); - if (!logfile) fprintf (stderr, "failed to open repair warp logfile %s\n", name); -} - -void RepairWarpMeasuresCloseLogfile () { - if (!logfile) return; - fclose (logfile); -} - -// XXX need to load_images first and generate some lookup tables -int RepairWarpMeasures (Catalog *catalog) { - - off_t Nimage; - Image *image = getimages (&Nimage, NULL); - - myAssert (!catalog->Nmeasure || catalog->measure, "programming error"); - - int NfixChipID = 0, NfixStackID = 0, NfixWarpID = 0, NfixWarpImageID = 0, NmissWarp = 0, NmissStack = 0, NbadWarp = 0; - - int onePercent = catalog->Nmeasure / 100; - - for (off_t j = 0; j < catalog->Nmeasure; j++) { - if (j % onePercent == 0) fprintf (stderr, "."); - - Measure *measure = &catalog->measure[j]; - - // get the associated average value: - int averef = measure->averef; - Average *average = &catalog->average[averef]; - - // double-check for consistency - myAssert (average->objID == measure->objID, "objID mismatch: %f %f : %d %d\n", average->R, average->D, average->objID, measure->objID); - - // repair extID for non-warps: - if (isGPC1chip(measure->photcode)) { - double mjd = ohana_sec_to_mjd (measure->t); - int ccdnum = measure->photcode % 100; - uint64_t extID = CreatePSPSDetectionID (mjd, ccdnum, measure->detID); - if (extID != measure->extID) { - measure->extID = extID; - NfixChipID ++; - } - continue; - } - - // repair extID for non-warps: - if (isGPC1stack(measure->photcode)) { - int im = getImageByID (measure->imageID); - if (im < 0) { - // if we don't have a matching image, we are probably assigned to an image on the wrong side of the sky - // reconstruct the imageID from the R,D,photcode & X,Y consistency - im = GetStackSeq (image, average->R, average->D, measure->photcode, measure->Xccd, measure->Yccd); - if (im < 0) { - // still failed, this is bad - if (NmissStack < 10) fprintf (stderr, "missing stack exposure: %f %f : %d %d\n", measure->R, measure->D, measure->imageID, measure->photcode); - NmissStack ++; - continue; - } - // the old imageID was wrong. replace with the new one: - if (logfile) fprintf (logfile, "fix stack %10.6f %10.6f %5d %6.1f %6.1f 0x%08x %d -> %d\n", average->R, average->D, measure->photcode, measure->Xccd, measure->Yccd, measure->detID, measure->imageID, image[im].imageID); - measure->imageID = image[im].imageID; - } else { - // we have a matching image, but check that it is correct by re-projecting the pixel position - double Xtst, Ytst; - RD_to_XY (&Xtst, &Ytst, average->R, average->D, &image[im].coords); - - // find the pixel offset - double dX = (Xtst - measure->Xccd); - double dY = (Ytst - measure->Yccd); - - // if dPos is small, we have the right image - double dPos = hypot(dX,dY); - if (dPos > 10.0) { - // if dPos is large, we don't the right imageID, reconstruct the imageID from - // the R,D,photcode & X,Y consistency - im = GetStackSeq (image, average->R, average->D, measure->photcode, measure->Xccd, measure->Yccd); - if (im < 0) { - // still failed, this is bad - if (NmissStack < 10) fprintf (stderr, "missing stack exposure: %f %f : %d %d\n", measure->R, measure->D, measure->imageID, measure->photcode); - NmissStack ++; - continue; - } - // the old imageID was wrong. replace with the new one: - if (logfile) fprintf (logfile, "fix stack %10.6f %10.6f %5d %6.1f %6.1f 0x%08x %d -> %d\n", average->R, average->D, measure->photcode, measure->Xccd, measure->Yccd, measure->detID, measure->imageID, image[im].imageID); - measure->imageID = image[im].imageID; - } - } - uint64_t extID = CreatePSPSStackDetectionID (35, image[im].externID, measure->detID); - if (extID != measure->extID) { - measure->extID = extID; - NfixStackID ++; - } - continue; - } - - // we are only going to repair warp detections - if (!isGPC1warp(measure->photcode)) continue; - - // warp coordinates to confirm warp - double X = measure->Xccd; - double Y = measure->Yccd; - - // check if this detection is far from its correct location - double Rwrp = measure->R; - double Dwrp = measure->D; - - double Rave = average->R; - double Dave = average->D; - - float csdec = cos(Dave * RAD_DEG); - - // find the ra,dec displacement in arcsec: - double dR = 3600.0*(Rwrp - Rave)*csdec; - double dD = 3600.0*(Dwrp - Dave); - - // skip detections which are within a small distance of the expected location - // NOTE: this actually works surprisingly well near the pole - double dPos = hypot(dR,dD); - if (dPos < 2.0) { - int im = getImageByID (measure->imageID); - uint64_t extID = CreatePSPSStackDetectionID (34, image[im].externID, measure->detID); - if (extID != measure->extID) { - measure->extID = extID; - NfixWarpID ++; - } - continue; - } - - // find the corrected warp ID: - int warpSeq = GetWarpSeq (image, measure->t, measure->photcode, Rave, Dave, X, Y); - if (warpSeq < 0) { - if (NmissWarp < 10) fprintf (stderr, "missing warp exposure: %f %f : %d %d\n", Rave, Dave, measure->t, measure->photcode); - NmissWarp ++; - continue; - } - - // assert on coords being TAN? - Coords *imcoords = &image[warpSeq].coords; - - double Rnew, Dnew; - XY_to_RD (&Rnew, &Dnew, X, Y, imcoords); - - Rnew = ohana_normalize_angle_to_midpoint(Rnew, 180.0); - Dnew = ohana_normalize_angle_to_midpoint(Dnew, 0.0); - - // find the ra,dec displacement in arcsec: - double dRnew = 3600.0*(Rnew - Rave)*csdec; - double dDnew = 3600.0*(Dnew - Dave); - - // detections should now be repaired; fail if not - dPos = hypot(dRnew,dDnew); - if (dPos > 5.0) { - fprintf (stderr, "measurement still far from average location: %f %f vs %f %f : %d %d\n", Rave, Dave, Rnew, Dnew, measure->t, measure->photcode); - NbadWarp ++; - } - - measure->imageID = image[warpSeq].imageID; - measure->R = Rnew; - measure->D = Dnew; - - uint64_t extID = CreatePSPSStackDetectionID (34, image[warpSeq].externID, measure->detID); - if (extID != measure->extID) { - measure->extID = extID; - NfixWarpID ++; - } - NfixWarpImageID ++; - } - - fprintf (stderr, "\n"); - fprintf (stderr, "repaired %s : warp image ID %d : ext ID chip %d stack %d warp %d : missed warp %d stack %d : bad warp: %d\n", catalog->filename, NfixWarpImageID, NfixChipID, NfixStackID, NfixWarpID, NmissWarp, NmissStack, NbadWarp); - - return (TRUE); -} Index: trunk/Ohana/src/relphot/src/ResetOps.c =================================================================== --- trunk/Ohana/src/relphot/src/ResetOps.c (revision 41647) +++ trunk/Ohana/src/relphot/src/ResetOps.c (revision 41647) @@ -0,0 +1,193 @@ +# include "relphot.h" + +static DVOAverageFlags photomBitsAverage = + ID_OBJ_EXT | // extended in our data (eg, PS) + ID_OBJ_EXT_ALT | // extended in external data (eg, 2MASS) + ID_OBJ_GOOD | // good-quality measurement in our data (eg,PS) + ID_OBJ_GOOD_ALT | // good-quality measurement in external data (eg, 2MASS) + ID_OBJ_GOOD_STACK | // good-quality object in the stack (> 1 good stack) + ID_OBJ_BEST_STACK | // the primary stack measurement are the best measurements + ID_OBJ_SUSPECT_STACK | // suspect object in the stack (> 1 good or suspect stack, < 2 good) + ID_OBJ_BAD_STACK; // good-quality object in the stack (> 1 good stack) + +// flags used by the photometry analysis (excluding UBERCAL) +static DVOMeasureFlags photomBitsMeasure = + ID_MEAS_NOCAL | // detection ignored for this analysis (photcode, time range) + ID_MEAS_POOR_PHOTOM | // detection is photometry outlier + ID_MEAS_SKIP_PHOTOM | // detection was ignored for photometry measurement + ID_MEAS_AREA | // detetion was outside acceptable area of device + ID_MEAS_SYNTH_MAG | // magnitude is synthetic + ID_MEAS_STACK_PRIMARY | // this stack measurement is in the primary skycell + ID_MEAS_STACK_PHOT_SRC | // this measurement supplied the stack photometry + ID_MEAS_PHOTOM_PSF | // this measurement is used for the mean psf mag + ID_MEAS_PHOTOM_APER | // this measurement is used for the mean ap mag + ID_MEAS_PHOTOM_KRON ; // this measurement is used for the mean kron mag + +// ID_MEAS_PHOTOM_UBERCAL -- externally-supplied zero point from ubercal analysis +// this is set by 'setphot', do not reset here + +// flags used by the photometry analysis (excluding UBERCAL) +// unsigned int secfiltFlags = +// ID_PHOTOM_PASS_0 | // average measured at pass 0 +// ID_PHOTOM_PASS_1 | // average measured at pass 1 +// ID_PHOTOM_PASS_2 | // average measured at pass 2 +// ID_PHOTOM_PASS_3 | // average measured at pass 3 +// ID_PHOTOM_PASS_4 | // average measured at pass 3 +// ID_SECF_USE_SYNTH | // average measured at pass 3 +// ID_SECF_USE_UBERCAL | // average measured at pass 3 +// ID_SECF_OBJ_EXT; // average measured at pass 3 + + +// Used in bcatalog.c, applied to the subset catalogs (bright, good stars used for calibration) +void ResetAverageActivePhotcodes (SecFilt *secfilt) { + + if (!RESET) return; + + int Ns; + + // only loop over the active photcodes + for (Ns = 0; Ns < Nphotcodes; Ns++) { + + int thisCode = photcodes[Ns][0].code; + int Nsec = GetPhotcodeNsec(thisCode); + + dvo_secfilt_init (&secfilt[Nsec], SECFILT_RESET_ALL); + } +} + +// Used in bcatalog.c, applied to the selected measurements which are already selected by photcode +void ResetMeasureZeroPoints (MeasureTiny *measure, off_t Nmeasure, off_t Ncat) { + + if (!RESET) return; + + // only reset Mcal for measures with a matching image + // do not reset Mcal for ubercal images unless explicitly requested + if (measure->dbFlags & ID_MEAS_PHOTOM_UBERCAL) { + if (!KEEP_UBERCAL) { + measure->McalPSF = 0.0; + measure->McalAPER = 0.0; + measure->dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL; + } + } else { + if (getImageEntry (Nmeasure, Ncat) >= 0) { + measure->McalPSF = 0.0; + measure->McalAPER = 0.0; + } + } + if (RESET_FLATCORR) { + measure->Mflat = 0.0; + } + measure->dbFlags &= ~photomBitsMeasure; +} + +// Used in load_images.c. Applied to the subset of selected images used for the calibration. +// This will be only those images in the selected region and only those with the active +// photcode *unless* -use-all-images is selected. +void ResetImages (Image *subset, off_t Nsubset) { + + // reset image values as needed. always allow 'few' images to succeed, if possible (new + // images / detections may have been added + + for (off_t i = 0; i < Nsubset; i++) { + // reset these bits regardless (we will re-determine) + subset[i].flags &= ~(ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW); + if (RESET) { + if (RESET_ZEROPTS) { + if (!KEEP_UBERCAL || !(subset[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) { + subset[i].McalPSF = 0.0; + subset[i].McalAPER = 0.0; + subset[i].dMcal = NAN; + subset[i].flags &= ~ID_IMAGE_PHOTOM_UBERCAL; + } + } + subset[i].flags &= ~(ID_IMAGE_MOSAIC_POOR | ID_IMAGE_NIGHT_POOR); + subset[i].flags &= ~(ID_IMAGE_IMAGE_PHOTCAL | ID_IMAGE_MOSAIC_PHOTCAL | ID_IMAGE_TGROUP_PHOTCAL); + subset[i].ubercalDist = 1000; + } + } +} + +// Used in setMrelFinal.c +void ResetAverageAndMeasure (Catalog *catalog) { + + int Nsecfilt = GetPhotcodeNsecfilt (); + + // XXX I should really deprecate the concept of applying the average + // calculation to a limited set of photcodes. + // for now, just do all photcodes here ( + // as it stands, only stacks are limited by photcode; mean exp and forced warp + // are applied to all Nsecfilt + + for (off_t i = 0; i < catalog->Naverage; i++) { + catalog->average[i].psfQF = NAN; // force recalculation in setMrelCatalog + catalog->average[i].psfQFperf = NAN; // force recalculation in setMrelCatalog + catalog->average[i].stargal = NAN; // force recalculation in setMrelCatalog + catalog->average[i].photFlagsUpper = 0; // reset (will be re-calculated) + catalog->average[i].photFlagsLower = 0; // reset (will be re-calculated) + catalog->average[i].NwarpOK = 0; // reset (will be re-calculated) + + // RESET_FLATCORR independent of RESET? + + if (!RESET) continue; + + catalog->average[i].flags &= ~photomBitsAverage; // reset all photometry bits (but not astrom) + + // Reset secfilt values only for the active photcodes + for (int Ns = 0; Ns < Nphotcodes; Ns++) { + + int thisCode = photcodes[Ns][0].code; + int Nsec = GetPhotcodeNsec(thisCode); + off_t N = Nsecfilt*i+Nsec; + + dvo_secfilt_init (&catalog->secfilt[N], SECFILT_RESET_ALL); + } + + // Reset measure values only for the active photcodes + + off_t m = catalog->average[i].measureOffset; + + for (off_t j = 0; j < catalog->average[i].Nmeasure; j++, m++) { + + /* select measurements for active photcodes */ + PhotCode *code = GetPhotcodebyCode (catalog->measure[m].photcode); + if (!code) continue; + // photcode exists, but is it one of the active ones? + int found = FALSE; + for (int Ns = 0; !found && (Ns < Nphotcodes); Ns++) { + if (code->equiv != photcodes[Ns][0].code) continue; + found = TRUE; + } + if (!found) continue; + + // reset the Mflat value if requested + if (RESET_FLATCORR) { + catalog->measure[m].Mflat = 0.0; + } + + /* select measurements by time */ + if (TimeSelect) { + if (catalog->measure[m].t < TSTART) continue; + if (catalog->measure[m].t > TSTOP) continue; + } + + // only reset Mcal for measures with a matching image + // do not reset Mcal for ubercal images unless explicitly requested + + if (catalog->measure[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) { + if (!KEEP_UBERCAL) { + catalog->measure[m].McalPSF = 0.0; + catalog->measure[m].McalAPER = 0.0; + catalog->measure[m].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL; + } + } else { + if (RESET_ZEROPTS && (getImageEntry (m, 0) >= 0)) { + catalog->measure[m].McalPSF = 0.0; + catalog->measure[m].McalAPER = 0.0; + } + } + catalog->measure[m].dbFlags &= ~photomBitsMeasure; + } + } +} + + Index: trunk/Ohana/src/relphot/src/Shutdown.c =================================================================== --- trunk/Ohana/src/relphot/src/Shutdown.c (revision 41601) +++ trunk/Ohana/src/relphot/src/Shutdown.c (revision 41647) @@ -10,13 +10,12 @@ int Shutdown (char *format, ...) { va_list argp; - char *formatplus; + char formatplus[1024]; - ALLOCATE (formatplus, char, strlen(format)); - strcpy (formatplus, format); - strcat (formatplus, "\n"); + // we cannot allocate this with the ohana allocation tools + // because it could deadlock + snprintf (formatplus, 1024, "%s\n", format); va_start (argp, format); vfprintf (stderr, formatplus, argp); - free (formatplus); va_end (argp); Index: trunk/Ohana/src/relphot/src/StackImageMaps.c =================================================================== --- trunk/Ohana/src/relphot/src/StackImageMaps.c (revision 41601) +++ (revision ) @@ -1,165 +1,0 @@ -# include "relphot.h" - -// generate an index of stack image ra,dec centers sorted indexes to the image IDs - -# define D_NSTACKS 1000 -# define NGROUPS 5 - -typedef struct { - unsigned short photcode; - double *Rcenter; - double *Dcenter; - int *imageID; - int *imageSeq; - int Nstacks; - int NSTACKS; -} StackGroupType; - -// one group for each photcode -StackGroupType *stackgroup = NULL; - -int MakeStackIndex (void) { - - off_t Nimage; - Image *image = getimages (&Nimage, NULL); - - ALLOCATE (stackgroup, StackGroupType, NGROUPS); - - for (int i = 0; i < NGROUPS; i++) { - stackgroup[i].photcode = 11000 + i*100; - - stackgroup[i].Nstacks = 0; - stackgroup[i].NSTACKS = D_NSTACKS; - ALLOCATE (stackgroup[i].Rcenter, double, stackgroup[i].NSTACKS); - ALLOCATE (stackgroup[i].Dcenter, double, stackgroup[i].NSTACKS); - ALLOCATE (stackgroup[i].imageID, int, stackgroup[i].NSTACKS); - ALLOCATE (stackgroup[i].imageSeq, int, stackgroup[i].NSTACKS); - } - - // first generate a list of obstime values - for (off_t i = 0; i < Nimage; i++) { - if (!isGPC1stack(image[i].photcode)) continue; - - int Ng = -1; - switch (image[i].photcode) { - case 11000: Ng = 0; break; - case 11100: Ng = 1; break; - case 11200: Ng = 2; break; - case 11300: Ng = 3; break; - case 11400: Ng = 4; break; - default: myAbort("impossible photcode"); - } - myAssert (Ng < NGROUPS, "oops"); - - int N = stackgroup[Ng].Nstacks; - XY_to_RD (&stackgroup[Ng].Rcenter[N], &stackgroup[Ng].Dcenter[N], 0.5*image[i].NX, 0.5*image[i].NY, &image[i].coords); - stackgroup[Ng].Rcenter[N] = ohana_normalize_angle (stackgroup[Ng].Rcenter[N]); - - stackgroup[Ng].imageSeq[N] = i; - stackgroup[Ng].imageID[N] = image[i].imageID; - - stackgroup[Ng].Nstacks ++; - - if (stackgroup[Ng].Nstacks >= stackgroup[Ng].NSTACKS) { - stackgroup[Ng].NSTACKS += D_NSTACKS; - REALLOCATE (stackgroup[Ng].Rcenter, double, stackgroup[Ng].NSTACKS); - REALLOCATE (stackgroup[Ng].Dcenter, double, stackgroup[Ng].NSTACKS); - REALLOCATE (stackgroup[Ng].imageID, int, stackgroup[Ng].NSTACKS); - REALLOCATE (stackgroup[Ng].imageSeq, int, stackgroup[Ng].NSTACKS); - } - } - - // sort the list: - for (int i = 0; i < NGROUPS; i++) { - sort_by_ra (stackgroup[i].Rcenter, stackgroup[i].Dcenter, stackgroup[i].imageID, stackgroup[i].imageSeq, stackgroup[i].Nstacks); - } - - return TRUE; -} - -// find the stack which yields the given Rstk, Dstk for the given X,Y (and photcode) -int GetStackSeq (Image *image, double Rstk, double Dstk, unsigned short photcode, float X, float Y) { - - if (!stackgroup) return -1; - - // we have the stack centers; find all stacks within 0.3 degrees of this point - Rstk = ohana_normalize_angle (Rstk); - double dD = 0.3; - double dR = dD / cos(RAD_DEG*Dstk); - double Rmin = Rstk - dR; - double Rmax = Rstk + dR; - double Dmin = Dstk - dD; - double Dmax = Dstk + dD; - - int Ng = (photcode / 100) % 10; - int N = ohana_bisection_double (stackgroup[Ng].Rcenter, stackgroup[Ng].Nstacks, Rmin); - - double dPosMin = NAN; - int nPosMin = -1; - - for (; N < stackgroup[Ng].Nstacks; N++) { - if (stackgroup[Ng].Dcenter[N] < Dmin) continue; - if (stackgroup[Ng].Dcenter[N] > Dmax) continue; - if (stackgroup[Ng].Rcenter[N] > Rmax) break; - - int im = stackgroup[Ng].imageSeq[N]; - - // project to image pixels - double Xtst, Ytst; - RD_to_XY (&Xtst, &Ytst, Rstk, Dstk, &image[im].coords); - - // find the pixel offset - double dX = (Xtst - X); - double dY = (Ytst - Y); - - // if dPos is small, we have the right image - double dPos = hypot(dX,dY); - if (dPos < 10.0) return im; // 10 pixels == 2.0 arcsec - - // check for multiple possible matches?? - - // record the min pos so we can see how things fail - if (isnan(dPosMin)) { - dPosMin = dPos; - nPosMin = im; - } else { - if (dPos < dPosMin) { - dPosMin = dPos; - nPosMin = im; - } - } - } - - // we did not find the right image, return what we found (-1 or closest match) - return nPosMin; -} - -void FreeStackGroups (void) { - - if (!stackgroup) return; - - for (int i = 0; i < NGROUPS; i++) { - free (stackgroup[i].Rcenter); - free (stackgroup[i].Dcenter); - free (stackgroup[i].imageID); - free (stackgroup[i].imageSeq); - } - free (stackgroup); -} - -void sort_by_ra (double *R, double *D, int *I, int *S, int N) { - -# define SWAPFUNC(A,B){ double dtmp; int itmp; \ - dtmp = R[A]; R[A] = R[B]; R[B] = dtmp; \ - dtmp = D[A]; D[A] = D[B]; D[B] = dtmp; \ - itmp = I[A]; I[A] = I[B]; I[B] = itmp; \ - itmp = S[A]; S[A] = S[B]; S[B] = itmp; \ - } -# define COMPARE(A,B)(R[A] < R[B]) - - OHANA_SORT (N, COMPARE, SWAPFUNC); - -# undef SWAPFUNC -# undef COMPARE - -} Index: trunk/Ohana/src/relphot/src/StarOps.c =================================================================== --- trunk/Ohana/src/relphot/src/StarOps.c (revision 41601) +++ trunk/Ohana/src/relphot/src/StarOps.c (revision 41647) @@ -1,3 +1,5 @@ # include "relphot.h" + +# define STAR_BAD (ID_SECF_STAR_FEW | ID_SECF_STAR_POOR) enum {THREAD_RUN, THREAD_DONE}; @@ -8,5 +10,4 @@ Catalog *catalog; int Ncatalog; - FlatCorrectionTable *flatcorr; SetMrelInfo summary; } ThreadInfo; @@ -14,6 +15,5 @@ // ** internal functions: void *setMrel_worker (void *data); -int setMrel_threaded (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr); -int setMrel_catalog (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt); +int setMrel_threaded (Catalog *catalog, int Ncatalog); int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure); @@ -292,16 +292,21 @@ // setMrel and setMrelOutput are extremely similar, but have slightly different implications: -// * setMrel uses the internal Tiny structures only -// * setMrelOutput skips stars for which there are too few good measurements -// * setMrelOutput is meant to be called repeatedly, relaxing the criteria for 'good' on each pass +// * both call setMrelCatalog, but setMrel set isSetMrelFinal to FALSE +// * setMrel only uses the internal Tiny structures +// * setMrel skips stars for which there are too few good measurements +// -> (setMrelOutput assigns an average mag for all objects, even those with only 1 measurement) +// * setMrelOutput sets psfQF, psfQFperf, extNsigma, rank info, etc +// * setMrelOutput sets average Map // * setMrelOutput updates 2MASS average flags // * setMrelOutput updates average EXT flags (PS1 and 2MASS) - -int setMrel (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) { +// * setMrelOutput sets average stack & warp photometry (if requested), setMrel always skips those +// * setMrelOutput is NOT threaded it reads/writes catalogs one at a time + +int setMrel (Catalog *catalog, int Ncatalog) { int i; if (NTHREADS) { - int status = setMrel_threaded (catalog, Ncatalog, flatcorr); + int status = setMrel_threaded (catalog, Ncatalog); return status; } @@ -314,5 +319,5 @@ for (i = 0; i < Ncatalog; i++) { - setMrel_catalog_alt (catalog, i, FALSE, flatcorr, &results, Nsecfilt); + setMrelCatalog (catalog, i, FALSE, &results, Nsecfilt); SetMrelInfoAccum (&summary, &results); } @@ -325,6 +330,5 @@ } -// setMrelOutput is NOT threaded because we read/write catalogs one at a time (for now) -int setMrelOutput (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) { +int setMrelOutput (Catalog *catalog, int Ncatalog) { int i; @@ -337,5 +341,5 @@ for (i = 0; i < Ncatalog; i++) { - setMrel_catalog_alt (catalog, i, TRUE, flatcorr, &results, Nsecfilt); + setMrelCatalog (catalog, i, TRUE, &results, Nsecfilt); SetMrelInfoAccum (&summary, &results); } @@ -347,5 +351,5 @@ } -int setMrel_threaded (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) { +int setMrel_threaded (Catalog *catalog, int Ncatalog) { int i; @@ -374,5 +378,4 @@ threadinfo[i].catalog = catalog; threadinfo[i].Ncatalog = Ncatalog; - threadinfo[i].flatcorr = flatcorr; SetMrelInfoInit (&threadinfo[i].summary, FALSE); pthread_create (&threads[i], NULL, setMrel_worker, &threadinfo[i]); @@ -435,7 +438,6 @@ Catalog *catalog = threadinfo->catalog; - FlatCorrectionTable *flatcorr = threadinfo->flatcorr; - - setMrel_catalog_alt (catalog, i, FALSE, flatcorr, &results, Nsecfilt); + + setMrelCatalog (catalog, i, FALSE, &results, Nsecfilt); SetMrelInfoAccum (&threadinfo->summary, &results); } @@ -444,13 +446,4 @@ return NULL; } - -# define SKIP_THIS_MEAS(REASON) { \ - catalog[Nc].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM; \ - if (catalog[Nc].measure) { \ - catalog[Nc].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM; \ - } \ - results->REASON ++; \ - continue; } - int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure) { @@ -473,13 +466,10 @@ } -/* set measure.Mcal for all measures except ID_MEAS_NOCAL and ID_IMAGE_PHOTOM_NOCAL */ -int setMcalOutput (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) { +/* set measure.Mcal for all measures except ID_MEAS_NOCAL */ +int setMcalOutput (Catalog *catalog, int Ncatalog) { int i; off_t j, k, m; - MEAS_BAD = ID_MEAS_NOCAL; - IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL; - for (i = 0; i < Ncatalog; i++) { for (j = 0; j < catalog[i].Naverage; j++) { @@ -487,14 +477,18 @@ m = catalog[i].averageT[j].measureOffset; for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) { - if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; + if (catalog[i].measureT[m].dbFlags & ID_MEAS_NOCAL) continue; float McalPSF = getMcal (m, i, MAG_CLASS_PSF); float McalAPER = getMcal (m, i, MAG_CLASS_APER); if (isnan(McalPSF)) continue; - float Mmos = getMmos (m, i); - if (isnan(Mmos)) continue; - float Mgrp = getMgrp (m, i, catalog[i].measureT[m].airmass, NULL); - if (isnan(Mgrp)) continue; - float Mgrid = getMgrid (m, i); - if (isnan(Mgrid)) continue; + + // Mmos, Mgrp have already been transfered to Mcal and should be zero here + // float Mmos = getMmos (m, i); + // if (isnan(Mmos)) continue; + // float Mgrp = getMgrp (m, i, catalog[i].measureT[m].airmass, NULL); + // if (isnan(Mgrp)) continue; + + // Mgrid is applied to Mflat and should be ignored here + // float Mgrid = getMgridTiny (&catalog[i].measureT[m]); + // if (isnan(Mgrid)) continue; // Note that this operation is setting measure->McalAPER to image->McalAPER only @@ -509,8 +503,6 @@ // set the output calibration - // float Moff = Mmos + Mgrp + Mgrid; - float Moff = Mmos + Mgrid; // Mgrp has already been applied to the image.Mcal - catalog[i].measure[m].McalPSF = McalPSF + Moff; - catalog[i].measure[m].McalAPER = useStackAper ? McalAPER + Moff : McalPSF + Moff; + catalog[i].measure[m].McalPSF = McalPSF; + catalog[i].measure[m].McalAPER = useStackAper ? McalAPER : McalPSF; if (catalog[i].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) { @@ -584,10 +576,16 @@ } +static float MinChiSqLim = NAN; +static float MaxChiSqLim = NAN; +static float MinScatterLim = NAN; +static float MaxScatterLim = NAN; + void clean_stars (Catalog *catalog, int Ncatalog) { - int i, j, Ndel, Nave, Ntot, mark, Ns, Nscat, Nchi, Nnan; - float dM; - double MaxScatter, MaxChisq; - double *xlist, *slist, *dlist; + int Ndel, Nave, Ntot, Nscat, Nchi, Nnan; + double *xlist, *slist; + + if (isnan (MaxChiSqLim)) MaxChiSqLim = STAR_CHISQ; + if (isnan (MaxScatterLim)) MaxScatterLim = STAR_SCATTER; StatType stats; @@ -597,241 +595,66 @@ /* find Mchisq median -> ChiSq lim must be > median */ - for (i = Ntot = 0; i < Ncatalog; i++) { + for (int i = Ntot = 0; i < Ncatalog; i++) { Ntot += catalog[i].Naverage; } ALLOCATE (xlist, double, Ntot); ALLOCATE (slist, double, Ntot); - ALLOCATE (dlist, double, Ntot); int Nsecfilt = GetPhotcodeNsecfilt (); // eliminate bad stars using the stats for a single secfilt at a time - for (Ns = 0; Ns < Nphotcodes; Ns ++) { + for (int Ns = 0; Ns < Nphotcodes; Ns ++) { int thisCode = photcodes[Ns][0].code; int Nsec = GetPhotcodeNsec(thisCode); - for (i = Ntot = 0; i < Ncatalog; i++) { - for (j = 0; j < catalog[i].Naverage; j++) { - if ( catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD ) continue; + for (int i = Ntot = 0; i < Ncatalog; i++) { + for (int j = 0; j < catalog[i].Naverage; j++) { + if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; float Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq; if (isnan(Mchisq)) continue; xlist[Ntot] = Mchisq; slist[Ntot] = catalog[i].secfilt[Nsecfilt*j+Nsec].dMpsfChp; - dlist[Ntot] = 1; Ntot ++; } } - - liststats (xlist, dlist, NULL, Ntot, &stats); - MaxChisq = MAX (STAR_CHISQ, 2*stats.median); - - liststats (slist, dlist, NULL, Ntot, &stats); - MaxScatter = MAX (STAR_SCATTER, 2*stats.median); - fprintf (stderr, "Max Scatter: %f, Max Chisq: %f\n", MaxScatter, MaxChisq); + + liststats (xlist, NULL, NULL, Ntot, &stats); + float ChiSqUpper90 = stats.Upper90; + if (isnan (MinChiSqLim)) MinChiSqLim = 2.0*stats.median; // chi-square cut cannot fall below this value (even if this is > MaxChiSqLim) + float ChiSqLimit = MAX(MinChiSqLim, MIN(MaxChiSqLim, ChiSqUpper90)); // chi-square cut should be between MinChiSqLim and MaxChiSqLim + + liststats (slist, NULL, NULL, Ntot, &stats); + float ScatterUpper90 = stats.Upper90; + if (isnan (MinScatterLim)) MinScatterLim = 2.0*stats.median; // scatter cut cannto fall below this value (even if this is > MaxScatterLim) + float ScatterLimit = MAX(MinScatterLim, MIN(MaxScatterLim, ScatterUpper90)); + + fprintf (stderr, "STARS: ChiSqLimit: %f, ScatterLimit: %f | ChiSquare Upper 90: %f, Scatter Upper 90: %f\n", ChiSqLimit, ScatterLimit, ChiSqUpper90, ScatterUpper90); Ndel = Nave = Nscat = Nnan = Nchi = 0; - for (i = 0; i < Ncatalog; i++) { - for (j = 0; j < catalog[i].Naverage; j++) { - dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dMpsfChp; + for (int i = 0; i < Ncatalog; i++) { + for (int j = 0; j < catalog[i].Naverage; j++) { + float dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dMpsfChp; float Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq; - mark = (dM > MaxScatter) || (isnan(Mchisq)) || (Mchisq > MaxChisq); + int mark = (dM > ScatterLimit) || (isnan(Mchisq)) || (Mchisq > ChiSqLimit); if (mark) { - catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_OBJ_POOR; + catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_STAR_POOR; Ndel ++; - if (dM > MaxScatter) { Nscat ++; } - if (isnan(Mchisq)) { Nnan ++; } - if (Mchisq > MaxChisq) { Nchi ++; } + if (dM > ScatterLimit) { Nscat ++; } + if (isnan(Mchisq)) { Nnan ++; } + if (Mchisq > ChiSqLimit) { Nchi ++; } } else { - catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_OBJ_POOR; + catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_SECF_STAR_POOR; } Nave ++; } } - fprintf (stderr, "%d stars marked variable (%d scat, %d nan, %d chi), %d total\n", Ndel, Nscat, Nnan, Nchi, Nave); + fprintf (stderr, "%d of %d stars marked variable (%d scat, %d nan, %d chi)\n", Ndel, Nave, Nscat, Nnan, Nchi); } free (xlist); free (slist); - free (dlist); -} - -// clean_measures examines the stats for a single star. It measures the INNER 50% mean -// and sigma, it then re-measures the mean and sigma using all stars with NSIGMA_CLIP (3) -// sigma of the INNER 50% mean. it then flags any measurements which are more than -// NSIGMA_REJECT (5) sigma of the mean - -// this function only operations on the PSF magnitudes - -// XXX EAM clean_measures is no longer used -- we rely on IRLS to down-weight outliers - -# define NSIGMA_CLIP 3.0 -# define NSIGMA_REJECT 5.0 -void clean_measures (Catalog *catalog, int Ncatalog, int final, FlatCorrectionTable *flatcorr) { - - off_t j, k, m, Nmax, Ndel, Nave; - int i, N, image_bad, TOOFEW; - off_t *ilist; - double *tlist, *list, *dlist; - float Msys, Mcal, Mmos, Mgrp, Mgrid; - int Ncal, Nmos, Ngrp, Ngrid, Nfew; - - int Nsecfilt = GetPhotcodeNsecfilt (); - - if (VERBOSE) fprintf (stderr, "marking poor measures\n"); - /* Nmeasure is now different, need to reallocate */ - Nmax = 0; - for (i = 0; i < Ncatalog; i++) { - for (j = 0; j < catalog[i].Naverage; j++) { - Nmax = MAX (Nmax, catalog[i].averageT[j].Nmeasure); - } - } - ALLOCATE (list, double, MAX (1, Nmax)); - ALLOCATE (dlist, double, MAX (1, Nmax)); - ALLOCATE (ilist, off_t, Nmax); - ALLOCATE (tlist, double, Nmax); - - /* it makes no sense to mark 3-sigma outliers with <5 measurements */ - TOOFEW = MAX (5, STAR_TOOFEW); - - // stats structures for inner and full stats - StatType instats, stats; - liststats_setmode (&instats, "INNER_MEAN"); - liststats_setmode (&stats, "MEAN"); - - Ndel = Nave = 0; - Ncal = Nmos = Ngrp = Ngrid = Nfew = 0; - for (i = 0; i < Ncatalog; i++) { - for (j = 0; j < catalog[i].Naverage; j++) { - - int Ns; - for (Ns = 0; Ns < Nphotcodes; Ns++) { - - int thisCode = photcodes[Ns][0].code; - int Nsec = GetPhotcodeNsec(thisCode); - - /* on final processing, skip stars already measured */ - if (final && catalog[i].found_t[Nsecfilt*j + Nsec]) continue; - - /* skip bad stars to prevent them from becoming good (on inner sample) */ - if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; - - /* accumulate list of valid measurements */ - m = catalog[i].averageT[j].measureOffset; - N = 0; - for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) { - // skip measurements that do not match the current photcode - int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode); - if (ecode != thisCode) { continue; } - - // SKIP gpc1 stack data - if (isGPC1stack(catalog[i].measureT[m].photcode)) continue; - - // SKIP gpc1 forced-warp data - if (isGPC1warp(catalog[i].measureT[m].photcode)) continue; - - // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers - Mcal = getMcal (m, i, MAG_CLASS_PSF); - if (isnan(Mcal)) { Ncal ++; continue; } - Mmos = getMmos (m, i); - if (isnan(Mmos)) { Nmos ++; continue; } - Mgrp = getMgrp (m, i, catalog[i].measureT[m].airmass, NULL); - if (isnan(Mgrp)) { Ngrp ++; continue; } - Mgrid = getMgrid (m, i); - if (isnan(Mgrid)) { Ngrid ++; continue; } - - Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt], MAG_CLASS_PSF); - list[N] = Msys - Mcal - Mmos - Mgrp - Mgrid; - dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR); - N++; - } - if (N <= TOOFEW) { - Nfew ++; - continue; - } - - /* 3-sigma clip based on stats of inner 50% */ - - // calculated mean of inner 50% - liststats (list, dlist, NULL, N, &instats); - instats.sigma = MAX (MIN_ERROR, instats.sigma); /* if measurements agree too well, sigma -> 0.0 */ - - // ignore entries > 3sigma from inner mean - for (k = m = 0; k < N; k++) { - if (fabs (list[k] - instats.median) < NSIGMA_CLIP*instats.sigma) { - list[m] = list[k]; - m++; - } - } - // recalculate the mean & sigma of the accepted measurements - liststats (list, dlist, NULL, m, &stats); - stats.sigma = MAX (MIN_ERROR, stats.sigma); - - /* apply to list of all relevant measurements, including IMAGE_POOR & IMAGE_FEW */ - image_bad = IMAGE_BAD; - IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL; - m = catalog[i].averageT[j].measureOffset; - N = 0; - for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) { - // skip measurements that do not match the current photcode - int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode); - if (ecode != thisCode) { continue; } - - // SKIP gpc1 stack data - if (isGPC1stack(catalog[i].measureT[m].photcode)) continue; - - // SKIP gpc1 forced-warp data - if (isGPC1warp(catalog[i].measureT[m].photcode)) continue; - - // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers - Mcal = getMcal (m, i, MAG_CLASS_PSF); - if (isnan(Mcal)) continue; - Mmos = getMmos (m, i); - if (isnan(Mmos)) continue; - Mgrid = getMgrid (m, i); - if (isnan(Mgrid)) continue; - - Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt], MAG_CLASS_PSF); - list[N] = Msys - Mcal - Mmos - Mgrid; - dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR); - ilist[N] = m; - N++; - Nave ++; - } - if (N < TOOFEW) continue; - - /* mark bad measures (> 5 sigma deviant) */ - for (k = 0; k < N; k++) { - // treat the scatter of the star as a systematic term. this is a bit of an - // over-estimage (a perfect Gauss distribution with perfect errors would have - // mySigma = sqrt(2) too large) - float mySigma = hypot (stats.sigma, dlist[k]); - if (fabs (list[k] - stats.median) > NSIGMA_REJECT*mySigma) { - catalog[i].measureT[ilist[k]].dbFlags |= ID_MEAS_POOR_PHOTOM; - if (final) { - // for the final pass, we have a duplicate set of values in measure and measureT - catalog[i].measure[ilist[k]].dbFlags |= ID_MEAS_POOR_PHOTOM; - } - Ndel ++; - } else { - catalog[i].measureT[ilist[k]].dbFlags &= ~ID_MEAS_POOR_PHOTOM; - if (final) { - // for the final pass, we have a duplicate set of values in measure and measureT - catalog[i].measure[ilist[k]].dbFlags &= ~ID_MEAS_POOR_PHOTOM; - } - } - } - IMAGE_BAD = image_bad; - } - } - } - if (VERBOSE) fprintf (stderr, OFF_T_FMT" measures marked poor, "OFF_T_FMT" total\n", Ndel, Nave); - - free (list); - free (dlist); - free (ilist); - free (tlist); -} - -StatType statsStarN (Catalog *catalog, int Ncatalog, int Nsec, int seccode, FlatCorrectionTable *flatcorr) { +} + +StatType statsStarN (Catalog *catalog, int Ncatalog, int Nsec, int seccode) { off_t j, k, m, Ntot; @@ -869,5 +692,5 @@ Mmos = getMmos (m, i); if (isnan(Mmos)) { continue; } - Mgrid = getMgrid (m, i); + Mgrid = getMgridTiny (&catalog[i].measureT[m]); if (isnan(Mgrid)) { continue;} N++; @@ -1105,6 +928,6 @@ float MrelPSF = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); // average magnitude - // float Mgrid = getMgrid (m, c); // camera offset (deprecated?) - // float Mflat = getMflat (m, c, flatcorr, catalog); // flat-field correction + // float Mgrid = getMgridTiny (&catalog[c].measureT[m]); // camera offset (deprecated?) + // float Mflat = getMflat (m, c, catalog); // flat-field correction off_t n = catalog[c].measureT[m].averef; @@ -1121,2 +944,11 @@ fclose (fout); } + +// maybe make this a macro +float getMflat (off_t meas, int cat, Catalog *catalog) { + float offset = catalog[cat].measureT[meas].Mflat; + if (!isfinite(offset)) { + offset = 0.0; + } + return (offset); +} Index: trunk/Ohana/src/relphot/src/TGroupOps.c =================================================================== --- trunk/Ohana/src/relphot/src/TGroupOps.c (revision 41601) +++ trunk/Ohana/src/relphot/src/TGroupOps.c (revision 41647) @@ -2,5 +2,5 @@ void assignMosaicsToTGroups (void); -int save_test_night_measures (FILE *fout, TGroup *myTGroup, Catalog *catalog, FlatCorrectionTable *flatcorr); +int save_test_night_measures (FILE *fout, TGroup *myTGroup, Catalog *catalog); // tgroupTimes carries the times of the tgroups (initially, just the photometric nights) @@ -647,11 +647,10 @@ Catalog *catalog; Image *image; - FlatCorrectionTable *flatcorr; SetMgrpInfo info; } ThreadInfo; -int setMgrp_tgroup (TGroup *tgroup, off_t Nmos, Image *image, Catalog *catalog, SetMgrpInfo *info, FlatCorrectionTable *flatcorr); +int setMgrp_tgroup (TGroup *tgroup, off_t Nmos, Image *image, Catalog *catalog, SetMgrpInfo *info); void *setMgrp_worker (void *data); -int setMgrp_threaded (Catalog *catalog, FlatCorrectionTable *flatcorr); +int setMgrp_threaded (Catalog *catalog); void SetMgrpInfoInit (SetMgrpInfo *info, off_t Nmax, int allocLists) { @@ -727,5 +726,5 @@ } -int setMgrp (Catalog *catalog, FlatCorrectionTable *flatcorr) { +int setMgrp (Catalog *catalog) { off_t N; @@ -738,5 +737,5 @@ // so do not run setMgrp in threaded mode if PLOTSTUFF is set if (NTHREADS && !PLOTSTUFF) { - int status = setMgrp_threaded (catalog, flatcorr); + int status = setMgrp_threaded (catalog); return status; } @@ -758,5 +757,5 @@ TGroup *tgroup = tgroupTimes[i][0].byCode; for (off_t j = 0; j < tgroupTimes[i][0].nCode; j++) { - setMgrp_tgroup (&tgroup[j], i, image, catalog, &info, flatcorr); + setMgrp_tgroup (&tgroup[j], i, image, catalog, &info); // fprintf (stderr, "TGROUP time %.0f photcode %d Mcal: %f, dK: %f, dMcal: %f, sigma: %f, chisq: %f, %d of %d, limiting negative clouds to %f\n", // ohana_sec_to_mjd(tgroupTimes[i][0].start), tgroup[j].photcode, tgroup[j].McalPSF, tgroup[j].dKlam, tgroup[j].dMcal, tgroup[j].stdev, tgroup[j].McalChiSq, tgroup[j].nFitPhotom, tgroup[j].nValPhotom, CLOUD_TOLERANCE); @@ -781,5 +780,5 @@ // 'tgroup' is a pointer to the current tgroup of interest (entry Ngrp) -int setMgrp_tgroup (TGroup *myTGroup, off_t Ngrp, Image *image, Catalog *catalog, SetMgrpInfo *info, FlatCorrectionTable *flatcorr) { +int setMgrp_tgroup (TGroup *myTGroup, off_t Ngrp, Image *image, Catalog *catalog, SetMgrpInfo *info) { off_t j; @@ -826,5 +825,5 @@ fout = fopen (filename, "w"); fprintf (fout, "# tgroup %s (%.0f)\n", ohana_sec_to_date(tgroup->start), mjdStart); - save_test_night_measures (fout, myTGroup, catalog, flatcorr); + save_test_night_measures (fout, myTGroup, catalog); fclose (fout); } @@ -853,5 +852,5 @@ continue; } - float Mgrid = getMgrid (m, c); + float Mgrid = getMgridTiny (&catalog[c].measureT[m]); if (isnan(Mgrid)) { info->Ngrid ++; @@ -875,7 +874,8 @@ // apply that offset as well here for this image (in other words, each detection is // being compared to the model, excluding the zero point, Mcal. The model includes - // the flat-correction. NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat) - - float Mflat = getMflat (m, c, flatcorr, catalog); + // the flat-correction. NOTE the sign of Mflat (Image.Mcal = Measure.Mcal + Mflat) + // this was inconsistent with PhotRel pre-r41606 + + float Mflat = getMflat (m, c, catalog); off_t n = catalog[c].measureT[m].averef; @@ -900,5 +900,5 @@ assert (Nbright >= 0); - float Moff = Mcal + Mmos + Mgrid - Mflat; + float Moff = Mcal + Mmos + Mgrid + Mflat; psfStars->alldata-> yVector[N] = MsysPSF - MrelPSF - Moff; @@ -1000,5 +1000,5 @@ } -int save_test_night_measures (FILE *fout, TGroup *myTGroup, Catalog *catalog, FlatCorrectionTable *flatcorr) { +int save_test_night_measures (FILE *fout, TGroup *myTGroup, Catalog *catalog) { int Nsecfilt = GetPhotcodeNsecfilt (); @@ -1013,7 +1013,7 @@ float Mmos = getMmos (m, c); // mosaic zero point - float Mgrid = getMgrid (m, c); // camera offset (deprecated?) + float Mgrid = getMgridTiny (&catalog[c].measureT[m]); // camera offset (deprecated?) float MrelPSF = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); // average magnitude - float Mflat = getMflat (m, c, flatcorr, catalog); // flat-field correction + float Mflat = getMflat (m, c, catalog); // flat-field correction off_t n = catalog[c].measureT[m].averef; @@ -1023,5 +1023,5 @@ // for historical reasons, Mflat is defined with the wrong sign - float delta = MsysPSF - MrelPSF - Mcal - Mmos - Mgrid + Mflat - myTGroup->McalPSF; + float delta = MsysPSF - MrelPSF - Mcal - Mmos - Mgrid - Mflat - myTGroup->McalPSF; int isBad = (catalog[c].measureT[m].dbFlags & MEAS_BAD); @@ -1037,5 +1037,5 @@ } -int setMgrp_threaded (Catalog *catalog, FlatCorrectionTable *flatcorr) { +int setMgrp_threaded (Catalog *catalog) { int i; @@ -1075,5 +1075,4 @@ threadinfo[i].catalog = catalog; threadinfo[i].image = image; - threadinfo[i].flatcorr = flatcorr; // we do NOT allocate the arrays here, we only supply basic info used in the threads @@ -1124,4 +1123,6 @@ free (threadinfo); + // XXX not allocated SetMgrpInfoFree (&summary); + npass_output ++; @@ -1145,10 +1146,9 @@ Catalog *catalog = threadinfo->catalog; - FlatCorrectionTable *flatcorr = threadinfo->flatcorr; Image *image = threadinfo->image; TGroup *tgroup = tgroupTimes[i][0].byCode; for (int j = 0; j < tgroupTimes[i][0].nCode; j++) { - setMgrp_tgroup (&tgroup[j], i, image, catalog, &results, flatcorr); + setMgrp_tgroup (&tgroup[j], i, image, catalog, &results); SetMgrpInfoAccum (&threadinfo->info, &results); } @@ -1172,11 +1172,9 @@ ALLOCATE (dlist, double, NtgroupTimes*Nphotcodes); - int NIGHT_BAD = ID_IMAGE_NIGHT_POOR | ID_IMAGE_PHOTOM_FEW; - int n = 0; for (int i = 0; i < NtgroupTimes; i++) { TGroup *tgroup = tgroupTimes[i][0].byCode; for (int j = 0; j < tgroupTimes[i][0].nCode; j++) { - if (tgroup[j].flags & NIGHT_BAD) continue; + if (tgroup[j].flags & (ID_IMAGE_NIGHT_POOR | ID_IMAGE_PHOTOM_FEW)) continue; list[n] = tgroup[j].McalPSF; dlist[n] = 1; @@ -1209,5 +1207,5 @@ TGroup *tgroup = tgroupTimes[i][0].byCode; for (int j = 0; j < tgroupTimes[i][0].nCode; j++) { - if (tgroup[j].flags & IMAGE_BAD) continue; + if (tgroup[j].flags & (ID_IMAGE_NIGHT_POOR | ID_IMAGE_PHOTOM_FEW)) continue; list[n] = tgroup[j].dMcal; dlist[n] = 1; @@ -1240,5 +1238,5 @@ TGroup *tgroup = tgroupTimes[i][0].byCode; for (int j = 0; j < tgroupTimes[i][0].nCode; j++) { - if (tgroup[j].flags & IMAGE_BAD) continue; + if (tgroup[j].flags & (ID_IMAGE_NIGHT_POOR | ID_IMAGE_PHOTOM_FEW)) continue; list[n] = tgroup[j].McalChiSq; dlist[n] = 1; @@ -1267,6 +1265,8 @@ */ -static float MinMaxChiSq = NAN; -static float MinMaxScatter = NAN; +static float MinChiSqLim = NAN; +static float MaxChiSqLim = NAN; +static float MinScatterLim = NAN; +static float MaxScatterLim = NAN; void clean_tgroups () { @@ -1277,4 +1277,7 @@ if (VERBOSE) fprintf (stderr, "marking poor tgroups\n"); + + if (isnan (MaxChiSqLim)) MaxChiSqLim = NIGHT_CHISQ; + if (isnan (MaxScatterLim)) MaxScatterLim = NIGHT_SCATTER; ALLOCATE (mlist, double, NtgroupTimes*Nphotcodes); @@ -1300,13 +1303,13 @@ liststats (mlist, NULL, NULL, N, &stats); float ChiSqUpper90 = stats.Upper90; - if (isnan (MinMaxChiSq)) MinMaxChiSq = 2.0*stats.median; - float MaxChiSq = MAX (MinMaxChiSq, ChiSqUpper90); + if (isnan (MinChiSqLim)) MinChiSqLim = 2.0*stats.median; // chi-square cut cannot fall below this value (even if this is > MaxChiSqLim) + float ChiSqLimit = MAX(MinChiSqLim, MIN(MaxChiSqLim, ChiSqUpper90)); // chi-square cut should be between MinChiSqLim and MaxChiSqLim liststats (slist, NULL, NULL, N, &stats); float ScatterUpper90 = stats.Upper90; - if (isnan (MinMaxScatter)) MinMaxScatter = 2.0*stats.median; - float MaxScatter = MAX (MinMaxScatter, ScatterUpper90); - - fprintf (stderr, "TGROUPS: Max ChiSq: %f, Max Scatter: %f | ChiSquare Upper 90: %f, Scatter Upper 90: %f\n", MaxChiSq, MaxScatter, ChiSqUpper90, ScatterUpper90); + if (isnan (MinScatterLim)) MinScatterLim = 2.0*stats.median; // scatter cut cannto fall below this value (even if this is > MaxScatterLim) + float ScatterLimit = MAX(MinScatterLim, MIN(MaxScatterLim, ScatterUpper90)); + + fprintf (stderr, "TGROUPS: ChiSqLimit: %f, ScatterLimit: %f | ChiSquare Upper 90: %f, Scatter Upper 90: %f\n", ChiSqLimit, ScatterLimit, ChiSqUpper90, ScatterUpper90); int Ntotal = 0, Npoor = 0, Nmark = 0, Nscatter = 0, Nchisq = 0, NfewNights = 0, NfewExp = 0; @@ -1337,10 +1340,10 @@ } // scatter too large - if (tgroup[j].stdev > MaxScatter) { + if (tgroup[j].stdev > ScatterLimit) { mark = TRUE; Nscatter ++; } // chisq too large - if (tgroup[j].McalChiSq > MaxChiSq) { + if (tgroup[j].McalChiSq > ChiSqLimit) { mark = TRUE; Nchisq ++; @@ -1394,5 +1397,5 @@ c = tgroup[i].catalog[j]; - if (catalog[c].measureT[m].dbFlags & (ID_MEAS_AREA | ID_MEAS_NOCAL)) continue; + if (catalog[c].measureT[m].dbFlags & MEAS_BAD) continue; // ave = catalog[c].measureT[m].averef; @@ -1501,6 +1504,6 @@ float MrelPSF = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); // average magnitude - // float Mgrid = getMgrid (m, c); // camera offset (deprecated?) - // float Mflat = getMflat (m, c, flatcorr, catalog); // flat-field correction + // float Mgrid = getMgridTiny (&catalog[c].measureT[m]); // camera offset (deprecated?) + float Mflat = getMflat (m, c, catalog); // flat-field correction off_t n = catalog[c].measureT[m].averef; @@ -1510,5 +1513,5 @@ // for historical reasons, Mflat is defined with the wrong sign - float delta = MsysPSF - MrelPSF - Mcal - Mmos - Mgrp; + float delta = MsysPSF - MrelPSF - Mcal - Mmos - Mgrp - Mflat; fprintf (fout, "%f %f : %f %f : %f %f %f : %f %f\n", catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysPSF, MrelPSF, Mcal, Mmos, Mgrp, catalog[c].measureT[m].airmass, delta); Index: trunk/Ohana/src/relphot/src/WarpImageMaps.c =================================================================== --- trunk/Ohana/src/relphot/src/WarpImageMaps.c (revision 41601) +++ (revision ) @@ -1,279 +1,0 @@ -# include "relphot.h" - -// we have a list of warp images, each with obstime & photcode. I need to make a list of the obstime groups - -# define D_NOBSTIMES 10000 - -typedef struct { - int obstime; - unsigned short photcode; - int Nwarps; - int NWARPS; - off_t *warps; - double *Rmin; - double *Rmax; - double *Dmin; - double *Dmax; - int *onBoundary; -} WarpGroup; - -WarpGroup *warpgroup = NULL; -int Nwarpgroup = 0; -int NWARPGROUP = 0; - -myIndexType *warpObstimeIndex = NULL; - -int FindWarpGroups (void) { - - off_t Nimage; - Image *image = getimages (&Nimage, NULL); - - int *obstimes = NULL; - int Nobstimes = 0; - int NOBSTIMES = D_NOBSTIMES; - ALLOCATE (obstimes, int, NOBSTIMES); - - // first generate a list of obstime values - for (off_t i = 0; i < Nimage; i++) { - if (!isGPC1warp(image[i].photcode)) continue; - - obstimes[Nobstimes] = image[i].tzero; - Nobstimes ++; - - CHECK_REALLOCATE (obstimes, int, NOBSTIMES, Nobstimes, D_NOBSTIMES); - } - if (!Nobstimes) { - free (obstimes); - return FALSE; - } - - // sort the list: - isort (obstimes, Nobstimes); - - // find the unique obstimes - int *uniqtimes = NULL; - int *uniqcount = NULL; - int Nuniqtimes = 0; - ALLOCATE (uniqtimes, int, Nobstimes); - ALLOCATE (uniqcount, int, Nobstimes); - - // generate a uniq set of obstimes - for (int i = 0; i < Nobstimes; Nuniqtimes ++) { - uniqtimes[Nuniqtimes] = obstimes[i]; - int Ndup = 0; - int lastValue = uniqtimes[Nuniqtimes]; - while ((i < Nobstimes) && (obstimes[i] == lastValue)) { - i++; - Ndup ++; - uniqcount[Nuniqtimes] = Ndup; - } - } - REALLOCATE (uniqtimes, int, Nuniqtimes); - REALLOCATE (uniqcount, int, Nuniqtimes); - - // now I need to assign all warp images to a warpgroup - // I have 2 options to go from warp obstime to warpgroup: - // a) use a bisection lookup [may be slow] - // b) use an index on obstime [require ~5*3.14*1e7*4 bytes ~ 600MB for index] - - // generate an index on obstime - warpObstimeIndex = myIndexInit (); - warpObstimeIndex->minID = obstimes[0]; - warpObstimeIndex->maxID = obstimes[Nobstimes-1]; - - myIndexSetRange (warpObstimeIndex); - - // generate the warp groups (create the index as we go) - Nwarpgroup = Nuniqtimes; - ALLOCATE (warpgroup, WarpGroup, Nwarpgroup); - - for (int i = 0; i < Nuniqtimes; i++) { - warpgroup[i].obstime = uniqtimes[i]; - warpgroup[i].photcode = 0; // not yet set - warpgroup[i].Nwarps = 0; - warpgroup[i].NWARPS = 100; - ALLOCATE (warpgroup[i].warps, off_t, warpgroup[i].NWARPS); - myIndexSetEntry (warpObstimeIndex, uniqtimes[i], i); - } - - // assign all warps to one of the warpgroups - for (off_t i = 0; i < Nimage; i++) { - if (!isGPC1warp(image[i].photcode)) continue; - - int seq = myIndexGetEntry (warpObstimeIndex, image[i].tzero); - myAssert (warpgroup[seq].obstime == image[i].tzero, "oops"); - - if (!warpgroup[seq].photcode) { - warpgroup[seq].photcode = image[i].photcode; - } else { - myAssert (warpgroup[seq].photcode == image[i].photcode, "oops"); - } - - int N = warpgroup[seq].Nwarps; - warpgroup[seq].warps[N] = i; - warpgroup[seq].Nwarps ++; - - CHECK_REALLOCATE (warpgroup[seq].warps, off_t, warpgroup[seq].NWARPS, warpgroup[seq].Nwarps, 100); - } - - // we now have the warps assigned to groups. now we need to generate the grid data to find the warp assignments - for (int i = 0; i < Nwarpgroup; i++) { - - myAssert (warpgroup[i].Nwarps == uniqcount[i], "failure"); - ALLOCATE (warpgroup[i].Rmin, double, warpgroup[i].Nwarps); - ALLOCATE (warpgroup[i].Rmax, double, warpgroup[i].Nwarps); - ALLOCATE (warpgroup[i].Dmin, double, warpgroup[i].Nwarps); - ALLOCATE (warpgroup[i].Dmax, double, warpgroup[i].Nwarps); - ALLOCATE (warpgroup[i].onBoundary, int, warpgroup[i].Nwarps); - - for (int j = 0; j < warpgroup[i].Nwarps; j++) { - - off_t N = warpgroup[i].warps[j]; - - int Nx = image[N].NX; - int Ny = image[N].NY; - - warpgroup[i].onBoundary[j] = FALSE; - - double R, D; - XY_to_RD (&R, &D, 0.0, 0.0, &image[N].coords); - R = ohana_normalize_angle_to_midpoint(R, 180.0); - warpgroup[i].Rmin[j] = R; - warpgroup[i].Rmax[j] = R; - warpgroup[i].Dmin[j] = D; - warpgroup[i].Dmax[j] = D; - - // XXX need to worry about 0,360 boundary - XY_to_RD (&R, &D, Nx, 0.0, &image[N].coords); - R = ohana_normalize_angle_to_midpoint(R, 180.0); - warpgroup[i].Rmin[j] = MIN (warpgroup[i].Rmin[j], R); - warpgroup[i].Rmax[j] = MAX (warpgroup[i].Rmax[j], R); - warpgroup[i].Dmin[j] = MIN (warpgroup[i].Dmin[j], D); - warpgroup[i].Dmax[j] = MAX (warpgroup[i].Dmax[j], D); - - // XXX need to worry about 0,360 boundary - XY_to_RD (&R, &D, 0.0, Ny, &image[N].coords); - R = ohana_normalize_angle_to_midpoint(R, 180.0); - warpgroup[i].Rmin[j] = MIN (warpgroup[i].Rmin[j], R); - warpgroup[i].Rmax[j] = MAX (warpgroup[i].Rmax[j], R); - warpgroup[i].Dmin[j] = MIN (warpgroup[i].Dmin[j], D); - warpgroup[i].Dmax[j] = MAX (warpgroup[i].Dmax[j], D); - - // XXX need to worry about 0,360 boundary - XY_to_RD (&R, &D, Nx, Ny, &image[N].coords); - R = ohana_normalize_angle_to_midpoint(R, 180.0); - warpgroup[i].Rmin[j] = MIN (warpgroup[i].Rmin[j], R); - warpgroup[i].Rmax[j] = MAX (warpgroup[i].Rmax[j], R); - warpgroup[i].Dmin[j] = MIN (warpgroup[i].Dmin[j], D); - warpgroup[i].Dmax[j] = MAX (warpgroup[i].Dmax[j], D); - - // bump Rmin,Rmax,Dmin,Dmax 10 arcsec worth of padding - double dR = 10.0/3600.0 / cos (RAD_DEG*warpgroup[i].Dmin[j]); - double dD = 10.0/3600.0; - - // at north pole, force test of all nearby skycells - if (warpgroup[i].Dmax[j] > 89.8) { - warpgroup[i].Rmin[j] = 0.0; - warpgroup[i].Rmax[j] = 360.0; - warpgroup[i].Dmax[j] = 90.0; - warpgroup[i].Dmin[j] -= dD; - continue; - } - - if (warpgroup[i].Rmax[j] - warpgroup[i].Rmin[j] > 270.0) { - // Rmin and Rmax are in the range 0 - 360. For images at the 0,360 boundary, - // "Rmax" is the lower edge, and "Rmin" is the upper edge. we need to flip them - // and then break the 0-360 range: - double tmp = warpgroup[i].Rmin[j]; - warpgroup[i].Rmin[j] = warpgroup[i].Rmax[j] - 360.0; - warpgroup[i].Rmax[j] = tmp; - warpgroup[i].onBoundary[j] = TRUE; - } - - warpgroup[i].Rmin[j] -= dR; - warpgroup[i].Rmax[j] += dR; - warpgroup[i].Dmin[j] -= dD; - warpgroup[i].Dmax[j] += dD; - } - } - - free (uniqtimes); - free (uniqcount); - free (obstimes); - return TRUE; -} - -int GetWarpSeq (Image *image, int obstime, unsigned short photcode, double Rave, double Dave, float X, float Y) { - - if (!warpObstimeIndex) return -1; - - int seq = myIndexGetEntry (warpObstimeIndex, obstime); - myAssert (seq > -1, "ooops"); - myAssert (warpgroup[seq].photcode == photcode, "oops"); - - double dPosMin = NAN; - int nPosMin = -1; - - // we now have the warp group, but which is the correct warp? - for (int i = 0; i < warpgroup[seq].Nwarps; i++) { - if (warpgroup[seq].onBoundary[i]) { - int inRange1 = (Rave >= warpgroup[seq].Rmin[i]) && (Rave <= warpgroup[seq].Rmax[i]); - int inRange2 = (Rave >= warpgroup[seq].Rmin[i] + 360.0) && (Rave <= warpgroup[seq].Rmax[i] + 360.0); - if (!inRange1 && !inRange2) continue; - } else { - if (Rave < warpgroup[seq].Rmin[i]) continue; - if (Rave > warpgroup[seq].Rmax[i]) continue; - } - if (Dave < warpgroup[seq].Dmin[i]) continue; - if (Dave > warpgroup[seq].Dmax[i]) continue; - - // this is a possible image: check the coordinates: - int N = warpgroup[seq].warps[i]; - - // project Rave,Dave to image pixels - double Xtst, Ytst; - RD_to_XY (&Xtst, &Ytst, Rave, Dave, &image[N].coords); - - // find the pixel offset - double dX = (Xtst - X); - double dY = (Ytst - Y); - - // skip detections which are within a small distance from the expected location - double dPos = hypot(dX,dY); - if (dPos < 20.0) return N; // 20.0 pixels = 5.0 arcsec - - // check for multiple possible matches?? - - // record the min pos so we can see how things fail - if (isnan(dPosMin)) { - dPosMin = dPos; - nPosMin = i; - } else { - if (dPos < dPosMin) { - dPosMin = dPos; - nPosMin = i; - } - } - } - - // if we did not find any matched, give up - if (nPosMin < 0) return -1; - - // return closest image - return warpgroup[seq].warps[nPosMin]; -} - -void FreeWarpGroups (void) { - - for (int i = 0; i < Nwarpgroup; i++) { - free (warpgroup[i].warps); - free (warpgroup[i].Rmin); - free (warpgroup[i].Rmax); - free (warpgroup[i].Dmin); - free (warpgroup[i].Dmax); - free (warpgroup[i].onBoundary); - } - - free (warpgroup); - myIndexFree (warpObstimeIndex); -} Index: trunk/Ohana/src/relphot/src/ZeroPointModes.c =================================================================== --- trunk/Ohana/src/relphot/src/ZeroPointModes.c (revision 41601) +++ trunk/Ohana/src/relphot/src/ZeroPointModes.c (revision 41647) @@ -22,5 +22,5 @@ } -void SetZeroPointModes (void) { +void SetZeroPointModes (Catalog *catalog, int Ncatalog) { if (TGROUP_ZEROPT && MOSAIC_ZEROPT) { @@ -49,5 +49,5 @@ clean_mosaics(); // do this on every pass } - if ((CurrentLoop >= 18) && (CurrentLoop <= 999)) { + if ((CurrentLoop >= 18) && (CurrentLoop < 21)) { // after iterating a few times on the TGroup & Mosaic zero points: // * identify the good / bad Mosaics @@ -59,4 +59,18 @@ clean_mosaics(); // do this on every pass clean_images(); // do this on every pass + clean_stars(catalog, Ncatalog); // do this on every pass + } + if ((CurrentLoop >= 21) && (CurrentLoop <= 999)) { + // after iterating a few times on the TGroup & Mosaic zero points: + // * identify the good / bad Mosaics + // * fit the chips from the bad Mosaics + TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_GOOD_NIGHT; // stop fitting the bad nights + MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_BAD_NIGHT_GOOD_MOSAIC; // only fit good mosaics on bad nights + IMAGE_ZPT_MODE = IMAGE_ZPT_MODE_BAD_NIGHT_BAD_MOSAIC; // only fit bad mosaics on bad nights + GRID_ZPT_MODE = GRID_ZPT_MODE_ALL; // only fit bad mosaics on bad nights + clean_tgroups(); // do this on every pass + clean_mosaics(); // do this on every pass + clean_images(); // do this on every pass + clean_stars(catalog, Ncatalog); // do this on every pass } } @@ -68,5 +82,5 @@ MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_NONE; } - if ((CurrentLoop > 4) && (CurrentLoop <= 999)) { + if ((CurrentLoop > 4) && (CurrentLoop <= 8)) { // after iterating a few times on the TGroup zero points: // * identify the photometric nights @@ -77,13 +91,23 @@ clean_tgroups(); // do this on every pass or just sometimes? } + if ((CurrentLoop > 8) && (CurrentLoop <= 999)) { + // after iterating a few times on the TGroup zero points: + // * identify the photometric nights + // * fit the images from the non-photometric nights + TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_GOOD_NIGHT; // stop fitting the bad nights + MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_NONE; + IMAGE_ZPT_MODE = IMAGE_ZPT_MODE_BAD_NIGHT; + GRID_ZPT_MODE = GRID_ZPT_MODE_ALL; // only fit bad mosaics on bad nights + clean_tgroups(); // do this on every pass or just sometimes? + } } if (!TGROUP_ZEROPT && MOSAIC_ZEROPT) { - if (CurrentLoop <= 4) { + if (CurrentLoop < 6) { TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_NONE; MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_ALL; IMAGE_ZPT_MODE = IMAGE_ZPT_MODE_NONE; } - if ((CurrentLoop > 4) && (CurrentLoop <= 999)) { + if ((CurrentLoop >= 6) && (CurrentLoop < 12)) { // after iterating a few times on the Mosaic zero points: // * identify the good / bad mosaics @@ -94,4 +118,27 @@ clean_mosaics(); // do this on every pass or just sometimes? } + if ((CurrentLoop >= 12) && (CurrentLoop < 18)) { + // after iterating a few times on the Mosaic zero points: + // * identify the good / bad mosaics + // * fit the images from the bad mosaics + TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_NONE; // stop fitting the bad nights + MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_GOOD_MOSAIC; + IMAGE_ZPT_MODE = IMAGE_ZPT_MODE_BAD_MOSAIC; + clean_mosaics(); // do this on every pass or just sometimes? + clean_images(); // do this on every pass + clean_stars(catalog, Ncatalog); // do this on every pass + } + if ((CurrentLoop >= 18) && (CurrentLoop < 999)) { + // after iterating a few times on the Mosaic zero points: + // * identify the good / bad mosaics + // * fit the images from the bad mosaics + TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_NONE; // stop fitting the bad nights + MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_GOOD_MOSAIC; + IMAGE_ZPT_MODE = IMAGE_ZPT_MODE_BAD_MOSAIC; + GRID_ZPT_MODE = GRID_ZPT_MODE_ALL; // only fit bad mosaics on bad nights + clean_mosaics(); // do this on every pass or just sometimes? + clean_images(); // do this on every pass + clean_stars(catalog, Ncatalog); // do this on every pass + } } Index: trunk/Ohana/src/relphot/src/args.c =================================================================== --- trunk/Ohana/src/relphot/src/args.c (revision 41601) +++ trunk/Ohana/src/relphot/src/args.c (revision 41647) @@ -198,11 +198,4 @@ } - NGRID = 8; - if ((N = get_argument (argc, argv, "-ngrid"))) { - remove_argument (N, &argc, argv); - NGRID = atof (argv[N]); - remove_argument (N, &argc, argv); - } - TEST_IMAGE1 = -1; if ((N = get_argument (argc, argv, "-test-image1"))) { @@ -227,5 +220,4 @@ RESET = TRUE; } - RESET_ZEROPTS = FALSE; if ((N = get_argument (argc, argv, "-reset-zpts"))) { @@ -233,4 +225,9 @@ RESET_ZEROPTS = TRUE; } + RESET_FLATCORR = FALSE; + if ((N = get_argument (argc, argv, "-reset-flat"))) { + remove_argument (N, &argc, argv); + RESET_FLATCORR = TRUE; + } REPAIR_WARPS = FALSE; @@ -244,4 +241,9 @@ remove_argument (N, &argc, argv); PRESERVE_PS1 = TRUE; + } + REQUIRE_PSFFIT = FALSE; + if ((N = get_argument (argc, argv, "-require-psffit"))) { + remove_argument (N, &argc, argv); + REQUIRE_PSFFIT = TRUE; } @@ -387,9 +389,42 @@ } - // USE_GRID is not valid for all cases, probably should be its own mode... - USE_GRID = FALSE; + // GRID_ZEROPT is not valid for all cases, probably should be its own mode... + GRID_ZEROPT = FALSE; + GRID_MEANFILE = NULL; // this is set by GridCorrectionSave() + GRID_ZPT_MODE = GRID_ZPT_MODE_NONE; // start with grid off if ((N = get_argument (argc, argv, "-grid"))) { remove_argument (N, &argc, argv); - USE_GRID = TRUE; + GRID_ZEROPT = TRUE; + } + if ((N = get_argument (argc, argv, "-grid-meanfile"))) { + GRID_ZEROPT = TRUE; + remove_argument (N, &argc, argv); + GRID_MEANFILE = strcreate (argv[N]); // used by relphot -averages and relphot -apply-offsets + remove_argument (N, &argc, argv); + } + + GRID_BIN_GPC1 = 16; + if ((N = get_argument (argc, argv, "-grid-bin-gpc1"))) { + remove_argument (N, &argc, argv); + GRID_BIN_GPC1 = atof (argv[N]); + remove_argument (N, &argc, argv); + } + GRID_BIN_GPC2 = 16; + if ((N = get_argument (argc, argv, "-grid-bin-gpc2"))) { + remove_argument (N, &argc, argv); + GRID_BIN_GPC2 = atof (argv[N]); + remove_argument (N, &argc, argv); + } + GRID_BIN_HSC = 8; + if ((N = get_argument (argc, argv, "-grid-bin-hsc"))) { + remove_argument (N, &argc, argv); + GRID_BIN_HSC = atof (argv[N]); + remove_argument (N, &argc, argv); + } + GRID_BIN_CFH = 8; + if ((N = get_argument (argc, argv, "-grid-bin-cfh"))) { + remove_argument (N, &argc, argv); + GRID_BIN_CFH = atof (argv[N]); + remove_argument (N, &argc, argv); } @@ -555,4 +590,6 @@ case SYNTH_PHOT: case UPDATE_AVERAGES: + // note: initialize sets PhotcodeList to the full set of average photcodes + // if the mode is UPDATE_AVERAGES if (argc != 1) relphot_usage(); break; @@ -608,4 +645,5 @@ FREE (REGION_FILE); FREE (SYNTH_ZERO_POINTS); + FREE (GRID_MEANFILE); FREE (IMAGE_TABLE); FREE (CATDIR); @@ -763,9 +801,13 @@ RESET = TRUE; } - RESET_ZEROPTS = FALSE; if ((N = get_argument (argc, argv, "-reset-zpts"))) { remove_argument (N, &argc, argv); RESET_ZEROPTS = TRUE; + } + RESET_FLATCORR = FALSE; + if ((N = get_argument (argc, argv, "-reset-flat"))) { + remove_argument (N, &argc, argv); + RESET_FLATCORR = TRUE; } @@ -864,4 +906,15 @@ } + // load the GridCorrection file + GRID_ZEROPT = FALSE; + GRID_MEANFILE = NULL; + if ((N = get_argument (argc, argv, "-grid"))) { + GRID_ZEROPT = TRUE; + remove_argument (N, &argc, argv); + GRID_MEANFILE = strcreate (argv[N]); + remove_argument (N, &argc, argv); + GRID_ZPT_MODE = GRID_ZPT_MODE_ALL; // since we are applying the grid, need to active this + } + ImagSelect = FALSE; if ((N = get_argument (argc, argv, "-instmag"))) { @@ -922,4 +975,5 @@ FREE (SYNTH_ZERO_POINTS); + FREE (GRID_MEANFILE); FREE (BOUNDARY_TREE); FREE (UPDATE_CATFORMAT); Index: trunk/Ohana/src/relphot/src/assign_images.c =================================================================== --- trunk/Ohana/src/relphot/src/assign_images.c (revision 41601) +++ trunk/Ohana/src/relphot/src/assign_images.c (revision 41647) @@ -36,6 +36,7 @@ if (MOSAIC_ZEROPT) { makeMosaics (image, Nimage, FALSE); - - // center coords and Mcal, dMcal, Mchisq for the mosaics + MARKTIME("set mosaic coordinates and Mcal values: %f sec\n", dtime); + + // center coords and zero Mcal, dMcal, Mchisq for the mosaics setMosaicCenters (image, Nimage); MARKTIME("set mosaic coordinates and Mcal values: %f sec\n", dtime); @@ -124,9 +125,9 @@ Mosaic *mosaic = getMosaicForImage (j); if (mosaic) { - Rc = mosaic->coords.crval1; - Dc = mosaic->coords.crval2; - // NOTE : have defined mosaic Rc,Dc to choose the side of 0,360 on which most of the - // chips are located. but, for host assignment, we rationalize to 0.0 - 360.0 - Rc = ohana_normalize_angle_to_midpoint (Rc, 180.0); + Rc = mosaic->coords.crval1; + Dc = mosaic->coords.crval2; + // NOTE : have defined mosaic Rc,Dc to choose the side of 0,360 on which most of the + // chips are located. but, for host assignment, we rationalize to 0.0 - 360.0 + Rc = ohana_normalize_angle_to_midpoint (Rc, 180.0); } } Index: trunk/Ohana/src/relphot/src/bcatalog.c =================================================================== --- trunk/Ohana/src/relphot/src/bcatalog.c (revision 41601) +++ trunk/Ohana/src/relphot/src/bcatalog.c (revision 41647) @@ -35,17 +35,4 @@ ALLOCATE (Nvalid, int, Nsecfilt); - // flags used by the photometry analysis (excluding UBERCAL) - unsigned int PHOTOM_FLAGS = - ID_MEAS_NOCAL | // detection ignored for this analysis (photcode, time range) - ID_MEAS_POOR_PHOTOM | // detection is photometry outlier - ID_MEAS_SKIP_PHOTOM | // detection was ignored for photometry measurement - ID_MEAS_AREA | // detetion was outside acceptable area of device - ID_MEAS_SYNTH_MAG | // magnitude is synthetic - ID_MEAS_STACK_PRIMARY | // this stack measurement is in the primary skycell - ID_MEAS_STACK_PHOT_SRC | // this measurement supplied the stack photometry - ID_MEAS_PHOTOM_PSF | // this measurement is used for the mean psf mag - ID_MEAS_PHOTOM_APER | // this measurement is used for the mean ap mag - ID_MEAS_PHOTOM_KRON ; // this measurement is used for the mean kron mag - /* exclude stars not in range or with too few measurements */ for (i = 0; i < catalog[0].Naverage; i++) { @@ -61,16 +48,6 @@ } - // reset the calculated average magnitudes (does not affect ubercal-tied measurements or images) - if (RESET) { - int Ns; - - for (Ns = 0; Ns < Nphotcodes; Ns++) { - - int thisCode = photcodes[Ns][0].code; - int Nsec = GetPhotcodeNsec(thisCode); - - dvo_secfilt_init (&subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec], SECFILT_RESET_ALL); - } - } + // reset the calculated average magnitudes for active photcode only + ResetAverageActivePhotcodes (&subcatalog[0].secfilt[Nsecfilt*Naverage]); Nm = 0; @@ -111,6 +88,6 @@ if (isnan(catalog[0].measure[offset].Mkron)) { Nnan ++; continue; } - // require 0x01 in photFlags (fitted with a PSF) - if ((catalog[0].measure[offset].photFlags & 0x01) == 0) { Nbad ++; continue; } + // require 0x01 in photFlags (fitted with a PSF) -- add to the test data + if (REQUIRE_PSFFIT && (catalog[0].measure[offset].photFlags & 0x01) == 0) { Nbad ++; continue; } // very loose cut on PSF - Kron @@ -156,21 +133,6 @@ subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_SKIP_PHOTOM; subcatalog[0].measureT[Nmeasure].averef = Naverage; - if (RESET) { - // only reset Mcal for measures with a matching image - // do not reset Mcal for ubercal images unless explicitly requested - if (subcatalog[0].measureT[Nmeasure].dbFlags & ID_MEAS_PHOTOM_UBERCAL) { - if (!KEEP_UBERCAL) { - subcatalog[0].measureT[Nmeasure].McalPSF = 0.0; - subcatalog[0].measureT[Nmeasure].McalAPER = 0.0; - subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL; - } - } else { - if (getImageEntry (Nmeasure, Ncat) >= 0) { - subcatalog[0].measureT[Nmeasure].McalPSF = 0.0; - subcatalog[0].measureT[Nmeasure].McalAPER = 0.0; - } - } - subcatalog[0].measureT[Nmeasure].dbFlags &= ~PHOTOM_FLAGS; - } + ResetMeasureZeroPoints (&subcatalog[0].measureT[Nmeasure], Nmeasure, Ncat); + Nmeasure ++; Nm ++; Index: trunk/Ohana/src/relphot/src/extra.c =================================================================== --- trunk/Ohana/src/relphot/src/extra.c (revision 41601) +++ trunk/Ohana/src/relphot/src/extra.c (revision 41647) @@ -1,3 +1,4 @@ # include "relphot.h" +# define DISABLE_CORNER_GPC1 0 int isMosaicChip (int photcode) { @@ -14,5 +15,5 @@ int whichGPC1filter (int photcode) { -# if (1) +# if (DISABLE_CORNER_GPC1) // disable the corner chips: if (((photcode > 10000) && (photcode < 10077)) || (photcode == 4100)) return PS1_g; // g-band @@ -37,5 +38,5 @@ int isGPC1chip (int photcode) { -# if (1) +# if (DISABLE_CORNER_GPC1) if (((photcode > 10000) && (photcode < 10077)) || (photcode == 4100)) return TRUE; // g-band if (((photcode > 10100) && (photcode < 10177)) || (photcode == 4200)) return TRUE; // r-band Index: trunk/Ohana/src/relphot/src/global_stats.c =================================================================== --- trunk/Ohana/src/relphot/src/global_stats.c (revision 41601) +++ trunk/Ohana/src/relphot/src/global_stats.c (revision 41647) @@ -1,5 +1,5 @@ # include "relphot.h" -void global_stats (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr, int nloop) { +void global_stats (Catalog *catalog, int Ncatalog, int nloop) { StatType stN, stX, stS, imN, imX, imM, imD, msM, msX, msN, msD, tgM, tgX, tgD; @@ -33,5 +33,5 @@ int seccode = photcodes[Ns][0].code; - stN = statsStarN (catalog, Ncatalog, Nsec, seccode, flatcorr); + stN = statsStarN (catalog, Ncatalog, Nsec, seccode); stX = statsStarX (catalog, Ncatalog, Nsec); stS = statsStarS (catalog, Ncatalog, Nsec); @@ -43,10 +43,10 @@ } - imN = statsImageN (catalog, flatcorr); + imN = statsImageN (catalog); imX = statsImageX (catalog); imM = statsImageM (catalog); imD = statsImagedM (catalog); - msN = statsMosaicN (catalog, flatcorr); + msN = statsMosaicN (catalog); msM = statsMosaicM (catalog); msD = statsMosaicdM (catalog); @@ -58,10 +58,10 @@ fprintf (stderr, "meas / image: %7.0f %7.0f %7.0f %7.0f %7.0f %6d\n", imN.median, imN.mean, imN.sigma, imN.min, imN.max, imN.Nmeas); - fprintf (stderr, "Mcal image: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n", imM.median, imM.mean, imM.sigma, imM.min, imM.max, imM.Nmeas); + fprintf (stderr, "Mcal image: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n", imM.median, imM.mean, imM.sigma, imM.min, imM.max, imM.Nmeas); fprintf (stderr, "dMcal image: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n", imD.median, imD.mean, imD.sigma, imD.min, imD.max, imD.Nmeas); fprintf (stderr, "chisq image: %7.1f %7.1f %7.1f %7.1f %7.1f %6d\n", imX.median, imX.mean, imX.sigma, imX.min, imX.max, imX.Nmeas); fprintf (stderr, "meas / mosaic: %7.0f %7.0f %7.0f %7.0f %7.0f %6d\n", msN.median, msN.mean, msN.sigma, msN.min, msN.max, msN.Nmeas); - fprintf (stderr, "Mcal mosaic: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n", msM.median, msM.mean, msM.sigma, msM.min, msM.max, msM.Nmeas); + fprintf (stderr, "Mcal mosaic: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n", msM.median, msM.mean, msM.sigma, msM.min, msM.max, msM.Nmeas); fprintf (stderr, "dMcal mosaic: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n", msD.median, msD.mean, msD.sigma, msD.min, msD.max, msD.Nmeas); fprintf (stderr, "chisq mosaic: %7.1f %7.1f %7.1f %7.1f %7.1f %6d\n", msX.median, msX.mean, msX.sigma, msX.min, msX.max, msX.Nmeas); Index: trunk/Ohana/src/relphot/src/initialize.c =================================================================== --- trunk/Ohana/src/relphot/src/initialize.c (revision 41601) +++ trunk/Ohana/src/relphot/src/initialize.c (revision 41647) @@ -27,13 +27,4 @@ } - if (USE_GRID && (Nphotcodes > 1)) { - fprintf (stderr, "grid correction analysis currently can only operate on a single photcode\n"); - exit (1); - } - - IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP; - STAR_BAD = ID_OBJ_POOR | ID_OBJ_FEW; - MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA; - if (SHOW_PARAMS) { int Ns; @@ -67,11 +58,14 @@ fprintf (stderr, "VERBOSE: %d, PLOTSTUFF: %d\n", VERBOSE, PLOTSTUFF); - fprintf (stderr, "GRID_X: %d, GRID_Y: %d, BINNING: %d == Nmx: %d, Nmy: %d\n", - RELPHOT_GRID_X, - RELPHOT_GRID_Y, - RELPHOT_GRID_BINNING, - (int)(RELPHOT_GRID_X/RELPHOT_GRID_BINNING), (int)(RELPHOT_GRID_Y/RELPHOT_GRID_BINNING)); fprintf (stderr, "STAR_SCATTER %lf\n", STAR_SCATTER); + fprintf (stderr, "STAR_CHISQ %lf\n", STAR_CHISQ); + + fprintf (stderr, "MOSAIC_SCATTER %lf\n", MOSAIC_SCATTER); + fprintf (stderr, "MOSAIC_CHISQ %lf\n", MOSAIC_CHISQ); + + fprintf (stderr, "NIGHT_SCATTER %lf\n", NIGHT_SCATTER); + fprintf (stderr, "NIGHT_CHISQ %lf\n", NIGHT_CHISQ); + fprintf (stderr, "IMAGE_SCATTER %lf\n", IMAGE_SCATTER); fprintf (stderr, "IMAGE_OFFSET %lf\n", IMAGE_OFFSET); Index: trunk/Ohana/src/relphot/src/launch_region_hosts.c =================================================================== --- trunk/Ohana/src/relphot/src/launch_region_hosts.c (revision 41601) +++ trunk/Ohana/src/relphot/src/launch_region_hosts.c (revision 41647) @@ -62,4 +62,5 @@ if (RESET) strextend (&command, "-reset"); if (RESET_ZEROPTS) strextend (&command, "-reset-zpts"); + if (RESET_FLATCORR) strextend (&command, "-reset-flat"); if (!KEEP_UBERCAL) strextend (&command, "-reset-ubercal"); if (DophotSelect) strextend (&command, "-dophot %d", DophotValue); Index: trunk/Ohana/src/relphot/src/load_catalogs.c =================================================================== --- trunk/Ohana/src/relphot/src/load_catalogs.c (revision 41601) +++ trunk/Ohana/src/relphot/src/load_catalogs.c (revision 41647) @@ -165,4 +165,5 @@ if (RESET) { strextend (&command, "-reset"); } if (RESET_ZEROPTS) { strextend (&command, "-reset-zpts"); } + if (RESET_FLATCORR) { strextend (&command, "-reset-flat"); } if (!KEEP_UBERCAL) { strextend (&command, "-reset-ubercal"); } if (DophotSelect) { strextend (&command, "-dophot %d", DophotValue); } Index: trunk/Ohana/src/relphot/src/load_images.c =================================================================== --- trunk/Ohana/src/relphot/src/load_images.c (revision 41601) +++ trunk/Ohana/src/relphot/src/load_images.c (revision 41647) @@ -57,21 +57,5 @@ // reset image values as needed. always allow 'few' images to succeed, if possible (new // images / detections may have been added - for (off_t i = 0; i < Nsubset; i++) { - // reset these two regardless (we will re-determine) - subset[i].flags &= ~(ID_IMAGE_PHOTOM_NOCAL | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP | ID_IMAGE_PHOTOM_FEW); - if (RESET) { - if (RESET_ZEROPTS) { - if (!KEEP_UBERCAL || !(subset[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) { - subset[i].McalPSF = 0.0; - subset[i].McalAPER = 0.0; - subset[i].dMcal = NAN; - subset[i].flags &= ~ID_IMAGE_PHOTOM_UBERCAL; - } - } - subset[i].flags &= ~(ID_IMAGE_MOSAIC_POOR | ID_IMAGE_NIGHT_POOR); - subset[i].flags &= ~(ID_IMAGE_IMAGE_PHOTCAL | ID_IMAGE_MOSAIC_PHOTCAL | ID_IMAGE_TGROUP_PHOTCAL); - subset[i].ubercalDist = 1000; - } - } + ResetImages (subset, Nsubset); // XXX consider allowing relphot to skip the AstroMap load (do we need precise astrometry here?) @@ -115,4 +99,11 @@ return TRUE; } + +// This function re-creates the image subset array using the array of LineNumbers saved +// initially by load_images. This function is only called by relphot_images.c after the +// loops have been completed (if any). I am not certain this function is doing anything +// at all : the values in the image subset are copied back to the master table in +// save_images_updates() so the values copied from the master to the subset should already +// be there. int reload_images (FITS_DB *db) { Index: trunk/Ohana/src/relphot/src/plot_scatter.c =================================================================== --- trunk/Ohana/src/relphot/src/plot_scatter.c (revision 41601) +++ trunk/Ohana/src/relphot/src/plot_scatter.c (revision 41647) @@ -1,5 +1,5 @@ # include "relphot.h" -void plot_scatter (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) { +void plot_scatter (Catalog *catalog, int Ncatalog) { off_t i, j, k, m, N, Ntot; @@ -28,32 +28,32 @@ N = 0; for (i = 0; i < Ncatalog; i++) { - for (j = 0; j < catalog[i].Naverage; j++) { + for (j = 0; j < catalog[i].Naverage; j++) { - /* calculate the average value for a single star */ - if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; - m = catalog[i].averageT[j].measureOffset; + /* calculate the average value for a single star */ + // if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; + m = catalog[i].averageT[j].measureOffset; - for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) { - // skip measurements that do not match the current photcode - int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode); - if (ecode != thisCode) { continue; } + for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) { + // skip measurements that do not match the current photcode + int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode); + if (ecode != thisCode) { continue; } - if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; - Mcal = getMcal (m, i, MAG_CLASS_PSF); - if (isnan(Mcal)) continue; - Mmos = getMmos (m, i); - if (isnan(Mmos)) continue; - Mgrid = getMgrid (m, i); - if (isnan(Mgrid)) continue; + if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; + Mcal = getMcal (m, i, MAG_CLASS_PSF); + if (isnan(Mcal)) continue; + Mmos = getMmos (m, i); + if (isnan(Mmos)) continue; + Mgrid = getMgridTiny (&catalog[i].measureT[m]); + if (isnan(Mgrid)) continue; - Mrel = catalog[i].secfilt[Nsecfilt*j+Nsec].MpsfChp; - if (isnan(Mrel)) continue; + Mrel = catalog[i].secfilt[Nsecfilt*j+Nsec].MpsfChp; + if (isnan(Mrel)) continue; - xlist[N] = Mrel; - ylist[N] = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt], MAG_CLASS_PSF) - Mcal - Mmos - Mgrid - Mrel; - ilist[N] = PhotInstTiny (&catalog[i].measureT[m], MAG_CLASS_PSF); - N++; - } + xlist[N] = Mrel; + ylist[N] = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt], MAG_CLASS_PSF) - Mcal - Mmos - Mgrid - Mrel; + ilist[N] = PhotInstTiny (&catalog[i].measureT[m], MAG_CLASS_PSF); + N++; } + } } Index: trunk/Ohana/src/relphot/src/reload_catalogs.c =================================================================== --- trunk/Ohana/src/relphot/src/reload_catalogs.c (revision 41601) +++ trunk/Ohana/src/relphot/src/reload_catalogs.c (revision 41647) @@ -9,5 +9,5 @@ gettimeofday (&start, (void *) NULL); -void reload_catalogs (SkyList *skylist, FlatCorrectionTable *flatcorr, int hostID, char *hostpath) { +void reload_catalogs (SkyList *skylist, int hostID, char *hostpath) { int i; @@ -17,13 +17,6 @@ struct timeval start, stop; - double dtime = 0.0; - double time1 = 0.0; - double time2 = 0.0; - double time3a = 0.0; - double time3b = 0.0; - double time4 = 0.0; - double time5 = 0.0; - double time6 = 0.0; - double time7 = 0.0; + double dtime = 0.0, time1 = 0.0, time2 = 0.0, time3 = 0.0; + double time4 = 0.0, time5 = 0.0, time6 = 0.0, time7 = 0.0, time8 = 0.0; // XXX need to decide how to determine PARALLEL mode... @@ -81,23 +74,17 @@ populate_tiny_values(&catalog, DVO_TV_MEASURE | DVO_TV_AVERAGE); - TIMESTAMP(time3a); + TIMESTAMP(time3); // XXX need to worry about the image subset data initImageBins (&catalog, 1, FALSE); - initMosaicBins (&catalog, 1, FALSE); - initTGroupBins (&catalog, 1); - - initGridBins (&catalog, 1); - TIMESTAMP(time3b); + TIMESTAMP(time4); findImages (&catalog, 1, FALSE); - findMosaics (&catalog, 1, FALSE); - findTGroups (&catalog, 1); - - TIMESTAMP(time4); + TIMESTAMP(time5); initMrel (&catalog, 1); - setMrelFinal (&catalog, flatcorr, FALSE); - TIMESTAMP(time5); + + setMrelFinal (&catalog, FALSE); + TIMESTAMP(time6); // modify the output format as desired (ignore current format on disk) @@ -119,13 +106,8 @@ free_tiny_values(&catalog); dvo_catalog_free (&catalog); - TIMESTAMP(time6); - - // XXX freeTgroupBins here: + TIMESTAMP(time7); freeImageBins (1, FALSE); - freeMosaicBins (1, FALSE); - freeTGroupBins (1); - freeGridBins (1); - TIMESTAMP(time7); + TIMESTAMP(time8); if (hostID) { @@ -134,12 +116,12 @@ } - fprintf (stderr, "time step 1 %10.3f sec : find catalog\n", time1); - fprintf (stderr, "time step 2 %10.3f sec : load catalog\n", time2); - fprintf (stderr, "time step 3 %10.3f sec : make tiny values\n", time3a); - fprintf (stderr, "time step 4 %10.3f sec : init imbins\n", time3b); - fprintf (stderr, "time step 5 %10.3f sec : find images\n", time4); - fprintf (stderr, "time step 6 %10.3f sec : set Mrel\n", time5); - fprintf (stderr, "time step 7 %10.3f sec : save catalog\n", time6); - fprintf (stderr, "time step 8 %10.3f sec : free catalog\n", time7); + fprintf (stderr, "time step 1 %10.3f sec : find catalog\n", time1); + fprintf (stderr, "time step 2 %10.3f sec : load catalog\n", time2); + fprintf (stderr, "time step 3 %10.3f sec : make tiny values\n", time3); + fprintf (stderr, "time step 4 %10.3f sec : init imbins\n", time4); + fprintf (stderr, "time step 5 %10.3f sec : find images\n", time5); + fprintf (stderr, "time step 6 %10.3f sec : set Mrel\n", time6); + fprintf (stderr, "time step 7 %10.3f sec : save catalog\n", time7); + fprintf (stderr, "time step 8 %10.3f sec : free catalog\n", time8); } @@ -231,32 +213,28 @@ char *command = NULL; strextend (&command, "relphot_client %s -update-catalogs %s", PhotcodeList, imageFile); - strextend (&command, "-hostID %d -D CATDIR %s -hostdir %s", group->hosts[i][0].hostID, CATDIR, group->hosts[i][0].pathname); + strextend (&command, "-hostID %d", group->hosts[i][0].hostID); + strextend (&command, "-D CATDIR %s", CATDIR); + strextend (&command, "-hostdir %s", group->hosts[i][0].pathname); strextend (&command, "-region %f %f %f %f", UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); strextend (&command, "-statmode %s", STATMODE); strextend (&command, "-D CAMERA %s", CAMERA); - strextend (&command, "-D STAR_TOOFEW %d -minerror %f", STAR_TOOFEW, MIN_ERROR); - - // options & configs which affect relphot_client -update-catalogs - // VERBOSE, VERBOSE2 - // RESET - // RESET_ZEROPTS - // TimeSelect - // AreaSelect - // STATMODE - // STAR_TOOFEW - // MIN_ERROR - + strextend (&command, "-D STAR_TOOFEW %d", STAR_TOOFEW); + strextend (&command, "-minerror %f", MIN_ERROR); + if (VERBOSE) { strextend (&command, "-v"); } if (VERBOSE2) { strextend (&command, "-vv"); } if (RESET) { strextend (&command, "-reset"); } if (RESET_ZEROPTS) { strextend (&command, "-reset-zpts"); } + if (RESET_FLATCORR) { strextend (&command, "-reset-flat"); } + if (!KEEP_UBERCAL) { strextend (&command, "-reset-ubercal"); } if (PRESERVE_PS1) { strextend (&command, "-preserve-ps1"); } + if (IS_DIFF_DB) { strextend (&command, "-is-diff-db"); } if (UPDATE) { strextend (&command, "-update"); } - if (IS_DIFF_DB) { strextend (&command, "-is-diff-db"); } - if (!KEEP_UBERCAL) { strextend (&command, "-reset-ubercal"); } if (UPDATE_CATFORMAT) { strextend (&command, "-update-catformat %s", UPDATE_CATFORMAT); } if (BOUNDARY_TREE) { strextend (&command, "-boundary-tree %s", BOUNDARY_TREE); } if (SYNTH_ZERO_POINTS) { strextend (&command, "-synthphot-zpts %s", SYNTH_ZERO_POINTS); } + if (GRID_ZEROPT) { strextend (&command, "-grid %s", GRID_MEANFILE); } if (USE_BASIC_CHECK) { strextend (&command, "-basic-image-search"); } + // if (USE_ALL_IMAGES) { strextend (&command, "-use-all-images"); } if (USE_MCAL_PSF_FOR_STACK_APER) { strextend (&command, "-use-mcal-psf-for-stack-aper"); } Index: trunk/Ohana/src/relphot/src/relphot.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot.c (revision 41601) +++ trunk/Ohana/src/relphot/src/relphot.c (revision 41647) @@ -34,5 +34,7 @@ switch (mode) { case UPDATE_IMAGES: + // calculate zero points for images (may group by exposure [mosaic], night [tgroups]; may measure flat-correction) // IF CALLED WITH NLOOP == 0, DOES NOT LOAD bcatalog, just loads image table and generates ImageSubset.dat + // If called with -update, calls reload_catalogs() just like apply_offsets relphot_images (skylist); relphot_free (sky, skylist); @@ -52,5 +54,6 @@ case PARALLEL_REGIONS: - // run image updates in parallel across multiple remote machines + // calculate zero points for images (may group my exposure [mosaic], night [tgroups]; may measure flat-correction) + // equivalent to relphot_images, but run in parallel across multiple remote machines relphot_parallel_regions (sky); relphot_free (sky, skylist); @@ -69,5 +72,7 @@ exit (2); } - reload_catalogs (skylist, NULL, 0, NULL); + // we do not need to load grid corrections here. the value of + // GRID_MEANFILE, if specified, is passed to the clients + reload_catalogs (skylist, 0, NULL); relphot_free (sky, skylist); exit (0); @@ -84,4 +89,9 @@ I would like to merge the functionality of relphot_images (Nloop == 0), relphot_objects, and apply_offsets + + relphot_objects vs reload_objects (apply offsets): + + relphot_objects + ***/ Index: trunk/Ohana/src/relphot/src/relphot_client.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_client.c (revision 41601) +++ trunk/Ohana/src/relphot/src/relphot_client.c (revision 41647) @@ -80,9 +80,11 @@ initImagesSubset (image, NULL, Nimage); - // load the flat-field correction table from CATDIR - FlatCorrectionTable *flatcorr = NULL; - reload_catalogs (skylist, flatcorr, HOST_ID, HOSTDIR); + // load grid corrections here + GridCorrectionLoad (GRID_MEANFILE); + + reload_catalogs (skylist, HOST_ID, HOSTDIR); freeImages ((char *)image); free (image); + freeGridBins (); client_logger_message ("updated catalogs\n"); relphot_client_free (sky, skylist); @@ -90,4 +92,5 @@ } +/* deprecated (absorbed into MODE_UPDATE) case MODE_UPDATE_OBJECTS: { // take the current set of detections and set the mean magnitudes @@ -97,4 +100,5 @@ break; } +*/ case MODE_SYNTH_PHOT: @@ -104,4 +108,5 @@ break; + case MODE_UPDATE_OBJECTS: default: fprintf (stderr, "impossible!"); Index: trunk/Ohana/src/relphot/src/relphot_images.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_images.c (revision 41601) +++ trunk/Ohana/src/relphot/src/relphot_images.c (revision 41647) @@ -38,5 +38,12 @@ MARKTIME("-- load image data: %f sec\n", dtime); - /* load regions and images based on specified sky patch and/or catalog */ + // load regions and images based on specified sky patch and/or catalog + // NOTE: load_images transfers zero points from images to mosaics (initMosaicMcal) + // and to tgroups (initTGroupMcal), resetting the image zero points to 0.0. + + // (if both tgroups and mosaics are used, the mosaics will hold the value and tgroups + // will get values of 0.0. This is not clearly the right way to go, but may not matter + // if we reset everything). + load_images (&db, skylist, &UserPatch, FALSE, USE_ALL_IMAGES); MARKTIME("-- load images: %f sec\n", dtime); @@ -45,14 +52,9 @@ if (!UPDATE) dvo_image_unlock (&db); - // load the flat correction table (if defined) - // we no longer need to do this: the flatcorrection is now saved in measure.Mflat -# if (0) - char flatcorrfile[256]; - sprintf (flatcorrfile, "%s/flatcorr.fits", CATDIR); - FlatCorrectionTable *flatcorr = FlatCorrectionLoad (flatcorrfile, VERBOSE); -# endif - FlatCorrectionTable *flatcorr = NULL; + initGridBins (); // allocates the empty array of corrections (elements are built below) if (CALIBRATE_STACKS_AND_WARPS || (NLOOP > 0)) { + initGrid (); // allocate grid correction entries for existing photcodes + /* load catalog data from region files (hostID is 0 since we are not a client */ catalog = load_catalogs (skylist, &Ncatalog, 0, NULL, NULL); @@ -68,5 +70,4 @@ initTGroupBins (catalog, Ncatalog); - initGridBins (catalog, Ncatalog); initMrel (catalog, Ncatalog); @@ -74,5 +75,5 @@ MARKTIME("-- set up image indexes: %f sec\n", dtime); - findMosaics (catalog, Ncatalog, TRUE); /* also sets Grid values */ + findMosaics (catalog, Ncatalog, TRUE); MARKTIME("-- set up mosaic indexes: %f sec\n", dtime); @@ -82,7 +83,8 @@ SAVEPLOT = FALSE; + // mark ID_MEAS_AREA (user-selected chip region) and ID_MEAS_NOCAL (time range & active photcode) setExclusions (catalog, Ncatalog, TRUE); - global_stats (catalog, Ncatalog, flatcorr, 0); + global_stats (catalog, Ncatalog, 0); if (PLOTSTUFF) { @@ -91,38 +93,19 @@ } - // 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_OBJ_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_OBJ_POOR; - setMrel (catalog, Ncatalog, flatcorr); - STAR_BAD = ID_OBJ_POOR | ID_OBJ_FEW; - setMgrid (catalog, flatcorr); - } - STAR_BAD = ID_OBJ_POOR | ID_OBJ_FEW; - STAR_TOOFEW = star_toofew; - } - /* determine fit values */ for (i = 0; i < NLOOP; i++) { SetZptIteration (i); - setMrel (catalog, Ncatalog, flatcorr); // threaded - - setMcal (catalog, flatcorr); - setMmos (catalog, flatcorr); - setMgrp (catalog, flatcorr); - - setMgrid (catalog, flatcorr); + setMrel (catalog, Ncatalog); // threaded (calls setMrelCatalog) + + setMcal (catalog); + setMmos (catalog); + setMgrp (catalog); + + setMgrid (catalog, Ncatalog); MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime); if (PLOTSTUFF) { - plot_scatter (catalog, Ncatalog, flatcorr); - plot_grid (catalog, flatcorr); + plot_scatter (catalog, Ncatalog); plot_mosaics (); plot_images (); @@ -131,6 +114,6 @@ } - global_stats (catalog, Ncatalog, flatcorr, i); - SetZeroPointModes (); + global_stats (catalog, Ncatalog, i); + SetZeroPointModes (catalog, Ncatalog); MARKTIME("-- finished loop %d: %f sec\n", i, dtime); @@ -139,6 +122,5 @@ if (PLOTSTUFF) { - plot_scatter (catalog, Ncatalog, flatcorr); - plot_grid (catalog, flatcorr); + plot_scatter (catalog, Ncatalog); plot_mosaics (); plot_images (); @@ -147,5 +129,5 @@ } - if (USE_GRID) dump_grid (); + // if (GRID_ZEROPT) dump_grid (); MARKTIME("-- finalize Mcal values: %f sec\n", dtime); @@ -169,5 +151,6 @@ freeMosaicBins (Ncatalog, TRUE); freeTGroupBins (Ncatalog); - freeGridBins (Ncatalog); + + GridCorrectionSave (); // end of if (NLOOP > 0) block : this loop determines the offsets per chip @@ -179,19 +162,33 @@ } + // once the zero points have been calculated, they are reassigned to the + // individual images. At this point, the mosaic and tgroup values should not + // be used. + freeMosaics (); + freeTGroups (); + MOSAIC_ZEROPT = FALSE; + TGROUP_ZEROPT = FALSE; + + // only change the real database files if -update is requested + if (!UPDATE) { + dvo_image_unlock (&db); + freeImages (db.ftable.buffer); + gfits_db_free (&db); + freeGridBins (); + return TRUE; + } + 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. In a parallel context, this function writes the image parameters as a subset table for the remote clients */ if (!CALIBRATE_STACKS_AND_WARPS) { - reload_catalogs (skylist, flatcorr, 0, NULL); + reload_catalogs (skylist, 0, NULL); // calls setMrelFinal which setMrelOutput which calls setMrelCatalog MARKTIME("-- updated all catalogs: %f sec\n", dtime); } if (UPDATE_CATFORMAT) { - // ensure the db format is updated + // ensure the Image db format is updated db.format = dvo_catalog_catformat (UPDATE_CATFORMAT); gfits_modify (&db.header, "FORMAT", "%s", 1, UPDATE_CATFORMAT); @@ -204,5 +201,6 @@ freeImages (db.ftable.buffer); - freeMosaics (); + freeGridBins (); + gfits_db_free (&db); Index: trunk/Ohana/src/relphot/src/relphot_objects.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_objects.c (revision 41601) +++ trunk/Ohana/src/relphot/src/relphot_objects.c (revision 41647) @@ -6,24 +6,10 @@ // has been correctly propagated to the measurement, the average will respect that value. -int relphot_objects_parallel (SkyList *sky); - int relphot_objects (SkyList *skylist, int hostID, char *hostpath) { - off_t i, j, k; - int Nsecfilt; - struct stat filestat; - - Catalog catalog; + FITS_DB db; INITTIME; - // XXX need to decide how to determine PARALLEL mode... - if (PARALLEL && !hostID) { - relphot_objects_parallel (skylist); - return TRUE; - } - - FITS_DB db; - set_db (&db); gfits_db_init (&db); @@ -41,268 +27,14 @@ MARKTIME("loaded images: %f sec\n", dtime); - // load the ZP corrections here - if (SYNTH_ZERO_POINTS) SynthZeroPointsLoad (SYNTH_ZERO_POINTS); + // load grid corrections here (specified by -grid-meanfile) + GridCorrectionLoad (GRID_MEANFILE); - if (REPAIR_WARPS) { - FindWarpGroups (); - MARKTIME("setup warp groups: %f sec\n", dtime); - MakeStackIndex (); - MARKTIME("setup stack index: %f sec\n", dtime); - RepairWarpMeasuresOpenLogfile (); - } - - // load data from each region file, only use bright stars - for (i = 0; i < skylist[0].Nregions; i++) { - - // does this host ID match the desired location for the table? - if (!HostTableTestHost(skylist[0].regions[i], hostID)) continue; - - // set up the basic catalog info - char hostfile[1024]; - snprintf (hostfile, 1024, "%s/%s.cpt", hostpath, skylist[0].regions[i]->name); - - dvo_catalog_init (&catalog, TRUE); - catalog.filename = hostID ? hostfile : skylist[0].filename[i]; - - // only update existing db tables - int status = stat (catalog.filename, &filestat); - if ((status == -1) && (errno == ENOENT)) { - if (VERBOSE) fprintf (stderr, "no file %s, skipping\n", catalog.filename); - continue; - } - - catalog.catflags = DVO_LOAD_AVERAGE | DVO_LOAD_MEASURE | DVO_LOAD_MISSING | DVO_LOAD_SECFILT; - // if we are going to update the format, we should update all tables - if (UPDATE_CATFORMAT) { - catalog.catflags |= DVO_LOAD_LENSING | DVO_LOAD_LENSOBJ | DVO_LOAD_STARPAR | DVO_LOAD_GALPHOT; - } - catalog.Nsecfilt = GetPhotcodeNsecfilt (); - - if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) { - fprintf (stderr, "ERROR: failure reading catalog %s, skipping\n", catalog.filename); - continue; - } - if (!catalog.Naverage_disk) { - if (VERBOSE) fprintf (stderr, "no data in %s, skipping\n", catalog.filename); - dvo_catalog_unlock (&catalog); - dvo_catalog_free (&catalog); - continue; - } - - // reset the calculated average magnitudes (does not affect ubercal-tied measurements or images) - if (RESET) { - Nsecfilt = catalog.Nsecfilt; - - DVOAverageFlags photomBits = - ID_OBJ_EXT | // extended in our data (eg, PS) - ID_OBJ_EXT_ALT | // extended in external data (eg, 2MASS) - ID_OBJ_GOOD | // good-quality measurement in our data (eg,PS) - ID_OBJ_GOOD_ALT | // good-quality measurement in external data (eg, 2MASS) - ID_OBJ_GOOD_STACK | // good-quality object in the stack (> 1 good stack) - ID_OBJ_BEST_STACK | // the primary stack measurement are the best measurements - ID_OBJ_SUSPECT_STACK | // suspect object in the stack (> 1 good or suspect stack, < 2 good) - ID_OBJ_BAD_STACK; // good-quality object in the stack (> 1 good stack) - - for (j = 0; j < catalog.Naverage; j++) { - catalog.average[j].flags &= ~photomBits; // reset all except astrometry bits - catalog.average[j].psfQF = NAN; // reset (will be re-calculated here) - catalog.average[j].psfQFperf = NAN; // reset (will be re-calculated here) - catalog.average[j].stargal = NAN; // reset (will be re-calculated here) - catalog.average[j].photFlagsUpper = 0; // reset (will be re-calculated here) - catalog.average[j].photFlagsLower = 0; // reset (will be re-calculated here) - catalog.average[j].NwarpOK = 0; // reset (will be re-calculated here) - for (k = 0; k < Nsecfilt; k++) { - dvo_secfilt_init (&catalog.secfilt[j*Nsecfilt + k], SECFILT_RESET_ALL); - } - } - } - - populate_tiny_values(&catalog, DVO_TV_MEASURE | DVO_TV_AVERAGE); - - // XXX need to worry about the image subset data - initImageBins (&catalog, 1, FALSE); - initMosaicBins (&catalog, 1, FALSE); - initGridBins (&catalog, 1); - - findImages (&catalog, 1, FALSE); - findMosaics (&catalog, 1, FALSE); - - // update the detection coordinates using the new image parameters - if (REPAIR_WARPS) RepairWarpMeasures (&catalog); - - initMrel (&catalog, 1); - setMrelFinal (&catalog, NULL, TRUE); - - // uncomment for extra verbosity - // MARKTIME("setMrelFinal for "OFF_T_FMT" average "OFF_T_FMT" measure : %f sec\n", catalog.Naverage, catalog.Nmeasure, dtime); - // gettimeofday (&startTimer, NULL); // reset timer - - // XXX if we want to have options for setting warp, chip, stack independently, we need to init only the desired ones - // NOTE flatcorr == NULL, but it should have been applied already by setphot - - if (!UPDATE) { - dvo_catalog_unlock (&catalog); - free_tiny_values(&catalog); - dvo_catalog_free (&catalog); - freeImageBins (1, FALSE); - freeMosaicBins (1, FALSE); - freeGridBins (1); - continue; - } - - if (VERBOSE) fprintf (stderr, "saving catalog %s\n", catalog.filename); - - char history[128]; - struct timeval now; - gettimeofday (&now, (void *) NULL); - char *moddate = ohana_sec_to_date (now.tv_sec); - gfits_modify (&catalog.header, "RELPHOT", "%s", 1, moddate); - if (REPAIR_WARPS) { - snprintf (history, 128, "repair warp measure.imageID: %s", moddate); - gfits_modify_alt (&catalog.header, "HISTORY", "%S", 0, history); // adds a new entry - } - free (moddate); - - // we can optionally convert output format here - // but it would be better to define a dvo crawler program to do this - // catalog.catformat = DVO_FORMAT_PS1_V1; - - // modify the output format as desired (ignore current format on disk) - if (UPDATE_CATFORMAT) { - catalog.catformat = dvo_catalog_catformat (UPDATE_CATFORMAT); - } - - SetProtect (TRUE); - if (!dvo_catalog_save (&catalog, VERBOSE)) { fprintf (stderr, "ERROR: failed to save %s\n", catalog.filename); exit (1); } - if (!dvo_catalog_unlock (&catalog)) { fprintf (stderr, "ERROR: failed to unlock %s\n", catalog.filename); exit (1); } - SetProtect (FALSE); - - free_tiny_values(&catalog); - dvo_catalog_free (&catalog); - - freeImageBins (1, FALSE); - freeMosaicBins (1, FALSE); - freeGridBins (1); - } - - if (REPAIR_WARPS) { - FreeWarpGroups (); - FreeStackGroups (); - RepairWarpMeasuresCloseLogfile (); - } + reload_catalogs (skylist, hostID, hostpath); freeImages(db.ftable.buffer); + freeGridBins (); + gfits_db_free (&db); return (TRUE); } - -// CATDIR is supplied globally -# define DEBUG 1 -int relphot_objects_parallel (SkyList *sky) { - - // launch the setphot_client jobs to the parallel hosts - - // load the list of hosts - HostTable *table = HostTableLoad (CATDIR, sky->hosts); - if (!table) { - fprintf (stderr, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR); - exit (1); - } - - if (BOUNDARY_TREE) { - char *tmppath = abspath(BOUNDARY_TREE, DVO_MAX_PATH); - free (BOUNDARY_TREE); - BOUNDARY_TREE = tmppath; - } - - int i; - for (i = 0; i < table->Nhosts; i++) { - - // ensure that the paths are absolute path names - char *tmppath = abspath (table->hosts[i].pathname, DVO_MAX_PATH); - free (table->hosts[i].pathname); - table->hosts[i].pathname = tmppath; - - // options / arguments that can affect relphot_client -update-objects: - // VERBOSE, VERBOSE2 - // KEEP_UBERCAL - // RESET (-reset) - // TimeSelect -time - // DophotSelect - // (note that psfQF is applied rigidly at 0.85, as is the galaxy test) - // MAG_LIM - // SIGMA_LIM - // ImagSelect, ImagMin, ImagMax - // MaxDensityUse, MaxDensityValue - - char *command = NULL; - strextend (&command, "relphot_client -update-objects"); - strextend (&command, "-hostID %d", table->hosts[i].hostID); - strextend (&command, "-D CATDIR %s", CATDIR); - strextend (&command, "-hostdir %s", table->hosts[i].pathname); - strextend (&command, "-region %f %f %f %f", UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); - strextend (&command, "-statmode %s", STATMODE); - strextend (&command, "-D CAMERA %s", CAMERA); - strextend (&command, "-D STAR_TOOFEW %d", STAR_TOOFEW); - strextend (&command, "-minerror %f", MIN_ERROR); - - if (VERBOSE) { strextend (&command, "-v"); } - if (VERBOSE2) { strextend (&command, "-vv"); } - if (RESET) { strextend (&command, "-reset"); } - if (RESET_ZEROPTS) { strextend (&command, "-reset-zpts"); } - if (REPAIR_WARPS) { strextend (&command, "-repair-warps"); } - if (PRESERVE_PS1) { strextend (&command, "-preserve-ps1"); } - if (IS_DIFF_DB) { strextend (&command, "-is-diff-db"); } - if (UPDATE) { strextend (&command, "-update"); } - if (!KEEP_UBERCAL) { strextend (&command, "-reset-ubercal"); } - if (UPDATE_CATFORMAT) { strextend (&command, "-update-catformat %s", UPDATE_CATFORMAT); } - if (BOUNDARY_TREE) { strextend (&command, "-boundary-tree %s", BOUNDARY_TREE); } - if (SYNTH_ZERO_POINTS) { strextend (&command, "-synthphot-zpts %s", SYNTH_ZERO_POINTS); } - if (USE_BASIC_CHECK) { strextend (&command, "-basic-image-search"); } - if (USE_ALL_IMAGES) { strextend (&command, "-use-all-images"); } - if (USE_MCAL_PSF_FOR_STACK_APER) { strextend (&command, "-use-mcal-psf-for-stack-aper"); } - - if (!(STAGES & STAGE_CHIP)) { strextend (&command, "-skip-chip"); } - if (!(STAGES & STAGE_WARP)) { strextend (&command, "-skip-warp"); } - if (!(STAGES & STAGE_STACK)) { strextend (&command, "-skip-stack"); } - - // if (SET_MREL_VERSION != 1) { strextend (&command, "-set-mrel-version %d", SET_MREL_VERSION); } - - fprintf (stderr, "command: %s\n", command); - - if (PARALLEL_MANUAL) { - free (command); - continue; - } - - if (PARALLEL_SERIAL) { - int status = system (command); - if (status) { - fprintf (stderr, "ERROR running relphot_client\n"); - exit (2); - } - } else { - // launch the job on the remote machine (no handshake) - int errorInfo = 0; - int pid = rconnect ("ssh", table->hosts[i].hostname, command, table->hosts[i].stdio, &errorInfo, FALSE); - if (!pid) { - if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", table->hosts[i].hostname, errorInfo); - exit (1); - } - table->hosts[i].pid = pid; // save for future reference - } - free (command); - } - - if (PARALLEL_MANUAL) { - fprintf (stderr, "run the relphot_client commands above. when these are done, hit return\n"); - getchar(); - } - if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) { - HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE); - } - - FreeHostTable (table); - return TRUE; -} Index: trunk/Ohana/src/relphot/src/relphot_parallel_images.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_parallel_images.c (revision 41601) +++ trunk/Ohana/src/relphot/src/relphot_parallel_images.c (revision 41647) @@ -8,4 +8,6 @@ * update the unowned detections for owned objects from neighbor regions * update the unowned objects for owned detections + + * TGROUP and GRID corrections are not calculated. */ @@ -39,13 +41,4 @@ initImages (image, NULL, Nimage); - - // load the flat correction table (if defined) - // we no longer need to do this: the flatcorrection is now saved in measure.Mflat -# if (0) - char flatcorrfile[256]; - sprintf (flatcorrfile, "%s/flatcorr.fits", CATDIR); - FlatCorrectionTable *flatcorr = FlatCorrectionLoad (flatcorrfile, VERBOSE); -# endif - FlatCorrectionTable *flatcorr = NULL; // UserPatch.Rmin,Rmax may have range from a few degrees < 0.0 to few degrees > 360.0. @@ -88,5 +81,5 @@ setExclusions (catalog, Ncatalog, TRUE); - global_stats (catalog, Ncatalog, flatcorr, 0); + global_stats (catalog, Ncatalog, 0); if (PLOTSTUFF) { @@ -101,5 +94,5 @@ // set the mean stellar mags given the measurements and the image calibrations - setMrel (catalog, Ncatalog, flatcorr); + setMrel (catalog, Ncatalog); // share mean mags for objects at the boundary (number of unowned meas > 0) @@ -112,7 +105,6 @@ // set the image (Mcal) and mosaic (Mmos) zero point offsets given the mean mags and measurements - setMcal (catalog, flatcorr); - setMmos (catalog, flatcorr); - // setMgrp (catalog, FALSE, flatcorr); XXX think this through: this may not make sense for tgroups + setMcal (catalog); + setMmos (catalog); MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime); @@ -124,7 +116,7 @@ client_logger_message ("slurped image data : loop %d \n", i); - global_stats (catalog, Ncatalog, flatcorr, i); + global_stats (catalog, Ncatalog, i); - SetZeroPointModes (); + SetZeroPointModes (catalog, Ncatalog); MARKTIME("-- finished loop %d: %f sec\n", i, dtime); } @@ -147,5 +139,4 @@ setMcalFromMosaics (); // copy per-mosaic calibrations to the images - // XXX keep? setMcalFromTGroups (); // copy per-tgroup calibrations to the images share_image_mags (regionHosts, -1); @@ -161,5 +152,4 @@ freeImageBins(Ncatalog, TRUE); freeMosaicBins (Ncatalog, TRUE); - freeGridBins (Ncatalog); freeImages((char *)image); free (image); Index: trunk/Ohana/src/relphot/src/relphot_parallel_regions.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_parallel_regions.c (revision 41601) +++ trunk/Ohana/src/relphot/src/relphot_parallel_regions.c (revision 41647) @@ -14,12 +14,4 @@ exit (2); } - - // load the flat correction table (if defined) -# if (0) - char flatcorrfile[256]; - sprintf (flatcorrfile, "%s/flatcorr.fits", CATDIR); - FlatCorrectionTable *flatcorr = FlatCorrectionLoad (flatcorrfile, VERBOSE); -# endif - FlatCorrectionTable *flatcorr = NULL; // register database handle with shutdown procedure @@ -81,5 +73,5 @@ /* update catalogs (in parallel) */ - reload_catalogs (skylist, flatcorr, 0, NULL); + reload_catalogs (skylist, 0, NULL); // save the changes to the image parameters Index: trunk/Ohana/src/relphot/src/select_images.c =================================================================== --- trunk/Ohana/src/relphot/src/select_images.c (revision 41601) +++ trunk/Ohana/src/relphot/src/select_images.c (revision 41647) @@ -103,8 +103,5 @@ if (!(i % 300000)) fprintf (stderr, "."); - // XXX I am just going for force this for the moment: - // goto found_it; - - /* exclude images by photcode (we are not doing this, but I'm not sure why) */ + // only include active photcodes in the analysis if (timage[i].photcode) { int Ns = GetActivePhotcodeIndex (timage[i].photcode); Index: trunk/Ohana/src/relphot/src/setMrelCatalog.c =================================================================== --- trunk/Ohana/src/relphot/src/setMrelCatalog.c (revision 41601) +++ trunk/Ohana/src/relphot/src/setMrelCatalog.c (revision 41647) @@ -60,5 +60,5 @@ int useOLS = TRUE; // chosen based on iteration below -int setMrel_catalog_alt (Catalog *catalog, int Nc, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) { +int setMrelCatalog (Catalog *catalog, int Nc, int isSetMrelFinal, SetMrelInfo *results, int Nsecfilt) { off_t j; @@ -75,15 +75,15 @@ if (STAGES & STAGE_CHIP) { - setMrelAverageExposure (&catalog[Nc], Nc, j, Nsecfilt, isSetMrelFinal, flatcorr, results); + setMrelAverageExposure (&catalog[Nc], Nc, j, Nsecfilt, isSetMrelFinal, results); } // only apply Stack operation on setMrelFinal in first pass if (isSetMrelFinal && (STAGES & STAGE_STACK) && !IS_DIFF_DB) { - setMrelAverageStack (&catalog[Nc], Nc, j, Nsecfilt, flatcorr); + setMrelAverageStack (&catalog[Nc], Nc, j, Nsecfilt); } // only measure force-warp mean values if issetMrelFinal (make it optional?) if (isSetMrelFinal && (STAGES & STAGE_WARP)) { - setMrelAverageForcedWarp (&catalog[Nc], Nc, j, Nsecfilt, flatcorr, results); + setMrelAverageForcedWarp (&catalog[Nc], Nc, j, Nsecfilt, results); } } @@ -113,5 +113,5 @@ // set mean of chip measurements (selected by photcode range for now): -int setMrelAverageExposure (Catalog *catalog, int cat, off_t ave, int Nsecfilt, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results) { +int setMrelAverageExposure (Catalog *catalog, int cat, off_t ave, int Nsecfilt, int isSetMrelFinal, SetMrelInfo *results) { int Nsec; @@ -225,4 +225,18 @@ } + /* some things to note here: + + Mcal, Mmos, Mgrid are only relevant (non-zero) during the relphot_images and + relphot_parallel_images analysis steps. For the final output steps (setMrelFinal), + these are turned off + + During the image analysis, Mflat is the prior stored value while Mgrid is the new + measurement. In the final output step, Mgrid is transferred to Mflat and zero + + On the final calculation of Mrel, the value of Mgrid has been applied to Mflat, but + it has not been removed from the GridOps.c structures. In order to avoid + double-counting, we need to skip Mgrid on the final calculation. + */ + // ** Choose the calibration (depends on the mode : do I have an image reference or not?) float Mcal = 0.0, Mmos = 0.0, Mgrid = 0.0, Mflat = 0.0, Mgrp = 0.0; @@ -233,17 +247,22 @@ // data for which the associated image has not been loaded (probably because of // overlaps). Msys + measure.Mcal is our best guess of the true magnitude - Mmos = Mgrid = 0; Mcal = measureT[k].McalPSF; // check that this is zero for loaded REF value } else { - // getMcal returns image[].Mcal; note the flat-field correction is stored in measure.Mflat + // getMcal returns image[].Mcal Mcal = getMcal (meas, cat, MAG_CLASS_PSF); if (isnan(Mcal)) SKIP_THIS_MEAS(Ncal); - Mmos = getMmos (meas, cat); - if (isnan(Mmos)) SKIP_THIS_MEAS(Nmos); - Mgrp = getMgrp (meas, cat, measureT[k].airmass, &dMgrp); - if (isnan(Mgrp)) SKIP_THIS_MEAS(Ngrp); - Mgrid = getMgrid (meas, cat); - if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid); + + // the flat-field correction is stored in measure.Mflat Mflat = isnan (measureT[k].Mflat) ? 0.0 : measureT[k].Mflat; + + // see note above re: final output vs image analysis + if (!isSetMrelFinal) { + Mmos = getMmos (meas, cat); + if (isnan(Mmos)) SKIP_THIS_MEAS(Nmos); + Mgrp = getMgrp (meas, cat, measureT[k].airmass, &dMgrp); + if (isnan(Mgrp)) SKIP_THIS_MEAS(Ngrp); + Mgrid = getMgridTiny (&measureT[k]); + if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid); + } } @@ -269,6 +288,6 @@ float dMap = NAN; - // XXX I am not convinced the sign on Mflat here is correct - float Moff = Mflat + Mcal + Mmos + Mgrp + Mgrid; + // This definition is consistent with PhotRel: Mrel = Msys - Mcal - Mflat + float Moff = Mcal + Mmos + Mgrp + Mgrid + Mflat; if (isSetMrelFinal) { @@ -349,5 +368,6 @@ for (Nsec = 0; Nsec < Nsecfilt; Nsec++) { - // if we detected this object in PS1, or do not request -preserve-ps1, keep the mean photometry values + // -preserve-ps1 means keep an existing average PS1 value + // if we did not detect this object in PS1, or do not request -preserve-ps1, keep the mean photometry values if (!PRESERVE_PS1 || !(secfilt[Nsec].flags & ID_SECF_HAS_PS1)) { dvo_secfilt_init (&secfilt[Nsec], SECFILT_RESET_CHIP); // this does not reset astrometry or STACK bits @@ -364,6 +384,7 @@ } + // if -preserve-ps1 is set and this object has PS1 data, skip the rest of the steps: + // NOTE: This test only applies if this bit was set by an earlier analysis if (PRESERVE_PS1 && (secfilt[Nsec].flags & ID_SECF_HAS_PS1)) continue; - // if -preserve-ps1 is set and this object has PS1 data, skip the rest of the steps: // XXX hardwired grizy = (01234) JHK = (567) w = (8) @@ -416,9 +437,10 @@ int NrankingPSF = useOLS ? magStatsByRanking (&results->psfData[Nsec], psfstats) : magStatsByRankingIRLS (&results->psfData[Nsec], psfstats); if (NrankingPSF < Nminmeas) { - secfilt[Nsec].flags |= ID_OBJ_FEW; + secfilt[Nsec].flags |= ID_SECF_STAR_FEW; } else { secfilt[Nsec].MpsfChp = psfstats->mean; secfilt[Nsec].dMpsfChp = psfstats->error; secfilt[Nsec].Mchisq = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN; + secfilt[Nsec].flags &= ~ID_SECF_STAR_FEW; } @@ -552,5 +574,5 @@ // measure->McalAPER to image->McalPSF for chips and warps, but not stacks -int setMrelAverageStack (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr) { +int setMrelAverageStack (Catalog *catalog, int cat, off_t ave, int Nsecfilt) { // we are guaranteed to have average, measure, secfilt @@ -562,5 +584,5 @@ SecFilt *secfilt = &catalog[0].secfilt[ave*Nsecfilt]; - float McalPSF = 0, McalAPER = 0, Mmos = 0, Mgrid = 0; + float McalPSF = 0, McalAPER = 0; // set the primary projection cell and skycell for this coordinate @@ -604,18 +626,5 @@ // reset all stack-related values for this secfilt: - secfilt[Nsec].stackBestOff = -1; - secfilt[Nsec].stackPrmryOff = -1; - secfilt[Nsec].FpsfStk = NAN; - secfilt[Nsec].dFpsfStk = NAN; - secfilt[Nsec].FkronStk = NAN; - secfilt[Nsec].dFkronStk = NAN; - secfilt[Nsec].FapStk = NAN; - secfilt[Nsec].dFapStk = NAN; - secfilt[Nsec].MpsfStk = NAN; - secfilt[Nsec].MkronStk = NAN; - secfilt[Nsec].MapStk = NAN; - secfilt[Nsec].Nstack = 0; - secfilt[Nsec].NstackDet = 0; - secfilt[Nsec].flags &= ~ID_SECF_STACK_FLAGS; // reset the stack flags + dvo_secfilt_init (&secfilt[Nsec], SECFILT_RESET_STACK); for (off_t k = 0; k < Nmeasure; k++) { @@ -638,6 +647,4 @@ haveStack = TRUE; haveStackObject = TRUE; - - // if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS_STACK(Nbad); int isPrimary = FALSE; @@ -726,12 +733,7 @@ McalPSF = measure[meas].McalPSF; // check that this is zero for loaded REF value McalAPER = McalPSF; // check that this is zero for loaded REF value - Mmos = 0.0; - Mgrid = 0.0; } else { McalPSF = getMcal (measSeq, cat, MAG_CLASS_PSF); McalAPER = USE_MCAL_PSF_FOR_STACK_APER ? getMcal (measSeq, cat, MAG_CLASS_PSF) : getMcal (measSeq, cat, MAG_CLASS_KRON); - Mmos = getMmos (measSeq, cat); - Mgrid = getMgrid (measSeq, cat); - // XXX can Mmos and Mgrid exist for stacks? } @@ -769,6 +771,6 @@ // get the zero point for the selected image // Use a different zero point for the PSF-like and APERTURE-like magnitudes - float zpPSF = PhotZeroPoint (&measure[meas], &average[0], &secfilt[0]) - (McalPSF + Mmos + Mgrid); - float zpAPER = PhotZeroPoint (&measure[meas], &average[0], &secfilt[0]) - (McalAPER + Mmos + Mgrid); + float zpPSF = PhotZeroPoint (&measure[meas], &average[0], &secfilt[0]) - McalPSF; + float zpAPER = PhotZeroPoint (&measure[meas], &average[0], &secfilt[0]) - McalAPER; // zpFactor to go from instrumental flux to Janskies @@ -884,6 +886,5 @@ // chips and warps, but not stacks -int setMrelAverageForcedWarp (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr, SetMrelInfo *results) { - OHANA_UNUSED_PARAM(flatcorr); +int setMrelAverageForcedWarp (Catalog *catalog, int cat, off_t ave, int Nsecfilt, SetMrelInfo *results) { // we are guaranteed to have average, measure, secfilt Index: trunk/Ohana/src/relphot/src/setMrelFinal.c =================================================================== --- trunk/Ohana/src/relphot/src/setMrelFinal.c (revision 41601) +++ trunk/Ohana/src/relphot/src/setMrelFinal.c (revision 41647) @@ -5,5 +5,5 @@ // output dbFlags values -void setMrelFinal (Catalog *catalog, FlatCorrectionTable *flatcorr, int simpleAverage) { +void setMrelFinal (Catalog *catalog, int simpleAverage) { off_t i; @@ -11,75 +11,8 @@ int Nsecfilt = GetPhotcodeNsecfilt (); - /*** RESET photometry flags and Mcal values ***/ - if (RESET) { - // flags used by the photometry analysis (excluding UBERCAL) - unsigned int PHOTOM_FLAGS = - ID_MEAS_NOCAL | // detection ignored for this analysis (photcode, time range) - ID_MEAS_POOR_PHOTOM | // detection is photometry outlier - ID_MEAS_SKIP_PHOTOM | // detection was ignored for photometry measurement - ID_MEAS_AREA | // detetion was outside acceptable area of device - ID_MEAS_SYNTH_MAG | // magnitude is synthetic - ID_MEAS_STACK_PRIMARY | // this stack measurement is in the primary skycell - ID_MEAS_STACK_PHOT_SRC; // this measurement supplied the stack photometry - - // ID_MEAS_PHOTOM_UBERCAL -- externally-supplied zero point from ubercal analysis - // this is set by 'setphot', do not reset here - - // flags used by the photometry analysis (excluding UBERCAL) - // unsigned int secfiltFlags = - // ID_PHOTOM_PASS_0 | // average measured at pass 0 - // ID_PHOTOM_PASS_1 | // average measured at pass 1 - // ID_PHOTOM_PASS_2 | // average measured at pass 2 - // ID_PHOTOM_PASS_3 | // average measured at pass 3 - // ID_PHOTOM_PASS_4 | // average measured at pass 3 - // ID_SECF_USE_SYNTH | // average measured at pass 3 - // ID_SECF_USE_UBERCAL | // average measured at pass 3 - // ID_SECF_OBJ_EXT; // average measured at pass 3 - - // XXX I should really deprecate the concept of applying the average - // calculation to a limited set of photcodes. - // for now, just do all photcodes here ( - // as it stands, only stacks are limited by photcode; mean exp and forced warp - // are applied to all Nsecfilt - - for (i = 0; i < catalog[0].Naverage; i++) { - int Ns; - for (Ns = 0; Ns < Nsecfilt; Ns++) { - - off_t N = Nsecfilt*i+Ns; - dvo_secfilt_init (&catalog[0].secfilt[N], SECFILT_RESET_ALL); - - off_t m = catalog[0].average[i].measureOffset; - off_t j; - for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) { - - /* select measurements by time */ - if (TimeSelect) { - if (catalog[0].measure[m].t < TSTART) continue; - if (catalog[0].measure[m].t > TSTOP) continue; - } - - // only reset Mcal for measures with a matching image - // do not reset Mcal for ubercal images unless explicitly requested - - if (catalog[0].measure[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) { - if (!KEEP_UBERCAL) { - catalog[0].measure[m].McalPSF = 0.0; - catalog[0].measure[m].McalAPER = 0.0; - catalog[0].measure[m].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL; - } - } else { - if (RESET_ZEROPTS && (getImageEntry (m, 0) >= 0)) { - catalog[0].measure[m].McalPSF = 0.0; - catalog[0].measure[m].McalAPER = 0.0; - } - } - catalog[0].measure[m].dbFlags &= ~PHOTOM_FLAGS; - } - } - } - } - - // this sets flags in the measureT element, not the measure element + ResetAverageAndMeasure (catalog); + + // This sets ID_MEAS flags in the measureT element, not the measure element + // These bits are thus ephemeral and not saved. setExclusions (catalog, 1, VERBOSE); /* mark by area */ @@ -92,53 +25,11 @@ } - for (i = 0; i < catalog[0].Naverage; i++) { - catalog[0].average[i].psfQF = NAN; // force recalculation below - catalog[0].average[i].psfQFperf = NAN; // force recalculation below - catalog[0].average[i].stargal = NAN; // force recalculation below - catalog[0].average[i].photFlagsUpper = 0; // reset (will be re-calculated here) - catalog[0].average[i].photFlagsLower = 0; // reset (will be re-calculated here) - catalog[0].average[i].NwarpOK = 0; // reset (will be re-calculated here) - } - - // XXX make this optional? (do not clean for -averages?) - // XXX : for the moment, disable clean_measures - if (!simpleAverage && FALSE) clean_measures (catalog, 1, TRUE, flatcorr); /* mark outliers ID_MEAS_POOR_PHOTOM */ - ALLOCATE (catalog[0].measureRank, char, catalog[0].Nmeasure); setMeasureRank (catalog); - setMrelOutput (catalog, 1, flatcorr); // sets the values secfilt.MpsfChp = - setMcalOutput (catalog, 1, flatcorr); // sets measure.Mcal = image.Mcal - /* clear ID_OBJ_POOR, ID_OBJ_FEW values before writing ??? */ - /* ID_MEAS_NOCAL is an internal bit, so it should be cleared */ - // XXX eventually: move these internal flags to an internal-only catalog->FOO element (like found) - for (i = 0; i < catalog[0].Naverage; i++) { - catalog[0].average[i].flags &= ~ID_OBJ_FEW; - catalog[0].average[i].flags &= ~ID_OBJ_POOR; - off_t j; - for (j = 0; j < Nsecfilt; j++) { - catalog[0].secfilt[i*Nsecfilt+j].flags &= ~ID_OBJ_FEW; - catalog[0].secfilt[i*Nsecfilt+j].flags &= ~ID_OBJ_POOR; - } - off_t m = catalog[0].average[i].measureOffset; - for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) { - catalog[0].measure[m].dbFlags &= ~ID_MEAS_NOCAL; - } - } + setMflatFromGrid (catalog); // Mgrid is used to set Mflat; Mgrid is the ignored in setMrelOutput / setMrelCatalog + setMrelOutput (catalog, 1); // sets the values secfilt.MpsfChp = + setMcalOutput (catalog, 1); // sets measure.Mcal = image.Mcal } - -/* ID_MEAS_SKIP marks measurements which were not used to calculate Mrel */ -// For each star & photcode, we are selecting the measurements to use in the average -// photometry value. This function is called several times, lowering the bar to -// acceptance on each pass. The rejections at the different passes are: - -// 0 : only use measurements thought to be GOOD (photflags not POOR or BAD) -// 1 : accept measurements thought to be POOR based on photflags -// 2 : accept the measurements marked as outliers; accept images marked as outliers -// 3 : accept measurements thought to be BAD based on photflags -// 4 : accept the measurements outside of the instrumental magnitude limits (eg, SAT) - -// XXX if I'm reading this correctly, passes 0, 1, 2 are identical & accept all -// non-outliers with Imag range // int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure) { @@ -161,126 +52,11 @@ // } -int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure); - -void skip_measurements (Catalog *catalog, int pass, FlatCorrectionTable *flatcorr) { - - off_t i, k, m; - off_t Ntot, Ntry, Nkeep, Nskip; - float mag; - - int Nsecfilt = GetPhotcodeNsecfilt (); - - Ntot = Ntry = Nskip = Nkeep = 0; - - /* allow measures from images marked POOR and FEW */ - if (pass >= 3) IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL; - - /* allow measures marked as outliers (POOR) and off image region (AREA) */ - if (pass >= 3) { - MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_SKIP_PHOTOM; - } else { - MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA; - } - - /* mark measures which should be ignored on this pass */ - for (i = 0; i < catalog[0].Naverage; i++) { - Ntot += catalog[0].average[i].Nmeasure; - - if (FALSE && (catalog[0].average[i].objID == 0x46a4) && (catalog[0].average[i].catID == 0xf40e)) { - fprintf (stderr, "test obj\n"); - print_measure_set (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], catalog[0].measure); - } - - // mark measurements for each secfilt separately - int Ns; - for (Ns = 0; Ns < Nphotcodes; Ns++) { - - int thisCode = photcodes[Ns][0].code; - int Nsec = GetPhotcodeNsec(thisCode); - - /* star/photcodes already calibrated */ - if (catalog[0].found_t[Nsecfilt*i+Nsec]) continue; - - m = catalog[0].average[i].measureOffset; - for (k = 0; k < catalog[0].average[i].Nmeasure; k++, m++) { - Ntry++; - - // skip measurements not related to this photcode - PhotCode *code = GetPhotcodebyCode (catalog[0].measure[m].photcode); - if (!code) continue; - if (code->equiv != thisCode) continue; - - // skip measurements by time range (mark as skipped for this secfilt) - // XXX note that this is a bit dangerous : some objects may never get calibrated - if (TimeSelect) { - if (catalog[0].measure[m].t < TSTART) goto skip; - if (catalog[0].measure[m].t > TSTOP) goto skip; - } - - // we now have a measurement relavant to this photcode (and time range) - - // clear SKIP for all measures at first - catalog[0].measureT[m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM; - catalog[0].measure [m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM; - - // skip measurements from BAD images and mosaics (not REF mags) - // do NOT skip measurements without a matching image (REF mags) - off_t Nim = getImageEntry (m, 0); - if (Nim > -1) { - if (isnan(getMcal (m, 0, MAG_CLASS_PSF))) goto skip; - if (isnan(getMmos (m, 0))) goto skip; - } - - // PASS 4 : skip measurements by inst mag limit (not REF mags) - if ((pass < 4) && ImagSelect) { - if (Nim > -1) { - mag = PhotInst (&catalog[0].measure[m], MAG_CLASS_PSF); - if (mag < ImagMin) goto skip; - if (mag > ImagMax) goto skip; - } - } - - // PASS 3 : accept bad measurements (eg, SAT, CR), internal outliers - if (pass < 3) { - if (catalog[0].measure[m].photFlags & code->photomBadMask) goto skip; - if ((catalog[0].measure[m].photcode > 10000) && (catalog[0].measure[m].photcode < 10500)) { - if (catalog[0].measure[m].psfQF < 0.85) goto skip; - } - } - - // PASS 2 : (no additional cut) - - // PASS 1 : accept poor measurements as well (eg, POOR FIT, etc) - if ((pass < 1) && (catalog[0].measure[m].photFlags & code->photomPoorMask)) goto skip; - - // PASS 0 : skip poor, outlier, bad measurements as well (eg, POOR FIT, etc) - - Nkeep ++; - - continue; - - skip: - catalog[0].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM; - catalog[0].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM; - Nskip ++; - } - } - } - if (VERBOSE) fprintf (stderr, "pass %d, Ntot: "OFF_T_FMT", Ntry: "OFF_T_FMT", Nskip: "OFF_T_FMT", Nkeep: "OFF_T_FMT"\n", - pass, Ntot, Ntry, Nskip, Nkeep); -} - +// This function is only called for the final output step. By this point, we have +// propagated the mosaic and tgroup flags to each image. void setMeasureRank (Catalog *catalog) { int i; - // make these global - int IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL; - int IMAGE_POOR = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP; - int MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_SKIP_PHOTOM; - int MEAS_POOR = ID_MEAS_POOR_PHOTOM | ID_MEAS_AREA; - Measure *measure = catalog[0].measure; - MeasureTiny *measureT = catalog[0].measureT; char *measureRank = catalog[0].measureRank; @@ -289,9 +65,5 @@ measureRank[i] = 11; // start at a low rank - // clear SKIP for all measures at first - measureT[i].dbFlags &= ~ID_MEAS_SKIP_PHOTOM; - measure [i].dbFlags &= ~ID_MEAS_SKIP_PHOTOM; - - // skip measurements without a valid photcode + // measurements without a valid photcode have lowest rank (should not be used anyway) PhotCode *code = GetPhotcodebyCode (measure[i].photcode); if (!code) continue; @@ -299,49 +71,50 @@ // measurements outside time range have poor rank if (TimeSelect) { - if (measure[i].t < TSTART) { measureRank[i] = 10; continue; } - if (measure[i].t > TSTOP) { measureRank[i] = 10; continue; } + if (measure[i].t < TSTART) { measureRank[i] = 10; continue; } + if (measure[i].t > TSTOP) { measureRank[i] = 10; continue; } } - // measurements from BAD images and mosaics (not REF mags) - off_t Nim = getImageEntry (i, 0); + // Nim < 0 for REF mags, imageFlags have bits for IMAGE, MOSAIC, NIGHT + off_t Nim = getImageEntry (i, 0); int imageFlags = getImageFlags (i, 0); - int mosaicFlags = getMosaicFlags (i, 0); if (Nim > -1) { - if (imageFlags & IMAGE_BAD) { measureRank[i] = 9; continue; } - - if (mosaicFlags & IMAGE_BAD) { measureRank[i] = 9; continue; } - // measurements ranked by inst mag limit (not REF mags) if (ImagSelect) { float mag = PhotInst (&measure[i], MAG_CLASS_PSF); - if (mag < ImagMin) { measureRank[i] = 8; continue; } - if (mag > ImagMax) { measureRank[i] = 8; continue; } + if (mag < ImagMin) { measureRank[i] = 9; continue; } + if (mag > ImagMax) { measureRank[i] = 9; continue; } } } - // RANK 1 : BAD photFlags - if (measure[i].dbFlags & MEAS_BAD) { measureRank[i] = 7; continue; } - - // RANK 3 : BAD photFlags (eg, SAT, CR), internal outliers - if (measure[i].photFlags & code->photomBadMask) { measureRank[i] = 6; continue; } - + // RANK 8 : Poor image if (Nim > -1) { - if (imageFlags & IMAGE_POOR) { measureRank[i] = 5; continue; } - if (mosaicFlags & IMAGE_POOR) { measureRank[i] = 5; continue; } + if (imageFlags & ID_IMAGE_PHOTOM_POOR) { measureRank[i] = 8; continue; } } - // RANK 2 : psfQF value - if (!isfinite(measure[i].psfQF) || measure[i].psfQF < 0.85) { measureRank[i] = 4; continue; } + // RANK 7 : BAD photFlags (eg, SAT, CR), internal outliers + if (measure[i].photFlags & code->photomBadMask) { measureRank[i] = 7; continue; } + + // RANK 6 : bad psfQF value + if (!isfinite(measure[i].psfQF) || measure[i].psfQF < 0.85) { measureRank[i] = 6; continue; } - // RANK 1 : POOR photFlags - if (measure[i].dbFlags & MEAS_POOR) { measureRank[i] = 3; continue; } + // RANK 5 : not in valid chip region + if (measure[i].dbFlags & ID_MEAS_AREA) { measureRank[i] = 5; continue; } - // RANK 1 : POOR photFlags - if (measure[i].photFlags & code->photomPoorMask) { measureRank[i] = 2; continue; } + // RANK 4 : POOR photFlags + if (measure[i].photFlags & code->photomPoorMask) { measureRank[i] = 4; continue; } - // RANK 1 : psfQFperfect value - if (!isfinite(measure[i].psfQFperf) || measure[i].psfQFperf < 0.85) { measureRank[i] = 1; continue; } + // RANK 3 : bad psfQFperfect value + if (!isfinite(measure[i].psfQFperf) || measure[i].psfQFperf < 0.85) { measureRank[i] = 3; continue; } + // RANK 2 : Poor mosaic + if (Nim > -1) { + if (imageFlags & ID_IMAGE_MOSAIC_POOR) { measureRank[i] = 2; continue; } + } + + // RANK 1 : Poor night + if (Nim > -1) { + if (imageFlags & ID_IMAGE_NIGHT_POOR) { measureRank[i] = 1; continue; } + } // RANK 0 : perfect measurement: measureRank[i] = 0; Index: trunk/Ohana/src/relphot/test/relphot.grid.dvo =================================================================== --- trunk/Ohana/src/relphot/test/relphot.grid.dvo (revision 41647) +++ trunk/Ohana/src/relphot/test/relphot.grid.dvo (revision 41647) @@ -0,0 +1,520 @@ +# -*-sh-*- + +$TAP_BREAK = 0 +$TAP_VERBOSE = 1 +input tap.dvo +input ptolemy.dvo + +# set globals +if (not($?PLOT)) set PLOT = 0 + +$USE_CELL_OFFSETS = 1; # turn on non-zero cell offsets +$SET_UC_RESID = 1; # write a non-zero resid value to the ubercal table (should not have an impact) +$DO_2MASS = 1 +$DO_SYNTH = 1 +$DO_TYCHO = 1 +$NC_CLOUDS = 0.20 +$UC_CLOUDS = 0.02 + +## full test set: generate a set of test files, insert into a catdir, check raw catalog, run relphot, check results +macro go + if ($0 != 3) + echo "go (fileroot) (catdir)" + break + end + + local i fileroot rootdir catdir + + $fileroot = $1 + $catdir = $2 + + dirname $1 -var rootdir + mkdir $rootdir + + break -auto off + exec rm -rf $catdir + exec rm -rf $fileroot* + break -auto on + + # define a set of photometric (ubercal) and non-photometric exposure times + init + + # generate the common list of stars (coords & mags) for all exposures + mkstars 1000 + + catdir $catdir + + tapPLAN {((mjd_uc[]*4 + mjd_nc[]*4) + (mjd_uc[] + mjd_nc[]) + (mjd_uc[] + mjd_nc[])) / 3} + + # XXX first, run the analysis with only ubercal data (photometric data) + # we should be able to recover the nightly zero points fitting only tgroups + + # generate the basic images and check they were correctly ingested by dvo + for i 0 mjd_uc[] + sprintf filerootout "%s.uc.%02d" $fileroot $i + mkexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] secz_uc[$i] $FILTER $filerootout {$RA_CENTER + dR_uc[$i]} {$DEC_CENTER + dD_uc[$i]} + ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] secz_uc[$i] $FILTER raw + end + + for i 0 mjd_nc[] + sprintf filerootout "%s.nc.%02d" $fileroot $i + mkexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] secz_nc[$i] $FILTER $filerootout $RA_CENTER $DEC_CENTER + ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] secz_nc[$i] $FILTER raw + end + + # run relphot on the db and check that the images now match the expected values + echo relphot -grid -tgroup-fit-airmass -images g -v -region $RA_MIN $RA_MAX $DEC_MIN $DEC_MAX -D CATDIR $catdir -D STAR_TOOFEW 1 -statmode WT_MEAN -cloud-limit 0.5 -update -nloop 25 + # tapEXEC relphot -tgroup-fit-airmass -tgroups tgroups.dat -images g -v -region $RA_MIN $RA_MAX $DEC_MIN $DEC_MAX -D CATDIR $catdir -D STAR_TOOFEW 1 -statmode WT_MEAN -cloud-limit 0.5 -update -nloop 12 + + break + + # reset the per-filter ZPT_OFF:$filter values + $ZPT_OFF:$FILTER = NAN + + for i 0 mjd_uc[] + ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $FILTER relphot_uc + end + for i 0 mjd_nc[] + ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $FILTER relphot_nc + end + tapDONE +end + +# set various global variables +macro init + + # generates ptolemy.rc, dvo.photcodes, $zpt_nominal:g,r $klam_nominal:g,r + mkptolemy.rc + + $RA_CENTER = 10.0 + $DEC_CENTER = 20.0 + $PLATE_SCALE = 0.25 + $NCHIP_X = 1 + $NCHIP_Y = 1 + $NCELL_X = 2 + $NCELL_Y = 2 + $CHIP_DX = 1000 + $CHIP_DY = 1000 + $DEFECT_FRAC = 0.00 + + # we have two sets of images: ubercal (photometric) and not-ubercaled + + # ********* UBERCAL IMAGES *********** + + # ubercal nights each have a distinct zero point and slope: + vlist night_uc 55000.0 55010.0 55030.0 55060.0 55080.0 55090.0 + vlist dz_nt_uc +0.0000 +0.0300 -0.0200 -0.0400 0.0100 0.0150 +# vlist dk_nt_uc +0.0000 +0.0005 -0.0005 -0.0002 + vlist dk_nt_uc +0.0000 +0.0200 -0.0200 +0.0300 0.0005 0.0050 + + write tgroups.dat night_uc + + $EXPTIME = 10.0 + $FILTER = g + $ZPT_OFF:$FILTER = NAN + + delete -q mjd_uc seq_uc zpt_uc klam_uc secz_uc exptime_uc dR_uc dD_uc + + # generate 5 exposures for each night + for i 0 16 + set secz_tmp = rnd(night_uc) + 1.0; # airmass in range 1.0 - 2.0 + set klam_tmp = $klam_nominal:$FILTER + dk_nt_uc + set exp_tmp = zero(night_uc) + $EXPTIME + set zpt_tmp = $zpt_nominal:$FILTER + dz_nt_uc + 2.5*log(exp_tmp) + klam_tmp*(secz_tmp - 1.0) + + set seq_tmp = ramp(night_uc) + set mjd_tmp = night_uc + 0.33 + 0.01*$i; # exposures are 14 minutes apart starting at 10pm HST + set dR_tmp = zero(night_uc) + (0.25*($i % 4) - 0.5) * $PLATE_SCALE * $CHIP_DX / 3600.0 ; # offset in arcsec + set dD_tmp = zero(night_uc) + (0.25*int($i / 4) - 0.5) * $PLATE_SCALE * $CHIP_DY / 3600.0 ; # offset in arcsec + + concat klam_tmp klam_uc + concat secz_tmp secz_uc + concat zpt_tmp zpt_uc + concat seq_tmp seq_uc + concat mjd_tmp mjd_uc + concat dR_tmp dR_uc + concat dD_tmp dD_uc + + concat exp_tmp exptime_uc + end + + # ********* NON-UBERCAL IMAGES *********** + + # non-ubercal nights each have a random zero points and fixed slope (klam_nominal) + # non-ubercal nights must not be the same set as the ubercal nights above + vlist night_nc 55005.0 55015.0 55035.0 55065.0 + + delete -q mjd_nc seq_nc zpt_nc klam_nc secz_nc + + # generate 5 exposures for each night + for i 0 5 + set secz_tmp = rnd(night_nc) + 1.0; # airmass in range 1.0 - 2.0 + set klam_tmp = zero(night_nc) + $klam_nominal:$FILTER + set exp_tmp = zero(night_nc) + $EXPTIME + set zpt_tmp = $zpt_nominal:$FILTER + 2.5*log(exp_tmp) + klam_tmp*(secz_tmp - 1.0) - $NC_CLOUDS*(0.1 + rnd(secz_tmp)) + # observed zero points for nc nights are 0.1 - 1.1 lower than the ubercal nights + + set seq_tmp = ramp(night_nc) + set mjd_tmp = night_nc + 0.33 + 0.01*$i; # exposures are 14 minutes apart starting at 10pm HST + + concat klam_tmp klam_nc + concat secz_tmp secz_nc + concat zpt_tmp zpt_nc + concat seq_tmp seq_nc + concat mjd_tmp mjd_nc + + concat exp_tmp exptime_nc + end +end + +macro try.relphot + if ($0 != 3) + echo "go (fileroot) (catdir)" + break + end + + local i fileroot rootdir catdir + + $fileroot = $1 + $catdir = $2 + + # run relphot on the db and check that the images now match the expected values + tapEXEC relphot -images g,r -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir -D STAR_TOOFEW 1 -statmode WT_MEAN -cloud-limit 0.5 -update -nloop 12 + + # reset the per-filter ZPT_OFF:$filter values + $ZPT_OFF:$FILTER = NAN + + for i 0 mjd_uc[] + # ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $FILTER relphot_uc + ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $FILTER relphot_uc + end +end + +macro ckexposure + if ($0 != 8) + echo "ckexposure (catdir) (mjd) (zpt) (exptime) (airmass) (filter) (mode)" + echo " mode == raw or corr" + break + end + + local CATDIR MJD_IMAGE ZPT_REAL ZPT_NOMINAL EXPTIME FILTER AIRMASS + + $CATDIR = $1 + $MJD_IMAGE = $2 + $ZPT_REAL = $3 + $EXPTIME = $4 + $AIRMASS = $5 + $FILTER = $6 + $FUNC_MODE = $7 + + # XXX need a function to extract the nominal zpt for a given filter / photcode from the db + $ZPT_NOMINAL = $zpt_nominal:$FILTER + $KLAM_NOMINAL = $klam_nominal:$FILTER + + $TIMEFORMAT = mjd + $TIMEREF = 0.0 + + catdir $CATDIR + + # assume we still have stars_ra, stars_dec, stars_mag in hand + if ($PLOT) + dev -n 0 + region $RA_CENTER $DEC_CENTER 0.2 + images + pmeasure -all -m 15 20 + end + + # RA_MIN, etc are defined by mkstars + skyregion $RA_MIN $RA_MAX $DEC_MIN $DEC_MAX + + ## calculate the average zero point offset + if (isnan($ZPT_OFF:$FILTER)) + avextract ra dec $FILTER + match2d -closest ra dec stars_ra stars_dec 0.001 -index1 index1 -index2 index2 + reindex stars_tru_m = stars_mag using index1 + set dm_tru = $FILTER - stars_tru_m + vstat -q dm_tru + $ZPT_OFF:$FILTER = $MEDIAN + end + + mextract ra dec mag xccd yccd dbflags photflags time mag:sys mag:rel mag:inst mag:ave where (abs(time - $MJD_IMAGE) < 0.0001) + + match2d -closest ra dec stars_ra stars_dec 0.001 -index1 index1 -index2 index2 + + reindex stars_ra_m = stars_ra using index1 + reindex stars_dec_m = stars_dec using index1 + reindex stars_mag_m = stars_mag using index1 + set dr = 3600*(ra - stars_ra_m) + set dd = 3600*(dec - stars_dec_m) + set dm = mag - stars_mag_m + + set good_photom = (not(photflags & 0x02)) + foreach field dm xccd yccd + subset $field\_good = $field if good_photom + end + + if ($PLOT) + dev -n 1 + lim mag dm; clear; box; plot mag dm; plot -c red -pt 7 mag dm where (photflags & 0x00000002) + dev -n 2 + lim xccd dm; clear; box; plot xccd dm; plot -c red -pt 7 xccd dm where (photflags & 0x00000002) + dev -n 3 + lim ra dm; clear; box; plot ra dm; plot -c red -pt 7 ra dm where (photflags & 0x00000002) + echo "red points are known bad photometry" + end + + $ZPT_REAL_NORM = $ZPT_REAL - 2.5*log($EXPTIME) - $KLAM_NOMINAL*($AIRMASS - 1.0) + + if ("$FUNC_MODE" == "raw") + set dm_adjust = dm_good - $ZPT_NOMINAL + $ZPT_REAL_NORM + + vstat -q dm_adjust + sprintf line "addstar raw ZP %8.2f $FILTER %6.3f %6.3f %6.3f (%7.4f): %7.4f %7.4f : %6.3f" $MJD_IMAGE $ZPT_REAL $ZPT_REAL_NORM $ZPT_NOMINAL {$ZPT_REAL_NORM - $ZPT_NOMINAL} $MEAN $SIGMA dm_adjust[] + tapOK {abs($MEAN) < 0.005} "RAW : $line" + return + end + + if (("$FUNC_MODE" == "relphot_uc") || ("$FUNC_MODE" == "relphot_nc")) + + ## relphot does not set an absolute zero point, so we need to find that and subtract + set dm_adjust = dm_good - $ZPT_OFF:$FILTER + + vstat -q dm_adjust + sprintf line "relphot ZP %8.2f $FILTER %6.3f %6.3f %6.3f (%7.4f) : %7.4f %7.4f : %3d" $MJD_IMAGE $ZPT_REAL $ZPT_REAL_NORM $ZPT_NOMINAL {$ZPT_REAL_NORM - $ZPT_NOMINAL} $MEAN $SIGMA dm_adjust[] + + if ("$FUNC_MODE" == "relphot_uc") + tapOK {abs($MEAN) < 0.001} "MEAN : $line" + tapOK {abs($SIGMA) < 0.001} "SIGMA: $line" + end + if ("$FUNC_MODE" == "relphot_nc") + # XXX NOTE relaxed condition : bad photometry points are partly clipped but still cause a bias + tapOK {abs($MEAN) < 0.020} "MEAN : $line" + tapOK {abs($SIGMA) < 0.020} "SIGMA: $line" + end + + if ($PLOT) + cursor 1 + if ("$KEY" == "Escape") break + end + return + end + + echo "ERROR: unknown mode for ckexposure : $FUNC_MODE" +end + +macro mkexposure + if ($0 != 10) + echo "mkexposure (catdir) (mjd) (zpt) (exptime) (airmass) (filter) (fileroot) (ra) (dec)" + break + end + + local ix iy date time datetime ra dec ROOT RAo DECo ZPT MJD FILTER EXPTIME AIRMASS CATDIR cmd + + $CATDIR = $1 + $MJD = $2 + $ZPT = $3 + $EXPTIME = $4 + $AIRMASS = $5 + $FILTER = $6 + $ROOT = $7 + $RAo = $8 + $DECo = $9 + + $TIMEFORMAT = mjd + $TIMEREF = 0.0 + + # XX ctimes -ref $MJD -var datetime + # XX substr $datetime 0 10 date + # XX substr $datetime 11 8 time + + for ix 0 $NCHIP_X + for iy 0 $NCHIP_Y + $dx = $CHIP_DX * $PLATE_SCALE * ($ix - 0.5*$NCHIP_X + 0.5) + $dy = $CHIP_DY * $PLATE_SCALE * ($iy - 0.5*$NCHIP_Y + 0.5) + $ra = $RAo - $dx / 3600.0 / dcos($DECo) + $dec = $DECo - $dy / 3600.0 + # echo $ra $dec $dx $dy + mkinput $ROOT.$ix.$iy.in.txt $ra $dec $ZPT + + # ra,dec is the center of this chip + local options + $options = -mjd $MJD + $options = $options -radec $ra $dec + + #$options = $options -type PS1_V2 + #$options = $options -type PS1_V3 + #$options = $options -type PS1_V4 + $options = $options -type PS1_V5 + + $options = $options -coords + $options = $options -photcode GPC1.$FILTER.XY$ix\$iy + $options = $options -no-noise + $options = $options -flags + $options = $options -bad-psfqf-frac $DEFECT_FRAC + $options = $options -size $CHIP_DX $CHIP_DY + $options = $options -crpix {0.5*$CHIP_DX} {0.5*$CHIP_DY} + $options = $options -airmass $AIRMASS + $options = $options -exptime $EXPTIME + + tapEXEC mkcmf $ROOT.$ix.$iy.in.txt $ROOT.$ix.$iy.cmf $options + + # the fake images have inconsistent ra,dec and airmass,sidtime values + tapEXEC addstar $ROOT.$ix.$iy.cmf -D CATDIR $CATDIR -D CAMERA gpc1 -quick-airmass + end + end +end + +# make a simple input file for mkcmf +macro mkinput + if ($0 != 5) + echo "mkinput (filename) (ra) (dec) (zpt)" + break + end + + # chip coordinate of the stars + # ra,dec is at the chip center, which corresponds to pixel (CHIP_DX,CHIP_DY)/2 + set stars_X = (stars_ra - $2) * 3600.0 * dcos($3) / $PLATE_SCALE + $CHIP_DX*0.5 + set stars_Y = (stars_dec - $3) * 3600.0 / $PLATE_SCALE + $CHIP_DY*0.5 + set stars_M = (stars_mag - $4) + + set valid = (stars_X > 0) && (stars_X < $CHIP_DX) && (stars_Y > 0) && (stars_Y < $CHIP_DX) + subset stars_x = stars_X if valid + subset stars_y = stars_Y if valid + subset stars_m = stars_M if valid + subset stars_r = stars_ra if valid + subset stars_d = stars_dec if valid + subset stars_mr = stars_mag if valid + + # if we want to modify the observed magnitude based on the camera locations, add code here: + + # cx,cy can have values of 0,1 (since NCELL_X, NCELL_Y = 2) + set stars_cx = int($NCELL_X * stars_x / $CHIP_DX) + set stars_cy = int($NCELL_Y * stars_y / $CHIP_DY) + + # if I've done this right the flat-field offsets are + # (0,0) = 0.02, (0,1) = -0.02, (1,0) = 0.03, (1,1) = -0.03 + set stars_cm = (stars_cx == 0) ? ((stars_cy == 0) ? zero(stars_cx) + 0.02 : zero(stars_cx) - 0.02) : ((stars_cy == 0) ? zero(stars_cx) + 0.03 : zero(stars_cx) - 0.03) + set stars_mi = stars_m + stars_cm + + # create bad detections / detections to ignore or avoid, save to a file? + + # create 3 types of things to ignore / avoid + # * badPhot (deviate the photometry) + # * badFlag (set a bad flag bit) + # * poorFlag (set a bad flag bit) + set badPhot = rnd(stars_r) < $DEFECT_FRAC + set badFlag = rnd(stars_r) < $DEFECT_FRAC + set poorFlag = rnd(stars_r) < $DEFECT_FRAC + + # create deviant photometry (1%) + # I'm setting the flag value 0x02 (EXTMODEL) if the object is deviated. This is a simple way to pass the info + # note that outlier rejection only occurs for Nmeas > 5 + set stars_mi = stars_mi + 0.2*badPhot + set stars_fl = zero(stars_mi) + 0x00000080*badFlag + 0x00000100*poorFlag + 0x00000002*badPhot + + # bad flags that we sould respect: + # PM_SOURCE_MODE_FAIL = 0x00000008, ///< Fit (non-linear) failed (non-converge, off-edge, run to zero) + # PM_SOURCE_MODE_SATSTAR = 0x00000080, ///< Source model peak is above saturation + # PM_SOURCE_MODE_BADPSF = 0x00000400, ///< Failed to get good estimate of object's PSF + # PM_SOURCE_MODE_DEFECT = 0x00000800, ///< Source is thought to be a defect + # PM_SOURCE_MODE_SATURATED = 0x00001000, ///< Source is thought to be saturated pixels (bleed trail) + # PM_SOURCE_MODE_CR_LIMIT = 0x00002000, ///< Source has crNsigma above limit + # PM_SOURCE_MODE_MOMENTS_FAILURE = 0x00008000, ///< could not measure the moments + # PM_SOURCE_MODE_SKY_FAILURE = 0x00010000, ///< could not measure the local sky + # PM_SOURCE_MODE_SKYVAR_FAILURE = 0x00020000, ///< could not measure the local sky variance + # PM_SOURCE_MODE_SIZE_SKIPPED = 0x10000000, ///< size could not be determined + + # poor flags that we should respect: + # PM_SOURCE_MODE_POOR = 0x00000010, ///< Fit succeeds, but low-SN, high-Chisq, or large (for PSF -- drop?) + # PM_SOURCE_MODE_PAIR = 0x00000020, ///< Source fitted with a double psf + # PM_SOURCE_MODE_BLEND = 0x00000100, ///< Source is a blend with other sources + # PM_SOURCE_MODE_BELOW_MOMENTS_SN = 0x00040000, ///< moments not measured due to low S/N + # PM_SOURCE_MODE_BLEND_FIT = 0x00400000, ///< source was fitted as a blend + # PM_SOURCE_MODE_ON_SPIKE = 0x20000000, ///< peak lands on diffraction spike + # PM_SOURCE_MODE_ON_GHOST = 0x40000000, ///< peak lands on ghost or glint + # PM_SOURCE_MODE_OFF_CHIP = 0x80000000, ///< peak lands off edge of chip + + write -f "%10.6f %10.6f %8.4f %6.1f %6.1f %8.4f 0x%x" $1 stars_r stars_d stars_mr stars_x stars_y stars_mi stars_fl +end + +# generate a set of stars with raw RA, DEC, MAG values to use for all of the output cmf files +macro mkstars + if ($0 != 2) + echo "USAGE: mkstars (Nstars)" + break + end + + local RA_RANGE DEC_RANGE + + # the images are oriented along N-S, E-W lines + # we make a dithered grid of images in each night + + # size of region of interest in linear degrees + $RA_RANGE = 2.5 * $CHIP_DX * $NCHIP_X * $PLATE_SCALE / 3600 + $DEC_RANGE = 2.5 * $CHIP_DY * $NCHIP_Y * $PLATE_SCALE / 3600 + + $RA_MIN = $RA_CENTER - 0.5*$RA_RANGE/dcos($DEC_CENTER) + $RA_MAX = $RA_CENTER + 0.5*$RA_RANGE/dcos($DEC_CENTER) + $DEC_MIN = $DEC_CENTER - 0.5*$DEC_RANGE + $DEC_MAX = $DEC_CENTER + 0.5*$DEC_RANGE + + create tmp 0 $1 + set stars_ra = $RA_CENTER + $RA_RANGE * (rnd(tmp) - 0.5) / dcos ($DEC_CENTER) + set stars_dec = $DEC_CENTER + $DEC_RANGE * (rnd(tmp) - 0.5) + set stars_mag = 10.0 + 10.0 * rnd(tmp) + set stars_dmag = 0.02 + 0.1 * rnd(tmp) +end + +macro plot.test.files + if ($0 != 2) + echo "USAGE: plot.test.files (N)" + break + end + + clear -s + resize 2000 1200 + + section a0 0.00 0.50 0.25 0.50 + section a1 0.25 0.50 0.25 0.50 + section a2 0.50 0.50 0.25 0.50 + section a3 0.75 0.50 0.25 0.50 + + section a4 0.00 0.00 0.25 0.50 + section a5 0.25 0.00 0.25 0.50 + section a6 0.50 0.00 0.25 0.50 + section a7 0.75 0.00 0.25 0.50 + + label -fn helvetica 12 + $BSTYLE = -tickpad 0.2 +xpad 0.2 +ypad 0.2 -ypad 5 + + for i 0 8 + sprintf name tgrp.%05d.%02d.dat $i $1 + data $name + exec ls -l --full-time $name + read Msys 4 Mrel 5 Mcal 7 Mmos 8 Mgrp 9 airmass 11 dMag 12 + + section a$i + lim Msys dMag; box $BSTYLE; plot Msys dMag -pt cir -op 0.2 -c blue + end +end + +macro plot.tcat.data + if ($0 != 2) + echo "USAGE: plot.tcat.data (Nlast)" + break + end + + $BSTYLE = -tickpad 0.2 +xpad 0.2 +ypad 0.2 -ypad 5 + + for i 0 $1 + sprintf name tcat.%02d.dat $i + data $name + dev -n $i + resize 900 500 + clear; + box $BSTYLE + read Msys 4 Mrel 5 Mcal 7 Mmos 8 Mgrp 9 airmass 11 dMag 12; + lim Msys dMag; plot Msys dMag -pt cir -op 0.2 -c blue + echo $name + end +end Index: trunk/Ohana/src/relphot/test/relphot.samples.cmd =================================================================== --- trunk/Ohana/src/relphot/test/relphot.samples.cmd (revision 41647) +++ trunk/Ohana/src/relphot/test/relphot.samples.cmd (revision 41647) @@ -0,0 +1,6 @@ +relphot -tgroup-fit-airmass -tgroups tgroups.dat -images g -v -region 9 11 19 21 -D CATDIR t0.catdir.v2 -D STAR_TOOFEW 1 -statmode WT_MEAN -cloud-limit 0.5 -update -nloop 12 +-images g -v -region 9 11 19 21 -D CATDIR t1.catdir.v4 -D STAR_TOOFEW 1 -statmode MEAN -cloud-limit 0.5 -update -nloop 12 + +relphot -tgroup-fit-airmass -tgroups tgroups.dat -mosaic -imfreeze -images g -v -region 9 11 19 21 -D CATDIR t0.catdir.v2 -D STAR_TOOFEW 1 -statmode WT_MEAN -cloud-limit 0.5 -update -nloop 12 +relphot -tgroup-fit-airmass -tgroups tgroups.dat -mosaic -images g -v -region 9 11 19 21 -D CATDIR t2.catdir.v0 -D STAR_TOOFEW 1 -statmode MEAN -cloud-limit 0.5 -update -nloop 12 +relphot -ref-weight -mosaic -images g -v -region 9 11 19 21 -D CATDIR t2.catdir.v0 -D STAR_TOOFEW 1 -statmode MEAN -cloud-limit 0.5 -nloop 12 Index: trunk/Ohana/src/relphot/test/sample.gdb =================================================================== --- trunk/Ohana/src/relphot/test/sample.gdb (revision 41647) +++ trunk/Ohana/src/relphot/test/sample.gdb (revision 41647) @@ -0,0 +1,9 @@ + +set $i = 0 + +while ($i < Nmosaic) +# printf "%f %f %f\n", dataset->alldata->xVector[$i], dataset->alldata->yVector[$i], dataset->alldata->dyVector[$i] +# printf "%f %f\n", results->psfData[Nsec].flxlist[$i], results->psfData[Nsec].errlist[$i] + printf "%d : %f %f\n", $i, mlist[$i], slist[$i] + set $i = $i + 1 +end