Index: trunk/Ohana/src/relphot/Makefile =================================================================== --- trunk/Ohana/src/relphot/Makefile (revision 41388) +++ trunk/Ohana/src/relphot/Makefile (revision 41390) @@ -32,4 +32,6 @@ $(SRC)/BrightCatalog.$(ARCH).o \ $(SRC)/MosaicOps.$(ARCH).o \ +$(SRC)/TGroupOps.$(ARCH).o \ +$(SRC)/ReadTGroupFile.$(ARCH).o \ $(SRC)/SetSignals.$(ARCH).o \ $(SRC)/Shutdown.$(ARCH).o \ @@ -89,4 +91,6 @@ $(SRC)/BrightCatalog.$(ARCH).o \ $(SRC)/MosaicOps.$(ARCH).o \ +$(SRC)/TGroupOps.$(ARCH).o \ +$(SRC)/ReadTGroupFile.$(ARCH).o \ $(SRC)/SetSignals.$(ARCH).o \ $(SRC)/Shutdown.$(ARCH).o \ Index: trunk/Ohana/src/relphot/include/relphot.h =================================================================== --- trunk/Ohana/src/relphot/include/relphot.h (revision 41388) +++ trunk/Ohana/src/relphot/include/relphot.h (revision 41390) @@ -90,4 +90,30 @@ } Mosaic; +typedef struct { + unsigned int start; + unsigned int stop; + short photcode; + float McalPSF; + float McalAPER; + float dKlam; + float dMcal; + float dMsys; + float McalChiSq; + unsigned short nFitPhotom; + unsigned int flags; + + off_t Nimage; + off_t NIMAGE; + off_t *image; + + off_t Nmeasure; + off_t NMEASURE; + off_t *measure; + + off_t Ncatalog; // XX not needed + off_t NCATALOG; // XX not needed + off_t *catalog; +} TGroup; + typedef enum { STATS_NONE, @@ -145,4 +171,5 @@ int Ncal; int Nmos; + int Ngrp; int Ngrid; @@ -305,4 +332,6 @@ int VERBOSE_IMAGE; int MOSAIC_ZEROPT; +int TGROUP_ZEROPT; +int TGROUP_FIT_AIRMASS; int TEST_IMAGE1; @@ -403,4 +432,5 @@ void clean_measures PROTO((Catalog *catalog, int Ncatalog, int final, FlatCorrectionTable *flatcorr)); void clean_mosaics PROTO((void)); +void clean_tgroups PROTO((void)); void clean_stars PROTO((Catalog *catalog, int Ncatalog)); int corner_check PROTO((double *x1, double *y1, double *x2, double *y2)); @@ -410,4 +440,5 @@ void findImages PROTO((Catalog *catalog, int Ncatalog, int doImageList)); int findMosaics PROTO((Catalog *catalog, int Ncatalog, int doMosaicList)); +int findTGroups PROTO((Catalog *catalog, int Ncatalog)); void clearImages (void); @@ -428,4 +459,6 @@ void freeImageBins PROTO((int Ncatalog, int doImageList)); void freeMosaicBins PROTO((int Ncatalog, int doMosaicList)); +void freeTGroupBins PROTO((int Ncatalog)); + void free_catalogs PROTO((Catalog *catalog, int Ncatalog)); int gcatalog PROTO((Catalog *catalog, int FINAL)); @@ -437,4 +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 airmass)); float getMrel PROTO((Catalog *catalog, off_t meas, int cat, dvoMagClassType class, dvoMagSourceType source)); short getUbercalDist PROTO((off_t meas, int cat)); @@ -452,4 +486,7 @@ void initMosaicGrid PROTO((Image *image, off_t Nimage)); 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)); RelphotMode initialize PROTO((int argc, char **argv)); @@ -459,4 +496,11 @@ 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)); +void loadTGroups PROTO((char *filename)); +void initTGroupsMcal PROTO((void)); +TGroup *getTGroupForImage PROTO((off_t im)); +off_t findTGroup PROTO((unsigned int start, int photcode)); Catalog *load_catalogs PROTO((SkyList *skylist, int *Ncatalog, int hostID, char *hostpath, char *syncfile)); @@ -470,4 +514,6 @@ void matchImage PROTO((Catalog *catalog, off_t meas, int cat, int doImageList)); void matchMosaics PROTO((Catalog *catalog, off_t meas, int cat, int doMosaicList)); +void matchTGroups PROTO((Catalog *catalog, off_t meas, int cat)); + double opening_angle PROTO((double x1, double y1, double x2, double y2, double x3, double y3)); void plot_chisq PROTO((Catalog *catalog, int Ncatalog)); @@ -491,8 +537,10 @@ int setExclusions PROTO((Catalog *catalog, int Ncatalog, int verbose)); void setMcal PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr)); -void setMcalFinal PROTO((void)); +void setMcalFromMosaics PROTO((void)); +void setMcalFromTGroups PROTO((void)); int setMcalOutput PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr)); void setMgrid PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr)); int setMmos PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr)); +int setMgrp PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr)); int setMrel PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr)); void setMrelFinal PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr, int simpleAverage)); @@ -677,2 +725,4 @@ void RepairWarpMeasuresOpenLogfile (); void RepairWarpMeasuresCloseLogfile (); + +int save_images_backup (FITS_DB *db); Index: trunk/Ohana/src/relphot/src/GridOps.c =================================================================== --- trunk/Ohana/src/relphot/src/GridOps.c (revision 41388) +++ trunk/Ohana/src/relphot/src/GridOps.c (revision 41390) @@ -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: trunk/Ohana/src/relphot/src/ImageOps.c =================================================================== --- trunk/Ohana/src/relphot/src/ImageOps.c (revision 41388) +++ trunk/Ohana/src/relphot/src/ImageOps.c (revision 41390) @@ -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: trunk/Ohana/src/relphot/src/MosaicOps.c =================================================================== --- trunk/Ohana/src/relphot/src/MosaicOps.c (revision 41388) +++ trunk/Ohana/src/relphot/src/MosaicOps.c (revision 41390) @@ -720,5 +720,5 @@ } -void setMcalFinal () { +void setMcalFromMosaics () { off_t i, j, im, Nimage; @@ -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: trunk/Ohana/src/relphot/src/ReadTGroupFile.c =================================================================== --- trunk/Ohana/src/relphot/src/ReadTGroupFile.c (revision 41390) +++ trunk/Ohana/src/relphot/src/ReadTGroupFile.c (revision 41390) @@ -0,0 +1,100 @@ +# include "relphot.h" + +# define NBUFFER 0x10000 +# define IsCSV FALSE +# define COLUMN_NUMBER 1 + +unsigned int *ReadTGroupFile (FILE *f, int *nelem) { + + // we allocate one extra byte into which we never read so there will always be a NULL terminating the string + ALLOCATE_PTR (buffer, char, NBUFFER + 1); + bzero (buffer, NBUFFER + 1); + + int Nstart = 0; // location of the last valid byte in the buffer (start filling here) + int EndOfFile = FALSE; + + // int TimeFormat = TIME_MJD; + // time_t TimeReference = 0; + + int Nelem = 0; + int NELEM = 1000; + ALLOCATE_PTR (elem, unsigned int, NELEM); + + while (!EndOfFile) { + + int Nbytes = NBUFFER - Nstart; + bzero (&buffer[Nstart], Nbytes + 1); + // ^-- we have allocated one extra byte into which we never read so there will always be a NULL terminating the string + + int nread = fread (&buffer[Nstart], 1, Nbytes, f); + if (ferror (f)) { + perror ("error reading data file"); + return NULL; + } + + // we still need to parse the rest of the buffer, but there might not be an EOL on the last line + if (nread == 0) { + EndOfFile = TRUE; + } + + int bufferStatus = TRUE; + + char *c0 = buffer; // c0 always marks the start of a line + while (bufferStatus) { + char *c1 = strchr (c0, '\n'); // find the end of this current line (also valid for a Mac: \r\n) + if (!c1) { + c1 = strchr (c0, '\r'); // try \r for Windows files + } + if (!c1) { + Nstart = strlen (c0); + if (EndOfFile) { + // if we have reached EOF, we need to do one last pass in case there is a line without a return + c1 = c0 + Nstart; + bufferStatus = FALSE; + if (Nstart == 0) continue; // if we have reached EOF and c0 points at the last valid character, we are done + } else { + // if we have not reached EOF, we need to shift the buffer to the start of this line and read more data + memmove (buffer, c0, Nstart); + bufferStatus = FALSE; + continue; + } + } + *c1 = 0; // mark the end of the line + + // skip to the next line (but if EOF, do not overrun buffer) + if (*c0 == '#') { if (!EndOfFile) { c0 = c1 + 1; } continue; } + if (*c0 == '!') { if (!EndOfFile) { c0 = c1 + 1; } continue; } + + // the times are written as MJD in the file: + double dvalue; + int readStatus = IsCSV ? dparse_csv (&dvalue, COLUMN_NUMBER, c0) : dparse (&dvalue, COLUMN_NUMBER, c0); + if (!readStatus) continue; // skip invalid entries + + // convert the (double) mjd value to (unsigned int) time + elem[Nelem] = ohana_mjd_to_sec (dvalue); + + Nelem ++; + if (Nelem == NELEM) { + NELEM += 1000; + REALLOCATE (elem, unsigned int, NELEM); + } + + if (!EndOfFile) { + c0 = c1 + 1; + } + } + } + FREE (buffer); + + *nelem = Nelem; + return elem; +} + +/* + time_t tvalue; + int readStatus = IsCSV ? tparse_csv (&tvalue, COLUMN_NUMBER, c0) : tparse (&tvalue, COLUMN_NUMBER, c0); + if (!readStatus) continue; // skip invalid entries + + dvalue = TimeValue (tvalue, TimeReference, TimeFormat); + elem[Nelem] = (int) dvalue; +*/ Index: trunk/Ohana/src/relphot/src/StarOps.c =================================================================== --- trunk/Ohana/src/relphot/src/StarOps.c (revision 41388) +++ trunk/Ohana/src/relphot/src/StarOps.c (revision 41390) @@ -105,4 +105,5 @@ results->Ncal = 0; results->Nmos = 0; + results->Ngrp = 0; results->Ngrid = 0; } @@ -264,4 +265,5 @@ summary->Ncal += results->Ncal ; summary->Nmos += results->Nmos ; + summary->Ngrp += results->Ngrp ; summary->Ngrid += results->Ngrid; } @@ -313,5 +315,5 @@ SetMrelInfoAccum (&summary, &results); } - if (VERBOSE2) fprintf (stderr, "%d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys); + if (VERBOSE2) fprintf (stderr, "%d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrp: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrp, summary.Ngrid, summary.Nsys); SetMrelInfoFree (&results); @@ -336,5 +338,5 @@ SetMrelInfoAccum (&summary, &results); } - if (VERBOSE2) fprintf (stderr, "%d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys); + if (VERBOSE2) fprintf (stderr, "%d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrp: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrp, summary.Ngrid, summary.Nsys); SetMrelInfoFree (&results); @@ -393,5 +395,5 @@ // report stats & summary from the threads for (i = 0; i < NTHREADS; i++) { - if (VERBOSE2) fprintf (stderr, "setMrel thread %d : %d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", + if (VERBOSE2) fprintf (stderr, "setMrel thread %d : %d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrp: %d, Ngrid: %d, Nsys: %d)\n", i, threadinfo[i].summary.Ncode, @@ -400,4 +402,5 @@ threadinfo[i].summary.Ncal, threadinfo[i].summary.Nmos, + threadinfo[i].summary.Ngrp, threadinfo[i].summary.Ngrid, threadinfo[i].summary.Nsys); @@ -405,5 +408,5 @@ SetMrelInfoFree (&threadinfo[i].summary); } - if (VERBOSE2) fprintf (stderr, "total : %d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys); + if (VERBOSE2) fprintf (stderr, "total : %d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrp: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrp, summary.Ngrid, summary.Nsys); free (threadinfo); @@ -488,4 +491,6 @@ float Mmos = getMmos (m, i); if (isnan(Mmos)) continue; + float Mgrp = getMgrp (m, i, catalog[i].measureT[m].airmass); + if (isnan(Mgrp)) continue; float Mgrid = getMgrid (m, i); if (isnan(Mgrid)) continue; @@ -502,6 +507,8 @@ // set the output calibration - catalog[i].measure[m].McalPSF = McalPSF + Mmos + Mgrid; - catalog[i].measure[m].McalAPER = useStackAper ? McalAPER + Mmos + Mgrid : McalPSF + Mmos + 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; if (catalog[i].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) { @@ -662,6 +669,6 @@ off_t *ilist; double *tlist, *list, *dlist; - float Msys, Mcal, Mmos, Mgrid; - int Ncal, Nmos, Ngrid, Nfew; + float Msys, Mcal, Mmos, Mgrp, Mgrid; + int Ncal, Nmos, Ngrp, Ngrid, Nfew; int Nsecfilt = GetPhotcodeNsecfilt (); @@ -689,5 +696,5 @@ Ndel = Nave = 0; - Ncal = Nmos = Ngrid = Nfew = 0; + Ncal = Nmos = Ngrp = Ngrid = Nfew = 0; for (i = 0; i < Ncatalog; i++) { for (j = 0; j < catalog[i].Naverage; j++) { @@ -724,9 +731,11 @@ Mmos = getMmos (m, i); if (isnan(Mmos)) { Nmos ++; continue; } + Mgrp = getMgrp (m, i, catalog[i].measureT[m].airmass); + if (isnan(Mgrp)) { Ngrp ++; continue; } Mgrid = getMgrid (m, i); if (isnan(Mgrid)) { Ngrid ++; continue; } 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: trunk/Ohana/src/relphot/src/TGroupOps.c =================================================================== --- trunk/Ohana/src/relphot/src/TGroupOps.c (revision 41390) +++ trunk/Ohana/src/relphot/src/TGroupOps.c (revision 41390) @@ -0,0 +1,1226 @@ +# include "relphot.h" + +// XXX is this needed? +void plot_setMcal (double *list, int Npts, StatType *stats, float clouds); + +// see discussion in ImagesOps.c re: IDX_T + +// array of tgroup definition structures +static off_t Ntgroup; +static TGroup *tgroup; + +// tgroup index for given image : ImageToTGroup[ImageIndex] = TGroupIndex (ImageIndex : 0 < Nimage) +static off_t *ImageToTGroup; + +// elsewhere, we have loaded a set of catalogs with measures (catalog[cat].measure[meas]) + +// each tgroup has N_onTGroup[TGroupIndex] measurements +// static off_t *N_onTGroup; // actual number of measurements on tgroup +// static off_t *N_ONTGROUP; // allocated number of measurements on tgroup + +// relationships between the measure,catalog set and the tgroups: +static off_t **MeasureToTGroup; // TGroup index from measure,catalog : MeasureToTGroup[cat][meas] = TGroupIndex +// static off_t **TGroupToCatalog; // catalog for given measure on tgroup : TGroupCatalog[TGroupIndex][i] = cat (i : 0 < NonTGroup[TGroupIndex]) +// static off_t **TGroupToMeasure; // measure for given measure on tgroup : TGroupMeasure[TGroupIndex][i] = cat (i : 0 < NonTGroup[TGroupIndex]) + +// tgroupTimes carries the times of the photometric nights +static int NtgroupTimes; +static unsigned int *tgroupTimes; + +void sort_times (unsigned int *T, int N); + +/* find tgroups (unique time periods) (NOTE : we do NOT require matching photcodes...) + we need to have a tgroup for each photcode available for a given photometric night + */ + +// 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) { + + if (!TGROUP_ZEROPT) return; + + // read start times from a file (currently a list of MJD values) + FILE *f = fopen (filename, "r"); + if (!f) { + fprintf (stderr, "tgroup file %s not found\n", filename); + exit (4); + } + tgroupTimes = ReadTGroupFile (f, &NtgroupTimes); + if (!tgroupTimes) { + fprintf (stderr, "error reading tgroup file %s\n", filename); + exit (4); + } + + // sort the times for quick bisection + sort_times (tgroupTimes, NtgroupTimes); + + // we define a tgroup for each photcode being considered + // below, if there are no images in this tgroup with a given + // photcode, it is left empty (no images) + Ntgroup = NtgroupTimes*Nphotcodes; + + ALLOCATE (tgroup, TGroup, Ntgroup); + + // the start of an MJD day is 14:00 HST, so the end time can be the start time + 23h 59m 59s (86399s) + // XXX NOTE that this is HST-centric + + for (int i = 0; i < NtgroupTimes; i++) { + for (int j = 0; j < Nphotcodes; j++) { + /* a new tgroup, define ranges */ + int Nt = i*Nphotcodes + j; + tgroup[Nt].start = tgroupTimes[i]; + tgroup[Nt].stop = tgroupTimes[i] + 86399; + 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; + tgroup[Nt].McalChiSq = 0.0; // NAN or 0.0? + tgroup[Nt].flags = 0; + tgroup[Nt].photcode = photcodes[j][0].code; + + tgroup[Nt].NIMAGE = 1000; + tgroup[Nt].Nimage = 0; + ALLOCATE (tgroup[Nt].image, off_t, tgroup[Nt].NIMAGE); + tgroup[Nt].image[0] = -1; + + tgroup[Nt].Nmeasure = 0; + tgroup[Nt].NMEASURE = 1000; + ALLOCATE (tgroup[Nt].measure, off_t, tgroup[Nt].NMEASURE); + ALLOCATE (tgroup[Nt].catalog, off_t, tgroup[Nt].NMEASURE); + } + } + return; +} + +// assign the images to the tgroups +void initTGroups (Image *subset, off_t *LineNumber, off_t Nsubset) { + + if (!TGROUP_ZEROPT) return; + + /* a 'tgroup' in relphot is a virtual concept: there is no + * entry in the image table that represents this tgroup. Instead, it is an + * internal construct that defines a group of related images + */ + + ALLOCATE (ImageToTGroup, off_t, Nsubset); // tgroup to which image belongs + + // assign each image to a tgroup + int Nsimple = 0; + for (off_t i = 0; i < Nsubset; i++) { + ImageToTGroup[i] = -1; + + // ignore non-GPC1 images + // if (!isGPC1chip(subset[i].photcode)) { + // Nsimple ++; + // continue; + // } + + // ignore mosaic images (photcode == 0) + if (!subset[i].photcode) { + continue; + } + + off_t j = findTGroup(subset[i].tzero, subset[i].photcode); + if (j == -1) { + // images which do not belong to a tgroup should be ignored in this analysis + // XXX mark the images and the measurements somehow + continue; + } + + // add reference from image to tgroup + ImageToTGroup[i] = j; + + // add image to tgroup image list + tgroup[j].image[tgroup[j].Nimage] = i; // should this be 'i' or LineNumber[i] ?? + tgroup[j].Nimage ++; + if (tgroup[j].Nimage == tgroup[j].NIMAGE) { + tgroup[j].NIMAGE += 1000; + REALLOCATE (tgroup[j].image, off_t, tgroup[j].NIMAGE); + } + } + + initTGroupsMcal (); + + 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; +} + +void initTGroupsMcal (void) { + + off_t Nimages; + off_t *LineNumber; + Image *images = getimages (&Nimages, &LineNumber); + + if (!TGROUP_ZEROPT) return; + + StatType stats; + liststats_setmode (&stats, "MEDIAN"); + + // init the tgroup array values + for (int i = 0; i < Ntgroup; i++) { + tgroup[i].McalPSF = 0.0; // note : at the end, tgroup.Mcal is added back to the input images + tgroup[i].McalAPER = 0.0; // note : tgroup stores only offsets relative to the original image values + tgroup[i].dMcal = 0.0; // note : at the end, tgroup.Mcal is added back to the input images + tgroup[i].dMsys = 0.0; + tgroup[i].McalChiSq = 0.0;// NAN or 0.0? + + if (tgroup[i].Nimage == 0) continue; // no images, ignore the tgroup + + // calculate the median of the Mcal values for the images in this tgroup: + + ALLOCATE_PTR (McalPSF, double, tgroup[i].Nimage); + ALLOCATE_PTR (McalAPER, double, tgroup[i].Nimage); + for (int j = 0; j < tgroup[i].Nimage; j++) { + int Nsub = LineNumber[tgroup[i].image[j]]; + McalPSF[j] = images[Nsub].McalPSF; + McalAPER[j] = images[Nsub].McalAPER; + + // set these back to 0.0 to have no future effect + images[Nsub].McalPSF = 0.0; + images[Nsub].McalAPER = 0.0; + } + + liststats (McalPSF, NULL, NULL, tgroup[i].Nimage, &stats); + tgroup[i].McalPSF = stats.median; + FREE (McalPSF); + + liststats (McalAPER, NULL, NULL, tgroup[i].Nimage, &stats); + tgroup[i].McalAPER = stats.median; + FREE (McalAPER); + } + return; +} + +TGroup *getTGroupForImage (off_t im) { + + if (im < 0) return NULL; + if (!ImageToTGroup) return NULL; + + off_t m = ImageToTGroup[im]; + if (m < 0) return NULL; + if (m >= Ntgroup) return NULL; + + return (&tgroup[m]); +} + +// use bisection to find the overlapping tgroup (returns exact match) +// tgroupTimes is a sorted, unique list of times +// assume a 1day range for now +off_t findTGroup (unsigned int start, int photcode) { + + // first find the matching photcode + int Ns = -1; + for (int i = 0; (Ns < 0) && (i < Nphotcodes); i++) { + int ecode = GetPhotcodeEquivCodebyCode (photcode); + if (photcodes[i][0].code != ecode) continue; + Ns = i; + } + if (Ns < 0) return -1; + + off_t Nlo, Nhi, N; + + // find the last tgroup before start + Nlo = 0; // first valid tgroupTimes value + Nhi = Ntgroup - 1; // last valid tgroupTimes value + + // if start is not in this range, return -1 + if (start < tgroupTimes[Nlo]) return (-1); + if (start > tgroupTimes[Nhi] + 86399) return (-1); // 1-day range for tgroup for now + + while (Nhi - Nlo > 4) { + N = 0.5*(Nlo + Nhi); + if (tgroupTimes[N] < start) { + Nlo = MAX(N, 0); + } else { + Nhi = MIN(N, Ntgroup - 1); + } + } + // we now have : tgroupTimes[Nlo] < start <= tgroupTimes[Nhi] + + // find a matched tgroup starting from Nlo, or return -1 + // tgroupTimes is a lower bound, the upper bound is tgroupTimes + 86399 (1d) + for (N = Nlo; N <= Nhi; N++) { + if (start < tgroupTimes[N]) continue; + if (start > tgroupTimes[N] + 86399) continue; // XXX since they are sorted, probably do not need this check + return (N*Nphotcodes + Ns); // we have the tgroupTime and photcodes, so we can construct the tgroup index + } + return (-1); +} + +void setMcalFromTGroups () { + + off_t i, j, im, Nimage; + Image *image; + + if (!TGROUP_ZEROPT) return; + + image = getimages (&Nimage, NULL); + + fprintf (stderr, "*** return Mcal from tgroup.Mcal to image.Mcal ***\n"); + + // copy the tgroup results to the images. set the tgroup Mcal to 0.0 since we have moved its + // impact to the images + for (i = 0; i < Ntgroup; i++) { + for (j = 0; j < tgroup[i].Nimage; j++) { + im = tgroup[i].image[j]; + 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; + image[im].dMagSys = tgroup[i].dMsys; + image[im].nFitPhotom = tgroup[i].nFitPhotom; + 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; + } +} + +void initTGroupBins (Catalog *catalog, int Ncatalog) { + + off_t i, j; + + /* measure -> tgroup */ + if (!TGROUP_ZEROPT) return; + + ALLOCATE (MeasureToTGroup, off_t *, Ncatalog); + for (i = 0; i < Ncatalog; i++) { + ALLOCATE (MeasureToTGroup[i], off_t, MAX (catalog[i].Nmeasure, 1)); + for (j = 0; j < catalog[i].Nmeasure; j++) MeasureToTGroup[i][j] = -1; + } +} + +void freeTGroupBins (int Ncatalog) { + + off_t i; + + /* measure -> tgroup */ + if (!TGROUP_ZEROPT) return; + + for (i = 0; i < Ncatalog; i++) { + free (MeasureToTGroup[i]); + } + free (MeasureToTGroup); +} + +int findTGroups (Catalog *catalog, int Ncatalog) { + + int i, ecode, found, Ns; + off_t j; + + if (!TGROUP_ZEROPT) return (FALSE); + // if we are calibrating by tgroup, redefine myDet == on one of my tgroups + + int Nmatch = 0; + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Nmeasure; j++) { + catalog[i].measureT[j].myDet = FALSE; // a detetion is not mine until proven otherwise + + if (TimeSelect) { + if (catalog[i].measureT[j].t < TSTART) continue; + if (catalog[i].measureT[j].t > TSTOP) continue; + } + + ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[j].photcode); + found = FALSE; + for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) { + if (ecode == photcodes[Ns][0].code) found = TRUE; + } + if (!found) continue; + + matchTGroups (catalog, j, i); + 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); +} + +void matchTGroups (Catalog *catalog, off_t meas, int cat) { + + MeasureTiny *measure = &catalog[cat].measureT[meas]; + + off_t ID = measure[0].imageID; // ID is a unique ID in DVO + off_t idx = getImageByID (ID); // idx is the sequence number in the subset list + if (idx == -1) { + if (VERBOSE2) fprintf (stderr, "missed measurement "OFF_T_FMT", %d\n", meas, cat); + return; + } + + off_t grpID = ImageToTGroup[idx]; + if (grpID < 0) { + // Image *image = getimage(idx); + // fprintf (stderr, "unmatched image %s\n", image[0].name); + // skip measurements from simple chips (not matched to a tgroup by definition) + return; + } + + // test to check we got the right match: + { + Image *image = getimage(idx); + // XXX we are now matching with just tzero. be careful for cameras with drift + // unsigned int imageStart = image[0].tzero - MAX(0.01*image[0].trate*image[0].NY, 1); + unsigned int imageStart = image[0].tzero; + if ((imageStart < tgroup[grpID].start) || (imageStart > tgroup[grpID].stop)) { + fprintf (stderr, "error in image to tgroup match\n"); + abort(); + } + } + + // this measurement is on one of my tgroups, mark it as mine. + catalog[cat].measureT[meas].myDet = TRUE; + MeasureToTGroup[cat][meas] = grpID; + + tgroup[grpID].catalog[tgroup[grpID].Nmeasure] = cat; + tgroup[grpID].measure[tgroup[grpID].Nmeasure] = meas; + tgroup[grpID].Nmeasure ++; + + if (tgroup[grpID].Nmeasure == tgroup[grpID].NMEASURE) { + tgroup[grpID].NMEASURE += 1000; + REALLOCATE (tgroup[grpID].catalog, off_t, tgroup[grpID].NMEASURE); + REALLOCATE (tgroup[grpID].measure, off_t, tgroup[grpID].NMEASURE); + } + return; +} + +float getMgrp (off_t meas, int cat, float airmass) { + + off_t i; + float value; + + if (!TGROUP_ZEROPT) return (0); + + // unassigned measurements belong to simple chips + i = MeasureToTGroup[cat][meas]; + if (i == -1) return (0.0); + + // if (tgroup[i].flags & IMAGE_BAD) return (NAN); + value = tgroup[i].McalPSF + tgroup[i].dKlam*(airmass - 1.0); + return (value); +} + +int getTGroupFlags (off_t meas, int cat) { + + if (!TGROUP_ZEROPT) return (0); + + // unassigned measurements belong to simple chips + int i = MeasureToTGroup[cat][meas]; + if (i == -1) return (0); + + return (tgroup[i].flags); +} + +typedef struct { + int Nfew; + int Nbad; + int Ncal; + int Nmos; + int Ngrid; + int Nrel; + int Nsys; + int Nskip; + off_t Nmax; + int PoorImages; + double *psfMagList; + double *psfAirList; + double *psfErrList; + double *brightMagList; + double *brightErrList; + double *kronMagList; + double *kronErrList; +} SetMgrpInfo; + +enum {THREAD_RUN, THREAD_DONE}; + +typedef struct { + int entry; + int state; + Catalog *catalog; + Image *image; + FlatCorrectionTable *flatcorr; + SetMgrpInfo info; +} ThreadInfo; + +int setMgrp_tgroup (TGroup *tgroup, off_t Nmos, Image *image, Catalog *catalog, SetMgrpInfo *info, FlatCorrectionTable *flatcorr); +void *setMgrp_worker (void *data); +int setMgrp_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr); + +void SetMgrpInfoInit (SetMgrpInfo *info, off_t Nmax, int allocLists, int PoorImages) { + info->Nfew = 0; + info->Nbad = 0; + info->Ncal = 0; + info->Nmos = 0; + info->Nrel = 0; + info->Ngrid = 0; + info->Nskip = 0; + info->Nsys = 0; + + info->Nmax = Nmax; + info->PoorImages = PoorImages; + + if (allocLists) { + ALLOCATE (info->psfMagList, double, Nmax); + ALLOCATE (info->psfAirList, double, Nmax); + ALLOCATE (info->psfErrList, double, Nmax); + ALLOCATE (info->kronMagList, double, Nmax); + ALLOCATE (info->kronErrList, double, Nmax); + ALLOCATE (info->brightMagList, double, Nmax); + ALLOCATE (info->brightErrList, double, Nmax); + } +} + +void SetMgrpInfoFree (SetMgrpInfo *info) { + free (info->psfMagList ); + free (info->psfAirList ); + free (info->psfErrList ); + free (info->kronMagList ); + free (info->kronErrList ); + free (info->brightMagList); + free (info->brightErrList); +} + +void SetMgrpInfoAccum (SetMgrpInfo *summary, SetMgrpInfo *results) { + summary->Nfew += results->Nfew ; + summary->Nsys += results->Nsys ; + summary->Nbad += results->Nbad ; + summary->Ncal += results->Ncal ; + summary->Nmos += results->Nmos ; + summary->Nrel += results->Nrel ; + summary->Ngrid += results->Ngrid; + summary->Nskip += results->Nskip; +} + +static int npass_output = 0; + +// mutex to lock setMgrp_worker operations +static pthread_mutex_t setMgrp_mutex = PTHREAD_MUTEX_INITIALIZER; +static int nextTGroup = 0; + +// we have an array of tgroups (tgroup, Ntgroup). we need to hand out tgroups one at a time to +// the worker threads as they need +off_t getNextTGroupForThread () { + + pthread_mutex_lock (&setMgrp_mutex); + if (nextTGroup >= Ntgroup) { + pthread_mutex_unlock (&setMgrp_mutex); + return (-1); + } + int thisTGroup = nextTGroup; + nextTGroup ++; + + pthread_mutex_unlock (&setMgrp_mutex); + return (thisTGroup); +} + +int setMgrp (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) { + + off_t i, N, Nmax; + Image *image; + + if (!TGROUP_ZEROPT) return (FALSE); + + // plots cannot be done in a threaded context (the plot commands collide) + // so do not run setMgrp in threaded mode if PLOTSTUFF is set + if (NTHREADS && !PLOTSTUFF) { + int status = setMgrp_threaded (catalog, PoorImages, flatcorr); + return status; + } + + image = getimages (&N, NULL); + + if (PoorImages) { + // XXX use bad stars and measurements for PoorImages? or not? + // IMAGE_BAD = STAR_BAD = MEAS_BAD = 0; + IMAGE_BAD = 0; + } + + Nmax = 0; + for (i = 0; i < Ntgroup; i++) { + Nmax = MAX (Nmax, tgroup[i].Nmeasure); + } + + SetMgrpInfo info; + SetMgrpInfoInit (&info, Nmax, TRUE, PoorImages); + + // int savePlotDelay = PLOTDELAY; + // if (PLOTSTUFF) PLOTDELAY = 0.0; + + 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); + + if (PLOTSTUFF) { + fprintf (stdout, "press return\n"); + if (fscanf (stdin, "%*c") != 1) fprintf (stderr, "\n"); + } + + npass_output ++; + + fprintf (stderr, "%d tgroups marked having too few measurements (Nbad: %d, Ncal: %d, Ngrid: %d, Nrel: %d, Nsys: %d)\n", info.Nfew, info.Nbad, info.Ncal, info.Ngrid, info.Nrel, info.Nsys); + + if (PoorImages) { + IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP; + STAR_BAD = ID_OBJ_POOR | ID_OBJ_FEW; + MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA; + } + return (TRUE); +} + +# define NIGHT_TOOFEW 100 +# define NIGHT_GOOD_FRACTION 0.1 + +// 'tgroup' is a pointer to the current tgroup of interest (Nmos) +int setMgrp_tgroup (TGroup *myTGroup, off_t Ngrp, Image *image, Catalog *catalog, SetMgrpInfo *info, FlatCorrectionTable *flatcorr) { + + off_t j, NimageReal; + + StatType stats; + liststats_setmode (&stats, "INNER_WTMEAN"); + + double *psfMagList = info->psfMagList; + double *psfAirList = info->psfAirList; + double *psfErrList = info->psfErrList; + double *kronMagList = info->kronMagList; + double *kronErrList = info->kronErrList; + double *brightMagList = info->brightMagList; + double *brightErrList = info->brightErrList; + + // Image *imageReal = getimages (&NimageReal, NULL); returned pointer is not used + getimages (&NimageReal, NULL); + + /* on PoorImages run, skip good images */ + if (info->PoorImages) { + int bad = myTGroup[0].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP); + if (!bad) return TRUE; + } + + int Nsecfilt = GetPhotcodeNsecfilt (); + + int testImage = FALSE; + testImage |= (abs(myTGroup[0].start - 1245283200) < 10); + // testImage |= (abs(myTGroup[0].start - 1324103823) < 10); + + FILE *fout = NULL; + if (testImage) { + char filename[64]; + snprintf (filename, 64, "test.%05d.%02d.dat", (int) Ngrp, npass_output); + fout = fopen (filename, "w"); + } + + // number of stars to measure the bright-end scatter + int Nbright = 0; + + int N = 0; + for (j = 0; j < myTGroup->Nmeasure; j++) { + + off_t m = myTGroup->measure[j]; + off_t c = myTGroup->catalog[j]; + + // NOTE : we are only using Mcal == McalPSF in this function; we set McalAPER to McalPSF + if (fout) { + + float Mcal = getMcal (m, c, MAG_CLASS_PSF); // image zero point + float Mmos = getMmos (m, c); // mosaic zero point + // XXX these two should both be 0.0 if we are fitting tgroups : skip those calls above? + + float Mgrid = getMgrid (m, c); // camera offset (deprecated?) + float MrelPSF = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); // average magnitude + float Mflat = getMflat (m, c, flatcorr, catalog); // flat-field correction + + off_t n = catalog[c].measureT[m].averef; + + // magnitude for this measurement + float MsysPSF = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF); + + // for historical reasons, Mflat is defined with the wrong sign + float delta = MsysPSF - MrelPSF - Mcal - Mmos - Mgrid + Mflat; + + int isBad = (catalog[c].measureT[m].dbFlags & MEAS_BAD); + + 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); + } + + if (catalog[c].measureT[m].dbFlags & MEAS_BAD) { + info->Nbad ++; + continue; + } + float Mcal = getMcal (m, c, MAG_CLASS_PSF); + if (isnan(Mcal)) { + info->Ncal++; + continue; + } + float Mmos = getMmos (m, c); + if (isnan(Mmos)) { + info->Nmos++; + continue; + } + float Mgrid = getMgrid (m, c); + if (isnan(Mgrid)) { + info->Ngrid ++; + continue; + } + + // Mrel* is the average magnitude for this star. For PS1 stacks, we have too much + // PSF variability. We need to calibrate the PSF magnitudes separately from the + // Aperture-like magnitues. (We have an option to use the kron magnitudes or the + // other apertures here). I basically need to do this analysis separately for each + // magnitude type + + float MrelPSF = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); + if (isnan(MrelPSF)) { + info->Nrel ++; + continue; + } + float MrelKron = getMrel (catalog, m, c, MAG_CLASS_KRON, MAG_SRC_CHP); + + // 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) + + float Mflat = getMflat (m, c, flatcorr, catalog); + + off_t n = catalog[c].measureT[m].averef; + float MsysPSF = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF); + if (isnan(MsysPSF)) { + info->Nsys++; + continue; + } + float MsysKron = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_KRON); + + PhotCode *code = GetPhotcodebyCode (catalog[c].measureT[m].photcode); + if (!code) goto skip; + if (code->equiv < 1) goto skip; + int Nsec = GetPhotcodeNsec (code->equiv); + if (Nsec == -1) goto skip; + + skip: + assert (N < info->Nmax); + assert (N >= 0); + assert (Nbright < info->Nmax); + assert (Nbright >= 0); + + 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; + kronErrList[N] = psfErrList[N]; + if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) { + brightMagList[Nbright] = psfMagList[N]; + brightErrList[Nbright] = psfErrList[N]; + Nbright ++; + } + N++; + } + if (fout) { fclose (fout); } + + /* too few good measurements or too many bad measurements (skip in PoorImages run) */ + if (!info->PoorImages) { + int mark = (N < NIGHT_TOOFEW) || (N < NIGHT_GOOD_FRACTION*myTGroup->Nmeasure); + if (mark) { + // XXX is this right (effect of subset?) : image[tgroup[Ngrp].image[0]].name + if (VERBOSE2) { fprintf (stderr, "marked tgroup %s ("OFF_T_FMT"), (%d < %d) || (%d < %f*"OFF_T_FMT")\n", image[tgroup[Ngrp].image[0]].name, Ngrp, N, NIGHT_TOOFEW, N, NIGHT_GOOD_FRACTION, myTGroup->Nmeasure); } + myTGroup[0].flags |= ID_IMAGE_PHOTOM_FEW; + info->Nfew ++; + if (testImage) { + fprintf (stderr, "NOTE: *** marked test image poor : %d %d %d***\n", (int) N, (int) NIGHT_TOOFEW, (int) (NIGHT_GOOD_FRACTION*myTGroup->Nmeasure)); + } + } else { + myTGroup[0].flags &= ~ID_IMAGE_PHOTOM_FEW; + } + } + + 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); + + // for now, I have no reason to measure these separately for camera-level images + myTGroup[0].McalPSF = stats.mean; + myTGroup[0].dMcal = stats.error; + myTGroup[0].McalChiSq = stats.chisq; + myTGroup[0].nFitPhotom = N; + myTGroup[0].dKlam = dk; + + // drop this for tgroups? + liststats (kronMagList, kronErrList, NULL, N, &stats); + myTGroup[0].McalAPER = stats.mean; + + if (testImage) { + fprintf (stderr, "test image %d (%d) %f %f %d ... ", (int) Ngrp, myTGroup[0].start, stats.mean, stats.error, myTGroup[0].nFitPhotom); + } + + if (PLOTSTUFF) { + fprintf (stderr, "Mgrp: %6.3f %6.3f +/- %6.3f %5d %5d | %s\n", stats.mean, stats.median, stats.sigma, stats.Nmeas, N, image[tgroup[Ngrp].image[0]].name); + plot_setMcal (psfMagList, N, &stats, CLOUD_TOLERANCE); + } + + // bright end scatter + liststats (brightMagList, brightErrList, NULL, Nbright, &stats); + myTGroup[0].dMsys = stats.sigma; + + // keep this?? + if (myTGroup[0].McalPSF < -CLOUD_TOLERANCE) { + myTGroup[0].McalPSF = 0.0; + } + + if (testImage) { + fprintf (stderr, "%f %f : %f\n", myTGroup[0].McalPSF, myTGroup[0].dMsys, myTGroup[0].McalChiSq); + } + + return TRUE; +} + +int setMgrp_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) { + + int i; + off_t N; + + Image *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? + // IMAGE_BAD = STAR_BAD = MEAS_BAD = 0; + IMAGE_BAD = 0; + } + + off_t Nmax = 0; + for (i = 0; i < Ntgroup; i++) { + Nmax = MAX (Nmax, tgroup[i].Nmeasure); + } + + SetMgrpInfo summary; + SetMgrpInfoInit (&summary, Nmax, FALSE, PoorImages); + + pthread_attr_t attr; + pthread_attr_init (&attr); + pthread_attr_setdetachstate (&attr, PTHREAD_CREATE_DETACHED); + + pthread_t *threads; + ALLOCATE (threads, pthread_t, NTHREADS); + + ThreadInfo *threadinfo; + ALLOCATE (threadinfo, ThreadInfo, NTHREADS); + + // each time this function is called, we cycle through the available tgroups + // make sure we start at 0 + nextTGroup = 0;; + + // launch N worker threads + for (i = 0; i < NTHREADS; i++) { + threadinfo[i].entry = i; + threadinfo[i].state = THREAD_RUN; + threadinfo[i].catalog = catalog; + threadinfo[i].image = image; + threadinfo[i].flatcorr = flatcorr; + SetMgrpInfoInit (&threadinfo[i].info, Nmax, FALSE, PoorImages); + pthread_create (&threads[i], NULL, setMgrp_worker, &threadinfo[i]); + } + pthread_attr_destroy (&attr); + + // wait until all threads have finished + while (1) { + int allDone = TRUE; + for (i = 0; i < NTHREADS; i++) { + if (threadinfo[i].state == THREAD_RUN) allDone = FALSE; + } + if (allDone) { + break; + } + usleep (500000); + } + + // all threads are done, free the threads array and grab the info + free (threads); + + // report stats & summary from the threads + for (i = 0; i < NTHREADS; i++) { + fprintf (stderr, "setMgrp thread %d : %d tgroups marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nrel: %d, Nsys: %d), %d partials skipped\n", + i, + threadinfo[i].info.Nfew, + threadinfo[i].info.Nbad, + threadinfo[i].info.Ncal, + threadinfo[i].info.Nmos, + threadinfo[i].info.Ngrid, + threadinfo[i].info.Nrel, + threadinfo[i].info.Nsys, + threadinfo[i].info.Nskip); + SetMgrpInfoAccum (&summary, &threadinfo[i].info); + } + fprintf (stderr, "total : %d tgroups marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nrel: %d, Nsys: %d), %d partials skipped\n", + summary.Nfew, + summary.Nbad, + summary.Ncal, + summary.Nmos, + summary.Ngrid, + summary.Nrel, + summary.Nsys, + summary.Nskip); + free (threadinfo); + + npass_output ++; + + if (PoorImages) { + IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP; + STAR_BAD = ID_OBJ_POOR | ID_OBJ_FEW; + MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA; + } + return TRUE; +} + +void *setMgrp_worker (void *data) { + + ThreadInfo *threadinfo = data; + + SetMgrpInfo results; + SetMgrpInfoInit (&results, threadinfo->info.Nmax, TRUE, threadinfo->info.PoorImages); // allocate list, dlist arrays here + + while (1) { + + off_t i = getNextTGroupForThread(); + if (i == -1) { + threadinfo->state = THREAD_DONE; + return NULL; + } + + Catalog *catalog = threadinfo->catalog; + FlatCorrectionTable *flatcorr = threadinfo->flatcorr; + Image *image = threadinfo->image; + + setMgrp_tgroup (&tgroup[i], i, image, catalog, &results, flatcorr); + SetMgrpInfoAccum (&threadinfo->info, &results); + } + + SetMgrpInfoFree (&results); + return NULL; +} + +StatType statsTGroupM (Catalog *catalog) { + OHANA_UNUSED_PARAM(catalog); + + off_t i; + int n; + double *list, *dlist; + StatType stats; + + bzero (&stats, sizeof (StatType)); + if (!TGROUP_ZEROPT) return (stats); + + ALLOCATE (list, double, Ntgroup); + ALLOCATE (dlist, double, Ntgroup); + + n = 0; + for (i = 0; i < Ntgroup; i++) { + if (tgroup[i].flags & IMAGE_BAD) continue; + list[n] = tgroup[i].McalPSF; + dlist[n] = 1; + n++; + } + + liststats_setmode (&stats, "MEAN"); + + liststats (list, dlist, NULL, n, &stats); + free (list); + free (dlist); + return (stats); +} + +StatType statsTGroupdM (Catalog *catalog) { + OHANA_UNUSED_PARAM(catalog); + + off_t i, n; + double *list, *dlist; + StatType stats; + + bzero (&stats, sizeof (StatType)); + if (!TGROUP_ZEROPT) return (stats); + + ALLOCATE (list, double, Ntgroup); + ALLOCATE (dlist, double, Ntgroup); + + n = 0; + for (i = 0; i < Ntgroup; i++) { + if (tgroup[i].flags & IMAGE_BAD) continue; + list[n] = tgroup[i].dMcal; + dlist[n] = 1; + n++; + } + + liststats_setmode (&stats, "MEAN"); + + liststats (list, dlist, NULL, n, &stats); + free (list); + free (dlist); + return (stats); +} + +StatType statsTGroupN (Catalog *catalog, FlatCorrectionTable *flatcorr) { + + off_t i, j, m, c, n, N; + double *list, *dlist; + float Mcal, Mgrid, Mrel, Mmos; + StatType stats; + + bzero (&stats, sizeof (StatType)); + if (!TGROUP_ZEROPT) return (stats); + + ALLOCATE (list, double, Ntgroup); + ALLOCATE (dlist, double, Ntgroup); + + n = 0; + for (i = 0; i < Ntgroup; i++) { + if (tgroup[i].flags & IMAGE_BAD) continue; + + N = 0; + for (j = 0; j < tgroup[i].Nmeasure; j++) { + + m = tgroup[i].measure[j]; + c = tgroup[i].catalog[j]; + + Mcal = getMcal (m, c, MAG_CLASS_PSF); + if (isnan(Mcal)) continue; + Mmos = getMmos (m, c); + if (isnan(Mmos)) continue; + Mgrid = getMgrid (m, c); + if (isnan(Mgrid)) continue; + Mrel = getMrel (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); + if (isnan(Mrel)) continue; + N++; + } + list[n] = N; + dlist[n] = 1; + n++; + } + // fprintf (stderr, "Ntgroup: "OFF_T_FMT", n: "OFF_T_FMT"\n", Ntgroup, n); + + liststats_setmode (&stats, "MEAN"); + + liststats (list, dlist, NULL, n, &stats); + free (list); + free (dlist); + return (stats); +} + +StatType statsTGroupX (Catalog *catalog) { + OHANA_UNUSED_PARAM(catalog); + + off_t i, n; + double *list, *dlist; + StatType stats; + + bzero (&stats, sizeof (StatType)); + if (!TGROUP_ZEROPT) return (stats); + + ALLOCATE (list, double, Ntgroup); + ALLOCATE (dlist, double, Ntgroup); + + n = 0; + for (i = 0; i < Ntgroup; i++) { + if (tgroup[i].flags & IMAGE_BAD) continue; + if (KEEP_UBERCAL && (tgroup[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue; + + list[n] = tgroup[i].McalChiSq; + dlist[n] = 1; + n++; + } + + liststats_setmode (&stats, "MEAN"); + + liststats (list, dlist, NULL, n, &stats); + free (list); + free (dlist); + return (stats); +} + +/* mark tgroup if: abs(Mcal - ) too large, dMcal too large */ +void clean_tgroups () { + + off_t i, N, mark, Nmark, Nscatter, Noffset; + double *mlist, *slist, *dlist; + double MaxOffset, MedOffset, MaxScatter; + + if (!TGROUP_ZEROPT) return; + + if (VERBOSE) fprintf (stderr, "marking poor tgroups\n"); + + ALLOCATE (mlist, double, Ntgroup); + ALLOCATE (slist, double, Ntgroup); + ALLOCATE (dlist, double, Ntgroup); + + for (i = N = 0; i < Ntgroup; i++) { + if (tgroup[i].flags & IMAGE_BAD) continue; + mlist[N] = tgroup[i].McalPSF; + slist[N] = tgroup[i].dMcal; + dlist[N] = 1; + N++; + } + + StatType stats; + liststats_setmode (&stats, "MEAN"); + + liststats (mlist, dlist, NULL, N, &stats); + MaxOffset = MAX (IMAGE_OFFSET, 2*stats.sigma); + MedOffset = stats.median; + + liststats (slist, dlist, NULL, N, &stats); + MaxScatter = MAX (IMAGE_SCATTER, 2*stats.median); + fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset); + + Nmark = Nscatter = Noffset = 0; + for (i = 0; i < Ntgroup; i++) { + // if we are keeping ubercal sacrosanct, then we should not be allowed to break them... + if (KEEP_UBERCAL && (tgroup[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue; + + if (tgroup[i].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP)) continue; + + mark = FALSE; + if (tgroup[i].dMcal > MaxScatter) { + mark = TRUE; + Nscatter ++; + } + if (fabs(tgroup[i].McalPSF - MedOffset) > MaxOffset) { + mark = TRUE; + Noffset ++; + } + if (mark) { + Nmark ++; + tgroup[i].flags |= ID_IMAGE_PHOTOM_POOR; + } else { + tgroup[i].flags &= ~ID_IMAGE_PHOTOM_POOR; + } + } + + fprintf (stderr, OFF_T_FMT" tgroups marked poor ("OFF_T_FMT" scatter, "OFF_T_FMT" offset)\n", Nmark, Nscatter, Noffset); + + free (mlist); + free (slist); + free (dlist); +} + +void plot_tgroup_fields (Catalog *catalog) { + + off_t i, j, m, c, N, Nimage; + double *xlist, *ylist; + char string[64]; + Graphdata graphdata; + + if (!TGROUP_ZEROPT) return; + + // Image *image = getimages (&Nimage, NULL); returned value ignored + getimages (&Nimage, NULL); + + N = 0; + for (i = 0; i < Ntgroup; i++) { + N = MAX (N, tgroup[i].Nmeasure); + } + + ALLOCATE (xlist, double, N); + ALLOCATE (ylist, double, N); + + for (i = 0; i < Ntgroup; i++) { + N = 0; + for (j = 0; j < tgroup[i].Nmeasure; j++) { + + m = tgroup[i].measure[j]; + c = tgroup[i].catalog[j]; + + if (catalog[c].measureT[m].dbFlags & (ID_MEAS_AREA | ID_MEAS_NOCAL)) continue; + + // ave = catalog[c].measureT[m].averef; + xlist[N] = catalog[c].measureT[m].R; + ylist[N] = catalog[c].measureT[m].D; + N++; + } + + sprintf (string, "TGroup "OFF_T_FMT, i); + plot_defaults (&graphdata); + plot_list (&graphdata, xlist, ylist, N, string, NULL); + } + + free (ylist); + free (xlist); +} + +void plot_tgroups () { + + off_t i, bin; + double *xlist, *Mlist, *dlist; + Graphdata graphdata; + + if (!TGROUP_ZEROPT) return; + + ALLOCATE (xlist, double, Ntgroup); + ALLOCATE (dlist, double, Ntgroup); + ALLOCATE (Mlist, double, Ntgroup); + + for (i = 0; i < Ntgroup; i++) { + Mlist[i] = tgroup[i].McalPSF; + dlist[i] = tgroup[i].dMcal; + xlist[i] = tgroup[i].start; + } + + plot_defaults (&graphdata); + graphdata.xmin = 0.95; + graphdata.xmax = 2.50; + graphdata.ymin = PlotdMmin; + graphdata.ymax = PlotdMmax; + plot_list (&graphdata, xlist, Mlist, Ntgroup, "airmass vs Mcal", "%s.airmass.png", OUTROOT); + plot_defaults (&graphdata); + graphdata.size = 1.5; + graphdata.ptype = 7; + plot_list (&graphdata, Mlist, dlist, Ntgroup, "Mcal vs dMcal", "%s.MdM.png", OUTROOT); + +# define NBIN 200 + REALLOCATE (xlist, double, NBIN); + REALLOCATE (Mlist, double, NBIN); + + /**** dMcal histgram ****/ + for (i = 0; i < NBIN; i++) xlist[i] = 0.00005*i; + bzero (Mlist, NBIN*sizeof(double)); + for (i = 0; i < Ntgroup; i++) { + bin = tgroup[i].dMcal / 0.00005; + bin = MAX (0, MIN (NBIN - 1, bin)); + Mlist[bin] += 1.0; + } + plot_defaults (&graphdata); + graphdata.style = 1; + plot_list (&graphdata, xlist, Mlist, NBIN, "dMcal hist", "%s.dMcalhist.png", OUTROOT); + + free (dlist); + free (xlist); + free (Mlist); +} + Index: trunk/Ohana/src/relphot/src/args.c =================================================================== --- trunk/Ohana/src/relphot/src/args.c (revision 41388) +++ trunk/Ohana/src/relphot/src/args.c (revision 41390) @@ -336,4 +336,35 @@ 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); + TGROUP_FILENAME = strcreate(argv[N]); + remove_argument (N, &argc, argv); + FREEZE_IMAGES = 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); } @@ -519,4 +550,9 @@ if (argc != 1) relphot_usage (); + if (TGROUP_ZEROPT) { + loadTGroups (TGROUP_FILENAME); + free (TGROUP_FILENAME); + } + return mode; } @@ -542,4 +578,6 @@ FreePhotcodeTable(); + freeTGroups(); + ohana_memcheck (VERBOSE); ohana_memdump (VERBOSE); Index: trunk/Ohana/src/relphot/src/liststats.c =================================================================== --- trunk/Ohana/src/relphot/src/liststats.c (revision 41388) +++ trunk/Ohana/src/relphot/src/liststats.c (revision 41390) @@ -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: trunk/Ohana/src/relphot/src/load_images.c =================================================================== --- trunk/Ohana/src/relphot/src/load_images.c (revision 41388) +++ trunk/Ohana/src/relphot/src/load_images.c (revision 41390) @@ -110,4 +110,6 @@ MARKTIME("init images: %f sec\n", dtime); + initTGroups (subset, LineNumber, Nsubset); + return TRUE; } Index: trunk/Ohana/src/relphot/src/reload_catalogs.c =================================================================== --- trunk/Ohana/src/relphot/src/reload_catalogs.c (revision 41388) +++ trunk/Ohana/src/relphot/src/reload_catalogs.c (revision 41390) @@ -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: trunk/Ohana/src/relphot/src/relphot_images.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_images.c (revision 41388) +++ trunk/Ohana/src/relphot/src/relphot_images.c (revision 41390) @@ -70,4 +70,6 @@ initMosaicBins (catalog, Ncatalog, TRUE); + initTGroupBins (catalog, Ncatalog); + initGridBins (catalog, Ncatalog); initMrel (catalog, Ncatalog); @@ -77,4 +79,7 @@ findMosaics (catalog, Ncatalog, TRUE); /* also sets Grid values */ + MARKTIME("-- set up mosaic indexes: %f sec\n", dtime); + + findTGroups (catalog, Ncatalog); MARKTIME("-- set up mosaic indexes: %f sec\n", dtime); @@ -129,4 +134,6 @@ setMmos (catalog, FALSE, flatcorr); + setMgrp (catalog, FALSE, flatcorr); + setMgrid (catalog, flatcorr); MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime); @@ -151,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 (); @@ -172,7 +180,9 @@ setMcal (catalog, onlyPoorImages, flatcorr); setMmos (catalog, onlyPoorImages, flatcorr); + setMgrp (catalog, onlyPoorImages, flatcorr); MARKTIME("-- finalize Mcal values: %f sec\n", dtime); - setMcalFinal (); // copy per-mosaic calibrations to the images + setMcalFromMosaics (); // copy per-mosaic calibrations to the images + setMcalFromTGroups (); // copy per-tgroup calibrations to the images if (SAVE_IMAGE_UPDATES) { @@ -249,4 +259,5 @@ freeImageBins (Ncatalog, TRUE); freeMosaicBins (Ncatalog, TRUE); + freeTGroupBins (Ncatalog); freeGridBins (Ncatalog); @@ -255,5 +266,6 @@ // If nloop == 0, the above pass is not performed, in which case // the Mcal values passed to the mosaics are not returned to the images... - setMcalFinal (); // copy per-mosaic calibrations to the images + setMcalFromMosaics (); // copy per-mosaic calibrations to the images + setMcalFromTGroups (); // copy per-tgroup calibrations to the images } Index: trunk/Ohana/src/relphot/src/relphot_parallel_images.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_parallel_images.c (revision 41388) +++ trunk/Ohana/src/relphot/src/relphot_parallel_images.c (revision 41390) @@ -116,4 +116,5 @@ setMcal (catalog, FALSE, flatcorr); setMmos (catalog, FALSE, flatcorr); + // setMgrp (catalog, FALSE, flatcorr); XXX think this through: this may not make sense for tgroups MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime); @@ -152,7 +153,9 @@ setMcal (catalog, onlyPoorImages, flatcorr); setMmos (catalog, onlyPoorImages, flatcorr); + // setMgrp (catalog, onlyPoorImages, flatcorr); XXX see above MARKTIME("-- finalize Mcal values: %f sec\n", dtime); - setMcalFinal (); // copy per-mosaic calibrations to the images + setMcalFromMosaics (); // copy per-mosaic calibrations to the images + // XXX keep? setMcalFromTGroups (); // copy per-tgroup calibrations to the images share_image_mags (regionHosts, -1); Index: trunk/Ohana/src/relphot/src/setMrelCatalog.c =================================================================== --- trunk/Ohana/src/relphot/src/setMrelCatalog.c (revision 41388) +++ trunk/Ohana/src/relphot/src/setMrelCatalog.c (revision 41390) @@ -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? } Index: trunk/Ohana/src/relphot/test/dvo.config =================================================================== --- trunk/Ohana/src/relphot/test/dvo.config (revision 41390) +++ trunk/Ohana/src/relphot/test/dvo.config (revision 41390) @@ -0,0 +1,17 @@ +IMAGE_SCATTER 0.075 # mark images POOR if stdev(Mcal) > IMAGE_SCATTER +STAR_SCATTER 0.005 # mark stars POOR if stdev(Mrel) > STAR_SCATTER +IMAGE_OFFSET 0.100 # mark images POOR if abs(delta(Mcal)) > IMAGE_OFFSET +STAR_CHISQ 10.000 # mark stars POOR if Xm > STAR_CHISQ +STAR_TOOFEW 2 # mark star FEW if N(good) < STAR_TOOFEW +GRID_TOOFEW 10 # keep grid FROZEN if N(good) < GRID_TOOFEW +IMAGE_TOOFEW 10 # mark image FEW if N(good) < IMAGE_TOOFEW +IMAGE_GOOD_FRACTION 0.050 # mark image FEW if N(good) < IMAGE_GOOD_FRACTION * Nstars + +SCATTER_LIM 15.000 +MAG_LIM 20.000 # select stars brighter than this for relphot analysis +SIGMA_LIM 0.015 # select measurements with dM < SIGMA_LIM for relphot analysis + +RELPHOT_GRID_BINNING 256 +RELPHOT_GRID_X 4 +RELPHOT_GRID_Y 4 + Index: trunk/Ohana/src/relphot/test/ptolemy.dvo =================================================================== --- trunk/Ohana/src/relphot/test/ptolemy.dvo (revision 41390) +++ trunk/Ohana/src/relphot/test/ptolemy.dvo (revision 41390) @@ -0,0 +1,62 @@ + +macro mkptolemy.rc + # generate a fake ptolemy.rc (to make this a contained test) + exec rm -f ptolemy.rc + output ptolemy.rc + echo "\#" + echo PHOTCODE_FILE dvo.photcodes + echo GSCFILE GSCregions.tbl + echo ZERO_PT 25.0 + echo ADDSTAR_RADIUS 1.0 + echo CATMODE SPLIT + echo CATFORMAT PS1_V5 + echo SKY_DEPTH 4 + echo CAMERA gpc1 + echo MOSAICNAME GPC1 + echo CAMERA_CONFIG dvo.layout + echo JD-KEYWORD NONE + echo MJD-KEYWORD MJD-OBS + echo CCDNUM-KEYWORD EXTNAME + echo EXPTIME-KEYWORD EXPTIME + echo AIRMASS-KEYWORD AIRMASS + echo "input dvo.config" + output stdout + + # we need to define a photcode table consistent with the fake construction + $zpt_nominal:g = 24.0 + $zpt_nominal:r = 24.5 + + # klam is negative, so klam*(airmass - 1) increases the zero point + $klam_nominal:g = -0.15 + $klam_nominal:r = -0.10 + + local found + file GSCregions.tbl found + if (not($found)) + file GSCregions.tbl.gz found + if (not($found)) + echo "missing GSCregions.tbl file or gzipped version" + break + end + exec gunzip -c GSCregions.tbl.gz > GSCregions.tbl + end + + # generate a fake dvo.photcode + exec rm -f dvo.photcodes + output dvo.photcodes + echo "\# code name type zero slope offset c1 c2 slope zero equiv sys scale scale sys poor bad poor bad" + echo "1 g sec 0.000 0.000 0.000 1 3 0.0000 0 1051 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo "2 r sec 0.000 0.000 0.000 2 3 0.0000 0 1052 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo "10000 GPC1.g.XY00 dep $zpt_nominal:g $klam_nominal:g 0.000 - - 0.0000 0 1 0.015 0.000 1.000 0.000 0x0000 0x0000 0xe0440130 0x1003bc88" + echo "10001 GPC1.g.XY01 dep $zpt_nominal:g $klam_nominal:g 0.000 - - 0.0000 0 1 0.015 0.000 1.000 0.000 0x0000 0x0000 0xe0440130 0x1003bc88" + echo "10010 GPC1.g.XY10 dep $zpt_nominal:g $klam_nominal:g 0.000 - - 0.0000 0 1 0.015 0.000 1.000 0.000 0x0000 0x0000 0xe0440130 0x1003bc88" + echo "10011 GPC1.g.XY11 dep $zpt_nominal:g $klam_nominal:g 0.000 - - 0.0000 0 1 0.015 0.000 1.000 0.000 0x0000 0x0000 0xe0440130 0x1003bc88" + echo "10100 GPC1.r.XY00 dep $zpt_nominal:r $klam_nominal:r 0.000 - - 0.0000 0 2 0.015 0.000 1.000 0.000 0x0000 0x0000 0xe0440130 0x1003bc88" + echo "10101 GPC1.r.XY01 dep $zpt_nominal:r $klam_nominal:r 0.000 - - 0.0000 0 2 0.015 0.000 1.000 0.000 0x0000 0x0000 0xe0440130 0x1003bc88" + echo "10110 GPC1.r.XY10 dep $zpt_nominal:r $klam_nominal:r 0.000 - - 0.0000 0 2 0.015 0.000 1.000 0.000 0x0000 0x0000 0xe0440130 0x1003bc88" + echo "10111 GPC1.r.XY11 dep $zpt_nominal:r $klam_nominal:r 0.000 - - 0.0000 0 2 0.015 0.000 1.000 0.000 0x0000 0x0000 0xe0440130 0x1003bc88" + echo "2011 2MASS_J ref 0.000 0.000 0.000 - - 0.0000 0 6 0.080 0.000 2.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo "3001 SYNTH.g ref 0.000 0.000 0.000 - - 0.0000 0 1 nan 0.000 2.000 0.600 0x0000 0x0000 0x0000 0x0000" + echo "2020 TYCHO_B ref 0.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + output stdout +end Index: trunk/Ohana/src/relphot/test/relphot.basic.dvo =================================================================== --- trunk/Ohana/src/relphot/test/relphot.basic.dvo (revision 41390) +++ trunk/Ohana/src/relphot/test/relphot.basic.dvo (revision 41390) @@ -0,0 +1,448 @@ +# -*-sh-*- + +$TAP_BREAK = 0 +$TAP_VERBOSE = 1 +input tap.dvo +input ptolemy.dvo + +# set globals +if (not($?PLOT)) set PLOT = 0 + +$USE_CELL_OFFSETS = 1; # turn on non-zero cell offsets +$SET_UC_RESID = 1; # write a non-zero resid value to the ubercal table (should not have an impact) +$DO_2MASS = 1 +$DO_SYNTH = 1 +$DO_TYCHO = 1 +$NC_CLOUDS = 0.20 +$UC_CLOUDS = 0.02 + +## full test set: generate a set of test files, insert into a catdir, check raw catalog, run relphot, check results +macro go + if ($0 != 3) + echo "go (fileroot) (catdir)" + break + end + + local i fileroot rootdir catdir + + $fileroot = $1 + $catdir = $2 + + dirname $1 -var rootdir + mkdir $rootdir + + break -auto off + exec rm -rf $catdir + exec rm -rf $fileroot* + break -auto on + + # define a set of photometric (ubercal) and non-photometric exposure times + init + + # generate the common list of stars (coords & mags) for all exposures + mkstars 1000 + + catdir $catdir + + tapPLAN {((mjd_uc[]*4 + mjd_nc[]*4) + (mjd_uc[] + mjd_nc[]) + (mjd_uc[] + mjd_nc[])) / 3} + + # generate the basic images and check they were correctly ingested by dvo + for i 0 mjd_uc[] + mkexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i $fileroot $RA_CENTER $DEC_CENTER + ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i raw + end + for i 0 mjd_nc[] + mkexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i $fileroot $RA_CENTER $DEC_CENTER + ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i raw + end + + # run relphot on the db and check that the images now match the expected values + tapEXEC relphot -images g,r -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir -D STAR_TOOFEW 1 -statmode WT_MEAN -cloud-limit 0.5 -update -nloop 12 + + # reset the per-filter ZPT_OFF:$filter values + for i 0 $filt_uc:n + $filter = $filt_uc:$i + $ZPT_OFF:$filter = NAN + end + for i 0 $filt_nc:n + $filter = $filt_nc:$i + $ZPT_OFF:$filter = NAN + end + + for i 0 mjd_uc[] + ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i relphot_uc + end + for i 0 mjd_nc[] + ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i relphot_nc + end + tapDONE +end + +# set various global variables +macro init + + # generates ptolemy.rc, dvo.photcodes, $zpt_nominal:g,r $klam_nominal:g,r + mkptolemy.rc + + $RA_CENTER = 10.0 + $DEC_CENTER = 20.0 + $PLATE_SCALE = 0.25 + $NSEASON = 3 + $NFILTER = 3 + $NCHIP_X = 2 + $NCHIP_Y = 2 + $NCELL_X = 2 + $NCELL_Y = 2 + $CHIP_DX = 1000 + $CHIP_DY = 1000 + $DEFECT_FRAC = 0.05 + $OFFSET_FRAC_UC = 0.0 + $OFFSET_FRAC_NC = 0.5 + + # we have two sets of images: ubercal (photometric) and not-ubercaled + + # ********* UBERCAL IMAGES *********** + + # the number of exposures is defined by filt_uc:n + list filt_uc -split g g g g r r r r + + # tmpseq is used to generate vectors in this function + create tmpseq 0 $filt_uc:n + + # exptime and airmasses for uc images + set exptime_uc = 10.0 + zero(tmpseq) + +# XXX range of airmass + set airmass_uc = 1.3 + zero(tmpseq) + set klam_uc = zero(tmpseq) + set zpt_uc = zero(tmpseq) + + # ubercal zero points are defined as ZP_nominal + 2.5log(exptime) + K*(airmass - 1.0) + # note that K is defined as a negative value (is this sensible?) + + # zero points and airmass slopes for these sequences + for i 0 $filt_uc:n + klam_uc[$i] = $klam_nominal:$filt_uc:$i + zpt_uc[$i] = $zpt_nominal:$filt_uc:$i + 2.5*log(exptime_uc[$i]) + klam_uc[$i]*(airmass_uc[$i] - 1.0) + $UC_CLOUDS*(rnd(klam_uc[$i]) - 0.5) + $ZPT_OFF:$filt_uc:$i = NAN + end + set mjd_uc = 55000.00 + 0.01*(tmpseq % 3) + 10.0*int(tmpseq / 3) + + # ********* NON-UBERCAL IMAGES *********** + + # mjd and zpt values for not-ubercal'ed images + list filt_nc -split g g g g r r r r + + # tmpseq is used to generate vectors in this function + create tmpseq 0 $filt_nc:n + + # place this within a valid season (55000.0 - 55010.0 - 55020.0 - 55030.0) + set exptime_nc = 10.0 + zero(tmpseq) + set airmass_nc = 1.6 + zero(tmpseq) + set klam_nc = zero(tmpseq) + set zpt_nc = zero(tmpseq) + + # airmass slopes for these sequnece + for i 0 $filt_nc:n + klam_nc[$i] = $klam_nominal:$filt_nc:$i + zpt_nc[$i] = $zpt_nominal:$filt_nc:$i + 2.5*log(exptime_nc[$i]) + klam_nc[$i]*(airmass_nc[$i] - 1.0) + $NC_CLOUDS*(rnd(klam_nc[$i]) - 0.5) + $ZPT_OFF:$filt_uc:$i = NAN + end + set mjd_nc = 55000.10 + 0.01*(tmpseq % 3) + 10.0*int(tmpseq / 3) + + delete tmpseq +end + +macro try.relphot + if ($0 != 3) + echo "go (fileroot) (catdir)" + break + end + + local i fileroot rootdir catdir + + $fileroot = $1 + $catdir = $2 + + # run relphot on the db and check that the images now match the expected values + tapEXEC relphot -images g,r -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir -D STAR_TOOFEW 1 -statmode WT_MEAN -cloud-limit 0.5 -update -nloop 12 + + # reset the per-filter ZPT_OFF:$filter values + for i 0 $filt_uc:n + $filter = $filt_uc:$i + $ZPT_OFF:$filter = NAN + end + for i 0 $filt_nc:n + $filter = $filt_nc:$i + $ZPT_OFF:$filter = NAN + end + + for i 0 mjd_uc[] + # ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i relphot_uc + ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i relphot_uc + end + for i 0 mjd_nc[] + # ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i relphot_nc + ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i relphot_nc + end +end + +macro ckexposure + if ($0 != 8) + echo "ckexposure (catdir) (mjd) (zpt) (exptime) (airmass) (filter) (mode)" + echo " mode == raw or corr" + break + end + + local CATDIR MJD_IMAGE ZPT_REAL ZPT_NOMINAL EXPTIME FILTER AIRMASS + + $CATDIR = $1 + $MJD_IMAGE = $2 + $ZPT_REAL = $3 + $EXPTIME = $4 + $AIRMASS = $5 + $FILTER = $6 + $FUNC_MODE = $7 + + # XXX need a function to extract the nominal zpt for a given filter / photcode from the db + $ZPT_NOMINAL = $zpt_nominal:$FILTER + $KLAM_NOMINAL = $klam_nominal:$FILTER + + $TIMEFORMAT = mjd + $TIMEREF = 0.0 + + catdir $CATDIR + + # assume we still have stars_ra, stars_dec, stars_mag in hand + if ($PLOT) + dev -n 0 + region $RA_CENTER $DEC_CENTER 0.2 + images + pmeasure -all -m 15 20 + end + + skyregion {$RA_CENTER - 0.2/dcos($DEC_CENTER)} {$RA_CENTER + 0.2/dcos($DEC_CENTER)} {$DEC_CENTER - 0.2} {$DEC_CENTER + 0.2} + + ## calculate the average zero point offset + if (isnan($ZPT_OFF:$FILTER)) + avextract ra dec $FILTER + match2d -closest ra dec stars_ra stars_dec 0.001 -index1 index1 -index2 index2 + reindex stars_tru_m = stars_mag using index1 + set dm_tru = $FILTER - stars_tru_m + vstat -q dm_tru + $ZPT_OFF:$FILTER = $MEDIAN + end + + mextract ra dec mag xccd yccd dbflags photflags time mag:sys mag:rel mag:inst mag:ave where (abs(time - $MJD_IMAGE) < 0.0001) + + match2d -closest ra dec stars_ra stars_dec 0.001 -index1 index1 -index2 index2 + + reindex stars_ra_m = stars_ra using index1 + reindex stars_dec_m = stars_dec using index1 + reindex stars_mag_m = stars_mag using index1 + set dr = 3600*(ra - stars_ra_m) + set dd = 3600*(dec - stars_dec_m) + set dm = mag - stars_mag_m + + set good_photom = (not(photflags & 0x02)) + foreach field dm xccd yccd + subset $field\_good = $field if good_photom + end + + if ($PLOT) + dev -n 1 + lim mag dm; clear; box; plot mag dm; plot -c red -pt 7 mag dm where (photflags & 0x00000002) + dev -n 2 + lim xccd dm; clear; box; plot xccd dm; plot -c red -pt 7 xccd dm where (photflags & 0x00000002) + dev -n 3 + lim ra dm; clear; box; plot ra dm; plot -c red -pt 7 ra dm where (photflags & 0x00000002) + echo "red points are known bad photometry" + end + + $ZPT_REAL_NORM = $ZPT_REAL - 2.5*log($EXPTIME) - $KLAM_NOMINAL*($AIRMASS - 1.0) + + if ("$FUNC_MODE" == "raw") + set dm_adjust = dm_good - $ZPT_NOMINAL + $ZPT_REAL_NORM + + vstat -q dm_adjust + sprintf line "addstar raw ZP %8.2f $FILTER %6.3f %6.3f %6.3f (%7.4f): %7.4f %7.4f : %6.3f" $MJD_IMAGE $ZPT_REAL $ZPT_REAL_NORM $ZPT_NOMINAL {$ZPT_REAL_NORM - $ZPT_NOMINAL} $MEAN $SIGMA dm_adjust[] + tapOK {abs($MEAN) < 0.005} "RAW : $line" + return + end + + if (("$FUNC_MODE" == "relphot_uc") || ("$FUNC_MODE" == "relphot_nc")) + + ## relphot does not set an absolute zero point, so we need to find that and subtract + set dm_adjust = dm_good - $ZPT_OFF:$FILTER + + vstat -q dm_adjust + sprintf line "relphot ZP %8.2f $FILTER %6.3f %6.3f %6.3f (%7.4f) : %7.4f %7.4f : %3d" $MJD_IMAGE $ZPT_REAL $ZPT_REAL_NORM $ZPT_NOMINAL {$ZPT_REAL_NORM - $ZPT_NOMINAL} $MEAN $SIGMA dm_adjust[] + + if ("$FUNC_MODE" == "relphot_uc") + tapOK {abs($MEAN) < 0.001} "MEAN : $line" + tapOK {abs($SIGMA) < 0.001} "SIGMA: $line" + end + if ("$FUNC_MODE" == "relphot_nc") + # XXX NOTE relaxed condition : bad photometry points are partly clipped but still cause a bias + tapOK {abs($MEAN) < 0.020} "MEAN : $line" + tapOK {abs($SIGMA) < 0.020} "SIGMA: $line" + end + + if ($PLOT) + cursor 1 + if ("$KEY" == "Escape") break + end + return + end + + echo "ERROR: unknown mode for ckexposure : $FUNC_MODE" +end + +macro mkexposure + if ($0 != 10) + echo "mkexposure (catdir) (mjd) (zpt) (exptime) (airmass) (filter) (fileroot) (ra) (dec)" + break + end + + local ix iy date time datetime ra dec ROOT RAo DECo ZPT MJD FILTER EXPTIME AIRMASS CATDIR cmd + + $CATDIR = $1 + $MJD = $2 + $ZPT = $3 + $EXPTIME = $4 + $AIRMASS = $5 + $FILTER = $6 + $ROOT = $7 + $RAo = $8 + $DECo = $9 + + $TIMEFORMAT = mjd + $TIMEREF = 0.0 + + # XX ctimes -ref $MJD -var datetime + # XX substr $datetime 0 10 date + # XX substr $datetime 11 8 time + + for ix 0 $NCHIP_X + for iy 0 $NCHIP_Y + $dx = $CHIP_DX * $PLATE_SCALE * ($ix - 0.5*$NCHIP_X + 0.5) + $dy = $CHIP_DY * $PLATE_SCALE * ($iy - 0.5*$NCHIP_Y + 0.5) + $ra = $RAo - $dx / 3600.0 / dcos($DECo) + $dec = $DECo - $dy / 3600.0 + # echo $ra $dec $dx $dy + mkinput test.in.txt $ra $dec $ZPT + + # ra,dec is the center of this chip + local options + $options = -mjd $MJD + $options = $options -radec $ra $dec + + #$options = $options -type PS1_V2 + #$options = $options -type PS1_V3 + #$options = $options -type PS1_V4 + $options = $options -type PS1_V5 + + $options = $options -coords + $options = $options -photcode GPC1.$FILTER.XY$ix\$iy + $options = $options -no-noise + $options = $options -flags + $options = $options -bad-psfqf-frac $DEFECT_FRAC + $options = $options -size $CHIP_DX $CHIP_DY + $options = $options -crpix {0.5*$CHIP_DX} {0.5*$CHIP_DY} + $options = $options -airmass $AIRMASS + $options = $options -exptime $EXPTIME + + tapEXEC mkcmf test.in.txt $ROOT.$ix.$iy.cmf $options + + # the fake images have inconsistent ra,dec and airmass,sidtime values + tapEXEC addstar $ROOT.$ix.$iy.cmf -D CATDIR $CATDIR -D CAMERA gpc1 -quick-airmass + end + end +end + +# make a simple input file for mkcmf +macro mkinput + if ($0 != 5) + echo "mkinput (filename) (ra) (dec) (zpt)" + break + end + + # chip coordinate of the stars + # ra,dec is at the chip center, which corresponds to pixel (CHIP_DX,CHIP_DY)/2 + set stars_X = (stars_ra - $2) * 3600.0 * dcos($3) / $PLATE_SCALE + $CHIP_DX*0.5 + set stars_Y = (stars_dec - $3) * 3600.0 / $PLATE_SCALE + $CHIP_DY*0.5 + set stars_M = (stars_mag - $4) + + set valid = (stars_X > 0) && (stars_X < $CHIP_DX) && (stars_Y > 0) && (stars_Y < $CHIP_DX) + subset stars_x = stars_X if valid + subset stars_y = stars_Y if valid + subset stars_m = stars_M if valid + subset stars_r = stars_ra if valid + subset stars_d = stars_dec if valid + subset stars_mr = stars_mag if valid + + # if we want to modify the observed magnitude based on the camera locations, add code here: + set stars_mi = stars_m + + # create bad detections / detections to ignore or avoid, save to a file? + + # create 3 types of things to ignore / avoid + # * badPhot (deviate the photometry) + # * badFlag (set a bad flag bit) + # * poorFlag (set a bad flag bit) + set badPhot = rnd(stars_r) < $DEFECT_FRAC + set badFlag = rnd(stars_r) < $DEFECT_FRAC + set poorFlag = rnd(stars_r) < $DEFECT_FRAC + + # create deviant photometry (1%) + # I'm setting the flag value 0x02 (EXTMODEL) if the object is deviated. This is a simple way to pass the info + # note that outlier rejection only occurs for Nmeas > 5 + set stars_mi = stars_mi + 0.2*badPhot + set stars_fl = zero(stars_mi) + 0x00000080*badFlag + 0x00000100*poorFlag + 0x00000002*badPhot + + # bad flags that we sould respect: + # PM_SOURCE_MODE_FAIL = 0x00000008, ///< Fit (non-linear) failed (non-converge, off-edge, run to zero) + # PM_SOURCE_MODE_SATSTAR = 0x00000080, ///< Source model peak is above saturation + # PM_SOURCE_MODE_BADPSF = 0x00000400, ///< Failed to get good estimate of object's PSF + # PM_SOURCE_MODE_DEFECT = 0x00000800, ///< Source is thought to be a defect + # PM_SOURCE_MODE_SATURATED = 0x00001000, ///< Source is thought to be saturated pixels (bleed trail) + # PM_SOURCE_MODE_CR_LIMIT = 0x00002000, ///< Source has crNsigma above limit + # PM_SOURCE_MODE_MOMENTS_FAILURE = 0x00008000, ///< could not measure the moments + # PM_SOURCE_MODE_SKY_FAILURE = 0x00010000, ///< could not measure the local sky + # PM_SOURCE_MODE_SKYVAR_FAILURE = 0x00020000, ///< could not measure the local sky variance + # PM_SOURCE_MODE_SIZE_SKIPPED = 0x10000000, ///< size could not be determined + + # poor flags that we should respect: + # PM_SOURCE_MODE_POOR = 0x00000010, ///< Fit succeeds, but low-SN, high-Chisq, or large (for PSF -- drop?) + # PM_SOURCE_MODE_PAIR = 0x00000020, ///< Source fitted with a double psf + # PM_SOURCE_MODE_BLEND = 0x00000100, ///< Source is a blend with other sources + # PM_SOURCE_MODE_BELOW_MOMENTS_SN = 0x00040000, ///< moments not measured due to low S/N + # PM_SOURCE_MODE_BLEND_FIT = 0x00400000, ///< source was fitted as a blend + # PM_SOURCE_MODE_ON_SPIKE = 0x20000000, ///< peak lands on diffraction spike + # PM_SOURCE_MODE_ON_GHOST = 0x40000000, ///< peak lands on ghost or glint + # PM_SOURCE_MODE_OFF_CHIP = 0x80000000, ///< peak lands off edge of chip + + write -f "%10.6f %10.6f %8.4f %6.1f %6.1f %8.4f 0x%x" $1 stars_r stars_d stars_mr stars_x stars_y stars_mi stars_fl +end + +# generate a set of stars with raw RA, DEC, MAG values to use for all of the output cmf files +macro mkstars + if ($0 != 2) + echo "USAGE: mkstars (Nstars)" + break + end + + local RA_RANGE DEC_RANGE + + # the images are oriented along N-S, E-W lines + + # size of region of interest in linear arcseconds + $RA_RANGE = 1.25 * $CHIP_DX * $NCHIP_X * $PLATE_SCALE + $DEC_RANGE = 1.25 * $CHIP_DY * $NCHIP_Y * $PLATE_SCALE + + create tmp 0 $1 + set stars_ra = $RA_CENTER + $RA_RANGE * (rnd(tmp) - 0.5) / 3600 / dcos ($DEC_CENTER) + set stars_dec = $DEC_CENTER + $DEC_RANGE * (rnd(tmp) - 0.5) / 3600 + set stars_mag = 10.0 + 10.0 * rnd(tmp) + set stars_dmag = 0.02 + 0.1 * rnd(tmp) +end Index: trunk/Ohana/src/relphot/test/relphot.tgroups.dvo =================================================================== --- trunk/Ohana/src/relphot/test/relphot.tgroups.dvo (revision 41390) +++ trunk/Ohana/src/relphot/test/relphot.tgroups.dvo (revision 41390) @@ -0,0 +1,446 @@ +# -*-sh-*- + +$TAP_BREAK = 0 +$TAP_VERBOSE = 1 +input tap.dvo +input ptolemy.dvo + +# set globals +if (not($?PLOT)) set PLOT = 0 + +$USE_CELL_OFFSETS = 1; # turn on non-zero cell offsets +$SET_UC_RESID = 1; # write a non-zero resid value to the ubercal table (should not have an impact) +$DO_2MASS = 1 +$DO_SYNTH = 1 +$DO_TYCHO = 1 +$NC_CLOUDS = 0.20 +$UC_CLOUDS = 0.02 + +## full test set: generate a set of test files, insert into a catdir, check raw catalog, run relphot, check results +macro go + if ($0 != 3) + echo "go (fileroot) (catdir)" + break + end + + local i fileroot rootdir catdir + + $fileroot = $1 + $catdir = $2 + + dirname $1 -var rootdir + mkdir $rootdir + + break -auto off + exec rm -rf $catdir + exec rm -rf $fileroot* + break -auto on + + # define a set of photometric (ubercal) and non-photometric exposure times + init + + # generate the common list of stars (coords & mags) for all exposures + mkstars 1000 + + catdir $catdir + + tapPLAN {((mjd_uc[]*4 + mjd_nc[]*4) + (mjd_uc[] + mjd_nc[]) + (mjd_uc[] + mjd_nc[])) / 3} + + # XXX first, run the analysis with only ubercal data (photometric data) + # we should be able to recover the nightly zero points fitting only tgroups + + # generate the basic images and check they were correctly ingested by dvo + for i 0 mjd_uc[] + sprintf filerootout "%s.%02d" $fileroot $i + mkexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] secz_uc[$i] $FILTER $filerootout $RA_CENTER $DEC_CENTER + ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] secz_uc[$i] $FILTER raw + end + + break + + # run relphot on the db and check that the images now match the expected values + tapEXEC relphot -tgroup-fit-airmass -tgroups tgroups.dat -images g -v -region $RA_MIN $RA_MAX $DEC_MIN $DEC_MAX -D CATDIR $catdir -D STAR_TOOFEW 1 -statmode WT_MEAN -cloud-limit 0.5 -update -nloop 12 + + break + + # reset the per-filter ZPT_OFF:$filter values + $ZPT_OFF:$FILTER = NAN + + for i 0 mjd_uc[] + ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $FILTER relphot_uc + end + for i 0 mjd_nc[] + ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $FILTER relphot_nc + end + tapDONE +end + +# set various global variables +macro init + + # generates ptolemy.rc, dvo.photcodes, $zpt_nominal:g,r $klam_nominal:g,r + mkptolemy.rc + + $RA_CENTER = 10.0 + $DEC_CENTER = 20.0 + $PLATE_SCALE = 0.25 + $NCHIP_X = 1 + $NCHIP_Y = 1 + $CHIP_DX = 1000 + $CHIP_DY = 1000 + $DEFECT_FRAC = 0.00 + + # we have two sets of images: ubercal (photometric) and not-ubercaled + + # ********* UBERCAL IMAGES *********** + + # ubercal nights each have a distinct zero point and slope: + vlist night_uc 55000.0 55010.0 55030.0 55060.0 + vlist dz_nt_uc +0.0000 +0.0300 -0.0200 -0.0400 +# vlist dk_nt_uc +0.0000 +0.0005 -0.0005 -0.0002 + vlist dk_nt_uc +0.0000 +0.0200 -0.0200 +0.0300 + + write tgroups.dat night_uc + + $EXPTIME = 10.0 + $FILTER = g + $ZPT_OFF:$FILTER = NAN + + delete -q mjd_uc seq_uc zpt_uc klam_uc secz_uc exptime_uc + + # generate 5 exposures for each night + for i 0 5 + set secz_tmp = rnd(night_uc) + 1.0; # airmass in range 1.0 - 2.0 + set klam_tmp = $klam_nominal:$FILTER + dk_nt_uc + set exp_tmp = zero(night_uc) + $EXPTIME + set zpt_tmp = $zpt_nominal:$FILTER + dz_nt_uc + 2.5*log(exp_tmp) + klam_tmp*(secz_tmp - 1.0) + + set seq_tmp = ramp(night_uc) + set mjd_tmp = night_uc + 0.33 + 0.01*$i; # exposures are 14 minutes apart starting at 10pm HST + + concat klam_tmp klam_uc + concat secz_tmp secz_uc + concat zpt_tmp zpt_uc + concat seq_tmp seq_uc + concat mjd_tmp mjd_uc + + concat exp_tmp exptime_uc + end + + # ********* NON-UBERCAL IMAGES *********** + + # non-ubercal nights each have a random zero points and fixed slope (klam_nominal) + # non-ubercal nights must not be the same set as the ubercal nights above + vlist night_nc 55005.0 55015.0 55035.0 55065.0 + + delete -q mjd_nc seq_nc zpt_nc klam_nc secz_nc + + # generate 5 exposures for each night + for i 0 5 + set secz_tmp = rnd(night_nc) + 1.0; # airmass in range 1.0 - 2.0 + set klam_tmp = zero(night_nc) + $klam_nominal:$FILTER + set exp_tmp = zero(night_nc) + $EXPTIME + set zpt_tmp = $zpt_nominal:$FILTER + 2.5*log(exp_tmp) + klam_tmp*(secz_tmp - 1.0) - $NC_CLOUDS*(0.1 + rnd(secz_tmp)) + # observed zero points for nc nights are 0.1 - 1.1 lower than the ubercal nights + + set seq_tmp = ramp(night_nc) + set mjd_tmp = night_nc + 0.33 + 0.01*$i; # exposures are 14 minutes apart starting at 10pm HST + + concat klam_tmp klam_nc + concat secz_tmp secz_nc + concat zpt_tmp zpt_nc + concat seq_tmp seq_nc + concat mjd_tmp mjd_nc + + concat exp_tmp exptime_nc + end +end + +macro try.relphot + if ($0 != 3) + echo "go (fileroot) (catdir)" + break + end + + local i fileroot rootdir catdir + + $fileroot = $1 + $catdir = $2 + + # run relphot on the db and check that the images now match the expected values + tapEXEC relphot -images g,r -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir -D STAR_TOOFEW 1 -statmode WT_MEAN -cloud-limit 0.5 -update -nloop 12 + + # reset the per-filter ZPT_OFF:$filter values + $ZPT_OFF:$FILTER = NAN + + for i 0 mjd_uc[] + # ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $FILTER relphot_uc + ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $FILTER relphot_uc + end +end + +macro ckexposure + if ($0 != 8) + echo "ckexposure (catdir) (mjd) (zpt) (exptime) (airmass) (filter) (mode)" + echo " mode == raw or corr" + break + end + + local CATDIR MJD_IMAGE ZPT_REAL ZPT_NOMINAL EXPTIME FILTER AIRMASS + + $CATDIR = $1 + $MJD_IMAGE = $2 + $ZPT_REAL = $3 + $EXPTIME = $4 + $AIRMASS = $5 + $FILTER = $6 + $FUNC_MODE = $7 + + # XXX need a function to extract the nominal zpt for a given filter / photcode from the db + $ZPT_NOMINAL = $zpt_nominal:$FILTER + $KLAM_NOMINAL = $klam_nominal:$FILTER + + $TIMEFORMAT = mjd + $TIMEREF = 0.0 + + catdir $CATDIR + + # assume we still have stars_ra, stars_dec, stars_mag in hand + if ($PLOT) + dev -n 0 + region $RA_CENTER $DEC_CENTER 0.2 + images + pmeasure -all -m 15 20 + end + + # RA_MIN, etc are defined by mkstars + skyregion $RA_MIN $RA_MAX $DEC_MIN $DEC_MAX + + ## calculate the average zero point offset + if (isnan($ZPT_OFF:$FILTER)) + avextract ra dec $FILTER + match2d -closest ra dec stars_ra stars_dec 0.001 -index1 index1 -index2 index2 + reindex stars_tru_m = stars_mag using index1 + set dm_tru = $FILTER - stars_tru_m + vstat -q dm_tru + $ZPT_OFF:$FILTER = $MEDIAN + end + + mextract ra dec mag xccd yccd dbflags photflags time mag:sys mag:rel mag:inst mag:ave where (abs(time - $MJD_IMAGE) < 0.0001) + + match2d -closest ra dec stars_ra stars_dec 0.001 -index1 index1 -index2 index2 + + reindex stars_ra_m = stars_ra using index1 + reindex stars_dec_m = stars_dec using index1 + reindex stars_mag_m = stars_mag using index1 + set dr = 3600*(ra - stars_ra_m) + set dd = 3600*(dec - stars_dec_m) + set dm = mag - stars_mag_m + + set good_photom = (not(photflags & 0x02)) + foreach field dm xccd yccd + subset $field\_good = $field if good_photom + end + + if ($PLOT) + dev -n 1 + lim mag dm; clear; box; plot mag dm; plot -c red -pt 7 mag dm where (photflags & 0x00000002) + dev -n 2 + lim xccd dm; clear; box; plot xccd dm; plot -c red -pt 7 xccd dm where (photflags & 0x00000002) + dev -n 3 + lim ra dm; clear; box; plot ra dm; plot -c red -pt 7 ra dm where (photflags & 0x00000002) + echo "red points are known bad photometry" + end + + $ZPT_REAL_NORM = $ZPT_REAL - 2.5*log($EXPTIME) - $KLAM_NOMINAL*($AIRMASS - 1.0) + + if ("$FUNC_MODE" == "raw") + set dm_adjust = dm_good - $ZPT_NOMINAL + $ZPT_REAL_NORM + + vstat -q dm_adjust + sprintf line "addstar raw ZP %8.2f $FILTER %6.3f %6.3f %6.3f (%7.4f): %7.4f %7.4f : %6.3f" $MJD_IMAGE $ZPT_REAL $ZPT_REAL_NORM $ZPT_NOMINAL {$ZPT_REAL_NORM - $ZPT_NOMINAL} $MEAN $SIGMA dm_adjust[] + tapOK {abs($MEAN) < 0.005} "RAW : $line" + return + end + + if (("$FUNC_MODE" == "relphot_uc") || ("$FUNC_MODE" == "relphot_nc")) + + ## relphot does not set an absolute zero point, so we need to find that and subtract + set dm_adjust = dm_good - $ZPT_OFF:$FILTER + + vstat -q dm_adjust + sprintf line "relphot ZP %8.2f $FILTER %6.3f %6.3f %6.3f (%7.4f) : %7.4f %7.4f : %3d" $MJD_IMAGE $ZPT_REAL $ZPT_REAL_NORM $ZPT_NOMINAL {$ZPT_REAL_NORM - $ZPT_NOMINAL} $MEAN $SIGMA dm_adjust[] + + if ("$FUNC_MODE" == "relphot_uc") + tapOK {abs($MEAN) < 0.001} "MEAN : $line" + tapOK {abs($SIGMA) < 0.001} "SIGMA: $line" + end + if ("$FUNC_MODE" == "relphot_nc") + # XXX NOTE relaxed condition : bad photometry points are partly clipped but still cause a bias + tapOK {abs($MEAN) < 0.020} "MEAN : $line" + tapOK {abs($SIGMA) < 0.020} "SIGMA: $line" + end + + if ($PLOT) + cursor 1 + if ("$KEY" == "Escape") break + end + return + end + + echo "ERROR: unknown mode for ckexposure : $FUNC_MODE" +end + +macro mkexposure + if ($0 != 10) + echo "mkexposure (catdir) (mjd) (zpt) (exptime) (airmass) (filter) (fileroot) (ra) (dec)" + break + end + + local ix iy date time datetime ra dec ROOT RAo DECo ZPT MJD FILTER EXPTIME AIRMASS CATDIR cmd + + $CATDIR = $1 + $MJD = $2 + $ZPT = $3 + $EXPTIME = $4 + $AIRMASS = $5 + $FILTER = $6 + $ROOT = $7 + $RAo = $8 + $DECo = $9 + + $TIMEFORMAT = mjd + $TIMEREF = 0.0 + + # XX ctimes -ref $MJD -var datetime + # XX substr $datetime 0 10 date + # XX substr $datetime 11 8 time + + for ix 0 $NCHIP_X + for iy 0 $NCHIP_Y + $dx = $CHIP_DX * $PLATE_SCALE * ($ix - 0.5*$NCHIP_X + 0.5) + $dy = $CHIP_DY * $PLATE_SCALE * ($iy - 0.5*$NCHIP_Y + 0.5) + $ra = $RAo - $dx / 3600.0 / dcos($DECo) + $dec = $DECo - $dy / 3600.0 + # echo $ra $dec $dx $dy + mkinput $ROOT.$ix.$iy.in.txt $ra $dec $ZPT + + # ra,dec is the center of this chip + local options + $options = -mjd $MJD + $options = $options -radec $ra $dec + + #$options = $options -type PS1_V2 + #$options = $options -type PS1_V3 + #$options = $options -type PS1_V4 + $options = $options -type PS1_V5 + + $options = $options -coords + $options = $options -photcode GPC1.$FILTER.XY$ix\$iy + $options = $options -no-noise + $options = $options -flags + $options = $options -bad-psfqf-frac $DEFECT_FRAC + $options = $options -size $CHIP_DX $CHIP_DY + $options = $options -crpix {0.5*$CHIP_DX} {0.5*$CHIP_DY} + $options = $options -airmass $AIRMASS + $options = $options -exptime $EXPTIME + + tapEXEC mkcmf $ROOT.$ix.$iy.in.txt $ROOT.$ix.$iy.cmf $options + + # the fake images have inconsistent ra,dec and airmass,sidtime values + tapEXEC addstar $ROOT.$ix.$iy.cmf -D CATDIR $CATDIR -D CAMERA gpc1 -quick-airmass + end + end +end + +# make a simple input file for mkcmf +macro mkinput + if ($0 != 5) + echo "mkinput (filename) (ra) (dec) (zpt)" + break + end + + # chip coordinate of the stars + # ra,dec is at the chip center, which corresponds to pixel (CHIP_DX,CHIP_DY)/2 + set stars_X = (stars_ra - $2) * 3600.0 * dcos($3) / $PLATE_SCALE + $CHIP_DX*0.5 + set stars_Y = (stars_dec - $3) * 3600.0 / $PLATE_SCALE + $CHIP_DY*0.5 + set stars_M = (stars_mag - $4) + + set valid = (stars_X > 0) && (stars_X < $CHIP_DX) && (stars_Y > 0) && (stars_Y < $CHIP_DX) + subset stars_x = stars_X if valid + subset stars_y = stars_Y if valid + subset stars_m = stars_M if valid + subset stars_r = stars_ra if valid + subset stars_d = stars_dec if valid + subset stars_mr = stars_mag if valid + + # if we want to modify the observed magnitude based on the camera locations, add code here: + set stars_mi = stars_m + + # create bad detections / detections to ignore or avoid, save to a file? + + # create 3 types of things to ignore / avoid + # * badPhot (deviate the photometry) + # * badFlag (set a bad flag bit) + # * poorFlag (set a bad flag bit) + set badPhot = rnd(stars_r) < $DEFECT_FRAC + set badFlag = rnd(stars_r) < $DEFECT_FRAC + set poorFlag = rnd(stars_r) < $DEFECT_FRAC + + # create deviant photometry (1%) + # I'm setting the flag value 0x02 (EXTMODEL) if the object is deviated. This is a simple way to pass the info + # note that outlier rejection only occurs for Nmeas > 5 + set stars_mi = stars_mi + 0.2*badPhot + set stars_fl = zero(stars_mi) + 0x00000080*badFlag + 0x00000100*poorFlag + 0x00000002*badPhot + + # bad flags that we sould respect: + # PM_SOURCE_MODE_FAIL = 0x00000008, ///< Fit (non-linear) failed (non-converge, off-edge, run to zero) + # PM_SOURCE_MODE_SATSTAR = 0x00000080, ///< Source model peak is above saturation + # PM_SOURCE_MODE_BADPSF = 0x00000400, ///< Failed to get good estimate of object's PSF + # PM_SOURCE_MODE_DEFECT = 0x00000800, ///< Source is thought to be a defect + # PM_SOURCE_MODE_SATURATED = 0x00001000, ///< Source is thought to be saturated pixels (bleed trail) + # PM_SOURCE_MODE_CR_LIMIT = 0x00002000, ///< Source has crNsigma above limit + # PM_SOURCE_MODE_MOMENTS_FAILURE = 0x00008000, ///< could not measure the moments + # PM_SOURCE_MODE_SKY_FAILURE = 0x00010000, ///< could not measure the local sky + # PM_SOURCE_MODE_SKYVAR_FAILURE = 0x00020000, ///< could not measure the local sky variance + # PM_SOURCE_MODE_SIZE_SKIPPED = 0x10000000, ///< size could not be determined + + # poor flags that we should respect: + # PM_SOURCE_MODE_POOR = 0x00000010, ///< Fit succeeds, but low-SN, high-Chisq, or large (for PSF -- drop?) + # PM_SOURCE_MODE_PAIR = 0x00000020, ///< Source fitted with a double psf + # PM_SOURCE_MODE_BLEND = 0x00000100, ///< Source is a blend with other sources + # PM_SOURCE_MODE_BELOW_MOMENTS_SN = 0x00040000, ///< moments not measured due to low S/N + # PM_SOURCE_MODE_BLEND_FIT = 0x00400000, ///< source was fitted as a blend + # PM_SOURCE_MODE_ON_SPIKE = 0x20000000, ///< peak lands on diffraction spike + # PM_SOURCE_MODE_ON_GHOST = 0x40000000, ///< peak lands on ghost or glint + # PM_SOURCE_MODE_OFF_CHIP = 0x80000000, ///< peak lands off edge of chip + + write -f "%10.6f %10.6f %8.4f %6.1f %6.1f %8.4f 0x%x" $1 stars_r stars_d stars_mr stars_x stars_y stars_mi stars_fl +end + +# generate a set of stars with raw RA, DEC, MAG values to use for all of the output cmf files +macro mkstars + if ($0 != 2) + echo "USAGE: mkstars (Nstars)" + break + end + + local RA_RANGE DEC_RANGE + + # the images are oriented along N-S, E-W lines + + # size of region of interest in linear degrees + $RA_RANGE = 1.25 * $CHIP_DX * $NCHIP_X * $PLATE_SCALE / 3600 + $DEC_RANGE = 1.25 * $CHIP_DY * $NCHIP_Y * $PLATE_SCALE / 3600 + + $RA_MIN = $RA_CENTER - 0.5*$RA_RANGE/dcos($DEC_CENTER) + $RA_MAX = $RA_CENTER + 0.5*$RA_RANGE/dcos($DEC_CENTER) + $DEC_MIN = $DEC_CENTER - 0.5*$DEC_RANGE + $DEC_MAX = $DEC_CENTER + 0.5*$DEC_RANGE + + create tmp 0 $1 + set stars_ra = $RA_CENTER + $RA_RANGE * (rnd(tmp) - 0.5) / dcos ($DEC_CENTER) + set stars_dec = $DEC_CENTER + $DEC_RANGE * (rnd(tmp) - 0.5) + set stars_mag = 10.0 + 10.0 * rnd(tmp) + set stars_dmag = 0.02 + 0.1 * rnd(tmp) +end Index: trunk/Ohana/src/relphot/test/tap.dvo =================================================================== --- trunk/Ohana/src/relphot/test/tap.dvo (revision 41390) +++ trunk/Ohana/src/relphot/test/tap.dvo (revision 41390) @@ -0,0 +1,107 @@ +# -*-sh-*- + +if (not($?TAP_BREAK)) set TAP_BREAK = 0 +if (not($?TAP_VERBOSE)) set TAP_VERBOSE = 0 + +macro tapEXEC + if ($0 < 3) + echo "USAGE: tapEXEC (command)" + break + end + + local i cmd result + + $cmd = "" + for i 1 $0 + $cmd = $cmd $$i + end + + if ($TAP_VERBOSE >= 2) echo $cmd + + break -auto off + exec $cmd >& tap.log + $result = $STATUS + break -auto on + + if (not($result)) + exec cat tap.log + echo "not ok : $cmd" + $TAP_NFAIL ++ + $TAP_LAST = 0 + if ($TAP_BREAK) + break + end + else + if ($TAP_VERBOSE >= 1) echo "ok : $cmd" + $TAP_LAST = 1 + end + $TAP_NDONE ++ +end + +macro tapOK + if ($0 != 3) + echo "USAGE: tapOK (condition) (message)" + break + end + + if ($1) + if ($TAP_VERBOSE >= 1) + echo "ok : $2" + end + $TAP_LAST = 1 + else + echo "not ok : $2" + $TAP_NFAIL ++ + $TAP_LAST = 0 + if ($TAP_BREAK) + break + end + end + $TAP_NDONE ++ +end + +macro tapPLAN + if ($0 != 2) + echo "USAGE: tapPLAN (Ntests)" + break + end + + $TAP_NTEST = $1 + $TAP_NFAIL = 0 + $TAP_NSKIP = 0 + $TAP_NDONE = 0 + if (not($?TAP_BREAK)) set TAP_BREAK = 0 +end + +macro tapSKIP + if ($0 != 2) + echo "USAGE: tapSKIP (nskip)" + break + end + + $TAP_NSKIP += $1 + $TAP_NDONE += $1 +end + +macro tapDONE + if ($0 != 1) + echo "USAGE: tapDONE" + break + end + + if ($TAP_NDONE != $TAP_NTEST) + echo "planned tests ($TAP_NTEST) not equal to done tests ($TAP_NDONE)" + end + + if ($TAP_NFAIL) + echo "failed $TAP_NFAIL of $TAP_NDONE" + end + + if ($TAP_NSKIP) + echo "skipped $TAP_NSKIP of $TAP_NDONE" + end + + if ($TAP_NFAIL + $TAP_NSKIP == 0) + echo "passed $TAP_NDONE tests" + end +end