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

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

@@ -1 @@
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-
 /******************************************************************************
-    one component, two slopes:
-    1 / (1 + z^Npow + St*z^Ntide)
+ * this file defines the SGAUSS 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
-    params->data.F32[8] = St
-*****************************************************************************/
-
-# define SQ(A)((A)*(A))
+   power-law with fitted slope and outer tidal radius
+   1 / (1 + z^N + kz^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)
+   * PM_PAR_8   8   - amplitude of the tidal cutoff (k)
+   *****************************************************************************/
+
+/***
+    XXXX the model in this file needs to be tested more carefully.
+    the code for guessing the power-law slope based on the radial profile 
+    is either too slow or does not work well.
+    fix up the code to follow conventions in the other model function files.
+***/
+
+XXX broken code
+
+# define PM_MODEL_FUNC       pmModelFunc_SGAUSS
+# define PM_MODEL_FLUX       pmModelFlux_SGAUSS
+# define PM_MODEL_GUESS      pmModelGuess_SGAUSS
+# define PM_MODEL_LIMITS     pmModelLimits_SGAUSS
+# define PM_MODEL_RADIUS     pmModelRadius_SGAUSS
+# define PM_MODEL_FROM_PSF   pmModelFromPSF_SGAUSS
+# define PM_MODEL_FIT_STATUS pmModelFitStatus_SGAUSS
+
 psF64 psImageEllipseContour (psEllipseAxes 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 PM_MODEL_FUNC (psVector *deriv,
+                     const psVector *params,
+                     const psVector *x)
 {
     psF32 *PAR = params->data.F32;
@@ -61 +78 @@
 }
 
-bool pmModelLimits_SGAUSS (psVector **beta_lim, psVector **params_min, psVector **params_max)
+bool PM_MODEL_LIMITS  (psVector **beta_lim, psVector **params_min, psVector **params_max)
 {
 
@@ -99 +116 @@
 
     return (TRUE);
+}
+
+bool PM_MODEL_GUESS  (pmModel *model, pmSource *source)
+{
+
+    pmMoments *sMoments = source->moments;
+    pmPeak    *peak     = source->peak;
+    psF32     *params   = model->params->data.F32;
+
+    psEllipseAxes axes;
+    psEllipseShape shape;
+    psEllipseMoments 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 = psEllipseMomentsToAxes(moments);
+    shape = psEllipseAxesToShape(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);
+}
+
+psF64 PM_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;
+    for (z = 0.005; z < 50; z += 0.01) {
+        psF32 pr = PAR[8]*z;
+        f = 1.0 / (1 + pow(z, PAR[7]) + PS_SQR(PS_SQR(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 PM_MODEL_RADIUS   (const psVector *params, psF64 flux)
+{
+    psF64 r, z = 0.0, 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 = PS_SQR(r);
+        pr = PAR[8]*z;
+        f = 1.0 / (1 + pow(z, PAR[7]) + PS_SQR(PS_SQR(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 PM_MODEL_FROM_PSF  (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];
+        // XXX: Verify this (from EAM change)
+        //out[i] = Polynomial2DEval_EAM(poly, out[2], out[3]);
+        out[i] = psPolynomial2DEval(poly, out[2], out[3]);
+    }
+    return(true);
+}
+
+bool PM_MODEL_FIT_STATUS  (pmModel *model)
+{
+
+    psF32 dP;
+    bool  status;
+    psEllipseAxes axes;
+    psEllipseShape shape;
+
+    psF32 *PAR  = model->params->data.F32;
+    psF32 *dPAR = model->dparams->data.F32;
+
+    shape.sx = 1.0 / PAR[4];
+    shape.sy = 1.0 / PAR[5];
+    shape.sxy = PAR[6];
+
+    axes = psEllipseShapeToAxes (shape);
+
+    dP = 0;
+    dP += PS_SQR(dPAR[4] / PAR[4]);
+    dP += PS_SQR(dPAR[5] / PAR[5]);
+    dP += PS_SQR(dPAR[7] / PAR[7]);
+    dP = sqrt (dP);
+
+    status = true;
+    status &= (dP < 0.5);
+    status &= (PAR[1] > 0);
+    status &= ((dPAR[1]/PAR[1]) < 0.5);
+    status &= (fabs(PAR[8]) < 0.5);
+    status &= (dPAR[8] < 0.1);
+    status &= (axes.major > 1.41);
+    status &= (axes.minor > 1.41);
+    status &= ((axes.major / axes.minor) < 5.0);
+    status &= (PAR[7] > 0.5);
+
+    if (status)
+        return true;
+    return false;
 }
 
@@ -142 +334 @@
 }
 
-bool pmModelGuess_SGAUSS (pmModel *model, pmSource *source)
+// XXX EAM : test version using flux contours to guess slope
+bool PM_MODEL_GUESS_HARD (pmModel *model, pmSource *source)
 {
 
@@ -148 +341 @@
     pmPeak    *peak     = source->peak;
     psF32     *params   = model->params->data.F32;
+    float f1, f2;
 
     psEllipseAxes axes;
@@ -169 +363 @@
     params[6] = shape.sxy;
 
-    # if (0)
-
-        f1 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels);
+    f1 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels);
     axes.major *= 2.0;
     axes.minor *= 2.0;
@@ -181 +373 @@
         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;
-
-    psEllipseAxes axes;
-    psEllipseShape shape;
-    psEllipseMoments 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 = psEllipseMomentsToAxes(moments);
-    shape = psEllipseAxesToShape(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[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;
-    for (z = 0.005; z < 50; z += 0.01) {
-        psF32 pr = PAR[8]*z;
-        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 pmModelRadius_SGAUSS  (const psVector *params, psF64 flux)
-{
-    psF64 r, z = 0.0, 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 = PAR[8]*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 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];
-        // XXX: Verify this (from EAM change)
-        //out[i] = Polynomial2DEval_EAM(poly, out[2], out[3]);
-        out[i] = psPolynomial2DEval(poly, out[2], out[3]);
-    }
-    return(true);
-}
-
-bool pmModelFitStatus_SGAUSS (pmModel *model)
-{
-
-    psF32 dP;
-    bool  status;
-    psEllipseAxes axes;
-    psEllipseShape shape;
-
-    psF32 *PAR  = model->params->data.F32;
-    psF32 *dPAR = model->dparams->data.F32;
-
-    shape.sx = 1.0 / PAR[4];
-    shape.sy = 1.0 / PAR[5];
-    shape.sxy = PAR[6];
-
-    axes = psEllipseShapeToAxes (shape);
-
-    dP = 0;
-    dP += PS_SQR(dPAR[4] / PAR[4]);
-    dP += PS_SQR(dPAR[5] / PAR[5]);
-    dP += PS_SQR(dPAR[7] / PAR[7]);
-    dP = sqrt (dP);
-
-    status = true;
-    status &= (dP < 0.5);
-    status &= (PAR[1] > 0);
-    status &= ((dPAR[1]/PAR[1]) < 0.5);
-    status &= (fabs(PAR[8]) < 0.5);
-    status &= (dPAR[8] < 0.1);
-    status &= (axes.major > 1.41);
-    status &= (axes.minor > 1.41);
-    status &= ((axes.major / axes.minor) < 5.0);
-    status &= (PAR[7] > 0.5);
-
-    if (status)
-        return true;
-    return false;
-}
+# 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