Changeset 9775 for trunk/psModules/src/objects/models/pmModel_WAUSS.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_WAUSS.c (modified) (4 diffs)

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

@@ -1 @@
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+/******************************************************************************
+ * this file defines the WAUSS 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:
+ 
+   power-law with fitted linear term
+   1 / (1 + Az + Bz^2 + z^3/6)
+ 
+   * 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 (SigmaX / sqrt(2))
+   * PM_PAR_SYY 5   - Y^2 term of elliptical contour (SigmaY / sqrt(2))
+   * PM_PAR_SXY 6   - X*Y term of elliptical contour
+   * PM_PAR_7   7   - amplitude of the linear component (A)
+   * PM_PAR_8   8   - amplitude of the quadratic component (B)
+   *****************************************************************************/
 
-/******************************************************************************
-    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;
-*****************************************************************************/
+/***
+    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.
+***/
 
-psF64 psModelFunc_WAUSS(psVector *deriv,
-                        const psVector *params,
-                        const psVector *x)
+XXX broken code
+
+# define PM_MODEL_FUNC       pmModelFunc_WAUSS
+# define PM_MODEL_FLUX       pmModelFlux_WAUSS
+# define PM_MODEL_GUESS      pmModelGuess_WAUSS
+# define PM_MODEL_LIMITS     pmModelLimits_WAUSS
+# define PM_MODEL_RADIUS     pmModelRadius_WAUSS
+# define PM_MODEL_FROM_PSF   pmModelFromPSF_WAUSS
+# define PM_MODEL_FIT_STATUS pmModelFitStatus_WAUSS
+
+psF64 PS_MODEL_FUNC (psVector *deriv,
+                     const psVector *params,
+                     const psVector *x)
 {
     psF32 X = x->data.F32[0] - params->data.F32[2];
@@ -43 +66 @@
 }
 
-// this is probably wrong since it uses the gauss integral 2 pi sigma^2
-psF64 psModelFlux_WAUSS(const psVector *params)
+bool PM_MODEL_LIMITS  (psVector **beta_lim, psVector **params_min, psVector **params_max)
 {
-    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;
+    *beta_lim   = psVectorAlloc (8, PS_TYPE_F32);
+    *params_min = psVectorAlloc (8, PS_TYPE_F32);
+    *params_max = psVectorAlloc (8, PS_TYPE_F32);
 
-    return(Flux);
+    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);
 }
 
-// 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);
-    if (flux >= 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, psSource *source)
+bool PS_MODEL_GUESS  (psModel *model, psSource *source)
 {
 
@@ -90 +121 @@
 }
 
-bool psModelFromPSF_WAUSS (psModel *modelPSF, psModel *modelFLT, pmPSF *psf)
+// this is probably wrong since it uses the gauss integral 2 pi sigma^2
+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);
+}
+
+// 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] <= 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 PS_MODEL_FROM_PSF  (psModel *modelPSF, psModel *modelFLT, pmPSF *psf)
 {
 
@@ -107 +167 @@
     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