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

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

@@ -1 @@
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
 /******************************************************************************
-    one component, two slopes:
-    1 / (1 + z^Npow + PAR8*z^4)
+ * this file defines the ZGAUSS 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] = 1 / SigmaX;
-    params->data.F32[5] = 1 / SigmaY;
-    params->data.F32[6] = Sxy;
-    params->data.F32[7] = Npow
-*****************************************************************************/
-
-# define SQ(A)((A)*(A))
+   power-law with fitted slope and tidal cutoff
+   1 / (1 + z^N + (Az)^4)
+ 
+   * 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   - slope of power-law component (N)
+   *****************************************************************************/
+
+/***
+    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_ZGAUSS
+# define PM_MODEL_FLUX       pmModelFlux_ZGAUSS
+# define PM_MODEL_GUESS      pmModelGuess_ZGAUSS
+# define PM_MODEL_LIMITS     pmModelLimits_ZGAUSS
+# define PM_MODEL_RADIUS     pmModelRadius_ZGAUSS
+# define PM_MODEL_FROM_PSF   pmModelFromPSF_ZGAUSS
+# define PM_MODEL_FIT_STATUS pmModelFitStatus_ZGAUSS
+
 # define PAR8 0.1
 
-psF64 psModelFunc_ZGAUSS(psVector *deriv,
-                         const psVector *params,
-                         const psVector *x)
+psF64 PS_MODEL_FUNC (psVector *deriv,
+                     const psVector *params,
+                     const psVector *x)
 {
     psF32 *PAR = params->data.F32;
@@ -54 +71 @@
 }
 
-psF64 psModelFlux_ZGAUSS(const psVector *params)
-{
-    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
-
-    // 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[7]) + SQ(SQ(pr)));
-        norm += f;
-    }
-    norm *= 0.01;
-
-    psF64 Flux = PAR[1] * 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;
-
-    psEllipseAxes axes;
-    psEllipseShape shape;
-
-    if (flux <= 0)
-        return (1.0);
-    if (PAR[1] <= 0)
-        return (1.0);
-    if (flux >= PAR[1])
-        return (1.0);
-
-    // convert Sx,Sy,Sxy to major/minor axes
-    shape.sx = 1.0 / PAR[4];
-    shape.sy = 1.0 / PAR[5];
-    shape.sxy = PAR[6];
-
-    axes = psEllipseShapeToAxes (shape);
-    psF64 dr = 1.0 / axes.major;
-    psF64 limit = flux / PAR[1];
-
-    // 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[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)
+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)
 {
 
@@ -149 +135 @@
 }
 
-bool psModelFromPSF_ZGAUSS (psModel *modelPSF, psModel *modelFLT, pmPSF *psf)
+psF64 PS_MODEL_FLUX (const psVector *params)
+{
+    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
+
+    // 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[7]) + SQ(SQ(pr)));
+        norm += f;
+    }
+    norm *= 0.01;
+
+    psF64 Flux = PAR[1] * 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 PS_MODEL_RADIUS   (const psVector *params, psF64 flux)
+{
+    psF64 r, z, pr, f;
+    psF32 *PAR = params->data.F32;
+
+    psEllipseAxes axes;
+    psEllipseShape shape;
+
+    if (flux <= 0)
+        return (1.0);
+    if (PAR[1] <= 0)
+        return (1.0);
+    if (flux >= PAR[1])
+        return (1.0);
+
+    // convert Sx,Sy,Sxy to major/minor axes
+    shape.sx = 1.0 / PAR[4];
+    shape.sy = 1.0 / PAR[5];
+    shape.sxy = PAR[6];
+
+    axes = psEllipseShapeToAxes (shape);
+    psF64 dr = 1.0 / axes.major;
+    psF64 limit = flux / PAR[1];
+
+    // 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[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 PS_MODEL_FROM_PSF  (psModel *modelPSF, psModel *modelFLT, pmPSF *psf)
 {
 
@@ -166 +220 @@
     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