Index: trunk/psLib/src/astro/psCoord.c
===================================================================
--- trunk/psLib/src/astro/psCoord.c	(revision 5588)
+++ trunk/psLib/src/astro/psCoord.c	(revision 5624)
@@ -10,6 +10,6 @@
 *  @author GLG, MHPCC
 *
-*  @version $Revision: 1.93 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-11-23 23:54:43 $
+*  @version $Revision: 1.94 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-11-30 02:00:00 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -61,12 +61,10 @@
 should rename this.
  
-To derive this transformation, start with the simple 2-by-2 matrix inversion
-based on discriminants.  This will invert the
-    (x_2, y_2) = matrix(a b c d) * vector (x, y)
-where there are no constant terms.  Then you substitute 
-    x_2 = x_1 - e
-    y_2 = y_1 - f
-for (x_2, y_2) to get the desired inverse.
+XXX: Use the ADD version which is based on determinants.
  *****************************************************************************/
+
+// XXX EAM : below is the code using the standard matrix representation.
+//           note that this inversion requires x->nX == 1, y->nY == 1 and
+//           x->nY <= 1, y->nX <= 1
 psPlaneTransform *p_psPlaneTransformLinearInvert(psPlaneTransform *transform)
 {
@@ -130,4 +128,32 @@
     //    printf("HMMM: out->y: (%f %f %f)\n", out->y->coeff[0][0], out->y->coeff[1][0], out->y->coeff[0][1]);
 
+
+    // unless the cross terms are available, set these matrix elements to 0
+    psF64 r12 = 0.0;
+    if (transform->x->nY == 1) {
+        r12 = transform->x->coeff[0][1];
+    }
+    psF64 r21 = 0.0;
+    if (transform->y->nX == 1) {
+        r21 = transform->y->coeff[1][0];
+    }
+    psF64 r11 = transform->x->coeff[1][0];
+    psF64 r22 = transform->y->coeff[0][1];
+    psF64 xo  = transform->x->coeff[0][0];
+    psF64 yo  = transform->y->coeff[0][0];
+
+    psF64 invDet = 1.0 / (r11 * r22 - r12 * r21);
+
+    // apply the results back to the polynomials
+    out->x->coeff[0][0] = -invDet * (r22 * xo - r12 * yo);
+    out->y->coeff[0][0] = -invDet * (r11 * yo - r21 * xo);
+    out->x->coeff[1][0] = +invDet * r22;
+    out->y->coeff[0][1] = +invDet * r11;
+    if (transform->x->nY == 1) {
+        out->x->coeff[0][1] = -invDet * r12;
+    }
+    if (transform->y->nX == 1) {
+        out->y->coeff[1][0] = -invDet * r21;
+    }
     return(out);
 }
@@ -327,6 +353,5 @@
 XXX: Private Function.
  
-piNormalize(): take an input angle in radians and convert it to the range
-0:2*PI.
+piNormalize(): take an input angle in radians and convert it to the range 0:2*PI.
  *****************************************************************************/
 psF32 piNormalize(psF32 angle)
@@ -380,6 +405,9 @@
     PS_ASSERT_PTR_NON_NULL(projection, NULL);
 
-    psF64   theta = 0.0;
-    psF64   phi   = 0.0;
+    psF64 phi, theta;
+    psF64 sinDp, cosDp, sinAlpha, cosAlpha, sinDelta, cosDelta;
+    psF64 sinTheta, cosPhiCT, sinPhiCT, zeta;
+
+    bool zenithal = (projection->type == PS_PROJ_TAN) ||(projection->type == PS_PROJ_SIN);
 
     // Allocate return value
@@ -391,29 +419,43 @@
     }
 
-    // Convert to projection spherical coordinate system
-    theta = asin( sin(coord->d)*sin(projection->D) +
-                  cos(coord->d)*cos(projection->D)*cos(coord->r-projection->R));
-    phi = atan2( -1.0*cos(coord->d)*sin(coord->r-projection->R),
-                 sin(coord->d)*cos(projection->D) - cos(coord->d)*sin(projection->D)*cos(coord->r-projection->R) );
+    if (zenithal) {
+        sinDp = sin(projection->D);
+        cosDp = cos(projection->D);
+        sinAlpha = sin(coord->r-projection->R);
+        cosAlpha = cos(coord->r-projection->R);
+        sinDelta = sin(coord->d);
+        cosDelta = cos(coord->d);
+
+        sinTheta =  sinDelta*sinDp + cosDelta*cosDp*cosAlpha;
+        cosPhiCT =  sinDelta*cosDp - cosDelta*sinDp*cosAlpha;
+        sinPhiCT = -cosDelta*sinAlpha;
+    } else {
+        phi = coord->r - projection->R;
+        theta = coord->d - projection->D;
+    }
 
     // Perform the specified projection
