Index: /trunk/ppSim/src/Makefile.am =================================================================== --- /trunk/ppSim/src/Makefile.am (revision 41172) +++ /trunk/ppSim/src/Makefile.am (revision 41173) @@ -15,4 +15,5 @@ ppSimLoadStars.c \ ppSimMakeStars.c \ + ppSimMakeStarCluster.c \ ppSimMakeStarGrid.c \ ppSimMakeGalaxies.c \ Index: /trunk/ppSim/src/ppSim.h =================================================================== --- /trunk/ppSim/src/ppSim.h (revision 41172) +++ /trunk/ppSim/src/ppSim.h (revision 41173) @@ -103,4 +103,5 @@ bool ppSimLoadStars (psArray *stars, pmFPA *fpa, pmConfig *config); bool ppSimMakeStars(psArray *stars, pmFPA *fpa, pmConfig *config, const psRandom *rng); +bool ppSimMakeStarCluster(psArray *stars, pmFPA *fpa, pmConfig *config, const psRandom *rng); bool ppSimMakeStarGrid(psArray *stars, pmFPA *fpa, pmConfig *config, const psRandom *rng); bool ppSimInsertStars (pmReadout *readout, psImage *expCorr, psArray *stars, pmConfig *config); Index: /trunk/ppSim/src/ppSimCreate.c =================================================================== --- /trunk/ppSim/src/ppSimCreate.c (revision 41172) +++ /trunk/ppSim/src/ppSimCreate.c (revision 41173) @@ -159,5 +159,8 @@ return false; } - simSources->save = true; + + // XXXX TEST this is causing trouble for unknown reasons -- output step is blowing up. + // simSources->save = true; + simSources->save = false; // if we have loaded an input image, we derive certain values from the image, if possible Index: /trunk/ppSim/src/ppSimLoadStars.c =================================================================== --- /trunk/ppSim/src/ppSimLoadStars.c (revision 41172) +++ /trunk/ppSim/src/ppSimLoadStars.c (revision 41173) @@ -24,10 +24,11 @@ float pa = psMetadataLookupF32(NULL, recipe, "PA"); // Position angle (radians) float seeing = psMetadataLookupF32(NULL, recipe, "SEEING"); // Seeing SIGMA (pixels) - float scale = psMetadataLookupF32(NULL, recipe, "PIXEL.SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel) + float scale = psMetadataLookupF32(NULL, recipe, "PIXEL.SCALE"); // Plate scale (arcsec/pixel) + float scaleRad = scale * M_PI / 3600.0 / 180.0; // Plate scale in radians/pixel for WCS below float expTime = psMetadataLookupF32(NULL, recipe, "EXPTIME"); // Exposure time (sec) // Size of FPA psRegion *bounds = ppSimFPABounds (fpa); - float radius = 0.5 * PS_MAX(bounds->x1 - bounds->x0, bounds->y1 - bounds->y0) * scale; + float radius = 0.5 * PS_MAX(bounds->x1 - bounds->x0, bounds->y1 - bounds->y0) * scaleRad; float x0fpa = 0.5*(bounds->x0 + bounds->x1); @@ -53,5 +54,5 @@ stars = psArrayRealloc (stars, refStars->n); - psProjection *proj = psProjectionAlloc(ra0, dec0, scale, scale, PS_PROJ_TAN); // Projection + psProjection *proj = psProjectionAlloc(ra0, dec0, scaleRad, scaleRad, PS_PROJ_TAN); // Projection // Conversion loop Index: /trunk/ppSim/src/ppSimLoop.c =================================================================== --- /trunk/ppSim/src/ppSimLoop.c (revision 41172) +++ /trunk/ppSim/src/ppSimLoop.c (revision 41173) @@ -51,5 +51,8 @@ // Add random stars - if (!ppSimMakeStarGrid (stars, fpa, config, rng)) ESCAPE (PS_ERR_UNKNOWN, "failed to make random stars"); + if (!ppSimMakeStarGrid (stars, fpa, config, rng)) ESCAPE (PS_ERR_UNKNOWN, "failed to make star grid"); + + // Add random stars + if (!ppSimMakeStarCluster (stars, fpa, config, rng)) ESCAPE (PS_ERR_UNKNOWN, "failed to make cluster stars"); // Add random galaxies Index: /trunk/ppSim/src/ppSimMakeGalaxies.c =================================================================== --- /trunk/ppSim/src/ppSimMakeGalaxies.c (revision 41172) +++ /trunk/ppSim/src/ppSimMakeGalaxies.c (revision 41173) @@ -45,9 +45,10 @@ float zp = psMetadataLookupF32(&mdok, recipe, "ZEROPOINT"); // Photometric zero point float seeing = psMetadataLookupF32(&mdok, recipe, "SEEING"); // Seeing sigma (pix) - float scale = psMetadataLookupF32(&mdok, recipe, "PIXEL.SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel) + float scale = psMetadataLookupF32(&mdok, recipe, "PIXEL.SCALE"); // Plate scale (arcsec/pixel) float skyRate = psMetadataLookupF32(&mdok, recipe, "SKY.RATE"); // Sky rate if (isnan(skyRate)) { float skyMags = psMetadataLookupF32(&mdok, recipe, "SKY.MAGS"); assert (mdok); skyRate = scale * scale * ppSimMagToFlux (skyMags, zp); + // skyMags is in mags / square arcsec so scale must be arcsec / pixel } @@ -94,5 +95,5 @@ // Normalisation, set by the specified stellar density at the specified bright magnitude - float refSum = galaxyDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI); + float refSum = galaxyDensity * xSize * ySize * PS_SQR(scale / 3600.0); float normLum = refSum / (galaxyAlpha * logf(10.0) * powf(10.0, (galaxyAlpha * brightMag))); float normScale = normLum * galaxyAlpha * logf(10.0); @@ -143,5 +144,4 @@ rndValue = galaxyGridRandom ? drand48() : i / (float) nTotal; float scale = (galaxyRmajorMin + rndValue * galaxyRmajorSlope); - galaxy->Rmaj = scale; Index: /trunk/ppSim/src/ppSimMakeStarCluster.c =================================================================== --- /trunk/ppSim/src/ppSimMakeStarCluster.c (revision 41173) +++ /trunk/ppSim/src/ppSimMakeStarCluster.c (revision 41173) @@ -0,0 +1,157 @@ +# include "ppSim.h" + +# define ESCAPE(MSG) { \ + psError(PS_ERR_BAD_PARAMETER_VALUE, true, MSG); \ + return false; } + +bool ppSimMakeStarCluster(psArray *stars, pmFPA *fpa, pmConfig *config, const psRandom *rng) { + + bool status; + assert (stars); + + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe + + bool starsCluster = psMetadataLookupBool(&status, recipe, "STARS.CLUSTER"); // add fake stars in a cluster? + if (!starsCluster) return true; + + // cluster stars are added with magnitudes distributed based on the given luminosity function: + float starsAlpha = psMetadataLookupF32(&status, recipe, "STARS.LUM"); // Star luminosity func slope + float brightMag = psMetadataLookupF32(&status, recipe, "STARS.MAG"); // Star brightest magnitude + float nSigmaLim = psMetadataLookupF32(&status, recipe, "STARS.SIGMA.LIM"); // significance of faintest stars + + float starsClusterDensity = psMetadataLookupF32(&status, recipe, "STARS.CLUSTER.DENSITY"); // Peak density of cluster (units?) + float starsClusterSigma = psMetadataLookupF32(&status, recipe, "STARS.CLUSTER.SIGMA"); // Sigma of cluster distribution + float starsClusterXcenter = psMetadataLookupF32(&status, recipe, "STARS.CLUSTER.XCENTER"); // Xo of cluster distribution + float starsClusterYcenter = psMetadataLookupF32(&status, recipe, "STARS.CLUSTER.YCENTER"); // Yo of cluster distribution + + 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, "PIXEL.SCALE"); // Plate scale (arcsec/pixel) + + if (isnan(expTime)) ESCAPE("EXPTIME is not defined"); + if (isnan(zp)) ESCAPE("ZEROPOINT is not defined"); + if (isnan(seeing)) ESCAPE("SEEING is not defined"); + if (isnan(scale)) ESCAPE("PIXEL.SCALE is not defined"); + + float skyRate = psMetadataLookupF32(&status, recipe, "SKY.RATE"); // Sky rate (counts / sec) + if (isnan(skyRate)) { + float skyMags = psMetadataLookupF32(&status, recipe, "SKY.MAGS"); assert (status); + skyRate = scale * scale * ppSimMagToFlux (skyMags, zp); + } + + // Size of FPA + psRegion *bounds = ppSimFPABounds (fpa); + // Size of focal plane (can't I get this from the config?) + int xSize = bounds->x1 - bounds->x0; + int ySize = bounds->y1 - bounds->y0; + psFree(bounds); + + // choose reference magnitude & density to set normalization + float refMag = 0; + float refSum = 0; + refMag = brightMag; + + // Integral of a 2D Gauss with sigma and Io: 2.0 * M_PI * PS_SQR(sigma) * Io + + // central density is in stars per square degree but + // starsClusterSigma is in pixels (or arcsec?) + refSum = starsClusterDensity * 2.0 * M_PI * PS_SQR(starsClusterSigma * scale / 3600.0); + 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(&status, recipe, "READNOISE"); // Default read noise + if (!status) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe."); + psFree(bounds); + return false; + } + + // Faintest and brightest (integrated) fluxes (actually fluence, since integrated) for random stars + + // faint limit is set by detection limit, in turn set by sky noise + // float skySigma = sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime); + // XXX add darkRate back in? + float skySigma = sqrtf(PS_SQR(readnoise) + (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."); + return true; + } + + // 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 nTotal = 0; + // if (flatLum) { + // nTotal = flatNum; + // } else { + // nTotal = normScale * powf (10.0, (starsAlpha * faintMag)); + // } + + nTotal = normScale * powf (10.0, (starsAlpha * faintMag)); + + + psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld stars between %f and %f mag\n", nTotal, brightMag, faintMag); + + long oldSize = stars->n; + psArrayRealloc (stars, stars->n + nTotal); + + for (long i = 0; i < nTotal; i++) { + ppSimStar *star = ppSimStarAlloc (); + + // make 10 tries to select a star in the image range + for (int j = 0; j < 10; j++) { + float Xo = ppSimRandomGaussian (rng, starsClusterXcenter, starsClusterSigma); + if (Xo < 0) continue; + if (Xo > xSize) continue; + star->x = Xo; + break; + } + for (int j = 0; j < 10; j++) { + float Yo = ppSimRandomGaussian (rng, starsClusterYcenter, starsClusterSigma); + if (Yo < 0) continue; + if (Yo > ySize) continue; + star->y = Yo; + break; + } + + float starMag = 0; + // if (flatLum) { + // starMag = psRandomUniform(rng) * (faintMag - brightMag) + brightMag; + // } else { + // starMag = PS_MIN (faintMag, PS_MAX (brightMag, (log10((i + 1.0) / normScale) / starsAlpha))); + // } + + starMag = PS_MIN (faintMag, PS_MAX (brightMag, (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 + nTotal; + + return true; +} Index: /trunk/ppSim/src/ppSimMakeStarGrid.c =================================================================== --- /trunk/ppSim/src/ppSimMakeStarGrid.c (revision 41172) +++ /trunk/ppSim/src/ppSimMakeStarGrid.c (revision 41173) @@ -22,5 +22,5 @@ float zp = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); // Photometric zero point float seeing = psMetadataLookupF32(&status, recipe, "SEEING"); // Seeing SIGMA (pixels) - float scale = psMetadataLookupF32(&status, recipe, "PIXEL.SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel) + float scale = psMetadataLookupF32(&status, recipe, "PIXEL.SCALE"); // Plate scale (arcsec/pixel) if (isnan(expTime)) ESCAPE("EXPTIME is not defined"); Index: /trunk/ppSim/src/ppSimMakeStars.c =================================================================== --- /trunk/ppSim/src/ppSimMakeStars.c (revision 41172) +++ /trunk/ppSim/src/ppSimMakeStars.c (revision 41173) @@ -27,5 +27,9 @@ float zp = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); // Photometric zero point float seeing = psMetadataLookupF32(&status, recipe, "SEEING"); // Seeing SIGMA (pixels) - float scale = psMetadataLookupF32(&status, recipe, "PIXEL.SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel) + float scale = psMetadataLookupF32(&status, recipe, "PIXEL.SCALE"); // Plate scale (arcsec/pixel) + // there has been some confusion over pixel scale: PIXEL.SCALE is supplied in arcsec / pixel. + // In some other places (ppSimLoadStars.c, ppSimUtils.c) it is needed in radius for WCS conversion. + // in the past, it was converted to radians on load (but applied inconsisently elsewhere) + // Now the radian version is carried as scaleRad to be more explicit if (isnan(darkRate)) darkRate = 0.0; @@ -42,4 +46,5 @@ float skyMags = psMetadataLookupF32(&status, recipe, "SKY.MAGS"); assert (status); skyRate = scale * scale * ppSimMagToFlux (skyMags, zp); + // skyMags is in mags / square arcsec, so scale must be in arcsec / pixel } @@ -61,9 +66,9 @@ if (!status) { refMag = brightMag; - refSum = starsDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI); + refSum = starsDensity * xSize * ySize * PS_SQR(scale / 3600.0); } } else { refMag = brightMag; - refSum = starsDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI); + refSum = starsDensity * xSize * ySize * PS_SQR(scale / 3600.0); } psTrace("ppSim", 6, "refMag: %f, refSum: %f\n", refMag, refSum); @@ -120,4 +125,9 @@ psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld stars between %f and %f mag\n", nTotal, brightMag, faintMag); + if (nTotal > 1e7) { + psLogMsg("ppSim", PS_LOG_INFO, "Trying to add %d stars (far more than 10M stars!), giving up\n", (int) nTotal); + exit (2); + } + long oldSize = stars->n; psArrayRealloc (stars, stars->n + nTotal); Index: /trunk/ppSim/src/ppSimUtils.c =================================================================== --- /trunk/ppSim/src/ppSimUtils.c (revision 41172) +++ /trunk/ppSim/src/ppSimUtils.c (revision 41173) @@ -16,5 +16,5 @@ float pa = psMetadataLookupF32(NULL, recipe, "PA"); // Position angle (radians) float scale = psMetadataLookupF32(NULL, recipe, "PIXEL.SCALE"); // plate scale in arcsec / pixel - scale *= M_PI / 3600.0 / 180.0; // convert plate scale to radians/pixel + float scaleRad = scale * M_PI / 3600.0 / 180.0; // convert plate scale to radians/pixel int binning = psMetadataLookupS32(NULL, recipe, "BINNING"); // Binning in x and y @@ -61,7 +61,7 @@ psMetadata *header = psMetadataAlloc(); // Header, to return pmAstromWCS *wcs = pmAstromWCSAlloc(1, 1); // WCS structure - wcs->toSky = psProjectionAlloc(ra0, dec0, scale * xParity, scale * yParity, PS_PROJ_TAN); - wcs->cdelt1 = scale * PM_DEG_RAD * xParity; - wcs->cdelt2 = scale * PM_DEG_RAD * yParity; + wcs->toSky = psProjectionAlloc(ra0, dec0, scaleRad * xParity, scaleRad * yParity, PS_PROJ_TAN); + wcs->cdelt1 = scaleRad * PM_DEG_RAD * xParity; + wcs->cdelt2 = scaleRad * PM_DEG_RAD * yParity; wcs->crpix1 = x0; wcs->crpix2 = y0;