Changeset 32159

Timestamp:

Aug 21, 2011, 10:46:36 AM (15 years ago)

Author:

eugene

Message:

set model flux for sersic-like models using analytical form; use GRID.PEAK for non-sersic like models; galaxy model flux should now be (nearly) exact

Location:

branches/eam_branches/ipp-20110710/ppSim/src

Files:

• ppSimInsertGalaxies.c (modified) (1 diff)
• ppSimInsertStars.c (modified) (1 diff)
• ppSimMakeGalaxies.c (modified) (7 diffs)

branches/eam_branches/ipp-20110710/ppSim/src/ppSimInsertGalaxies.c

@@ -137 +137 @@
         axes.minor       = galaxy->Rmin;
         axes.theta       = galaxy->theta;
-        pmPSF_AxesToModel (PAR, axes);
+        pmPSF_AxesToModel (PAR, axes, type);
         psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
 

branches/eam_branches/ipp-20110710/ppSim/src/ppSimInsertStars.c

@@ -117 +117 @@
         // XXX set the mag & err values (should this be done in pmSourceFromModel?)
         // XXX i should be applying the gain and the correct effective area
-        psEllipseAxes axes = pmPSF_ModelToAxes (model->params->data.F32, 20.0);
+        psEllipseAxes axes = pmPSF_ModelToAxes (model->params->data.F32, 20.0, model->type);
         float Area = 2.0 * M_PI * axes.major * axes.minor;
 

branches/eam_branches/ipp-20110710/ppSim/src/ppSimMakeGalaxies.c

@@ -20 +20 @@
     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
@@ -46 +50 @@
     }
 
+    // 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;
+    }
+
     if (galaxyDensity <= 0) return true;
 
@@ -87 +99 @@
     long num = expf(logf(norm) + galaxyLum * logf(faint)) + 0.5;
 
+    fprintf (stderr, "galaxy grid inst mag: %f\n", galaxyGridMag - zp - 2.5*log10(expTime));
+
     if (galaxyGrid) {
         num = (int)(xSize / galaxyGridDX) * (int)(ySize / galaxyGridDY);
@@ -98 +112 @@
 
         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 +120 @@
 
                 // 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))
@@ -122 +138 @@
                 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;
 
@@ -135 +146 @@
                 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);