Changeset 6481 for trunk/psphot/src/psSparse.c

Timestamp:

Feb 23, 2006, 6:16:11 PM (20 years ago)

Author:

eugene

Message:

fixed some errors with double sources, added second-pass linear PSF fit, multiple-depths

File:

• trunk/psphot/src/psSparse.c (modified) (4 diffs)

trunk/psphot/src/psSparse.c

@@ -34 +34 @@
     sparse->Bfj->data.F32[2] = B->data.F32[2];
 
-    x = psSparseSolve (x, sparse, 0);
-    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
-
-    x = psSparseSolve (x, sparse, 1);
-    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
-
-    x = psSparseSolve (x, sparse, 2);
-    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
-
-    x = psSparseSolve (x, sparse, 3);
-    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
-
-    x = psSparseSolve (x, sparse, 4);
+    psSparseConstraint constraint;
+    constraint.paramMin   = -1e8;
+    constraint.paramMax   = +1e8;
+    constraint.paramDelta = +1e8;
+
+    x = psSparseSolve (x, constraint, sparse, 0);
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    x = psSparseSolve (x, constraint, sparse, 1);
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    x = psSparseSolve (x, constraint, sparse, 2);
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    x = psSparseSolve (x, constraint, sparse, 3);
+    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
+
+    x = psSparseSolve (x, constraint, sparse, 4);
     fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
     return;
@@ -166 +171 @@
 }
 
-psVector *psSparseSolve (psVector *guess, psSparse *sparse, int Niter) {
+psVector *psSparseSolve (psVector *guess, psSparseConstraint constraint, psSparse *sparse, int Niter) {
 
     psF32 dG;
@@ -182 +187 @@
         for (int i = 0; i < dQ->n; i++) {
             dG = (dQ->data.F32[i] - Bfj->data.F32[i]) / Qii->data.F32[i];
+            if (fabs (dG) > constraint.paramDelta) {
+                if (dG > 0) {
+                    dG = +constraint.paramDelta;
+                } else {
+                    dG = -constraint.paramDelta;
+                }
+            }
             guess->data.F32[i] -= dG;
+            guess->data.F32[i] = PS_MAX (guess->data.F32[i], constraint.paramMin);
+            guess->data.F32[i] = PS_MIN (guess->data.F32[i], constraint.paramMax);
         }
     }
@@ -220 +234 @@
 }
 
-