Index: /branches/eam_branches/relphot.20200519/include/relphot.h =================================================================== --- /branches/eam_branches/relphot.20200519/include/relphot.h (revision 41386) +++ /branches/eam_branches/relphot.20200519/include/relphot.h (revision 41387) @@ -96,4 +96,5 @@ float McalPSF; float McalAPER; + float dKlam; float dMcal; float dMsys; @@ -110,6 +111,6 @@ off_t *measure; - off_t Ncatalog; - off_t NCATALOG; + off_t Ncatalog; // XX not needed + off_t NCATALOG; // XX not needed off_t *catalog; } TGroup; @@ -332,4 +333,5 @@ int MOSAIC_ZEROPT; int TGROUP_ZEROPT; +int TGROUP_FIT_AIRMASS; int TEST_IMAGE1; @@ -457,5 +459,5 @@ void freeImageBins PROTO((int Ncatalog, int doImageList)); void freeMosaicBins PROTO((int Ncatalog, int doMosaicList)); -void freeTGroupBins PROTO((int Ncatalog, int doTGroupList)); +void freeTGroupBins PROTO((int Ncatalog)); void free_catalogs PROTO((Catalog *catalog, int Ncatalog)); @@ -468,5 +470,5 @@ float getMgrid PROTO((off_t meas, int cat)); float getMmos PROTO((off_t meas, int cat)); -float getMgrp PROTO((off_t meas, int cat)); +float getMgrp PROTO((off_t meas, int cat, float airmass)); float getMrel PROTO((Catalog *catalog, off_t meas, int cat, dvoMagClassType class, dvoMagSourceType source)); short getUbercalDist PROTO((off_t meas, int cat)); @@ -485,4 +487,5 @@ void initMosaics PROTO((Image *subset, off_t Nsubset, Image *image, char *inSubset, off_t Nimage)); void initTGroups PROTO((Image *subset, off_t *LineNumber, off_t Nsubset)); +void freeTGroups PROTO((void)); void initTGroupBins PROTO((Catalog *catalog, int Ncatalog)); void initMrel PROTO((Catalog *catalog, int Ncatalog)); @@ -493,4 +496,5 @@ int liststats_init PROTO((StatType *stats)); int liststats_irls PROTO((StatDataSet *dataset, int Npoints, StatType *stats)); +int liststats_fit1d PROTO((double *value, double *err, double *x, int Npts, StatType *stats, double *dk)); unsigned int *ReadTGroupFile PROTO((FILE *f, int *nelem)); @@ -721,2 +725,4 @@ void RepairWarpMeasuresOpenLogfile (); void RepairWarpMeasuresCloseLogfile (); + +int save_images_backup (FITS_DB *db); Index: /branches/eam_branches/relphot.20200519/src/GridOps.c =================================================================== --- /branches/eam_branches/relphot.20200519/src/GridOps.c (revision 41386) +++ /branches/eam_branches/relphot.20200519/src/GridOps.c (revision 41387) @@ -289,7 +289,7 @@ int **gotstar, **gridmeas; - int i, j, k, Ngood, Nbad, Nmos, Ncal, Nrel, Nsys; + int i, j, k, Ngood, Nbad, Nmos, Ncal, Nrel, Nsys, Ngrp; double **A, **B, *Mjx, *Wjx; - float Msys, Mcal, Mmos, Merr, Wsys; + float Msys, Mcal, Mmos, Mgrp, Merr, Wsys; double Mj, Wj; @@ -305,5 +305,5 @@ } - Ngood = Nbad = Ncal = Nmos = Nrel = Nsys = 0; + Ngood = Nbad = Ncal = Nmos = Nrel = Nsys = Ngrp = 0; ALLOCATE (gotstar, int *, Ncatalog); @@ -411,4 +411,11 @@ if (isnan(Mmos)) { Nmos ++; + continue; + } + + // skip mosaics marked as BAD + Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass); + if (isnan(Mgrp)) { + Ngrp ++; continue; } @@ -526,7 +533,7 @@ void setMgrid (Catalog *catalog, FlatCorrectionTable *flatcorr) { - int i, j, m, c, n, N, Nmax, Nbad, Nmos, Ncal, Nrel, Nsys, Nfit; + int i, j, m, c, n, N, Nmax, Nbad, Nmos, Ngrp, Ncal, Nrel, Nsys, Nfit; double *list, *dlist; - float Msys, Mrel, Mcal, Mmos; + float Msys, Mrel, Mcal, Mmos, Mgrp; StatType stats; @@ -544,5 +551,5 @@ ALLOCATE (dlist, double, Nmax); - Nbad = Ncal = Nmos = Nrel = Nsys = Nfit = 0; + Nbad = Ncal = Nmos = Ngrp = Nrel = Nsys = Nfit = 0; for (i = 0; i < Ngrid; i++) { @@ -568,4 +575,9 @@ continue; } + Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass); + if (isnan(Mgrp)) { + Ngrp ++; + continue; + } Mrel = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); if (isnan(Mrel)) { @@ -580,5 +592,5 @@ continue; } - list[N] = Msys - Mrel - Mcal - Mmos; + list[N] = Msys - Mrel - Mcal - Mmos - Mgrp; dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); N++; @@ -616,5 +628,5 @@ int i, j, m, c, n, N, Narea; - float Msys, Mrel, Mcal, Mmos; + float Msys, Mrel, Mcal, Mmos, Mgrp; double *xlist, *Mlist, *dlist, *ylist; Graphdata graphdata; @@ -649,4 +661,6 @@ Mmos = getMmos (m, c); if (isnan(Mmos)) continue; + Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass); + if (isnan(Mgrp)) continue; Mrel = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); if (isnan(Mrel)) continue; @@ -657,6 +671,6 @@ xlist[N] = Xmeas[c][m]; ylist[N] = Ymeas[c][m]; - Mlist[N] = Msys - Mrel - Mcal - Mmos; - dlist[N] = Msys - Mrel - Mcal - Mmos - gridM[i]; + Mlist[N] = Msys - Mrel - Mcal - Mmos - Mgrp; + dlist[N] = Msys - Mrel - Mcal - Mmos - Mgrp - gridM[i]; N++; } Index: /branches/eam_branches/relphot.20200519/src/ImageOps.c =================================================================== --- /branches/eam_branches/relphot.20200519/src/ImageOps.c (revision 41386) +++ /branches/eam_branches/relphot.20200519/src/ImageOps.c (revision 41387) @@ -586,5 +586,5 @@ off_t i, j, m, c, n; - int mark, bad, Nfew, Nbad, Nmos, Nrel, Ngrid, Nsys; + int mark, bad, Nfew, Nbad, Nmos, Ngrp, Nrel, Ngrid, Nsys; StatType stats; @@ -612,5 +612,5 @@ StatDataSet *brightStars = StatDataSetAlloc (1, Nmax); - Nfew = Nbad = Nmos = Ngrid = Nrel = Nsys = 0; + Nfew = Nbad = Nmos = Ngrp = Ngrid = Nrel = Nsys = 0; int Ncalibrated = 0; @@ -631,5 +631,7 @@ // only freeze the gpc1 chips - if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue; + // XXX : when do we FREEZE_IMAGES in a mix of GPC1, GPC2, etc? + // if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue; + if (FREEZE_IMAGES) continue; if (image[i].photcode == 0) continue; // skip the PHU images @@ -666,4 +668,9 @@ continue; } + float Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass); + if (isnan(Mgrp)) { + Ngrp ++; + continue; + } float Mgrid = getMgrid (m, c); if (isnan(Mgrid)) { @@ -706,5 +713,5 @@ skip: - psfStars->flxlist[Nref] = MsysPSF - MrelPSF - Mmos - Mgrid + Mflat; + psfStars->flxlist[Nref] = MsysPSF - MrelPSF - Mmos - Mgrp - Mgrid + Mflat; psfStars->errlist[Nref] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); psfStars->wgtlist[Nref] = 1; @@ -718,5 +725,5 @@ if (isfinite(MrelKron) && isfinite(MsysKron)) { - kronStars->flxlist[Nkron] = MsysKron - MrelKron - Mmos - Mgrid + Mflat; + kronStars->flxlist[Nkron] = MsysKron - MrelKron - Mmos - Mgrp - Mgrid + Mflat; kronStars->errlist[Nkron] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); kronStars->wgtlist[Nkron] = 1; @@ -1023,5 +1030,5 @@ off_t i, j, m, c, n, N; double *list, *dlist; - float Mcal, Mmos, Mgrid; + float Mcal, Mmos, Mgrp, Mgrid; StatType stats; @@ -1050,4 +1057,6 @@ Mmos = getMmos (m, c); if (isnan(Mmos)) continue; + Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass); + if (isnan(Mgrp)) continue; Mgrid = getMgrid (m, c); if (isnan(Mgrid)) continue; Index: /branches/eam_branches/relphot.20200519/src/MosaicOps.c =================================================================== --- /branches/eam_branches/relphot.20200519/src/MosaicOps.c (revision 41386) +++ /branches/eam_branches/relphot.20200519/src/MosaicOps.c (revision 41387) @@ -919,4 +919,5 @@ int Nbad; int Ncal; + int Ngrp; int Ngrid; int Nrel; @@ -1154,4 +1155,5 @@ if (fout) { float Mcal = getMcal (m, c, MAG_CLASS_PSF); + // float Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass); float Mgrid = getMgrid (m, c); float MrelPSF = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); @@ -1175,4 +1177,9 @@ if (isnan(Mcal)) { info->Ncal++; + continue; + } + float Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass); + if (isnan(Mgrp)) { + info->Ngrp++; continue; } @@ -1692,5 +1699,5 @@ off_t i, j, m, c, n, N; double *list, *dlist; - float Mcal, Mgrid, Mrel; + float Mcal, Mgrp, Mgrid, Mrel; StatType stats; @@ -1716,4 +1723,6 @@ Mcal = getMcal (m, c, MAG_CLASS_PSF); if (isnan(Mcal)) continue; + Mgrp = getMgrp (m, c, catalog[c].measureT[m].airmass); + if (isnan(Mgrp)) continue; Mgrid = getMgrid (m, c); if (isnan(Mgrid)) continue; Index: /branches/eam_branches/relphot.20200519/src/StarOps.c =================================================================== --- /branches/eam_branches/relphot.20200519/src/StarOps.c (revision 41386) +++ /branches/eam_branches/relphot.20200519/src/StarOps.c (revision 41387) @@ -491,5 +491,5 @@ float Mmos = getMmos (m, i); if (isnan(Mmos)) continue; - float Mgrp = getMgrp (m, i); + float Mgrp = getMgrp (m, i, catalog[i].measureT[m].airmass); if (isnan(Mgrp)) continue; float Mgrid = getMgrid (m, i); @@ -507,5 +507,6 @@ // set the output calibration - float Moff = Mmos + Mgrp + Mgrid; + // float Moff = Mmos + Mgrp + Mgrid; + float Moff = Mmos + Mgrid; // Mgrp has already been applied to the image.Mcal catalog[i].measure[m].McalPSF = McalPSF + Moff; catalog[i].measure[m].McalAPER = useStackAper ? McalAPER + Moff : McalPSF + Moff; @@ -730,5 +731,5 @@ Mmos = getMmos (m, i); if (isnan(Mmos)) { Nmos ++; continue; } - Mgrp = getMgrp (m, i); + Mgrp = getMgrp (m, i, catalog[i].measureT[m].airmass); if (isnan(Mgrp)) { Ngrp ++; continue; } Mgrid = getMgrid (m, i); @@ -736,5 +737,5 @@ Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt], MAG_CLASS_PSF); - list[N] = Msys - Mcal - Mmos - Mgrid; + list[N] = Msys - Mcal - Mmos - Mgrp - Mgrid; dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR); N++; Index: /branches/eam_branches/relphot.20200519/src/TGroupOps.c =================================================================== --- /branches/eam_branches/relphot.20200519/src/TGroupOps.c (revision 41386) +++ /branches/eam_branches/relphot.20200519/src/TGroupOps.c (revision 41387) @@ -36,4 +36,5 @@ // load the photometric nights from a file : these will be used to define the tgroups // XXX call this in args.c where TGROUP_ZEROPT is set +// XXX cannot call this in args.c : the Nphotcodes have not yet been loaded. void loadTGroups (char *filename) { @@ -73,4 +74,5 @@ tgroup[Nt].McalPSF = 0.0; // note : at the end, tgroup.Mcal is added back to the input images tgroup[Nt].McalAPER = 0.0; // note : tgroup stores only offsets relative to the original image values + tgroup[Nt].dKlam = 0.0; tgroup[Nt].dMcal = 0.0; // note : at the end, tgroup.Mcal is added back to the input images tgroup[Nt].dMsys = 0.0; @@ -111,6 +113,11 @@ // ignore non-GPC1 images - if (!isGPC1chip(subset[i].photcode)) { - Nsimple ++; + // if (!isGPC1chip(subset[i].photcode)) { + // Nsimple ++; + // continue; + // } + + // ignore mosaic images (photcode == 0) + if (!subset[i].photcode) { continue; } @@ -138,4 +145,20 @@ fprintf (stderr, "matched %d images to %d tgroups, %d simple chips not matched to tgroups\n", (int) (Nsubset - Nsimple), (int) Ntgroup, (int) Nsimple); + return; +} + +void freeTGroups (void) { + + if (!TGROUP_ZEROPT) return; + + free (ImageToTGroup); + free (tgroupTimes); + + for (int i = 0; i < Ntgroup; i++) { + free (tgroup[i].image); + free (tgroup[i].measure); + free (tgroup[i].catalog); + } + free (tgroup); return; } @@ -221,5 +244,5 @@ // if start is not in this range, return -1 if (start < tgroupTimes[Nlo]) return (-1); - if (start > tgroupTimes[Nhi]) return (-1); + if (start > tgroupTimes[Nhi] + 86399) return (-1); // 1-day range for tgroup for now while (Nhi - Nlo > 4) { @@ -259,6 +282,7 @@ for (j = 0; j < tgroup[i].Nimage; j++) { im = tgroup[i].image[j]; - image[im].McalPSF += tgroup[i].McalPSF; - image[im].McalAPER += tgroup[i].McalAPER; + double Mgrp = (TGROUP_FIT_AIRMASS) ? tgroup[i].McalPSF + tgroup[i].dKlam*(image[im].secz - 1.0) : tgroup[i].McalPSF; + image[im].McalPSF += Mgrp; + image[im].McalAPER += Mgrp; image[im].dMcal = tgroup[i].dMcal; image[im].McalChiSq = tgroup[i].McalChiSq; @@ -267,7 +291,10 @@ image[im].flags |= (tgroup[i].flags & ID_IMAGE_PHOTOM_FEW); image[im].flags |= (tgroup[i].flags & ID_IMAGE_PHOTOM_POOR); + + fprintf (stderr, "TG to IMAGE: %f -> %f (%d)\n", tgroup[i].McalPSF, image[im].McalPSF, tgroup[i].nFitPhotom); } tgroup[i].McalPSF = 0.0; tgroup[i].McalAPER = 0.0; + tgroup[i].dKlam = 0.0; } } @@ -287,5 +314,5 @@ } -void freeTGroupBins (int Ncatalog, int doTGroupList) { +void freeTGroupBins (int Ncatalog) { off_t i; @@ -329,5 +356,8 @@ } } - // fprintf (stderr, "matched %d detections to tgroups\n", Nmatch); + fprintf (stderr, "Matched %d detections to tgroups\n", Nmatch); + for (i = 0; i < Ntgroup; i++) { + fprintf (stderr, " TGROUP time %d, Nmeasure: %d, Nimage: %d\n", tgroup[i].start, (int) tgroup[i].Nmeasure, (int) tgroup[i].Nimage); + } return (TRUE); } @@ -380,5 +410,5 @@ } -float getMgrp (off_t meas, int cat) { +float getMgrp (off_t meas, int cat, float airmass) { off_t i; @@ -392,5 +422,5 @@ // if (tgroup[i].flags & IMAGE_BAD) return (NAN); - value = tgroup[i].McalPSF; + value = tgroup[i].McalPSF + tgroup[i].dKlam*(airmass - 1.0); return (value); } @@ -419,4 +449,5 @@ int PoorImages; double *psfMagList; + double *psfAirList; double *psfErrList; double *brightMagList; @@ -456,4 +487,5 @@ if (allocLists) { ALLOCATE (info->psfMagList, double, Nmax); + ALLOCATE (info->psfAirList, double, Nmax); ALLOCATE (info->psfErrList, double, Nmax); ALLOCATE (info->kronMagList, double, Nmax); @@ -466,4 +498,5 @@ void SetMgrpInfoFree (SetMgrpInfo *info) { free (info->psfMagList ); + free (info->psfAirList ); free (info->psfErrList ); free (info->kronMagList ); @@ -522,6 +555,4 @@ image = getimages (&N, NULL); - fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE); - if (PoorImages) { // XXX use bad stars and measurements for PoorImages? or not? @@ -543,4 +574,5 @@ for (i = 0; i < Ntgroup; i++) { setMgrp_tgroup (&tgroup[i], i, image, catalog, &info, flatcorr); + fprintf (stderr, "TGROUP Mcal: %f, dK: %f, limiting negative clouds to %f\n", tgroup[i].McalPSF, tgroup[i].dKlam, CLOUD_TOLERANCE); } SetMgrpInfoFree (&info); @@ -575,4 +607,5 @@ double *psfMagList = info->psfMagList; + double *psfAirList = info->psfAirList; double *psfErrList = info->psfErrList; double *kronMagList = info->kronMagList; @@ -593,5 +626,5 @@ int testImage = FALSE; - // testImage |= (abs(myTGroup[0].start - 1324104046) < 10); + testImage |= (abs(myTGroup[0].start - 1245283200) < 10); // testImage |= (abs(myTGroup[0].start - 1324103823) < 10); @@ -633,5 +666,5 @@ int isBad = (catalog[c].measureT[m].dbFlags & MEAS_BAD); - fprintf (fout, "%f %f : %f %f %f %f %f : %f %d\n", catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysPSF, MrelPSF, Mcal, Mgrid, Mflat, delta, isBad); + fprintf (fout, "%f %f : %f %f %f %f %f : %f %f %d\n", catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysPSF, MrelPSF, Mcal, Mgrid, Mflat, catalog[c].measureT[m].airmass, delta, isBad); } @@ -697,4 +730,5 @@ psfMagList[N] = MsysPSF - MrelPSF - Mcal - Mgrid + Mflat; + psfAirList[N] = (catalog[c].measureT[m].airmass - 1.0); psfErrList[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); kronMagList[N] = MsysKron - MrelKron - Mcal - Mgrid + Mflat; @@ -725,5 +759,10 @@ } - liststats (psfMagList, psfErrList, NULL, N, &stats); + double dk = 0.0; + if (TGROUP_FIT_AIRMASS) { + liststats_fit1d (psfMagList, psfErrList, psfAirList, N, &stats, &dk); + } else { + liststats (psfMagList, psfErrList, NULL, N, &stats); + } if (VERBOSE2 && info->PoorImages) fprintf (stderr, "Mgrp: %f %f %d %d\n", stats.mean, stats.sigma, stats.Nmeas, N); @@ -733,4 +772,5 @@ myTGroup[0].McalChiSq = stats.chisq; myTGroup[0].nFitPhotom = N; + myTGroup[0].dKlam = dk; // drop this for tgroups? Index: /branches/eam_branches/relphot.20200519/src/args.c =================================================================== --- /branches/eam_branches/relphot.20200519/src/args.c (revision 41386) +++ /branches/eam_branches/relphot.20200519/src/args.c (revision 41387) @@ -317,32 +317,54 @@ } + FREEZE_IMAGES = FALSE; + if ((N = get_argument (argc, argv, "-imfreeze"))) { + remove_argument (N, &argc, argv); + FREEZE_IMAGES = TRUE; + } + FREEZE_MOSAICS = FALSE; + if ((N = get_argument (argc, argv, "-mosfreeze"))) { + remove_argument (N, &argc, argv); + FREEZE_MOSAICS = TRUE; + } + CALIBRATE_STACKS_AND_WARPS = FALSE; + if ((N = get_argument (argc, argv, "-only-stacks-and-warps"))) { + remove_argument (N, &argc, argv); + CALIBRATE_STACKS_AND_WARPS = TRUE; + } + USE_MCAL_PSF_FOR_STACK_APER = FALSE; + if ((N = get_argument (argc, argv, "-use-mcal-psf-for-stack-aper"))) { + remove_argument (N, &argc, argv); + USE_MCAL_PSF_FOR_STACK_APER = TRUE; + } + // calibrate the photometry by night / time period TGROUP_ZEROPT = FALSE; + char *TGROUP_FILENAME = NULL; if ((N = get_argument (argc, argv, "-tgroup"))) { TGROUP_ZEROPT = TRUE; remove_argument (N, &argc, argv); - loadTGroups (argv[N]); - remove_argument (N, &argc, argv); - } - - FREEZE_IMAGES = FALSE; - if ((N = get_argument (argc, argv, "-imfreeze"))) { + TGROUP_FILENAME = strcreate(argv[N]); remove_argument (N, &argc, argv); FREEZE_IMAGES = TRUE; } - FREEZE_MOSAICS = FALSE; - if ((N = get_argument (argc, argv, "-mosfreeze"))) { - remove_argument (N, &argc, argv); - FREEZE_MOSAICS = TRUE; - } - CALIBRATE_STACKS_AND_WARPS = FALSE; - if ((N = get_argument (argc, argv, "-only-stacks-and-warps"))) { - remove_argument (N, &argc, argv); - CALIBRATE_STACKS_AND_WARPS = TRUE; - } - USE_MCAL_PSF_FOR_STACK_APER = FALSE; - if ((N = get_argument (argc, argv, "-use-mcal-psf-for-stack-aper"))) { - remove_argument (N, &argc, argv); - USE_MCAL_PSF_FOR_STACK_APER = TRUE; + if ((N = get_argument (argc, argv, "-tgroups"))) { + if (TGROUP_ZEROPT) { + fprintf (stderr, "use only one of -tgroups and -tgroup (same meaning)\n"); + exit (2); + } + TGROUP_ZEROPT = TRUE; + remove_argument (N, &argc, argv); + TGROUP_FILENAME = strcreate(argv[N]); + remove_argument (N, &argc, argv); + FREEZE_IMAGES = TRUE; + } + TGROUP_FIT_AIRMASS = FALSE; + if ((N = get_argument (argc, argv, "-tgroup-fit-airmass"))) { + if (!TGROUP_ZEROPT) { + fprintf (stderr, "-tgroup-fit-airmass requires -tgroup (filename) option\n"); + exit (2); + } + TGROUP_FIT_AIRMASS = TRUE; + remove_argument (N, &argc, argv); } @@ -528,4 +550,9 @@ if (argc != 1) relphot_usage (); + if (TGROUP_ZEROPT) { + loadTGroups (TGROUP_FILENAME); + free (TGROUP_FILENAME); + } + return mode; } @@ -551,4 +578,6 @@ FreePhotcodeTable(); + freeTGroups(); + ohana_memcheck (VERBOSE); ohana_memdump (VERBOSE); Index: /branches/eam_branches/relphot.20200519/src/liststats.c =================================================================== --- /branches/eam_branches/relphot.20200519/src/liststats.c (revision 41386) +++ /branches/eam_branches/relphot.20200519/src/liststats.c (revision 41387) @@ -172,9 +172,11 @@ for (i = ks; i < ke; i++) { M = SQ (value[i] - Mo); - dM = SQ (dvalue[i]); - X2 += M / dM; + if (dvalue) { + dM = SQ (dvalue[i]); + X2 += M / dM; + } dS += M; } - X2 = X2 / (Nm - 1); + X2 = dvalue ? X2 / (Nm - 1) : NAN; dS = sqrt (dS / (Nm - 1)); @@ -187,4 +189,31 @@ return (TRUE); } + +int liststats_fit1d (double *value, double *err, double *x, int Npts, StatType *stats, double *dk) { + + double M = 0.0, MZ = 0.0, Z1 = 0.0, Z2 = 0.0, R = 0.0; + + for (int i = 0; i < Npts; i++) { + double wt = 1.0 / (err[i]*err[i]); + double V = value[i]; + double Z = x[i]; + M += V * wt; + MZ += V*Z * wt; + Z1 += Z * wt; + Z2 += Z*Z * wt; + R += 1.0 * wt; + } + + double DetInv = (R*Z2 - Z1*Z1); + double Det = (fabs(DetInv) < -1e6) ? 0.0 : 1.0 / DetInv; // do not allow unstable solutions + + double zp = (M*Z2 - Z1*MZ)*Det; + *dk = (MZ*R - M*Z1)*Det; + + stats->mean = zp; + stats->chisq = NAN; + return TRUE; +} + // These should probably be tunable: Index: /branches/eam_branches/relphot.20200519/src/reload_catalogs.c =================================================================== --- /branches/eam_branches/relphot.20200519/src/reload_catalogs.c (revision 41386) +++ /branches/eam_branches/relphot.20200519/src/reload_catalogs.c (revision 41387) @@ -86,4 +86,6 @@ initImageBins (&catalog, 1, FALSE); initMosaicBins (&catalog, 1, FALSE); + initTGroupBins (&catalog, 1); + initGridBins (&catalog, 1); TIMESTAMP(time3b); @@ -91,4 +93,6 @@ findImages (&catalog, 1, FALSE); findMosaics (&catalog, 1, FALSE); + findTGroups (&catalog, 1); + TIMESTAMP(time4); @@ -117,6 +121,9 @@ TIMESTAMP(time6); + // XXX freeTgroupBins here: + freeImageBins (1, FALSE); freeMosaicBins (1, FALSE); + freeTGroupBins (1); freeGridBins (1); TIMESTAMP(time7); Index: /branches/eam_branches/relphot.20200519/src/relphot_images.c =================================================================== --- /branches/eam_branches/relphot.20200519/src/relphot_images.c (revision 41386) +++ /branches/eam_branches/relphot.20200519/src/relphot_images.c (revision 41387) @@ -81,5 +81,5 @@ MARKTIME("-- set up mosaic indexes: %f sec\n", dtime); - findTGroups (catalog, Ncatalog); /* also sets Grid values */ + findTGroups (catalog, Ncatalog); MARKTIME("-- set up mosaic indexes: %f sec\n", dtime); @@ -158,4 +158,5 @@ if ((i > 8) && (i % 8 == 3)) clean_stars (catalog, Ncatalog); if ((i > 8) && (i % 8 == 5)) clean_mosaics (); + if ((i > 8) && (i % 8 == 5)) clean_tgroups (); if ((i > 8) && (i % 8 == 5)) clean_images (); @@ -258,4 +259,5 @@ freeImageBins (Ncatalog, TRUE); freeMosaicBins (Ncatalog, TRUE); + freeTGroupBins (Ncatalog); freeGridBins (Ncatalog); Index: /branches/eam_branches/relphot.20200519/src/setMrelCatalog.c =================================================================== --- /branches/eam_branches/relphot.20200519/src/setMrelCatalog.c (revision 41386) +++ /branches/eam_branches/relphot.20200519/src/setMrelCatalog.c (revision 41387) @@ -221,5 +221,5 @@ // ** Choose the calibration (depends on the mode : do I have an image reference or not?) - float Mcal = 0.0, Mmos = 0.0, Mgrid = 0.0, Mflat = 0.0; + float Mcal = 0.0, Mmos = 0.0, Mgrid = 0.0, Mflat = 0.0, Mgrp = 0.0; off_t meas = measureOffset + k; if (getImageEntry (meas, cat) < 0) { @@ -235,4 +235,6 @@ Mmos = getMmos (meas, cat); if (isnan(Mmos)) SKIP_THIS_MEAS(Nmos); + Mgrp = getMgrp (meas, cat, measureT[k].airmass); + if (isnan(Mgrp)) SKIP_THIS_MEAS(Ngrp); Mgrid = getMgrid (meas, cat); if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid); @@ -261,5 +263,5 @@ if (CHECK_VALID_MAG(Map, dMap)) { int Nap = results->aperData[Nsec].Nlist; - results->aperData[Nsec].flxlist[Nap] = Map - Mflat - Mcal - Mmos - Mgrid; // this is consistent with PhotRel + results->aperData[Nsec].flxlist[Nap] = Map - Mflat - Mcal - Mmos - Mgrp - Mgrid; // this is consistent with PhotRel results->aperData[Nsec].errlist[Nap] = dMap; results->aperData[Nsec].wgtlist[Nap] = useUbercalWeight ? UBERCAL_WEIGHT : 1.0; @@ -281,5 +283,5 @@ if (CHECK_VALID_MAG(Mkron, dMkron)) { int Nkron = results->kronData[Nsec].Nlist; - results->kronData[Nsec].flxlist[Nkron] = Mkron - Mflat - Mcal - Mmos - Mgrid; // this is consistent with PhotRel + results->kronData[Nsec].flxlist[Nkron] = Mkron - Mflat - Mcal - Mmos - Mgrp - Mgrid; // this is consistent with PhotRel results->kronData[Nsec].errlist[Nkron] = dMkron; results->kronData[Nsec].wgtlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0; @@ -294,5 +296,5 @@ if (CHECK_VALID_MAG(Mpsf, dMpsf)) { int Npsf = results->psfData[Nsec].Nlist; - results->psfData[Nsec].flxlist[Npsf] = Mpsf - Mflat - Mcal - Mmos - Mgrid; // this is consistent with PhotRel + results->psfData[Nsec].flxlist[Npsf] = Mpsf - Mflat - Mcal - Mmos - Mgrp - Mgrid; // this is consistent with PhotRel results->psfData[Nsec].errlist[Npsf] = dMpsf; results->psfData[Nsec].wgtlist[Npsf] = isUbercal ? UBERCAL_WEIGHT : 1.0; @@ -711,4 +713,5 @@ Mmos = getMmos (measSeq, cat); Mgrid = getMgrid (measSeq, cat); + // XXX can Mmos and Mgrid exist for stacks? }