Index: /trunk/Ohana/src/libdvo/src/coordops.c
===================================================================
--- /trunk/Ohana/src/libdvo/src/coordops.c	(revision 8293)
+++ /trunk/Ohana/src/libdvo/src/coordops.c	(revision 8294)
@@ -41,13 +41,11 @@
 
   /** extra polynomial terms **/
-  if (Polynomial) {
-    if (coords[0].Npolyterms > 1) {
-      L += X*X*coords[0].polyterms[0][0] + X*Y*coords[0].polyterms[1][0] + Y*Y*coords[0].polyterms[2][0];
-      M += X*X*coords[0].polyterms[0][1] + X*Y*coords[0].polyterms[1][1] + Y*Y*coords[0].polyterms[2][1];
-    }
-    if (coords[0].Npolyterms > 2) {
-      L += X*X*X*coords[0].polyterms[3][0] + X*X*Y*coords[0].polyterms[4][0] + X*Y*Y*coords[0].polyterms[5][0] + Y*Y*Y*coords[0].polyterms[6][0];
-      M += X*X*X*coords[0].polyterms[3][1] + X*X*Y*coords[0].polyterms[4][1] + X*Y*Y*coords[0].polyterms[5][1] + Y*Y*Y*coords[0].polyterms[6][1];
-    }
+  if (coords[0].Npolyterms > 1) {
+    L += X*X*coords[0].polyterms[0][0] + X*Y*coords[0].polyterms[1][0] + Y*Y*coords[0].polyterms[2][0];
+    M += X*X*coords[0].polyterms[0][1] + X*Y*coords[0].polyterms[1][1] + Y*Y*coords[0].polyterms[2][1];
+  }
+  if (coords[0].Npolyterms > 2) {
+    L += X*X*X*coords[0].polyterms[3][0] + X*X*Y*coords[0].polyterms[4][0] + X*Y*Y*coords[0].polyterms[5][0] + Y*Y*Y*coords[0].polyterms[6][0];
+    M += X*X*X*coords[0].polyterms[3][1] + X*X*Y*coords[0].polyterms[4][1] + X*Y*Y*coords[0].polyterms[5][1] + Y*Y*Y*coords[0].polyterms[6][1];
   }
 
@@ -253,5 +251,5 @@
 
   /** extra polynomial terms **/
-  if (Polynomial) {
+  if (coords[0].Npolyterms > 1) {
     for (i = 0; i < 10; i++) {
       Lo = (X*coords[0].pc1_1 + Y*coords[0].pc1_2);
@@ -310,5 +308,5 @@
 int GetCoords (Coords *coords, Header *header) {
   
-  int status, itmp, Polynomial;
+  int status, itmp, Polynomial, Polyterm;
   double Lambda, rotate, scale;
   double equinox;
@@ -345,34 +343,35 @@
       }
 
+      /* set NPLYTERM based on header.  if NPLYTERM is missing, it should have a 
+	 value of 0, unless the projection type is one of PLY, DIS, WRP, in which
+	 case it should be set to 3 */
       ctype = &coords[0].ctype[4];
       Polynomial = !strcmp (ctype, "-PLY") || !strcmp (ctype, "-DIS") || !strcmp (ctype, "-WRP");
-
-      if (Polynomial) {
-	if (gfits_scan (header, "NPLYTERM", "%d", 1, &itmp)) {
-	  coords[0].Npolyterms = itmp;
-	} else {
-	  coords[0].Npolyterms = 3;
-	}
-	switch (coords[0].Npolyterms) {
-	  case 3:
-	    status &= gfits_scan (header, "PCA1X3Y0", "%f", 1, &coords[0].polyterms[3][0]);
-	    status &= gfits_scan (header, "PCA1X2Y1", "%f", 1, &coords[0].polyterms[4][0]);
-	    status &= gfits_scan (header, "PCA1X1Y2", "%f", 1, &coords[0].polyterms[5][0]);
-	    status &= gfits_scan (header, "PCA1X0Y3", "%f", 1, &coords[0].polyterms[6][0]);
-	    status &= gfits_scan (header, "PCA2X3Y0", "%f", 1, &coords[0].polyterms[3][1]);
-	    status &= gfits_scan (header, "PCA2X2Y1", "%f", 1, &coords[0].polyterms[4][1]);
-	    status &= gfits_scan (header, "PCA2X1Y2", "%f", 1, &coords[0].polyterms[5][1]);
-	    status &= gfits_scan (header, "PCA2X0Y3", "%f", 1, &coords[0].polyterms[6][1]);
-	  case 2:
-	    status &= gfits_scan (header, "PCA1X2Y0", "%f", 1, &coords[0].polyterms[0][0]);
-	    status &= gfits_scan (header, "PCA1X1Y1", "%f", 1, &coords[0].polyterms[1][0]);
-	    status &= gfits_scan (header, "PCA1X0Y2", "%f", 1, &coords[0].polyterms[2][0]);
-	    status &= gfits_scan (header, "PCA2X2Y0", "%f", 1, &coords[0].polyterms[0][1]);
-	    status &= gfits_scan (header, "PCA2X1Y1", "%f", 1, &coords[0].polyterms[1][1]);
-	    status &= gfits_scan (header, "PCA2X0Y2", "%f", 1, &coords[0].polyterms[2][1]);
-	  case 0:
-	  case 1:
-	    break;
-	}
+      Polyterm = gfits_scan (header, "NPLYTERM", "%d", 1, &itmp);
+
+      coords[0].Npolyterms = 0;
+      if (Polynomial && !Polyterm) coords[0].Npolyterms = 3;
+      if (Polyterm) coords[0].Npolyterms = itmp;
+
+      switch (coords[0].Npolyterms) {
+	case 3:
+	  status &= gfits_scan (header, "PCA1X3Y0", "%f", 1, &coords[0].polyterms[3][0]);
+	  status &= gfits_scan (header, "PCA1X2Y1", "%f", 1, &coords[0].polyterms[4][0]);
+	  status &= gfits_scan (header, "PCA1X1Y2", "%f", 1, &coords[0].polyterms[5][0]);
+	  status &= gfits_scan (header, "PCA1X0Y3", "%f", 1, &coords[0].polyterms[6][0]);
+	  status &= gfits_scan (header, "PCA2X3Y0", "%f", 1, &coords[0].polyterms[3][1]);
+	  status &= gfits_scan (header, "PCA2X2Y1", "%f", 1, &coords[0].polyterms[4][1]);
+	  status &= gfits_scan (header, "PCA2X1Y2", "%f", 1, &coords[0].polyterms[5][1]);
+	  status &= gfits_scan (header, "PCA2X0Y3", "%f", 1, &coords[0].polyterms[6][1]);
+	case 2:
+	  status &= gfits_scan (header, "PCA1X2Y0", "%f", 1, &coords[0].polyterms[0][0]);
+	  status &= gfits_scan (header, "PCA1X1Y1", "%f", 1, &coords[0].polyterms[1][0]);
+	  status &= gfits_scan (header, "PCA1X0Y2", "%f", 1, &coords[0].polyterms[2][0]);
+	  status &= gfits_scan (header, "PCA2X2Y0", "%f", 1, &coords[0].polyterms[0][1]);
+	  status &= gfits_scan (header, "PCA2X1Y1", "%f", 1, &coords[0].polyterms[1][1]);
+	  status &= gfits_scan (header, "PCA2X0Y2", "%f", 1, &coords[0].polyterms[2][1]);
+	case 0:
+	case 1:
+	  break;
       }
     } else {
@@ -574,2 +573,23 @@
 
 # endif
+
+/*
+
+Projections and Transformations
+
+The Coords structure is used to represent the standard FITS header representations of the WCS terms.
+The FITS WCS encapsulates several concepts in one system: coordinate frame, projection, and
+transformations.  The ctype variable defines both the frame (ie, RA/DEC, GLON/GLAT, etc) and the
+projection (ie, AIT, SIN, TAN, etc).  The associated terms (PC00i00j) define transformations from
+the projection system to another linear coordinate system.  I have extended the basic WCS
+transformation terms to include higher-order polynomial terms.  The presence of higher-order terms
+is indicated in the header by the NPLYTERM keyword, with a value greater than 1 (a value of 0 or 1
+implies no higher-order coefficients).  The coefficients have keywords of the form PCAnXiYj where
+'n' is the axis corresponding to the dependent variable, 'i' is the order of the X component and 'j'
+is the order of the Y component.  A value of 2 indicates that the second-order coefficients are
+defined (PCAnX2Y0, PCAnX1Y1, PCAnX0Y2); A value of 3 indicates that the third-order coefficients are
+also defined (PCAnX3Y0, PCAnX2Y1, PCAnX1Y2, PCAnX0Y3).  Some headers in the past were generated
+without the NPLYTERM keyword.  a value of PLY, DIS, or WRP for the projection without a
+corresponding value for NPLYTERM implies that the value should be interpreted as '3'.
+
+*/