Index: /branches/czw_branch/20160809/psLib/src/astro/psCoord.c
===================================================================
--- /branches/czw_branch/20160809/psLib/src/astro/psCoord.c	(revision 39683)
+++ /branches/czw_branch/20160809/psLib/src/astro/psCoord.c	(revision 39684)
@@ -234,4 +234,7 @@
     const psPlane* coords)
 {
+  if (!transform) {
+    fprintf (stderr, "problem\n");
+  }
     PS_ASSERT_PTR_NON_NULL(transform, NULL);
     PS_ASSERT_PTR_NON_NULL(transform->x, NULL);
@@ -305,5 +308,7 @@
 void projectionFree(psProjection *p)
 {
-    // There are no dynamically allocated items
+  if (!p) return;
+  if (!p->radial) return;
+  psFree (p->radial);
 }
 
@@ -321,4 +326,5 @@
     p->Ys = Ys;
     p->type = type;
+    p->radial = NULL;
 
     psMemSetDeallocator(p, (psFreeFunc) projectionFree);
@@ -368,4 +374,5 @@
       case PS_PROJ_ZEA:
       case PS_PROJ_ZPL:
+      case PS_PROJ_ZPN:
         class = PS_PROJECTION_CLASS_ZENITHAL;
         break;
@@ -442,4 +449,35 @@
 	    out->y = -zeta * cosPhiCT;
 	    break;
+
+	  case PS_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?
+
+	    double cosTheta = hypot(sinPhiCT, cosPhiCT);
+	    double theta = atan2 (sinTheta, cosTheta);
+	    
+	    double gamma = M_PI_2 - theta;
+	    
+	    // i = 0 .. N - 1 (1 <= Npolyterms <= 21)
+	    double Ro = 0.0;
+	    for (int i = projection->radial->n - 1; i > 0; i--) {
+	      double Pi = projection->radial->data.F64[i];
+	      Ro = (Ro + Pi)*gamma;
+	    }
+	    Ro += projection->radial->data.F64[0];
+	    
+	    out->x = (cosTheta == 0.0) ? 0.0 : +Ro * sinPhiCT / cosTheta ;
+	    out->y = (cosTheta == 0.0) ? 0.0 : -Ro * cosPhiCT / cosTheta ;
+	    break;
+	  }
+
 	  default:
 	    psAbort("invalid projection");
@@ -525,4 +563,5 @@
       case PS_PROJ_ZEA:
       case PS_PROJ_ZPL:
+      case PS_PROJ_ZPN:
         class = PS_PROJECTION_CLASS_ZENITHAL;
         break;
@@ -574,9 +613,10 @@
 	  case PS_PROJ_SIN:
             // Orhtographic deprojection
+            cosTheta = R;
             sinTheta = sqrt (1 - R*R);
-            cosTheta = R;
             break;
 	  case PS_PROJ_STG:
-            psAbort("STG not defined");
+	    sinTheta = (4 - R) / (4 + R);
+	    cosTheta = sqrt (1 - sinTheta*sinTheta);
             break;
 	  case PS_PROJ_ZEA:
@@ -587,4 +627,48 @@
             cosTheta = sqrt (1 - PS_SQR(sinTheta));
             break;
+
+      case PS_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
+
+	{
+	  // find the roots of f(gamma) - R = 0
+	  // starting guess for gamma is (R - P0) / P1
+	  double gamma = (R - projection->radial->data.F64[0]) / projection->radial->data.F64[1];
+	  
+	  for (int iter = 0; iter < 5; iter++) {
+	    
+	    double Rc = 0.0; // this will hold the ander 
+	    double dR = 0.0;
+	    for (int i = projection->radial->n - 1; i > 1; i--) {
+	      double Pi = projection->radial->data.F64[i];
+	      Rc = (Rc + Pi)*gamma;
+	      dR = (dR + i*Pi)*gamma;
+	    }
+	    double P0 = projection->radial->data.F64[0];
+	    double P1 = projection->radial->data.F64[1];
+	    Rc = (Rc + P1)*gamma + P0;
+	    dR = (dR + P1);
+
+	    double gamma_new = gamma - (Rc - R) / dR;
+	    gamma = gamma_new;
+	  }
+	  
+	  double theta = M_PI_2 - gamma ;
+	  cosTheta = cos (theta);
+	  sinTheta = sin (theta);
+	  break;
+	}
+	
 	  default:
             psAbort("invalid projection");
