Changeset 8882 for trunk/psphot/src/models/pmModel_SGAUSS.c

Timestamp:

Sep 22, 2006, 2:29:31 AM (20 years ago)

Author:

rhl

Message:

[Start to] provide symbolic names for parameters

File:

• trunk/psphot/src/models/pmModel_SGAUSS.c (modified) (5 diffs)

trunk/psphot/src/models/pmModel_SGAUSS.c

@@ -25 +25 @@
     psF32 *PAR = params->data.F32;
 
-    psF32 X  = x->data.F32[0] - PAR[2];
-    psF32 Y  = x->data.F32[1] - PAR[3];
-    psF32 px = PAR[4]*X;
-    psF32 py = PAR[5]*Y;
-    psF32 z  = PS_MAX((0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[6]*X*Y), 1e-8);
-    // note that if z -> 0, dPAR[7] -> -inf
-    // also z^(PAR[7]-1) -> Inf 
-
-    psF32 pr = z*PAR[8];
+    psF32 X  = x->data.F32[0] - PAR[PM_PAR_XPOS];
+    psF32 Y  = x->data.F32[1] - PAR[PM_PAR_YPOS];
+    psF32 px = PAR[PM_PAR_SXX]*X;
+    psF32 py = PAR[PM_PAR_SYY]*Y;
+    psF32 z  = PS_MAX((0.5*PS_SQR(px) + 0.5*PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y), 1e-8);
+    // note that if z -> 0, dPAR[PM_PAR_7] -> -inf
+    // also z^(PAR[PM_PAR_7]-1) -> Inf 
+
+    psF32 pr = z*PAR[PM_PAR_8];
     psF32 pr3 = pr*pr*pr;
-    psF32 p  = pow(z, PAR[7] - 1.0);
+    psF32 p  = pow(z, PAR[PM_PAR_7] - 1.0);
     psF32 r  = 1.0 / (1 + z*p + pr*pr3);
-    psF32 f  = PAR[1]*r + PAR[0];
+    psF32 f  = PAR[PM_PAR_FLUX]*r + PAR[PM_PAR_SKY];
 
     if (deriv != NULL) {
         // note difference from a pure gaussian: q = params->data.F32[1]*r
-        psF32 t = PAR[1]*r*r;
-        psF32 q = t*(PAR[7]*p + 4*PAR[8]*pr3);
+        psF32 t = PAR[PM_PAR_FLUX]*r*r;
+        psF32 q = t*(PAR[PM_PAR_7]*p + 4*PAR[PM_PAR_8]*pr3);
 
         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[2] = q*(2.0*px*PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]*Y);
+        deriv->data.F32[3] = q*(2.0*py*PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]*X);
         deriv->data.F32[4] = -2.0*q*px*X;
         deriv->data.F32[5] = -2.0*q*py*Y;
@@ -230 +230 @@
     psF32 *PAR = params->data.F32;
 
-    psF64 A1   = PS_SQR(PAR[4]);
-    psF64 A2   = PS_SQR(PAR[5]);
-    psF64 A3   = PS_SQR(PAR[6]);
+    psF64 A1   = PS_SQR(PAR[PM_PAR_SXX]);
+    psF64 A2   = PS_SQR(PAR[PM_PAR_SYY]);
+    psF64 A3   = PS_SQR(PAR[PM_PAR_SXY]);
     psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3);
     // Area is equivalent to 2 pi sigma^2
@@ -239 +239 @@
     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)));
+      psF32 pr = PAR[PM_PAR_8]*z;
+        f = 1.0 / (1 + pow(z, PAR[PM_PAR_7]) + SQ(SQ(pr)));
         norm += f;
     }
     norm *= 0.01;
     
-    psF64 Flux = PAR[1] * Area * norm;
+    psF64 Flux = PAR[PM_PAR_FLUX] * Area * norm;
 
     return(Flux);
@@ -262 +262 @@
 
     if (flux <= 0) return (1.0);
-    if (PAR[1] <= 0) return (1.0);
-    if (flux >= PAR[1]) return (1.0);
+    if (PAR[PM_PAR_FLUX] <= 0) return (1.0);
+    if (flux >= PAR[PM_PAR_FLUX]) 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];
+    shape.sx = 1.0 / PAR[PM_PAR_SXX];
+    shape.sy = 1.0 / PAR[PM_PAR_SYY];
+    shape.sxy = PAR[PM_PAR_SXY];
 
     axes = EllipseShapeToAxes (shape);
     psF64 dr = 1.0 / axes.major;
-    psF64 limit = flux / PAR[1];
+    psF64 limit = flux / PAR[PM_PAR_FLUX];
 
     // 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)));
+        pr = PAR[PM_PAR_8]*z;
+        f = 1.0 / (1 + pow(z, PAR[PM_PAR_7]) + SQ(SQ(pr)));
         if (f < limit) break;
     }
@@ -315 +315 @@
     psF32 *dPAR = model->dparams->data.F32;
 
-    shape.sx = 1.0 / PAR[4];
-    shape.sy = 1.0 / PAR[5];
-    shape.sxy = PAR[6];
+    shape.sx = 1.0 / PAR[PM_PAR_SXX];
+    shape.sy = 1.0 / PAR[PM_PAR_SYY];
+    shape.sxy = PAR[PM_PAR_SXY];
 
     axes = EllipseShapeToAxes (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 += PS_SQR(dPAR[PM_PAR_SXX] / PAR[PM_PAR_SXX]);
+    dP += PS_SQR(dPAR[PM_PAR_SYY] / PAR[PM_PAR_SYY]);
+    dP += PS_SQR(dPAR[PM_PAR_7] / PAR[PM_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 &= (PAR[PM_PAR_FLUX] > 0);
+    status &= ((dPAR[PM_PAR_FLUX]/PAR[PM_PAR_FLUX]) < 0.5);
+    status &= (fabs(PAR[PM_PAR_8]) < 0.5);
+    status &= (dPAR[PM_PAR_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);
+    status &= (PAR[PM_PAR_7] > 0.5);
 
     if (status) return true;