Index: trunk/stac/src/stacRead.c
===================================================================
--- trunk/stac/src/stacRead.c	(revision 2764)
+++ trunk/stac/src/stacRead.c	(revision 2897)
@@ -3,6 +3,4 @@
 #include "pslib.h"
 #include "stac.h"
-
-#define ORDER 2
 
 psArray *stacReadImages(stacConfig *config)
@@ -54,57 +52,75 @@
 	// Read the file
 	psTrace("stac.read.maps",5,"Reading map file %s....\n",mapfile);
-	psPlaneTransform *map = psPlaneTransformAlloc(ORDER,ORDER);
 
-	// Format is as following:
-	// A0  A1
-	// B00 B01
-	// B10 B11
-	//
-	// where:
-	// x_sky = A0 + B00*x + B10*y
-	// y_sky = A1 + B01*x + B11*y
+	// Format is now:
+	// order
+	// x coefficients
+	// y coefficients
+	// 
+	// where the order is 1 for linear, 2 for quadratic, 3 for cubic.
+	// and the coefficients are read by the following pseudo-code:
+	// 
+	// 	for (int k = 0; k < order + 1; k++)
+	// 	    for (int j = 0; j < k; j++)
+	// 		int i = k - j;
+	//              read coefficient of x^i y^j
+	// 
+	// This produces the ordering:
+	// x^0 y^0, x^1 y^0, x^0 y^1, x^2 y^0, x^1 y^1, x^0 y^2, x^3 y^0, x^2 y^1, x^1 y^2, x^0 y^3
+	// 
+	// This is, of course, for third order polynomials.
+	// For lower orders, the list is truncated at the appropriate level.
 
-	double a0, a1, b00, b10, b01, b11;
-
-	if (fscanf(mapfp,"%lf %lf %lf %lf %lf %lf", &a0, &a1, &b00, &b01, &b10, &b11) != 6) {
-	    psLogMsg("stac.read.maps", PS_LOG_ERROR, "Couldn't read all 6 transformation values from %s\n",
-		     mapfile);
+	int order = 0;			// Polynomial order
+	if (fscanf(mapfp, "%d", &order) != 1) {
+	    psLogMsg("stac.read.maps", PS_LOG_ERROR, "Unable to read map file %s\n", mapfile);
 	    exit(EXIT_FAILURE);
 	}
-	fclose(mapfp);
-	map->x->coeff[0][0] = a0;
-	map->x->coeff[1][0] = b00;
-	map->x->coeff[0][1] = b10;
-	map->y->coeff[0][0] = a1;
-	map->y->coeff[1][0] = b01;
-	map->y->coeff[0][1] = b11;
 
-	// Mask out the cross-terms
-	map->x->mask[1][1] = 1;
-	map->y->mask[1][1] = 1;
-	// Set the cross-terms to zero, just in case...
-	map->x->coeff[1][1] = 0.0;
-	map->y->coeff[1][1] = 0.0;
+	psTrace("stac.read.maps", 5, "Polynomial order: %d\n", order);
+
+	psPlaneTransform *map = psPlaneTransformAlloc(order + 1, order + 1); // The transformation
+	// Set coefficient masks
+	for (int i = 0; i < order + 1; i++) {
+	    for (int j = 0; j < order + 1; j++) {
+		if (i + j > order + 1) {
+		    map->x->mask[i][j] = 1;
+		    map->y->mask[i][j] = 1;
+		} else {
+		    map->x->mask[i][j] = 0;
+		    map->y->mask[i][j] = 0;
+		}
+	    }
+	}
+
+	// Read x coefficients
+	psTrace("stac.read.maps", 7, "x' = \n");
+	for (int k = 0; k < order + 1; k++) {
+	    for (int j = 0; j < k + 1; j++) {
+		int i = k - j;
+		if (fscanf(mapfp, "%lf", &map->x->coeff[i][j]) != 1) {
+		    psLogMsg("stac.read.maps", PS_LOG_ERROR, "Unable to read map file %s\n", mapfile);
+		    exit(EXIT_FAILURE);
+		}
+		psTrace("stac.read.maps", 7, "      %f x^%d y^%d\n", map->x->coeff[i][j], i, j);
+	    }
+	}
+	// Read y coefficients
+	psTrace("stac.read.maps", 7, "y' = \n");
+	for (int k = 0; k < order + 1; k++) {
+	    for (int j = 0; j < k + 1; j++) {
+		int i = k - j;
+		if (fscanf(mapfp, "%lf", &map->y->coeff[i][j]) != 1) {
+		    psLogMsg("stac.read.maps", PS_LOG_ERROR, "Unable to read map file %s\n", mapfile);
+		    exit(EXIT_FAILURE);
+		}
+		psTrace("stac.read.maps", 7, "      %f x^%d y^%d\n", map->y->coeff[i][j], i, j);
+	    }
+	}
+
+   	fclose(mapfp);
 
 	// Plug it into the array
 	maps->data[i] = map;
-
-#ifdef TESTING
-	// Check psLib's answer is the same as I would get with old version
-	double detx = 256.0, dety = 256.0;
-	// My code
-	double myx = a0 + b00*detx + b10*dety;
-	double myy = a1 + b01*detx + b11*dety;
-	// psLib
-	psPlane *det = psAlloc(sizeof(psPlane));
-	det->x = 256.0;
-	det->y = 256.0;
-	psPlane *pslib = psPlaneTransformApply(NULL, map, det);
-
-	psTrace("stac.read.maps.test",0,"me: %lf,%lf\n",myx,myy);
-	psTrace("stac.read.maps.test",0,"ps: %lf,%lf\n",pslib->x,pslib->y);
-	psFree(det);
-	psFree(pslib);
-#endif
 
     }