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

Timestamp:

Oct 24, 2006, 2:00:17 PM (20 years ago)

Author:

Paul Price

Message:

Changed definition of psVectorAlloc.

File:

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

trunk/psphot/src/models/pmModel_SGAUSS.c

@@ -20 +20 @@
 
 psF32 pmModelFunc_SGAUSS(psVector *deriv,
-                         const psVector *params,
-                         const psVector *x)
+                         const psVector *params,
+                         const psVector *x)
 {
     psF32 *PAR = params->data.F32;
@@ -31 +31 @@
     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 
+    // also z^(PAR[PM_PAR_7]-1) -> Inf
 
     psF32 pr = z*PAR[PM_PAR_8];
@@ -42 +42 @@
         // note difference from a pure gaussian: q = params->data.F32[1]*r
         psF32 t = PAR[PM_PAR_I0]*r*r;
-        psF32 q = t*(PAR[PM_PAR_7]*p + 4*PAR[PM_PAR_8]*pr3);
+        psF32 q = t*(PAR[PM_PAR_7]*p + 4*PAR[PM_PAR_8]*pr3);
 
         deriv->data.F32[0] = +1.0;
@@ -111 +111 @@
     xo = xc - image->col0;
     yo = yc - image->row0;
-    
-    psVector *contour = psVectorAlloc (Nt, PS_TYPE_F32);
+
+    psVector *contour = psVectorAllocEmpty (Nt, PS_TYPE_F32);
     for (t = 0, N = 0; (t < 2*M_PI) && (N < Nt); t += dt) {
-        x = ct*axes.major*cos(t) + st*axes.minor*sin(t) + xo;
-        y = ct*axes.minor*sin(t) + st*axes.major*cos(t) + yo;
-        value = p_psImageGetElementF64(image, x, y);
-        if (isfinite(value)) {
-            contour->data.F32[N] = value;
-            N++;
-        }
+        x = ct*axes.major*cos(t) + st*axes.minor*sin(t) + xo;
+        y = ct*axes.minor*sin(t) + st*axes.major*cos(t) + yo;
+        value = p_psImageGetElementF64(image, x, y);
+        if (isfinite(value)) {
+            contour->data.F32[N] = value;
+            N++;
+        }
     }
     contour->n = N;
@@ -134 +134 @@
     return (value);
 }
-    
+
 bool pmModelGuess_SGAUSS (pmModel *model, pmSource *source) {
 
@@ -144 +144 @@
     EllipseShape shape;
     EllipseMoments moments;
-    
+
     moments.x2 = PS_SQR(sMoments->Sx);
     moments.y2 = PS_SQR(sMoments->Sy);
@@ -161 +161 @@
     params[6] = shape.sxy;
 
-# if (0)    
+# if (0)
     f1 = psImageEllipseContour (axes, peak->x, peak->y, source->pixels);
     axes.major *= 2.0;
@@ -168 +168 @@
 
     if (f1 > f2) {
-        params[7] = PS_MIN (3.0, PS_MAX (0.5, log(2.0*(f1/f2) - 1.0) / log(2.0)));
+        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[7] = 0.5;
     }
 # endif
@@ -192 +192 @@
     EllipseShape shape;
     EllipseMoments moments;
-    
+
     moments.x2 = PS_SQR(sMoments->Sx);
     moments.y2 = PS_SQR(sMoments->Sy);
@@ -215 +215 @@
 
     if (f1 > f2) {
-        params[7] = PS_MIN (3.0, PS_MAX (0.5, log(2.0*(f1/f2) - 1.0) / log(2.0)));
+        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[7] = 0.5;
     }
     params[8] = 0.1;
@@ -240 +240 @@
     for (z = 0.005; z < 50; z += 0.01) {
       psF32 pr = PAR[PM_PAR_8]*z;
-        f = 1.0 / (1 + pow(z, PAR[PM_PAR_7]) + SQ(SQ(pr)));
-        norm += f;
+        f = 1.0 / (1 + pow(z, PAR[PM_PAR_7]) + SQ(SQ(pr)));
+        norm += f;
     }
     norm *= 0.01;
-    
+
     psF64 Flux = PAR[PM_PAR_I0] * Area * norm;
 
@@ -276 +276 @@
     // 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[PM_PAR_8]*z;
-        f = 1.0 / (1 + pow(z, PAR[PM_PAR_7]) + SQ(SQ(pr)));
-        if (f < limit) break;
+        z = SQ(r);
+        pr = PAR[PM_PAR_8]*z;
+        f = 1.0 / (1 + pow(z, PAR[PM_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");
+        fprintf (stderr, "error in code\n");
     }
     return (radius);
@@ -292 +292 @@
     psF32 *out = modelPSF->params->data.F32;
     psF32 *in  = modelFLT->params->data.F32;
-    
+
     out[0] = in[0];
     out[1] = in[1];
@@ -299 +299 @@
 
     for (int i = 4; i < 9; i++) {
-        psPolynomial2D *poly = psf->params->data[i-4];
-        out[i] = psPolynomial2DEval (poly, out[2], out[3]);
+        psPolynomial2D *poly = psf->params->data[i-4];
+        out[i] = psPolynomial2DEval (poly, out[2], out[3]);
     }
     return(true);
@@ -311 +311 @@
     EllipseAxes axes;
     EllipseShape shape;
-  
+
     psF32 *PAR  = model->params->data.F32;
     psF32 *dPAR = model->dparams->data.F32;