Index: unk/psphot/src/models/pmModel_GAUSS.c
===================================================================
--- /trunk/psphot/src/models/pmModel_GAUSS.c	(revision 9805)
+++ 	(revision )
@@ -1,152 +1,0 @@
-
-/******************************************************************************
-    params->data.F32[0] = So;
-    params->data.F32[1] = Zo;
-    params->data.F32[2] = Xo;
-    params->data.F32[3] = Yo;
-    params->data.F32[4] = sqrt(2.0) / SigmaX;
-    params->data.F32[5] = sqrt(2.0) / SigmaY;
-    params->data.F32[6] = Sxy;
-*****************************************************************************/
-
-psF32 pmModelFunc_GAUSS(psVector *deriv,
-			const psVector *params,
-			const psVector *x)
-{
-    psF32 X  = x->data.F32[0] - params->data.F32[2];
-    psF32 Y  = x->data.F32[1] - params->data.F32[3];
-    psF32 px = params->data.F32[4]*X;
-    psF32 py = params->data.F32[5]*Y;
-    psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
-    psF32 r  = exp(-z);
-    psF32 q  = params->data.F32[1]*r;
-    psF32 f  = q + params->data.F32[0];
-
-    if (deriv != NULL) {
-        deriv->data.F32[0] = +1.0;
-        deriv->data.F32[1] = +r;
-        deriv->data.F32[2] = q*(2*px*params->data.F32[4] + params->data.F32[6]*Y);
-        deriv->data.F32[3] = q*(2*py*params->data.F32[5] + params->data.F32[6]*X);
-        deriv->data.F32[4] = -2.0*q*px*X;
-        deriv->data.F32[5] = -2.0*q*py*Y;
-        deriv->data.F32[6] = -q*X*Y;
-    }
-    return(f);
-}
-
-psF64 pmModelFlux_GAUSS(const psVector *params)
-{
-    psF64 A1   = PS_SQR(params->data.F32[4]);
-    psF64 A2   = PS_SQR(params->data.F32[5]);
-    psF64 A3   = PS_SQR(params->data.F32[6]);
-    psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3);
-    // Area is equivalent to 2 pi sigma^2
-
-    psF64 Flux = params->data.F32[1] * Area;
-
-    return(Flux);
-}
-
-// return the radius which yields the requested flux
-psF64 pmModelRadius_GAUSS  (const psVector *params, psF64 flux)
-{
-    
-    if (flux <= 0) return (1.0);
-    if (params->data.F32[1] <= 0) return (1.0);
-    if (flux >= params->data.F32[1]) return (1.0);
-
-    psF64 sigma  = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]);
-    psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux));
-    return (radius);
-}
-
-// define the parameter limits
-bool pmModelLimits_GAUSS (psVector **beta_lim, psVector **params_min, psVector **params_max) {
-
-    *beta_lim   = psVectorAlloc (7, PS_TYPE_F32);
-    *params_min = psVectorAlloc (7, PS_TYPE_F32);
-    *params_max = psVectorAlloc (7, PS_TYPE_F32);
-
-    beta_lim[0][0].data.F32[0] = 1000;
-    beta_lim[0][0].data.F32[1] = 10000;
-    beta_lim[0][0].data.F32[2] = 5;
-    beta_lim[0][0].data.F32[3] = 5;
-    beta_lim[0][0].data.F32[4] = 0.5;
-    beta_lim[0][0].data.F32[5] = 0.5;
-    beta_lim[0][0].data.F32[6] = 0.5;
-
-    params_min[0][0].data.F32[0] = -1000;
-    params_min[0][0].data.F32[1] = 0;
-    params_min[0][0].data.F32[2] = -100;
-    params_min[0][0].data.F32[3] = -100;
-    params_min[0][0].data.F32[4] = 0.01;
-    params_min[0][0].data.F32[5] = 0.01;
-    params_min[0][0].data.F32[6] = -5.0;
-
-    params_max[0][0].data.F32[0] = 1e5;
-    params_max[0][0].data.F32[1] = 1e6;
-    params_max[0][0].data.F32[2] = 1e4;  // this should be set by image dimensions!
-    params_max[0][0].data.F32[3] = 1e4;  // this should be set by image dimensions!
-    params_max[0][0].data.F32[4] = 2.0;
-    params_max[0][0].data.F32[5] = 2.0;
-    params_max[0][0].data.F32[6] = +5.0;
-
-    return (TRUE);
-}
-
-// make an initial guess for parameters
-bool pmModelGuess_GAUSS (pmModel *model, pmSource *source) {
-
-    pmMoments *moments = source->moments;
-    psF32     *params  = model->params->data.F32;
-
-    params[0] = moments->Sky;
-    params[1] = moments->Peak - moments->Sky;
-    params[2] = moments->x;
-    params[3] = moments->y;
-    params[4] = 1.2 / moments->Sx;
-    params[5] = 1.2 / moments->Sy;
-    params[6] = 0.0;
-    return(true);
-}
-
-// construct the PSF model from the FLT model and the psf
-bool pmModelFromPSF_GAUSS (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) {
-
-    psF32 *out = modelPSF->params->data.F32;
-    psF32 *in  = modelFLT->params->data.F32;
-    
-    out[0] = in[0];
-    out[1] = in[1];
-    out[2] = in[2];
-    out[3] = in[3];
-
-    for (int i = 4; i < 7; i++) {
-	psPolynomial2D *poly = psf->params->data[i-4];
-	out[i] = psPolynomial2DEval (poly, out[2], out[3]);
-    }
-    return(true);
-}
-
-// check the status of the fitted model
-bool pmModelFitStatus_GAUSS (pmModel *model) {
-
-    psF32 dP;
-    bool  status;
-
-    psF32 *PAR  = model->params->data.F32;
-    psF32 *dPAR = model->dparams->data.F32;
-
-    dP = 0;
-    dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
-    dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
-    dP = sqrt (dP);
-
-    status = true;
-    status &= (dP < 0.5);
-    status &= (PAR[PM_PAR_I0] > 0);
-    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);
-
-    if (status) return true;
-    return false;
-}
Index: unk/psphot/src/models/pmModel_PGAUSS.c
===================================================================
--- /trunk/psphot/src/models/pmModel_PGAUSS.c	(revision 9805)
+++ 	(revision )
@@ -1,167 +1,0 @@
-
-/******************************************************************************
-    params->data.F32[0] = So;
-    params->data.F32[1] = Zo;
-    params->data.F32[2] = Xo;
-    params->data.F32[3] = Yo;
-    params->data.F32[4] = sqrt(2.0) / SigmaX;
-    params->data.F32[5] = sqrt(2.0) / SigmaY;
-    params->data.F32[6] = Sxy;
-*****************************************************************************/
-
-psF32 pmModelFunc_PGAUSS(psVector *deriv,
-			 const psVector *params,
-			 const psVector *x)
-{
-    psF32 *PAR = params->data.F32;
-
-    psF32 X  = x->data.F32[0] - PAR[PM_PAR_XPOS];
-    psF32 Y  = x->data.F32[1] - PAR[PM_PAR_YPOS];
-    psF32 px = PAR[PM_PAR_SXX]*X;
-    psF32 py = PAR[PM_PAR_SYY]*Y;
-    psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y;
-    psF32 t  = 1 + z + z*z/2.0;
-    psF32 r  = 1.0 / (t + z*z*z/6.0); /* exp (-Z) */
-    psF32 f  = PAR[PM_PAR_I0]*r + PAR[PM_PAR_SKY];
-
-    if (deriv != NULL) {
-        // note difference from a pure gaussian: q = PAR[PM_PAR_I0]*r
-        psF32 q = PAR[PM_PAR_I0]*r*r*t;
-        deriv->data.F32[0] = +1.0;
-        deriv->data.F32[1] = +r;
-        deriv->data.F32[2] = q*(2.0*px*PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]*Y);
-        deriv->data.F32[3] = q*(2.0*py*PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]*X);
-        deriv->data.F32[4] = -2.0*q*px*X;
-        deriv->data.F32[5] = -2.0*q*py*Y;
-        deriv->data.F32[6] = -q*X*Y;
-    }
-    return(f);
-}
-
-bool pmModelLimits_PGAUSS (psVector **beta_lim, psVector **params_min, psVector **params_max) {
-
-    *beta_lim   = psVectorAlloc (7, PS_TYPE_F32);
-    *params_min = psVectorAlloc (7, PS_TYPE_F32);
-    *params_max = psVectorAlloc (7, PS_TYPE_F32);
-
-    beta_lim[0][0].data.F32[0] = 1000;
-    beta_lim[0][0].data.F32[1] = 10000;
-    beta_lim[0][0].data.F32[2] = 5;
-    beta_lim[0][0].data.F32[3] = 5;
-    beta_lim[0][0].data.F32[4] = 0.5;
-    beta_lim[0][0].data.F32[5] = 0.5;
-    beta_lim[0][0].data.F32[6] = 0.5;
-
-    params_min[0][0].data.F32[0] = -1000;
-    params_min[0][0].data.F32[1] = 0;
-    params_min[0][0].data.F32[2] = -100;
-    params_min[0][0].data.F32[3] = -100;
-    params_min[0][0].data.F32[4] = 0.01;
-    params_min[0][0].data.F32[5] = 0.01;
-    params_min[0][0].data.F32[6] = -5.0;
-
-    params_max[0][0].data.F32[0] = 1e5;
-    params_max[0][0].data.F32[1] = 1e6;
-    params_max[0][0].data.F32[2] = 1e4;  // this should be set by image dimensions!
-    params_max[0][0].data.F32[3] = 1e4;  // this should be set by image dimensions!
-    params_max[0][0].data.F32[4] = 2.0;
-    params_max[0][0].data.F32[5] = 2.0;
-    params_max[0][0].data.F32[6] = +5.0;
-
-    return (TRUE);
-}
-
-psF64 pmModelFlux_PGAUSS(const psVector *params)
-{
-    float f, norm, z;
-
-    psF32 *PAR = params->data.F32;
-
-    psF64 A1   = PS_SQR(PAR[PM_PAR_SXX]);
-    psF64 A2   = PS_SQR(PAR[PM_PAR_SYY]);
-    psF64 A3   = PS_SQR(PAR[PM_PAR_SXY]);
-    psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3);
-    // Area is equivalent to 2 pi sigma^2
-
-    // the area needs to be multiplied by the integral of f(z)
-    norm = 0.0;
-    for (z = 0.005; z < 50; z += 0.01) {
-	f = 1.0 / (1 + z + z*z/2 + z*z*z/6);
-	norm += f;
-    }
-    norm *= 0.01;
-    
-    psF64 Flux = params->data.F32[PM_PAR_I0] * Area * norm;
-
-    return(Flux);
-}
-
-// define this function so it never returns Inf or NaN
-// return the radius which yields the requested flux
-psF64 pmModelRadius_PGAUSS  (const psVector *params, psF64 flux)
-{
-    if (flux <= 0) return (1.0);
-    if (params->data.F32[1] <= 0) return (1.0);
-    if (flux >= params->data.F32[1]) return (1.0);
-
-    psF64 sigma  = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]);
-    psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux));
-    if (isnan(radius)) {
-      fprintf (stderr, "error in code\n");
-    }
-    return (radius);
-}
-
-bool pmModelGuess_PGAUSS (pmModel *model, pmSource *source) {
-
-    pmMoments *moments = source->moments;
-    psF32     *params  = model->params->data.F32;
-
-    params[0] = moments->Sky;
-    params[1] = moments->Peak - moments->Sky;
-    params[2] = moments->x;
-    params[3] = moments->y;
-    params[4] = 1.2 / moments->Sx;
-    params[5] = 1.2 / moments->Sy;
-    params[6] = 0.0;
-    return(true);
-}
-
-bool pmModelFromPSF_PGAUSS (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) {
-
-    psF32 *out = modelPSF->params->data.F32;
-    psF32 *in  = modelFLT->params->data.F32;
-    
-    out[0] = in[0];
-    out[1] = in[1];
-    out[2] = in[2];
-    out[3] = in[3];
-
-    for (int i = 4; i < 7; i++) {
-	psPolynomial2D *poly = psf->params->data[i-4];
-	out[i] = psPolynomial2DEval (poly, out[2], out[3]);
-    }
-    return(true);
-}
-
-bool pmModelFitStatus_PGAUSS (pmModel *model) {
-
-    psF32 dP;
-    bool  status;
-
-    psF32 *PAR  = model->params->data.F32;
-    psF32 *dPAR = model->dparams->data.F32;
-
-    dP = 0;
-    dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
-    dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
-    dP = sqrt (dP);
-
-    status = true;
-    status &= (dP < 0.5);
-    status &= (PAR[PM_PAR_I0] > 0);
-    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);
-
-    if (status) return true;
-    return false;
-}
Index: unk/psphot/src/models/pmModel_QGAUSS.c
===================================================================
--- /trunk/psphot/src/models/pmModel_QGAUSS.c	(revision 9805)
+++ 	(revision )
@@ -1,197 +1,0 @@
-
-/******************************************************************************
-    one component, two slopes:
-    1 / (1 + z^M + z^N)
-
-    params->data.F32[0] = So;
-    params->data.F32[1] = Zo;
-    params->data.F32[2] = Xo;
-    params->data.F32[3] = Yo;
-    params->data.F32[4] = sqrt(2.0) / SigmaX;
-    params->data.F32[5] = sqrt(2.0) / SigmaY;
-    params->data.F32[6] = Sxy;
-    params->data.F32[7] = 
-    params->data.F32[8] = 
-*****************************************************************************/
-
-/* XXX EAM : we need a way to have user-set values for fix parameters */
-# define QG_S1  2.5
-# define dQG_S1 1.5
-/* dQG_S1 is QG_S1 - 1.0 */
-
-psF32 pmModelFunc_QGAUSS(psVector *deriv,
-			 const psVector *params,
-			 const psVector *x)
-{
-    psF32 *PAR = params->data.F32;
-
-    psF32 X  = x->data.F32[0] - PAR[PM_PAR_XPOS];
-    psF32 Y  = x->data.F32[1] - PAR[PM_PAR_YPOS];
-    psF32 px = PAR[PM_PAR_SXX]*X;
-    psF32 py = PAR[PM_PAR_SYY]*Y;
-    psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y;
-
-    psF32 r  = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, QG_S1));
-    psF32 f  = PAR[PM_PAR_I0]*r + PAR[PM_PAR_SKY];
-
-    if (deriv != NULL) {
-        // note difference from a pure gaussian: q = params->data.F32[1]*r
-        psF32 t = PAR[PM_PAR_I0]*r*r;
-	psF32 q = t*(PAR[PM_PAR_7] + QG_S1*pow(z, dQG_S1));
-
-        deriv->data.F32[0] = +1.0;
-        deriv->data.F32[1] = +r;
-        deriv->data.F32[2] = q*(2.0*px*PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]*Y);
-        deriv->data.F32[3] = q*(2.0*py*PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]*X);
-        deriv->data.F32[4] = -2.0*q*px*X;
-        deriv->data.F32[5] = -2.0*q*py*Y;
-        deriv->data.F32[6] = -q*X*Y;
-        deriv->data.F32[7] = -t*z;
-    }
-    return(f);
-}
-
-bool pmModelLimits_QGAUSS (psVector **beta_lim, psVector **params_min, psVector **params_max) {
-
-    *beta_lim   = psVectorAlloc (8, PS_TYPE_F32);
-    *params_min = psVectorAlloc (8, PS_TYPE_F32);
-    *params_max = psVectorAlloc (8, PS_TYPE_F32);
-
-    beta_lim[0][0].data.F32[0] = 1000;
-    beta_lim[0][0].data.F32[1] = 10000;
-    beta_lim[0][0].data.F32[2] = 5;
-    beta_lim[0][0].data.F32[3] = 5;
-    beta_lim[0][0].data.F32[4] = 0.5;
-    beta_lim[0][0].data.F32[5] = 0.5;
-    beta_lim[0][0].data.F32[6] = 0.5;
-    beta_lim[0][0].data.F32[7] = 0.5;
-
-    params_min[0][0].data.F32[0] = -1000;
-    params_min[0][0].data.F32[1] = 0;
-    params_min[0][0].data.F32[2] = -100;
-    params_min[0][0].data.F32[3] = -100;
-    params_min[0][0].data.F32[4] = 0.01;
-    params_min[0][0].data.F32[5] = 0.01;
-    params_min[0][0].data.F32[6] = -5.0;
-    params_min[0][0].data.F32[7] = 0.1;
-
-    params_max[0][0].data.F32[0] = 1e5;
-    params_max[0][0].data.F32[1] = 1e6;
-    params_max[0][0].data.F32[2] = 1e4;  // this should be set by image dimensions!
-    params_max[0][0].data.F32[3] = 1e4;  // this should be set by image dimensions!
-    params_max[0][0].data.F32[4] = 2.0;
-    params_max[0][0].data.F32[5] = 2.0;
-    params_max[0][0].data.F32[6] = +5.0;
-    params_max[0][0].data.F32[7] = 10.0;
-
-    return (TRUE);
-}
-
-bool pmModelGuess_QGAUSS (pmModel *model, pmSource *source) {
-
-    pmMoments *moments = source->moments;
-    pmPeak    *peak    = source->peak;
-    psF32     *params  = model->params->data.F32;
-
-    params[0] = moments->Sky;
-    params[1] = moments->Peak - moments->Sky;
-    params[2] = peak->x;
-    params[3] = peak->y;
-    params[4] = 1.2 / moments->Sx;
-    params[5] = 1.2 / moments->Sy;
-    params[6] = 0.0;
-    params[7] = 1.0;
-    return(true);
-}
-
-psF64 pmModelFlux_QGAUSS(const psVector *params)
-{
-    float f, norm, z;
-
-    psF32 *PAR = params->data.F32;
-
-    psF64 A1   = PS_SQR(PAR[PM_PAR_SXX]);
-    psF64 A2   = PS_SQR(PAR[PM_PAR_SYY]);
-    psF64 A3   = PS_SQR(PAR[PM_PAR_SXY]);
-    psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3);
-    // Area is equivalent to 2 pi sigma^2
-
-    // the area needs to be multiplied by the integral of f(z)
-    norm = 0.0;
-    for (z = 0.005; z < 50; z += 0.01) {
-	f = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, QG_S1));
-	norm += f;
-    }
-    norm *= 0.01;
-    
-    psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
-
-    return(Flux);
-}
-
-// define this function so it never returns Inf or NaN
-// return the radius which yields the requested flux
-psF64 pmModelRadius_QGAUSS  (const psVector *params, psF64 flux)
-{
-    psF64 z, f;
-    psF32 *PAR = params->data.F32;
-
-    if (flux <= 0) return (1.0);
-    if (PAR[PM_PAR_I0] <= 0) return (1.0);
-    if (flux >= PAR[PM_PAR_I0]) return (1.0);
-
-    // if Sx == Sy, sigma = Sx == Sy
-    psF64 sigma = hypot (1.0 / PAR[PM_PAR_SXX], 1.0 / PAR[PM_PAR_SYY]) / sqrt(2.0);
-    psF64 dz = 1.0 / (2.0 * sigma*sigma);
-    psF64 limit = flux / PAR[PM_PAR_I0];
-
-    // we can do this much better with intelligent choices here
-    for (z = 0.0; z < 20.0; z += dz) {
-	f = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, QG_S1));
-	if (f < limit) break;
-    }
-    psF64 radius = sigma * sqrt (2.0 * z);
-    if (isnan(radius)) {
-      fprintf (stderr, "error in code\n");
-    }
-    return (radius);
-}
-
-bool pmModelFromPSF_QGAUSS (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) {
-
-    psF32 *out = modelPSF->params->data.F32;
-    psF32 *in  = modelFLT->params->data.F32;
-    
-    out[0] = in[0];
-    out[1] = in[1];
-    out[2] = in[2];
-    out[3] = in[3];
-
-    for (int i = 4; i < 8; i++) {
-	psPolynomial2D *poly = psf->params->data[i-4];
-	out[i] = psPolynomial2DEval (poly, out[2], out[3]);
-    }
-    return(true);
-}
-
-bool pmModelFitStatus_QGAUSS (pmModel *model) {
-
-    psF32 dP;
-    bool  status;
-
-    psF32 *PAR  = model->params->data.F32;
-    psF32 *dPAR = model->dparams->data.F32;
-
-    dP = 0;
-    dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
-    dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
-    dP = sqrt (dP);
-
-    status = true;
-    status &= (dP < 0.5);
-    status &= (PAR[PM_PAR_I0] > 0);
-    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);
-
-    if (status) return true;
-    return false;
-}
Index: unk/psphot/src/models/pmModel_RGAUSS.c
===================================================================
--- /trunk/psphot/src/models/pmModel_RGAUSS.c	(revision 9805)
+++ 	(revision )
@@ -1,151 +1,0 @@
-
-/******************************************************************************
-    one component, two slopes:
-    1 / (1 + z + z^Npow)
-
-    params->data.F32[0] = So;
-    params->data.F32[1] = Zo;
-    params->data.F32[2] = Xo;
-    params->data.F32[3] = Yo;
-    params->data.F32[4] = 1 / SigmaX;
-    params->data.F32[5] = 1 / SigmaY;
-    params->data.F32[6] = Sxy;
-    params->data.F32[7] = Npow
-*****************************************************************************/
-
-psF64 psModelFunc_RGAUSS(psVector *deriv,
-			 const psVector *params,
-			 const psVector *x)
-{
-    psF32 *PAR = params->data.F32;
-
-    psF32 X  = x->data.F32[0] - PAR[PM_PAR_XPOS];
-    psF32 Y  = x->data.F32[1] - PAR[PM_PAR_YPOS];
-    psF32 px = PAR[PM_PAR_SXX]*X;
-    psF32 py = PAR[PM_PAR_SYY]*Y;
-    psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y;
-
-    // psF32 FACT = 1 + 5*exp(-5*PAR[7]);
-    
-    psF32 p  = pow(z, PAR[PM_PAR_7] - 1.0);
-    psF32 r  = 1.0 / (1 + z + z*p);
-    psF32 f  = PAR[PM_PAR_I0]*r + PAR[PM_PAR_SKY];
-
-    if (deriv != NULL) {
-        // note difference from a pure gaussian: q = params->data.F32[1]*r
-        psF32 t = PAR[PM_PAR_I0]*r*r;
-	psF32 q = t*(1 + PAR[PM_PAR_7]*p);
-
-        deriv->data.F32[0] = +1.0;
-        deriv->data.F32[1] = +r;
-        deriv->data.F32[2] = q*(2.0*px*PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]*Y);
-        deriv->data.F32[3] = q*(2.0*py*PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]*X);
-        deriv->data.F32[4] = -q*px*X*2;
-        deriv->data.F32[5] = -q*py*Y*2;
-        deriv->data.F32[6] = -q*X*Y;
-        deriv->data.F32[7] = -5.0*t*log(z)*p*z;
-
-        // deriv->data.F32[4] = -1.8*q*px*X*2;
-        // deriv->data.F32[5] = -1.8*q*py*Y*2;
-        // deriv->data.F32[6] = -1.5*q*X*Y;
-        // deriv->data.F32[7] = -5.0*t*log(z)*p*z;
-    }
-    return(f);
-}
-
-psF64 psModelFlux_RGAUSS(const psVector *params)
-{
-    float f, norm, z;
-
-    psF32 *PAR = params->data.F32;
-
-    psF64 A1   = PS_SQR(PAR[PM_PAR_SXX]);
-    psF64 A2   = PS_SQR(PAR[PM_PAR_SYY]);
-    psF64 A3   = PS_SQR(PAR[PM_PAR_SXY]);
-    psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3);
-    // Area is equivalent to 2 pi sigma^2
-
-    // the area needs to be multiplied by the integral of f(z)
-    norm = 0.0;
-    for (z = 0.005; z < 50; z += 0.01) {
-	f = 1.0 / (1 + z + pow(z, PAR[PM_PAR_7]));
-	norm += f;
-    }
-    norm *= 0.01;
-    
-    psF64 Flux = params->data.F32[1] * Area * norm;
-
-    return(Flux);
-}
-
-// define this function so it never returns Inf or NaN
-// return the radius which yields the requested flux
-psF64 psModelRadius_RGAUSS  (const psVector *params, psF64 flux)
-{
-    psF64 z, f, p;
-    psF32 *PAR = params->data.F32;
-
-    if (flux <= 0) return (1.0);
-    if (PAR[PM_PAR_I0] <= 0) return (1.0);
-    if (flux >= PAR[PM_PAR_I0]) return (1.0);
-
-    // if Sx == Sy, sigma = Sx == Sy
-    psF64 sigma = hypot (1.0 / PAR[PM_PAR_SXX], 1.0 / PAR[PM_PAR_SYY]) / sqrt(2.0);
-    psF64 dz = 1.0 / (2.0 * sigma*sigma);
-    psF64 limit = flux / PAR[PM_PAR_I0];
-
-    // we can do this much better with intelligent choices here
-    for (z = 0.0; z < 20.0; z += dz) {
-	p = pow(z, PAR[PM_PAR_7]);
-	f = 1.0 / (1 + z + p);
-	if (f < limit) break;
-    }
-    psF64 radius = sigma * sqrt (2.0 * z);
-    if (isnan(radius)) {
-      fprintf (stderr, "error in code\n");
-    }
-    return (radius);
-}
-
-bool psModelGuess_RGAUSS (psModel *model, psSource *source) {
-
-    psVector *params = model->params;
-
-    EllipseAxes axes;
-    EllipseShape shape;
-    EllipseMoments moments;
-    
-    moments.x2 = PS_SQR(source->moments->Sx);
-    moments.y2 = PS_SQR(source->moments->Sy);
-    moments.xy = source->moments->Sxy;
-
-    axes = EllipseMomentsToAxes(moments);
-    shape = EllipseAxesToShape(axes);
-
-    params->data.F32[0] = source->moments->Sky;
-    params->data.F32[1] = source->peak->counts - source->moments->Sky;
-    params->data.F32[2] = source->moments->x;
-    params->data.F32[3] = source->moments->y;
-    params->data.F32[4] = 1.0 / shape.sx;
-    params->data.F32[5] = 1.0 / shape.sy;
-    params->data.F32[6] = shape.sxy;
-    params->data.F32[7] = 2.0;
-    return(true);
-}
-
-bool psModelFromPSF_RGAUSS (psModel *modelPSF, psModel *modelFLT, pmPSF *psf) {
-
-    psF32 *out = modelPSF->params->data.F32;
-    psF32 *in  = modelFLT->params->data.F32;
-    
-    out[0] = in[0];
-    out[1] = in[1];
-    out[2] = in[2];
-    out[3] = in[3];
-
-    for (int i = 4; i < 8; i++) {
-	psPolynomial2D *poly = psf->params->data[i-4];
-	out[i] = Polynomial2DEval (poly, out[2], out[3]);
-    }
-    return(true);
-}
Index: unk/psphot/src/models/pmModel_SGAUSS.c
===================================================================
--- /trunk/psphot/src/models/pmModel_SGAUSS.c	(revision 9805)
+++ 	(revision )
@@ -1,342 +1,0 @@
-
-/******************************************************************************
-    one component, two slopes:
-    1 / (1 + z^Npow + St*z^Ntide)
-
-    params->data.F32[0] = So;
-    params->data.F32[1] = Zo;
-    params->data.F32[2] = Xo;
-    params->data.F32[3] = Yo;
-    params->data.F32[4] = 1 / SigmaX;
-    params->data.F32[5] = 1 / SigmaY;
-    params->data.F32[6] = Sxy;
-    params->data.F32[7] = Npow
-    params->data.F32[8] = St
-*****************************************************************************/
-
-# define SQ(A)((A)*(A))
-psF64 psImageEllipseContour (EllipseAxes axes, double xc, double yc, psImage *image);
-psF64 p_psImageGetElementF64(psImage *a, int i, int j);
-
-psF32 pmModelFunc_SGAUSS(psVector *deriv,
-                         const psVector *params,
-                         const psVector *x)
-{
-    psF32 *PAR = params->data.F32;
-
-    psF32 X  = x->data.F32[0] - PAR[PM_PAR_XPOS];
-    psF32 Y  = x->data.F32[1] - PAR[PM_PAR_YPOS];
-    psF32 px = PAR[PM_PAR_SXX]*X;
-    psF32 py = PAR[PM_PAR_SYY]*Y;
-    psF32 z  = PS_MAX((0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y), 1e-8);
-    // note that if z -> 0, dPAR[PM_PAR_7] -> -inf
-    // also z^(PAR[PM_PAR_7]-1) -> Inf
-
-    psF32 pr = z*PAR[PM_PAR_8];
-    psF32 pr3 = pr*pr*pr;
-    psF32 p  = pow(z, PAR[PM_PAR_7] - 1.0);
-    psF32 r  = 1.0 / (1 + z*p + pr*pr3);
-    psF32 f  = PAR[PM_PAR_I0]*r + PAR[PM_PAR_SKY];
-
-    if (deriv != NULL) {
-        // note difference from a pure gaussian: q = params->data.F32[1]*r
-        psF32 t = PAR[PM_PAR_I0]*r*r;
-        psF32 q = t*(PAR[PM_PAR_7]*p + 4*PAR[PM_PAR_8]*pr3);
-
-        deriv->data.F32[0] = +1.0;
-        deriv->data.F32[1] = +r;
-        deriv->data.F32[2] = q*(2.0*px*PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]*Y);
-        deriv->data.F32[3] = q*(2.0*py*PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]*X);
-        deriv->data.F32[4] = -2.0*q*px*X;
-        deriv->data.F32[5] = -2.0*q*py*Y;
-        deriv->data.F32[6] = -q*X*Y;
-        deriv->data.F32[7] = -2*t*log(z)*z*p;
-        deriv->data.F32[8] = -2*t*4*z*pr3;
-    }
-    return(f);
-}
-
-bool pmModelLimits_SGAUSS (psVector **beta_lim, psVector **params_min, psVector **params_max) {
-
-    *beta_lim   = psVectorAlloc (9, PS_TYPE_F32);
-    *params_min = psVectorAlloc (9, PS_TYPE_F32);
-    *params_max = psVectorAlloc (9, PS_TYPE_F32);
-
-    beta_lim[0][0].data.F32[0] = 1000;
-    beta_lim[0][0].data.F32[1] = 10000;
-    beta_lim[0][0].data.F32[2] = 5;
-    beta_lim[0][0].data.F32[3] = 5;
-    beta_lim[0][0].data.F32[4] = 0.5;
-    beta_lim[0][0].data.F32[5] = 0.5;
-    beta_lim[0][0].data.F32[6] = 0.5;
-    beta_lim[0][0].data.F32[7] = 0.5;
-    beta_lim[0][0].data.F32[8] = 0.05;
-
-    params_min[0][0].data.F32[0] = -1000;
-    params_min[0][0].data.F32[1] = 0;
-    params_min[0][0].data.F32[2] = -100;
-    params_min[0][0].data.F32[3] = -100;
-    params_min[0][0].data.F32[4] = 0.01;
-    params_min[0][0].data.F32[5] = 0.01;
-    params_min[0][0].data.F32[6] = -5.0;
-    params_min[0][0].data.F32[7] = 0.5;
-    params_min[0][0].data.F32[8] = 0.05;
-
-    params_max[0][0].data.F32[0] = 1e5;
-    params_max[0][0].data.F32[1] = 1e6;
-    params_max[0][0].data.F32[2] = 1e4;  // this should be set by image dimensions!
-    params_max[0][0].data.F32[3] = 1e4;  // this should be set by image dimensions!
-    params_max[0][0].data.F32[4] = 2.0;
-    params_max[0][0].data.F32[5] = 2.0;
-    params_max[0][0].data.F32[6] = +3.0;
-    params_max[0][0].data.F32[7] = 5.0;
-    params_max[0][0].data.F32[8] = 0.5;
-
-    return (TRUE);
-}
-
-// measure the flux for the elliptical contour
-psF64 psImageEllipseContour (EllipseAxes axes, double xc, double yc, psImage *image) {
-
-    double t, dt, ct, st, xo, yo, value;
-    int N, Nt, x, y;
-
-    // choose dt to uniformly divide contour, with ~1 pix spacing at most
-    dt = asin (1 / axes.minor);
-    Nt = (int)(2*M_PI / dt) + 1;
-    dt = 2*M_PI / Nt;
-
-    ct = cos(axes.theta);
-    st = sin(axes.theta);
-    xo = xc - image->col0;
-    yo = yc - image->row0;
-
-    psVector *contour = psVectorAllocEmpty (Nt, PS_TYPE_F32);
-    for (t = 0, N = 0; (t < 2*M_PI) && (N < Nt); t += dt) {
-        x = ct*axes.major*cos(t) + st*axes.minor*sin(t) + xo;
-        y = ct*axes.minor*sin(t) + st*axes.major*cos(t) + yo;
-        value = p_psImageGetElementF64(image, x, y);
-        if (isfinite(value)) {
-            contour->data.F32[N] = value;
-            N++;
-        }
-    }
-    contour->n = N;
-    // accept every pixel: double counting is not so problematic here...
-
-    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
-    psVectorStats (stats, contour, NULL, NULL, 0);
-    value = stats->sampleMedian;
-
-    psFree (stats);
-    psFree (contour);
-
-    return (value);
-}
-
-bool pmModelGuess_SGAUSS (pmModel *model, pmSource *source) {
-
-    pmMoments *sMoments = source->moments;
-    pmPeak    *peak     = source->peak;
-    psF32     *params   = model->params->data.F32;
-
-    EllipseAxes axes;
-    EllipseShape shape;
-    EllipseMoments moments;
-
-    moments.x2 = PS_SQR(sMoments->Sx);
-    moments.y2 = PS_SQR(sMoments->Sy);
-    moments.xy = sMoments->Sxy;
-
-    // solve the math to go from Moments To Shape
-    axes = EllipseMomentsToAxes(moments);
-    shape = EllipseAxesToShape(axes);
-
-    params[0] = sMoments->Sky;
-    params[1] = sMoments->Peak - sMoments->Sky;
-    params[2] = peak->x;
-    params[3] = peak->y;
-    params[4] = 1.0 / shape.sx;
-    params[5] = 1.0 / shape.sy;
-    params[6] = shape.sxy;
-
-# if (0)
-    f1 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels);
-    axes.major *= 2.0;
-    axes.minor *= 2.0;
-    f2 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels);
-
-    if (f1 > f2) {
-        params[7] = PS_MIN (3.0, PS_MAX (0.5, log(2.0*(f1/f2) - 1.0) / log(2.0)));
-    } else {
-        params[7] = 0.5;
-    }
-# endif
-
-    params[7] = 1.8;
-    params[8] = 0.1;
-
-
-    return(true);
-}
-
-// XXX EAM : test version using flux contours to guess slope
-bool pmModelGuess_SGAUSS_HARD (pmModel *model, pmSource *source) {
-
-    pmMoments *sMoments = source->moments;
-    pmPeak    *peak     = source->peak;
-    psF32     *params   = model->params->data.F32;
-    float f1, f2;
-
-    EllipseAxes axes;
-    EllipseShape shape;
-    EllipseMoments moments;
-
-    moments.x2 = PS_SQR(sMoments->Sx);
-    moments.y2 = PS_SQR(sMoments->Sy);
-    moments.xy = sMoments->Sxy;
-
-    // solve the math to go from Moments To Shape
-    axes = EllipseMomentsToAxes(moments);
-    shape = EllipseAxesToShape(axes);
-
-    params[0] = sMoments->Sky;
-    params[1] = sMoments->Peak - sMoments->Sky;
-    params[2] = peak->x;
-    params[3] = peak->y;
-    params[4] = 1.0 / shape.sx;
-    params[5] = 1.0 / shape.sy;
-    params[6] = shape.sxy;
-
-    f1 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels);
-    axes.major *= 2.0;
-    axes.minor *= 2.0;
-    f2 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels);
-
-    if (f1 > f2) {
-        params[7] = PS_MIN (3.0, PS_MAX (0.5, log(2.0*(f1/f2) - 1.0) / log(2.0)));
-    } else {
-        params[7] = 0.5;
-    }
-    params[8] = 0.1;
-
-    return(true);
-}
-
-psF64 pmModelFlux_SGAUSS(const psVector *params)
-{
-    float f, norm, z;
-
-    psF32 *PAR = params->data.F32;
-
-    psF64 A1   = PS_SQR(PAR[PM_PAR_SXX]);
-    psF64 A2   = PS_SQR(PAR[PM_PAR_SYY]);
-    psF64 A3   = PS_SQR(PAR[PM_PAR_SXY]);
-    psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3);
-    // Area is equivalent to 2 pi sigma^2
-
-    // the area needs to be multiplied by the integral of f(z)
-    norm = 0.0;
-    for (z = 0.005; z < 50; z += 0.01) {
-      psF32 pr = PAR[PM_PAR_8]*z;
-        f = 1.0 / (1 + pow(z, PAR[PM_PAR_7]) + SQ(SQ(pr)));
-        norm += f;
-    }
-    norm *= 0.01;
-
-    psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
-
-    return(Flux);
-}
-
-// XXX need to define the radius along the major axis
-// define this function so it never returns Inf or NaN
-// return the radius which yields the requested flux
-psF64 pmModelRadius_SGAUSS  (const psVector *params, psF64 flux)
-{
-    psF64 r, z, pr, f;
-    psF32 *PAR = params->data.F32;
-
-    EllipseAxes axes;
-    EllipseShape shape;
-
-    if (flux <= 0) return (1.0);
-    if (PAR[PM_PAR_I0] <= 0) return (1.0);
-    if (flux >= PAR[PM_PAR_I0]) return (1.0);
-
-    // convert Sx,Sy,Sxy to major/minor axes
-    shape.sx = 1.0 / PAR[PM_PAR_SXX];
-    shape.sy = 1.0 / PAR[PM_PAR_SYY];
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    axes = EllipseShapeToAxes (shape);
-    psF64 dr = 1.0 / axes.major;
-    psF64 limit = flux / PAR[PM_PAR_I0];
-
-    // XXX : we can do this faster with an intelligent starting choice
-    for (r = 0.0; r < 20.0; r += dr) {
-        z = SQ(r);
-        pr = PAR[PM_PAR_8]*z;
-        f = 1.0 / (1 + pow(z, PAR[PM_PAR_7]) + SQ(SQ(pr)));
-        if (f < limit) break;
-    }
-    psF64 radius = 2.0 * axes.major * sqrt (z);
-    if (isnan(radius)) {
-        fprintf (stderr, "error in code\n");
-    }
-    return (radius);
-}
-
-bool pmModelFromPSF_SGAUSS (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) {
-
-    psF32 *out = modelPSF->params->data.F32;
-    psF32 *in  = modelFLT->params->data.F32;
-
-    out[0] = in[0];
-    out[1] = in[1];
-    out[2] = in[2];
-    out[3] = in[3];
-
-    for (int i = 4; i < 9; i++) {
-        psPolynomial2D *poly = psf->params->data[i-4];
-        out[i] = psPolynomial2DEval (poly, out[2], out[3]);
-    }
-    return(true);
-}
-
-bool pmModelFitStatus_SGAUSS (pmModel *model) {
-
-    psF32 dP;
-    bool  status;
-    EllipseAxes axes;
-    EllipseShape shape;
-
-    psF32 *PAR  = model->params->data.F32;
-    psF32 *dPAR = model->dparams->data.F32;
-
-    shape.sx = 1.0 / PAR[PM_PAR_SXX];
-    shape.sy = 1.0 / PAR[PM_PAR_SYY];
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    axes = EllipseShapeToAxes (shape);
-
-    dP = 0;
-    dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
-    dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
-    dP += PS_SQR(dPAR[PM_PAR_7] / PAR[PM_PAR_7]);
-    dP = sqrt (dP);
-
-    status = true;
-    status &= (dP < 0.5);
-    status &= (PAR[PM_PAR_I0] > 0);
-    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);
-    status &= (fabs(PAR[PM_PAR_8]) < 0.5);
-    status &= (dPAR[PM_PAR_8] < 0.1);
-    status &= (axes.major > 1.41);
-    status &= (axes.minor > 1.41);
-    status &= ((axes.major / axes.minor) < 5.0);
-    status &= (PAR[PM_PAR_7] > 0.5);
-
-    if (status) return true;
-    return false;
-}
Index: unk/psphot/src/models/pmModel_TAUSS.c
===================================================================
--- /trunk/psphot/src/models/pmModel_TAUSS.c	(revision 9805)
+++ 	(revision )
@@ -1,163 +1,0 @@
-/******************************************************************************
-    params->data.F32[0] = So;
-    params->data.F32[1] = Zo;
-    params->data.F32[2] = Xo;
-    params->data.F32[3] = Yo;
-    params->data.F32[4] = sqrt(2.0) / SigmaX;
-    params->data.F32[5] = sqrt(2.0) / SigmaY;
-    params->data.F32[6] = Sxy;
-*****************************************************************************/
-
-psF32 pmModelFunc_TAUSS(psVector *deriv,
-                        const psVector *params,
-                        const psVector *x)
-{
-    psF32 X  = x->data.F32[0] - params->data.F32[2];
-    psF32 Y  = x->data.F32[1] - params->data.F32[3];
-    psF32 px = params->data.F32[4]*X;
-    psF32 py = params->data.F32[5]*Y;
-    psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
-    psF32 r  = exp(-z);
-    psF32 q  = params->data.F32[1]*r;
-    psF32 f  = q + params->data.F32[0];
-
-    if (deriv != NULL) {
-        deriv->data.F32[0] = +1.0;
-        deriv->data.F32[1] = +r;
-        deriv->data.F32[2] = q*(2*px*params->data.F32[4] + params->data.F32[6]*Y);
-        deriv->data.F32[3] = q*(2*py*params->data.F32[5] + params->data.F32[6]*X);
-        deriv->data.F32[4] = -2.0*q*px*X;
-        deriv->data.F32[5] = -2.0*q*py*Y;
-        deriv->data.F32[6] = -q*X*Y;
-    }
-    return(f);
-}
-
-psF64 pmModelFlux_TAUSS(const psVector *params)
-{
-    psF64 A1   = PS_SQR(params->data.F32[4]);
-    psF64 A2   = PS_SQR(params->data.F32[5]);
-    psF64 A3   = PS_SQR(params->data.F32[6]);
-    psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3);
-    // Area is equivalent to 2 pi sigma^2
-
-    psF64 Flux = params->data.F32[1] * Area;
-
-    return(Flux);
-}
-
-// return the radius which yields the requested flux
-psF64 pmModelRadius_TAUSS  (const psVector *params, psF64 flux)
-{
-
-    if (flux <= 0)
-        return (1.0);
-    if (params->data.F32[1] <= 0)
-        return (1.0);
-    if (flux >= params->data.F32[1])
-        return (1.0);
-
-    psF64 sigma  = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]);
-    psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux));
-    return (radius);
-}
-
-// define the parameter limits
-bool pmModelLimits_TAUSS (psVector **beta_lim, psVector **params_min, psVector **params_max)
-{
-
-    *beta_lim   = psVectorAlloc (7, PS_TYPE_F32);
-    *params_min = psVectorAlloc (7, PS_TYPE_F32);
-    *params_max = psVectorAlloc (7, PS_TYPE_F32);
-
-    beta_lim[0][0].data.F32[0] = 1000;
-    beta_lim[0][0].data.F32[1] = 10000;
-    beta_lim[0][0].data.F32[2] = 5;
-    beta_lim[0][0].data.F32[3] = 5;
-    beta_lim[0][0].data.F32[4] = 0.5;
-    beta_lim[0][0].data.F32[5] = 0.5;
-    beta_lim[0][0].data.F32[6] = 0.5;
-
-    params_min[0][0].data.F32[0] = -1000;
-    params_min[0][0].data.F32[1] = 0;
-    params_min[0][0].data.F32[2] = -100;
-    params_min[0][0].data.F32[3] = -100;
-    params_min[0][0].data.F32[4] = 0.01;
-    params_min[0][0].data.F32[5] = 0.01;
-    params_min[0][0].data.F32[6] = -5.0;
-
-    params_max[0][0].data.F32[0] = 1e5;
-    params_max[0][0].data.F32[1] = 1e6;
-    params_max[0][0].data.F32[2] = 1e4;  // this should be set by image dimensions!
-    params_max[0][0].data.F32[3] = 1e4;  // this should be set by image dimensions!
-    params_max[0][0].data.F32[4] = 2.0;
-    params_max[0][0].data.F32[5] = 2.0;
-    params_max[0][0].data.F32[6] = +5.0;
-
-    return (TRUE);
-}
-
-// make an initial guess for parameters
-bool pmModelGuess_TAUSS (pmModel *model, pmSource *source)
-{
-
-    pmMoments *moments = source->moments;
-    psF32     *params  = model->params->data.F32;
-
-    params[0] = moments->Sky;
-    params[1] = moments->Peak - moments->Sky;
-    params[2] = moments->x;
-    params[3] = moments->y;
-    params[4] = 1.2 / moments->Sx;
-    params[5] = 1.2 / moments->Sy;
-    params[6] = 0.0;
-    return(true);
-}
-
-// construct the PSF model from the FLT model and the psf
-bool pmModelFromPSF_TAUSS (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf)
-{
-
-    psF32 *out = modelPSF->params->data.F32;
-    psF32 *in  = modelFLT->params->data.F32;
-
-    out[0] = in[0];
-    out[1] = in[1];
-    out[2] = in[2];
-    out[3] = in[3];
-
-    for (int i = 4; i < 7; i++) {
-        psPolynomial2D *poly = psf->params->data[i-4];
-        // XXX: be wary of a bug here.  EAM had his own version of psPolynomial2DEval().
-        // I think the reason was that hewas using nX=xOrder.  Not sure.  I changed this
-        // without verifying.
-        // out[i] = Polynomial2DEval_EAM (poly, out[2], out[3]);
-        out[i] = psPolynomial2DEval(poly, out[2], out[3]);
-    }
-    return(true);
-}
-
-// check the status of the fitted model
-bool pmModelFitStatus_TAUSS (pmModel *model)
-{
-
-    psF32 dP;
-    bool  status;
-
-    psF32 *PAR  = model->params->data.F32;
-    psF32 *dPAR = model->dparams->data.F32;
-
-    dP = 0;
-    dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
-    dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
-    dP = sqrt (dP);
-
-    status = true;
-    status &= (dP < 0.5);
-    status &= (PAR[PM_PAR_I0] > 0);
-    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);
-
-    if (status)
-        return true;
-    return false;
-}
Index: unk/psphot/src/models/pmModel_TGAUSS.c
===================================================================
--- /trunk/psphot/src/models/pmModel_TGAUSS.c	(revision 9805)
+++ 	(revision )
@@ -1,200 +1,0 @@
-
-/******************************************************************************
-    one component, two slopes:
-    1 / (1 + z^M + z^N)
-
-    params->data.F32[0] = So;
-    params->data.F32[1] = Zo;
-    params->data.F32[2] = Xo;
-    params->data.F32[3] = Yo;
-    params->data.F32[4] = sqrt(2.0) / SigmaX;
-    params->data.F32[5] = sqrt(2.0) / SigmaY;
-    params->data.F32[6] = Sxy;
-    params->data.F32[7] = 
-    params->data.F32[8] = 
-*****************************************************************************/
-
-/* XXX EAM : we need a way to have user-set values for fix parameters */
-# define TG_S  2.60
-# define dTG_S 1.60
-# define TRF 0.50
-/* dTG_S is TG_S - 1.0 */
-
-psF32 pmModelFunc_TGAUSS(psVector *deriv,
-			 const psVector *params,
-			 const psVector *x)
-{
-    psF32 *PAR = params->data.F32;
-
-    psF32 X  = x->data.F32[0] - PAR[PM_PAR_XPOS];
-    psF32 Y  = x->data.F32[1] - PAR[PM_PAR_YPOS];
-    psF32 px = PAR[PM_PAR_SXX]*X;
-    psF32 py = PAR[PM_PAR_SYY]*Y;
-    psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y;
-    psF32 er = exp(TRF*z);
-
-    psF32 r  = 1.0 / (er + PAR[PM_PAR_7]*z + pow(z, TG_S));    // (1/R)
-    psF32 f  = PAR[PM_PAR_I0]*r + PAR[PM_PAR_SKY];
-
-    if (deriv != NULL) {
-	psF32 t = PAR[PM_PAR_I0]*r*r;	// df/dR
-	psF32 q = t*(TRF*er + PAR[PM_PAR_7] + TG_S*pow(z, dTG_S));  // (df/dR)(dR/dz)
-
-        deriv->data.F32[0] = +1.0;
-        deriv->data.F32[1] = +r;
-        deriv->data.F32[2] = q*(2.0*px*PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]*Y);
-	deriv->data.F32[3] = q*(2.0*py*PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]*X);
-        deriv->data.F32[4] = -2.0*q*px*X;
-        deriv->data.F32[5] = -2.0*q*py*Y;
-        deriv->data.F32[6] = -q*X*Y;
-        deriv->data.F32[7] = -t*z;
-    }
-    return(f);
-}
-
-bool pmModelLimits_TGAUSS (psVector **beta_lim, psVector **params_min, psVector **params_max) {
-
-    *beta_lim   = psVectorAlloc (8, PS_TYPE_F32);
-    *params_min = psVectorAlloc (8, PS_TYPE_F32);
-    *params_max = psVectorAlloc (8, PS_TYPE_F32);
-
-    beta_lim[0][0].data.F32[0] = 1000;
-    beta_lim[0][0].data.F32[1] = 10000;
-    beta_lim[0][0].data.F32[2] = 5;
-    beta_lim[0][0].data.F32[3] = 5;
-    beta_lim[0][0].data.F32[4] = 0.5;
-    beta_lim[0][0].data.F32[5] = 0.5;
-    beta_lim[0][0].data.F32[6] = 0.5;
-    beta_lim[0][0].data.F32[7] = 0.5;
-
-    params_min[0][0].data.F32[0] = -1000;
-    params_min[0][0].data.F32[1] = 0;
-    params_min[0][0].data.F32[2] = -100;
-    params_min[0][0].data.F32[3] = -100;
-    params_min[0][0].data.F32[4] = 0.01;
-    params_min[0][0].data.F32[5] = 0.01;
-    params_min[0][0].data.F32[6] = -5.0;
-    params_min[0][0].data.F32[7] = 0.1;
-
-    params_max[0][0].data.F32[0] = 1e5;
-    params_max[0][0].data.F32[1] = 1e6;
-    params_max[0][0].data.F32[2] = 1e4;  // this should be set by image dimensions!
-    params_max[0][0].data.F32[3] = 1e4;  // this should be set by image dimensions!
-    params_max[0][0].data.F32[4] = 2.0;
-    params_max[0][0].data.F32[5] = 2.0;
-    params_max[0][0].data.F32[6] = +5.0;
-    params_max[0][0].data.F32[7] = 10.0;
-
-    return (TRUE);
-}
-
-bool pmModelGuess_TGAUSS (pmModel *model, pmSource *source) {
-
-    pmMoments *moments = source->moments;
-    pmPeak    *peak    = source->peak;
-    psF32     *params  = model->params->data.F32;
-
-    params[0] = moments->Sky;
-    params[1] = moments->Peak - moments->Sky;
-    params[2] = peak->x;
-    params[3] = peak->y;
-    params[4] = 1.2 / moments->Sx;
-    params[5] = 1.2 / moments->Sy;
-    params[6] = 0.0;
-    params[7] = 1.0;
-    return(true);
-}
-
-psF64 pmModelFlux_TGAUSS(const psVector *params)
-{
-    float f, norm, z;
-
-    psF32 *PAR = params->data.F32;
-
-    psF64 A1   = PS_SQR(PAR[PM_PAR_SXX]);
-    psF64 A2   = PS_SQR(PAR[PM_PAR_SYY]);
-    psF64 A3   = PS_SQR(PAR[PM_PAR_SXY]);
-    psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3);
-    // Area is equivalent to 2 pi sigma^2
-
-    // the area needs to be multiplied by the integral of f(z)
-    norm = 0.0;
-    for (z = 0.005; z < 50; z += 0.01) {
-        psF32 er = exp(TRF*z);
-	f = 1.0 / (er + PAR[PM_PAR_7]*z + pow(z, TG_S));
-	norm += f;
-    }
-    norm *= 0.01;
-    
-    psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
-
-    return(Flux);
-}
-
-// define this function so it never returns Inf or NaN
-// return the radius which yields the requested flux
-psF64 pmModelRadius_TGAUSS  (const psVector *params, psF64 flux)
-{
-    psF64 z, f;
-    psF32 *PAR = params->data.F32;
-
-    if (flux <= 0) return (1.0);
-    if (PAR[PM_PAR_I0] <= 0) return (1.0);
-    if (flux >= PAR[PM_PAR_I0]) return (1.0);
-
-    // if Sx == Sy, sigma = Sx == Sy
-    psF64 sigma = hypot (1.0 / PAR[PM_PAR_SXX], 1.0 / PAR[PM_PAR_SYY]) / sqrt(2.0);
-    psF64 dz = 1.0 / (2.0 * sigma*sigma);
-    psF64 limit = flux / PAR[PM_PAR_I0];
-
-    // we can do this much better with intelligent choices here
-    for (z = 0.0; z < 20.0; z += dz) {
-        psF32 er = exp(TRF*z);
-	f = 1.0 / (er + PAR[PM_PAR_7]*z + pow(z, TG_S));
-	if (f < limit) break;
-    }
-    psF64 radius = sigma * sqrt (2.0 * z);
-    if (isnan(radius)) {
-      fprintf (stderr, "error in code\n");
-    }
-    return (radius);
-}
-
-bool pmModelFromPSF_TGAUSS (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) {
-
-    psF32 *out = modelPSF->params->data.F32;
-    psF32 *in  = modelFLT->params->data.F32;
-    
-    out[0] = in[0];
-    out[1] = in[1];
-    out[2] = in[2];
-    out[3] = in[3];
-
-    for (int i = 4; i < 8; i++) {
-	psPolynomial2D *poly = psf->params->data[i-4];
-	out[i] = psPolynomial2DEval (poly, out[2], out[3]);
-    }
-    return(true);
-}
-
-bool pmModelFitStatus_TGAUSS (pmModel *model) {
-
-    psF32 dP;
-    bool  status;
-
-    psF32 *PAR  = model->params->data.F32;
-    psF32 *dPAR = model->dparams->data.F32;
-
-    dP = 0;
-    dP += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
-    dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
-    dP = sqrt (dP);
-
-    status = true;
-    status &= (dP < 0.5);
-    status &= (PAR[PM_PAR_I0] > 0);
-    status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5);
-
-    if (status) return true;
-    return false;
-}
Index: unk/psphot/src/models/pmModel_WAUSS.c
===================================================================
--- /trunk/psphot/src/models/pmModel_WAUSS.c	(revision 9805)
+++ 	(revision )
@@ -1,99 +1,0 @@
-
-/******************************************************************************
-    params->data.F32[0] = So;
-    params->data.F32[1] = Zo;
-    params->data.F32[2] = Xo;
-    params->data.F32[3] = Yo;
-    params->data.F32[4] = sqrt(2.0) / SigmaX;
-    params->data.F32[5] = sqrt(2.0) / SigmaY;
-    params->data.F32[6] = Sxy;
-*****************************************************************************/
-
-psF64 psModelFunc_WAUSS(psVector *deriv,
-			const psVector *params,
-			const psVector *x)
-{
-    psF32 X = x->data.F32[0] - params->data.F32[2];
-    psF32 Y = x->data.F32[1] - params->data.F32[2];
-    psF32 px = params->data.F32[4]*X;
-    psF32 py = params->data.F32[5]*Y;
-    psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
-    psF32 t = 0.5*z*z*(1.0 + params->data.F32[8]*z/3.0);
-    psF32 r = 1.0 / (1.0 + z + params->data.F32[7]*t); /* exp (-Z) */
-    psF32 f = params->data.F32[1]*r + params->data.F32[0];
-
-    if (deriv != NULL) {
-        // note difference from gaussian: q = params->data.F32[1]*r
-        psF32 q = params->data.F32[1]*r*r*(1.0 + params->data.F32[7]*z*(1.0 + params->data.F32[8]*z/2.0));
-        deriv->data.F32[0] = +1.0;
-        deriv->data.F32[1] = +r;
-        deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y);
-        deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X);
-        deriv->data.F32[4] = -2.0*q*px*X;
-        deriv->data.F32[5] = -2.0*q*py*Y;
-        deriv->data.F32[6] = -q*X*Y;
-        deriv->data.F32[7] = -100.0*params->data.F32[1]*r*r*t;
-        deriv->data.F32[8] = -100.0*params->data.F32[1]*r*r*params->data.F32[7]*(z*z*z)/6.0;
-        // The values of 100 dampen the swing of params->data.F32[7,8] */
-    }
-    return(f);
-}
-
-// this is probably wrong since it uses the gauss integral 2 pi sigma^2
-psF64 psModelFlux_WAUSS(const psVector *params)
-{
-    psF64 A1   = 1 / PS_SQR(params->data.F32[4]);
-    psF64 A2   = 1 / PS_SQR(params->data.F32[5]);
-    psF64 A3   = params->data.F32[6];
-    psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - PS_SQR(A3));
-    // Area is equivalent to 2 pi sigma^2
-
-    psF64 Flux = params->data.F32[1] * Area;
-
-    return(Flux);
-}
-
-// return the radius which yields the requested flux
-psF64 psModelRadius_WAUSS  (const psVector *params, psF64 flux)
-{
-    if (flux <= 0) return (1.0);
-    if (params->data.F32[1] <= 0) return (1.0);
-
-    psF64 sigma  = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]);
-    psF64 radius = sigma * sqrt (2.0 * log(params->data.F32[1] / flux));
-    return (radius);
-}
-
-bool psModelGuess_WAUSS (psModel *model, pmSource *source) {
-
-    psVector *params = model->params;
-
-    params->data.F32[0] = source->moments->Sky;
-    params->data.F32[1] = source->peak->counts - source->moments->Sky;
-    params->data.F32[2] = source->moments->x;
-    params->data.F32[3] = source->moments->y;
-    params->data.F32[4] = sqrt(2.0) / source->moments->Sx;
-    params->data.F32[5] = sqrt(2.0) / source->moments->Sy;
-    params->data.F32[6] = source->moments->Sxy;
-    // XXX: What are these?
-    // params->data.F32[7] = B2;
-    // params->data.F32[8] = B3;
-    return(true);
-}
-
-bool psModelFromPSF_WAUSS (psModel *modelPSF, psModel *modelFLT, pmPSF *psf) {
-
-    psF32 *out = modelPSF->params->data.F32;
-    psF32 *in  = modelFLT->params->data.F32;
-    
-    out[0] = in[0];
-    out[1] = in[1];
-    out[2] = in[2];
-    out[3] = in[3];
-
-    for (int i = 4; i < 9; i++) {
-	psPolynomial2D *poly = psf->params->data[i-4];
-	out[i] = Polynomial2DEval (poly, out[2], out[3]);
-    }
-    return(true);
-}
Index: unk/psphot/src/models/pmModel_ZGAUSS.c
===================================================================
--- /trunk/psphot/src/models/pmModel_ZGAUSS.c	(revision 9805)
+++ 	(revision )
@@ -1,158 +1,0 @@
-
-/******************************************************************************
-    one component, two slopes:
-    1 / (1 + z^Npow + PAR8*z^4)
-
-    params->data.F32[0] = So;
-    params->data.F32[1] = Zo;
-    params->data.F32[2] = Xo;
-    params->data.F32[3] = Yo;
-    params->data.F32[4] = 1 / SigmaX;
-    params->data.F32[5] = 1 / SigmaY;
-    params->data.F32[6] = Sxy;
-    params->data.F32[7] = Npow
-*****************************************************************************/
-
-# define SQ(A)((A)*(A))
-# define PAR8 0.1
-
-psF64 psModelFunc_ZGAUSS(psVector *deriv,
-			 const psVector *params,
-			 const psVector *x)
-{
-    psF32 *PAR = params->data.F32;
-
-    psF32 X  = x->data.F32[0] - PAR[PM_PAR_XPOS];
-    psF32 Y  = x->data.F32[1] - PAR[PM_PAR_YPOS];
-    psF32 px = PAR[PM_PAR_SXX]*X;
-    psF32 py = PAR[PM_PAR_SYY]*Y;
-    psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y;
-
-    psF32 pr = PAR8*z;
-    psF32 p  = pow(z, PAR[PM_PAR_7] - 1.0);
-    psF32 r  = 1.0 / (1 + z*p + SQ(SQ(pr)));
-    psF32 f  = PAR[PM_PAR_I0]*r + PAR[PM_PAR_SKY];
-
-    if (deriv != NULL) {
-        // note difference from a pure gaussian: q = params->data.F32[1]*r
-        psF32 t = PAR[PM_PAR_I0]*r*r;
-	psF32 q = t*(PAR[PM_PAR_7]*p + 4*PAR8*pr*pr*pr);
-
-        deriv->data.F32[0] = +1.0;
-        deriv->data.F32[1] = +r;
-        deriv->data.F32[2] = q*(2.0*px*PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]*Y);
-        deriv->data.F32[3] = q*(2.0*py*PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]*X);
-        deriv->data.F32[4] = -q*px*X;
-        deriv->data.F32[5] = -q*py*Y;
-        deriv->data.F32[6] = -q*X*Y;
-        deriv->data.F32[7] = -t*log(z)*z*p;
-    }
-    return(f);
-}
-
-psF64 psModelFlux_ZGAUSS(const psVector *params)
-{
-    float f, norm, z;
-
-    psF32 *PAR = params->data.F32;
-
-    psF64 A1   = PS_SQR(PAR[PM_PAR_SXX]);
-    psF64 A2   = PS_SQR(PAR[PM_PAR_SYY]);
-    psF64 A3   = PS_SQR(PAR[PM_PAR_SXY]);
-    psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3);
-    // Area is equivalent to 2 pi sigma^2
-
-    // the area needs to be multiplied by the integral of f(z)
-    norm = 0.0;
-    psF32 pr = PAR8*z;
-    for (z = 0.005; z < 50; z += 0.01) {
-	f = 1.0 / (1 + pow(z, PAR[PM_PAR_7]) + SQ(SQ(pr)));
-	norm += f;
-    }
-    norm *= 0.01;
-    
-    psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
-
-    return(Flux);
-}
-
-// XXX need to define the radius along the major axis
-// define this function so it never returns Inf or NaN
-// return the radius which yields the requested flux
-psF64 psModelRadius_ZGAUSS  (const psVector *params, psF64 flux)
-{
-    psF64 r, z, pr, f;
-    psF32 *PAR = params->data.F32;
-
-    EllipseAxes axes;
-    EllipseShape shape;
-
-    if (flux <= 0) return (1.0);
-    if (PAR[PM_PAR_I0] <= 0) return (1.0);
-    if (flux >= PAR[PM_PAR_I0]) return (1.0);
-
-    // convert Sx,Sy,Sxy to major/minor axes
-    shape.sx = 1.0 / PAR[PM_PAR_SXX];
-    shape.sy = 1.0 / PAR[PM_PAR_SYY];
-    shape.sxy = PAR[PM_PAR_SXY];
-
-    axes = EllipseShapeToAxes (shape);
-    psF64 dr = 1.0 / axes.major;
-    psF64 limit = flux / PAR[PM_PAR_I0];
-
-    // XXX : we can do this faster with an intelligent starting choice
-    for (r = 0.0; r < 20.0; r += dr) {
-	z = SQ(r);
-	pr = PAR8*z;
-	f = 1.0 / (1 + pow(z, PAR[PM_PAR_7]) + SQ(SQ(pr)));
-	if (f < limit) break;
-    }
-    psF64 radius = 2.0 * axes.major * sqrt (z);
-    if (isnan(radius)) {
-	fprintf (stderr, "error in code\n");
-    }
-    return (radius);
-}
-
-bool psModelGuess_ZGAUSS (psModel *model, psSource *source) {
-
-    psVector *params = model->params;
-
-    EllipseAxes axes;
-    EllipseShape shape;
-    EllipseMoments moments;
-    
-    moments.x2 = PS_SQR(source->moments->Sx);
-    moments.y2 = PS_SQR(source->moments->Sy);
-    moments.xy = source->moments->Sxy;
-
-    axes = EllipseMomentsToAxes(moments);
-    shape = EllipseAxesToShape(axes);
-
-    params->data.F32[0] = source->moments->Sky;
-    params->data.F32[1] = source->peak->counts - source->moments->Sky;
-    params->data.F32[2] = source->moments->x;
-    params->data.F32[3] = source->moments->y;
-    params->data.F32[4] = 1.0 / shape.sx;
-    params->data.F32[5] = 1.0 / shape.sy;
-    params->data.F32[6] = shape.sxy;
-    params->data.F32[7] = 1.9;
-    return(true);
-}
-
-bool psModelFromPSF_ZGAUSS (psModel *modelPSF, psModel *modelFLT, pmPSF *psf) {
-
-    psF32 *out = modelPSF->params->data.F32;
-    psF32 *in  = modelFLT->params->data.F32;
-    
-    out[0] = in[0];
-    out[1] = in[1];
-    out[2] = in[2];
-    out[3] = in[3];
-
-    for (int i = 4; i < 8; i++) {
-	psPolynomial2D *poly = psf->params->data[i-4];
-	out[i] = Polynomial2DEval (poly, out[2], out[3]);
-    }
-    return(true);
-}