Index: trunk/Ohana/src/dvomerge/src/dvo_image_merge_dbs.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvo_image_merge_dbs.c (revision 31425) +++ trunk/Ohana/src/dvomerge/src/dvo_image_merge_dbs.c (revision 31450) @@ -202,5 +202,5 @@ if (newID == 0) { fprintf (stderr, "cannot find image ID "OFF_T_FMT"\n", oldID); - exit (2); + // exit (2); } catalog[0].measure[i].imageID = newID; Index: trunk/Ohana/src/dvomerge/src/dvoverify.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvoverify.c (revision 31425) +++ trunk/Ohana/src/dvomerge/src/dvoverify.c (revision 31450) @@ -148,5 +148,5 @@ int VerifyTableFile (char *filename) { - int status; + int status, Next; off_t Nbytes; Header header; @@ -205,4 +205,5 @@ // scan all extentions Nbytes = 0; + Next = -1; if (DEBUG) fprintf (stderr, "sizes: ("OFF_T_FMT" vs "OFF_T_FMT")\n", Nbytes, fileStats.st_size); while (Nbytes < fileStats.st_size) { @@ -210,5 +211,9 @@ // Check on the PHU if (!gfits_fread_header (file, &header)) { - fprintf (stderr, "unable to read PHU header for %s\n", filename); + if (Next == -1) { + fprintf (stderr, "unable to read PHU header for %s\n", filename); + } else { + fprintf (stderr, "unable to read header for %s, extension %d (or file has excess bytes)\n", filename, Next); + } fclose(file); return (FALSE); @@ -243,4 +248,5 @@ } } + Next ++; } if (DEBUG) fprintf (stderr, "file is good: %s\n", filename); Index: trunk/Ohana/src/libdvo/include/dvo.h =================================================================== --- trunk/Ohana/src/libdvo/include/dvo.h (revision 31425) +++ trunk/Ohana/src/libdvo/include/dvo.h (revision 31450) @@ -20,7 +20,7 @@ DVO_FORMAT_PS1_DEV_2, DVO_FORMAT_PS1_DEV_3, + DVO_FORMAT_PS1_REF, DVO_FORMAT_PS1_V1, DVO_FORMAT_PS1_V2, - DVO_FORMAT_PS1_REF } DVOTableFormat; @@ -210,4 +210,31 @@ } PhotCodeData; +// a reduced-subset structure for relphot +typedef struct { + double R; + double D; + unsigned short Nmeasure; + int measureOffset; + uint32_t flags; +} AverageTiny; + +// a reduced-subset structure for relphot +typedef struct { + float dR; + float dD; + float M; + float Mcal; + float dM; + float airmass; + float Xccd; + float Yccd; + float dt; + int t; + unsigned int averef; + unsigned int imageID; + unsigned int dbFlags; + unsigned short photcode; +} MeasureTiny; + /* a catalog contains this data */ typedef struct Catalog { @@ -226,4 +253,8 @@ off_t Naves_disk, Nmeas_disk, Nmiss_disk, Nsecf_disk; /* current number of each component on disk */ off_t Naves_off, Nmeas_off, Nmiss_off, Nsecf_off; /* index of first loaded data value */ + + // note that we use these for the full-sky relphot analysis + AverageTiny *averageT; + MeasureTiny *measureT; /* the Nsecf_* values above are number of table rows (eg, Naverage*Nsecfilt) */ @@ -330,4 +361,19 @@ int PhotColor (Average *average, SecFilt *secfilt, Measure *measure, int c1, int c2, double *color); +float PhotInstTiny (MeasureTiny *measure); +float PhotCatTiny (MeasureTiny *measure); +float PhotAperTiny (MeasureTiny *measure); +float PhotSysTiny (MeasureTiny *measure, AverageTiny *average, SecFilt *secfilt); +float PhotRelTiny (MeasureTiny *measure, AverageTiny *average, SecFilt *secfilt); +float PhotCalTiny (MeasureTiny *thisone, AverageTiny *average, SecFilt *secfilt, MeasureTiny *measure, PhotCode *code); +float PhotAveTiny (PhotCode *code, AverageTiny *average, SecFilt *secfilt); +float PhotRefTiny (PhotCode *code, AverageTiny *average, SecFilt *secfilt, MeasureTiny *measure); +float PhotXmTiny (PhotCode *code, AverageTiny *average, SecFilt *secfilt); +float PhotdMTiny (PhotCode *code, AverageTiny *average, SecFilt *secfilt); + +float PhotColorForCodeTiny (AverageTiny *average, SecFilt *secfilt, MeasureTiny *measure, PhotCode *code); +int PhotColorTiny (AverageTiny *average, SecFilt *secfilt, MeasureTiny *measure, int c1, int c2, double *color); + + PhotCodeData *GetPhotcodeTable (void); void SetPhotcodeTable (PhotCodeData *); Index: trunk/Ohana/src/libdvo/src/dvo_catalog_split.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_catalog_split.c (revision 31425) +++ trunk/Ohana/src/libdvo/src/dvo_catalog_split.c (revision 31450) @@ -448,5 +448,6 @@ catalog[0].measure = FtableToMeasure (&ftable, &Nmeasure, &catalog[0].catformat); if (Nmeasure != Nrows) { - fprintf (stderr, "Warning: mismatch between Nmeasure in PHU and Table headers ("OFF_T_FMT" vs "OFF_T_FMT")\n", Nmeasure, Nrows); + // XXX this condition denotes the eof has been reached; not an error or a warning + // fprintf (stderr, "Warning: mismatch between Nmeasure in PHU and Table headers ("OFF_T_FMT" vs "OFF_T_FMT")\n", Nmeasure, Nrows); } gfits_free_header (&header); Index: trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c (revision 31425) +++ trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c (revision 31450) @@ -536,4 +536,324 @@ return (TRUE); } + +/******** alternate photometry conversion functions using MeasureTiny and AverageTiny *********/ +float PhotInstTiny (MeasureTiny *measure) { + + int Np; + float M; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + M = measure[0].M; + return (M); + } + + M = measure[0].M - measure[0].dt - ZERO_POINT; + + return (M); + +} + +float PhotCatTiny (MeasureTiny *measure) { + + int Np; + float Mcat; + PhotCode *code; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Mcat = measure[0].M; + return (Mcat); + } + code = &photcodes[0].code[Np]; + Mcat = measure[0].M - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C; + + return (Mcat); +} + +# if (0) +float PhotAperTiny (MeasureTiny *measure) { + + int Np; + float Mcat; + PhotCode *code; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Mcat = measure[0].Map; + return (Mcat); + } + code = &photcodes[0].code[Np]; + Mcat = measure[0].Map - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C; + + return (Mcat); +} +# endif + +float PhotSysTiny (MeasureTiny *measure, AverageTiny *average, SecFilt *secfilt) { + + int i, Np; + float Mcat, Mcol, Msys, mc, Mc; + PhotCode *code; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Msys = measure[0].M; + return (Msys); + } + code = &photcodes[0].code[Np]; + Mcat = measure[0].M - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C; + + /* for DEP, color must be made of PRI/SEC */ + mc = PhotColorForCodeTiny (average, secfilt, NULL, code); + if (isnan(mc)) return (Mcat); + mc = mc - SCALE*code[0].dX; + + Mc = mc; + Mcol = 0; + for (i = 0; i < code[0].Nc; i++) { + Mcol += code[0].X[i]*Mc; + Mc *= mc; + } + Msys = Mcat + Mcol; + return (Msys); +} + +float PhotRelTiny (MeasureTiny *measure, AverageTiny *average, SecFilt *secfilt) { + + int i, Np; + float Mcat, Mcol, Mrel, mc, Mc; + PhotCode *code; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Mcat = measure[0].M; + return (Mcat); + } + code = &photcodes[0].code[Np]; + Mrel = measure[0].M - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C - measure[0].Mcal; + + /* for DEP, color must be made of PRI/SEC */ + mc = PhotColorForCodeTiny (average, secfilt, NULL, code); + if (isnan(mc)) return (Mrel); + mc = mc - SCALE*code[0].dX; + + Mc = mc; + Mcol = 0; + for (i = 0; i < code[0].Nc; i++) { + Mcol += code[0].X[i]*Mc; + Mc *= mc; /* the 0.001 is needed for higher order terms to keep the units mag = mag^n */ + } + Mrel += Mcol; + return (Mrel); +} + +/* return calibrated magnitude from measure for given photcode */ +float PhotCalTiny (MeasureTiny *thisone, AverageTiny *average, SecFilt *secfilt, MeasureTiny *measure, PhotCode *code) { + + int i, Np; + float Mcal, Mrel, Mcol, mc, Mc; + + if (code == NULL) return NAN; + + /* code must be the matching PRI/SEC code for this measurement or an equivalent ALT */ + Np = photcodes[0].hashcode[thisone[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Mrel = thisone[0].M; + return (Mrel); + } + if (code[0].code != photcodes[0].code[Np].equiv) return (NAN); + + Mcal = PhotRelTiny (thisone, average, secfilt) + SCALE*code[0].C; + + mc = PhotColorForCodeTiny (average, secfilt, measure, code); + if (isnan(mc)) return (Mcal); + mc = mc - SCALE*code[0].dX; + + Mc = mc; + Mcol = 0; + for (i = 0; i < code[0].Nc; i++) { + Mcol += code[0].X[i]*Mc; + Mc *= mc; + } + Mcal += Mcol; + return (Mcal); +} + +/* color term may not use DEP magnitude */ +float PhotColorForCodeTiny (AverageTiny *average, SecFilt *secfilt, MeasureTiny *measure, PhotCode *code) { + + int i, Ns1, Ns2, Ns; + float m1, m2, mc; + PhotCode *color; + + m1 = m2 = NAN; + + if (measure == NULL) { + Ns1 = photcodes[0].hashNsec[code[0].c1]; + Ns2 = photcodes[0].hashNsec[code[0].c2]; + + m1 = (Ns1 == -1) ? NAN : secfilt[Ns1].M; + m2 = (Ns2 == -1) ? NAN : secfilt[Ns2].M; + mc = (isnan(m1) || isnan(m2)) ? NAN : (m1 - m2); + return (mc); + } + + /* find magnitude matching first color term */ + color = GetPhotcodebyCode (code[0].c1); + if (color == NULL) return (NAN); + if (color[0].type == PHOT_REF) { + for (i = 0; (i < average[0].Nmeasure) && (isnan(m1)); i++) { + if (measure[i].photcode == color[0].code) { + m1 = measure[i].M; + } + } + } else { + Ns = photcodes[0].hashNsec[color[0].code]; + m1 = (Ns == -1) ? NAN : secfilt[Ns].M; + } + + /* find magnitude matching second color term */ + color = GetPhotcodebyCode (code[0].c2); + if (color == NULL) return (NAN); + if (color[0].type == PHOT_REF) { + for (i = 0; (i < average[0].Nmeasure) && (isnan(m2)); i++) { + if (measure[i].photcode == color[0].code) { + m2 = measure[i].M; + } + } + } else { + Ns = photcodes[0].hashNsec[color[0].code]; + m2 = (Ns == -1) ? NAN : secfilt[Ns].M; + } + mc = (isnan(m1) || isnan(m2)) ? NAN : (m1 - m2); + return (mc); +} + +/* return calibrated magnitude from average/secfilt for given photcode */ +float PhotRefTiny (PhotCode *code, AverageTiny *average, SecFilt *secfilt, MeasureTiny *measure) { + + int i, Ns; + float Mave, Mref, Mcol, mc; + double Mc; + + if (code == NULL) return NAN; + + Ns = photcodes[0].hashNsec[code[0].code]; + Mave = (Ns == -1) ? NAN : secfilt[Ns].M; + Mref = Mave + SCALE*code[0].C; + + mc = PhotColorForCodeTiny (average, secfilt, measure, code); + if (isnan(mc)) return (Mref); + mc = mc - SCALE*code[0].dX; + + Mc = mc; + Mcol = 0; + for (i = 0; i < code[0].Nc; i++) { + Mcol += code[0].X[i]*Mc; + Mc *= mc; /* the 0.001 is needed for higher order terms to keep the units mag = mag^n */ + } + Mref += Mcol; + return (Mref); +} + +/***/ +float PhotAveTiny (PhotCode *code, AverageTiny *average, SecFilt *secfilt) { + + int Ns; + float Mave; + + if (code == NULL) return NAN; + + Ns = photcodes[0].hashNsec[code[0].code]; + Mave = (Ns == -1) ? NAN : secfilt[Ns].M; + return (Mave); +} + +float PhotdMTiny (PhotCode *code, AverageTiny *average, SecFilt *secfilt) { + + int Ns; + float dM; + + if (code == NULL) return NAN; + + Ns = photcodes[0].hashNsec[code[0].code]; + dM = (Ns == -1) ? NAN : secfilt[Ns].dM; + return (dM); +} + +// XXX return NAN or NAN_S_SHORT? (secfilt->Xm is short) +float PhotXmTiny (PhotCode *code, AverageTiny *average, SecFilt *secfilt) { + + int Ns; + short Mi; + float Xm; + + if (code == NULL) return NAN; + + Ns = photcodes[0].hashNsec[code[0].code]; + Mi = (Ns == -1) ? NAN : secfilt[Ns].Xm; + Xm = (isnan(Mi)) ? -1.0 : pow (10.0, 0.01*Mi); + return (Xm); +} + +/* given a photcode pair c1 & c2, return the color of this star (NaN if not found) */ +int PhotColorTiny (AverageTiny *average, SecFilt *secfilt, MeasureTiny *measure, int c1, int c2, double *color) { + + int i, Ns; + double M1, M2, dM; + PhotCode *code; + + code = GetPhotcodebyCode (c1); + if (code == NULL) return (FALSE); + if (code[0].type == PHOT_REF) { + for (i = 0; i < average[0].Nmeasure; i++) { + if (measure[i].photcode == c1) { + M1 = measure[i].M; + goto filter1; + } + } + return (FALSE); + } else { + Ns = photcodes[0].hashNsec[code[0].code]; + M1 = (Ns == -1) ? NAN : secfilt[Ns].M; + } + +filter1: + code = GetPhotcodebyCode (c2); + if (code == NULL) return (FALSE); + if (code[0].type == PHOT_REF) { + for (i = 0; i < average[0].Nmeasure; i++) { + if (measure[i].photcode == c2) { + M2 = measure[i].M; + goto filter2; + } + } + return (FALSE); + } else { + Ns = photcodes[0].hashNsec[code[0].code]; + M2 = (Ns == -1) ? NAN : secfilt[Ns].M; + } + +filter2: + + dM = M1 - M2; + *color = dM; + + return (TRUE); +} +/***********************************************/ // Create a map between the secfilt values from one photcode table to another Index: trunk/Ohana/src/opihi/cmd.data/Makefile =================================================================== --- trunk/Ohana/src/opihi/cmd.data/Makefile (revision 31425) +++ trunk/Ohana/src/opihi/cmd.data/Makefile (revision 31450) @@ -141,4 +141,5 @@ $(SRC)/vpop.$(ARCH).o \ $(SRC)/vroll.$(ARCH).o \ +$(SRC)/vshift.$(ARCH).o \ $(SRC)/vsmooth.$(ARCH).o \ $(SRC)/vstats.$(ARCH).o \ Index: trunk/Ohana/src/opihi/cmd.data/gridify.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/gridify.c (revision 31425) +++ trunk/Ohana/src/opihi/cmd.data/gridify.c (revision 31450) @@ -4,5 +4,5 @@ int i, Nx, Ny, Xb, Yb, Normalize, N; - float Xmin, Xmax, dX, Ymin, Ymax, dY; + float Xmin, Xmax, dX, Ymin, Ymax, dY, initValue; float *buf, *val; int *Nval; @@ -15,4 +15,11 @@ remove_argument (N, &argc, argv); Normalize = FALSE; + } + + initValue = 0.0; + if ((N = get_argument (argc, argv, "-init-value"))) { + remove_argument (N, &argc, argv); + initValue = atof(argv[N]); + remove_argument (N, &argc, argv); } @@ -69,7 +76,7 @@ buf = (float *) bf[0].matrix.buffer; for (i = 0; i < Nx*Ny; i++) { + buf[i] = initValue; if (Normalize) { if (Nval[i] == 0) { - buf[i] = 0; continue; } Index: trunk/Ohana/src/opihi/cmd.data/init.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/init.c (revision 31425) +++ trunk/Ohana/src/opihi/cmd.data/init.c (revision 31450) @@ -130,4 +130,5 @@ int vstats PROTO((int, char **)); int vroll PROTO((int, char **)); +int vshift PROTO((int, char **)); int vpop PROTO((int, char **)); int vsmooth PROTO((int, char **)); @@ -274,5 +275,6 @@ {1, "vmaxwell", vmaxwell, "fit a Maxwellian to a vector"}, {1, "vpop", vpop, "remove first element of a vector"}, - {1, "vroll", vroll, "roll vector elements"}, + {1, "vroll", vroll, "roll vector elements by 1 entry"}, + {1, "vshift", vshift, "shift vector elements by arbitrary amount"}, {1, "vsmooth", vsmooth, "Gaussian smooth of a vector"}, {1, "vstats", vstats, "statistics on a vector"}, Index: trunk/Ohana/src/opihi/cmd.data/vshift.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/vshift.c (revision 31450) +++ trunk/Ohana/src/opihi/cmd.data/vshift.c (revision 31450) @@ -0,0 +1,73 @@ +# include "data.h" + +int vshift (int argc, char **argv) { + + int isPos, N, Npix, Nshift, delta, ROLL; + Vector *ivec, *ovec; + + ROLL = FALSE; + if ((N = get_argument (argc, argv, "-roll"))) { + remove_argument (N, &argc, argv); + ROLL = TRUE; + } + + if (argc != 4) { + gprint (GP_ERR, "USAGE: vshift (input) (output) (delta) [-roll]\n"); + gprint (GP_ERR, " shift vector by (delta) elements\n"); + gprint (GP_ERR, " a positive value move element (i) to (i+delta)\n"); + gprint (GP_ERR, " -roll : move dropped values to the other size (no elements are lost)\n"); + return (FALSE); + } + + if ((ivec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((ovec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) return (FALSE); + delta = atoi(argv[3]); + isPos = delta > 0; + delta = abs(delta); + + ResetVector (ovec, ivec->type, ivec->Nelements); + + Npix = ivec[0].Nelements; + if (delta > Npix) { + if (!ROLL) { + if (ivec[0].type == OPIHI_FLT) { + memset (ovec[0].elements.Flt, 0, Npix*sizeof(opihi_flt)); + } else { + memset (ovec[0].elements.Flt, 0, Npix*sizeof(opihi_flt)); + } + return TRUE; + } + delta = delta % Npix; + } + + if (ivec[0].type == OPIHI_FLT) { + memset (ovec[0].elements.Flt, 0, Npix*sizeof(opihi_flt)); + Nshift = Npix - delta; + if (isPos) { + memcpy (&ovec[0].elements.Flt[delta], &ivec[0].elements.Flt[0], Nshift*sizeof(opihi_flt)); + if (ROLL) { + memcpy (&ovec[0].elements.Flt[0], &ivec[0].elements.Flt[Nshift], delta*sizeof(opihi_flt)); + } + } else { + memcpy (&ovec[0].elements.Flt[0], &ivec[0].elements.Flt[delta], Nshift*sizeof(opihi_flt)); + if (ROLL) { + memcpy (&ovec[0].elements.Flt[Nshift], &ivec[0].elements.Flt[0], delta*sizeof(opihi_flt)); + } + } + } else { + memset (ovec[0].elements.Int, 0, Npix*sizeof(opihi_int)); + Nshift = Npix - delta; + if (isPos) { + memcpy (&ovec[0].elements.Int[delta], &ivec[0].elements.Int[0], Nshift*sizeof(opihi_int)); + if (ROLL) { + memcpy (&ovec[0].elements.Int[0], &ivec[0].elements.Int[Nshift], delta*sizeof(opihi_int)); + } + } else { + memcpy (&ovec[0].elements.Int[0], &ivec[0].elements.Int[delta], Nshift*sizeof(opihi_int)); + if (ROLL) { + memcpy (&ovec[0].elements.Int[Nshift], &ivec[0].elements.Int[0], delta*sizeof(opihi_int)); + } + } + } + return (TRUE); +} Index: trunk/Ohana/src/opihi/dvo/Makefile =================================================================== --- trunk/Ohana/src/opihi/dvo/Makefile (revision 31425) +++ trunk/Ohana/src/opihi/dvo/Makefile (revision 31450) @@ -80,4 +80,5 @@ $(SRC)/mextract.$(ARCH).o \ $(SRC)/mmextract.$(ARCH).o \ +$(SRC)/objectcoverage.$(ARCH).o \ $(SRC)/photcodes.$(ARCH).o \ $(SRC)/pmeasure.$(ARCH).o \ Index: trunk/Ohana/src/opihi/dvo/avextract.c =================================================================== --- trunk/Ohana/src/opihi/dvo/avextract.c (revision 31425) +++ trunk/Ohana/src/opihi/dvo/avextract.c (revision 31450) @@ -5,5 +5,5 @@ off_t i, j, n, m; int N, Npts, NPTS, last, next, state, Nfields, Nreturn, Ncstack, Nstack; - int Nsecfilt, mode, VERBOSE; + int Nsecfilt, mode, VERBOSE, needMeasures; char **cstack, name[1024]; void *Signal; @@ -101,4 +101,13 @@ } + // check the requested fields : are all average/secfilt entries, or do we need measures? + needMeasures = FALSE; + for (i = 0; !needMeasures && (i < Nfields); i++) { + if (fields[i].magMode == MAG_NONE) continue; + if (fields[i].photcode == NULL) continue; // assert this? + if (fields[i].photcode[0].type == PHOT_REF) needMeasures = TRUE; + if (fields[i].photcode[0].type == PHOT_DEP) needMeasures = TRUE; + } + // grab data from all selected sky regions Signal = signal (SIGINT, handle_interrupt); @@ -107,5 +116,8 @@ /* lock, load, unlock catalog */ catalog.filename = skylist[0].filename[i]; - catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_SECF; + catalog.catflags = LOAD_AVES | LOAD_SECF; + if (needMeasures) { + catalog.catflags |= LOAD_MEAS; + } catalog.Nsecfilt = 0; Index: trunk/Ohana/src/opihi/dvo/init.c =================================================================== --- trunk/Ohana/src/opihi/dvo/init.c (revision 31425) +++ trunk/Ohana/src/opihi/dvo/init.c (revision 31450) @@ -40,5 +40,6 @@ int lightcurve PROTO((int, char **)); int mextract PROTO((int, char **)); -int mmextract PROTO((int, char **)); +int mmextract PROTO((int, char **)); +int objectcoverage PROTO((int, char **)); int pcat PROTO((int, char **)); int photcodes PROTO((int, char **)); @@ -92,4 +93,5 @@ {1, "mextract", mextract, "extract measure data values"}, {1, "mmextract", mmextract, "extract joined measurements"}, + {1, "objectcoverage", objectcoverage, "plot catalog boundaries"}, {1, "pcat", skycat, "plot catalog boundaries"}, {1, "photcodes", photcodes, "list photometry codes"}, Index: trunk/Ohana/src/opihi/dvo/objectcoverage.c =================================================================== --- trunk/Ohana/src/opihi/dvo/objectcoverage.c (revision 31450) +++ trunk/Ohana/src/opihi/dvo/objectcoverage.c (revision 31450) @@ -0,0 +1,235 @@ +# include "dvoshell.h" + +int wordhash (char *word); + +int objectcoverage (int argc, char **argv) { + + void *Signal; + int ShowDensity; + int N, status, TimeSelect, xs, ys; + time_t tzero, tend; + double pixscale, r, d, Xs, Ys, trange, RaCenter, DecCenter; + char projection[16]; + float *V; + int k, j, invalid, Nx, Ny; + Buffer *buf; + Coords coords; + int Nsecfilt; + SkyList *skylist; + SkyRegionSelection *selection; + + Catalog catalog; + catalog.average = NULL; + catalog.secfilt = NULL; + catalog.measure = NULL; + skylist = NULL; + selection = NULL; + RaCenter = 0.0; + DecCenter = 0.0; + if ((selection = SetRegionSelection (&argc, argv)) == NULL) goto escape; + if ((N = get_argument (argc, argv, "-center"))) { + remove_argument (N, &argc, argv); + RaCenter = atof (argv[N]); + remove_argument (N, &argc, argv); + DecCenter = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + Nx = 0; + Ny = 0; + if ((N = get_argument (argc, argv, "-size"))) { + remove_argument (N, &argc, argv); + Nx = atof (argv[N]); + remove_argument (N, &argc, argv); + Ny = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + pixscale = 1.0; + if ((N = get_argument (argc, argv, "-scale"))) { + remove_argument (N, &argc, argv); + pixscale = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + ShowDensity = FALSE; + if ((N = get_argument (argc, argv, "-density"))) { + remove_argument (N, &argc, argv); + ShowDensity = TRUE; + } + + strcpy (projection, "DEC--AIT"); + if ((N = get_argument (argc, argv, "-proj"))) { + remove_argument (N, &argc, argv); + if (!strcasecmp(argv[N], "TAN")) { + strcpy (projection, "DEC--TAN"); + } + if (!strcasecmp(argv[N], "SIN")) { + strcpy (projection, "DEC--SIN"); + } + if (!strcasecmp(argv[N], "GLS")) { + strcpy (projection, "DEC--GLS"); + } + if (!strcasecmp(argv[N], "PAR")) { + strcpy (projection, "DEC--PAR"); + } + remove_argument (N, &argc, argv); + } + + TimeSelect = FALSE; + if ((N = get_argument (argc, argv, "-time"))) { + remove_argument (N, &argc, argv); + if (!ohana_str_to_time (argv[N], &tzero)) { + gprint (GP_ERR, "syntax error\n"); + return (FALSE); + } + remove_argument (N, &argc, argv); + if (!ohana_str_to_dtime (argv[N], &trange)) { + gprint (GP_ERR, "syntax error\n"); + return (FALSE); + } + remove_argument (N, &argc, argv); + if (trange < 0) { + trange = fabs (trange); + tzero -= trange; + } + TimeSelect = TRUE; + } + if ((N = get_argument (argc, argv, "-trange"))) { + remove_argument (N, &argc, argv); + if (!ohana_str_to_time (argv[N], &tzero)) { + gprint (GP_ERR, "syntax error\n"); + return (FALSE); + } + remove_argument (N, &argc, argv); + if (!ohana_str_to_time (argv[N], &tend)) { + gprint (GP_ERR, "syntax error\n"); + return (FALSE); + } + remove_argument (N, &argc, argv); + trange = tend - tzero; + if (trange < 0) { + trange = fabs (trange); + tzero -= trange; + } + TimeSelect = TRUE; + } + + if (argc != 3) { + gprint (GP_ERR, "USAGE: objectcoverage (buffer) (photcode)\n"); + gprint (GP_ERR, " options: [-scale pixscale] [-center ra dec] [-size Nx Nx] [-proj projection] [-time start range] [-trange start stop] [-name name] [+mosaic] [-mosaic] [-density]\n"); + gprint (GP_ERR, " (buffer) saves bitmapped image\n"); + gprint (GP_ERR, " (photcode) ..........\n"); + gprint (GP_ERR, " -scale (pixscale) : specifies the pixel size in degrees [1.0]\n"); + gprint (GP_ERR, " -center (ra) (dec) : specifies the center of the field [0.0, 0.0]\n"); + gprint (GP_ERR, " -size (Nx) (Ny) : specifies the size of the image [360/scale, 180/scale]\n"); + gprint (GP_ERR, " -proj (projection) : specifies the projection choice [AIT]\n"); + gprint (GP_ERR, " -density : create image with relative density (else binary on/off)\n"); + gprint (GP_ERR, " note: we need 64800 / (pixscale)^2 pixels to represent the sky\n"); + return (FALSE); + } + + if ((buf = SelectBuffer (argv[1], ANYBUFFER, TRUE)) == NULL) return (FALSE); + int myCode = GetPhotcodeCodebyName (argv[2]); + if (myCode == 0) { + gprint (GP_ERR, "invalid photcode\n"); + return (FALSE); + } + + if (!Nx || !Ny) { + Nx = 360/pixscale; + Ny = 180/pixscale; + } + + gfits_free_matrix (&buf[0].matrix); + gfits_free_header (&buf[0].header); + CreateBuffer (buf, Nx, Ny, -32, 0.0, 1.0); + strcpy (buf[0].file, "(empty)"); + + coords.crval1 = RaCenter; + coords.crval2 = DecCenter; + coords.crpix1 = 0.5*Nx; + coords.crpix2 = 0.5*Ny; + strcpy (coords.ctype, projection); + coords.pc1_1 = -1.0; + coords.pc2_2 = +1.0; + coords.pc1_2 = 0.0; + coords.pc2_1 = 0.0; + coords.cdelt1 = coords.cdelt2 = pixscale; + coords.Npolyterms = 0; + + PutCoords (&coords, &buf[0].header); + + V = (float *)buf[0].matrix.buffer; + + for (ys = 0; ys < Ny; ys++) { + for (xs = 0; xs < Nx; xs++) { + status = XY_to_RD (&r, &d, (double)(xs), (double)(ys), &coords); + status &= (r >= 0); + status &= (r <= 360); + if (status) { + V[ys*Nx + xs] = ShowDensity ? 0 : 2; + } else { + V[ys*Nx + xs] = ShowDensity ? -1 : 0; + } + } + } + + Nsecfilt = GetPhotcodeNsecfilt(); + // grab data from all selected sky regions + Signal = signal (SIGINT, handle_interrupt); + interrupt = FALSE; + + /* load region corresponding to selection above */ + if ((skylist = SelectRegions (selection)) == NULL) goto escape; + + /* loop over regions, extract data for each region */ + for (k = 0; (k < skylist[0].Nregions) && !interrupt; k++) { + /* lock, load, unlock catalog */ + catalog.filename = skylist[0].filename[k]; + catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_SECF; + catalog.Nsecfilt = 0; + + // an error exit status here is a significant error + if (!dvo_catalog_open (&catalog, NULL, FALSE, "r")) { + fprintf (stderr, "ERROR: failure to open catalog file %s\n", catalog.filename); + exit (2); + } + dvo_catalog_unlock (&catalog); + + int Nsec = GetPhotcodeNsec(myCode); + + for (j = 0; j < catalog.Naverage; j++) { + if (catalog.average[j].Nmeasure < 2) { continue; } + + + if (catalog.secfilt[j*Nsecfilt+Nsec].Ncode < 2) { continue; } + + invalid = ((catalog.secfilt[j*Nsecfilt + Nsec].M < 1.0) || (isnan(catalog.secfilt[j*Nsecfilt + Nsec].M))); + if (invalid) continue; + + double r = catalog.average[j].R; + double d = catalog.average[j].D; + status = RD_to_XY (&Xs, &Ys, r, d, &coords); + if (Xs < 0) continue; + if (Ys < 0) continue; + if (Xs >= Nx) continue; + if (Ys >= Ny) continue; + if (status) { + xs = (int)Xs; + ys = (int)Ys; + if (ShowDensity) { + V[ys*Nx + xs] += 1; + } else { + V[ys*Nx + xs] = 1; + } + } + } + } + return (TRUE); + +escape: + return (FALSE); +} + + Index: trunk/Ohana/src/photdbc/src/make_subcatalog.c =================================================================== --- trunk/Ohana/src/photdbc/src/make_subcatalog.c (revision 31425) +++ trunk/Ohana/src/photdbc/src/make_subcatalog.c (revision 31450) @@ -71,4 +71,9 @@ } + // if the input catalog is an old type, generate the catID entries: + if (catalog[0].catformat < DVO_FORMAT_PS1_V1) { + subcatalog[0].average[Naverage].catID = catalog[0].catID; + } + minMag = 32; minSigma = 32; @@ -127,4 +132,10 @@ subcatalog[0].measure[Nmeasure] = catalog[0].measure[offset]; subcatalog[0].measure[Nmeasure].averef = Naverage; + + // if the input catalog is an old type, generate the catID entries: + if (catalog[0].catformat < DVO_FORMAT_PS1_V1) { + subcatalog[0].measure[Nmeasure].catID = catalog[0].catID; + } + Nmeasure ++; Nm ++; Index: trunk/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 31425) +++ trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 31450) @@ -189,5 +189,10 @@ coords.crval1 = R[0]; coords.crval2 = D[0]; - Tmean /= (float) N; + + if (FIT_TARGET == TARGET_HIGH_SPEED) { + Tmean = 0.5*(Tmax - Tmin); + } else { + Tmean /= (float) N; + } XVERB = FALSE && (catalog[i].measure[m].dM < 0.01) && (N == 6) && (mode == FIT_PM_ONLY); Index: trunk/Ohana/src/relphot/doc/allsky.txt =================================================================== --- trunk/Ohana/src/relphot/doc/allsky.txt (revision 31450) +++ trunk/Ohana/src/relphot/doc/allsky.txt (revision 31450) @@ -0,0 +1,50 @@ + +2011.04.12 + + I have modified the code to operate on multiple filter at a time, + and to have a limited density, with the goal of running the analysis + on a huge area. However, we are somewhat limited by the fact that the + 'measure' and 'average' structures are huge and have far more than + is needed for just the relphot (or even relastro) analysis. + + I am exploring the idea of using a stripped-down structure in the + relphot analysis -- it would be assigned in bcatalog. + + Here are the elements we actually need in relphot proper: + +average: + FIELD R, RA, double, RA, decimal degrees + FIELD D, DEC, double, DEC, decimal degrees + FIELD Nmeasure, NMEASURE, unsigned short, number of psf measurements + FIELD measureOffset, OFF_MEASURE, int, offset to first psf measurement + FIELD flags, FLAGS, uint32_t, average object flags (star; ghost; etc) + = 2*8 + 2 + 4 + 4 = 26 bytes + +measure: + FIELD dR, D_RA, float, RA offset, arcsec + FIELD dD, D_DEC, float, DEC offset, arcsec + FIELD M, MAG, float, catalog mag, mag + FIELD Mcal, M_CAL, float, image cal mag, mag + + FIELD dM, MAG_ERR, float, mag error, mag + FIELD airmass, AIRMASS, float, (airmass - 1), airmass + FIELD Xccd, X_CCD, float, X coord on chip, pixels + FIELD Yccd, Y_CCD, float, Y coord on chip, pixels + + FIELD t, TIME, int, time in seconds (UNIX) + FIELD averef, AVE_REF, unsigned int, reference to average entry + FIELD imageID, IMAGE_ID, unsigned int, reference to DVO image ID + FIELD dbFlags, DB_FLAGS, unsigned int, flags supplied by analysis in database + FIELD photcode, PHOTCODE, unsigned short, photcode + = 8*4 + 4*4 + 2 = 50 bytes + + (note that psfQual and photFlags are only used in bcatalog to select measurements) + +eg, 1 year 3 pi has 1 average and 20 measures per object (and 8 secfilt entries) + +* current: 112 + 8*24 + 20*148 = 3264 bytes / star +* propose: 26 + 8*24 + 20*50 = 1218 bytes / star + +for the full 3pi area, in 30GB, we could use 820 stars per square +degree (5740 per exposure), as opposed to the current limit of 300 +(2100 per exposure). Index: trunk/Ohana/src/relphot/include/relphot.h =================================================================== --- trunk/Ohana/src/relphot/include/relphot.h (revision 31425) +++ trunk/Ohana/src/relphot/include/relphot.h (revision 31450) @@ -11,4 +11,5 @@ unsigned int start; unsigned int stop; + short photcode; float Mcal; float dMcal; @@ -85,9 +86,14 @@ int RELPHOT_GRID_BINNING; -PhotCode *photcode; -int PhotNsec; -int PhotSec; +int *photseclist; +int Nphotcodes; +PhotCode **photcodes; +// int PhotSec; +// int PhotNsec; PhotCode *refPhotcode; + +int MaxDensityUse; +double MaxDensityValue; int AreaSelect; @@ -213,7 +219,7 @@ StatType statsMosaicX PROTO((Catalog *catalog)); StatType statsMosaicdM PROTO((Catalog *catalog)); -StatType statsStarN PROTO((Catalog *catalog, int Ncatalog)); -StatType statsStarS PROTO((Catalog *catalog, int Ncatalog)); -StatType statsStarX PROTO((Catalog *catalog, int Ncatalog)); +StatType statsStarN PROTO((Catalog *catalog, int Ncatalog, int Nsec, int seccode)); +StatType statsStarS PROTO((Catalog *catalog, int Ncatalog, int Nsec)); +StatType statsStarX PROTO((Catalog *catalog, int Ncatalog, int Nsec)); void wcatalog PROTO((Catalog *catalog)); void wimages PROTO((void)); @@ -223,2 +229,10 @@ void relphot_usage (void); void relphot_help (int argc, char **argv); + +off_t getImageByID (off_t ID); + +int rationalize_mosaics (); +int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog); + +int populate_tiny_values (Catalog *catalog); +int free_tiny_values (Catalog *catalog); Index: trunk/Ohana/src/relphot/src/GridOps.c =================================================================== --- trunk/Ohana/src/relphot/src/GridOps.c (revision 31425) +++ trunk/Ohana/src/relphot/src/GridOps.c (revision 31450) @@ -348,4 +348,8 @@ // sums for each star which touches as cell on both bases. + int Nsecfilt = GetPhotcodeNsecfilt (); + int thisCode = photcodes[0][0].code; + int Nsec = GetPhotcodeNsec(thisCode); + for (i = 0; i < Ngrid; i++) { @@ -362,5 +366,5 @@ mx = mlist[i][j]; c = clist[i][j]; - n = catalog[c].measure[mx].averef; + n = catalog[c].measureT[mx].averef; // if we have already visited this star, skip the stuff below @@ -369,5 +373,5 @@ // skip stars marked as BAD - if (catalog[c].average[n].flags & STAR_BAD) { + if (catalog[c].secfilt[n*Nsecfilt+Nsec].flags & STAR_BAD) { Nrel ++; continue; @@ -389,5 +393,5 @@ // skip measurements marked as BAD - if (catalog[c].measure[m].dbFlags & MEAS_BAD) { + if (catalog[c].measureT[m].dbFlags & MEAS_BAD) { Nbad ++; continue; @@ -410,6 +414,5 @@ // select the color- and airmass-corrected observed magnitude for this star // XXX need to be able to turn off the color-correction until initial average mags are found - // Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]); - Msys = PhotCat (&catalog[c].measure[m]); + Msys = PhotCatTiny (&catalog[c].measureT[m]); if (isnan(Msys)) { Nsys++; @@ -418,5 +421,5 @@ // mag-error for this measurement - Merr = MAX (catalog[c].measure[m].dM, MIN_ERROR); + Merr = MAX (catalog[c].measureT[m].dM, MIN_ERROR); // Wsys = 1.0 / SQ(Merr); @@ -528,4 +531,6 @@ if (!USE_GRID) return; + int Nsecfilt = GetPhotcodeNsecfilt (); + Nmax = Nlist[0]; for (i = 0; i < Ngrid; i++) { @@ -545,5 +550,5 @@ c = clist[i][j]; - if (catalog[c].measure[m].dbFlags & MEAS_BAD) { + if (catalog[c].measureT[m].dbFlags & MEAS_BAD) { Nbad ++; continue; @@ -565,6 +570,6 @@ } - n = catalog[c].measure[m].averef; - Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]); + n = catalog[c].measureT[m].averef; + Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]); if (isnan(Msys)) { Nsys++; @@ -572,5 +577,5 @@ } list[N] = Msys - Mrel - Mcal - Mmos; - dlist[N] = MAX (catalog[c].measure[m].dM, MIN_ERROR); + dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); N++; } @@ -613,4 +618,6 @@ if (!USE_GRID) return; + int Nsecfilt = GetPhotcodeNsecfilt (); + N = 0; for (i = 0; i < Ngrid; i++) @@ -630,5 +637,5 @@ c = clist[i][j]; - if (catalog[c].measure[m].dbFlags & MEAS_BAD) { + if (catalog[c].measureT[m].dbFlags & MEAS_BAD) { Narea ++; continue; @@ -641,6 +648,6 @@ if (isnan(Mrel)) continue; - n = catalog[c].measure[m].averef; - Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]); + n = catalog[c].measureT[m].averef; + Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]); xlist[N] = Xmeas[c][m]; @@ -704,5 +711,5 @@ gfits_create_matrix (&header, &matrix); gfits_modify (&header, "NEXTEND", OFF_T_FMT, 1, Nimage + 3); - gfits_modify (&header, "FILTER", "%s", 1, photcode[0].name); + gfits_modify (&header, "FILTER", "%s", 1, photcodes[0][0].name); // XXXX note that this expects a single photcode, enforced in initialize.d gfits_modify_alt (&header, "COMMENT", "%S", 1, "Mosaic Photometry Grid Analysis"); @@ -727,5 +734,5 @@ theader.bitpix = -32; gfits_create_Theader (&theader, "IMAGE"); - gfits_modify (&theader, "FILTER", "%s", 1, photcode[0].name); + gfits_modify (&theader, "FILTER", "%s", 1, photcodes[0][0].name); gfits_modify (&theader, "EXTNAME", "%s", 1, "MAG_OFFSET"); gfits_create_matrix (&theader, &matrix); @@ -742,5 +749,5 @@ /* save grid Nmeas values */ gfits_modify (&theader, "EXTNAME", "%s", 1, "NMEAS"); - gfits_modify (&theader, "FILTER", "%s", 1, photcode[0].name); + gfits_modify (&theader, "FILTER", "%s", 1, photcodes[0][0].name); gfits_create_matrix (&theader, &matrix); for (i = 0; i < gridX; i++) { @@ -756,5 +763,5 @@ /* save grid sigma values */ gfits_modify (&theader, "EXTNAME", "%s", 1, "SIGMA"); - gfits_modify (&theader, "FILTER", "%s", 1, photcode[0].name); + gfits_modify (&theader, "FILTER", "%s", 1, photcodes[0][0].name); gfits_create_matrix (&theader, &matrix); for (i = 0; i < gridX; i++) { @@ -776,5 +783,5 @@ gfits_modify (&theader, "EXTNAME", "%s", 1, camera.ccdname[N]); - gfits_modify (&theader, "FILTER", "%s", 1, photcode[0].name); + gfits_modify (&theader, "FILTER", "%s", 1, photcodes[0][0].name); gfits_modify (&theader, "NX", "%d", 1, camera.Nx); gfits_modify (&theader, "NY", "%d", 1, camera.Ny); Index: trunk/Ohana/src/relphot/src/ImageOps.c =================================================================== --- trunk/Ohana/src/relphot/src/ImageOps.c (revision 31425) +++ trunk/Ohana/src/relphot/src/ImageOps.c (revision 31450) @@ -1,5 +1,3 @@ # include "relphot.h" - -# define USE_IMAGE_ID 1 static off_t **bin; // link from catalog,measure to image @@ -12,14 +10,7 @@ static off_t Nimage; // number of available images -// if we search by image ID, we sort (imageIDs, imageIdx) by imageIDs to get a sorted -// index - -# if USE_IMAGE_ID +// to search by image ID, we sort (imageIDs, imageIdx) by imageIDs to get a sorted index static off_t *imageIDs; // list of all image IDs static off_t *imageIdx; // list of index for image IDs -# else -static unsigned int *start; -static unsigned int *stop; -# endif Image *getimages (off_t *N) { @@ -40,5 +31,4 @@ Nimage = N; -# if USE_IMAGE_ID ALLOCATE (imageIDs, off_t, Nimage); ALLOCATE (imageIdx, off_t, Nimage); @@ -49,13 +39,4 @@ } llsortpair (imageIDs, imageIdx, Nimage); -# else - ALLOCATE (start, unsigned, Nimage); - ALLOCATE (stop, unsigned, Nimage); - - for (i = 0; i < Nimage; i++) { - start[i] = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1); - stop[i] = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1); - } -# endif } @@ -65,5 +46,4 @@ // use bisection to find the specified image ID -# if USE_IMAGE_ID off_t Nlo, Nhi, N; @@ -82,5 +62,4 @@ return (imageIdx[N]); } -# endif return (-1); @@ -128,22 +107,29 @@ } -/* select all image equivalent to the current photcode */ +/* select all image equivalent to the active photcode set */ void findImages (Catalog *catalog, int Ncatalog) { off_t j; - int i, ecode; - - for (i = 0; i < Ncatalog; i++) { + int i, ecode, Ns, found; + + int Nmatch = 0; + for (i = 0; i < Ncatalog; i++) { for (j = 0; j < catalog[i].Nmeasure; j++) { - ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[j].photcode); - if (photcode[0].code != ecode) 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; matchImage (catalog, j, i); - } - } -} - -int findCCD (off_t idx, off_t meas, int cat, Measure *measure) { - - int ccdnum; + Nmatch ++; + } + } + // fprintf (stderr, "matched %d detections to images\n", Nmatch); +} + +int findCCD (off_t idx, off_t meas, int cat, MeasureTiny *measure) { + + int ccdnum, found, Ns, ecode; double X, Y; char *pname, *filter, *p, base[256]; @@ -151,5 +137,14 @@ /* identify the ccd on the basis of the photcode name */ pname = GetPhotcodeNamebyCode (image[idx].photcode); - filter = photcode[0].name; + + // XXX this seems quite terrible... + ecode = GetPhotcodeEquivCodebyCode (measure[0].photcode); + found = FALSE; + for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) { + if (ecode == photcodes[Ns][0].code) found = TRUE; + } + if (!found) return (FALSE); + + filter = photcodes[Ns][0].name; sprintf (base, "%s.%s.", MOSAICNAME, filter); if (strncmp (pname, base, strlen (base))) return (FALSE); @@ -170,7 +165,7 @@ // old code to add this measurement to the grid cell for this chip - // ave = measure[0].averef; - // ra = catalog[cat].average[ave].R - measure[0].dR / 3600.0; - // dec = catalog[cat].average[ave].D - measure[0].dD / 3600.0; + // ave = measureT[0].averef; + // ra = catalog[cat].averageT[ave].R - measureT[0].dR / 3600.0; + // dec = catalog[cat].averageT[ave].D - measureT[0].dD / 3600.0; // RD_to_XY (&X, &Y, ra, dec, &image[i].coords); @@ -183,12 +178,11 @@ } -# if USE_IMAGE_ID void matchImage (Catalog *catalog, off_t meas, int cat) { off_t idx, ID; int status; - Measure *measure; + MeasureTiny *measure; - measure = &catalog[cat].measure[meas]; + measure = &catalog[cat].measureT[meas]; ID = measure[0].imageID; @@ -226,45 +220,4 @@ } -# else -// this is the time-based match -void matchImage (Catalog *catalog, off_t meas, int cat) { - - off_t i; - int status; - Measure *measure; - - measure = &catalog[cat].measure[meas]; - for (i = 0; i < Nimage; i++) { - if (image[0].photcode == -1) continue; - if (measure[0].photcode != image[i].photcode) continue; - if (measure[0].t < start[i]) continue; - if (measure[0].t > stop[i]) continue; - - if (USE_GRID) { - status = findCCD(i, meas, cat); - if (!status) continue; - } - - // index for (catalog, measure) -> image - bin[cat][meas] = i; - - // index for image, Nentry -> catalog - clist[i][Nlist[i]] = cat; - - // index for image, Nentry -> measure - mlist[i][Nlist[i]] = meas; - Nlist[i] ++; - - if (Nlist[i] == NLIST[i]) { - NLIST[i] += 100; - REALLOCATE (clist[i], off_t, NLIST[i]); - REALLOCATE (mlist[i], off_t, NLIST[i]); - } - return; - } - /* fprintf (stderr, "can't find source image for this measurement: %d (%d)\n", measure[0].t, measure[0].photcode); */ -} -# endif - off_t getImageEntry (off_t meas, int cat) { @@ -307,4 +260,6 @@ if (FREEZE_IMAGES) return; + + int Nsecfilt = GetPhotcodeNsecfilt (); if (PoorImages) { @@ -339,5 +294,5 @@ c = clist[i][j]; - if (catalog[c].measure[m].dbFlags & MEAS_BAD) { + if (catalog[c].measureT[m].dbFlags & MEAS_BAD) { Nbad++; continue; @@ -359,6 +314,6 @@ } - n = catalog[c].measure[m].averef; - Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]); + n = catalog[c].measureT[m].averef; + Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]); if (isnan(Msys)) { Nsys++; @@ -366,6 +321,6 @@ } list[N] = Msys - Mrel - Mmos - Mgrid; - dlist[N] = MAX (catalog[c].measure[m].dM, MIN_ERROR); - if (catalog[c].measure[m].dM < IMFIT_SYS_SIGMA_LIM) { + dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); + if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) { McalBright += list[N]; McalBright2 += SQ(list[N]); @@ -632,3 +587,2 @@ return (stats); } - Index: trunk/Ohana/src/relphot/src/MosaicOps.c =================================================================== --- trunk/Ohana/src/relphot/src/MosaicOps.c (revision 31425) +++ trunk/Ohana/src/relphot/src/MosaicOps.c (revision 31450) @@ -9,4 +9,7 @@ static off_t **imlist; /* mosaic -> image[] */ static off_t **bin; /* catalog, measure -> mosaic */ + +// list of mosaics associated with an image +static off_t *mosimage; // list of mosaic associated with each image @@ -33,4 +36,6 @@ ALLOCATE (NIMLIST, off_t, NMOSAIC); + ALLOCATE (mosimage, off_t, Nimage); // mosaic to which image belongs + /* a 'mosaic' in relphot is (unlike relastro) a virtual concept: there is no * entry in the image table that represents this mosaic. Instead, it is an @@ -40,4 +45,5 @@ /* generate list of unique mosaics */ for (i = 0; i < Nimage; i++) { + mosimage[i] = -1; /* select valid mosaic images by photcode */ @@ -56,4 +62,7 @@ if (start > mosaic[j].stop) continue; found = TRUE; + + // add reference from image to mosaic + mosimage[i] = j; /* add image to mosaic image list */ @@ -77,4 +86,5 @@ mosaic[Nmosaic].flags = image[i].flags; mosaic[Nmosaic].secz = image[i].secz; + mosaic[Nmosaic].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode); /* add image to mosaic image list */ @@ -84,4 +94,7 @@ imlist[Nmosaic][0] = i; + // add reference from image to mosaic + mosimage[i] = Nmosaic; + Nmosaic ++; if (Nmosaic == NMOSAIC) { @@ -250,20 +263,27 @@ int findMosaics (Catalog *catalog, int Ncatalog) { - int i, ecode; + int i, ecode, found, Ns; off_t j; if (!MOSAIC_ZEROPT) return (FALSE); + int Nmatch = 0; for (i = 0; i < Ncatalog; i++) { for (j = 0; j < catalog[i].Nmeasure; j++) { if (TimeSelect) { - if (catalog[i].measure[j].t < TSTART) continue; - if (catalog[i].measure[j].t > TSTOP) continue; - } - ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[j].photcode); - if (photcode[0].code != ecode) continue; + 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; matchMosaics (catalog, j, i); - } - } + Nmatch ++; + } + } + // fprintf (stderr, "matched %d detections to mosaics\n", Nmatch); return (TRUE); } @@ -271,36 +291,44 @@ void matchMosaics (Catalog *catalog, off_t meas, int cat) { - int i; - - for (i = 0; i < Nmosaic; i++) { - if (catalog[cat].measure[meas].t < mosaic[i].start) continue; - if (catalog[cat].measure[meas].t > mosaic[i].stop) continue; + off_t idx, ID, mosID; + MeasureTiny *measure; + + measure = &catalog[cat].measureT[meas]; + + ID = measure[0].imageID; + idx = getImageByID (ID); + if (idx == -1) { + if (VERBOSE2) fprintf (stderr, "missed measurement "OFF_T_FMT", %d\n", meas, cat); + return; + } + + mosID = mosimage[idx]; + if (mosID < 0) { + Image *image = getimage(idx); + fprintf (stderr, "unmatched image %s\n", image[0].name); + return; + } + + // test to check we got the right match: + { + Image *image = getimage(idx); + unsigned int imageStart = image[0].tzero - MAX(0.01*image[0].trate*image[0].NY, 1); + if (imageStart != mosaic[mosID].start) { + fprintf (stderr, "error in image to mosaic match\n"); + abort(); + } + } + + bin[cat][meas] = mosID; + + clist[mosID][Nlist[mosID]] = cat; + mlist[mosID][Nlist[mosID]] = meas; + Nlist[mosID] ++; -# ifdef GRID_V1 - if (USE_GRID) { - ave = catalog[cat].measure[meas].averef; - ra = catalog[cat].average[ave].R - catalog[cat].measure[meas].dR / 3600.0; - dec = catalog[cat].average[ave].D - catalog[cat].measure[meas].dD / 3600.0; - - /* X,Y always positive-definite in range 0,0 - dX, dY */ - RD_to_XY (&X, &Y, ra, dec, &mosaic[i].coords); - setGridMeasure (meas, cat, X, Y); - } -# endif - - bin[cat][meas] = i; - - clist[i][Nlist[i]] = cat; - mlist[i][Nlist[i]] = meas; - Nlist[i] ++; - - if (Nlist[i] == NLIST[i]) { - NLIST[i] += 100; - REALLOCATE (clist[i], int, NLIST[i]); - REALLOCATE (mlist[i], off_t, NLIST[i]); - } - return; - } - fprintf (stderr, "missed measurement\n"); + if (Nlist[mosID] == NLIST[mosID]) { + NLIST[mosID] += 100; + REALLOCATE (clist[mosID], int, NLIST[mosID]); + REALLOCATE (mlist[mosID], off_t, NLIST[mosID]); + } return; } @@ -324,5 +352,5 @@ off_t i, j, m, c, N, Nmax; - int n, mark, bad, Nfew, Nbad, Ncal, Nrel, Ngrid, Nsys; + int n, mark, bad, Nfew, Nbad, Ncal, Nrel, Ngrid, Nsys, Nsecfilt; float Msys, Mrel, Mcal, Mgrid; double *list, *dlist, *Mlist, *dMlist; @@ -334,4 +362,6 @@ image = getimages (&N); + + Nsecfilt = GetPhotcodeNsecfilt (); if (PoorImages) { @@ -364,5 +394,5 @@ c = clist[i][j]; - if (catalog[c].measure[m].dbFlags & MEAS_BAD) { + if (catalog[c].measureT[m].dbFlags & MEAS_BAD) { Nbad ++; continue; @@ -384,6 +414,6 @@ } - n = catalog[c].measure[m].averef; - Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]); + n = catalog[c].measureT[m].averef; + Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]); if (isnan(Msys)) { Nsys++; @@ -391,5 +421,5 @@ } list[N] = Msys - Mrel - Mcal - Mgrid; - dlist[N] = MAX (catalog[c].measure[m].dM, MIN_ERROR); + dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); Mlist[N] = Msys; dMlist[N] = list[N]; @@ -433,4 +463,198 @@ } +// When we rationalize the images/mosaics, we are driving the negative cloud images back +// to 0.0. At the same time, we make a guess to the effective impact on all other images, +// driven by the coupling of common stars. +int rationalize_mosaics (Catalog *catalog, int Ncatalog) { + + double *imageOffset, **starOffset; + int **starNcount, *seclist; + int **Slist, *NSlist, *NSLIST; + int i, j, k, m, nMos, Ns, found; + + off_t Nimage; + Image *image; + + // set a test value for now + float CLOUD_TOLERANCE = 0.025; + + if (!MOSAIC_ZEROPT) return (FALSE); + if (FREEZE_MOSAICS) return (FALSE); + + image = getimages (&Nimage); + + ALLOCATE (imageOffset, double, Nmosaic); + + ALLOCATE ( Slist, int *, Nmosaic); // array of calibrated star indexes on this mosaic + ALLOCATE (NSlist, int, Nmosaic); // number of stars on mosaic + ALLOCATE (NSLIST, int, Nmosaic); // number of Slist entries allocated + memset (Slist, 0, Nmosaic*sizeof(int *)); + + // find the images / mosaics with negative clouds and save their offset + for (i = 0; i < Nmosaic; i++) { + + NSlist[i] = 0; + NSLIST[i] = 100; + ALLOCATE (Slist[i], int, NSLIST[i]); + + imageOffset[i] = 0.0; + + if (VERBOSE2 && (fabs(mosaic[i].Mcal) < CLOUD_TOLERANCE)) { + fprintf (stderr, "cloud-free: %s : %f\n", image[imlist[i][0]].name, mosaic[i].Mcal); + } + if (VERBOSE2 && (mosaic[i].Mcal < -CLOUD_TOLERANCE)) { + imageOffset[i] = -mosaic[i].Mcal; + // NOTE the negative sign: down below, we are going to add in the negative of Mcal + // to this image, and the propagated mean values for other images + fprintf (stderr, "anti-clouds: %s : %f\n", image[imlist[i][0]].name, mosaic[i].Mcal); + } + if (VERBOSE2 && (mosaic[i].Mcal > CLOUD_TOLERANCE)) { + fprintf (stderr, "cloudy : %s : %f\n", image[imlist[i][0]].name, mosaic[i].Mcal); + } + } + + int Nsecfilt = GetPhotcodeNsecfilt (); + ALLOCATE (seclist, int, Nphotcodes); + for (Ns = 0; Ns < Nphotcodes; Ns ++) { + int thisCode = photcodes[Ns][0].code; + seclist[Ns] = GetPhotcodeNsec(thisCode); + } + + // allocate an array for star offsets + int Nstars = 0; + int NSTARS = 1000; + ALLOCATE (starOffset, double *, NSTARS); + ALLOCATE (starNcount, int *, NSTARS); + memset (starOffset, 0, NSTARS*sizeof(double *)); + memset (starNcount, 0, NSTARS*sizeof(int *)); + + // find the mean offset for each star + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + ALLOCATE (starOffset[Nstars], double, Nphotcodes); + ALLOCATE (starNcount[Nstars], int, Nphotcodes); + memset (starOffset[Nstars], 0, Nphotcodes*sizeof(double)); + memset (starNcount[Nstars], 0, Nphotcodes*sizeof(int)); + + m = catalog[i].averageT[j].measureOffset; + + // determine the mosaic for each measurement + for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) { + + // skip unused measurements + if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; + + // skip unused measurements + int Nsec; + int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode); + found = FALSE; + for (Ns = 0; !found && (Ns < Nphotcodes); Ns ++) { + if (ecode == photcodes[Ns][0].code) { + found = TRUE; + break; + } + } + if (!found) continue; + Nsec = seclist[Ns]; + + // bad stars for this secfilt + if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; + + // skip REF measurements (not tied to an image) + if (getImageEntry (m, i) < 0) continue; + + // find the source of this measurement (skip unassigned measurements) + nMos = bin[i][m]; + if (nMos == -1) continue; + + if (mosaic[nMos].photcode != ecode) { + fprintf (stderr, "*"); + } + + assert (Ns >= 0); + assert (Ns < Nphotcodes); + + // accumulate the offsets from the negative cloud images (others have 0.0 value) + starOffset[Nstars][Ns] += imageOffset[nMos]; + starNcount[Nstars][Ns] ++; + + // record the mosaic->star reference + Slist[nMos][NSlist[nMos]] = Nstars; + NSlist[nMos] ++; + if (NSlist[nMos] == NSLIST[nMos]) { + NSLIST[nMos] += 100; + REALLOCATE (Slist[nMos], int, NSLIST[nMos]); + } + } + Nstars ++; + if (Nstars == NSTARS) { + NSTARS += 1000; + REALLOCATE (starOffset, double *, NSTARS); + REALLOCATE (starNcount, int *, NSTARS); + memset (&starOffset[NSTARS-1000], 0, 1000*sizeof(double *)); + memset (&starNcount[NSTARS-1000], 0, 1000*sizeof(int *)); + } + } + } + + // find the mean offset of the images without negative clouds + for (i = 0; i < Nmosaic; i++) { + + found = FALSE; + for (Ns = 0; !found && (Ns < Nphotcodes); Ns ++) { + if (mosaic[i].photcode == photcodes[Ns][0].code) { + found = TRUE; + break; + } + } + if (!found) { + fprintf (stderr, "invalid photcode for mosaic?\n"); + abort(); + } + + // a negative cloud image (cloud: Mcal > 0; anti-clouds: Mcal < 0; imageOffset = -Mcal) + if (imageOffset[i] > 0.0) continue; + + // we need to actually have cross-references to count + if (NSlist[i] < 2) continue; + + float dM = 0.0; + for (j = 0; j < NSlist[i]; j++) { + Nstars = Slist[i][j]; + if (starNcount[Nstars][Ns] > 1) { + dM += (starOffset[Nstars][Ns] / starNcount[Nstars][Ns]); + } + } + imageOffset[i] = dM / NSlist[i]; + // fprintf (stderr, "adjust image: %s : (%f %d) : %f\n", image[imlist[i][0]].name, dM, NSlist[i], imageOffset[i]); + } + + // for (i = 0; i < Nmosaic; i++) { + // fprintf (stderr, "correction: %s : %f\n", image[imlist[i][0]].name, imageOffset[i]); + // } + + // apply all offset values to the mosaics + // find the images / mosaics with negative clouds and save their offset + for (i = 0; i < Nmosaic; i++) { + mosaic[i].Mcal += imageOffset[i]; + } + + for (i = 0; i < Nstars; i++) { + free (starOffset[i]); + free (starNcount[i]); + } + free (starOffset); + free (starNcount); + + for (i = 0; i < Nmosaic; i++){ + free (Slist[i]); + } + free (NSlist); + free (NSLIST); + free (imageOffset); + free (seclist); + return (TRUE); +} + StatType statsMosaicM (Catalog *catalog) { @@ -524,5 +748,5 @@ n++; } - fprintf (stderr, "Nmosaic: "OFF_T_FMT", n: "OFF_T_FMT"\n", Nmosaic, n); + // fprintf (stderr, "Nmosaic: "OFF_T_FMT", n: "OFF_T_FMT"\n", Nmosaic, n); liststats (list, dlist, n, &stats); @@ -563,5 +787,5 @@ void clean_mosaics () { - off_t i, N, mark, Nmark; + off_t i, N, mark, Nmark, Nscatter, Noffset; double *mlist, *slist, *dlist; double MaxOffset, MedOffset, MaxScatter; @@ -592,8 +816,15 @@ fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset); - Nmark = 0; + Nmark = Nscatter = Noffset = 0; for (i = 0; i < Nmosaic; i++) { mark = FALSE; - mark = (mosaic[i].dMcal > MaxScatter) || (fabs(mosaic[i].Mcal - MedOffset) > MaxOffset); + if (mosaic[i].dMcal > MaxScatter) { + mark = TRUE; + Nscatter ++; + } + if (fabs(mosaic[i].Mcal - MedOffset) > MaxOffset) { + mark = TRUE; + Noffset ++; + } if (mark) { Nmark ++; @@ -604,5 +835,5 @@ } - fprintf (stderr, OFF_T_FMT" mosaics marked poor\n", Nmark); + fprintf (stderr, OFF_T_FMT" mosaics marked poor ("OFF_T_FMT" scatter, "OFF_T_FMT" offset)\n", Nmark, Nscatter, Noffset); initstats (STATMODE); free (mlist); @@ -640,9 +871,9 @@ c = clist[i][j]; - if (catalog[c].measure[m].dbFlags & (ID_MEAS_AREA | ID_MEAS_NOCAL)) continue; - - ave = catalog[c].measure[m].averef; - xlist[N] = catalog[c].average[ave].R - catalog[c].measure[m].dR / 3600.0; - ylist[N] = catalog[c].average[ave].D - catalog[c].measure[m].dD / 3600.0; + if (catalog[c].measureT[m].dbFlags & (ID_MEAS_AREA | ID_MEAS_NOCAL)) continue; + + ave = catalog[c].measureT[m].averef; + xlist[N] = catalog[c].averageT[ave].R - catalog[c].measureT[m].dR / 3600.0; + ylist[N] = catalog[c].averageT[ave].D - catalog[c].measureT[m].dD / 3600.0; N++; } Index: trunk/Ohana/src/relphot/src/StarOps.c =================================================================== --- trunk/Ohana/src/relphot/src/StarOps.c (revision 31425) +++ trunk/Ohana/src/relphot/src/StarOps.c (revision 31450) @@ -5,4 +5,10 @@ static double *dlist; +// When we rationalize the images/mosaics, we are driving the negative cloud images back +// to 0.0. At the same time, we make a guess to the effective impact on all other images, +// driven by the coupling of common stars. This array carries the impact of those offsets +// on each star +static double *Moffset; + void initMrel (Catalog *catalog, int Ncatalog) { @@ -12,21 +18,30 @@ for (i = 0; i < Ncatalog; i++) { for (j = 0; j < catalog[i].Naverage; j++) { - Nmax = MAX (Nmax, catalog[i].average[j].Nmeasure); - } - } - - ALLOCATE (list, double, MAX (1, Nmax)); - ALLOCATE (dlist, double, MAX (1, Nmax)); + Nmax = MAX (Nmax, catalog[i].averageT[j].Nmeasure); + } + } + + ALLOCATE (list, double, MAX (1, Nmax)); + ALLOCATE (dlist, double, MAX (1, Nmax)); + ALLOCATE (Moffset, double, MAX (1, Nmax)); } float getMrel (Catalog *catalog, off_t meas, int cat) { + int Nsec, Nsecfilt, photcode; int ave; float value; - ave = catalog[cat].measure[meas].averef; - if (catalog[cat].average[ave].flags & STAR_BAD) return (NAN); + ave = catalog[cat].measureT[meas].averef; + photcode = catalog[cat].measureT[meas].photcode; + + int ecode = GetPhotcodeEquivCodebyCode (photcode); + Nsec = GetPhotcodeNsec(ecode); + Nsecfilt = GetPhotcodeNsecfilt (); + + // is this star OK? + if (catalog[cat].secfilt[Nsecfilt*ave+Nsec].flags & STAR_BAD) return (NAN); - value = catalog[cat].secfilt[PhotNsec*ave+PhotSec].M; + value = catalog[cat].secfilt[Nsecfilt*ave+Nsec].M; return (value); } @@ -39,4 +54,5 @@ StatType stats; + int Nsecfilt = GetPhotcodeNsecfilt (); Nfew = Nsys = Nbad = Ncal = Nmos = Ngrid = 0; @@ -44,68 +60,81 @@ for (j = 0; j < catalog[i].Naverage; j++) { - /* calculate the average value for a single star */ - if (catalog[i].average[j].flags & STAR_BAD) continue; - m = catalog[i].average[j].measureOffset; - - N = 0; - for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) { - if (catalog[i].measure[m].dbFlags & MEAS_BAD) { - Nbad ++; - continue; - } - // XXX allow REF stars (no Image Entry) to be included in the calculation this - // should be optionally set, and should allow for REF stars to be downweighted by - // more than their reported errors. how such info is carried is unclear... - if (getImageEntry (m, i) < 0) { - Mcal = Mmos = Mgrid = 0; - } else { - Mcal = getMcal (m, i); - if (isnan(Mcal)) { - Ncal ++; + int Ns; + for (Ns = 0; Ns < Nphotcodes; Ns++) { + + int thisCode = photcodes[Ns][0].code; + int Nsec = GetPhotcodeNsec(thisCode); + + /* calculate the average value for a single star */ + + // skip bad stars + if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; + m = catalog[i].averageT[j].measureOffset; + + N = 0; + for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) { + // skip measurements that do not match the current photcode + int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode); + if (ecode != thisCode) { continue; } + + if (catalog[i].measureT[m].dbFlags & MEAS_BAD) { + Nbad ++; continue; } - Mmos = getMmos (m, i); - if (isnan(Mmos)) { - Nmos ++; + // XXX allow REF stars (no Image Entry) to be included in the calculation this + // should be optionally set, and should allow for REF stars to be downweighted by + // more than their reported errors. how such info is carried is unclear... + if (getImageEntry (m, i) < 0) { + Mcal = Mmos = Mgrid = 0; + } else { + Mcal = getMcal (m, i); + if (isnan(Mcal)) { + Ncal ++; + continue; + } + Mmos = getMmos (m, i); + if (isnan(Mmos)) { + Nmos ++; + 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]); + if (isnan(Msys)) { + Nsys++; continue; } - Mgrid = getMgrid (m, i); - if (isnan(Mgrid)) { - Ngrid++; - continue; + list[N] = Msys - Mcal - Mmos - Mgrid; + dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR); + + // tie down reference photometry if the -refcode (code) option is selected + if (refPhotcode) { + if (GetPhotcodeEquivCodebyCode(catalog[i].measureT[m].photcode) == refPhotcode[0].equiv) { + // increase the weight by a factor of 100: + dlist[N] = 0.01*catalog[i].measureT[m].dM; + } } - } - - Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]); - if (isnan(Msys)) { - Nsys++; - continue; - } - list[N] = Msys - Mcal - Mmos - Mgrid; - dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR); - - if (refPhotcode) { - if (GetPhotcodeEquivCodebyCode(catalog[i].measure[m].photcode) == refPhotcode[0].equiv) { - // increase the weight by a factor of 100: - dlist[N] = 0.01*catalog[i].measure[m].dM; - } - } - N++; - } - - // when performing the grid analysis, STAR_TOOFEW will be set to 1; - - if (N <= STAR_TOOFEW) { /* too few measurements */ - catalog[i].average[j].flags |= ID_STAR_FEW; - Nfew ++; - } else { - catalog[i].average[j].flags &= ~ID_STAR_FEW; - } - - liststats (list, dlist, N, &stats); - - catalog[i].secfilt[PhotNsec*j+PhotSec].M = stats.mean; - catalog[i].secfilt[PhotNsec*j+PhotSec].dM = stats.sigma; - catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT; + N++; + } + + // when performing the grid analysis, STAR_TOOFEW will be set to 1; + if (N <= STAR_TOOFEW) { /* too few measurements */ + catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_FEW; + Nfew ++; + } else { + catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_FEW; + } + + liststats (list, dlist, N, &stats); + + catalog[i].secfilt[Nsecfilt*j+Nsec].M = stats.mean; + catalog[i].secfilt[Nsecfilt*j+Nsec].dM = stats.sigma; + catalog[i].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT; + } } } @@ -122,4 +151,7 @@ double *list, *dlist; StatType stats; + int Nsec, Nsecfilt, ecode; + + Nsecfilt = GetPhotcodeNsecfilt (); /* Nmeasure is now different, need to reallocate */ @@ -127,5 +159,5 @@ for (i = 0; i < Ncatalog; i++) { for (j = 0; j < catalog[i].Naverage; j++) { - Nmax = MAX (Nmax, catalog[i].average[j].Nmeasure); + Nmax = MAX (Nmax, catalog[i].averageT[j].Nmeasure); } } @@ -135,40 +167,50 @@ for (i = 0; i < Ncatalog; i++) { for (j = 0; j < catalog[i].Naverage; j++) { - /* skip stars already calibrated */ if (catalog[i].found[j]) continue; - N = 0; - m = catalog[i].average[j].measureOffset; - for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) { - if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue; - // XXX allow REF stars (no Image Entry) to be included in the calculation this - // should be optionally set, and should allow for REF stars to be downweighted by - // more than their reported errors. how such info is carried is unclear... - if (getImageEntry (m, i) < 0) { - Mcal = Mmos = Mgrid = 0; - } else { - Mcal = getMcal (m, i); - if (isnan(Mcal)) continue; - Mmos = getMmos (m, i); - if (isnan(Mmos)) continue; - Mgrid = getMgrid (m, i); - if (isnan(Mgrid)) continue; - } - - Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]); - list[N] = Msys - Mcal - Mmos - Mgrid; - dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR); - N++; - } - if (N < 1) continue; - - liststats (list, dlist, N, &stats); - if (mark) catalog[i].found[j] = TRUE; - - /* use sigma or error in dM for output? */ - catalog[i].secfilt[PhotNsec*j+PhotSec].M = stats.mean; - catalog[i].secfilt[PhotNsec*j+PhotSec].dM = MAX (stats.error, stats.sigma); - catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT; + int Ns; + for (Ns = 0; Ns < Nphotcodes; Ns++) { + int thisCode = photcodes[Ns][0].code; + Nsec = GetPhotcodeNsec(thisCode); + + N = 0; + m = catalog[i].averageT[j].measureOffset; + for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) { + // skip measurements that do not match the current photcode + ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode); + if (ecode != thisCode) { continue; } + + if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; + + // XXX allow REF stars (no Image Entry) to be included in the calculation this + // should be optionally set, and should allow for REF stars to be downweighted by + // more than their reported errors. how such info is carried is unclear... + if (getImageEntry (m, i) < 0) { + Mcal = Mmos = Mgrid = 0; + } else { + Mcal = getMcal (m, i); + if (isnan(Mcal)) continue; + Mmos = getMmos (m, i); + if (isnan(Mmos)) continue; + Mgrid = getMgrid (m, i); + if (isnan(Mgrid)) continue; + } + + Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt]); + list[N] = Msys - Mcal - Mmos - Mgrid; + dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR); + N++; + } + if (N < 1) continue; + + liststats (list, dlist, N, &stats); + if (mark) catalog[i].found[j] = TRUE; + + /* use sigma or error in dM for output? */ + catalog[i].secfilt[Nsecfilt*j+Nsec].M = stats.mean; + catalog[i].secfilt[Nsecfilt*j+Nsec].dM = MAX (stats.error, stats.sigma); + catalog[i].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT; + } } } @@ -179,5 +221,7 @@ } -// for each average object, set the average mags based on existing equiv photometry +// For each average object, set the average mags based on existing equiv photometry. +// NOTE: this function operates on the real Measure & Average structures, not the +// MeasureTiny & AverageTiny structures int setMave (Catalog *catalog, int Ncatalog) { @@ -202,5 +246,5 @@ ALLOCATE (dlist, double, MAX (1, Nmax)); - Nsecfilt = catalog[0].Nsecfilt; + Nsecfilt = GetPhotcodeNsecfilt (); # define PSFQUALSTATS 1 @@ -220,5 +264,5 @@ if (GetPhotcodeEquivCodebyCode (catalog[i].measure[m].photcode) != Nc) continue; - Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]); + Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); if (isnan(Msys)) continue; @@ -321,7 +365,7 @@ for (j = 0; j < catalog[i].Naverage; j++) { - m = catalog[i].average[j].measureOffset; - for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) { - if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue; + m = catalog[i].averageT[j].measureOffset; + for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) { + if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; Mcal = getMcal (m, i); if (isnan(Mcal)) continue; @@ -330,4 +374,6 @@ Mgrid = getMgrid (m, i); if (isnan(Mgrid)) continue; + + // set the output calibration catalog[i].measure[m].Mcal = Mcal + Mmos + Mgrid; } @@ -339,5 +385,5 @@ void clean_stars (Catalog *catalog, int Ncatalog) { - int i, j, Ndel, Nave, Ntot, mark; + int i, j, Ndel, Nave, Ntot, mark, Ns; float dM, Xm; double Chisq, MaxScatter, MaxChisq; @@ -354,42 +400,53 @@ ALLOCATE (slist, double, Ntot); ALLOCATE (dlist, double, Ntot); - for (i = Ntot = 0; i < Ncatalog; i++) { - for (j = 0; j < catalog[i].Naverage; j++) { - if (catalog[i].average[j].flags & STAR_BAD) continue; - Xm = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm; - if (Xm == -1) continue; - Chisq = pow (10.0, 0.01*Xm); - xlist[Ntot] = Chisq; - slist[Ntot] = catalog[i].secfilt[PhotNsec*j+PhotSec].dM; - dlist[Ntot] = 1; - Ntot ++; - } - } + + int Nsecfilt = GetPhotcodeNsecfilt (); + + // eliminate bad stars using the stats for a single secfilt at a time + // XXX DEP replace average.flags with secfilt flags + for (Ns = 0; Ns < Nphotcodes; Ns ++) { + + int thisCode = photcodes[Ns][0].code; + int Nsec = GetPhotcodeNsec(thisCode); + + for (i = Ntot = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + if ( catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD ) continue; + Xm = catalog[i].secfilt[Nsecfilt*j+Nsec].Xm; + if (Xm == -1) continue; + Chisq = pow (10.0, 0.01*Xm); + xlist[Ntot] = Chisq; + slist[Ntot] = catalog[i].secfilt[Nsecfilt*j+Nsec].dM; + dlist[Ntot] = 1; + Ntot ++; + } + } - initstats ("MEAN"); - liststats (xlist, dlist, Ntot, &stats); - MaxChisq = MAX (STAR_CHISQ, 2*stats.median); - liststats (slist, dlist, Ntot, &stats); - MaxScatter = MAX (STAR_SCATTER, 2*stats.median); - fprintf (stderr, "Max Scatter: %f, Max Chisq: %f\n", MaxScatter, MaxChisq); - - Ndel = Nave = 0; - for (i = 0; i < Ncatalog; i++) { - for (j = 0; j < catalog[i].Naverage; j++) { - dM = catalog[i].secfilt[PhotNsec*j+PhotSec].dM; - Xm = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm; - Chisq = pow (10.0, 0.01*Xm); - mark = (dM > MaxScatter) || (Xm == NAN_S_SHORT) || (Chisq > MaxChisq); - if (mark) { - catalog[i].average[j].flags |= ID_STAR_POOR; - Ndel ++; - } else { - catalog[i].average[j].flags &= ~ID_STAR_POOR; - } - Nave ++; - } - } - fprintf (stderr, "%d stars marked variable, %d total\n", Ndel, Nave); - initstats (STATMODE); + initstats ("MEAN"); + liststats (xlist, dlist, Ntot, &stats); + MaxChisq = MAX (STAR_CHISQ, 2*stats.median); + liststats (slist, dlist, Ntot, &stats); + MaxScatter = MAX (STAR_SCATTER, 2*stats.median); + fprintf (stderr, "Max Scatter: %f, Max Chisq: %f\n", MaxScatter, MaxChisq); + + Ndel = Nave = 0; + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dM; + Xm = catalog[i].secfilt[Nsecfilt*j+Nsec].Xm; + Chisq = pow (10.0, 0.01*Xm); + mark = (dM > MaxScatter) || (Xm == NAN_S_SHORT) || (Chisq > MaxChisq); + if (mark) { + catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_POOR; + Ndel ++; + } else { + catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_POOR; + } + Nave ++; + } + } + fprintf (stderr, "%d stars marked variable, %d total\n", Ndel, Nave); + initstats (STATMODE); + } free (xlist); free (slist); @@ -409,4 +466,6 @@ int Ncal, Nmos, Ngrid, Nfew; + int Nsecfilt = GetPhotcodeNsecfilt (); + if (VERBOSE) fprintf (stderr, "marking poor measures\n"); /* Nmeasure is now different, need to reallocate */ @@ -414,5 +473,5 @@ for (i = 0; i < Ncatalog; i++) { for (j = 0; j < catalog[i].Naverage; j++) { - Nmax = MAX (Nmax, catalog[i].average[j].Nmeasure); + Nmax = MAX (Nmax, catalog[i].averageT[j].Nmeasure); } } @@ -430,79 +489,98 @@ for (j = 0; j < catalog[i].Naverage; j++) { - /* skip bad stars to prevent them from becoming good (on inner sample) */ - if (catalog[i].average[j].flags & STAR_BAD) continue; - - /* on final processing, skip stars already measured */ - if (final && catalog[i].found[j]) continue; - - /* accumulate list of valid measurements */ - m = catalog[i].average[j].measureOffset; - N = 0; - for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) { - /* if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue; */ - Mcal = getMcal (m, i); - if (isnan(Mcal)) { Ncal ++; continue; } - Mmos = getMmos (m, i); - if (isnan(Mmos)) { Nmos ++; continue; } - Mgrid = getMgrid (m, i); - if (isnan(Mgrid)) { Ngrid ++; continue; } - - Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]); - list[N] = Msys - Mcal - Mmos - Mgrid; - dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR); - N++; - } - if (N <= TOOFEW) { Nfew ++; continue; } - - /* 3-sigma clip based on stats of inner 50% */ - - // calculated mean of inner 50% - initstats ("INNER_MEAN"); - liststats (list, dlist, N, &stats); - stats.sigma = MAX (MIN_ERROR, stats.sigma); /* if measurements agree too well, sigma -> 0.0 */ - - // ignore entries > 3sigma from inner mean - for (k = m = 0; k < N; k++) { - if (fabs (list[k] - stats.median) < NSIGMA_CLIP*stats.sigma) { - list[m] = list[k]; - m++; - } - } - // recalculate the mean & sigma of the accepted measurements - initstats ("MEAN"); - liststats (list, dlist, m, &stats); - stats.sigma = MAX (MIN_ERROR, stats.sigma); - - /* apply to list of all relevant measurements, including IMAGE_POOR & IMAGE_FEW */ - image_bad = IMAGE_BAD; - IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL; - m = catalog[i].average[j].measureOffset; - N = 0; - for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) { - /* if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue; */ - Mcal = getMcal (m, i); - if (isnan(Mcal)) continue; - Mmos = getMmos (m, i); - if (isnan(Mmos)) continue; - Mgrid = getMgrid (m, i); - if (isnan(Mgrid)) continue; - - Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]); - list[N] = Msys - Mcal - Mmos - Mgrid; - dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR); - ilist[N] = m; - N++; - Nave ++; - } - if (N < TOOFEW) continue; - - /* mark bad measures (> 3 sigma deviant) */ - for (k = 0; k < N; k++) { - if (fabs (list[k] - stats.median) > NSIGMA_REJECT*stats.sigma) { - catalog[i].measure[ilist[k]].dbFlags |= ID_MEAS_POOR_PHOTOM; - Ndel ++; - } - } - IMAGE_BAD = image_bad; + int Ns; + for (Ns = 0; Ns < Nphotcodes; Ns++) { + + /* on final processing, skip stars already measured */ + if (final && catalog[i].found[j]) continue; + + int thisCode = photcodes[Ns][0].code; + int Nsec = GetPhotcodeNsec(thisCode); + + /* skip bad stars to prevent them from becoming good (on inner sample) */ + if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; + + /* accumulate list of valid measurements */ + m = catalog[i].averageT[j].measureOffset; + N = 0; + for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) { + // skip measurements that do not match the current photcode + int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode); + if (ecode != thisCode) { continue; } + + /* if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; */ + Mcal = getMcal (m, i); + if (isnan(Mcal)) { Ncal ++; continue; } + Mmos = getMmos (m, i); + if (isnan(Mmos)) { Nmos ++; 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]); + list[N] = Msys - Mcal - Mmos - Mgrid; + dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR); + N++; + } + if (N <= TOOFEW) { Nfew ++; continue; } + + /* 3-sigma clip based on stats of inner 50% */ + + // calculated mean of inner 50% + initstats ("INNER_MEAN"); + liststats (list, dlist, N, &stats); + stats.sigma = MAX (MIN_ERROR, stats.sigma); /* if measurements agree too well, sigma -> 0.0 */ + + // ignore entries > 3sigma from inner mean + for (k = m = 0; k < N; k++) { + if (fabs (list[k] - stats.median) < NSIGMA_CLIP*stats.sigma) { + list[m] = list[k]; + m++; + } + } + // recalculate the mean & sigma of the accepted measurements + initstats ("MEAN"); + liststats (list, dlist, m, &stats); + stats.sigma = MAX (MIN_ERROR, stats.sigma); + + /* apply to list of all relevant measurements, including IMAGE_POOR & IMAGE_FEW */ + image_bad = IMAGE_BAD; + IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL; + m = catalog[i].averageT[j].measureOffset; + N = 0; + for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) { + // skip measurements that do not match the current photcode + int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode); + if (ecode != thisCode) { continue; } + + /* if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; */ + Mcal = getMcal (m, i); + if (isnan(Mcal)) continue; + Mmos = getMmos (m, i); + if (isnan(Mmos)) continue; + Mgrid = getMgrid (m, i); + if (isnan(Mgrid)) continue; + + Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt]); + list[N] = Msys - Mcal - Mmos - Mgrid; + dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR); + ilist[N] = m; + N++; + Nave ++; + } + if (N < TOOFEW) continue; + + /* mark bad measures (> 3 sigma deviant) */ + for (k = 0; k < N; k++) { + if (fabs (list[k] - stats.median) > NSIGMA_REJECT*stats.sigma) { + catalog[i].measureT[ilist[k]].dbFlags |= ID_MEAS_POOR_PHOTOM; + if (final) { + // for the final pass, we have a duplicate set of values in measure and measureT + catalog[i].measure[ilist[k]].dbFlags |= ID_MEAS_POOR_PHOTOM; + } + Ndel ++; + } + } + IMAGE_BAD = image_bad; + } } } @@ -513,5 +591,5 @@ } -StatType statsStarN (Catalog *catalog, int Ncatalog) { +StatType statsStarN (Catalog *catalog, int Ncatalog, int Nsec, int seccode) { off_t j, k, m, Ntot; @@ -520,4 +598,8 @@ float Mcal, Mmos, Mgrid; StatType stats; + // int N1, N2, N3, N4, N0; + // N1 = N2 = N3 = N4 = N0 = 0; + + int Nsecfilt = GetPhotcodeNsecfilt (); Ntot = 0; @@ -534,15 +616,17 @@ /* calculate the average value for a single star */ - if (catalog[i].average[j].flags & STAR_BAD) continue; - m = catalog[i].average[j].measureOffset; + if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) { continue; } + m = catalog[i].averageT[j].measureOffset; N = 0; - for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) { + for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) { + int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode); + if (ecode != seccode) { continue;} Mcal = getMcal (m, i); - if (isnan(Mcal)) continue; + if (isnan(Mcal)) { continue;} Mmos = getMmos (m, i); - if (isnan(Mmos)) continue; + if (isnan(Mmos)) { continue; } Mgrid = getMgrid (m, i); - if (isnan(Mgrid)) continue; + if (isnan(Mgrid)) { continue;} N++; } @@ -554,4 +638,5 @@ } + // fprintf (stderr, "N1: %d, N2: %d, N3: %d, N4: %d, N0: %d\n", N1, N2, N3, N4, N0); liststats (list, dlist, n, &stats); free (list); @@ -560,5 +645,6 @@ } -StatType statsStarX (Catalog *catalog, int Ncatalog) { +// stats for a single secfilt at a time +StatType statsStarX (Catalog *catalog, int Ncatalog, int Nsec) { off_t j, Ntot; @@ -567,4 +653,6 @@ StatType stats; + int Nsecfilt = GetPhotcodeNsecfilt (); + Ntot = 0; for (i = 0; i < Ncatalog; i++) { @@ -580,7 +668,7 @@ /* calculate the average value for a single star */ - if (catalog[i].average[j].flags & STAR_BAD) continue; - - Xm = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm; + if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; + + Xm = catalog[i].secfilt[Nsecfilt*j+Nsec].Xm; if (Xm == NAN_S_SHORT) continue; list[n] = pow (10.0, 0.01*Xm); @@ -596,5 +684,5 @@ } -StatType statsStarS (Catalog *catalog, int Ncatalog) { +StatType statsStarS (Catalog *catalog, int Ncatalog, int Nsec) { int i, n; @@ -604,4 +692,6 @@ StatType stats; + int Nsecfilt = GetPhotcodeNsecfilt (); + Ntot = 0; for (i = 0; i < Ncatalog; i++) { @@ -617,7 +707,7 @@ /* calculate the average value for a single star */ - if (catalog[i].average[j].flags & STAR_BAD) continue; - - dM = catalog[i].secfilt[PhotNsec*j+PhotSec].dM; + if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; + + dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dM; list[n] = dM; dlist[n] = 1; @@ -644,20 +734,28 @@ ALLOCATE (Mlist, double, NBIN); - for (i = 0; i < NBIN; i++) xlist[i] = 0.00025*i; - bzero (Mlist, NBIN*sizeof(double)); - for (i = 0; i < Ncatalog; i++) { - for (j = 0; j < catalog[i].Naverage; j++) { - if (catalog[i].average[j].flags & STAR_BAD) continue; - dMrel = catalog[i].secfilt[PhotNsec*j+PhotSec].dM; - bin = dMrel / 0.00025; - bin = MAX (0, MIN (NBIN-1, bin)); - Mlist[bin] += 1.0; - } - } - - plot_defaults (&graphdata); - graphdata.style = 1; - plot_list (&graphdata, xlist, Mlist, NBIN, "dMrel hist", "%s.dMhist.png", OUTROOT); - + int Nsecfilt = GetPhotcodeNsecfilt (); + + int Ns; + for (Ns = 0; Ns < Nphotcodes; Ns++) { + + int thisCode = photcodes[Ns][0].code; + int Nsec = GetPhotcodeNsec(thisCode); + + for (i = 0; i < NBIN; i++) xlist[i] = 0.00025*i; + bzero (Mlist, NBIN*sizeof(double)); + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; + dMrel = catalog[i].secfilt[Nsecfilt*j+Nsec].dM; + bin = dMrel / 0.00025; + bin = MAX (0, MIN (NBIN-1, bin)); + Mlist[bin] += 1.0; + } + } + + plot_defaults (&graphdata); + graphdata.style = 1; + plot_list (&graphdata, xlist, Mlist, NBIN, "dMrel hist", "%s.dMhist.png", OUTROOT); + } free (xlist); free (Mlist); @@ -671,4 +769,6 @@ Graphdata graphdata; + int Nsecfilt = GetPhotcodeNsecfilt (); + Ntotal = 0; for (i = 0; i < Ncatalog; i++) Ntotal += catalog[i].Naverage; @@ -677,19 +777,27 @@ ALLOCATE (ylist, double, Ntotal); - N = 0; - for (i = 0; i < Ncatalog; i++) { - for (j = 0; j < catalog[i].Naverage; j++) { - if (catalog[i].average[j].flags & STAR_BAD) continue; - xlist[N] = catalog[i].secfilt[PhotNsec*j+PhotSec].M; - value = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm; - if (value == NAN_S_SHORT) continue; - ylist[N] = 0.01*value; - N++; - } - } - - plot_defaults (&graphdata); - graphdata.ymin = -3.0; - plot_list (&graphdata, xlist, ylist, N, "chisq", "%s.chisq.png", OUTROOT); + int Ns; + for (Ns = 0; Ns < Nphotcodes; Ns++) { + + int thisCode = photcodes[Ns][0].code; + int Nsec = GetPhotcodeNsec(thisCode); + + N = 0; + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; + xlist[N] = catalog[i].secfilt[Nsecfilt*j+Nsec].M; + value = catalog[i].secfilt[Nsecfilt*j+Nsec].Xm; + if (value == NAN_S_SHORT) continue; + ylist[N] = 0.01*value; + N++; + } + } + + plot_defaults (&graphdata); + graphdata.ymin = -3.0; + plot_list (&graphdata, xlist, ylist, N, "chisq", "%s.chisq.png", OUTROOT); + } + free (xlist); free (ylist); @@ -716,6 +824,6 @@ for (i = 0; i < Ncatalog; i++) { for (j = 0; j < catalog[i].Naverage; j++) { - xlist[N] = catalog[i].average[j].R; - ylist[N] = catalog[i].average[j].D; + xlist[N] = catalog[i].averageT[j].R; + ylist[N] = catalog[i].averageT[j].D; N++; } Index: trunk/Ohana/src/relphot/src/args.c =================================================================== --- trunk/Ohana/src/relphot/src/args.c (revision 31425) +++ trunk/Ohana/src/relphot/src/args.c (revision 31450) @@ -128,4 +128,12 @@ } + MaxDensityUse = FALSE; + if ((N = get_argument (argc, argv, "-max-density"))) { + remove_argument (N, &argc, argv); + MaxDensityValue = atof(argv[N]); + remove_argument (N, &argc, argv); + MaxDensityUse = TRUE; + } + SHOW_PARAMS = FALSE; if ((N = get_argument (argc, argv, "-params"))) { Index: trunk/Ohana/src/relphot/src/bcatalog.c =================================================================== --- trunk/Ohana/src/relphot/src/bcatalog.c (revision 31425) +++ trunk/Ohana/src/relphot/src/bcatalog.c (revision 31450) @@ -1,16 +1,49 @@ # include "relphot.h" + +extern double drand48(); + +int CopyAverageTiny (AverageTiny *averageT, Average *average) { + + averageT[0].R = average[0].R; + averageT[0].D = average[0].D; + averageT[0].flags = average[0].flags; + averageT[0].Nmeasure = average[0].Nmeasure; + averageT[0].measureOffset = average[0].measureOffset; + + // make Nmeasure & measureOffset optional? + + return (TRUE); +} + +int CopyMeasureTiny (MeasureTiny *measureT, Measure *measure) { + + measureT[0].dR = measure[0].dR; + measureT[0].dD = measure[0].dD; + measureT[0].M = measure[0].M; + measureT[0].Mcal = measure[0].Mcal; + measureT[0].dM = measure[0].dM; + measureT[0].airmass = measure[0].airmass; + measureT[0].Xccd = measure[0].Xccd; + measureT[0].Yccd = measure[0].Yccd; + measureT[0].t = measure[0].t; + measureT[0].dt = measure[0].dt; + measureT[0].averef = measure[0].averef; + measureT[0].imageID = measure[0].imageID; + measureT[0].dbFlags = measure[0].dbFlags; + measureT[0].photcode = measure[0].photcode; + + return (TRUE); +} int bcatalog (Catalog *subcatalog, Catalog *catalog) { off_t i, j, offset; - int ecode; + int ecode, found, Ns; off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm; float mag; int Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Ngalaxy, Npsfqf; - // XXX PhotNsec as a global is a bad idea; either get it from catalog - // or get it from: - // Nsecfilt = GetPhotcodeNsecfilt (); - // assert (catalog[0].Nsecfilt == Nsecfilt); + int Nsecfilt = GetPhotcodeNsecfilt (); + assert (Nsecfilt == catalog[0].Nsecfilt); /* we are moving only the subset of measurements from catalog[0] to subcatalog[0] */ @@ -18,7 +51,7 @@ NAVERAGE = 50; NMEASURE = 1000; - ALLOCATE (subcatalog[0].average, Average, NAVERAGE); - ALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*PhotNsec); - ALLOCATE (subcatalog[0].measure, Measure, NMEASURE); + ALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE); + ALLOCATE (subcatalog[0].averageT, AverageTiny, NAVERAGE); + ALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt); Nmeasure = Naverage = 0; @@ -30,15 +63,23 @@ /* start with all stars good */ - subcatalog[0].average[Naverage] = catalog[0].average[i]; - subcatalog[0].average[Naverage].measureOffset = Nmeasure; - for (j = 0; j < PhotNsec; j++) { - subcatalog[0].secfilt[PhotNsec*Naverage+j] = catalog[0].secfilt[PhotNsec*i+j]; + CopyAverageTiny (&subcatalog[0].averageT[Naverage], &catalog[0].average[i]); + subcatalog[0].averageT[Naverage].measureOffset = Nmeasure; + + for (j = 0; j < Nsecfilt; j++) { + subcatalog[0].secfilt[Nsecfilt*Naverage+j] = catalog[0].secfilt[Nsecfilt*i+j]; } if (RESET) { - subcatalog[0].secfilt[PhotNsec*Naverage+PhotSec].M = NAN; - subcatalog[0].secfilt[PhotNsec*Naverage+PhotSec].dM = NAN; - subcatalog[0].average[Naverage].flags &= ~ID_STAR_FEW; - subcatalog[0].average[Naverage].flags &= ~ID_STAR_POOR; + int Ns; + for (Ns = 0; Ns < Nphotcodes; Ns++) { + + int thisCode = photcodes[Ns][0].code; + int Nsec = GetPhotcodeNsec(thisCode); + + subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M = NAN; + subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].dM = NAN; + subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].flags &= ~ID_STAR_FEW; + subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].flags &= ~ID_STAR_POOR; + } } @@ -52,5 +93,12 @@ /* select measurements by photcode */ ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[offset].photcode); - if (ecode != photcode[0].code) { Ncode ++; continue; } + found = FALSE; + for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) { + if (ecode == photcodes[Ns][0].code) found = TRUE; + } + if (!found) { + Ncode ++; + continue; + } /* select measurements by time */ @@ -61,6 +109,4 @@ /* select measurements by quality */ - // XXX ignore this criterion for REF measurements? - // XXX chnage this to select by bitflags if (DophotSelect && ((catalog[0].measure[offset].photFlags >> 16) != DophotValue)) { Ndophot ++; continue; } @@ -70,9 +116,9 @@ // check for galaxies if (!isnan(catalog[0].measure[offset].Map)) { - if (catalog[0].measure[offset].M - catalog[0].measure[offset].Map > 0.15) { - nEXT ++; - } else { - nPSF ++; - } + if (catalog[0].measure[offset].M - catalog[0].measure[offset].Map > 0.15) { + nEXT ++; + } else { + nPSF ++; + } } @@ -91,13 +137,13 @@ } - subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_SKIP_PHOTOM; - subcatalog[0].measure[Nmeasure] = catalog[0].measure[offset]; - subcatalog[0].measure[Nmeasure].averef = Naverage; + CopyMeasureTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]); + subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_SKIP_PHOTOM; + subcatalog[0].measureT[Nmeasure].averef = Naverage; if (RESET) { - subcatalog[0].measure[Nmeasure].Mcal = 0; - subcatalog[0].measure[Nmeasure].dbFlags &= 0xff00; - subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_POOR_PHOTOM; - subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_AREA; - subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_NOCAL; + subcatalog[0].measureT[Nmeasure].Mcal = 0; + subcatalog[0].measureT[Nmeasure].dbFlags &= 0xff00; + subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_POOR_PHOTOM; + subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_AREA; + subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_NOCAL; } Nmeasure ++; @@ -105,5 +151,5 @@ if (Nmeasure == NMEASURE) { NMEASURE += 1000; - REALLOCATE (subcatalog[0].measure, Measure, NMEASURE); + REALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE); } } @@ -122,20 +168,24 @@ continue; } - subcatalog[0].average[Naverage].Nmeasure = Nm; + subcatalog[0].averageT[Naverage].Nmeasure = Nm; Naverage ++; if (Naverage == NAVERAGE) { NAVERAGE += 50; - REALLOCATE (subcatalog[0].average, Average, NAVERAGE); - REALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*PhotNsec); - } - } - REALLOCATE (subcatalog[0].average, Average, MAX (Naverage, 1)); - REALLOCATE (subcatalog[0].measure, Measure, MAX (Nmeasure, 1)); - REALLOCATE (subcatalog[0].secfilt, SecFilt, PhotNsec*MAX (Naverage, 1)); + REALLOCATE (subcatalog[0].averageT, AverageTiny, NAVERAGE); + REALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt); + } + } + REALLOCATE (subcatalog[0].averageT, AverageTiny, MAX (Naverage, 1)); + REALLOCATE (subcatalog[0].measureT, MeasureTiny, MAX (Nmeasure, 1)); + REALLOCATE (subcatalog[0].secfilt, SecFilt, Nsecfilt*MAX (Naverage, 1)); subcatalog[0].Naverage = Naverage; subcatalog[0].Nmeasure = Nmeasure; subcatalog[0].Nsecfilt = catalog[0].Nsecfilt; subcatalog[0].Nsecf_mem = Naverage * catalog[0].Nsecfilt; - assert (PhotNsec == catalog[0].Nsecfilt); + + // limit the total number of stars in the catalog + if (MaxDensityUse) { + LimitDensityCatalog (subcatalog, catalog); + } if (VERBOSE) { @@ -147,2 +197,119 @@ return (TRUE); } + +int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog) { + + Catalog tmpcatalog; + + double Rmin, Rmax, Dmin, Dmax; + + int Nsecfilt = GetPhotcodeNsecfilt (); + + gfits_scan (&catalog[0].header, "RA0", "%lf", 1, &Rmin); + gfits_scan (&catalog[0].header, "DEC0", "%lf", 1, &Dmin); + gfits_scan (&catalog[0].header, "RA1", "%lf", 1, &Rmax); + gfits_scan (&catalog[0].header, "DEC1", "%lf", 1, &Dmax); + + if (VERBOSE2) fprintf (stderr, "extracting from catalog covering region %f,%f to %f,%f\n", Rmin, Dmin, Rmax, Dmax); + + float AREA = fabs(Dmax - Dmin) * fabs(Rmax - Rmin) * cos (0.5*RAD_DEG*(Dmax + Dmin)); + assert (AREA > 0); + + off_t Nmax = MaxDensityValue * AREA; + if (subcatalog[0].Naverage <= Nmax) { + if (VERBOSE) { + fprintf (stderr, "subcatalog has less than the max density\n"); + } + return (TRUE); + } + + off_t Naverage = subcatalog[0].Naverage; + + // select a random subset of Nmax stars from subcatalog using Fisher-Yates + + // we are going to select Nmax entries by generating a random-sorted index list + off_t *index, tmp, i, j, ave; + ALLOCATE (index, off_t, Naverage); + for (i = 0; i < Naverage; i++) { + index[i] = i; + } + for (i = 0; i < Naverage; i++) { + j = (Naverage - i) * drand48() + i; // a number between i and Naverage + tmp = index[j]; + index[j] = index[i]; + index[i] = tmp; + } + + // count the number of measurements this selection will yield + off_t NMEASURE = 0; + for (i = 0; i < Nmax; i++) { + ave = index[i]; + NMEASURE += subcatalog[0].averageT[ave].Nmeasure; + } + + // allocate the output data + ALLOCATE (tmpcatalog.averageT, AverageTiny, Nmax); + ALLOCATE (tmpcatalog.measureT, MeasureTiny, NMEASURE); + ALLOCATE (tmpcatalog.secfilt, SecFilt, Nmax * Nsecfilt); + + off_t Nmeasure = 0; + + // copy the Nmax selected entries from subcatalog to tmpcatalog (adjusting links) + for (i = 0; i < Nmax; i++) { + ave = index[i]; + tmpcatalog.averageT[i] = subcatalog[0].averageT[ave]; + tmpcatalog.averageT[i].measureOffset = Nmeasure; + for (j = 0; j < tmpcatalog.averageT[i].Nmeasure; j++) { + off_t offset = subcatalog[0].averageT[ave].measureOffset + j; + tmpcatalog.measureT[Nmeasure] = subcatalog[0].measureT[offset]; + tmpcatalog.measureT[Nmeasure].averef = i; + Nmeasure ++; + } + } + + if (VERBOSE) { + fprintf (stderr, "limited to "OFF_T_FMT" of "OFF_T_FMT" stars ("OFF_T_FMT" of "OFF_T_FMT" measures) for catalog %s\n", + Nmax, subcatalog[0].Naverage, Nmeasure, subcatalog[0].Nmeasure, catalog[0].filename); + } + + free (subcatalog[0].averageT); + free (subcatalog[0].measureT); + free (subcatalog[0].secfilt); + + subcatalog[0].averageT = tmpcatalog.averageT; + subcatalog[0].measureT = tmpcatalog.measureT; + subcatalog[0].secfilt = tmpcatalog.secfilt; + subcatalog[0].Naverage = Nmax; + subcatalog[0].Nmeasure = Nmeasure; + subcatalog[0].Nsecfilt = catalog[0].Nsecfilt; + subcatalog[0].Nsecf_mem = Naverage * catalog[0].Nsecfilt; + + return (TRUE); +} + +// for the cases where we are not using a subset of the data, we still need to have a copy of these fields +int populate_tiny_values (Catalog *catalog) { + + off_t i; + + ALLOCATE (catalog[0].measureT, MeasureTiny, catalog[0].Nmeasure); + ALLOCATE (catalog[0].averageT, AverageTiny, catalog[0].Naverage); + + for (i = 0; i < catalog[0].Naverage; i++) { + CopyAverageTiny (&catalog[0].averageT[i], &catalog[0].average[i]); + } + + for (i = 0; i < catalog[0].Nmeasure; i++) { + CopyMeasureTiny (&catalog[0].measureT[i], &catalog[0].measure[i]); + } + + return (TRUE); +} + +int free_tiny_values (Catalog *catalog) { + + free (catalog[0].averageT); + free (catalog[0].measureT); + return (TRUE); +} + Index: trunk/Ohana/src/relphot/src/global_stats.c =================================================================== --- trunk/Ohana/src/relphot/src/global_stats.c (revision 31425) +++ trunk/Ohana/src/relphot/src/global_stats.c (revision 31450) @@ -7,34 +7,44 @@ initstats ("MEAN"); - stN = statsStarN (catalog, Ncatalog); - stX = statsStarX (catalog, Ncatalog); - stS = statsStarS (catalog, Ncatalog); + fprintf (stderr, "\n"); + fprintf (stderr, "STATS median mean sigma min max Nmeas\n"); + int Ns; + for (Ns = 0; Ns < Nphotcodes; Ns++) { + + int thisCode = photcodes[Ns][0].code; + int Nsec = GetPhotcodeNsec(thisCode); + int seccode = photcodes[Ns][0].code; + + stN = statsStarN (catalog, Ncatalog, Nsec, seccode); + stX = statsStarX (catalog, Ncatalog, Nsec); + stS = statsStarS (catalog, Ncatalog, Nsec); + + fprintf (stderr, " --- stats for %s ---\n", photcodes[Ns][0].name); + fprintf (stderr, "meas / star: %7.0f %7.1f %7.1f %7.0f %7.0f %6d\n", stN.median, stN.mean, stN.sigma, stN.min, stN.max, stN.Nmeas); + fprintf (stderr, "dMrel star: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n", stS.median, stS.mean, stS.sigma, stS.min, stS.max, stS.Nmeas); + fprintf (stderr, "chisq star: %7.1f %7.1f %7.1f %7.1f %7.1f %6d\n", stX.median, stX.mean, stX.sigma, stX.min, stX.max, stX.Nmeas); + } + imN = statsImageN (catalog); imX = statsImageX (catalog); imM = statsImageM (catalog); imD = statsImagedM (catalog); - + msN = statsMosaicN (catalog); msM = statsMosaicM (catalog); msD = statsMosaicdM (catalog); msX = statsMosaicX (catalog); - - fprintf (stderr, "STATS median mean sigma min max Nmeas\n"); + fprintf (stderr, "meas / image: %7.0f %7.0f %7.0f %7.0f %7.0f %6d\n", imN.median, imN.mean, imN.sigma, imN.min, imN.max, imN.Nmeas); - fprintf (stderr, "meas / mosaic: %7.0f %7.0f %7.0f %7.0f %7.0f %6d\n", msN.median, msN.mean, msN.sigma, msN.min, msN.max, msN.Nmeas); - fprintf (stderr, "meas / star: %7.0f %7.1f %7.1f %7.0f %7.0f %6d\n", stN.median, stN.mean, stN.sigma, stN.min, stN.max, stN.Nmeas); - - fprintf (stderr, "chisq image: %7.1f %7.1f %7.1f %7.1f %7.1f %6d\n", imX.median, imX.mean, imX.sigma, imX.min, imX.max, imX.Nmeas); - fprintf (stderr, "chisq mosaic: %7.1f %7.1f %7.1f %7.1f %7.1f %6d\n", msX.median, msX.mean, msX.sigma, msX.min, msX.max, msX.Nmeas); - fprintf (stderr, "chisq star: %7.1f %7.1f %7.1f %7.1f %7.1f %6d\n", stX.median, stX.mean, stX.sigma, stX.min, stX.max, stX.Nmeas); - fprintf (stderr, "Mcal image: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n", imM.median, imM.mean, imM.sigma, imM.min, imM.max, imM.Nmeas); fprintf (stderr, "dMcal image: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n", imD.median, imD.mean, imD.sigma, imD.min, imD.max, imD.Nmeas); + fprintf (stderr, "chisq image: %7.1f %7.1f %7.1f %7.1f %7.1f %6d\n", imX.median, imX.mean, imX.sigma, imX.min, imX.max, imX.Nmeas); + fprintf (stderr, "meas / mosaic: %7.0f %7.0f %7.0f %7.0f %7.0f %6d\n", msN.median, msN.mean, msN.sigma, msN.min, msN.max, msN.Nmeas); fprintf (stderr, "Mcal mosaic: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n", msM.median, msM.mean, msM.sigma, msM.min, msM.max, msM.Nmeas); fprintf (stderr, "dMcal mosaic: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n", msD.median, msD.mean, msD.sigma, msD.min, msD.max, msD.Nmeas); - - fprintf (stderr, "dMrel star: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n\n", stS.median, stS.mean, stS.sigma, stS.min, stS.max, stS.Nmeas); + fprintf (stderr, "chisq mosaic: %7.1f %7.1f %7.1f %7.1f %7.1f %6d\n", msX.median, msX.mean, msX.sigma, msX.min, msX.max, msX.Nmeas); + initstats (STATMODE); Index: trunk/Ohana/src/relphot/src/initialize.c =================================================================== --- trunk/Ohana/src/relphot/src/initialize.c (revision 31425) +++ trunk/Ohana/src/relphot/src/initialize.c (revision 31450) @@ -13,4 +13,6 @@ N = UserPatchSelect ? 1 : 2; +# if (0) + // XXX DEP if ((photcode = GetPhotcodebyName (argv[N])) == NULL) { fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", argv[N]); @@ -21,8 +23,39 @@ exit (1); } + // PhotSec is used to select the single average photcode being processed PhotSec = GetPhotcodeNsec (photcode[0].code); +# endif + + Nphotcodes = 0; + photcodes = NULL; + int NPHOTCODES = 10; + ALLOCATE (photcodes, PhotCode *, NPHOTCODES); + + /* parse the comma-separated list of photcodesKeep */ + char *myList = strcreate(argv[N]); + char *list = myList; + char *codename = NULL; + char *ptr = NULL; + while ((codename = strtok_r (list, ",", &ptr)) != NULL) { + list = NULL; // pass NULL on successive strtok_r calls + fprintf (stderr, "PHOTCODE LIST: %s\n", myList); + fprintf (stderr, "codename: %s\n", codename); + if ((photcodes[Nphotcodes] = GetPhotcodebyName (codename)) == NULL) { + fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename); + exit (1); + } + if (photcodes[Nphotcodes][0].type != PHOT_SEC) { + fprintf (stderr, "photcode %s is not an filter type (SEC)\n", codename); + exit (1); + } + Nphotcodes ++; + CHECK_REALLOCATE (photcodes, PhotCode *, NPHOTCODES, Nphotcodes, 10); + } } - - PhotNsec = GetPhotcodeNsecfilt (); + // XXX DEP PhotNsec = GetPhotcodeNsecfilt (); + if (USE_GRID && (Nphotcodes > 1)) { + fprintf (stderr, "grid correction analysis currently can only operate on a single photcode\n"); + exit (1); + } initstats (STATMODE); @@ -55,4 +88,9 @@ exit (0); } + + // init the random seed + long A, B; + A = time(NULL); + for (B = 0; A == time(NULL); B++); + srand48(B); } - Index: trunk/Ohana/src/relphot/src/load_images.c =================================================================== --- trunk/Ohana/src/relphot/src/load_images.c (revision 31425) +++ trunk/Ohana/src/relphot/src/load_images.c (revision 31450) @@ -44,5 +44,5 @@ // select the images which overlap the selected sky regions subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset); - MARKTIME("selected images: %f sec\n", dtime); + MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime); // generate db->vtable from db->ftable based on the selection Index: trunk/Ohana/src/relphot/src/plot_scatter.c =================================================================== --- trunk/Ohana/src/relphot/src/plot_scatter.c (revision 31425) +++ trunk/Ohana/src/relphot/src/plot_scatter.c (revision 31450) @@ -18,41 +18,54 @@ ALLOCATE (ilist, double, Ntot); - N = 0; - for (i = 0; i < Ncatalog; i++) { - for (j = 0; j < catalog[i].Naverage; j++) { + int Nsecfilt = GetPhotcodeNsecfilt (); - /* calculate the average value for a single star */ - if (catalog[i].average[j].flags & STAR_BAD) continue; - m = catalog[i].average[j].measureOffset; + int Ns; + for (Ns = 0; Ns < Nphotcodes; Ns++) { - for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) { - if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue; - Mcal = getMcal (m, i); - if (isnan(Mcal)) continue; - Mmos = getMmos (m, i); - if (isnan(Mmos)) continue; - Mgrid = getMgrid (m, i); - if (isnan(Mgrid)) continue; + int thisCode = photcodes[Ns][0].code; + int Nsec = GetPhotcodeNsec(thisCode); - Mrel = catalog[i].secfilt[PhotNsec*j+PhotSec].M; - xlist[N] = Mrel; - ylist[N] = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]) - Mcal - Mmos - Mgrid - Mrel; - ilist[N] = PhotInst (&catalog[i].measure[m]); - N++; - } + N = 0; + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + + /* calculate the average value for a single star */ + if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; + m = catalog[i].average[j].measureOffset; + + for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) { + // skip measurements that do not match the current photcode + int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode); + if (ecode != thisCode) { continue; } + + if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; + Mcal = getMcal (m, i); + if (isnan(Mcal)) continue; + Mmos = getMmos (m, i); + if (isnan(Mmos)) continue; + Mgrid = getMgrid (m, i); + if (isnan(Mgrid)) continue; + + Mrel = catalog[i].secfilt[Nsecfilt*j+Nsec].M; + xlist[N] = Mrel; + ylist[N] = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt]) - Mcal - Mmos - Mgrid - Mrel; + ilist[N] = PhotInstTiny (&catalog[i].measureT[m]); + N++; + } + } } + + plot_defaults (&graphdata); + graphdata.xmin = PlotMmin; + graphdata.xmax = PlotMmax; + graphdata.ymin = PlotdMmin; + graphdata.ymax = PlotdMmax; + plot_list (&graphdata, xlist, ylist, N, "mag vs dmag", "%s.Mag.png", OUTROOT); + + plot_defaults (&graphdata); + graphdata.ymin = PlotdMmin; + graphdata.ymax = PlotdMmax; + plot_list (&graphdata, ilist, ylist, N, "imag vs dmag", "%s.iMag.png", OUTROOT); } - - plot_defaults (&graphdata); - graphdata.xmin = PlotMmin; - graphdata.xmax = PlotMmax; - graphdata.ymin = PlotdMmin; - graphdata.ymax = PlotdMmax; - plot_list (&graphdata, xlist, ylist, N, "mag vs dmag", "%s.Mag.png", OUTROOT); - - plot_defaults (&graphdata); - graphdata.ymin = PlotdMmin; - graphdata.ymax = PlotdMmax; - plot_list (&graphdata, ilist, ylist, N, "imag vs dmag", "%s.iMag.png", OUTROOT); free (xlist); free (ylist); Index: trunk/Ohana/src/relphot/src/reload_catalogs.c =================================================================== --- trunk/Ohana/src/relphot/src/reload_catalogs.c (revision 31425) +++ trunk/Ohana/src/relphot/src/reload_catalogs.c (revision 31450) @@ -1,3 +1,9 @@ # include "relphot.h" + +# define TIMESTAMP(TIME) \ + gettimeofday (&stop, (void *) NULL); \ + dtime = DTIME (stop, start); \ + TIME += dtime; \ + gettimeofday (&start, (void *) NULL); void reload_catalogs (SkyList *skylist) { @@ -8,8 +14,19 @@ Catalog catalog; + struct timeval start, stop; + double dtime = 0.0; + double time1 = 0.0; + double time2 = 0.0; + double time3 = 0.0; + double time4 = 0.0; + double time5 = 0.0; + double time6 = 0.0; + double time7 = 0.0; + if (VERBOSE) fprintf (stderr, "re-loading catalog data\n"); /* load data from each region file */ for (i = 0; i < skylist[0].Nregions; i++) { + gettimeofday (&start, (void *) NULL); catalog.filename = skylist[0].filename[i]; @@ -20,4 +37,5 @@ continue; } + TIMESTAMP(time1); catalog.catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data @@ -36,20 +54,39 @@ continue; } + TIMESTAMP(time2); + + populate_tiny_values(&catalog); initImageBins (&catalog, 1); initMosaicBins (&catalog, 1); - initGridBins (&catalog, 1); + initGridBins (&catalog, 1); + TIMESTAMP(time3); - findImages (&catalog, 1); - findMosaics (&catalog, 1); + findImages (&catalog, 1); // FX + findMosaics (&catalog, 1); // + TIMESTAMP(time4); setMrelFinal (&catalog); + TIMESTAMP(time5); + dvo_catalog_save (&catalog, VERBOSE); dvo_catalog_unlock (&catalog); + + free_tiny_values(&catalog); dvo_catalog_free (&catalog); + TIMESTAMP(time6); freeImageBins (1); freeMosaicBins (1); freeGridBins (1); + TIMESTAMP(time7); } + + fprintf (stderr, "time1 %f : find catalog\n", time1); + fprintf (stderr, "time2 %f : load catalog\n", time2); + fprintf (stderr, "time3 %f : init imbins\n", time3); + fprintf (stderr, "time4 %f : find images\n", time4); + fprintf (stderr, "time5 %f : set Mrel\n", time5); + fprintf (stderr, "time6 %f : save catalog\n", time6); + fprintf (stderr, "time7 %f : free catalog\n", time7); } Index: trunk/Ohana/src/relphot/src/relphot.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot.c (revision 31425) +++ trunk/Ohana/src/relphot/src/relphot.c (revision 31450) @@ -47,5 +47,5 @@ if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename); } - MARKTIME("load image data: %f sec\n", dtime); + MARKTIME("-- load image data: %f sec\n", dtime); /* load regions and images based on specified sky patch */ @@ -54,5 +54,5 @@ // XXX this is fairly lame: argv[1] is photcode if UserPatchSelect is true skylist = load_images (&db, argv[1], &UserPatch, UserPatchSelect); - MARKTIME("load images: %f sec\n", dtime); + MARKTIME("-- load images: %f sec\n", dtime); /* unlock, if we can (else, unlocked below) */ @@ -61,5 +61,5 @@ /* load catalog data from region files */ catalog = load_catalogs (skylist, &Ncatalog); - MARKTIME("load catalog data: %f sec\n", dtime); + MARKTIME("-- load catalog data: %f sec\n", dtime); /* add in a loop over the catalogs calling dvo_catalog_chipcoords */ @@ -67,5 +67,5 @@ /* match measurements with images, mosaics */ initImageBins (catalog, Ncatalog); - MARKTIME("make image bins: %f sec\n", dtime); + MARKTIME("-- make image bins: %f sec\n", dtime); initMosaicBins (catalog, Ncatalog); @@ -74,12 +74,14 @@ findImages (catalog, Ncatalog); - MARKTIME("set up image indexes: %f sec\n", dtime); + MARKTIME("-- set up image indexes: %f sec\n", dtime); findMosaics (catalog, Ncatalog); /* also sets Grid values */ - MARKTIME("set up mosaic indexes: %f sec\n", dtime); + MARKTIME("-- set up mosaic indexes: %f sec\n", dtime); SAVEPLOT = FALSE; setExclusions (catalog, Ncatalog); + + global_stats (catalog, Ncatalog); if (PLOTSTUFF) { @@ -91,22 +93,4 @@ if (USE_GRID) { int star_toofew; - -# if (USE_DIRECT) - // until we finish the grid analysis, do not reject stars out-of-hand based on ID_STAR_FEW - // XXX this is kind of poor: need to have a better distinctions about STAR_BAD in setMrel vs getMrel - star_toofew = STAR_TOOFEW; - STAR_TOOFEW = 0; - STAR_BAD = ID_STAR_POOR; - - showGridCount (); - setMgridDirect (catalog, Ncatalog); - - STAR_BAD = ID_STAR_POOR | ID_STAR_FEW; - STAR_TOOFEW = star_toofew; - - dump_grid (); - exit (0); - -# else // until we finish the grid analysis, do not reject stars out-of-hand based on ID_STAR_FEW @@ -122,5 +106,4 @@ STAR_BAD = ID_STAR_POOR | ID_STAR_FEW; STAR_TOOFEW = star_toofew; -# endif } @@ -140,19 +123,28 @@ plot_chisq (catalog, Ncatalog); } - if ((i == 1) || (i == 5) || (i == 9) || (i == 13)) clean_measures (catalog, Ncatalog, FALSE); - if ((i == 2) || (i == 6) || (i == 10) || (i == 14)) clean_stars (catalog, Ncatalog); - if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_mosaics (); - if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_images (); + if (i % 6 == 1) rationalize_mosaics (catalog, Ncatalog); + // if (i % 6 == 1) rationalize_images (); + if (i % 6 == 2) clean_measures (catalog, Ncatalog, FALSE); + if (i % 6 == 3) clean_stars (catalog, Ncatalog); + if (i % 6 == 5) clean_mosaics (); + if (i % 6 == 5) clean_images (); + + // if ((i == 1) || (i == 5) || (i == 9) || (i == 13)) clean_measures (catalog, Ncatalog, FALSE); + // if ((i == 2) || (i == 6) || (i == 10) || (i == 14)) clean_stars (catalog, Ncatalog); + // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_mosaics (); + // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_images (); global_stats (catalog, Ncatalog); + MARKTIME("-- finished loop %d: %f sec\n", i, dtime); } - SAVEPLOT = TRUE; - plot_scatter (catalog, Ncatalog); - plot_grid (catalog); - plot_mosaics (); - plot_images (); - plot_stars (catalog, Ncatalog); - plot_chisq (catalog, Ncatalog); - + if (PLOTSTUFF) { + plot_scatter (catalog, Ncatalog); + plot_grid (catalog); + plot_mosaics (); + plot_images (); + plot_stars (catalog, Ncatalog); + plot_chisq (catalog, Ncatalog); + } + if (USE_GRID) dump_grid (); if (!UPDATE) exit (0); @@ -161,7 +153,11 @@ setMcal (catalog, TRUE); setMmos (catalog, TRUE); + MARKTIME("-- finalize Mcal values: %f sec\n", dtime); /* at this point, we have correct cal coeffs in the image/mosaic structures */ - for (i = 0; i < Ncatalog; i++) dvo_catalog_free (&catalog[i]); + for (i = 0; i < Ncatalog; i++) { + free_tiny_values (&catalog[i]); + dvo_catalog_free (&catalog[i]); + } freeImageBins (Ncatalog); freeMosaicBins (Ncatalog); @@ -170,6 +166,7 @@ /* load catalog data from region files, update Mrel include all data */ reload_catalogs (skylist); + MARKTIME("-- updated all catalogs: %f sec\n", dtime); + setMcalFinal (); - reload_images (&db); Index: trunk/Ohana/src/relphot/src/relphot_objects.c =================================================================== --- trunk/Ohana/src/relphot/src/relphot_objects.c (revision 31425) +++ trunk/Ohana/src/relphot/src/relphot_objects.c (revision 31450) @@ -3,5 +3,5 @@ int relphot_objects () { - off_t i, j, k, m; + off_t i, j, k; int Nsecfilt; @@ -38,5 +38,5 @@ } - // XXX consider what gets reset (only PHOTOM flags) + // reset if (RESET) { Nsecfilt = catalog.Nsecfilt; @@ -51,9 +51,5 @@ catalog.secfilt[j*Nsecfilt + k].Ncode = 0; catalog.secfilt[j*Nsecfilt + k].Nused = 0; - } - m = catalog.average[j].measureOffset; - for (k = 0; k < catalog.average[j].Nmeasure; k++) { - catalog.measure[m+k].dbFlags = 0; - catalog.measure[m+k].Mcal = 0; + // XXX reset the photometry flags for secfilt entries? } } Index: trunk/Ohana/src/relphot/src/select_images.c =================================================================== --- trunk/Ohana/src/relphot/src/select_images.c (revision 31425) +++ trunk/Ohana/src/relphot/src/select_images.c (revision 31450) @@ -64,4 +64,6 @@ DmaxSkyRegion = -90.0; + // FILE *ftest = fopen ("relphot.dump.dat", "w"); + /* compare with each region file */ for (i = 0; i < skylist[0].Nregions; i++) { @@ -117,5 +119,10 @@ /* exclude images by photcode */ ecode = GetPhotcodeEquivCodebyCode (timage[i].photcode); - if (ecode != photcode[0].code) continue; + found = FALSE; + int Ns; + for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) { + if (ecode == photcodes[Ns][0].code) found = TRUE; + } + if (!found) continue; /* exclude images by time */ @@ -130,4 +137,7 @@ continue; } + + // XXX temporary test : record center coords for each accepted and rejected chip/mosaic + // double RAo, DECo; /* define image corners - note the DIS images (mosaic phu) are special */ @@ -138,4 +148,5 @@ Xi[3] = -0.5*timage[i].NX; Yi[3] = +0.5*timage[i].NY; Xi[4] = -0.5*timage[i].NX; Yi[4] = -0.5*timage[i].NY; + // XY_to_RD(&RAo, &DECo, 0.0, 0.0, &timage[i].coords); } else { Xi[0] = 0; Yi[0] = 0; @@ -144,10 +155,11 @@ Xi[3] = 0; Yi[3] = timage[i].NY; Xi[4] = 0; Yi[4] = 0; + // XY_to_RD(&RAo, &DECo, 0.5*timage[i].NX, 0.5*timage[i].NY, &timage[i].coords); } found = FALSE; /* transform corners to ra,dec -- costs ~3sec for 3M images (pikake) */ - double RminImage = 360.0; - double RmaxImage = 0.0; + double RminImage = RmidSkyRegion + 180.0; + double RmaxImage = RmidSkyRegion - 180.0; double DminImage = +90.0; double DmaxImage = -90.0; @@ -218,4 +230,7 @@ found_it: + // XXX We claim this is a good image: write to a test file + // fprintf (ftest, "%s : %lf %lf : %f %f %d %x\n", timage[i].name, RAo, DECo, timage[i].Mcal, timage[i].dMcal, timage[i].nstar, timage[i].flags); + image[nimage] = timage[i]; /* always allow 'few' images to succeed, if possible */ @@ -239,4 +254,6 @@ } MARKTIME("finish image selection: %f sec\n", dtime); + + // fclose (ftest); if (VERBOSE) fprintf (stderr, "found "OFF_T_FMT" images\n", nimage); Index: trunk/Ohana/src/relphot/src/setExclusions.c =================================================================== --- trunk/Ohana/src/relphot/src/setExclusions.c (revision 31425) +++ trunk/Ohana/src/relphot/src/setExclusions.c (revision 31450) @@ -1,8 +1,13 @@ # include "relphot.h" + +// this function sets the NOCAL and AREA dbFlags bits for the MeasureTiny elements these +// are used elsewhere (StarOps.c, ImageOps.c, MosaicOps.c, GridOps.c, etc) to skip bad +// measurements. The only exception is 'setMave' which is called by 'relphot_objects', +// and uses the bits read from disk as the test int setExclusions (Catalog *catalog, int Ncatalog) { off_t i, j, k, m, Narea, Nnocal, Ngood; - int ecode; + int ecode, found, Ns; Coords *coords; double r, d, x, y; @@ -11,21 +16,25 @@ for (i = 0; i < Ncatalog; i++) { for (j = 0; j < catalog[i].Naverage; j++) { - m = catalog[i].average[j].measureOffset; - for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) { + m = catalog[i].averageT[j].measureOffset; + for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) { /* select measurements by photcode */ - ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[m].photcode); - if (ecode != photcode[0].code) goto mark_nocal; + ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode); + found = FALSE; + for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) { + if (ecode == photcodes[Ns][0].code) found = TRUE; + } + if (!found) goto mark_nocal; /* select measurements by time */ if (TimeSelect) { - if (catalog[i].measure[m].t < TSTART) goto mark_nocal; - if (catalog[i].measure[m].t > TSTOP) goto mark_nocal; + if (catalog[i].measureT[m].t < TSTART) goto mark_nocal; + if (catalog[i].measureT[m].t > TSTOP) goto mark_nocal; } /* select measurements by mag limit */ if (AreaSelect) { - r = catalog[i].average[j].R + catalog[i].measure[m].dR / 3600.0; - d = catalog[i].average[j].D + catalog[i].measure[m].dD / 3600.0; + r = catalog[i].averageT[j].R + catalog[i].measureT[m].dR / 3600.0; + d = catalog[i].averageT[j].D + catalog[i].measureT[m].dD / 3600.0; if ((coords = getCoords (m, i)) == NULL) goto markbad; RD_to_XY (&x, &y, r, d, coords); @@ -39,10 +48,10 @@ markbad: - catalog[i].measure[m].dbFlags |= ID_MEAS_AREA; + catalog[i].measureT[m].dbFlags |= ID_MEAS_AREA; Narea ++; continue; mark_nocal: - catalog[i].measure[m].dbFlags |= ID_MEAS_NOCAL; + catalog[i].measureT[m].dbFlags |= ID_MEAS_NOCAL; Nnocal ++; continue; Index: trunk/Ohana/src/relphot/src/setMrelFinal.c =================================================================== --- trunk/Ohana/src/relphot/src/setMrelFinal.c (revision 31425) +++ trunk/Ohana/src/relphot/src/setMrelFinal.c (revision 31450) @@ -1,3 +1,7 @@ # include "relphot.h" + +// we've just reloaded the data from disk; we now need to apply the Image/Mosaic/Grid +// calibrations determined by the rest of the program. We also need to set the final +// output dbFlags values void setMrelFinal (Catalog *catalog) { @@ -9,32 +13,42 @@ if (RESET) { - for (i = 0; i < catalog[0].Naverage; i++) { - catalog[0].secfilt[PhotNsec*i+PhotSec].M = NAN; - catalog[0].secfilt[PhotNsec*i+PhotSec].dM = NAN; - catalog[0].secfilt[PhotNsec*i+PhotSec].Xm = NAN_S_SHORT; + int Nsecfilt = GetPhotcodeNsecfilt (); - m = catalog[0].average[i].measureOffset; - for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) { + int Ns; + for (Ns = 0; Ns < Nphotcodes; Ns++) { - /* select measurements by photcode */ - ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[m].photcode); - if (ecode != photcode[0].code) continue; + int thisCode = photcodes[Ns][0].code; + int Nsec = GetPhotcodeNsec(thisCode); + + for (i = 0; i < catalog[0].Naverage; i++) { + catalog[0].secfilt[Nsecfilt*i+Nsec].M = NAN; + catalog[0].secfilt[Nsecfilt*i+Nsec].dM = NAN; + catalog[0].secfilt[Nsecfilt*i+Nsec].Xm = NAN_S_SHORT; + + m = catalog[0].average[i].measureOffset; + for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) { - /* select measurements by time */ - if (TimeSelect) { - if (catalog[0].measure[m].t < TSTART) continue; - if (catalog[0].measure[m].t > TSTOP) continue; + // skip measurements that do not match the current photcode + ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[m].photcode); + if (ecode != thisCode) { continue; } + + /* select measurements by time */ + if (TimeSelect) { + if (catalog[0].measure[m].t < TSTART) continue; + if (catalog[0].measure[m].t > TSTOP) continue; + } + + catalog[0].measure[m].Mcal = 0; + catalog[0].measure[m].dbFlags &= 0xff00; + catalog[0].measure[m].dbFlags &= ~ID_MEAS_POOR_PHOTOM; + catalog[0].measure[m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM; + catalog[0].measure[m].dbFlags &= ~ID_MEAS_AREA; + catalog[0].measure[m].dbFlags &= ~ID_MEAS_NOCAL; } - - catalog[0].measure[m].Mcal = 0; - catalog[0].measure[m].dbFlags &= 0xff00; - catalog[0].measure[m].dbFlags &= ~ID_MEAS_POOR_PHOTOM; - catalog[0].measure[m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM; - catalog[0].measure[m].dbFlags &= ~ID_MEAS_AREA; - catalog[0].measure[m].dbFlags &= ~ID_MEAS_NOCAL; } } } + // this sets flags in the measureT element, not the measure element setExclusions (catalog, 1); /* mark by area */ @@ -53,8 +67,15 @@ setMcalOutput (catalog, 1); - /* clear ID_STAR_POOR, ID_STAR_FEW, ID_MEAS_NOCAL values before writing ??? */ + int Nsecfilt = GetPhotcodeNsecfilt (); + + /* clear ID_STAR_POOR, ID_STAR_FEW values before writing ??? */ + /* ID_MEAS_NOCAL is an internal bit, so it should be cleared */ for (i = 0; i < catalog[0].Naverage; i++) { catalog[0].average[i].flags &= ~ID_STAR_FEW; catalog[0].average[i].flags &= ~ID_STAR_POOR; + for (j = 0; j < Nsecfilt; j++) { + catalog[0].secfilt[i*Nsecfilt+j].flags &= ~ID_STAR_FEW; + catalog[0].secfilt[i*Nsecfilt+j].flags &= ~ID_STAR_POOR; + } m = catalog[0].average[i].measureOffset; for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) { @@ -68,5 +89,5 @@ off_t i, k, m; - int ecode; + int ecode, found, Ns; off_t Ntot, Ntry, Nkeep, Nskip; float mag; @@ -90,10 +111,14 @@ /* clear SKIP for all measures at first */ - catalog[0].measure[m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM; + catalog[0].measureT[m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM; /** never use these measurements (wrong photcode, bad time range) */ /* skipped via NOCAL, don't mark as skipped */ ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[m].photcode); - if (ecode != photcode[0].code) continue; + found = FALSE; + for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) { + if (ecode == photcodes[Ns][0].code) found = TRUE; + } + if (!found) continue; /* skip measurements by time range */ @@ -134,5 +159,6 @@ skip: - catalog[0].measure[m].dbFlags |= ID_MEAS_SKIP_PHOTOM; + catalog[0].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM; + catalog[0].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM; Nskip ++; }