Index: trunk/psModules/src/objects/Makefile.am
===================================================================
--- trunk/psModules/src/objects/Makefile.am	(revision 14530)
+++ trunk/psModules/src/objects/Makefile.am	(revision 14652)
@@ -7,9 +7,12 @@
      pmMoments.c \
      pmModel.c \
-     pmModelGroup.c \
+     pmModelClass.c \
+     pmModelUtils.c \
      pmSource.c \
+     pmSourceUtils.c \
      pmSourceSky.c \
      pmSourceContour.c \
      pmSourceFitModel.c \
+     pmSourceFitSet.c \
      pmSourcePhotometry.c \
      pmSourceIO.c \
@@ -31,10 +34,10 @@
 
 EXTRA_DIST = \
-	models/pmModel_GAUSS.c \
-	models/pmModel_PGAUSS.c \
-	models/pmModel_QGAUSS.c \
-	models/pmModel_SGAUSS.c \
-	models/pmModel_RGAUSS.c \
-	models/pmModel_SERSIC.c
+     models/pmModel_GAUSS.c \
+     models/pmModel_PGAUSS.c \
+     models/pmModel_QGAUSS.c \
+     models/pmModel_SGAUSS.c \
+     models/pmModel_RGAUSS.c \
+     models/pmModel_SERSIC.c
 
 pkginclude_HEADERS = \
@@ -42,9 +45,12 @@
      pmMoments.h \
      pmModel.h \
-     pmModelGroup.h \
+     pmModelClass.h \
+     pmModelUtils.h \
      pmSource.h \
+     pmSourceUtils.h \
      pmSourceSky.h \
      pmSourceContour.h \
      pmSourceFitModel.h \
+     pmSourceFitSet.h \
      pmSourcePhotometry.h \
      pmSourceIO.h \
Index: trunk/psModules/src/objects/models/pmModel_GAUSS.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_GAUSS.c	(revision 14530)
+++ trunk/psModules/src/objects/models/pmModel_GAUSS.c	(revision 14652)
@@ -19,12 +19,14 @@
  *****************************************************************************/
 
