Index: /trunk/psphot/src/psphotImageQuality.c
===================================================================
--- /trunk/psphot/src/psphotImageQuality.c	(revision 19896)
+++ /trunk/psphot/src/psphotImageQuality.c	(revision 19896)
@@ -0,0 +1,74 @@
+# include "psphotInternal.h"
+
+// selecting the 'good' stars (likely to be psf stars), measure the M_cn, M_sn terms for n = 2,3,4
+bool psphotImageQuality (psMetadata *recipe, psArray *sources) {
+
+    psVector *M2 = psVectorAllocEmpty (100);
+    psVector *M3 = psVectorAllocEmpty (100);
+    psVector *M4 = psVectorAllocEmpty (100);
+
+    for (int i = 0; i < sources->n; i++) {
+
+	pmSource *source = sources->data[i];
+	if (!source) continue;
+
+	// select by S/N?
+	// select by PSFSTAR mode?
+	// ??
+        if (source->type != PM_SOURCE_TYPE_STAR) continue;
+	if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
+
+	pmMoments *moments = source->moments;
+	if (!moments) continue;
+
+	// if (source->moments->SN < XXX) continue;
+
+	// M_cn = \sum (f r^2 cos(n t)) / \sum (f)
+	// M_sn = \sum (f r^2 sin(n t)) / \sum (f)
+
+	// r^2 cos(2t) = r^2 cos^2 t - r^2 sin^2 t = x^2 - y^2
+	// r^2 sin(2t) = r^2 2 cos t sin t = 2 x y
+
+	// r^2 cos(3t) = r^2 cos^3 t - r^2 3 cos t sin^2 t = (x^3 - 3 x y^2) / r
+	// r^2 sin(3t) = r^2 3 cos^2 t sin t - sin^3 t = (3 x^2 y - y^3) / r
+
+	// r^2 cos(4t) = r^2 cos^4 t - r^2 6 cos^2 t sin^2 t + r^2 sin^4 t = (x^4 - 6 x^2 y^2 + y^4) / r^2
+	// r^2 sin(4t) = r^2 4 cos^3 t sin t - 4 sin^3 t cos t = (4 x^3 y - 4 y^3 x) / r^2
+
+	// the moments->Mnnn terms contain the straight sums: for example: moments->Mxxx contains \sum x^3 / r
+	M_c2 = moments->Mxx - moments->Myy;
+	M_s2 = 2*moments->Mxy;
+	psVectorAppend (M2, sqrt (M_c2^2 + M_s2^2));
+	
+	M_c3 = moments->Mxxx - 3.0*moments->Mxyy;
+	M_s3 = 3.0*moments->Mxxy - moments->Myyy;
+	psVectorAppend (M3, sqrt (M_c3^2 + M_s3^2));
+	
+	M_c4 = moments->Mxxxx - 6.0*moments->Mxxyy + moments->Myyyy;
+	M_s4 = 4.0*moments->Mxxxy - 4.0*moments->Mxyyy;
+	psVectorAppend (M4, sqrt (M_c2^2 + M_s2^2));
+    }	
+
+    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN, PS_STAT_ROBUST_STDEV);
+
+    psVectorStats (stats, M2, NULL, NULL, 0);
+    float vM2 = stats->robustMedian;
+    float dM2 = stats->robustStdev;
+    psMetadataAddF32 (header, PS_LIST_TAIL, "IQ_M2",   PS_META_REPLACE, "M_2 = sqrt (M_c2^2 + M_s2^2)", vM2);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "IQ_D_M2", PS_META_REPLACE, "Stdev of M_2", dM2);
+
+    psVectorStats (stats, M3, NULL, NULL, 0);
+    float vM3 = stats->robustMedian;
+    float dM3 = stats->robustStdev;
+    psMetadataAddF32 (header, PS_LIST_TAIL, "IQ_M3",   PS_META_REPLACE, "M_3 = sqrt (M_c3^2 + M_s3^2)", vM3);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "IQ_D_M3", PS_META_REPLACE, "Stdev of M_3", dM3);
+
+    psVectorStats (stats, M3, NULL, NULL, 0);
+    float vM4 = stats->robustMedian;
+    float dM4 = stats->robustStdev;
+    psMetadataAddF32 (header, PS_LIST_TAIL, "IQ_M4",   PS_META_REPLACE, "M_4 = sqrt (M_c4^2 + M_s4^2)", vM4);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "IQ_D_M4", PS_META_REPLACE, "Stdev of M_4", dM4);
+
+    psLogMsg ("psphot", PS_LOG_INFO, "Image Quality Stats : M_2 : %f +/- %f, M_3 : %f +/- %f, M_4 : %f +/- %f\n", vM2, dM2, vM3, dM3, vM4, dM4);
+    return true;
+}