Index: /branches/eam_rel8_b0/psLib/src/astro/psCoord.c
===================================================================
--- /branches/eam_rel8_b0/psLib/src/astro/psCoord.c	(revision 5553)
+++ /branches/eam_rel8_b0/psLib/src/astro/psCoord.c	(revision 5554)
@@ -10,6 +10,6 @@
 *  @author GLG, MHPCC
 *
-*  @version $Revision: 1.88.4.1 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-11-18 20:22:35 $
+*  @version $Revision: 1.88.4.2 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-11-21 19:04:02 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -76,5 +76,5 @@
 XXX: Use the ADD version which is based on determinants.
  *****************************************************************************/
-psPlaneTransform *p_psPlaneTransformLinearInvert(psPlaneTransform *transform)
+psPlaneTransform *p_psPlaneTransformLinearInvert_MHPCC(psPlaneTransform *transform)
 {
     PS_ASSERT_PTR_NON_NULL(transform, 0);
@@ -104,26 +104,26 @@
     }
 
-    psPlaneTransform *out = psPlaneTransformAlloc(2, 2);
+    psPlaneTransform *out = psPlaneTransformAlloc(1, 1);
 
     /* This is sample code from IfA.  It didn't work initially, and I did not
        spend any time debugging it.
-
-        psF64 a = transform->x->coeff[1][0];
-        psF64 b = transform->x->coeff[0][1];
-        psF64 c = transform->y->coeff[1][0];
-        psF64 d = transform->y->coeff[0][1];
-        psF64 e = transform->x->coeff[0][0];
-        psF64 f = transform->y->coeff[0][0];
-
-        psF64 invDet = 1.0 / (a * d - b * c); // Inverse of the determinant
-
-        // Not entirely sure why this works, but it appears to do so....................................!
-        out->x->coeff[1][0] = invDet * a;
-        out->x->coeff[0][1] = - invDet * b;
-        out->y->coeff[1][0] = - invDet * c;
-        out->y->coeff[0][1] = invDet * d;
-
-        out->x->coeff[0][0] = - invDet * (d * e + c * f);
-        out->y->coeff[0][0] = - invDet * (b * e + a * f);
+       XXX EAM : here is the correct matrix inversion code
+
+        psF64 r11 = transform->x->coeff[1][0];
+        psF64 r12 = transform->x->coeff[0][1];
+        psF64 r21 = transform->y->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); // Inverse of the determinant
+
+        out->x->coeff[1][0] = +invDet * r22;
+        out->x->coeff[0][1] = -invDet * r12;
+        out->y->coeff[1][0] = -invDet * r21;
+        out->y->coeff[0][1] = +invDet * r11;
+
+        out->x->coeff[0][0] = - invDet * (r22 * xo - r12 * yo);
+        out->y->coeff[0][0] = - invDet * (r11 * yo - r21 * xo);
     */
     out->x->coeff[0][0] = (-D + ((F*A)/C)) / (E - ((F*B)/C));
@@ -134,4 +134,55 @@
     out->y->coeff[0][1] =  1.0 / (F - ((C*E)/B));
 
+    return(out);
+}
+
+// XXX EAM : the above code yielded NaNs for the y coeffs.  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)
+{
+    PS_ASSERT_PTR_NON_NULL(transform, 0);
+    PS_ASSERT_PTR_NON_NULL(transform->x, 0);
+    PS_ASSERT_PTR_NON_NULL(transform->y, 0);
+
+    if (transform->x->nX != 1)
+        return NULL;
+    if (transform->y->nY != 1)
+        return NULL;
+    if (transform->x->nY > 1)
+        return NULL;
+    if (transform->y->nX > 1)
+        return NULL;
+
+    // this choice is consistent with the nOrder form for the polynomials
+    psPlaneTransform *out = psPlaneTransformAlloc(1, 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);
 }
@@ -389,35 +440,52 @@
     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
     psPlane* out = psPlaneAlloc();
 
-    // 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,
@@ -441,4 +509,10 @@
     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;
     psF64  phi   = 0.0;
@@ -450,26 +524,40 @@
     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,
@@ -479,10 +567,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