Index: trunk/ppSim/src/ppSimMakeGalaxies.c =================================================================== --- trunk/ppSim/src/ppSimMakeGalaxies.c (revision 31157) +++ trunk/ppSim/src/ppSimMakeGalaxies.c (revision 32350) @@ -11,6 +11,6 @@ if (!galaxyFake) return true; - float galaxyLum = psMetadataLookupF32(&mdok, recipe, "GALAXY.LUM"); // Galaxy luminosity func slope - float galaxyMag = psMetadataLookupF32(&mdok, recipe, "GALAXY.MAG"); // Galaxy brightest magnitude + float galaxyAlpha = psMetadataLookupF32(&mdok, recipe, "GALAXY.LUM"); // Galaxy luminosity func slope + float brightMag = psMetadataLookupF32(&mdok, recipe, "GALAXY.MAG"); // Galaxy brightest magnitude float galaxyDensity = psMetadataLookupF32(&mdok, recipe, "GALAXY.DENSITY"); // Density of fakes @@ -20,5 +20,9 @@ bool galaxyGridRandom = psMetadataLookupBool(&mdok, recipe, "GALAXY.GRID.RANDOM"); // Density of fakes float galaxyGridPeak = psMetadataLookupF32 (&mdok, recipe, "GALAXY.GRID.PEAK"); // peak flux of fakes + float galaxyGridMag = psMetadataLookupF32 (&mdok, recipe, "GALAXY.GRID.MAG"); // peak flux of fakes + float galaxyGridXoff = psMetadataLookupF32 (&mdok, recipe, "GALAXY.GRID.XOFF"); // centroid offset + float galaxyGridYoff = psMetadataLookupF32 (&mdok, recipe, "GALAXY.GRID.YOFF"); // centroid offset + float galaxyRmajorMax = psMetadataLookupF32(&mdok, recipe, "GALAXY.RMAJOR.MAX"); // Density of fakes float galaxyRmajorMin = psMetadataLookupF32(&mdok, recipe, "GALAXY.RMAJOR.MIN"); // Density of fakes @@ -35,8 +39,8 @@ float galaxyIndexMax = psMetadataLookupF32(&mdok, recipe, "GALAXY.INDEX.MAX"); // Density of fakes - float darkRate = psMetadataLookupF32(&mdok, recipe, "DARK.RATE"); // Dark 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 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 skyRate = psMetadataLookupF32(&mdok, recipe, "SKY.RATE"); // Sky rate @@ -44,4 +48,12 @@ float skyMags = psMetadataLookupF32(&mdok, recipe, "SKY.MAGS"); assert (mdok); skyRate = scale * scale * ppSimMagToFlux (skyMags, zp); + } + + // determine the galaxy model + char *modelName = psMetadataLookupStr(&mdok, recipe, "GALAXY.MODEL"); // galaxy model name + pmModelType type = pmModelClassGetType (modelName); + if (type == -1) { + psError (PS_ERR_UNKNOWN, false, "invalid model name"); + return false; } @@ -64,31 +76,37 @@ // 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) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe."); - psFree(bounds); - return false; - } - - // Peak fluxes: faintest and brightest levels for random galaxy - float faint = 0.1 * 10.0 * sqrt(2.0*M_PI) * seeing * sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime); - float bright = powf(10.0, -0.4 * (galaxyMag - zp)) / sqrt(2.0*M_PI) / seeing * expTime; - if (bright < faint) { - psLogMsg("ppSim", PS_LOG_INFO, - "Image noise is above brightest random galaxy --- no random galaxy added."); + // float readnoise = psMetadataLookupF32(&mdok, recipe, "READNOISE"); // Default read noise + // if (!mdok) { + // psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe."); + // psFree(bounds); + // return false; + // } + + // faintest and brightest magnitudes for random galaxies + // float brightFlux = ppSimMagToFlux (brightMag, zp) * expTime; + float faintMag = brightMag + 5.0; + if (brightMag > faintMag) { + psLogMsg("ppSim", PS_LOG_INFO, "Image noise is above brightest random galaxy --- no random galaxy added."); return true; } // Normalisation, set by the specified stellar density at the specified bright magnitude - float norm = galaxyDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI) / - powf(bright, galaxyLum); + float refSum = galaxyDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI); + float normLum = refSum / (galaxyAlpha * logf(10.0) * powf(10.0, (galaxyAlpha * brightMag))); + float normScale = normLum * galaxyAlpha * logf(10.0); // Total number of galaxy down to the faint flux end - long num = expf(logf(norm) + galaxyLum * logf(faint)) + 0.5; + long nTotal = 0; + if (galaxyAlpha < 0.1) { + nTotal = 100; + } else { + nTotal = normScale * powf (10.0, (galaxyAlpha * faintMag)); + } if (galaxyGrid) { - num = (int)(xSize / galaxyGridDX) * (int)(ySize / galaxyGridDY); + fprintf (stderr, "galaxy grid inst mag: %f\n", galaxyGridMag - zp - 2.5*log10(expTime)); + nTotal = (int)(xSize / galaxyGridDX) * (int)(ySize / galaxyGridDY); - psLogMsg("ppSim", PS_LOG_INFO, "Adding grid of %ld galaxies\n", num); + psLogMsg("ppSim", PS_LOG_INFO, "Adding grid of %ld galaxies\n", nTotal); float galaxyIndexSlope = (galaxyIndexMax - galaxyIndexMin); @@ -98,6 +116,6 @@ int i = 0; - float refNorm = 1.0; - float ourNorm = 1.0; + // float refNorm = 1.0; + // float ourNorm = 1.0; for (long iy = 0.5*galaxyGridDY; iy < ySize; iy += galaxyGridDY) { @@ -106,7 +124,9 @@ // make fpa center of distribution - galaxy->x = ix; - galaxy->y = iy; - galaxy->peak = galaxyGridPeak; + // center the galaxy peak on the center of a pixel + galaxy->x = ix + 0.5 + galaxyGridXoff; + galaxy->y = iy + 0.5 + galaxyGridYoff; + + float galaxyGridFlux = expTime * powf(10.0, -0.4 * (galaxyGridMag - zp)); // galaxyIndex from user should be for function of this form: exp(-r^(1/n)) @@ -115,33 +135,45 @@ float rndValue; - rndValue = galaxyGridRandom ? drand48() : i / (float) num; + rndValue = galaxyGridRandom ? drand48() : i / (float) nTotal; float index = (galaxyIndexMin + rndValue * galaxyIndexSlope); galaxy->index = 0.5/index; // factor of 0.5 because the Sersic model creates exp(-z^n), not exp(-r^n) - rndValue = galaxyGridRandom ? drand48() : i / (float) num; + rndValue = galaxyGridRandom ? drand48() : i / (float) nTotal; float scale = (galaxyRmajorMin + rndValue * galaxyRmajorSlope); - // for a sersic model, - float bn = 1.9992*index - 0.3271; - float fR = 1.0 / (sqrt(2.0) * pow (bn, index)); - float Io = exp(bn); - galaxy->Rmaj = scale; - rndValue = galaxyGridRandom ? drand48() : i / (float) num; + rndValue = galaxyGridRandom ? drand48() : i / (float) nTotal; galaxy->Rmin = (galaxyARatioMin + rndValue * galaxyARatioSlope) * galaxy->Rmaj; - rndValue = galaxyGridRandom ? drand48() : i / (float) num; + rndValue = galaxyGridRandom ? drand48() : i / (float) nTotal; galaxy->theta = (galaxyThetaMin + rndValue * galaxyThetaSlope); - if (i == 0) { - refNorm = Io*scale*scale; - } - ourNorm = refNorm / (Io*scale*scale); - galaxy->peak *= ourNorm; - - if (0) { - fprintf (stderr, "Rmaj: %f, scale: %f, index: %f, bn: %f, Ro: %f, Io: %f, theta: %f\n", galaxy->Rmaj, scale, index, bn, fR, Io, galaxy->theta); + bool sersicType = false; + sersicType |= (type == pmModelClassGetType ("PS_MODEL_SERSIC")); + sersicType |= (type == pmModelClassGetType ("PS_MODEL_EXP")); + sersicType |= (type == pmModelClassGetType ("PS_MODEL_DEV")); + if (sersicType) { + if (type == pmModelClassGetType ("PS_MODEL_DEV")) index = 4.0; + if (type == pmModelClassGetType ("PS_MODEL_EXP")) index = 1.0; + + // for a sersic model, + float bn = 1.9992*index - 0.3271; + float Io = exp(bn); + + // the integral of a Sersic has an analytical form as follows: + float logGamma = lgamma(2.0*index); + float bnFactor = pow(bn, 2.0*index); + float norm = 2.0 * M_PI * PS_SQR(galaxy->Rmaj) * index * Io * exp(logGamma) / bnFactor; + + // XXX probably should limit the allowed thinness of a galaxy + galaxy->peak = galaxyGridFlux / norm / (galaxy->Rmin / galaxy->Rmaj); + } else { + galaxy->peak = galaxyGridPeak; } + + // if (0) { + // fprintf (stderr, "Rmaj: %f, scale: %f, index: %f, bn: %f, Ro: %f, Io: %f, theta: %f\n", galaxy->Rmaj, scale, index, bn, fR, Io, galaxy->theta); + // } psArrayAdd (galaxies, 100, galaxy); @@ -151,4 +183,5 @@ } else { + fprintf (stderr, "galaxy ref inst mag: %f\n", brightMag - zp - 2.5*log10(expTime)); float galaxyIndexSlope = (galaxyIndexMax - galaxyIndexMin); float galaxyThetaSlope = (galaxyThetaMax - galaxyThetaMin); @@ -156,8 +189,7 @@ float galaxyRmajorSlope = (galaxyRmajorMax - galaxyRmajorMin); - psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld galaxies down to %f mag\n", num, - -2.5 * log10(faint * sqrt(2.0*M_PI) * seeing / expTime) + zp); - - for (long i = 0; i < num; i++) { + psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld galaxies down to %f mag\n", nTotal, faintMag); + + for (long i = 0; i < nTotal; i++) { ppSimGalaxy *galaxy = ppSimGalaxyAlloc (); @@ -170,23 +202,34 @@ float scale = (psRandomUniform(rng) * galaxyRmajorSlope + galaxyRmajorMin); - - // for a sersic model, - float bn = 1.9992*index - 0.3271; - float fR = 1.0 / (sqrt(2.0) * pow (bn, index)); - float Io = exp(bn); - galaxy->Rmaj = scale; - galaxy->Rmin = (PS_SQR(psRandomUniform(rng)) * galaxyARatioSlope + galaxyARatioMin) * galaxy->Rmaj; - galaxy->theta = (psRandomUniform(rng) * galaxyThetaSlope + galaxyThetaMin); - // the flux is only approximately correct (scaling of the peak is wrong) - // elsewhere (ppSimInsertGalaxies.c) we calculate the true galaxy flux - galaxy->flux = expf((logf(i + 1) - logf(norm)) / galaxyLum); // Peak flux - galaxy->peak = galaxy->flux / (2.0*M_PI*PS_SQR(seeing)); - - if (0) { - fprintf (stderr, "Rmaj: %f, scale: %f, index: %f, bn: %f, Ro: %f, Io: %f\n", galaxy->Rmaj, scale, index, bn, fR, Io); + bool sersicType = false; + sersicType |= (type == pmModelClassGetType ("PS_MODEL_SERSIC")); + sersicType |= (type == pmModelClassGetType ("PS_MODEL_EXP")); + sersicType |= (type == pmModelClassGetType ("PS_MODEL_DEV")); + if (sersicType) { + if (type == pmModelClassGetType ("PS_MODEL_DEV")) index = 4.0; + if (type == pmModelClassGetType ("PS_MODEL_EXP")) index = 1.0; + + // for a sersic model, + float bn = 1.9992*index - 0.3271; + float Io = exp(bn); + + // the integral of a Sersic has an analytical form as follows: + float logGamma = lgamma(2.0*index); + float bnFactor = pow(bn, 2.0*index); + float norm = 2.0 * M_PI * PS_SQR(galaxy->Rmaj) * index * Io * exp(logGamma) / bnFactor; + + // XXX probably should limit the allowed thinness of a galaxy + float galaxyMag = PS_MIN (faintMag, PS_MAX (brightMag, (log10((i + 1.0) / normScale) / galaxyAlpha))); + + galaxy->flux = ppSimMagToFlux (galaxyMag, zp) * expTime; + galaxy->peak = galaxy->flux / norm / (galaxy->Rmin / galaxy->Rmaj); + } + + if (1) { + fprintf (stderr, "Rmaj: %f, Rmin: %f, theta: %f, scale: %f, index: %f, peak: %f, mag: %f\n", galaxy->Rmaj, galaxy->Rmin, galaxy->theta*180.0/M_PI, scale, index, galaxy->peak, -2.5*log10(galaxy->flux)); }