Index: trunk/psphot/src/models/pmModel_TGAUSS.c
===================================================================
--- trunk/psphot/src/models/pmModel_TGAUSS.c	(revision 4954)
+++ trunk/psphot/src/models/pmModel_TGAUSS.c	(revision 5130)
@@ -12,7 +12,14 @@
     params->data.F32[6] = Sxy;
     params->data.F32[7] = 
+    params->data.F32[8] = 
 *****************************************************************************/
 
-psF64 psModelFunc_TGAUSS(psVector *deriv,
+/* 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)
@@ -25,35 +32,100 @@
     psF32 py = PAR[5]*Y;
     psF32 z  = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y;
-    
-    psF32 p  = pow(z, 1.2);
-    psF32 r  = 1.0 / (1 + z + PAR[7]*z*p);
+    psF32 er = exp(TRF*z);
+
+    psF32 r  = 1.0 / (er + PAR[7]*z + pow(z, TG_S));    // (1/R)
     psF32 f  = PAR[1]*r + PAR[0];
 
     if (deriv != NULL) {
-        // note difference from a pure gaussian: q = params->data.F32[1]*r
-	psF32 t = PAR[1]*r*r;
-	psF32 q = t*(1 + PAR[7]*2.2*p);
+	psF32 t = PAR[1]*r*r;	                         // df/dR
+	psF32 q = t*(TRF*er + 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[4] + PAR[6]*Y);
-        deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[6]*X);
+	deriv->data.F32[3] = q*(2.0*py*PAR[5] + PAR[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] = -t*z*p;
+        deriv->data.F32[7] = -t*z;
     }
     return(f);
 }
 
-psF64 psModelFlux_TGAUSS(const psVector *params)
+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)
 {
-    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));
+    float f, norm, z;
+
+    psF32 *PAR = params->data.F32;
+
+    psF64 A1   = PS_SQR(PAR[4]);
+    psF64 A2   = PS_SQR(PAR[5]);
+    psF64 A3   = PS_SQR(PAR[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;
+    // 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[7]*z + pow(z, TG_S));
+	norm += f;
+    }
+    norm *= 0.01;
+    
+    psF64 Flux = PAR[1] * Area * norm;
 
     return(Flux);
@@ -62,12 +134,25 @@
 // define this function so it never returns Inf or NaN
 // return the radius which yields the requested flux
-psF64 psModelRadius_TGAUSS  (const psVector *params, psF64 flux)
+psF64 pmModelRadius_TGAUSS  (const psVector *params, psF64 flux)
 {
+    psF64 z, f;
+    psF32 *PAR = params->data.F32;
+
     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 (PAR[1] <= 0) return (1.0);
+    if (flux >= PAR[1]) return (1.0);
+
+    // if Sx == Sy, sigma = Sx == Sy
+    psF64 sigma = hypot (1.0 / PAR[4], 1.0 / PAR[5]) / sqrt(2.0);
+    psF64 dz = 1.0 / (2.0 * sigma*sigma);
+    psF64 limit = flux / PAR[1];
+
+    // 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[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");
@@ -76,21 +161,5 @@
 }
 
-bool psModelGuess_TGAUSS (psModel *model, psSource *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] = 1.0/source->moments->Sx;
-    params->data.F32[5] = 1.0/source->moments->Sy;
-    // params->data.F32[6] = source->moments->Sxy;
-    params->data.F32[6] = 0.0;
-    params->data.F32[7] = 5.0;
-    return(true);
-}
-
-bool psModelFromPSF_TGAUSS (psModel *modelPSF, psModel *modelFLT, pmPSF *psf) {
+bool pmModelFromPSF_TGAUSS (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) {
 
     psF32 *out = modelPSF->params->data.F32;
@@ -104,6 +173,28 @@
     for (int i = 4; i < 8; i++) {
 	psPolynomial2D *poly = psf->params->data[i-4];
-	out[i] = Polynomial2DEval (poly, out[2], out[3]);
+	out[i] = Polynomial2DEval_EAM (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[4] / PAR[4]);
+    dP += PS_SQR(dPAR[5] / PAR[5]);
+    dP = sqrt (dP);
+
+    status = true;
+    status &= (dP < 0.5);
+    status &= (PAR[1] > 0);
+    status &= ((dPAR[1]/PAR[1]) < 0.5);
+
+    if (status) return true;
+    return false;
+}