Changeset 25062 for branches/pap/psphot/src/psphotFake.c

Timestamp:

Aug 13, 2009, 1:23:29 PM (17 years ago)

Author:

Paul Price

Message:

Working on fake sources in psphot. Still not done. Being put on back burner for now.

File:

• branches/pap/psphot/src/psphotFake.c (modified) (3 diffs)

branches/pap/psphot/src/psphotFake.c

@@ -9 +9 @@
 
 
-// Calculate the limiting flux for an image
+// Calculate the limiting magnitude for an image
 //
 // We limit ourselves to calculating the peak flux in the smoothed image, which should be close (modulo
@@ -62 +62 @@
     }
 
-    return thresh * sqrtf(meanVar * factor); // Limiting peak flux in smoothed image
+    return -2.5 * log10(thresh * sqrtf(meanVar * factor));
 }
+
+static bool fakeGenerate(psImage **xSrc, psImage **ySrc, // Positions of sources
+                         const pmReadout *ro,            // Readout of interest
+                         const psMetadata *recipe, // psphot recipe
+                         float refMag              // Reference magnitude
+                         )
+{
+    PM_ASSERT_READOUT_NON_NULL(ro, NULL);
+    PM_ASSERT_READOUT_IMAGE(ro, NULL);
+    PS_ASSERT_METADATA_NON_NULL(recipe, NULL);
+
+    psMetadata *fake = psMetadataLookupMetadata(NULL, recipe, "FAKE"); // Fake information
+    if (!fake) {
+        psError(PSPHOT_ERR_CONFIG, false, "Unable to find FAKE information in recipe");
+        return NULL;
+    }
+
+    psVector *magOffsets = psMetadataLookupVector(NULL, recipe, "FAKE.MAG"); // Magnitude offsets
+    if (!magOffset || magOffset->type.type != PS_TYPE_F32) {
+        psError(PSPHOT_ERR_CONFIG, false, "Unable to find FAKE.MAG F32 vector in recipe");
+        return NULL;
+    }
+    int numBins = magOffset->n;                                     // Number of bins
+    int numSources = psMetadataLookupS32(NULL, recipe, "FAKE.NUM"); // Number of sources for each bin
+
+    int numCols = ro->image->numCols, numRows = ro->image->numRows; // Size of image
+
+
+
+    *xSrc = psImageRecycle(*xSrc, numSources, numBins, PS_TYPE_F32);
+    *ySrc = psImageRecycle(*ySrc, numSources, numBins, PS_TYPE_F32);
+
+    // Master list, for image creation
+    psVector *xAll = psVectorAlloc(numBins * numSources, PS_TYPE_F32);
+    psVector *yAll = psVectorAlloc(numBins * numSources, PS_TYPE_F32);
+    psVector *magAll = psVectorAlloc(numBins * numSources, PS_TYPE_F32);
+
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
+    for (int i = 0; i <= numBins; i++) {
+        psArray *data = psArrayAlloc(3); // Sources for bin
+
+        psVector *x = psVectorAlloc(numSources, PS_TYPE_F32);
+        psVector *y = psVectorAlloc(numSources, PS_TYPE_F32);
+
+        float mag = refMag + magOffset->data.F32[i]; // Instrumental magnitude of sources
+
+        for (int j = 0; j <= numSources; j++) {
+            x->data.F32[j] = psRandomUniform(rng) * numCols;
+            y->data.F32[j] = psRandomUniform(rng) * numRows;
+            magAll->data.F32[i * numSources + j] = mag;
+        }
+        memcpy(xSrc->data.F32[i], x->data.F32, numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
+        memcpy(ySrc->data.F32[i], y->data.F32, numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
+        memcpy(&xAll->data.F32[i * numSources], x->data.F32, numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
+        memcpy(&yAll->data.F32[i * numSources], y->data.F32, numSources * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
+    }
+    psFree(rng);
+
+    pmReadout *fakeRO = pmReadoutAlloc(); // Fake readout
+    pmReadoutFakeFromVectors(fakeRO, xAll, yAll, mag,
+
 
 
@@ -80 +141 @@
 
     // Generate array of fake sources and define the pixels with pmSourceDefinePixels
-    psArray *fakeOrig = psphotFakeGenerate(config, readout, recipe);
+    psArray *fakeOrig = fakeGenerate(config, readout, recipe);
 
     psArray *fakeMeasured = psArrayAlloc(fakeOrig->n); // Measured values for the sources in fakeOrig