Changeset 26674 for branches/eam_branches/20091201/ppStack/src/ppStackSources.c

Timestamp:

Jan 22, 2010, 5:32:32 PM (16 years ago)

Author:

Paul Price

Message:

Get running with scaling of kernel parameters.

File:

• branches/eam_branches/20091201/ppStack/src/ppStackSources.c (modified) (3 diffs)

branches/eam_branches/20091201/ppStack/src/ppStackSources.c

@@ -61 +61 @@
     PS_ASSERT_PTR_NON_NULL(config, false);
 
-    if (!options->matchZPs) {
-        int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
-        options->norm = psVectorAlloc(num, PS_TYPE_F32);
-        psVectorInit (options->norm, 0.0);
-
-        // XXX do I need to set this?
-        // options->sumExposure = sumExpTime;
-
-        return true;
+    if (!options->matchZPs && !options->photometry) {
+        int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
+        options->norm = psVectorAlloc(num, PS_TYPE_F32);
+        psVectorInit (options->norm, 0.0);
+
+        // XXX do I need to set this?
+        // options->sumExposure = sumExpTime;
+
+        return true;
     }
 
@@ -191 +191 @@
 #endif
 
-    psVector *trans = pmSourceMatchRelphot(matches, zp, tol, iter1, starRej1, starSys1,
-                                           iter2, starRej2, starSys2, starLimit,
-                                           transIter, transRej, transThresh); // Transparencies for each image
-    if (!trans) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to measure transparencies");
-        return false;
-    }
+    if (options->matchZPs) {
+        psVector *trans = pmSourceMatchRelphot(matches, zp, tol, iter1, starRej1, starSys1,
+                                               iter2, starRej2, starSys2, starLimit,
+                                               transIter, transRej, transThresh); // Transparencies per image
+        if (!trans) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to measure transparencies");
+            return false;
+        }
+        for (int i = 0; i < trans->n; i++) {
+            if (!isfinite(trans->data.F32[i])) {
+                inputMask->data.U8[i] |= PPSTACK_MASK_CAL;
+            }
+        }
+        // M = m + c0 + c1 * airmass - 2.5log(t) + transparency
+        // Want sources to have m corresponding to airmass = 1 and t = sumExpTime and transparency = 0
+        // m_0 + c1 * airmass_0 + 2.5log(t_0) - trans_0 = m_1 + c1 * airmass_1 + 2.5log(t_1) - trans_1
+        // m_0 = m_1 + zp_1 + trans_1 - c1 * airmass_0 - 2.5log(t_0)
+        // We don't need to know the magnitude zero point for the filter, since it cancels out
+        if (options->matchZPs) {
+            options->norm = psVectorAlloc(num, PS_TYPE_F32);
+            for (int i = 0; i < num; i++) {
+                if (!isfinite(trans->data.F32[i])) {
+                    continue;
+                }
+                psArray *sources = sourceLists->data[i]; // Sources of interest
+                float magCorr = airmassTerm - 2.5*log10(sumExpTime) - zp->data.F32[i] - trans->data.F32[i];
+                options->norm->data.F32[i] = magCorr;
+                psLogMsg("ppStack", PS_LOG_INFO, "Applying magnitude correction to image %d: %f\n",
+                         i, magCorr);
+
+                for (int j = 0; j < sources->n; j++) {
+                    pmSource *source = sources->data[j]; // Source of interest
+                    if (!source) {
+                        continue;
+                    }
+                    source->psfMag += magCorr;
+                }
+            }
+        }
+
+#ifdef TESTING
+        dumpMatches("source_mags.dat", num, matches, zp, trans);
+#endif
+        psFree(trans);
+
+#ifdef TESTING
+        // Double check: all transparencies should be zero
+        {
+            psArray *matches = pmSourceMatchSources(sourceLists, radius); // List of matches
+            if (!matches) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to match sources");
+                psFree(zp);
+                return false;
+            }
+            psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej,
+                                                   transThresh, starRej, starSys);
+            for (int i = 0; i < num; i++) {
+                fprintf(stderr, "Transparency of image %d: %f\n", i, trans->data.F32[i]);
+            }
+            psFree(trans);
+        }
+#endif
+    }
+
 
 #if 0
@@ -216 +273 @@
 #endif
 
