Index: /trunk/Ohana/src/relphot/Makefile =================================================================== --- /trunk/Ohana/src/relphot/Makefile (revision 37036) +++ /trunk/Ohana/src/relphot/Makefile (revision 37037) @@ -34,4 +34,5 @@ $(SRC)/args.$(ARCH).o \ $(SRC)/help.$(ARCH).o \ +$(SRC)/extra.$(ARCH).o \ $(SRC)/bcatalog.$(ARCH).o \ $(SRC)/global_stats.$(ARCH).o \ @@ -82,8 +83,10 @@ $(SRC)/args.$(ARCH).o \ $(SRC)/help.$(ARCH).o \ +$(SRC)/extra.$(ARCH).o \ $(SRC)/synthetic_mags.$(ARCH).o \ $(SRC)/plotstuff.$(ARCH).o \ $(SRC)/liststats.$(ARCH).o \ $(SRC)/initialize.$(ARCH).o \ +$(SRC)/syncfile.$(ARCH).o \ $(SRC)/load_catalogs.$(ARCH).o \ $(SRC)/reload_catalogs.$(ARCH).o \ Index: /trunk/Ohana/src/relphot/doc/mosaic.txt =================================================================== --- /trunk/Ohana/src/relphot/doc/mosaic.txt (revision 37036) +++ /trunk/Ohana/src/relphot/doc/mosaic.txt (revision 37037) @@ -1,2 +1,26 @@ + + +2014.07.09 + +* trying to make sense of mosaic.Mcal vs image.Mcal in the context of UBERCAL: + + * relphot_images: + * load_images (convert raw FITS table to Image structure, select subset matching selection) + * initMosaics (associate mosaics to images (gpc1 photcodes only)) + * initMosaicGrid (define spatial range of mosaics [only used by grid analysis]; set Mcal to and image.Mcal to 0.0) + --- calculate image or mosaic Mcal values + * setMcalFinal (set image.Mcal = mosaicMcal) + + * relphot_parallel_regions: + * assign_images + * makeMosaics (equivalent to initMosaics, but works on full image table, not subset) + * setMosaicCenters (set Mcal to and image.Mcal to 0.0) <--- this is wrong! + (sends images ONLY, not mosaics, to remote machines) + (slurps back new image values, applies to db) + + * relphot_parallel_images: + * makeMosaics (associate mosaics to images, gpc1 photcodes only) + * initMosaicGrid (needed in this one) +--- For 'mosaic' zero points, I define the mosaics by grouping the images Index: /trunk/Ohana/src/relphot/doc/pv2.txt =================================================================== --- /trunk/Ohana/src/relphot/doc/pv2.txt (revision 37037) +++ /trunk/Ohana/src/relphot/doc/pv2.txt (revision 37037) @@ -0,0 +1,53 @@ + +2014.07.03 + +NOTE: 2 questions on these measurements: chisq and ubercal dist? + +I think I have worked out all of the issues in the mag / flux +calibrations for stacks and forced warps. It took me a few passes to +get the right combinations of the 3.56 or the 8.9 offset to the mags. +But now, the stack psf mags, warp psf mags, stack flux, warp flux, and +mean psf mags all agree at the 0.1 mmag level. + +I still have some concerns regarding the flags and astrometry. I am +not set on my laptop to test the primary cell stuff (well, I could +probably reproduce it, but easier to wait until I'm home). I need to +review the 'waterfall' and check that the choices make sense (and +write it down for review). + +2014.07.02 + +I am trying to run a (rather artificial) test of relphot +specifically, I have taken a single smf and an associated warp cmf and +modified the headers so I can load them in DVO and run relphot. Some +issues I've hit: + +* needed to specify PHOTCODE for the stack and fforce fake cmfs +* needed to modify IMAGEID and MJD-OBS so the chips associated + correctly with different gpc1 images +* modified relphot/MosaicOps.c to explicitly skip stack and warp by + photcode when matching chips into mosaics (kind of a hack as it + won't work with other cameras) + +2014.07.01 + +thinking about the issues with relphot for the new PV2 / PV3 schema. +I now have 3 separate sets of mean photometry for each of the stages +cam/staticsky/warp. I need to make sure I am handling everthing +needed. + +Calibration sequence: + +* relphot_images / relphot_parallel_images : calculate the zero points + +Loop around with the standard sequence, but only fitting for the mean +camera mags and camera zero points. At the end of the iteration, fit +the stacks and warps to the new mean object mags. + +Some things I've started to worry about: + +1) QSOs & other variables. Am I detecting and avoiding these well +enough? + +2) what is the appropriate ubercal weight (10, 100, 1000?) + Index: /trunk/Ohana/src/relphot/include/relphot.h =================================================================== --- /trunk/Ohana/src/relphot/include/relphot.h (revision 37036) +++ /trunk/Ohana/src/relphot/include/relphot.h (revision 37037) @@ -4,13 +4,4 @@ # include # include - -# define MARKTIME(MSG,...) { \ - gettimeofday (&stopTimer, (void *) NULL); \ - float dtime = DTIME (stopTimer, startTimer); \ - fprintf (stderr, MSG, __VA_ARGS__); } - -# define INITTIME \ - struct timeval startTimer, stopTimer; \ - gettimeofday (&startTimer, (void *) NULL); /* # define GRID_V1 */ @@ -119,5 +110,5 @@ float M; float dM; - float Xm; + float Mchisq; int Nsec; unsigned int objID; @@ -130,5 +121,5 @@ float dMcal; float dMagSys; - float Xm; + short Xm; int nFitPhotom; int flags; @@ -317,5 +308,5 @@ int findMosaics PROTO((Catalog *catalog, int Ncatalog, int doMosaicList)); -void makeMosaics (Image *image, off_t Nimage); +void makeMosaics (Image *image, off_t Nimage, int mergeMcal); Mosaic *getMosaicForImage (off_t im); void setMosaicCenters (Image *image, off_t Nimage); @@ -360,6 +351,6 @@ void liststats_setmode PROTO((StatType *stats, char *strmode)); int liststats PROTO((double *value, double *dvalue, double *wvalue, int N, StatType *stats)); -Catalog *load_catalogs PROTO((SkyList *skylist, int *Ncatalog, int hostID, char *hostpath)); -Catalog *load_catalogs_parallel PROTO((SkyList *sky, int *Ncatalog)); +Catalog *load_catalogs PROTO((SkyList *skylist, int *Ncatalog, int hostID, char *hostpath, char *syncfile)); +Catalog *load_catalogs_parallel PROTO((SkyList *sky, int *Ncatalog, char *syncfile)); SkyList *load_images PROTO((FITS_DB *db, char *regionName, SkyRegion *region)); @@ -509,2 +500,9 @@ int relphot_parallel_images (); int relphot_parallel_regions (); + +// in extra.c +int isGPC1chip (int photcode); +int isGPC1stack (int photcode); +int isGPC1warp (int photcode); + + Index: /trunk/Ohana/src/relphot/src/BrightCatalog.c =================================================================== --- /trunk/Ohana/src/relphot/src/BrightCatalog.c (revision 37036) +++ /trunk/Ohana/src/relphot/src/BrightCatalog.c (revision 37037) @@ -63,6 +63,6 @@ // need to create and assign to flat-field correction - GET_COLUMN(dR, "RA_OFF", float); - GET_COLUMN(dD, "DEC_OFF", float); + GET_COLUMN(R, "RA", double); + GET_COLUMN(D, "DEC", double); GET_COLUMN(M, "MAG_SYS", float); GET_COLUMN(Mcal, "MAG_CAL", float); @@ -85,6 +85,6 @@ ALLOCATE (measure, MeasureTiny, Nrow); for (i = 0; i < Nrow; i++) { - measure[i].dR = dR[i]; - measure[i].dD = dD[i]; + measure[i].R = R[i]; + measure[i].D = D[i]; measure[i].M = M[i]; measure[i].Mcal = Mcal[i]; @@ -105,6 +105,6 @@ fprintf (stderr, "loaded data for %lld measure\n", (long long) Nrow); - free (dR ); - free (dD ); + free (R ); + free (D ); free (M ); free (Mcal ); @@ -136,6 +136,6 @@ // need to create and assign to flat-field correction - GET_COLUMN(R, "RA", double); - GET_COLUMN(D, "DEC", double); + GET_COLUMN(R, "RA", double); + GET_COLUMN(D, "DEC", double); GET_COLUMN(Nmeasure, "NMEAS", int); GET_COLUMN(measureOffset, "MEASURE_OFF", int); @@ -182,12 +182,12 @@ // need to create and assign to flat-field correction - GET_COLUMN(M, "MAG", float); - GET_COLUMN(dM, "MAG_ERR", float); - GET_COLUMN(Xm, "MAG_CHI", float); - GET_COLUMN(flags, "FLAGS", int); - GET_COLUMN(Ncode, "NCODE", short); - GET_COLUMN(Nused, "NUSED", short); - GET_COLUMN(M_20, "MAG_20", short); - GET_COLUMN(M_80, "MAG_80", short); + GET_COLUMN(M, "MAG", float); + GET_COLUMN(dM, "MAG_ERR", float); + GET_COLUMN(Mchisq, "MAG_CHI", float); + GET_COLUMN(flags, "FLAGS", int); + GET_COLUMN(Ncode, "NCODE", short); + GET_COLUMN(Nused, "NUSED", short); + GET_COLUMN(Mmin, "MAG_MIN", float); + GET_COLUMN(Mmax, "MAG_MAX", float); gfits_free_header (&theader); gfits_free_table (&ftable); @@ -196,23 +196,23 @@ ALLOCATE (secfilt, SecFilt, Nrow); for (i = 0; i < Nrow; i++) { - secfilt[i].M = M[i]; - secfilt[i].dM = dM[i]; - secfilt[i].Xm = Xm[i]; - secfilt[i].flags = flags[i]; - secfilt[i].Ncode = Ncode[i]; - secfilt[i].Nused = Nused[i]; - secfilt[i].M_20 = M_20[i]; - secfilt[i].M_80 = M_80[i]; + secfilt[i].M = M[i]; + secfilt[i].dM = dM[i]; + secfilt[i].Mchisq = Mchisq[i]; + secfilt[i].flags = flags[i]; + secfilt[i].Ncode = Ncode[i]; + secfilt[i].Nused = Nused[i]; + secfilt[i].Mmin = Mmin[i]; + secfilt[i].Mmax = Mmax[i]; } fprintf (stderr, "loaded data for %lld secfilt\n", (long long) Nrow); - free (M ); - free (dM ); - free (Xm ); - free (flags); - free (Ncode); - free (Nused); - free (M_20 ); - free (M_80 ); + free (M ); + free (dM ); + free (Mchisq); + free (flags ); + free (Ncode ); + free (Nused ); + free (Mmin ); + free (Mmax ); catalog->secfilt = secfilt; // assert Nsecfilt * Naverage = Nrow? @@ -265,6 +265,6 @@ gfits_create_table_header (&theader, "BINTABLE", "MEASURE_TINY"); - gfits_define_bintable_column (&theader, "E", "RA_OFF", "ra offset", "arcsec", 1.0, 0.0); - gfits_define_bintable_column (&theader, "E", "DEC_OFF", "dec offset", "arcsec", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "RA", "ra", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "DEC", "dec", "degree", 1.0, 0.0); gfits_define_bintable_column (&theader, "E", "MAG_SYS", "magnitude (sys)", NULL, 1.0, 0.0); gfits_define_bintable_column (&theader, "E", "MAG_CAL", "magnitude (cal)", NULL, 1.0, 0.0); @@ -274,11 +274,11 @@ gfits_define_bintable_column (&theader, "E", "Y_CCD", "ccd y coord", "pix", 1.0, 0.0); gfits_define_bintable_column (&theader, "E", "EXPTIME", "-2.5 * log (exposure time)", "sec", 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "TIME", "time of exp", "sec", 1.0, 1.0*0x8000); - gfits_define_bintable_column (&theader, "J", "AVE_REF", "pointer to average table", NULL, 1.0, 1.0*0x8000); - gfits_define_bintable_column (&theader, "J", "IMAGE_ID", "image", NULL, 1.0, 1.0*0x8000); - gfits_define_bintable_column (&theader, "J", "DB_FLAGS", "flags", NULL, 1.0, 1.0*0x8000); - gfits_define_bintable_column (&theader, "J", "PHOT_FLAGS", "photflags", NULL, 1.0, 1.0*0x8000); - gfits_define_bintable_column (&theader, "J", "CAT_ID", "catalog", NULL, 1.0, 1.0*0x8000); - gfits_define_bintable_column (&theader, "I", "PHOTCODE", "photcode", NULL, 1.0, 1.0*0x80); + gfits_define_bintable_column (&theader, "J", "TIME", "time of exp", "sec", 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "J", "AVE_REF", "pointer to average table", NULL, 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "J", "IMAGE_ID", "image", NULL, 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "J", "DB_FLAGS", "flags", NULL, 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "J", "PHOT_FLAGS", "photflags", NULL, 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "J", "CAT_ID", "catalog", NULL, 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "I", "PHOTCODE", "photcode", NULL, 1.0, FT_BZERO_INT16); // generate the output array that carries the data @@ -286,26 +286,26 @@ // create intermediate storage arrays - float *dR ; ALLOCATE (dR , float, catalog->Nmeasure); - float *dD ; ALLOCATE (dD , float, catalog->Nmeasure); - float *M ; ALLOCATE (M , float, catalog->Nmeasure); - float *Mcal ; ALLOCATE (Mcal , float, catalog->Nmeasure); - float *dM ; ALLOCATE (dM , float, catalog->Nmeasure); - float *airmass ; ALLOCATE (airmass , float, catalog->Nmeasure); - float *Xccd ; ALLOCATE (Xccd , float, catalog->Nmeasure); - float *Yccd ; ALLOCATE (Yccd , float, catalog->Nmeasure); - float *dt ; ALLOCATE (dt , float, catalog->Nmeasure); - int *t ; ALLOCATE (t , int , catalog->Nmeasure); - int *averef ; ALLOCATE (averef , int , catalog->Nmeasure); - int *imageID ; ALLOCATE (imageID , int , catalog->Nmeasure); - int *dbFlags ; ALLOCATE (dbFlags , int , catalog->Nmeasure); - int *photFlags ; ALLOCATE (photFlags, int , catalog->Nmeasure); - int *catID ; ALLOCATE (catID , int , catalog->Nmeasure); - short *photcode ; ALLOCATE (photcode , short, catalog->Nmeasure); + double *R ; ALLOCATE (R , double, catalog->Nmeasure); + double *D ; ALLOCATE (D , double, catalog->Nmeasure); + float *M ; ALLOCATE (M , float, catalog->Nmeasure); + float *Mcal ; ALLOCATE (Mcal , float, catalog->Nmeasure); + float *dM ; ALLOCATE (dM , float, catalog->Nmeasure); + float *airmass ; ALLOCATE (airmass , float, catalog->Nmeasure); + float *Xccd ; ALLOCATE (Xccd , float, catalog->Nmeasure); + float *Yccd ; ALLOCATE (Yccd , float, catalog->Nmeasure); + float *dt ; ALLOCATE (dt , float, catalog->Nmeasure); + int *t ; ALLOCATE (t , int , catalog->Nmeasure); + int *averef ; ALLOCATE (averef , int , catalog->Nmeasure); + int *imageID ; ALLOCATE (imageID , int , catalog->Nmeasure); + int *dbFlags ; ALLOCATE (dbFlags , int , catalog->Nmeasure); + int *photFlags ; ALLOCATE (photFlags, int , catalog->Nmeasure); + int *catID ; ALLOCATE (catID , int , catalog->Nmeasure); + short *photcode ; ALLOCATE (photcode , short, catalog->Nmeasure); // assign the storage arrays MeasureTiny *measure = catalog->measure; for (i = 0; i < catalog->Nmeasure; i++) { - dR[i] = measure[i].dR ; - dD[i] = measure[i].dD ; + R[i] = measure[i].R ; + D[i] = measure[i].D ; M[i] = measure[i].M ; Mcal[i] = measure[i].Mcal ; @@ -325,6 +325,6 @@ // add the columns to the output array - gfits_set_bintable_column (&theader, &ftable, "RA_OFF", dR, catalog->Nmeasure); - gfits_set_bintable_column (&theader, &ftable, "DEC_OFF", dD, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "RA", R, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "DEC", D, catalog->Nmeasure); gfits_set_bintable_column (&theader, &ftable, "MAG_SYS", M, catalog->Nmeasure); gfits_set_bintable_column (&theader, &ftable, "MAG_CAL", Mcal, catalog->Nmeasure); @@ -342,6 +342,6 @@ gfits_set_bintable_column (&theader, &ftable, "PHOTCODE", photcode, catalog->Nmeasure); - free (dR ); - free (dD ); + free (R ); + free (D ); free (M ); free (Mcal ); @@ -428,12 +428,12 @@ gfits_create_table_header (&theader, "BINTABLE", "SECFILT"); - gfits_define_bintable_column (&theader, "E", "MAG", "ra offset", "arcsec", 1.0, 0.0); - gfits_define_bintable_column (&theader, "E", "MAG_ERR", "dec offset", "arcsec", 1.0, 0.0); - gfits_define_bintable_column (&theader, "E", "MAG_CHI", "magnitude (sys)", NULL, 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "FLAGS", "magnitude (cal)", NULL, 1.0, 0.0); - gfits_define_bintable_column (&theader, "I", "NCODE", "magnitude (err)", NULL, 1.0, 0.0); - gfits_define_bintable_column (&theader, "I", "NUSED", "airmass", NULL, 1.0, 0.0); - gfits_define_bintable_column (&theader, "I", "MAG_20", "ccd x coord", "pix", 1.0, 0.0); - gfits_define_bintable_column (&theader, "I", "MAG_80", "ccd y coord", "pix", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG", "", "arcsec", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_ERR", "", "arcsec", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_CHI", "", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "FLAGS", "", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "I", "NCODE", "", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "I", "NUSED", "", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_MIN", "min valid mag", "mag", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_MAX", "max valid mag", "mag", 1.0, 0.0); // generate the output array that carries the data @@ -447,42 +447,42 @@ float *M ; ALLOCATE (M , float, Nsec); float *dM ; ALLOCATE (dM , float, Nsec); - float *Xm ; ALLOCATE (Xm , float, Nsec); + float *Mchisq ; ALLOCATE (Mchisq , float, Nsec); int *flags ; ALLOCATE (flags , int, Nsec); short *Ncode ; ALLOCATE (Ncode , short, Nsec); short *Nused ; ALLOCATE (Nused , short, Nsec); - short *M_20 ; ALLOCATE (M_20 , short, Nsec); - short *M_80 ; ALLOCATE (M_80 , short, Nsec); + float *Mmin ; ALLOCATE (Mmin , float, Nsec); + float *Mmax ; ALLOCATE (Mmax , float, Nsec); // assign the storage arrays SecFilt *secfilt = catalog->secfilt; for (i = 0; i < Nsec; i++) { - M [i] = secfilt[i]. M ; - dM [i] = secfilt[i]. dM ; - Xm [i] = secfilt[i]. Xm ; - flags [i] = secfilt[i]. flags ; - Ncode [i] = secfilt[i]. Ncode ; - Nused [i] = secfilt[i]. Nused ; - M_20 [i] = secfilt[i]. M_20 ; - M_80 [i] = secfilt[i]. M_80 ; + M [i] = secfilt[i].M ; + dM [i] = secfilt[i].dM ; + Mchisq[i] = secfilt[i].Mchisq ; + flags [i] = secfilt[i].flags ; + Ncode [i] = secfilt[i].Ncode ; + Nused [i] = secfilt[i].Nused ; + Mmin [i] = secfilt[i].Mmin ; + Mmax [i] = secfilt[i].Mmax ; } // add the columns to the output array - gfits_set_bintable_column (&theader, &ftable, "MAG", M , Nsec); - gfits_set_bintable_column (&theader, &ftable, "MAG_ERR", dM , Nsec); - gfits_set_bintable_column (&theader, &ftable, "MAG_CHI", Xm , Nsec); - gfits_set_bintable_column (&theader, &ftable, "FLAGS", flags, Nsec); - gfits_set_bintable_column (&theader, &ftable, "NCODE", Ncode, Nsec); - gfits_set_bintable_column (&theader, &ftable, "NUSED", Nused, Nsec); - gfits_set_bintable_column (&theader, &ftable, "MAG_20", M_20 , Nsec); - gfits_set_bintable_column (&theader, &ftable, "MAG_80", M_80 , Nsec); + gfits_set_bintable_column (&theader, &ftable, "MAG", M , Nsec); + gfits_set_bintable_column (&theader, &ftable, "MAG_ERR", dM , Nsec); + gfits_set_bintable_column (&theader, &ftable, "MAG_CHI", Mchisq, Nsec); + gfits_set_bintable_column (&theader, &ftable, "FLAGS", flags , Nsec); + gfits_set_bintable_column (&theader, &ftable, "NCODE", Ncode , Nsec); + gfits_set_bintable_column (&theader, &ftable, "NUSED", Nused , Nsec); + gfits_set_bintable_column (&theader, &ftable, "MAG_MIN", Mmin , Nsec); + gfits_set_bintable_column (&theader, &ftable, "MAG_MAX", Mmax , Nsec); free (M ); free (dM ); - free (Xm ); + free (Mchisq ); free (flags ); free (Ncode ); free (Nused ); - free (M_20 ); - free (M_80 ); + free (Mmin ); + free (Mmax ); gfits_fwrite_Theader (f, &theader); Index: /trunk/Ohana/src/relphot/src/GridOps.c =================================================================== --- /trunk/Ohana/src/relphot/src/GridOps.c (revision 37036) +++ /trunk/Ohana/src/relphot/src/GridOps.c (revision 37037) @@ -414,5 +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 = PhotCatTiny (&catalog[c].measureT[m]); + Msys = PhotCatTiny (&catalog[c].measureT[m], MAG_CLASS_PSF); if (isnan(Msys)) { Nsys++; @@ -573,5 +573,5 @@ n = catalog[c].measureT[m].averef; - Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]); + Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF); if (isnan(Msys)) { Nsys++; @@ -651,5 +651,5 @@ n = catalog[c].measureT[m].averef; - Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]); + Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF); xlist[N] = Xmeas[c][m]; Index: /trunk/Ohana/src/relphot/src/ImageMagIO.c =================================================================== --- /trunk/Ohana/src/relphot/src/ImageMagIO.c (revision 37036) +++ /trunk/Ohana/src/relphot/src/ImageMagIO.c (revision 37037) @@ -61,9 +61,9 @@ GET_COLUMN (dMcal, "MCAL_ERR", float); GET_COLUMN (dMagSys, "MCAL_SYSERR", float); - GET_COLUMN (Xm, "CHISQ", float); GET_COLUMN (nFitPhotom, "NFIT", int); GET_COLUMN (flags, "FLAGS", int); GET_COLUMN (ubercalDist, "UDIST", int); GET_COLUMN (imageID, "ID", int); + GET_COLUMN (Xm, "CHISQ", short); // free the memory associated with the FITS files @@ -121,12 +121,12 @@ gfits_create_table_header (&theader, "BINTABLE", "IMAGE_MAGS"); - gfits_define_bintable_column (&theader, "E", "MCAL", "cal offset", "magnitudes", 1.0, 0.0); - gfits_define_bintable_column (&theader, "E", "MCAL_ERR", "cal error", "magnitudes", 1.0, 0.0); - gfits_define_bintable_column (&theader, "E", "MCAL_SYSERR", "systematic error", "magnitudes", 1.0, 0.0); - gfits_define_bintable_column (&theader, "E", "CHISQ", "cal chisq", "unitless", 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "NFIT", "number of fitted stars", "unitless", 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "FLAGS", "analysis flags", "unitless", 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "UDIST", "distance to ubercal images", "images", 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "ID", "image ID", "unitless", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MCAL", "cal offset", "magnitudes", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MCAL_ERR", "cal error", "magnitudes", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MCAL_SYSERR", "systematic error", "magnitudes", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "NFIT", "number of fitted stars", "unitless", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "FLAGS", "analysis flags", "unitless", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "UDIST", "distance to ubercal images", "images", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "ID", "image ID", "unitless", 1.0, 0.0); + gfits_define_bintable_column (&theader, "I", "CHISQ", "cal chisq", "unitless", 1.0, FT_BZERO_INT16); // generate the output array that carries the data Index: /trunk/Ohana/src/relphot/src/ImageOps.c =================================================================== --- /trunk/Ohana/src/relphot/src/ImageOps.c (revision 37036) +++ /trunk/Ohana/src/relphot/src/ImageOps.c (revision 37037) @@ -531,5 +531,5 @@ liststats_setmode (&stats, STATMODE); - if (FREEZE_IMAGES) return; + // FREEZE_IMAGES only applies to mosaic data (eg, gpc1) fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE); @@ -561,4 +561,6 @@ if (!bad) continue; } + + if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue; // UBERCAL image: if this is an ubercal image, set minUbercalDist to 0: @@ -608,5 +610,5 @@ n = catalog[c].measureT[m].averef; - Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]); + Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF); if (isnan(Msys)) { Nsys++; @@ -695,5 +697,5 @@ double MaxOffset, MaxScatter, MedOffset; - if (FREEZE_IMAGES) return; + // if (FREEZE_IMAGES) return; if (VERBOSE) fprintf (stderr, "marking poor images\n"); @@ -706,4 +708,7 @@ for (i = N = 0; i < Nimage; i++) { if (image[i].flags & IMAGE_BAD) continue; + + if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue; + mlist[N] = image[i].Mcal; slist[N] = image[i].dMcal; @@ -797,5 +802,5 @@ Graphdata graphdata; - if (FREEZE_IMAGES) return; + // if (FREEZE_IMAGES) return; ALLOCATE (xlist, double, Nimage); @@ -810,8 +815,14 @@ float mindMcal = +100.0; float maxdMcal = -100.0; - for (i = 0; i < Nimage; i++) { - Mlist[i] = image[i].Mcal; - dlist[i] = image[i].dMcal; - xlist[i] = image[i].secz; + + int Nplot = 0; + + for (i = 0; i < Nimage; i++) { + + if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue; + + Mlist[Nplot] = image[i].Mcal; + dlist[Nplot] = image[i].dMcal; + xlist[Nplot] = image[i].secz; minAirmass = MIN (image[i].secz, minAirmass); maxAirmass = MAX (image[i].secz, maxAirmass); @@ -820,4 +831,6 @@ mindMcal = MIN (image[i].dMcal, mindMcal); maxdMcal = MAX (image[i].dMcal, maxdMcal); + + Nplot ++; } @@ -852,5 +865,8 @@ for (i = 0; i < NBIN; i++) xlist[i] = 0.00025*i; bzero (Mlist, NBIN*sizeof(double)); - for (i = 0; i < Nimage; i++) { + + for (i = 0; i < Nimage; i++) { + if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue; + bin = image[i].dMcal / 0.00025; bin = MAX (0, MIN (NBIN - 1, bin)); @@ -877,5 +893,7 @@ StatType stats; bzero (&stats, sizeof (StatType)); - if (FREEZE_IMAGES) return (stats); + + // we no longer blindly apply FREEZE_IMAGES to all images, only to mosaics + // if (FREEZE_IMAGES) return (stats); ALLOCATE (list, double, Nimage); @@ -885,4 +903,6 @@ for (i = 0; i < Nimage; i++) { if (image[i].flags & IMAGE_BAD) continue; + + if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue; N = 0; @@ -920,5 +940,5 @@ bzero (&stats, sizeof (StatType)); - if (FREEZE_IMAGES) return (stats); + // if (FREEZE_IMAGES) return (stats); ALLOCATE (list, double, Nimage); @@ -929,4 +949,6 @@ if (image[i].flags & IMAGE_BAD) continue; + + if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue; list[n] = pow (10.0, 0.01*image[i].Xm); @@ -950,5 +972,5 @@ bzero (&stats, sizeof (StatType)); - if (FREEZE_IMAGES) return (stats); + // if (FREEZE_IMAGES) return (stats); ALLOCATE (list, double, Nimage); @@ -959,4 +981,6 @@ if (image[i].flags & IMAGE_BAD) continue; + + if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue; list[n] = image[i].Mcal; @@ -980,5 +1004,5 @@ bzero (&stats, sizeof (StatType)); - if (FREEZE_IMAGES) return (stats); + // if (FREEZE_IMAGES) return (stats); ALLOCATE (list, double, Nimage); @@ -989,4 +1013,6 @@ if (image[i].flags & IMAGE_BAD) continue; + + if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue; list[n] = image[i].dMcal; Index: /trunk/Ohana/src/relphot/src/ImageSubset.c =================================================================== --- /trunk/Ohana/src/relphot/src/ImageSubset.c (revision 37036) +++ /trunk/Ohana/src/relphot/src/ImageSubset.c (revision 37037) @@ -126,19 +126,15 @@ gfits_create_table_header (&theader, "BINTABLE", "IMAGE_SUBSET"); - gfits_define_bintable_column (&theader, "E", "MCAL", "zero point offset", "magnitudes", 1.0, 0.0); - gfits_define_bintable_column (&theader, "E", "MCAL_ERR", "zero point error", "magnitudes", 1.0, 0.0); - - gfits_define_bintable_column (&theader, "J", "IMAGE_ID", "image ID", NULL, 1.0, 1.0*0x8000); - gfits_define_bintable_column (&theader, "J", "PHOTOM_MAP", "map", NULL, 1.0, 1.0*0x8000); - gfits_define_bintable_column (&theader, "J", "FLAGS", "flags", NULL, 1.0, 1.0*0x8000); - - gfits_define_bintable_column (&theader, "J", "TESS_ID", "ID", NULL, 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "PROJ_ID", "ID", NULL, 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "SKYCELL_ID", "ID", NULL, 1.0, 0.0); - - gfits_define_bintable_column (&theader, "J", "TZERO", "exposure start", NULL, 1.0, 1.0*0x8000); - gfits_define_bintable_column (&theader, "I", "TRATE", "tti rate", NULL, 1.0, 0.0); - - gfits_define_bintable_column (&theader, "I", "UBERCAL_DIST", "ubercal distance", NULL, 1.0, 1.0*0x80); + gfits_define_bintable_column (&theader, "E", "MCAL", "zero point offset", "magnitudes", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MCAL_ERR", "zero point error", "magnitudes", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "IMAGE_ID", "image ID", NULL, 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "J", "PHOTOM_MAP", "map", NULL, 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "J", "FLAGS", "flags", NULL, 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "J", "TESS_ID", "ID", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "PROJ_ID", "ID", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "SKYCELL_ID", "ID", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "TZERO", "exposure start", NULL, 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "I", "TRATE", "tti rate", NULL, 1.0, FT_BZERO_INT16); + gfits_define_bintable_column (&theader, "I", "UBERCAL_DIST", "ubercal distance", NULL, 1.0, 0.0); // generate the output array that carries the data @@ -152,6 +148,6 @@ // create intermediate storage arrays - ALLOCATE (Mcal, float, Nimage); - ALLOCATE (dMcal, float, Nimage); + ALLOCATE (Mcal, float, Nimage); + ALLOCATE (dMcal, float, Nimage); ALLOCATE (imageID, unsigned int, Nimage); ALLOCATE (map, unsigned int, Nimage); Index: /trunk/Ohana/src/relphot/src/MeanMagIO.c =================================================================== --- /trunk/Ohana/src/relphot/src/MeanMagIO.c (revision 37036) +++ /trunk/Ohana/src/relphot/src/MeanMagIO.c (revision 37037) @@ -60,5 +60,5 @@ GET_COLUMN (M, "MAG", float); GET_COLUMN (dM, "MAG_ERR", float); - GET_COLUMN (Xm, "MAG_CHISQ", float); + GET_COLUMN (Mchisq, "MAG_CHISQ", float); GET_COLUMN (Nsec, "NSEC", int); GET_COLUMN (objID, "OBJ_ID", int); @@ -73,19 +73,19 @@ ALLOCATE (meanmags, MeanMag, Nrow); for (i = 0; i < Nrow; i++) { - meanmags[i].M = M [i]; - meanmags[i].dM = dM [i]; - meanmags[i].Xm = Xm [i]; - meanmags[i].Nsec = Nsec [i]; - meanmags[i].objID = objID[i]; - meanmags[i].catID = catID[i]; + meanmags[i].M = M [i]; + meanmags[i].dM = dM [i]; + meanmags[i].Mchisq = Mchisq [i]; + meanmags[i].Nsec = Nsec [i]; + meanmags[i].objID = objID [i]; + meanmags[i].catID = catID [i]; } fprintf (stderr, "loaded data for %lld objects (* filters)\n", (long long) Nrow); - free (M ); - free (dM ); - free (Xm ); - free (Nsec ); - free (objID); - free (catID); + free (M ); + free (dM ); + free (Mchisq); + free (Nsec ); + free (objID ); + free (catID ); *nmeanmags = Nrow; @@ -118,12 +118,12 @@ gfits_define_bintable_column (&theader, "E", "MAG_ERR", "mean magnitude error", "magnitudes", 1.0, 0.0); gfits_define_bintable_column (&theader, "E", "MAG_CHISQ", "mean magnitude chisq", "unitless", 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "NSEC", "secfilt sequence", NULL, 1.0, 1.0*0x8000); - gfits_define_bintable_column (&theader, "J", "OBJ_ID", "object ID", NULL, 1.0, 1.0*0x8000); - gfits_define_bintable_column (&theader, "J", "CAT_ID", "catalog ID", NULL, 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "NSEC", "secfilt sequence", NULL, 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "J", "OBJ_ID", "object ID", NULL, 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "J", "CAT_ID", "catalog ID", NULL, 1.0, FT_BZERO_INT32); // generate the output array that carries the data gfits_create_table (&theader, &ftable); - float *M, *dM, *Xm; + float *M, *dM, *Mchisq; int *Nsec; unsigned int *objID, *catID; @@ -132,5 +132,5 @@ ALLOCATE (M, float, Nmeanmags); ALLOCATE (dM, float, Nmeanmags); - ALLOCATE (Xm, float, Nmeanmags); + ALLOCATE (Mchisq, float, Nmeanmags); ALLOCATE (Nsec, int, Nmeanmags); ALLOCATE (objID, unsigned int, Nmeanmags); @@ -139,10 +139,10 @@ // assign the storage arrays for (i = 0; i < Nmeanmags; i++) { - M [i] = meanmags[i].M ; - dM [i] = meanmags[i].dM ; - Xm [i] = meanmags[i].Xm ; - Nsec [i] = meanmags[i].Nsec ; - objID[i] = meanmags[i].objID; - catID[i] = meanmags[i].catID; + M [i] = meanmags[i].M ; + dM [i] = meanmags[i].dM ; + Mchisq[i] = meanmags[i].Mchisq; + Nsec [i] = meanmags[i].Nsec ; + objID [i] = meanmags[i].objID ; + catID [i] = meanmags[i].catID ; } @@ -150,15 +150,15 @@ gfits_set_bintable_column (&theader, &ftable, "MAG", M, Nmeanmags); gfits_set_bintable_column (&theader, &ftable, "MAG_ERR", dM, Nmeanmags); - gfits_set_bintable_column (&theader, &ftable, "MAG_CHISQ", Xm, Nmeanmags); + gfits_set_bintable_column (&theader, &ftable, "MAG_CHISQ", Mchisq, Nmeanmags); gfits_set_bintable_column (&theader, &ftable, "NSEC", Nsec, Nmeanmags); gfits_set_bintable_column (&theader, &ftable, "OBJ_ID", objID, Nmeanmags); gfits_set_bintable_column (&theader, &ftable, "CAT_ID", catID, Nmeanmags); - free (M ); - free (dM ); - free (Xm ); - free (Nsec ); - free (objID); - free (catID); + free (M ); + free (dM ); + free (Mchisq); + free (Nsec ); + free (objID ); + free (catID ); FILE *f = fopen (filename, "w"); Index: /trunk/Ohana/src/relphot/src/MosaicOps.c =================================================================== --- /trunk/Ohana/src/relphot/src/MosaicOps.c (revision 37036) +++ /trunk/Ohana/src/relphot/src/MosaicOps.c (revision 37037) @@ -1,4 +1,5 @@ # include "relphot.h" void plot_setMcal (double *list, int Npts, StatType *stats, float clouds); +off_t findMosaic (unsigned int *startTimes, off_t Nmosaic, unsigned int start); // see discussion in ImagesOps.c re: IDX_T @@ -36,4 +37,18 @@ // MosaicN_image was 'Nimlist' // MosaicToImage was 'imlist' + +void sort_times (unsigned int *T, int N) { + +# define SWAPFUNC(A,B){ unsigned int tmp; \ + tmp = T[A]; T[A] = T[B]; T[B] = tmp; \ +} +# define COMPARE(A,B)(T[A] < T[B]) + + OHANA_SORT (N, COMPARE, SWAPFUNC); + +# undef SWAPFUNC +# undef COMPARE + +} /* find mosaic frames (unique time periods & photcode name matches mosaic) */ @@ -96,12 +111,12 @@ /* a new mosaic, define ranges -- preserve the original values incase this image is not used */ - mosaic[Nmosaic].start = start; - mosaic[Nmosaic].stop = stop; - mosaic[Nmosaic].Mcal = 0.0; // note : mosaic stores only offsets relative to the original image values - mosaic[Nmosaic].dMcal = 0.0; // note : at the end, mosaic.Mcal is added back to the input images - mosaic[Nmosaic].dMsys = 0.0; - mosaic[Nmosaic].Xm = 0.0; + mosaic[Nmosaic].start = start; + mosaic[Nmosaic].stop = stop; + mosaic[Nmosaic].Mcal = 0.0; // note : mosaic stores only offsets relative to the original image values + mosaic[Nmosaic].dMcal = 0.0; // note : at the end, mosaic.Mcal is added back to the input images + mosaic[Nmosaic].dMsys = 0.0; + mosaic[Nmosaic].Xm = 0.0; mosaic[Nmosaic].flags = image[i].flags; - mosaic[Nmosaic].secz = image[i].secz; + mosaic[Nmosaic].secz = image[i].secz; mosaic[Nmosaic].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode); @@ -136,20 +151,4 @@ } -off_t findMosaic (unsigned int *startTimes, off_t Nmosaic, unsigned int start); - -void sort_times (unsigned int *T, int N) { - -# define SWAPFUNC(A,B){ unsigned int tmp; \ - tmp = T[A]; T[A] = T[B]; T[B] = tmp; \ -} -# define COMPARE(A,B)(T[A] < T[B]) - - OHANA_SORT (N, COMPARE, SWAPFUNC); - -# undef SWAPFUNC -# undef COMPARE - -} - /* find mosaic frames (unique time periods) (NOTE : we do NOT require matching photcodes...) this function will also identify the images NOT in the subset which belong to a selected mosaic @@ -157,7 +156,6 @@ void initMosaics (Image *subset, off_t Nsubset, Image *image, char *inSubset, off_t Nimage) { - off_t i, j, status, found, NMOSAIC, *MosaicN_IMAGE; + off_t i, j, found, NMOSAIC, *MosaicN_IMAGE; unsigned int start, stop, *startTimes, *startTimesMosaic; - char *pname; if (!MOSAIC_ZEROPT) return; @@ -172,8 +170,11 @@ // generate a list of all subset image start times ALLOCATE (startTimes, unsigned int, Nsubset); + int Nmoschip = 0; for (i = 0; i < Nsubset; i++) { - startTimes[i] = subset[i].tzero; - } - sort_times (startTimes, Nsubset); + if (!isGPC1chip(subset[i].photcode)) continue; + startTimes[Nmoschip] = subset[i].tzero; + Nmoschip ++; + } + sort_times (startTimes, Nmoschip); MARKTIME("create array of all image obstimes: %f sec\n", dtime); @@ -184,9 +185,6 @@ // generate a list of the unique start times (these define the mosaics) - for (i = 0; i < Nsubset; i++) { - if (startTimes[i] < startTimesMosaic[Nmosaic]) { - fprintf (stderr, "error?\n"); - abort(); - } + for (i = 0; i < Nmoschip; i++) { + myAssert (startTimes[i] >= startTimesMosaic[Nmosaic], "times out of order?"); if (startTimes[i] == startTimesMosaic[Nmosaic]) continue; Nmosaic ++; @@ -210,14 +208,14 @@ for (i = 0; i < Nmosaic; i++) { /* a new mosaic, define ranges */ - mosaic[i].start = startTimesMosaic[i]; - mosaic[i].stop = 0; - mosaic[i].Mcal = 0.0; - mosaic[i].dMcal = 0.0; - mosaic[i].dMsys = 0.0; - mosaic[i].Xm = 0.0; - mosaic[i].flags = 0; - mosaic[i].secz = NAN; + mosaic[i].start = startTimesMosaic[i]; + mosaic[i].stop = 0; + mosaic[i].Mcal = 0.0; + mosaic[i].dMcal = 0.0; + mosaic[i].dMsys = 0.0; + mosaic[i].Xm = 0.0; + mosaic[i].flags = 0; + mosaic[i].secz = NAN; mosaic[i].photcode = 0; - mosaic[i].skipCal = FALSE; + mosaic[i].skipCal = FALSE; memset (&mosaic[i].coords, 0, sizeof(Coords)); @@ -238,9 +236,5 @@ } - /* select valid mosaic images by photcode */ - pname = GetPhotcodeNamebyCode (image[i].photcode); - if (!pname) continue; - status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME)); - if (status) continue; + if (!isGPC1chip(image[i].photcode)) continue; /* set image time range */ @@ -262,11 +256,12 @@ // assign each image to a mosaic + int Nsimple = 0; for (i = 0; i < Nsubset; i++) { ImageToMosaic[i] = -1; - /* select valid mosaic images by photcode */ - pname = GetPhotcodeNamebyCode (subset[i].photcode); - status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME)); - if (status) continue; + if (!isGPC1chip(subset[i].photcode)) { + Nsimple ++; + continue; + } start = subset[i].tzero; @@ -299,11 +294,11 @@ abort(); } - mosaic[j].stop = stop; - mosaic[j].Mcal = 0.0; - mosaic[j].dMcal = 0.0; - mosaic[j].Xm = 0.0; - mosaic[j].dMsys = subset[i].flags; - mosaic[j].flags = subset[i].flags; - mosaic[j].secz = subset[i].secz; + mosaic[j].stop = stop; + mosaic[j].Mcal = 0.0; + mosaic[j].dMcal = 0.0; + mosaic[j].Xm = 0.0; + mosaic[j].dMsys = subset[i].flags; + mosaic[j].flags = subset[i].flags; + mosaic[j].secz = subset[i].secz; mosaic[j].photcode = GetPhotcodeEquivCodebyCode (subset[i].photcode); } @@ -316,17 +311,15 @@ initMosaicGrid (subset, Nsubset); - fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nsubset, (int) Nmosaic); + fprintf (stderr, "matched %d images to %d mosaics, %d simple chips not matched to mosaics\n", (int) (Nsubset - Nsimple), (int) Nmosaic, (int) Nsimple); return; } -/* find mosaic frames (unique time periods) (NOTE : we do NOT require matching photcodes but we - match images by MOSAICNAME. this last point is weak: it forces a single camera at a time. - we can extend the logic to multiple cameras if we make list of MOSAICNAMES (better to assign a camera ID) +/* find mosaic frames (unique time periods) (NOTE : require gpc1 chips, which is pretty limiting) + if mergeMcal is TRUE, values will be saved on Mosaic.Mcal */ -void makeMosaics (Image *image, off_t Nimage) { - - off_t i, j, status, found, NMOSAIC, *MosaicN_IMAGE; +void makeMosaics (Image *image, off_t Nimage, int mergeMcal) { + + off_t i, j, found, NMOSAIC, *MosaicN_IMAGE; unsigned int start, stop, *startTimes, *startTimesMosaic; - char *pname; if (!MOSAIC_ZEROPT) return; @@ -341,8 +334,11 @@ // generate a list of all image start times ALLOCATE (startTimes, unsigned int, Nimage); + int Nmoschip = 0; for (i = 0; i < Nimage; i++) { - startTimes[i] = image[i].tzero; - } - sort_times (startTimes, Nimage); + if (!isGPC1chip(image[i].photcode)) continue; + startTimes[Nmoschip] = image[i].tzero; + Nmoschip ++; + } + sort_times (startTimes, Nmoschip); MARKTIME("create array of all image obstimes: %f sec\n", dtime); @@ -353,9 +349,6 @@ // generate a list of the unique start times (these define the mosaics) - for (i = 0; i < Nimage; i++) { - if (startTimes[i] < startTimesMosaic[Nmosaic]) { - fprintf (stderr, "error?\n"); - abort(); - } + for (i = 0; i < Nmoschip; i++) { + myAssert (startTimes[i] >= startTimesMosaic[Nmosaic], "times out of order?"); if (startTimes[i] == startTimesMosaic[Nmosaic]) continue; Nmosaic ++; @@ -379,14 +372,14 @@ // init the mosaic array values for (i = 0; i < Nmosaic; i++) { - mosaic[i].start = startTimesMosaic[i]; - mosaic[i].stop = 0; - mosaic[i].Mcal = 0.0; - mosaic[i].dMcal = 0.0; - mosaic[i].dMsys = 0.0; - mosaic[i].Xm = 0.0; - mosaic[i].flags = 0; - mosaic[i].secz = NAN; + mosaic[i].start = startTimesMosaic[i]; + mosaic[i].stop = 0; + mosaic[i].Mcal = 0.0; + mosaic[i].dMcal = 0.0; + mosaic[i].dMsys = 0.0; + mosaic[i].Xm = 0.0; + mosaic[i].flags = 0; + mosaic[i].secz = NAN; mosaic[i].photcode = 0; - mosaic[i].skipCal = FALSE; + mosaic[i].skipCal = FALSE; memset (&mosaic[i].coords, 0, sizeof(Coords)); @@ -401,13 +394,12 @@ // assign each image to a mosaic + int Nsimple = 0; for (i = 0; i < Nimage; i++) { ImageToMosaic[i] = -1; - /* select valid mosaic images by photcode */ - pname = GetPhotcodeNamebyCode (image[i].photcode); - if (!pname) continue; - - status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME)); - if (status) continue; + if (!isGPC1chip(image[i].photcode)) { + Nsimple ++; + continue; + } start = image[i].tzero; @@ -440,11 +432,11 @@ abort(); } - mosaic[j].stop = stop; - mosaic[j].Mcal = 0.0; - mosaic[j].dMcal = 0.0; - mosaic[j].Xm = 0.0; - mosaic[j].dMsys = image[i].flags; - mosaic[j].flags = image[i].flags; - mosaic[j].secz = image[i].secz; + mosaic[j].stop = stop; + mosaic[j].Mcal = 0.0; + mosaic[j].dMcal = 0.0; + mosaic[j].Xm = 0.0; + mosaic[j].dMsys = image[i].flags; + mosaic[j].flags = image[i].flags; + mosaic[j].secz = image[i].secz; mosaic[j].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode); } @@ -456,6 +448,9 @@ free (startTimesMosaic); + if (mergeMcal) { + initMosaicGrid (image, Nimage); + } + fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nimage, (int) Nmosaic); - return; } @@ -513,5 +508,4 @@ off_t i, j, m, NX, NY, NC, Nc; double R, D, Rmid, Dmid; - double Mcal, dMcal, Xm; double *Rc, *Dc; @@ -523,5 +517,4 @@ Nc = 0; Rmid = Dmid = NAN; - Mcal = dMcal = Xm = 0; for (j = 0; j < MosaicN_Image[i]; j++) { m = MosaicToImage[i][j]; @@ -563,15 +556,4 @@ REALLOCATE (Dc, double, NC); } - - Mcal += image[m].Mcal; - dMcal += image[m].dMcal; - Xm += image[m].Xm; - - // for ubercal images, we (elsewhere) keep Mcal frozen - - /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */ - image[m].Mcal = 0.0; - image[m].dMcal = NAN; - image[m].Xm = NAN_S_SHORT; } @@ -617,7 +599,7 @@ mosaic[i].coords.pc2_1 = 0.0; - mosaic[i].Mcal = Mcal / MosaicN_Image[i]; - mosaic[i].dMcal = dMcal / MosaicN_Image[i]; - mosaic[i].Xm = Xm / MosaicN_Image[i]; + mosaic[i].Mcal = 0.0; + mosaic[i].dMcal = 0.0; + mosaic[i].Xm = 0.0; } return; @@ -635,4 +617,6 @@ double R, D, Rmin, Rmax, Dmin, Dmax; double Mcal, dMcal, Xm; + + fprintf (stderr, "*** moving Mcal from image.Mcal to mosaic.Mcal ***\n"); dXmax = dYmax = 0.0; @@ -667,19 +651,19 @@ Dmin = MIN (Dmin, D); Dmax = MAX (Dmax, D); + + /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */ + + // XXX: how does this work with UBERCAL? We want to keep the Mcal values supplied by ubercal, but + // solve for a single offset for each exposure (Mosaic.Mcal). + // we also want to keep the flat-field terms for each exposure (regardless of ubercal or not) + // if it helps, note that ubercal uses a single zp per exposure, so the mean of those values is the same as the value + Mcal += image[m].Mcal; dMcal += image[m].dMcal; Xm += image[m].Xm; - // XXX: how does this work with UBERCAL? We want to keep the Mcal values supplied by ubercal, but - // solve for a single offset for each exposure (Mosaic.Mcal). - - // we also want to keep the flat-field terms for each exposure (regardless of ubercal or not) - - // if it helps, note that ubercal uses a single zp per exposure, so the mean of those values is the same as the value - - /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */ - image[m].Mcal = 0.0; - image[m].dMcal = NAN; - image[m].Xm = NAN_S_SHORT; + image[m].Mcal = 0.0; + image[m].dMcal = NAN; + image[m].Xm = NAN_S_SHORT; } dS /= MosaicN_Image[i]; @@ -697,7 +681,7 @@ RD_to_XY (&dX, &dY, Rmax, Dmax, &mosaic[i].coords); - mosaic[i].Mcal = Mcal / MosaicN_Image[i]; - mosaic[i].dMcal = dMcal / MosaicN_Image[i]; - mosaic[i].Xm = Xm / MosaicN_Image[i]; + mosaic[i].Mcal = Mcal / MosaicN_Image[i]; + mosaic[i].dMcal = dMcal / MosaicN_Image[i]; + mosaic[i].Xm = Xm / MosaicN_Image[i]; } if (!USE_GRID) return; @@ -709,5 +693,4 @@ } -// XXX : what about mosaics with skipCal == TRUE? void setMcalFinal () { @@ -719,6 +702,5 @@ image = getimages (&Nimage, NULL); - // XXX I think this is OK in the ubercal context, but probably need to skip UBERCAL - // images? (no need to update them) + fprintf (stderr, "*** return Mcal from mosaic.Mcal to image.Mcal ***\n"); // copy the mosaic results to the images. set the mosaic Mcal to 0.0 since we have moved its @@ -840,6 +822,7 @@ mosID = ImageToMosaic[idx]; if (mosID < 0) { - Image *image = getimage(idx); - fprintf (stderr, "unmatched image %s\n", image[0].name); + // Image *image = getimage(idx); + // fprintf (stderr, "unmatched image %s\n", image[0].name); + // skip measurements from simple chips (not matched to a mosaic by definition) return; } @@ -882,6 +865,7 @@ if (!MOSAIC_ZEROPT) return (0); + // unassigned measurements belong to simple chips i = MeasureToMosaic[cat][meas]; - if (i == -1) return (NAN); + if (i == -1) return (0.0); if (mosaic[i].flags & IMAGE_BAD) return (NAN); @@ -1127,5 +1111,5 @@ off_t n = catalog[c].measureT[m].averef; - Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]); + Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF); float delta = Msys - Mrel - Mcal - Mgrid + Mflat; @@ -1164,5 +1148,5 @@ off_t n = catalog[c].measureT[m].averef; - Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]); + Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF); if (isnan(Msys)) { info->Nsys++; @@ -1239,5 +1223,5 @@ if (testImage) { - fprintf (stderr, "%f %f : %d %f\n", myMosaic[0].Mcal, myMosaic[0].dMsys, myMosaic[0].Xm, pow(10.0, 0.01*myMosaic[0].Xm)); + fprintf (stderr, "%f %f : %f\n", myMosaic[0].Mcal, myMosaic[0].dMsys, pow (10.0, 0.01*myMosaic[0].Xm)); } @@ -1592,4 +1576,6 @@ if (mosaic[i].flags & IMAGE_BAD) continue; if (mosaic[i].skipCal) continue; + if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue; + list[n] = mosaic[i].Mcal; dlist[n] = 1; @@ -1622,4 +1608,6 @@ if (mosaic[i].flags & IMAGE_BAD) continue; if (mosaic[i].skipCal) continue; + if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue; + list[n] = mosaic[i].dMcal; dlist[n] = 1; @@ -1653,4 +1641,5 @@ if (mosaic[i].flags & IMAGE_BAD) continue; if (mosaic[i].skipCal) continue; + if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue; N = 0; @@ -1699,5 +1688,7 @@ if (mosaic[i].flags & IMAGE_BAD) continue; if (mosaic[i].skipCal) continue; - list[n] = pow(10.0, 0.01*mosaic[i].Xm); + if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue; + + list[n] = pow (10.0, 0.01*mosaic[i].Xm); dlist[n] = 1; n++; @@ -1782,7 +1773,6 @@ void plot_mosaic_fields (Catalog *catalog) { - off_t i, j, m, c, N, ave, Nimage; + off_t i, j, m, c, N, Nimage; double *xlist, *ylist; - // double Xmin, Xmax, Ymin, Ymax; char string[64]; Graphdata graphdata; @@ -1803,6 +1793,4 @@ for (i = 0; i < Nmosaic; i++) { N = 0; - // Xmin = Ymin = +360.0; - // Xmax = Ymax = -360.0; for (j = 0; j < N_onMosaic[i]; j++) { @@ -1812,7 +1800,7 @@ 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; + // ave = catalog[c].measureT[m].averef; + xlist[N] = catalog[c].measureT[m].R; + ylist[N] = catalog[c].measureT[m].D; N++; } Index: /trunk/Ohana/src/relphot/src/Shutdown.c =================================================================== --- /trunk/Ohana/src/relphot/src/Shutdown.c (revision 37036) +++ /trunk/Ohana/src/relphot/src/Shutdown.c (revision 37037) @@ -23,5 +23,5 @@ SetProtect (TRUE); if (db) gfits_db_close (db); - fprintf (stderr, "ERROR: addstar halted\n"); + fprintf (stderr, "ERROR: relphot halted\n"); exit (1); } Index: /trunk/Ohana/src/relphot/src/StarOps.c =================================================================== --- /trunk/Ohana/src/relphot/src/StarOps.c (revision 37036) +++ /trunk/Ohana/src/relphot/src/StarOps.c (revision 37037) @@ -153,4 +153,7 @@ int i; + // fprintf (stderr, "this version of setMrel is invalid for now\n"); + // exit (1); + int Nsecfilt = GetPhotcodeNsecfilt (); @@ -407,5 +410,5 @@ /* star/photcodes already calibrated */ - if ( isSetMrelFinal && catalog[Nc].found[Nsecfilt*j+Nsec]) continue; + if ( isSetMrelFinal && catalog[Nc].found_t[Nsecfilt*j+Nsec]) continue; // skip bad stars @@ -467,5 +470,5 @@ // skip some absurd values NAN, < 0.0, > 30.0 - Msys = PhotSysTiny (&catalog[Nc].measureT[m], &catalog[Nc].averageT[j], &catalog[Nc].secfilt[j*Nsecfilt]); + Msys = PhotSysTiny (&catalog[Nc].measureT[m], &catalog[Nc].averageT[j], &catalog[Nc].secfilt[j*Nsecfilt], MAG_CLASS_PSF); if (isnan(Msys)) SKIP_THIS_MEAS(Nsys); if (Msys < 0.0) SKIP_THIS_MEAS(Nsys); @@ -670,7 +673,7 @@ liststats (list, dlist, wlist, N, &stats); - catalog[Nc].secfilt[Nsecfilt*j+Nsec].M = stats.mean; - catalog[Nc].secfilt[Nsecfilt*j+Nsec].dM = stats.error; - catalog[Nc].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq + 1e-4) : NAN_S_SHORT; + catalog[Nc].secfilt[Nsecfilt*j+Nsec].M = stats.mean; + catalog[Nc].secfilt[Nsecfilt*j+Nsec].dM = stats.error; + catalog[Nc].secfilt[Nsecfilt*j+Nsec].Mchisq = (stats.Nmeas > 1) ? stats.chisq : NAN; // when running -averages, we have no information about the images, so we cannot set this @@ -680,5 +683,5 @@ if (isSetMrelFinal) { - catalog[Nc].found[Nsecfilt*j+Nsec] = TRUE; + catalog[Nc].found_t[Nsecfilt*j+Nsec] = TRUE; catalog[Nc].secfilt[Nsecfilt*j+Nsec].Mstdev = 1000.0*stats.sigma; // Mstdev is in millimags (not enough space for more precision) @@ -950,6 +953,6 @@ int i, j, Ndel, Nave, Ntot, mark, Ns, Nscat, Nchi, Nnan; - float dM, Xm; - double Chisq, MaxScatter, MaxChisq; + float dM; + double MaxScatter, MaxChisq; double *xlist, *slist, *dlist; @@ -978,8 +981,7 @@ 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; + float Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq; + if (isnan(Mchisq)) continue; + xlist[Ntot] = Mchisq; slist[Ntot] = catalog[i].secfilt[Nsecfilt*j+Nsec].dM; dlist[Ntot] = 1; @@ -999,13 +1001,12 @@ 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); + float Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq; + mark = (dM > MaxScatter) || (isnan(Mchisq)) || (Mchisq > MaxChisq); if (mark) { catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_POOR; Ndel ++; - if (dM > MaxScatter) { Nscat ++; } - if (Xm == NAN_S_SHORT) { Nnan ++; } - if (Chisq > MaxChisq) { Nchi ++; } + if (dM > MaxScatter) { Nscat ++; } + if (isnan(Mchisq)) { Nnan ++; } + if (Mchisq > MaxChisq) { Nchi ++; } } else { catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_POOR; @@ -1072,5 +1073,5 @@ /* on final processing, skip stars already measured */ - if (final && catalog[i].found[Nsecfilt*j + Nsec]) continue; + if (final && catalog[i].found_t[Nsecfilt*j + Nsec]) continue; /* skip bad stars to prevent them from becoming good (on inner sample) */ @@ -1093,5 +1094,5 @@ if (isnan(Mgrid)) { Ngrid ++; continue; } - Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt]); + Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt], MAG_CLASS_PSF); list[N] = Msys - Mcal - Mmos - Mgrid; dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR); @@ -1138,5 +1139,5 @@ if (isnan(Mgrid)) continue; - Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt]); + Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt], MAG_CLASS_PSF); list[N] = Msys - Mcal - Mmos - Mgrid; dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR); @@ -1233,5 +1234,5 @@ off_t j, Ntot; - int i, n, Xm; + int i, n; double *list, *dlist; StatType stats; @@ -1254,7 +1255,7 @@ 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); + float Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq; + if (isnan(Mchisq)) continue; + list[n] = Mchisq; dlist[n] = 1; n++; @@ -1377,7 +1378,7 @@ 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; + value = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq; + if (isnan(value)) continue; + ylist[N] = value; N++; } Index: /trunk/Ohana/src/relphot/src/args.c =================================================================== --- /trunk/Ohana/src/relphot/src/args.c (revision 37036) +++ /trunk/Ohana/src/relphot/src/args.c (revision 37037) @@ -227,8 +227,8 @@ } - SHOW_PARAMS = FALSE; + SHOW_PARAMS = TRUE; if ((N = get_argument (argc, argv, "-params"))) { remove_argument (N, &argc, argv); - SHOW_PARAMS = TRUE; + SHOW_PARAMS = FALSE; } Index: /trunk/Ohana/src/relphot/src/assign_images.c =================================================================== --- /trunk/Ohana/src/relphot/src/assign_images.c (revision 37036) +++ /trunk/Ohana/src/relphot/src/assign_images.c (revision 37037) @@ -25,7 +25,7 @@ if (MOSAIC_ZEROPT) { - makeMosaics (image, Nimage); - - // center coords and Mcal, dMcal, Xm for the mosaics + makeMosaics (image, Nimage, FALSE); + + // center coords and Mcal, dMcal, Mchisq for the mosaics setMosaicCenters (image, Nimage); MARKTIME("set mosaic coordinates and Mcal values: %f sec\n", dtime); Index: /trunk/Ohana/src/relphot/src/bcatalog.c =================================================================== --- /trunk/Ohana/src/relphot/src/bcatalog.c (revision 37036) +++ /trunk/Ohana/src/relphot/src/bcatalog.c (revision 37037) @@ -23,4 +23,7 @@ Nmeasure = Naverage = 0; + int NmMax = 0; + int NmMin = catalog[0].Nmeasure; + Ncode = Ntime = Ndophot = Nmag = Nsigma = Nimag = Nfew = Npsfqf = Ngalaxy = Nnan = Nbad = Npoor = 0; @@ -55,4 +58,6 @@ if (RESET) { int Ns; + +# if (0) DVOSecfiltFlags secfiltBits = ID_SECF_STAR_FEW | @@ -66,4 +71,6 @@ ID_SECF_USE_UBERCAL | ID_SECF_OBJ_EXT; +# endif + for (Ns = 0; Ns < Nphotcodes; Ns++) { @@ -71,15 +78,18 @@ int Nsec = GetPhotcodeNsec(thisCode); - subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M = NAN; - subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Map = NAN; - subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].dM = NAN; - subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Mstdev = NAN_S_SHORT; - subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Xm = NAN; - subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M_20 = NAN_S_SHORT; - subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M_80 = NAN_S_SHORT; - subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Ncode = 0; - subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Nused = 0; - subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].flags &= ~secfiltBits; - subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].ubercalDist = 1000; + fprintf (stderr, "really use dvo_secfilt_init?"); + dvo_secfilt_init (&subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec]); + + // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M = NAN; + // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Map = NAN; + // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].dM = NAN; + // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Mstdev = NAN_S_SHORT; + // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Mchisq = NAN; + // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M_20 = NAN_S_SHORT; + // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M_80 = NAN_S_SHORT; + // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Ncode = 0; + // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Nused = 0; + // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].flags &= ~secfiltBits; + // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].ubercalDist = 1000; } } @@ -136,5 +146,5 @@ /* select measurements by mag limit */ - mag = PhotCat (&catalog[0].measure[offset]); + mag = PhotCat (&catalog[0].measure[offset], MAG_CLASS_PSF); if (mag > MAG_LIM) { Nmag ++; continue; } @@ -144,5 +154,5 @@ /* select measurements by mag limit */ if (ImagSelect) { - mag = PhotInst (&catalog[0].measure[offset]); + mag = PhotInst (&catalog[0].measure[offset], MAG_CLASS_PSF); if (mag < ImagMin) { Nimag ++; continue; } if (mag > ImagMax) { Nimag ++; continue; } @@ -180,5 +190,5 @@ } myAssert (Nmeasure < NMEASURE, "realloc failure"); - } + } // end of catalog.average.Nmeasure loop // skip object if it is likely to be a galaxy @@ -205,4 +215,6 @@ continue; } + NmMin = MIN (Nm, NmMin); + NmMax = MAX (Nm, NmMax); // for w-band photometry (& other cases?) convert gri photometry @@ -235,4 +247,11 @@ subcatalog[0].Nsecf_mem = Naverage * catalog[0].Nsecfilt; + if (VERBOSE) { + fprintf (stderr, "using "OFF_T_FMT" stars ("OFF_T_FMT" measures) of "OFF_T_FMT" ("OFF_T_FMT" measures) for catalog %s, %d < Nm < %d\n", + subcatalog[0].Naverage, subcatalog[0].Nmeasure, catalog[0].Naverage, catalog[0].Nmeasure, catalog[0].filename, NmMin, NmMax); + fprintf (stderr, "rejections: %d stars have too few measures:\n %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d psfqf, %d Nnan, %d galaxies, %d bad, %d poor\n", + Nfew, Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Npsfqf, Nnan, Ngalaxy, Nbad, Npoor); + } + // limit the total number of stars in the catalog if (MaxDensityUse) { @@ -240,10 +259,4 @@ } - if (VERBOSE) { - fprintf (stderr, "using "OFF_T_FMT" stars ("OFF_T_FMT" measures) of "OFF_T_FMT" for catalog %s\n", - subcatalog[0].Naverage, subcatalog[0].Nmeasure, i, catalog[0].filename); - fprintf (stderr, "rejections: %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d few, %d psfqf, %d Nnan, %d galaxies, %d bad, %d poor\n", - Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Npsfqf, Nnan, Ngalaxy, Nbad, Npoor); - } return (TRUE); } Index: /trunk/Ohana/src/relphot/src/extra.c =================================================================== --- /trunk/Ohana/src/relphot/src/extra.c (revision 37037) +++ /trunk/Ohana/src/relphot/src/extra.c (revision 37037) @@ -0,0 +1,40 @@ +# include "relphot.h" + +// for now (20140710) I need to identify gpc1 chips explicitly. generalize in the future +int isGPC1chip (int photcode) { + + if ((photcode > 10000) && (photcode < 10077)) return TRUE; // g-band + if ((photcode > 10100) && (photcode < 10177)) return TRUE; // r-band + if ((photcode > 10200) && (photcode < 10277)) return TRUE; // i-band + if ((photcode > 10300) && (photcode < 10377)) return TRUE; // z-band + if ((photcode > 10400) && (photcode < 10477)) return TRUE; // y-band + if ((photcode > 10500) && (photcode < 10577)) return TRUE; // w-band + + return FALSE; +} + +// for now (20140710) I need to identify gpc1 stacks explicitly. generalize in the future +int isGPC1stack (int photcode) { + + if (photcode == 11000) return TRUE; // g-band + if (photcode == 11100) return TRUE; // r-band + if (photcode == 11200) return TRUE; // i-band + if (photcode == 11300) return TRUE; // z-band + if (photcode == 11400) return TRUE; // y-band + if (photcode == 11500) return TRUE; // w-band + + return FALSE; +} + +// for now (20140710) I need to identify gpc1 stacks explicitly. generalize in the future +int isGPC1warp (int photcode) { + + if (photcode == 12000) return TRUE; // g-band + if (photcode == 12100) return TRUE; // r-band + if (photcode == 12200) return TRUE; // i-band + if (photcode == 12300) return TRUE; // z-band + if (photcode == 12400) return TRUE; // y-band + if (photcode == 12500) return TRUE; // w-band + + return FALSE; +} Index: /trunk/Ohana/src/relphot/src/initialize.c =================================================================== --- /trunk/Ohana/src/relphot/src/initialize.c (revision 37036) +++ /trunk/Ohana/src/relphot/src/initialize.c (revision 37037) @@ -36,5 +36,14 @@ if (SHOW_PARAMS) { - fprintf (stderr, "current parameter settings:\n"); + int Ns; + fprintf (stderr, "subset selection criteria:\n"); + fprintf (stderr, " photcodes "); + for (Ns = 0; Ns < Nphotcodes; Ns++) { + if (Ns == Nphotcodes - 1) { + fprintf (stderr, "%s\n", photcodes[Ns][0].name); + } else { + fprintf (stderr, "%s, ", photcodes[Ns][0].name); + } + } if (TimeSelect) { fprintf (stderr, "TimeSelect: TRUE (%s - %s)\n", ohana_sec_to_date (TSTART), ohana_sec_to_date (TSTOP)); @@ -42,4 +51,18 @@ fprintf (stderr, "TimeSelect: FALSE\n"); } + if (DophotSelect) { + fprintf (stderr, "DophotSelect: TRUE (%d)\n", DophotValue); + } else { + fprintf (stderr, "DophotSelect: FALSE\n"); + } + fprintf (stderr, "PSF_QF limit: 0.85 (hardwired)\n"); + + // fprintf (stderr, "Photom Bad Mask: 0x%08x, Photom Poor Mask: 0x%08x\n"); + + fprintf (stderr, "MAG_LIM: %f, SIGMA_LIM: %f\n", MAG_LIM, SIGMA_LIM); + fprintf (stderr, "INST_MAG_MIN: %f, INST_MAG_MAX: %f\n", ImagMin, ImagMax); + + fprintf (stderr, "STAR_TOOFEW: %d\n", STAR_TOOFEW); + fprintf (stderr, "VERBOSE: %d, PLOTSTUFF: %d\n", VERBOSE, PLOTSTUFF); fprintf (stderr, "GRID_X: %d, GRID_Y: %d, BINNING: %d == Nmx: %d, Nmy: %d\n", @@ -49,5 +72,4 @@ (int)(RELPHOT_GRID_X/RELPHOT_GRID_BINNING), (int)(RELPHOT_GRID_Y/RELPHOT_GRID_BINNING)); - fprintf (stderr, "MAG_LIM %lf\n", MAG_LIM); fprintf (stderr, "STAR_SCATTER %lf\n", STAR_SCATTER); fprintf (stderr, "IMAGE_SCATTER %lf\n", IMAGE_SCATTER); @@ -56,5 +78,4 @@ fprintf (stderr, "GSCFILE %s\n", GSCFILE); fprintf (stderr, "CATDIR %s\n", CATDIR); - exit (0); } Index: /trunk/Ohana/src/relphot/src/launch_region_hosts.c =================================================================== --- /trunk/Ohana/src/relphot/src/launch_region_hosts.c (revision 37036) +++ /trunk/Ohana/src/relphot/src/launch_region_hosts.c (revision 37037) @@ -55,6 +55,15 @@ char command[1024]; - snprintf (command, 1024, "relphot %s -parallel-images %s -region-hosts %s -region-hostID %d -D CATDIR %s -region %f %f %f %f -statmode %s -D CAMERA %s -D STAR_TOOFEW %d -minerror %f", - PhotcodeList, filename, REGION_FILE, host->hostID, CATDIR, host->RminCat, host->RmaxCat, host->DminCat, host->DmaxCat, STATMODE, CAMERA, STAR_TOOFEW, MIN_ERROR); + snprintf (command, 1024, "relphot %s", PhotcodeList); + strextend (command, "-parallel-images %s", filename); + strextend (command, "-region-hosts %s", REGION_FILE); + strextend (command, "-region-hostID %d", host->hostID); + strextend (command, "-D CATDIR %s", CATDIR); + strextend (command, "-region %f %f %f %f", host->RminCat, host->RmaxCat, host->DminCat, host->DmaxCat); + strextend (command, "-statmode %s", STATMODE); + strextend (command, "-D CAMERA %s", CAMERA); + strextend (command, "-D STAR_TOOFEW %d", STAR_TOOFEW); + strextend (command, "-minerror %f", MIN_ERROR); + strextend (command, "-cloud-limit %f", CLOUD_TOLERANCE); if (VERBOSE) strextend (command, "-v"); Index: /trunk/Ohana/src/relphot/src/load_catalogs.c =================================================================== --- /trunk/Ohana/src/relphot/src/load_catalogs.c (revision 37036) +++ /trunk/Ohana/src/relphot/src/load_catalogs.c (revision 37037) @@ -12,5 +12,5 @@ // if this function is called with a specified hostID, then only the fraction of the // database hosted by that hostID is loaded -Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int hostID, char *hostpath) { +Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int hostID, char *hostpath, char *syncfile) { off_t i, Nmeas, Nstar, Nmeas_total, Nstar_total; @@ -19,5 +19,5 @@ // XXX need to decide how to determine PARALLEL mode... if (PARALLEL && !hostID) { - catalog = load_catalogs_parallel (skylist, Ncatalog); + catalog = load_catalogs_parallel (skylist, Ncatalog, syncfile); return catalog; } @@ -92,4 +92,10 @@ // in regular relphot, we shutdown here; in relphot_client, we generate and return an empty table (for consistency) + // if we are running with parallel_images but not a parallel database, we need to + // release the lock so the next image host can proceed + if (!hostID && syncfile) { + update_sync_file (syncfile, 1); + } + // XXX consider only returning the populated catalogs *Ncatalog = skylist[0].Nregions; @@ -105,5 +111,5 @@ // CATDIR is supplied globally # define DEBUG 1 -Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog) { +Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog, char *syncfile) { char uniquer[12]; @@ -206,4 +212,17 @@ } + // update syncfile here (save lots of I/O time): + + // at this point, the remote relastro_client jobs are done loading their data. in a + // parallel_images mode, the next image host can be launched while this image host now + // reads that + + // NOTE: if I let all hosts load blindly, I saturate the data clients with too many + // relphot_client requests. I need to have the master mediate this. the master + // will not launch the next remote job until this one says it is done + if (syncfile) { + update_sync_file (syncfile, 1); + } + // each host generates a BrightCatalog structure, with the measure, average, etc value // loaded into a single set of arrays (of MeasureTiny, AverageTiny, Secfilt). I need to Index: /trunk/Ohana/src/relphot/src/plot_scatter.c =================================================================== --- /trunk/Ohana/src/relphot/src/plot_scatter.c (revision 37036) +++ /trunk/Ohana/src/relphot/src/plot_scatter.c (revision 37037) @@ -51,6 +51,6 @@ 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]); + ylist[N] = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt], MAG_CLASS_PSF) - Mcal - Mmos - Mgrid - Mrel; + ilist[N] = PhotInstTiny (&catalog[i].measureT[m], MAG_CLASS_PSF); N++; } Index: /trunk/Ohana/src/relphot/src/relphot_client.c =================================================================== --- /trunk/Ohana/src/relphot/src/relphot_client.c (revision 37036) +++ /trunk/Ohana/src/relphot/src/relphot_client.c (revision 37037) @@ -39,5 +39,5 @@ case MODE_LOAD: { int Ncatalog; - Catalog *catalog = load_catalogs (skylist, &Ncatalog, HOST_ID, HOSTDIR); + Catalog *catalog = load_catalogs (skylist, &Ncatalog, HOST_ID, HOSTDIR, NULL); if (!catalog) { fprintf (stderr, "ERROR loading catalogs from %s\n", CATDIR); Index: /trunk/Ohana/src/relphot/src/relphot_images.c =================================================================== --- /trunk/Ohana/src/relphot/src/relphot_images.c (revision 37036) +++ /trunk/Ohana/src/relphot/src/relphot_images.c (revision 37037) @@ -52,5 +52,5 @@ if (NLOOP > 0) { /* load catalog data from region files (hostID is 0 since we are not a client */ - catalog = load_catalogs (skylist, &Ncatalog, 0, NULL); + catalog = load_catalogs (skylist, &Ncatalog, 0, NULL, NULL); MARKTIME("-- load catalog data: %f sec\n", dtime); Index: /trunk/Ohana/src/relphot/src/relphot_parallel_images.c =================================================================== --- /trunk/Ohana/src/relphot/src/relphot_parallel_images.c (revision 37036) +++ /trunk/Ohana/src/relphot/src/relphot_parallel_images.c (revision 37037) @@ -38,6 +38,5 @@ // unlink (IMAGE_TABLE); // XXX a bit risky, add some protection? - // XXX need to deal with mosaic vs image... - makeMosaics (image, Nimage); + makeMosaics (image, Nimage, TRUE); initImages (image, NULL, Nimage); @@ -56,14 +55,11 @@ /* load catalog data from region files (hostID is 0 since we are not a client */ - catalog = load_catalogs (skylist, &Ncatalog, 0, NULL); + char *syncfile = make_filename (CATDIR, regionHosts->hosts[myHost].hostname, REGION_HOST_ID, "loadcat.sync"); + catalog = load_catalogs (skylist, &Ncatalog, 0, NULL, syncfile); MARKTIME("-- load catalog data: %f sec\n", dtime); + free (syncfile); + client_logger_message ("loaded catalog data\n"); - // NOTE: if I let all hosts load blindly, I saturate the data clients with too many - // relphot_client requests. I need to have the master mediate this. the master - // will not launch the next remote job until this one says it is done - char *syncfile = make_filename (CATDIR, regionHosts->hosts[myHost].hostname, REGION_HOST_ID, "loadcat.sync"); - update_sync_file (syncfile, 1); - // generate tables go from catID,objID -> catSeq,objSeq indexCatalogs (catalog, Ncatalog); Index: /trunk/Ohana/src/relphot/src/relphot_parallel_regions.c =================================================================== --- /trunk/Ohana/src/relphot/src/relphot_parallel_regions.c (revision 37036) +++ /trunk/Ohana/src/relphot/src/relphot_parallel_regions.c (revision 37037) @@ -1,3 +1,4 @@ # include "relphot.h" +int save_images_backup (FITS_DB *db); int relphot_parallel_regions () { @@ -63,4 +64,7 @@ SkyList *skylist = SkyListByBounds (sky, -1, regionHosts->Rmin, regionHosts->Rmax, regionHosts->Dmin, regionHosts->Dmax); + // I have to save a copy because dvo_image_save and _unlock swap and free the data + save_images_backup (&db); + /* update catalogs (in parallel) */ reload_catalogs (skylist, flatcorr, 0, NULL); @@ -74,2 +78,37 @@ } +int save_images_backup (FITS_DB *db) { + + FITS_DB dbX; + char filename[1024]; + + INITTIME; + + gfits_db_init (&dbX); + dbX.lockstate = LCK_XCLD; + dbX.timeout = 60.0; + dbX.mode = db->mode; + dbX.format = db->format; + + snprintf (filename, 1024, "%s.bck", ImageCat); + if (!gfits_db_lock (&dbX, filename)) { + fprintf (stderr, "can't lock backup image image catalog\n"); + return (FALSE); + } + + // copy Images.dat data (all or only the subset / vtable elements?) I think I need to dump + // the entire Image table, but with the updates in place I think this says: re-work the + // ftable/vtable usage in this program to be more sensible... + gfits_copy_header (&db->header, &dbX.header); + gfits_copy_matrix (&db->matrix, &dbX.matrix); + gfits_copy_header (&db->theader, &dbX.theader); + gfits_copy_ftable (&db->ftable, &dbX.ftable); + + dbX.ftable.header = &dbX.theader; + + dvo_image_save (&dbX, VERBOSE); + dvo_image_unlock (&dbX); + MARKTIME("-- save Image.dat.bck: %f sec\n", dtime); + + return TRUE; +} Index: /trunk/Ohana/src/relphot/src/setExclusions.c =================================================================== --- /trunk/Ohana/src/relphot/src/setExclusions.c (revision 37036) +++ /trunk/Ohana/src/relphot/src/setExclusions.c (revision 37037) @@ -35,6 +35,6 @@ /* select measurements by mag limit */ if (AreaSelect) { - 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; + r = catalog[i].measureT[m].R; + d = catalog[i].measureT[m].D; if ((coords = getCoords (m, i)) == NULL) goto markbad; RD_to_XY (&x, &y, r, d, coords); Index: /trunk/Ohana/src/relphot/src/setMrelCatalog.c =================================================================== --- /trunk/Ohana/src/relphot/src/setMrelCatalog.c (revision 37036) +++ /trunk/Ohana/src/relphot/src/setMrelCatalog.c (revision 37037) @@ -1,3 +1,4 @@ # include "relphot.h" +# define UBERCAL_WEIGHT 100.0 # define SKIP_THIS_MEAS(REASON) { \ @@ -13,4 +14,9 @@ continue; } +static float MagToFlux (float Mag) { + float Flux = pow(10.0, -0.4*(Mag)); + return (Flux); +} + int print_measure_set_alt (Average *average, SecFilt *secfilt, Measure *measure) { @@ -29,4 +35,6 @@ return (TRUE); } + +int setMrelAverageForcedWarp (off_t measureOffset, int cat, int pass, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t *found); int setMrel_catalog_alt (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) { @@ -52,5 +60,5 @@ Average *average = catalog[Nc].average ? &catalog[Nc].average[j] : NULL; Measure *measure = catalog[Nc].measure ? &catalog[Nc].measure[m] : NULL; - setMrelAverageExposure (m, Nc, pass, flatcorr, results, average, &catalog[Nc].averageT[j], &catalog[Nc].secfilt[j*Nsecfilt], measure, &catalog[Nc].measureT[m], &catalog[Nc].found[Nsecfilt*j]); + setMrelAverageExposure (m, Nc, pass, flatcorr, results, average, &catalog[Nc].averageT[j], &catalog[Nc].secfilt[j*Nsecfilt], measure, &catalog[Nc].measureT[m], &catalog[Nc].found_t[Nsecfilt*j]); // only apply Stack operation on setMrelFinal in first pass @@ -59,4 +67,10 @@ setGlobalObjStats (&catalog[Nc].average[j], &catalog[Nc].measure[m]); } + + // only measure force-warp mean values if issetMrelFinal (make it optional?) + if (isSetMrelFinal) { + setMrelAverageForcedWarp (m, Nc, pass, results, &catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], &catalog[Nc].measure[m], &catalog[Nc].foundWarp_t[Nsecfilt*j]); + setGlobalObjStats (&catalog[Nc].average[j], &catalog[Nc].measure[m]); + } } if (primaryCell) free (primaryCell); @@ -64,4 +78,5 @@ } +// set mean of chip measurements (selected by photcode range for now): int setMrelAverageExposure (off_t measureOffset, int cat, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT, off_t *found) { @@ -144,7 +159,11 @@ // are we a PS1 exposure photcode? - if ((measureT[k].photcode >= 10000) && (measureT[k].photcode <= 10500)) { - NexpPS1 ++; - } + if (isGPC1chip(measureT[k].photcode)) NexpPS1 ++; + + // SKIP gpc1 stack data + if (isGPC1stack(measureT[k].photcode)) continue; + + // SKIP gpc1 forced-warp data + if (isGPC1warp(measureT[k].photcode)) continue; if (measureT[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad); @@ -166,5 +185,5 @@ // skip some absurd values NAN, < 0.0, > 30.0 - Msys = PhotSysTiny (&measureT[k], &averageT[0], &secfilt[0]); + Msys = PhotSysTiny (&measureT[k], &averageT[0], &secfilt[0], MAG_CLASS_PSF); if (isnan(Msys)) SKIP_THIS_MEAS(Nsys); if (Msys < 0.0) SKIP_THIS_MEAS(Nsys); @@ -197,9 +216,4 @@ } havePS1 = TRUE; - } - - // gpc1 stack data - if ((measure[k].photcode >= 11000) && (measure[k].photcode <= 11400)) { - continue; } @@ -255,17 +269,17 @@ // Map (Maplist) and Mkron (Mkronlist,dkronlist) are used to calculate mean mags per filter - float Map = PhotAperCat (&measure[k]); + float Map = PhotCat (&measure[k], MAG_CLASS_APER); if (!isnan(Map)) { Maplist[Nap] = Map - Mcal - Mmos - Mgrid; - daplist[Nap] = dMpsf; - waplist[Nap] = isUbercal ? 10.0 : 1.0; + daplist[Nap] = dMpsf;// XXX check on this... + waplist[Nap] = isUbercal ? UBERCAL_WEIGHT : 1.0; Nap ++; } - float Mkron = PhotKronCat (&measure[k]); + float Mkron = PhotCat (&measure[k], MAG_CLASS_KRON); if (!isnan(Mkron)) { Mkronlist[Nkron] = Mkron - Mcal - Mmos - Mgrid; dkronlist[Nkron] = measure[k].dMkron; - wkronlist[Nkron] = isUbercal ? 10.0 : 1.0; + wkronlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0; Nkron ++; } @@ -280,5 +294,5 @@ Mpsflist[Npsf] = Msys - Mcal - Mmos - Mgrid; dpsflist[Npsf] = dMpsf; - wpsflist[Npsf] = isUbercal ? 10.0 : 1.0; + wpsflist[Npsf] = isUbercal ? UBERCAL_WEIGHT : 1.0; // NOTE: @@ -293,5 +307,5 @@ if (refPhotcode) { if (code->code == refPhotcode->code) { - wpsflist[Npsf] = 100.0; + wpsflist[Npsf] = UBERCAL_WEIGHT; } } @@ -304,5 +318,5 @@ if (isSetMrelFinal && (pass == 0)) { if ((thisCode < 6) || (thisCode == 9)) { - secfilt[Nsec].Ncode = NexpPS1; + secfilt[Nsec].Ncode = NexpPS1; } else { secfilt[Nsec].Ncode = Ncode; // 2MASS data if it exists @@ -334,7 +348,7 @@ liststats (Mpsflist, dpsflist, wpsflist, Npsf, psfstats); - secfilt[Nsec].M = psfstats->mean; - secfilt[Nsec].dM = psfstats->error; - secfilt[Nsec].Xm = (psfstats->Nmeas > 1) ? 100.0*log10(psfstats->chisq + 1e-4) : NAN_S_SHORT; + secfilt[Nsec].M = psfstats->mean; + secfilt[Nsec].dM = psfstats->error; + secfilt[Nsec].Mchisq = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN; // when running -averages, we have no information about the images, so we cannot set this @@ -346,18 +360,23 @@ found[Nsec] = TRUE; - secfilt[Nsec].Mstdev = 1000.0*psfstats->sigma; // Mstdev is in millimags (not enough space for more precision) + secfilt[Nsec].Mstdev = psfstats->sigma; // Mstdev is in millimags (not enough space for more precision) // secfilt[Nsec].Ncode = Ncode; secfilt[Nsec].Nused = psfstats->Nmeas; - secfilt[Nsec].M_80 = 1000 * psfstats->Upper80; - secfilt[Nsec].M_20 = 1000 * psfstats->Lower20; + secfilt[Nsec].Mmax = psfstats->max; + secfilt[Nsec].Mmin = psfstats->min; // NOTE : use the modified weight for apmags as well as psf mags liststats (Maplist, daplist, waplist, Nap, apstats); secfilt[Nsec].Map = Nap > 0 ? apstats->mean : NAN; + secfilt[Nsec].dMap = Nap > 0 ? apstats->error : NAN; + secfilt[Nsec].sMap = Nap > 0 ? apstats->sigma : NAN; + secfilt[Nsec].NusedAp = Nap; liststats (Mkronlist, dkronlist, wkronlist, Nkron, kronstats); secfilt[Nsec].Mkron = Nkron > 0 ? kronstats->mean : NAN; secfilt[Nsec].dMkron = Nkron > 0 ? kronstats->error : NAN; + secfilt[Nsec].sMkron = Nkron > 0 ? kronstats->sigma : NAN; + secfilt[Nsec].NusedKron = Nkron; // NOTE: for 2MASS measurements, Next should be 1, as should N @@ -472,7 +491,5 @@ // only examine gpc1 stack data - // XXX this is absurdly hardwired (along with several photcode tests) - if (measure[k].photcode < 11000) continue; - if (measure[k].photcode > 11400) continue; + if (!isGPC1stack(measure[k].photcode)) continue; haveStack = TRUE; @@ -493,5 +510,5 @@ if (isnan(Mgrid)) SKIP_THIS_MEAS_STACK(Ngrid); - Msys = PhotSys (&measure[k], &average[0], &secfilt[0]); + Msys = PhotSys (&measure[k], &average[0], &secfilt[0], MAG_CLASS_PSF); if (isnan(Msys)) SKIP_THIS_MEAS_STACK(Nsys); @@ -581,8 +598,14 @@ // need to put in AB mag factor to get to Janskies (or uJy?) - secfilt[Nsec].FluxPSF = zpFactor * measure[k].FluxPSF; - secfilt[Nsec].dFluxPSF = zpFactor * measure[k].dFluxPSF; - secfilt[Nsec].FluxKron = zpFactor * measure[k].FluxKron; - secfilt[Nsec].dFluxKron = zpFactor * measure[k].dFluxKron; + secfilt[Nsec].FpsfStk = zpFactor * measure[k].FluxPSF; + secfilt[Nsec].dFpsfStk = zpFactor * measure[k].dFluxPSF; + secfilt[Nsec].FkronStk = zpFactor * measure[k].FluxKron; + secfilt[Nsec].dFkronStk = zpFactor * measure[k].dFluxKron; + secfilt[Nsec].FapStk = zpFactor * measure[k].FluxAp; + secfilt[Nsec].dFapStk = zpFactor * measure[k].dFluxAp; + + secfilt[Nsec].MpsfStk = (measure[k].FluxPSF > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FpsfStk) : NAN; + secfilt[Nsec].MkronStk = (measure[k].FluxKron > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FkronStk) : NAN; + secfilt[Nsec].MapStk = (measure[k].FluxAp > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FapStk) : NAN; secfilt[Nsec].stackDetectID = ID; @@ -607,4 +630,153 @@ } +# undef SKIP_THIS_MEAS +# define SKIP_THIS_MEAS(REASON) { \ + measure[k].dbFlags |= ID_MEAS_SKIP_PHOTOM; \ + results->REASON ++; \ + continue; } + +// set mean of forced-warp measurements (selected by photcode range for now): +// somewhat simplified relative to chip-photometry: +// * no grid, no mosaic, no 2MASS, no SYNTH, no Ubercal, no flatcorr +int setMrelAverageForcedWarp (off_t measureOffset, int cat, int pass, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t *found) { + + off_t k; + float Fsys = 0, Mcal= 0; + + // we are measuring means for 3 types of FLUXes: psf, ap, kron. I am using the psf mag + // error for the ap mags, but krons have their own errors. it is an open question if I + // should be doing weighted or unweighted fits (this is a user option) + double *Fpsflist = results->Mpsflist; + double *dpsflist = results->dpsflist; + double *wpsflist = results->wpsflist; + + double *Faplist = results->Maplist; + double *daplist = results->daplist; + double *waplist = results->waplist; + + double *Fkronlist = results->Mkronlist; + double *dkronlist = results->dkronlist; + double *wkronlist = results->wkronlist; + + StatType *psfstats = &results->psfstats; + StatType *apstats = &results->apstats; + StatType *kronstats = &results->kronstats; + + // option for a test print + if (FALSE && (average[0].objID == 0x7146) && (average[0].catID == 0x49d8)) { + fprintf (stderr, "test obj\n"); + print_measure_set_alt (average, secfilt, measure); + } + + int Ns; + for (Ns = 0; Ns < Nphotcodes; Ns++) { + + int thisCode = photcodes[Ns][0].code; + int Nsec = GetPhotcodeNsec(thisCode); + + /* calculate the average mag in this SEC photcode for a single star */ + + /* star/photcodes already calibrated */ + if (found[Nsec]) continue; + + off_t meas = measureOffset; + + int Nap = 0; + int Npsf = 0; + int Nkron = 0; + for (k = 0; k < average[0].Nmeasure; k++, meas++) { + + // only examine gpc1 forced-warp data + if (!isGPC1warp(measure[k].photcode)) continue; + + // skip measurements that do not match the current photcode + PhotCode *code = GetPhotcodebyCode (measure[k].photcode); + if (!code) continue; + if (code->equiv != thisCode) { continue; } + + if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad); + + if (getImageEntry (meas, cat) < 0) { + // measurements without an image are either external reference photometry or + // data for which the associated image has not been loaded (probably because of + // overlaps). Msys + measure.Mcal is our best guess of the true magnitude + Mcal = measure[k].Mcal; // check that this is zero for loaded REF value + } else { + // use getMcal not getMcal_alt? + Mcal = getMcal_alt (meas, cat, NULL, measure[k].Xccd, measure[k].Yccd); + // Mcal = getMcal (meas, cat); + if (isnan(Mcal)) SKIP_THIS_MEAS(Ncal); + } + float Fcal = MagToFlux(-Mcal); + + // in the calculations below, + // ...list gives the error per measurement, wlist gives the weight + // we can modify the error and weight in a few ways: + // 1) MIN_ERROR guarantees a floor + // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode + + // skip some absurd values NAN + // NOTE : I am using PhotCat not PhotSys for now since GPC1 chip-to-chip color terms + // are small (and not measured) + float Fpsf = PhotFluxCat (&measure[k], MAG_CLASS_PSF); + if (isnan(Fsys)) SKIP_THIS_MEAS(Nsys); + // if (Msys < 0.0) SKIP_THIS_MEAS(Nsys); + // if (Msys > 30.0) SKIP_THIS_MEAS(Nsys); + + float dFpsf = PhotFluxCatErr (&measure[k], MAG_CLASS_PSF); + dFpsf = MAX (dFpsf, MIN_ERROR*Fpsf); // MIN_ERROR is a fractional error + Fpsflist[Npsf] = Fpsf * Fcal; + dpsflist[Npsf] = dFpsf * Fcal; + wpsflist[Npsf] = 1.0; + Npsf ++; + + float Fap = PhotFluxCat (&measure[k], MAG_CLASS_APER); + float dFap = PhotFluxCatErr (&measure[k], MAG_CLASS_APER); + if (!isnan(Fap)) { + Faplist[Nap] = Fap * Fcal; + daplist[Nap] = dFap * Fcal; + waplist[Nap] = 1.0; // drop weight lists? + Nap ++; + } + + float Fkron = PhotFluxCat (&measure[k], MAG_CLASS_KRON); + float dFkron = PhotFluxCatErr (&measure[k], MAG_CLASS_KRON); + if (!isnan(Fkron)) { + Fkronlist[Nkron] = Fkron * Fcal; + dkronlist[Nkron] = dFkron * Fcal; + wkronlist[Nkron] = 1.0; + Nkron ++; + } + } + if (Npsf < 1) continue; + + found[Nsec] = TRUE; + + liststats (Fpsflist, dpsflist, wpsflist, Npsf, psfstats); + + secfilt[Nsec].FpsfWrp = psfstats->mean; + secfilt[Nsec].dFpsfWrp = psfstats->error; + secfilt[Nsec].sFpsfWrp = psfstats->sigma; // Mstdev is in millimags (not enough space for more precision) + secfilt[Nsec].NusedWrp = psfstats->Nmeas; + secfilt[Nsec].MpsfWrp = isnan(secfilt[Nsec].FpsfWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FpsfWrp); // 8.9 since flux is in Jy + + // NOTE : use the modified weight for apmags as well as psf mags + liststats (Faplist, daplist, waplist, Nap, apstats); + secfilt[Nsec].FapWrp = Nap > 0 ? apstats->mean : NAN; + secfilt[Nsec].dFapWrp = Nap > 0 ? apstats->error : NAN; + secfilt[Nsec].sFapWrp = Nap > 0 ? apstats->sigma : NAN; + secfilt[Nsec].NusedApWrp = Nap; + secfilt[Nsec].MapWrp = isnan(secfilt[Nsec].FapWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FapWrp); // 8.9 since flux is in Jy + + liststats (Fkronlist, dkronlist, wkronlist, Nkron, kronstats); + secfilt[Nsec].FkronWrp = Nkron > 0 ? kronstats->mean : NAN; + secfilt[Nsec].dFkronWrp = Nkron > 0 ? kronstats->error : NAN; + secfilt[Nsec].sFkronWrp = Nkron > 0 ? kronstats->sigma : NAN; + secfilt[Nsec].NusedKronWrp = Nkron; + secfilt[Nsec].MkronWrp = isnan(secfilt[Nsec].FkronWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FkronWrp); // 8.9 since flux is in Jy + } + return (TRUE); +} + int setGlobalObjStats (Average *average, Measure *measure) { Index: /trunk/Ohana/src/relphot/src/setMrelFinal.c =================================================================== --- /trunk/Ohana/src/relphot/src/setMrelFinal.c (revision 37036) +++ /trunk/Ohana/src/relphot/src/setMrelFinal.c (revision 37037) @@ -85,7 +85,9 @@ /* set catalog[0].found[i] = FALSE */ - ALLOCATE (catalog[0].found, off_t, MAX (1, Nsecfilt*catalog[0].Naverage)); + ALLOCATE (catalog[0].found_t, off_t, MAX (1, Nsecfilt*catalog[0].Naverage)); + ALLOCATE (catalog[0].foundWarp_t, off_t, MAX (1, Nsecfilt*catalog[0].Naverage)); for (i = 0; i < Nsecfilt*catalog[0].Naverage; i++) { - catalog[0].found[i] = FALSE; + catalog[0].found_t[i] = FALSE; + catalog[0].foundWarp_t[i] = FALSE; } @@ -194,5 +196,5 @@ /* star/photcodes already calibrated */ - if (catalog[0].found[Nsecfilt*i+Nsec]) continue; + if (catalog[0].found_t[Nsecfilt*i+Nsec]) continue; m = catalog[0].average[i].measureOffset; @@ -229,5 +231,5 @@ if ((pass < 4) && ImagSelect) { if (Nim > -1) { - mag = PhotInst (&catalog[0].measure[m]); + mag = PhotInst (&catalog[0].measure[m], MAG_CLASS_PSF); if (mag < ImagMin) goto skip; if (mag > ImagMax) goto skip; Index: /trunk/Ohana/src/relphot/src/share_mean_mags.c =================================================================== --- /trunk/Ohana/src/relphot/src/share_mean_mags.c (revision 37036) +++ /trunk/Ohana/src/relphot/src/share_mean_mags.c (revision 37037) @@ -123,5 +123,5 @@ meanmags->M = secfilt->M; meanmags->dM = secfilt->dM; - meanmags->Xm = secfilt->Xm; + meanmags->Mchisq = secfilt->Mchisq; meanmags->Nsec = Nsec; // key to secfilt entry