Index: /branches/eam_branches/ipp-20211108/Ohana/src/relphot/Makefile =================================================================== --- /branches/eam_branches/ipp-20211108/Ohana/src/relphot/Makefile (revision 42072) +++ /branches/eam_branches/ipp-20211108/Ohana/src/relphot/Makefile (revision 42073) @@ -54,4 +54,6 @@ $(SRC)/save_images.$(ARCH).o \ $(SRC)/MagResidSave.$(ARCH).o \ +$(SRC)/setXradAverages.$(ARCH).o \ +$(SRC)/lensing.$(ARCH).o \ $(SRC)/synthetic_mags.$(ARCH).o \ $(SRC)/plot_scatter.$(ARCH).o \ @@ -117,4 +119,6 @@ $(SRC)/save_images.$(ARCH).o \ $(SRC)/MagResidSave.$(ARCH).o \ +$(SRC)/setXradAverages.$(ARCH).o \ +$(SRC)/lensing.$(ARCH).o \ $(SRC)/load_catalogs.$(ARCH).o \ $(SRC)/reload_catalogs.$(ARCH).o \ Index: /branches/eam_branches/ipp-20211108/Ohana/src/relphot/include/lensing.h =================================================================== --- /branches/eam_branches/ipp-20211108/Ohana/src/relphot/include/lensing.h (revision 42073) +++ /branches/eam_branches/ipp-20211108/Ohana/src/relphot/include/lensing.h (revision 42073) @@ -0,0 +1,33 @@ + +// As a temporary hack (for UNIONS DR3), I am going to overload the following fields with +// radial aperture values from the convolved images: + +# define F_ApR5_C1 X11_sm_obj +# define dF_ApR5_C1 E1_sm_obj +# define F_ApR6_C1 X22_sm_obj +# define dF_ApR6_C1 E2_sm_obj +# define F_ApR7_C1 X11_sh_obj +# define dF_ApR7_C1 E2_sh_obj + +# define F_ApR5_C2 X11_sm_psf +# define dF_ApR5_C2 E1_sm_psf +# define F_ApR6_C2 X22_sm_psf +# define dF_ApR6_C2 E2_sm_psf +# define F_ApR7_C2 X11_sh_psf +# define dF_ApR7_C2 E2_sh_psf + + +// counter to track number of valid measurements for each radius +typedef struct { + int N5_C0; int N5_C1; int N5_C2; + int N6_C0; int N6_C1; int N6_C2; + int N7_C0; int N7_C1; int N7_C2; + int Nmeas; +} Lensctr; + +Lensctr *dvo_lensctr_init (int Nsec); +int dvo_lensctr_reset (Lensctr *lensctr, int Nsec); +int dvo_lensing_accum (Lensobj *lensobj, Lensctr *lensctr, Lensing *lensing, float Fcal); +int dvo_lensctr_has_values (Lensctr *lensctr); +int dvo_lensobj_stat (float *mean, float *error, float *stdev, float *fill, int count); +int dvo_lensobj_aves (Lensobj *lensobj, Lensctr *lensctr); Index: /branches/eam_branches/ipp-20211108/Ohana/src/relphot/include/relphot.h =================================================================== --- /branches/eam_branches/ipp-20211108/Ohana/src/relphot/include/relphot.h (revision 42072) +++ /branches/eam_branches/ipp-20211108/Ohana/src/relphot/include/relphot.h (revision 42073) @@ -486,4 +486,5 @@ char *OUTROOT; char *UPDATE_CATFORMAT; +int UPDATE_XRAD; int PLOTDELAY; int UpdateAverages; @@ -877,4 +878,5 @@ int MagResidSave(char *filename, Catalog *catalog); +int setXradAverages (Catalog *catalog); void ResetAverageAndMeasure (Catalog *catalog); Index: /branches/eam_branches/ipp-20211108/Ohana/src/relphot/src/args.c =================================================================== --- /branches/eam_branches/ipp-20211108/Ohana/src/relphot/src/args.c (revision 42072) +++ /branches/eam_branches/ipp-20211108/Ohana/src/relphot/src/args.c (revision 42073) @@ -259,4 +259,10 @@ UPDATE_CATFORMAT = strcreate (argv[N]); remove_argument (N, &argc, argv); + } + + UPDATE_XRAD = FALSE; + if ((N = get_argument (argc, argv, "-update-xrad-average"))) { + remove_argument (N, &argc, argv); + UPDATE_XRAD = TRUE; } Index: /branches/eam_branches/ipp-20211108/Ohana/src/relphot/src/lensing.c =================================================================== --- /branches/eam_branches/ipp-20211108/Ohana/src/relphot/src/lensing.c (revision 42073) +++ /branches/eam_branches/ipp-20211108/Ohana/src/relphot/src/lensing.c (revision 42073) @@ -0,0 +1,151 @@ +# include "relphot.h" +# include "lensing.h" + +Lensctr *dvo_lensctr_init (int Nsec) { + + ALLOCATE_PTR (lensctr, Lensctr, Nsec); + return (lensctr); +} + +int dvo_lensctr_reset (Lensctr *lensctr, int Nsec) { + + for (int i = 0; i < Nsec; i++) { + memset (&lensctr[i], 0, sizeof (Lensctr)); + } + return TRUE; +} + +int dvo_lensing_accum (Lensobj *lensobj, Lensctr *lensctr, Lensing *lensing, float Fcal) { + + int hasLensing = FALSE; + + if (isfinite(lensing-> F_ApR5)) { + lensobj-> F_ApR5 += Fcal * lensing-> F_ApR5; + lensobj->dF_ApR5 += SQ(Fcal * lensing->dF_ApR5); + lensobj->sF_ApR5 += SQ(Fcal * lensing-> F_ApR5); + lensobj->fF_ApR5 += lensing->fF_ApR5; + lensctr->N5_C0 ++; + hasLensing = TRUE; + } + + if (isfinite(lensing-> F_ApR6)) { + lensobj-> F_ApR6 += Fcal * lensing-> F_ApR6; + lensobj->dF_ApR6 += SQ(Fcal * lensing->dF_ApR6); + lensobj->sF_ApR6 += SQ(Fcal * lensing-> F_ApR6); + lensobj->fF_ApR6 += lensing->fF_ApR6; + lensctr->N6_C0 ++; + hasLensing = TRUE; + } + + if (isfinite(lensing-> F_ApR7)) { + lensobj-> F_ApR7 += Fcal * lensing-> F_ApR7; + lensobj->dF_ApR7 += SQ(Fcal * lensing->dF_ApR7); + lensobj->sF_ApR7 += SQ(Fcal * lensing-> F_ApR7); + lensobj->fF_ApR7 += lensing->fF_ApR7; + lensctr->N7_C0 ++; + hasLensing = TRUE; + } + + if (isfinite(lensing-> F_ApR5_C1)) { + lensobj-> F_ApR5_C1 += Fcal * lensing-> F_ApR5_C1; + lensobj->dF_ApR5_C1 += SQ(Fcal * lensing->dF_ApR5_C1); + lensctr->N5_C1 ++; + hasLensing = TRUE; + } + if (isfinite(lensing-> F_ApR6_C1)) { + lensobj-> F_ApR6_C1 += Fcal * lensing-> F_ApR6_C1; + lensobj->dF_ApR6_C1 += SQ(Fcal * lensing->dF_ApR6_C1); + lensctr->N6_C1 ++; + hasLensing = TRUE; + } + if (isfinite(lensing-> F_ApR7_C1)) { + lensobj-> F_ApR7_C1 += Fcal * lensing-> F_ApR7_C1; + lensobj->dF_ApR7_C1 += SQ(Fcal * lensing->dF_ApR7_C1); + lensctr->N7_C1 ++; + hasLensing = TRUE; + } + + if (isfinite(lensing-> F_ApR5_C2)) { + lensobj-> F_ApR5_C2 += Fcal * lensing-> F_ApR5_C2; + lensobj->dF_ApR5_C2 += SQ(Fcal * lensing->dF_ApR5_C2); + lensctr->N5_C2 ++; + hasLensing = TRUE; + } + if (isfinite(lensing-> F_ApR6_C2)) { + lensobj-> F_ApR6_C2 += Fcal * lensing-> F_ApR6_C2; + lensobj->dF_ApR6_C2 += SQ(Fcal * lensing->dF_ApR6_C2); + lensctr->N6_C2 ++; + hasLensing = TRUE; + } + if (isfinite(lensing-> F_ApR7_C2)) { + lensobj-> F_ApR7_C2 += Fcal * lensing-> F_ApR7_C2; + lensobj->dF_ApR7_C2 += SQ(Fcal * lensing->dF_ApR7_C2); + lensctr->N7_C2 ++; + hasLensing = TRUE; + } + if (hasLensing) { lensctr->Nmeas ++; } + return TRUE; +} + +// return TRUE if any counter is non-zero: +int dvo_lensctr_has_values (Lensctr *lensctr) { + if (lensctr->N5_C0) return TRUE; + if (lensctr->N6_C0) return TRUE; + if (lensctr->N7_C0) return TRUE; + + if (lensctr->N5_C1) return TRUE; + if (lensctr->N6_C1) return TRUE; + if (lensctr->N7_C1) return TRUE; + + if (lensctr->N5_C2) return TRUE; + if (lensctr->N6_C2) return TRUE; + if (lensctr->N7_C2) return TRUE; + + return FALSE; +} + +int dvo_lensobj_stat (float *mean, float *error, float *stdev, float *fill, int count) { + + if (count) { + *mean /= (float) count; + *error = sqrt(*error / (float) count); + + if (fill) { *fill /= (float) count; } + + if (stdev) { + if (count < 2) { + *stdev = NAN; + } else { + double S1 = SQ(*mean); // ^2 + double S2 = *stdev / (float) count; // sum(f^2) / N + *stdev = sqrt(S2 - S1) * (count / (count - 1.0)); // correct to sample stdev + } + } + } else { + *mean = NAN; + *error = NAN; + if (stdev) { *stdev = NAN; } + if (fill) { *fill = NAN; } + } + return TRUE; +} + +int dvo_lensobj_aves (Lensobj *lensobj, Lensctr *lensctr) { + + dvo_lensobj_stat (&lensobj-> F_ApR5, &lensobj->dF_ApR5, &lensobj->sF_ApR5, &lensobj->fF_ApR5, lensctr->N5_C0); + dvo_lensobj_stat (&lensobj-> F_ApR6, &lensobj->dF_ApR6, &lensobj->sF_ApR6, &lensobj->fF_ApR6, lensctr->N6_C0); + dvo_lensobj_stat (&lensobj-> F_ApR7, &lensobj->dF_ApR7, &lensobj->sF_ApR7, &lensobj->fF_ApR7, lensctr->N7_C0); + + dvo_lensobj_stat (&lensobj-> F_ApR5_C1, &lensobj->dF_ApR5_C1, NULL, NULL, lensctr->N5_C1); + dvo_lensobj_stat (&lensobj-> F_ApR6_C1, &lensobj->dF_ApR6_C1, NULL, NULL, lensctr->N6_C1); + dvo_lensobj_stat (&lensobj-> F_ApR7_C1, &lensobj->dF_ApR7_C1, NULL, NULL, lensctr->N7_C1); + + dvo_lensobj_stat (&lensobj-> F_ApR5_C2, &lensobj->dF_ApR5_C2, NULL, NULL, lensctr->N5_C2); + dvo_lensobj_stat (&lensobj-> F_ApR6_C2, &lensobj->dF_ApR6_C2, NULL, NULL, lensctr->N6_C2); + dvo_lensobj_stat (&lensobj-> F_ApR7_C2, &lensobj->dF_ApR7_C2, NULL, NULL, lensctr->N7_C2); + + lensobj->Nmeas = lensctr->Nmeas; + + return TRUE; +} + Index: /branches/eam_branches/ipp-20211108/Ohana/src/relphot/src/reload_catalogs.c =================================================================== --- /branches/eam_branches/ipp-20211108/Ohana/src/relphot/src/reload_catalogs.c (revision 42072) +++ /branches/eam_branches/ipp-20211108/Ohana/src/relphot/src/reload_catalogs.c (revision 42073) @@ -55,7 +55,11 @@ catalog.catflags = DVO_LOAD_AVERAGE | DVO_LOAD_MEASURE | DVO_LOAD_MISSING | DVO_LOAD_SECFILT; + // if we want to update the average XRAD values in lensobj, need to load lensobj and lensing + if (UPDATE_XRAD || UPDATE_CATFORMAT) { + catalog.catflags |= DVO_LOAD_LENSING | DVO_LOAD_LENSOBJ; + } // if we are going to update the format, we should update all tables if (UPDATE_CATFORMAT) { - catalog.catflags |= DVO_LOAD_LENSING | DVO_LOAD_LENSOBJ | DVO_LOAD_STARPAR | DVO_LOAD_GALPHOT; + catalog.catflags |= DVO_LOAD_STARPAR | DVO_LOAD_GALPHOT; } catalog.Nsecfilt = GetPhotcodeNsecfilt (); // set the desired number in case we need to create the catalog @@ -87,4 +91,8 @@ setMrelFinal (&catalog, FALSE); TIMESTAMP(time6); + + if (UPDATE_XRAD) { + setXradAverages (&catalog); + } // modify the output format as desired (ignore current format on disk) Index: /branches/eam_branches/ipp-20211108/Ohana/src/relphot/src/relphot_images.c =================================================================== --- /branches/eam_branches/ipp-20211108/Ohana/src/relphot/src/relphot_images.c (revision 42072) +++ /branches/eam_branches/ipp-20211108/Ohana/src/relphot/src/relphot_images.c (revision 42073) @@ -133,4 +133,7 @@ // if (GRID_ZEROPT) dump_grid (); + setMcal (catalog); + setMmos (catalog); + setMgrp (catalog); MARKTIME("-- finalize Mcal values: %f sec\n", dtime); Index: /branches/eam_branches/ipp-20211108/Ohana/src/relphot/src/setXradAverages.c =================================================================== --- /branches/eam_branches/ipp-20211108/Ohana/src/relphot/src/setXradAverages.c (revision 42073) +++ /branches/eam_branches/ipp-20211108/Ohana/src/relphot/src/setXradAverages.c (revision 42073) @@ -0,0 +1,153 @@ +# include "relphot.h" +# include "lensing.h" + +# define SCALE 0.001 +static float MagToFlux (float Mag) { + float Flux = pow(10.0, -0.4*(Mag)); + return (Flux); +} + +// This function simultaneously generates the lensobj table entries for each object +// with lensing measurements and calculates the average values +// NOTE: 'lensing' and 'lensobj' names are archaic : they really hold xrad values + +int setXradAverages (Catalog *catalog) { + + if (VERBOSE2) fprintf (stderr, "lensobj catalog %d : "OFF_T_FMT" ave, "OFF_T_FMT" meas, "OFF_T_FMT" lensing\n", catalog->catID, catalog->Naverage, catalog->Nmeasure, catalog->Nlensing); + + int Nsecfilt = GetPhotcodeNsecfilt (); + + // save the photcodes for the Nsec values (assigned to lensobj below) + ALLOCATE_PTR (photcodeVals, int, Nsecfilt); + for (int Nsec = 0; Nsec < Nsecfilt; Nsec++) { + PhotCode *code = GetPhotcodebyNsec (Nsec); + photcodeVals[Nsec] = code->code; + } + + // myLensobj[] is an intermediate accumulation structure + ALLOCATE_PTR (myLensobj, Lensobj, Nsecfilt); + + // myLensctr holds counters for each of the radial apertures / types + Lensctr *myLensctr = dvo_lensctr_init (Nsecfilt); + + // I think we have already allocated lensobj + int NLENSOBJ = 10000; + REALLOCATE (catalog->lensobj, Lensobj, NLENSOBJ); + + // Nlensobj tracks how many lensobj entries have I actually generated so far for this + // catalog. An object may have 0, 1, or more lensobj entries (max = Nsecfilt). In + // UNIONS DR3, I have lensobj values only for the i-band, but Nsecfilt = 6 + + int Nlensobj = 0; + + // In a sorted database, the lensing and measure values for a single object are all in the range: + // average->measureOffset : average->measureOffset + average->Nmeasure + // average->lensingOffset : average->lensingOffset + average->Nlensing + // NOTE: all lensing entries must have a matching measurement entry (same detID, same imageID) + // but not all measurement entries must have a lensing entry + + // For each object: + // 1) I need to examine the lensing entries and match them to their corresponding measure entries + // 2) determine the set of unique photcodes in the lensing set (XXX in future, overload photcode in Lensing). + + for (off_t i = 0; i < catalog->Naverage; i++) { + + Average *average = &catalog->average[i]; + if (average->Nlensing == 0) continue; + + // start of the measure and lensing sequences + off_t Loff = average->lensingOffset; + off_t Moff = average->measureOffset; + + // reset the myLensobj accumulators + for (int Nsec = 0; Nsec < Nsecfilt; Nsec ++) { + dvo_lensobj_init (&myLensobj[Nsec], TRUE); // init accumulated values to 0 + myLensobj[Nsec].photcode = photcodeVals[Nsec]; // set the photcodes for the accumulators + } + dvo_lensctr_reset (myLensctr, Nsecfilt); // init counters to 0 + + // assign the lensing values to the appropriate lensobj + for (int Lj = 0; Lj < average->Nlensing; Lj ++) { + Lensing *lensing = &catalog->lensing[Loff + Lj]; + if (lensing->detID < 0) continue; // XXX some invalid lensing entries? + + // need to find the corresponding Measure to get the photcode + int foundMeasure = FALSE; + for (int Mj = 0; !foundMeasure && (Mj < average->Nmeasure); Mj ++) { + Measure *measure = &catalog->measure[Moff + Mj]; + + // skip the mismatched entries + if (lensing->detID != measure->detID) continue; + if (lensing->imageID != measure->imageID) continue; + + // ** this is the matched entry + foundMeasure = TRUE; + lensing->oldImID = measure->photcode; // XXX save the photcode on the lensing structure (unused element) + + PhotCode *code = GetPhotcodebyCode (measure->photcode); + myAssert (code, "missing photcode?"); + myAssert (code->equiv > -1, "photcode not equivalent to secfilt?"); + + int Nsec = GetPhotcodeNsec (code->equiv); + myAssert (Nsec > -1, "cannot find Nsec?"); + + // relphot sets measure->Mcal (setMcalOutput.c, called by setMrelFinal.c) + // XXX : I'm using McalAPER since these lens measurements are aperture-like, right? + float Mcal = code[0].K*(measure->airmass - 1.000) + SCALE*code->C - measure->McalAPER; + float Fcal = 3630.8 * MagToFlux(Mcal); + + // lensing->F_ApR5, etc are in units of DN/sec + // Fcal * lensing->F_ApR5 is in Jy + + // F_ApR5 is , the mean of the lensing entries. we will save stdev of F_ApR5_i, + // not the r.m.s. of sF_ApR5 + + // stdev = sqrt(F_Ap_R5_i^2 / N - ) + + // accumulate valid values + dvo_lensing_accum (&myLensobj[Nsec], &myLensctr[Nsec], lensing, Fcal); + + // XXX TEST: + if (average->objID == 1702) { + fprintf (stderr, "%e : %e : %e\n", Fcal, lensing->F_ApR5, myLensobj[Nsec].F_ApR5); + } + } + myAssert (foundMeasure, "oops, unmatched lensing entry"); + } + + average->Nlensobj = 0; + + // now loop over the Nsec values and calculate averages for each + for (int Nsec = 0; Nsec < Nsecfilt; Nsec ++) { + if (!dvo_lensctr_has_values (&myLensctr[Nsec])) continue; + + // calculate averages + dvo_lensobj_aves (&myLensobj[Nsec], &myLensctr[Nsec]); + + // copy into the output array + catalog->lensobj[Nlensobj] = myLensobj[Nsec]; + catalog->lensobj[Nlensobj].objID = average->objID; + catalog->lensobj[Nlensobj].catID = average->catID; + + if (!average->Nlensobj) { + average->lensobjOffset = Nlensobj; + } + average->Nlensobj ++; + Nlensobj ++; + + CHECK_REALLOCATE (catalog->lensobj, Lensobj, NLENSOBJ, Nlensobj, 1000); + } + } + + catalog->Nlensobj = Nlensobj; + catalog->Nlensobj_disk = 0; + catalog->Nlensobj_off = 0; + REALLOCATE (catalog->lensobj, Lensobj, Nlensobj); + + free (myLensobj); + free (myLensctr); + free (photcodeVals); + + return TRUE; +} +