Changeset 25433 for branches/eam_branches/20090715/psphot/src/psphotEllipticalContour.c

Timestamp:

Sep 17, 2009, 2:26:32 PM (17 years ago)

Author:

eugene

Message:

change radius for extended sources to use footprint; clean up some of the visualizations; plug some leaks

File:

• branches/eam_branches/20090715/psphot/src/psphotEllipticalContour.c (modified) (6 diffs)

branches/eam_branches/20090715/psphot/src/psphotEllipticalContour.c

@@ -1 +1 @@
 # include "psphotInternal.h"
-
-// XXX consistency : theta in radians here and in calling functions
 
 // model parameters
@@ -7 +5 @@
 psF32 psphotEllipticalContourFunc (psVector *deriv, const psVector *params, const psVector *coord);
 
-bool psphotEllipticalContour (pmPetrosian *petrosian) {
+bool psphotEllipticalContour (pmSource *source, pmPetrosian *petrosian) {
 
     // use LMM to fit theta vs radius to an ellipse
@@ -19 +17 @@
     // arrays to hold the data to be fitted
     // we fit x and y vs theta in separate passes.
-    psArray *x = psArrayAlloc(2*radius->n);
-    psVector *y = psVectorAlloc(2*radius->n, PS_TYPE_F32);
-    psVector *yErr = psVectorAlloc(2*radius->n, PS_TYPE_F32);
+    psArray *x = psArrayAllocEmpty(2*radius->n);
+    psVector *y = psVectorAllocEmpty(2*radius->n, PS_TYPE_F32);
+    psVector *yErr = psVectorAllocEmpty(2*radius->n, PS_TYPE_F32);
 
     int n = 0;
     for (int i = 0; i < radius->n; i++) {
+        if (!isfinite(radius->data.F32[i])) continue;
 
         psVector *coord = NULL;
@@ -50 +49 @@
         Rmax = MAX (Rmax, radius->data.F32[i]);
     }   
-    assert (x->n == n);
-    assert (y->n == n);
+    x->n = n;
+    y->n = n;
+    yErr->n = n;
+
+    if (n < 4) {
+        psFree (x);
+        psFree (y);
+        psFree (yErr);
+        return false;
+    }
 
     psVector *params = psVectorAlloc (3, PS_TYPE_F32);
@@ -61 +68 @@
 
     // XXX for now, no parameter masks, skip checkLimits
+    // XXX might help to add a limit to the angle or re-parameterize the ellipse in terms of Rxx, Ryy, Rxy
     // constraint->checkLimits = psastroModelBoresiteLimits;
 
@@ -84 +92 @@
     }
 
-    fprintf (stderr, "# fitted values:\n");
-    fprintf (stderr, "Po:  %f\n", params->data.F32[PAR_PHI]*PS_DEG_RAD);
-    fprintf (stderr, "Ep:  %f\n", params->data.F32[PAR_EPSILON]);
-    fprintf (stderr, "Rm:  %f\n", params->data.F32[PAR_RMIN]);
-
+    psTrace ("psphot", 4, "# fitted values:\n");
+    psTrace ("psphot", 4, "Phi:   %f\n", petrosian->axes.theta*PS_DEG_RAD);
+    psTrace ("psphot", 4, "Rmaj:  %f\n", petrosian->axes.major);
+    psTrace ("psphot", 4, "Rmin:  %f\n", petrosian->axes.minor);
+    
     // show the results
-    psphotPetrosianVisualEllipticalContour (petrosian);
+    // psphotPetrosianVisualEllipticalContour (petrosian);
 
     psFree (x);
     psFree (y);
     psFree (yErr);
+    psFree (params);
+    psFree (covar);
+    psFree (myMin);
+    psFree (constraint);
 
     return true;