Index: /branches/eam_branches/ipp-20150419/Ohana/src/relphot/include/relphot.h =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/relphot/include/relphot.h (revision 38313) +++ /branches/eam_branches/ipp-20150419/Ohana/src/relphot/include/relphot.h (revision 38314) @@ -137,4 +137,5 @@ 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 Index: /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/ImageOps.c =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/ImageOps.c (revision 38313) +++ /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/ImageOps.c (revision 38314) @@ -692,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]; @@ -722,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: /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/StarOps.c =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/StarOps.c (revision 38313) +++ /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/StarOps.c (revision 38314) @@ -112,4 +112,5 @@ results->Nmeas[i] = 0; + results->NmeasGood[i] = 0; results->Next[i] = 0; results->NexpPS1[i] = 0; @@ -149,4 +150,5 @@ ALLOCATE (results->Nmeas, int, results->Nsecfilt); + ALLOCATE (results->NmeasGood, int, results->Nsecfilt); ALLOCATE (results->Next, int, results->Nsecfilt); ALLOCATE (results->NexpPS1, int, results->Nsecfilt); @@ -1240,4 +1242,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: /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/bcatalog.c =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/bcatalog.c (revision 38313) +++ /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/bcatalog.c (revision 38314) @@ -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++) { Index: /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/plotstuff.c =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/plotstuff.c (revision 38313) +++ /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/plotstuff.c (revision 38314) @@ -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: /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/setMrelCatalog.c =================================================================== --- /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/setMrelCatalog.c (revision 38313) +++ /branches/eam_branches/ipp-20150419/Ohana/src/relphot/src/setMrelCatalog.c (revision 38314) @@ -1,3 +1,5 @@ # include "relphot.h" +int markMeasureByRanking (StatDataSet *dataset, Measure *measure, int minrank, DVOMeasureFlags flags); + # define UBERCAL_WEIGHT 100.0 void sort_StatDataSet (StatDataSet *dataset); @@ -11,4 +13,5 @@ # define CHECK_VALID_MAG(MAG) (!isnan(MAG) && (MAG > 0.0) && (MAG < 30.0)) +# define CHECK_VALID_FLUX(MAG) (!isnan(FLUX)) # if (0) @@ -318,4 +321,11 @@ if (isSetMrelFinal) { + if ((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) { @@ -341,4 +351,5 @@ secfilt[Nsec].NusedAp = Nranking; } + markMeasureByRanking (&results->aperData[Nsec], measure, minRank, ID_MEAS_PHOTOM_APER); StatType *kronstats = &results->kronstats; @@ -350,4 +361,5 @@ secfilt[Nsec].NusedKron = Nranking; } + markMeasureByRanking (&results->kronData[Nsec], measure, minRank, ID_MEAS_PHOTOM_KRON); // does this object appear extended in > 50% of measurements? @@ -629,10 +641,4 @@ } -# 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: @@ -640,6 +646,4 @@ int setMrelAverageForcedWarp (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr, SetMrelInfo *results) { -# if (0) - off_t measureOffset = average->measureOffset; @@ -652,23 +656,4 @@ 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; - // option for a test print if ((average[0].objID == 0xd6e) && (average[0].catID == 0x4984)) { @@ -683,178 +668,157 @@ } - 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; - + for (Nsec = 0; Nsec < Nsecfilt; Nsec++) { // 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; - + 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; - } + if (!FindImageSkycellID (secfilt[Nsec].stackBestOff, cat, &stkTessID[Nsec], &stkProjID[Nsec], &stkSkycellID[Nsec])) continue; + } + } + + // assign measurements to the photcode lists + for (k = 0; k < average[0].Nmeasure; k++) { + + // 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->NwarpGood[Nsec] ++; + + measure[k].dbFlags &= ~ID_MEAS_WARP_USED; + + // use primary skycell for DIFF, stack skycell for WARP + if (IS_DIFF_DB) { + if (MatchImageSkycellID (meas, cat, priTessID, priProjID, priSkycellID)) { + measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY; } 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; - continue; - } - - found[Nsec] = TRUE; - - liststats (Fpsflist, dpsflist, wpsflist, Npsf, psfstats); + continue; + } + } else { + if (!MatchImageSkycellID (meas, cat, stkTessID, stkProjID, stkSkycellID)) continue; + } + + // assign the max psfQF values for each Nsec value + if (IS_DIFF_DB) { + 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); + + // XX measure[k].dbFlags &= ~ID_MEAS_WARP_USED; + // XXX measure[k].dbFlags |= ID_MEAS_WARP_USED; + + // 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)) { + 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)) { + 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)) { + int Nkron = results->kronData[Nsec].Nlist; + results->kronData[Nsec].flxlist[Nkron] = Fap * Fcal; + results->kronData[Nsec].errlist[Nkron] = dFap * 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 ++; + } + } + + // now calculate the mean stats for the Nsec bands. + for (Nsec = 0; Nsec < Nsecfilt; Nsec++) { + // XXX cannot do this -- resets psf mags + // dvo_secfilt_init (&secfilt[Nsec]); // XXX does this reset astrometry bits? no, right? 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 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, IID_MEAS_WARP_USED); + + // if too few valid measurements meet the minimum criteria, go to the next entry + StatType *aperstats = &results->aperstats; + Nranking = magStatsByRanking (&results->aperData[Nsec], aperstats); + if (Nranking) { + secfilt[Nsec].FapWrp = aperstats->mean; + secfilt[Nsec].dFapWrp = aperstats->error; + secfilt[Nsec].sFapWrp = (aperstats->Nmeas > 1) ? aperstats->chisq : NAN; + secfilt[Nsec].NusedApWrp = aperstats->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]; } setGlobalObjStats (average, measure); - -# endif return (TRUE); @@ -909,4 +873,18 @@ 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; }