Changeset 9775 for trunk/psModules/src/objects/models/pmModel_TGAUSS.c

Timestamp:

Oct 29, 2006, 5:02:59 PM (20 years ago)

Author:

magnier

Message:

major cleanup of model code, fixes to use new param defs

File:

• trunk/psModules/src/objects/models/pmModel_TGAUSS.c (modified) (4 diffs)

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

@@ -1 @@
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
 /******************************************************************************
-    one component, two slopes:
-    1 / (1 + z^M + z^N)
+ * this file defines the TGAUSS source shape model (XXX need a better name!).  Note that these
+ * model functions are loaded by pmModelGroup.c using 'include', and thus need no 'include'
+ * statements of their own.  The models use a psVector to represent the set of parameters, with
+ * the sequence used to specify the meaning of the parameter.  The meaning of the parameters
+ * may thus vary depending on the specifics of the model.  All models which are used a PSF
+ * representations share a few parameters, for which # define names are listed in pmModel.h:
  
-    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] =
-*****************************************************************************/
-
-psF64 psModelFunc_TGAUSS(psVector *deriv,
-                         const psVector *params,
-                         const psVector *x)
+   power-law with fixed slope and fitted amplitude
+   1 / (1 + z + kz^2.2)
+ 
+   * PM_PAR_SKY 0   - local sky : note that this is unused and may be dropped in the future
+   * PM_PAR_I0 1    - central intensity
+   * PM_PAR_XPOS 2  - X center of object
+   * PM_PAR_YPOS 3  - Y center of object
+   * PM_PAR_SXX 4   - X^2 term of elliptical contour (sqrt(2) / SigmaX)
+   * PM_PAR_SYY 5   - Y^2 term of elliptical contour (sqrt(2) / SigmaY)
+   * PM_PAR_SXY 6   - X*Y term of elliptical contour
+   * PM_PAR_7   7   - amplitude of high-order component (k)
+   *****************************************************************************/
+
+/***
+    XXXX the model in this file needs to be tested more carefully.
+    fix up the code to follow conventions in the other model function files.
+***/
+
+XXX broken code
+
+# define PM_MODEL_FUNC       pmModelFunc_TGAUSS
+# define PM_MODEL_FLUX       pmModelFlux_TGAUSS
+# define PM_MODEL_GUESS      pmModelGuess_TGAUSS
+# define PM_MODEL_LIMITS     pmModelLimits_TGAUSS
+# define PM_MODEL_RADIUS     pmModelRadius_TGAUSS
+# define PM_MODEL_FROM_PSF   pmModelFromPSF_TGAUSS
+# define PM_MODEL_FIT_STATUS pmModelFitStatus_TGAUSS
+
+psF64 PS_MODEL_FUNC (psVector *deriv,
+                     const psVector *params,
+                     const psVector *x)
 {
     psF32 *PAR = params->data.F32;
@@ -50 +68 @@
 }
 
-psF64 psModelFlux_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));
-    // Area is equivalent to 2 pi sigma^2
-
-    psF64 Flux = params->data.F32[1] * Area;
-
-    return(Flux);
-}
-
-// 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)
-{
-    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 psModelGuess_TGAUSS (psModel *model, psSource *source)
+bool PM_MODEL_LIMITS  (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[PM_PAR_SKY] = 1000;
+    beta_lim[0][0].data.F32[PM_PAR_I0] = 3e6;
+    beta_lim[0][0].data.F32[PM_PAR_XPOS] = 5;
+    beta_lim[0][0].data.F32[PM_PAR_YPOS] = 5;
+    beta_lim[0][0].data.F32[PM_PAR_SXX] = 0.5;
+    beta_lim[0][0].data.F32[PM_PAR_SYY] = 0.5;
+    beta_lim[0][0].data.F32[PM_PAR_SXY] = 0.5;
+    beta_lim[0][0].data.F32[PM_PAR_7] = 0.5;
+
+    params_min[0][0].data.F32[PM_PAR_SKY] = -1000;
+    params_min[0][0].data.F32[PM_PAR_I0] = 0;
+    params_min[0][0].data.F32[PM_PAR_XPOS] = -100;
+    params_min[0][0].data.F32[PM_PAR_YPOS] = -100;
+    params_min[0][0].data.F32[PM_PAR_SXX] = 0.5;
+    params_min[0][0].data.F32[PM_PAR_SYY] = 0.5;
+    params_min[0][0].data.F32[PM_PAR_SXY] = -5.0;
+    params_min[0][0].data.F32[PM_PAR_7] = 0.1;
+
+    params_max[0][0].data.F32[PM_PAR_SKY] = 1e5;
+    params_max[0][0].data.F32[PM_PAR_I0] = 1e8;
+    params_max[0][0].data.F32[PM_PAR_XPOS] = 1e4;  // this should be set by image dimensions!
+    params_max[0][0].data.F32[PM_PAR_YPOS] = 1e4;  // this should be set by image dimensions!
+    params_max[0][0].data.F32[PM_PAR_SXX] = 100.0;
+    params_max[0][0].data.F32[PM_PAR_SYY] = 100.0;
+    params_max[0][0].data.F32[PM_PAR_SXY] = +5.0;
+    params_max[0][0].data.F32[PM_PAR_7] = 10.0;
+
+    return (TRUE);
+}
+
+bool PS_MODEL_GUESS  (psModel *model, psSource *source)
 {
 
@@ -99 +122 @@
 }
 
-bool psModelFromPSF_TGAUSS (psModel *modelPSF, psModel *modelFLT, pmPSF *psf)
+psF64 PS_MODEL_FLUX (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);
+}
+
+// define this function so it never returns Inf or NaN
+// return the radius which yields the requested flux
+psF64 PS_MODEL_RADIUS   (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 PS_MODEL_FROM_PSF  (psModel *modelPSF, psModel *modelFLT, pmPSF *psf)
 {
 
@@ -116 +171 @@
     return(true);
 }
+
+bool PM_MODEL_FIT_STATUS (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);
+
+    return status;
+}
+
+# undef PM_MODEL_FUNC
+# undef PM_MODEL_FLUX
+# undef PM_MODEL_GUESS
+# undef PM_MODEL_LIMITS
+# undef PM_MODEL_RADIUS
+# undef PM_MODEL_FROM_PSF
+# undef PM_MODEL_FIT_STATUS