Index: /branches/eam_branch_20080511/ppSim/src/Makefile.am =================================================================== --- /branches/eam_branch_20080511/ppSim/src/Makefile.am (revision 17704) +++ /branches/eam_branch_20080511/ppSim/src/Makefile.am (revision 17705) @@ -28,4 +28,5 @@ ppSimMergeReadouts.c \ ppSimMosaicChip.c \ + ppSimRandomGaussian.c \ ppSimBadPixels.c Index: /branches/eam_branch_20080511/ppSim/src/ppSim.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSim.c (revision 17704) +++ /branches/eam_branch_20080511/ppSim/src/ppSim.c (revision 17705) @@ -2,4 +2,5 @@ int failure (pmConfig *config, psExit status, char *message) { + ppSimRandomGaussianFree (); psErrorStackPrint(stderr, message); psFree(config); @@ -31,4 +32,5 @@ } + ppSimRandomGaussianFree (); psFree(config); pmModelClassCleanup(); Index: /branches/eam_branch_20080511/ppSim/src/ppSim.h =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSim.h (revision 17704) +++ /branches/eam_branch_20080511/ppSim/src/ppSim.h (revision 17705) @@ -176,3 +176,7 @@ bool ppSimMergeReadouts (pmReadout *outReadout, pmReadout *inReadout); +double ppSimRandomGaussian (const psRandom *rnd, double mean, double sigma); +double ppSimRandomGaussianNorm (const psRandom *rnd); +void ppSimRandomGaussianFree(); + #endif Index: /branches/eam_branch_20080511/ppSim/src/ppSimAddNoise.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimAddNoise.c (revision 17704) +++ /branches/eam_branch_20080511/ppSim/src/ppSimAddNoise.c (revision 17705) @@ -32,5 +32,5 @@ // (generate a static array with the cumulative distribution, use the // random number to select a bin from the histogram) - signal->data.F32[y][x] += sqrtf(variance->data.F32[y][x]) * psRandomGaussian(rng); + signal->data.F32[y][x] += sqrtf(variance->data.F32[y][x]) * ppSimRandomGaussianNorm(rng); signal->data.F32[y][x] /= gain; // Converting to ADU } @@ -40,2 +40,3 @@ return true; } + Index: /branches/eam_branch_20080511/ppSim/src/ppSimAddOverscan.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimAddOverscan.c (revision 17704) +++ /branches/eam_branch_20080511/ppSim/src/ppSimAddOverscan.c (revision 17705) @@ -7,4 +7,5 @@ // skip this step if we did not generate a bias level if (biasCols == NULL) return true; + if (biasRows == NULL) return true; psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe Index: /branches/eam_branch_20080511/ppSim/src/ppSimArguments.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimArguments.c (revision 17704) +++ /branches/eam_branch_20080511/ppSim/src/ppSimArguments.c (revision 17705) @@ -22,5 +22,4 @@ { bool status; - bool mdok; // Status of MD lookup assert(config); @@ -154,5 +153,5 @@ ppSimType type = ppSimTypeFromString (typeStr); - char *filter = ppSimArgToRecipeStr(&status, options, "FILTER", arguments, "-filter"); // Filter name + ppSimArgToRecipeStr(&status, options, "FILTER", arguments, "-filter"); // Filter name // set the exposure time @@ -178,7 +177,4 @@ ppSimArgToRecipeS32(&status, options, "BINNING", arguments, "-bin"); - // XX don't do this: recipe is not yet fully realized (camera options not interpolated) - psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe - if (type == PPSIM_TYPE_OBJECT) { // Load values required for adding stars @@ -195,5 +191,5 @@ } - float scale = ppSimArgToRecipeF32(&status, options, "PIXEL.SCALE", arguments, "-scale"); // Plate scale + ppSimArgToRecipeF32(&status, options, "PIXEL.SCALE", arguments, "-scale"); // Plate scale (arcsec / pixel) ppSimArgToRecipeF32(&status, options, "SKY.MAGS", arguments, "-skymags"); // Plate scale ppSimArgToRecipeStr(&status, options, "PSF.MODEL", arguments, "-psfclass"); // PSF model class @@ -201,16 +197,10 @@ // 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, "SEEING", 0, "Seeing FWHM (arcsec)", seeing); 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); - assert (filter != NULL); - ppSimArgToRecipeF32(&status, options, "ZEROPOINT", arguments, "-zp"); // Zero point - if (!status) { - float zp = ppSimGetZeroPoint (recipe, filter); - psMetadataAddF32(options, PS_LIST_TAIL, "ZEROPOINT", 0, "Photometric zeropoint", zp); - } } Index: /branches/eam_branch_20080511/ppSim/src/ppSimCreate.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimCreate.c (revision 17704) +++ /branches/eam_branch_20080511/ppSim/src/ppSimCreate.c (revision 17705) @@ -53,5 +53,36 @@ config->formatName = psStringCopy (output->formatName); + // the recipe is now fully realized for the desired camera psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe + + char *typeStr = psMetadataLookupStr(NULL, recipe, "IMAGE.TYPE"); // Type of image to simulate + ppSimType type = ppSimTypeFromString (typeStr); // Type of image to simulate + + if (type == PPSIM_TYPE_OBJECT) { + // adjust the seeing by the scale + float seeing = psMetadataLookupF32(&status, recipe, "SEEING"); + if (isnan(seeing)) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "seeing is not defined"); + psFree(fpa); + return NULL; + } + float scale = psMetadataLookupF32(&status, recipe, "PIXEL.SCALE"); + if (isnan(scale)) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "pixel scale is not defined"); + psFree(fpa); + return NULL; + } + psMetadataAddF32(recipe, PS_LIST_TAIL, "SEEING", PS_META_REPLACE, "Seeing SIGMA (pixels)", seeing / 2.0 / sqrt(2.0 * log(2.0)) / scale); + } + + if ((type == PPSIM_TYPE_OBJECT) || (type == PPSIM_TYPE_FLAT)) { + // determine the zeropoint from the filter + float zp = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); + if (isnan(zp)) { + char *filter = psMetadataLookupStr(&status, recipe, "FILTER"); + float zp = ppSimGetZeroPoint (recipe, filter); + psMetadataAddF32(recipe, PS_LIST_TAIL, "ZEROPOINT", PS_META_REPLACE, "Photometric zeropoint", zp); + } + } // For photometry, we operate on the chip-mosaicked image. we create a copy of the mosaicked @@ -115,6 +146,6 @@ psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe - ppSimArgToRecipeF32(&status, recipe, "EXPTIME", fpa->concepts, "FPA.EXPOSURE"); - char *filter = ppSimArgToRecipeStr(&status, recipe, "FILTER", fpa->concepts, "FPA.FILTER"); + ppSimArgToRecipeF32(&status, recipe, "EXPTIME", input->fpa->concepts, "FPA.EXPOSURE"); + char *filter = ppSimArgToRecipeStr(&status, recipe, "FILTER", input->fpa->concepts, "FPA.FILTER"); float zp = ppSimGetZeroPoint (recipe, filter); Index: /branches/eam_branch_20080511/ppSim/src/ppSimInsertStars.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimInsertStars.c (revision 17704) +++ /branches/eam_branch_20080511/ppSim/src/ppSimInsertStars.c (revision 17705) @@ -121,4 +121,5 @@ // insert the source flux in the image pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY); + pmSourceAddWithOffset (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, 0xff, dX, dY); psArrayAdd (sources, 100,source); psFree(source); // Drop reference Index: /branches/eam_branch_20080511/ppSim/src/ppSimLoop.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimLoop.c (revision 17704) +++ /branches/eam_branch_20080511/ppSim/src/ppSimLoop.c (revision 17705) @@ -170,5 +170,12 @@ if (inputCell && inputCell->readouts->data[i]) { - if (!ppSimMergeReadouts (readout, inputCell->readouts->data[i])) ESCAPE (PS_ERR_UNKNOWN, "problem merging input image with simulated image"); + pmReadout *inReadout = inputCell->readouts->data[i]; + if (!inReadout->weight) { + if (!pmReadoutGenerateWeight(inReadout, true)) { + psError (PS_ERR_UNKNOWN, false, "trouble creating weight"); + return false; + } + } + if (!ppSimMergeReadouts (readout, inReadout)) ESCAPE (PS_ERR_UNKNOWN, "problem merging input image with simulated image"); } } Index: /branches/eam_branch_20080511/ppSim/src/ppSimMakeBiassec.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimMakeBiassec.c (revision 17704) +++ /branches/eam_branch_20080511/ppSim/src/ppSimMakeBiassec.c (revision 17705) @@ -6,4 +6,7 @@ psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe + + bool bias = psMetadataLookupBool(&mdok, recipe, "BIAS"); // Generate a Bias? + if (!bias) return NULL; psList *biassec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.BIASSEC"); // Bias regions Index: /branches/eam_branch_20080511/ppSim/src/ppSimMakeSky.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimMakeSky.c (revision 17704) +++ /branches/eam_branch_20080511/ppSim/src/ppSimMakeSky.c (revision 17705) @@ -17,5 +17,6 @@ psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe - bool sky = psMetadataLookupBool(&status, recipe, "SKY"); // Generate a SKY flux? + bool sky = psMetadataLookupBool(&status, recipe, "SKY"); // Generate a SKY flux? + bool flat = psMetadataLookupBool(&status, recipe, "FLAT"); // Apply flat-field term? float expTime = psMetadataLookupF32(&status, recipe, "EXPTIME"); // Exposure time @@ -24,10 +25,10 @@ float flatRate = psMetadataLookupF32(&status, recipe, "FLAT.RATE"); // Flat-field rate float shutterTime = psMetadataLookupF32(&status, recipe, "SHUTTER.TIME"); // Shutter time - float scatterScale = psMetadataLookupF32(&status, recipe, "SCATTER.SCALE"); // scattered light term + float scatterFrac = psMetadataLookupF32(&status, recipe, "SCATTER.FRAC"); // scattered light fraction (max) float skyRate = psMetadataLookupF32(&status, recipe, "SKY.RATE"); // Sky rate - if (isnan(skyRate)) { + float skyMags = psMetadataLookupF32(&status, recipe, "SKY.MAGS"); assert (status); + if (!isnan(skyMags)) { float zp = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); assert (status); float scale = psMetadataLookupF32(&status, recipe, "PIXEL.SCALE"); assert (status); - float skyMags = psMetadataLookupF32(&status, recipe, "SKY.MAGS"); assert (status); skyRate = scale * scale * ppSimMagToFlux (skyMags, zp); } @@ -69,6 +70,10 @@ float realExpTime = expTime + shutterTime * (xFPA + yFPA + 2.0) / 4.0; - // Gaussian flat-field over the FPA - float flatValue = expf(-0.5 / PS_SQR(flatSigma) * (PS_SQR(yFPA) + PS_SQR(xFPA))) / flatSigma / sqrtf(2.0 * M_PI); + // Gaussian flat-field over the FPA with flatValue = 1.0 at the field center + float flatValue = 1.0; + if (flat) { + // we make the flat-field have a response of 1.0 at the field center (like a vignetting) + flatValue = expf(-0.5 / PS_SQR(flatSigma) * (PS_SQR(yFPA) + PS_SQR(xFPA))); + } float scatterRate = 0.0; @@ -78,9 +83,9 @@ if (type == PPSIM_TYPE_FLAT) { float xF = 2.0*(xFPA / dXfpa) - 1.0; - scatterRate = scatterScale * PS_SQR(xF); + scatterRate = scatterFrac * PS_SQR(xF); } // Sky background - float skyFlux = (skyRate * flatValue + scatterRate) * realExpTime; // Flux from sky + float skyFlux = (skyRate * (flatValue + scatterRate)) * realExpTime; // Flux from sky signal->data.F32[y][x] += skyFlux; variance->data.F32[y][x] += skyFlux; Index: /branches/eam_branch_20080511/ppSim/src/ppSimRandomGaussian.c =================================================================== --- /branches/eam_branch_20080511/ppSim/src/ppSimRandomGaussian.c (revision 17705) +++ /branches/eam_branch_20080511/ppSim/src/ppSimRandomGaussian.c (revision 17705) @@ -0,0 +1,102 @@ +# include "ppSim.h" + +static int Ngaussint = 0; +static double *gaussint = NULL; + +extern double drand48(); + +double p_ppSimGaussian (double x, double mean, double sigma) { + + double f; + + f = exp (-0.5 * PS_SQR(x - mean) / PS_SQR(sigma)) / sqrt(2 * M_PI * PS_SQR(sigma)); + + return (f); + +} + +void ppSimRandomGaussianFree() +{ + psFree (gaussint); + return; +} + +void ppSimRandomGaussianAlloc (int Nbin) { + + gaussint = (double *) psAlloc(Nbin*sizeof(double)); + return; +} + +/* integrate a gaussian from -5 sigma to +5 sigma */ +void p_ppSimRandomGaussianInit () { + + int i; + long A, B; + double val, x, dx, dx1, dx2, dx3, df; + double mean, sigma; + + /* no need to generate this if it already exists */ + if (gaussint) return; + + A = time(NULL); + for (B = 0; A == time(NULL); B++); + srand48(B); + + Ngaussint = 0x1000; + ppSimRandomGaussianAlloc (Ngaussint + 1); + + val = 0; + dx = 1.0 / Ngaussint; + dx1 = dx / 3.0; + dx2 = 2.0*dx/3.0; + dx3 = dx; + mean = 0.0; + sigma = 1.0; + + for (i = 0, x = -7.0; (i < Ngaussint) && (x < 7.0); x += dx) { + df = (3.0*p_ppSimGaussian(x , mean, sigma) + + 9.0*p_ppSimGaussian(x+dx1, mean, sigma) + + 9.0*p_ppSimGaussian(x+dx2, mean, sigma) + + 3.0*p_ppSimGaussian(x+dx3, mean, sigma)) * (dx1/8.0); + val += df; + if (val > (i + 0.5) / (double) Ngaussint) { + gaussint[i] = x + dx / 2.0; + i++; + } + } +} + +// XXX we are using drand48() rather than the random var supplied by rnd +double ppSimRandomGaussian (const psRandom *rnd, double mean, double sigma) { + + int i; + double y; + + if (gaussint == NULL) { + p_ppSimRandomGaussianInit (); + } + + y = drand48(); + i = Ngaussint*y; + y = gaussint[i]*sigma + mean; + + return (y); + +} + +// XXX we are using drand48() rather than the random var supplied by rnd +double ppSimRandomGaussianNorm (const psRandom *rnd) { + + int i; + double y; + + if (gaussint == NULL) { + p_ppSimRandomGaussianInit (); + } + + y = drand48(); + i = Ngaussint*y; + y = gaussint[i]; + + return (y); +}