Changeset 35768

Timestamp:

Jul 3, 2013, 2:37:22 PM (13 years ago)

Author:

eugene

Message:

deprecate KiiOpen,KiiClose (now KapaOpen,etc); major rework of psEllipse translations : use common functions pmModelAxesToParams and pmModelParamsToAxes ; use new convergence method in pmPCM_MinimizeChisq; add convergence crerition options to psMinimization; threaded versions of pmPSFtryFitEXT and pmPSFtryFitPSF

Location:

trunk/psModules

Files:

• . (modified) (1 prop)
• src/astrom/pmAstrometryUtils.c (modified) (1 diff)
• src/astrom/pmAstrometryVisual.c (modified) (1 diff)
• src/camera/pmReadoutFake.c (modified) (2 diffs)
• src/extras/Makefile.am (modified) (2 diffs)
• src/extras/pmThreadTools.c (copied) (copied from branches/eam_branches/ipp-20130509/psModules/src/extras/pmThreadTools.c )
• src/extras/pmThreadTools.h (copied) (copied from branches/eam_branches/ipp-20130509/psModules/src/extras/pmThreadTools.h )
• src/imcombine/pmStackVisual.c (modified) (1 diff)
• src/imcombine/pmSubtractionVisual.c (modified) (1 diff)
• src/objects/models/pmModel_DEV.c (modified) (9 diffs)
• src/objects/models/pmModel_EXP.c (modified) (8 diffs)
• src/objects/models/pmModel_GAUSS.c (modified) (7 diffs)
• src/objects/models/pmModel_PGAUSS.c (modified) (7 diffs)
• src/objects/models/pmModel_PS1_V1.c (modified) (8 diffs)
• src/objects/models/pmModel_QGAUSS.c (modified) (12 diffs)
• src/objects/models/pmModel_RGAUSS.c (modified) (7 diffs)
• src/objects/models/pmModel_SERSIC.c (modified) (9 diffs)
• src/objects/models/pmModel_TRAIL.c (modified) (2 diffs)
• src/objects/pmModelUtils.c (modified) (2 diffs)
• src/objects/pmModelUtils.h (modified) (1 diff)
• src/objects/pmPCM_MinimizeChisq.c (modified) (6 diffs)
• src/objects/pmPCMdata.c (modified) (2 diffs)
• src/objects/pmPSF.c (modified) (6 diffs)
• src/objects/pmPSF.h (modified) (1 diff)
• src/objects/pmPSFtry.h (modified) (1 diff)
• src/objects/pmPSFtryFitEXT.c (modified) (3 diffs)
• src/objects/pmPSFtryFitPSF.c (modified) (5 diffs)
• src/objects/pmPSFtryMakePSF.c (modified) (2 diffs)
• src/objects/pmSource.c (modified) (2 diffs)
• src/objects/pmSourceFitModel.c (modified) (4 diffs)
• src/objects/pmSourceFitModel.h (modified) (1 diff)
• src/objects/pmSourceFitPCM.c (modified) (4 diffs)
• src/objects/pmSourceFitSet.c (modified) (1 diff)
• src/objects/pmSourceIO_CMF.c.in (modified) (1 diff)
• src/objects/pmSourceIO_CMP.c (modified) (1 diff)
• src/objects/pmSourceIO_OBJ.c (modified) (1 diff)
• src/objects/pmSourceIO_PS1_CAL_0.c (modified) (2 diffs)
• src/objects/pmSourceIO_PS1_DEV_0.c (modified) (1 diff)
• src/objects/pmSourceIO_PS1_DEV_1.c (modified) (2 diffs)
• src/objects/pmSourceIO_SMPDATA.c (modified) (1 diff)
• src/objects/pmSourceIO_SX.c (modified) (1 diff)
• src/objects/pmSourceOutputs.c (modified) (1 diff)
• src/objects/pmSourcePlotPSFModel.c (modified) (1 diff)
• src/objects/pmSourceVisual.c (modified) (1 diff)
• src/psmodules.h (modified) (1 diff)

trunk/psModules/src/astrom/pmAstrometryUtils.c

@@ -78 +78 @@
             Alpha->data.F32[1][1] = psPolynomial2DEval (YdY, Xo, Yo);
 
+            // Beta = (-L,-M) for the current guess
             Beta->data.F32[0] = -1.0 * psPolynomial2DEval (trans->x, Xo, Yo);
             Beta->data.F32[1] = -1.0 * psPolynomial2DEval (trans->y, Xo, Yo);
-
+            // fprintf (stderr, " Beta: %f %f\n", Beta->data.F32[0], Beta->data.F32[1]);
+
+            // since we want (L,M) = (0,0), Beta is also the offset
             if (!psMatrixGJSolve (Alpha, Beta)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to solve for center.");

trunk/psModules/src/astrom/pmAstrometryVisual.c

@@ -51 +51 @@
 bool pmAstromVisualClose(void)
 {
-    if (kapa1 != -1) KiiClose(kapa1);
-    if (kapa2 != -1) KiiClose(kapa2);
+    if (kapa1 != -1) KapaClose(kapa1);
+    if (kapa2 != -1) KapaClose(kapa2);
     return true;
 }

trunk/psModules/src/camera/pmReadoutFake.c

@@ -35 +35 @@
 #include "pmReadoutFake.h"
 
-#define MAX_AXIS_RATIO 20.0             // Maximum axis ratio for PSF model
+// XXX this is now hard-wired in pmModelParamsToAxes
+// #define MAX_AXIS_RATIO 20.0             // Maximum axis ratio for PSF model
+
 #define MODEL_MASK (PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE | \
                     PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_LIMITS) // Mask to apply to models
@@ -52 +54 @@
 
     psF32 *params = model->params->data.F32; // Model parameters
-    psEllipseAxes axes = pmPSF_ModelToAxes(params, MAX_AXIS_RATIO, model->type); // Ellipse axes
+    psEllipseAxes axes = pmPSF_ModelToAxes(params, model->type); // Ellipse axes
     // Curiously, the minor axis can be larger than the major axis, so need to check.
     if (axes.major >= axes.minor) {

trunk/psModules/src/extras/Makefile.am

@@ -10 +10 @@
         pmVisual.c \
         pmVisualUtils.c \
+        pmThreadTools.c \
         ippStages.c \
         pmCensor.c
@@ -20 +21 @@
         pmVisual.h \
         pmVisualUtils.h \
+        pmThreadTools.h \
         ippDiffMode.h \
         ippStages.h \

trunk/psModules/src/imcombine/pmStackVisual.c

@@ -36 +36 @@
 {
     if(kapa != -1)
-        KiiClose(kapa);
+        KapaClose(kapa);
     if(kapa2 != -1)
-        KiiClose(kapa2);
+        KapaClose(kapa2);
     return true;
 }

trunk/psModules/src/imcombine/pmSubtractionVisual.c

@@ -56 +56 @@
 bool pmSubtractionVisualClose(void)
 {
-    if(kapa1 != -1) KiiClose(kapa1);
-    if(kapa2 != -1) KiiClose(kapa2);
+    if(kapa1 != -1) KapaClose(kapa1);
+    if(kapa2 != -1) KapaClose(kapa2);
     return true;
 }

trunk/psModules/src/objects/models/pmModel_DEV.c

@@ -123 +123 @@
 
         // first, use Rmajor and index to find the central pixel flux (fraction of total flux)
-        psEllipseShape shape;
-
-        shape.sx  = PAR[PM_PAR_SXX];
-        shape.sy  = PAR[PM_PAR_SYY];
-        shape.sxy = PAR[PM_PAR_SXY];
-
-        // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio
-        psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+        psEllipseAxes axes;
+        pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
 
         // get the central pixel flux from the lookup table
@@ -238 +232 @@
         // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
         // angle and let f2,f1 fight it out
-        q2 = 0.5*sqrtf(q1);
+        // NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy
+        q2 = 2.0*0.5*sqrtf(q1);
     }
 
@@ -303 +298 @@
 
     // set the shape parameters
-    // XXX adjust this?
-    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
+    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, true)) {
       return false;
     }
@@ -330 +324 @@
 psF64 PM_MODEL_FLUX (const psVector *params)
 {
-    psEllipseShape shape;
-
     psF32 *PAR = params->data.F32;
 
-    shape.sx  = PAR[PM_PAR_SXX];
-    shape.sy  = PAR[PM_PAR_SYY];
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    // for a non-circular DeVaucouleur, the flux of the Rmajor equivalent is scaled by the AspectRatio
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
     float AspectRatio = axes.minor / axes.major;
 
@@ -359 +347 @@
 psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
 {
-    psEllipseShape shape;
-
     psF32 *PAR = params->data.F32;
 
@@ -370 +356 @@
         return (1.0);
 
-    shape.sx  = PAR[PM_PAR_SXX];
-    shape.sy  = PAR[PM_PAR_SYY];
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
 
     // f = Io exp(-z^n) -> z^n = ln(Io/f)
@@ -382 +365 @@
     psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", 
               PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY]);
-
     return (radius);
 }
@@ -407 +389 @@
     // the 2D PSF model fits polarization terms (E0,E1,E2)
     // convert to shape terms (SXX,SYY,SXY)
