Index: trunk/psModules/src/objects/pmSourceMatch.c
===================================================================
--- trunk/psModules/src/objects/pmSourceMatch.c	(revision 20948)
+++ trunk/psModules/src/objects/pmSourceMatch.c	(revision 20949)
@@ -282,29 +282,25 @@
 
 
-
-
-
-psVector *pmSourceMatchRelphot(const psArray *matches, // Array of matches
-                               const psVector *zp, // Zero points for each image (including airmass term)
-                               float cloudClip // Clipping for clouds
-                               )
-{
-    PS_ASSERT_ARRAY_NON_NULL(matches, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(zp, NULL);
-    PS_ASSERT_VECTOR_TYPE(zp, PS_TYPE_F32, NULL);
-    PS_ASSERT_FLOAT_LARGER_THAN(cloudClip, 0.0, NULL);
-
+// Iterate on the star magnitudes and image transparencies
+// Returns the solution chi^2
+float sourceMatchRelphotIterate(psVector *trans, // Transparencies
+                                psVector *stars, // Star magnitudes
+                                const psArray *matches, // Array of matches
+                                const psVector *zp // Zero points for each image (including airmass term)
+                                )
+{
     int numImages = zp->n;              // Number of images
     int numStars = matches->n;          // Number of stars
-    psVector *trans = psVectorAlloc(numImages, PS_TYPE_F32); // Transparencies for each image, magnitudes
+
+    psAssert(trans && trans->type.type == PS_TYPE_F32, "Need transparencies");
+    psAssert(trans->n == numImages, "Not enough transparencies: %ld\n", trans->n);
+    psAssert(stars && stars->type.type == PS_TYPE_F32, "Need star magnitudes");
+    psAssert(stars->n == numStars, "Not enough stars: %ld\n", stars->n);
+    psAssert(matches, "Need list of matches");
+    psAssert(zp && zp->type.type == PS_TYPE_F32, "Need zero points");
+    psAssert(zp->n == numImages, "Not enough ZPs: %ld", zp->n);
+
     psVector *transErr = psVectorAlloc(numImages, PS_TYPE_F32); // Errors in the transparency
-    psVector *stars = psVectorAlloc(numStars, PS_TYPE_F32); // Magnitudes for each star
     psVector *starsErr = psVectorAlloc(numStars, PS_TYPE_F32); // Errors in the star magnitudes
-
-
-    // XXX Need to apply provided ZPs
-
-    psVectorInit(trans, 0.0);
-    psVectorInit(transErr, 0.0);
 
     // Solve the star magnitudes
@@ -318,5 +314,5 @@
             float magErr = match->magErr->data.F32[j]; // Error in measured magnitude
             float invErr2 = 1.0 / PS_SQR(magErr); // Inverse square error
-            stars->data.F32[i] += (mag - trans->data.F32[index]) * invErr2;
+            stars->data.F32[i] += (mag + zp->data.F32[index] - trans->data.F32[index]) * invErr2;
             starsErr->data.F32[i] += invErr2;
         }
@@ -347,12 +343,57 @@
         transErr->data.F32[i] = sqrtf(inverse);
 
-        psTrace("psModules.objects", 1, "Transparency for image %d: %f +/- %f\n",
+        psTrace("psModules.objects", 3, "Transparency for image %d: %f +/- %f\n",
                 i, trans->data.F32[i], transErr->data.F32[i]);
     }
 
+    // Once more through to evaluate chi^2
+    float chi2 = 0.0;                   // chi^2 for iteration
+    for (int i = 0; i < numStars; i++) {
+        pmSourceMatch *match = matches->data[i]; // Matched stars
+        for (int j = 0; j < match->num; j++) {
+            int index = match->image->data.U32[j]; // Image index
+            float mag = match->mag->data.F32[j]; // Measured magnitude
+            float magErr = match->magErr->data.F32[j]; // Error in measured magnitude
+            chi2 += PS_SQR(mag - stars->data.F32[index] - trans->data.F32[index] + zp->data.F32[index]) /
+                PS_SQR(magErr);
+        }
+    }
+    chi2 /= numStars + numImages;
+
     psFree(transErr);
-    psFree(stars);
     psFree(starsErr);
 
+    return chi2;
+}
+
+
+psVector *pmSourceMatchRelphot(const psArray *matches, // Array of matches
+                               const psVector *zp, // Zero points for each image (including airmass term)
+                               int maxIter, // Maximum number of iterations
+                               float tol, // Relative tolerance for convergence
+                               float cloudClip, // Clipping for clouds
+                               float starClip // Clipping for stars
+                               )
+{
+    PS_ASSERT_ARRAY_NON_NULL(matches, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(zp, NULL);
+    PS_ASSERT_VECTOR_TYPE(zp, PS_TYPE_F32, NULL);
+    PS_ASSERT_FLOAT_LARGER_THAN(cloudClip, 0.0, NULL);
+
+    int numImages = zp->n;              // Number of images
+    int numStars = matches->n;          // Number of stars
+    psVector *trans = psVectorAlloc(numImages, PS_TYPE_F32); // Transparencies for each image, magnitudes
+    psVectorInit(trans, 0.0);
+    psVector *stars = psVectorAlloc(numStars, PS_TYPE_F32); // Magnitudes for each star
+
+    float chi2 = sourceMatchRelphotIterate(trans, stars, matches, zp); // chi^2 for solution
+    psTrace("psModules.objects", 1, "Initial: chi^2 = %f\n", chi2);
+    float lastChi2 = INFINITY;          // chi^2 on last iteration
+    for (int i = 0; i < maxIter && (lastChi2 - chi2) / chi2 > tol; i++) {
+        lastChi2 = chi2;
+        chi2 = sourceMatchRelphotIterate(trans, stars, matches, zp); // chi^2 for solution
+        psTrace("psModules.objects", 1, "iter = %d: chi^2 = %f\n", i, chi2);
+    }
+
     return trans;
 }