Index: /trunk/Ohana/src/libohana/src/coordops.c
===================================================================
--- /trunk/Ohana/src/libohana/src/coordops.c	(revision 3295)
+++ /trunk/Ohana/src/libohana/src/coordops.c	(revision 3296)
@@ -4,4 +4,6 @@
 int XY_to_RD (double *ra, double *dec, double x, double y, Coords *coords) {
 
+  int Polynomial, Zenithal, Cartesian, PseudoCyl;
+  char *type;
   double L, M, X, Y, T, Z;
   double R, sphi, cphi, stht, ctht;
@@ -11,46 +13,54 @@
   *dec = 0;
   stht = ctht = 1;
+  type = &coords[0].ctype[4];
+  
+  Polynomial = !strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-WRP");
+  Zenithal   = !strcmp(type, "-PLY") || !strcmp(type, "-TAN") || !strcmp(type, "-SIN") || !strcmp(type, "-ZEA") || !strcmp(&coords[0].ctype[0], "MM");
+  Cartesian  = !strcmp(type, "-LIN") || !strcmp(type, "-WRP") || !strcmp(&coords[0].ctype[0], "GENE");
+  PseudoCyl  = !strcmp(type, "-AIT") || !strcmp(type, "-GLS") || !strcmp(type, "-PAR");
+  if (!Zenithal && !Cartesian && !PseudoCyl) return (FALSE);
 
   /** convert pixel coordinates to cartesian system **/
   X = coords[0].cdelt1*(x - coords[0].crpix1);
   Y = coords[0].cdelt2*(y - coords[0].crpix2);
-  if (!strcmp(&coords[0].ctype[4], "-PLY")) {
-    if (coords[0].Npolyterms > 2) {
-      X += coords[0].cdelt1*(x*x*coords[0].polyterms[0][0] + x*y*coords[0].polyterms[1][0] + y*y*coords[0].polyterms[2][0]);
-      Y += coords[0].cdelt2*(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) {
-      X += coords[0].cdelt1*(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]);
-      Y += coords[0].cdelt2*(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]);
-    }
-  }
 
   L = (X*coords[0].pc1_1 + Y*coords[0].pc1_2);
   M = (X*coords[0].pc2_1 + Y*coords[0].pc2_2);
-  /* in FITS ref, L,M = x, y  alpha, delta = phi, theta */
-
+
+  /** 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];
+    }
+  }
+
+  /**** Locally Cartesian Projections ****/
+  if (Cartesian) {
+    *ra  = L + coords[0].crval1;
+    *dec = M + coords[0].crval2;
+    return (TRUE);
+  }
+  
   /**** Zenithal Projections ****/