-#ifdef TESTING
-    dumpMatches("source_mags.dat", num, matches, zp, trans);
-#endif
-
-    for (int i = 0; i < trans->n; i++) {
-        if (!isfinite(trans->data.F32[i])) {
-            inputMask->data.U8[i] |= PPSTACK_MASK_CAL;
-        }
-    }
-
-    // Save best matches SOMEWHERE for future photometry
-    // XXX this is a really poor output location; clean up the pmFPAfiles used in ppStack
-    pmCell *sourcesCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT");
-    psArray *sourcesBest = psArrayAllocEmpty(matches->n);
-
-    // XXX something of a hack: require at least 2 detections or the nominated fraction of the max possible
-    int minMatches = PS_MAX(2, fracMatch * num);// Minimum number of matches required
-    for (int i = 0; i < matches->n; i++) {
-        pmSourceMatch *match = matches->data[i]; // Match of interest
-        if (match->num < minMatches) {
-            continue;
-        }
-
-        // We need to grab a single instance of this source: just take the first available
-        int image = match->image->data.S32[0]; // Index of image
-        int index = match->index->data.S32[0]; // Index of source within image
-        psArray *sources = sourceLists->data[image]; // Sources for image
-        pmSource *source = sources->data[index]; // Source of interest
-
-        psArrayAdd(sourcesBest, sourcesBest->n, source);
-    }
-    psMetadataAdd(sourcesCell->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY | PS_META_REPLACE,
-                  "psphot sources", sourcesBest);
-    psLogMsg("ppStack", PS_LOG_INFO, "Selected %ld sources for photometry analysis", sourcesBest->n);
-    psFree(sourcesBest);
+    if (options->photometry) {
+        // Save best matches for future photometry
+        psArray *sourcesBest = options->sources = psArrayAllocEmpty(matches->n);
+
+        // XXX something of a hack: require at least 2 detections or the nominated fraction of the max possible
+        int minMatches = PS_MAX(2, fracMatch * num);// Minimum number of matches required
+        for (int i = 0; i < matches->n; i++) {
+            pmSourceMatch *match = matches->data[i]; // Match of interest
+            if (match->num < minMatches) {
+                continue;
+            }
+
+            // We need to grab a single instance of this source: just take the first available
+            int image = match->image->data.S32[0]; // Index of image
+            int index = match->index->data.S32[0]; // Index of source within image
+            psArray *sources = sourceLists->data[image]; // Sources for image
+            pmSource *source = sources->data[index]; // Source of interest
+
+            psArrayAdd(sourcesBest, sourcesBest->n, source);
+        }
+        psLogMsg("ppStack", PS_LOG_INFO, "Selected %ld sources for photometry analysis", sourcesBest->n);
+    }
 
     psFree(matches);
-
-    // M = m + c0 + c1 * airmass - 2.5log(t) + transparency
-    // Want sources to have m corresponding to airmass = 1 and t = sumExpTime and transparency = 0
-    // m_0 + c1 * airmass_0 + 2.5log(t_0) - trans_0 = m_1 + c1 * airmass_1 + 2.5log(t_1) - trans_1
-    // m_0 = m_1 + zp_1 + trans_1 - c1 * airmass_0 - 2.5log(t_0)
-    // We don't need to know the magnitude zero point for the filter, since it cancels out
-    options->norm = psVectorAlloc(num, PS_TYPE_F32);
-    for (int i = 0; i < num; i++) {
-        if (!isfinite(trans->data.F32[i])) {
-            continue;
-        }
-        psArray *sources = sourceLists->data[i]; // Sources of interest
-        float magCorr = airmassTerm - 2.5*log10(sumExpTime) - zp->data.F32[i] - trans->data.F32[i];
-        options->norm->data.F32[i] = magCorr;
-        psLogMsg("ppStack", PS_LOG_INFO, "Applying magnitude correction to image %d: %f\n", i, magCorr);
-
-        for (int j = 0; j < sources->n; j++) {
-            pmSource *source = sources->data[j]; // Source of interest
-            if (!source) {
-                continue;
-            }
-            source->psfMag += magCorr;
-        }
-    }
-    psFree(trans);
-
-#ifdef TESTING
-    // Double check: all transparencies should be zero
-    {
-        psArray *matches = pmSourceMatchSources(sourceLists, radius); // List of matches
-        if (!matches) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to match sources");
-            psFree(zp);
-            return false;
-        }
-        psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej,
-                                               transThresh, starRej, starSys);
-        for (int i = 0; i < num; i++) {
-            fprintf(stderr, "Transparency of image %d: %f\n", i, trans->data.F32[i]);
-        }
-        psFree(trans);
-    }
-#endif
 
     return true;