Changeset 19898

Timestamp:

Oct 4, 2008, 6:14:24 AM (18 years ago)

Author:

eugene

Message:

third and fourth moments are f r^{2 cos(nt) (divide by r and r}2

File:

• trunk/psModules/src/objects/pmSourceMoments.c (modified) (5 diffs)

trunk/psModules/src/objects/pmSourceMoments.c

@@ -6 +6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-10-04 02:29:19 $
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-10-04 16:14:24 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -204 +204 @@
 
             // radius is just a function of (xDiff, yDiff)
-            if (!VALID_RADIUS(xDiff, yDiff, R2)) continue;
+            psF32 r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+            psF32 r  = sqrt(r2);
+            if (r > radius) continue;
 
             psF32 pDiff = *vPix - sky;
@@ -210 +212 @@
 
             // XXX EAM : should this limit be user-defined?
+
             if (PS_SQR(pDiff) < wDiff) continue;
             if (pDiff < 0) continue;
@@ -222 +225 @@
             psF32 yy = yDiff * y;
 
-            psF32 xxx = xDiff * xx;
-            psF32 xxy = xDiff * xy;
-            psF32 xyy = xDiff * yy;
-            psF32 yyy = yDiff * yy;
-
-            psF32 xxxx = xDiff * xxx;
-            psF32 xxxy = xDiff * xxy;
-            psF32 xxyy = xDiff * xyy;
-            psF32 xyyy = xDiff * yyy;
-            psF32 yyyy = yDiff * yyy;
+            psF32 xxx = xDiff * xx / r;
+            psF32 xxy = xDiff * xy / r;
+            psF32 xyy = xDiff * yy / r;
+            psF32 yyy = yDiff * yy / r;
+
+            psF32 xxxx = xDiff * xxx / r2;
+            psF32 xxxy = xDiff * xxy / r2;
+            psF32 xxyy = xDiff * xyy / r2;
+            psF32 xyyy = xDiff * yyy / r2;
+            psF32 yyyy = yDiff * yyy / r2;
 
             XX  += xx;
@@ -264 +267 @@
     source->moments->Mxyyy = XYYY/Sum;
     source->moments->Myyyy = YYYY/Sum;
+
+    if (source->moments->Mxx < 0) {
+      fprintf (stderr, "error: neg second moment??\n");
+    }
+    if (source->moments->Myy < 0) {
+      fprintf (stderr, "error: neg second moment??\n");
+    }
 
     psTrace ("psModules.objects", 4, "Mxx: %f  Mxy: %f  Myy: %f  Mxxx: %f  Mxxy: %f  Mxyy: %f  Myyy: %f  Mxxxx: %f  Mxxxy: %f  Mxxyy: %f  Mxyyy: %f  Mxyyy: %f\n",