-    if (!pmPSF_FitToModel (out, 0.1)) {
+    bool useReff = pmModelUseReff (modelPSF->type);
+    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
         psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
         return false;
@@ -460 +443 @@
     // convert to shape terms (SXX,SYY,SXY)
     // XXX user-defined value for limit?
-    if (!pmPSF_FitToModel (PAR, 0.1)) {
+    bool useReff = pmModelUseReff (model->type);
+    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
         psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
         return false;

trunk/psModules/src/objects/models/pmModel_EXP.c

@@ -115 +115 @@
 
         // first, use Rmajor and index to find the central pixel flux (fraction of total flux)
-        psEllipseShape shape;
-
-        shape.sx  = PAR[PM_PAR_SXX];
-        shape.sy  = PAR[PM_PAR_SYY];
-        shape.sxy = PAR[PM_PAR_SXY];
-
-        // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio
-        psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+        psEllipseAxes axes;
+        pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
 
         // get the central pixel flux from the lookup table
@@ -230 +224 @@
         // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
         // angle and let f2,f1 fight it out
-        q2 = 0.5*sqrtf(q1);
+        // NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy
+        q2 = 2.0*0.5*sqrtf(q1);
     }
 
@@ -295 +290 @@
 
     // set the shape parameters
-    // XXX adjust this?
-    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
+    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, true)) {
       return false;
     }
