Index: /trunk/Ohana/src/relphot/include/relphot.h =================================================================== --- /trunk/Ohana/src/relphot/include/relphot.h (revision 38465) +++ /trunk/Ohana/src/relphot/include/relphot.h (revision 38466) @@ -96,4 +96,13 @@ typedef struct { + double *flxlist; // list of measure.mag values for a given star + double *errlist; // mag errors for a star + double *wgtlist; // weights to use for mean mags + int *ranking; // weights to use for mean mags + int *measSeq; // weights to use for mean mags + int Nlist; +} StatDataSet; + +typedef struct { int Nfew; int Ncode; @@ -105,22 +114,37 @@ // NOTE: the following arrays are (possibly) pre-allocated and carried down to each - // thread. The first 3 (list, dlist, wlist) are used in all relphot analyses; the - // others are only used on the final output steps. - - double *Mpsflist; // list of measure.mag values for a given star - double *dpsflist; // mag errors for a star - double *wpsflist; // weights to use for mean mags - - double *Maplist; // ap mags for a star - double *daplist; // ap mags for a star - double *waplist; // ap mags for a star - - double *Mkronlist; // kron mags for a star - double *dkronlist; // kron mag errors - double *wkronlist; // kron mag errors + // thread. The psfData are used in all relphot analyses; the others are only used on + // the final output steps. + + int Nsecfilt; + + StatDataSet *psfData; // one is allocated for each primary (average) photcode + StatDataSet *aperData; + StatDataSet *kronData; double *psfqf_list; // psfqf for all filters double *psfqfperf_list; // psfqfperf for all filters double *stargal_list; // stargal for all filters + + int *havePS1; // this secfilt has synthetic mags + int *haveSYN; // this secfilt has synthetic mags + int *measSYN; // this measurement is the synthetic mag for this secfilt + int *needSYN; // this secfilt mag should use synthetic mags + float *minSYN; // minimum synthetic mag below which synthetic should be forced + + float *psfQfMax; // max psfQf value for this secfilt + float *psfQfPerfMax; // max psfQfperf value for this secfilt + + int *Nmeas; // count of PS1 exposure (chip) measurements for this secfilt + int *NmeasGood; // count of PS1 exposure (chip) measurements for this secfilt + int *Next; // count of PS1 exposure (chip) measurements for this secfilt + int *NexpPS1; // count of PS1 exposure (chip) measurements for this secfilt + int *haveUbercal; // does this secfilt have any ubercal data? + + int *tessID; // tess,proj,skycell to use for warp and diff analysis + int *projID; + int *skycellID; + + float *minUbercalDist; StatType psfstats; @@ -375,4 +399,5 @@ void liststats_setmode PROTO((StatType *stats, char *strmode)); int liststats PROTO((double *value, double *dvalue, double *wvalue, int N, StatType *stats)); +int liststats_init PROTO((StatType *stats)); Catalog *load_catalogs PROTO((SkyList *skylist, int *Ncatalog, int hostID, char *hostpath, char *syncfile)); Catalog *load_catalogs_parallel PROTO((SkyList *sky, int *Ncatalog, char *syncfile)); @@ -412,5 +437,5 @@ 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, int pass, FlatCorrectionTable *flatcorr)); +int setMrelOutput PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr)); int setMave PROTO((Catalog *catalog, int Ncatalog)); void set_ZP PROTO((double ZERO)); @@ -474,12 +499,18 @@ int print_measure_set_alt (Average *average, SecFilt *secfilt, Measure *measure); -int setMrel_catalog_alt (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt); -int setMrelAverageExposure (off_t meas, int cat, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT, off_t *found); -int setMrelAverageStack (off_t measureOffset, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t Nmeasure); -int setMrelAverageForcedWarp (off_t measureOffset, int cat, int pass, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t *found, off_t Nmeasure); +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); int setGlobalObjStats (Average *average, Measure *measure); void SetMrelInfoInit (SetMrelInfo *results, int allocLists); +void SetMrelInfoFree (SetMrelInfo *results); +void SetMrelInfoReset (SetMrelInfo *results); +void SetMrelInfoResetObject (SetMrelInfo *results); + +void StatDataSetFree (StatDataSet *dataset, int Nsecfilt); +StatDataSet *StatDataSetAlloc (int Nsecfilt, int Nmax); int init_synthetic_mags (); @@ -533,7 +564,11 @@ int isGPC1synth (int photcode); int whichGPC1filter (int photcode); - - -SynthZeroPoints *SynthZeroPointsLoad (char *filename); +int is2MASS (int photcode); +int isTYCHO (int photcode); + +int magStatsByRanking (StatDataSet *dataset, StatType *stats); + +int SynthZeroPointsLoad (char *filename); +SynthZeroPoints *SynthZeroPointsGet (); int relphot_synthphot (int hostID, char *hostpath); @@ -542,2 +577,7 @@ int relphot_synthphot_catalog (Catalog *catalog, SynthZeroPoints *zpts); int relphot_synthphot_average (Average *average, SecFilt *secfilt, Measure *measure, SynthZeroPoints *zpts); + +void setMeasureRank (Catalog *catalog); +int getImageFlags (off_t meas, int cat); +int getMosaicFlags (off_t meas, int cat); + Index: /trunk/Ohana/src/relphot/src/ImageOps.c =================================================================== --- /trunk/Ohana/src/relphot/src/ImageOps.c (revision 38465) +++ /trunk/Ohana/src/relphot/src/ImageOps.c (revision 38466) @@ -397,4 +397,15 @@ i = MeasureToImage[cat][meas]; return (i); +} + +int getImageFlags (off_t meas, int cat) { + + off_t i; + + if (!MeasureToImage) return 0; + + i = MeasureToImage[cat][meas]; + + return (image[i].flags); } @@ -681,4 +692,8 @@ list[N] = Msys - Mrel - Mmos - Mgrid + Mflat; dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); + if (fabs(list[N]) > 0.03) { + // fprintf (stderr, "deviant\n"); + } + if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) { Mlist[Nbright] = list[N]; @@ -711,5 +726,5 @@ if (PLOTSTUFF) { - fprintf (stderr, "Mcal for : %s\n", image[i].name); + fprintf (stderr, "Mcal for : %s : %7.4f %7.4f\n", image[i].name, image[i].Mcal, image[i].dMcal); plot_setMcal (list, N, &stats, CLOUD_TOLERANCE); } Index: /trunk/Ohana/src/relphot/src/MosaicOps.c =================================================================== --- /trunk/Ohana/src/relphot/src/MosaicOps.c (revision 38465) +++ /trunk/Ohana/src/relphot/src/MosaicOps.c (revision 38466) @@ -882,4 +882,15 @@ } +int getMosaicFlags (off_t meas, int cat) { + + if (!MOSAIC_ZEROPT) return (0); + + // unassigned measurements belong to simple chips + int i = MeasureToMosaic[cat][meas]; + if (i == -1) return (0); + + return (mosaic[i].flags); +} + typedef struct { int Nfew; Index: /trunk/Ohana/src/relphot/src/StarOps.c =================================================================== --- /trunk/Ohana/src/relphot/src/StarOps.c (revision 38465) +++ /trunk/Ohana/src/relphot/src/StarOps.c (revision 38466) @@ -1,5 +1,3 @@ # include "relphot.h" - -static int Nmax; enum {THREAD_RUN, THREAD_DONE}; @@ -20,6 +18,8 @@ int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure); -// we want to allocate the stats list,dlist arrays only once (or once per thread). -// this function finds the largest array so we can allocate that max size when needed +// we want to allocate the StatDataSet arrays only once (or once per thread). this +// function finds the largest value of Nmeasure so we can allocate that max size when +// needed +static int Nmax; void initMrel (Catalog *catalog, int Ncatalog) { @@ -54,5 +54,35 @@ } -void SetMrelInfoInit (SetMrelInfo *results, int allocLists) { +StatDataSet *StatDataSetAlloc (int Nsecfilt, int Nmax) { + + int i; + + StatDataSet *dataset = NULL; + ALLOCATE (dataset, StatDataSet, Nsecfilt); + for (i = 0; i < Nsecfilt; i++) { + ALLOCATE (dataset[i].flxlist, double, Nmax); + ALLOCATE (dataset[i].wgtlist, double, Nmax); + ALLOCATE (dataset[i].errlist, double, Nmax); + ALLOCATE (dataset[i].ranking, int, Nmax); + ALLOCATE (dataset[i].measSeq, int, Nmax); + } + return dataset; +} + +void StatDataSetFree (StatDataSet *dataset, int Nsecfilt) { + + int i; + + for (i = 0; i < Nsecfilt; i++) { + FREE (dataset[i].flxlist); + FREE (dataset[i].wgtlist); + FREE (dataset[i].errlist); + FREE (dataset[i].ranking); + FREE (dataset[i].measSeq); + } + FREE (dataset); +} + +void SetMrelInfoReset (SetMrelInfo *results) { results->Nfew = 0; results->Ncode = 0; @@ -62,15 +92,147 @@ results->Nmos = 0; results->Ngrid = 0; +} + +static int NsynthNames = 6; +static float synthMaxMags[] = {13.5, 13.5, 13.5, 13.0, 12.0, 13.5}; +static char *synthNames[] = { "g", "r", "i", "z", "y", "w"}; + +void SetMrelInfoResetObject (SetMrelInfo *results) { + + int i; + for (i = 0; i < results->Nsecfilt; i++) { + results->havePS1[i] = FALSE; + results->haveSYN[i] = FALSE; + results->needSYN[i] = FALSE; + results->measSYN[i] = -1; + + results->psfQfMax[i] = NAN; + results->psfQfPerfMax[i] = NAN; + + results->Nmeas[i] = 0; + results->NmeasGood[i] = 0; + results->Next[i] = 0; + results->NexpPS1[i] = 0; + results->havePS1[i] = 0; + results->haveUbercal[i] = 0; + + results->tessID[i] = 0; + results->projID[i] = 0; + results->skycellID[i] = 0; + + results->minUbercalDist[i] = 1000; + + results->psfData[i].Nlist = 0; + results->aperData[i].Nlist = 0; + results->kronData[i].Nlist = 0; + } +} + +void SetMrelInfoInit (SetMrelInfo *results, int allocLists) { + + int i; + SetMrelInfoReset (results); if (allocLists) { - ALLOCATE (results->Mpsflist, double, Nmax); - ALLOCATE (results->dpsflist, double, Nmax); - ALLOCATE (results->wpsflist, double, Nmax); + results->Nsecfilt = GetPhotcodeNsecfilt (); + + results->psfData = StatDataSetAlloc (results->Nsecfilt, Nmax); + results->aperData = StatDataSetAlloc (results->Nsecfilt, Nmax); + results->kronData = StatDataSetAlloc (results->Nsecfilt, Nmax); + + ALLOCATE (results->psfqf_list, double, Nmax); + ALLOCATE (results->psfqfperf_list, double, Nmax); + ALLOCATE (results->stargal_list, double, Nmax); + + ALLOCATE (results->havePS1, int, results->Nsecfilt); + ALLOCATE (results->haveSYN, int, results->Nsecfilt); + ALLOCATE (results->needSYN, int, results->Nsecfilt); + ALLOCATE (results->measSYN, int, results->Nsecfilt); + ALLOCATE (results->minSYN, float, results->Nsecfilt); + + ALLOCATE (results->psfQfMax, float, results->Nsecfilt); + ALLOCATE (results->psfQfPerfMax, float, results->Nsecfilt); + + ALLOCATE (results->Nmeas, int, results->Nsecfilt); + ALLOCATE (results->NmeasGood, int, results->Nsecfilt); + ALLOCATE (results->Next, int, results->Nsecfilt); + ALLOCATE (results->NexpPS1, int, results->Nsecfilt); + ALLOCATE (results->haveUbercal, int, results->Nsecfilt); + + ALLOCATE (results->tessID, int, results->Nsecfilt); + ALLOCATE (results->projID, int, results->Nsecfilt); + ALLOCATE (results->skycellID, int, results->Nsecfilt); + + ALLOCATE (results->minUbercalDist, float, results->Nsecfilt); + + SetMrelInfoResetObject (results); + + for (i = 0; i < NsynthNames; i++) { + int photcode = GetPhotcodeCodebyName (synthNames[i]); + if (!photcode) continue; + int Nsec = GetPhotcodeNsec (photcode); + if (Nsec < 0) continue; + results->minSYN[Nsec] = synthMaxMags[i]; + } + } else { + results->Nsecfilt = 0; + results->psfData = NULL; + results->aperData = NULL; + results->kronData = NULL; + + results->psfqf_list = NULL; + results->psfqfperf_list = NULL; + results->stargal_list = NULL; + + results->havePS1 = NULL; + results->haveSYN = NULL; + results->needSYN = NULL; + results->measSYN = NULL; + results->minSYN = NULL; + + results->psfQfMax = NULL; + results->psfQfPerfMax = NULL; + + results->Nmeas = NULL; + results->NmeasGood = NULL; + results->Next = NULL; + results->NexpPS1 = NULL; + results->haveUbercal = NULL; + + results->tessID = NULL; + results->projID = NULL; + results->skycellID = NULL; + + results->minUbercalDist = NULL; } } void SetMrelInfoFree (SetMrelInfo *results) { - free (results->Mpsflist); - free (results->dpsflist); - free (results->wpsflist); + StatDataSetFree (results->psfData, results->Nsecfilt); + StatDataSetFree (results->aperData, results->Nsecfilt); + StatDataSetFree (results->kronData, results->Nsecfilt); + FREE (results->psfqf_list); + FREE (results->psfqfperf_list); + FREE (results->stargal_list); + + FREE (results->havePS1); + FREE (results->haveSYN); + FREE (results->needSYN); + FREE (results->measSYN); + FREE (results->minSYN ); + + FREE (results->psfQfMax); + FREE (results->psfQfPerfMax); + + FREE (results->Nmeas); + FREE (results->NmeasGood); + FREE (results->Next); + FREE (results->NexpPS1); + FREE (results->haveUbercal); + + FREE (results->tessID); + FREE (results->projID); + FREE (results->skycellID); + + FREE (results->minUbercalDist); } @@ -128,5 +290,5 @@ for (i = 0; i < Ncatalog; i++) { - setMrel_catalog_alt (catalog, i, -1, flatcorr, &results, Nsecfilt); + setMrel_catalog_alt (catalog, i, FALSE, flatcorr, &results, Nsecfilt); SetMrelInfoAccum (&summary, &results); } @@ -138,11 +300,9 @@ } -int setMrelOutput (Catalog *catalog, int Ncatalog, int pass, FlatCorrectionTable *flatcorr) { +// setMrelOutput is NOT threaded because we read/write catalogs one at a time (for now) +int setMrelOutput (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) { int i; - // fprintf (stderr, "this version of setMrel is invalid for now\n"); - // exit (1); - int Nsecfilt = GetPhotcodeNsecfilt (); @@ -151,18 +311,6 @@ SetMrelInfoInit (&results, TRUE); // allocates results->list,dlist,wlist - ALLOCATE (results.Maplist, double, Nmax); - ALLOCATE (results.daplist, double, Nmax); - ALLOCATE (results.waplist, double, Nmax); - - ALLOCATE (results.Mkronlist, double, Nmax); - ALLOCATE (results.dkronlist, double, Nmax); - ALLOCATE (results.wkronlist, double, Nmax); - - ALLOCATE (results.psfqf_list, double, Nmax); - ALLOCATE (results.psfqfperf_list, double, Nmax); - ALLOCATE (results.stargal_list, double, Nmax); - for (i = 0; i < Ncatalog; i++) { - setMrel_catalog_alt (catalog, i, pass, flatcorr, &results, Nsecfilt); + setMrel_catalog_alt (catalog, i, TRUE, flatcorr, &results, Nsecfilt); SetMrelInfoAccum (&summary, &results); } @@ -170,15 +318,4 @@ SetMrelInfoFree (&results); - free (results.Maplist); - free (results.daplist); - free (results.waplist); - - free (results.Mkronlist); - free (results.dkronlist); - free (results.wkronlist); - - free (results.psfqf_list); - free (results.psfqfperf_list); - free (results.stargal_list); return (TRUE); } @@ -261,5 +398,4 @@ while (1) { - off_t i = getNextCatalogForThread(threadinfo->Ncatalog); if (i == -1) { @@ -271,6 +407,5 @@ FlatCorrectionTable *flatcorr = threadinfo->flatcorr; - // pass == -1 for anything other than the final pass - setMrel_catalog_alt (catalog, i, -1, flatcorr, &results, Nsecfilt); + setMrel_catalog_alt (catalog, i, FALSE, flatcorr, &results, Nsecfilt); SetMrelInfoAccum (&threadinfo->summary, &results); } @@ -1142,4 +1277,10 @@ } 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; + } } } Index: /trunk/Ohana/src/relphot/src/args.c =================================================================== --- /trunk/Ohana/src/relphot/src/args.c (revision 38465) +++ /trunk/Ohana/src/relphot/src/args.c (revision 38466) @@ -324,4 +324,6 @@ } + // if -synthphot is chosen, an artificial w-band measurement is generated for each object + // with r and i photometry. this was used to apply the ubercal photometry to w-band images. SyntheticPhotometry = FALSE; if ((N = get_argument (argc, argv, "-synthphot"))) { @@ -365,5 +367,15 @@ mode = UPDATE_AVERAGES; } + + // if -synthphot-zpts is supplied, the map of corrections needed for the synthetic photometry is supplied + // and used to tie down the synthetic magnitudes + // NOTE: SYNTH_ZERO_POINTS is used for the same file in both -synthphot-zpts and -synthphot_means SYNTH_ZERO_POINTS = NULL; + if ((N = get_argument (argc, argv, "-synthphot-zpts"))) { + remove_argument (N, &argc, argv); + myAssert (N < argc, "missing argument to -synthphot-zpts"); + SYNTH_ZERO_POINTS = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } if ((N = get_argument (argc, argv, "-synthphot_means"))) { mode = SYNTH_PHOT; @@ -476,4 +488,10 @@ } SYNTH_ZERO_POINTS = NULL; + if ((N = get_argument (argc, argv, "-synthphot-zpts"))) { + remove_argument (N, &argc, argv); + myAssert (N < argc, "missing argument to -synthphot-zpts"); + SYNTH_ZERO_POINTS = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } if ((N = get_argument (argc, argv, "-synthphot_means"))) { MODE = MODE_SYNTH_PHOT; Index: /trunk/Ohana/src/relphot/src/bcatalog.c =================================================================== --- /trunk/Ohana/src/relphot/src/bcatalog.c (revision 38465) +++ /trunk/Ohana/src/relphot/src/bcatalog.c (revision 38466) @@ -37,9 +37,15 @@ // 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_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++) { @@ -59,18 +65,4 @@ int Ns; -# if (0) - DVOSecfiltFlags secfiltBits = - ID_SECF_STAR_FEW | - ID_SECF_STAR_POOR | - ID_PHOTOM_PASS_0 | - ID_PHOTOM_PASS_1 | - ID_PHOTOM_PASS_2 | - ID_PHOTOM_PASS_3 | - ID_PHOTOM_PASS_4 | - ID_SECF_USE_SYNTH | - ID_SECF_USE_UBERCAL | - ID_SECF_OBJ_EXT; -# endif - for (Ns = 0; Ns < Nphotcodes; Ns++) { @@ -78,18 +70,5 @@ int Nsec = GetPhotcodeNsec(thisCode); - // fprintf (stderr, "really use dvo_secfilt_init?"); - dvo_secfilt_init (&subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec]); - - // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M = NAN; - // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Map = NAN; - // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].dM = NAN; - // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Mstdev = NAN_S_SHORT; - // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Mchisq = NAN; - // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M_20 = NAN_S_SHORT; - // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M_80 = NAN_S_SHORT; - // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Ncode = 0; - // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Nused = 0; - // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].flags &= ~secfiltBits; - // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].ubercalDist = 1000; + dvo_secfilt_init (&subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec], SECFILT_RESET_ALL); } } Index: /trunk/Ohana/src/relphot/src/extra.c =================================================================== --- /trunk/Ohana/src/relphot/src/extra.c (revision 38465) +++ /trunk/Ohana/src/relphot/src/extra.c (revision 38466) @@ -70,2 +70,15 @@ return FALSE; } + +int is2MASS (int photcode) { + + if ((photcode >= 2011) && (photcode <= 2013)) return TRUE; + return FALSE; +} + +int isTYCHO (int photcode) { + + if ((photcode == 2020) || (photcode == 2021)) return TRUE; + return FALSE; +} + Index: /trunk/Ohana/src/relphot/src/global_stats.c =================================================================== --- /trunk/Ohana/src/relphot/src/global_stats.c (revision 38465) +++ /trunk/Ohana/src/relphot/src/global_stats.c (revision 38466) @@ -4,4 +4,16 @@ StatType stN, stX, stS, imN, imX, imM, imD, msM, msX, msN, msD; + + liststats_init (&stN); + liststats_init (&stX); + liststats_init (&stS); + liststats_init (&imN); + liststats_init (&imX); + liststats_init (&imM); + liststats_init (&imD); + liststats_init (&msN); + liststats_init (&msX); + liststats_init (&msM); + liststats_init (&msD); // INITTIME; @@ -26,7 +38,7 @@ fprintf (stderr, " --- stats for %s ---\n", photcodes[Ns][0].name); - fprintf (stderr, "meas / star: %7.0f %7.1f %7.1f %7.0f %7.0f %6d\n", stN.median, stN.mean, stN.sigma, stN.min, stN.max, stN.Nmeas); + fprintf (stderr, "meas / star: %7.0f %7.1f %7.1f %7.0f %7.0f %6d\n", stN.median, stN.mean, stN.sigma, stN.min, stN.max, stN.Nmeas); fprintf (stderr, "dMrel star: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n", stS.median, stS.mean, stS.sigma, stS.min, stS.max, stS.Nmeas); - fprintf (stderr, "chisq star: %7.1f %7.1f %7.1f %7.1f %7.1f %6d\n", stX.median, stX.mean, stX.sigma, stX.min, stX.max, stX.Nmeas); + fprintf (stderr, "chisq star: %7.1f %7.1f %7.1f %7.1f %7.1f %6d\n", stX.median, stX.mean, stX.sigma, stX.min, stX.max, stX.Nmeas); } Index: /trunk/Ohana/src/relphot/src/liststats.c =================================================================== --- /trunk/Ohana/src/relphot/src/liststats.c (revision 38465) +++ /trunk/Ohana/src/relphot/src/liststats.c (revision 38466) @@ -16,4 +16,19 @@ } +int liststats_init (StatType *stats) { + stats->median = NAN; + stats->mean = NAN; + stats->sigma = NAN; + stats->error = NAN; + stats->chisq = NAN; + stats->min = NAN; + stats->max = NAN; + stats->Upper80 = NAN; + stats->Lower20 = NAN; + stats->total = NAN; + stats->Nmeas = 0; + return TRUE; +} + int liststats (double *value, double *dvalue, double *weight, int N, StatType *stats) { @@ -22,4 +37,6 @@ myAssert (stats->statmode != STATS_NONE, "programming error, liststats mode not set"); + + liststats_init (stats); ke = ks = dMo = 0; Index: /trunk/Ohana/src/relphot/src/plotstuff.c =================================================================== --- /trunk/Ohana/src/relphot/src/plotstuff.c (revision 38465) +++ /trunk/Ohana/src/relphot/src/plotstuff.c (revision 38466) @@ -109,4 +109,13 @@ } +void plot_wait () { + if (PLOTDELAY > 500000) { + fprintf (stdout, "press return\n"); + if (fscanf (stdin, "%*c") != 1) fprintf (stderr, "\n"); + } else { + usleep (PLOTDELAY); + } +} + // plot the vector pair to a file with name defined by the varargs format void plot_list (Graphdata *graphdata, double *xlist, double *ylist, int N, char *label, char *format, ...) { @@ -151,10 +160,5 @@ if (SAVEPLOT) JpegPlot (graphdata, 0, filename); - if (PLOTDELAY > 500000) { - fprintf (stdout, "press return\n"); - if (fscanf (stdin, "%*c") != 1) fprintf (stderr, "\n"); - } else { - usleep (PLOTDELAY); - } + plot_wait(); free (filename); } @@ -166,4 +170,5 @@ PlotVector (Npts, xlist, 0, 0, "x"); PlotVector (Npts, ylist, 1, 0, "y"); + plot_wait(); } Index: /trunk/Ohana/src/relphot/src/reload_catalogs.c =================================================================== --- /trunk/Ohana/src/relphot/src/reload_catalogs.c (revision 38465) +++ /trunk/Ohana/src/relphot/src/reload_catalogs.c (revision 38466) @@ -33,4 +33,7 @@ } + // load the ZP corrections here + if (SYNTH_ZERO_POINTS) SynthZeroPointsLoad (SYNTH_ZERO_POINTS); + if (VERBOSE) fprintf (stderr, "re-loading catalog data\n"); @@ -63,6 +66,7 @@ if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) { - fprintf (stderr, "ERROR: failure reading catalog %s\n", catalog.filename); - exit (1); + fprintf (stderr, "ERROR: failure reading catalog %s, skipping\n", catalog.filename); + continue; + // exit (1); } if (VERBOSE && (catalog.Naverage_disk == 0)) { @@ -168,4 +172,9 @@ free (BOUNDARY_TREE); BOUNDARY_TREE = tmppath; + } + if (SYNTH_ZERO_POINTS) { + char *tmppath = abspath(SYNTH_ZERO_POINTS, DVO_MAX_PATH); + free (SYNTH_ZERO_POINTS); + SYNTH_ZERO_POINTS = tmppath; } @@ -234,4 +243,5 @@ 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); } // deprecate Index: /trunk/Ohana/src/relphot/src/relphot_objects.c =================================================================== --- /trunk/Ohana/src/relphot/src/relphot_objects.c (revision 38465) +++ /trunk/Ohana/src/relphot/src/relphot_objects.c (revision 38466) @@ -30,4 +30,7 @@ } + // load the ZP corrections here + if (SYNTH_ZERO_POINTS) SynthZeroPointsLoad (SYNTH_ZERO_POINTS); + // load data from each region file, only use bright stars for (i = 0; i < skylist[0].Nregions; i++) { @@ -60,18 +63,5 @@ if (RESET) { Nsecfilt = catalog.Nsecfilt; - DVOSecfiltFlags secfiltBits = - ID_SECF_STAR_FEW | - ID_SECF_STAR_POOR | - ID_PHOTOM_PASS_0 | - ID_PHOTOM_PASS_1 | - ID_PHOTOM_PASS_2 | - ID_PHOTOM_PASS_3 | - ID_PHOTOM_PASS_4 | - ID_SECF_USE_SYNTH | - ID_SECF_USE_UBERCAL | - ID_SECF_HAS_STACK | - ID_SECF_HAS_PS1 | - ID_SECF_OBJ_EXT ; - + DVOAverageFlags astromBits = ID_STAR_FIT_AVE | @@ -89,6 +79,5 @@ catalog.average[j].stargal = NAN; // reset (will be re-calculated here) for (k = 0; k < Nsecfilt; k++) { - dvo_secfilt_init (&catalog.secfilt[j*Nsecfilt + k]); - catalog.secfilt[j*Nsecfilt + k].flags &= ~secfiltBits; // XXX: only reset photom bits: + dvo_secfilt_init (&catalog.secfilt[j*Nsecfilt + k], SECFILT_RESET_ALL); } } @@ -99,4 +88,5 @@ initMrel (&catalog, 1); setMrelFinal (&catalog, NULL, TRUE); + // XXX if we want to have options for setting warp, chip, stack independently, we need to init only the desired ones if (!UPDATE) { @@ -164,4 +154,5 @@ if (UPDATE) { strextend (&command, "-update"); } if (!KEEP_UBERCAL) { strextend (&command, "-reset-ubercal"); } + if (SYNTH_ZERO_POINTS) { snprintf (tmpline, 1024, "%s -synthphot-zpts %s", command, SYNTH_ZERO_POINTS); strcpy (command, tmpline); } // if (SET_MREL_VERSION != 1) { strextend (&command, "-set-mrel-version %d", SET_MREL_VERSION); } Index: /trunk/Ohana/src/relphot/src/relphot_synthphot.c =================================================================== --- /trunk/Ohana/src/relphot/src/relphot_synthphot.c (revision 38465) +++ /trunk/Ohana/src/relphot/src/relphot_synthphot.c (revision 38466) @@ -28,5 +28,6 @@ // load the ZP corrections here - SynthZeroPoints *zpts = SynthZeroPointsLoad (SYNTH_ZERO_POINTS); + SynthZeroPointsLoad (SYNTH_ZERO_POINTS); + SynthZeroPoints *zpts = SynthZeroPointsGet (); if (!zpts) exit (2); Index: /trunk/Ohana/src/relphot/src/setMrelCatalog.c =================================================================== --- /trunk/Ohana/src/relphot/src/setMrelCatalog.c (revision 38465) +++ /trunk/Ohana/src/relphot/src/setMrelCatalog.c (revision 38466) @@ -1,16 +1,17 @@ # include "relphot.h" +int markMeasureByRanking (StatDataSet *dataset, Measure *measure, int minrank, DVOMeasureFlags flags); + # define UBERCAL_WEIGHT 100.0 - -# define SKIP_THIS_MEAS(REASON) { \ - measureT[k].dbFlags |= ID_MEAS_SKIP_PHOTOM; \ - if (measure) { \ - measure[k].dbFlags |= ID_MEAS_SKIP_PHOTOM; \ - } \ - results->REASON ++; \ +void sort_StatDataSet (StatDataSet *dataset); + +# define SKIP_THIS_MEAS(REASON) { \ + results->REASON ++; \ continue; } # define SKIP_THIS_MEAS_STACK(REASON) { \ - measure[k].dbFlags |= ID_MEAS_SKIP_PHOTOM; \ continue; } + +# define CHECK_VALID_MAG(MAG,D_MAG) (isfinite(MAG) && isfinite(D_MAG) && (MAG > -5.0) && (MAG < 30.0)) +# define CHECK_VALID_FLUX(FLUX,D_FLUX) (isfinite(FLUX) && isfinite(D_FLUX)) static float MagToFlux (float Mag) { @@ -36,402 +37,390 @@ } -int setMrel_catalog_alt (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) { +int setMrel_catalog_alt (Catalog *catalog, int Nc, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) { off_t j; - liststats_setmode (&results->psfstats, STATMODE); - liststats_setmode (&results->apstats, STATMODE); + liststats_setmode (&results->psfstats, STATMODE); + liststats_setmode (&results->apstats, STATMODE); liststats_setmode (&results->kronstats, STATMODE); - SetMrelInfoInit (results, FALSE); // do not allocate list,dlist,wlist arrays - - int isSetMrelFinal = (pass >= 0); - - // XX char *primaryCell = NULL; - // XX if (isSetMrelFinal) { - // XX ALLOCATE (primaryCell, char, DVO_MAX_PATH); - // XX } + SetMrelInfoReset (results); // reset the countesrs for (j = 0; j < catalog[Nc].Naverage; j++) { - // we are guaranteed to have averageT, but not average - off_t m = catalog[Nc].averageT[j].measureOffset; - Average *average = catalog[Nc].average ? &catalog[Nc].average[j] : NULL; - Measure *measure = catalog[Nc].measure ? &catalog[Nc].measure[m] : NULL; - SecFilt *secfilt = catalog[Nc].secfilt ? &catalog[Nc].secfilt[j*Nsecfilt] : NULL; - setMrelAverageExposure (m, Nc, pass, flatcorr, results, average, &catalog[Nc].averageT[j], secfilt, measure, &catalog[Nc].measureT[m], &catalog[Nc].found_t[Nsecfilt*j]); + + setMrelAverageExposure (&catalog[Nc], Nc, j, Nsecfilt, isSetMrelFinal, flatcorr, results); // only apply Stack operation on setMrelFinal in first pass - if (isSetMrelFinal && (pass == 0) && !IS_DIFF_DB) { - setMrelAverageStack (m, Nc, flatcorr, results, average, secfilt, &catalog[Nc].measure[m], catalog[Nc].Nmeasure); - setGlobalObjStats (&catalog[Nc].average[j], &catalog[Nc].measure[m]); + if (isSetMrelFinal && !IS_DIFF_DB) { + setMrelAverageStack (&catalog[Nc], Nc, j, Nsecfilt, flatcorr); } // only measure force-warp mean values if issetMrelFinal (make it optional?) if (isSetMrelFinal) { - setMrelAverageForcedWarp (m, Nc, pass, results, average, secfilt, measure, &catalog[Nc].foundWarp_t[Nsecfilt*j], catalog[Nc].Nmeasure); - setGlobalObjStats (average, measure); - } - } - // if (primaryCell) free (primaryCell); + setMrelAverageForcedWarp (&catalog[Nc], Nc, j, Nsecfilt, flatcorr, results); + } + } + return (TRUE); } +// NOTE: +// Msys is measure[i].M + zp corrections +// Mcal is image[j].Mcal +// Mmos and Mgrid are offsets for mosaic and grid + +// tie down reference photometry if the -refcode (code) option is selected +// eg, -refcode g_SDSS +// this probably makes no sense in the context of multifilter analysis +// XXX probably need to use the photocde table to assign reference mag weights. + +// dlist gives the error per measurement, wlist gives the weight +// we can modify the error and weight in a few ways: +// 1) MIN_ERROR guarantees a floor +// 2) photomErrSys is added in quadrature as a sytematic error, set per photcode +// 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages +// 4) some reference photcode of some kind can be specified as fixed and have a high weight + // set mean of chip measurements (selected by photcode range for now): -int setMrelAverageExposure (off_t measureOffset, int cat, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT, off_t *found) { - +int setMrelAverageExposure (Catalog *catalog, int cat, off_t ave, int Nsecfilt, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results) { + + int Nsec; off_t k; - float Msys = 0, Mcal= 0, Mmos = 0, Mgrid = 0; - - // we are measuring means for 3 types of mags: psf, ap, kron. I am using the psf mag - // error for the ap mags, but krons have their own errors. it is an open question if I - // should be doing weighted or unweighted fits (this is a user option) - double *Mpsflist = results->Mpsflist; - double *dpsflist = results->dpsflist; - double *wpsflist = results->wpsflist; - - double *Maplist = results->Maplist; - double *daplist = results->daplist; - double *waplist = results->waplist; - - double *Mkronlist = results->Mkronlist; - double *dkronlist = results->dkronlist; - double *wkronlist = results->wkronlist; - - StatType *psfstats = &results->psfstats; - StatType *apstats = &results->apstats; - StatType *kronstats = &results->kronstats; - - int isSetMrelFinal = (pass >= 0); + + // we are guaranteed to have averageT, but not average + AverageTiny *averageT = &catalog[0].averageT[ave]; + + off_t measureOffset = averageT->measureOffset; + MeasureTiny *measureT = &catalog[0].measureT[measureOffset]; + + // we are NOT guaranteed to have average, measure, secfilt + Average *average = catalog[0].average ? &catalog[0].average[ave] : NULL; + Measure *measure = catalog[0].measure ? &catalog[0].measure[measureOffset] : NULL; + SecFilt *secfilt = catalog[0].secfilt ? &catalog[0].secfilt[ave*Nsecfilt] : NULL; + char *measureRank = catalog[0].measureRank ? &catalog[0].measureRank[measureOffset] : NULL; + + // we are measuring means for 3 types of mags: psf, ap, kron. + + // in the final assignment, set the mean mags even if only 1 measurements exists + int Nminmeas = isSetMrelFinal ? 1 : STAR_TOOFEW + 1; + + // isSetMrelFinal : in the final pass, we set psf, kron, ap mag values and extra stats; other passes only do psf mags // option for a test print - if (FALSE && (average[0].objID == 0x7146) && (average[0].catID == 0x49d8)) { + if (FALSE && (averageT[0].objID == 0x7146) && (averageT[0].catID == 0x49d8)) { fprintf (stderr, "test obj\n"); print_measure_set_alt (average, secfilt, measure); } + SetMrelInfoResetObject (results); // reset the per-object arrays + + int NextPS1 = 0; + int NpsfPS1 = 0; int GoodPS1 = FALSE; int Good2MASS = FALSE; int Galaxy2MASS = FALSE; - - int NextPS1 = 0; - int NpsfPS1 = 0; - - int Ns; - for (Ns = 0; Ns < Nphotcodes; Ns++) { - - int thisCode = photcodes[Ns][0].code; - int Nsec = GetPhotcodeNsec(thisCode); - - /* calculate the average mag in this SEC photcode for a single star */ - - /* star/photcodes already calibrated */ - if ( isSetMrelFinal && found[Nsec]) continue; - - // skip bad stars - if (!isSetMrelFinal && (secfilt[Nsec].flags & STAR_BAD)) continue; - - int NexpPS1 = 0; - int Ncode = 0; - int Next = 0; - int haveSynth = FALSE; - int havePS1 = FALSE; - - int forceSynth = FALSE; - int forceSynthEntry = -1; - - int haveUbercal = FALSE; - - int minUbercalDist = 1000; + int haveTYCHO = FALSE; + + float stargalmax = 0.0; + + // assign measurements to the photcode lists + for (k = 0; k < averageT[0].Nmeasure; k++) { + + // these bits should not be set unless we use them in this pass + measure[k].dbFlags &= ~(ID_MEAS_PHOTOM_PSF | ID_MEAS_PHOTOM_APER | ID_MEAS_PHOTOM_KRON); + + // skip measurements that do not have a valid photcode (raise exception?) + PhotCode *code = GetPhotcodebyCode (measureT[k].photcode); + myAssert (code, "invalid photcode??"); + + // SKIP gpc1 stack data (hard-wired photcodes) + if (isGPC1stack(measureT[k].photcode)) continue; - off_t meas = measureOffset; - - float psfQfMax = 0.0; - float psfQfPerfMax = 0.0; - - int Nap = 0; - int Npsf = 0; - int Nkron = 0; - for (k = 0; k < averageT[0].Nmeasure; k++, meas++) { - - // skip measurements that do not match the current photcode - PhotCode *code = GetPhotcodebyCode (measureT[k].photcode); - if (!code) continue; - if (code->equiv != thisCode) { continue; } - Ncode ++; - - // are we a PS1 exposure photcode? - if (isGPC1chip(measureT[k].photcode)) NexpPS1 ++; - - // SKIP gpc1 stack data - if (isGPC1stack(measureT[k].photcode)) continue; - - // SKIP gpc1 forced-warp data - if (isGPC1warp(measureT[k].photcode)) continue; - - if (isSetMrelFinal) { - if (measure[k].psfQF > psfQfMax) psfQfMax = measure[k].psfQF; - if (measure[k].psfQFperf > psfQfPerfMax) psfQfPerfMax = measure[k].psfQFperf; - } - - if (measureT[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad); - - if (getImageEntry (meas, cat) < 0) { - // measurements without an image are either external reference photometry or - // 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].Mcal; // check that this is zero for loaded REF value + // SKIP gpc1 forced-warp data (hard-wired photcodes) + if (isGPC1warp(measureT[k].photcode)) continue; + + if (isTYCHO(measureT[k].photcode)) { haveTYCHO = TRUE; } + + if (is2MASS(measureT[k].photcode)) { + if (measureT[k].photFlags & 0x00c00000) { + Galaxy2MASS = TRUE; // per object value + } + if (measureT[k].photFlags & 0x00000007) { + Good2MASS = TRUE; // per object value + } + } + + // assign the Nsec value so we can assign to the right lists + int Nsec = GetPhotcodeNsec (code->equiv); + if (Nsec < 0) continue; // skip measurements which do not have an equiv average photcode + + results->Nmeas[Nsec] ++; + + // various steps only relevant to the final pass + if (isSetMrelFinal) { + if (isfinite(measure[k].psfQF) && (measure[k].psfQF > results->psfQfMax[Nsec])) results->psfQfMax[Nsec] = measure[k].psfQF; + if (isfinite(measure[k].psfQFperf) && (measure[k].psfQFperf > results->psfQfPerfMax[Nsec])) results->psfQfPerfMax[Nsec] = measure[k].psfQFperf; + + // are we a PS1 exposure photcode? (hard-wired photcodes) + if (isGPC1chip(measure[k].photcode)) { + results->NexpPS1[Nsec] ++; + results->havePS1[Nsec] = TRUE; + + // only count psfQF > 0.85 here + if (isfinite(measure[k].extNsigma) && isfinite(measure[k].psfQF) && (measure[k].psfQF > 0.85)) { + stargalmax = MAX (stargalmax, measure[k].extNsigma); + } + } + + // force the use of SYN even if we have PS1 mags? + if (isGPC1synth(measure[k].photcode)) { + results->haveSYN[Nsec] = TRUE; + results->measSYN[Nsec] = k; + measureT[k].dbFlags |= ID_MEAS_SYNTH_MAG; // redundant with photcode + // note that synthetic mags have the real mags in measure.M (no zero point offset) + if (measure[k].M < results->minSYN[Nsec]) results->needSYN[Nsec] = TRUE; + continue; + } + } + + // ** Choose the calibration (depends on the mode : do I have an image reference or not?) + float Mcal = 0, Mmos = 0, Mgrid = 0; + off_t meas = measureOffset + k; + if (getImageEntry (meas, cat) < 0) { + // measurements without an image are either external reference photometry or + // 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].Mcal; // check that this is zero for loaded REF value + } else { + // getMcal_alt returns image[].Mcal modified by flatcorr(image.photom_map_id,x,y) + Mcal = getMcal_alt (meas, cat, flatcorr, measureT[k].Xccd, measureT[k].Yccd); + if (isnan(Mcal)) SKIP_THIS_MEAS(Ncal); + Mmos = getMmos (meas, cat); + if (isnan(Mmos)) SKIP_THIS_MEAS(Nmos); + Mgrid = getMgrid (meas, cat); + if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid); + } + + int myUbercalDist = getUbercalDist(meas, cat); + results->minUbercalDist[Nsec] = MIN(results->minUbercalDist[Nsec], myUbercalDist); + + int isUbercal = (measureT[k].dbFlags & ID_MEAS_PHOTOM_UBERCAL); + + if (isUbercal) results->haveUbercal[Nsec] = TRUE; // haveUbercal is set per secfilt, isUbercal is per measure XXX define this array + + int useUbercalWeight = isUbercal || (refPhotcode && (code->code == refPhotcode->code)); + + float Map = NAN; + float dMap = NAN; + if (isSetMrelFinal) { + Map = PhotCat (&measure[k], MAG_CLASS_APER); + dMap = MAX (hypot(measure[k].dMap, code->photomErrSys), MIN_ERROR); + if (CHECK_VALID_MAG(Map, dMap)) { + int Nap = results->aperData[Nsec].Nlist; + results->aperData[Nsec].flxlist[Nap] = Map - Mcal - Mmos - Mgrid; + results->aperData[Nsec].errlist[Nap] = dMap; + results->aperData[Nsec].wgtlist[Nap] = useUbercalWeight ? UBERCAL_WEIGHT : 1.0; + results->aperData[Nsec].ranking[Nap] = measureRank[k]; + results->aperData[Nsec].measSeq[Nap] = k; + results->aperData[Nsec].Nlist ++; + } + } + + float Mkron, dMkron; + if (isSetMrelFinal) { + Mkron = PhotCat (&measure[k], MAG_CLASS_KRON); + dMkron = MAX (hypot(measure[k].dMkron, code->photomErrSys), MIN_ERROR); + if (CHECK_VALID_MAG(Mkron, dMkron)) { + int Nkron = results->kronData[Nsec].Nlist; + results->kronData[Nsec].flxlist[Nkron] = Mkron - Mcal - Mmos - Mgrid; + results->kronData[Nsec].errlist[Nkron] = dMkron; + results->kronData[Nsec].wgtlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0; + results->kronData[Nsec].ranking[Nkron] = measureRank[k]; + results->kronData[Nsec].measSeq[Nkron] = k; + results->kronData[Nsec].Nlist ++; + } + } + + float Mpsf = PhotSysTiny (&measureT[k], &averageT[0], &secfilt[0], MAG_CLASS_PSF); + float dMpsf = MAX (hypot(measureT[k].dM, code->photomErrSys), MIN_ERROR); + if (CHECK_VALID_MAG(Mpsf, dMpsf)) { + int Npsf = results->psfData[Nsec].Nlist; + results->psfData[Nsec].flxlist[Npsf] = Mpsf - Mcal - Mmos - Mgrid; + results->psfData[Nsec].errlist[Npsf] = dMpsf; + results->psfData[Nsec].wgtlist[Npsf] = isUbercal ? UBERCAL_WEIGHT : 1.0; + results->psfData[Nsec].ranking[Npsf] = measureRank ? measureRank[k] : 0; + results->psfData[Nsec].measSeq[Npsf] = k; + results->psfData[Nsec].Nlist ++; + } + + // count the extended detections (all PS1 bands) + if (isSetMrelFinal && isGPC1chip(measureT[k].photcode) && !isnan(Map) && !isnan(Mpsf)) { + float dMagAp = Mpsf - Map; + float SigmaAp = hypot(0.1, 2.5*dMpsf); + // XXX this is still quite ad hoc, but at least it: + // (a) converges to 0.1 mag offset at the bright end + // (b) converges to 0.5 mag offset at the faint end (dM = 0.2) + if (dMagAp > SigmaAp) { + results->Next[Nsec] ++; + NextPS1 ++; } else { - Mcal = getMcal_alt (meas, cat, flatcorr, measureT[k].Xccd, measureT[k].Yccd); - if (isnan(Mcal)) SKIP_THIS_MEAS(Ncal); - Mmos = getMmos (meas, cat); - if (isnan(Mmos)) SKIP_THIS_MEAS(Nmos); - Mgrid = getMgrid (meas, cat); - if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid); - } - - // skip some absurd values NAN, < 0.0, > 30.0 - Msys = PhotSysTiny (&measureT[k], &averageT[0], &secfilt[0], MAG_CLASS_PSF); - if (isnan(Msys)) SKIP_THIS_MEAS(Nsys); - if (Msys < 0.0) SKIP_THIS_MEAS(Nsys); - if (Msys > 30.0) SKIP_THIS_MEAS(Nsys); - - int myUbercalDist = getUbercalDist(meas, cat); - minUbercalDist = MIN(minUbercalDist, myUbercalDist); - - int isUbercal = (measureT[k].dbFlags & ID_MEAS_PHOTOM_UBERCAL); - float dMpsf = MAX (hypot(measureT[k].dM, code->photomErrSys), MIN_ERROR); - - if (isUbercal) haveUbercal = TRUE; // haveUbercal is set per secfilt, isUbercal is per measure - - if (isSetMrelFinal) { - // special options for PS1 data - if ((measure[k].photcode >= 10000) && (measure[k].photcode <= 10500)) { - // count the extended detections - if (!isnan(measure[k].Map)) { - float dMagAp = measure[k].M - measure[k].Map; - float SigmaAp = hypot(0.1, 2.5*measure[k].dM); - // XXX this is still quite ad hoc, but at least it: - // (a) converges to 0.1 mag offset at the bright end - // (b) converges to 0.5 mag offset at the faint end (dM = 0.2) - if (dMagAp > SigmaAp) { - Next ++; - NextPS1 ++; - } else { - NpsfPS1 ++; - } - } - havePS1 = TRUE; + NpsfPS1 ++; + } + } + } + + SynthZeroPoints *synthzpts = SynthZeroPointsGet(); + + // find max values across filters + float psfQfMax = 0.0; + float psfQfPerfMax = 0.0; + + // now calculate the mean stats for the Nsec bands. + for (Nsec = 0; Nsec < Nsecfilt; Nsec++) { + dvo_secfilt_init (&secfilt[Nsec], SECFILT_RESET_CHIP); // this does not reset astrometry or STACK bits + + // XXX hardwired + if ((Nsec < 4) || (Nsec == 8)) { + secfilt[Nsec].Ncode = results->NexpPS1[Nsec]; + } else { + secfilt[Nsec].Ncode = results->Nmeas[Nsec]; // 2MASS data if it exists + } + + if (haveTYCHO) { + secfilt[Nsec].flags |= ID_SECF_HAS_TYCHO; + } + + // force the synthetic magnitude (apply zpt calibration if no tycho photometry exists) + if (results->needSYN[Nsec]) { + if (!isSetMrelFinal) continue; + + // use the single SYNTH value instead of the other mags here + float Mpsf = measure[results->measSYN[Nsec]].M; + float ZP = 0.0; + + // if we have loaded the synth zpt correction, and do NOT have a tycho measurement, apply the correction + if (!haveTYCHO && synthzpts) { + // need to look up the (X,Y) coords from this (R,D) location + double X, Y; + double R = ohana_normalize_angle_to_midpoint (measure->R, 0.0); // XXX 0.0 or 180.0 for center? + RD_to_XY (&X, &Y, R, measure->D, &synthzpts->coords); + if (X < 0) return FALSE; + if (Y < 0) return FALSE; + if (X >= synthzpts->Nx) return FALSE; + if (Y >= synthzpts->Ny) return FALSE; + + int Xpix = X; + int Ypix = Y; + int Npix = Xpix + Ypix*synthzpts->Nx; + + float *value = (float *) synthzpts->matrix[Nsec].buffer; + ZP = !isnan(value[Npix]) ? value[Npix] : 0.0; + if (isfinite(value[Npix])) { + secfilt[Nsec].flags |= ID_SECF_FIX_SYNTH; } - - // count extended detections for 2MASS (XXX NOTE hardwired photcodes 2011, 2012, 2013) - if ((measure[k].photcode >= 2011) && (measure[k].photcode <= 2013)) { - if (measure[k].photFlags & 0x00c00000) { - Next ++; - Galaxy2MASS = TRUE; - } - if (pass == 0) { - if (measure[k].photFlags & 0x00000007) { - Good2MASS = TRUE; - } else { - // detections without one of these bits should only be used in PASS_1 - SKIP_THIS_MEAS(Nbad); - } - } - } - - // Blindly accepth the SYNTH mags if we are above saturation, otherwise, - // ignore SYNTH photcodes until PASS == 4 (where we also accept saturated stars) - if ((measure[k].photcode >= 3001) && (measure[k].photcode <= 3005)) { - // something of a hack: force object to use synth values if synth mags >> - // saturation (3pi instrumental mags < -15) - float MaxMagForceSynth = NAN; - switch (measure[k].photcode) { - case 3001: - MaxMagForceSynth = 13.64; - break; - case 3002: - MaxMagForceSynth = 13.76; - break; - case 3003: - MaxMagForceSynth = 13.74; - break; - case 3004: - MaxMagForceSynth = 12.94; - break; - case 3005: - MaxMagForceSynth = 12.01; - break; - } - if (measureT[k].M < MaxMagForceSynth) { - forceSynth = TRUE; - forceSynthEntry = Npsf; - } else { - if (pass < 4) { - SKIP_THIS_MEAS(Nbad); - } - haveSynth = TRUE; - } - } - - // Map (Maplist) and Mkron (Mkronlist,dkronlist) are used to calculate mean mags per filter - float Map = PhotCat (&measure[k], MAG_CLASS_APER); - if (!isnan(Map)) { - Maplist[Nap] = Map - Mcal - Mmos - Mgrid; - daplist[Nap] = dMpsf;// XXX check on this... - waplist[Nap] = isUbercal ? UBERCAL_WEIGHT : 1.0; - Nap ++; - } - - float Mkron = PhotCat (&measure[k], MAG_CLASS_KRON); - if (!isnan(Mkron)) { - Mkronlist[Nkron] = Mkron - Mcal - Mmos - Mgrid; - dkronlist[Nkron] = measure[k].dMkron; - wkronlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0; - Nkron ++; - } - } // if (isSetMrelFinal) - - // dlist gives the error per measurement, wlist gives the weight - // we can modify the error and weight in a few ways: - // 1) MIN_ERROR guarantees a floor - // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode - // 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages - // 4) some reference photcode of some kind can be specified as fixed and have a high weight - Mpsflist[Npsf] = Msys - Mcal - Mmos - Mgrid; - dpsflist[Npsf] = dMpsf; - wpsflist[Npsf] = isUbercal ? UBERCAL_WEIGHT : 1.0; - - // NOTE: - // Msys is measure[i].M + zp corrections - // Mcal is image[j].Mcal - // Mmos and Mgrid are offsets for mosaic and grid - - // tie down reference photometry if the -refcode (code) option is selected - // eg, -refcode g_SDSS - // this probably makes no sense in the context of multifilter analysis - // XXX probably need to use the photocde table to assign reference mag weights. - if (refPhotcode) { - if (code->code == refPhotcode->code) { - wpsflist[Npsf] = UBERCAL_WEIGHT; - } - } - Npsf ++; - } - - int Nminmeas = isSetMrelFinal ? 1 : STAR_TOOFEW + 1; - - // XXX : ugh : another hard-wired photcode entry... - if (isSetMrelFinal && (pass == 0)) { - if ((thisCode < 6) || (thisCode == 9)) { - secfilt[Nsec].Ncode = NexpPS1; - } else { - secfilt[Nsec].Ncode = Ncode; // 2MASS data if it exists - } - } - - // when performing the grid analysis, STAR_TOOFEW should be set to 1; - if (Npsf < Nminmeas) { /* too few measurements */ - // fprintf (f, "%10.6f %10.6f %d %d %d\n", averageT[0].R, averageT[0].D, measureT[0].imageID, Npsf, STAR_TOOFEW); + } + secfilt[Nsec].M = Mpsf + ZP; + secfilt[Nsec].dM = 0.6; + secfilt[Nsec].Mchisq = 0.0; + secfilt[Nsec].flags |= ID_SECF_USE_SYNTH; + continue; + } + + // if too few valid measurements meet the minimum criteria, go to the next entry + StatType *psfstats = &results->psfstats; + int Nranking = magStatsByRanking (&results->psfData[Nsec], psfstats); + if (Nranking < Nminmeas) { secfilt[Nsec].flags |= ID_STAR_FEW; - if (Ncode == 0) { - results->Ncode ++; - } else { - results->Nfew ++; - } - continue; } else { - secfilt[Nsec].flags &= ~ID_STAR_FEW; - } - - if (forceSynth) { - // use the single SYNTH value instead of the other mags here - myAssert ((forceSynthEntry < Npsf) && (forceSynthEntry >= 0), "programming error"); - Mpsflist[0] = Mpsflist[forceSynthEntry]; - dpsflist[0] = dpsflist[forceSynthEntry]; - wpsflist[0] = wpsflist[forceSynthEntry]; - Npsf = 1; - } - liststats (Mpsflist, dpsflist, wpsflist, Npsf, psfstats); - - secfilt[Nsec].M = psfstats->mean; - secfilt[Nsec].dM = psfstats->error; - secfilt[Nsec].Mchisq = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN; - - // when running -averages, we have no information about the images, so we cannot set this - if (minUbercalDist > -1) { - secfilt[Nsec].ubercalDist = minUbercalDist; - } + secfilt[Nsec].M = psfstats->mean; + secfilt[Nsec].dM = psfstats->error; + secfilt[Nsec].Mchisq = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN; + } + int minRank = (Nranking > 0) ? results->psfData[Nsec].ranking[0] : 10; if (isSetMrelFinal) { - found[Nsec] = TRUE; - - secfilt[Nsec].Mstdev = psfstats->sigma; // Mstdev is in millimags (not enough space for more precision) - // secfilt[Nsec].Ncode = Ncode; - secfilt[Nsec].Nused = psfstats->Nmeas; - - secfilt[Nsec].Mmax = psfstats->max; - secfilt[Nsec].Mmin = psfstats->min; - - secfilt[Nsec].psfQfMax = psfQfMax; - secfilt[Nsec].psfQfPerfMax = psfQfPerfMax; - - // NOTE : use the modified weight for apmags as well as psf mags - liststats (Maplist, daplist, waplist, Nap, apstats); - secfilt[Nsec].Map = Nap > 0 ? apstats->mean : NAN; - secfilt[Nsec].dMap = Nap > 0 ? apstats->error : NAN; - secfilt[Nsec].sMap = Nap > 0 ? apstats->sigma : NAN; - secfilt[Nsec].NusedAp = Nap; - - liststats (Mkronlist, dkronlist, wkronlist, Nkron, kronstats); - secfilt[Nsec].Mkron = Nkron > 0 ? kronstats->mean : NAN; - secfilt[Nsec].dMkron = Nkron > 0 ? kronstats->error : NAN; - secfilt[Nsec].sMkron = Nkron > 0 ? kronstats->sigma : NAN; - secfilt[Nsec].NusedKron = Nkron; - - // NOTE: for 2MASS measurements, Next should be 1, as should N - if ((Next > 0) && (Next > 0.5*Npsf)) { + if (FALSE && (average->objID == 0x1) && (average->catID)) { + fprintf (stderr, "test obj\n"); + } + + // mark the measurements matching this ranking + markMeasureByRanking (&results->psfData[Nsec], measure, minRank, ID_MEAS_PHOTOM_PSF); + + // when running -averages, we have no information about the images, so we cannot set this + if (results->minUbercalDist[Nsec] > -1) { + secfilt[Nsec].ubercalDist = results->minUbercalDist[Nsec]; + } + + if (Nranking) { + secfilt[Nsec].Mstdev = psfstats->sigma; + secfilt[Nsec].Nused = psfstats->Nmeas; + secfilt[Nsec].Mmax = psfstats->max; + secfilt[Nsec].Mmin = psfstats->min; + } + + secfilt[Nsec].psfQfMax = results->psfQfMax[Nsec]; + secfilt[Nsec].psfQfPerfMax = results->psfQfPerfMax[Nsec]; + + if (isfinite(secfilt[Nsec].psfQfMax) && (secfilt[Nsec].psfQfMax > psfQfMax)) psfQfMax = secfilt[Nsec].psfQfMax; + if (isfinite(secfilt[Nsec].psfQfPerfMax) && (secfilt[Nsec].psfQfPerfMax > psfQfPerfMax)) psfQfPerfMax = secfilt[Nsec].psfQfPerfMax; + + StatType *apstats = &results->apstats; + Nranking = magStatsByRanking (&results->aperData[Nsec], apstats); + if (Nranking) { + secfilt[Nsec].Map = apstats->mean; + secfilt[Nsec].dMap = apstats->error; + secfilt[Nsec].sMap = apstats->sigma; + secfilt[Nsec].NusedAp = Nranking; + } + markMeasureByRanking (&results->aperData[Nsec], measure, minRank, ID_MEAS_PHOTOM_APER); + + StatType *kronstats = &results->kronstats; + Nranking = magStatsByRanking (&results->kronData[Nsec], kronstats); + if (Nranking) { + secfilt[Nsec].Mkron = kronstats->mean; + secfilt[Nsec].dMkron = kronstats->error; + secfilt[Nsec].sMkron = kronstats->sigma; + secfilt[Nsec].NusedKron = Nranking; + } + markMeasureByRanking (&results->kronData[Nsec], measure, minRank, ID_MEAS_PHOTOM_KRON); + + // does this object appear extended in > 50% of measurements? + if ((results->Next[Nsec] > 0) && (results->Next[Nsec] > 0.5*psfstats->Nmeas)) { secfilt[Nsec].flags |= ID_SECF_OBJ_EXT; } - switch (pass) { + switch (minRank) { case 0: - secfilt[Nsec].flags |= ID_PHOTOM_PASS_0; - if (havePS1) GoodPS1 = TRUE; + secfilt[Nsec].flags |= ID_SECF_RANK_0; + if (results->havePS1[Nsec]) GoodPS1 = TRUE; break; case 1: - secfilt[Nsec].flags |= ID_PHOTOM_PASS_1; - if (havePS1) GoodPS1 = TRUE; + secfilt[Nsec].flags |= ID_SECF_RANK_1; + if (results->havePS1[Nsec]) GoodPS1 = TRUE; break; case 2: - secfilt[Nsec].flags |= ID_PHOTOM_PASS_2; - if (havePS1) GoodPS1 = TRUE; + secfilt[Nsec].flags |= ID_SECF_RANK_2; + if (results->havePS1[Nsec]) GoodPS1 = TRUE; break; case 3: - secfilt[Nsec].flags |= ID_PHOTOM_PASS_3; + secfilt[Nsec].flags |= ID_SECF_RANK_3; break; case 4: - secfilt[Nsec].flags |= ID_PHOTOM_PASS_4; + secfilt[Nsec].flags |= ID_SECF_RANK_4; break; - } - if (haveSynth) { - secfilt[Nsec].flags |= ID_SECF_USE_SYNTH; - } - if (havePS1) { + default: + break; + } + if (results->havePS1[Nsec]) { secfilt[Nsec].flags |= ID_SECF_HAS_PS1; } - if (haveUbercal) { + if (results->haveUbercal[Nsec]) { secfilt[Nsec].flags |= ID_SECF_USE_UBERCAL; - } + } } } if (isSetMrelFinal) { - if (pass == 0) { - DVOAverageFlags flagBits = ID_OBJ_EXT | ID_OBJ_EXT_ALT | ID_OBJ_GOOD | ID_OBJ_GOOD_ALT; - // we attempt to set a few flags here; reset those bits before trying: - average[0].flags &= ~flagBits; - } - if (NextPS1 && (NextPS1 > NpsfPS1)) { average[0].flags |= ID_OBJ_EXT; @@ -446,4 +435,6 @@ average[0].flags |= ID_OBJ_GOOD_ALT; } + average[0].psfQF = psfQfMax; + average[0].psfQFperf = psfQfPerfMax; } return (TRUE); @@ -455,9 +446,17 @@ // 2) select the BEST detections per filter (regardless of PRIMARY) // 3) apply the zero point and AB->Jy transformations -int setMrelAverageStack (off_t measureOffset, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t Nmeasure) { +int setMrelAverageStack (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr) { + + // we are guaranteed to have average, measure, secfilt + Average *average = &catalog[0].average[ave]; + off_t measureOffset = average->measureOffset; + int Nmeasure = average->Nmeasure; + + Measure *measure = &catalog[0].measure[measureOffset]; + SecFilt *secfilt = &catalog[0].secfilt[ave*Nsecfilt]; off_t k; - float Mcal= 0, Mmos = 0, Mgrid = 0, Finst = 0; + float Mcal = 0, Mmos = 0, Mgrid = 0, Finst = 0; // set the primary projection cell and skycell for this coordinate @@ -500,5 +499,5 @@ off_t meas = measureOffset; - for (k = 0; k < average[0].Nmeasure; k++, meas++) { + for (k = 0; k < Nmeasure; k++, meas++) { // only examine gpc1 stack data @@ -528,5 +527,5 @@ measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY; secfilt[Nsec].stackPrmryOff = meas; - myAssert (secfilt[Nsec].stackPrmryOff <= Nmeasure, "stackPrmryOff out of range"); + myAssert (secfilt[Nsec].stackPrmryOff <= catalog[0].Nmeasure, "stackPrmryOff out of range"); } @@ -555,10 +554,15 @@ } - // NOTE: negative and insignificant vlues values are allowed, but not NAN flux values + // NOTE: negative and insignificant flux values are allowed, but not NAN flux values Finst = PhotFluxInst (&measure[k], MAG_CLASS_PSF); if (isnan(Finst)) SKIP_THIS_MEAS_STACK(Ninst); + // choose the best psfQFperf value for the BEST measurement if (measure[k].psfQFperf > psfQFbest) { psfQFbest = measure[k].psfQFperf; + stackBestMeasure = k; + } + // ... UNLESS psfQFperf > 0.98 for the primary, in which case just use the primary. + if ((measure[k].dbFlags & ID_MEAS_STACK_PRIMARY) && (measure[k].psfQFperf > 0.98)) { stackBestMeasure = k; } @@ -640,5 +644,5 @@ secfilt[Nsec].stackBestOff = k + measureOffset; - myAssert (secfilt[Nsec].stackBestOff <= Nmeasure, "stackBestOff out of range"); + myAssert (secfilt[Nsec].stackBestOff <= catalog[0].Nmeasure, "stackBestOff out of range"); secfilt[Nsec].Nstack = Nstack; @@ -661,251 +665,246 @@ average[0].flags |= ID_OBJ_BAD_STACK; } + return (TRUE); } - -# undef SKIP_THIS_MEAS -# define SKIP_THIS_MEAS(REASON) { \ - measure[k].dbFlags |= ID_MEAS_SKIP_PHOTOM; \ - results->REASON ++; \ - continue; } // set mean of forced-warp measurements (selected by photcode range for now): // somewhat simplified relative to chip-photometry: // * no grid, no mosaic, no 2MASS, no SYNTH, no Ubercal, no flatcorr -int setMrelAverageForcedWarp (off_t measureOffset, int cat, int pass, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t *found, off_t Nmeasure) { - - off_t k; - float Mcal= 0; - - // we are measuring means for 3 types of FLUXes: psf, ap, kron. I am using the psf mag - // error for the ap mags, but krons have their own errors. it is an open question if I - // should be doing weighted or unweighted fits (this is a user option) - double *Fpsflist = results->Mpsflist; - double *dpsflist = results->dpsflist; - double *wpsflist = results->wpsflist; - - double *Faplist = results->Maplist; - double *daplist = results->daplist; - double *waplist = results->waplist; - - double *Fkronlist = results->Mkronlist; - double *dkronlist = results->dkronlist; - double *wkronlist = results->wkronlist; - - StatType *psfstats = &results->psfstats; - StatType *apstats = &results->apstats; - StatType *kronstats = &results->kronstats; +// analysis is done on flux, not mags (as the faintest objects will be nearly insignificant) +int setMrelAverageForcedWarp (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr, SetMrelInfo *results) { + + // we are guaranteed to have average, measure, secfilt + Average *average = &catalog[0].average[ave]; + off_t measureOffset = average->measureOffset; + + Measure *measure = &catalog[0].measure[measureOffset]; + SecFilt *secfilt = &catalog[0].secfilt[ave*Nsecfilt]; + char *measureRank = &catalog[0].measureRank[measureOffset]; + + off_t k, Nsec; + float Mcal = 0; + + SetMrelInfoResetObject (results); // reset the per-object arrays // option for a test print - if ((average[0].objID == 0xd6e) && (average[0].catID == 0x4984)) { + if (FALSE && (average[0].objID == 0x0000020d) && (average[0].catID == 0x000049e1)) { fprintf (stderr, "test obj\n"); print_measure_set_alt (average, secfilt, measure); } - // for a diff db, we just want primary, if available or not - int priTessID, priProjID, priSkycellID; + // for a diff db, we just want primary, if available or not. for warp we want to use the + // skycell which provided the stack best entry + int tessID, projID, skycellID; if (IS_DIFF_DB) { - get_tess_ids(&priTessID, &priProjID, &priSkycellID, average[0].R, average[0].D); - } - - int Ns; - for (Ns = 0; Ns < Nphotcodes; Ns++) { - - int thisCode = photcodes[Ns][0].code; - int Nsec = GetPhotcodeNsec(thisCode); - - /* calculate the average mag in this SEC photcode for a single star */ - - /* star/photcodes already calibrated */ - if (found[Nsec]) continue; - - // I need to restrict the measurements coming only from the same skycell as the secfilt values - // if we cannot make that association, skip this set of warps - int stkTessID, stkProjID, stkSkycellID; - + get_tess_ids(&tessID, &projID, &skycellID, average[0].R, average[0].D); + } + for (Nsec = 0; Nsec < Nsecfilt; Nsec++) { if (!IS_DIFF_DB) { - myAssert (secfilt[Nsec].stackBestOff <= Nmeasure, "stackBestOff out of range"); - if (!FindImageSkycellID (secfilt[Nsec].stackBestOff, cat, &stkTessID, &stkProjID, &stkSkycellID)) continue; - } - - float psfQfMax = 0.0; - float psfQfPerfMax = 0.0; - - off_t meas = measureOffset; - - int Nap = 0; - int Npsf = 0; - int Nkron = 0; - - int Nwarp = 0; - int NwarpGood = 0; - for (k = 0; k < average[0].Nmeasure; k++, meas++) { - - // only examine gpc1 forced-warp data - if (!isGPC1warp(measure[k].photcode)) continue; - - // skip measurements that do not match the current photcode - PhotCode *code = GetPhotcodebyCode (measure[k].photcode); - if (!code) continue; - if (code->equiv != thisCode) { continue; } - - Nwarp ++; - if (measure[k].psfQFperf > 0.85) NwarpGood ++; - - measure[k].dbFlags &= ~ID_MEAS_WARP_USED; - - if (IS_DIFF_DB) { - if (MatchImageSkycellID (meas, cat, priTessID, priProjID, priSkycellID)) { - measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY; - } else { - continue; - } - } else { - if (!MatchImageSkycellID (meas, cat, stkTessID, stkProjID, stkSkycellID)) continue; - } - - if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad); - - if (getImageEntry (meas, cat) < 0) { - // measurements without an image are either external reference photometry or - // 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 - Mcal = measure[k].Mcal; // check that this is zero for loaded REF value - } else { - // use getMcal not getMcal_alt? - Mcal = getMcal_alt (meas, cat, NULL, measure[k].Xccd, measure[k].Yccd); - // Mcal = getMcal (meas, cat); - if (isnan(Mcal)) SKIP_THIS_MEAS(Ncal); - } - float Fcal = MagToFlux(-Mcal); - - // these are only used below if IS_DIFF_DB - if (IS_DIFF_DB) { - if (measure[k].psfQF > psfQfMax) psfQfMax = measure[k].psfQF; - if (measure[k].psfQFperf > psfQfPerfMax) psfQfPerfMax = measure[k].psfQFperf; - } - - // in the calculations below, - // ...list gives the error per measurement, wlist gives the weight - // we can modify the error and weight in a few ways: - // 1) MIN_ERROR guarantees a floor - // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode - - // skip some absurd values NAN - // NOTE : I am using PhotCat not PhotSys for now since GPC1 chip-to-chip color terms - // are small (and not measured) - float Fpsf = PhotFluxCat (&measure[k], MAG_CLASS_PSF); - if (isnan(Fpsf)) SKIP_THIS_MEAS(Nsys); - // if (Msys < 0.0) SKIP_THIS_MEAS(Nsys); - // if (Msys > 30.0) SKIP_THIS_MEAS(Nsys); - - float dFpsf = PhotFluxCatErr (&measure[k], MAG_CLASS_PSF); - if (isnan(dFpsf)) SKIP_THIS_MEAS (Nsys); - - int isBad = (measure[k].psfQF < 0.85); - isBad |= isnan(measure[k].psfQF); - if (isBad) SKIP_THIS_MEAS (Nsys); - - int isPoor = (measure[k].psfQFperf < 0.85); - if ((pass == 0) && isPoor) SKIP_THIS_MEAS (Nsys); - - dFpsf = MAX (dFpsf, MIN_ERROR*Fpsf); // MIN_ERROR is a fractional error - Fpsflist[Npsf] = Fpsf * Fcal; - dpsflist[Npsf] = dFpsf * Fcal; - wpsflist[Npsf] = 1.0; - Npsf ++; - - measure[k].dbFlags |= ID_MEAS_WARP_USED; - - float Fap = PhotFluxCat (&measure[k], MAG_CLASS_APER); - float dFap = PhotFluxCatErr (&measure[k], MAG_CLASS_APER); - if (!isnan(Fap)) { - Faplist[Nap] = Fap * Fcal; - daplist[Nap] = dFap * Fcal; - waplist[Nap] = 1.0; // drop weight lists? - Nap ++; - } - - float Fkron = PhotFluxCat (&measure[k], MAG_CLASS_KRON); - float dFkron = PhotFluxCatErr (&measure[k], MAG_CLASS_KRON); - if (!isnan(Fkron)) { - Fkronlist[Nkron] = Fkron * Fcal; - dkronlist[Nkron] = dFkron * Fcal; - wkronlist[Nkron] = 1.0; - Nkron ++; - } - } - if (Npsf < 1) { - secfilt[Nsec].Nwarp = Nwarp; - secfilt[Nsec].NwarpGood = NwarpGood; - secfilt[Nsec].NusedWrp = Npsf; - secfilt[Nsec].NusedApWrp = Nap; - secfilt[Nsec].NusedKronWrp = Nkron; + myAssert (secfilt[Nsec].stackBestOff <= catalog->Nmeasure, "stackBestOff out of range"); + if (!FindImageSkycellID (secfilt[Nsec].stackBestOff, cat, &tessID, &projID, &skycellID)) continue; + } + results->tessID[Nsec] = tessID; + results->projID[Nsec] = projID; + results->skycellID[Nsec] = skycellID; + } + + // assign measurements to the photcode lists + for (k = 0; k < average->Nmeasure; k++) { + off_t meas = measureOffset + k; + + measure[k].dbFlags &= ~ID_MEAS_WARP_USED; // this should not be set unless we use it in this pass + + // skip measurements that do not have a valid photcode (raise exception?) + PhotCode *code = GetPhotcodebyCode (measure[k].photcode); + myAssert (code, "invalid photcode??"); + + // only examine gpc1 forced-warp data + if (!isGPC1warp(measure[k].photcode)) continue; + + // assign the Nsec value so we can assign to the right lists + int Nsec = GetPhotcodeNsec (code->equiv); + if (Nsec < 0) continue; // skip measurements which do not have an equiv average photcode + + results->Nmeas[Nsec] ++; + if (measure[k].psfQFperf > 0.85) results->NmeasGood[Nsec] ++; + + // use primary skycell for DIFF, stack skycell for WARP + if (MatchImageSkycellID (meas, cat, results->tessID[Nsec], results->projID[Nsec], results->skycellID[Nsec])) { continue; } - found[Nsec] = TRUE; - - liststats (Fpsflist, dpsflist, wpsflist, Npsf, psfstats); + // assign the max psfQF values for each Nsec value + if (IS_DIFF_DB) { + measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY; + if (measure[k].psfQF > results->psfQfMax[Nsec]) results->psfQfMax[Nsec] = measure[k].psfQF; + if (measure[k].psfQFperf > results->psfQfPerfMax[Nsec]) results->psfQfPerfMax[Nsec] = measure[k].psfQFperf; + } + + if (getImageEntry (meas, cat) < 0) { + // measurements without an image are either external reference photometry or + // 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 + Mcal = measure[k].Mcal; // check that this is zero for loaded REF value + } else { + // use getMcal not getMcal_alt? + Mcal = getMcal_alt (meas, cat, NULL, measure[k].Xccd, measure[k].Yccd); + // Mcal = getMcal (meas, cat); + if (isnan(Mcal)) SKIP_THIS_MEAS(Ncal); + } + float Fcal = MagToFlux(-Mcal); + + // NOTE : I am using PhotCat not PhotSys for now since GPC1 chip-to-chip color terms + // are small (and not measured) + float Fpsf = PhotFluxCat (&measure[k], MAG_CLASS_PSF); + float dFpsf = PhotFluxCatErr (&measure[k], MAG_CLASS_PSF); + if (CHECK_VALID_FLUX(Fpsf, dFpsf)) { + dFpsf = MAX (hypot(dFpsf, code->photomErrSys*Fpsf), MIN_ERROR*Fpsf); // bump up the error by a systematic floor + int Npsf = results->psfData[Nsec].Nlist; + results->psfData[Nsec].flxlist[Npsf] = Fpsf * Fcal; + results->psfData[Nsec].errlist[Npsf] = dFpsf * Fcal; + results->psfData[Nsec].wgtlist[Npsf] = 1.0; + results->psfData[Nsec].ranking[Npsf] = measureRank ? measureRank[k] : 0; + results->psfData[Nsec].measSeq[Npsf] = k; + results->psfData[Nsec].Nlist ++; + } + + float Fap = PhotFluxCat (&measure[k], MAG_CLASS_APER); + float dFap = PhotFluxCatErr (&measure[k], MAG_CLASS_APER); + if (CHECK_VALID_FLUX(Fap, dFap)) { + dFap = MAX (hypot(dFap, code->photomErrSys*Fap), MIN_ERROR*Fap); // bump up the error by a systematic floor + int Naper = results->aperData[Nsec].Nlist; + results->aperData[Nsec].flxlist[Naper] = Fap * Fcal; + results->aperData[Nsec].errlist[Naper] = dFap * Fcal; + results->aperData[Nsec].wgtlist[Naper] = 1.0; + results->aperData[Nsec].ranking[Naper] = measureRank ? measureRank[k] : 0; + results->aperData[Nsec].measSeq[Naper] = k; + results->aperData[Nsec].Nlist ++; + } + + float Fkron = PhotFluxCat (&measure[k], MAG_CLASS_KRON); + float dFkron = PhotFluxCatErr (&measure[k], MAG_CLASS_KRON); + if (CHECK_VALID_FLUX(Fkron, dFkron)) { + dFkron = MAX (hypot(dFkron, code->photomErrSys*Fkron), MIN_ERROR*Fkron); // bump up the error by a systematic floor + int Nkron = results->kronData[Nsec].Nlist; + results->kronData[Nsec].flxlist[Nkron] = Fkron * Fcal; + results->kronData[Nsec].errlist[Nkron] = dFkron * Fcal; + results->kronData[Nsec].wgtlist[Nkron] = 1.0; + results->kronData[Nsec].ranking[Nkron] = measureRank ? measureRank[k] : 0; + results->kronData[Nsec].measSeq[Nkron] = k; + results->kronData[Nsec].Nlist ++; + } + } + + // find max values across filters + float psfQfMax = 0.0; + float psfQfPerfMax = 0.0; + + // now calculate the mean stats for the Nsec bands. + for (Nsec = 0; Nsec < Nsecfilt; Nsec++) { + dvo_secfilt_init (&secfilt[Nsec], SECFILT_RESET_WARP); if (IS_DIFF_DB) { // for non DIFF_DB, these are set in setMrelAverageExposure - secfilt[Nsec].psfQfMax = psfQfMax; - secfilt[Nsec].psfQfPerfMax = psfQfPerfMax; - } - - secfilt[Nsec].FpsfWrp = psfstats->mean; - secfilt[Nsec].dFpsfWrp = psfstats->error; - secfilt[Nsec].sFpsfWrp = psfstats->sigma; // Mstdev is in millimags (not enough space for more precision) - secfilt[Nsec].NusedWrp = psfstats->Nmeas; - secfilt[Nsec].MpsfWrp = isnan(secfilt[Nsec].FpsfWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FpsfWrp); // 8.9 since flux is in Jy - - // NOTE : use the modified weight for apmags as well as psf mags - liststats (Faplist, daplist, waplist, Nap, apstats); - secfilt[Nsec].FapWrp = Nap > 0 ? apstats->mean : NAN; - secfilt[Nsec].dFapWrp = Nap > 0 ? apstats->error : NAN; - secfilt[Nsec].sFapWrp = Nap > 0 ? apstats->sigma : NAN; - secfilt[Nsec].NusedApWrp = Nap; - secfilt[Nsec].MapWrp = isnan(secfilt[Nsec].FapWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FapWrp); // 8.9 since flux is in Jy - - liststats (Fkronlist, dkronlist, wkronlist, Nkron, kronstats); - secfilt[Nsec].FkronWrp = Nkron > 0 ? kronstats->mean : NAN; - secfilt[Nsec].dFkronWrp = Nkron > 0 ? kronstats->error : NAN; - secfilt[Nsec].sFkronWrp = Nkron > 0 ? kronstats->sigma : NAN; - secfilt[Nsec].NusedKronWrp = Nkron; - secfilt[Nsec].MkronWrp = isnan(secfilt[Nsec].FkronWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FkronWrp); // 8.9 since flux is in Jy - - secfilt[Nsec].Nwarp = Nwarp; - secfilt[Nsec].NwarpGood = NwarpGood; + secfilt[Nsec].psfQfMax = results->psfQfMax[Nsec]; + secfilt[Nsec].psfQfPerfMax = results->psfQfPerfMax[Nsec]; + + if (isfinite(secfilt[Nsec].psfQfMax) && (secfilt[Nsec].psfQfMax > psfQfMax)) psfQfMax = secfilt[Nsec].psfQfMax; + if (isfinite(secfilt[Nsec].psfQfPerfMax) && (secfilt[Nsec].psfQfPerfMax > psfQfPerfMax)) psfQfPerfMax = secfilt[Nsec].psfQfPerfMax; + } + + // if too few valid measurements meet the minimum criteria, go to the next entry + StatType *psfstats = &results->psfstats; + int Nranking = magStatsByRanking (&results->psfData[Nsec], psfstats); + if (Nranking) { + secfilt[Nsec].FpsfWrp = psfstats->mean; + secfilt[Nsec].dFpsfWrp = psfstats->error; + secfilt[Nsec].sFpsfWrp = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN; + secfilt[Nsec].NusedWrp = psfstats->Nmeas; + secfilt[Nsec].MpsfWrp = isnan(secfilt[Nsec].FpsfWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FpsfWrp); // 8.9 since flux is in Jy + } + int minRank = (Nranking > 0) ? results->psfData[Nsec].ranking[0] : 10; + markMeasureByRanking (&results->psfData[Nsec], measure, minRank, ID_MEAS_WARP_USED); + + // if too few valid measurements meet the minimum criteria, go to the next entry + StatType *apstats = &results->apstats; + Nranking = magStatsByRanking (&results->aperData[Nsec], apstats); + if (Nranking) { + secfilt[Nsec].FapWrp = apstats->mean; + secfilt[Nsec].dFapWrp = apstats->error; + secfilt[Nsec].sFapWrp = (apstats->Nmeas > 1) ? apstats->chisq : NAN; + secfilt[Nsec].NusedApWrp = apstats->Nmeas; + secfilt[Nsec].MapWrp = isnan(secfilt[Nsec].FapWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FapWrp); // 8.9 since flux is in Jy + } + + // if too few valid measurements meet the minimum criteria, go to the next entry + StatType *kronstats = &results->kronstats; + Nranking = magStatsByRanking (&results->kronData[Nsec], kronstats); + if (Nranking) { + secfilt[Nsec].FkronWrp = kronstats->mean; + secfilt[Nsec].dFkronWrp = kronstats->error; + secfilt[Nsec].sFkronWrp = (kronstats->Nmeas > 1) ? kronstats->chisq : NAN; + secfilt[Nsec].NusedKronWrp = kronstats->Nmeas; + secfilt[Nsec].MkronWrp = isnan(secfilt[Nsec].FkronWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FkronWrp); // 8.9 since flux is in Jy + } + + secfilt[Nsec].Nwarp = results->Nmeas[Nsec]; + secfilt[Nsec].NwarpGood = results->NmeasGood[Nsec]; } return (TRUE); } -int setGlobalObjStats (Average *average, Measure *measure) { - - off_t k; - - float psfQFmax = 0.0; - float psfQFperfmax = 0.0; - float stargalmax = 0.0; - - for (k = 0; k < average[0].Nmeasure; k++) { - - if (isfinite(measure[k].psfQF)) { - psfQFmax = MAX (psfQFmax, measure[k].psfQF); - } - if (isfinite(measure[k].psfQFperf)) { - psfQFperfmax = MAX (psfQFperfmax, measure[k].psfQFperf); - } - - // only count psfQF > 0.85 here - if (isfinite(measure[k].extNsigma) && (measure[k].psfQF > 0.85)) { - stargalmax = MAX (stargalmax, measure[k].extNsigma); - } - } - average[0].psfQF = psfQFmax; - average[0].psfQFperf = psfQFperfmax; - average[0].stargal = stargalmax; - - return (TRUE); +int magStatsByRanking (StatDataSet *dataset, StatType *stats) { + + liststats_init (stats); + + if (dataset->Nlist == 0) return 0; + + // we have a list of measurements for this Nsec value, along with errors, weights, and their rank + // first, sort by the rank (increasing) + sort_StatDataSet (dataset); + + // find the values with the same minimum rank: + int Nranking = 0; + int minRank = dataset->ranking[0]; // MIN (5, result->psfData[Nsec].ranking[0]); -- only allow rank 5 or < + + int i; + for (i = 0; (i < dataset->Nlist) && (dataset->ranking[i] == minRank); i++, Nranking++); + + liststats (dataset->flxlist, dataset->errlist, dataset->wgtlist, Nranking, stats); + return (Nranking); } + +int markMeasureByRanking (StatDataSet *dataset, Measure *measure, int minrank, DVOMeasureFlags flags) { + int i; + + for (i = 0; i < dataset->Nlist; i++) { + int k = dataset->measSeq[i]; + if (dataset->ranking[i] > minrank) { + measure[k].dbFlags &= ~flags; + } else { + measure[k].dbFlags |= flags; + } + } + return TRUE; +} + +/* sort a coordinate pair (X,Y) and the associated index (S) */ +void sort_StatDataSet (StatDataSet *dataset) { + +# define SWAPFUNC(A,B){ double dtmp; int itmp; \ + dtmp = dataset->flxlist[A]; dataset->flxlist[A] = dataset->flxlist[B]; dataset->flxlist[B] = dtmp; \ + dtmp = dataset->errlist[A]; dataset->errlist[A] = dataset->errlist[B]; dataset->errlist[B] = dtmp; \ + dtmp = dataset->wgtlist[A]; dataset->wgtlist[A] = dataset->wgtlist[B]; dataset->wgtlist[B] = dtmp; \ + itmp = dataset->ranking[A]; dataset->ranking[A] = dataset->ranking[B]; dataset->ranking[B] = itmp; \ + itmp = dataset->measSeq[A]; dataset->measSeq[A] = dataset->measSeq[B]; dataset->measSeq[B] = itmp; \ + } +# define COMPARE(A,B)(dataset->ranking[A] < dataset->ranking[B]) + + OHANA_SORT (dataset->Nlist, COMPARE, SWAPFUNC); + +# undef SWAPFUNC +# undef COMPARE + +} + Index: /trunk/Ohana/src/relphot/src/setMrelFinal.c =================================================================== --- /trunk/Ohana/src/relphot/src/setMrelFinal.c (revision 38465) +++ /trunk/Ohana/src/relphot/src/setMrelFinal.c (revision 38466) @@ -12,5 +12,5 @@ int Nsecfilt = GetPhotcodeNsecfilt (); - /* if we reset the catalog, reset all the current measurements */ + /*** RESET photometry flags and Mcal values ***/ if (RESET) { // flags used by the photometry analysis (excluding UBERCAL) @@ -46,5 +46,5 @@ for (i = 0; i < catalog[0].Naverage; i++) { off_t N = Nsecfilt*i+Nsec; - dvo_secfilt_init (&catalog[0].secfilt[N]); + dvo_secfilt_init (&catalog[0].secfilt[N], SECFILT_RESET_ALL); off_t m = catalog[0].average[i].measureOffset; @@ -100,8 +100,8 @@ // XXX make this optional? (do not clean for -averages?) if (!simpleAverage) clean_measures (catalog, 1, TRUE, flatcorr); /* mark outliers ID_MEAS_POOR_PHOTOM */ - for (i = 0; i < 5; i++) { - skip_measurements (catalog, i, flatcorr); /* set ID_MEAS_SKIP for measures to be skipped */ - setMrelOutput (catalog, 1, i, flatcorr); /* set Mrel using remaining measures */ - } + + ALLOCATE (catalog[0].measureRank, char, catalog[0].Nmeasure); + setMeasureRank (catalog); + setMrelOutput (catalog, 1, flatcorr); if (!simpleAverage) setMcalOutput (catalog, 1, flatcorr); @@ -267,14 +267,79 @@ } -/* dophot types: - -1 - star -2 - galaxy -3 - star -4 - weak -5 - edge -7 - weak -9 - weak -10 - sat -*/ - +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; + + /* set measureRank[] based on various quality measurements */ + for (i = 0; i < catalog[0].Nmeasure; i++) { + 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 + PhotCode *code = GetPhotcodebyCode (measure[i].photcode); + if (!code) continue; + + // 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; } + } + + // measurements from BAD images and mosaics (not REF mags) + 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; } + } + } + + // 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; } + + if (Nim > -1) { + if (imageFlags & IMAGE_POOR) { measureRank[i] = 5; continue; } + if (mosaicFlags & IMAGE_POOR) { measureRank[i] = 5; continue; } + } + + // RANK 2 : psfQF value + if (!isfinite(measure[i].psfQF) || measure[i].psfQF < 0.85) { measureRank[i] = 4; continue; } + + // RANK 1 : POOR photFlags + if (measure[i].dbFlags & MEAS_POOR) { measureRank[i] = 3; continue; } + + // RANK 1 : POOR photFlags + if (measure[i].photFlags & code->photomPoorMask) { measureRank[i] = 2; continue; } + + // RANK 1 : psfQFperfect value + if (!isfinite(measure[i].psfQFperf) || measure[i].psfQFperf < 0.85) { measureRank[i] = 1; continue; } + + // RANK 0 : perfect measurement: + measureRank[i] = 0; + } +} Index: /trunk/Ohana/src/relphot/src/synthetic_zpts.c =================================================================== --- /trunk/Ohana/src/relphot/src/synthetic_zpts.c (revision 38465) +++ /trunk/Ohana/src/relphot/src/synthetic_zpts.c (revision 38466) @@ -7,10 +7,12 @@ static char *extname[5] = {"map_g", "map_r", "map_i", "map_z", "map_y"}; -SynthZeroPoints *SynthZeroPointsLoad (char *filename) { +static SynthZeroPoints *zpts = NULL; + +int SynthZeroPointsLoad (char *filename) { FILE *f = fopen (filename, "r"); if (!f) { fprintf (stderr, "ERROR: cannot open file %s\n", filename); - return NULL; + return FALSE; } @@ -23,5 +25,4 @@ // map_y - SynthZeroPoints *zpts = NULL; ALLOCATE (zpts, SynthZeroPoints, 1); @@ -29,5 +30,5 @@ if (VERBOSE) fprintf (stderr, "can't read header\n"); fclose (f); - return NULL; + return FALSE; } @@ -49,5 +50,8 @@ zpts->Ny = zpts->matrix[0].Naxis[1]; + return TRUE; +} + +SynthZeroPoints *SynthZeroPointsGet () { return zpts; } -