Index: /branches/czw_branch/20160809/Ohana/src/libdvo/src/coordops.c
===================================================================
--- /branches/czw_branch/20160809/Ohana/src/libdvo/src/coordops.c	(revision 39687)
+++ /branches/czw_branch/20160809/Ohana/src/libdvo/src/coordops.c	(revision 39688)
@@ -84,4 +84,5 @@
 
   /** extra polynomial terms **/
+  // for ZPN, these are used to modify the radial distance and not the X,Y coords
   if ((coords[0].Npolyterms > 1) && (proj != PROJ_ZPN)) {
     X2 = X*X;
@@ -157,8 +158,17 @@
 	  ctht = 0.0;
 	} else {
-	  T = DEG_RAD / R;
-	  stht =   T / sqrt ( 1.0 + T*T);
-	  ctht = 1.0 / sqrt ( 1.0 + T*T);
+	  // T = DEG_RAD / R; // T in 1/radians
+	  // stht =   T / sqrt ( 1.0 + T*T);
+	  // ctht = 1.0 / sqrt ( 1.0 + T*T);
+
+	  T = RAD_DEG * R;
+	  stht = 1.0 / sqrt ( 1.0 + T*T);
+	  ctht =   T / sqrt ( 1.0 + T*T);
 	}
+	break;
+      case PROJ_SIN:
+	// R = (180/pi) cos (theta)
+	ctht = RAD_DEG * R;
+	stht = sqrt (1 - ctht*ctht);
 	break;
       case PROJ_STG:
@@ -167,9 +177,4 @@
 	ctht = sqrt (1 - stht*stht);
 	break;
-      case PROJ_SIN:
-	// R = (180/pi) cos (theta)
-	ctht = RAD_DEG * R;
-	stht = sqrt (1 - ctht*ctht);
-	break;
       case PROJ_ARC:
 	// R = 90 - theta (degrees)
@@ -177,36 +182,4 @@
 	stht = cos (RAD_DEG * R);
 	break;
-
-      case PROJ_ZPN:
-	// R = 90 - theta (degrees)
-	// this is wrong because we are ignoring the distortion
-	// XXX For now, just solve for terms up to n = 3
-
-	// Ro = (pi/180)(90 - theta)
-	// R = (180/pi)sum (P_i R^i)
-
-	if (UKIRT_ONLY) {
-	  double Ro = RAD_DEG * R;
-	  double P1 = coords[0].polyterms[0][1];
-	  double P3 = coords[0].polyterms[0][3];
-
-	  // find the roots of f(gamma) = P1 gamma + P3 gamma^3 - Ro
-	  // starting guess for gamma is Ro / P1
-	  double gamma = Ro / P1;
-
-	  int i;
-	  for (i = 0; i < 5; i++) {
-	    double F = P1*gamma + P3*gamma*gamma*gamma - Ro;
-	    double dFdgamma = P1 + 3.0*P3*gamma*gamma;
-
-	    double gamma_new = gamma - F / dFdgamma;
-	    gamma = gamma_new;
-	  }
-	  
-	  double theta = 90.0 - gamma * DEG_RAD ;
-	  ctht = cos (RAD_DEG * theta);
-	  stht = sin (RAD_DEG * theta);
-	  break;
-	}
 
       case PROJ_ZEA:
@@ -220,4 +193,51 @@
 	ctht = sqrt (1 - stht*stht);
 	break;