-    // Gnomonic projection
-    if (projection->type == PS_PROJ_TAN) {
-        out->x = (cos(theta)*sin(phi))/sin(theta);
-        out->y = (-1.0*cos(theta)*cos(phi))/sin(theta);
+    switch (projection->type) {
+    case PS_PROJ_TAN:
+        // Gnomonic projection
+        out->x = +sinPhiCT / sinTheta;
+        out->y = -cosPhiCT / sinTheta;
+        break;
+    case PS_PROJ_SIN:
         // Othrographic projection
-    } else if (projection->type == PS_PROJ_SIN) {
-        out->x = cos(theta)*sin(phi);
-        out->y = -1.0*cos(theta)*cos(phi);
+        out->x = +sinPhiCT;
+        out->y = -cosPhiCT;
+        break;
+    case PS_PROJ_AIT:
         // Hammer-Aitoff projection
-    } else if ( projection->type == PS_PROJ_AIT) {
-        psF64 zeta = 1.0/sqrt(0.5*(1.0+cos(theta)*cos(phi/2.0)));
+        zeta = 1.0/sqrt(0.5*(1.0+cos(theta)*cos(phi/2.0)));
         out->x = 2.0*zeta*cos(theta)*sin(phi/2.0);
         out->y = zeta*sin(theta);
+        break;
+    case PS_PROJ_PAR:
         // Parabolic projection
-    } else if ( projection->type == PS_PROJ_PAR) {
         out->x = phi*(2.0*cos(2.0*theta/3.0) - 1.0);
         out->y = M_PI*sin(theta/3.0);
-    } else {
+    default:
         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                 PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNKNOWN,
@@ -424,6 +466,6 @@
 
     // Apply plate scales
-    out->x *= projection->Xs;
-    out->y *= projection->Ys;
+    out->x /= projection->Xs;
+    out->y /= projection->Ys;
 
     // Return output
@@ -444,4 +486,10 @@
     PS_ASSERT_PTR_NON_NULL(coord, NULL);
     PS_ASSERT_PTR_NON_NULL(projection, NULL);
+
+    psF64 rho      = 0.0;
+    psF64 sinTheta = 0.0;
+    psF64 cosTheta = 0.0;
+    psF64 sinPhi   = 0.0;
+    psF64 cosPhi   = 0.0;
 
     psF64  theta = 0.0;
@@ -458,30 +506,42 @@
 
     // Remove plate scales
-    // XXX: Verify this.  EAM suggested we do a multiply, however that does
-    // not make sense if we also do the multiply in psProject().
-    psF64  x = coord->x/projection->Xs;
-    psF64  y = coord->y/projection->Ys;
+    psF64  x = coord->x*projection->Xs;
+    psF64  y = coord->y*projection->Ys;
+    psF64  R = sqrt(x*x + y*y);
+
+    bool zenithal = (projection->type == PS_PROJ_TAN) ||(projection->type == PS_PROJ_SIN);
 
     // Perform inverse projection
-    // Gnonomic deprojection
-    if ( projection->type == PS_PROJ_TAN) {
-        phi = atan(-1.0*x/y);
-        theta = atan(1.0/sqrt(x*x+y*y));
+    switch (projection->type) {
+    case PS_PROJ_TAN:
+        // Gnonomic deprojection
+        rho      = sqrt (1 + R*R);
+        sinTheta = 1 / rho;
+        cosTheta = R / rho;
+        sinPhi   = (R == 0) ? 0.0 : +x / R;
+        cosPhi   = (R == 0) ? 1.0 : -y / R;
+        break;
+    case PS_PROJ_SIN:
         // Orhtographic deprojection
-    } else if ( projection->type == PS_PROJ_SIN) {
-        phi = atan((-1.0*x)/y);
-        theta = atan( sqrt(1.0-(x*x+y*y)) / sqrt(x*x+y*y));
+        sinTheta = sqrt (1 - R*R);
+        cosTheta = R;
+        sinPhi   = (R == 0) ? 0.0 : +x / R;
+        cosPhi   = (R == 0) ? 1.0 : -y / R;
+        break;
+    case PS_PROJ_AIT:
         // Hammer-Aitoff deprojection
-    } else if ( projection->type == PS_PROJ_AIT) {
-        psF64 z = sqrt(1.0 - ((x/4.0)*(x/4.0)) - ((y/2.0)*(y/2.0)));
-        phi = 2.0*atan((z*x) / (2.0*(2.0*z*z-1.0)) );
-        theta = asin(y*z);
+        // XXX EAM : need range check on z^2 : must be > 0
+        // XXX EAM : old code, ADD, and elixir code are discrepant re x/4, y/2
+        rho = sqrt(1.0 - PS_SQR(x/4.0) - PS_SQR(y/2.0));
+        phi = 2.0*atan2((2.0*rho*rho-1.0), x*rho);
+        theta = asin(y*rho);
+        break;
+    case PS_PROJ_PAR:
         // Parabolic deprojection
-    } else if ( projection->type == PS_PROJ_PAR) {
-        psF64 rho = y/M_PI;
+        rho = y/M_PI;
         phi = x/(1.0 - 4.0*rho*rho);
         theta = 3.0*asin(rho);
-        // Invalid deprojection type
-    } else {
+        break;
+    default:
         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                 PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNKNOWN,
@@ -491,10 +551,19 @@
     }
 
-    // Convert from projection spherical coordinates
-    out->d = asin( sin(theta)*sin(projection->D) +
-                   cos(theta)*cos(projection->D)*cos(phi) );
-    out->r = projection->R + atan2( -1.0*cos(theta)*sin(phi),
-                                    sin(theta)*cos(projection->D) -
-                                    cos(theta)*sin(projection->D)*cos(phi) );
+    if (zenithal) {
+        psF64 sinDp = sin(projection->D);
+        psF64 cosDp = cos(projection->D);
+
+        // Convert from projection spherical coordinates
+        psF64 delta = asin(sinTheta*sinDp + cosTheta*cosDp*cosPhi);
+        psF64 sinAlphaF = -cosTheta*sinPhi;
+        psF64 cosAlphaF = -cosTheta*cosPhi*sinDp + sinTheta*cosDp;
+
+        out->d = delta;
+        out->r = atan2(sinAlphaF, cosAlphaF) + projection->R;
+    } else {
+        out->r = phi   + projection->R;
+        out->d = theta + projection->D;
+    }
 
     // Return sphere coordinate