Changeset 34134

Timestamp:

Jul 10, 2012, 2:34:42 PM (14 years ago)

Author:

eugene

Message:

a bit of cleanup on the galaxy creation; add TRAIL model for galaxies

Location:

branches/eam_branches/ipp-20120627/ppSim/src

Files:

• ppSimInsertGalaxies.c (modified) (4 diffs)
• ppSimMakeGalaxies.c (modified) (6 diffs)

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

@@ -76 +76 @@
     if (type == -1) {
         psError (PS_ERR_UNKNOWN, false, "invalid model name");
-        return false;
+        return false;
     }
 
@@ -136 +136 @@
         PAR[PM_PAR_YPOS] = yChip;
 
-        psEllipseAxes axes;
-        axes.major       = galaxy->Rmaj;
-        axes.minor       = galaxy->Rmin;
-        axes.theta       = galaxy->theta;
-        pmPSF_AxesToModel (PAR, axes, type);
-        psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
+        if (type == pmModelClassGetType ("PS_MODEL_TRAIL")) {
+            PAR[PM_PAR_LENGTH] = galaxy->Rmaj;
+            PAR[PM_PAR_SIGMA]  = galaxy->Rmin;
+            PAR[PM_PAR_THETA]  = galaxy->theta;
+        } else {
+            psEllipseAxes axes;
+            axes.major       = galaxy->Rmaj;
+            axes.minor       = galaxy->Rmin;
+            axes.theta       = galaxy->theta;
+            pmPSF_AxesToModel (PAR, axes, type);
+        }
+        psF64 Area = 2.0 * M_PI * galaxy->Rmaj * galaxy->Rmin;
 
         if (nParam == 8) {
@@ -147 +153 @@
         }
 
-        // set the normalization to get the desired flux
-        // XXX for now, use peak: pmModelSetFlux (model, galaxy->flux);
-
         // XXX let the flux limit be a user-defined number of sky sigmas (not just 1.0)
         float radius = model->modelRadius (model->params, 0.1*roughNoise);
         radius = PS_MAX (radius, 1.0);
+
         // XXX the exp(-r^0.25) models can go way out if allowed...
         radius = PS_MIN (radius, 300.0); 
@@ -176 +180 @@
 
         float par8 = (model->params->n == 8) ? model->params->data.F32[7] : 0.0;
-        fprintf (outfile, "%8.3f %8.3f %10.2f  %2d  %7.3f %5.3f  %5.3f %5.3f %5.3f %5.3f\n", galaxy->x, galaxy->y, flux, 1, source->psfMag, source->psfMagErr, axes.major, axes.minor, axes.theta, par8);
+        fprintf (outfile, "%8.3f %8.3f %10.2f  %2d  %7.3f %5.3f  %5.3f %5.3f %5.3f %5.3f\n", galaxy->x, galaxy->y, flux, 1, source->psfMag, source->psfMagErr, galaxy->Rmaj, galaxy->Rmin, galaxy->theta, par8);
 
         psArrayAdd (sources, 100,source);

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

@@ -1 +1 @@
 # include "ppSim.h"
+
+bool ppSimSetGalaxyPeak (ppSimGalaxy *galaxy, pmModelType type, float index, float seeing);
 
 bool ppSimMakeGalaxies(psArray *galaxies, pmFPA *fpa, pmConfig *config, const psRandom *rng) {
@@ -19 +21 @@
     int galaxyGridDY      = psMetadataLookupS32 (&mdok, recipe, "GALAXY.GRID.DY"); // Density of fakes
     bool galaxyGridRandom = psMetadataLookupBool(&mdok, recipe, "GALAXY.GRID.RANDOM"); // Density of fakes
-    float galaxyGridPeak  = psMetadataLookupF32 (&mdok, recipe, "GALAXY.GRID.PEAK"); // peak flux 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
     
@@ -42 +44 @@
     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
@@ -150 +152 @@
                 galaxy->theta = (galaxyThetaMin  + rndValue * galaxyThetaSlope);
                 
-                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;
-                }
+                galaxy->flux = galaxyGridFlux;
+                ppSimSetGalaxyPeak (galaxy, type, index, seeing);
 
                 // if (0) {
@@ -206 +188 @@
             galaxy->theta = (psRandomUniform(rng) * galaxyThetaSlope  + galaxyThetaMin);
 
-            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);
-            }
+            // 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;
+            ppSimSetGalaxyPeak (galaxy, type, index, seeing);
 
             if (1) {
@@ -239 +203 @@
     return true;
 }
+
+bool ppSimSetGalaxyPeak (ppSimGalaxy *galaxy, pmModelType type, float index, float seeing) {
+
+    bool isSersicType = false;
+    if (type == pmModelClassGetType ("PS_MODEL_DEV")) {
+        index = 4.0;
+        isSersicType = true;
+    }
+    if (type == pmModelClassGetType ("PS_MODEL_EXP")) {
+        index = 1.0;
+        isSersicType = true;
+    }
+    if (type == pmModelClassGetType ("PS_MODEL_SERSIC")) {
+        isSersicType = true;
+    }
+
+    if (isSersicType) {
+        // 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 = galaxy->flux / norm / (galaxy->Rmin / galaxy->Rmaj);
+        return true;
+    }
+
+    if (type == pmModelClassGetType ("PS_MODEL_TRAIL")) {
+        // galaxy->Rmaj -> PM_PAR_LENGTH
+        // seeing -> PM_PAR_SIGMA
+        galaxy->Rmin = seeing / 2.35; // Rmin is used for sigma
+        galaxy->peak = galaxy->flux / (galaxy->Rmaj * galaxy->Rmin * 2.0 * sqrt(2.0 * M_PI));
+        return true;
+    }
+
+    galaxy->peak = galaxy->flux / (2 * M_PI * galaxy->Rmaj * galaxy->Rmin);
+    return true;
+}