Index: trunk/ppSim/src/ppSim.h =================================================================== --- trunk/ppSim/src/ppSim.h (revision 17032) +++ trunk/ppSim/src/ppSim.h (revision 17557) @@ -133,4 +133,9 @@ ); +float ppSimStarSkyNoise (float skySigma, float seeingSigma); +float ppSimStarPeakToFlux (float peak, float seeingSigma); +float ppSimStarFluxToPeak (float flux, float seeingSigma); +float ppSimFluxToMag (float flux, float zp); +float ppSimMagToFlux (float mag, float zp); #endif Index: trunk/ppSim/src/ppSimArguments.c =================================================================== --- trunk/ppSim/src/ppSimArguments.c (revision 17032) +++ trunk/ppSim/src/ppSimArguments.c (revision 17557) @@ -18,31 +18,26 @@ } -// Get a value from the command-line arguments or the recipe, and add it to the target -bool valueArgRecipe(pmConfig *config, // Configuration - psMetadata *arguments, // Command-line arguments - const char *argName, // Argument name in the command-line arguments - const psMetadata *recipe, // Recipe +// Get a value from the command-line arguments and add it to recipe options +bool valueArgRecipe(psMetadata *options, // Target to which to add value const char *recipeName, // Name for value in the recipe - psMetadata *target // Target to which to add value + psMetadata *arguments, // Command-line arguments + const char *argName // Argument name in the command-line arguments ) { - float value = psMetadataLookupF32(NULL, arguments, argName); // Value of interest - if (isnan(value)) { - bool mdok; // Status of MD lookup - value = psMetadataLookupF32(&mdok, recipe, recipeName); - if (!mdok) { - psErrorStackPrint(stderr, "Unable to find %s in recipe %s", recipeName, PPSIM_RECIPE); - psFree((psPtr)arguments); - psFree(config); - exit(PS_EXIT_CONFIG_ERROR); - } - } - return psMetadataAddF32(target, PS_LIST_TAIL, recipeName, 0, NULL, value); + bool status; // Status of MD lookup + float value = psMetadataLookupF32(&status, arguments, argName); // Value of interest + if (isnan(value)) return true; + status = psMetadataAddF32(options, PS_LIST_TAIL, recipeName, 0, NULL, value); + return status; } +// this function supplements the RECIPE:OPTIONS folder with command-line options void ppSimArguments(int argc, char *argv[], pmConfig *config) { + bool mdok; // Status of MD lookup + assert(config); + // save the following command-line options in the arguments structure psMetadata *arguments = psMetadataAlloc(); // Command-line arguments psMetadataAddStr(arguments, PS_LIST_TAIL, "-format", 0, "Camera format name", NULL); @@ -68,4 +63,5 @@ psMetadataAddF32(arguments, PS_LIST_TAIL, "-starsdensity", 0, "Density of fake stars at magnitude", NAN); psMetadataAddStr(arguments, PS_LIST_TAIL, "-psfclass", 0, "Type of PSF model", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-galmodel", 0, "Type of Galaxy model", NULL); psMetadataAddS32(arguments, PS_LIST_TAIL, "-bin", 0, "Binning in x and y", 1); @@ -92,4 +88,5 @@ } + // apply an alternate camera format psString formatName = psMetadataLookupStr(NULL, arguments, "-format"); // Name of format if (formatName) { @@ -114,4 +111,6 @@ } + // specify the type of simulated image to produce + // XXX this should not be required if we supplied INPUT const char *typeStr = psMetadataLookupStr(NULL, arguments, "-type"); // Exposure type if (!typeStr) { @@ -138,4 +137,14 @@ psMetadataAddStr(config->arguments, PS_LIST_TAIL, "FILTER", 0, "Filter name", filter); + // save the following additional recipe values based on command-line options + // these options override the PPSIM recipe values loaded from recipe files + psMetadata *options = pmConfigRecipeOptions (config, PPSIM_RECIPE); + if (!options) { + psErrorStackPrint(stderr, "Unable to find recipe options for %s", PPSIM_RECIPE); + psFree(arguments); + psFree(config); + exit(PS_EXIT_CONFIG_ERROR); + } + float expTime; if (type == PPSIM_TYPE_BIAS) { @@ -148,32 +157,25 @@ } } - psMetadataAddF32(config->arguments, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (s)", expTime); - - psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe for ppSim - if (!recipe) { - psErrorStackPrint(stderr, "Unable to find recipe %s", PPSIM_RECIPE); - psFree(arguments); - psFree(config); - exit(PS_EXIT_CONFIG_ERROR); - } - - valueArgRecipe(config, arguments, "-biaslevel", recipe, "BIAS.LEVEL", config->arguments); - valueArgRecipe(config, arguments, "-biasrange", recipe, "BIAS.RANGE", config->arguments); - valueArgRecipe(config, arguments, "-darkrate", recipe, "DARK.RATE", config->arguments); - valueArgRecipe(config, arguments, "-flatsigma", recipe, "FLAT.SIGMA", config->arguments); - valueArgRecipe(config, arguments, "-flatrate", recipe, "FLAT.RATE", config->arguments); - valueArgRecipe(config, arguments, "-shuttertime", recipe, "SHUTTER.TIME", config->arguments); - valueArgRecipe(config, arguments, "-skyrate", recipe, "SKY.RATE", config->arguments); - valueArgRecipe(config, arguments, "-starslum", recipe, "STARS.LUM", config->arguments); - valueArgRecipe(config, arguments, "-starsmag", recipe, "STARS.MAG", config->arguments); - valueArgRecipe(config, arguments, "-starsdensity", recipe, "STARS.DENSITY", config->arguments); - - bool mdok; // Status of MD lookup + psMetadataAddF32(options, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (s)", expTime); + + // these values all get moved from arguments to RECIPE:OPTIONS + valueArgRecipe(options, "BIAS.LEVEL", arguments, "-biaslevel"); + valueArgRecipe(options, "BIAS.RANGE", arguments, "-biasrange"); + valueArgRecipe(options, "DARK.RATE", arguments, "-darkrate"); + valueArgRecipe(options, "FLAT.SIGMA", arguments, "-flatsigma"); + valueArgRecipe(options, "FLAT.RATE", arguments, "-flatrate"); + valueArgRecipe(options, "SHUTTER.TIME", arguments, "-shuttertime"); + valueArgRecipe(options, "SKY.RATE", arguments, "-skyrate"); + valueArgRecipe(options, "STARS.LUM", arguments, "-starslum"); + valueArgRecipe(options, "STARS.MAG", arguments, "-starsmag"); + valueArgRecipe(options, "STARS.DENSITY", arguments, "-starsdensity"); + + psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe + int biasOrder = psMetadataLookupS32(&mdok, recipe, "BIAS.ORDER"); // Overscan polynomial order if (!mdok) { psWarning("BIAS.ORDER(S32) is not set in the recipe %s --- assuming %d", PPSIM_RECIPE, biasOrder); } - psMetadataAddS32(config->arguments, PS_LIST_TAIL, "BIAS.ORDER", 0, - "Overscan polynomial order", biasOrder); + psMetadataAddS32(options, PS_LIST_TAIL, "BIAS.ORDER", 0, "Overscan polynomial order", biasOrder); int binning = psMetadataLookupS32(NULL, arguments, "-bin"); // Binning in x and y @@ -190,6 +192,6 @@ float ra0 = psMetadataLookupF32(NULL, arguments, "-ra"); // Right Ascension of boresight float dec0 = psMetadataLookupF32(NULL, arguments, "-dec"); // Declination of boresight - float pa = psMetadataLookupF32(NULL, arguments, "-pa"); // Position angle - float seeing = psMetadataLookupF32(NULL, arguments, "-seeing"); // Zero point + float pa = psMetadataLookupF32(NULL, arguments, "-pa"); // Position angle + float seeing = psMetadataLookupF32(NULL, arguments, "-seeing"); // seeing (FWHM in arcsec) // XXX scale and zp should be supplied by the config file (allow override, but this is camera-dependent) @@ -199,5 +201,5 @@ } - float zp = psMetadataLookupF32(NULL, arguments, "-zp"); // Zero point + float zp = psMetadataLookupF32(NULL, arguments, "-zp"); // Zero point if (isnan(zp)) { // use the filter to get the zeropoint from the recipe @@ -234,15 +236,20 @@ } - psMetadataAddF32(config->arguments, PS_LIST_TAIL, "RA", 0, "Boresight RA (radians)", ra0 * M_PI / 180.0); - psMetadataAddF32(config->arguments, PS_LIST_TAIL, "DEC", 0, "Boresight Declination (radians)", dec0 * M_PI / 180.0); - psMetadataAddF32(config->arguments, PS_LIST_TAIL, "PA", 0, "Boresight position angle (radians)",pa * M_PI / 180.0); - psMetadataAddF32(config->arguments, PS_LIST_TAIL, "SEEING", 0, "Seeing sigma (pix)", seeing / 2.0 / sqrt(2.0 * log(2.0)) / scale); - - psMetadataAddF32(config->arguments, PS_LIST_TAIL, "SCALE", 0, "Plate scale (arcsec/pix)", scale); - psMetadataAddF32(config->arguments, PS_LIST_TAIL, "ZEROPOINT", 0, "Photometric zeropoint", zp); - psMetadataAddF32(config->arguments, PS_LIST_TAIL, "SKY.MAGS", 0, "sky surface brightness", skymags); - - const char *psfClass = psMetadataLookupStr(NULL, arguments, "-psfclass"); // Filter name + psMetadataAddF32(options, PS_LIST_TAIL, "RA", 0, "Boresight RA (radians)", ra0 * M_PI / 180.0); + psMetadataAddF32(options, PS_LIST_TAIL, "DEC", 0, "Boresight Declination (radians)", dec0 * M_PI / 180.0); + psMetadataAddF32(options, PS_LIST_TAIL, "PA", 0, "Boresight position angle (radians)",pa * M_PI / 180.0); + + // the user supplies FWHM in arcsec; here we convert to Sigma in pixels + psMetadataAddF32(options, PS_LIST_TAIL, "SEEING", 0, "Seeing SIGMA (pixels)", seeing / 2.0 / sqrt(2.0 * log(2.0)) / scale); + + psMetadataAddF32(options, PS_LIST_TAIL, "SCALE", 0, "Plate scale (arcsec/pix)", scale); + psMetadataAddF32(options, PS_LIST_TAIL, "ZEROPOINT", 0, "Photometric zeropoint", zp); + psMetadataAddF32(options, PS_LIST_TAIL, "SKY.MAGS", 0, "sky surface brightness", skymags); + + const char *psfClass = psMetadataLookupStr(NULL, arguments, "-psfclass"); // PSF model class psMetadataAddStr(config->arguments, PS_LIST_TAIL, "PSF.MODEL", 0, "PSF model class", psfClass); + + const char *galModel = psMetadataLookupStr(NULL, arguments, "-galmodel"); // Galaxy model name + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "GALAXY.MODEL", 0, "Galaxy model", galModel); } @@ -250,2 +257,17 @@ return; } + +/* the following elements come from the config->arguments: + + PSPHOT.PSF + INPUT + TYPE + FILTER + BIAS.ORDER + BINNING + OUTPUT + PSF.MODEL + GALAXY.MODEL + + all othr values should come from the recipe +*/ Index: trunk/ppSim/src/ppSimCreate.c =================================================================== --- trunk/ppSim/src/ppSimCreate.c (revision 17032) +++ trunk/ppSim/src/ppSimCreate.c (revision 17557) @@ -18,4 +18,10 @@ pmFPAfile *input = pmFPAfileDefineFromArgs (&status, config, "PPIMAGE.INPUT", "INPUT"); if (!input) { + // if we have not specified the camera already, we need to interpolate the recipes associated with this camera, and read other command-line recipes + if (!pmConfigReadRecipes(config, PM_RECIPE_SOURCE_CL)) { + psError(PS_ERR_IO, false, "Error merging recipes from camera config for %s", config->cameraName); + return NULL; + } + simImage = true; fpa = pmFPAConstruct(config->camera); // FPA to contain the observation @@ -24,4 +30,5 @@ return NULL; } + } else { simImage = false; @@ -33,4 +40,5 @@ } + // define the output image file pmFPAfile *file = pmFPAfileDefineOutput(config, fpa, OUTPUT_FILE); if (!file) { Index: trunk/ppSim/src/ppSimInsertGalaxies.c =================================================================== --- trunk/ppSim/src/ppSimInsertGalaxies.c (revision 17032) +++ trunk/ppSim/src/ppSimInsertGalaxies.c (revision 17557) @@ -22,6 +22,7 @@ psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe - float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time - float darkRate = psMetadataLookupF32(NULL, config->arguments, "DARK.RATE"); // Dark rate + float expTime = psMetadataLookupF32(NULL, recipe, "EXPTIME"); // Exposure time + float darkRate = psMetadataLookupF32(NULL, recipe, "DARK.RATE"); // Dark rate + float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE");// CCD read noise, e if (isnan(readnoise)) { @@ -34,10 +35,10 @@ } - float skyRate = psMetadataLookupF32(NULL, config->arguments, "SKY.RATE"); // Sky rate + float skyRate = psMetadataLookupF32(NULL, recipe, "SKY.RATE"); // Sky rate if (isnan(skyRate)) { - float zp = psMetadataLookupF32(&mdok, config->arguments, "ZEROPOINT"); assert (mdok); - float scale = psMetadataLookupF32(&mdok, config->arguments, "SCALE"); assert (mdok); - float skyMags = psMetadataLookupF32(&mdok, config->arguments, "SKY.MAGS"); assert (mdok); - skyRate = scale * scale * pow (10.0, -0.4*(skyMags - zp)); + float zp = psMetadataLookupF32(&mdok, recipe, "ZEROPOINT"); assert (mdok); + float scale = psMetadataLookupF32(&mdok, recipe, "SCALE"); assert (mdok); + float skyMags = psMetadataLookupF32(&mdok, recipe, "SKY.MAGS"); assert (mdok); + skyRate = scale * scale * ppSimMagToFlux (skyMags, zp); } @@ -58,5 +59,5 @@ // determine the galaxy model - char *modelName = psMetadataLookupStr(&mdok, config->arguments, "GALAXY.MODEL"); // Seeing sigma (pix) + char *modelName = psMetadataLookupStr(&mdok, config->arguments, "GALAXY.MODEL"); // galaxy model name if (modelName == NULL) { modelName = defaultModel; Index: trunk/ppSim/src/ppSimInsertStars.c =================================================================== --- trunk/ppSim/src/ppSimInsertStars.c (revision 17032) +++ trunk/ppSim/src/ppSimInsertStars.c (revision 17557) @@ -27,6 +27,7 @@ psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe - float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time - float darkRate = psMetadataLookupF32(NULL, config->arguments, "DARK.RATE"); // Dark rate + float expTime = psMetadataLookupF32(NULL, recipe, "EXPTIME"); // Exposure time + float darkRate = psMetadataLookupF32(NULL, recipe, "DARK.RATE"); // Dark rate + float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE");// CCD read noise, e if (isnan(readnoise)) { @@ -39,10 +40,10 @@ } - float skyRate = psMetadataLookupF32(NULL, config->arguments, "SKY.RATE"); // Sky rate + float skyRate = psMetadataLookupF32(NULL, recipe, "SKY.RATE"); // Sky rate if (isnan(skyRate)) { - float zp = psMetadataLookupF32(&mdok, config->arguments, "ZEROPOINT"); assert (mdok); - float scale = psMetadataLookupF32(&mdok, config->arguments, "SCALE"); assert (mdok); - float skyMags = psMetadataLookupF32(&mdok, config->arguments, "SKY.MAGS"); assert (mdok); - skyRate = scale * scale * pow (10.0, -0.4*(skyMags - zp)); + float zp = psMetadataLookupF32(&mdok, recipe, "ZEROPOINT"); assert (mdok); + float scale = psMetadataLookupF32(&mdok, recipe, "SCALE"); assert (mdok); + float skyMags = psMetadataLookupF32(&mdok, recipe, "SKY.MAGS"); assert (mdok); + skyRate = scale * scale * ppSimMagToFlux (skyMags, zp); } Index: trunk/ppSim/src/ppSimLoadStars.c =================================================================== --- trunk/ppSim/src/ppSimLoadStars.c (revision 17032) +++ trunk/ppSim/src/ppSimLoadStars.c (revision 17557) @@ -18,11 +18,12 @@ } - float zp = psMetadataLookupF32(NULL, config->arguments, "ZEROPOINT"); // Photometric zero point - float ra0 = psMetadataLookupF32(NULL, config->arguments, "RA"); // Boresight RA (radians) - float dec0 = psMetadataLookupF32(NULL, config->arguments, "DEC"); // Boresight Dec (radians) - float pa = psMetadataLookupF32(NULL, config->arguments, "PA"); // Position angle (radians) - float seeing = psMetadataLookupF32(NULL, config->arguments, "SEEING"); // Seeing sigma (pix) - float scale = psMetadataLookupF32(NULL, config->arguments, "SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel) - float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time + // XXX push these into the recipe in ppSimArguments() + float zp = psMetadataLookupF32(NULL, recipe, "ZEROPOINT"); // Photometric zero point + float ra0 = psMetadataLookupF32(NULL, recipe, "RA"); // Boresight RA (radians) + float dec0 = psMetadataLookupF32(NULL, recipe, "DEC"); // Boresight Dec (radians) + float pa = psMetadataLookupF32(NULL, recipe, "PA"); // Position angle (radians) + float seeing = psMetadataLookupF32(NULL, recipe, "SEEING"); // Seeing SIGMA (pixels) + float scale = psMetadataLookupF32(NULL, recipe, "SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel) + float expTime = psMetadataLookupF32(NULL, recipe, "EXPTIME"); // Exposure time (sec) // Size of FPA @@ -68,6 +69,7 @@ // Convert magnitude to peak flux - star->flux = powf(10.0, -0.4 * (star->mag - zp)) * expTime; - star->peak = star->flux / (2.0*M_PI * PS_SQR(seeing)); + star->flux = ppSimMagToFlux (star->mag, zp) * expTime; + star->peak = ppSimStarFluxToPeak (star->flux, seeing); + stars->data[oldSize + i] = star; @@ -75,8 +77,21 @@ } stars->n = oldSize + refStars->n; + + pmLumFunc *lumfunc = psastroLuminosityFunction (refStars); psFree(refStars); psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "FPA.RA", PS_META_REPLACE, "Right ascension", ra0); psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "FPA.DEC", PS_META_REPLACE, "Declination", dec0); + + if (lumfunc) { + psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.MAG.MIN", PS_META_REPLACE, "min valid magnitude", lumfunc->mMin); + psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.MAG.MAX", PS_META_REPLACE, "max valid magnitude", lumfunc->mMax); + psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.LF.SLOPE", PS_META_REPLACE, "log-mag histogram slope", lumfunc->slope); + psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.LF.OFFSET", PS_META_REPLACE, "log-mag histogram offset", lumfunc->offset); + psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.MAG.PEAK", PS_META_REPLACE, "magnitude of peak bin", lumfunc->mPeak); + psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.NUM.PEAK", PS_META_REPLACE, "number of stars in peak bin", lumfunc->nPeak); + psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "STARS.REAL.SUM.PEAK", PS_META_REPLACE, "sum of stars up to peak bin", lumfunc->sPeak); + psFree (lumfunc); + } return stars; Index: trunk/ppSim/src/ppSimMakeBias.c =================================================================== --- trunk/ppSim/src/ppSimMakeBias.c (revision 17032) +++ trunk/ppSim/src/ppSimMakeBias.c (revision 17557) @@ -9,7 +9,7 @@ psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe - float biasLevel = psMetadataLookupF32(NULL, config->arguments, "BIAS.LEVEL"); // Bias level - float biasRange = psMetadataLookupF32(NULL, config->arguments, "BIAS.RANGE"); // Bias range - int biasOrder = psMetadataLookupS32(NULL, config->arguments, "BIAS.ORDER"); // Bias order + float biasLevel = psMetadataLookupF32(NULL, recipe, "BIAS.LEVEL"); // Bias level + float biasRange = psMetadataLookupF32(NULL, recipe, "BIAS.RANGE"); // Bias range + int biasOrder = psMetadataLookupS32(NULL, recipe, "BIAS.ORDER"); // Bias order float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE");// CCD read noise, e Index: trunk/ppSim/src/ppSimMakeDark.c =================================================================== --- trunk/ppSim/src/ppSimMakeDark.c (revision 17032) +++ trunk/ppSim/src/ppSimMakeDark.c (revision 17557) @@ -4,9 +4,13 @@ bool ppSimMakeDark (pmReadout *readout, pmConfig *config) { + bool mdok; + psImage *signal = readout->image; psImage *variance = readout->weight; - float darkRate = psMetadataLookupF32(NULL, config->arguments, "DARK.RATE"); // Dark rate - float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time + psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe + + float darkRate = psMetadataLookupF32(NULL, recipe, "DARK.RATE"); // Dark rate + float expTime = psMetadataLookupF32(NULL, recipe, "EXPTIME"); // Exposure time for (int y = 0; y < signal->numRows; y++) { Index: trunk/ppSim/src/ppSimMakeGalaxies.c =================================================================== --- trunk/ppSim/src/ppSimMakeGalaxies.c (revision 17032) +++ trunk/ppSim/src/ppSimMakeGalaxies.c (revision 17557) @@ -11,11 +11,11 @@ if (!galaxyFake) return NULL; - float galaxyLum = psMetadataLookupF32(&mdok, recipe, "GALAXY.LUM"); // Galaxy luminosity func slope - float galaxyMag = psMetadataLookupF32(&mdok, recipe, "GALAXY.MAG"); // Galaxy brightest magnitude - float galaxyDensity = psMetadataLookupF32(&mdok, recipe, "GALAXY.DENSITY"); // Density of fakes + float galaxyLum = psMetadataLookupF32(&mdok, recipe, "GALAXY.LUM"); // Galaxy luminosity func slope + float galaxyMag = psMetadataLookupF32(&mdok, recipe, "GALAXY.MAG"); // Galaxy brightest magnitude + float galaxyDensity = psMetadataLookupF32(&mdok, recipe, "GALAXY.DENSITY"); // Density of fakes - bool galaxyGrid = psMetadataLookupBool(&mdok, recipe, "GALAXY.GRID"); // Density of fakes - int galaxyGridDX = psMetadataLookupS32(&mdok, recipe, "GALAXY.GRID.DX"); // Density of fakes - int galaxyGridDY = psMetadataLookupS32(&mdok, recipe, "GALAXY.GRID.DY"); // Density of fakes + bool galaxyGrid = psMetadataLookupBool(&mdok, recipe, "GALAXY.GRID"); // Density of fakes + int galaxyGridDX = psMetadataLookupS32(&mdok, recipe, "GALAXY.GRID.DX"); // Density of fakes + int galaxyGridDY = psMetadataLookupS32(&mdok, recipe, "GALAXY.GRID.DY"); // Density of fakes float galaxyRmajorMax = psMetadataLookupF32(&mdok, recipe, "GALAXY.RMAJOR.MAX"); // Density of fakes @@ -25,22 +25,19 @@ float galaxyARatioMin = psMetadataLookupF32(&mdok, recipe, "GALAXY.ARATIO.MIN"); // Density of fakes - float galaxyThetaMax = psMetadataLookupF32(&mdok, recipe, "GALAXY.THETA.MAX"); // Density of fakes - float galaxyThetaMin = psMetadataLookupF32(&mdok, recipe, "GALAXY.THETA.MIN"); // Density of fakes + float galaxyThetaMax = psMetadataLookupF32(&mdok, recipe, "GALAXY.THETA.MAX"); // Density of fakes + float galaxyThetaMin = psMetadataLookupF32(&mdok, recipe, "GALAXY.THETA.MIN"); // Density of fakes float galaxyIndexMin = psMetadataLookupF32(&mdok, recipe, "GALAXY.INDEX.MIN"); // Density of fakes float galaxyIndexMax = psMetadataLookupF32(&mdok, recipe, "GALAXY.INDEX.MAX"); // Density of fakes - // XXX push these into the recipe... - float darkRate = psMetadataLookupF32(&mdok, config->arguments, "DARK.RATE"); // Dark rate - float expTime = psMetadataLookupF32(&mdok, config->arguments, "EXPTIME"); // Exposure time - - float zp = psMetadataLookupF32(&mdok, config->arguments, "ZEROPOINT"); // Photometric zero point - float seeing = psMetadataLookupF32(&mdok, config->arguments, "SEEING"); // Seeing sigma (pix) - float scale = psMetadataLookupF32(&mdok, config->arguments, "SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel) - - float skyRate = psMetadataLookupF32(&mdok, config->arguments, "SKY.RATE"); // Sky rate + float darkRate = psMetadataLookupF32(&mdok, recipe, "DARK.RATE"); // Dark rate + float expTime = psMetadataLookupF32(&mdok, recipe, "EXPTIME"); // Exposure time + float zp = psMetadataLookupF32(&mdok, recipe, "ZEROPOINT"); // Photometric zero point + float seeing = psMetadataLookupF32(&mdok, recipe, "SEEING"); // Seeing sigma (pix) + float scale = psMetadataLookupF32(&mdok, recipe, "SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel) + float skyRate = psMetadataLookupF32(&mdok, recipe, "SKY.RATE"); // Sky rate if (isnan(skyRate)) { - float skyMags = psMetadataLookupF32(&mdok, config->arguments, "SKY.MAGS"); assert (mdok); - skyRate = scale * scale * pow (10.0, -0.4*(skyMags - zp)); + float skyMags = psMetadataLookupF32(&mdok, recipe, "SKY.MAGS"); assert (mdok); + skyRate = scale * scale * ppSimMagToFlux (skyMags, zp); } Index: trunk/ppSim/src/ppSimMakeSky.c =================================================================== --- trunk/ppSim/src/ppSimMakeSky.c (revision 17032) +++ trunk/ppSim/src/ppSimMakeSky.c (revision 17557) @@ -13,15 +13,16 @@ pmFPA *fpa = chip->parent; - float expTime = psMetadataLookupF32(&status, config->arguments, "EXPTIME"); // Exposure time - float flatSigma = psMetadataLookupF32(&status, config->arguments, "FLAT.SIGMA"); // Flat-field coefficient - float flatRate = psMetadataLookupF32(&status, config->arguments, "FLAT.RATE"); // Flat-field rate - float shutterTime = psMetadataLookupF32(&status, config->arguments, "SHUTTER.TIME"); // Shutter time + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe - float skyRate = psMetadataLookupF32(&status, config->arguments, "SKY.RATE"); // Sky rate + float expTime = psMetadataLookupF32(&status, recipe, "EXPTIME"); // Exposure time + float flatSigma = psMetadataLookupF32(&status, recipe, "FLAT.SIGMA"); // Flat-field coefficient + float flatRate = psMetadataLookupF32(&status, recipe, "FLAT.RATE"); // Flat-field rate + float shutterTime = psMetadataLookupF32(&status, recipe, "SHUTTER.TIME"); // Shutter time + float skyRate = psMetadataLookupF32(&status, recipe, "SKY.RATE"); // Sky rate if (isnan(skyRate)) { - float zp = psMetadataLookupF32(&status, config->arguments, "ZEROPOINT"); assert (status); - float scale = psMetadataLookupF32(&status, config->arguments, "SCALE"); assert (status); - float skyMags = psMetadataLookupF32(&status, config->arguments, "SKY.MAGS"); assert (status); - skyRate = scale * scale * pow (10.0, -0.4*(skyMags - zp)); + float zp = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); assert (status); + float scale = psMetadataLookupF32(&status, recipe, "SCALE"); assert (status); + float skyMags = psMetadataLookupF32(&status, recipe, "SKY.MAGS"); assert (status); + skyRate = scale * scale * ppSimMagToFlux (skyMags, zp); } Index: trunk/ppSim/src/ppSimMakeStars.c =================================================================== --- trunk/ppSim/src/ppSimMakeStars.c (revision 17032) +++ trunk/ppSim/src/ppSimMakeStars.c (revision 17557) @@ -1,35 +1,33 @@ # include "ppSim.h" - -#define FAINT_FUDGE_FACTOR 0.4 // Fraction of the noise in a (boxcar) PSF for the faint limit - bool ppSimMakeStars(psArray *stars, pmFPA *fpa, pmConfig *config, const psRandom *rng) { - bool mdok; + bool status; assert (stars); - psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe - bool starsFake = psMetadataLookupBool(&mdok, recipe, "STARS.FAKE"); // Density of fakes + bool starsFake = psMetadataLookupBool(&status, recipe, "STARS.FAKE"); // Density of fakes if (!starsFake) return true; - float starsLum = psMetadataLookupF32(NULL, config->arguments, "STARS.LUM"); // Star luminosity func slope - float starsMag = psMetadataLookupF32(NULL, config->arguments, "STARS.MAG"); // Star brightest magnitude - float starsDensity = psMetadataLookupF32(NULL, config->arguments, "STARS.DENSITY"); // Density of fakes + bool starsReal = psMetadataLookupBool(&status, recipe, "STARS.REAL"); // were real stars generated? + bool matchDensity = psMetadataLookupBool(&status, recipe, "MATCH.DENSITY"); // match real star density? - float darkRate = psMetadataLookupF32(NULL, config->arguments, "DARK.RATE"); // Dark rate - float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time + float starsAlpha = psMetadataLookupF32(&status, recipe, "STARS.LUM"); // Star luminosity func slope + float starsDensity = psMetadataLookupF32(&status, recipe, "STARS.DENSITY"); // Density of fakes + float brightMag = psMetadataLookupF32(&status, recipe, "STARS.MAG"); // Star brightest magnitude + float nSigmaLim = psMetadataLookupF32(&status, recipe, "STARS.SIGMA.LIM"); // significance of faintest stars - float zp = psMetadataLookupF32(NULL, config->arguments, "ZEROPOINT"); // Photometric zero point - float seeing = psMetadataLookupF32(NULL, config->arguments, "SEEING"); // Seeing sigma (pix) - float scale = psMetadataLookupF32(NULL, config->arguments, "SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel) + float darkRate = psMetadataLookupF32(&status, recipe, "DARK.RATE"); // Dark rate + float expTime = psMetadataLookupF32(&status, recipe, "EXPTIME"); // Exposure time + float zp = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); // Photometric zero point + float seeing = psMetadataLookupF32(&status, recipe, "SEEING"); // Seeing SIGMA (pixels) + float scale = psMetadataLookupF32(&status, recipe, "SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel) - float skyRate = psMetadataLookupF32(&mdok, config->arguments, "SKY.RATE"); // Sky rate + float skyRate = psMetadataLookupF32(&status, recipe, "SKY.RATE"); // Sky rate if (isnan(skyRate)) { - float skyMags = psMetadataLookupF32(&mdok, config->arguments, "SKY.MAGS"); assert (mdok); - skyRate = scale * scale * pow (10.0, -0.4*(skyMags - zp)); + float skyMags = psMetadataLookupF32(&status, recipe, "SKY.MAGS"); assert (status); + skyRate = scale * scale * ppSimMagToFlux (skyMags, zp); } - - if (starsDensity <= 0) return true; // Size of FPA @@ -40,8 +38,23 @@ psFree(bounds); + // choose reference magnitude & density to set normalization + float refMag = 0; + float refSum = 0; + if (starsReal && matchDensity) { + refMag = psMetadataLookupF32(&status, fpa->concepts, "STARS.REAL.MAG.PEAK"); // Star brightest magnitude + refSum = psMetadataLookupF32(&status, fpa->concepts, "STARS.REAL.SUM.PEAK"); // Star brightest magnitude + assert (status); + } else { + refMag = brightMag; + refSum = starsDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI); + } + psTrace("ppSim", 6, "refMag: %f, refSum: %f\n", refMag, refSum); + + if (refSum <= 0) return true; + // Grabbing read noise from the recipe rather than the cell, which is a potential danger, but it // shouldn't be too bad. - float readnoise = psMetadataLookupF32(&mdok, recipe, "READNOISE"); // Default read noise - if (!mdok) { + float readnoise = psMetadataLookupF32(&status, recipe, "READNOISE"); // Default read noise + if (!status) { psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe."); psFree(bounds); @@ -50,10 +63,18 @@ // Faintest and brightest (integrated) fluxes (actually fluence, since integrated) for random stars - float faint = FAINT_FUDGE_FACTOR * sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime) * - (2.0 * sqrt(2.0 * log(2.0)) * seeing); // Faint limit is related to the noise in a (boxcar) PSF - float bright = powf(10.0, -0.4 * (starsMag - zp)) * expTime; // Bright limit is specified by user as mag - psTrace("ppSim", 6, "Faint limit: %f\n", faint); - psTrace("ppSim", 6, "Bright limit: %f\n", bright); - if (bright < faint) { + + // faint limit is set by detection limit, in turn set by sky noise + float skySigma = sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime); + float skyNoise = ppSimStarSkyNoise (skySigma, seeing); + + // total flux of faintest star: + float faintCounts = nSigmaLim * skyNoise; + float faintMag = ppSimFluxToMag ((faintCounts / expTime), zp); + + float brightCounts = ppSimMagToFlux (brightMag, zp) * expTime; // Bright limit is specified by user as mag + + psTrace("ppSim", 6, "Faint limit: %f counts = %f mags\n", faintCounts, faintMag); + psTrace("ppSim", 6, "Bright limit: %f counts = %f mags\n", brightCounts, brightMag); + if (brightCounts < faintCounts) { psLogMsg("ppSim", PS_LOG_INFO, "Image noise is above brightest random star --- no random stars added."); @@ -61,18 +82,23 @@ } - // Normalisation, set by the specified stellar density at the specified bright magnitude - float norm = starsDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI) / - powf(bright, starsLum); + // given log_10 (dN / dmag) = alpha mag + beta + // or dN / dmag = No 10^(alpha mag), then: + // N(m < Mo) = alpha * No * log(10.0) * 10^(alpha*Mo) + + // XXX this needs to handle the case of a flat distribution for efficiency testing + + // Normalization, set by the specified stellar density at the specified bright magnitude + float norm = refSum / (starsAlpha * logf(10.0) * powf(10.0, (starsAlpha * refMag))); + float normScale = norm * starsAlpha * logf(10.0); // Total number of stars down to the faint flux end - long num = expf(logf(norm) + starsLum * logf(faint)) + 0.5; + long nTotal = normScale * powf (10.0, (starsAlpha * faintMag)); - psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld stars down to %f mag\n", - num, -2.5 * log10(faint / expTime) + zp); + psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld stars down to %f mag\n", nTotal, faintMag); long oldSize = stars->n; - psArrayRealloc (stars, stars->n + num); + psArrayRealloc (stars, stars->n + nTotal); - for (long i = 0; i < num; i++) { + for (long i = 0; i < nTotal; i++) { ppSimStar *star = ppSimStarAlloc (); @@ -81,10 +107,13 @@ star->y = psRandomUniform(rng) * ySize; // y position - star->flux = expf((logf(i + 1) - logf(norm)) / starsLum); // Peak flux - star->peak = star->flux / (2.0*M_PI*PS_SQR(seeing)); + // XXX this needs to respect the bright limit + float starMag = log10((i + 1.0) / normScale) / starsAlpha; + + star->flux = ppSimMagToFlux (starMag, zp) * expTime; + star->peak = ppSimStarFluxToPeak (star->flux, seeing); stars->data[oldSize + i] = star; } - stars->n = oldSize + num; + stars->n = oldSize + nTotal; return true; Index: trunk/ppSim/src/ppSimSequence.c =================================================================== --- trunk/ppSim/src/ppSimSequence.c (revision 17032) +++ trunk/ppSim/src/ppSimSequence.c (revision 17557) @@ -1,8 +1,10 @@ # include "ppSimSequence.h" +# include int main (int argc, char **argv) { + bool status; int argNum; - bool status; + char line[1024]; unsigned int nFail; @@ -17,12 +19,14 @@ char *path = NULL; - if (psArgumentGet (argc, argv, "-path")) { + if ((argNum = psArgumentGet (argc, argv, "-path"))) { psArgumentRemove(argNum, &argc, argv); path = psStringCopy (argv[argNum]); psArgumentRemove(argNum, &argc, argv); + sprintf (line, "mkdir -p %s", path); + system (line); } char *workdir = NULL; - if (psArgumentGet (argc, argv, "-workdir")) { + if ((argNum = psArgumentGet (argc, argv, "-workdir"))) { psArgumentRemove(argNum, &argc, argv); workdir = psStringCopy (argv[argNum]); @@ -31,5 +35,5 @@ char *basename = NULL; - if (psArgumentGet (argc, argv, "-basename")) { + if ((argNum = psArgumentGet (argc, argv, "-basename"))) { psArgumentRemove(argNum, &argc, argv); basename = psStringCopy (argv[argNum]); @@ -130,4 +134,5 @@ exit (1); } + exit (0); } Index: trunk/ppSim/src/ppSimSetPSF.c =================================================================== --- trunk/ppSim/src/ppSimSetPSF.c (revision 17032) +++ trunk/ppSim/src/ppSimSetPSF.c (revision 17557) @@ -1,10 +1,12 @@ # include "ppSim.h" -static char *defaultModel = "PS_MODEL_GAUSS"; +static char *defaultModel = "PS_MODEL_QGAUSS"; bool ppSimSetPSF (pmChip *chip, pmConfig *config) { - bool status; + bool status, mdok; pmPSF *psf = NULL; pmTrend2D *param = NULL; + + psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe // the pmPSF IO functions stores the PSF on the chip->analysis @@ -14,7 +16,8 @@ } - // no supplied PSF, build one using supplied value for seeing - // XXX need to correct for the pixel scale - float seeing = psMetadataLookupF32(&status, config->arguments, "SEEING"); // Seeing sigma (pix) + // no supplied PSF, build one using supplied value for seeing seeing is already corrected + // for the pixel scale, and is converted from FWHM to SIGMA (this is done in + // ppSimArguments) + float seeing = psMetadataLookupF32(&status, recipe, "SEEING"); // Seeing SIGMA (pixels) char *psfModelName = psMetadataLookupStr(&status, config->arguments, "PSF.MODEL"); // Name of PSF model @@ -50,5 +53,5 @@ psEllipsePol pol; - // supply the semi-major axis + // supply the semi-major axis (these are SIGMA values in PIXELS) axes.major = seeing; axes.minor = seeing; Index: trunk/ppSim/src/ppSimStars.c =================================================================== --- trunk/ppSim/src/ppSimStars.c (revision 17032) +++ trunk/ppSim/src/ppSimStars.c (revision 17557) @@ -26,2 +26,34 @@ return galaxy; } + +float ppSimStarSkyNoise (float skySigma, float seeingSigma) { + + float skyNoise = skySigma * sqrt(4*M_PI*PS_SQR(seeingSigma)); + return skyNoise; +} + +float ppSimStarPeakToFlux (float peak, float seeingSigma) { + + float psfArea = 2.0*M_PI*PS_SQR(seeingSigma); + float flux = peak * psfArea; + return flux; +} + +float ppSimStarFluxToPeak (float flux, float seeingSigma) { + + float psfArea = 2.0*M_PI*PS_SQR(seeingSigma); + float peak = flux / psfArea; + return peak; +} + +float ppSimFluxToMag (float flux, float zp) { + + float mag = -2.5*log10(flux) + zp; + return mag; +} + +float ppSimMagToFlux (float mag, float zp) { + + float flux = powf (10.0, -0.4*(mag - zp)); + return flux; +} Index: trunk/ppSim/src/ppSimUtils.c =================================================================== --- trunk/ppSim/src/ppSimUtils.c (revision 17032) +++ trunk/ppSim/src/ppSimUtils.c (revision 17557) @@ -8,9 +8,12 @@ pmCell *cell) { + bool mdok; - float ra0 = psMetadataLookupF32(NULL, config->arguments, "RA"); // Boresight RA (radians) - float dec0 = psMetadataLookupF32(NULL, config->arguments, "DEC"); // Boresight Dec (radians) - float pa = psMetadataLookupF32(NULL, config->arguments, "PA"); // Position angle (radians) - float scale = psMetadataLookupF32(NULL, config->arguments, "SCALE"); // plate scale in arcsec / pixel + psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe + + float ra0 = psMetadataLookupF32(NULL, recipe, "RA"); // Boresight RA (radians) + float dec0 = psMetadataLookupF32(NULL, recipe, "DEC"); // Boresight Dec (radians) + float pa = psMetadataLookupF32(NULL, recipe, "PA"); // Position angle (radians) + float scale = psMetadataLookupF32(NULL, recipe, "SCALE"); // plate scale in arcsec / pixel scale *= M_PI / 3600.0 / 180.0; // convert plate scale to radians/pixel @@ -94,5 +97,9 @@ bool ppSimUpdateConceptsFPA (pmFPA *fpa, pmConfig *config) { - float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time + bool mdok; + + psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe + + float expTime = psMetadataLookupF32(NULL, recipe, "EXPTIME"); // Exposure time const char *filter = psMetadataLookupStr(NULL, config->arguments, "FILTER"); // Filter name @@ -128,6 +135,10 @@ bool ppSimUpdateConceptsCell (pmCell *cell, pmConfig *config) { + bool mdok; + + psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe + int binning = psMetadataLookupS32(NULL, config->arguments, "BINNING"); // Binning in x and y - float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time + float expTime = psMetadataLookupF32(NULL, recipe, "EXPTIME"); // Exposure time psMetadataAddF32(cell->concepts, PS_LIST_TAIL, "CELL.EXPOSURE", PS_META_REPLACE,