Index: trunk/psLib/src/astro/psCoord.c
===================================================================
--- trunk/psLib/src/astro/psCoord.c	(revision 2058)
+++ trunk/psLib/src/astro/psCoord.c	(revision 2077)
@@ -10,6 +10,6 @@
 *  @author George Gusciora, MHPCC
 *
-*  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2004-10-12 20:53:02 $
+*  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2004-10-13 19:40:42 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -61,4 +61,5 @@
 XXX: Verify this arc tan function.
 XXX: macro
+XXX: I don't think this arg() is correct.
  *****************************************************************************/
 
@@ -66,18 +67,24 @@
               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 (p_psArg(%f, %f).\n", x, y);
+    if (y!= 0.0) {
+        return(atan(y/x));
+    } else {
+        return(0.5 * M_PI);
+    }
     return (0.0);
+    /*
+        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 (p_psArg(%f, %f).\n", x, y);
+    */
 }
 
@@ -92,12 +99,6 @@
         out = (psPlane* ) psAlloc(sizeof(psPlane));
     }
-    out->x = transform->x->coeff[0][0] +
-             (transform->x->coeff[1][0] * coords->x) +
-             (transform->x->coeff[0][1] * coords->y);
-
-    out->y = transform->y->coeff[0][0] +
-             (transform->y->coeff[1][0] * coords->x) +
-             (transform->y->coeff[0][1] * coords->y);
-
+    out->x = psDPolynomial2DEval(coords->x, coords->y, transform->x);
+    out->y = psDPolynomial2DEval(coords->x, coords->y, transform->y);
     return (out);
 }
@@ -116,17 +117,6 @@
         out = (psPlane* ) psAlloc(sizeof(psPlane));
     }
-
-    out->x = transform->x->coeff[0][0][0][0] +
-             (transform->x->coeff[1][0][0][0] * coords->x) +
-             (transform->x->coeff[0][1][0][0] * coords->y) +
-             (transform->x->coeff[0][0][1][0] * color) +
-             (transform->x->coeff[0][0][0][1] * magnitude);
-
-    out->y = transform->y->coeff[0][0][0][0] +
-             (transform->y->coeff[1][0][0][0] * coords->x) +
-             (transform->y->coeff[0][1][0][0] * coords->y) +
-             (transform->y->coeff[0][0][1][0] * color) +
-             (transform->y->coeff[0][0][0][1] * magnitude);
-
+    out->x = psDPolynomial4DEval(coords->x, coords->y, color, magnitude, transform->x);
+    out->y = psDPolynomial4DEval(coords->x, coords->y, color, magnitude, transform->y);
     return (out);
 }
@@ -196,4 +186,5 @@
     double xo = 0.0;
 
+    psFree(tmTime);
     return (psSphereTransformAlloc(phi, Xo, xo));
 }
@@ -208,4 +199,5 @@
     double xo = 0.0;
 
+    psFree(tmTime);
     return (psSphereTransformAlloc(phi, Xo, xo));
 }
@@ -251,5 +243,6 @@
 
     } 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)));
+        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);
@@ -325,4 +318,9 @@
  
 XXX: Do I need to somehow scale this projection?
+ 
+XXX: Does PS_LINEAR mode make sense?  The result must be returned in psSphere
+     regardless of the mode.
+ 
+XXX: How to compound errors?
  *****************************************************************************/
 psSphere* psSphereGetOffset(const psSphere* restrict position1,
@@ -351,16 +349,24 @@
         return (tmp);
     } else if (mode == PS_SPHERICAL) {
-        tmpR = position2->r - position1->r;
-        tmpD = position2->d - position1->d;
+        tmp = (psSphere* ) psAlloc(sizeof(psSphere));
+        tmp->r = position2->r - position1->r;
+        tmp->d = position2->d - position1->d;
+
+        // XXX: Wrap these to an acceptable range.  This assumes that all
+        // angles are in radians.
+        tmp->r = fmod(tmp->r, 2*M_PI);
+        tmp->d = fmod(tmp->d, 2*M_PI);
+        tmp->rErr = 0.0;
+        tmp->dErr = 0.0;
 
         if (unit == PS_ARCSEC) {
-            tmpR = (tmpR * 180.0 * 60.0 * 60.0) / M_PI;
-            tmpD = (tmpR * 180.0 * 60.0 * 60.0) / M_PI;
+            tmp->r = (tmp->r * 180.0 * 60.0 * 60.0) / M_PI;
+            tmp->d = (tmp->d * 180.0 * 60.0 * 60.0) / M_PI;
         } else if (unit == PS_ARCMIN) {
-            tmpR = (tmpR * 180.0 * 60.0) / M_PI;
-            tmpD = (tmpR * 180.0 * 60.0) / M_PI;
+            tmp->r = (tmp->r * 180.0 * 60.0) / M_PI;
+            tmp->d = (tmp->d * 180.0 * 60.0) / M_PI;
         } else if (unit == PS_DEGREE) {
-            tmpR = (tmpR * 180.0) / M_PI;
-            tmpD = (tmpR * 180.0) / M_PI;
+            tmp->r = (tmp->r * 180.0) / M_PI;
+            tmp->d = (tmp->d * 180.0) / M_PI;
         } else if (unit == PS_RADIAN) {}
         else {
@@ -368,11 +374,4 @@
         }
 
-        tmp = (psSphere* ) psAlloc(sizeof(psSphere));
-        // XXX: Wrap these to an acceptable range.  Is this correct?
-        tmp->r = fmod(tmpR, 2*M_PI);
-        tmp->d = fmod(tmpD, 2*M_PI);
-        tmp->rErr = 0.0;
-        tmp->dErr = 0.0;
-
         return (tmp);
     }
@@ -388,4 +387,6 @@
  
 XXX: I copied the algorithm from the ADD exactly.
+ 
+XXX: Should we compound errors?
  *****************************************************************************/
 
@@ -432,5 +433,4 @@
 
         tmp = (psSphere* ) psAlloc(sizeof(psSphere));
-        // XXX: Is this an acceptable way to wrap the angular output?
         tmp->r = position->r + tmpR;
         tmp->r = fmod(tmp->r, 2.0*M_PI);