@@ -316 +310 @@
 psF64 PM_MODEL_FLUX (const psVector *params)
 {
-    psEllipseShape shape;
-
     psF32 *PAR = params->data.F32;
 
-    shape.sx  = PAR[PM_PAR_SXX];
-    shape.sy  = PAR[PM_PAR_SYY];
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    // for a non-circular Exponential, the flux of the Rmajor equivalent is scaled by the AspectRatio
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
     float AspectRatio = axes.minor / axes.major;
 
@@ -345 +333 @@
 psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
 {
-    psEllipseShape shape;
-
     psF32 *PAR = params->data.F32;
 
@@ -356 +342 @@
         return (1.0);
 
-    shape.sx  = PAR[PM_PAR_SXX];
-    shape.sy  = PAR[PM_PAR_SYY];
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
 
     // f = Io exp(-sqrt(z)) -> sqrt(z) = ln(Io/f)
@@ -392 +375 @@
     // the 2D PSF model fits polarization terms (E0,E1,E2)
     // convert to shape terms (SXX,SYY,SXY)
-    if (!pmPSF_FitToModel (out, 0.1)) {
+    bool useReff = pmModelUseReff (modelPSF->type);
+    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
         psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
         return false;
@@ -445 +429 @@
     // convert to shape terms (SXX,SYY,SXY)
     // XXX user-defined value for limit?
-    if (!pmPSF_FitToModel (PAR, 0.1)) {
+    bool useReff = pmModelUseReff (model->type);
+    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
         psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
         return false;

trunk/psModules/src/objects/models/pmModel_GAUSS.c

@@ -129 +129 @@
     float q2 = NAN;
     if (nParam == PM_PAR_SXY) {
-        float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
-        float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
+        // NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy
+        float f1 = 2.0 / PS_SQR(params[PM_PAR_SYY]) + 2.0 / PS_SQR(params[PM_PAR_SXX]);
+        float f2 = 2.0 / PS_SQR(params[PM_PAR_SYY]) - 2.0 / PS_SQR(params[PM_PAR_SXX]);
         float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
         q1 = (q1 < 0.0) ? 0.0 : q1;
@@ -200 +201 @@
 
     // set the shape parameters
-    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
+    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, false)) {
       return false;
     }
@@ -219 +220 @@
 psF64 PM_MODEL_FLUX (const psVector *params)
 {
-
-    psEllipseShape shape;
-
     psF32 *PAR = params->data.F32;
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
-    shape.sxy = PAR[PM_PAR_SXY];
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
 
     // Area is equivalent to 2 pi sigma^2
-    // axes ratio < 20
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
     psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
 
@@ -242 +237 @@
 psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
 {
-    psEllipseShape shape;
-
     psF32 *PAR = params->data.F32;
 
@@ -253 +246 @@
         return (1.0);
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
+
     psF64 radius = axes.major * sqrt (2.0 * log(PAR[PM_PAR_I0] / flux));
     psAssert (isfinite(radius), "fix this code: radius should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", 
@@ -285 +276 @@
     }
 
-    // the OLD 2D model for SXY actually fitted SXY / (SXX^-2 + SYY^-2); correct here
-    // out[PM_PAR_SXY] = pmPSF_SXYtoModel (out);
-
     // the 2D PSF model fits polarization terms (E0,E1,E2)
     // convert to shape terms (SXX,SYY,SXY)
-    // XXX user-defined value for limit?
-    if (!pmPSF_FitToModel (out, 0.1)) {
-        // psError(PM_ERR_PSF, false, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
+    bool useReff = pmModelUseReff (modelPSF->type);
+    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
         psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
         return false;
@@ -343 +330 @@
     // the 2D PSF model fits polarization terms (E0,E1,E2)
     // convert to shape terms (SXX,SYY,SXY)
-    // XXX user-defined value for limit?
-    if (!pmPSF_FitToModel (PAR, 0.1)) {
+    bool useReff = pmModelUseReff (model->type);
+    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
         psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
         return false;

trunk/psModules/src/objects/models/pmModel_PGAUSS.c

@@ -129 +129 @@
     float q2 = NAN;
     if (nParam == PM_PAR_SXY) {
-        float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
-        float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
+        // NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy
+        float f1 = 2.0 / PS_SQR(params[PM_PAR_SYY]) + 2.0 / PS_SQR(params[PM_PAR_SXX]);
+        float f2 = 2.0 / PS_SQR(params[PM_PAR_SYY]) - 2.0 / PS_SQR(params[PM_PAR_SXX]);
         float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
         q1 = (q1 < 0.0) ? 0.0 : q1;
@@ -201 +202 @@
 
     // set the shape parameters
-    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
+    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, false)) {
       return false;
     }
@@ -222 +223 @@
 {
     float z, norm;
-    psEllipseShape shape;
 
     psF32 *PAR = params->data.F32;
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
+
     float AspectRatio = axes.minor / axes.major;
 
@@ -262 +260 @@
 {
     psF64 z;
-    psEllipseShape shape;
 
     psF32 *PAR = params->data.F32;
@@ -273 +270 @@
         return (1.0);
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    // this estimates the radius assuming f(z) is roughly exp(-z)
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
+
     psF64 sigma = axes.major;
 
@@ -347 +341 @@
     }
 
-    // the OLD 2D model for SXY actually fitted SXY / (SXX^-2 + SYY^-2); correct here
-    // out[PM_PAR_SXY] = pmPSF_SXYtoModel (out);
-
     // the 2D PSF model fits polarization terms (E0,E1,E2)
     // convert to shape terms (SXX,SYY,SXY)
-    if (!pmPSF_FitToModel (out, 0.1)) {
+    bool useReff = pmModelUseReff (modelPSF->type);
+    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
         psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
         return false;
@@ -403 +395 @@
     // the 2D PSF model fits polarization terms (E0,E1,E2)
     // convert to shape terms (SXX,SYY,SXY)
-    // XXX user-defined value for limit?
-    if (!pmPSF_FitToModel (PAR, 0.1)) {
+    bool useReff = pmModelUseReff (model->type);
+    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
         psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
         return false;

trunk/psModules/src/objects/models/pmModel_PS1_V1.c

@@ -1 +1 @@
 /******************************************************************************
- * this file defines the PS1_V1 source shape model.  Note that these model functions are loaded
- * by pmModelClass.c using 'include', and thus need no 'include' statements of their own.  The
- * models use a psVector to represent the set of parameters, with the sequence used to specify
- * the meaning of the parameter.  The meaning of the parameters may thus vary depending on the
- * specifics of the model.  All models which are used as a PSF representations share a few
- * parameters, for which # define names are listed in pmModel.h:
+ * this file defines the PS1_V1 source shape model.  Note that these model functions are
+ * loaded by pmModelClass.c using 'include', and thus need no 'include' statements of
+ * their own.  The models use a psVector to represent the set of parameters, with the
+ * sequence used to specify the meaning of the parameter.  The meaning of the parameters
+ * may thus vary depending on the specifics of the model.  All models which are used as a
+ * PSF representations share a few parameters, for which # define names are listed in
+ * pmModel.h:
 
    power-law with fitted linear term
@@ -148 +149 @@
     float q2 = NAN;
     if (nParam == PM_PAR_SXY) {
-        float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
-        float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
+        float f1 = 2.0 / PS_SQR(params[PM_PAR_SYY]) + 2.0 / PS_SQR(params[PM_PAR_SXX]);
+        float f2 = 2.0 / PS_SQR(params[PM_PAR_SYY]) - 2.0 / PS_SQR(params[PM_PAR_SXX]);
         float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
         q1 = (q1 < 0.0) ? 0.0 : q1;
@@ -220 +221 @@
 
     // set the shape parameters
-    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
+    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, false)) {
       return false;
     }
@@ -244 +245 @@
 {
     float z, norm;
-    psEllipseShape shape;
 
     psF32 *PAR = params->data.F32;
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
     float AspectRatio = axes.minor / axes.major;
 
@@ -284 +281 @@
 {
     psF64 z;
-    psEllipseShape shape;
 
     psF32 *PAR = params->data.F32;
@@ -292 +288 @@
     if (flux >= PAR[PM_PAR_I0]) return 1.0;
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
     psF64 sigma = axes.major;
 
@@ -363 +356 @@
     // the 2D PSF model fits polarization terms (E0,E1,E2)
     // convert to shape terms (SXX,SYY,SXY)
-    if (!pmPSF_FitToModel (out, 0.1)) {
+    bool useReff = pmModelUseReff (modelPSF->type);
+    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
         psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
         return false;
@@ -416 +410 @@
     // convert to shape terms (SXX,SYY,SXY)
     // XXX user-defined value for limit?
-    if (!pmPSF_FitToModel (PAR, 0.1)) {
+    bool useReff = pmModelUseReff (model->type);
+    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
         psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
         return false;

trunk/psModules/src/objects/models/pmModel_QGAUSS.c

@@ -1 +1 @@
 /******************************************************************************
- * this file defines the QGAUSS source shape model (XXX need a better name!).  Note that these
- * model functions are loaded by pmModelClass.c using 'include', and thus need no 'include'
- * statements of their own.  The models use a psVector to represent the set of parameters, with
- * the sequence used to specify the meaning of the parameter.  The meaning of the parameters
- * may thus vary depending on the specifics of the model.  All models which are used a PSF
- * representations share a few parameters, for which # define names are listed in pmModel.h:
+ * this file defines the QGAUSS source shape model.  Note that these model functions are
+ * loaded by pmModelClass.c using 'include', and thus need no 'include' statements of
+ * their own.  The models use a psVector to represent the set of parameters, with the
+ * sequence used to specify the meaning of the parameter.  The meaning of the parameters
+ * may thus vary depending on the specifics of the model.  All models which are used as a
+ * PSF representations share a few parameters, for which # define names are listed in
+ * pmModel.h:
 
    power-law with fitted linear term
@@ -14 +15 @@
    * PM_PAR_XPOS 2  - X center of object
    * PM_PAR_YPOS 3  - Y center of object
-   * PM_PAR_SXX 4   - X^2 term of elliptical contour (sqrt(2) / SigmaX)
-   * PM_PAR_SYY 5   - Y^2 term of elliptical contour (sqrt(2) / SigmaY)
+   * PM_PAR_SXX 4   - X^2 term of elliptical contour (SigmaX / sqrt(2))
+   * PM_PAR_SYY 5   - Y^2 term of elliptical contour (SigmaY / sqrt(2))
    * PM_PAR_SXY 6   - X*Y term of elliptical contour
    * PM_PAR_7   7   - amplitude of the linear component (k)
@@ -138 +139 @@
 # define AR_MAX 20.0
 # define AR_RATIO 0.99
-
 bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
 {
@@ -149 +149 @@
     float q2 = NAN;
     if (nParam == PM_PAR_SXY) {
-        float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
-        float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
+        // NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy
+        float f1 = 2.0 / PS_SQR(params[PM_PAR_SYY]) + 2.0 / PS_SQR(params[PM_PAR_SXX]);
+        float f2 = 2.0 / PS_SQR(params[PM_PAR_SYY]) - 2.0 / PS_SQR(params[PM_PAR_SXX]);
         float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
         q1 = (q1 < 0.0) ? 0.0 : q1;
@@ -203 +204 @@
           return true;
       }
-    default:
+      default:
         psAbort("invalid choice for limits");
     }
@@ -221 +222 @@
 
     // set the shape parameters
-    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
+    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, false)) {
       return false;
     }
@@ -245 +246 @@
 {
     float z, norm;
-    psEllipseShape shape;
 
     psF32 *PAR = params->data.F32;
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
     float AspectRatio = axes.minor / axes.major;
 
@@ -285 +282 @@
 {
     psF64 z;
-    psEllipseShape shape;
 
     psF32 *PAR = params->data.F32;
@@ -293 +289 @@
     if (flux >= PAR[PM_PAR_I0]) return 1.0;
 
-    // if (PAR[PM_PAR_7] == 0.0) return powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA);
-
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
     psF64 sigma = axes.major;
 
@@ -307 +298 @@
         return ( sigma * sqrt (2.0 * z) );
     }
-
     psF64 limit = flux / PAR[PM_PAR_I0];
 
@@ -367 +357 @@
     // the 2D PSF model fits polarization terms (E0,E1,E2)
     // convert to shape terms (SXX,SYY,SXY)
-    if (!pmPSF_FitToModel (out, 0.1)) {
+    bool useReff = pmModelUseReff (modelPSF->type);
+    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
         psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
         return false;
@@ -424 +415 @@
     // the 2D PSF model fits polarization terms (E0,E1,E2)
     // convert to shape terms (SXX,SYY,SXY)
-    // XXX user-defined value for limit?
-    if (!pmPSF_FitToModel (PAR, 0.1)) {
+    bool useReff = pmModelUseReff (model->type);
+    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
         psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
         return false;

trunk/psModules/src/objects/models/pmModel_RGAUSS.c

@@ -138 +138 @@
     float q2 = NAN;
     if (nParam == PM_PAR_SXY) {
-        float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);
-        float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);
+        // NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy
+        float f1 = 2.0 / PS_SQR(params[PM_PAR_SYY]) + 2.0 / PS_SQR(params[PM_PAR_SXX]);
+        float f2 = 2.0 / PS_SQR(params[PM_PAR_SYY]) - 2.0 / PS_SQR(params[PM_PAR_SXX]);
         float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);
         q1 = (q1 < 0.0) ? 0.0 : q1;
@@ -210 +211 @@
 
     // set the shape parameters
-    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments)) {
+    if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, false)) {
       return false;
     }
@@ -234 +235 @@
 {
     float z, norm;
-    psEllipseShape shape;
 
     psF32 *PAR = params->data.F32;
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
     float AspectRatio = axes.minor / axes.major;
 
@@ -274 +271 @@
 {
     psF64 z;
-    psEllipseShape shape;
 
     psF32 *PAR = params->data.F32;
@@ -285 +281 @@
         return (1.0);
 
-    shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false);
     psF64 sigma = axes.major;
 
@@ -357 +350 @@
     // the 2D PSF model fits polarization terms (E0,E1,E2)
     // convert to shape terms (SXX,SYY,SXY)
-    if (!pmPSF_FitToModel (out, 0.1)) {
+    bool useReff = pmModelUseReff (modelPSF->type);
+    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
         psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
         return false;
@@ -409 +403 @@
     // the 2D PSF model fits polarization terms (E0,E1,E2)
     // convert to shape terms (SXX,SYY,SXY)
-    // XXX user-defined value for limit?
-    if (!pmPSF_FitToModel (PAR, 0.1)) {
+    bool useReff = pmModelUseReff (model->type);
+    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
         psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
         return false;

trunk/psModules/src/objects/models/pmModel_SERSIC.c

@@ -125 +125 @@
 
         // first, use Rmajor and index to find the central pixel flux (fraction of total flux)
-        psEllipseShape shape;
-
-        shape.sx  = PAR[PM_PAR_SXX];
-        shape.sy  = PAR[PM_PAR_SYY];
-        shape.sxy = PAR[PM_PAR_SXY];
-
-        // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio
-        psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+        psEllipseAxes axes;
+        pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
 
         // get the central pixel flux from the lookup table
@@ -250 +244 @@
         // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg..  Saturate at that
         // angle and let f2,f1 fight it out
-        q2 = 0.5*sqrtf(q1);
+        // NOTE: the factor of 2 is needed to convert par[SXX,SYY] to shape.sx,sy
+        q2 = 2.0*0.5*sqrtf(q1);
     }
 
@@ -347 +342 @@
     axes.major = Rmajor;
     axes.minor = Rminor;
-    psEllipseShape shape = psEllipseAxesToShape (axes);
-
-    if (!isfinite(shape.sx))  return false;
-    if (!isfinite(shape.sy))  return false;
-    if (!isfinite(shape.sxy)) return false;
+
+    pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], axes, true);
 
     float bn = 1.9992*index - 0.3271;
@@ -357 +349 @@
     float Io = exp(0.5*bn);
 
-    float Sxx = PS_MAX(0.5, shape.sx);
-    float Syy = PS_MAX(0.5, shape.sy);
-
-    PAR[PM_PAR_SXX]  = Sxx;
-    PAR[PM_PAR_SYY]  = Syy;
-    PAR[PM_PAR_SXY]  = shape.sxy;
-
     // set the model normalization (adjust for Sersic best guess)
     if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) {
@@ -381 +366 @@
 psF64 PM_MODEL_FLUX (const psVector *params)
 {
-    psEllipseShape shape;
-
     psF32 *PAR = params->data.F32;
 
-    shape.sx  = PAR[PM_PAR_SXX];
-    shape.sy  = PAR[PM_PAR_SYY];
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    // for a non-circular Sersic, the flux of the Rmajor equivalent is scaled by the AspectRatio
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
     float AspectRatio = axes.minor / axes.major;
 
@@ -410 +389 @@
 psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
 {
-    psEllipseShape shape;
-
     psF32 *PAR = params->data.F32;
 
@@ -421 +398 @@
         return (1.0);
 
-    shape.sx  = PAR[PM_PAR_SXX];
-    shape.sy  = PAR[PM_PAR_SYY];
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
 
     // f = Io exp(-z^n) -> z^n = ln(Io/f)
@@ -457 +431 @@
     // the 2D PSF model fits polarization terms (E0,E1,E2)
     // convert to shape terms (SXX,SYY,SXY)
-    if (!pmPSF_FitToModel (out, 0.1)) {
+    bool useReff = pmModelUseReff (modelPSF->type);
+    if (!pmPSF_FitToModel (out, 0.1, useReff)) {
         psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
         return false;
@@ -510 +485 @@
     // convert to shape terms (SXX,SYY,SXY)
     // XXX user-defined value for limit?
-    if (!pmPSF_FitToModel (PAR, 0.1)) {
+    bool useReff = pmModelUseReff (model->type);
+    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
         psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
         return false;

trunk/psModules/src/objects/models/pmModel_TRAIL.c

@@ -350 +350 @@
     PAR[PM_PAR_SKY]  = 0.0;
 
-    // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms:
-    psEllipseShape psfShape;
-    psfShape.sx  = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2;
-    psfShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY];
-    psfShape.sy  = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2;
-    psEllipseAxes psfAxes = psEllipseShapeToAxes (psfShape, 20.0);
+    psF32 *psfPAR  = source->modelPSF->params->data.F32;
+    bool useReff = pmModelUseReff (source->modelPSF->type);
+
+    psEllipseAxes psfAxes;
+    pmModelParamsToAxes (&psfAxes, psfPAR[PM_PAR_SXX], psfPAR[PM_PAR_SXY], psfPAR[PM_PAR_SYY], useReff);
 
     psEllipseMoments emoments;
@@ -369 +368 @@
     if (!isfinite(axes.theta)) return false;
 
-    float size = (axes.major > sqrt(source->moments->Mrf)) ? axes.major : sqrt(source->moments->Mrf);
-    //    if (size > psfAxes.major) { size -= psfAxes.major; }
-    //else { size = psfAxes.major; }
+    float size = NAN;
+    if (!isfinite(source->moments->Mrf)) {
+      size = axes.major;
+    } else {
+      size = (axes.major > sqrt(source->moments->Mrf)) ? axes.major : sqrt(source->moments->Mrf);
+    }
 
     float theta, peak;

trunk/psModules/src/objects/pmModelUtils.c

@@ -118 +118 @@
 }
 
-bool pmModelSetShape (float *Sxx, float *Sxy, float *Syy, pmMoments *moments) {
+bool pmModelUseReff (pmModelType type) {
+    bool useReff = false;
+    useReff |= (type == pmModelClassGetType ("PS_MODEL_SERSIC"));
+    useReff |= (type == pmModelClassGetType ("PS_MODEL_DEV"));
+    useReff |= (type == pmModelClassGetType ("PS_MODEL_EXP"));
+    return useReff;
+}
+
+// this function and the one below handle the two cases, where the model shape is uses R_eff or Sigma
+bool pmModelAxesToParams (float *Sxx, float *Sxy, float *Syy, psEllipseAxes axes, bool useReff)  {
+
+    psEllipseShape shape = psEllipseAxesToShape (axes);
+
+    if (!isfinite(shape.sx))  return false;
+    if (!isfinite(shape.sy))  return false;
+    if (!isfinite(shape.sxy)) return false;
+
+    // set the shape parameters
+    if (useReff) {
+        *Sxx  = PS_MAX(0.5, shape.sx);
+        *Syy  = PS_MAX(0.5, shape.sy);
+        *Sxy  = shape.sxy * 2.0;
+    } else {
+        *Sxx  = PS_MAX(0.5, M_SQRT2*shape.sx);
+        *Syy  = PS_MAX(0.5, M_SQRT2*shape.sy);
+        *Sxy  = shape.sxy;
+    }
+
+    return true;
+}
+
+bool pmModelParamsToAxes (psEllipseAxes *axes, float Sxx, float Sxy, float Syy, bool useReff)  {
+
+    psEllipseShape shape;
+
+    // set the shape parameters
+    if (useReff) {
+        shape.sx  = Sxx;
+        shape.sy  = Syy;
+        shape.sxy = Sxy / 2.0;
+    } else {
+        shape.sx  = Sxx / M_SQRT2;
+        shape.sy  = Syy / M_SQRT2;
+        shape.sxy = Sxy;
+    }
+
+    if ((shape.sx == 0) || (shape.sy == 0)) {
+        axes->major = 0.0;
+        axes->minor = 0.0;
+        axes->theta = 0.0;
+    } else {
+        // axes ratio < 20
+        // replace with maxAR argument?
+        *axes = psEllipseShapeToAxes (shape, 20.0);
+    }
+
+    return true;
+}
+
+// Reff says if this is a model which uses R_eff (like exp or dev) instead of Sigma
+// set the parameter values SXX, SXY, SYY
+bool pmModelSetShape (float *Sxx, float *Sxy, float *Syy, pmMoments *moments, bool useReff) {
 
     psEllipseMoments emoments;
@@ -137 +198 @@
     axes.minor *= scale;
 
-    psEllipseShape shape = psEllipseAxesToShape (axes);
-
-    if (!isfinite(shape.sx))  return false;
-    if (!isfinite(shape.sy))  return false;
-    if (!isfinite(shape.sxy)) return false;
-
-    // set the shape parameters
-    *Sxx  = PS_MAX(0.5, M_SQRT2*shape.sx);
-    *Syy  = PS_MAX(0.5, M_SQRT2*shape.sy);
-    *Sxy  = shape.sxy;
+    pmModelAxesToParams (Sxx, Sxy, Syy, axes, useReff);
 
     return true;

trunk/psModules/src/objects/pmModelUtils.h

@@ -44 +44 @@
 bool pmModelSetPosition (float *Xo, float *Yo, pmSource *source);
 bool pmModelSetNorm (float *Io, pmSource *source);
-bool pmModelSetShape (float *Sxx, float *Sxy, float *Syy, pmMoments *moments);
+bool pmModelSetShape (float *Sxx, float *Sxy, float *Syy, pmMoments *moments, bool useReff);
+
+bool pmModelUseReff (pmModelType type);
+bool pmModelAxesToParams (float *Sxx, float *Sxy, float *Syy, psEllipseAxes axes, bool useReff);
+bool pmModelParamsToAxes (psEllipseAxes *axes, float Sxx, float Sxy, float Syy, bool useReff);
 
 // XXX void pmModelSetModelVarOption (bool option);

trunk/psModules/src/objects/pmPCM_MinimizeChisq.c

@@ -81 +81 @@
         psAbort ("programming error: no unmasked parameters to be fit\n");
     }
+    psAssert (pcm->nPar == Beta->n, "did we set the masked parameters correctly??");
 
     // allocate internal arrays (current vs Guess)
-    psImage *alpha   = psImageAlloc(Alpha->numCols, Alpha->numRows, PS_TYPE_F32);
-    psVector *beta   = psVectorAlloc(Beta->n, PS_TYPE_F32);
+    psImage *alpha   = psImageAlloc(pcm->nPar, pcm->nPar, PS_TYPE_F32);
+    psVector *beta   = psVectorAlloc(pcm->nPar, PS_TYPE_F32);
     psVector *Params = psVectorAlloc(params->n, PS_TYPE_F32);
 
@@ -90 +91 @@
     psF32 lambda = 0.001;
     psF32 dLinear = 0.0;
+    psF32 nu = 2.0;
 
 # if (USE_FFT && PRE_CONVOLVE)
@@ -120 +122 @@
     bool done = (min->iter >= min->maxIter);
     while (!done) {
-        psTrace("psModules.objects", 5, "Iteration number %d.  (max iterations is %d).\n", min->iter, min->maxIter);
-        psTrace("psModules.objects", 5, "Last delta is %f.  stop if < %f, accept if < %f\n", min->lastDelta, min->minTol, min->maxTol);
+        psTrace(FACILITY, 5, "Iteration number %d.  (max iterations is %d).\n", min->iter, min->maxIter);
+
+        if (min->chisqConvergence) {
+            psTrace(FACILITY, 5, "Last delta is %f.  stop if < %f, accept if < %f\n", min->lastDelta, min->minTol, min->maxTol);
+        } else {
+            psTrace(FACILITY, 5, "Last delta is %f.  stop if < %f, accept if < %f\n", min->rParSigma, min->minTol*pcm->nPar, min->maxTol*pcm->nPar);
+        }
 
         // set a new guess for Alpha, Beta, Params
@@ -140 +147 @@
             p_psVectorPrint(psTraceGetDestination(), Params, "params guess (1)");
         }
+
+        // calculate the parameter change (rParDelta) and error radius (rParSigma)
+        //    rParDelta : radius of parameter change;
+        //    rParSigma : radius of parameter error 
+        
+        // note that (before SetABX) Alpha[i][i] is the covariance matrix and
+        // Beta is the actual parameter change for this pass
+
+        // note that Alpha & Beta only represent unmasked parameters, while params and Params have all 
+
+        // dParSigma = Alpha[i][i] : error (squared) on parameter i
+        // dParDelta = Params->data.F32[i] - params->data.F32[i]     : change on parameter i
+        float rParSigma = 0.0;
+        for (int j = 0, J = 0; j < Params->n; j++) {
+            if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[j])) {
+                continue;
+            }
+            rParSigma += PS_SQR(Params->data.F32[j] - params->data.F32[j]) / Alpha->data.F32[J][J];
+            J++;
+        }
+        rParSigma = sqrt(rParSigma);
+        psTrace(FACILITY, 5, "rParSigma: %f, Niter: %d\n", rParSigma, min->iter);
+        // fprintf (stderr, "rParSigma: %f, Niter: %d\n", rParSigma, min->iter);
 
         // calculate Chisq for new guess, update Alpha & Beta
@@ -164 +194 @@
         }
 
-        /* if (Chisq < min->value) {  */
+        // change in chisq/nDOF since last minimum
+        min->lastDelta = (min->value - Chisq) / pcm->nDOF;
+
+        // rho is positive if the new chisq is smaller; allow for some insignificant change (slight negative rho)
+
+        // XXX the old version of lambda changes:
+        // XXX : Madsen gives suggestion for better use of rho
+        // rho is positive if the new chisq is smaller
         if (rho >= -1e-6) {
-            min->lastDelta = (min->value - Chisq) / (source->pixels->numCols*source->pixels->numRows - params->n);
             min->value = Chisq;
             alpha  = psImageCopy(alpha, Alpha, PS_TYPE_F32);
             beta   = psVectorCopy(beta, Beta, PS_TYPE_F32);
             params = psVectorCopy(params, Params, PS_TYPE_F32);
-            lambda *= 0.25;
 
             // save the new convolved model image
             psFree (source->modelFlux);
             source->modelFlux = pmPCMdataSaveImage(pcm);
-        } else {
-            lambda *= 10.0;
-        }
+        } 
+        switch (min->gainFactorMode) {
+          case 0:
+            if (rho >= -1e-6) {
+                lambda *= 0.25;
+            } else {
+                lambda *= 10.0;
+            }
+            break;
+
+          case 1:
+            // adjust the gain ratio (lambda) based on rho
+            if (rho < 0.25) {
+                lambda *= 2.0;
+            } 
+            if (rho > 0.75) {
+                lambda *= 0.333;
+            }
+            break;
+
+          case 2:
+            if (rho > 0.0) {
+                lambda *= PS_MAX(0.33, (1.0 - pow(2.0*rho - 1.0, 3.0)));
+                nu = 2.0;
+            } else {
+                lambda *= nu;
+                nu *= 2.0;
+            }
+            break;
+        }
         min->iter++;
 
+        // ending conditions:
+        // 1) hard limit : too many iterations
         done = (min->iter >= min->maxIter);
         
-        // check for convergence:
+        // 2) require deltaChi > 1e-6 (ie, chisq is decreasing, but accept an insignificant change)
+        if (min->lastDelta < -1e-6) {
+            continue;
+        }
+
+        // save this value in case we stop iterating
+        min->rParSigma = rParSigma;
+
+        // 2) require chisqDOF < maxChisqDOF (if maxChisqDOF is not NAN)
+        // keep iterating regardless of rParSigma in this case
         float chisqDOF = Chisq / pcm->nDOF;
-        if (!isfinite(min->maxChisqDOF) || ((chisqDOF < min->maxChisqDOF) && isfinite(min->lastDelta))) {
+        if (isfinite(min->maxChisqDOF) && (chisqDOF > min->maxChisqDOF)) {
+            continue;
+        }
+
+        // delta-chisq or rParSigma ?
+        if (min->chisqConvergence) {
             done |= (min->lastDelta < min->minTol);
+        } else {
+            done |= (rParSigma < min->minTol*pcm->nPar);
         }
     }
@@ -220 +300 @@
 
     // if the last improvement was at least as good as maxTol, accept the fit:
-    if (min->lastDelta <= min->maxTol) {
-        psTrace(FACILITY, 6, "---- end (true) ----\n");
-        return(true);
+    if (min->chisqConvergence) {
+        if (min->lastDelta <= min->maxTol) {
+            psTrace(FACILITY, 6, "---- end (true) ----\n");
+            return(true);
+        }
+    } else {
+        if (min->rParSigma <= min->maxTol*pcm->nPar) {
+            psTrace(FACILITY, 6, "---- end (true) ----\n");
+            return(true);
+        }
     }
     psTrace(FACILITY, 6, "---- end (false) ----\n");

trunk/psModules/src/objects/pmPCMdata.c

@@ -291 +291 @@
     psAssert (modelPSF, "psf model must be defined");
     
-    psEllipseShape shape;
     psEllipseAxes axes;
-
-    shape.sx  = modelPSF->params->data.F32[PM_PAR_SXX];
-    shape.sy  = modelPSF->params->data.F32[PM_PAR_SYY];
-    shape.sxy = modelPSF->params->data.F32[PM_PAR_SXY];
-    axes = psEllipseShapeToAxes (shape, 20.0);
+    bool useReff = pmModelUseReff (modelPSF->type);
+    psF32 *PAR = modelPSF->params->data.F32;
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff);
     
-    float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*modelPSF->params->data.F32[PM_PAR_I0]);
+    float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]);
     float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
 
@@ -451 +448 @@
         psAssert (modelPSF, "psf model must be defined");
     
-        psEllipseShape shape;
         psEllipseAxes axes;
-
-        shape.sx  = modelPSF->params->data.F32[PM_PAR_SXX];
-        shape.sy  = modelPSF->params->data.F32[PM_PAR_SYY];
-        shape.sxy = modelPSF->params->data.F32[PM_PAR_SXY];
-        axes = psEllipseShapeToAxes (shape, 20.0);
+        bool useReff = pmModelUseReff (modelPSF->type);
+        psF32 *PAR = modelPSF->params->data.F32;
+        pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff);
     
-        float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*modelPSF->params->data.F32[PM_PAR_I0]);
+        float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5*PAR[PM_PAR_I0]);
         float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
 

trunk/psModules/src/objects/pmPSF.c

@@ -249 +249 @@
 // Mxy = SXY * (SXX^-4 + SYY^-4 - 2 SXY ^2)
 
+// XXX deprecated
 // input: model->param, output: psf->param[PM_PAR_SXY]
-double pmPSF_SXYfromModel (psF32 *modelPar)
-{
-    PS_ASSERT_PTR_NON_NULL(modelPar, NAN);
-
-    double SXX = modelPar[PM_PAR_SXX];
-    double SYY = modelPar[PM_PAR_SYY];
-    double SXY = modelPar[PM_PAR_SXY];
-
-    double par = SXY / PS_SQR(1.0 / PS_SQR(SXX) + 1.0 / PS_SQR(SYY));
-    return (par);
-}
-
+// XXX double pmPSF_SXYfromModel (psF32 *modelPar)
+// XXX {
+// XXX     PS_ASSERT_PTR_NON_NULL(modelPar, NAN);
+// XXX 
+// XXX     double SXX = modelPar[PM_PAR_SXX];
+// XXX     double SYY = modelPar[PM_PAR_SYY];
+// XXX     double SXY = modelPar[PM_PAR_SXY];
+// XXX 
+// XXX     double par = SXY / PS_SQR(1.0 / PS_SQR(SXX) + 1.0 / PS_SQR(SYY));
+// XXX     return (par);
+// XXX }
+
+// XXX deprecated
 // input: fitted psf->param, output: model->param[PM_PAR_SXY]
-double pmPSF_SXYtoModel (psF32 *fittedPar)
-{
-    PS_ASSERT_PTR_NON_NULL(fittedPar, NAN);
-
-    double SXX = fittedPar[PM_PAR_SXX];
-    double SYY = fittedPar[PM_PAR_SYY];
-    double fit = fittedPar[PM_PAR_SXY];
-
-    double SXY = fit * PS_SQR(1.0 / PS_SQR(SXX) + 1.0 / PS_SQR(SYY));
-
-    assert (!isnan(SXY));
-
-    return SXY;
-}
-
-// New Concept: the PSF modelling function fits the polarization terms e0, e1, e2:
-
-// convert the parameters used in the fitted source model
-// to the parameters used in the 2D PSF model
-// XXX this function may be invalid for SERSIC, DEV, EXP models (SQRT2 not used?)
-bool pmPSF_FitToModel (psF32 *fittedPar, float minMinorAxis)
+// XXX double pmPSF_SXYtoModel (psF32 *fittedPar)
+// XXX {
+// XXX     PS_ASSERT_PTR_NON_NULL(fittedPar, NAN);
+// XXX 
+// XXX     double SXX = fittedPar[PM_PAR_SXX];
+// XXX     double SYY = fittedPar[PM_PAR_SYY];
+// XXX     double fit = fittedPar[PM_PAR_SXY];
+// XXX 
+// XXX     double SXY = fit * PS_SQR(1.0 / PS_SQR(SXX) + 1.0 / PS_SQR(SYY));
+// XXX 
+// XXX     assert (!isnan(SXY));
+// XXX 
+// XXX     return SXY;
+// XXX }
+
+// The PSF modelling function fits the polarization terms e0, e1, e2:
+
+// the FIT is the 2D representation of the shape using polarization parameters for the elliptical contour
+// the MODEL is the realized psf model for a given location
+
+// convert the parameters (in situ) used in the fitted source model to the parameters used in
+// the 2D PSF model
+bool pmPSF_FitToModel (psF32 *fittedPar, float minMinorAxis, bool useReff)
 {
     PS_ASSERT_PTR_NON_NULL(fittedPar, false);
@@ -298 +302 @@
         return false;
     }
-    psEllipseShape shape = psEllipseAxesToShape (axes);
-
-    fittedPar[PM_PAR_SXX] = shape.sx * M_SQRT2;
-    fittedPar[PM_PAR_SYY] = shape.sy * M_SQRT2;
-    fittedPar[PM_PAR_SXY] = shape.sxy;
-
+
+    pmModelAxesToParams (&fittedPar[PM_PAR_SXX], &fittedPar[PM_PAR_SXY], &fittedPar[PM_PAR_SYY], axes, useReff);
     return true;
 }
 
-// convert the PSF parameters used in the 2D PSF model fit into the
-// parameters used in the source model
-// XXX this function may be invalid for SERSIC, DEV, EXP models (SQRT2 not used?)
-psEllipsePol pmPSF_ModelToFit (psF32 *modelPar)
+// convert the parameters (in situ) used in the 2D PSF model fit into the parameters used in
+// the source model
+psEllipsePol pmPSF_ModelToFit (psF32 *modelPar, bool useReff)
 {
     // must assert non-NULL input parameter
@@ -319 +318 @@
     PS_ASSERT_PTR_NON_NULL(modelPar, pol);
 
-    psEllipseShape shape;
-
-    shape.sx  = modelPar[PM_PAR_SXX] / M_SQRT2;
-    shape.sy  = modelPar[PM_PAR_SYY] / M_SQRT2;
-    shape.sxy = modelPar[PM_PAR_SXY];
-
-    pol = psEllipseShapeToPol (shape);
+    psEllipseAxes axes;
+    pmModelParamsToAxes (&axes, modelPar[PM_PAR_SXX], modelPar[PM_PAR_SXY], modelPar[PM_PAR_SYY], useReff);
+
+    pol = psEllipseAxesToPol (axes);
 
     return pol;
@@ -332 +328 @@
 // convert the parameters used in the fitted source model to the psEllipseAxes representation
 // (major,minor,theta)
-psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, double maxAR, pmModelType type)
-{
-    psEllipseShape shape;
+psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, pmModelType type)
+{
     psEllipseAxes axes;
     axes.major = NAN;
     axes.minor = NAN;
     axes.theta = NAN;
-    //   XXX: must assert non-NULL input parameter
+
     PS_ASSERT_PTR_NON_NULL(modelPar, axes);
 
-    bool useReff = false;
-    useReff |= (type == pmModelClassGetType ("PS_MODEL_SERSIC"));
-    useReff |= (type == pmModelClassGetType ("PS_MODEL_DEV"));
-    useReff |= (type == pmModelClassGetType ("PS_MODEL_EXP"));
-
-    if (useReff) {
-        shape.sx  = modelPar[PM_PAR_SXX];
-        shape.sy  = modelPar[PM_PAR_SYY];
-        shape.sxy = modelPar[PM_PAR_SXY] / 2.0;
-        // XXX I *think* dividing by 2.0 is the right direction, but this 
-        // needs to be checked with a real test
-    } else {
-        shape.sx  = modelPar[PM_PAR_SXX] / M_SQRT2;
-        shape.sy  = modelPar[PM_PAR_SYY] / M_SQRT2;
-        shape.sxy = modelPar[PM_PAR_SXY];
-    }
-
-    if ((shape.sx == 0) || (shape.sy == 0)) {
-        axes.major = 0.0;
-        axes.minor = 0.0;
-        axes.theta = 0.0;
-    } else {
-        // XXX this is not really consistent with the model fit range above
-        axes = psEllipseShapeToAxes (shape, maxAR);
-    }
-
+    bool useReff = pmModelUseReff (type);
+    pmModelParamsToAxes (&axes, modelPar[PM_PAR_SXX], modelPar[PM_PAR_SXY], modelPar[PM_PAR_SYY], useReff);
     return axes;
 }
@@ -377 +348 @@
     PS_ASSERT_PTR_NON_NULL(modelPar, false);
 
+    modelPar[PM_PAR_SXX] = 0.0;
+    modelPar[PM_PAR_SYY] = 0.0;
+    modelPar[PM_PAR_SXY] = 0.0;
+    
     if ((axes.major <= 0) || (axes.minor <= 0)) {
-        modelPar[PM_PAR_SXX] = 0.0;
-        modelPar[PM_PAR_SYY] = 0.0;
-        modelPar[PM_PAR_SXY] = 0.0;
         return true;
     }
-
-    psEllipseShape shape = psEllipseAxesToShape (axes);
-
-    bool useReff = false;
-    useReff |= ( type == pmModelClassGetType ("PS_MODEL_SERSIC"));
-    useReff |= ( type == pmModelClassGetType ("PS_MODEL_DEV"));
-    useReff |= ( type == pmModelClassGetType ("PS_MODEL_EXP"));
-
-    if (useReff) {
-        modelPar[PM_PAR_SXX] = shape.sx;
-        modelPar[PM_PAR_SYY] = shape.sy;
-        modelPar[PM_PAR_SXY] = shape.sxy * 2.0; // XXX NEED factor of 2 here for correct angle conversion
-    } else {
-        modelPar[PM_PAR_SXX] = shape.sx * M_SQRT2;
-        modelPar[PM_PAR_SYY] = shape.sy * M_SQRT2;
-        modelPar[PM_PAR_SXY] = shape.sxy;
-    }
+    
+    bool useReff = pmModelUseReff (type);
+    pmModelAxesToParams (&modelPar[PM_PAR_SXX], &modelPar[PM_PAR_SXY], &modelPar[PM_PAR_SYY], axes, useReff);
     return true;
 }
