Index: trunk/psLib/src/astronomy/psCoord.c
===================================================================
--- trunk/psLib/src/astronomy/psCoord.c	(revision 1287)
+++ trunk/psLib/src/astronomy/psCoord.c	(revision 1293)
@@ -10,6 +10,6 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-07-23 03:13:39 $
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-07-24 02:42:59 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -22,5 +22,7 @@
 #include "psCoord.h"
 #include "psMemory.h"
-
+#include "psAbort.h"
+#include <math.h>
+#include <float.h>
 
 psPlane *psPlaneTransformApply(psPlane *out,
@@ -127,4 +129,117 @@
 }
 
+// XXX: Is this the correct way to calculate this?
+float cot(float x)
+{
+    return(1.0 / atan(x));
+}
+
+float arg(float x, float y)
+{
+    if (x > 0) {
+        return(atan(y/x));
+    } else if ((x==0) && (y==0)) {
+        return(0.5 * M_PI);
+    } else if ((x==0) && (y==0)) {
+        return(-0.5 * M_PI);
+    } else if ((x==0) && (y==0)) {
+        return(M_PI + atan(y/x));
+    } else if ((x==0) && (y==0)) {
+        return(-M_PI + atan(y/x));
+    }
+
+    psAbort(__func__, "Unacceptable range for (arg(%f, %f).\n", x, y);
+    return(0.0);
+}
+
+
+psPlane *psProject(const psSphere *coord,
+                   const psProjection *projection)
+{
+    float R = 0.0;
+    float alpha = 0.0;
+    psPlane *tmp= (psPlane *) psAlloc(sizeof(psPlane));
+
+    if (projection->type == PS_PROJ_TAN) {
+        R = cot(coord->r) * (180.0 / M_PI);
+        tmp->x = R * sin(coord->d);
+        tmp->y = R * cos(coord->d);
+
+    } else if (projection->type == PS_PROJ_SIN) {
+        R = cos(coord->r) * (180.0 / M_PI);
+        tmp->x = R * sin(coord->d);
+        tmp->y = R * cos(coord->d);
+
+    } else if (projection->type == PS_PROJ_CAR) {
+        tmp->x = coord->d;
+        tmp->y = coord->r;
+
+    } else if (projection->type == PS_PROJ_MER) {
+        tmp->x = coord->d;
+        tmp->y = log(tan(45.0 + (0.5 * coord->r))) * 180.0/M_PI;
+
+    } else if (projection->type == PS_PROJ_AIT) {
+        alpha = 1.0 / ((180.0 / M_PI) *
+                       sqrt(1.0 + (cos(coord->r) * cos(0.5 * coord->d) * 0.5)));
+
+        tmp->x = 2.0 * alpha * cos(coord->r) * sin(0.5 * coord->d);
+        tmp->y = alpha * sin(coord->d);
+
+    } else if (projection->type == PS_PROJ_PAR) {
+        psAbort(__func__, "The projection type PS_PROJ_PAR is undefined.\n");
+
+    } else if (projection->type == PS_PROJ_GLS) {
+        psAbort(__func__, "The projection type PS_PROJ_GLG is undefined.\n");
+    }
+
+    return(tmp);
+}
+
+psSphere *psDeproject(const psPlane *coord,
+                      const psProjection *projection)
+{
+    float R = 0.0;
+    float chu = 0.0;
+    float chu1 = 0.0;
+    float chu2 = 0.0;
+    psSphere *tmp= (psSphere *) psAlloc(sizeof(psSphere));
+
+    if (projection->type == PS_PROJ_TAN) {
+        R = sqrt((coord->x * coord->x) + (coord->y * coord->y));
+        tmp->d = arg(-coord->y, coord->x);
+        tmp->r = atan(180.0 / (R * M_PI));
+
+    } else if (projection->type == PS_PROJ_SIN) {
+        R = sqrt((coord->x * coord->x) + (coord->y * coord->y));
+        tmp->d = arg(-coord->y, coord->x);
+        tmp->r = acos((R * M_PI) / 180.0);
+
+    } else if (projection->type == PS_PROJ_CAR) {
+        tmp->d = coord->x;
+        tmp->r = coord->y;
+
+    } else if (projection->type == PS_PROJ_MER) {
+        tmp->d = coord->x;
+        tmp->r = (2.0 * atan(exp((coord->y * M_PI/180.0)))) - 180.0;
+
+    } else if (projection->type == PS_PROJ_AIT) {
+        chu1 = (coord->x * M_PI) / 720.0;
+        chu1*= chu1;
+        chu2 = (coord->y * M_PI) / 360.0;
+        chu2*= chu2;
+        chu = sqrt(1.0 - chu1 - chu2);
+        tmp->d = 2.0 * arg((2.0 * chu * chu) - 1.0,
+                           (coord->x * chu * M_PI)/360.0);
+        tmp->r = asin((coord->y * chu * M_PI) / 180.0);
+
+    } else if (projection->type == PS_PROJ_PAR) {
+        psAbort(__func__, "The projection type PS_PROJ_PAR is undefined.\n");
+
+    } else if (projection->type == PS_PROJ_GLS) {
+        psAbort(__func__, "The projection type PS_PROJ_GLG is undefined.\n");
+    }
+
+    return(tmp);
+}
 
 psSphere *psSphereGetOffset(const psSphere *restrict position1,
@@ -133,8 +248,11 @@
                             psSphereOffsetUnit unit)
 {
-    return(NULL);
-}
-
-
+    psAbort(__func__, "This function has not be dedfined.\n");
+    return(NULL);
+}
+
+
+// XXX: I copied the algorithm from the ADD exactly.  Arguments mode and unit
+// are ignored.
 psSphere *psSphereSetOffset(const psSphere *restrict position,
                             const psSphere *restrict offset,
@@ -142,4 +260,24 @@
                             psSphereOffsetUnit unit)
 {
-    return(NULL);
-}
+    psPlane lin;
+    psSphere *tmp;
+    psProjection proj;
+
+    if (mode == PS_LINEAR) {
+        proj.R = position->r;
+        proj.D = position->d;
+        proj.Xs = 0.0;
+        proj.Ys = 0.0;
+        proj.type = PS_PROJ_TAN;
+
+        lin.x = offset->r;
+        lin.y = offset->d;
+
+        tmp = psDeproject(&lin, &proj);
+        return(tmp);
+    } else if (mode == PS_SPHERICAL) {
+        psAbort(__func__, "Sperical offset modes are not defined.\n");
+    }
+    psAbort(__func__, "Unrecognized offset mode\n");
+    return(NULL);
+}