Changeset 36031

Timestamp:

Aug 25, 2013, 12:30:55 PM (13 years ago)

Author:

eugene

Message:

loop over Reff and Sindex

File:

• branches/eam_branches/ipp-20130711/psphot/src/psphotSourceFits.c (modified) (3 diffs)

branches/eam_branches/ipp-20130711/psphot/src/psphotSourceFits.c

@@ -760 +760 @@
 # define N_INDEX_GUESS_INV 4
 
+float reffGuess[] = {3.0, 10.0, 30.0};
+# define N_REFF_GUESS 3
+
 // A sersic model is very sensitive to the index.  attempt to find the index first by grid search in just the index
 // for a sersic model, attempt to fit just the index and normalization with a modest number of iterations
@@ -787 +790 @@
     psEllipseAxes momentAxes = psEllipseMomentsToAxes (moments, 20.0);
 
-    psEllipseAxes guessAxes;
-    guessAxes.major = source->moments->Mrf;
-    guessAxes.minor = (momentAxes.minor / momentAxes.minor) * guessAxes.major;
-    guessAxes.theta = momentAxes.theta;
-
-    if (!isfinite(guessAxes.major)) return false;
-    if (!isfinite(guessAxes.minor)) return false;
-    if (!isfinite(guessAxes.theta)) return false;
-
-    // convert the major,minor,theta to shape parameters for an Reff-like model
-    pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true);
+    if (0) {
+        psEllipseAxes guessAxes;
+        guessAxes.major = source->moments->Mrf;
+        guessAxes.minor = (momentAxes.minor / momentAxes.minor) * guessAxes.major;
+        guessAxes.theta = momentAxes.theta;
+
+        if (!isfinite(guessAxes.major)) return false;
+        if (!isfinite(guessAxes.minor)) return false;
+        if (!isfinite(guessAxes.theta)) return false;
+
+        // convert the major,minor,theta to shape parameters for an Reff-like model
+        pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true);
+    }
 
     // set the model position
@@ -813 +818 @@
     float iMin = NAN;
     float sMin = NAN;
+    float rMin = NAN;
 
     // loop over index and Reff, keeping the ARatio and Theta constant?
     // loop over index guesses and find the best fit
-    for (int i = 0; i < N_INDEX_GUESS_INV; i++) {
-        PAR[PM_PAR_7] = indexGuessInv[i];
-
-        // generated the modelFlux
-        pmPCMMakeModel (source, pcm->modelConv, maskVal, psfSize);
+    for (int j = 0; j < N_REFF_GUESS; j++) {
+
+        psEllipseAxes guessAxes;
+        guessAxes.major = reffGuess[j];
+        guessAxes.minor = (momentAxes.minor / momentAxes.minor) * guessAxes.major;
+        guessAxes.theta = momentAxes.theta;
+
+        if (!isfinite(guessAxes.major)) return false;
+        if (!isfinite(guessAxes.minor)) return false;
+        if (!isfinite(guessAxes.theta)) return false;
+
+        // convert the major,minor,theta to shape parameters for an Reff-like model
+        pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true);
+
+        for (int i = 0; i < N_INDEX_GUESS_INV; i++) {
+            PAR[PM_PAR_7] = indexGuessInv[i];
+
+            // generated the modelFlux
+            pmPCMMakeModel (source, pcm->modelConv, maskVal, psfSize);
         
-        float YY = 0.0;
-        float YM = 0.0;
-        float MM = 0.0;
-        bool usePoisson = false;
-
-        for (int iy = 0; iy < source->pixels->numRows; iy++) {
-            for (int ix = 0; ix < source->pixels->numCols; ix++) {
-                // skip masked points
-                if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) {
-                    continue;
+            float YY = 0.0;
+            float YM = 0.0;
+            float MM = 0.0;
+            bool usePoisson = false;
+
+            for (int iy = 0; iy < source->pixels->numRows; iy++) {
+                for (int ix = 0; ix < source->pixels->numCols; ix++) {
+                    // skip masked points
+                    if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) {
+                        continue;
+                    }
+                    // skip zero-variance points
+                    if (source->variance->data.F32[iy][ix] == 0) {
+                        continue;
+                    }
+                    // skip nan value points
+                    if (!isfinite(source->pixels->data.F32[iy][ix])) {
+                        continue;
+                    }
+
+                    float fy = source->pixels->data.F32[iy][ix];
+                    float fm = source->modelFlux->data.F32[iy][ix];
+                    float wt = (usePoisson) ? 1.0 / source->variance->data.F32[iy][ix] : 1.0;
+
+                    YY += PS_SQR(fy) * wt;
+                    YM += fm * fy * wt;
+                    MM += PS_SQR(fm) * wt;
                 }
-                // skip zero-variance points
-                if (source->variance->data.F32[iy][ix] == 0) {
-                    continue;
-                }
-                // skip nan value points
-                if (!isfinite(source->pixels->data.F32[iy][ix])) {
-                    continue;
-                }
-
-                float fy = source->pixels->data.F32[iy][ix];
-                float fm = source->modelFlux->data.F32[iy][ix];
-                float wt = (usePoisson) ? 1.0 / source->variance->data.F32[iy][ix] : 1.0;
-
-                YY += PS_SQR(fy) * wt;
-                YM += fm * fy * wt;
-                MM += PS_SQR(fm) * wt;
             }
-        }
-
-        float Io = YM / MM;
-        float Chisq = YY - 2 * Io * YM + Io * Io * MM;
-        if ((i == 0) || (Chisq < xMin)) {
-            xMin = Chisq;
-            iMin = Io;
-            sMin = indexGuessInv[i];
-        }
-        fprintf (stderr, "%d | %f %f %f | %f %f %f", i, indexGuessInv[i], Io, Chisq, sMin, iMin, xMin);
-        fprintf (stderr, "\n");
+
+            float Io = YM / MM;
+            float Chisq = YY - 2 * Io * YM + Io * Io * MM;
+            if (isnan(xMin) || (Chisq < xMin)) {
+                xMin = Chisq;
+                iMin = Io;
+                sMin = indexGuessInv[i];
+                rMin = reffGuess[j];
+            }
+            fprintf (stderr, "%d | %f %f %f | %f %f %f", i, indexGuessInv[i], reffGuess[j], Io, Chisq, sMin, rMin, iMin, xMin);
+            fprintf (stderr, "\n");
+        }
+    }
+
+    {
+        psEllipseAxes guessAxes;
+        guessAxes.major = rMin;
+        guessAxes.minor = (momentAxes.minor / momentAxes.minor) * guessAxes.major;
+        guessAxes.theta = momentAxes.theta;
+
+        if (!isfinite(guessAxes.major)) return false;
+        if (!isfinite(guessAxes.minor)) return false;
+        if (!isfinite(guessAxes.theta)) return false;
+
+        // convert the major,minor,theta to shape parameters for an Reff-like model
+        pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true);
     }