Changeset 36325 for branches/eam_branches/ipp-20130904/psphot/src/psphotSourceFits.c

Timestamp:

Nov 28, 2013, 1:25:15 PM (13 years ago)

Author:

eugene

Message:

fix up grid search for sersic index

File:

• branches/eam_branches/ipp-20130904/psphot/src/psphotSourceFits.c (modified) (9 diffs)

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

@@ -653 +653 @@
     if (TIMING) { t4 = psTimerMark ("psphotFitPCM"); }
 
-    // psphotPCMfitCheckSize (pcm, source, maskVal, psfSize);
+    if (options.isInteractive) psphotPCMfitCheckSize (pcm, source, maskVal, psfSize);
     // if (pcm->modelConv->nIter == fitOptions->nIter) {
     //  psphotPCMfitRetry (pcm, source, &options, maskVal, markVal, psfSize);
@@ -791 +791 @@
 
 // 0.5 / n for (0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0)
-float indexGuessInv[] = {1.00, 0.50, 0.333, 0.25, 1.666, 0.125, 0.10, 0.0833};
+float indexGuessInv[] = {1.00, 0.50, 0.333, 0.25, 0.166, 0.125, 0.10, 0.0833};
 float indexGuessR1q[] = {1.06, 1.19, 1.335, 1.48, 1.840, 2.290, 2.84, 3.5300};
 # define N_INDEX_GUESS_INV 8
@@ -902 +902 @@
                 rMin = reffGuess[j] / indexGuessR1q[i];
             }
-            // fprintf (stderr, "%d | %f %f %f %f | %f %f %f %f", i, indexGuessInv[i], reffGuess[j], Io, Chisq, sMin, rMin, iMin, xMin);
-            // fprintf (stderr, "\n");
+            fprintf (stderr, "%d | %f %f %f %f | %f %f %f %f", i, indexGuessInv[i], reffGuess[j], Io, Chisq, sMin, rMin, iMin, xMin);
+            fprintf (stderr, "\n");
         }
     }
@@ -1006 +1006 @@
     psVector *Sidx = psVectorAllocEmpty (16, PS_TYPE_F32);
 
-    PAR[PM_PAR_7] = indexGuessInv[nStart];
+    float Sm = NAN, Sp = NAN, So = NAN;
+    if (nStart == 0) {
+        Sm = indexGuessInv[nStart];
+        So = 0.5*(indexGuessInv[nStart + 1] + indexGuessInv[nStart]);
+        Sp = indexGuessInv[nStart + 1];
+    } else if (nStart == N_INDEX_GUESS_INV - 1) {
+        Sp = indexGuessInv[nStart];
+        So = 0.5*(indexGuessInv[nStart - 1] + indexGuessInv[nStart]);
+        Sm = indexGuessInv[nStart - 1];
+    } else {
+        Sm = 0.5*(indexGuessInv[nStart - 1] + indexGuessInv[nStart]);
+        So = indexGuessInv[nStart];
+        Sp = 0.5*(indexGuessInv[nStart + 1] + indexGuessInv[nStart]);
+    }
+    
+    PAR[PM_PAR_7] = Sm;
     pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
     if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]);
@@ -1012 +1027 @@
     psVectorAppend (chi2, model->chisqNorm);
 
-    PAR[PM_PAR_7] = (nStart < N_INDEX_GUESS_INV - 1) ? 0.5*(indexGuessInv[nStart + 1] + indexGuessInv[nStart]) : indexGuessInv[nStart - 1];
+    PAR[PM_PAR_7] = So;
     pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
     if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]);
@@ -1018 +1033 @@
     psVectorAppend (chi2, model->chisqNorm);
 
-    PAR[PM_PAR_7] = (nStart > 0) ? 0.5*(indexGuessInv[nStart - 1] + indexGuessInv[nStart]) : indexGuessInv[nStart + 1];
+    PAR[PM_PAR_7] = Sp;
     pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
     if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]);
@@ -1036 +1051 @@
     float Smin = -0.5 * poly->coeff[1] / poly->coeff[2] / 100.0;
 
-    // constrain Smin to be in a valid range (1.0 - 0.1, corresponding to 0.5 (Gauss) to 5.0 (slightly peakier than Dev)
-    Smin = PS_MAX(PS_MIN(Smin, 1.0), 0.1);
+    // constrain Smin to be in a valid range: allow the fitted range to go a bit beyond the 3 trial points, but no further
+    float Smx = Sm - 0.25*(So - Sm);
+    float Spx = Sp + 0.25*(Sp - So);
+    Smin = PS_MAX(PS_MIN(Smin, Smx), Spx);
     PAR[PM_PAR_7] = Smin;
 
@@ -1249 +1266 @@
 
     // loop over Reff, keeping the ARatio and Theta constant
-    for (int j = -20; j < 21; j++) {
-
-        float dref = j * 0.02;
+    for (int j = -4; j <= 4; j++) {
+
+        float dref = j * 0.01;
 
         psEllipseAxes guessAxes;
@@ -1306 +1323 @@
             rMin = dref;
         }
-        fprintf (stderr, "%d | %f %f %f | %f %f %f\n", j, dref, Io, Chisq, rMin, iMin, xMin);
+        // fprintf (stderr, "%d | %f %f %f | %f %f %f\n", j, dref, Io, Chisq, rMin, iMin, xMin);
     }