Changeset 5844 for trunk/psModules/src/objects/pmPSFtry.c

Timestamp:

Dec 23, 2005, 3:24:38 PM (20 years ago)

Author:

desonia

Message:

changes from eam_rel9_b1

File:

• trunk/psModules/src/objects/pmPSFtry.c (modified) (4 diffs)

trunk/psModules/src/objects/pmPSFtry.c

@@ -5 +5 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-12 21:14:38 $
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-24 01:24:32 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -176 +176 @@
 
     // XXX this function wants aperture radius for pmSourcePhotometry
-    pmPSFtryMetric (psfTry, RADIUS);
+    pmPSFtryMetric_Alt (psfTry, RADIUS);
     psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n",
               modelName,
@@ -288 +288 @@
 
     // 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
@@ -314 +316 @@
     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);
+    fprintf (stderr, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
+
+    try
+        ->psf->ApResid = poly->coeff[0];
+    try
+        ->psf->dApResid = stats->sampleStdev;
+    try
+        ->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
+
+    fprintf (stderr, "*******************************************************************************\n");
+
+    FILE *f;
+    f = fopen ("apresid.dat", "w");
+    if (f == NULL)
+        psAbort ("pmPSFtry", "can't open output file");
+
+    for (int i = 0; i < try
+                ->sources->n; i++) {
+            fprintf (f, "%3d %8.4f %12.5e %8.4f %3d\n", i, try
+                         ->fitMag->data.F64[i], rflux->data.F64[i], try
+                             ->metric->data.F64[i], try
+                                 ->mask->data.U8[i]);
+        }
+    fclose (f);
+
+    psFree (rflux);
+    psFree (poly);
+    psFree (stats);
+
+    // psFree (daFit);
+    // psFree (daResid);
+
+    return true;
+}