Index: /branches/eam_branches/ipp-20150419/Ohana/src/libdvo/include/libdvo_astro.h
===================================================================
--- /branches/eam_branches/ipp-20150419/Ohana/src/libdvo/include/libdvo_astro.h	(revision 38182)
+++ /branches/eam_branches/ipp-20150419/Ohana/src/libdvo/include/libdvo_astro.h	(revision 38183)
@@ -26,4 +26,5 @@
   PROJ_GLS, // pseudocyl
   PROJ_PAR, // pseudocyl
+  PROJ_ZPN, // zenithal
 } OhanaProjection;
 
Index: /branches/eam_branches/ipp-20150419/Ohana/src/libdvo/src/coordops.c
===================================================================
--- /branches/eam_branches/ipp-20150419/Ohana/src/libdvo/src/coordops.c	(revision 38182)
+++ /branches/eam_branches/ipp-20150419/Ohana/src/libdvo/src/coordops.c	(revision 38183)
@@ -1,3 +1,4 @@
 # include <dvo.h>
+# define UKIRT_ONLY 1
 
 /* note that Coords.ctype carries the DEC (ctype2) value */
@@ -81,5 +82,5 @@
 
   /** extra polynomial terms **/
-  if (coords[0].Npolyterms > 1) {
+  if ((coords[0].Npolyterms > 1) && (proj != PROJ_ZPN)) {
     X2 = X*X;
     Y2 = Y*Y;
@@ -171,4 +172,37 @@
 	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:
       case PROJ_ZPL:
@@ -238,4 +272,5 @@
 int RD_to_LM (double *L, double *M, double ra, double dec, Coords *coords) {
 
+  int i;
   double phi, theta;
   double sphi, cphi, stht, ctht;
@@ -303,4 +338,5 @@
 	*M = -DEG_RAD * cphi;
 	return (stht > 0);
+
       case PROJ_ARC:
 	// R = 90 - theta
@@ -311,4 +347,39 @@
 	*M = (ctht == 0.0) ? 0.0 : -Rc * cphi / ctht ;
 	return (TRUE);
+
+      case PROJ_ZPN:
+	// Ro = (pi/180)(90 - theta)
+	// R = (180/pi)sum (P_i R^i)
+	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];
+	}
+	Rc = DEG_RAD * Rc;
+
+	*L = (ctht == 0.0) ? 0.0 : +Rc * sphi / ctht ;
+	*M = (ctht == 0.0) ? 0.0 : -Rc * cphi / ctht ;
+	return (TRUE);
+
       case PROJ_ZEA:
       case PROJ_ZPL:
@@ -386,5 +457,6 @@
 
   /** extra polynomial terms **/
-  if (coords[0].Npolyterms > 1) {
+  // polyterm values are only valid for some projections.  for PROJ_ZPN, they are applied to R, not X,Y
+  if ((coords[0].Npolyterms > 1) && (proj != PROJ_ZPN)) {
     for (i = 0; i < 10; i++) {
       // find derivatives at Xo, Yo
@@ -531,5 +603,5 @@
 int GetCoords (Coords *coords, Header *header) {
   
-  int status, status1, status2, itmp, Polynomial, Polyterm;
+  int i, status, status1, status2, itmp, Polynomial, Polyterm;
   double Lambda, rotate, rotate1, rotate2, scale;
   double equinox;
@@ -663,4 +735,29 @@
       coords[0].pc2_1 /= scale;
       coords[0].pc2_2 /= scale;
+
+      if (!strcmp (&coords[0].ctype[4], "-ZPN")) {
+	int found;
+	for (i = 0; i < 14; i++) {
+	  char name[64];
+	  sprintf (name, "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;
+	  }
+	}
+      }
       break;
 
@@ -895,4 +992,5 @@
   if (!strcmp(&ctype[4], "-ZEA")) return PROJ_ZEA;
   if (!strcmp(&ctype[4], "-ZPL")) return PROJ_ZPL;
+  if (!strcmp(&ctype[4], "-ZPN")) return PROJ_ZPN;
   if (!strcmp(&ctype[4], "-ARC")) return PROJ_ARC;
   if (!strcmp(&ctype[4], "-STG")) return PROJ_STG;
@@ -916,4 +1014,5 @@
     case PROJ_ZEA: strcpy(&ctype[4], "-ZEA"); return TRUE;
     case PROJ_ZPL: strcpy(&ctype[4], "-ZPL"); return TRUE;
+    case PROJ_ZPN: strcpy(&ctype[4], "-ZPN"); return TRUE;
     case PROJ_ARC: strcpy(&ctype[4], "-ARC"); return TRUE;
     case PROJ_STG: strcpy(&ctype[4], "-STG"); return TRUE;
@@ -937,4 +1036,5 @@
     case PROJ_ZEA:
     case PROJ_ZPL:
+    case PROJ_ZPN:
     case PROJ_ARC:
     case PROJ_STG: