Index: /trunk/Ohana/src/relphot/Makefile =================================================================== --- /trunk/Ohana/src/relphot/Makefile (revision 42787) +++ /trunk/Ohana/src/relphot/Makefile (revision 42788) @@ -82,6 +82,7 @@ $(SRC)/setExclusions.$(ARCH).o \ $(SRC)/setMrelFinal.$(ARCH).o \ -$(SRC)/setMrelCatalog.$(ARCH).o \ -$(SRC)/BoundaryTreeOps.$(ARCH).o \ +$(SRC)/setMrelCatalog.$(ARCH).o \ +$(SRC)/BoundaryTreeOps.$(ARCH).o \ +$(SRC)/DumpAllMags.$(ARCH).o \ $(SRC)/write_coords.$(ARCH).o Index: /trunk/Ohana/src/relphot/include/relphot.h =================================================================== --- /trunk/Ohana/src/relphot/include/relphot.h (revision 42787) +++ /trunk/Ohana/src/relphot/include/relphot.h (revision 42788) @@ -160,4 +160,7 @@ float dY; // bin = int(Ychip * dY), dY = Ny / NyChip // NxChip, NyChip = 4900,4900 for now + + int **nAlloc; // allocated vector length + double ***dMval; // values used to calculate corrections } GridCorrectionType; @@ -607,4 +610,5 @@ void initGridBins PROTO((void)); GridCorrectionType *getGridCorrByCode PROTO((int code)); +GridCorrectionType *newGridCorrByCode PROTO((int code)); GridCorrectionType *getGridCorrNext PROTO((int *Nlast)); int GridCorrectionSave PROTO((void)); @@ -894,2 +898,4 @@ double get_median_zpt_tgroups (short photcode); void set_median_zpt_tgroups (short photcode, double zpt); + +int DumpAllMags(char *filename, Catalog *catalog, int Ncatalog); Index: /trunk/Ohana/src/relphot/src/DumpAllMags.c =================================================================== --- /trunk/Ohana/src/relphot/src/DumpAllMags.c (revision 42788) +++ /trunk/Ohana/src/relphot/src/DumpAllMags.c (revision 42788) @@ -0,0 +1,231 @@ +# include "relphot.h" + +// STATUS is value expected for success +# define CHECK_STATUS(STATUS,MSG,...) \ + if (!(STATUS)) { \ + fprintf (stderr, MSG, __VA_ARGS__); \ + return FALSE; \ + } + +int DumpAllMags(char *filename, Catalog *catalog, int Ncatalog) { + + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + gfits_init_header (&header); + header.extend = TRUE; + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + + gfits_create_table_header (&theader, "BINTABLE", "DATA"); + + 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", "MCAL", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MGRP", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MMOS", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MFLT", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MGRD", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MREL", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MSYS", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MOFF", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "AIRMASS", "", "unitless", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "MJD", "", "unitless", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "EXPTIME", "", "unitless", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "PHOTCODE", "", "unitless", 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "J", "IX", "", "unitless", 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "J", "IY", "", "unitless", 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "J", "XCCD", "", "unitless", 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "J", "YCCD", "", "unitless", 1.0, FT_BZERO_INT32); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + // count the max number of measurements: + int NptsMax = 0; + for (int nc = 0; nc < Ncatalog; nc++) { + for (int na = 0; na < catalog[nc].Naverage; na++) { + + int nm = catalog[nc].averageT[na].measureOffset; + for (int k = 0; k < catalog[nc].averageT[na].Nmeasure; k++, nm++) { + NptsMax++; + } + } + } + + // create intermediate storage arrays + ALLOCATE_PTR (Rs, double, NptsMax); + ALLOCATE_PTR (Ds, double, NptsMax); + ALLOCATE_PTR (Mcal, float, NptsMax); + ALLOCATE_PTR (Mgrp, float, NptsMax); + ALLOCATE_PTR (Mmos, float, NptsMax); + ALLOCATE_PTR (Mflt, float, NptsMax); + ALLOCATE_PTR (Mgrd, float, NptsMax); + ALLOCATE_PTR (Mrel, float, NptsMax); + ALLOCATE_PTR (Msys, float, NptsMax); + ALLOCATE_PTR (Moff, float, NptsMax); + ALLOCATE_PTR (airmass, float, NptsMax); + ALLOCATE_PTR (mjd, double, NptsMax); + ALLOCATE_PTR (exptime, float, NptsMax); + ALLOCATE_PTR (photcode, int, NptsMax); + ALLOCATE_PTR (iX, int, NptsMax); + ALLOCATE_PTR (iY, int, NptsMax); + ALLOCATE_PTR (Xccd, int, NptsMax); + ALLOCATE_PTR (Yccd, int, NptsMax); + + // XXXX fill in the vectors here + int Nsecfilt = GetPhotcodeNsecfilt (); + + int Npts = 0; + for (int nc = 0; nc < Ncatalog; nc++) { + for (int na = 0; na < catalog[nc].Naverage; na++) { + + int nm = catalog[nc].averageT[na].measureOffset; + for (int k = 0; k < catalog[nc].averageT[na].Nmeasure; k++, nm++) { + + // skip measurements marked by AREA or TIME + if (catalog[nc].measureT[nm].dbFlags & MEAS_BAD) continue; + + Mcal[Npts] = getMcal (nm, nc, MAG_CLASS_PSF); + if (isnan(Mcal[Npts])) continue; + + Mgrp[Npts] = getMgrp (nm, nc, catalog[nc].measureT[nm].airmass, NULL); + if (isnan(Mgrp[Npts])) continue; + + Mmos[Npts] = getMmos (nm, nc); + if (isnan(Mmos[Npts])) continue; + + Mflt[Npts] = getMflat (nm, nc, catalog); + if (isnan(Mflt[Npts])) continue; + + // Mrel* is the average magnitude for this star. For PS1 stacks, we have too much + // PSF variability. We need to calibrate the PSF magnitudes separately from the + // Aperture-like magnitues. (We have an option to use the kron magnitudes or the + // other apertures here). I basically need to do this analysis separately for each + // magnitude type + + // XXX works if we stop loop here + // Npts ++; + // continue; + + Mrel[Npts] = getMrel (catalog, nm, nc, MAG_CLASS_PSF, MAG_SRC_CHP); + if (isnan(Mrel[Npts])) continue; + + Msys[Npts] = PhotSysTiny (&catalog[nc].measureT[nm], &catalog[nc].averageT[na], &catalog[nc].secfilt[na*Nsecfilt], MAG_CLASS_PSF); + if (isnan(Msys[Npts])) continue; + + Mgrd[Npts] = getMgridTiny (&catalog[nc].measureT[nm]); + if (isnan(Mgrd[Npts])) continue; + + Moff[Npts] = Mcal[Npts] + Mgrp[Npts] + Mmos[Npts] + Mflt[Npts] + Mgrd[Npts]; + + photcode[Npts] = catalog[nc].measureT[nm].photcode; + + if (GRID_ZEROPT) { + GridCorrectionType *grid = getGridCorrByCode(photcode[Npts]); + if (!grid) continue; // does not match one of our images, skip + + // edge effects could cause some positions to be slightly out of range + // probably should trap extreme outliers + iX[Npts] = MIN(MAX(0, (int)(catalog[nc].measureT[nm].Xccd * grid->dX)), grid->Nx - 1); + iY[Npts] = MIN(MAX(0, (int)(catalog[nc].measureT[nm].Yccd * grid->dY)), grid->Ny - 1); + } else { + iX[Npts] = 0; + iY[Npts] = 0; + } + + Xccd[Npts] = catalog[nc].measureT[nm].Xccd; + Yccd[Npts] = catalog[nc].measureT[nm].Yccd; + + Rs[Npts] = catalog[nc].measureT[nm].R; + Ds[Npts] = catalog[nc].measureT[nm].D; + + airmass[Npts] = catalog[nc].measureT[nm].airmass; + + exptime[Npts] = 0.0; // catalog[nc].measureT[nm].dt; + mjd[Npts] = 0.0; // ohana_sec_to_mjd (catalog[nc].measureT[nm].tzero); + + Npts ++; + } + } + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "RA", Rs, Npts); + gfits_set_bintable_column (&theader, &ftable, "DEC", Ds, Npts); + gfits_set_bintable_column (&theader, &ftable, "MCAL", Mcal, Npts); + gfits_set_bintable_column (&theader, &ftable, "MGRP", Mgrp, Npts); + gfits_set_bintable_column (&theader, &ftable, "MMOS", Mmos, Npts); + gfits_set_bintable_column (&theader, &ftable, "MFLT", Mflt, Npts); + gfits_set_bintable_column (&theader, &ftable, "MGRD", Mgrd, Npts); + gfits_set_bintable_column (&theader, &ftable, "MREL", Mrel, Npts); + gfits_set_bintable_column (&theader, &ftable, "MSYS", Msys, Npts); + gfits_set_bintable_column (&theader, &ftable, "MOFF", Moff, Npts); + gfits_set_bintable_column (&theader, &ftable, "AIRMASS", airmass, Npts); + gfits_set_bintable_column (&theader, &ftable, "MJD", mjd, Npts); + gfits_set_bintable_column (&theader, &ftable, "EXPTIME", exptime, Npts); + gfits_set_bintable_column (&theader, &ftable, "PHOTCODE", photcode, Npts); + gfits_set_bintable_column (&theader, &ftable, "IX", iX, Npts); + gfits_set_bintable_column (&theader, &ftable, "IY", iY, Npts); + gfits_set_bintable_column (&theader, &ftable, "XCCD", Xccd, Npts); + gfits_set_bintable_column (&theader, &ftable, "YCCD", Yccd, Npts); + + free (Rs); + free (Ds); + free (Mcal); + free (Mgrp); + free (Mmos); + free (Mflt); + free (Mgrd); + free (Mrel); + free (Msys); + free (Moff); + free (airmass); + free (mjd); + free (exptime); + free (photcode); + free (iX); + free (iY); + free (Xccd); + free (Yccd); + + FILE *f = fopen (filename, "w"); + if (!f) { + fprintf (stderr, "ERROR: cannot open meanmag file for output %s\n", filename); + return FALSE; + } + + int status; + status = gfits_fwrite_header (f, &header); + CHECK_STATUS (status, "ERROR: cannot write header for meanmags %s\n", filename); + + status = gfits_fwrite_matrix (f, &matrix); + CHECK_STATUS (status, "ERROR: cannot write matrix for meanmags %s\n", filename); + + status = gfits_fwrite_Theader (f, &theader); + CHECK_STATUS (status, "ERROR: cannot write table header for meanmags %s\n", filename); + + status = gfits_fwrite_table (f, &ftable); + CHECK_STATUS (status, "ERROR: cannot write table data for meanmags %s\n", filename); + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + int fd = fileno (f); + + status = fflush (f); + CHECK_STATUS (!status, "ERROR: cannot flush file mags %s\n", filename); + + status = fsync (fd); + CHECK_STATUS (!status, "ERROR: cannot flush file mags %s\n", filename); + + status = fclose (f); + CHECK_STATUS (!status, "ERROR: problem closing mags file %s\n", filename); + + return TRUE; +} + Index: /trunk/Ohana/src/relphot/src/GridIO.c =================================================================== --- /trunk/Ohana/src/relphot/src/GridIO.c (revision 42787) +++ /trunk/Ohana/src/relphot/src/GridIO.c (revision 42788) @@ -151,5 +151,5 @@ if (!gfits_scan (&header, "PHOTCODE", "%d", 1, &photcode)) myAbort ("failed to get PHOTCODE"); - GridCorrectionType *GridCorr = getGridCorrByCode (photcode); + GridCorrectionType *GridCorr = newGridCorrByCode (photcode); GridCorr->photcode = photcode; Index: /trunk/Ohana/src/relphot/src/GridOps.c =================================================================== --- /trunk/Ohana/src/relphot/src/GridOps.c (revision 42787) +++ /trunk/Ohana/src/relphot/src/GridOps.c (revision 42788) @@ -30,6 +30,6 @@ # define NX_CHIP_DEFAULT 1000 # define NY_CHIP_DEFAULT 1000 -# define NX_BIN_DEFAULT 2 -# define NY_BIN_DEFAULT 2 +# define NX_BIN_DEFAULT 4 +# define NY_BIN_DEFAULT 4 void initGridBins (void) { @@ -76,4 +76,11 @@ if (GridCorr[code]->dMgrid) FREE (GridCorr[code]->dMgrid[ix]); if (GridCorr[code]->nMgrid) FREE (GridCorr[code]->nMgrid[ix]); + if (GridCorr[code]->nAlloc) FREE (GridCorr[code]->nAlloc[ix]); + + // dMval is a 2D array of vectors + for (int iy = 0; iy < GridCorr[code]->Ny; iy++) { + FREE (GridCorr[code]->dMval[ix][iy]); + } + FREE (GridCorr[code]->dMval[ix]); } @@ -81,4 +88,6 @@ FREE (GridCorr[code]->dMgrid); FREE (GridCorr[code]->nMgrid); + FREE (GridCorr[code]->nAlloc); + FREE (GridCorr[code]-> dMval); FREE(GridCorr[code]); @@ -109,5 +118,6 @@ } -GridCorrectionType *getGridCorrByCode (int code) { +// generate a new GridCorr structure in the array +GridCorrectionType *newGridCorrByCode (int code) { if (GridCorr == NULL) return NULL; @@ -116,5 +126,18 @@ if (code < 0) return NULL; + myAssert (GridCorr[code] == NULL, "oops, should not already be allocated!"); ALLOCATE(GridCorr[code], GridCorrectionType, 1) + GridCorrectionType *result = GridCorr[code]; + return result; +} + +// return the GridCode structure corresponding to 'code' +GridCorrectionType *getGridCorrByCode (int code) { + + if (GridCorr == NULL) return NULL; + + if (code >= NGridCorr) return NULL; + if (code < 0) return NULL; + GridCorrectionType *result = GridCorr[code]; return result; @@ -163,11 +186,18 @@ // we are normally accessing this array randomly, so there is no advantage to // doing Mgrid[y][x] vs Mgrid[x][y] + ALLOCATE (GridCorr[code]-> dMval, double **, NX_BIN); ALLOCATE (GridCorr[code]-> Mgrid, float *, NX_BIN); ALLOCATE (GridCorr[code]->dMgrid, float *, NX_BIN); ALLOCATE (GridCorr[code]->nMgrid, int *, NX_BIN); + ALLOCATE (GridCorr[code]->nAlloc, int *, NX_BIN); for (int ix = 0; ix < NX_BIN; ix++) { ALLOCATE (GridCorr[code]-> Mgrid[ix], float, NY_BIN); ALLOCATE (GridCorr[code]->dMgrid[ix], float, NY_BIN); ALLOCATE (GridCorr[code]->nMgrid[ix], int, NY_BIN); + ALLOCATE (GridCorr[code]->nAlloc[ix], int, NY_BIN); + + // when we first allocate this array, set the elements to 0 (NULL) + // these arrays will be allocated or reallocated in resetMgrid + ALLOCATE_ZERO (GridCorr[code]-> dMval[ix], double *, NX_BIN); } } @@ -190,4 +220,12 @@ GridCorr[code]->dMgrid[ix][iy] = 0.0; GridCorr[code]->nMgrid[ix][iy] = 0; + GridCorr[code]->nAlloc[ix][iy] = 100; + + // allocate or reset vector length to default + if (GridCorr[code]->dMval[ix][iy]) { + REALLOCATE (GridCorr[code]->dMval[ix][iy], double, GridCorr[code]->nAlloc[ix][iy]); + } else { + ALLOCATE (GridCorr[code]->dMval[ix][iy], double, GridCorr[code]->nAlloc[ix][iy]); + } } } @@ -263,10 +301,49 @@ if (fabs(dM) > 0.5) continue; - grid-> Mgrid[ix][iy] += dM; - grid->dMgrid[ix][iy] += dM*dM; + int ival = grid->nMgrid[ix][iy]; + grid-> dMval[ix][iy][ival] = dM; grid->nMgrid[ix][iy] ++; + CHECK_REALLOCATE (grid-> dMval[ix][iy], double, grid->nAlloc[ix][iy], grid->nMgrid[ix][iy], 100); + + /* XXX old code to calculate simple means: + grid-> Mgrid[ix][iy] += dM; + grid->dMgrid[ix][iy] += dM*dM; + grid->nMgrid[ix][iy] ++; + */ } } } + + StatType stats; + liststats_setmode (&stats, "MEDIAN"); + + // now calculate Mgrid, dMgrid, nMgrid from Sum, Sum, Npt + for (int code = 0; code < NGridCorr; code++) { + if (!GridCorr[code]) continue; + + for (int ix = 0; ix < GridCorr[code]->Nx; ix++) { + for (int iy = 0; iy < GridCorr[code]->Ny; iy++) { + + // cells without sufficient coverage stay at 0.0 + if (GridCorr[code]->nMgrid[ix][iy] < 5) { + GridCorr[code]-> Mgrid[ix][iy] = 0.0; + GridCorr[code]->dMgrid[ix][iy] = 0.0; + continue; + } + + liststats_init (&stats); + + liststats (GridCorr[code]->dMval[ix][iy], NULL, NULL, GridCorr[code]->nMgrid[ix][iy], &stats); + double altSigma = (stats.Upper80 - stats.Lower20) / 1.6; // 20% to 80% encompasses 60% of the values, corresponds to the range (-0.85 sigma : +0.85 sigma) + + GridCorr[code]-> Mgrid[ix][iy] = stats.median; + GridCorr[code]->dMgrid[ix][iy] = altSigma; // sample stdev + } + } + } + return; +} + +# if (0) // now calculate Mgrid, dMgrid, nMgrid from Sum, Sum, Npt @@ -300,7 +377,5 @@ // float GridAve = GridSum / GridCnt; // fprintf (stderr, "grid average: %f\n", GridAve); - } - return; -} +# endif float getMgrid (Measure *measure) { Index: /trunk/Ohana/src/relphot/src/relphot_images.c =================================================================== --- /trunk/Ohana/src/relphot/src/relphot_images.c (revision 42787) +++ /trunk/Ohana/src/relphot/src/relphot_images.c (revision 42788) @@ -147,4 +147,6 @@ MagResidSave ("mag.resid.fits", catalog); + // DumpAllMags ("all.mags.fits", catalog, Ncatalog); + save_images_updates (&db);