Index: /branches/eam_branches/ipp-20150625/Ohana/src/addstar/Makefile =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/addstar/Makefile (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/addstar/Makefile (revision 38983) @@ -39,5 +39,5 @@ # programs in 'SERVER' use the client-server concept and are out of date -INSTALL = addstar sedstar loadstarpar loadstarpar_client loadICRF loadICRF_client load2mass loadwise loadtycho loadsupercos skycells mkcmf dumpskycells findskycell +INSTALL = addstar sedstar loadgalphot loadstarpar loadstarpar_client loadICRF loadICRF_client skycells mkcmf dumpskycells findskycell load2mass loadwise loadtycho loadsupercos SERVER = addstar_client addstarc addstard addstart @@ -259,4 +259,5 @@ $(SRC)/loadgalphot_join.$(ARCH).o \ $(SRC)/loadgalphot_table.$(ARCH).o \ +$(SRC)/resort_catalog.$(ARCH).o \ $(SRC)/strhash.$(ARCH).o \ $(SRC)/sortIDs.$(ARCH).o \ @@ -284,4 +285,5 @@ $(SRC)/loadstarpar_table.$(ARCH).o \ $(SRC)/loadstarpar_extras.$(ARCH).o \ +$(SRC)/resort_catalog.$(ARCH).o \ $(SRC)/psps_ids.$(ARCH).o @@ -302,5 +304,8 @@ $(SRC)/loadstarpar_table.$(ARCH).o \ $(SRC)/loadstarpar_extras.$(ARCH).o \ -$(SRC)/psps_ids.$(ARCH).o +$(SRC)/resort_catalog.$(ARCH).o \ +$(SRC)/psps_ids.$(ARCH).o + +# drop resort_catalog? LOAD-ICRF = \ Index: /branches/eam_branches/ipp-20150625/Ohana/src/addstar/include/addstar.h =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/addstar/include/addstar.h (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/addstar/include/addstar.h (revision 38983) @@ -333,4 +333,9 @@ Catalog *addstar_catalog_init (int Nstars); +void resort_catalog_measure (Catalog *catalog); +void resort_catalog_lensing (Catalog *catalog); +void resort_catalog_starpar (Catalog *catalog); +void resort_catalog_galphot (Catalog *catalog); + /** there is an inconsistency to be resolved: fixed structures (like Image) Index: /branches/eam_branches/ipp-20150625/Ohana/src/addstar/include/loadgalphot.h =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/addstar/include/loadgalphot.h (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/addstar/include/loadgalphot.h (revision 38983) @@ -36,9 +36,9 @@ } GalPhotIDset; -AddstarClientOptions args_loadgalphot (int argc, char **argv, AddstarClientOptions options); +AddstarClientOptions args_loadgalphot (int *argc, char **argv, AddstarClientOptions options); int loadgalphot_table (SkyList *skylistInput, HostTable *hosts, char *filename, AddstarClientOptions *options); -GalPhot_Stars *loadgalphot_readstars (char *filename, int *nstars); +GalPhot_Stars *loadgalphot_readstars (char *filename, int *nstars, AddstarClientOptions *options); int loadgalphot_sortStars (GalPhot_Stars *stars, int Nstars); Index: /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/ReadXradFITS.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/ReadXradFITS.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/ReadXradFITS.c (revision 38983) @@ -14,6 +14,11 @@ /* load the table data */ if (!gfits_fread_ftable_data (f, &table, FALSE)) { - fprintf (stderr, "ERROR: can't read table header\n"); - exit (1); + fprintf (stderr, "WARNING: no xrad data, skipping\n"); + return TRUE; + } + + if (theader[0].Naxis[1] == 0) { + fprintf (stderr, "WARNING: no xrad data, skipping\n"); + return TRUE; } @@ -65,9 +70,10 @@ for (i = 0; i < catalog->Nmeasure; i++) { if (catalog->measure[i].detID < RadID[Nap]) { - Nap ++; continue; + // this is a psf measurement which does not have a radial aperture } if (catalog->measure[i].detID > RadID[Nap]) { myAbort("radial apertures out of order? seems like a bug\n"); + // this would be a radial aperture which does not have a radial aperture } @@ -95,5 +101,5 @@ Nap ++; } - catalog->Nlensing = Nap; + myAssert (Nap == Nrow, "did we go too far???"); gfits_free_table (&table); Index: /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/args_loadgalphot.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/args_loadgalphot.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/args_loadgalphot.c (revision 38983) @@ -4,11 +4,11 @@ static void help (void); -AddstarClientOptions args_loadgalphot (int argc, char **argv, AddstarClientOptions options) { +AddstarClientOptions args_loadgalphot (int *argc, char **argv, AddstarClientOptions options) { int N; /* check for help request */ - if (get_argument (argc, argv, "-help") || - get_argument (argc, argv, "-h")) { + if (get_argument (*argc, argv, "-help") || + get_argument (*argc, argv, "-h")) { help (); } @@ -27,20 +27,20 @@ HOST_ID = 0; PARALLEL = FALSE; - if ((N = get_argument (argc, argv, "-parallel"))) { + if ((N = get_argument (*argc, argv, "-parallel"))) { PARALLEL = TRUE; - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); } // this is a test mode : rather than launching the remote jobs and waiting for completion, // relphot will simply list the remote command and wait for the user to signal completion PARALLEL_MANUAL = FALSE; - if ((N = get_argument (argc, argv, "-parallel-manual"))) { + if ((N = get_argument (*argc, argv, "-parallel-manual"))) { PARALLEL = TRUE; // -parallel-manual implies -parallel PARALLEL_MANUAL = TRUE; - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); } // this is a test mode : rather than launching the relphot_client jobs remotely, they are // run in serial via 'system' PARALLEL_SERIAL = FALSE; - if ((N = get_argument (argc, argv, "-parallel-serial"))) { + if ((N = get_argument (*argc, argv, "-parallel-serial"))) { if (PARALLEL_MANUAL) { fprintf (stderr, "ERROR: cannot mix -parallel-manual and -parallel-serial\n"); @@ -49,35 +49,46 @@ PARALLEL = TRUE; // -parallel-serial implies -parallel PARALLEL_SERIAL = TRUE; - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); } /* only add to existing regions */ options.existing_regions = FALSE; - if ((N = get_argument (argc, argv, "-existing-regions"))) { + if ((N = get_argument (*argc, argv, "-existing-regions"))) { options.existing_regions = TRUE; - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); } /* only add to existing objects */ options.only_match = FALSE; - if ((N = get_argument (argc, argv, "-only-match"))) { + if ((N = get_argument (*argc, argv, "-only-match"))) { options.only_match = TRUE; - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); } /* replace measurement, don't duplicate (ref/cat only) */ options.replace = FALSE; - if ((N = get_argument (argc, argv, "-replace"))) { + if ((N = get_argument (*argc, argv, "-replace"))) { options.replace = TRUE; - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); + } + /* override any header PHOTCODE values */ + options.photcode = 0; + if ((N = get_argument (*argc, argv, "-p"))) { + remove_argument (N, argc, argv); + options.photcode = GetPhotcodeCodebyName (argv[N]); + remove_argument (N, argc, argv); + } + if ((N = get_argument (*argc, argv, "-photcode"))) { + remove_argument (N, argc, argv); + options.photcode = GetPhotcodeCodebyName (argv[N]); + remove_argument (N, argc, argv); } /* extra error messages */ VERBOSE = FALSE; - if ((N = get_argument (argc, argv, "-v"))) { + if ((N = get_argument (*argc, argv, "-v"))) { VERBOSE = TRUE; - remove_argument (N, &argc, argv); + remove_argument (N, argc, argv); } /* other addstar options which cannot be used in loadgalphot */ - options.photcode = 0; options.timeref = 0; options.mosaic = FALSE; @@ -95,5 +106,5 @@ DUMP = NULL; - if (argc < 2) { + if (*argc < 2) { fprintf (stderr, "USAGE: loadgalphot [options] (fitsfile) [..more files]\n"); exit (2); Index: /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/find_matches_galphot.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/find_matches_galphot.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/find_matches_galphot.c (revision 38983) @@ -199,6 +199,4 @@ REALLOCATE (catalog[0].galphot, GalPhot, Ngalphot); - catalog[0].sorted = FALSE; - /* check if the catalog has changed? if no change, no need to write */ catalog[0].objID = objID; // new max value, save on catalog close @@ -206,4 +204,10 @@ catalog[0].Ngalphot = Ngalphot; catalog[0].Nsecfilt_mem = Nave*Nsecfilt; + + // we need to for + catalog[0].sorted = FALSE; + resort_catalog_galphot (catalog); + catalog[0].sorted = TRUE; + if (VERBOSE) fprintf (stderr, "Nstars, Nave, Ngalphot: %d "OFF_T_FMT" "OFF_T_FMT" ("OFF_T_FMT" matches)\n", Nstars, Nave, Ngalphot, Nmatch); Index: /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/find_matches_starpar.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/find_matches_starpar.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/find_matches_starpar.c (revision 38983) @@ -200,5 +200,6 @@ REALLOCATE (catalog[0].starpar, StarPar, Nstarpar); - catalog[0].sorted = FALSE; + resort_catalog_starpar (catalog); + catalog[0].sorted = TRUE; /* check if the catalog has changed? if no change, no need to write */ Index: /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/loadgalphot.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/loadgalphot.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/loadgalphot.c (revision 38983) @@ -15,5 +15,5 @@ SetSignals (); options = ConfigInit (&argc, argv); - options = args_loadgalphot (argc, argv, options); + options = args_loadgalphot (&argc, argv, options); // load the full sky description table (dvodb must exist) Index: /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/loadgalphot_fit2d.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/loadgalphot_fit2d.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/loadgalphot_fit2d.c (revision 38983) @@ -86,5 +86,7 @@ int N, i, j, nx, ny, K, k, n, Npt, status; - double mean, sigma, maxsigma; + double mean = 0.0; + double sigma = 0.0; + double maxsigma = 0.0; float *x, *y, *z, *zf; Index: /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/loadgalphot_readstars.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/loadgalphot_readstars.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/loadgalphot_readstars.c (revision 38983) @@ -17,5 +17,5 @@ static float ZeroPt = 0; -GalPhot_Stars *loadgalphot_readstars (char *filename, int *nstars) { +GalPhot_Stars *loadgalphot_readstars (char *filename, int *nstars, AddstarClientOptions *options) { // read in the full FITS files ('cause I don't have a partial read option) Index: /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/loadgalphot_table.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/loadgalphot_table.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/loadgalphot_table.c (revision 38983) @@ -6,5 +6,5 @@ int i, Nstars; - GalPhot_Stars *stars = loadgalphot_readstars (filename, &Nstars); + GalPhot_Stars *stars = loadgalphot_readstars (filename, &Nstars, options); // sort the stars by RA Index: /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/mkcmf.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/mkcmf.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/mkcmf.c (revision 38983) @@ -25,4 +25,9 @@ void writeStars_PS1_DEV_0 (FTable *ftable, double *X, double *Y, double *M, int Nstars); +int WriteXSRCtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars, float zeroPt, float exptime); +int WriteXFITtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars); +int WriteXRADtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars, int Nrad); +int WriteDETFtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars); + int ADDNOISE = TRUE; float BAD_PSFQF_FRAC = 0.0; @@ -54,4 +59,21 @@ } + int isStack = FALSE; + int isForcedWarp = FALSE; + int WriteXSRC = FALSE; + int WriteXFIT = FALSE; + int WriteXRAD = FALSE; + int WriteDETF = FALSE; + if ((N = get_argument (argc, argv, "-stack"))) { + remove_argument (N, &argc, argv); + WriteXSRC = WriteXFIT = WriteXRAD = WriteDETF = TRUE; + isStack = TRUE; + } + if ((N = get_argument (argc, argv, "-forcedwarp"))) { + remove_argument (N, &argc, argv); + WriteXSRC = WriteXFIT = WriteXRAD = WriteDETF = TRUE; + isForcedWarp = TRUE; + } + static char *photcode = "SIMTEST.r.Chip"; if ((N = get_argument (argc, argv, "-photcode"))) { @@ -76,4 +98,18 @@ } + static char *tess_id = "RINGS.V3"; + if ((N = get_argument (argc, argv, "-tess_id"))) { + remove_argument (N, &argc, argv); + tess_id = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + + static char *skycell = "skycell.1133.081"; + if ((N = get_argument (argc, argv, "-skycell"))) { + remove_argument (N, &argc, argv); + skycell = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + static char *time = "00:00:00"; if ((N = get_argument (argc, argv, "-time"))) { @@ -83,9 +119,11 @@ } - char extname[80], exthead[80]; + char extname[80], exthead[80], extroot[80]; + strcpy (extroot, "Chip"); strcpy (extname, "Chip.psf"); strcpy (exthead, "Chip.hdr"); - if ((N = get_argument (argc, argv, "-extname"))) { - remove_argument (N, &argc, argv); + if ((N = get_argument (argc, argv, "-extroot"))) { + remove_argument (N, &argc, argv); + snprintf (extroot, 80, "%s", argv[N]); snprintf (extname, 80, "%s.psf", argv[N]); snprintf (exthead, 80, "%s.hdr", argv[N]); @@ -206,6 +244,11 @@ coords.crpix2 = CRPIX2; - coords.cdelt1 = 0.25/3600.0; - coords.cdelt2 = 0.25/3600.0; + if (isStack || isForcedWarp) { + coords.cdelt1 = 0.25/3600.0; + coords.cdelt2 = 0.25/3600.0; + } else { + coords.cdelt1 = 0.257/3600.0; + coords.cdelt2 = 0.257/3600.0; + } // load stars and generate complete output fields @@ -280,4 +323,6 @@ gfits_print (&header, "IMNAXIS1", "%d", 1, (int)(Xmax)); gfits_print (&header, "IMNAXIS2", "%d", 1, (int)(Ymax)); + gfits_print (&header, "NAXIS1", "%d", 1, (int)(Xmax)); + gfits_print (&header, "NAXIS2", "%d", 1, (int)(Ymax)); gfits_modify (&header, "NSTARS", "%d", 1, Nstars); @@ -289,5 +334,5 @@ gfits_modify (&header, "UTC-OBS", "%s", 1, time); } - gfits_modify (&header, "ZERO_PT", "%lf", 1, ZERO_POINT); + gfits_modify (&header, "ZERO_PT", "%lf", 1, ZERO_POINT); // this is now not used gfits_modify (&header, "EXPTIME", "%lf", 1, exptime); gfits_modify (&header, "AIRMASS", "%lf", 1, airmass); @@ -301,4 +346,53 @@ PutCoords (&coords, &header); gfits_modify (&header, "EXTNAME", "%s", 1, exthead); + + // various other fields needed to ipptopsps + gfits_modify (&header, "FWHM_MAJ", "%f", 1, FX*2.8); + gfits_modify (&header, "FWHM_MIN", "%f", 1, FY*2.8); + gfits_modify (&header, "MSKY_MN", "%f", 1, SKY); + gfits_modify (&header, "MSKY_SIG", "%f", 1, DSKY); + gfits_modify (&header, "ZPT_ERR", "%f", 1, 0.03); + gfits_modify (&header, "ANGLE", "%f", 1, 0.1); + gfits_modify (&header, "IQ_FW1", "%f", 1, FX*2.8); + gfits_modify (&header, "IQ_FW2", "%f", 1, FX*2.8); + gfits_modify (&header, "IQ_M2C", "%f", 1, FX*2.8); + gfits_modify (&header, "IQ_M2S", "%f", 1, FX*2.8); + gfits_modify (&header, "IQ_M3", "%f", 1, FX*2.8); + gfits_modify (&header, "IQ_M4", "%f", 1, FX*2.8); + gfits_modify (&header, "APMIFIT", "%f", 1, 0.1); + gfits_modify (&header, "DAPMIFIT", "%f", 1, 0.03); + gfits_modify (&header, "DETECTOR", "%s", 1, "CCID58-1-02a2"); + + gfits_modify (&header, "HIERARCH DETEFF.MAGREF", "%f", 1, -5.0); + + gfits_modify (&header, "PSFMODEL", "%s", 1, "PS_MODEL_PS1_V1"); + gfits_modify (&header, "AST_CDX", "%f", 1, 0.05); + gfits_modify (&header, "AST_CDY", "%f", 1, 0.04); + + float zeroPt = 25.0 + 0.3*(drand48() - 0.5); + + if (!isStack && !isForcedWarp) { + gfits_modify (&header, "HIERARCH DETREND.MASK", "%s", 1, "detref615.XY33.fits"); + gfits_modify (&header, "HIERARCH DETREND.DARK", "%s", 1, "GPC1.DARKTEST.norm.856.0.XY33.fits"); + gfits_modify (&header, "HIERARCH DETREND.FLAT", "%s", 1, "GPC1.FLATTEST.300.XY33.co.fits"); + gfits_modify (&header, "HIERARCH DETREND.NOISEMAP", "%s", 1, "GPC1.noisemap.norm.965.0.XY33.fits"); + gfits_modify (&header, "HIERARCH DETREND.NONLIN", "%s", 1, "linearity_data.XY33.fits"); + gfits_modify (&header, "HIERARCH DETREND.VIDEODARK", "%s", 1, "GPC1.VIDEODARK.979.0.XY33.fits"); + gfits_modify (&header, "ZPT_OBS", "%f", 1, zeroPt); + } + + if (isStack) { + gfits_modify (&header, "HIERARCH FPA.ZP", "%f", 1, zeroPt); + gfits_modify (&header, "TESS_ID", "%s", 1, tess_id); + gfits_modify (&header, "SKYCELL", "%s", 1, skycell); + gfits_modify (&header, "ZPT_OBS", "%f", 1, zeroPt); + } + if (isForcedWarp) { + gfits_modify (&header, "HIERARCH FPA.FILTERID", "%s", 1, "r.00000"); // does this affect anything? I don't actually think so... + gfits_modify (&header, "HIERARCH FPA.ZP", "%f", 1, zeroPt); + gfits_modify (&header, "TESS_ID", "%s", 1, tess_id); + gfits_modify (&header, "SKYCELL", "%s", 1, skycell); + gfits_modify (&header, "PHOT_V", "%s", 1, "38100"); + } int i; @@ -380,4 +474,13 @@ gfits_fwrite_Theader (fits, ftable.header); gfits_fwrite_table (fits, &ftable); + + if (WriteXSRC) WriteXSRCtable (fits, extroot, X, Y, M, Flag, Nstars, zeroPt, exptime); + if (WriteXFIT) WriteXFITtable (fits, extroot, X, Y, M, Flag, Nstars); + + if (WriteXRAD && isStack) WriteXRADtable (fits, extroot, X, Y, M, Flag, Nstars, 3); + if (WriteXRAD && isForcedWarp) WriteXRADtable (fits, extroot, X, Y, M, Flag, Nstars, 1); + + if (WriteDETF) WriteDETFtable (fits, extroot, X, Y, M, Flag, Nstars); + fclose (fits); @@ -475,13 +578,15 @@ fSN = 1.0 / sqrt(flux); + if (ADDNOISE) { + X[i] += FX * fSN * rnd_gauss(0.0, 1.0); + Y[i] += FY * fSN * rnd_gauss(0.0, 1.0); + M[i] += fSN*rnd_gauss(0.0, 1.0); + flux = pow (10.0, -0.4*M[i]); + fSN = 1.0 / sqrt(flux); + } + stars[i].X = X[i]; stars[i].Y = Y[i]; stars[i].M = M[i]; - - if (ADDNOISE) { - stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0); - stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0); - stars[i].M += fSN*rnd_gauss(0.0, 1.0); - } stars[i].dX = FX * fSN; @@ -519,13 +624,15 @@ fSN = 1.0 / sqrt(flux); + if (ADDNOISE) { + X[i] += FX * fSN * rnd_gauss(0.0, 1.0); + Y[i] += FY * fSN * rnd_gauss(0.0, 1.0); + M[i] += fSN*rnd_gauss(0.0, 1.0); + flux = pow (10.0, -0.4*M[i]); + fSN = 1.0 / sqrt(flux); + } + stars[i].X = X[i]; stars[i].Y = Y[i]; stars[i].M = M[i]; - - if (ADDNOISE) { - stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0); - stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0); - stars[i].M += fSN*rnd_gauss(0.0, 1.0); - } stars[i].dX = FX * fSN; @@ -566,13 +673,15 @@ fSN = 1.0 / sqrt(flux); + if (ADDNOISE) { + X[i] += FX * fSN * rnd_gauss(0.0, 1.0); + Y[i] += FY * fSN * rnd_gauss(0.0, 1.0); + M[i] += fSN*rnd_gauss(0.0, 1.0); + flux = pow (10.0, -0.4*M[i]); + fSN = 1.0 / sqrt(flux); + } + stars[i].X = X[i]; stars[i].Y = Y[i]; stars[i].M = M[i]; - - if (ADDNOISE) { - stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0); - stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0); - stars[i].M += fSN*rnd_gauss(0.0, 1.0); - } stars[i].dX = FX * fSN; @@ -615,13 +724,15 @@ fSN = 1.0 / sqrt(flux); + if (ADDNOISE) { + X[i] += FX * fSN * rnd_gauss(0.0, 1.0); + Y[i] += FY * fSN * rnd_gauss(0.0, 1.0); + M[i] += fSN*rnd_gauss(0.0, 1.0); + flux = pow (10.0, -0.4*M[i]); + fSN = 1.0 / sqrt(flux); + } + stars[i].X = X[i]; stars[i].Y = Y[i]; stars[i].M = M[i]; - - if (ADDNOISE) { - stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0); - stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0); - stars[i].M += fSN*rnd_gauss(0.0, 1.0); - } // randomly give poor PSFQF values @@ -670,14 +781,16 @@ fSN = 1.0 / sqrt(flux); + if (ADDNOISE) { + X[i] += FX * fSN * rnd_gauss(0.0, 1.0); + Y[i] += FY * fSN * rnd_gauss(0.0, 1.0); + M[i] += fSN*rnd_gauss(0.0, 1.0); + flux = pow (10.0, -0.4*M[i]); + fSN = 1.0 / sqrt(flux); + } + stars[i].X = X[i]; stars[i].Y = Y[i]; stars[i].M = M[i]; stars[i].Map = M[i] - 0.05; - - if (ADDNOISE) { - stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0); - stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0); - stars[i].M += fSN*rnd_gauss(0.0, 1.0); - } // randomly give poor PSFQF values @@ -729,4 +842,12 @@ fSN = 1.0 / sqrt(flux); + if (ADDNOISE) { + X[i] += FX * fSN * rnd_gauss(0.0, 1.0); + Y[i] += FY * fSN * rnd_gauss(0.0, 1.0); + M[i] += fSN*rnd_gauss(0.0, 1.0); + flux = pow (10.0, -0.4*M[i]); + fSN = 1.0 / sqrt(flux); + } + stars[i].X = X[i]; stars[i].Y = Y[i]; @@ -740,10 +861,4 @@ stars[i].apFlux = pow(10.0, -0.4*stars[i].Map); stars[i].apFluxErr = stars[i].apFlux * fSN; - - if (ADDNOISE) { - stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0); - stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0); - stars[i].M += fSN*rnd_gauss(0.0, 1.0); - } // randomly give poor PSFQF values @@ -796,4 +911,12 @@ fSN = 1.0 / sqrt(flux); + if (ADDNOISE) { + X[i] += FX * fSN * rnd_gauss(0.0, 1.0); + Y[i] += FY * fSN * rnd_gauss(0.0, 1.0); + M[i] += fSN*rnd_gauss(0.0, 1.0); + flux = pow (10.0, -0.4*M[i]); + fSN = 1.0 / sqrt(flux); + } + stars[i].detID = i; stars[i].X = X[i]; @@ -803,5 +926,5 @@ stars[i].posangle = 10.0; - stars[i].pltscale = 0.25; + stars[i].pltscale = coords.cdelt1; stars[i].M = M[i]; @@ -875,14 +998,4 @@ stars[i].nFrames = 1; - - if (ADDNOISE) { - stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0); - stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0); - float Moff = fSN*rnd_gauss(0.0, 1.0); - stars[i].M += Moff; - stars[i].Map += Moff; - stars[i].MapRaw += Moff; - stars[i].Mcalib += Moff; - } } @@ -904,4 +1017,12 @@ flux = pow (10.0, -0.4*M[i]); fSN = 1.0 / sqrt(flux); + + if (ADDNOISE) { + X[i] += FX * fSN * rnd_gauss(0.0, 1.0); + Y[i] += FY * fSN * rnd_gauss(0.0, 1.0); + M[i] += fSN*rnd_gauss(0.0, 1.0); + flux = pow (10.0, -0.4*M[i]); + fSN = 1.0 / sqrt(flux); + } stars[i].detID = i; @@ -984,14 +1105,4 @@ stars[i].nFrames = 1; - - if (ADDNOISE) { - stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0); - stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0); - float Moff = fSN*rnd_gauss(0.0, 1.0); - stars[i].M += Moff; - stars[i].Map += Moff; - stars[i].MapRaw += Moff; - stars[i].Mcalib += Moff; - } } @@ -1005,4 +1116,12 @@ CMF_PS1_V5_R2_Lensing *stars; float flux, fSN; + + // We want some fraction of the stars to have zero psfqfperfect so that we can test the + // ForcedWarpMasked table construction, but we don't want things to be random so we + // change every Nstars / modulo to zero + int modulo = 0; + if (BAD_PSFQF_FRAC > 0) { + modulo = (int) (1.0 / BAD_PSFQF_FRAC); + } // XXX add gaussian-distributed noise based on counts @@ -1013,4 +1132,12 @@ flux = pow (10.0, -0.4*M[i]); fSN = 1.0 / sqrt(flux); + + if (ADDNOISE) { + X[i] += FX * fSN * rnd_gauss(0.0, 1.0); + Y[i] += FY * fSN * rnd_gauss(0.0, 1.0); + M[i] += fSN*rnd_gauss(0.0, 1.0); + flux = pow (10.0, -0.4*M[i]); + fSN = 1.0 / sqrt(flux); + } stars[i].detID = i; @@ -1057,8 +1184,7 @@ stars[i].fwhmMin = FY*2.8; - // randomly give poor PSFQF values - if ((BAD_PSFQF_FRAC > 0.0) && (drand48() < BAD_PSFQF_FRAC)) { - stars[i].psfQF = 0.25; - stars[i].psfQFperf = 0.24; + if (modulo && ((i+1) % modulo == 0)) { + stars[i].psfQF = 0.0; + stars[i].psfQFperf = 0.0; } else { stars[i].psfQF = PSFQUAL; @@ -1103,4 +1229,8 @@ stars[i].E2_sh_psf = FX; + stars[i].srcChipNum = 23; + stars[i].srcChipX = X[i] + 10.0; + stars[i].srcChipY = Y[i] + 10.0; + stars[i].kronFlux = flux * 1.25; stars[i].kronFluxErr = fSN * flux * 1.25; @@ -1117,14 +1247,4 @@ stars[i].nFrames = 1; - - if (ADDNOISE) { - stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0); - stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0); - float Moff = fSN*rnd_gauss(0.0, 1.0); - stars[i].M += Moff; - stars[i].Map += Moff; - stars[i].MapRaw += Moff; - stars[i].Mcalib += Moff; - } } @@ -1133,2 +1253,427 @@ } +int WriteXSRCtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars, float zeroPt, float exptime) { + + int i, j; + Header header; + FTable ftable; + + char extdata[80]; + snprintf (extdata, 80, "%s.xsrc", extroot); + + /* Example Code to create a bintable */ + { + gfits_create_table_header (&header, "BINTABLE", extdata); + + gfits_define_bintable_column (&header, "1J", "IPP_IDET" , "no comment", NULL, 1.0, FT_BZERO_INT32); // unsigned + gfits_define_bintable_column (&header, "1E", "X_EXT" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "1E", "Y_EXT" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "1E", "X_EXT_SIG" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "1E", "Y_EXT_SIG" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "1E", "F25_ARATIO" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "1E", "F25_THETA" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "1E", "PETRO_MAG" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "1E", "PETRO_MAG_ERR" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "1E", "PETRO_RADIUS" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "1E", "PETRO_RADIUS_ERR" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "1E", "PETRO_RADIUS_50" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "1E", "PETRO_RADIUS_50_ERR", "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "1E", "PETRO_RADIUS_90" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "1E", "PETRO_RADIUS_90_ERR", "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "1E", "PETRO_FILL" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "1E", "HALF_LIGHT_RADIUS" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "1E", "G_RT" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "1E", "G_RA" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "1E", "G_S2" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "1E", "G_A" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "1E", "G_BUMPY" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "9E", "PROF_SB" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "9E", "PROF_FLUX" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "9E", "PROF_FILL" , "no comment", NULL, 1.0, 0.0); + + // generate the output array that carries the data + gfits_create_table (&header, &ftable); + + // create intermediate storage arrays + int *IPP_IDET ; ALLOCATE (IPP_IDET , int , Nstars); + float *X_EXT ; ALLOCATE (X_EXT , float , Nstars); + float *Y_EXT ; ALLOCATE (Y_EXT , float , Nstars); + float *X_EXT_SIG ; ALLOCATE (X_EXT_SIG , float , Nstars); + float *Y_EXT_SIG ; ALLOCATE (Y_EXT_SIG , float , Nstars); + float *F25_ARATIO ; ALLOCATE (F25_ARATIO , float , Nstars); + float *F25_THETA ; ALLOCATE (F25_THETA , float , Nstars); + float *PETRO_MAG ; ALLOCATE (PETRO_MAG , float , Nstars); + float *PETRO_MAG_ERR ; ALLOCATE (PETRO_MAG_ERR , float , Nstars); + float *PETRO_RADIUS ; ALLOCATE (PETRO_RADIUS , float , Nstars); + float *PETRO_RADIUS_ERR ; ALLOCATE (PETRO_RADIUS_ERR , float , Nstars); + float *PETRO_RADIUS_50 ; ALLOCATE (PETRO_RADIUS_50 , float , Nstars); + float *PETRO_RADIUS_50_ERR ; ALLOCATE (PETRO_RADIUS_50_ERR , float , Nstars); + float *PETRO_RADIUS_90 ; ALLOCATE (PETRO_RADIUS_90 , float , Nstars); + float *PETRO_RADIUS_90_ERR ; ALLOCATE (PETRO_RADIUS_90_ERR , float , Nstars); + float *PETRO_FILL ; ALLOCATE (PETRO_FILL , float , Nstars); + float *HALF_LIGHT_RADIUS ; ALLOCATE (HALF_LIGHT_RADIUS , float , Nstars); + float *G_RT ; ALLOCATE (G_RT , float , Nstars); + float *G_RA ; ALLOCATE (G_RA , float , Nstars); + float *G_S2 ; ALLOCATE (G_S2 , float , Nstars); + float *G_A ; ALLOCATE (G_A , float , Nstars); + float *G_BUMPY ; ALLOCATE (G_BUMPY , float , Nstars); + float *PROF_SB ; ALLOCATE (PROF_SB , float , 9*Nstars); + float *PROF_FLUX ; ALLOCATE (PROF_FLUX , float , 9*Nstars); + float *PROF_FILL ; ALLOCATE (PROF_FILL , float , 9*Nstars); + + float magtime = 2.5*log10(exptime); + + // assign the storage arrays + for (i = 0; i < Nstars; i++) { + float flux = pow (10.0, -0.4*M[i]); + float fSN = 1.0 / sqrt(flux); + + IPP_IDET [i] = i; + X_EXT [i] = X[i]; + Y_EXT [i] = Y[i]; + X_EXT_SIG [i] = FX * fSN; + Y_EXT_SIG [i] = FY * fSN; + F25_ARATIO [i] = 0.9; + F25_THETA [i] = 10.0; // seems to be in degrees in the cmfs + PETRO_MAG [i] = M[i] + zeroPt + magtime; + PETRO_MAG_ERR [i] = 1.0 / fSN; // note that the real cmfs have PETRO_MAG_ERR inverted + PETRO_RADIUS [i] = 8.0; + if (i == 2) { + // This case is here to exercise a case for ipp2psps. + PETRO_RADIUS_ERR [i] = NAN; + } + else { + PETRO_RADIUS_ERR [i] = 0.1; + } + PETRO_RADIUS_50 [i] = 4.0; + PETRO_RADIUS_50_ERR [i] = 0.1; + PETRO_RADIUS_90 [i] = 12.0; + PETRO_RADIUS_90_ERR [i] = 0.1; + PETRO_FILL [i] = 0.95; + HALF_LIGHT_RADIUS [i] = 6.0; + G_RT [i] = 1.0; + G_RA [i] = 1.0; + G_S2 [i] = 1.0; + G_A [i] = 1.0; + G_BUMPY [i] = 1.0; + for (j = 0; j < 9; j++) { + PROF_SB [j + 9*i] = flux / 10.0; + PROF_FLUX [j + 9*i] = flux; + PROF_FILL [j + 9*i] = 0.9; + } + } + + // set the table data + gfits_set_bintable_column (&header, &ftable, "IPP_IDET" , IPP_IDET , Nstars); + gfits_set_bintable_column (&header, &ftable, "X_EXT" , X_EXT , Nstars); + gfits_set_bintable_column (&header, &ftable, "Y_EXT" , Y_EXT , Nstars); + gfits_set_bintable_column (&header, &ftable, "X_EXT_SIG" , X_EXT_SIG , Nstars); + gfits_set_bintable_column (&header, &ftable, "Y_EXT_SIG" , Y_EXT_SIG , Nstars); + gfits_set_bintable_column (&header, &ftable, "F25_ARATIO" , F25_ARATIO , Nstars); + gfits_set_bintable_column (&header, &ftable, "F25_THETA" , F25_THETA , Nstars); + gfits_set_bintable_column (&header, &ftable, "PETRO_MAG" , PETRO_MAG , Nstars); + gfits_set_bintable_column (&header, &ftable, "PETRO_MAG_ERR" , PETRO_MAG_ERR , Nstars); + gfits_set_bintable_column (&header, &ftable, "PETRO_RADIUS" , PETRO_RADIUS , Nstars); + gfits_set_bintable_column (&header, &ftable, "PETRO_RADIUS_ERR" , PETRO_RADIUS_ERR , Nstars); + gfits_set_bintable_column (&header, &ftable, "PETRO_RADIUS_50" , PETRO_RADIUS_50 , Nstars); + gfits_set_bintable_column (&header, &ftable, "PETRO_RADIUS_50_ERR", PETRO_RADIUS_50_ERR , Nstars); + gfits_set_bintable_column (&header, &ftable, "PETRO_RADIUS_90" , PETRO_RADIUS_90 , Nstars); + gfits_set_bintable_column (&header, &ftable, "PETRO_RADIUS_90_ERR", PETRO_RADIUS_90_ERR , Nstars); + gfits_set_bintable_column (&header, &ftable, "PETRO_FILL" , PETRO_FILL , Nstars); + gfits_set_bintable_column (&header, &ftable, "HALF_LIGHT_RADIUS" , HALF_LIGHT_RADIUS , Nstars); + gfits_set_bintable_column (&header, &ftable, "G_RT" , G_RT , Nstars); + gfits_set_bintable_column (&header, &ftable, "G_RA" , G_RA , Nstars); + gfits_set_bintable_column (&header, &ftable, "G_S2" , G_S2 , Nstars); + gfits_set_bintable_column (&header, &ftable, "G_A" , G_A , Nstars); + gfits_set_bintable_column (&header, &ftable, "G_BUMPY" , G_BUMPY , Nstars); + gfits_set_bintable_column (&header, &ftable, "PROF_SB" , PROF_SB , Nstars); + gfits_set_bintable_column (&header, &ftable, "PROF_FLUX" , PROF_FLUX , Nstars); + gfits_set_bintable_column (&header, &ftable, "PROF_FILL" , PROF_FILL , Nstars); + } + + gfits_fwrite_Theader (fits, &header); + gfits_fwrite_table (fits, &ftable); + gfits_free_header (&header); + gfits_free_table (&ftable); + + return TRUE; +} + +int WriteXFITtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars) { + + int i; + Header header; + FTable ftable; + + char extdata[80]; + snprintf (extdata, 80, "%s.xfit", extroot); + + /* Example Code to create a bintable */ + { + gfits_create_table_header (&header, "BINTABLE", extdata); + + // define the table layout + gfits_define_bintable_column (&header, "J", "IPP_IDET" , "no comment", NULL, 1.0, FT_BZERO_INT32); // unsigned + gfits_define_bintable_column (&header, "E", "X_EXT" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "E", "Y_EXT" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "E", "X_EXT_SIG" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "E", "Y_EXT_SIG" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "E", "SKY_EXT" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "E", "RA_EXT" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "E", "DEC_EXT" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "E", "EXT_INST_MAG" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "E", "EXT_INST_MAG_SIG" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "E", "EXT_CHISQ" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "J", "EXT_NDOF" , "no comment", NULL, 1.0, FT_BZERO_INT32); // unsigned + gfits_define_bintable_column (&header, "J", "EXT_MODEL_TYPE" , "no comment", NULL, 1.0, FT_BZERO_INT32); // unsigned + gfits_define_bintable_column (&header, "I", "EXT_FLAGS" , "no comment", NULL, 1.0, FT_BZERO_INT16); // unsigned + gfits_define_bintable_column (&header, "E", "PSF_INST_MAG" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "E", "AP_MAG" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "E", "KRON_MAG" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "J", "NPARAMS" , "no comment", NULL, 1.0, FT_BZERO_INT32); // unsigned + gfits_define_bintable_column (&header, "15A", "MODEL_TYPE" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "E", "EXT_WIDTH_MAJ" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "E", "EXT_WIDTH_MIN" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "E", "EXT_THETA" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "E", "EXT_WIDTH_MAJ_ERR", "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "E", "EXT_WIDTH_MIN_ERR", "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "E", "EXT_THETA_ERR" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "E", "EXT_PAR_07" , "no comment", NULL, 1.0, 0.0); + + // generate the output array that carries the data + gfits_create_table (&header, &ftable); + + // create intermediate storage arrays + int *IPP_IDET ; ALLOCATE (IPP_IDET , int , Nstars); + float *X_EXT ; ALLOCATE (X_EXT , float , Nstars); + float *Y_EXT ; ALLOCATE (Y_EXT , float , Nstars); + float *X_EXT_SIG ; ALLOCATE (X_EXT_SIG , float , Nstars); + float *Y_EXT_SIG ; ALLOCATE (Y_EXT_SIG , float , Nstars); + float *SKY_EXT ; ALLOCATE (SKY_EXT , float , Nstars); + float *RA_EXT ; ALLOCATE (RA_EXT , float , Nstars); + float *DEC_EXT ; ALLOCATE (DEC_EXT , float , Nstars); + float *EXT_INST_MAG ; ALLOCATE (EXT_INST_MAG , float , Nstars); + float *EXT_INST_MAG_SIG ; ALLOCATE (EXT_INST_MAG_SIG , float , Nstars); + float *EXT_CHISQ ; ALLOCATE (EXT_CHISQ , float , Nstars); + int *EXT_NDOF ; ALLOCATE (EXT_NDOF , int , Nstars); + int *EXT_MODEL_TYPE ; ALLOCATE (EXT_MODEL_TYPE , int , Nstars); + short *EXT_FLAGS ; ALLOCATE (EXT_FLAGS , short , Nstars); + float *PSF_INST_MAG ; ALLOCATE (PSF_INST_MAG , float , Nstars); + float *AP_MAG ; ALLOCATE (AP_MAG , float , Nstars); + float *KRON_MAG ; ALLOCATE (KRON_MAG , float , Nstars); + int *NPARAMS ; ALLOCATE (NPARAMS , int , Nstars); + char *MODEL_TYPE ; ALLOCATE (MODEL_TYPE , char , Nstars*15); + float *EXT_WIDTH_MAJ ; ALLOCATE (EXT_WIDTH_MAJ , float , Nstars); + float *EXT_WIDTH_MIN ; ALLOCATE (EXT_WIDTH_MIN , float , Nstars); + float *EXT_THETA ; ALLOCATE (EXT_THETA , float , Nstars); + float *EXT_WIDTH_MAJ_ERR ; ALLOCATE (EXT_WIDTH_MAJ_ERR, float , Nstars); + float *EXT_WIDTH_MIN_ERR ; ALLOCATE (EXT_WIDTH_MIN_ERR, float , Nstars); + float *EXT_THETA_ERR ; ALLOCATE (EXT_THETA_ERR , float , Nstars); + float *EXT_PAR_07 ; ALLOCATE (EXT_PAR_07 , float , Nstars); + + // assign the storage arrays + for (i = 0; i < Nstars; i++) { + float flux = pow (10.0, -0.4*M[i]); + float fSN = 1.0 / sqrt(flux); + + double ra, dec; + XY_to_RD (&ra, &dec, X[i], Y[i], &coords); + + IPP_IDET [i] = i; + X_EXT [i] = X[i]; + Y_EXT [i] = Y[i]; + X_EXT_SIG [i] = FX * fSN; + Y_EXT_SIG [i] = FY * fSN; + SKY_EXT [i] = 0.1; + RA_EXT [i] = ra; + DEC_EXT [i] = dec; + EXT_INST_MAG [i] = M[i]; + EXT_INST_MAG_SIG [i] = fSN; + EXT_CHISQ [i] = 1.0; + EXT_NDOF [i] = 10; + EXT_MODEL_TYPE [i] = -1; + EXT_FLAGS [i] = 0x0010; + PSF_INST_MAG [i] = M[i]; + AP_MAG [i] = M[i]; + KRON_MAG [i] = M[i]; + NPARAMS [i] = 7; + EXT_WIDTH_MAJ [i] = FX*2.8; + EXT_WIDTH_MIN [i] = FY*2.8; + EXT_THETA [i] = 0.57; + EXT_WIDTH_MAJ_ERR[i] = FX*0.28; + EXT_WIDTH_MIN_ERR[i] = FY*0.28; + EXT_THETA_ERR [i] = 0.0057; + EXT_PAR_07 [i] = NAN; + + strcpy (&MODEL_TYPE[15*i], "PS_MODEL_DEV"); + } + + // set the table data + gfits_set_bintable_column (&header, &ftable, "IPP_IDET" , IPP_IDET , Nstars); + gfits_set_bintable_column (&header, &ftable, "X_EXT" , X_EXT , Nstars); + gfits_set_bintable_column (&header, &ftable, "Y_EXT" , Y_EXT , Nstars); + gfits_set_bintable_column (&header, &ftable, "X_EXT_SIG" , X_EXT_SIG , Nstars); + gfits_set_bintable_column (&header, &ftable, "Y_EXT_SIG" , Y_EXT_SIG , Nstars); + gfits_set_bintable_column (&header, &ftable, "SKY_EXT" , SKY_EXT , Nstars); + gfits_set_bintable_column (&header, &ftable, "RA_EXT" , RA_EXT , Nstars); + gfits_set_bintable_column (&header, &ftable, "DEC_EXT" , DEC_EXT , Nstars); + gfits_set_bintable_column (&header, &ftable, "EXT_INST_MAG" , EXT_INST_MAG , Nstars); + gfits_set_bintable_column (&header, &ftable, "EXT_INST_MAG_SIG" , EXT_INST_MAG_SIG , Nstars); + gfits_set_bintable_column (&header, &ftable, "EXT_CHISQ" , EXT_CHISQ , Nstars); + gfits_set_bintable_column (&header, &ftable, "EXT_NDOF" , EXT_NDOF , Nstars); + gfits_set_bintable_column (&header, &ftable, "EXT_MODEL_TYPE" , EXT_MODEL_TYPE , Nstars); + gfits_set_bintable_column (&header, &ftable, "EXT_FLAGS" , EXT_FLAGS , Nstars); + gfits_set_bintable_column (&header, &ftable, "PSF_INST_MAG" , PSF_INST_MAG , Nstars); + gfits_set_bintable_column (&header, &ftable, "AP_MAG" , AP_MAG , Nstars); + gfits_set_bintable_column (&header, &ftable, "KRON_MAG" , KRON_MAG , Nstars); + gfits_set_bintable_column (&header, &ftable, "NPARAMS" , NPARAMS , Nstars); + gfits_set_bintable_column (&header, &ftable, "MODEL_TYPE" , MODEL_TYPE , Nstars); + gfits_set_bintable_column (&header, &ftable, "EXT_WIDTH_MAJ" , EXT_WIDTH_MAJ , Nstars); + gfits_set_bintable_column (&header, &ftable, "EXT_WIDTH_MIN" , EXT_WIDTH_MIN , Nstars); + gfits_set_bintable_column (&header, &ftable, "EXT_THETA" , EXT_THETA , Nstars); + gfits_set_bintable_column (&header, &ftable, "EXT_WIDTH_MAJ_ERR", EXT_WIDTH_MAJ_ERR, Nstars); + gfits_set_bintable_column (&header, &ftable, "EXT_WIDTH_MIN_ERR", EXT_WIDTH_MIN_ERR, Nstars); + gfits_set_bintable_column (&header, &ftable, "EXT_THETA_ERR" , EXT_THETA_ERR , Nstars); + gfits_set_bintable_column (&header, &ftable, "EXT_PAR_07" , EXT_PAR_07 , Nstars); + } + /**/ + + gfits_fwrite_Theader (fits, &header); + gfits_fwrite_table (fits, &ftable); + gfits_free_header (&header); + gfits_free_table (&ftable); + + return TRUE; +} + +static float psfFWHM[3] = {4.8, 6.0, 8.0}; + +int WriteXRADtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars, int Nrad) { + + int i, j, k; + Header header; + FTable ftable; + + char extdata[80]; + snprintf (extdata, 80, "%s.xrad", extroot); + + /* Example Code to create a bintable */ + { + gfits_create_table_header (&header, "BINTABLE", extdata); + + // define the table layout + gfits_define_bintable_column (&header, "J", "IPP_IDET" , "no comment", NULL, 1.0, FT_BZERO_INT32); // unsigned + gfits_define_bintable_column (&header, "E", "X_APER" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "E", "Y_APER" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "E", "PSF_FWHM" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "9E", "APER_FLUX" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "9E", "APER_FLUX_ERR" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "9E", "APER_FLUX_STDEV" , "no comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "9E", "APER_FILL" , "no comment", NULL, 1.0, 0.0); + + // generate the output array that carries the data + gfits_create_table (&header, &ftable); + + // create intermediate storage arrays + int *IPP_IDET ; ALLOCATE (IPP_IDET , int , Nrad*Nstars); + float *X_APER ; ALLOCATE (X_APER , float , Nrad*Nstars); + float *Y_APER ; ALLOCATE (Y_APER , float , Nrad*Nstars); + float *PSF_FWHM ; ALLOCATE (PSF_FWHM , float , Nrad*Nstars); + float *APER_FLUX ; ALLOCATE (APER_FLUX , float , Nrad*Nstars*9); + float *APER_FLUX_ERR ; ALLOCATE (APER_FLUX_ERR , float , Nrad*Nstars*9); + float *APER_FLUX_STDEV ; ALLOCATE (APER_FLUX_STDEV , float , Nrad*Nstars*9); + float *APER_FILL ; ALLOCATE (APER_FILL , float , Nrad*Nstars*9); + + // assign the storage arrays + for (i = 0; i < Nstars; i++) { + float flux = pow (10.0, -0.4*M[i]); + // float fSN = 1.0 / sqrt(flux); + + for (k = 0; k < Nrad; k++) { + int N = Nrad*i + k; + IPP_IDET [N] = i; + X_APER [N] = X[i]; + Y_APER [N] = Y[i]; + PSF_FWHM [N] = psfFWHM[k]; + for (j = 0; j < 9; j ++) { + float Flux = flux*(1.0 + 0.02*j - 0.2); + APER_FLUX [j + 9*N] = Flux; + APER_FLUX_ERR [j + 9*N] = sqrt(Flux); + APER_FLUX_STDEV[j + 9*N] = sqrt(Flux*1.1); + APER_FILL [j + 9*N] = 0.9; + } + } + } + + // set the table data + gfits_set_bintable_column (&header, &ftable, "IPP_IDET" , IPP_IDET , Nrad*Nstars); + gfits_set_bintable_column (&header, &ftable, "X_APER" , X_APER , Nrad*Nstars); + gfits_set_bintable_column (&header, &ftable, "Y_APER" , Y_APER , Nrad*Nstars); + gfits_set_bintable_column (&header, &ftable, "PSF_FWHM" , PSF_FWHM , Nrad*Nstars); + gfits_set_bintable_column (&header, &ftable, "APER_FLUX" , APER_FLUX , Nrad*Nstars); + gfits_set_bintable_column (&header, &ftable, "APER_FLUX_ERR" , APER_FLUX_ERR , Nrad*Nstars); + gfits_set_bintable_column (&header, &ftable, "APER_FLUX_STDEV" , APER_FLUX_STDEV , Nrad*Nstars); + gfits_set_bintable_column (&header, &ftable, "APER_FILL" , APER_FILL , Nrad*Nstars); + } + /**/ + + gfits_fwrite_Theader (fits, &header); + gfits_fwrite_table (fits, &ftable); + gfits_free_header (&header); + gfits_free_table (&ftable); + + return TRUE; +} + +int WriteDETFtable (FILE *fits, char *extroot, double *X, double *Y, double *M, unsigned int *Flag, int Nstars) { + return TRUE; +} + +int WriteDummyTable (FILE *fits) { + + Header header; + FTable ftable; + + /* Example Code to create a bintable + { + gfits_create_table_header (&header, "BINTABLE", extroot); + + // define the table layout + gfits_define_bintable_column (&header, "D", "NAME", "comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "E", "NAME", "comment", NULL, 1.0, 0.0); + gfits_define_bintable_column (&header, "J", "NAME", "comment", NULL, 1.0, FT_BZERO_INT32); // unsigned + gfits_define_bintable_column (&header, "I", "NAME", "comment", NULL, 1.0, FT_BZERO_INT16); // unsigned + + // generate the output array that carries the data + gfits_create_table (&header, &ftable); + + // create intermediate storage arrays + double *R; ALLOCATE (R, double, Nvalue); + float *M; ALLOCATE (M, float , Nvalue); + int *I; ALLOCATE (I, int , Nvalue); + short *P; ALLOCATE (P, short , Nvalue); + + // assign the storage arrays + for (i = 0; Nvalue; i++) { + R[i] = value; + M[i] = value; + I[i] = value; + P[i] = value; + } + + // set the table data + gfits_set_bintable_column (&header, &ftable, "NAME", R, catalog->Nmeasure); + gfits_set_bintable_column (&header, &ftable, "NAME", M, catalog->Nmeasure); + gfits_set_bintable_column (&header, &ftable, "NAME", I, catalog->Nmeasure); + gfits_set_bintable_column (&header, &ftable, "NAME", P, catalog->Nmeasure); + } + */ + + gfits_fwrite_Theader (fits, &header); + gfits_fwrite_table (fits, &ftable); + gfits_free_header (&header); + gfits_free_table (&ftable); + + return TRUE; +} + Index: /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/resort_catalog.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/resort_catalog.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/resort_catalog.c (revision 38983) @@ -2,34 +2,4 @@ void SortAveMatch (off_t *MEAS, off_t *AVE, off_t N); -void resort_catalog_measure (Catalog *catalog); -void resort_catalog_lensing (Catalog *catalog); -void resort_catalog_starpar (Catalog *catalog); - -void resort_catalog_old (Catalog *catalog) { - - off_t *next_meas, *next_lens; - off_t Naves, Nmeas, Nlens; - - if (catalog[0].sorted == TRUE) return; - - INITTIME; - - /* internal counters */ - Nmeas = catalog[0].Nmeasure; - Nlens = catalog[0].Nlensing; - Naves = catalog[0].Naverage; - - /* set up pointers for linked list of measure, missing */ - next_meas = build_measure_links (catalog[0].average, Naves, catalog[0].measure, Nmeas); - next_lens = build_lensing_links (catalog[0].average, Naves, catalog[0].lensing, Nlens); - - catalog[0].sorted = TRUE; - catalog[0].measure = sort_measure (catalog[0].average, Naves, catalog[0].measure, Nmeas, next_meas); - catalog[0].lensing = sort_lensing (catalog[0].average, Naves, catalog[0].lensing, Nlens, next_lens); - - MARKTIME (" match time %9.4f sec for %7lld measures, %6lld average\n", dtime, (long long) Nmeas, (long long) Naves); - - return; -} void resort_catalog (Catalog *catalog) { @@ -40,4 +10,6 @@ resort_catalog_lensing (catalog); resort_catalog_starpar (catalog); + resort_catalog_galphot (catalog); + catalog[0].sorted = TRUE; } @@ -83,4 +55,6 @@ // earlier formats did not carry the objID or catID, so they are not available (we could assign on load, but we don't) myAssert(average[averageSeq[i]].catID == measure[measureSeq[i]].catID, "object / detection mismatch"); + myAssert (average[averageSeq[i]].objID == measure[measureSeq[i]].objID, "object ID mismatch?"); +# if (0) // myAssert(average[averageSeq[i]].objID == measure[measureSeq[i]].objID, "object / detection mismatch"); @@ -107,4 +81,5 @@ } } +# endif # if (0) // for reasons I do not understand, the mini dvodbs generated on stsci1X had a handful of detections with an inconsistency between averef and objID. @@ -193,6 +168,4 @@ // MARKTIME(" match time %9.4f sec for %7lld measures, %6lld average\n", dtime, (long long) Nmeasure, (long long) Naverage); - catalog[0].sorted = TRUE; - FREE (measureSeq); FREE (averageSeq); @@ -246,22 +219,5 @@ } - // check that averageSeq is now in order - // for (i = 1; i < Nlensing; i++) { - // if (averageSeq[i] < averageSeq[i-1]) { - // fprintf (stderr, "%d ", (int) i); - // } - // } - // fprintf (stderr, "\n"); - SortAveMatch(lensingSeq, averageSeq, Nlensing); - // MARKTIME("sort : %f sec\n", dtime); - - // check that averageSeq is now in order - // for (i = 1; i < Nlensing; i++) { - // if (averageSeq[i] < averageSeq[i-1]) { - // fprintf (stderr, "%d ", (int) i); - // } - // } - // fprintf (stderr, "\n"); // copy the lensing entries in the sorted order @@ -271,5 +227,4 @@ lensingTMP[i] = lensing[j]; } - // MARKTIME("assign lensing : %f sec\n", dtime); // update the values of average.lensingOffset and average.Nlensing @@ -290,7 +245,5 @@ N++; } - // N++; average[currentAve].Nlensing = N; - // MARKTIME("update Nlensing : %f sec\n", dtime); int NlensingTotal = 0; @@ -318,8 +271,4 @@ fprintf (stderr, "ERROR: catalog %s has an invalid Nlensing\n", catalog[0].filename); } - - // MARKTIME(" match time %9.4f sec for %7lld lensings, %6lld average\n", dtime, (long long) Nlensing, (long long) Naverage); - - catalog[0].sorted = TRUE; FREE (lensingSeq); @@ -374,22 +323,5 @@ } - // check that averageSeq is now in order - // for (i = 1; i < Nstarpar; i++) { - // if (averageSeq[i] < averageSeq[i-1]) { - // fprintf (stderr, "%d ", (int) i); - // } - // } - // fprintf (stderr, "\n"); - SortAveMatch(starparSeq, averageSeq, Nstarpar); - // MARKTIME("sort : %f sec\n", dtime); - - // check that averageSeq is now in order - // for (i = 1; i < Nstarpar; i++) { - // if (averageSeq[i] < averageSeq[i-1]) { - // fprintf (stderr, "%d ", (int) i); - // } - // } - // fprintf (stderr, "\n"); // copy the starpar entries in the sorted order @@ -399,5 +331,4 @@ starparTMP[i] = starpar[j]; } - // MARKTIME("assign starpar : %f sec\n", dtime); // update the values of average.starparOffset and average.Nstarpar @@ -418,7 +349,5 @@ N++; } - // N++; average[currentAve].Nstarpar = N; - // MARKTIME("update Nstarpar : %f sec\n", dtime); int NstarparTotal = 0; @@ -447,9 +376,109 @@ } - // MARKTIME(" match time %9.4f sec for %7lld starpars, %6lld average\n", dtime, (long long) Nstarpar, (long long) Naverage); - - catalog[0].sorted = TRUE; - FREE (starparSeq); + FREE (averageSeq); + + return; +} + +void resort_catalog_galphot (Catalog *catalog) { + + off_t Naverage, Ngalphot; + GalPhot *galphot; + Average *average; + off_t i, j, N, currentAve; + + off_t *galphotSeq = NULL; + off_t *averageSeq = NULL; + GalPhot *galphotTMP = NULL; + + // struct timeval start, stop; + // gettimeofday (&start, NULL); + + /* internal counters */ + Ngalphot = catalog[0].Ngalphot; + Naverage = catalog[0].Naverage; + + if (!Ngalphot) return; + + galphot = catalog[0].galphot; + average = catalog[0].average; + + // we have a table of average objects and an unsorted table of measurements. each measurement + // has a reference to the average object sequence (as well as an ID) + // galphot[i].averef -> average[averef] + // galphot[i].objID = average[averef].objID + // galphot[i].catID = average[averef].catID + + // we want a sorted galphot array with all averef entries in sequence + + ALLOCATE (galphotSeq, off_t, Ngalphot); + ALLOCATE (averageSeq, off_t, Ngalphot); + + for (i = 0; i < Ngalphot; i++) { + galphotSeq[i] = i; + averageSeq[i] = galphot[i].averef; + + if (catalog[0].catformat >= DVO_FORMAT_PS1_V1) { + // earlier formats did not carry the objID or catID, so they are not available (we could assign on load, but we don't) + myAssert(average[averageSeq[i]].catID == galphot[galphotSeq[i]].catID, "object / detection mismatch"); + myAssert(average[averageSeq[i]].objID == galphot[galphotSeq[i]].objID, "object / detection mismatch"); + } + } + + SortAveMatch(galphotSeq, averageSeq, Ngalphot); + + // copy the galphot entries in the sorted order + ALLOCATE (galphotTMP, GalPhot, Ngalphot); + for (i = 0; i < Ngalphot; i++) { + j = galphotSeq[i]; + galphotTMP[i] = galphot[j]; + } + + // update the values of average.galphotOffset and average.Ngalphot + FREE(galphot); + catalog[0].galphot = galphotTMP; + + N = 0; + currentAve = averageSeq[0]; + average[currentAve].galphotOffset = 0; + for (i = 0; i < Ngalphot; i++) { + if (averageSeq[i] != currentAve) { + // we have hit the next entry in the list + average[currentAve].Ngalphot = N; + N = 0; + currentAve = averageSeq[i]; + average[currentAve].galphotOffset = i; + } + N++; + } + average[currentAve].Ngalphot = N; + + int NgalphotTotal = 0; + int galphotOffsetOK = TRUE; + for (i = 0; i < Naverage; i++) { + NgalphotTotal += catalog[0].average[i].Ngalphot; + if (VERBOSE && !(NgalphotTotal <= catalog[0].Ngalphot)) { + fprintf (stderr, "too few galphot: %d %d %d\n", (int) i, NgalphotTotal, (int) catalog[0].Ngalphot); + } + galphotOffsetOK &= (catalog[0].average[i].galphotOffset < catalog[0].Ngalphot); + if (VERBOSE && !(catalog[0].average[i].galphotOffset < catalog[0].Ngalphot)) { + fprintf (stderr, "offset too large: %d %d %d\n", (int) i, catalog[0].average[i].Ngalphot, (int) catalog[0].Ngalphot); + } + galphotOffsetOK &= (catalog[0].average[i].galphotOffset + catalog[0].average[i].Ngalphot <= catalog[0].Ngalphot); + if (VERBOSE && !(catalog[0].average[i].galphotOffset + catalog[0].average[i].Ngalphot <= catalog[0].Ngalphot)) { + fprintf (stderr, "orrset + Ngalphot too large: %d + %d > %d %d\n", (int) i, catalog[0].average[i].galphotOffset, catalog[0].average[i].Ngalphot, (int) catalog[0].Ngalphot); + } + } + + if (!galphotOffsetOK) { + fprintf (stderr, "ERROR: catalog %s has an invalid galphotOffset\n", catalog[0].filename); + } + + if (NgalphotTotal != catalog[0].Ngalphot) { + fprintf (stderr, "ERROR: catalog %s has an invalid Ngalphot\n", catalog[0].filename); + } + + FREE (galphotSeq); FREE (averageSeq); Index: /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/resort_threaded.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/resort_threaded.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/resort_threaded.c (revision 38983) @@ -72,10 +72,5 @@ // this is an overloaded value to mean 'force sort' if (threadData->forcesort) catalog.sorted = FALSE; - - if (OLD_RESORT) { - resort_catalog_old (&catalog); - } else { - resort_catalog (&catalog); - } + resort_catalog (&catalog); /* report total updated values */ Index: /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/resort_unthreaded.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/resort_unthreaded.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/addstar/src/resort_unthreaded.c (revision 38983) @@ -42,12 +42,7 @@ } - // this is an overloaded value to mean 'force sort' + // to force the sort, we need to set sorted 'false' if (ForceSort) catalog.sorted = FALSE; - - if (OLD_RESORT) { - resort_catalog_old (&catalog); - } else { - resort_catalog (&catalog); - } + resort_catalog (&catalog); // report total updated values Index: /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/Makefile =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/Makefile (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/Makefile (revision 38983) @@ -32,4 +32,6 @@ $(SRC)/insert_FWobjects_dvopsps.$(ARCH).o \ $(SRC)/insert_FWobjects_dvopsps_catalog.$(ARCH).o \ +$(SRC)/insert_FGshape_dvopsps.$(ARCH).o \ +$(SRC)/insert_FGshape_dvopsps_catalog.$(ARCH).o \ $(SRC)/insert_diffobj_dvopsps.$(ARCH).o \ $(SRC)/insert_diffobj_dvopsps_catalog.$(ARCH).o \ @@ -52,4 +54,6 @@ $(SRC)/insert_FWobjects_dvopsps.$(ARCH).o \ $(SRC)/insert_FWobjects_dvopsps_catalog.$(ARCH).o \ +$(SRC)/insert_FGshape_dvopsps.$(ARCH).o \ +$(SRC)/insert_FGshape_dvopsps_catalog.$(ARCH).o \ $(SRC)/insert_diffobj_dvopsps.$(ARCH).o \ $(SRC)/insert_diffobj_dvopsps_catalog.$(ARCH).o \ Index: /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/include/dvopsps.h =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/include/dvopsps.h (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/include/dvopsps.h (revision 38983) @@ -20,4 +20,5 @@ float decErr; float zp; + float zpFactor; float telluricExt; float airmass; @@ -115,4 +116,13 @@ int insert_FWobjects_mysql_commit PROTO((IOBuffer *ave_buffer, IOBuffer *sec_buffer, IOBuffer *cpy_buffer, MYSQL *mysql)); +int insert_FGshape_dvopsps PROTO((void)); +int insert_FGshape_dvopsps_parallel PROTO((SkyList *sky)); +int insert_FGshape_dvopsps_catalog PROTO((Catalog *catalog, char *basename, MYSQL *mysql)); + +int insert_FGshape_mysql_create_tables PROTO((char *basename, MYSQL *mysql)); +int insert_FGshape_mysql_init PROTO((IOBuffer *ave_buffer, IOBuffer *cpy_buffer, char *basename)); +int insert_FGshape_mysql_value PROTO((IOBuffer *ave_buffer, IOBuffer *cpy_buffer, Average *average, GalPhot *galphot, int Ngalphot)); +int insert_FGshape_mysql_commit PROTO((IOBuffer *ave_buffer, IOBuffer *cpy_buffer, MYSQL *mysql)); + int insert_skytable PROTO((void)); Index: /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/src/dvopsps.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/src/dvopsps.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/src/dvopsps.c (revision 38983) @@ -21,4 +21,7 @@ status = insert_FWobjects_dvopsps (); } + if (!strcasecmp (argv[1], "forced_galaxy_shape")) { + status = insert_FGshape_dvopsps (); + } if (!strcasecmp (argv[1], "skytable")) { status = insert_skytable (); @@ -30,5 +33,9 @@ } - if (!status) exit (1); + if (!status) { + fprintf (stdout, "ERROR running dvopsps!\n"); + exit (1); + } + fprintf (stdout, "SUCCESS!\n"); exit (0); } Index: /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/src/dvopsps_client.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/src/dvopsps_client.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/src/dvopsps_client.c (revision 38983) @@ -28,4 +28,7 @@ status = insert_FWobjects_dvopsps (); } + if (!strcasecmp (argv[1], "forced_galaxy_shape")) { + status = insert_FGshape_dvopsps (); + } if (status == -1) { fprintf (stderr, "invalid mode, should be : detections, skytable\n"); Index: /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/src/initialize_dvopsps.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/src/initialize_dvopsps.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/src/initialize_dvopsps.c (revision 38983) @@ -21,5 +21,5 @@ fprintf (stderr, "\n"); - fprintf (stderr, " (mode) is one of detections, objects, forced_warp_objects, skytable\n"); + fprintf (stderr, " (mode) is one of detections, objects, forced_warp_objects, forced_galaxy_shaoe, skytable\n"); fprintf (stderr, "\n"); Index: /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/src/insert_FGshape_dvopsps.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/src/insert_FGshape_dvopsps.c (revision 38983) +++ /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/src/insert_FGshape_dvopsps.c (revision 38983) @@ -0,0 +1,173 @@ +# include "dvopsps.h" +# define USE_MYSQL 1 + +// determine the relevant catalogs, launch parallel clients if desired +int insert_FGshape_dvopsps () { + + SkyTable *sky = NULL; + SkyList *skylist = NULL; + Catalog catalog; + off_t i; + + // load the current sky table (layout of all SkyRegions) + sky = SkyTableLoadOptimal (CATDIR, NULL, NULL, FALSE, -1, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + + // determine the populated SkyRegions overlapping the requested area (default depth) + if (SINGLE_CPT) { + skylist = SkyRegionByCPT (sky, SINGLE_CPT); + } else { + skylist = SkyListByPatch (sky, -1, &UserPatch); + } + myAssert (skylist, "ooops!"); + + // pass the skylist to the parallel version, and only launch on relevant hosts + if (PARALLEL && !HOST_ID) { + int status = insert_FGshape_dvopsps_parallel (skylist); + return status; + } + +# if (USE_MYSQL) + // NOTE: mysql connection happens here since each dvopsps_client makes its own connection + MYSQL mysqlBase; + MYSQL *mysqlReal = mysql_dvopsps_connect (&mysqlBase); + if (!mysqlReal) { + fprintf (stderr, "failed to connect to mysql\n"); + exit (1); + } +# else + MYSQL *mysqlReal = NULL; +# endif + + // select measurements for each populated catalog + for (i = 0; i < skylist[0].Nregions; i++) { + + // does this host ID match the desired location for the table? + if (!HostTableTestHost(skylist[0].regions[i], HOST_ID)) continue; + + char localFilename[1024]; + snprintf (localFilename, 1024, "%s/%s.cpt", HOSTDIR, skylist->regions[i]->name); + + // set up the basic catalog info + dvo_catalog_init (&catalog, TRUE); + catalog.filename = HOST_ID ? localFilename : skylist[0].filename[i]; + catalog.catflags = DVO_LOAD_AVERAGE | DVO_LOAD_GALPHOT; + + if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "r")) { + fprintf (stderr, "ERROR: failure reading catalog %s\n", catalog.filename); + exit (1); + } + if (!catalog.Naverage_disk) { + if (VERBOSE) fprintf (stderr, "no data in %s, skipping\n", catalog.filename); + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + continue; + } + + // NOTE: this is where the real action happens + insert_FGshape_dvopsps_catalog (&catalog, skylist->regions[i]->name, mysqlReal); + + // NOTE : unlike setastrom or relphot, this program is read-only wrt dvo + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + } + + if (0) { + char query[256]; + MYSQL_RES *result; + + sprintf (query, "commit;"); + if (mysql_query (mysqlReal, query)) { + fprintf (stderr, "failed to turn off autocommit\n"); + fprintf (stderr, "%s\n", mysql_error (mysqlReal)); + return FALSE; + } + result = mysql_store_result (mysqlReal); + mysql_free_result (result); + } + + return (TRUE); +} + +# define DEBUG 1 + +int insert_FGshape_dvopsps_parallel (SkyList *sky) { + + // launch the dvopsps_client jobs to the parallel hosts + + // load the list of hosts + HostTable *table = HostTableLoad (CATDIR, sky->hosts); + + // if sky->Nregions < 0.5*Nhosts, check if the host used before launching... + + int i, j; + for (i = 0; i < table->Nhosts; i++) { + + if (sky->Nregions < table->Nhosts) { + // do any of the regions want this host? + int wantThisHost = FALSE; + for (j = 0; j < sky->Nregions; j++) { + if (HostTableTestHost (sky->regions[j], table->hosts[i].hostID)) { + wantThisHost = TRUE; + break; + } + } + if (!wantThisHost) { + // fprintf (stderr, "skip host %s\n", table->hosts[i].hostname); + continue; + } + } + + // ensure that the paths are absolute path names + char *tmppath = abspath (table->hosts[i].pathname, DVO_MAX_PATH); + free (table->hosts[i].pathname); + table->hosts[i].pathname = tmppath; + + char command[1024]; + snprintf (command, 1024, "dvopsps_client forced_galaxy_shape -hostID %d -catdir %s -hostdir %s -region %f %f %f %f", + table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, + UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); + + char tmpline[1024]; + snprintf (tmpline, 1024, "%s -dbhost %s", command, DATABASE_HOST); strcpy (command, tmpline); + snprintf (tmpline, 1024, "%s -dbuser %s", command, DATABASE_USER); strcpy (command, tmpline); + snprintf (tmpline, 1024, "%s -dbpass %s", command, DATABASE_PASS); strcpy (command, tmpline); + snprintf (tmpline, 1024, "%s -dbname %s", command, DATABASE_NAME); strcpy (command, tmpline); + + if (VERBOSE) { snprintf (tmpline, 1024, "%s -v", command); strcpy (command, tmpline); } + if (SINGLE_CPT) { snprintf (tmpline, 1024, "%s -cpt %s", command, SINGLE_CPT); strcpy (command, tmpline); } + + fprintf (stderr, "command: %s\n", command); + + if (PARALLEL_MANUAL) continue; + + if (PARALLEL_SERIAL) { + int status = system (command); + if (status) { + fprintf (stderr, "ERROR running dvopsps_client\n"); + exit (2); + } + } else { + // launch the job on the remote machine (no handshake) + int errorInfo = 0; + int pid = rconnect ("ssh", table->hosts[i].hostname, command, table->hosts[i].stdio, &errorInfo, FALSE); + if (!pid) { + if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", table->hosts[i].hostname, errorInfo); + exit (1); + } + table->hosts[i].pid = pid; // save for future reference + } + } + + // wait for the remote jobs to be completed + if (PARALLEL_MANUAL) { + fprintf (stderr, "run the dvopsps_client commands above. when these are done, hit return\n"); + getchar(); + } + if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) { + int status = HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE); + if (!status) return FALSE; + } + + return (TRUE); +} Index: /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/src/insert_FGshape_dvopsps_catalog.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/src/insert_FGshape_dvopsps_catalog.c (revision 38983) +++ /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/src/insert_FGshape_dvopsps_catalog.c (revision 38983) @@ -0,0 +1,289 @@ +# include "dvopsps.h" +# define DEBUG 0 + +int insert_FGshape_dvopsps_catalog (Catalog *catalog, char *basename, MYSQL *mysql) { + + off_t i; + int missingID = 0; + int noGalphot = 0; + + IOBuffer ave_buffer, gal_buffer; + ave_buffer.Nalloc = 0; + gal_buffer.Nalloc = 0; + + Average *average = catalog->average; + GalPhot *galphot = catalog->galphot; + + off_t found = 0; + + char *cleanname = strcreate (basename); + for (i = 0; i < strlen(cleanname); i++) { + if (cleanname[i] == '.') cleanname[i] = '_'; + if (cleanname[i] == '/') cleanname[i] = '_'; + } + insert_FGshape_mysql_create_tables (cleanname, mysql); + + INITTIME; + insert_FGshape_mysql_init (&ave_buffer, &gal_buffer, cleanname); + int Ninsert = 0; + + // NOTE for testing, just do a few objects + if (!mysql) { + catalog[0].Naverage = 5; + } + for (i = 0; i < catalog[0].Naverage; i++) { + + if (average[i].extID == 0) { + missingID ++; + } + + off_t m = average[i].galphotOffset; + off_t Ngalphot = average[i].Ngalphot; + + // skip objects with no lensing object + if (!Ngalphot) { + noGalphot ++; + continue; + } + + // XXX check return status + insert_FGshape_mysql_value (&ave_buffer, &gal_buffer, &average[i], &galphot[m], Ngalphot); + + // Average is a big table, but Ngalphot*galphot might be bigger, so check both + int fullBuffer = FALSE; + fullBuffer = fullBuffer || (ave_buffer.Nbuffer > MAX_BUFFER); + fullBuffer = fullBuffer || (gal_buffer.Nbuffer > MAX_BUFFER); + + if (fullBuffer) { + insert_FGshape_mysql_commit (&ave_buffer, &gal_buffer, mysql); + if (DEBUG) fprintf (stderr, "inserted %d rows\n", Ninsert); + Ninsert = 0; + // empty buffers but do not use FlushBuffer as that resets the allocation + ave_buffer.Nbuffer = 0; + bzero (ave_buffer.buffer, ave_buffer.Nalloc); + gal_buffer.Nbuffer = 0; + bzero (gal_buffer.buffer, gal_buffer.Nalloc); + insert_FGshape_mysql_init (&ave_buffer, &gal_buffer, cleanname); + } + Ninsert ++; + found ++; + } + if (Ninsert) { + insert_FGshape_mysql_commit (&ave_buffer, &gal_buffer, mysql); + } + if (VERBOSE) fprintf (stderr, "inserted "OFF_T_FMT" average, "OFF_T_FMT" galphot\n", + catalog[0].Naverage, catalog[0].Ngalphot); + FreeIOBuffer (&ave_buffer); + FreeIOBuffer (&gal_buffer); + free (cleanname); + + MARKTIME("-- inserted "OFF_T_FMT" objects in %f sec, skipped %d without IDs, %d with no galphot\n", found, dtime, missingID, noGalphot); + return (TRUE); +} + +int insert_FGshape_mysql_create_tables (char *basename, MYSQL *mysql) { + + int status; + IOBuffer buffer; + InitIOBuffer (&buffer, 1024); + // MYSQL_RES *result; + + // *** Average / cpt *** + PrintIOBuffer (&buffer, "DROP TABLE %s_cpt\n", basename); + status = mysql_query(mysql, buffer.buffer); + if (status) { + fprintf (stderr, "failed to drop table:\n"); + fprintf (stderr, "%s\n", mysql_error(mysql)); + } + buffer.Nbuffer = 0; + bzero (buffer.buffer, buffer.Nalloc); + + // Only send the necessary fields (eg, do not sent parallax and pm) + PrintIOBuffer (&buffer, "CREATE TABLE %s_cpt (" + "EXT_ID BIGINT, " + "OBJ_ID INT, " + "CAT_ID INT " + ")\n", basename); + + if (DEBUG) fprintf (stderr, "%s\n", buffer.buffer); + status = mysql_query(mysql, buffer.buffer); + if (status) { + fprintf (stderr, "failed to create table:\n"); + fprintf (stderr, "%s\n", mysql_error(mysql)); + fprintf (stderr, "Nbuffer: %d\n", buffer.Nbuffer); + } + buffer.Nbuffer = 0; + bzero (buffer.buffer, buffer.Nalloc); + + // *** Galphot / cpq *** + PrintIOBuffer (&buffer, "DROP TABLE %s_cpq\n", basename); + status = mysql_query(mysql, buffer.buffer); + if (status) { + fprintf (stderr, "failed to drop table:\n"); + fprintf (stderr, "%s\n", mysql_error(mysql)); + } + buffer.Nbuffer = 0; + bzero (buffer.buffer, buffer.Nalloc); + + PrintIOBuffer (&buffer, "CREATE TABLE %s_cpq (" + "OBJ_ID INT, " + "CAT_ID INT, " + "DET_ID INT, " + "PHOTCODE SMALLINT, " + "MODEL_TYPE SMALLINT, " + + "MAG FLOAT, " + "MAG_ERR FLOAT, " + "MAJOR_AXIS FLOAT, " + "MINOR_AXIS FLOAT, " + "MAJOR_AXIS_ERR FLOAT, " + "MINOR_AXIS_ERR FLOAT, " + "THETA FLOAT, " + "THETA_ERR FLOAT, " + "INDEX FLOAT, " + "CHISQ FLOAT " + + "NPIX INT " + ")\n", basename); + + if (DEBUG) fprintf (stderr, "%s\n", buffer.buffer); + status = mysql_query(mysql, buffer.buffer); + if (status) { + fprintf (stderr, "failed to create table:\n"); + fprintf (stderr, "%s\n", mysql_error(mysql)); + fprintf (stderr, "Nbuffer: %d\n", buffer.Nbuffer); + } + + FreeIOBuffer (&buffer); + return TRUE; +} + +int insert_FGshape_mysql_init (IOBuffer *ave_buffer, IOBuffer *gal_buffer, char *basename) { + + // *** Average / cpt *** + if (ave_buffer->Nalloc == 0) { + InitIOBuffer (ave_buffer, 1024); + } + PrintIOBuffer (ave_buffer, "INSERT INTO %s_cpt (" + "EXT_ID, " + "OBJ_ID, " + "CAT_ID " + ") VALUES \n", basename); + + // *** Galphot / cpq *** + if (gal_buffer->Nalloc == 0) { + InitIOBuffer (gal_buffer, 1024); + } + PrintIOBuffer (gal_buffer, "INSERT INTO %s_cpq (" + "OBJ_ID, " + "CAT_ID, " + "DET_ID, " + "NPIX, " + + "PHOTCODE, " + "MODEL_TYPE, " + "MAG, " + "MAG_ERR, " + + "MAJOR_AXIS, " + "MINOR_AXIS, " + "MAJOR_AXIS_ERR, " + "MINOR_AXIS_ERR, " + + "THETA, " + "THETA_ERR, " + "INDEX, " + "CHISQ" + ") VALUES \n", basename); + + return TRUE; +} + +# define PRINT_FLOAT(BUFFER,FIELD,FORMAT) \ + if (isinf(FIELD) || isnan(FIELD)) PrintIOBuffer (BUFFER, "NULL, "); \ + else PrintIOBuffer (BUFFER, FORMAT, FIELD); + +int insert_FGshape_mysql_value (IOBuffer *ave_buffer, IOBuffer *gal_buffer, Average *average, GalPhot *galphot, int Ngalphot) { + + int i; + + PrintIOBuffer (ave_buffer, " ("); + PrintIOBuffer (ave_buffer, "%lu, ", average->extID); + PrintIOBuffer (ave_buffer, "%u, ", average->objID); + PrintIOBuffer (ave_buffer, "%u ", average->catID); + PrintIOBuffer (ave_buffer, "),\n"); + + for (i = 0; i < Ngalphot; i++) { + PrintIOBuffer (gal_buffer, " ("); + PrintIOBuffer (gal_buffer, "%u, ", galphot->objID); + PrintIOBuffer (gal_buffer, "%u, ", galphot->catID); + PrintIOBuffer (gal_buffer, "%u, ", galphot->detID); + + PrintIOBuffer (gal_buffer, "%hd, ", galphot->photcode); + PrintIOBuffer (gal_buffer, "%hd, ", galphot->modelType); + + PRINT_FLOAT(gal_buffer, galphot->mag, "%.6e, "); + PRINT_FLOAT(gal_buffer, galphot->magErr, "%.6e, "); + PRINT_FLOAT(gal_buffer, galphot->majorAxis, "%.6e, "); + PRINT_FLOAT(gal_buffer, galphot->minorAxis, "%.6e, "); + PRINT_FLOAT(gal_buffer, galphot->majorAxisErr, "%.6e, "); + PRINT_FLOAT(gal_buffer, galphot->minorAxisErr, "%.6e, "); + PRINT_FLOAT(gal_buffer, galphot->theta, "%.6e, "); + PRINT_FLOAT(gal_buffer, galphot->thetaErr, "%.6e, "); + PRINT_FLOAT(gal_buffer, galphot->index, "%.6e, "); + PRINT_FLOAT(gal_buffer, galphot->chisq, "%.6e, "); + + PrintIOBuffer (gal_buffer, "%d ", (int) galphot->Npix); // NOTE: PRINT_FLOAT always adds a trailing comma -- need to use a different method on the last entry + PrintIOBuffer (gal_buffer, "),\n"); + galphot ++; + } + + return TRUE; +} + +int insert_FGshape_mysql_commit (IOBuffer *ave_buffer, IOBuffer *gal_buffer, MYSQL *mysql) { + + MYSQL_RES *result; + + // check that the last two chars are ,\n and replace with ;\n + if (!strcmp(&ave_buffer->buffer[ave_buffer->Nbuffer-2], ",\n")) { + ave_buffer->buffer[ave_buffer->Nbuffer-2] = ';'; + } else { + fprintf (stderr, "invalid sql?\n"); + return FALSE; + } + + // check that the last two chars are ,\n and replace with ;\n + if (!strcmp(&gal_buffer->buffer[gal_buffer->Nbuffer-2], ",\n")) { + gal_buffer->buffer[gal_buffer->Nbuffer-2] = ';'; + } else { + fprintf (stderr, "invalid sql?\n"); + return FALSE; + } + + // XXX check return status + if (mysql) { + int status; + if (DEBUG) fprintf (stderr, "%s\n", ave_buffer->buffer); + status = mysql_query(mysql, ave_buffer->buffer); + if (status) { + fprintf (stderr, "error with insert:\n"); + fprintf (stderr, "%s\n", mysql_error(mysql)); + fprintf (stderr, "Nbuffer: %d\n", ave_buffer->Nbuffer); + } + result = mysql_store_result (mysql); + if (result) mysql_free_result (result); + + if (DEBUG) fprintf (stderr, "%s\n", gal_buffer->buffer); + status = mysql_query(mysql, gal_buffer->buffer); + if (status) { + fprintf (stderr, "error with insert:\n"); + fprintf (stderr, "%s\n", mysql_error(mysql)); + fprintf (stderr, "Nbuffer: %d\n", gal_buffer->Nbuffer); + } + result = mysql_store_result (mysql); + if (result) mysql_free_result (result); + } + + return TRUE; +} Index: /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/src/insert_detections_dvopsps_catalog.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/src/insert_detections_dvopsps_catalog.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/src/insert_detections_dvopsps_catalog.c (revision 38983) @@ -183,5 +183,5 @@ PrintIOBuffer (buffer, "INSERT INTO dvoDetectionFull (objID, detectID, ippObjID, ippDetectID, imageID, catID, "); - PrintIOBuffer (buffer, "ra, dec_, raErr, decErr, zp, telluricExt, airmass, expTime, "); + PrintIOBuffer (buffer, "ra, dec_, raErr, decErr, zp, zpFactor, telluricExt, airmass, expTime, "); PrintIOBuffer (buffer, "Mpsf, dMpsf, Mkron, dMkron, Map, dMap, flags) VALUES \n"); @@ -253,4 +253,5 @@ PRINT_FLOAT(buffer, detection->decErr, "%.6f, "); PRINT_FLOAT(buffer, detection->zp, "%.6f, "); + PRINT_FLOAT(buffer, detection->zpFactor, "%.6e, "); PRINT_FLOAT(buffer, detection->telluricExt, "%.6f, "); PRINT_FLOAT(buffer, detection->airmass, "%.6f, "); @@ -307,4 +308,5 @@ float nominalZP = code->C * 0.001 + code->K * (measure->airmass - 1); float zp = nominalZP - measure->Mcal; + float zpFactor = pow(10.0, -0.4*zp + 3.56); float telluricExt = - measure->Mcal; float expTime = pow(10.0, 0.4 * measure->dt); @@ -324,4 +326,5 @@ detection->zp = zp; + detection->zpFactor = zpFactor; detection->telluricExt = telluricExt; detection->airmass = airmass; Index: /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/src/insert_objects_dvopsps_catalog.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/src/insert_objects_dvopsps_catalog.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/src/insert_objects_dvopsps_catalog.c (revision 38983) @@ -43,18 +43,13 @@ // fprintf (stderr, "%f : %d %d %d\n", average[i].ChiSqPM, isinf(average[i].ChiSqPM), isnan(average[i].ChiSqPM), isfinite(average[i].ChiSqPM)); - // Average is a bigger table, but Nsecfilt*secfilt might be bigger, so check both + // Average is a big table, but Nsecfilt*secfilt might be bigger, so check both if ((ave_buffer.Nbuffer > MAX_BUFFER) || (sec_buffer.Nbuffer > MAX_BUFFER)) { insert_objects_mysql_commit (&ave_buffer, &sec_buffer, mysql); if (DEBUG) fprintf (stderr, "inserted %d rows\n", Ninsert); Ninsert = 0; - if (0) { - FlushIOBuffer (&ave_buffer); - FlushIOBuffer (&sec_buffer); - } else { - ave_buffer.Nbuffer = 0; - bzero (ave_buffer.buffer, ave_buffer.Nalloc); - sec_buffer.Nbuffer = 0; - bzero (sec_buffer.buffer, sec_buffer.Nalloc); - } + ave_buffer.Nbuffer = 0; + bzero (ave_buffer.buffer, ave_buffer.Nalloc); + sec_buffer.Nbuffer = 0; + bzero (sec_buffer.buffer, sec_buffer.Nalloc); insert_objects_mysql_init (&ave_buffer, &sec_buffer, cleanname); } Index: /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/src/insert_skytable.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/src/insert_skytable.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/dvopsps/src/insert_skytable.c (revision 38983) @@ -64,5 +64,5 @@ } - PrintIOBuffer (buffer, "INSERT INTO dvoSkyTable (R_MIN, R_MAX, D_MIN, D_MAX, INDEX_, NAME) VALUES \n"); + PrintIOBuffer (buffer, "INSERT INTO dvoSkyTable (R_MIN, R_MAX, D_MIN, D_MAX, REGION_ID, NAME) VALUES \n"); return TRUE; Index: /branches/eam_branches/ipp-20150625/Ohana/src/getstar/src/MatchImages.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/getstar/src/MatchImages.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/getstar/src/MatchImages.c (revision 38983) @@ -24,9 +24,10 @@ /* for mosaic astrometry, the grid should be w.r.t. the tangent-plane, not chip coords */ - InitCoords (&tcoords, "DEC--ARC"); + InitCoords (&tcoords, "DEC--TAN"); tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; tcoords.crval1 = image[0].coords.crval1; tcoords.crval2 = image[0].coords.crval2; + if (!strcmp (&image[0].coords.ctype[4], "-WRP")) { myAssert (image[0].coords.mosaic, "should already have tested this"); Index: /branches/eam_branches/ipp-20150625/Ohana/src/libautocode/def/average.d =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/libautocode/def/average.d (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/libautocode/def/average.d (revision 38983) @@ -41,5 +41,5 @@ FIELD Nlensing, NLENSING, unsigned short, number of lensing measurements FIELD Nlensobj, NLENSOBJ, unsigned short, number of lensing measurements -FIELD Ngalphot, NGALPHOT, unsigned short, number of galphot measurements +FIELD Ngalphot, NGALPHOT, unsigned short, number of galphot measurements FIELD measureOffset, OFF_MEASURE, int, offset to first psf measurement Index: /branches/eam_branches/ipp-20150625/Ohana/src/libautocode/def/cmf-ps1-v5-lensing-alt.d =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/libautocode/def/cmf-ps1-v5-lensing-alt.d (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/libautocode/def/cmf-ps1-v5-lensing-alt.d (revision 38983) @@ -57,27 +57,24 @@ FIELD M4s, MOMENTS_M4S, float, fourth moment sin(t), pixels^4 -FIELD X11_sm_obj, X11_SM_OBJ, float, lensing smear -FIELD X12_sm_obj, X12_SM_OBJ, float, lensing smear -FIELD X22_sm_obj, X22_SM_OBJ, float, lensing smear -FIELD E1_sm_obj, E1_SM_OBJ, float, lensing smear -FIELD E2_sm_obj, E2_SM_OBJ, float, lensing smear - -FIELD X11_sh_obj, X11_SH_OBJ, float, lensing shear -FIELD X12_sh_obj, X12_SH_OBJ, float, lensing shear -FIELD X22_sh_obj, X22_SH_OBJ, float, lensing shear -FIELD E1_sh_obj, E1_SH_OBJ, float, lensing shear -FIELD E2_sh_obj, E2_SH_OBJ, float, lensing shear - -FIELD X11_sm_psf, X11_SM_PSF, float, lensing smear -FIELD X12_sm_psf, X12_SM_PSF, float, lensing smear -FIELD X22_sm_psf, X22_SM_PSF, float, lensing smear -FIELD E1_sm_psf, E1_SM_PSF, float, lensing smear -FIELD E2_sm_psf, E2_SM_PSF, float, lensing smear - -FIELD X11_sh_psf, X11_SH_PSF, float, lensing shear -FIELD X12_sh_psf, X12_SH_PSF, float, lensing shear -FIELD X22_sh_psf, X22_SH_PSF, float, lensing shear -FIELD E1_sh_psf, E1_SH_PSF, float, lensing shear -FIELD E2_sh_psf, E2_SH_PSF, float, lensing shear +FIELD X11_sm_obj, LENS_X11_SM_OBJ, float, lensing smear +FIELD X12_sm_obj, LENS_X12_SM_OBJ, float, lensing smear +FIELD X22_sm_obj, LENS_X22_SM_OBJ, float, lensing smear +FIELD E1_sm_obj, LENS_E1_SM_OBJ, float, lensing smear +FIELD E2_sm_obj, LENS_E2_SM_OBJ, float, lensing smear +FIELD X11_sh_obj, LENS_X11_SH_OBJ, float, lensing shear +FIELD X12_sh_obj, LENS_X12_SH_OBJ, float, lensing shear +FIELD X22_sh_obj, LENS_X22_SH_OBJ, float, lensing shear +FIELD E1_sh_obj, LENS_E1_SH_OBJ, float, lensing shear +FIELD E2_sh_obj, LENS_E2_SH_OBJ, float, lensing shear +FIELD X11_sm_psf, LENS_X11_SM_PSF, float, lensing smear +FIELD X12_sm_psf, LENS_X12_SM_PSF, float, lensing smear +FIELD X22_sm_psf, LENS_X22_SM_PSF, float, lensing smear +FIELD E1_sm_psf, LENS_E1_SM_PSF, float, lensing smear +FIELD E2_sm_psf, LENS_E2_SM_PSF, float, lensing smear +FIELD X11_sh_psf, LENS_X11_SH_PSF, float, lensing shear +FIELD X12_sh_psf, LENS_X12_SH_PSF, float, lensing shear +FIELD X22_sh_psf, LENS_X22_SH_PSF, float, lensing shear +FIELD E1_sh_psf, LENS_E1_SH_PSF, float, lensing shear +FIELD E2_sh_psf, LENS_E2_SH_PSF, float, lensing shear FIELD Mr1, MOMENTS_R1, float, first radial moment, pixels Index: /branches/eam_branches/ipp-20150625/Ohana/src/libautocode/def/cmf-ps1-v5-lensing.d =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/libautocode/def/cmf-ps1-v5-lensing.d (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/libautocode/def/cmf-ps1-v5-lensing.d (revision 38983) @@ -57,27 +57,24 @@ FIELD M4s, MOMENTS_M4S, float, fourth moment sin(t), pixels^4 -FIELD X11_sm_obj, X11_SM_OBJ, float, lensing smear -FIELD X12_sm_obj, X12_SM_OBJ, float, lensing smear -FIELD X22_sm_obj, X22_SM_OBJ, float, lensing smear -FIELD E1_sm_obj, E1_SM_OBJ, float, lensing smear -FIELD E2_sm_obj, E2_SM_OBJ, float, lensing smear - -FIELD X11_sh_obj, X11_SH_OBJ, float, lensing shear -FIELD X12_sh_obj, X12_SH_OBJ, float, lensing shear -FIELD X22_sh_obj, X22_SH_OBJ, float, lensing shear -FIELD E1_sh_obj, E1_SH_OBJ, float, lensing shear -FIELD E2_sh_obj, E2_SH_OBJ, float, lensing shear - -FIELD X11_sm_psf, X11_SM_PSF, float, lensing smear -FIELD X12_sm_psf, X12_SM_PSF, float, lensing smear -FIELD X22_sm_psf, X22_SM_PSF, float, lensing smear -FIELD E1_sm_psf, E1_SM_PSF, float, lensing smear -FIELD E2_sm_psf, E2_SM_PSF, float, lensing smear - -FIELD X11_sh_psf, X11_SH_PSF, float, lensing shear -FIELD X12_sh_psf, X12_SH_PSF, float, lensing shear -FIELD X22_sh_psf, X22_SH_PSF, float, lensing shear -FIELD E1_sh_psf, E1_SH_PSF, float, lensing shear -FIELD E2_sh_psf, E2_SH_PSF, float, lensing shear +FIELD X11_sm_obj, LENS_X11_SM_OBJ, float, lensing smear +FIELD X12_sm_obj, LENS_X12_SM_OBJ, float, lensing smear +FIELD X22_sm_obj, LENS_X22_SM_OBJ, float, lensing smear +FIELD E1_sm_obj, LENS_E1_SM_OBJ, float, lensing smear +FIELD E2_sm_obj, LENS_E2_SM_OBJ, float, lensing smear +FIELD X11_sh_obj, LENS_X11_SH_OBJ, float, lensing shear +FIELD X12_sh_obj, LENS_X12_SH_OBJ, float, lensing shear +FIELD X22_sh_obj, LENS_X22_SH_OBJ, float, lensing shear +FIELD E1_sh_obj, LENS_E1_SH_OBJ, float, lensing shear +FIELD E2_sh_obj, LENS_E2_SH_OBJ, float, lensing shear +FIELD X11_sm_psf, LENS_X11_SM_PSF, float, lensing smear +FIELD X12_sm_psf, LENS_X12_SM_PSF, float, lensing smear +FIELD X22_sm_psf, LENS_X22_SM_PSF, float, lensing smear +FIELD E1_sm_psf, LENS_E1_SM_PSF, float, lensing smear +FIELD E2_sm_psf, LENS_E2_SM_PSF, float, lensing smear +FIELD X11_sh_psf, LENS_X11_SH_PSF, float, lensing shear +FIELD X12_sh_psf, LENS_X12_SH_PSF, float, lensing shear +FIELD X22_sh_psf, LENS_X22_SH_PSF, float, lensing shear +FIELD E1_sh_psf, LENS_E1_SH_PSF, float, lensing shear +FIELD E2_sh_psf, LENS_E2_SH_PSF, float, lensing shear FIELD srcChipNum, SRC_CHIP_NUM, short, source chip in warp Index: /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/include/dvo.h =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/include/dvo.h (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/include/dvo.h (revision 38983) @@ -1217,3 +1217,20 @@ # undef LENSFIELD +# define GALPHOT_FIELD(NAME, VALUE) float GalphotValue_##NAME (PhotCode *code, dvoMagClassType class, GalPhot *galphot, int Ngalphot) + +GALPHOT_FIELD(GAL_MAG, mag); +GALPHOT_FIELD(GAL_MAG_ERR, magErr); +GALPHOT_FIELD(GAL_MAJ, majorAxis); +GALPHOT_FIELD(GAL_MAJ_ERR, majorAxisErr); +GALPHOT_FIELD(GAL_MIN, minorAxis); +GALPHOT_FIELD(GAL_MIN_ERR, minorAxisErr); +GALPHOT_FIELD(GAL_THETA, theta); +GALPHOT_FIELD(GAL_THETA_ERR, thetaErr); +GALPHOT_FIELD(GAL_INDEX, index); +GALPHOT_FIELD(GAL_CHISQ, chisq); +GALPHOT_FIELD(GAL_NPIX, Npix); +GALPHOT_FIELD(GAL_TYPE, modelType); + +# undef GALPHOT_FIELD + # endif // DVO_H Index: /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/include/dvodb.h =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/include/dvodb.h (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/include/dvodb.h (revision 38983) @@ -97,4 +97,17 @@ MAG_OPTION_E1, MAG_OPTION_E2, + + MAG_OPTION_GAL_MAG, + MAG_OPTION_GAL_MAG_ERR, + MAG_OPTION_GAL_MAJ, + MAG_OPTION_GAL_MAJ_ERR, + MAG_OPTION_GAL_MIN, + MAG_OPTION_GAL_MIN_ERR, + MAG_OPTION_GAL_THETA, + MAG_OPTION_GAL_THETA_ERR, + MAG_OPTION_GAL_INDEX, + MAG_OPTION_GAL_CHISQ, + MAG_OPTION_GAL_NPIX, + MAG_OPTION_GAL_TYPE, } dvoMagOptionType; @@ -107,4 +120,7 @@ MAG_CLASS_KRON, MAG_CLASS_APER, + MAG_CLASS_DEV, // DeVaucouleur Model (only for galphot) + MAG_CLASS_EXP, // Exponential Model (only for galphot) + MAG_CLASS_SER, // Sersic Model (only for galphot) } dvoMagClassType; @@ -484,5 +500,5 @@ int ParseImageField PROTO((dbField *field, char *fieldName)); -dbValue dbExtractAverages PROTO((Average *average, SecFilt *secfilt, Measure *measure, Lensobj *lensobj, StarPar *starpar, dbField *field)); +dbValue dbExtractAverages PROTO((Average *average, SecFilt *secfilt, Measure *measure, Lensobj *lensobj, StarPar *starpar, GalPhot *galphot, dbField *field)); dbValue dbExtractMeasures PROTO((Average *average, SecFilt *secfilt, Measure *measure, Lensing *lensing, StarPar *starpar, dbField *field)); dbValue dbExtractImages PROTO((Image *image, off_t Nimage, off_t N, dbField *field)); @@ -497,4 +513,5 @@ int dbFieldNeedLensing (dbField *fields, int Nfields); int dbFieldNeedStarpar (dbField *fields, int Nfields, int isAverage); +int dbFieldNeedGalphot (dbField *fields, int Nfields); void FreeImageSelection (void); Index: /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/cmf-ps1-v5-r0-lensing.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/cmf-ps1-v5-r0-lensing.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/cmf-ps1-v5-r0-lensing.c (revision 38983) @@ -269,24 +269,24 @@ gfits_define_bintable_column (header, "E", "MOMENTS_M4C", "fourth moment cos(t)", "pixels^4", 1.0, 0.0); gfits_define_bintable_column (header, "E", "MOMENTS_M4S", "fourth moment sin(t)", "pixels^4", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X11_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X12_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X22_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E1_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E2_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X11_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X12_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X22_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E1_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E2_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X11_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X12_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X22_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E1_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E2_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X11_SH_PSF", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X12_SH_PSF", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X22_SH_PSF", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E1_SH_PSF", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E2_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X11_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X12_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X22_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E1_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E2_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X11_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X12_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X22_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E1_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E2_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X11_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X12_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X22_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E1_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E2_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X11_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X12_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X22_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E1_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E2_SH_PSF", "lensing shear", "", 1.0, 0.0); gfits_define_bintable_column (header, "E", "MOMENTS_R1", "first radial moment", "pixels", 1.0, 0.0); gfits_define_bintable_column (header, "E", "MOMENTS_RH", "half radial moment", "pixels^1/2", 1.0, 0.0); @@ -365,24 +365,24 @@ gfits_define_bintable_column (header, "E", "MOMENTS_M4C", "fourth moment cos(t)", "pixels^4", 1.0, 0.0); gfits_define_bintable_column (header, "E", "MOMENTS_M4S", "fourth moment sin(t)", "pixels^4", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X11_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X12_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X22_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E1_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E2_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X11_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X12_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X22_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E1_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E2_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X11_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X12_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X22_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E1_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E2_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X11_SH_PSF", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X12_SH_PSF", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X22_SH_PSF", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E1_SH_PSF", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E2_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X11_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X12_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X22_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E1_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E2_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X11_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X12_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X22_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E1_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E2_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X11_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X12_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X22_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E1_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E2_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X11_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X12_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X22_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E1_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E2_SH_PSF", "lensing shear", "", 1.0, 0.0); gfits_define_bintable_column (header, "E", "MOMENTS_R1", "first radial moment", "pixels", 1.0, 0.0); gfits_define_bintable_column (header, "E", "MOMENTS_RH", "half radial moment", "pixels^1/2", 1.0, 0.0); Index: /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/cmf-ps1-v5-r1-lensing.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/cmf-ps1-v5-r1-lensing.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/cmf-ps1-v5-r1-lensing.c (revision 38983) @@ -275,24 +275,24 @@ gfits_define_bintable_column (header, "E", "MOMENTS_M4C", "fourth moment cos(t)", "pixels^4", 1.0, 0.0); gfits_define_bintable_column (header, "E", "MOMENTS_M4S", "fourth moment sin(t)", "pixels^4", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X11_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X12_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X22_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E1_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E2_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X11_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X12_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X22_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E1_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E2_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X11_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X12_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X22_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E1_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E2_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X11_SH_PSF", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X12_SH_PSF", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X22_SH_PSF", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E1_SH_PSF", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E2_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X11_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X12_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X22_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E1_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E2_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X11_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X12_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X22_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E1_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E2_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X11_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X12_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X22_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E1_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E2_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X11_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X12_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X22_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E1_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E2_SH_PSF", "lensing shear", "", 1.0, 0.0); gfits_define_bintable_column (header, "I", "SRC_CHIP_NUM", "number of source chip for warp", "", 1.0, 0.0); @@ -377,24 +377,24 @@ gfits_define_bintable_column (header, "E", "MOMENTS_M4C", "fourth moment cos(t)", "pixels^4", 1.0, 0.0); gfits_define_bintable_column (header, "E", "MOMENTS_M4S", "fourth moment sin(t)", "pixels^4", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X11_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X12_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X22_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E1_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E2_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X11_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X12_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X22_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E1_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E2_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X11_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X12_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X22_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E1_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E2_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X11_SH_PSF", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X12_SH_PSF", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X22_SH_PSF", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E1_SH_PSF", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E2_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X11_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X12_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X22_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E1_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E2_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X11_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X12_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X22_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E1_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E2_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X11_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X12_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X22_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E1_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E2_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X11_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X12_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X22_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E1_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E2_SH_PSF", "lensing shear", "", 1.0, 0.0); gfits_define_bintable_column (header, "I", "SRC_CHIP_NUM", "number of source chip for warp", "", 1.0, 0.0); Index: /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/cmf-ps1-v5-r2-lensing.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/cmf-ps1-v5-r2-lensing.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/cmf-ps1-v5-r2-lensing.c (revision 38983) @@ -278,26 +278,26 @@ gfits_define_bintable_column (header, "E", "MOMENTS_M4C", "fourth moment cos(t)", "pixels^4", 1.0, 0.0); gfits_define_bintable_column (header, "E", "MOMENTS_M4S", "fourth moment sin(t)", "pixels^4", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X11_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X12_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X22_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E1_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E2_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X11_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X12_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X22_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E1_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E2_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X11_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X12_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X22_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E1_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E2_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X11_SH_PSF", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X12_SH_PSF", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X22_SH_PSF", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E1_SH_PSF", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E2_SH_PSF", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E1_PSF", "psf polarization", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E2_PSF", "psf polarization", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X11_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X12_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X22_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E1_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E2_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X11_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X12_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X22_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E1_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E2_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X11_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X12_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X22_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E1_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E2_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X11_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X12_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X22_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E1_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E2_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E1_PSF", "psf polarization", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E2_PSF", "psf polarization", "", 1.0, 0.0); gfits_define_bintable_column (header, "I", "SRC_CHIP_NUM", "number of source chip for warp", "", 1.0, 0.0); @@ -382,26 +382,26 @@ gfits_define_bintable_column (header, "E", "MOMENTS_M4C", "fourth moment cos(t)", "pixels^4", 1.0, 0.0); gfits_define_bintable_column (header, "E", "MOMENTS_M4S", "fourth moment sin(t)", "pixels^4", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X11_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X12_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X22_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E1_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E2_SM_OBJ", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X11_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X12_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X22_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E1_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E2_SH_OBJ", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X11_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X12_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X22_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E1_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E2_SM_PSF", "lensing smear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X11_SH_PSF", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X12_SH_PSF", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "X22_SH_PSF", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E1_SH_PSF", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E2_SH_PSF", "lensing shear", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E1_PSF", "psf polarization", "", 1.0, 0.0); - gfits_define_bintable_column (header, "E", "E2_PSF", "psf polarization", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X11_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X12_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X22_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E1_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E2_SM_OBJ", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X11_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X12_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X22_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E1_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E2_SH_OBJ", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X11_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X12_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X22_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E1_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E2_SM_PSF", "lensing smear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X11_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X12_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_X22_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E1_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E2_SH_PSF", "lensing shear", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E1_PSF", "psf polarization", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LENS_E2_PSF", "psf polarization", "", 1.0, 0.0); gfits_define_bintable_column (header, "I", "SRC_CHIP_NUM", "number of source chip for warp", "", 1.0, 0.0); Index: /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/dbExtractAverages.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/dbExtractAverages.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/dbExtractAverages.c (revision 38983) @@ -50,5 +50,5 @@ /* return average.field based on the selection */ -dbValue dbExtractAverages (Average *average, SecFilt *secfilt, Measure *measure, Lensobj *lensobj, StarPar *starpar, dbField *field) { +dbValue dbExtractAverages (Average *average, SecFilt *secfilt, Measure *measure, Lensobj *lensobj, StarPar *starpar, GalPhot *galphot, dbField *field) { // off_t i; @@ -408,4 +408,17 @@ break; + case MAG_OPTION_GAL_MAG: { value.Flt = GalphotValue_GAL_MAG (field->photcode, field->magClass, galphot, average->Ngalphot); break; } + case MAG_OPTION_GAL_MAG_ERR: { value.Flt = GalphotValue_GAL_MAG_ERR (field->photcode, field->magClass, galphot, average->Ngalphot); break; } + case MAG_OPTION_GAL_MAJ: { value.Flt = GalphotValue_GAL_MAJ (field->photcode, field->magClass, galphot, average->Ngalphot); break; } + case MAG_OPTION_GAL_MAJ_ERR: { value.Flt = GalphotValue_GAL_MAJ_ERR (field->photcode, field->magClass, galphot, average->Ngalphot); break; } + case MAG_OPTION_GAL_MIN: { value.Flt = GalphotValue_GAL_MIN (field->photcode, field->magClass, galphot, average->Ngalphot); break; } + case MAG_OPTION_GAL_MIN_ERR: { value.Flt = GalphotValue_GAL_MIN_ERR (field->photcode, field->magClass, galphot, average->Ngalphot); break; } + case MAG_OPTION_GAL_THETA: { value.Flt = GalphotValue_GAL_THETA (field->photcode, field->magClass, galphot, average->Ngalphot); break; } + case MAG_OPTION_GAL_THETA_ERR: { value.Flt = GalphotValue_GAL_THETA_ERR(field->photcode, field->magClass, galphot, average->Ngalphot); break; } + case MAG_OPTION_GAL_INDEX: { value.Flt = GalphotValue_GAL_INDEX (field->photcode, field->magClass, galphot, average->Ngalphot); break; } + case MAG_OPTION_GAL_CHISQ: { value.Flt = GalphotValue_GAL_CHISQ (field->photcode, field->magClass, galphot, average->Ngalphot); break; } + case MAG_OPTION_GAL_NPIX: { value.Flt = GalphotValue_GAL_NPIX (field->photcode, field->magClass, galphot, average->Ngalphot); break; } + case MAG_OPTION_GAL_TYPE: { value.Flt = GalphotValue_GAL_TYPE (field->photcode, field->magClass, galphot, average->Ngalphot); break; } + default: break; @@ -413,4 +426,5 @@ break; + // this section assumes a 1-to-1 match between average and starpar case AVE_E_BV: { value.Flt = starpar->Ebv; break; } case AVE_E_BV_ERR: { value.Flt = starpar->dEbv; break; } Index: /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/dbExtractMeasures.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/dbExtractMeasures.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/dbExtractMeasures.c (revision 38983) @@ -287,4 +287,17 @@ case MAG_OPTION_E2: + case MAG_OPTION_GAL_MAG: + case MAG_OPTION_GAL_MAG_ERR: + case MAG_OPTION_GAL_MAJ: + case MAG_OPTION_GAL_MAJ_ERR: + case MAG_OPTION_GAL_MIN: + case MAG_OPTION_GAL_MIN_ERR: + case MAG_OPTION_GAL_THETA: + case MAG_OPTION_GAL_THETA_ERR: + case MAG_OPTION_GAL_INDEX: + case MAG_OPTION_GAL_CHISQ: + case MAG_OPTION_GAL_NPIX: + case MAG_OPTION_GAL_TYPE: + case MAG_OPTION_NONE: break; Index: /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/dbFields.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/dbFields.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/dbFields.c (revision 38983) @@ -161,4 +161,25 @@ } +int dbFieldNeedGalphot (dbField *fields, int Nfields) { + int i; + for (i = 0; i < Nfields; i++) { + if (fields[i].magOption == MAG_OPTION_NONE) continue; // non-measure fields do not have a photcode + + if (fields[i].magOption == MAG_OPTION_GAL_MAG ) return TRUE; + if (fields[i].magOption == MAG_OPTION_GAL_MAG_ERR ) return TRUE; + if (fields[i].magOption == MAG_OPTION_GAL_MAJ ) return TRUE; + if (fields[i].magOption == MAG_OPTION_GAL_MAJ_ERR ) return TRUE; + if (fields[i].magOption == MAG_OPTION_GAL_MIN ) return TRUE; + if (fields[i].magOption == MAG_OPTION_GAL_MIN_ERR ) return TRUE; + if (fields[i].magOption == MAG_OPTION_GAL_THETA ) return TRUE; + if (fields[i].magOption == MAG_OPTION_GAL_THETA_ERR ) return TRUE; + if (fields[i].magOption == MAG_OPTION_GAL_INDEX ) return TRUE; + if (fields[i].magOption == MAG_OPTION_GAL_CHISQ ) return TRUE; + if (fields[i].magOption == MAG_OPTION_GAL_NPIX ) return TRUE; + if (fields[i].magOption == MAG_OPTION_GAL_TYPE ) return TRUE; + } + return FALSE; +} + dvoMagSourceType GetMagSource (char *string) { @@ -193,4 +214,10 @@ if (!strcasecmp (string, "aper")) return (MAG_CLASS_APER); if (!strcasecmp (string, "ap")) return (MAG_CLASS_APER); + + // these model classes are only valid for galphot mag options + if (!strcasecmp (string, "dev")) return (MAG_CLASS_DEV); + if (!strcasecmp (string, "exp")) return (MAG_CLASS_EXP); + if (!strcasecmp (string, "ser")) return (MAG_CLASS_SER); + return MAG_CLASS_NONE; @@ -247,4 +274,17 @@ if (!strcasecmp (string, "E2")) return MAG_OPTION_E2; + if (!strcasecmp (string, "GAL_MAG" )) return MAG_OPTION_GAL_MAG ; + if (!strcasecmp (string, "GAL_MAG_ERR")) return MAG_OPTION_GAL_MAG_ERR ; + if (!strcasecmp (string, "GAL_MAJ" )) return MAG_OPTION_GAL_MAJ ; + if (!strcasecmp (string, "GAL_MAJ_ERR")) return MAG_OPTION_GAL_MAJ_ERR ; + if (!strcasecmp (string, "GAL_MIN" )) return MAG_OPTION_GAL_MIN ; + if (!strcasecmp (string, "GAL_MIN_ERR")) return MAG_OPTION_GAL_MIN_ERR ; + if (!strcasecmp (string, "GAL_THETA" )) return MAG_OPTION_GAL_THETA ; + if (!strcasecmp (string, "GAL_THETA_ERR")) return MAG_OPTION_GAL_THETA_ERR ; + if (!strcasecmp (string, "GAL_INDEX" )) return MAG_OPTION_GAL_INDEX ; + if (!strcasecmp (string, "GAL_CHISQ" )) return MAG_OPTION_GAL_CHISQ ; + if (!strcasecmp (string, "GAL_NPIX" )) return MAG_OPTION_GAL_NPIX ; + if (!strcasecmp (string, "GAL_TYPE" )) return MAG_OPTION_GAL_TYPE ; + return MAG_OPTION_NONE; } @@ -373,4 +413,8 @@ free (word); return FALSE; + } + + if (dbFieldNeedGalphot (field, 1)) { + field->magClass = MAG_CLASS_SER; // for galphot, we need to set a default model class (SERSIC) } Index: /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/dvo_photcode_ops.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/dvo_photcode_ops.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/dvo_photcode_ops.c (revision 38983) @@ -1009,5 +1009,7 @@ // mag_AB = -2.5 log (flux_cgs) - 48.6 (by definition) [~Vega flux in V-band] // flux_Jy = flux_cgs * 10^23 + // -2.5 log (flux_cgs) = -2.5 log (flux_Jy * 10^-23) = -2.5 log (flux_Jy) + 2.5*23 + // -2.5 log (flux_cgs) = -2.5 log (flux_Jy) + 57.5 // from these we get: @@ -1017,5 +1019,5 @@ // mag_AB = mag_inst + ZP - // log(flux_Jy) = -0.4*mag_inst -0.4*ZP + 3.56 + // log(flux_Jy) = -0.4*mag_inst - 0.4*ZP + 3.56 // mag_inst = -2.5 log (flux_Jy) - ZP + 8.9 @@ -1023,7 +1025,13 @@ // -2.5 log (flux_inst) = -2.5 log (flux_Jy) - ZP + 8.9 + // log (flux_inst) = log (flux_Jy) + 0.4 ZP - 3.56 + // flux_inst = flux_Jy * 10^(0.4*ZP - 3.56) + // flux_Jy = flux_inst * 10^(3.56 - 0.4*ZP) = flux_inst * zpFactor + + // flux_Jy = flux_inst * 10^(-0.4*(ZP - 8.9)) = flux_inst * zpFactor + // flux_Jy = flux_inst * 3630.8 * 10^(-0.4*ZP) @@ -1186,11 +1194,11 @@ switch (class) { case MAG_CLASS_PSF: - Fave = MagToFlux(secfilt[Ns].M); + Fave = MagToFlux(secfilt[Ns].M - 8.9); break; case MAG_CLASS_KRON: - Fave = MagToFlux(secfilt[Ns].Mkron); + Fave = MagToFlux(secfilt[Ns].Mkron - 8.9); break; case MAG_CLASS_APER: - Fave = MagToFlux(secfilt[Ns].Map); + Fave = MagToFlux(secfilt[Ns].Map - 8.9); break; default: @@ -1201,11 +1209,11 @@ switch (class) { case MAG_CLASS_PSF: - Fave = isnan (secfilt[Ns].FpsfWrp) ? MagToFlux(secfilt[Ns].MpsfWrp) : secfilt[Ns].FpsfWrp; + Fave = isnan (secfilt[Ns].FpsfWrp) ? MagToFlux(secfilt[Ns].MpsfWrp - 8.9) : secfilt[Ns].FpsfWrp; break; case MAG_CLASS_KRON: - Fave = isnan (secfilt[Ns].FkronWrp) ? MagToFlux(secfilt[Ns].MkronWrp) : secfilt[Ns].FkronWrp; + Fave = isnan (secfilt[Ns].FkronWrp) ? MagToFlux(secfilt[Ns].MkronWrp - 8.9) : secfilt[Ns].FkronWrp; break; case MAG_CLASS_APER: - Fave = isnan (secfilt[Ns].FapWrp) ? MagToFlux(secfilt[Ns].MapWrp) : secfilt[Ns].FapWrp; + Fave = isnan (secfilt[Ns].FapWrp) ? MagToFlux(secfilt[Ns].MapWrp - 8.9) : secfilt[Ns].FapWrp; break; default: @@ -1216,11 +1224,11 @@ switch (class) { case MAG_CLASS_PSF: - Fave = isnan (secfilt[Ns].FpsfStk) ? MagToFlux(secfilt[Ns].MpsfStk) : secfilt[Ns].FpsfStk; + Fave = isnan (secfilt[Ns].FpsfStk) ? MagToFlux(secfilt[Ns].MpsfStk - 8.9) : secfilt[Ns].FpsfStk; break; case MAG_CLASS_KRON: - Fave = isnan (secfilt[Ns].FkronStk) ? MagToFlux(secfilt[Ns].MkronStk) : secfilt[Ns].FkronStk; + Fave = isnan (secfilt[Ns].FkronStk) ? MagToFlux(secfilt[Ns].MkronStk - 8.9) : secfilt[Ns].FkronStk; break; case MAG_CLASS_APER: - Fave = isnan (secfilt[Ns].FapStk) ? MagToFlux(secfilt[Ns].MapStk) : secfilt[Ns].FapStk; + Fave = isnan (secfilt[Ns].FapStk) ? MagToFlux(secfilt[Ns].MapStk - 8.9) : secfilt[Ns].FapStk; break; default: @@ -1246,11 +1254,11 @@ switch (class) { case MAG_CLASS_PSF: - dFave = secfilt[Ns].dM * MagToFlux(secfilt[Ns].M); + dFave = secfilt[Ns].dM * MagToFlux(secfilt[Ns].M - 8.9); break; case MAG_CLASS_KRON: - dFave = secfilt[Ns].dMkron * MagToFlux(secfilt[Ns].Mkron); + dFave = secfilt[Ns].dMkron * MagToFlux(secfilt[Ns].Mkron - 8.9); break; case MAG_CLASS_APER: - dFave = secfilt[Ns].dMap * MagToFlux(secfilt[Ns].Map); + dFave = secfilt[Ns].dMap * MagToFlux(secfilt[Ns].Map - 8.9); break; default: @@ -1880,2 +1888,70 @@ LENSFIELD(E1); LENSFIELD(E2); + +# if (0) +float GalphotValue_GAL_MAG (PhotCode *code, dvoMagClassType class, GalPhot *galphot, int Ngalphot) { + int n; + if (code == NULL) return NAN; + for (n = 0; n < Ngalphot; n++) { + short equivCode = GetPhotcodeEquivCodebyCode (galphot[n].photcode); + if (!equivCode) continue; + if (equivCode != code->code) continue; + switch (class) { + case MAG_CLASS_SER: + if (galphot[n].modelType != 5) continue; + break; + case MAG_CLASS_EXP: + if (galphot[n].modelType != 6) continue; + break; + case MAG_CLASS_DEV: + if (galphot[n].modelType != 7) continue; + break; + default: + return NAN; + } + float value = galphot[n].mag; + return (value); + } + return NAN; +} +# endif + +# define GALPHOT_FIELD(NAME, VALUE) \ + float GalphotValue_##NAME (PhotCode *code, dvoMagClassType class, GalPhot *galphot, int Ngalphot) { \ + int n; \ + if (code == NULL) return NAN; \ + for (n = 0; n < Ngalphot; n++) { \ + short equivCode = GetPhotcodeEquivCodebyCode (galphot[n].photcode); \ + if (!equivCode) continue; \ + if (equivCode != code->code) continue; \ + switch (class) { \ + case MAG_CLASS_SER: \ + if (galphot[n].modelType != 5) continue; \ + break; \ + case MAG_CLASS_EXP: \ + if (galphot[n].modelType != 6) continue; \ + break; \ + case MAG_CLASS_DEV: \ + if (galphot[n].modelType != 7) continue; \ + break; \ + default: \ + return NAN; \ + } \ + float value = galphot[n].VALUE; \ + return (value); \ + } \ + return NAN; \ + } + +GALPHOT_FIELD(GAL_MAG, mag) +GALPHOT_FIELD(GAL_MAG_ERR, magErr) +GALPHOT_FIELD(GAL_MAJ, majorAxis) +GALPHOT_FIELD(GAL_MAJ_ERR, majorAxisErr) +GALPHOT_FIELD(GAL_MIN, minorAxis) +GALPHOT_FIELD(GAL_MIN_ERR, minorAxisErr) +GALPHOT_FIELD(GAL_THETA, theta) +GALPHOT_FIELD(GAL_THETA_ERR, thetaErr) +GALPHOT_FIELD(GAL_INDEX, index) +GALPHOT_FIELD(GAL_CHISQ, chisq) +GALPHOT_FIELD(GAL_NPIX, Npix) +GALPHOT_FIELD(GAL_TYPE, modelType) Index: /branches/eam_branches/ipp-20150625/Ohana/src/opihi/cmd.astro/fitpm_irls.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/opihi/cmd.astro/fitpm_irls.c (revision 38983) +++ /branches/eam_branches/ipp-20150625/Ohana/src/opihi/cmd.astro/fitpm_irls.c (revision 38983) @@ -0,0 +1,492 @@ +# include "astro.h" +# define J2000 51544.5 /* Modified Julian date at standard epoch J2000 */ + +# define ESCAPE(MSG,...) { \ + gprint (GP_ERR, MSG, __VA_ARGS__); \ + return FALSE; } + +typedef struct { + double Ro, dRo; + double Do, dDo; + + double uR, duR; + double uD, duD; + + double chisq; + int Nfit; +} PMFit_IRLS; + +int FitPMonly_IRLS (PMFit_IRLS *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE); +int IRLS_converged (PMFit_IRLS *fit); + +int fitpm_irls (int argc, char **argv) { + + int i, N; + + Vector *rvec, *dvec, *tvec, *dRvec, *dDvec; + + Vector *mvec = NULL; // mask vector + if ((N = get_argument (argc, argv, "-mask"))) { + remove_argument (N, &argc, argv); + if ((mvec = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + CastVector (mvec, OPIHI_INT); + } + + int VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + remove_argument (N, &argc, argv); + VERBOSE = TRUE; + } + if ((N = get_argument (argc, argv, "-vv"))) { + remove_argument (N, &argc, argv); + VERBOSE = 2; + } + + if (argc != 6) { + gprint (GP_ERR, "USAGE: fitplx_irls (ra) (dR) (dec) (dD) (mjd) [-mask mask]\n"); + // what about the errors? + return (FALSE); + } + + /* select input / output buffers */ + // if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE); + // Nx = buf[0].header.Naxis[0]; + // Ny = buf[0].header.Naxis[1]; + + if ((rvec = SelectVector (argv[1], ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", argv[1]); + if ((dvec = SelectVector (argv[3], ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", argv[3]); + if ((tvec = SelectVector (argv[5], ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", argv[5]); + + if ((dRvec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", argv[2]); + if ((dDvec = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", argv[4]); + + double *R = rvec->elements.Flt; + double *D = dvec->elements.Flt; + double *T = tvec->elements.Flt; + + double *dR = dRvec->elements.Flt; + double *dD = dDvec->elements.Flt; + + int *mask = NULL; + if (mvec) { + mask = mvec->elements.Int; + } + + N = tvec->Nelements; // XXX check other lengths + + // find mean values to remove + double Npts = 0; + double Tmean = 0; + double Rmean = 0; + double Dmean = 0; + double Tmin = +1000000; + double Tmax = -1000000; + for (i = 0; i < N; i++) { + if (mask && !mask[i]) continue; + Rmean += R[i]; + Dmean += D[i]; + Tmean += T[i]; + Tmin = MIN(Tmin, T[i]); + Tmax = MAX(Tmax, T[i]); + Npts += 1.0; + } + Rmean /= Npts; + Dmean /= Npts; + Tmean /= Npts; + + float Trange = Tmax - Tmin; + // fprintf (stderr, "R,D : %f,%f, T: %f, Trange: %f, Tmin: %f, Tmax: %f\n", Rmean, Dmean, Tmean, Trange, Tmin, Tmax); + + /* project coordinates to a plane centered on the object with units of arcsec */ + Coords coords; + InitCoords (&coords, "DEC--SIN"); + coords.crval1 = Rmean; + coords.crval2 = Dmean; + coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0; + + double *X, *Y, *t, *dX, *dY; + ALLOCATE (X, double, N); + ALLOCATE (Y, double, N); + ALLOCATE (dX, double, N); + ALLOCATE (dY, double, N); + ALLOCATE (t, double, N); + + int n = 0; + for (i = 0; i < N; i++) { + if (mask && !mask[i]) continue; + RD_to_XY (&X[n], &Y[n], R[i], D[i], &coords); + dX[n] = dR[i]; + dY[n] = dD[i]; + t[n] = (T[i] - Tmean) / 365.25; + n++; + } + + PMFit_IRLS fit; + if (!FitPMonly_IRLS (&fit, X, dX, Y, dY, t, n, VERBOSE)) { + return FALSE; + } + + // fprintf (stderr, "Roff, Doff: %f, %f; dRo, dDo: %f, %f\n", fit.Ro, fit.Do, fit.dRo, fit.dDo); + + XY_to_RD (&Rmean, &Dmean, fit.Ro, fit.Do, &coords); + if (VERBOSE) { + fprintf (stderr, "Ro, Do: %f, %f +/- %f, %f\n", Rmean, Dmean, fit.dRo, fit.dDo); + fprintf (stderr, "uR, uD: %f, %f; duR, duD: %f, %f\n", fit.uR, fit.uD, fit.duR, fit.duD); + fprintf (stderr, "chisq: %f Nfit %d\n", fit.chisq, fit.Nfit); + } + + set_variable ("RA", Rmean); + set_variable ("DEC", Dmean); + set_variable ("dR", fit.dRo); + set_variable ("dD", fit.dDo); + set_variable ("uR", fit.uR); + set_variable ("uD", fit.uD); + set_variable ("duR", fit.duR); + set_variable ("duD", fit.duD); + set_variable ("plx", 0.0); + set_variable ("dplx", 0.0); + + set_variable ("Tmean", Tmean); + set_variable ("Trange", Trange); + set_variable ("Prange", 0.0); + + set_variable ("chisq", fit.chisq); + set_variable ("Nfit", fit.Nfit); + + return (TRUE); +} + +/* do we want an init function which does the alloc and a clear function to free? */ +int FitPMonly_IRLS (PMFit_IRLS *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE) { + + int i; + + static double **A, **B; + + + double chisq, Xf, Yf; + + double **Cov; + double *Beta, *Beta_prev; + + double sigma_robust, sigma_ols, sigma_final, sigma_hat; + double *Wx, *Wy; + double *rx, *ry; + double *ux, *uy; + int dof = 2 * Npts - 4; + int p = 4; + int n = 2 * Npts; + double tolerance; + int convergence; + int iterations; + + /* do I need to do this as 2 2x2 matrix equations? */ + if (A == NULL) { + ALLOCATE (A, double *, 4); + ALLOCATE (B, double *, 4); + for (i = 0; i < 4; i++) { + ALLOCATE (A[i], double, 4); + ALLOCATE (B[i], double, 1); + memset (A[i], 0, 4*sizeof(double)); + memset (B[i], 0, 1*sizeof(double)); + } + } + + // things we need + ALLOCATE (Cov, double *, 4); + for (i = 0; i < 4; i++) { + ALLOCATE ( Cov[i], double, 4); + } + + ALLOCATE(Beta, double, 4); + ALLOCATE(Beta_prev, double, 4); + ALLOCATE(Wx, double, Npts); + ALLOCATE(Wy, double, Npts); + ALLOCATE(rx, double, Npts); + ALLOCATE(ry, double, Npts); + ALLOCATE(u, double, Npts); + + // Convert the measurement errors into initial weights. + for (i = 0; i < Npts; i++) { + Wx[i] = 1 / dX[i]; + Wy[i] = 1 / dY[i]; + } + + // Solve OLS equation + if (!weighted_LS(T,X,Wx,Y,Wy,Npts, + A,B,VERBOSE)) { + // Handle fail case + return(FALSE); + } + + // Calculate r vector of residuals and least squares sigma + sigma_ols = 0.0; + for (i = 0; i < Npts; i++) { + rx[i] = X[i] - (T[i] * B[0][0] + B[1][0]); + ry[i] = Y[i] - (T[i] * B[2][0] + B[3][0]); + // u[i] = r[i] / + sigma_ols += SQ(rx[i]) + SQ(ry[i]); + + } + sigma_ols = sqrt(sigma_ols / dof); + + // Save OLS covariance; + for (i = 0; i < 4; i++) { + for (j = 0; j < 4; j++) { + Cov[i][j] = A[i][j]; + } + } + + // Save Beta + for (i = 0; i < 4; i++) { + Beta[i] = B[i][0]; + } + + // Iterately reweight and solve + converged = FALSE; + iterations = 0; + do { + // Save Beta. + for (i = 0; i < 4; i ++) { + Beta_prev[i] = Beta[i]; + } + + // Assign W + for (i = 0; i < Npts; i++) { + Wx[i] = weight_cauchy(rx[i] / dX[i]); + Wy[i] = weight_cauchy(ry[i] / dY[i]); + } + + // Solve + if (!weighted_LS(T,X,Wx,Y,Wy,Npts, + A,B,VERBOSE)) { + // Handle fail case + return(FALSE); + } + + for (i = 0; i < 4; i++) { + Beta[i] = B[i][0]; + } + + // r + sigma_hat = 0.0; + for (i = 0; i < Npts; i++) { + rx[i] = X[i] - (T[i] * B[0][0] + B[1][0]); + ry[i] = Y[i] - (T[i] * B[2][0] + B[3][0]); + u[i] = sqrt(SQ(rx[i] / dX[i]) + SQ(ry[i] / dY[i])); + } + sigma_hat = MAD(u,Npts) / 0.6745; + + // Check convergence + converged = TRUE; + tolerance = 1e-4; // This should probably be tunable. + for (i = 0; i < 4; i++) { + if (fabs(Beta[i] - Beta_prev[i]) > tolerance * abs(Beta[i])) { + converged = FALSE; + } + } + + iterations++; + if (iterations >= 10) { + converged = TRUE; + // Throw a warning or something here. + } + + } while (!converged); + + double ax, ay; + double bx, by; + double lambda; + double sigma_robust_x, sigma_robust_y; + double sigma_final_x, sigma_final_y; + double Sum_Wx, Sum_Wy; + + ax = 0.0; ay = 0.0; + bx = 0.0; by = 0.0; + lambda = 0.0; + for (i = 0; i < Npts; i++) { + Wx[i] = weight_cauchy(rx[i] / dX[i]); + Wy[i] = weight_cauchy(ry[i] / dY[i]); + + ax += dpsi_cauchy(rx[i] / dX[i]); + ay += dpsi_cauchy(ry[i] / dY[i]); + + bx += SQ(Wx[i]); + by += SQ(Wy[i]); + + Sum_Wx += Wx[i]; + Sum_Wy += Wy[i]; + } + ax /= 1.0 * Npts; // mean(psi_dot(r)) + ay /= 1.0 * Npts; + bx /= 1.0 * (Npts - p); // mean(psi^2(r)) * (N / (N-p)) + by /= 1.0 * (Npts - p); + + sigma_robust_x = lamba * sqrt(bx) * sigma_hat * 2.385 / ax; + sigma_robust_y = lamba * sqrt(by) * sigma_hat * 2.385 / ay; + + // This is actually sigma^2, as that's the factor in the covariance (dumouchel 4.1) + sigma_final_x = fmax(SQ(sigma_robust_x), (n * SQ(sigma_robust_x) + SQ(p * sigma_OLS)) / (n + SQ(p))); + sigma_final_y = fmax(SQ(sigma_robust_y), (n * SQ(sigma_robust_y) + SQ(p * sigma_OLS)) / (n + SQ(p))); + + for (i = 0; i < 4; i++) { + for (j = 0; j < 4; j++) { + // This uses the original OLS covariance. + if ((i < 2)&&(j < 2)) { // Upper portion + Cov[i][j] *= sigma_final_x; + } + elsif ((i > 1)&&(j > 1)) { // Lower portion + Cov[i][j] *= sigma_final_y; + } + else { // Cross term + Cov[i][j] *= sqrt(sigma_final_x * sigma_final_y); + } + } + } + + // Finish. + fit[0].Ro = Beta[0]; + fit[0].uR = Beta[1]; + fit[0].Do = Beta[2]; + fit[0].uD = Beta[3]; + + fit[0].dRo = sqrt(Cov[0][0]); + fit[0].duR = sqrt(Cov[1][1]); + fit[0].dDo = sqrt(Cov[2][2]); + fit[0].duD = sqrt(Cov[3][3]); + + // Sort out the final weight threshold. + + // add up the chi square for the fit + chisq = 0.0; + fit[0].Nfit = 0; + for (i = 0; i < Npts; i++) { + if ((Wx[i] > 0.1 * Sum_Wx / (1.0 * Npts))|| + (Wy[i] > 0.1 * Sum_Wy / (1.0 * Npts))) { + Xf = fit[0].Ro + fit[0].uR*T[i]; + Yf = fit[0].Do + fit[0].uD*T[i]; + chisq += SQ(X[i] - Xf) / SQ(dX[i]); + chisq += SQ(Y[i] - Yf) / SQ(dY[i]); + fit[0].Nfit += 1; + } + // if (VERBOSE) fprintf (stderr, "chisq contrib : %f %f : %f %f : %f %f : %f %f : %f\n", Xf, Yf, X[i] - Xf, Y[i] - Yf, dX[i], dY[i], (X[i] - Xf) / dX[i], (Y[i] - Yf) / dY[i], chisq); + } + // fit[0].Nfit = Npts; + + // the reduced chisq is divided by (Ndof = 2*Npts - 4) + fit[0].chisq = chisq / (2.0*Npts - 4.0); + return (TRUE); +} + + +double weight_cauchy (double x) { + double r = x / 2.385; + return (1.0 / (1.0 + SQ(r))); +} + +// dpsi = (d/dx) (x * weight(x)) +double dpsi_cauchy (double x) { + double r2 = SQ(x / 2.385); + return ((1.0 - r2) / (SQ(1 + r2))); +} + + +// median absolute deviation +// MAD = median(abs(x - median(x))) +double MAD(double *in, int N) { + double *x; + double median = 0.0; + int i; + + ALLOCATE(x,double,N); + for (i = 0; i < N; i++) { + x[i] = in[i]; + } + + dsort(x,N); + + if (N % 2) { + median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]); + } else { + median = x[(int)(0.5*N)]; + } + + for (i = 0; i < N; i++ ) { + x[i] = fabs(x[i] - median); + } + + dsort(x,N); + + if (N % 2) { + median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]); + } else { + median = x[(int)(0.5*N)]; + } + + return(median); +} + + + +int weighted_LS (double *T, double *X, double *WX, double *Y, double *WY, int Npts, + double **A, double **B, int VERBOSE) { + + int i,j; + double Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT; + Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0; + for (i = 0; i < Npts; i++) { + Wx += WX[i]; + Wy += WY[i]; + + Tx += T[i]*WX[i]; + Ty += T[i]*WY[i]; + + Tx2 += SQ(T[i])*WX[i]; + Ty2 += SQ(T[i])*WY[i]; + + Xs += X[i]*WX[i]; + Ys += Y[i]*WY[i]; + + XT += X[i]*T[i]*WX[i]; + YT += Y[i]*T[i]*WY[i]; + } + + // X^T W X + A[0][0] = Wx; + A[0][1] = Tx; + + A[1][0] = Tx; + A[1][1] = Tx2; + + A[2][2] = Wy; + A[2][3] = Ty; + + A[3][2] = Ty; + A[3][3] = Ty2; + + // X^T W Y + B[0][0] = Xs; + B[1][0] = XT; + B[2][0] = Ys; + B[3][0] = YT; + + if (!dgaussjordan ((double **)A, (double **)B, 4, 1)) { + if (VERBOSE) fprintf (stderr, "error in fit\n"); + if (VERBOSE == 2) { + int j; + for (i = 0; i < 4; i++) { + for (j = 0; j < 4; j++) { + fprintf (stderr, "%e ", A[i][j]); + } + fprintf (stderr, " : %e\n", A[i][0]); + } + } + return FALSE; + } + + // A => (X^T W X)^{-1} + // B => beta + + return TRUE; +} Index: /branches/eam_branches/ipp-20150625/Ohana/src/opihi/cmd.basic/list.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/opihi/cmd.basic/list.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/opihi/cmd.basic/list.c (revision 38983) @@ -1,22 +1,34 @@ # include "basic.h" +# include # define D_NLINES 100 static char prompt[] = ">> "; +int SplitByCharToList (int argc, char **argv, int EXCEL_STYLE); +int GlobToList (char *listname, char *Glob, int EXCEL_STYLE); +int CommandToList (char *listname, char *Command, int EXCEL_STYLE); + int list (int argc, char **argv) { int ThisList, depth, i, done, found; char *input, line[1024]; - int N, Nbytes, NBYTES, Nread, status; - int RunCommand; - char *A, *B, *val, *Cmd; - FILE *f; - - Cmd = NULL; - RunCommand = FALSE; + int N; + + char *Command = NULL; if ((N = get_argument (argc, argv, "-x"))) { remove_argument (N, &argc, argv); - Cmd = strcreate (argv[N]); - remove_argument (N, &argc, argv); - RunCommand = TRUE; + Command = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-exec"))) { + remove_argument (N, &argc, argv); + Command = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + + char *Glob = NULL; + if ((N = get_argument (argc, argv, "-glob"))) { + remove_argument (N, &argc, argv); + Glob = strcreate (argv[N]); + remove_argument (N, &argc, argv); } @@ -37,5 +49,5 @@ return (FALSE); } - ListVectorsToList (argv[1]); + ListVectorsToList (argv[1]); // <-- generate the list return TRUE; } @@ -47,6 +59,25 @@ return (FALSE); } - ListBuffersToList (argv[1]); + ListBuffersToList (argv[1]); // <-- generate the list return TRUE; + } + + if (Command) { + if (argc != 2) { + gprint (GP_ERR, "USAGE: list (root) -exec (command)\n"); + gprint (GP_ERR, " OR: list (root) -x (command)\n"); + return (FALSE); + } + int status = CommandToList (argv[1], Command, EXCEL_STYLE); // <-- generate the list + return status; + } + + if (Glob) { + if (argc != 2) { + gprint (GP_ERR, "USAGE: list (root) -glob (fileglob)\n"); + return (FALSE); + } + int status = GlobToList (argv[1], Glob, EXCEL_STYLE); // <-- generate the list + return status; } @@ -74,37 +105,6 @@ return (FALSE); } - - int j, nWords; - char splitter; - char *new, *word, *ptr; - - nWords = 0; - splitter = argv[3][0]; - - for (i = 0; i < argc - 3; i++) { - new = strcreate (argv[i+3]); - for (j = 0; new[j]; j++) { - if (new[j] == splitter) new[j] = ' '; - } - - ptr = new; - while (ptr) { - word = thisword (ptr); - if (!word) break; - - // sprintf (line, "%s:%d", argv[1], nWords); - set_list_varname (line, argv[1], nWords, EXCEL_STYLE); - - set_str_variable (line, word); - FREE (word); - ptr = nextword (ptr); - nWords ++; - } - FREE (new); - } - sprintf (line, "%s:n", argv[1]); - set_int_variable (line, nWords); - - return (TRUE); + int status = SplitByCharToList (argc, argv, EXCEL_STYLE); // <-- generate the list + return status; } @@ -200,4 +200,5 @@ gprint (GP_ERR, "USAGE: list (root) : supply list data, terminate with 'END'\n"); gprint (GP_ERR, "USAGE: list (root) -x (command) : create list from shell output\n"); + gprint (GP_ERR, "USAGE: list (root) -exec (command) : create list from shell output\n"); gprint (GP_ERR, "USAGE: list (root) -vectors : create list from vector names\n"); gprint (GP_ERR, "USAGE: list (root) -buffers : create list from buffer names\n"); @@ -208,59 +209,4 @@ gprint (GP_ERR, "USAGE: list (root) -del (word) : delete the entry by value\n"); return (FALSE); - } - - if (RunCommand) { - - /* val will hold the result */ - NBYTES = 1024; - ALLOCATE (val, char, NBYTES); - - /* need to loop until command produces no more output, - REALLOCATING as needed. */ - f = popen (Cmd, "r"); - done = FALSE; - Nbytes = 0; - while (!done) { - Nread = fread (&val[Nbytes], 1, 1023, f); - if (Nread < 0) { - gprint (GP_ERR, "error reading from command\n"); - done = TRUE; - } - if (Nread > 0) { - Nbytes += Nread; - NBYTES = 1024 + Nbytes; - REALLOCATE (val, char, NBYTES); - } - if (Nread == 0) { - done = TRUE; - } - - } - val[Nbytes] = 0; - status = pclose (f); - free (Cmd); - - if (status) { - gprint (GP_ERR, "warning: exit status of command %d\n", status); - } - - A = B = val; - for (i = 0; B != (char *) NULL;) { - while (isspace (*A) && (*A != 0)) A++; - B = strchr (A, '\n'); - if (B != (char *) NULL) { *B = 0; } - if (*A != 0) { - // sprintf (line, "%s:%d", argv[1], i); - set_list_varname (line, argv[1], i, EXCEL_STYLE); - set_str_variable (line, A); - A = B + 1; - i++; - } - } - free (val); - - sprintf (line, "%s:n", argv[1]); - set_int_variable (line, i); - return (TRUE); } @@ -309,2 +255,123 @@ return (TRUE); } + +int CommandToList (char *listname, char *Command, int EXCEL_STYLE) { + + char line[1024]; + + /* val will hold the result */ + int NBYTES = 1024; + char *val; + ALLOCATE (val, char, NBYTES); + + /* need to loop until command produces no more output, + REALLOCATING as needed. */ + FILE *f = popen (Command, "r"); + + int done = FALSE; + int Nbytes = 0; + while (!done) { + int Nread = fread (&val[Nbytes], 1, 1023, f); + if (Nread < 0) { + gprint (GP_ERR, "error reading from command\n"); + done = TRUE; + } + if (Nread > 0) { + Nbytes += Nread; + NBYTES = 1024 + Nbytes; + REALLOCATE (val, char, NBYTES); + } + if (Nread == 0) { + done = TRUE; + } + } + val[Nbytes] = 0; + int status = pclose (f); + free (Command); + + if (status) { + gprint (GP_ERR, "warning: exit status of command %d\n", status); + } + + int i; + char *A = val; + char *B = val; + for (i = 0; B != (char *) NULL;) { + while (isspace (*A) && (*A != 0)) A++; + B = strchr (A, '\n'); + if (B != (char *) NULL) { *B = 0; } + if (*A != 0) { + set_list_varname (line, listname, i, EXCEL_STYLE); + set_str_variable (line, A); + A = B + 1; + i++; + } + } + free (val); + + sprintf (line, "%s:n", listname); + set_int_variable (line, i); + return (TRUE); +} + +int GlobToList (char *listname, char *Glob, int EXCEL_STYLE) { + + int i; + char line[1024]; + glob_t globList; + + // parse the filename as a glob + globList.gl_offs = 0; + glob (Glob, 0, NULL, &globList); + + // if the glob does not match, save the literal word: + // otherwise save all glob matches + if (globList.gl_pathc == 0) { + sprintf (line, "%s:n", listname); + set_int_variable (line, 0); + return TRUE; + } + + int Nfile = globList.gl_pathc; + for (i = 0; i < Nfile; i++) { + set_list_varname (line, listname, i, EXCEL_STYLE); + set_str_variable (line, globList.gl_pathv[i]); + } + sprintf (line, "%s:n", listname); + set_int_variable (line, Nfile); + return (TRUE); +} + +int SplitByCharToList (int argc, char **argv, int EXCEL_STYLE) { + + int i, j; + char line[1024]; + + int nWords = 0; + char splitter = argv[3][0]; + + for (i = 0; i < argc - 3; i++) { + char *new = strcreate (argv[i+3]); + for (j = 0; new[j]; j++) { + if (new[j] == splitter) new[j] = ' '; + } + + char *ptr = new; + while (ptr) { + char *word = thisword (ptr); + if (!word) break; + + set_list_varname (line, argv[1], nWords, EXCEL_STYLE); + + set_str_variable (line, word); + FREE (word); + ptr = nextword (ptr); + nWords ++; + } + FREE (new); + } + sprintf (line, "%s:n", argv[1]); + set_int_variable (line, nWords); + + return (TRUE); +} Index: /branches/eam_branches/ipp-20150625/Ohana/src/opihi/dvo/avextract.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/opihi/dvo/avextract.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/opihi/dvo/avextract.c (revision 38983) @@ -167,4 +167,5 @@ int needLensobj = dbFieldNeedLensobj (fields, Nfields); int needStarpar = dbFieldNeedStarpar (fields, Nfields, TRUE); + int needGalphot = dbFieldNeedGalphot (fields, Nfields); // grab data from all selected sky regions @@ -191,4 +192,5 @@ catalog.catflags |= needLensobj ? DVO_LOAD_LENSOBJ : DVO_SKIP_LENSOBJ; catalog.catflags |= needStarpar ? DVO_LOAD_STARPAR : DVO_SKIP_STARPAR; + catalog.catflags |= needGalphot ? DVO_LOAD_GALPHOT : DVO_SKIP_GALPHOT; catalog.Nsecfilt = 0; @@ -219,4 +221,7 @@ StarPar *starpar = needStarpar ? &catalog.starpar[m] : NULL; + m = average->galphotOffset; + GalPhot *galphot = needGalphot ? &catalog.galphot[m] : NULL; + m = j*Nsecfilt; SecFilt *secfilt = &catalog.secfilt[m]; @@ -224,5 +229,5 @@ for (n = 0; n < Nfields; n++) { // we are passing in the *first* measure, but average->Nmeasure gives the count - values[n] = dbExtractAverages (average, secfilt, measure, lensobj, starpar, &fields[n]); + values[n] = dbExtractAverages (average, secfilt, measure, lensobj, starpar, galphot, &fields[n]); } Index: /branches/eam_branches/ipp-20150625/Ohana/src/opihi/dvo/avmatch.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/opihi/dvo/avmatch.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/opihi/dvo/avmatch.c (revision 38983) @@ -238,5 +238,5 @@ dbExtractAveragesInitAve (); for (n = 0; n < Nfields; n++) { - values[n] = dbExtractAverages (average, secfilt, measure, lensobj, starpar, &fields[n]); + values[n] = dbExtractAverages (average, secfilt, measure, lensobj, starpar, NULL, &fields[n]); } Index: /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/high_speed_utils.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/high_speed_utils.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/high_speed_utils.c (revision 38983) @@ -141,5 +141,5 @@ off_t n; for (n = 0; n < NfieldsA; n++) { - valuesA[n] = dbExtractAverages (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], &catalog[0].measure[m], NULL, NULL, &fieldsA[n]); + valuesA[n] = dbExtractAverages (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], &catalog[0].measure[m], NULL, NULL, NULL, &fieldsA[n]); } return dbBooleanCond(stackA, NstackA, valuesA); @@ -155,5 +155,5 @@ off_t n; for (n = 0; n < NfieldsB; n++) { - valuesB[n] = dbExtractAverages (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], &catalog[0].measure[m], NULL, NULL, &fieldsB[n]); + valuesB[n] = dbExtractAverages (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], &catalog[0].measure[m], NULL, NULL, NULL, &fieldsB[n]); } return dbBooleanCond(stackB, NstackB, valuesB); Index: /branches/eam_branches/ipp-20150625/Ohana/src/relphot/src/setMrelCatalog.c =================================================================== --- /branches/eam_branches/ipp-20150625/Ohana/src/relphot/src/setMrelCatalog.c (revision 38982) +++ /branches/eam_branches/ipp-20150625/Ohana/src/relphot/src/setMrelCatalog.c (revision 38983) @@ -14,4 +14,7 @@ # define CHECK_VALID_MAG(MAG,D_MAG) (isfinite(MAG) && isfinite(D_MAG) && (MAG > -5.0) && (MAG < 30.0)) # define CHECK_VALID_FLUX(FLUX,D_FLUX) (isfinite(FLUX) && isfinite(D_FLUX)) + +// set A to B if B > A or if A is NAN, as long as B is not NAN +# define MAX_NOT_NAN(A,B) { if (isfinite(B) && (!isfinite(A) || (B > A))) A = B; } static float MagToFlux (float Mag) { @@ -162,6 +165,6 @@ // various steps only relevant to the final pass if (isSetMrelFinal) { - if (isfinite(measure[k].psfQF) && (measure[k].psfQF > results->psfQfMax[Nsec])) results->psfQfMax[Nsec] = measure[k].psfQF; - if (isfinite(measure[k].psfQFperf) && (measure[k].psfQFperf > results->psfQfPerfMax[Nsec])) results->psfQfPerfMax[Nsec] = measure[k].psfQFperf; + MAX_NOT_NAN (results->psfQfMax[Nsec], measure[k].psfQF); + MAX_NOT_NAN (results->psfQfPerfMax[Nsec], measure[k].psfQFperf); // are we a PS1 exposure photcode? (hard-wired photcodes) @@ -219,5 +222,10 @@ if (isSetMrelFinal) { Map = PhotCat (&measure[k], MAG_CLASS_APER); - dMap = MAX (hypot(measure[k].dMap, code->photomErrSys), MIN_ERROR); + + // NOTE: for PV3, apFluxErr is broken (see pmSourcePhotometry.c:245). we are going to + // use PSF flux error instead here: + // dMap = MAX (hypot(measure[k].dMap, code->photomErrSys), MIN_ERROR); + dMap = MAX (hypot(measure[k].dM, code->photomErrSys), MIN_ERROR); + if (CHECK_VALID_MAG(Map, dMap)) { int Nap = results->aperData[Nsec].Nlist; @@ -367,6 +375,6 @@ secfilt[Nsec].psfQfPerfMax = results->psfQfPerfMax[Nsec]; - if (isfinite(secfilt[Nsec].psfQfMax) && (secfilt[Nsec].psfQfMax > psfQfMax)) psfQfMax = secfilt[Nsec].psfQfMax; - if (isfinite(secfilt[Nsec].psfQfPerfMax) && (secfilt[Nsec].psfQfPerfMax > psfQfPerfMax)) psfQfPerfMax = secfilt[Nsec].psfQfPerfMax; + MAX_NOT_NAN (psfQfMax, secfilt[Nsec].psfQfMax); + MAX_NOT_NAN (psfQfPerfMax, secfilt[Nsec].psfQfPerfMax); StatType *apstats = &results->apstats; @@ -640,5 +648,9 @@ secfilt[Nsec].dFkronStk = zpFactor * measure[k].dFluxKron; secfilt[Nsec].FapStk = zpFactor * measure[k].FluxAp; - secfilt[Nsec].dFapStk = zpFactor * measure[k].dFluxAp; + + // NOTE: for PV3, apFluxErr is broken (see pmSourcePhotometry.c:245). we are going to + // use PSF flux error instead here: + // secfilt[Nsec].dFapStk = zpFactor * measure[k].dFluxAp; + secfilt[Nsec].dFapStk = zpFactor * measure[k].dFluxPSF; // Jy to AB mags @@ -742,6 +754,6 @@ if (IS_DIFF_DB) { measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY; - if (measure[k].psfQF > results->psfQfMax[Nsec]) results->psfQfMax[Nsec] = measure[k].psfQF; - if (measure[k].psfQFperf > results->psfQfPerfMax[Nsec]) results->psfQfPerfMax[Nsec] = measure[k].psfQFperf; + MAX_NOT_NAN (results->psfQfMax[Nsec], measure[k].psfQF); + MAX_NOT_NAN (results->psfQfPerfMax[Nsec], measure[k].psfQFperf); } @@ -775,5 +787,10 @@ float Fap = PhotFluxCat (&measure[k], MAG_CLASS_APER); - float dFap = PhotFluxCatErr (&measure[k], MAG_CLASS_APER); + + // NOTE: for PV3, apFluxErr is broken (see pmSourcePhotometry.c:245). we are going to + // use PSF flux error instead here: + // float dFap = PhotFluxCatErr (&measure[k], MAG_CLASS_APER); + float dFap = PhotFluxCatErr (&measure[k], MAG_CLASS_PSF); + if (CHECK_VALID_FLUX(Fap, dFap)) { dFap = MAX (hypot(dFap, code->photomErrSys*Fap), MIN_ERROR*Fap); // bump up the error by a systematic floor @@ -814,6 +831,6 @@ secfilt[Nsec].psfQfPerfMax = results->psfQfPerfMax[Nsec]; - if (isfinite(secfilt[Nsec].psfQfMax) && (secfilt[Nsec].psfQfMax > psfQfMax)) psfQfMax = secfilt[Nsec].psfQfMax; - if (isfinite(secfilt[Nsec].psfQfPerfMax) && (secfilt[Nsec].psfQfPerfMax > psfQfPerfMax)) psfQfPerfMax = secfilt[Nsec].psfQfPerfMax; + MAX_NOT_NAN (psfQfMax, secfilt[Nsec].psfQfMax); + MAX_NOT_NAN (psfQfPerfMax, secfilt[Nsec].psfQfPerfMax); }