+
+      case PROJ_ZPN:
+
+	// the forward projection is:
+	// theta = atan2(stht, ctht)
+	// gamma = (pi/2 - theta) : theta in radians
+	// Ro = sum (P_i gamma^i)
+	// R  = (180/pi) Ro
+
+	// given R, we need to find theta:
+	// Ro = R * (pi / 180) = sum (P_i gamma^i)
+	// solve sum (P_i gamma^i) - Ro = 0 using Newton-Raphson
+
+	// use Ro to get a guess for gamma and iterate
+
+	{
+	  double Ro = RAD_DEG * R;
+	  
+	  // find the roots of f(gamma) - Ro = 0
+	  // starting guess for gamma is (Ro - P0) / P1
+	  double gamma = (Ro - coords[0].polyterms[0][0]) / coords[0].polyterms[1][0];
+	  
+	  int iter;
+	  for (iter = 0; iter < 5; iter++) {
+	    
+	    double Rc = 0.0; // this will hold the ander 
+	    double dR = 0.0;
+	    for (int i = coords[0].Npolyterms - 1; i > 1; i--) {
+	      double Pi = (i < 7) ? coords[0].polyterms[i][0] : coords[0].polyterms[i-7][1];
+	      Rc = (Rc + Pi)*gamma;
+	      dR = (dR + i*Pi)*gamma;
+	    }
+	    double P0 = coords[0].polyterms[0][0];
+	    double P1 = coords[0].polyterms[1][0];
+	    Rc = (Rc + P1)*gamma + P0;
+	    dR = (dR + P1);
+
+	    double gamma_new = gamma - (Rc - Ro) / dR;
+	    gamma = gamma_new;
+	  }
+	  
+	  double theta = 90.0 - gamma * DEG_RAD ;
+	  ctht = cos (RAD_DEG * theta);
+	  stht = sin (RAD_DEG * theta);
+	  break;
+	}
+
       default:
 	return (FALSE);
@@ -354,32 +374,30 @@
 
       case PROJ_ZPN:
+	// the forward projection is:
+	// theta = atan2(stht, ctht)
+	// gamma = (pi/2 - theta) : theta in radians
+	// Ro = sum (P_i gamma^i)
+	// R  = (180/pi) Ro
+
 	// Ro = (pi/180)(90 - theta)
 	// R = (180/pi)sum (P_i R^i)
+
+	// is ZPN defined for Npolyterms = 0 or 1?
+
 	ctht = hypot(sphi, cphi);
 	theta = atan2 (stht, ctht);
 