-# define PM_MODEL_FUNC       pmModelFunc_GAUSS
-# define PM_MODEL_FLUX       pmModelFlux_GAUSS
-# define PM_MODEL_GUESS      pmModelGuess_GAUSS
-# define PM_MODEL_LIMITS     pmModelLimits_GAUSS
-# define PM_MODEL_RADIUS     pmModelRadius_GAUSS
-# define PM_MODEL_FROM_PSF   pmModelFromPSF_GAUSS
-# define PM_MODEL_FIT_STATUS pmModelFitStatus_GAUSS
-
+# define PM_MODEL_FUNC       	  pmModelFunc_GAUSS
+# define PM_MODEL_FLUX       	  pmModelFlux_GAUSS
+# define PM_MODEL_GUESS      	  pmModelGuess_GAUSS
+# define PM_MODEL_LIMITS     	  pmModelLimits_GAUSS
+# define PM_MODEL_RADIUS     	  pmModelRadius_GAUSS
+# define PM_MODEL_FROM_PSF   	  pmModelFromPSF_GAUSS
+# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_GAUSS
+# define PM_MODEL_FIT_STATUS      pmModelFitStatus_GAUSS
+
+// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y)
 psF32 PM_MODEL_FUNC(psVector *deriv,
                     const psVector *params,
@@ -38,4 +40,6 @@
     psF32 py = Y / PAR[PM_PAR_SYY];
     psF32 z  = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y;
+    assert (z >= 0.0);
+
     psF32 r  = exp(-z);
     psF32 q  = PAR[PM_PAR_I0]*r;
@@ -61,4 +65,5 @@
 # define AR_MAX 20.0
 # define AR_RATIO 0.99
+
 bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
 {
@@ -93,8 +98,8 @@
             break;
         case PM_PAR_SXX:
-            beta_lim = 0.5;
+            beta_lim = 2.0;
             break;
         case PM_PAR_SYY:
-            beta_lim = 0.5;
+            beta_lim = 2.0;
             break;
         case PM_PAR_SXY:
@@ -257,18 +262,17 @@
 bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf)
 {
-
     psF32 *out = modelPSF->params->data.F32;
     psF32 *in  = modelFLT->params->data.F32;
 
     // we require these two parameters to exist
-    assert (psf->params_NEW->n > PM_PAR_YPOS);
-    assert (psf->params_NEW->n > PM_PAR_XPOS);
+    assert (psf->params->n > PM_PAR_YPOS);
+    assert (psf->params->n > PM_PAR_XPOS);
 
     // supply the model-fitted parameters, or copy from the input
-    for (int i = 0; i < psf->params_NEW->n; i++) {
-        if (psf->params_NEW->data[i] == NULL) {
+    for (int i = 0; i < psf->params->n; i++) {
+        if (psf->params->data[i] == NULL) {
             out[i] = in[i];
         } else {
-            psPolynomial2D *poly = psf->params_NEW->data[i];
+            psPolynomial2D *poly = psf->params->data[i];
             out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);
         }
@@ -289,7 +293,7 @@
     // apply the model limits here: this truncates excessive extrapolation
     // XXX do we need to do this still?  should we put in asserts to test?
-    for (int i = 0; i < psf->params_NEW->n; i++) {
+    for (int i = 0; i < psf->params->n; i++) {
         // apply the limits to all components or just the psf-model parameters?
-        if (psf->params_NEW->data[i] == NULL)
+        if (psf->params->data[i] == NULL)
             continue;
 
@@ -306,4 +310,56 @@
 }
 
+// construct the PSF model from the FLT model and the psf
+// XXX is this sufficiently general do be a global function, not a pmModelClass function?
+bool PM_MODEL_PARAMS_FROM_PSF (pmModel *model, pmPSF *psf, float Xo, float Yo, float Io)
+{
+    psF32 *PAR = model->params->data.F32;
+
+    // we require these two parameters to exist
+    assert (psf->params->n > PM_PAR_YPOS);
+    assert (psf->params->n > PM_PAR_XPOS);
+
+    PAR[PM_PAR_SKY]  = 0.0;
+    PAR[PM_PAR_I0]   = Io;
+    PAR[PM_PAR_XPOS] = Xo;
+    PAR[PM_PAR_YPOS] = Yo;
+    
+    // supply the model-fitted parameters, or copy from the input
+    for (int i = 0; i < psf->params->n; i++) {
+	if (i == PM_PAR_SKY) continue;
+	if (i == PM_PAR_I0) continue;
+	if (i == PM_PAR_XPOS) continue;
+	if (i == PM_PAR_YPOS) continue;
+	psPolynomial2D *poly = psf->params->data[i];
+	assert (poly);
+	PAR[i] = psPolynomial2DEval(poly, Xo, Yo);
+    }
+
+    // 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)) {
+	psError(PM_ERR_PSF, false, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
+	return false;
+    }
+
+    // apply the model limits here: this truncates excessive extrapolation
+    // XXX do we need to do this still?  should we put in asserts to test?
+    for (int i = 0; i < psf->params->n; i++) {
+        // apply the limits to all components or just the psf-model parameters?
+        if (psf->params->data[i] == NULL)
+            continue;
+
+	bool status = true;
+        status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MIN, i, PAR, NULL);
+        status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MAX, i, PAR, NULL);
+	if (!status) {
+	    psTrace ("psModules.objects", 5, "Hitting parameter limits at (r,c) = (%.1f, %.1f)", Xo, Yo);
+	    model->flags |= PM_MODEL_STATUS_LIMITS;
+	}
+    }
+    return(true);
+}
+
 // check the status of the fitted model
 // this test is invalid if the parameters are derived
@@ -311,5 +367,4 @@
 bool PM_MODEL_FIT_STATUS (pmModel *model)
 {
-
     psF32 dP;
     bool  status;
@@ -339,3 +394,4 @@
 # undef PM_MODEL_RADIUS
 # undef PM_MODEL_FROM_PSF
+# undef PM_MODEL_PARAMS_FROM_PSF
 # undef PM_MODEL_FIT_STATUS
Index: trunk/psModules/src/objects/models/pmModel_PGAUSS.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_PGAUSS.c	(revision 14530)
+++ trunk/psModules/src/objects/models/pmModel_PGAUSS.c	(revision 14652)
@@ -19,11 +19,12 @@
  *****************************************************************************/
 
-# define PM_MODEL_FUNC       pmModelFunc_PGAUSS
-# define PM_MODEL_FLUX       pmModelFlux_PGAUSS
-# define PM_MODEL_GUESS      pmModelGuess_PGAUSS
-# define PM_MODEL_LIMITS     pmModelLimits_PGAUSS
-# define PM_MODEL_RADIUS     pmModelRadius_PGAUSS
-# define PM_MODEL_FROM_PSF   pmModelFromPSF_PGAUSS
-# define PM_MODEL_FIT_STATUS pmModelFitStatus_PGAUSS
+# define PM_MODEL_FUNC       	  pmModelFunc_PGAUSS
+# define PM_MODEL_FLUX       	  pmModelFlux_PGAUSS
+# define PM_MODEL_GUESS      	  pmModelGuess_PGAUSS
+# define PM_MODEL_LIMITS     	  pmModelLimits_PGAUSS
+# define PM_MODEL_RADIUS     	  pmModelRadius_PGAUSS
+# define PM_MODEL_FROM_PSF   	  pmModelFromPSF_PGAUSS
+# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_PGAUSS
+# define PM_MODEL_FIT_STATUS      pmModelFitStatus_PGAUSS
 
 // the model is a function of the pixel coordinate (pixcoord[0,1] = x,y)
@@ -66,12 +67,4 @@
 # define AR_RATIO 0.99
 
-// we constraint limits by the min valid minor axis:
-# define MIN_MINOR_AXIS 0.5
-
-// f3 = (s_b^-2 - s_a^-2); F3_SQ_MAX is MIN_MINOR_AXIS^-4
-# define F3_SQ_MAX 16.0 
-
-static float saveParams[8];
-
 bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
 {
@@ -150,5 +143,4 @@
             psAbort("invalid parameter %d for param min test", nParam);
         }
-	saveParams[nParam] = params[nParam];
 	if (params[nParam] < params_min) {
             params[nParam] = params_min;
@@ -184,5 +176,4 @@
             psAbort("invalid parameter %d for param max test", nParam);
         }
-	saveParams[nParam] = params[nParam];
         if (params[nParam] > params_max) {
             params[nParam] = params_max;
@@ -263,4 +254,6 @@
 psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
 {
+    psF64 z, f;
+    int Nstep = 0;
     psEllipseShape shape;
 
@@ -280,5 +273,36 @@
     // this estimates the radius assuming f(z) is roughly exp(-z)
     psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0);
-    psF64 radius = axes.major * sqrt (2.0 * log(PAR[PM_PAR_I0] / flux));
+    psF64 sigma = axes.major;
+
+    psF64 limit = flux / PAR[PM_PAR_I0];
+
+    // use the fact that f is monotonically decreasing
+    z = 0;
+    Nstep = 0;
+
+    // choose a z value guaranteed to be beyond our limit
+    float z0 = pow((1.0 / limit), (1.0 / 3.0));
+    float z1 = (1.0 / limit);
+    z1 = PS_MAX (z0, z1);
+    z0 = 0.0;
+
+    // perform a type of bisection to find the value
+    float f0 = 1.0 / (1 + z0 + z0*z0/2.0 + z0*z0*z0/6.0);
+    float f1 = 1.0 / (1 + z1 + z1*z1/2.0 + z1*z1*z1/6.0);
+    while ((Nstep < 10) && (fabs(z1 - z0) > 0.5)) {
+        z = 0.5*(z0 + z1);
+	f = 1.0 / (1 + z + z*z/2.0 + z*z*z/6.0);
+        if (f > limit) {
+            z0 = z;
+            f0 = f;
+        } else {
+            z1 = z;
+            f1 = f;
+        }
+        Nstep ++;
+    }
+    psF64 radius = sigma * sqrt (2.0 * z);
+
+    // psF64 radius = axes.major * sqrt (2.0 * log(PAR[PM_PAR_I0] / flux));
 
     if (isnan(radius))
@@ -297,12 +321,12 @@
 
     // we require these two parameters to exist
-    assert (psf->params_NEW->n > PM_PAR_YPOS);
-    assert (psf->params_NEW->n > PM_PAR_XPOS);
-
-    for (int i = 0; i < psf->params_NEW->n; i++) {
-        if (psf->params_NEW->data[i] == NULL) {
+    assert (psf->params->n > PM_PAR_YPOS);
+    assert (psf->params->n > PM_PAR_XPOS);
+
+    for (int i = 0; i < psf->params->n; i++) {
+        if (psf->params->data[i] == NULL) {
             out[i] = in[i];
         } else {
-            psPolynomial2D *poly = psf->params_NEW->data[i];
+            psPolynomial2D *poly = psf->params->data[i];
             out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);
         }
@@ -322,7 +346,7 @@
     // apply the model limits here: this truncates excessive extrapolation
     // XXX do we need to do this still?  should we put in asserts to test?
-    for (int i = 0; i < psf->params_NEW->n; i++) {
+    for (int i = 0; i < psf->params->n; i++) {
         // apply the limits to all components or just the psf-model parameters?
-        if (psf->params_NEW->data[i] == NULL)
+        if (psf->params->data[i] == NULL)
             continue;
 
@@ -339,4 +363,56 @@
 }
 
+// construct the PSF model from the FLT model and the psf
+// XXX is this sufficiently general do be a global function, not a pmModelClass function?
+bool PM_MODEL_PARAMS_FROM_PSF (pmModel *model, pmPSF *psf, float Xo, float Yo, float Io)
+{
+    psF32 *PAR = model->params->data.F32;
+
+    // we require these two parameters to exist
+    assert (psf->params->n > PM_PAR_YPOS);
+    assert (psf->params->n > PM_PAR_XPOS);
+
+    PAR[PM_PAR_SKY]  = 0.0;
+    PAR[PM_PAR_I0]   = Io;
+    PAR[PM_PAR_XPOS] = Xo;
+    PAR[PM_PAR_YPOS] = Yo;
+    
+    // supply the model-fitted parameters, or copy from the input
+    for (int i = 0; i < psf->params->n; i++) {
+	if (i == PM_PAR_SKY) continue;
+	if (i == PM_PAR_I0) continue;
+	if (i == PM_PAR_XPOS) continue;
+	if (i == PM_PAR_YPOS) continue;
+	psPolynomial2D *poly = psf->params->data[i];
+	assert (poly);
+	PAR[i] = psPolynomial2DEval(poly, Xo, Yo);
+    }
+
+    // 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)) {
+	psError(PM_ERR_PSF, false, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
+	return false;
+    }
+
+    // apply the model limits here: this truncates excessive extrapolation
+    // XXX do we need to do this still?  should we put in asserts to test?
+    for (int i = 0; i < psf->params->n; i++) {
+        // apply the limits to all components or just the psf-model parameters?
+        if (psf->params->data[i] == NULL)
+            continue;
+
+	bool status = true;
+        status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MIN, i, PAR, NULL);
+        status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MAX, i, PAR, NULL);
+	if (!status) {
+	    psTrace ("psModules.objects", 5, "Hitting parameter limits at (r,c) = (%.1f, %.1f)", Xo, Yo);
+	    model->flags |= PM_MODEL_STATUS_LIMITS;
+	}
+    }
+    return(true);
+}
+
 bool PM_MODEL_FIT_STATUS (pmModel *model)
 {
@@ -367,3 +443,4 @@
 # undef PM_MODEL_RADIUS
 # undef PM_MODEL_FROM_PSF
+# undef PM_MODEL_PARAMS_FROM_PSF
 # undef PM_MODEL_FIT_STATUS
Index: trunk/psModules/src/objects/models/pmModel_QGAUSS.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_QGAUSS.c	(revision 14530)
+++ trunk/psModules/src/objects/models/pmModel_QGAUSS.c	(revision 14652)
@@ -20,11 +20,12 @@
    *****************************************************************************/
 
-# define PM_MODEL_FUNC       pmModelFunc_QGAUSS
-# define PM_MODEL_FLUX       pmModelFlux_QGAUSS
-# define PM_MODEL_GUESS      pmModelGuess_QGAUSS
-# define PM_MODEL_LIMITS     pmModelLimits_QGAUSS
-# define PM_MODEL_RADIUS     pmModelRadius_QGAUSS
-# define PM_MODEL_FROM_PSF   pmModelFromPSF_QGAUSS
-# define PM_MODEL_FIT_STATUS pmModelFitStatus_QGAUSS
+# define PM_MODEL_FUNC       	  pmModelFunc_QGAUSS
+# define PM_MODEL_FLUX       	  pmModelFlux_QGAUSS
+# define PM_MODEL_GUESS      	  pmModelGuess_QGAUSS
+# define PM_MODEL_LIMITS     	  pmModelLimits_QGAUSS
+# define PM_MODEL_RADIUS     	  pmModelRadius_QGAUSS
+# define PM_MODEL_FROM_PSF   	  pmModelFromPSF_QGAUSS
+# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_QGAUSS
+# define PM_MODEL_FIT_STATUS      pmModelFitStatus_QGAUSS
 
 psF32 PM_MODEL_FUNC (psVector *deriv,
@@ -350,12 +351,12 @@
 
     // we require these two parameters to exist
-    assert (psf->params_NEW->n > PM_PAR_YPOS);
-    assert (psf->params_NEW->n > PM_PAR_XPOS);
-
-    for (int i = 0; i < psf->params_NEW->n; i++) {
-        if (psf->params_NEW->data[i] == NULL) {
+    assert (psf->params->n > PM_PAR_YPOS);
+    assert (psf->params->n > PM_PAR_XPOS);
+
+    for (int i = 0; i < psf->params->n; i++) {
+        if (psf->params->data[i] == NULL) {
             out[i] = in[i];
         } else {
-            psPolynomial2D *poly = psf->params_NEW->data[i];
+            psPolynomial2D *poly = psf->params->data[i];
             out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);
         }
@@ -372,7 +373,7 @@
     // apply the model limits here: this truncates excessive extrapolation
     // XXX do we need to do this still?  should we put in asserts to test?
-    for (int i = 0; i < psf->params_NEW->n; i++) {
+    for (int i = 0; i < psf->params->n; i++) {
         // apply the limits to all components or just the psf-model parameters?
-        if (psf->params_NEW->data[i] == NULL)
+        if (psf->params->data[i] == NULL)
             continue;
 
@@ -390,4 +391,56 @@
 }
 
+// construct the PSF model from the FLT model and the psf
+// XXX is this sufficiently general do be a global function, not a pmModelClass function?
+bool PM_MODEL_PARAMS_FROM_PSF (pmModel *model, pmPSF *psf, float Xo, float Yo, float Io)
+{
+    psF32 *PAR = model->params->data.F32;
+
+    // we require these two parameters to exist
+    assert (psf->params->n > PM_PAR_YPOS);
+    assert (psf->params->n > PM_PAR_XPOS);
+
+    PAR[PM_PAR_SKY]  = 0.0;
+    PAR[PM_PAR_I0]   = Io;
+    PAR[PM_PAR_XPOS] = Xo;
+    PAR[PM_PAR_YPOS] = Yo;
+    
+    // supply the model-fitted parameters, or copy from the input
+    for (int i = 0; i < psf->params->n; i++) {
+	if (i == PM_PAR_SKY) continue;
+	if (i == PM_PAR_I0) continue;
+	if (i == PM_PAR_XPOS) continue;
+	if (i == PM_PAR_YPOS) continue;
+	psPolynomial2D *poly = psf->params->data[i];
+	assert (poly);
+	PAR[i] = psPolynomial2DEval(poly, Xo, Yo);
+    }
+
+    // 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)) {
+	psError(PM_ERR_PSF, false, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
+	return false;
+    }
+
+    // apply the model limits here: this truncates excessive extrapolation
+    // XXX do we need to do this still?  should we put in asserts to test?
+    for (int i = 0; i < psf->params->n; i++) {
+        // apply the limits to all components or just the psf-model parameters?
+        if (psf->params->data[i] == NULL)
+            continue;
+
+	bool status = true;
+        status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MIN, i, PAR, NULL);
+        status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MAX, i, PAR, NULL);
+	if (!status) {
+	    psTrace ("psModules.objects", 5, "Hitting parameter limits at (r,c) = (%.1f, %.1f)", Xo, Yo);
+	    model->flags |= PM_MODEL_STATUS_LIMITS;
+	}
+    }
+    return(true);
+}
+
 bool PM_MODEL_FIT_STATUS (pmModel *model)
 {
@@ -418,3 +471,4 @@
 # undef PM_MODEL_RADIUS
 # undef PM_MODEL_FROM_PSF
+# undef PM_MODEL_PARAMS_FROM_PSF
 # undef PM_MODEL_FIT_STATUS
Index: trunk/psModules/src/objects/models/pmModel_RGAUSS.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_RGAUSS.c	(revision 14530)
+++ trunk/psModules/src/objects/models/pmModel_RGAUSS.c	(revision 14652)
@@ -20,11 +20,12 @@
  *****************************************************************************/
 
-# define PM_MODEL_FUNC       pmModelFunc_RGAUSS
-# define PM_MODEL_FLUX       pmModelFlux_RGAUSS
-# define PM_MODEL_GUESS      pmModelGuess_RGAUSS
-# define PM_MODEL_LIMITS     pmModelLimits_RGAUSS
-# define PM_MODEL_RADIUS     pmModelRadius_RGAUSS
-# define PM_MODEL_FROM_PSF   pmModelFromPSF_RGAUSS
-# define PM_MODEL_FIT_STATUS pmModelFitStatus_RGAUSS
+# define PM_MODEL_FUNC       	  pmModelFunc_RGAUSS
+# define PM_MODEL_FLUX       	  pmModelFlux_RGAUSS
+# define PM_MODEL_GUESS      	  pmModelGuess_RGAUSS
+# define PM_MODEL_LIMITS     	  pmModelLimits_RGAUSS
+# define PM_MODEL_RADIUS     	  pmModelRadius_RGAUSS
+# define PM_MODEL_FROM_PSF   	  pmModelFromPSF_RGAUSS
+# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_RGAUSS
+# define PM_MODEL_FIT_STATUS      pmModelFitStatus_RGAUSS
 
 psF32 PM_MODEL_FUNC (psVector *deriv,
@@ -343,12 +344,12 @@
 
     // we require these two parameters to exist
-    assert (psf->params_NEW->n > PM_PAR_YPOS);
-    assert (psf->params_NEW->n > PM_PAR_XPOS);
-
-    for (int i = 0; i < psf->params_NEW->n; i++) {
-        if (psf->params_NEW->data[i] == NULL) {
+    assert (psf->params->n > PM_PAR_YPOS);
+    assert (psf->params->n > PM_PAR_XPOS);
+
+    for (int i = 0; i < psf->params->n; i++) {
+        if (psf->params->data[i] == NULL) {
             out[i] = in[i];
         } else {
-            psPolynomial2D *poly = psf->params_NEW->data[i];
+            psPolynomial2D *poly = psf->params->data[i];
             out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);
         }
@@ -365,7 +366,7 @@
     // apply the model limits here: this truncates excessive extrapolation
     // XXX do we need to do this still?  should we put in asserts to test?
-    for (int i = 0; i < psf->params_NEW->n; i++) {
+    for (int i = 0; i < psf->params->n; i++) {
         // apply the limits to all components or just the psf-model parameters?
-        if (psf->params_NEW->data[i] == NULL)
+        if (psf->params->data[i] == NULL)
             continue;
 
@@ -383,4 +384,56 @@
 }
 
+// construct the PSF model from the FLT model and the psf
+// XXX is this sufficiently general do be a global function, not a pmModelClass function?
+bool PM_MODEL_PARAMS_FROM_PSF (pmModel *model, pmPSF *psf, float Xo, float Yo, float Io)
+{
+    psF32 *PAR = model->params->data.F32;
+
+    // we require these two parameters to exist
+    assert (psf->params->n > PM_PAR_YPOS);
+    assert (psf->params->n > PM_PAR_XPOS);
+
+    PAR[PM_PAR_SKY]  = 0.0;
+    PAR[PM_PAR_I0]   = Io;
+    PAR[PM_PAR_XPOS] = Xo;
+    PAR[PM_PAR_YPOS] = Yo;
+    
+    // supply the model-fitted parameters, or copy from the input
+    for (int i = 0; i < psf->params->n; i++) {
+	if (i == PM_PAR_SKY) continue;
+	if (i == PM_PAR_I0) continue;
+	if (i == PM_PAR_XPOS) continue;
+	if (i == PM_PAR_YPOS) continue;
+	psPolynomial2D *poly = psf->params->data[i];
+	assert (poly);
+	PAR[i] = psPolynomial2DEval(poly, Xo, Yo);
+    }
+
+    // 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)) {
+	psError(PM_ERR_PSF, false, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
+	return false;
+    }
+
+    // apply the model limits here: this truncates excessive extrapolation
+    // XXX do we need to do this still?  should we put in asserts to test?
+    for (int i = 0; i < psf->params->n; i++) {
+        // apply the limits to all components or just the psf-model parameters?
+        if (psf->params->data[i] == NULL)
+            continue;
+
+	bool status = true;
+        status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MIN, i, PAR, NULL);
+        status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MAX, i, PAR, NULL);
+	if (!status) {
+	    psTrace ("psModules.objects", 5, "Hitting parameter limits at (r,c) = (%.1f, %.1f)", Xo, Yo);
+	    model->flags |= PM_MODEL_STATUS_LIMITS;
+	}
+    }
+    return(true);
+}
+
 bool PM_MODEL_FIT_STATUS (pmModel *model)
 {
@@ -411,3 +464,4 @@
 # undef PM_MODEL_RADIUS
 # undef PM_MODEL_FROM_PSF
+# undef PM_MODEL_PARAMS_FROM_PSF
 # undef PM_MODEL_FIT_STATUS
Index: trunk/psModules/src/objects/models/pmModel_SERSIC.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_SERSIC.c	(revision 14530)
+++ trunk/psModules/src/objects/models/pmModel_SERSIC.c	(revision 14652)
@@ -23,11 +23,12 @@
    *****************************************************************************/
 
-# define PM_MODEL_FUNC       pmModelFunc_SERSIC
-# define PM_MODEL_FLUX       pmModelFlux_SERSIC
-# define PM_MODEL_GUESS      pmModelGuess_SERSIC
-# define PM_MODEL_LIMITS     pmModelLimits_SERSIC
-# define PM_MODEL_RADIUS     pmModelRadius_SERSIC
-# define PM_MODEL_FROM_PSF   pmModelFromPSF_SERSIC
-# define PM_MODEL_FIT_STATUS pmModelFitStatus_SERSIC
+# define PM_MODEL_FUNC       	  pmModelFunc_SERSIC
+# define PM_MODEL_FLUX       	  pmModelFlux_SERSIC
+# define PM_MODEL_GUESS      	  pmModelGuess_SERSIC
+# define PM_MODEL_LIMITS     	  pmModelLimits_SERSIC
+# define PM_MODEL_RADIUS     	  pmModelRadius_SERSIC
+# define PM_MODEL_FROM_PSF   	  pmModelFromPSF_SERSIC
+# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_SERSIC
+# define PM_MODEL_FIT_STATUS      pmModelFitStatus_SERSIC
 
 psF32 PM_MODEL_FUNC (psVector *deriv,
@@ -360,12 +361,12 @@
 
     // we require these two parameters to exist
-    assert (psf->params_NEW->n > PM_PAR_YPOS);
-    assert (psf->params_NEW->n > PM_PAR_XPOS);
-
-    for (int i = 0; i < psf->params_NEW->n; i++) {
-        if (psf->params_NEW->data[i] == NULL) {
+    assert (psf->params->n > PM_PAR_YPOS);
+    assert (psf->params->n > PM_PAR_XPOS);
+
+    for (int i = 0; i < psf->params->n; i++) {
+        if (psf->params->data[i] == NULL) {
             out[i] = in[i];
         } else {
-            psPolynomial2D *poly = psf->params_NEW->data[i];
+            psPolynomial2D *poly = psf->params->data[i];
             out[i] = psPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);
         }
@@ -382,7 +383,7 @@
     // apply the model limits here: this truncates excessive extrapolation
     // XXX do we need to do this still?  should we put in asserts to test?
-    for (int i = 0; i < psf->params_NEW->n; i++) {
+    for (int i = 0; i < psf->params->n; i++) {
         // apply the limits to all components or just the psf-model parameters?
-        if (psf->params_NEW->data[i] == NULL)
+        if (psf->params->data[i] == NULL)
             continue;
 
@@ -400,4 +401,56 @@
 }
 
+// construct the PSF model from the FLT model and the psf
+// XXX is this sufficiently general do be a global function, not a pmModelClass function?
+bool PM_MODEL_PARAMS_FROM_PSF (pmModel *model, pmPSF *psf, float Xo, float Yo, float Io)
+{
+    psF32 *PAR = model->params->data.F32;
+
+    // we require these two parameters to exist
+    assert (psf->params->n > PM_PAR_YPOS);
+    assert (psf->params->n > PM_PAR_XPOS);
+
+    PAR[PM_PAR_SKY]  = 0.0;
+    PAR[PM_PAR_I0]   = Io;
+    PAR[PM_PAR_XPOS] = Xo;
+    PAR[PM_PAR_YPOS] = Yo;
+    
+    // supply the model-fitted parameters, or copy from the input
+    for (int i = 0; i < psf->params->n; i++) {
+	if (i == PM_PAR_SKY) continue;
+	if (i == PM_PAR_I0) continue;
+	if (i == PM_PAR_XPOS) continue;
+	if (i == PM_PAR_YPOS) continue;
+	psPolynomial2D *poly = psf->params->data[i];
+	assert (poly);
+	PAR[i] = psPolynomial2DEval(poly, Xo, Yo);
+    }
+
+    // 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)) {
+	psError(PM_ERR_PSF, false, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
+	return false;
+    }
+
+    // apply the model limits here: this truncates excessive extrapolation
+    // XXX do we need to do this still?  should we put in asserts to test?
+    for (int i = 0; i < psf->params->n; i++) {
+        // apply the limits to all components or just the psf-model parameters?
+        if (psf->params->data[i] == NULL)
+            continue;
+
+	bool status = true;
+        status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MIN, i, PAR, NULL);
+        status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MAX, i, PAR, NULL);
+	if (!status) {
+	    psTrace ("psModules.objects", 5, "Hitting parameter limits at (r,c) = (%.1f, %.1f)", Xo, Yo);
+	    model->flags |= PM_MODEL_STATUS_LIMITS;
+	}
+    }
+    return(true);
+}
+
 bool PM_MODEL_FIT_STATUS (pmModel *model)
 {
@@ -428,3 +481,4 @@
 # undef PM_MODEL_RADIUS
 # undef PM_MODEL_FROM_PSF
+# undef PM_MODEL_PARAMS_FROM_PSF
 # undef PM_MODEL_FIT_STATUS
Index: trunk/psModules/src/objects/pmGrowthCurve.c
===================================================================
--- trunk/psModules/src/objects/pmGrowthCurve.c	(revision 14530)
+++ trunk/psModules/src/objects/pmGrowthCurve.c	(revision 14652)
@@ -5,6 +5,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -23,8 +23,8 @@
 #include "pmMoments.h"
 #include "pmResiduals.h"
+#include "pmGrowthCurve.h"
+#include "pmPSF.h"
 #include "pmModel.h"
 #include "pmSource.h"
-#include "pmGrowthCurve.h"
-#include "pmPSF.h"
 #include "psVectorBracket.h"
 
Index: trunk/psModules/src/objects/pmModel.c
===================================================================
--- trunk/psModules/src/objects/pmModel.c	(revision 14530)
+++ trunk/psModules/src/objects/pmModel.c	(revision 14652)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-06-20 02:22:26 $
+ *  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -21,16 +21,11 @@
 #include <string.h>
 #include <pslib.h>
+#include "pmHDU.h"
+#include "pmFPA.h"
 #include "pmResiduals.h"
+#include "pmGrowthCurve.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-
-/* The following types and functions need to be known by pmModel.c and are defined in
-   pmModelGroup.h.  The use of pmSource and other types in pmModelGroup.h prevent us from
-   directly including it here  ***/
-
-typedef psMinimizeLMChi2Func pmModelFunc;
-typedef bool (*pmModelFitStatusFunc)(pmModel *model);
-pmModelFunc pmModelFunc_GetFunction (pmModelType type);
-pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction (pmModelType type);
-int pmModelParameterCount(pmModelType type);
+#include "pmModelClass.h"
 
 static void modelFree(pmModel *tmp)
@@ -45,5 +40,4 @@
 pmModelAlloc(): Allocate the pmModel structure, along with its parameters,
 and initialize the type member.  Initialize the params to 0.0.
-XXX EAM: simplifying code with pmModelParameterCount
 *****************************************************************************/
 pmModel *pmModelAlloc(pmModelType type)
@@ -53,4 +47,10 @@
     pmModel *tmp = (pmModel *) psAlloc(sizeof(pmModel));
     psMemSetDeallocator(tmp, (psFreeFunc) modelFree);
+
+    pmModelClass *class = pmModelClassSelect (type);
+    if (class == NULL) {
+        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
+        return(NULL);
+    }
 
     tmp->type = type;
@@ -62,12 +62,11 @@
     tmp->residuals = NULL;              // XXX should the model own this memory?
 
-    psS32 Nparams = pmModelParameterCount(type);
-    if (Nparams == 0) {
-        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
-        return(NULL);
-    }
+    psS32 Nparams = pmModelClassParameterCount(type);
+    assert (Nparams);
 
     tmp->params  = psVectorAlloc(Nparams, PS_TYPE_F32);
     tmp->dparams = psVectorAlloc(Nparams, PS_TYPE_F32);
+    assert (tmp->params);
+    assert (tmp->dparams);
 
     for (psS32 i = 0; i < tmp->params->n; i++) {
@@ -75,4 +74,13 @@
         tmp->dparams->data.F32[i] = 0.0;
     }
+
+    tmp->modelFunc   	    = class->modelFunc;
+    tmp->modelFlux   	    = class->modelFlux;
+    tmp->modelRadius 	    = class->modelRadius;
+    tmp->modelLimits 	    = class->modelLimits;
+    tmp->modelGuess         = class->modelGuess;
+    tmp->modelFromPSF       = class->modelFromPSF;
+    tmp->modelParamsFromPSF = class->modelParamsFromPSF;
+    tmp->modelFitStatus     = class->modelFitStatus;
 
     psTrace("psModules.objects", 3, "---- %s() end ----\n", __func__);
@@ -92,11 +100,11 @@
     new->radiusFit = model->radiusFit;
 
-    // XXX note that model->residuals is just a reference
-    new->residuals = model->residuals;
-
     for (int i = 0; i < new->params->n; i++) {
         new->params->data.F32[i]  = model->params->data.F32[i];
         new->dparams->data.F32[i] = model->dparams->data.F32[i];
     }
+
+    // note that model->residuals is just a reference
+    new->residuals = model->residuals;
 
     return (new);
@@ -123,8 +131,26 @@
     x->data.F32[1] = (psF32) (row + image->row0);
     psF32 tmpF;
-    pmModelFunc modelFunc;
-
-    modelFunc = pmModelFunc_GetFunction (model->type);
-    tmpF = modelFunc (NULL, model->params, x);
+
+    tmpF = model->modelFunc (NULL, model->params, x);
+    psFree(x);
+    psTrace("psModules.objects", 3, "---- %s() end ----\n", __func__);
+    return(tmpF);
+}
+
+psF32 pmModelEvalWithOffset(pmModel *model, psImage *image, psS32 col, psS32 row, int dx, int dy)
+{
+    psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
+    PS_ASSERT_PTR_NON_NULL(image, false);
+    PS_ASSERT_PTR_NON_NULL(model, false);
+    PS_ASSERT_PTR_NON_NULL(model->params, false);
+
+    // Allocate the x coordinate structure and convert row/col to image space.
+    //
+    psVector *x = psVectorAlloc(2, PS_TYPE_F32);
+    x->data.F32[0] = (psF32) (col + image->col0 + dx);
+    x->data.F32[1] = (psF32) (row + image->row0 + dy);
+    psF32 tmpF;
+
+    tmpF = model->modelFunc (NULL, model->params, x);
     psFree(x);
     psTrace("psModules.objects", 3, "---- %s() end ----\n", __func__);
@@ -137,5 +163,7 @@
                           pmModelOpMode mode,
                           bool add,
-                          psMaskType maskVal
+                          psMaskType maskVal,
+			  int dx,
+			  int dy
     )
 {
@@ -148,5 +176,5 @@
     psVector *x = psVectorAlloc(2, PS_TYPE_F32);
     psVector *params = model->params;
-    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
+
     psS32 imageCol;
     psS32 imageRow;
@@ -209,6 +237,6 @@
             // Convert i/j to image coord space:
             // XXX should we use using 0.5 pixel offset?
-            imageCol = ix + image->col0;
-            imageRow = iy + image->row0;
+            imageCol = ix + image->col0 + dx;
+            imageRow = iy + image->row0 + dy;
 
             x->data.F32[0] = (float) imageCol;
@@ -219,5 +247,5 @@
             // add in the desired components for this coordinate
             if (mode & PM_MODEL_OP_FUNC) {
-                pixelValue += modelFunc (NULL, params, x);
+                pixelValue += model->modelFunc (NULL, params, x);
             }
 
@@ -280,5 +308,5 @@
 {
     psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
-    psBool rc = AddOrSubModel(image, mask, model, mode, true, maskVal);
+    psBool rc = AddOrSubModel(image, mask, model, mode, true, maskVal, 0.0, 0.0);
     psTrace("psModules.objects", 3, "---- %s(%d) end ----\n", __func__, rc);
     return(rc);
@@ -294,19 +322,38 @@
 {
     psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
-    psBool rc = AddOrSubModel(image, mask, model, mode, false, maskVal);
+    psBool rc = AddOrSubModel(image, mask, model, mode, false, maskVal, 0.0, 0.0);
     psTrace("psModules.objects", 3, "---- %s(%d) end ----\n", __func__, rc);
     return(rc);
 }
 
-// check for success of model fit
-bool pmModelFitStatus (pmModel *model)
-{
-
-    bool status;
-
-    pmModelFitStatusFunc statusFunc = pmModelFitStatusFunc_GetFunction (model->type);
-    status = statusFunc (model);
-
-    return (status);
-}
-
+/******************************************************************************
+ *****************************************************************************/
+bool pmModelAddWithOffset(psImage *image,
+			  psImage *mask,
+			  pmModel *model,
+			  pmModelOpMode mode,
+			  psMaskType maskVal,
+			  int dx,
+			  int dy)
+{
+    psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
+    psBool rc = AddOrSubModel(image, mask, model, mode, true, maskVal, dx, dy);
+    psTrace("psModules.objects", 3, "---- %s(%d) end ----\n", __func__, rc);
+    return(rc);
+}
+
+/******************************************************************************
+ *****************************************************************************/
+bool pmModelSubWithOffset(psImage *image,
+			  psImage *mask,
+			  pmModel *model,
+			  pmModelOpMode mode,
+			  psMaskType maskVal,
+			  int dx,
+			  int dy)
+{
+    psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
+    psBool rc = AddOrSubModel(image, mask, model, mode, false, maskVal, dx, dy);
+    psTrace("psModules.objects", 3, "---- %s(%d) end ----\n", __func__, rc);
+    return(rc);
+}
Index: trunk/psModules/src/objects/pmModel.h
===================================================================
--- trunk/psModules/src/objects/pmModel.h	(revision 14530)
+++ trunk/psModules/src/objects/pmModel.h	(revision 14652)
@@ -1,3 +1,3 @@
-/* @file  pmObjects.h
+/* @file  pmModel.h
  * @brief Functions to define and manipulate object models
  *
@@ -5,6 +5,6 @@
  * @author EAM, IfA
  *
- * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-06-20 02:22:26 $
+ * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-08-24 00:11:02 $
  *
  * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -17,6 +17,6 @@
 /// @{
 
-// type of model carried by the pmModel structure
-typedef int pmModelType;
+/* pointers for the functions types below are supplied to each pmModel, and can be used by
+   the programmer without needing to know the model class */
 
 typedef enum {
@@ -41,4 +41,35 @@
 } pmModelOpMode;
 
+typedef struct pmModel pmModel;
+typedef struct pmSource pmSource;
+
+//  This function is the model chi-square minimization function for this model.
+typedef psMinimizeLMChi2Func pmModelFunc;
+
+//  This function sets the parameter limits for this model.
+typedef psMinimizeLMLimitFunc pmModelLimits;
+
+// This function returns the integrated flux for the given model parameters.
+typedef psF64 (*pmModelFlux)(const psVector *params);
+
+// This function returns the radius at which the given model and parameters
+// achieves the given flux.
+typedef psF64 (*pmModelRadius)(const psVector *params, double flux);
+
+//  This function provides the model guess parameters based on the details of
+//  the given source.
+typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source);
+
+//  This function constructs the PSF model for the given source based on the
+//  supplied psf and the EXT model for the object.
+typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelEXT, pmPSF *psf);
+
+//  This function sets the model parameters based on the PSF for a given coordinate and central
+//  intensity
+typedef bool (*pmModelParamsFromPSF)(pmModel *model, pmPSF *psf, float Xo, float Yo, float Io);
+
+//  This function returns the success / failure status of the given model fit
+typedef bool (*pmModelFitStatusFunc)(pmModel *model);
+
 /** pmModel data structure
  *
@@ -50,6 +81,5 @@
  *
  */
-typedef struct
-{
+struct pmModel {
     pmModelType type;                   ///< Model to be used.
     psVector *params;                   ///< Paramater values.
@@ -62,11 +92,15 @@
     float radiusFit;                    ///< fit radius actually used
     pmResiduals *residuals;             ///< normalized PSF residuals
-}
-pmModel;
 
-/* XXX we are currently saving the residuals with the pmModel.  It might be better to save this
- * in the pmSource.  we may want an API to Add/Sub a pmModel (analytical model only) or a
- * pmSource (analytical + residuals).
- */
+    // functions for this model which depend on the model class
+    pmModelFunc          modelFunc;
+    pmModelFlux          modelFlux;
+    pmModelRadius        modelRadius;
+    pmModelLimits        modelLimits;
+    pmModelGuessFunc     modelGuess;
+    pmModelFromPSFFunc   modelFromPSF;
+    pmModelParamsFromPSF modelParamsFromPSF;
+    pmModelFitStatusFunc modelFitStatus;
+};
 
 /** Symbolic names for the elements of [d]params
@@ -130,4 +164,20 @@
 );
 
+bool pmModelAddWithOffset(psImage *image,
+			  psImage *mask,
+			  pmModel *model,
+			  pmModelOpMode mode,
+			  psMaskType maskVal,
+			  int dx,
+			  int dy);
+
+bool pmModelSubWithOffset(psImage *image,
+			  psImage *mask,
+			  pmModel *model,
+			  pmModelOpMode mode,
+			  psMaskType maskVal,
+			  int dx,
+			  int dy);
+
 /** pmModelFitStatus()
  *
Index: trunk/psModules/src/objects/pmModelClass.c
===================================================================
--- trunk/psModules/src/objects/pmModelClass.c	(revision 14652)
+++ trunk/psModules/src/objects/pmModelClass.c	(revision 14652)
@@ -0,0 +1,141 @@
+/** @file  pmModelClass.c
+ *
+ *  Functions to define and manipulate object model attributes
+ *
+ *  @author GLG, MHPCC
+ *  @author EAM, IfA
+ *
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
+ *
+ *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <pslib.h>
+#include "pmHDU.h"
+#include "pmFPA.h"
+#include "pmPeaks.h"
+#include "pmMoments.h"
+#include "pmResiduals.h"
+#include "pmGrowthCurve.h"
+#include "pmPSF.h"
+#include "pmModel.h"
+#include "pmSource.h"
+#include "pmModelClass.h"
+#include "pmErrorCodes.h"
+
+// XXX shouldn't these be defined for us in pslib.h ???
+double hypot(double x, double y);
+double sqrt (double x);
+
+# include "models/pmModel_GAUSS.c"
+# include "models/pmModel_PGAUSS.c"
+# include "models/pmModel_QGAUSS.c"
+# include "models/pmModel_RGAUSS.c"
+# include "models/pmModel_SERSIC.c"
+
+static pmModelClass defaultModels[] = {
+    {"PS_MODEL_GAUSS",        7, pmModelFunc_GAUSS,   pmModelFlux_GAUSS,   pmModelRadius_GAUSS,   pmModelLimits_GAUSS,   pmModelGuess_GAUSS,  pmModelFromPSF_GAUSS,  pmModelParamsFromPSF_GAUSS,  pmModelFitStatus_GAUSS},
+    {"PS_MODEL_PGAUSS",       7, pmModelFunc_PGAUSS,  pmModelFlux_PGAUSS,  pmModelRadius_PGAUSS,  pmModelLimits_PGAUSS,  pmModelGuess_PGAUSS, pmModelFromPSF_PGAUSS, pmModelParamsFromPSF_PGAUSS, pmModelFitStatus_PGAUSS},
+    {"PS_MODEL_QGAUSS",       8, pmModelFunc_QGAUSS,  pmModelFlux_QGAUSS,  pmModelRadius_QGAUSS,  pmModelLimits_QGAUSS,  pmModelGuess_QGAUSS, pmModelFromPSF_QGAUSS, pmModelParamsFromPSF_QGAUSS, pmModelFitStatus_QGAUSS},
+    {"PS_MODEL_RGAUSS",       8, pmModelFunc_RGAUSS,  pmModelFlux_RGAUSS,  pmModelRadius_RGAUSS,  pmModelLimits_RGAUSS,  pmModelGuess_RGAUSS, pmModelFromPSF_RGAUSS, pmModelParamsFromPSF_RGAUSS, pmModelFitStatus_RGAUSS},
+    {"PS_MODEL_SERSIC",       8, pmModelFunc_SERSIC,  pmModelFlux_SERSIC,  pmModelRadius_SERSIC,  pmModelLimits_SERSIC,  pmModelGuess_SERSIC, pmModelFromPSF_SERSIC, pmModelParamsFromPSF_SERSIC, pmModelFitStatus_SERSIC}
+};
+
+static pmModelClass *models = NULL;
+static int Nmodels = 0;
+
+static void ModelClassFree (pmModelClass *modelClass)
+{
+
+    if (modelClass == NULL)
+        return;
+    return;
+}
+
+pmModelClass *pmModelClassAlloc (int nModels)
+{
+    pmModelClass *modelClass = (pmModelClass *) psAlloc (nModels * sizeof(pmModelClass));
+    psMemSetDeallocator(modelClass, (psFreeFunc) ModelClassFree);
+    return (modelClass);
+}
+
+void pmModelClassAdd (pmModelClass *model)
+{
+    if (models == NULL) {
+        pmModelClassInit ();
+    }
+
+    Nmodels ++;
+    models = (pmModelClass *) psRealloc (models, Nmodels*sizeof(pmModelClass));
+    models[Nmodels-1] = model[0];
+    return;
+}
+
+bool pmModelClassInit (void)
+{
+    // if we do not need to init, return false;
+    if (models != NULL)
+        return false;
+
+    int Nnew = sizeof (defaultModels) / sizeof (pmModelClass);
+
+    models = pmModelClassAlloc (Nnew);
+    for (int i = 0; i < Nnew; i++) {
+        models[i] = defaultModels[i];
+    }
+    Nmodels = Nnew;
+    return true;
+}
+
+pmModelClass *pmModelClassSelect (pmModelType type)
+{
+    if ((type < 0) || (type >= Nmodels)) {
+        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
+        return (NULL);
+    }
+    return (&models[type]);
+}
+
+void pmModelClassCleanup (void)
+{
+    psFree (models);
+    models = NULL;
+    return;
+}
+
+psS32 pmModelClassParameterCount (pmModelType type)
+{
+    if ((type < 0) || (type >= Nmodels)) {
+        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
+        return (0);
+    }
+    return (models[type].nParams);
+}
+
+psS32 pmModelClassGetType (char *name)
+{
+    for (int i = 0; i < Nmodels; i++) {
+        if (!strcmp(models[i].name, name)) {
+            return (i);
+        }
+    }
+    return (-1);
+}
+
+char *pmModelClassGetName (pmModelType type)
+{
+    if ((type < 0) || (type >= Nmodels)) {
+        psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
+        return (NULL);
+    }
+    return (models[type].name);
+}
Index: trunk/psModules/src/objects/pmModelClass.h
===================================================================
--- trunk/psModules/src/objects/pmModelClass.h	(revision 14652)
+++ trunk/psModules/src/objects/pmModelClass.h	(revision 14652)
@@ -0,0 +1,73 @@
+/* @file  pmModelClass.h
+ *
+ * The object model function types are desined to allow for the flexible addition of new
+ * object models. Every object model, with parameters represented by pmModel, has an
+ * associated set of functions which provide necessary support operations.  A These
+ * functions allow the programmer to select the approriate function or property for a
+ * specific object model class.
+ *
+ * Every model instance belongs to a class of models, defined by the value of the
+ * pmModelType type entry. Various functions need access to information about each of the
+ * models. Some of this information varies from model to model, and may depend on the
+ * current parameter values or other data quantities. In order to keep the code from
+ * requiring the information about each model to be coded into the low-level fitting
+ * routines, we define a collection of functions which allow us to abstract this type of
+ * model-dependent information. These generic functions take the model type and return the
+ * corresponding function pointer for the specified model. Each model is defined by
+ * creating this collection of specific functions, and placing them in a single file for
+ * each model. We define the following structure to carry the collection of information
+ * about the models.
+ *
+ * @author EAM, IfA
+ *
+ * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-08-24 00:11:02 $
+ * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ */
+
+# ifndef PM_MODEL_CLASS_H
+# define PM_MODEL_CLASS_H
+
+/// @addtogroup Objects Object Detection / Analysis Functions
+/// @{
+
+typedef struct
+{
+    char *name;
+    int nParams;
+    pmModelFunc          modelFunc;
+    pmModelFlux          modelFlux;
+    pmModelRadius        modelRadius;
+    pmModelLimits        modelLimits;
+    pmModelGuessFunc     modelGuess;
+    pmModelFromPSFFunc   modelFromPSF;
+    pmModelParamsFromPSF modelParamsFromPSF;
+    pmModelFitStatusFunc modelFitStatus;
+} pmModelClass;
+
+// allocate a pmModelClass to hold nModels entries
+pmModelClass *pmModelClassAlloc (int nModels);
+
+// initialize the internal (static) model class with the default models
+bool pmModelClassInit (void);
+
+// free the internal (static) model class
+void pmModelClassCleanup (void);
+
+// add a new model class to the collection of model classes
+void pmModelClassAdd (pmModelClass *modelClass);
+
+// get the specified model class
+pmModelClass *pmModelClassSelect (pmModelType type);
+
+// This function returns the number of parameters used by the listed function.
+int pmModelClassParameterCount (pmModelType type);
+
+// This function returns the user-space model names for the specified model type.
+char *pmModelClassGetName (pmModelType type);
+
+// This function returns the internal model type code for the user-space model names.
+pmModelType pmModelClassGetType (char *name);
+
+/// @}
+# endif /* PM_MODEL_CLASS_H */
Index: trunk/psModules/src/objects/pmModelGroup.c
===================================================================
--- trunk/psModules/src/objects/pmModelGroup.c	(revision 14530)
+++ trunk/psModules/src/objects/pmModelGroup.c	(revision 14652)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-07-20 00:31:43 $
+ *  @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -28,7 +28,4 @@
 #include "pmResiduals.h"
 #include "pmModel.h"
-#include "pmResiduals.h"
-#include "pmPSF.h"
-#include "pmSource.h"
 #include "pmModelGroup.h"
 #include "pmErrorCodes.h"
@@ -198,33 +195,2 @@
     return (models[type].name);
 }
-
-/******************************************************************************
-    pmSourceModelGuess(source, model): This function allocates a new
-    pmModel structure based on the given modelType specified in the argument list.
-    The corresponding pmModelGuess function is returned, and used to
-    supply the values of the params array in the pmModel structure.
- 
-    XXX: Many parameters are based on the src->moments structure, which is in
-    image, not subImage coords.  Therefore, the calls to the model evaluation
-    functions will be in image, not subImage coords.  Remember this.
-*****************************************************************************/
-pmModel *pmSourceModelGuess(pmSource *source,
-                            pmModelType modelType)
-{
-    psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_PTR_NON_NULL(source, NULL);
-    PS_ASSERT_PTR_NON_NULL(source->moments, NULL);
-    PS_ASSERT_PTR_NON_NULL(source->peak, NULL);
-
-    pmModel *model = pmModelAlloc(modelType);
-
-    pmModelGuessFunc modelGuessFunc = pmModelGuessFunc_GetFunction(modelType);
-    if (!modelGuessFunc(model, source)) {
-	psFree (model);
-	return NULL;
-    }
-
-    psTrace("psModules.objects", 3, "---- %s() end ----\n", __func__);
-    return(model);
-}
-
Index: trunk/psModules/src/objects/pmModelGroup.h
===================================================================
--- trunk/psModules/src/objects/pmModelGroup.h	(revision 14530)
+++ trunk/psModules/src/objects/pmModelGroup.h	(revision 14652)
@@ -1,14 +1,26 @@
 /* @file  pmModelGroup.h
  *
- * The object model function types are defined to allow for the flexible addition
- * of new object models. Every object model, with parameters represented by
- * pmModel, has an associated set of functions which provide necessary support
- * operations. A set of abstract functions allow the programmer to select the
- * approriate function or property for a specific named object model.
+ * The object model function types are desined to allow for the flexible addition of new object
+ * models. Every object model, with parameters represented by pmModel, has an associated set of
+ * functions which provide necessary support operations. A set of abstract functions allow the
+ * programmer to select the approriate function or property for a specific named object model.
+ *
+ * Every model instance belongs to a class of models, defined by the value of
+ * the pmModelType type entry. Various functions need access to information about
+ * each of the models. Some of this information varies from model to model, and
+ * may depend on the current parameter values or other data quantities. In order
+ * to keep the code from requiring the information about each model to be coded
+ * into the low-level fitting routines, we define a collection of functions which
+ * allow us to abstract this type of model-dependent information. These generic
+ * functions take the model type and return the corresponding function pointer
+ * for the specified model. Each model is defined by creating this collection of
+ * specific functions, and placing them in a single file for each model. We
+ * define the following structure to carry the collection of information about
+ * the models.
  *
  * @author EAM, IfA
  *
- * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-01-24 02:54:15 $
+ * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-08-24 00:11:02 $
  * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
@@ -23,5 +35,5 @@
 typedef psMinimizeLMChi2Func pmModelFunc;
 
-//  This function is the model chi-square minimization function for this model.
+//  This function sets the parameter limits for this model.
 typedef psMinimizeLMLimitFunc pmModelLimits;
 
@@ -29,41 +41,23 @@
 typedef psF64 (*pmModelFlux)(const psVector *params);
 
-
 // This function returns the radius at which the given model and parameters
 // achieves the given flux.
 typedef psF64 (*pmModelRadius)(const psVector *params, double flux);
 
-/*  This function sets the model parameter limits vectors for the given model
- */
-// typedef bool (*pmModelLimits)(psVector **beta_lim, psVector **params_min, psVector **params_max);
-
-/*  This function provides the model guess parameters based on the details of
- *   the given source.
- */
+//  This function provides the model guess parameters based on the details of
+//  the given source.
 typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source);
 
-
-/*  This function constructs the PSF model for the given source based on the
- *  supplied psf and the EXT model for the object.
- */
+//  This function constructs the PSF model for the given source based on the
+//  supplied psf and the EXT model for the object.
 typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelEXT, pmPSF *psf);
 
-/*  This function returns the success / failure status of the given model fit
- */
+//  This function sets the model parameters based on the PSF for a given coordinate and central
+//  intensity
+typedef bool (*pmModelParamsFromPSF)(pmModel *model, pmPSF *psf, float Xo, float Yo, float Io);
+
+//  This function returns the success / failure status of the given model fit
 typedef bool (*pmModelFitStatusFunc)(pmModel *model);
 
-/* Every model instance belongs to a class of models, defined by the value of
- * the pmModelType type entry. Various functions need access to information about
- * each of the models. Some of this information varies from model to model, and
- * may depend on the current parameter values or other data quantities. In order
- * to keep the code from requiring the information about each model to be coded
- * into the low-level fitting routines, we define a collection of functions which
- * allow us to abstract this type of model-dependent information. These generic
- * functions take the model type and return the corresponding function pointer
- * for the specified model. Each model is defined by creating this collection of
- * specific functions, and placing them in a single file for each model. We
- * define the following structure to carry the collection of information about
- * the models.
- */
 typedef struct
 {
@@ -76,4 +70,5 @@
     pmModelGuessFunc     modelGuessFunc;
     pmModelFromPSFFunc   modelFromPSFFunc;
+    pmModelParamsFromPSF modelParamsFromPSF;
     pmModelFitStatusFunc modelFitStatusFunc;
 }
Index: trunk/psModules/src/objects/pmModelUtils.c
===================================================================
--- trunk/psModules/src/objects/pmModelUtils.c	(revision 14530)
+++ trunk/psModules/src/objects/pmModelUtils.c	(revision 14652)
@@ -5,6 +5,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-08-16 18:33:37 $
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -20,38 +20,24 @@
 #include <string.h>
 #include <pslib.h>
+#include "pmHDU.h"
+#include "pmFPA.h"
 #include "pmResiduals.h"
+#include "pmGrowthCurve.h"
+#include "pmPSF.h"
 #include "pmModel.h"
+#include "pmErrorCodes.h"
+#include "pmModelUtils.h"
 
-// instantiate a model for the PSF at this location (normalized or not?)
-// NOTE: psf and (x,y) are defined wrt chip coordinates
-pmModel *pmModelFromPSFforXY (pmPSF *psf, float x, float y, float flux) {
-
-    assert (psf);
-
-    // allocate a new pmModel to hold the PSF version
-    pmModel *modelEXT = pmModelAlloc (psf->type);
-
-    modelEXT->params->data.F32[PM_PAR_SKY]  = 0.0;
-    modelEXT->params->data.F32[PM_PAR_I0]   = 1.0;
-    modelEXT->params->data.F32[PM_PAR_XPOS] = x;
-    modelEXT->params->data.F32[PM_PAR_YPOS] = y;
-
-    // find function used to set the model parameters
-    pmModelFromPSFFunc modelFromPSFFunc = pmModelFromPSFFunc_GetFunction (psf->type);
-
+/*****************************************************************************
+pmModelFromPSF (*modelEXT, *psf):  use the model position parameters to
+construct a realization of the PSF model at the object coordinates
+ *****************************************************************************/
+pmModel *pmModelFromPSF (pmModel *modelEXT, pmPSF *psf)
+{
     // allocate a new pmModel to hold the PSF version
     pmModel *modelPSF = pmModelAlloc (psf->type);
 
-    // adjust the normalization: 
-    pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type);
-    normFlux = modelFluxFunc (model->params);
-    assert (isfinite(normFlux));
-    assert (normFlux > 0);
-
-    // set the correct normalization
-    modelEXT->params->data.F32[PM_PAR_I0] = flux / normFlux;
-
     // set model parameters for this source based on PSF information
-    if (!modelFromPSFFunc (modelPSF, modelEXT, psf)) {
+    if (!modelEXT->modelFromPSF (modelPSF, modelEXT, psf)) {
         psError(PM_ERR_PSF, false, "Failed to set model params from PSF");
         psFree(modelPSF);
@@ -60,35 +46,44 @@
     // XXX note that model->residuals is just a reference
     modelPSF->residuals = psf->residuals;
-    psFree (modelEXT);
 
     return (modelPSF);
 }
 
-pmSource *pmSourceFromModel (pmModel *model, pmReadout *readout, pmSourceType type) {
+// instantiate a model for the PSF at this location with peak flux
+// NOTE: psf and (Xo,Yo) are defined wrt chip coordinates
+pmModel *pmModelFromPSFforXY (pmPSF *psf, float Xo, float Yo, float Io) {
 
-    pmSource *source = pmSourceAlloc ();
+    assert (psf);
 
-    // use the model centroid for the peak
-    switch (type) {
-      case PM_SOURCE_TYPE_STAR:
-	source->modelPSF = model;
-	break;
-      case PM_SOURCE_TYPE_EXTENDED:
-	source->modelEXT = model;
-	break;
-      default:
-	psAbort ("psphot", "error");
+    // allocate a new pmModel to hold the PSF version
+    pmModel *modelPSF = pmModelAlloc (psf->type);
+
+    // set model parameters for this source based on PSF information
+    if (!modelPSF->modelParamsFromPSF (modelPSF, psf, Xo, Yo, Io)) {
+        psError(PM_ERR_PSF, false, "Failed to set model params from PSF");
+        psFree(modelPSF);
+        return NULL;
     }
 
-    source->peak = pmPeakAlloc ();
+    // XXX note that model->residuals is just a reference
+    modelPSF->residuals = psf->residuals;
 
-    float x = model->params->data.F32[PM_PAR_XPOS];
-    float y = model->params->data.F32[PM_PAR_YPOS];
+    return (modelPSF);
+}
 
-    // XXX need to define the radius in some rational way
-    // XXX x,y are defined wrt readout->image parent, but the model
-    // parameters are defined wrt chip coordinates
-    pmSourceDefinePixels (source, readout, x, y, radius);
+// set this model to have the requested flux
+bool pmModelSetFlux (pmModel *model, float flux) {
 
-    return (source);
+    // set Io to be 1.0
+    model->params->data.F32[PM_PAR_I0] = 1.0;
+
+    // determine the normalized flux
+    float normFlux = model->modelFlux (model->params);
+    assert (isfinite(normFlux));
+    assert (normFlux > 0);
+
+    // set the desired normalization
+    model->params->data.F32[PM_PAR_I0] = flux / normFlux;
+
+    return true;
 }
Index: trunk/psModules/src/objects/pmModelUtils.h
===================================================================
--- trunk/psModules/src/objects/pmModelUtils.h	(revision 14652)
+++ trunk/psModules/src/objects/pmModelUtils.h	(revision 14652)
@@ -0,0 +1,47 @@
+/* @file  pmModelUtils.h
+ *
+ * Utility functions for working with pmSources
+ *
+ * @author EAM, IfA
+ *
+ * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-08-24 00:11:02 $
+ * Copyright 2007 IfA, University of Hawaii
+ */
+
+# ifndef PM_MODEL_UTILS_H
+# define PM_MODEL_UTILS_H
+
+/// @addtogroup Objects Object Detection / Analysis Functions
+/// @{
+
+/**
+ *
+ * This function constructs a pmModel instance based on the pmPSF description
+ * of the PSF. The input is a pmModel with at least the values of the centroid
+ * coordinates (possibly normalization if this is needed) defined. The values of
+ * the PSF-dependent parameters are specified for the specific realization based
+ * on the coordinates of the object.
+ *
+ */
+pmModel *pmModelFromPSF(
+    pmModel *model,                     ///< Add comment
+    pmPSF *psf                          ///< Add comment
+);
+
+pmModel *pmModelFromPSFforXY (
+    pmPSF *psf, 
+    float Xo, 
+    float Yo, 
+    float Io
+    );
+
+bool pmModelSetFlux (
+    pmModel *model, 
+    float flux
+    );
+
+
+
+/// @}
+# endif /* PM_MODEL_UTILS_H */
Index: trunk/psModules/src/objects/pmPSF.c
===================================================================
--- trunk/psModules/src/objects/pmPSF.c	(revision 14530)
+++ trunk/psModules/src/objects/pmPSF.c	(revision 14652)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-06-20 02:22:26 $
+ *  @version $Revision: 1.26 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -28,9 +28,10 @@
 #include "pmMoments.h"
 #include "pmResiduals.h"
+#include "pmGrowthCurve.h"
+#include "pmPSF.h"
 #include "pmModel.h"
 #include "pmSource.h"
-#include "pmGrowthCurve.h"
-#include "pmPSF.h"
-#include "pmModelGroup.h"
+#include "pmModelClass.h"
+#include "pmModelUtils.h"
 #include "pmSourcePhotometry.h"
 #include "pmFPAMaskWeight.h"
@@ -74,5 +75,5 @@
     psFree (psf->ApTrend);
     psFree (psf->growth);
-    psFree (psf->params_NEW);
+    psFree (psf->params);
     psFree (psf->residuals);
     return;
@@ -103,4 +104,6 @@
     psf->skyBias  = 0.0;
     psf->skySat   = 0.0;
+    psf->nPSFstars  = 0;
+    psf->nApResid   = 0;
     psf->poissonErrors = poissonErrors;
 
@@ -124,10 +127,10 @@
     psf->residuals = NULL;
 
-    Nparams = pmModelParameterCount (type);
+    Nparams = pmModelClassParameterCount (type);
     if (!Nparams) {
         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
         return(NULL);
     }
-    psf->params_NEW = psArrayAlloc(Nparams);
+    psf->params = psArrayAlloc(Nparams);
 
     // the order of the PSF parameter (X,Y) fits is determined by the psfTrendMask polynomial
@@ -139,5 +142,5 @@
 
     if (psfTrendMask) {
-        for (int i = 0; i < psf->params_NEW->n; i++) {
+        for (int i = 0; i < psf->params->n; i++) {
             if (i == PM_PAR_SKY)
                 continue;
@@ -155,5 +158,5 @@
                 }
             }
-            psf->params_NEW->data[i] = param;
+            psf->params->data[i] = param;
         }
     }
@@ -161,31 +164,4 @@
     psMemSetDeallocator(psf, (psFreeFunc) pmPSFFree);
     return(psf);
-}
-
-/*****************************************************************************
-pmModelFromPSF (*modelEXT, *psf):  use the model position parameters to
-construct a realization of the PSF model at the object coordinates
- *****************************************************************************/
-pmModel *pmModelFromPSF (pmModel *modelEXT, pmPSF *psf)
-{
-
-    // need to define the relationship between the modelEXT and modelPSF ?
-
-    // find function used to set the model parameters
-    pmModelFromPSFFunc modelFromPSFFunc = pmModelFromPSFFunc_GetFunction (psf->type);
-
-    // allocate a new pmModel to hold the PSF version
-    pmModel *modelPSF = pmModelAlloc (psf->type);
-
-    // set model parameters for this source based on PSF information
-    if (!modelFromPSFFunc (modelPSF, modelEXT, psf)) {
-        psError(PM_ERR_PSF, false, "Failed to set model params from PSF");
-        psFree(modelPSF);
-        return NULL;
-    }
-    // XXX note that model->residuals is just a reference
-    modelPSF->residuals = psf->residuals;
-
-    return (modelPSF);
 }
 
@@ -335,9 +311,9 @@
     va_start(ap, sxy);
 
-    pmModelType type = pmModelSetType (typeName);
+    pmModelType type = pmModelClassGetType (typeName);
     psPolynomial2D *psfTrend = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 0, 0);
     pmPSF *psf = pmPSFAlloc (type, true, psfTrend);
 
-    psVector *par = psVectorAlloc (psf->params_NEW->n, PS_TYPE_F32);
+    psVector *par = psVectorAlloc (psf->params->n, PS_TYPE_F32);
     par->data.F32[PM_PAR_SXX] = sxx;
     par->data.F32[PM_PAR_SYY] = syy;
@@ -348,15 +324,15 @@
     // set the psf shape parameters
     psPolynomial2D *poly = NULL;
-    poly = psf->params_NEW->data[PM_PAR_E0];
+    poly = psf->params->data[PM_PAR_E0];
     poly->coeff[0][0] = pol.e0;
 
-    poly = psf->params_NEW->data[PM_PAR_E1];
+    poly = psf->params->data[PM_PAR_E1];
     poly->coeff[0][0] = pol.e1;
 
-    poly = psf->params_NEW->data[PM_PAR_E2];
+    poly = psf->params->data[PM_PAR_E2];
     poly->coeff[0][0] = pol.e2;
 
-    for (int i = PM_PAR_SXY + 1; i < psf->params_NEW->n; i++) {
-        poly = psf->params_NEW->data[i];
+    for (int i = PM_PAR_SXY + 1; i < psf->params->n; i++) {
+        poly = psf->params->data[i];
         poly->coeff[0][0] = (psF32)va_arg(ap, psF64);
     }
Index: trunk/psModules/src/objects/pmPSF.h
===================================================================
--- trunk/psModules/src/objects/pmPSF.h	(revision 14530)
+++ trunk/psModules/src/objects/pmPSF.h	(revision 14652)
@@ -6,6 +6,6 @@
  * @author EAM, IfA
  *
- * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-06-20 02:22:26 $
+ * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-08-24 00:11:02 $
  * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
@@ -16,4 +16,7 @@
 /// @addtogroup Objects Object Detection / Analysis Functions
 /// @{
+
+// type of model carried by the pmModel structure
+typedef int pmModelType;
 
 typedef enum {
@@ -44,9 +47,6 @@
 typedef struct
 {
-    pmModelType type;   ///< PSF Model in use
-    psArray *params_NEW;   ///< Model parameters (psPolynomial2D)
-    // XXXXX I am changing params: we will allocate elements for the
-    // unfitted elements (So, Io, Xo, Yo) and leave them as NULL
-    // I am using a new name to catch all refs to params with gcc
+    pmModelType type;			///< PSF Model in use
+    psArray *params;			///< Model parameters (psPolynomial2D)
     psPolynomial1D *ChiTrend;         ///< Chisq vs flux fit (correction for systematic errors)
     psPolynomial4D *ApTrend;            ///< ApResid vs (x,y,rflux) (rflux = ten(0.4*mInst))
@@ -79,18 +79,4 @@
 );
 
-/**
- *
- * This function constructs a pmModel instance based on the pmPSF description
- * of the PSF. The input is a pmModel with at least the values of the centroid
- * coordinates (possibly normalization if this is needed) defined. The values of
- * the PSF-dependent parameters are specified for the specific realization based
- * on the coordinates of the object.
- *
- */
-pmModel *pmModelFromPSF(
-    pmModel *model,                     ///< Add comment
-    pmPSF *psf                          ///< Add comment
-);
-
 bool pmPSFMaskApTrend (psPolynomial4D *trend, pmPSFApTrendOptions option);
 
Index: trunk/psModules/src/objects/pmPSF_IO.c
===================================================================
--- trunk/psModules/src/objects/pmPSF_IO.c	(revision 14530)
+++ trunk/psModules/src/objects/pmPSF_IO.c	(revision 14652)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-07-14 03:20:22 $
+ *  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -40,17 +40,17 @@
 #include "pmGrowthCurve.h"
 #include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmPSF_IO.h"
 #include "pmSource.h"
-#include "pmModelGroup.h"
+#include "pmModelClass.h"
 #include "pmSourceIO.h"
 
-bool pmFPAviewCheckDataStatusForPSFmodel (const pmFPAview *view, const pmFPAfile *file)
+bool pmPSFmodelCheckDataStatusForView (const pmFPAview *view, const pmFPAfile *file)
 {
     pmFPA *fpa = file->fpa;
 
     if (view->chip == -1) {
-        bool exists = pmFPACheckDataStatusForPSFmodel (fpa);
+        bool exists = pmPSFmodelCheckDataStatusForFPA (fpa);
         return exists;
     }
@@ -62,5 +62,5 @@
 
     if (view->cell == -1) {
-        bool exists = pmChipCheckDataStatusForPSFmodel (chip);
+        bool exists = pmPSFmodelCheckDataStatusForChip (chip);
         return exists;
     }
@@ -69,60 +69,28 @@
         return false;
     }
-    pmCell *cell = chip->cells->data[view->cell];
-
-    if (view->readout == -1) {
-        bool exists = pmCellCheckDataStatusForPSFmodel (cell);
-        return exists;
-    }
-
-    if (view->readout >= cell->readouts->n) {
-        psError(PS_ERR_IO, true, "Requested readout == %d >= cell->readouds->n == %ld", view->readout, cell->readouts->n);
-        return false;
-    }
-    pmReadout *readout = cell->readouts->data[view->readout];
-
-    return pmReadoutCheckDataStatusForPSFmodel (readout);
-}
-
-bool pmFPACheckDataStatusForPSFmodel (const pmFPA *fpa) {
+    psError(PS_ERR_IO, false, "PSF only valid at the chip level");
+    return false;
+}
+
+bool pmPSFmodelCheckDataStatusForFPA (const pmFPA *fpa) {
 
     for (int i = 0; i < fpa->chips->n; i++) {
         pmChip *chip = fpa->chips->data[i];
         if (!chip) continue;
-        if (pmChipCheckDataStatusForPSFmodel (chip)) return true;
+        if (pmPSFmodelCheckDataStatusForChip (chip)) return true;
     }
     return false;
 }
 
-bool pmChipCheckDataStatusForPSFmodel (const pmChip *chip) {
-
-    for (int i = 0; i < chip->cells->n; i++) {
-        pmCell *cell = chip->cells->data[i];
-        if (!cell) continue;
-        if (pmCellCheckDataStatusForPSFmodel (cell)) return true;
-    }
-    return false;
-}
-
-bool pmCellCheckDataStatusForPSFmodel (const pmCell *cell) {
-
-    for (int i = 0; i < cell->readouts->n; i++) {
-        pmReadout *readout = cell->readouts->data[i];
-        if (!readout) continue;
-        if (pmReadoutCheckDataStatusForPSFmodel (readout)) return true;
-    }
-    return false;
-}
-
-bool pmReadoutCheckDataStatusForPSFmodel (const pmReadout *readout) {
+bool pmPSFmodelCheckDataStatusForChip (const pmChip *chip) {
 
     bool status;
 
     // select the psf of interest
-    pmPSF *psf = psMetadataLookupPtr(&status, readout->analysis, "PSPHOT.PSF");
+    pmPSF *psf = psMetadataLookupPtr(&status, chip->analysis, "PSPHOT.PSF");
     return psf ? true : false;
 }
 
-bool pmFPAviewWritePSFmodel (const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+bool pmPSFmodelWriteForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
 {
 
@@ -130,5 +98,5 @@
 
     if (view->chip == -1) {
-        if (!pmFPAWritePSFmodel(fpa, view, file, config)) {
+        if (!pmPSFmodelWriteFPA(fpa, view, file, config)) {
             psError(PS_ERR_IO, false, "Failed to write PSF for fpa");
             return false;
@@ -143,5 +111,5 @@
 
     if (view->cell == -1) {
-        if (!pmChipWritePSFmodel (chip, view, file, config)) {
+        if (!pmPSFmodelWriteChip (chip, view, file, config)) {
             psError(PS_ERR_IO, false, "Failed to write PSF for chip");
             return false;
@@ -150,36 +118,15 @@
     }
 
-    if (view->cell >= chip->cells->n) {
-        return false;
-    }
-    pmCell *cell = chip->cells->data[view->cell];
-
-    if (view->readout == -1) {
-        if (!pmCellWritePSFmodel (cell, view, file, config)) {
-            psError(PS_ERR_IO, false, "Failed to write PSF for cell");
-            return false;
-        }
-        return true;
-    }
-
-    if (view->readout >= cell->readouts->n) {
-        return false;
-    }
-    pmReadout *readout = cell->readouts->data[view->readout];
-
-    if (!pmReadoutWritePSFmodel (readout, view, file, config)) {
-        psError(PS_ERR_IO, false, "Failed to write PSF for readout");
-        return false;
-    }
-    return true;
+    psError(PS_ERR_IO, false, "PSF must be written at the chip level");
+    return false;
 }
 
 // read in all chip-level PSFmodel files for this FPA
-bool pmFPAWritePSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+bool pmPSFmodelWriteFPA (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
 {
 
     for (int i = 0; i < fpa->chips->n; i++) {
         pmChip *chip = fpa->chips->data[i];
-        if (!pmChipWritePSFmodel (chip, view, file, config)) {
+        if (!pmPSFmodelWriteChip (chip, view, file, config)) {
             psError(PS_ERR_IO, false, "Failed to write PSF for %dth chip", i);
             return false;
@@ -190,39 +137,25 @@
 
 // read in all cell-level PSFmodel files for this chip
-bool pmChipWritePSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-    for (int i = 0; i < chip->cells->n; i++) {
-        pmCell *cell = chip->cells->data[i];
-        if (!pmCellWritePSFmodel (cell, view, file, config)) {
-            psError(PS_ERR_IO, false, "Failed to write PSF for %dth cell", i);
-            return false;
-        }
+bool pmPSFmodelWriteChip (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+{
+    if (!pmPSFmodelWrite (chip->analysis, view, file, config)) {
+	psError(PS_ERR_IO, false, "Failed to write PSF for chip");
+	return false;
     }
     return true;
 }
 
-// read in all readout-level PSFmodel files for this cell
-bool pmCellWritePSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-    for (int i = 0; i < cell->readouts->n; i++) {
-        pmReadout *readout = cell->readouts->data[i];
-        if (!pmReadoutWritePSFmodel (readout, view, file, config)) {
-            psError(PS_ERR_IO, false, "Failed to write PSF for %dth readout", i);
-            return false;
-        }
-    }
-    return true;
-}
-
-// for each Readout (ie, analysed image), we write out
+// for a pmPSF supplied on the analysis metadata, we write out
 // - image header        : FITS Image NAXIS = 0
 // - psf table (+header) : FITS Table
 // - psf resid (+header) : FITS Image
 // if needed, we also write out a PHU blank header
-bool pmReadoutWritePSFmodel (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+bool pmPSFmodelWrite (psMetadata *analysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
 {
     bool status;
     pmHDU *hdu;
     char *headName, *tableName, *residName;
+
+    if (!analysis) return false;
 
     // write a PHU? (only if input image is MEF)
@@ -300,7 +233,7 @@
 
     // select the psf of interest
-    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
+    pmPSF *psf = psMetadataLookupPtr (&status, analysis, "PSPHOT.PSF");
     if (!psf) {
-	psError(PS_ERR_UNKNOWN, true, "missing PSF for this readout");
+	psError(PS_ERR_UNKNOWN, true, "missing PSF for this analysis metadata");
 	psFree (tableName);
 	psFree (residName);
@@ -313,10 +246,10 @@
 	psMetadata *header = psMetadataAlloc();
 
-	char *modelName = pmModelGetType (psf->type);
+	char *modelName = pmModelClassGetName (psf->type);
 	psMetadataAddStr (header, PS_LIST_TAIL, "PSF_NAME", 0, "PSF model name", modelName);
 
 	psMetadataAddBool (header, PS_LIST_TAIL, "POISSON",  0, "Use Poisson errors in fits?", psf->poissonErrors);
 
-	int nPar = pmModelParameterCount (psf->type)    ;
+	int nPar = pmModelClassParameterCount (psf->type)    ;
 	psMetadataAdd (header, PS_LIST_TAIL, "PSF_NPAR", PS_DATA_S32, "PSF model parameter count", nPar);
 
@@ -324,5 +257,5 @@
 	for (int i = 0; i < nPar; i++) {
 	    char name[9];
-	    psPolynomial2D *poly = psf->params_NEW->data[i];
+	    psPolynomial2D *poly = psf->params->data[i];
 	    if (poly == NULL) continue;
 	    snprintf (name, 9, "PAR%02d_NX", i);
@@ -344,5 +277,5 @@
 	psArray *psfTable = psArrayAllocEmpty (100);
 	for (int i = 0; i < nPar; i++) {
-	    psPolynomial2D *poly = psf->params_NEW->data[i];
+	    psPolynomial2D *poly = psf->params->data[i];
 	    if (poly == NULL) continue; // skip unset parameters (eg, XPOS)
 	    for (int ix = 0; ix <= poly->nX; ix++) {
@@ -445,5 +378,5 @@
 
 // if this file needs to have a PHU written out, write one
-bool pmPSF_WritePHU (const pmFPAview *view, pmFPAfile *file, const pmConfig *config) {
+bool pmPSFmodelWritePHU (const pmFPAview *view, pmFPAfile *file, const pmConfig *config) {
 
     // not needed if already written
@@ -487,5 +420,5 @@
 }
 
-bool pmFPAviewReadPSFmodel (const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+bool pmPSFmodelReadForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
 {
 
@@ -493,5 +426,5 @@
 
     if (view->chip == -1) {
-        return pmFPAReadPSFmodel(fpa, view, file, config);
+        return pmPSFmodelReadFPA(fpa, view, file, config);
     }
 
@@ -502,31 +435,18 @@
 
     if (view->cell == -1) {
-        return pmChipReadPSFmodel(chip, view, file, config);
-    }
-
-    if (view->cell >= chip->cells->n) {
-        psAbort("Programming error: view does not apply to FPA.");
-    }
-    pmCell *cell = chip->cells->data[view->cell];
-
-    if (view->readout == -1) {
-        return pmCellReadPSFmodel(cell, view, file, config);
-    }
-
-    if (view->readout >= cell->readouts->n) {
-        psAbort("Programming error: view does not apply to FPA.");
-    }
-    pmReadout *readout = cell->readouts->data[view->readout];
-
-    return pmReadoutReadPSFmodel(readout, view, file, config);
+        return pmPSFmodelReadChip(chip, view, file, config);
+    }
+
+    psError(PS_ERR_IO, false, "PSF must be read at the chip level");
+    return false;
 }
 
 // read in all chip-level PSFmodel files for this FPA
-bool pmFPAReadPSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+bool pmPSFmodelReadFPA (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
 {
     bool success = true;                // Was everything successful?
     for (int i = 0; i < fpa->chips->n; i++) {
         pmChip *chip = fpa->chips->data[i];
-        success &= pmChipReadPSFmodel(chip, view, file, config);
+        success &= pmPSFmodelReadChip(chip, view, file, config);
     }
     return success;
@@ -534,28 +454,16 @@
 
 // read in all cell-level PSFmodel files for this chip
-bool pmChipReadPSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-    bool success = true;                // Was everything successful?
-    for (int i = 0; i < chip->cells->n; i++) {
-        pmCell *cell = chip->cells->data[i];
-        success &= pmCellReadPSFmodel (cell, view, file, config);
-    }
-    return success;
-}
-
-// read in all readout-level PSFmodel files for this cell
-bool pmCellReadPSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-    bool success = true;                // Was everything successful?
-    for (int i = 0; i < cell->readouts->n; i++) {
-        pmReadout *readout = cell->readouts->data[i];
-        success &= pmReadoutReadPSFmodel(readout, view, file, config);
-    }
-    return success;
+bool pmPSFmodelReadChip (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+{
+    if (!pmPSFmodelRead (chip->analysis, view, file, config)) {
+	psError(PS_ERR_IO, false, "Failed to write PSF for chip");
+	return false;
+    }
+    return true;
 }
 
 // for each Readout (ie, analysed image), we write out: header + table with PSF model parameters,
 // and header + image for the PSF residual images
-bool pmReadoutReadPSFmodel(pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+bool pmPSFmodelRead (psMetadata *analysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
 {
     bool status;
@@ -599,5 +507,5 @@
     // load the PSF model parameters from the FITS table
     char *modelName = psMetadataLookupStr (&status, header, "PSF_NAME");
-    pmModelType type = pmModelSetType (modelName);
+    pmModelType type = pmModelClassGetType (modelName);
     if (type == -1) {
         psError(PS_ERR_UNKNOWN, true, "invalid model name %s in psf file %s", modelName, file->filename);
@@ -612,5 +520,5 @@
     // check the number of expected parameters
     int nPar = psMetadataLookupS32 (&status, header, "PSF_NPAR");
-    if (nPar != pmModelParameterCount (psf->type))
+    if (nPar != pmModelClassParameterCount (psf->type))
         psAbort("mismatch model par count");
 
@@ -627,5 +535,5 @@
 	    return false;
 	}
-        psf->params_NEW->data[i] = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXorder, nYorder);
+        psf->params->data[i] = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXorder, nYorder);
     }
 
@@ -650,5 +558,5 @@
         // XXX sanity check here
 
-        psPolynomial2D *poly = psf->params_NEW->data[iPar];
+        psPolynomial2D *poly = psf->params->data[iPar];
 	if (poly == NULL) {
 	    psError(PS_ERR_UNKNOWN, true, "values for parameter %d, but missing NX", iPar);
@@ -698,5 +606,5 @@
     }
 
-    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.PSF",     PS_DATA_UNKNOWN,  "psphot psf", psf);
+    psMetadataAdd (analysis, PS_LIST_TAIL, "PSPHOT.PSF",     PS_DATA_UNKNOWN,  "psphot psf", psf);
     psFree (psf);
 
@@ -708,6 +616,5 @@
 }
 
-/************ old support functions, deprecate? **************/
-
+// create a psMetadata representation (human-readable) of a psf model
 psMetadata *pmPSFtoMetadata (psMetadata *metadata, pmPSF *psf)
 {
@@ -717,12 +624,12 @@
     }
 
-    char *modelName = pmModelGetType (psf->type);
+    char *modelName = pmModelClassGetName (psf->type);
     psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NAME", PS_DATA_STRING, "PSF model name", modelName);
 
-    int nPar = pmModelParameterCount (psf->type)    ;
+    int nPar = pmModelClassParameterCount (psf->type)    ;
     psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NPAR", PS_DATA_S32, "PSF model parameter count", nPar);
 
     for (int i = 0; i < nPar; i++) {
-        psPolynomial2D *poly = psf->params_NEW->data[i];
+        psPolynomial2D *poly = psf->params->data[i];
         if (poly == NULL)
             continue;
@@ -742,4 +649,5 @@
 }
 
+// parse a psMetadata representation (human-readable) of a psf model
 pmPSF *pmPSFfromMetadata (psMetadata *metadata)
 {
@@ -749,5 +657,5 @@
 
     char *modelName = psMetadataLookupPtr (&status, metadata, "PSF_MODEL_NAME");
-    pmModelType type = pmModelSetType (modelName);
+    pmModelType type = pmModelClassGetType (modelName);
 
     bool poissonErrors = psMetadataLookupPtr (&status, metadata, "PSF_POISSON_ERRORS");
@@ -759,5 +667,5 @@
 
     int nPar = psMetadataLookupS32 (&status, metadata, "PSF_MODEL_NPAR");
-    if (nPar != pmModelParameterCount (psf->type))
+    if (nPar != pmModelClassParameterCount (psf->type))
         psAbort("mismatch model par count");
 
@@ -770,6 +678,6 @@
             continue;
         psPolynomial2D *poly = psPolynomial2DfromMetadata (folder);
-        psFree (psf->params_NEW->data[i]);
-        psf->params_NEW->data[i] = poly;
+        psFree (psf->params->data[i]);
+        psf->params->data[i] = poly;
     }
     sprintf (keyword, "APTREND");
@@ -789,64 +697,2 @@
     return (psf);
 }
-
-// read in all readout-level Objects files for this cell
-bool pmReadoutWritePSFmodel_Config (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-    bool status;
-    char *filename;
-    char *realname;
-
-    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
-
-    switch (file->type) {
-      case PM_FPA_FILE_PSF:
-        filename = pmFPAfileNameFromRule (file->filerule, file, view);
-        bool create = file->mode == PM_FPA_MODE_WRITE ? true : false;
-        realname = pmConfigConvertFilename (filename, config, create);
-
-        psMetadata *psfData = pmPSFtoMetadata (NULL, psf);
-        psMetadataConfigWrite (psfData, realname);
-        psFree (psfData);
-        psFree (realname);
-        psFree (filename);
-        return true;
-
-      default:
-        fprintf (stderr, "warning: type mismatch\n");
-        break;
-    }
-    return false;
-}
-
-// read in all readout-level Objects files for this cell
-bool pmReadoutReadPSFmodel_Config (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-
-    unsigned int Nfail;
-    char *filename;
-    char *realname;
-
-    switch (file->type) {
-      case PM_FPA_FILE_PSF:
-        filename = pmFPAfileNameFromRule (file->filerule, file, view);
-        bool create = file->mode == PM_FPA_MODE_WRITE ? true : false;
-        realname = pmConfigConvertFilename (filename, config, create);
-
-        psMetadata *psfData = psMetadataConfigRead(NULL, &Nfail, realname, FALSE);
-        pmPSF *psf = pmPSFfromMetadata (psfData);
-        psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN, "psphot psf", psf);
-
-        psFree (psf);
-        psFree (psfData);
-        psFree (realname);
-        psFree (filename);
-
-        return true;
-
-      default:
-        fprintf (stderr, "warning: type mismatch\n");
-        break;
-    }
-    return false;
-}
-
Index: trunk/psModules/src/objects/pmPSF_IO.h
===================================================================
--- trunk/psModules/src/objects/pmPSF_IO.h	(revision 14530)
+++ trunk/psModules/src/objects/pmPSF_IO.h	(revision 14652)
@@ -6,6 +6,6 @@
  * @author EAM, IfA
  *
- * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-07-14 03:20:22 $
+ * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-08-24 00:11:02 $
  * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
@@ -17,25 +17,22 @@
 /// @{
 
+bool pmPSFmodelWriteForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
+bool pmPSFmodelWriteFPA (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
+bool pmPSFmodelWriteChip (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
+bool pmPSFmodelWrite (psMetadata *analysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
+
+bool pmPSFmodelWritePHU (const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
+
+bool pmPSFmodelReadForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
+bool pmPSFmodelReadFPA (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
+bool pmPSFmodelReadChip (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
+bool pmPSFmodelRead (psMetadata *analysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
+
+bool pmPSFmodelCheckDataStatusForView (const pmFPAview *view, const pmFPAfile *file);
+bool pmPSFmodelCheckDataStatusForFPA (const pmFPA *fpa);
+bool pmPSFmodelCheckDataStatusForChip (const pmChip *chip);
+
 psMetadata *pmPSFtoMetadata (psMetadata *metadata, pmPSF *psf);
 pmPSF *pmPSFfromMetadata (psMetadata *metadata);
-bool pmFPAviewWritePSFmodel (const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
-bool pmFPAWritePSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
-bool pmChipWritePSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
-bool pmCellWritePSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
-bool pmReadoutWritePSFmodel (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
-
-bool pmPSF_WritePHU (const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
-
-bool pmFPAviewReadPSFmodel (const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
-bool pmFPAReadPSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
-bool pmChipReadPSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
-bool pmCellReadPSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
-bool pmReadoutReadPSFmodel (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
-
-bool pmFPAviewCheckDataStatusForPSFmodel (const pmFPAview *view, const pmFPAfile *file);
-bool pmFPACheckDataStatusForPSFmodel (const pmFPA *fpa);
-bool pmChipCheckDataStatusForPSFmodel (const pmChip *chip);
-bool pmCellCheckDataStatusForPSFmodel (const pmCell *cell);
-bool pmReadoutCheckDataStatusForPSFmodel (const pmReadout *readout);
 
 /// @}
Index: trunk/psModules/src/objects/pmPSFtry.c
===================================================================
--- trunk/psModules/src/objects/pmPSFtry.c	(revision 14530)
+++ trunk/psModules/src/objects/pmPSFtry.c	(revision 14652)
@@ -5,6 +5,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.43 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-06-20 02:22:26 $
+ *  @version $Revision: 1.44 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -23,10 +23,12 @@
 #include "pmMoments.h"
 #include "pmResiduals.h"
+#include "pmGrowthCurve.h"
+#include "pmPSF.h"
 #include "pmModel.h"
 #include "pmSource.h"
-#include "pmGrowthCurve.h"
-#include "pmPSF.h"
+#include "pmSourceUtils.h"
 #include "pmPSFtry.h"
-#include "pmModelGroup.h"
+#include "pmModelClass.h"
+#include "pmModelUtils.h"
 #include "pmSourceFitModel.h"
 #include "pmSourcePhotometry.h"
@@ -59,5 +61,5 @@
 
     // validate the requested model name
-    pmModelType type = pmModelSetType (modelName);
+    pmModelType type = pmModelClassGetType (modelName);
     if (type == -1) {
         psError (PS_ERR_UNKNOWN, true, "invalid model name %s", modelName);
@@ -412,9 +414,9 @@
     // This way, the parameters masked by one of the fits will be applied to the others
     for (int i = 0; i < stats->clipIter; i++) {
-        psVectorClipFitPolynomial2D (psf->params_NEW->data[PM_PAR_E0], stats, psfTry->mask, 0xff, e0, dz, x, y);
+        psVectorClipFitPolynomial2D (psf->params->data[PM_PAR_E0], stats, psfTry->mask, 0xff, e0, dz, x, y);
         psTrace ("psModules.pmPSFtry", 4, "clipped E0 : keeping %ld of %ld\n", stats->clippedNvalues, e0->n);
-        psVectorClipFitPolynomial2D (psf->params_NEW->data[PM_PAR_E1], stats, psfTry->mask, 0xff, e1, dz, x, y);
+        psVectorClipFitPolynomial2D (psf->params->data[PM_PAR_E1], stats, psfTry->mask, 0xff, e1, dz, x, y);
         psTrace ("psModules.pmPSFtry", 4, "clipped E1 : keeping %ld of %ld\n", stats->clippedNvalues, e1->n);
-        psVectorClipFitPolynomial2D (psf->params_NEW->data[PM_PAR_E2], stats, psfTry->mask, 0xff, e2, dz, x, y);
+        psVectorClipFitPolynomial2D (psf->params->data[PM_PAR_E2], stats, psfTry->mask, 0xff, e2, dz, x, y);
         psTrace ("psModules.pmPSFtry", 4, "clipped E2 : keeping %ld of %ld\n", stats->clippedNvalues, e2->n);
     }
@@ -436,5 +438,5 @@
 
     // skip the unfitted parameters (X, Y, Io, Sky) and the shape parameters (SXX, SYY, SXY)
-    for (int i = 0; i < psf->params_NEW->n; i++) {
+    for (int i = 0; i < psf->params->n; i++) {
         switch (i) {
           case PM_PAR_SKY:
@@ -461,5 +463,5 @@
         // the mask is carried from previous steps and updated with this operation
         // the weight is either the flux error or NULL, depending on 'applyWeights'
-        if (!psVectorClipFitPolynomial2D(psf->params_NEW->data[i], stats, psfTry->mask, 0xff, z, NULL, x, y)) {
+        if (!psVectorClipFitPolynomial2D(psf->params->data[i], stats, psfTry->mask, 0xff, z, NULL, x, y)) {
             psError(PS_ERR_UNKNOWN, false, "failed to build psf model for parameter %d", i);
             psFree(stats);
@@ -489,6 +491,6 @@
             fprintf (f, "%f %f : ", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
 
-            for (int i = 0; i < psf->params_NEW->n; i++) {
-                if (psf->params_NEW->data[i] == NULL)
+            for (int i = 0; i < psf->params->n; i++) {
+                if (psf->params->data[i] == NULL)
                     continue;
                 fprintf (f, "%f %f : ", model->params->data.F32[i], modelPSF->params->data.F32[i]);
Index: trunk/psModules/src/objects/pmPeaks.c
===================================================================
--- trunk/psModules/src/objects/pmPeaks.c	(revision 14530)
+++ trunk/psModules/src/objects/pmPeaks.c	(revision 14652)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-26 01:20:29 $
+ *  @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -35,5 +35,5 @@
                         pmPeakType type)
 {
-    psTrace(__func__, 5, "---- begin ----\n");
+    psTrace("psModules.objects", 5, "---- begin ----\n");
 
     if (peaks == NULL) {
@@ -67,5 +67,5 @@
     psFree (peak);
 
-    psTrace(__func__, 5, "---- end ----\n");
+    psTrace("psModules.objects", 5, "---- end ----\n");
     return(peaks);
 }
@@ -118,5 +118,5 @@
 *****************************************************************************/
 static void peakFree(pmPeak *tmp)
-{} // used by pmIsPeak()
+{} // used by pmPeakTest()
 
 pmPeak *pmPeakAlloc(psS32 x,
@@ -145,5 +145,5 @@
 }
 
-bool pmIsPeak(const psPtr ptr)
+bool pmPeakTest(const psPtr ptr)
 {
     return (psMemGetDeallocator(ptr) == (psFreeFunc)peakFree);
@@ -193,5 +193,5 @@
 
 /******************************************************************************
-pmFindVectorPeaks(vector, threshold): Find all local peaks in the given vector
+pmPeaksInVector(vector, threshold): Find all local peaks in the given vector
 above the given threshold.  Returns a vector of type PS_TYPE_U32 containing
 the location (x value) of all peaks.
@@ -203,6 +203,6 @@
 Depending upon actual use, this may need to be optimized.
 *****************************************************************************/
-psVector *pmFindVectorPeaks(const psVector *vector,
-                            psF32 threshold)
+psVector *pmPeaksInVector(const psVector *vector,
+			 psF32 threshold)
 {
     psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
@@ -295,5 +295,5 @@
 
 /******************************************************************************
-pmFindImagePeaks(image, threshold): Find all local peaks in the given psImage
+pmPeaksInImage(image, threshold): Find all local peaks in the given psImage
 above the given threshold.  Returns a psArray containing location (x/y value)
 of all peaks.
@@ -309,5 +309,5 @@
 
 *****************************************************************************/
-psArray *pmFindImagePeaks(const psImage *image, psF32 threshold)
+psArray *pmPeaksInImage(const psImage *image, psF32 threshold)
 {
     psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
@@ -327,6 +327,6 @@
     row = 0;
     tmpRow = getRowVectorFromImage((psImage *) image, row);
-    psVector *row1 = pmFindVectorPeaks(tmpRow, threshold);
-    // pmFindVectorPeaks returns coords in the vector, not corrected for col0
+    psVector *row1 = pmPeaksInVector(tmpRow, threshold);
+    // pmPeaksInVector returns coords in the vector, not corrected for col0
 
     for (psU32 i = 0 ; i < row1->n ; i++ ) {
@@ -382,5 +382,5 @@
     for (row = 1 ; row < (image->numRows - 1) ; row++) {
         tmpRow = getRowVectorFromImage((psImage *) image, row);
-        row1 = pmFindVectorPeaks(tmpRow, threshold);
+        row1 = pmPeaksInVector(tmpRow, threshold);
 
         // Step through all local peaks in this row.
@@ -461,5 +461,5 @@
     row = image->numRows - 1;
     tmpRow = getRowVectorFromImage((psImage *) image, row);
-    row1 = pmFindVectorPeaks(tmpRow, threshold);
+    row1 = pmPeaksInVector(tmpRow, threshold);
     for (psU32 i = 0 ; i < row1->n ; i++ ) {
         col = row1->data.U32[i];
@@ -501,48 +501,6 @@
 }
 
-
-/******************************************************************************
-psCullPeaks(peaks, maxValue, valid): eliminate peaks from the psArray that have
-a peak value above the given maximum, or fall outside the valid region.
- 
-XXX: Should the sky value be used when comparing the maximum?
- 
-XXX: warning message if valid is NULL?
- 
-XXX: changed API to create a NEW output psArray (should change name as well)
- 
-XXX: Do we free the psList elements of those culled peaks?
- 
-XXX EAM : do we still need pmCullPeaks, or only pmPeaksSubset?
-*****************************************************************************/
-psList *pmCullPeaks(psList *peaks,
-                    psF32 maxValue,
-                    const psRegion valid)
-{
-    psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_PTR_NON_NULL(peaks, NULL);
-
-    psListElem *tmpListElem = (psListElem *) peaks->head;
-    psS32 indexNum = 0;
-
-    //    printf("pmCullPeaks(): list size is %d\n", peaks->size);
-    while (tmpListElem != NULL) {
-        pmPeak *tmpPeak = (pmPeak *) tmpListElem->data;
-        if ((tmpPeak->value > maxValue) ||
-                (true == isItInThisRegion(valid, tmpPeak->x, tmpPeak->y))) {
-            psListRemoveData(peaks, (psPtr) tmpPeak);
-        }
-
-        indexNum++;
-        tmpListElem = tmpListElem->next;
-    }
-
-    psTrace("psModules.objects", 3, "---- %s() end ----\n", __func__);
-    return(peaks);
-}
-
-// XXX EAM: I changed this to return a new, subset array
-//          rather than alter the existing one
-// XXX: Fix the *valid pointer.
+// return a new array of peaks which are in the valid region and below threshold
+// XXX this function is unused and probably could be dropped
 psArray *pmPeaksSubset(
     psArray *peaks,
@@ -555,5 +513,5 @@
     psArray *output = psArrayAllocEmpty (200);
 
-    psTrace (".pmObjects.pmCullPeaks", 3, "list size is %ld\n", peaks->n);
+    psTrace ("psModules.objects", 3, "list size is %ld\n", peaks->n);
 
     for (int i = 0; i < peaks->n; i++) {
Index: trunk/psModules/src/objects/pmPeaks.h
===================================================================
--- trunk/psModules/src/objects/pmPeaks.h	(revision 14530)
+++ trunk/psModules/src/objects/pmPeaks.h	(revision 14652)
@@ -10,6 +10,6 @@
  * @author GLG, MHPCC
  *
- * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-01-24 02:54:15 $
+ * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-08-24 00:11:02 $
  * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
@@ -72,7 +72,7 @@
 );
 
-bool pmIsPeak(const psPtr ptr);
+bool pmPeakTest(const psPtr ptr);
 
-/** pmFindVectorPeaks()
+/** pmPeaksInVector()
  *
  * Find all local peaks in the given vector above the given threshold. A peak
@@ -91,5 +91,5 @@
  *
  */
-psVector *pmFindVectorPeaks(
+psVector *pmPeaksInVector(
     const psVector *vector,  ///< The input vector (float)
     float threshold   ///< Threshold above which to find a peak
@@ -97,5 +97,5 @@
 
 
-/** pmFindImagePeaks()
+/** pmPeaksInImage()
  *
  * Find all local peaks in the given image above the given threshold. This
@@ -111,5 +111,5 @@
  *
  */
-psArray *pmFindImagePeaks(
+psArray *pmPeaksInImage(
     const psImage *image,  ///< The input image where peaks will be found (float)
     float threshold   ///< Threshold above which to find a peak
Index: trunk/psModules/src/objects/pmSource.c
===================================================================
--- trunk/psModules/src/objects/pmSource.c	(revision 14530)
+++ trunk/psModules/src/objects/pmSource.c	(revision 14652)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-08-16 18:33:00 $
+ *  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -27,4 +27,6 @@
 #include "pmMoments.h"
 #include "pmResiduals.h"
+#include "pmGrowthCurve.h"
+#include "pmPSF.h"
 #include "pmModel.h"
 #include "pmSource.h"
@@ -113,5 +115,5 @@
 }
 
-bool pmIsSource(const psPtr ptr)
+bool pmSourceTest(const psPtr ptr)
 {
     return (psMemGetDeallocator(ptr) == (psFreeFunc)sourceFree);
@@ -337,5 +339,5 @@
             return emptyClump;
         }
-        peaks = pmFindImagePeaks (splane, stats[0].max / 2);
+        peaks = pmPeaksInImage (splane, stats[0].max / 2);
         psTrace ("psModules.objects", 2, "clump threshold is %f\n", stats[0].max/2);
 
@@ -797,5 +799,5 @@
 
 // should we call pmSourceCacheModel if it does not exist?
-bool pmSourceOp (pmSource *source, pmModelOpMode mode, bool add, psMaskType maskVal) {
+bool pmSourceOp (pmSource *source, pmModelOpMode mode, bool add, psMaskType maskVal, int dx, int dy) {
 
     bool status;
@@ -856,8 +858,9 @@
         target = source->weight;
     }
+
     if (add) {
-        status = pmModelAdd (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal);
+        status = pmModelAddWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
     } else {
-        status = pmModelSub (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal);
+        status = pmModelSubWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
     }
 
@@ -866,9 +869,17 @@
 
 bool pmSourceAdd (pmSource *source, pmModelOpMode mode, psMaskType maskVal) {
-    return pmSourceOp (source, mode, true, maskVal);
+    return pmSourceOp (source, mode, true, maskVal, 0, 0);
 }
 
 bool pmSourceSub (pmSource *source, pmModelOpMode mode, psMaskType maskVal) {
-    return pmSourceOp (source, mode, false, maskVal);
+    return pmSourceOp (source, mode, false, maskVal, 0, 0);
+}
+
+bool pmSourceAddWithOffset (pmSource *source, pmModelOpMode mode, psMaskType maskVal, int dx, int dy) {
+    return pmSourceOp (source, mode, true, maskVal, dx, dy);
+}
+
+bool pmSourceSubWithOffset (pmSource *source, pmModelOpMode mode, psMaskType maskVal, int dx, int dy) {
+    return pmSourceOp (source, mode, false, maskVal, dx, dy);
 }
 
Index: trunk/psModules/src/objects/pmSource.h
===================================================================
--- trunk/psModules/src/objects/pmSource.h	(revision 14530)
+++ trunk/psModules/src/objects/pmSource.h	(revision 14652)
@@ -3,6 +3,6 @@
  * @author EAM, IfA; GLG, MHPCC
  *
- * @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-08-15 20:21:18 $
+ * @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-08-24 00:11:02 $
  * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
@@ -57,6 +57,5 @@
  *
  */
-typedef struct
-{
+struct pmSource {
     const int id;                       ///< Unique ID for object
     pmPeak *peak;                       ///< Description of peak pixel.
@@ -81,6 +80,5 @@
     psRegion region;                    ///< area on image covered by selected pixels
     float sky, skyErr;                  ///< The sky and its error at the center of the object
-}
-pmSource;
+};
 
 /** pmPSFClump data structure
@@ -114,5 +112,5 @@
 void pmSourceFreePixels(pmSource *source);
 
-bool pmIsSource(const psPtr ptr);
+bool pmSourceTest(const psPtr ptr);
 
 /** pmSourceDefinePixels()
@@ -218,6 +216,8 @@
 bool pmSourceAdd (pmSource *source, pmModelOpMode mode, psMaskType maskVal);
 bool pmSourceSub (pmSource *source, pmModelOpMode mode, psMaskType maskVal);
-
-bool pmSourceOp (pmSource *source, pmModelOpMode mode, bool add, psMaskType maskVal);
+bool pmSourceAddWithOffset (pmSource *source, pmModelOpMode mode, psMaskType maskVal, int dx, int dy);
+bool pmSourceSubWithOffset (pmSource *source, pmModelOpMode mode, psMaskType maskVal, int dx, int dy);
+
+bool pmSourceOp (pmSource *source, pmModelOpMode mode, bool add, psMaskType maskVal, int dx, int dy);
 bool pmSourceCacheModel (pmSource *source, psMaskType maskVal);
 bool pmSourceCachePSF (pmSource *source, psMaskType maskVal);
Index: trunk/psModules/src/objects/pmSourceContour.c
===================================================================
--- trunk/psModules/src/objects/pmSourceContour.c	(revision 14530)
+++ trunk/psModules/src/objects/pmSourceContour.c	(revision 14652)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -26,4 +26,6 @@
 #include "pmMoments.h"
 #include "pmResiduals.h"
+#include "pmGrowthCurve.h"
+#include "pmPSF.h"
 #include "pmModel.h"
 #include "pmSource.h"
Index: trunk/psModules/src/objects/pmSourceFitModel.c
===================================================================
--- trunk/psModules/src/objects/pmSourceFitModel.c	(revision 14530)
+++ trunk/psModules/src/objects/pmSourceFitModel.c	(revision 14652)
@@ -6,6 +6,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-06-21 22:58:11 $
+ *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -27,8 +27,8 @@
 #include "pmGrowthCurve.h"
 #include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmSource.h"
-#include "pmModelGroup.h"
+#include "pmModelClass.h"
 #include "pmSourceFitModel.h"
 
@@ -116,15 +116,8 @@
     psVector *dparams = model->dparams;
 
-    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
-    if (!modelFunc)
-        psAbort("invalid model function");
-    pmModelLimits checkLimits = pmModelLimits_GetFunction (model->type);
-    if (!checkLimits)
-        psAbort("invalid model limits function");
-
     // create the minimization constraints
     psMinConstraint *constraint = psMinConstraintAlloc();
     constraint->paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
-    constraint->checkLimits = checkLimits;
+    constraint->checkLimits = model->modelLimits;
 
     // set parameter mask based on fitting mode
@@ -156,6 +149,6 @@
     // force the floating parameters to fall within the contraint ranges
     for (int i = 0; i < params->n; i++) {
-        checkLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL);
-        checkLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL);
+        model->modelLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL);
+        model->modelLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL);
     }
 
@@ -174,5 +167,5 @@
     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
 
-    fitStatus = psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, modelFunc);
+    fitStatus = psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, model->modelFunc);
     for (int i = 0; i < dparams->n; i++) {
         if (psTraceGetLevel("psModules.objects") >= 4) {
@@ -201,5 +194,5 @@
             altmask->data.U8[i] = (constraint->paramMask->data.U8[i]) ? 0 : 1;
         }
-        psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, modelFunc);
+        psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, model->modelFunc);
         for (int i = 0; i < dparams->n; i++) {
             if (!constraint->paramMask->data.U8[i])
@@ -237,332 +230,2 @@
 }
 
-/********************* Source Model Set Functions ***************************/
-
-// these static variables are used by one pass of pmSourceFitSet to store
-// data for a model set.  If we are going to make psphot thread-safe, these
-// will have to go in a structure of their own or be allocated once per thread
-// sky, p1.1, p1.2, p1.3,... p1.n, p2.1, p2.2,
-// nPar = nSrc*(nOnePar - 1) + 1
-static pmModelFunc oneModelFunc;
-static pmModelLimits oneCheckLimits;
-static psVector *onePar;
-static psVector *oneDeriv;
-static int nOnePar;
-
-bool pmSourceFitSetInit (pmModelType type)
-{
-
-    oneModelFunc = pmModelFunc_GetFunction (type);
-    oneCheckLimits = pmModelLimits_GetFunction (type);
-    nOnePar = pmModelParameterCount (type);
-
-    onePar = psVectorAlloc (nOnePar, PS_DATA_F32);
-    oneDeriv = psVectorAlloc (nOnePar, PS_DATA_F32);
-
-    return true;
-}
-
-void pmSourceFitSetClear (void)
-{
-
-    psFree (onePar);
-    psFree (oneDeriv);
-    return;
-}
-
-bool pmSourceFitSet_CheckLimits (psMinConstraintMode mode, int nParam, float *params, float *betas)
-{
-    // convert the value of nParam into corresponding single model parameter entry
-    // convert params into corresponding single model parameter pointer
-
-    int nParamSingle = (nParam - 1) % (nOnePar - 1) + 1;
-    float *paramSingle = params + nParam - nParamSingle;
-    float *betaSingle = betas + nParam - nParamSingle;
-    bool status = oneCheckLimits (mode, nParamSingle, paramSingle, betaSingle);
-    return status;
-}
-
-psF32 pmSourceFitSet_Function(psVector *deriv,
-                              const psVector *params,
-                              const psVector *x)
-{
-
-    psF32 value;
-    psF32 model;
-
-    psF32 *PAR = onePar->data.F32;
-    psF32 *dPAR = oneDeriv->data.F32;
-
-    psF32 *pars = params->data.F32;
-    psF32 *dpars = (deriv == NULL) ? NULL : deriv->data.F32;
-
-    int nSrc = (params->n - 1) / (nOnePar - 1);
-
-    PAR[0] = model = pars[0];
-    for (int i = 0; i < nSrc; i++) {
-        int nOff = i*nOnePar - i;
-        for (int n = 1; n < nOnePar; n++) {
-            PAR[n] = pars[nOff + n];
-        }
-        if (deriv == NULL) {
-            value = oneModelFunc (NULL, onePar, x);
-        } else {
-            value = oneModelFunc (oneDeriv, onePar, x);
-            for (int n = 1; n < nOnePar; n++) {
-                dpars[nOff + n] = dPAR[n];
-            }
-        }
-        model += value;
-    }
-    if (deriv != NULL) {
-        dpars[0] = dPAR[0]*2.0;
-    }
-    return (model);
-}
-
-/*
-i:       0                   1                 2
-n:         1  2  3  4  5  6  1  2  3  4  5  6  1  2  3  4  5  6
-i*6 + n: 0 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18
-*/
-
-bool pmSourceFitSet (pmSource *source,
-                     psArray *modelSet,
-                     pmSourceFitMode mode,
-                     psMaskType maskVal)
-{
-    psTrace("psModules.objects", 3, "---- %s begin ----\n", __func__);
-    PS_ASSERT_PTR_NON_NULL(source, false);
-    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
-    PS_ASSERT_PTR_NON_NULL(source->maskObj, false);
-    PS_ASSERT_PTR_NON_NULL(source->weight, false);
-
-    psBool fitStatus = true;
-    psBool onPic     = true;
-    psBool rc        = true;
-
-    // maximum number of valid pixels
-    psS32 nPix = source->pixels->numRows * source->pixels->numCols;
-
-    // construct the coordinate and value entries
-    psArray *x = psArrayAllocEmpty(nPix);
-    psVector *y = psVectorAllocEmpty(nPix, PS_TYPE_F32);
-    psVector *yErr = psVectorAllocEmpty(nPix, PS_TYPE_F32);
-
-    // fill in the coordinate and value entries
-    nPix = 0;
-    for (psS32 i = 0; i < source->pixels->numRows; i++) {
-        for (psS32 j = 0; j < source->pixels->numCols; j++) {
-            // skip masked points
-            if (source->maskObj->data.U8[i][j] & maskVal) {
-                continue;
-            }
-            // skip zero-weight points
-            if (source->weight->data.F32[i][j] == 0) {
-                continue;
-            }
-            // skip nan values in image
-            if (!isfinite(source->pixels->data.F32[i][j])) {
-                continue;
-            }
-
-            psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
-
-            // Convert i/j to image space:
-            coord->data.F32[0] = (psF32) (j + source->pixels->col0);
-            coord->data.F32[1] = (psF32) (i + source->pixels->row0);
-            x->data[nPix] = (psPtr *) coord;
-            y->data.F32[nPix] = source->pixels->data.F32[i][j];
-
-            // psMinimizeLMChi2 takes wt = 1/dY^2.  suggestion from RHL is to use the local sky
-            // as weight to avoid the bias from systematic errors here we would just use the
-            // source sky variance
-            if (PM_SOURCE_FIT_MODEL_PIX_WEIGHTS) {
-                yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j];
-            } else {
-                yErr->data.F32[nPix] = 1.0 / PM_SOURCE_FIT_MODEL_WEIGHT;
-            }
-            nPix++;
-        }
-    }
-    x->n = nPix;
-    y->n = nPix;
-    yErr->n = nPix;
-
-    // base values on first model (all models must be identical)
-    pmModel *model = modelSet->data[0];
-
-    // determine number of model parameters
-    int nSrc = modelSet->n;
-    int nPar = model->params->n - 1;  // number of object parameters (excluding sky)
-
-    // define parameter vectors for source set
-    psVector *params = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
-    psVector *dparams = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
-
-    // set the static variables
-    pmSourceFitSetInit (model->type);
-
-    // create the minimization constraints
-    psMinConstraint *constraint = psMinConstraintAlloc();
-    constraint->paramMask = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_U8);
-    constraint->checkLimits = pmSourceFitSet_CheckLimits;
-
-    // set the parameter guesses for the multiple models
-    // first the value for the single sky parameter
-    params->data.F32[0] = model->params->data.F32[0];
-
-    // next, the values for the source parameters
-    for (int i = 0; i < nSrc; i++) {
-        model = modelSet->data[i];
-        for (int n = 1; n < nPar + 1; n++) {
-            params->data.F32[i*nPar + n] = model->params->data.F32[n];
-        }
-    }
-
-    if (psTraceGetLevel("psModules.objects") >= 5) {
-        for (int i = 0; i < params->n; i++) {
-            fprintf (stderr, "%d %f %d\n", i, params->data.F32[i], constraint->paramMask->data.U8[i]);
-        }
-    }
-
-    // set the parameter masks based on the fitting mode
-    int nParams = 0;
-    switch (mode) {
-    case PM_SOURCE_FIT_NORM:
-        // NORM-only model fits only source normalization (Io)
-        nParams = nSrc;
-        psVectorInit (constraint->paramMask, 1);
-        for (int i = 0; i < nSrc; i++) {
-            constraint->paramMask->data.U8[1 + i*nPar] = 0;
-        }
-        break;
-    case PM_SOURCE_FIT_PSF:
-        // PSF model only fits x,y,Io
-        nParams = 3*nSrc;
-        psVectorInit (constraint->paramMask, 1);
-        for (int i = 0; i < nSrc; i++) {
-            constraint->paramMask->data.U8[1 + i*nPar] = 0;
-            constraint->paramMask->data.U8[2 + i*nPar] = 0;
-            constraint->paramMask->data.U8[3 + i*nPar] = 0;
-        }
-        break;
-    case PM_SOURCE_FIT_EXT:
-        // EXT model fits all params (except sky)
-        nParams = nPar*nSrc;
-        psVectorInit (constraint->paramMask, 0);
-        constraint->paramMask->data.U8[0] = 1;
-        break;
-    default:
-        psAbort("invalid fitting mode");
-    }
-
-    // force the floating parameters to fall within the contraint ranges
-    for (int i = 0; i < params->n; i++) {
-        pmSourceFitSet_CheckLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL);
-        pmSourceFitSet_CheckLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL);
-    }
-
-    if (nPix <  nParams + 1) {
-        psTrace (__func__, 4, "insufficient valid pixels\n");
-        psTrace("psModules.objects", 3, "---- %s() end : fail pixels ----\n", __func__);
-        model->flags |= PM_MODEL_STATUS_BADARGS;
-        psFree (x);
-        psFree (y);
-        psFree (yErr);
-        psFree (params);
-        psFree (dparams);
-        psFree(constraint);
-        pmSourceFitSetClear ();
-        return(false);
-    }
-
-    psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_TOLERANCE);
-
-    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
-
-    fitStatus = psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, pmSourceFitSet_Function);
-    if (!fitStatus) {
-        psTrace("psModules.objects", 4, "Failed to fit model (%d)\n", nSrc);
-    }
-
-    // parameter errors from the covariance matrix
-    for (int i = 0; i < dparams->n; i++) {
-        if ((constraint->paramMask != NULL) && constraint->paramMask->data.U8[i])
-            continue;
-        dparams->data.F32[i] = sqrt(covar->data.F32[i][i]);
-    }
-
-    // get the Gauss-Newton distance for fixed model parameters
-    if (constraint->paramMask != NULL) {
-        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F32);
-        psVector *altmask = psVectorAlloc (params->n, PS_TYPE_U8);
-        altmask->data.U8[0] = 1;
-        for (int i = 1; i < dparams->n; i++) {
-            altmask->data.U8[i] = (constraint->paramMask->data.U8[i]) ? 0 : 1;
-        }
-        psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, pmSourceFitSet_Function);
-        for (int i = 0; i < dparams->n; i++) {
-            if (!constraint->paramMask->data.U8[i])
-                continue;
-            // note that delta is the value *subtracted* from the parameter
-            // to get the new guess.  for dparams to represent the direction
-            // of motion, we need to take -delta
-            dparams->data.F32[i] = -delta->data.F32[i];
-        }
-        psFree (delta);
-        psFree (altmask);
-    }
-
-    // assign back the parameters to the models
-    for (int i = 0; i < nSrc; i++) {
-        model = modelSet->data[i];
-        model->params->data.F32[0] = params->data.F32[0];
-        for (int n = 1; n < nPar + 1; n++) {
-            if (psTraceGetLevel("psModules.objects") >= 4) {
-                fprintf (stderr, "%f ", params->data.F32[i*nPar + n]);
-            }
-            model->params->data.F32[n] = params->data.F32[i*nPar + n];
-            model->dparams->data.F32[n] = dparams->data.F32[i*nPar + n];
-        }
-        psTrace ("psModules.objects", 4, " src %d", i);
-
-        // save the resulting chisq, nDOF, nIter
-        // these are not unique for any one source
-        model->chisq = myMin->value;
-        model->nIter = myMin->iter;
-        model->nDOF  = y->n - nParams;
-
-        // set the model success or failure status
-        model->flags |= PM_MODEL_STATUS_FITTED;
-        if (!fitStatus) model->flags |= PM_MODEL_STATUS_NONCONVERGE;
-
-        // models can go insane: reject these
-        onPic &= (model->params->data.F32[PM_PAR_XPOS] >= source->pixels->col0);
-        onPic &= (model->params->data.F32[PM_PAR_XPOS] <  source->pixels->col0 + source->pixels->numCols);
-        onPic &= (model->params->data.F32[PM_PAR_XPOS] >= source->pixels->row0);
-        onPic &= (model->params->data.F32[PM_PAR_XPOS] <  source->pixels->row0 + source->pixels->numRows);
-        if (!onPic) model->flags |= PM_MODEL_STATUS_OFFIMAGE;
-    }
-    psTrace ("psModules.objects", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value);
-
-    source->mode |= PM_SOURCE_MODE_FITTED;
-
-    psFree(x);
-    psFree(y);
-    psFree(yErr);
-    psFree(myMin);
-    psFree(covar);
-    psFree(constraint);
-    psFree(params);
-    psFree(dparams);
-
-    // free static memory used by pmSourceFitSet
-    pmSourceFitSetClear ();
-
-    rc = (onPic && fitStatus);
-    psTrace (__func__, 5, "onPic: %d, fitStatus: %d, nIter: %d, chisq: %f, nDof: %d\n", onPic, fitStatus, model->nIter, model->chisq, model->nDOF);
-    psTrace("psModules.objects", 5, "---- %s end (%d) ----\n", __func__, rc);
-    return(rc);
-}
-
Index: trunk/psModules/src/objects/pmSourceFitSet.c
===================================================================
--- trunk/psModules/src/objects/pmSourceFitSet.c	(revision 14652)
+++ trunk/psModules/src/objects/pmSourceFitSet.c	(revision 14652)
@@ -0,0 +1,417 @@
+/** @file  pmSourceFitModel.c
+ *
+ *  fit single source models to image pixels
+ *
+ *  @author EAM, IfA
+ *  @author GLG, MHPCC
+ *
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
+ *
+ *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <pslib.h>
+#include "pmHDU.h"
+#include "pmFPA.h"
+#include "pmPeaks.h"
+#include "pmMoments.h"
+#include "pmGrowthCurve.h"
+#include "pmResiduals.h"
+#include "pmPSF.h"
+#include "pmModel.h"
+#include "pmSource.h"
+#include "pmModelClass.h"
+#include "pmSourceFitModel.h"
+#include "pmSourceFitSet.h"
+
+// save as static values so they may be set externally
+static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15;
+static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1;
+static psF32 PM_SOURCE_FIT_MODEL_WEIGHT = 1.0;
+static bool  PM_SOURCE_FIT_MODEL_PIX_WEIGHTS = true;
+
+/********************* Source Model Set Functions ***************************/
+
+static pmSourceFitSetData *thisSet = NULL;
+
+static void pmSourceFitSetDataFree (pmSourceFitSetData *set) {
+    if (!set) return;
+
+    psFree (set->modelSet);
+    psFree (set->paramSet);
+    psFree (set->derivSet);
+    return;
+}
+
+pmSourceFitSetData *pmSourceFitSetDataAlloc (psArray *modelSet) {
+
+    pmSourceFitSetData *set = (pmSourceFitSetData *) psAlloc(sizeof(pmSourceFitSetData));
+    psMemSetDeallocator(set, (psFreeFunc) pmSourceFitSetDataFree);
+
+    set->modelSet = psMemIncrRefCounter (modelSet);
+    set->paramSet = psArrayAlloc (modelSet->n);
+    set->derivSet = psArrayAlloc (modelSet->n);
+    set->nParamSet = 0;
+
+    for (int i = 0; i < modelSet->n; i++) {
+	pmModel *model = modelSet->data[i];
+	
+	int nParams = pmModelClassParameterCount (model->type);
+
+	set->paramSet->data[i] = psVectorAlloc (nParams, PS_DATA_F32);
+	set->derivSet->data[i] = psVectorAlloc (nParams, PS_DATA_F32);
+
+	set->nParamSet += nParams;
+    }
+
+    return set;
+}
+
+// this function is called with the full set of parameters and the beta values in a single vector
+bool pmSourceFitSetCheckLimits (psMinConstraintMode mode, int nParam, float *params, float *betas)
+{
+    assert (thisSet);
+
+    // nParam is the parameter in the full sequence.  determine which single model this comes from
+    int nModel = -1;
+    int nParamOne = -1;
+    int nParamBase = 0;
+    for (int i = 0; i < thisSet->modelSet->n; i++) {
+	psVector *param = thisSet->paramSet->data[i];
+	if ((nParamBase <= nParam) && (nParam < nParamBase + param->n)) {
+	    nModel = i;
+	    nParamOne = nParam - nParamBase;
+	    break;
+	}
+	nParamBase += param->n;
+    }
+    assert (nModel > -1);
+
+    pmModel *model = thisSet->modelSet->data[nModel];
+
+    // pass the single model function a pointer to the start of that model's sequence
+    float *paramOne = params + nParamBase;
+    float *betaOne = betas + nParamBase;
+    bool status = model->modelLimits (mode, nParamOne, paramOne, betaOne);
+    return status;
+}
+
+bool pmSourceFitSetJoin (psVector *deriv, psVector *param, pmSourceFitSetData *set) {
+
+    assert (set);
+    assert (set->paramSet->n == set->derivSet->n);
+
+    int n = 0;
+    for (int i = 0; i < set->paramSet->n; i++) {
+	
+	psVector *paramOne = set->paramSet->data[i];
+	psVector *derivOne = set->derivSet->data[i];
+	assert ((deriv == NULL) || (paramOne->n == derivOne->n));
+
+	for (int j = 0; j < paramOne->n; j++, n++) {
+	    param->data.F32[n] = paramOne->data.F32[j];
+	    if (deriv) {
+		deriv->data.F32[n] = derivOne->data.F32[j];
+	    }
+	}
+    }
+    return true;
+}
+
+bool pmSourceFitSetSplit (pmSourceFitSetData *set, const psVector *deriv, const psVector *param) {
+
+    assert (set);
+    assert (set->paramSet->n == set->derivSet->n);
+
+    int n = 0;
+    for (int i = 0; i < set->paramSet->n; i++) {
+	
+	psVector *paramOne = set->paramSet->data[i];
+	psVector *derivOne = set->derivSet->data[i];
+	assert ((deriv == NULL) || (paramOne->n == derivOne->n));
+
+	for (int j = 0; j < paramOne->n; j++, n++) {
+	    paramOne->data.F32[j] = param->data.F32[n];
+	    if (deriv) {
+		derivOne->data.F32[j] = deriv->data.F32[n];
+	    }
+	}
+    }
+    return true;
+}
+
+bool pmSourceFitSetValues (pmSourceFitSetData *set, const psVector *dparam, const psVector *param, pmSource *source, psMinimization *myMin, int nPix, bool fitStatus) {
+
+    bool onPic = true;
+
+    assert (set);
+
+    int n = 0;
+    for (int i = 0; i < set->paramSet->n; i++) {
+	
+	pmModel *model = set->modelSet->data[i];
+
+	for (int j = 0; j < model->params->n; j++, n++) {
+            if (psTraceGetLevel("psModules.objects") >= 4) {
+                fprintf (stderr, "%f ", param->data.F32[n]);
+            }
+	    model->params->data.F32[j] = param->data.F32[n];
+	    model->dparams->data.F32[j] = dparam->data.F32[n];
+	}
+        psTrace ("psModules.objects", 4, " src %d", i);
+
+        // save the resulting chisq, nDOF, nIter
+        // these are not unique for any one source
+        model->chisq = myMin->value;
+        model->nIter = myMin->iter;
+        model->nDOF  = nPix - model->params->n;
+
+        // set the model success or failure status
+        model->flags |= PM_MODEL_STATUS_FITTED;
+        if (!fitStatus) model->flags |= PM_MODEL_STATUS_NONCONVERGE;
+
+        // models can go insane: reject these
+        onPic &= (model->params->data.F32[PM_PAR_XPOS] >= source->pixels->col0);
+        onPic &= (model->params->data.F32[PM_PAR_XPOS] <  source->pixels->col0 + source->pixels->numCols);
+        onPic &= (model->params->data.F32[PM_PAR_XPOS] >= source->pixels->row0);
+        onPic &= (model->params->data.F32[PM_PAR_XPOS] <  source->pixels->row0 + source->pixels->numRows);
+        if (!onPic) model->flags |= PM_MODEL_STATUS_OFFIMAGE;
+    }
+    return true;
+}
+
+psF32 pmSourceFitSetFunction(psVector *deriv, const psVector *param, const psVector *x)
+{
+    float chisqSum = 0.0;
+    float chisqOne = 0.0;
+
+    assert (thisSet);
+    pmSourceFitSetSplit (thisSet, deriv, param);
+
+    for (int i = 0; i < thisSet->modelSet->n; i++) {
+	
+	pmModel *model = thisSet->modelSet->data[i];
+
+	psVector *paramOne = thisSet->paramSet->data[i];
+	psVector *derivOne = thisSet->derivSet->data[i];
+
+	chisqOne = model->modelFunc (derivOne, paramOne, x);
+	chisqSum += chisqOne;
+    }
+
+    return (chisqSum);
+}
+
+// XXX allow the mode to be a function of the object (eg, S/N)
+bool pmSourceFitSetMasks (psMinConstraint *constraint, pmSourceFitSetData *set, pmSourceFitMode mode) {
+
+    // unmask everyone
+    psVectorInit (constraint->paramMask, 0);
+
+    int n = 0;
+    for (int i = 0; i < set->paramSet->n; i++) {
+	psVector *paramOne = set->paramSet->data[i];
+
+	switch (mode) {
+	  case PM_SOURCE_FIT_NORM:
+	    // mask all but Xo,Yo,Io
+	    for (int j = 0; j < paramOne->n; j++) {
+		if (j == PM_PAR_I0) continue;
+		constraint->paramMask->data.U8[n + j] = 1;
+	    }
+	    break;
+	  case PM_SOURCE_FIT_PSF:
+	    // mask all but Xo,Yo,Io
+	    for (int j = 0; j < paramOne->n; j++) {
+		if (j == PM_PAR_XPOS) continue;
+		if (j == PM_PAR_YPOS) continue;
+		if (j == PM_PAR_I0) continue;
+		constraint->paramMask->data.U8[n + j] = 1;
+	    }
+	    break;
+	  case PM_SOURCE_FIT_EXT:
+	    // EXT model fits all params (except sky)
+	    constraint->paramMask->data.U8[n + PM_PAR_SKY] = 1;
+	    break;
+	  default:
+	    psAbort("invalid fitting mode");
+	}
+	n += paramOne->n;
+    }
+    return true;
+}
+
+bool pmSourceFitSet (pmSource *source,
+                     psArray *modelSet,
+                     pmSourceFitMode mode,
+                     psMaskType maskVal)
+{
+    psTrace("psModules.objects", 3, "---- %s begin ----\n", __func__);
+    PS_ASSERT_PTR_NON_NULL(source, false);
+    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
+    PS_ASSERT_PTR_NON_NULL(source->maskObj, false);
+    PS_ASSERT_PTR_NON_NULL(source->weight, false);
+
+    bool fitStatus = true;
+    bool onPic     = true;
+
+    // maximum number of valid pixels
+    int nPix = source->pixels->numRows * source->pixels->numCols;
+
+    // construct the coordinate and value entries
+    psArray *x = psArrayAllocEmpty(nPix);
+    psVector *y = psVectorAllocEmpty(nPix, PS_TYPE_F32);
+    psVector *yErr = psVectorAllocEmpty(nPix, PS_TYPE_F32);
+
+    // fill in the coordinate and value entries
+    nPix = 0;
+    for (psS32 i = 0; i < source->pixels->numRows; i++) {
+        for (psS32 j = 0; j < source->pixels->numCols; j++) {
+            // skip masked points
+            if (source->maskObj->data.U8[i][j] & maskVal) {
+                continue;
+            }
+            // skip zero-weight points
+            if (source->weight->data.F32[i][j] == 0) {
+                continue;
+            }
+            // skip nan values in image
+            if (!isfinite(source->pixels->data.F32[i][j])) {
+                continue;
+            }
+
+            psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
+
+            // Convert i/j to image space:
+            coord->data.F32[0] = (psF32) (j + source->pixels->col0);
+            coord->data.F32[1] = (psF32) (i + source->pixels->row0);
+            x->data[nPix] = (psPtr *) coord;
+            y->data.F32[nPix] = source->pixels->data.F32[i][j];
+
+            // psMinimizeLMChi2 takes wt = 1/dY^2.  suggestion from RHL is to use the local sky
+            // as weight to avoid the bias from systematic errors here we would just use the
+            // source sky variance
+            if (PM_SOURCE_FIT_MODEL_PIX_WEIGHTS) {
+                yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j];
+            } else {
+                yErr->data.F32[nPix] = 1.0 / PM_SOURCE_FIT_MODEL_WEIGHT;
+            }
+            nPix++;
+        }
+    }
+    x->n = nPix;
+    y->n = nPix;
+    yErr->n = nPix;
+
+    thisSet = pmSourceFitSetDataAlloc (modelSet);
+
+    // define param and deriv vectors for complete set of parameters
+    psVector *params = psVectorAlloc (thisSet->nParamSet, PS_TYPE_F32);
+
+    // set the param and deriv vectors based on the curent values
+    pmSourceFitSetJoin (NULL, params, thisSet);
+
+    // create the minimization constraints
+    psMinConstraint *constraint = psMinConstraintAlloc();
+    constraint->paramMask = psVectorAlloc (thisSet->nParamSet, PS_TYPE_U8);
+    constraint->checkLimits = pmSourceFitSetCheckLimits;
+
+    pmSourceFitSetMasks (constraint, thisSet, mode);
+
+    // force the floating parameters to fall within the contraint ranges
+    for (int i = 0; i < params->n; i++) {
+        pmSourceFitSetCheckLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL);
+        pmSourceFitSetCheckLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL);
+    }
+
+    if (psTraceGetLevel("psModules.objects") >= 5) {
+        for (int i = 0; i < params->n; i++) {
+            fprintf (stderr, "%d %f %d\n", i, params->data.F32[i], constraint->paramMask->data.U8[i]);
+        }
+    }
+
+    if (nPix <  thisSet->nParamSet + 1) {
+        psTrace (__func__, 4, "insufficient valid pixels\n");
+        psTrace("psModules.objects", 3, "---- %s() end : fail pixels ----\n", __func__);
+	for (int i = 0; i < modelSet->n; i++) {
+	    pmModel *model = modelSet->data[i];
+	    model->flags |= PM_MODEL_STATUS_BADARGS;
+	}
+        psFree (x);
+        psFree (y);
+        psFree (yErr);
+        psFree (params);
+        psFree(constraint);
+        psFree (thisSet);
+	thisSet = NULL;
+        return(false);
+    }
+
+    psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_TOLERANCE);
+
+    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
+
+    fitStatus = psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, pmSourceFitSetFunction);
+    if (!fitStatus) {
+        psTrace("psModules.objects", 4, "Failed to fit model (%ld components)\n", modelSet->n);
+    }
+
+    // parameter errors from the covariance matrix
+    psVector *dparams = psVectorAlloc (thisSet->nParamSet, PS_TYPE_F32);
+    for (int i = 0; i < dparams->n; i++) {
+        if ((constraint->paramMask != NULL) && constraint->paramMask->data.U8[i])
+            continue;
+        dparams->data.F32[i] = sqrt(covar->data.F32[i][i]);
+    }
+
+    // get the Gauss-Newton distance for fixed model parameters
+    if (constraint->paramMask != NULL) {
+        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F32);
+        psVector *altmask = psVectorAlloc (params->n, PS_TYPE_U8);
+        altmask->data.U8[0] = 1;
+        for (int i = 1; i < dparams->n; i++) {
+            altmask->data.U8[i] = (constraint->paramMask->data.U8[i]) ? 0 : 1;
+        }
+        psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, pmSourceFitSetFunction);
+        for (int i = 0; i < dparams->n; i++) {
+            if (!constraint->paramMask->data.U8[i])
+                continue;
+            // note that delta is the value *subtracted* from the parameter
+            // to get the new guess.  for dparams to represent the direction
+            // of motion, we need to take -delta
+            dparams->data.F32[i] = -delta->data.F32[i];
+        }
+        psFree (delta);
+        psFree (altmask);
+    }
+
+    pmSourceFitSetValues (thisSet, dparams, params, source, myMin, y->n, fitStatus);
+    psTrace ("psModules.objects", 5, "onPic: %d, fitStatus: %d, nIter: %d, chisq: %f, nPix: %ld\n", onPic, fitStatus, myMin->iter, myMin->value, y->n);
+
+    source->mode |= PM_SOURCE_MODE_FITTED;
+
+    psFree(x);
+    psFree(y);
+    psFree(yErr);
+    psFree(myMin);
+    psFree(covar);
+    psFree(constraint);
+    psFree(params);
+    psFree(dparams);
+    psFree(thisSet);
+
+    thisSet = NULL;
+
+    bool rc = (onPic && fitStatus);
+    psTrace("psModules.objects", 5, "---- %s end (%d) ----\n", __func__, rc);
+    return(rc);
+}
Index: trunk/psModules/src/objects/pmSourceFitSet.h
===================================================================
--- trunk/psModules/src/objects/pmSourceFitSet.h	(revision 14652)
+++ trunk/psModules/src/objects/pmSourceFitSet.h	(revision 14652)
@@ -0,0 +1,53 @@
+/* @file  pmSourceFitSet.h
+ *
+ * @author EAM, IfA; GLG, MHPCC
+ *
+ * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-08-24 00:11:02 $
+ * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ */
+
+# ifndef PM_SOURCE_FIT_SET_H
+# define PM_SOURCE_FIT_SET_H
+
+/// @addtogroup Objects Object Detection / Analysis Functions
+/// @{
+
+typedef struct {
+    psArray *modelSet;
+    psArray *paramSet;
+    psArray *derivSet;
+    int nParamSet;
+} pmSourceFitSetData;
+
+// initialize data for a group of object models
+pmSourceFitSetData *pmSourceFitSetDataAlloc (psArray *modelSet);
+
+// function used to set limits for a group of models
+bool pmSourceFitSetCheckLimits (psMinConstraintMode mode, int nParam, float *params, float *betas);
+
+bool pmSourceFitSetJoin (psVector *deriv, psVector *param, pmSourceFitSetData *set);
+bool pmSourceFitSetSplit (pmSourceFitSetData *set, const psVector *deriv, const psVector *param);
+bool pmSourceFitSetValues (pmSourceFitSetData *set, const psVector *dparam, const psVector *param, pmSource *source, psMinimization *myMin, int nPix, bool fitStatus);
+
+psF32 pmSourceFitSetFunction(psVector *deriv, const psVector *param, const psVector *x);
+bool pmSourceFitSetMasks (psMinConstraint *constraint, pmSourceFitSetData *set, pmSourceFitMode mode);
+
+/** pmSourceFitSet()
+ *
+ * Fit the requested model to the specified source. The starting guess for the model is given
+ * by the input source.model parameter values. The pixels of interest are specified by the
+ * source.pixels and source.mask entries. This function calls psMinimizeLMChi2() on the image
+ * data. The function returns TRUE on success or FALSE on failure.
+ *
+ */
+bool pmSourceFitSet(
+    pmSource *source,   ///< The input pmSource
+    psArray *modelSet,   ///< model to be fitted
+    pmSourceFitMode mode,  ///< define parameters to be fitted
+    psMaskType maskVal                  ///< Vale to mask
+
+);
+
+/// @}
+# endif /* PM_SOURCE_FIT_MODEL_H */
Index: trunk/psModules/src/objects/pmSourceIO.c
===================================================================
--- trunk/psModules/src/objects/pmSourceIO.c	(revision 14530)
+++ trunk/psModules/src/objects/pmSourceIO.c	(revision 14652)
@@ -3,6 +3,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.45 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-07-14 03:20:44 $
+ *  @version $Revision: 1.46 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -34,8 +34,8 @@
 #include "pmGrowthCurve.h"
 #include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmSource.h"
-#include "pmModelGroup.h"
+#include "pmModelClass.h"
 #include "pmSourceIO.h"
 
Index: trunk/psModules/src/objects/pmSourceIO_CMP.c
===================================================================
--- trunk/psModules/src/objects/pmSourceIO_CMP.c	(revision 14530)
+++ trunk/psModules/src/objects/pmSourceIO_CMP.c	(revision 14652)
@@ -3,6 +3,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-05-03 00:13:03 $
+ *  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -32,8 +32,8 @@
 #include "pmGrowthCurve.h"
 #include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmSource.h"
-#include "pmModelGroup.h"
+#include "pmModelClass.h"
 #include "pmSourceIO.h"
 
@@ -154,9 +154,9 @@
 
     // define PSF model type
-    int modelType = pmModelSetType ("PS_MODEL_GAUSS");
+    int modelType = pmModelClassGetType ("PS_MODEL_GAUSS");
 
     char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME");
     if (PSF_NAME != NULL) {
-        modelType = pmModelSetType (PSF_NAME);
+        modelType = pmModelClassGetType (PSF_NAME);
     }
 
Index: trunk/psModules/src/objects/pmSourceIO_OBJ.c
===================================================================
--- trunk/psModules/src/objects/pmSourceIO_OBJ.c	(revision 14530)
+++ trunk/psModules/src/objects/pmSourceIO_OBJ.c	(revision 14652)
@@ -3,6 +3,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-05-03 00:13:03 $
+ *  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -32,8 +32,8 @@
 #include "pmGrowthCurve.h"
 #include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmSource.h"
-#include "pmModelGroup.h"
+#include "pmModelClass.h"
 #include "pmSourceIO.h"
 
Index: trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c
===================================================================
--- trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c	(revision 14530)
+++ trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c	(revision 14652)
@@ -3,6 +3,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-05-03 00:13:42 $
+ *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -32,8 +32,8 @@
 #include "pmGrowthCurve.h"
 #include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmSource.h"
-#include "pmModelGroup.h"
+#include "pmModelClass.h"
 #include "pmSourceIO.h"
 
@@ -152,9 +152,9 @@
 
     // define PSF model type
-    int modelType = pmModelSetType ("PS_MODEL_GAUSS");
+    int modelType = pmModelClassGetType ("PS_MODEL_GAUSS");
 
     char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME");
     if (PSF_NAME != NULL) {
-        modelType = pmModelSetType (PSF_NAME);
+        modelType = pmModelClassGetType (PSF_NAME);
     }
 
Index: trunk/psModules/src/objects/pmSourceIO_RAW.c
===================================================================
--- trunk/psModules/src/objects/pmSourceIO_RAW.c	(revision 14530)
+++ trunk/psModules/src/objects/pmSourceIO_RAW.c	(revision 14652)
@@ -3,6 +3,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -32,8 +32,8 @@
 #include "pmGrowthCurve.h"
 #include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmSource.h"
-#include "pmModelGroup.h"
+#include "pmModelClass.h"
 #include "pmSourcePhotometry.h"
 #include "pmSourceIO.h"
Index: trunk/psModules/src/objects/pmSourceIO_SMPDATA.c
===================================================================
--- trunk/psModules/src/objects/pmSourceIO_SMPDATA.c	(revision 14530)
+++ trunk/psModules/src/objects/pmSourceIO_SMPDATA.c	(revision 14652)
@@ -3,6 +3,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-05-03 00:13:42 $
+ *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -30,10 +30,10 @@
 #include "pmPeaks.h"
 #include "pmMoments.h"
+#include "pmResiduals.h"
 #include "pmGrowthCurve.h"
-#include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmSource.h"
-#include "pmModelGroup.h"
+#include "pmModelClass.h"
 #include "pmSourceIO.h"
 
@@ -132,9 +132,9 @@
 
     // define PSF model type
-    int modelType = pmModelSetType ("PS_MODEL_GAUSS");
+    int modelType = pmModelClassGetType ("PS_MODEL_GAUSS");
 
     char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME");
     if (PSF_NAME != NULL) {
-        modelType = pmModelSetType (PSF_NAME);
+        modelType = pmModelClassGetType (PSF_NAME);
     }
 
Index: trunk/psModules/src/objects/pmSourceIO_SX.c
===================================================================
--- trunk/psModules/src/objects/pmSourceIO_SX.c	(revision 14530)
+++ trunk/psModules/src/objects/pmSourceIO_SX.c	(revision 14652)
@@ -3,6 +3,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-27 22:14:08 $
+ *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -32,8 +32,8 @@
 #include "pmGrowthCurve.h"
 #include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmSource.h"
-#include "pmModelGroup.h"
+#include "pmModelClass.h"
 #include "pmSourceIO.h"
 
Index: trunk/psModules/src/objects/pmSourcePhotometry.c
===================================================================
--- trunk/psModules/src/objects/pmSourcePhotometry.c	(revision 14530)
+++ trunk/psModules/src/objects/pmSourcePhotometry.c	(revision 14652)
@@ -3,6 +3,6 @@
  *  @author EAM, IfA; GLG, MHPCC
  *
- *  @version $Revision: 1.28 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-06-20 02:22:26 $
+ *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -26,8 +26,8 @@
 #include "pmGrowthCurve.h"
 #include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmSource.h"
-#include "pmModelGroup.h"
+#include "pmModelClass.h"
 #include "pmSourcePhotometry.h"
 
@@ -250,6 +250,5 @@
 
     // measure fitMag
-    pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type);
-    fitSum = modelFluxFunc (model->params);
+    fitSum = model->modelFlux (model->params);
     if (fitSum <= 0)
         return false;
@@ -324,5 +323,4 @@
 
     // the model function returns the source flux at a position
-    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
     psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
 
@@ -354,5 +352,5 @@
 
             // for the full model, add all points
-            value = modelFunc (NULL, params, coord) - sky;
+            value = model->modelFunc (NULL, params, coord) - sky;
             modelSum += value;
 
Index: trunk/psModules/src/objects/pmSourcePlotApResid.c
===================================================================
--- trunk/psModules/src/objects/pmSourcePlotApResid.c	(revision 14530)
+++ trunk/psModules/src/objects/pmSourcePlotApResid.c	(revision 14652)
@@ -4,6 +4,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-05-22 21:45:48 $
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *  Copyright 2006 IfA, University of Hawaii
  */
@@ -27,8 +27,8 @@
 #include "pmPeaks.h"
 #include "pmMoments.h"
+#include "pmResiduals.h"
 #include "pmGrowthCurve.h"
-#include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmSource.h"
 #include "pmSourcePlots.h"
Index: trunk/psModules/src/objects/pmSourcePlotMoments.c
===================================================================
--- trunk/psModules/src/objects/pmSourcePlotMoments.c	(revision 14530)
+++ trunk/psModules/src/objects/pmSourcePlotMoments.c	(revision 14652)
@@ -5,6 +5,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-05-22 21:45:48 $
+ *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2006 IfA, University of Hawaii
@@ -32,6 +32,6 @@
 #include "pmGrowthCurve.h"
 #include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmSource.h"
 #include "pmSourcePlots.h"
Index: trunk/psModules/src/objects/pmSourcePlotPSFModel.c
===================================================================
--- trunk/psModules/src/objects/pmSourcePlotPSFModel.c	(revision 14530)
+++ trunk/psModules/src/objects/pmSourcePlotPSFModel.c	(revision 14652)
@@ -5,6 +5,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-05-22 21:45:48 $
+ *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2006 IfA, University of Hawaii
@@ -32,6 +32,6 @@
 #include "pmGrowthCurve.h"
 #include "pmResiduals.h"
+#include "pmPSF.h"
 #include "pmModel.h"
-#include "pmPSF.h"
 #include "pmSource.h"
 #include "pmSourcePlots.h"
Index: trunk/psModules/src/objects/pmSourceSky.c
===================================================================
--- trunk/psModules/src/objects/pmSourceSky.c	(revision 14530)
+++ trunk/psModules/src/objects/pmSourceSky.c	(revision 14652)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-06-20 02:22:26 $
+ *  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -27,4 +27,6 @@
 #include "pmMoments.h"
 #include "pmResiduals.h"
+#include "pmGrowthCurve.h"
+#include "pmPSF.h"
 #include "pmModel.h"
 #include "pmSource.h"
Index: trunk/psModules/src/objects/pmSourceUtils.c
===================================================================
--- trunk/psModules/src/objects/pmSourceUtils.c	(revision 14652)
+++ trunk/psModules/src/objects/pmSourceUtils.c	(revision 14652)
@@ -0,0 +1,112 @@
+/** @file  pmSource.c
+ *
+ *  Functions to define and manipulate sources on images
+ *
+ *  @author GLG, MHPCC
+ *  @author EAM, IfA: significant modifications.
+ *
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-24 00:11:02 $
+ *
+ *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <pslib.h>
+#include "pmHDU.h"
+#include "pmFPA.h"
+#include "pmFPAMaskWeight.h"
+#include "pmPeaks.h"
+#include "pmMoments.h"
+#include "pmResiduals.h"
+#include "pmGrowthCurve.h"
+#include "pmPSF.h"
+#include "pmModel.h"
+#include "pmSource.h"
+#include "pmSourceUtils.h"
+
+/******************************************************************************
+    pmSourceModelGuess(source, model): This function allocates a new
+    pmModel structure based on the given modelType specified in the argument list.
+    The corresponding pmModelGuess function is returned, and used to
+    supply the values of the params array in the pmModel structure.
+ 
+    XXX: Many parameters are based on the src->moments structure, which is in
+    image, not subImage coords.  Therefore, the calls to the model evaluation
+    functions will be in image, not subImage coords.  Remember this.
+*****************************************************************************/
+pmModel *pmSourceModelGuess(pmSource *source,
+                            pmModelType modelType)
+{
+    psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
+    PS_ASSERT_PTR_NON_NULL(source, NULL);
+    PS_ASSERT_PTR_NON_NULL(source->moments, NULL);
+    PS_ASSERT_PTR_NON_NULL(source->peak, NULL);
+
+    pmModel *model = pmModelAlloc(modelType);
+
+    if (!model->modelGuess(model, source)) {
+	psFree (model);
+	return NULL;
+    }
+
+    psTrace("psModules.objects", 3, "---- %s() end ----\n", __func__);
+    return(model);
+}
+
+pmSource *pmSourceFromModel (pmModel *model, pmReadout *readout, float radius, pmSourceType type) {
+
+    pmSource *source = pmSourceAlloc ();
+
+    // use the model centroid for the peak
+    switch (type) {
+      case PM_SOURCE_TYPE_STAR:
+	source->modelPSF = model;
+	break;
+      case PM_SOURCE_TYPE_EXTENDED:
+	source->modelEXT = model;
+	break;
+      default:
+	psAbort ("invalid source type");
+    }
+    source->type = type;
+
+    pmCell *cell = readout->parent;
+
+    float Io    = model->params->data.F32[PM_PAR_I0];
+    float xChip = model->params->data.F32[PM_PAR_XPOS];
+    float yChip = model->params->data.F32[PM_PAR_YPOS];
+
+    source->peak = pmPeakAlloc (xChip, yChip, Io, PM_PEAK_LONE);
+
+    int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0");
+    int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0");
+    int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY");
+    int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY");
+
+    // XXX fix the binning : currently not selected from concepts
+    // int xBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); // Binning in x and y
+    // int yBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN"); // Binning in x and y
+    int xBin = 1;
+    int yBin = 1;
+
+    // Position on the cell 
+    float xCell = PM_CHIP_TO_CELL(xChip, x0Cell, xParityCell, xBin);
+    float yCell = PM_CHIP_TO_CELL(yChip, y0Cell, yParityCell, yBin);
+
+    // Position on the readout
+    // float xReadout = CELL_TO_READOUT(xCell, x0Readout);
+    // float yReadout = CELL_TO_READOUT(yCell, y0Readout);
+    
+    pmSourceDefinePixels (source, readout, xCell, yCell, radius);
+
+    return (source);
+}
+
Index: trunk/psModules/src/objects/pmSourceUtils.h
===================================================================
--- trunk/psModules/src/objects/pmSourceUtils.h	(revision 14652)
+++ trunk/psModules/src/objects/pmSourceUtils.h	(revision 14652)
@@ -0,0 +1,40 @@
+/* @file  pmSourceUtils.h
+ *
+ * Utility functions for working with pmSources
+ *
+ * @author EAM, IfA
+ *
+ * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-08-24 00:11:02 $
+ * Copyright 2007 IfA, University of Hawaii
+ */
+
+# ifndef PM_SOURCE_UTILS_H
+# define PM_SOURCE_UTILS_H
+
+/// @addtogroup Objects Object Detection / Analysis Functions
+/// @{
+
+/** pmSourceModelGuess()
+ *
+ * Convert available data to an initial guess for the given model. This
+ * function allocates a pmModel entry for the pmSource structure based on the
+ * provided model selection. The method of defining the model parameter guesses
+ * are specified for each model below. The guess values are placed in the model
+ * parameters. The function returns TRUE on success or FALSE on failure.
+ *
+ */
+pmModel *pmSourceModelGuess(
+    pmSource *source,   ///< The input pmSource
+    pmModelType model   ///< The type of model to be created.
+);
+
+pmSource *pmSourceFromModel (
+  pmModel *model, 
+  pmReadout *readout, 
+  float radius,
+  pmSourceType type
+  );
+
+/// @}
+# endif /* PM_MODEL_CLASS_H */