Index: trunk/psLib/src/astro/psCoord.c
===================================================================
--- trunk/psLib/src/astro/psCoord.c	(revision 6484)
+++ trunk/psLib/src/astro/psCoord.c	(revision 7079)
@@ -10,6 +10,6 @@
 *  @author GLG, MHPCC
 *
-*  @version $Revision: 1.118 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-02-24 23:43:15 $
+*  @version $Revision: 1.119 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-05-06 03:34:21 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -32,4 +32,5 @@
 #include <float.h>
 
+# define ELIXIR_CODE 1
 
 static void planeFree(psPlane *p)
@@ -375,14 +376,23 @@
         cosDelta = cos(coord->d);
 
+        # if ELIXIR_CODE
+
+        sinTheta =  cosDelta*cosAlpha*cosDp + sinDelta*sinDp;
+    sinPhiCT =  cosDelta*sinAlpha;
+    cosPhiCT =  cosDelta*cosAlpha*sinDp - sinDelta*cosDp;
+    # else
+
         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
-    switch (projection->type) {
+    cosPhiCT =  sinDelta*cosDp - cosDelta*sinDp*cosAlpha;
+    sinPhiCT = -cosDelta*sinAlpha;
+    # endif
+
+} else {
+    phi = coord->r - projection->R;
+    theta = coord->d - projection->D;
+}
+
+// Perform the specified projection
+switch (projection->type) {
     case PS_PROJ_TAN:
         // Gnomonic projection
@@ -468,9 +478,21 @@
     case PS_PROJ_TAN:
         // Gnonomic deprojection
+        // the MHPCC version here was fine
+        # if 0
+
+        sinPhi   = (R == 0) ? 0.0 : +x / R;
+        cosPhi   = (R == 0) ? 1.0 : -y / R;
         rho      = sqrt (1 + R*R);
+        sinTheta = 1 / rho;
+        cosTheta = R / rho;
+        # else
+
+            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;
+        # endif
+
         break;
     case PS_PROJ_SIN:
@@ -508,10 +530,19 @@
 
         // Convert from projection spherical coordinates
-        psF64 delta = asin(sinTheta*sinDp + cosTheta*cosDp*cosPhi);
-        psF64 sinAlphaF = -cosTheta*sinPhi;
-        psF64 cosAlphaF = -cosTheta*cosPhi*sinDp + sinTheta*cosDp;
+        // psLib versions:
+        # if ELIXIR_CODE
+        // XXX the elixir version : does the ADD have a sign error
+        psF64 delta = asin(sinTheta*sinDp - cosTheta*cosPhi*cosDp);
+        psF64 sinAlpha = cosTheta*sinPhi;
+        psF64 cosAlpha = sinTheta*cosDp + cosTheta*cosPhi*sinDp;
+        # else
+
+            psF64 delta = asin(sinTheta*sinDp + cosTheta*cosDp*cosPhi);
+        psF64 sinAlpha = -cosTheta*sinPhi;
+        psF64 cosAlpha = -cosTheta*cosPhi*sinDp + sinTheta*cosDp;
+        # endif
 
         out->d = delta;
-        out->r = atan2(sinAlphaF, cosAlphaF) + projection->R;
+        out->r = atan2(sinAlpha, cosAlpha) + projection->R;
     } else {
         out->r = phi   + projection->R;