-	double Rc;
-	if (UKIRT_ONLY) {
-	  double P1 = coords[0].polyterms[0][1];
-	  double P3 = coords[0].polyterms[0][3];
-	  double gamma = RAD_DEG * (90 - DEG_RAD * theta);
-	  Rc = P1*gamma + P3*gamma*gamma*gamma;
-	} else {
-	  double Ro = RAD_DEG * (90 - DEG_RAD * theta);
-
-	  // i = 0 .. Npolyterms - 1 (1 <= Npolyterms <= 21)
-	  i = coords[0].Npolyterms - 1;
-	  Rc = 0.0;
-	  while (i > 0) {
-	    if (i < 7) {
-	      Rc = (Rc + coords[0].polyterms[0][i])*Ro;
-	    } else {
-	      Rc = (Rc + coords[0].polyterms[1][i-7])*Ro;
-	    }
-	    i --;
-	  }
-	  Rc += coords[0].polyterms[0][i];
+	double Ro;
+	double gamma = M_PI_2 - theta;
+
+	// i = 0 .. Npolyterms - 1 (1 <= Npolyterms <= 21)
+	Ro = 0.0;
+	for (i = coords[0].Npolyterms - 1; i > 0; i --) {
+	  double Pi = (i < 7) ? coords[0].polyterms[i][0] : coords[0].polyterms[i-7][1];
+	  Ro = (Ro + Pi)*gamma;
 	}
-	Rc = DEG_RAD * Rc;
+	Ro += coords[0].polyterms[0][0];
+
+	Rc = DEG_RAD * Ro;
 
 	*L = (ctht == 0.0) ? 0.0 : +Rc * sphi / ctht ;
@@ -609,7 +627,9 @@
 enum {COORD_TYPE_NONE, COORD_TYPE_PC, COORD_TYPE_ROT, COORD_TYPE_CD, COORD_TYPE_LIN};
 
+int GetRadialZPN (Coords *coords, Header *header);
+
 int GetCoords (Coords *coords, Header *header) {
   
-  int i, status, status1, status2, itmp, Polynomial, Polyterm;
+  int status, status1, status2, itmp, Polynomial, Polyterm;
   double Lambda, rotate, rotate1, rotate2, scale;
   double equinox;
@@ -666,8 +686,16 @@
       status &= gfits_scan (header, "PC002002", "%f",  1, &coords[0].pc2_2);
 
+      ctype = &coords[0].ctype[4];
+
+      // read the ZPN coeffients PV2_i (i = 0 < 14)
+      // ZPN is inconsistent with the other Polynomial types
+      if (!strcmp (ctype, "-ZPN")) { 
+	GetRadialZPN (coords, header);
+	break;
+      }
+
       /* 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");
       Polyterm = gfits_scan (header, "NPLYTERM", "%d", 1, &itmp);
@@ -723,4 +751,7 @@
       coords[0].pc2_1 =  sin(rotate*RAD_DEG) / Lambda;
       coords[0].pc2_2 =  cos(rotate*RAD_DEG);
+
+      // read the ZPN coeffients PV2_i (i = 0 < 14)
+      if (!strcmp (&coords[0].ctype[4], "-ZPN")) GetRadialZPN (coords, header);
       break;
 
@@ -744,28 +775,6 @@
       coords[0].pc2_2 /= scale;
 
-      if (!strcmp (&coords[0].ctype[4], "-ZPN")) {
-	int found;
-	for (i = 0; i < 14; i++) {
-	  char name[64];
-	  snprintf (name, 64, "PV2_%d", i);
-	  if (i < 7) {
-	    found = gfits_scan (header, name, "%f", 1, &coords[0].polyterms[0][i]);
-	  } else {
-	    found = gfits_scan (header, name, "%f", 1, &coords[0].polyterms[1][i-7]);
-	  }
-	  if ((i == 0) && !found) {
-	    coords[0].polyterms[0][0] = 0.0;
-	    continue;
-	  }
-	  if ((i == 1) && !found) {
-	    coords[0].polyterms[0][1] = 1.0;
-	    continue;
-	  }
-	  if (!found) {
-	    coords[0].Npolyterms = i;
-	    break;
-	  }
-	}
-      }
+      // read the ZPN coeffients PV2_i (i = 0 < 14)
+      if (!strcmp (&coords[0].ctype[4], "-ZPN")) GetRadialZPN (coords, header);
       break;
 
@@ -819,4 +828,37 @@
   }
   return (status);
+}
+
+int GetRadialZPN (Coords *coords, Header *header) {
+
+  // RA---ZPN can have up to 14 radial polynomial terms.  these are stored in
+  // polyterms[0][0] - [6][0] for the first 7 and [0][1] - [6][1] for the rest
+  // these terms are coeffients of a polynomial of the radial distances
+
+  // read the ZPN coeffients PV2_i (i = 0 < 14)
+  int found;
+  int Nmax = 0;
+  for (int i = 0; i < 14; i++) {
+    char name[64];
+    snprintf (name, 64, "PV2_%d", i);
+    if (i < 7) {
+      coords[0].polyterms[i][0] = 0.0;
+      found = gfits_scan (header, name, "%f", 1, &coords[0].polyterms[i][0]);
+    } else {
+      coords[0].polyterms[i-7][1] = 0.0;
+      found = gfits_scan (header, name, "%f", 1, &coords[0].polyterms[i-7][1]);
+    }
+    // PV2_1 is implicit if not present
+    if ((i == 1) && !found) {
+      coords[0].polyterms[1][0] = 1.0;
+      continue;
+    }
+    // set Npolyterms based on the largest coefficient found
+    if (found) {
+      Nmax = i;
+    }
+  }
+  coords[0].Npolyterms = Nmax + 1;
+  return TRUE;
 }