-  if (!strcmp(&coords[0].ctype[4], "-PLY") || 
-      !strcmp(&coords[0].ctype[4], "-TAN") || 
-      !strcmp(&coords[0].ctype[4], "-SIN") || 
-      !strcmp(&coords[0].ctype[4], "-ZEA") || 
-      !strcmp(&coords[0].ctype[0], "MM")) {
+  if (Zenithal) {
     R = hypot (L,M);
     if ((L == 0) && (M == 0)) {
       sphi = 0;
       cphi = 1;
-    }
-    else {
+    } else {
       sphi =  L / R;
       cphi = -M / R;
     }
 
-    if (!strcmp(&coords[0].ctype[4], "-PLY") || 
-	!strcmp(&coords[0].ctype[4], "-TAN")) {
+    if (!strcmp(type, "-PLY") || !strcmp(type, "-TAN") || !strcmp(type, "-DIS")) {
       if (R == 0) {
 	stht = 1.0;
 	ctht = 0.0;
-      }
-      else {
+      } else {
 	T = DEG_RAD / R;
 	stht =   T / sqrt ( 1.0 + T*T);
@@ -58,10 +68,9 @@
       }
     }
-    if (!strcmp(&coords[0].ctype[4], "-SIN") || 
-	!strcmp(&coords[0].ctype[0], "MM")) {
+    if (!strcmp(type, "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) {
       ctht = RAD_DEG * R;
       stht = sqrt (1 - ctht*ctht);
     }
-    if (!strcmp(&coords[0].ctype[4], "-ZEA")) {
+    if (!strcmp(type, "-ZEA")) {
       stht = 1 - 0.5*SQ(R*RAD_DEG);
       ctht = sqrt (1 - stht*stht);
@@ -79,68 +88,66 @@
     *ra  = DEG_RAD*alpha + coords[0].crval1;
     *dec = DEG_RAD*delta;
-  }
-  
-  /**** Locally Cartesian Projections ****/
-  if (!strcmp(&coords[0].ctype[4], "-LIN") || 
-      !strcmp(&coords[0].ctype[0], "GENE")) {
-    *ra  = L + coords[0].crval1;
-    *dec = M + coords[0].crval2;
+    return (TRUE);
   }
   
   /**** Other Conventional Projections ****/
-  if (!strcmp(&coords[0].ctype[4], "-AIT")) {
-    Z = sqrt (1.0 - SQ(RAD_DEG*0.25*L) - SQ(RAD_DEG*0.5*M));
-    alpha = 2.0 * DEG_RAD * atan2 (RAD_DEG*0.5*Z*L, 2.0*SQ(Z) - 1.0);
-    delta = DEG_RAD * asin (RAD_DEG*M*Z);
-    *ra  = alpha + coords[0].crval1;
-    *dec = delta + coords[0].crval2;
-  }
-  if (!strcmp(&coords[0].ctype[4], "-GLS")) {
-    /* L,M in degrees, alpha,delta in degrees */
-    alpha = L / cos (RAD_DEG * M);
-    delta = M;
-    *ra  = alpha + coords[0].crval1;
-    *dec = delta + coords[0].crval2;
-  }
-  if (!strcmp(&coords[0].ctype[4], "-PAR")) {
-    /* L,M in degrees, alpha,delta in degrees */
-    alpha = L / (1.0 - SQ(2.0*M/180));
-    delta = 3 * DEG_RAD * asin (M/180.0);
-    *ra  = alpha + coords[0].crval1;
-    *dec = delta + coords[0].crval2;
-  }
-  return (TRUE);
+  if (PseudoCyl) {
+    if (!strcmp(type, "-AIT")) {
+      Z = sqrt (1.0 - SQ(RAD_DEG*0.25*L) - SQ(RAD_DEG*0.5*M));
+      alpha = 2.0 * DEG_RAD * atan2 (RAD_DEG*0.5*Z*L, 2.0*SQ(Z) - 1.0);
+      delta = DEG_RAD * asin (RAD_DEG*M*Z);
+      *ra  = alpha + coords[0].crval1;
+      *dec = delta + coords[0].crval2;
+    }
+    if (!strcmp(type, "-GLS")) {
+      /* L,M in degrees, alpha,delta in degrees */
+      alpha = L / cos (RAD_DEG * M);
+      delta = M;
+      *ra  = alpha + coords[0].crval1;
+      *dec = delta + coords[0].crval2;
+    }
+    if (!strcmp(type, "-PAR")) {
+      /* L,M in degrees, alpha,delta in degrees */
+      alpha = L / (1.0 - SQ(2.0*M/180));
+      delta = 3 * DEG_RAD * asin (M/180.0);
+      *ra  = alpha + coords[0].crval1;
+      *dec = delta + coords[0].crval2;
+    }
+    return (TRUE);
+  }
+  return (FALSE);
 }
 
 int RD_to_XY (double *x, double *y, double ra, double dec, Coords *coords) {
 
+  char *type;
+  int i, status, Polynomial, Zenithal, Cartesian, PseudoCyl;
   double phi, theta;
-  double tmp_d;
-  double X, Y, sphi, cphi, stht;
+  double determ;
+  double X, Y, L, M, Lo, Mo, dL, dM;
+  double sphi, cphi, stht;
   double salp, calp, sdel, cdel, sdp, cdp;
   double P, A, Rc;
-  int status;
-
-  if (!strcmp(&coords[0].ctype[4], "-PLY")) {
-    /* fprintf (stderr, "approximate to polynomial TAN plane fit\n");  */
-  }
-
-  X = Y = 1;
+
   status = TRUE;
   *x = 0;
   *y = 0;
+  type = &coords[0].ctype[4];
+  L = M = 0;
+
+  Polynomial = !strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-WRP");
+  Zenithal   = !strcmp(type, "-PLY") || !strcmp(type, "-TAN") || !strcmp(type, "-SIN") || !strcmp(type, "-ZEA") || !strcmp(&coords[0].ctype[0], "MM");
+  Cartesian  = !strcmp(type, "-LIN") || !strcmp(type, "-WRP") || !strcmp(&coords[0].ctype[0], "GENE");
+  PseudoCyl  = !strcmp(type, "-AIT") || !strcmp(type, "-GLS") || !strcmp(type, "-PAR");
+  if (!Zenithal && !Cartesian && !PseudoCyl) return (FALSE);
 
   /**** Locally Cartesian Projections ****/
-  if (!strcmp(&coords[0].ctype[0], "GENE") || !strcmp(&coords[0].ctype[4], "-LIN")) {
-    X = (ra  - coords[0].crval1);
-    Y = (dec - coords[0].crval2);
+  if (Cartesian) {
+    L = (ra  - coords[0].crval1);
+    M = (dec - coords[0].crval2);
   }
   
   /**** Zenithal Projections ****/
-  if (!strcmp(&coords[0].ctype[4], "-PLY") || 
-      !strcmp(&coords[0].ctype[4], "-TAN") || 
-      !strcmp(&coords[0].ctype[4], "-ZEA") || 
-      !strcmp(&coords[0].ctype[4], "-SIN") || 
-      !strcmp(&coords[0].ctype[0], "MM")) {
+  if (Zenithal)  {
     sdp  = sin(RAD_DEG*coords[0].crval2);
     cdp  = cos(RAD_DEG*coords[0].crval2);
@@ -155,16 +162,16 @@
     if (stht < 0) status = FALSE;
 
-    if (!strcmp(&coords[0].ctype[4], "-PLY") || !strcmp(&coords[0].ctype[4], "-TAN") ) {
-      X =  DEG_RAD * sphi / stht;
-      Y = -DEG_RAD * cphi / stht;
-    }
-    if (!strcmp(&coords[0].ctype[4], "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) {
-      X =  DEG_RAD * sphi;
-      Y = -DEG_RAD * cphi;
-    }
-    if (!strcmp(&coords[0].ctype[4], "-ZEA")) {
+    if (!strcmp(type, "-PLY") || !strcmp(type, "-TAN") || !strcmp(type, "-DIS")) {
+      L =  DEG_RAD * sphi / stht;
+      M = -DEG_RAD * cphi / stht;
+    }
+    if (!strcmp(type, "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) {
+      L =  DEG_RAD * sphi;
+      M = -DEG_RAD * cphi;
+    }
+    if (!strcmp(type, "-ZEA")) {
       Rc = DEG_RAD * M_SQRT2 / sqrt (1 + stht);
-      X =  Rc * sphi;
-      Y = -Rc * cphi;
+      L =  Rc * sphi;
+      M = -Rc * cphi;
       status = TRUE;
     }
@@ -172,33 +179,61 @@
 
   /**** Other Standard Projections ****/
-  if (!strcmp(&coords[0].ctype[4], "-AIT")) {
-    phi = RAD_DEG*(ra - coords[0].crval1);
-    theta = RAD_DEG*(dec - coords[0].crval2);
-    P = 1.0 + cos (theta) * cos (0.5*phi);
-    if (P != 0.0) {
-      A =  DEG_RAD * sqrt (2.0 / P);
-      X =  2.0 * A * cos (theta) * sin (0.5*phi);
-      Y =  A * sin (theta);
-    } else { 
-      X =  0.0;
-      Y =  0.0;
-    }	
-  }
-  if (!strcmp(&coords[0].ctype[4], "-GLS")) {
-    phi = ra - coords[0].crval1;
-    theta = dec - coords[0].crval2;
-    X = phi * cos(RAD_DEG * theta);
-    Y = theta;
-  }
-  if (!strcmp(&coords[0].ctype[4], "-PAR")) {
-    phi = ra - coords[0].crval1;
-    theta = dec - coords[0].crval2;
-    X = phi * (2.0*cos(2*RAD_DEG*theta/3.0) - 1);
-    Y = 180.0 * sin (RAD_DEG*theta/3.0);
-  }
-    
-  tmp_d = 1.0 / (coords[0].pc1_1*coords[0].pc2_2 - coords[0].pc1_2*coords[0].pc2_1);
-  *x = tmp_d * (coords[0].pc2_2*X - coords[0].pc1_2*Y) / coords[0].cdelt1 + coords[0].crpix1;
-  *y = tmp_d * (coords[0].pc1_1*Y - coords[0].pc2_1*X) / coords[0].cdelt2 + coords[0].crpix2;
+  if (PseudoCyl) {
+    if (!strcmp(type, "-AIT")) {
+      phi = RAD_DEG*(ra - coords[0].crval1);
+      theta = RAD_DEG*(dec - coords[0].crval2);
+      P = 1.0 + cos (theta) * cos (0.5*phi);
+      if (P != 0.0) {
+	A =  DEG_RAD * sqrt (2.0 / P);
+	L =  2.0 * A * cos (theta) * sin (0.5*phi);
+	M =  A * sin (theta);
+      } else { 
+	L =  0.0;
+	M =  0.0;
+      }	
+    }
+    if (!strcmp(type, "-GLS")) {
+      phi = ra - coords[0].crval1;
+      theta = dec - coords[0].crval2;
+      L = phi * cos(RAD_DEG * theta);
+      M = theta;
+    }
+    if (!strcmp(type, "-PAR")) {
+      phi = ra - coords[0].crval1;
+      theta = dec - coords[0].crval2;
+      L = phi * (2.0*cos(2*RAD_DEG*theta/3.0) - 1);
+      M = 180.0 * sin (RAD_DEG*theta/3.0);
+    }
+  }
+
+  /* convert L,M to X,Y */
+  determ = 1.0 / (coords[0].pc1_1*coords[0].pc2_2 - coords[0].pc1_2*coords[0].pc2_1);
+  X = determ * (coords[0].pc2_2*L - coords[0].pc1_2*M);
+  Y = determ * (coords[0].pc1_1*M - coords[0].pc2_1*L);
+
+  /** extra polynomial terms **/
+  if (Polynomial) {
+    for (i = 0; i < 3; i++) {
+      Lo = (X*coords[0].pc1_1 + Y*coords[0].pc1_2);
+      Mo = (X*coords[0].pc2_1 + Y*coords[0].pc2_2);
+      if (coords[0].Npolyterms > 1) {
+	Lo += X*X*coords[0].polyterms[0][0] + X*Y*coords[0].polyterms[1][0] + Y*Y*coords[0].polyterms[2][0];
+	Mo += 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) {
+	Lo += 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];
+	Mo += 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];
+      }
+      dL = (L - Lo);
+      dM = (M - Mo);
+
+      X += determ * (coords[0].pc2_2*dL - coords[0].pc1_2*dM);
+      Y += determ * (coords[0].pc1_1*dM - coords[0].pc2_1*dL);
+    }
+  }
+  /* check for correct size (iterate?) */
+
+  *x = X / coords[0].cdelt1 + coords[0].crpix1;
+  *y = Y / coords[0].cdelt2 + coords[0].crpix2;
 
   return (status);
@@ -391,2 +426,28 @@
    PCA2X0Y3 = coords.polyterm[6][1] = y^2                                               
 */
+
+# if (0)
+
+  /** convert pixel coordinates to cartesian system **/
+  X = coords[0].cdelt1*(x - coords[0].crpix1);
+  Y = coords[0].cdelt2*(y - coords[0].crpix2);
+  if (Polynomi) {
+    if (coords[0].Npolyterms > 2) {
+      X += coords[0].cdelt1*(x*x*coords[0].polyterms[0][0] + x*y*coords[0].polyterms[1][0] + y*y*coords[0].polyterms[2][0]);
+      Y += coords[0].cdelt2*(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) {
+      X += coords[0].cdelt1*(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]);
+      Y += coords[0].cdelt2*(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]);
+    }
+  }
+
+  L = (X*coords[0].pc1_1 + Y*coords[0].pc1_2);
+  M = (X*coords[0].pc2_1 + Y*coords[0].pc2_2);
+/** this code is the old method used for higher order terms.  they
+    are essentially 6th order, with weird coupled terms.
+    I don't think any real data used these terms, but they should 
+    be re-calculated, I would think 
+**/
+
+# endif