Index: trunk/psLib/src/astronomy/psCoord.c
===================================================================
--- trunk/psLib/src/astronomy/psCoord.c	(revision 2941)
+++ trunk/psLib/src/astronomy/psCoord.c	(revision 2983)
@@ -10,6 +10,6 @@
 *  @author GLG, MHPCC
 *
-*  @version $Revision: 1.44 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-01-10 19:47:10 $
+*  @version $Revision: 1.45 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-01-13 22:45:28 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -24,4 +24,5 @@
 #include "psConstants.h"
 #include "psError.h"
+#include "psLogMsg.h"
 #include "psAstronomyErrors.h"
 #include <math.h>
@@ -183,4 +184,23 @@
 
 /******************************************************************************
+XXX: Private Function.
+ 
+piNormalize(): take an input angle in radians and convert it to the range 0:PI.
+ *****************************************************************************/
+psF32 piNormalize(psF32 angle)
+{
+    while (angle < FLT_EPSILON) {
+        angle+=M_PI;
+    }
+
+    while (angle >= M_PI) {
+        angle-=M_PI;
+    }
+    return(angle);
+}
+
+/******************************************************************************
+XXX: We convert Right Ascension angles to the range 0:PI.  Is that acceptable?
+XXX: Should we do something for Declination as well?
  *****************************************************************************/
 psSphere* psSphereTransformApply(psSphere* out,
@@ -210,5 +230,6 @@
     psF64 phi = atan2(eq56, eq57) + transform->alphaP;
 
-    out->r = phi;
+
+    out->r = piNormalize(phi);
     out->d = theta;
 
@@ -281,4 +302,5 @@
 XXX: Waiting for the definition of the PS_PROJ_PAR projection.
 XXX: Waiting for the definition of the PS_PROJ_GLS projection.
+XXX: Must apply scaling at the end.
  *****************************************************************************/
 psPlane* psProject(const psSphere* coord,
@@ -288,13 +310,44 @@
     PS_PTR_CHECK_NULL(projection, NULL);
 
-    //    float R = 0.0;
-    //    float alpha = 0.0;
     psPlane* tmp = (psPlane* ) psAlloc(sizeof(psPlane));
-
     if (projection->type == PS_PROJ_TAN) {
-        tmp->x = cos(coord->r) * sin(coord->d) / sin(coord->r);
-        tmp->y = -cos(coord->r) * cos(coord->d) / sin(coord->r);
+
+        /*
+                // ********************************************
+                // From the ADD
+                // ********************************************
+                // psLogMsg(__func__, PS_LOG_WARN, "WARNING: psProject(): These projections don't work.");
+                // return(NULL);
+                psF32 sinTheta;
+                psF32 cosThetaCosPhi;
+                psF32 cosThetaSinPhi;
+         
+                sinTheta = (sin(coord->d) * sin(projection->D)) + 
+                           (cos(coord->d) * cos(projection->D) * (cos(coord->r - projection->R)));
+                cosThetaCosPhi = (sin(coord->d) * cos(projection->D)) - 
+                                 (cos(coord->d) * sin(projection->D) * (cos(coord->r - projection->R)));
+                cosThetaSinPhi = - cos(coord->d) * sin(coord->r - projection->R);
+         
+                tmp->x =  -cosThetaSinPhi / sinTheta;
+                tmp->y = cosThetaCosPhi / sinTheta;
+         
+        */
+
+        // ********************************************
+        // From the mathworks: http://mathworld.wolfram.com/GnomonicProjection.html
+        // delta_0 and phi_1 are the projection centers, not the point being projected.
+        // (delta_0, phi_1) == (projection->R, projection->D)
+        // ********************************************
+        psF32 cosC = (sin(projection->D) * sin(coord->d)) +
+                     (cos(projection->D) * cos(coord->d) * cos(coord->r - projection->R));
+        tmp->x = (cos(coord->d) * sin(coord->r - projection->R)) / cosC;
+        tmp->y = ((cos(projection->D) * sin(coord->d)) -
+                  (sin(projection->D) * cos(coord->d) * cos(coord->r - projection->R))) /
+                 cosC;
 
     } else if (projection->type == PS_PROJ_SIN) {
+        // psLogMsg(__func__, PS_LOG_WARN, "WARNING: psProject(): These projections don't work.");
+        // return(NULL);
+
         tmp->x = cos(coord->d) * sin(coord->d);
         tmp->y = -cos(coord->d) * cos(coord->d);
@@ -332,6 +385,8 @@
 
 /******************************************************************************
+XXX: Waiting for algorithms.
 XXX: Waiting for the definition of the PS_PROJ_PAR projection.
 XXX: Waiting for the definition of the PS_PROJ_GLS projection.
+XXX: Must apply scaling at the beginning.
  *****************************************************************************/
 psSphere* psDeproject(const psPlane* coord,
@@ -340,5 +395,6 @@
     PS_PTR_CHECK_NULL(coord, NULL);
     PS_PTR_CHECK_NULL(projection, NULL);
-
+    //    psLogMsg(__func__, PS_LOG_WARN, "WARNING: psDeproject(): These projections don't work.");
+    //    return(NULL);
     float R = 0.0;
     float chu = 0.0;
@@ -348,7 +404,39 @@
 
     if (projection->type == PS_PROJ_TAN) {
-        tmp->r = atan2(-coord->y, coord->x);
-        R = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y));
-        tmp->d = atan(1.0 / R);
+        /*
+                // ********************************************
+                // From the ADD
+                // ********************************************
+                psF32 phi;
+                psF32 theta;
+                psF32 sinDelta;
+                psF32 cosDeltaCosAlphaMinusAlphaP;
+                psF32 cosDeltaSinAlphaMinusAlphaP;
+         
+                phi = atan2(-coord->y, coord->x);
+                R = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y));
+                theta = atan(1.0 / R);
+         
+                sinDelta = (sin(theta) * sin(projection->D)) + 
+                           (cos(theta) * cos(projection->D) * cos(phi));
+                cosDeltaCosAlphaMinusAlphaP = 
+                           (sin(theta) * cos(projection->D)) - 
+                           (cos(theta) * sin(projection->D) * cos(phi));
+                cosDeltaSinAlphaMinusAlphaP = -cos(theta) * sin(phi);
+                tmp->d = asin(sinDelta);
+                tmp->r = atan2(cosDeltaSinAlphaMinusAlphaP, cosDeltaCosAlphaMinusAlphaP) + projection->R;
+        */
+        // ********************************************
+        // From the mathworks: http://mathworld.wolfram.com/GnomonicProjection.html
+        // delta_0 and phi_1 are the projection centers, not the point being projected.
+        // (delta_0, phi_1) == (projection->R, projection->D)
+        // ********************************************
+        psF32 row = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y));
+        psF32 C = atan(row);
+        tmp->d = asin((cos(C) * sin(projection->D)) + ((coord->y * sin(C) * cos(projection->D)) / row));
+        psF32 tmpAtan = atan((coord->x * sin(C)) /
+                             ((row * cos(projection->D) * cos(C)) - (coord->y * sin(projection->D) * sin(C))));
+        tmp->r = projection->R + tmpAtan;
+
 
     } else if (projection->type == PS_PROJ_SIN) {
@@ -415,4 +503,13 @@
     PS_PTR_CHECK_NULL(position2, NULL);
 
+    if (position1->d >= DEG_TO_RAD(90.0)) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: psDeproject(): position1->d is larger than 90 degrees.  Returning NULL.\n");
+        return(NULL);
+    }
+    if (position2->d >= DEG_TO_RAD(90.0)) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: psDeproject(): position2->d is larger than 90 degrees.  Returning NULL.\n");
+        return(NULL);
+    }
+
     psPlane* lin;
     psProjection proj;
@@ -426,15 +523,8 @@
         proj.type = PS_PROJ_TAN;
 
-        printf("------------ position 1 is (%f, %f)\n", position1->r, position1->d);
-        printf("------------ position 2 is (%f, %f)\n", position2->r, position2->d);
         lin = psProject(position1, &proj);
-        printf("------------ projected position 1 is (%f, %f)\n", lin->y, lin->x);
         lin = psProject(position2, &proj);
-        printf("------------ projected position 2 is (%f, %f)\n", lin->y, lin->x);
         tmp = psDeproject(lin, &proj);
-        printf("------------ deprojected position 2 is (%f, %f)\n", tmp->r, tmp->d);
-
         psFree(lin);
-
         // XXX: Do we need to convert units in tmp?
         return (tmp);