@@ -423 +381 @@
     par->data.F32[PM_PAR_SXY] = sxy;
 
-    psEllipsePol pol = pmPSF_ModelToFit(par->data.F32);
+    bool useReff = pmModelUseReff (options->type);
+    psEllipsePol pol = pmPSF_ModelToFit(par->data.F32, useReff);
 
     pmTrend2D *trend = NULL;

trunk/psModules/src/objects/pmPSF.h

@@ -107 +107 @@
 
 bool pmPSF_AxesToModel (psF32 *modelPar, psEllipseAxes axes, pmModelType type);
-bool pmPSF_FitToModel (psF32 *fittedPar, float minMinorAxis);
+bool pmPSF_FitToModel (psF32 *fittedPar, float minMinorAxis, bool useReff);
 
-psEllipsePol pmPSF_ModelToFit (psF32 *modelPar);
-psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, double maxAR, pmModelType type);
+psEllipsePol pmPSF_ModelToFit (psF32 *modelPar, bool useReff);
+psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, pmModelType type);
 
 /// Calculate FWHM value from a PSF

trunk/psModules/src/objects/pmPSFtry.h

@@ -99 +99 @@
 /** fit EXT models to all possible psf sources */
 bool pmPSFtryFitEXT (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal);
+bool pmPSFtryFitEXT_Threaded (psThreadJob *job);
 
 bool pmPSFtryMakePSF (bool *pGoodFit, pmPSFtry *psfTry);
 
 bool pmPSFtryFitPSF (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal);
