Changeset 5725 for branches/eam_rel9_b0/psModules/src/objects/pmPSFtry.c

Timestamp:

Dec 7, 2005, 9:16:37 AM (21 years ago)

Author:

magnier

Message:

fixed order for psPolynomial1DAlloc arguments
called correct for of psVectorClipFitPolynomial1D

added pmPSFtryMetric_Alt, an alternative (preferred) measurement of
the PSF metric.

File:

• branches/eam_rel9_b0/psModules/src/objects/pmPSFtry.c (modified) (3 diffs)

branches/eam_rel9_b0/psModules/src/objects/pmPSFtry.c

@@ -183 +183 @@
 
     // XXX this function wants aperture radius for pmSourcePhotometry
-    pmPSFtryMetric (try
-                    , RADIUS);
+    pmPSFtryMetric_Alt (try
+                        , RADIUS);
     psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n",
               modelName, try
@@ -305 +305 @@
 
     // linear fit to rfBin, daBin
-    psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1);
+    psPolynomial1D *poly = psPolynomial1DAlloc(1, PS_POLYNOMIAL_ORD);
+    psStats *fitstat = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+    poly = psVectorClipFitPolynomial1D (poly, fitstat, maskB, 1, daBin, NULL, rfBin);
 
     // XXX EAM : this is the intended API (cycle 7? cycle 8?)
-    poly = psVectorFitPolynomial1D(poly, maskB, 1, daBin, NULL, rfBin);
+    // poly = psVectorFitPolynomial1D(poly, maskB, 1, daBin, NULL, rfBin);
 
     // XXX EAM : replace this when the above version is implemented
@@ -334 +336 @@
     psFree (poly);
     psFree (stats);
+    psFree (fitstat);
 
     return true;
 }
+
+bool pmPSFtryMetric_Alt (pmPSFtry *try
+                         , float RADIUS)
+{
+
+    // the measured (aperture - fit) magnitudes (dA == try->metric)
+    //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
+    //     dA = dAo + dsky/flux
+    //   where flux is the flux of the star
+    // we fit this trend to find the infinite flux aperture correction (dAo),
+    //   the nominal sky bias (dsky), and the error on dAo
+    // the values of dA are contaminated by stars with close neighbors in the aperture
+    //   we use an outlier rejection to avoid this bias
+
+    // rflux = ten(0.4*fitMag);
+    psVector *rflux = psVectorAlloc (try
+                                     ->sources->n, PS_TYPE_F64);
+    for (int i = 0; i < try
+                ->sources->n; i++) {
+            if (try
+                    ->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
+            rflux->data.F64[i] = pow(10.0, 0.4*try
+                                     ->fitMag->data.F64[i]);
+        }
+
+    // XXX EAM : try 3hi/1lo sigma clipping on the rflux vs metric fit
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+
+    // XXX EAM
+    stats->min = 1.0;
+    stats->max = 3.0;
+    stats->clipIter = 3;
+
+    // linear clipped fit to rfBin, daBin
+    psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
+    poly = psVectorClipFitPolynomial1D (poly, stats, try
+                                            ->mask, PSFTRY_MASK_ALL, try
+                                                ->metric, NULL, rflux);
+
+    try
+        ->psf->ApResid = poly->coeff[0];
+    try
+        ->psf->dApResid = stats->sampleStdev;
+    try
+        ->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
+
+    psFree (rflux);
+    psFree (poly);
+    psFree (stats);
+
+    // psFree (daFit);
+    // psFree (daResid);
+
+    return true;
+}