Changeset 31438 for branches/eam_branches/ipp-20110404/psModules/src/objects/pmGrowthCurveGenerate.c

Timestamp:

May 4, 2011, 4:44:43 PM (15 years ago)

Author:

eugene

Message:

explicitly include the reference aperture in the growth curve; add some debugging output to pmGrowthCurveGenerate.c; use the reference aperture magnitude when using sources for the growth curve; extend the model integration to 150 sigmas^{2 (not just 50); re-generate modelFlux pixels if at all deviant from source}

File:

• branches/eam_branches/ipp-20110404/psModules/src/objects/pmGrowthCurveGenerate.c (modified) (5 diffs)

branches/eam_branches/ipp-20110404/psModules/src/objects/pmGrowthCurveGenerate.c

@@ -226 +226 @@
 }
 
+# define DEBUG 0
+# if (DEBUG)
+static FILE *fgr = NULL;
+# endif
+
 // we generate the growth curve for the center of the image with the specified psf model
 bool pmGrowthCurveGenerateFromSources (pmReadout *readout, pmPSF *psf, psArray *sources, bool INTERPOLATE_AP, psImageMaskType maskVal, psImageMaskType markVal)
@@ -240 +245 @@
     psArray *growths = psArrayAllocEmpty (100);
 
+# if (DEBUG)
+    fgr = fopen ("growth.mags.dat", "w");
+# endif
+
     for (int i = 0; i < sources->n; i++) {
 
@@ -254 +263 @@
     psAssert (growths->n, "cannot build growth curve (no valid PSF stars?)");
 
-# if (0)
-    FILE *f = fopen ("growth.dat", "w");
-    assert (f);
-
-    for (int j = 0; j < psf->growth->radius->n; j++) {
-
-        fprintf (f, "%f : ", psf->growth->radius->data.F32[j]);
-
-        // XXX dump the growth curves:
-        for (int i = 0; i < growths->n; i++) {
-            
-            pmGrowthCurve *growth = growths->data[i];
-            if (!growth) continue;
-
-            fprintf (f, "%8.4f ", growth->apMag->data.F32[j]);
-        }
-        fprintf (f, "\n");
-    }
-    fclose (f);
+# if (DEBUG)
+    fclose (fgr);
 # endif
 
@@ -290 +282 @@
             pmGrowthCurve *growth = growths->data[j];
             if (!isfinite(growth->apMag->data.F32[i])) continue;
-            psVectorAppend (values, growth->apMag->data.F32[i] - growth->fitMag);
+            psVectorAppend (values, growth->apMag->data.F32[i] - growth->refMag);
         }
         if (values->n == 0) {
@@ -354 +346 @@
         growth->apMag->data.F32[i] = source->apMag;
     }
+    psAssert(growth->refBin >= 0, "invalid growth reference bin");
+    psAssert(growth->refBin < growth->apMag->n, "invalid growth reference bin");
+    growth->refMag = growth->apMag->data.F32[growth->refBin];
+
+    // Loop over the range of radii
+# if (DEBUG)
+    for (int i = 0; i < growth->radius->n; i++) {
+        fprintf (fgr, "%f %f  %f %f %f %f\n", xc, yc, growth->radius->data.F32[i], growth->apMag->data.F32[i], growth->fitMag, growth->refMag);
+    }
+# endif
+
     return growth;
 }
-
-# if (0)
-    // solve for the best growth and the offsets per star to minimize the scatter
-
-    // generate the per-star and per-radius sums
-    int Nradius = psf->growth->radius->n;
-    int Nstar = growths->n;
-
-    psVector *Mstar   = psVectorAlloc (Nstar, PS_DATA_F32);
-    psVector *Mradius = psVectorAlloc (Nradius, PS_DATA_F32);
-
-    psVectorInit (Mstar, 0.0);
-    psVectorInit (Mradius, 0.0);
-
-    for (int i = 0; i < Nstar; i++) {
-        pmGrowthCurve *growth = growths->data[i];
-        for (int j = 0; j < Nradius; j++) {
-            Mstar->data.F32[i] += growth->apMag->data.F32[j];
-            Mradius->data.F32[j] += growth->apMag->data.F32[j];
-        }
-    }
-
-    psImage *A = psImageAlloc(Nstar + Nradius, Nstar + Nradius, PS_DATA_F32);
-    psVector *B = psVectorAlloc(Nstar + Nradius, PS_DATA_F32);
-
-    psImageInit (A, 0.0);
-    psVectorInit (B, 0.0);
-
-    for (int i = 0; i < Nstar; i++) {
-        A->data.F32[i][i] = Nradius;
-        B->data.F32[i] = Mstar->data.F32[i];
-    }
-    for (int i = 0; i < Nradius; i++) {
-        int j = i + Nstar;
-        A->data.F32[j][j] = Nstar;
-        B->data.F32[j] = Mradius->data.F32[i];
-    }
-    
-    if (!psMatrixGJSolve(A, B)) {
-        psAbort("GJ failure");
-    }
-
-    for (int i = 0; i < Nradius; i++) {
-        psf->growth->apMag->data.F32[i] = B->data.F32[Nstar+i];
-    }
-
-    // XXX this analysis produces a biased growth curve: points with small radius are more
-    // likely to err on the side of being too faint (bacause of mis-aligned aperture), while on
-    // the bright side the tend to be too bright (because of contaminating neighbors).  Bright
-    // stars are less likely to be biased, 
-
-    // XXX at this point, we could iterate and exclude some of the stars with low data quality
-    // XXX we could also weight the above by flux or something
-
-# endif