+bool pmPSFtryFitPSF_Threaded (psThreadJob *job);
+
+bool pmPSFThreads (void);
 
 /** pmPSFtryMetric()

trunk/psModules/src/objects/pmPSFtryFitEXT.c

@@ -46 +46 @@
 #include "pmSourceVisual.h"
 
+bool pmPSFThreads (void) {
+
+    psThreadTask *task = NULL;
+
+    task = psThreadTaskAlloc("PSF_TRY_FIT_EXT", 6);
+    task->function = &pmPSFtryFitEXT_Threaded;
+    psThreadTaskAdd(task);
+    psFree(task);
+
+    task = psThreadTaskAlloc("PSF_TRY_FIT_PSF", 6);
+    task->function = &pmPSFtryFitPSF_Threaded;
+    psThreadTaskAdd(task);
+    psFree(task);
+
+    return true;
+}
+
+static int Next = 0;
+
 // Fit an EXT model to all candidates PSF sources.
 // Note: this is independent of the modeled 2D variations in the PSF.
 bool pmPSFtryFitEXT (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal) {
-
-    bool status;
 
     psTimerStart ("psf.fit");
@@ -60 +77 @@
     maskVal |= markVal;
 
-    int Next = 0;
+    Next = 0;
     for (int i = 0; i < psfTry->sources->n; i++) {
 
         pmSource *source = psfTry->sources->data[i];
-        if (!source->moments) {
-            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
-            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : no moments\n", i, source->peak->x, source->peak->y);
-            continue;
-        }
-        if (!source->moments->nPixels) {
-            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : no pixels\n", i, source->peak->x, source->peak->y);
-            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
-            continue;
-        }
-        // If mask object does not exist, mark the source as bad.
-        // We cannot proceed with it because psImageMaskPixels leaves an uncleared error code last which causes
-        // psphot to exit with a fault. 
-        if (source->maskObj == NULL) {
-            psTrace ("psModules.objects", 4, "source %d (%d,%d) : null maskObj\n", i, source->peak->x, source->peak->y);
-            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
-            continue;
-        }
 
-        source->modelEXT = pmSourceModelGuess (source, options->type, maskVal, markVal);
-        if (source->modelEXT == NULL) {
-            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : failed to generate model guess\n", i, source->peak->x, source->peak->y);
-            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
-            continue;
-        }
+        if (!source->moments) {
+            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
+            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : no moments\n", i, source->peak->x, source->peak->y);
+            continue;
+        }
+        if (!source->moments->nPixels) {
+            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : no pixels\n", i, source->peak->x, source->peak->y);
+            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
+            continue;
+        }
+        // If mask object does not exist, mark the source as bad.
+        // We cannot proceed with it because psImageMaskPixels leaves an uncleared error code last which causes
+        // psphot to exit with a fault. 
+        if (source->maskObj == NULL) {
+            psTrace ("psModules.objects", 4, "source %d (%d,%d) : null maskObj\n", i, source->peak->x, source->peak->y);
+            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
+            continue;
+        }
 
-        // set object mask to define valid pixels
-        // XXX 0.5 PIX: is the circle symmetric about the peak coordinate (given 0.5,0.5 center)?
-        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "OR", markVal);
+        source->modelEXT = pmSourceModelGuess (source, options->type, maskVal, markVal);
+        if (source->modelEXT == NULL) {
+            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : failed to generate model guess\n", i, source->peak->x, source->peak->y);
+            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
+            continue;
+        }
 
-        // fit model as EXT, not PSF
-        status = pmSourceFitModel (source, source->modelEXT, options->fitOptions, maskVal);
+        // do some actual work on this source
+        psThreadJob *job = psThreadJobAlloc ("PSF_TRY_FIT_EXT");
+        psArrayAdd(job->args, 1, source);
+        psArrayAdd(job->args, 1, psfTry);
+        psArrayAdd(job->args, 1, options);
+        
+        PS_ARRAY_ADD_SCALAR(job->args, i,        PS_TYPE_S32);
 
-        // clear object mask to define valid pixels
-        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
+        PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
+        PS_ARRAY_ADD_SCALAR(job->args, markVal,  PS_TYPE_IMAGE_MASK);
 
-        // exclude the poor fits
-        if (!status) {
-            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : status is poor\n", i, source->peak->x, source->peak->y);
-            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
-            continue;
-        }
-        Next ++;
+# if (1)
+        if (!psThreadJobAddPending(job)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to create psf model.");
+            return false;
+        }
+# else
+        if (!pmPSFtryFitEXT_Threaded(job)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to create psf model.");
+            return false;
+        }
+        psFree(job);
+# endif
     }
+
+    // wait for the threads to finish and manage results
+    if (!psThreadPoolWait (false, true)) {
+        psError(PS_ERR_UNKNOWN, false, "failure to model psf");
+        return false;
+    }
+
+    // we have only supplied one type of job, so we can assume the types here
+    psThreadJob *job = NULL;
+    while ((job = psThreadJobGetDone()) != NULL) {
+        // we have no returned data from this operation
+        if (job->args->n < 1) fprintf (stderr, "error with job\n");
+        psFree(job);
+    }
+
     psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit ext:   %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Next, psfTry->sources->n);
     psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (EXT)\n", Next, psfTry->sources->n);
@@ -118 +157 @@
     return true;
 }
+
+bool pmPSFtryFitEXT_Threaded (psThreadJob *job) {
+
+    pmSource *source =      job->args->data[0];
+    pmPSFtry *psfTry =      job->args->data[1];
+    pmPSFOptions *options = job->args->data[2];
+
+    int i = PS_SCALAR_VALUE(job->args->data[3], S32);
+
+    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
+    psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
+
+    // set object mask to define valid pixels
+    // XXX 0.5 PIX: is the circle symmetric about the peak coordinate (given 0.5,0.5 center)?
+    psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "OR", markVal);
+
+    // fit model as EXT, not PSF
+    bool status = pmSourceFitModel (source, source->modelEXT, options->fitOptions, maskVal);
+
+    // clear object mask to define valid pixels
+    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
+
+    // exclude the poor fits
+    if (!status) {
+        psTrace ("psModules.objects", 4, "masking %d (%d,%d) : status is poor\n", i, source->peak->x, source->peak->y);
+        psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
+        return true;
+    }
+    Next ++;
+    
+    return true;
+}

trunk/psModules/src/objects/pmPSFtryFitPSF.c

@@ -43 +43 @@
 #include "pmSourceVisual.h"
 
+static int Npsf = 0;
+
 // stage 3: Refit with fixed shape parameters.  This function uses the LMM fitting, but could
 // be re-written to use the simultaneous linear fitting (see psphotFitSourcesLinear.c)
 bool pmPSFtryFitPSF (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal) {
 
-    bool status;
-
     psTimerStart ("psf.fit");
 
@@ -57 +57 @@
     maskVal |= markVal;
 
-    // DEBUG code: save the PSF model fit data in detail
-# ifdef DEBUG
-    char filename[64];
-    snprintf (filename, 64, "psffit.%dx%d.dat", psfTry->psf->trendNx, psfTry->psf->trendNy);
-    FILE *f = fopen (filename, "w");
-    psAssert (f, "failed open");
-# endif
-
-    int Npsf = 0;
+    Npsf = 0;
     for (int i = 0; i < psfTry->sources->n; i++) {
 
@@ -77 +69 @@
         }
 
-        // set shape for this model based on PSF
+        // set shape for this model based on PSF
         psFree (source->modelPSF);
-        source->modelPSF = pmModelFromPSF (source->modelEXT, psfTry->psf);
-        if (source->modelPSF == NULL) {
-            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL;
-            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : bad PSF fit\n", i, source->peak->x, source->peak->y);
-            continue;
-        }
-        // PSF fit and aperture mags use different radii
-        source->modelPSF->fitRadius = options->fitRadius;
-        source->apRadius = options->apRadius;
-
-        // set object mask to define valid pixels for PSF model fit
-        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "OR", markVal);
-
-        // fit the PSF model to the source
-        status = pmSourceFitModel (source, source->modelPSF, options->fitOptions, maskVal);
-
-        // skip poor fits
-        if (!status) {
-            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
-            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_PSF_FAIL;
-            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y);
-            continue;
-        }
-
-        // set object mask to define valid pixels for APERTURE magnitude
-        if (options->fitRadius != options->apRadius) {
-            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
-            psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->apRadius, "OR", markVal);
+        source->modelPSF = pmModelFromPSF (source->modelEXT, psfTry->psf);
+        if (source->modelPSF == NULL) {
+            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL;
+            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : bad PSF fit\n", i, source->peak->x, source->peak->y);
+            return false;
         }
 
-        // This function calculates the psf and aperture magnitudes
-        status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal, markVal, options->apRadius); // raw PSF mag, AP mag
-        if (!status || isnan(source->apMag)) {
-            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
-            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_PHOT;
-            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y);
-            continue;
-        }
-
-        // clear object mask to define valid pixels
-        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
-
-        psfTry->fitMag->data.F32[i] = source->psfMag;
-        psfTry->metric->data.F32[i] = source->apMag - source->psfMag;
-        psfTry->metricErr->data.F32[i] = source->psfMagErr;
-
-        // XXX this did not work: modifies shape of psf too much
-        // psfTry->metric->data.F32[i] = -2.5*log10(source->moments->Sum) - source->psfMag;
-
+        // do some actual work on this source
+        psThreadJob *job = psThreadJobAlloc ("PSF_TRY_FIT_PSF");
+        psArrayAdd(job->args, 1, source);
+        psArrayAdd(job->args, 1, psfTry);
+        psArrayAdd(job->args, 1, options);
+        
+        PS_ARRAY_ADD_SCALAR(job->args, i,        PS_TYPE_S32);
+
+        PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
+        PS_ARRAY_ADD_SCALAR(job->args, markVal,  PS_TYPE_IMAGE_MASK);
+
+# if (1)
+        if (!psThreadJobAddPending(job)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to create psf model.");
+            return false;
+        }
+# else
+        if (!pmPSFtryFitPSF_Threaded(job)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to create psf model.");
+            return false;
+        }
+        psFree(job);
+# endif
+    }
+
+    // wait for the threads to finish and manage results
+    if (!psThreadPoolWait (false, true)) {
+        psError(PS_ERR_UNKNOWN, false, "failure to model psf");
+        return false;
+    }
+
+    // we have only supplied one type of job, so we can assume the types here
+    psThreadJob *job = NULL;
+    while ((job = psThreadJobGetDone()) != NULL) {
+        // we have no returned data from this operation
+        if (job->args->n < 1) fprintf (stderr, "error with job\n");
+        psFree(job);
+    }
+    psfTry->psf->nPSFstars = Npsf;
+
+    // DEBUG code: save the PSF model fit data in detail
 # ifdef DEBUG
+
+    char filename[64];
+    snprintf (filename, 64, "psffit.%dx%d.dat", psfTry->psf->trendNx, psfTry->psf->trendNy);
+    FILE *f = fopen (filename, "w");
+    psAssert (f, "failed open");
+
+    for (int i = 0; i < psfTry->sources->n; i++) {
+
+        // skip masked sources
+        if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) continue;
+
+        pmSource *source = psfTry->sources->data[i];
+
         fprintf (f, "%6.1f %6.1f : %6.1f %6.1f : %8.3f %8.3f %8.3f : %f : %f %f %f : %f\n",
                  source->peak->xf, source->peak->yf, 
@@ -136 +140 @@
                  source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY], 
                  source->modelPSF->params->data.F32[PM_PAR_SYY], source->modelPSF->params->data.F32[PM_PAR_7]);
-# endif
-
-        psTrace ("psModules.object", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n);
-        Npsf ++;
-    }
-    psfTry->psf->nPSFstars = Npsf;
-
-# ifdef DEBUG
+    }
     fclose (f);
 # endif
@@ -159 +156 @@
     return true;
 }
+
+bool pmPSFtryFitPSF_Threaded (psThreadJob *job) {
+
+    pmSource *source =      job->args->data[0];
+    pmPSFtry *psfTry =      job->args->data[1];
+    pmPSFOptions *options = job->args->data[2];
+
+    int i = PS_SCALAR_VALUE(job->args->data[3], S32);
+
+    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
+    psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
+
+    // PSF fit and aperture mags use different radii
+    source->modelPSF->fitRadius = options->fitRadius;
+    source->apRadius            = options->apRadius;
+
+    // set object mask to define valid pixels for PSF model fit
+    psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "OR", markVal);
+
+    // fit the PSF model to the source
+    bool status = pmSourceFitModel (source, source->modelPSF, options->fitOptions, maskVal);
+
+    // skip poor fits
+    if (!status) {
+        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
+        psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_PSF_FAIL;
+        psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y);
+        return true;
+    }
+
+    // set object mask to define valid pixels for APERTURE magnitude
+    if (options->fitRadius != options->apRadius) {
+        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
+        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->apRadius, "OR", markVal);
+    }
+
+    // This function calculates the psf and aperture magnitudes
+    status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal, markVal, options->apRadius); // raw PSF mag, AP mag
+    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
+
+    if (!status || isnan(source->apMag)) {
+        psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_PHOT;
+        psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y);
+        return true;
+    }
+
+    psfTry->fitMag->data.F32[i] = source->psfMag;
+    psfTry->metric->data.F32[i] = source->apMag - source->psfMag;
+    psfTry->metricErr->data.F32[i] = source->psfMagErr;
+
+    psTrace ("psModules.object", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n);
+    Npsf ++;
+    return true;
+}

trunk/psModules/src/objects/pmPSFtryMakePSF.c

@@ -165 +165 @@
 
             for (int i = 0; i < psf->params->n; i++) {
-                if (psf->params->data[i] == NULL) {
-                    psFree(modelPSF);
-                    continue;
-                }
+                if (psf->params->data[i] == NULL) continue;
                 fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[i], modelPSF->params->data.F32[i]);
             }
@@ -214 +211 @@
         assert (source->modelEXT); // all unmasked sources should have modelEXT
 
-        psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32);
+        bool useReff = pmModelUseReff (source->modelEXT->type);
+        psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32, useReff);
 
         e0->data.F32[i] = pol.e0;

trunk/psModules/src/objects/pmSource.c

@@ -1145 +1145 @@
     bool status;
     psEllipseShape oldshape;
-    psEllipseShape newshape;
     psEllipseAxes axes;
 
@@ -1166 +1165 @@
     if (!isfinite(oldI0)) return false;
 
+    bool useReff = pmModelUseReff (model->type);
+    pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], useReff);
+
     // increase size and height of source
-    axes = psEllipseShapeToAxes (oldshape, 20.0);
     axes.major *= SIZE;
     axes.minor *= SIZE;
-    newshape = psEllipseAxesToShape (axes);
+
+    pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], axes, useReff);
     PAR[PM_PAR_I0]  = FACTOR*oldI0;
-    PAR[PM_PAR_SXX] = newshape.sx;
-    PAR[PM_PAR_SYY] = newshape.sy;
-    PAR[PM_PAR_SXY] = newshape.sxy;
 
     psImage *target = source->variance;

trunk/psModules/src/objects/pmSourceFitModel.c

@@ -62 +62 @@
     opt->poissonErrors = true;
     opt->saveCovariance = false;
+
+    // we default to the old algorithm
+    opt->gainFactorMode = 0;
+    opt->chisqConvergence = true;
 
     return opt;
@@ -241 +245 @@
 
     psMinimization *myMin = psMinimizationAlloc (options->nIter, options->minTol, options->maxTol);
+    myMin->gainFactorMode = options->gainFactorMode;
+    myMin->chisqConvergence = options->chisqConvergence;
 
     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
@@ -252 +258 @@
     }
     if (options->saveCovariance) {
+        psFree (model->covar);
         model->covar = psMemIncrRefCounter(covar);
     }
@@ -273 +280 @@
     model->flags |= PM_MODEL_STATUS_FITTED;
     if (!fitStatus) model->flags |= PM_MODEL_STATUS_NONCONVERGE;
-    if (myMin->lastDelta > myMin->minTol) model->flags |= PM_MODEL_STATUS_WEAK_FIT;
+
+    if (myMin->chisqConvergence) {
+      if (myMin->lastDelta > myMin->minTol) model->flags |= PM_MODEL_STATUS_WEAK_FIT;
+    } else {
+      if (myMin->rParSigma > myMin->minTol*nParams) model->flags |= PM_MODEL_STATUS_WEAK_FIT;
+    }
 
     // get the Gauss-Newton distance for fixed model parameters

trunk/psModules/src/objects/pmSourceFitModel.h

@@ -35 +35 @@
     bool poissonErrors;                 ///< use poisson errors for fits?
     bool saveCovariance;
+    int gainFactorMode;
+    bool chisqConvergence; 
 } pmSourceFitOptions;
 

trunk/psModules/src/objects/pmSourceFitPCM.c

@@ -66 +66 @@
     // set up the minimization process
     psMinimization *myMin = psMinimizationAlloc (fitOptions->nIter, fitOptions->minTol, fitOptions->maxTol);
+    myMin->chisqConvergence = fitOptions->chisqConvergence;
+    myMin->gainFactorMode = fitOptions->gainFactorMode;
 
     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
@@ -118 +120 @@
     if (!fitStatus) pcm->modelConv->flags |= PM_MODEL_STATUS_NONCONVERGE;
 
+    if (myMin->chisqConvergence) {
+      if (myMin->lastDelta > myMin->minTol) pcm->modelConv->flags |= PM_MODEL_STATUS_WEAK_FIT;
+    } else {
+      if (myMin->rParSigma > myMin->minTol*pcm->nPar) pcm->modelConv->flags |= PM_MODEL_STATUS_WEAK_FIT;
+    }
+
     // once we have fitted a model, we need to record that this model is a PCM model:
     pcm->modelConv->isPCM = true;
@@ -145 +153 @@
 bool pmSourceModelGuessPCM (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal) {
 
-  if (!pcm->modelConv->modelGuess(pcm->modelConv, source, maskVal, markVal)) {
+    if (!pcm->modelConv->modelGuess(pcm->modelConv, source, maskVal, markVal)) {
         return false;
     }
     return true;
 
+    // the following was an attempt to make analytical modifications to the shape terms based on the psf
+    // this has been replaced with a more empirical approach
+# if (0)
     // generate copy of the model
     // XXX we could modify the parameter values or even the model 
@@ -196 +207 @@
 
     return true;
+# endif
 }

trunk/psModules/src/objects/pmSourceFitSet.c

@@ -568 +568 @@
 
     psMinimization *myMin = psMinimizationAlloc (options->nIter, options->minTol, options->maxTol);
+    myMin->gainFactorMode = options->gainFactorMode;
+    myMin->chisqConvergence = options->chisqConvergence;
 
     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);

trunk/psModules/src/objects/pmSourceIO_CMF.c.in

@@ -925 +925 @@
                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR",     0, "EXT angle err (SXY, isnan)", dPAR[PM_PAR_SXY]);
                 } else {
-                    psEllipseAxes axes = pmPSF_ModelToAxes (PAR, 20.0, model->type);
+                    psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type);
                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width (major axis), length for trail", axes.major);
                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width (minor axis), sigma for trail",  axes.minor);

trunk/psModules/src/objects/pmSourceIO_CMP.c

@@ -135 +135 @@
         lsky = (source->sky < 1.0) ? 0.0 : log10(source->sky);
 
-        axes = pmPSF_ModelToAxes (PAR, 20.0, model->type);
+        axes = pmPSF_ModelToAxes (PAR, model->type);
 
         float psfMagErr = isfinite(source->psfMagErr) ? source->psfMagErr : 999;

trunk/psModules/src/objects/pmSourceIO_OBJ.c

@@ -92 +92 @@
         }
 
-        axes = pmPSF_ModelToAxes (PAR, 20.0, model->type);
+        axes = pmPSF_ModelToAxes (PAR, model->type);
 
         psLineInit (line);

trunk/psModules/src/objects/pmSourceIO_PS1_CAL_0.c

@@ -114 +114 @@
             yErr = dPAR[PM_PAR_YPOS];
             if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
-                axes = pmPSF_ModelToAxes (PAR, 20.0, model->type);
+                axes = pmPSF_ModelToAxes (PAR, model->type);
             } else {
                 axes.major = NAN;
@@ -623 +623 @@
             yErr = dPAR[PM_PAR_YPOS];
 
-            axes = pmPSF_ModelToAxes (PAR, 20.0, model->type);
+            axes = pmPSF_ModelToAxes (PAR, model->type);
 
             // generate RA,DEC

trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c

@@ -90 +90 @@
             yErr = dPAR[PM_PAR_YPOS];
 
-            axes = pmPSF_ModelToAxes (PAR, 20.0, model->type);
+            axes = pmPSF_ModelToAxes (PAR, model->type);
         } else {
             // XXX: This code seg faults if source->peak is NULL.

trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c

@@ -96 +96 @@
             yErr = dPAR[PM_PAR_YPOS];
             if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
-                axes = pmPSF_ModelToAxes (PAR, 20.0, model->type);
+                axes = pmPSF_ModelToAxes (PAR, model->type);
             } else {
                 axes.major = NAN;
@@ -523 +523 @@
             yErr = dPAR[PM_PAR_YPOS];
 
-            axes = pmPSF_ModelToAxes (PAR, 20.0, model->type);
+            axes = pmPSF_ModelToAxes (PAR, model->type);
 
             row = psMetadataAlloc ();

trunk/psModules/src/objects/pmSourceIO_SMPDATA.c

@@ -92 +92 @@
             lsky = (source->sky < 1.0) ? 0.0 : log10(source->sky);
 
-            axes = pmPSF_ModelToAxes (PAR, 20.0, model->type);
+            axes = pmPSF_ModelToAxes (PAR, model->type);
 
         } else {

trunk/psModules/src/objects/pmSourceIO_SX.c

@@ -81 +81 @@
         // pmSourceSextractType (source, &type, &flags);
 
-        axes = pmPSF_ModelToAxes (PAR, 20.0, model->type);
+        axes = pmPSF_ModelToAxes (PAR, model->type);
 
         psLineInit (line);

trunk/psModules/src/objects/pmSourceOutputs.c

@@ -107 +107 @@
         }
         if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXY]) && isfinite(PAR[PM_PAR_SYY])) {
-            axes = pmPSF_ModelToAxes (PAR, 20.0, model->type);
+            axes = pmPSF_ModelToAxes (PAR, model->type);
             outputs->psfMajor = axes.major;
             outputs->psfMinor = axes.minor;

trunk/psModules/src/objects/pmSourcePlotPSFModel.c

@@ -146 +146 @@
         // force the axis ratio to be < 20.0
         psEllipseAxes axes_mnt = psEllipseMomentsToAxes (moments, 20.0);
-        psEllipseAxes axes_psf = pmPSF_ModelToAxes (PAR, 20.0, model->type);
+        psEllipseAxes axes_psf = pmPSF_ModelToAxes (PAR, model->type);
 
         // moments major axis

trunk/psModules/src/objects/pmSourceVisual.c

@@ -46 +46 @@
 bool pmSourceVisualClose() {
     if (kapa1 != -1)
-        KiiClose(kapa1);
+        KapaClose(kapa1);
     return true;
 }

trunk/psModules/src/psmodules.h

@@ -158 +158 @@
 #include <pmReadoutFake.h>
 #include <pmPSFEnvelope.h>
+#include <pmThreadTools.h>
 
 #endif