Index: /branches/eam_branches/ipp-20111122/Ohana/src/relphot/include/relphot.h =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relphot/include/relphot.h (revision 33156) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relphot/include/relphot.h (revision 33157) @@ -79,4 +79,5 @@ int FREEZE_MOSAICS; int USE_GRID; +int KEEP_UBERCAL; char *OUTROOT; int PLOTDELAY; @@ -157,4 +158,5 @@ off_t getImageEntry PROTO((off_t meas, int cat)); float getMcal PROTO((off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog)); +float getMflat PROTO((off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog)); float getMgrid PROTO((off_t meas, int cat)); float getMmos PROTO((off_t meas, int cat)); @@ -196,5 +198,5 @@ int reload_images PROTO((FITS_DB *db)); int setExclusions PROTO((Catalog *catalog, int Ncatalog)); -void setMcal PROTO((Catalog *catalog, int Poor)); +void setMcal PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr)); void setMcalFinal PROTO((void)); int setMcalOutput PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr)); Index: /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/ImageOps.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/ImageOps.c (revision 33156) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/ImageOps.c (revision 33157) @@ -1,3 +1,4 @@ # include "relphot.h" +void plot_setMcal (double *list, int Npts, StatType *stats, float clouds); // we have an array of images: image[ImageIndex] (ImageIndex : 0 < Nimage) @@ -253,4 +254,5 @@ } +// returns image.Mcal - ff(x,y) float getMcal (off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog) { @@ -268,9 +270,25 @@ int flat_id = image[i].photom_map_id; if (flat_id) { - offset = FlatCorrectionOffset (flatcorr, flat_id, catalog[cat].measure[meas].Xccd, catalog[cat].measure[meas].Yccd); - } - value += offset; + offset = FlatCorrectionOffset (flatcorr, flat_id, catalog[cat].measureT[meas].Xccd, catalog[cat].measureT[meas].Yccd); + } + value -= offset; return (value); +} + +float getMflat (off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog) { + + off_t i = MeasureToImage[cat][meas]; + if (i == -1) return (NAN); + + float offset = 0.0; + + // to do this, I need to pass in the catalog and flatcorr pointers + int flat_id = image[i].photom_map_id; + if (flat_id) { + offset = FlatCorrectionOffset (flatcorr, flat_id, catalog[cat].measureT[meas].Xccd, catalog[cat].measureT[meas].Yccd); + } + + return (offset); } @@ -285,9 +303,9 @@ /* determine Mcal values for all images */ -void setMcal (Catalog *catalog, int PoorImages) { +void setMcal (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) { off_t i, j, m, c, n, N, Nmax, mark, bad; int Nfew, Nbad, Nmos, Ngrid, Nrel, Nsys, Nbright; - float Msys, Mrel, Mmos, Mgrid, McalBright, McalBright2; + float Msys, Mrel, Mmos, Mgrid, Mflat, McalBright, McalBright2; double *list, *dlist; StatType stats; @@ -321,4 +339,7 @@ if (!bad) continue; } + + /* we optionally do not recalibrate images with UBERCAL zero points */ + if (KEEP_UBERCAL && (image[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue; N = 0; @@ -348,4 +369,11 @@ } + // image.Mcal is not supposed to include the flat-field correction, so we need to + // apply that offset as well here for this image (in other words, each detection is + // being compared to the model, excluding the zero point, Mcal. The model includes + // the flat-correction. NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat) + + Mflat = getMflat (m, c, flatcorr, catalog); + n = catalog[c].measureT[m].averef; Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]); @@ -354,5 +382,8 @@ continue; } - list[N] = Msys - Mrel - Mmos - Mgrid; + + fprintf (stderr, "%1d, %3d : %3d, %3d : %10.6f %10.6f : %6.3f %6.3f %6.3f %6.3f %6.3f\n", (int) i, (int) j, (int) c, (int) m, catalog[c].averageT[n].R, catalog[c].averageT[n].D, Msys, Mrel, Mmos, Mgrid, Mflat); + + list[N] = Msys - Mrel - Mmos - Mgrid + Mflat; dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) { @@ -385,4 +416,6 @@ image[i].Xm = 100.0*log10(stats.chisq); + plot_setMcal (list, N, &stats, CLOUD_TOLERANCE); + if (image[i].Mcal < -CLOUD_TOLERANCE) { image[i].Mcal = 0.0; @@ -454,4 +487,36 @@ } +void plot_setMcal (double *list, int Npts, StatType *stats, float clouds) { + + off_t i; + double *xlist; + Graphdata graphdata; + + if (!PLOTSTUFF) return; + + ALLOCATE (xlist, double, Npts); + + /**** Mcal vs seq ****/ + float minMcal = +100.0; + float maxMcal = -100.0; + for (i = 0; i < Npts; i++) { + xlist[i] = i; + minMcal = MIN (list[i], minMcal); + maxMcal = MAX (list[i], maxMcal); + } + + float McalRange = MAX(1.2*(maxMcal - minMcal), 0.21); + float McalCenter = 0.5*(maxMcal + minMcal); + + plot_defaults (&graphdata); + graphdata.xmin = -1; + graphdata.xmax = Npts + 1; + graphdata.ymin = McalCenter - 0.5*McalRange; + graphdata.ymax = McalCenter + 0.5*McalRange; + plot_list (&graphdata, xlist, list, Npts, "sequence vs Mcal", "%s.seq.png", OUTROOT); + + free (xlist); +} + void plot_images () { @@ -467,15 +532,43 @@ /**** dMcal vs airmass ****/ + float minAirmass = 1000.0; + float maxAirmass = 0.0; + float minMcal = +100.0; + float maxMcal = -100.0; + float mindMcal = +100.0; + float maxdMcal = -100.0; for (i = 0; i < Nimage; i++) { Mlist[i] = image[i].Mcal; dlist[i] = image[i].dMcal; xlist[i] = image[i].secz; - } + minAirmass = MIN (image[i].secz, minAirmass); + maxAirmass = MAX (image[i].secz, maxAirmass); + minMcal = MIN (image[i].Mcal, minMcal); + maxMcal = MAX (image[i].Mcal, maxMcal); + mindMcal = MIN (image[i].dMcal, mindMcal); + maxdMcal = MAX (image[i].dMcal, maxdMcal); + } + + float AirmassRange = MAX(1.2*(maxAirmass - minAirmass), 0.25); + float AirmassCenter = 0.5*(maxAirmass + minAirmass); + + float McalRange = MAX(1.2*(maxMcal - minMcal), 0.21); + float McalCenter = 0.5*(maxMcal + minMcal); + + float dMcalRange = MAX(1.2*(maxdMcal - mindMcal), 0.21); + float dMcalCenter = 0.5*(maxdMcal + mindMcal); plot_defaults (&graphdata); - graphdata.ymin = PlotdMmin; - graphdata.ymax = PlotdMmax; + graphdata.xmin = AirmassCenter - 0.5*AirmassRange; + graphdata.xmax = AirmassCenter + 0.5*AirmassRange; + graphdata.ymin = McalCenter - 0.5*McalRange; + graphdata.ymax = McalCenter + 0.5*McalRange; plot_list (&graphdata, xlist, Mlist, Nimage, "airmass vs Mcal", "%s.airmass.png", OUTROOT); + plot_defaults (&graphdata); + graphdata.xmin = McalCenter - 0.5*McalRange; + graphdata.xmax = McalCenter + 0.5*McalRange; + graphdata.ymin = dMcalCenter - 0.5*dMcalRange; + graphdata.ymax = dMcalCenter + 0.5*dMcalRange; plot_list (&graphdata, Mlist, dlist, Nimage, "Mcal vs dMcal", "%s.Mcal.dMcal.png", OUTROOT); @@ -495,4 +588,6 @@ plot_defaults (&graphdata); graphdata.style = 1; + graphdata.xmin = dMcalCenter - 0.5*dMcalRange; + graphdata.xmax = dMcalCenter + 0.5*dMcalRange; plot_list (&graphdata, xlist, Mlist, NBIN, "dMcal hist", "%s.dMcalhist.png", OUTROOT); Index: /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/MosaicOps.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/MosaicOps.c (revision 33156) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/MosaicOps.c (revision 33157) @@ -102,4 +102,6 @@ mosaic[Nmosaic].secz = image[i].secz; mosaic[Nmosaic].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode); + + // XXX do we need to do something about flag consistency across a mosaic? /* add image to mosaic image list */ @@ -175,4 +177,12 @@ dMcal += image[m].dMcal; Xm += image[m].Xm; + + // XXX: how does this work with UBERCAL? We want to keep the Mcal values supplied by ubercal, but + // solve for a single offset for each exposure (Mosaic.Mcal). + + // we also want to keep the flat-field terms for each exposure (regardless of ubercal or not) + + // if it helps, note that ubercal uses a single zp per exposure, so the mean of those values is the same as the value + /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */ image[m].Mcal = 0.0; @@ -215,4 +225,7 @@ image = getimages (&Nimage, NULL); + // XXX I think this is OK in the ubercal context, but probably need to skip UBERCAL + // images? (no need to update them) + // copy the mosaic results to the images. set the mosaic Mcal to 0.0 since we have moved its // impact to the images @@ -380,5 +393,5 @@ off_t i, j, m, c, N, Nmax; int n, mark, bad, Nfew, Nbad, Ncal, Nrel, Ngrid, Nsys, Nsecfilt; - float Msys, Mrel, Mcal, Mgrid; + float Msys, Mrel, Mcal, Mgrid, Mflat; double *list, *dlist, *Mlist, *dMlist; StatType stats; @@ -414,4 +427,7 @@ if (!bad) continue; } + + /* we optionally do not recalibrate images with UBERCAL zero points */ + if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue; N = 0; @@ -441,4 +457,11 @@ } + // image.Mcal is not supposed to include the flat-field correction, so we need to + // apply that offset as well here for this image (in other words, each detection is + // being compared to the model, excluding the zero point, Mcal. The model includes + // the flat-correction. NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat) + + Mflat = getMflat (m, c, flatcorr, catalog); + n = catalog[c].measureT[m].averef; Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]); @@ -447,5 +470,5 @@ continue; } - list[N] = Msys - Mrel - Mcal - Mgrid; + list[N] = Msys - Mrel - Mcal - Mgrid + Mflat; dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); Mlist[N] = Msys; @@ -499,4 +522,5 @@ // to 0.0. At the same time, we make a guess to the effective impact on all other images, // driven by the coupling of common stars. +// XXX this function is no longer used because we force significantly negative clouds to 0.0 int rationalize_mosaics (Catalog *catalog, int Ncatalog) { Index: /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/StarOps.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/StarOps.c (revision 33156) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/StarOps.c (revision 33157) @@ -60,4 +60,6 @@ for (j = 0; j < catalog[i].Naverage; j++) { + // XXX accumulate all secfilt values in a single pass? + int Ns; for (Ns = 0; Ns < Nphotcodes; Ns++) { @@ -82,4 +84,7 @@ continue; } + + // XXX include the ubercal measurements? yes, of course! these are the well-determined values, supposedly! + // XXX allow REF stars (no Image Entry) to be included in the calculation this // should be optionally set, and should allow for REF stars to be downweighted by @@ -144,4 +149,8 @@ } +// setMrel and setMrelOutput are extremely similar, but have slightly different implications: +// * setMrel uses the internal Tiny structures only +// * setMrelOutput skips stars for which there are too few good measurements +// * setMrelOutput is meant to be called repeatedly, relaxing the criteria for 'good' on each pass int setMrelOutput (Catalog *catalog, int Ncatalog, int mark, FlatCorrectionTable *flatcorr) { @@ -167,4 +176,5 @@ for (i = 0; i < Ncatalog; i++) { for (j = 0; j < catalog[i].Naverage; j++) { + /* skip stars already calibrated */ if (catalog[i].found[j]) continue; @@ -172,4 +182,5 @@ int Ns; for (Ns = 0; Ns < Nphotcodes; Ns++) { + int thisCode = photcodes[Ns][0].code; Nsec = GetPhotcodeNsec(thisCode); Index: /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/args.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/args.c (revision 33156) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/args.c (revision 33157) @@ -189,4 +189,10 @@ } + KEEP_UBERCAL = FALSE; + if ((N = get_argument (argc, argv, "-keep-ubercal"))) { + remove_argument (N, &argc, argv); + KEEP_UBERCAL = TRUE; + } + MIN_ERROR = 0.001; if ((N = get_argument (argc, argv, "-minerror"))) { Index: /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/help.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/help.c (revision 33156) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/help.c (revision 33157) @@ -37,4 +37,5 @@ fprintf (stderr, " -imfreeze\n"); fprintf (stderr, " -grid\n"); + fprintf (stderr, " -keep-ubercal : do not change zero point for ubercal-ed images\n"); fprintf (stderr, " -area Xmin Xmax Ymin Ymax\n"); fprintf (stderr, " -instmag min max\n"); Index: /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/plot_scatter.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/plot_scatter.c (revision 33156) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/plot_scatter.c (revision 33157) @@ -11,5 +11,5 @@ for (i = 0; i < Ncatalog; i++) { for (j = 0; j < catalog[i].Naverage; j++) { - Ntot += catalog[i].average[j].Nmeasure; + Ntot += catalog[i].averageT[j].Nmeasure; } } @@ -32,7 +32,7 @@ /* calculate the average value for a single star */ if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; - m = catalog[i].average[j].measureOffset; + m = catalog[i].averageT[j].measureOffset; - for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) { + for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) { // skip measurements that do not match the current photcode int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode); Index: /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/plotstuff.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/plotstuff.c (revision 33156) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/plotstuff.c (revision 33157) @@ -165,5 +165,5 @@ graphdata[0].color = black; graphdata[0].lweight = 0; - graphdata[0].size = 0.5; + graphdata[0].size = 1.0; graphdata[0].xmin = dUNDEF; Index: /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/relphot.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/relphot.c (revision 33156) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/relphot.c (revision 33157) @@ -120,6 +120,12 @@ /* determine fit values */ for (i = 0; i < NLOOP; i++) { + if (PLOTSTUFF) { + plot_scatter (catalog, Ncatalog, flatcorr); + } setMrel (catalog, Ncatalog, flatcorr); - setMcal (catalog, FALSE); + if (PLOTSTUFF) { + plot_scatter (catalog, Ncatalog, flatcorr); + } + setMcal (catalog, FALSE, flatcorr); setMmos (catalog, FALSE, flatcorr); setMgrid (catalog, flatcorr); @@ -160,5 +166,5 @@ /* set Mcal & Mmos for bad images */ - setMcal (catalog, TRUE); + setMcal (catalog, TRUE, flatcorr); setMmos (catalog, TRUE, flatcorr); MARKTIME("-- finalize Mcal values: %f sec\n", dtime);