@@ -1183,43 +1267,23 @@
 
     switch (type) {
-      case PS_PROJ_LIN:
-        psStringAppend (&name, "%s-LIN", prefix);
-        return name;
-      case PS_PROJ_PLY:
-        psStringAppend (&name, "%s-PLY", prefix);
-        return name;
-      case PS_PROJ_WRP:
-        psStringAppend (&name, "%s-WRP", prefix);
-        return name;
-      case PS_PROJ_TAN:
-        psStringAppend (&name, "%s-TAN", prefix);
-        return name;
-      case PS_PROJ_TNX:
-        psStringAppend (&name, "%s-TNX", prefix);
-        return name;
-      case PS_PROJ_DIS:
-        psStringAppend (&name, "%s-DIS", prefix);
-        return name;
-      case PS_PROJ_SIN:
-        psStringAppend (&name, "%s-SIN", prefix);
-        return name;
-      case PS_PROJ_STG:
-        psStringAppend (&name, "%s-STG", prefix);
-        return name;
-      case PS_PROJ_AIT:
-        psStringAppend (&name, "%s-AIT", prefix);
-        return name;
-      case PS_PROJ_PAR:
-        psStringAppend (&name, "%s-PAR", prefix);
-        return name;
-      case PS_PROJ_GLS:
-        psStringAppend (&name, "%s-GLS", prefix);
-        return name;
-      case PS_PROJ_CAR:
-        psStringAppend (&name, "%s-CAR", prefix);
-        return name;
-      case PS_PROJ_MER:
-        psStringAppend (&name, "%s-MER", prefix);
-        return name;
+      case PS_PROJ_LIN: psStringAppend (&name, "%s-LIN", prefix); return name;
+      case PS_PROJ_PLY: psStringAppend (&name, "%s-PLY", prefix); return name;
+      case PS_PROJ_WRP: psStringAppend (&name, "%s-WRP", prefix); return name;
+      case PS_PROJ_TAN: psStringAppend (&name, "%s-TAN", prefix); return name;
+
+      case PS_PROJ_DIS: psStringAppend (&name, "%s-DIS", prefix); return name;
+      case PS_PROJ_SIN: psStringAppend (&name, "%s-SIN", prefix); return name;
+      case PS_PROJ_STG: psStringAppend (&name, "%s-STG", prefix); return name;
+      case PS_PROJ_TNX: psStringAppend (&name, "%s-TNX", prefix); return name;
+
+      case PS_PROJ_ZEA: psStringAppend (&name, "%s-ZEA", prefix); return name;
+      case PS_PROJ_ZPL: psStringAppend (&name, "%s-ZPL", prefix); return name;
+      case PS_PROJ_ZPN: psStringAppend (&name, "%s-ZPN", prefix); return name;
+      case PS_PROJ_AIT: psStringAppend (&name, "%s-AIT", prefix); return name;
+
+      case PS_PROJ_PAR: psStringAppend (&name, "%s-PAR", prefix); return name;
+      case PS_PROJ_GLS: psStringAppend (&name, "%s-GLS", prefix); return name;
+      case PS_PROJ_CAR: psStringAppend (&name, "%s-CAR", prefix); return name;
+      case PS_PROJ_MER: psStringAppend (&name, "%s-MER", prefix); return name;
 
       default:
@@ -1238,8 +1302,11 @@
     if (!strcmp (&name[4], "-WRP")) return PS_PROJ_WRP;
     if (!strcmp (&name[4], "-TAN")) return PS_PROJ_TAN;
-    if (!strcmp (&name[4], "-TNX")) return PS_PROJ_TNX;
     if (!strcmp (&name[4], "-DIS")) return PS_PROJ_DIS;
     if (!strcmp (&name[4], "-SIN")) return PS_PROJ_SIN;
     if (!strcmp (&name[4], "-STG")) return PS_PROJ_STG;
+    if (!strcmp (&name[4], "-TNX")) return PS_PROJ_TNX;
+    if (!strcmp (&name[4], "-ZEA")) return PS_PROJ_ZEA;
+    if (!strcmp (&name[4], "-ZPL")) return PS_PROJ_ZPL;
+    if (!strcmp (&name[4], "-ZPN")) return PS_PROJ_ZPN;
     if (!strcmp (&name[4], "-AIT")) return PS_PROJ_AIT;
     if (!strcmp (&name[4], "-PAR")) return PS_PROJ_PAR;
Index: /branches/czw_branch/20160809/psLib/src/astro/psCoord.h
===================================================================
--- /branches/czw_branch/20160809/psLib/src/astro/psCoord.h	(revision 39683)
+++ /branches/czw_branch/20160809/psLib/src/astro/psCoord.h	(revision 39684)
@@ -130,4 +130,5 @@
     PS_PROJ_ZEA,                        ///< Sine projection
     PS_PROJ_ZPL,                        ///< Sine projection
+    PS_PROJ_ZPN,                        ///< Sine projection
     PS_PROJ_AIT,                        ///< Aitoff projection
     PS_PROJ_PAR,                        ///< Par projection
@@ -150,6 +151,6 @@
     double Ys;                         ///< plate-scale in Y direction
     psProjectionType type;             ///< Projection type
-}
-psProjection;
+    psVector *radial;	       ///< radial distortion terms 
+} psProjection;
 
 /** Allocates a psPlane