Changeset 34261 for trunk/ppSim/src/ppSimMakeGalaxies.c

Timestamp:

Jul 31, 2012, 4:03:41 PM (14 years ago)

Author:

eugene

Message:

add trails to test images from ppSim

File:

• trunk/ppSim/src/ppSimMakeGalaxies.c (modified) (6 diffs)

trunk/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; // 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;
+}