Index: trunk/psLib/src/astro/psCoord.c
===================================================================
--- trunk/psLib/src/astro/psCoord.c	(revision 1406)
+++ trunk/psLib/src/astro/psCoord.c	(revision 1407)
@@ -1,2 +1,3 @@
+
 /** @file  psCoord.c
 *
@@ -10,6 +11,6 @@
 *  @author George Gusciora, MHPCC
 *
-*  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2004-08-06 22:34:05 $
+*  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2004-08-07 00:06:06 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -27,68 +28,56 @@
 #include <float.h>
 
-static float cot( float x );
-static float arg( float x, float y );
+static float cot(float x);
+static float arg(float x, float y);
 
 // This is the only function in this file which I understand.
-psPlane *psPlaneTransformApply( psPlane *out,
-                                const psPlaneTransform *transform,
-                                const psPlane *coords )
-{
-    if ( out == NULL ) {
-        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 );
-
-    return ( out );
-}
-
+psPlane *psPlaneTransformApply(psPlane * out, const psPlaneTransform * transform, const psPlane * coords)
+{
+    if (out == NULL) {
+        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);
+
+    return (out);
+}
 
 // This transformation takes into account parameters beyond an objects
 // spatial coordinates: term3 and term4.
-psPlane *psPlaneDistortApply( psPlane *out,
-                              const psPlaneDistort *transform,
-                              const psPlane *coords,
-                              float term3,
-                              float term4 )
-{
-    if ( out == NULL ) {
-        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 ] * term3 ) +
-             ( transform->x->coeff[ 0 ][ 0 ][ 0 ][ 1 ] * term4 );
-
-    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 ] * term3 ) +
-             ( transform->y->coeff[ 0 ][ 0 ][ 0 ][ 1 ] * term4 );
-
-    return ( out );
-}
-
+psPlane *psPlaneDistortApply(psPlane * out,
+                             const psPlaneDistort * transform,
+                             const psPlane * coords, float term3, float term4)
+{
+    if (out == NULL) {
+        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] * term3) + (transform->x->coeff[0][0][0][1] * term4);
+
+    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] * term3) + (transform->y->coeff[0][0][0][1] * term4);
+
+    return (out);
+}
 
 // This function prototype has been modified since the SDRS.
-psSphereTransform *psSphereTransformAlloc( double NPlat,
-        double Xo,
-        double xo )
-{
-    psSphereTransform * tmp = ( psSphereTransform * ) psAlloc( sizeof( psSphereTransform ) );
-
-    tmp->sinPhi = sin( NPlat );
-    tmp->cosPhi = cos( NPlat );
+psSphereTransform *psSphereTransformAlloc(double NPlat, double Xo, double xo)
+{
+    psSphereTransform *tmp = (psSphereTransform *) psAlloc(sizeof(psSphereTransform));
+
+    tmp->sinPhi = sin(NPlat);
+    tmp->cosPhi = cos(NPlat);
     tmp->Xo = Xo;
     tmp->xo = xo;
 
-    return ( tmp );
+    return (tmp);
 }
 
@@ -100,7 +89,5 @@
 // there are no typo's.
 
-psSphere *psSphereTransformApply( psSphere *out,
-                                  const psSphereTransform *transform,
-                                  const psSphere *coord )
+psSphere *psSphereTransformApply(psSphere * out, const psSphereTransform * transform, const psSphere * coord)
 {
     double sinY = 0.0;
@@ -112,6 +99,6 @@
     double dx = 0.0;
 
-    if ( out == NULL ) {
-        out = ( psSphere * ) psAlloc( sizeof( psSphere ) );
+    if (out == NULL) {
+        out = (psSphere *) psAlloc(sizeof(psSphere));
     }
 
@@ -119,20 +106,19 @@
     y = coord->d;
     dx = x - transform->xo;
-    sinY = cos( y ) * sin( dx ) * transform->sinPhi + sin( y ) * transform->cosPhi;
-    cosY = sqrt( 1.0 - sinY * sinY );
-    sinX = ( cos( y ) * sin( dx ) * transform->cosPhi - sin( y ) * transform->sinPhi ) /
-           cos( y );
-    cosX = cos( y ) * cos( dx ) / cos( y );
-
-    out->r = atan2( sinX, cosX ) + transform->Xo;
-    out->d = atan2( sinY, cosY );
-
-    return ( out );
-}
-
-psSphereTransform *psSphereTransformICRStoEcliptic( psTime time )
-{
-    struct tm * tmTime = psTimeToTM( time );
-    double year = ( double ) ( 1900 + tmTime->tm_year );
+    sinY = cos(y) * sin(dx) * transform->sinPhi + sin(y) * transform->cosPhi;
+    cosY = sqrt(1.0 - sinY * sinY);
+    sinX = (cos(y) * sin(dx) * transform->cosPhi - sin(y) * transform->sinPhi) / cos(y);
+    cosX = cos(y) * cos(dx) / cos(y);
+
+    out->r = atan2(sinX, cosX) + transform->Xo;
+    out->d = atan2(sinY, cosY);
+
+    return (out);
+}
+
+psSphereTransform *psSphereTransformICRStoEcliptic(psTime time)
+{
+    struct tm *tmTime = psTimeToTM(time);
+    double year = (double)(1900 + tmTime->tm_year);
     double T = year / 100.0;
     double phi = -23.452294 + 0.013013 * T + 0.000001639 * T * T - 0.000000503 * T * T * T;
@@ -140,112 +126,99 @@
     double xo = 0.0;
 
-    return ( psSphereTransformAlloc( phi, Xo, xo ) );
-}
-
-psSphereTransform *psSphereTransformEcliptictoICRS( psTime time )
-{
-    struct tm * tmTime = psTimeToTM( time );
-    double year = ( double ) ( 1900 + tmTime->tm_year );
+    return (psSphereTransformAlloc(phi, Xo, xo));
+}
+
+psSphereTransform *psSphereTransformEcliptictoICRS(psTime time)
+{
+    struct tm *tmTime = psTimeToTM(time);
+    double year = (double)(1900 + tmTime->tm_year);
     double T = year / 100.0;
-    double phi = + 23.452294 - 0.013013 * T - 0.000001639 * T * T + 0.000000503 * T * T * T;
+    double phi = +23.452294 - 0.013013 * T - 0.000001639 * T * T + 0.000000503 * T * T * T;
     double Xo = 0.0;
     double xo = 0.0;
 
-    return ( psSphereTransformAlloc( phi, Xo, xo ) );
-}
-
-psSphereTransform *psSphereTransformICRStoGalatic( void )
-{
-    return ( psSphereTransformAlloc( 62.6, 282.25, 33.0 ) );
-}
-
-psSphereTransform *psSphereTransformGalatictoICRS( void )
-{
-    return ( psSphereTransformAlloc( -62.6, 33.0, 282.25 ) );
+    return (psSphereTransformAlloc(phi, Xo, xo));
+}
+
+psSphereTransform *psSphereTransformICRStoGalatic(void)
+{
+    return (psSphereTransformAlloc(62.6, 282.25, 33.0));
+}
+
+psSphereTransform *psSphereTransformGalatictoICRS(void)
+{
+    return (psSphereTransformAlloc(-62.6, 33.0, 282.25));
 }
 
 // XXX: Is this the correct way to calculate this?
-float cot( float x )
-{
-    return ( 1.0 / atan( x ) );
+float cot(float x)
+{
+    return (1.0 / atan(x));
 }
 
 // This is some kind of arc tan function.
-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 );
+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);
 }
 
 // XXX: Waiting for the definition of the PS_PROJ_PAR projection.
 // XXX: Waiting for the definition of the PS_PROJ_GLS projection.
-psPlane *psProject( const psSphere *coord,
-                    const psProjection *projection )
+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_GLS is undefined.\n" );
-                            }
-
-    return ( tmp );
+    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_GLS is undefined.\n");
+    }
+
+    return (tmp);
 }
 
 // XXX: Waiting for the definition of the PS_PROJ_PAR projection.
 // XXX: Waiting for the definition of the PS_PROJ_GLS projection.
-psSphere *psDeproject( const psPlane *coord,
-                       const psProjection *projection )
+psSphere *psDeproject(const psPlane * coord, const psProjection * projection)
 {
     float R = 0.0;
@@ -253,58 +226,50 @@
     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 *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);
 }
 
 // XXX: Do I need to check for unacceptable transformation parameters?
 // Maybe, if the points are on the North/South Pole, etc?
-psSphere *psSphereGetOffset( const psSphere *restrict position1,
-                             const psSphere *restrict position2,
-                             psSphereOffsetMode mode,
-                             psSphereOffsetUnit unit )
-{
-    //    psPlane *lin;
+psSphere *psSphereGetOffset(const psSphere * restrict position1,
+                            const psSphere * restrict position2,
+                            psSphereOffsetMode mode, psSphereOffsetUnit unit)
+{
+    // psPlane *lin;
     psProjection proj;
     psSphere *tmp;
@@ -312,5 +277,5 @@
     double tmpD = 0.0;
 
-    if ( mode == PS_LINEAR ) {
+    if (mode == PS_LINEAR) {
         // XXX: I have no idea how to construct this.  Maybe project both
         // sperical positions onto the plane, set the origin at one of the
@@ -318,5 +283,5 @@
 
         // XXX: Do I need to somehow scale this projection?
-        // project position1?  Will it project to (0.0, 0.0)?
+        // project position1? Will it project to (0.0, 0.0)?
         proj.R = position1->r;
         proj.D = position1->d;
@@ -325,51 +290,45 @@
         proj.type = PS_PROJ_TAN;
 
-        //        lin = psProject(position2, proj);
-        //        tmp = psDeproject(lin, proj);
+        // lin = psProject(position2, proj);
+        // tmp = psDeproject(lin, proj);
 
         // XXX: Do we need to convert units in tmp?
-        return ( tmp );
-    } else
-        if ( mode == PS_SPHERICAL ) {
-            tmpR = position2->r - position1->r;
-            tmpD = position2->d - position1->d;
-
-            if ( unit == PS_ARCSEC ) {
-                tmpR = ( tmpR * 180.0 * 60.0 * 60.0 ) / M_PI;
-                tmpD = ( tmpR * 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;
-                } else
-                    if ( unit == PS_DEGREE ) {
-                        tmpR = ( tmpR * 180.0 ) / M_PI;
-                        tmpD = ( tmpR * 180.0 ) / M_PI;
-                    } else
-                    if ( unit == PS_RADIAN ) {}
-                        else {
-                            psAbort( __func__, "Unknown offset unit: 0x%x\n", unit );
-                        }
-
-            tmp = ( psSphere * ) psAlloc( sizeof( psSphere ) );
-            tmp->r = tmpR;
-            tmp->d = tmpD;
-            tmp->rErr = 0.0;
-            tmp->dErr = 0.0;
-            // XXX: Do we need to wrap these to an acceptable range?
-            return ( tmp );
+        return (tmp);
+    } else if (mode == PS_SPHERICAL) {
+        tmpR = position2->r - position1->r;
+        tmpD = position2->d - position1->d;
+
+        if (unit == PS_ARCSEC) {
+            tmpR = (tmpR * 180.0 * 60.0 * 60.0) / M_PI;
+            tmpD = (tmpR * 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;
+        } else if (unit == PS_DEGREE) {
+            tmpR = (tmpR * 180.0) / M_PI;
+            tmpD = (tmpR * 180.0) / M_PI;
+        } else if (unit == PS_RADIAN) {}
+        else {
+            psAbort(__func__, "Unknown offset unit: 0x%x\n", unit);
         }
-    psAbort( __func__, "Unrecognized offset mode\n" );
-    return ( NULL );
-}
-
+
+        tmp = (psSphere *) psAlloc(sizeof(psSphere));
+        tmp->r = tmpR;
+        tmp->d = tmpD;
+        tmp->rErr = 0.0;
+        tmp->dErr = 0.0;
+        // XXX: Do we need to wrap these to an acceptable range?
+        return (tmp);
+    }
+    psAbort(__func__, "Unrecognized offset mode\n");
+    return (NULL);
+}
 
 // XXX: Do I need to check for unacceptable transformation parameters?
 // Maybe, if the points are on the North/South Pole, etc?
 // XXX: I copied the algorithm from the ADD exactly.
-psSphere *psSphereSetOffset( const psSphere *restrict position,
-                             const psSphere *restrict offset,
-                             psSphereOffsetMode mode,
-                             psSphereOffsetUnit unit )
+psSphere *psSphereSetOffset(const psSphere * restrict position,
+                            const psSphere * restrict offset,
+                            psSphereOffsetMode mode, psSphereOffsetUnit unit)
 {
     psPlane lin;
@@ -379,5 +338,5 @@
     double tmpD = 0.0;
 
-    if ( mode == PS_LINEAR ) {
+    if (mode == PS_LINEAR) {
         proj.R = position->r;
         proj.D = position->d;
@@ -389,39 +348,35 @@
         lin.y = offset->d;
 
-        tmp = psDeproject( &lin, &proj );
-        return ( tmp );
-
-    } else
-        if ( mode == PS_SPHERICAL ) {
-            if ( unit == PS_ARCSEC ) {
-                tmpR = ( M_PI * offset->r ) / ( 180.0 * 60.0 * 60.0 );
-                tmpD = ( M_PI * offset->d ) / ( 180.0 * 60.0 * 60.0 );
-            } else
-                if ( unit == PS_ARCMIN ) {
-                    tmpR = ( M_PI * offset->r ) / ( 180.0 * 60.0 );
-                    tmpD = ( M_PI * offset->d ) / ( 180.0 * 60.0 );
-                } else
-                    if ( unit == PS_DEGREE ) {
-                        tmpR = ( M_PI * offset->r ) / ( 180.0 );
-                        tmpD = ( M_PI * offset->d ) / ( 180.0 );
-                    } else
-                        if ( unit == PS_RADIAN ) {
-                            tmpR = offset->r;
-                            tmpD = offset->d;
-                        } else {
-                            psAbort( __func__, "Unknown offset unit: 0x%x\n", unit );
-                        }
-
-            tmp = ( psSphere * ) psAlloc( sizeof( psSphere ) );
-            tmp->r = position->r + tmpR;
-            tmp->r = position->d + tmpD;
-            tmp->rErr = 0.0;
-            tmp->dErr = 0.0;
-
-            // XXX: wrap tmp->r and tmp->d to the allowed range (-PI to PI)
-            // and (0 to 2*PI).
-            return ( tmp );
+        tmp = psDeproject(&lin, &proj);
+        return (tmp);
+
+    } else if (mode == PS_SPHERICAL) {
+        if (unit == PS_ARCSEC) {
+            tmpR = (M_PI * offset->r) / (180.0 * 60.0 * 60.0);
+            tmpD = (M_PI * offset->d) / (180.0 * 60.0 * 60.0);
+        } else if (unit == PS_ARCMIN) {
+            tmpR = (M_PI * offset->r) / (180.0 * 60.0);
+            tmpD = (M_PI * offset->d) / (180.0 * 60.0);
+        } else if (unit == PS_DEGREE) {
+            tmpR = (M_PI * offset->r) / (180.0);
+            tmpD = (M_PI * offset->d) / (180.0);
+        } else if (unit == PS_RADIAN) {
+            tmpR = offset->r;
+            tmpD = offset->d;
+        } else {
+            psAbort(__func__, "Unknown offset unit: 0x%x\n", unit);
         }
-    psAbort( __func__, "Unrecognized offset mode\n" );
-    return ( NULL );
-}
+
+        tmp = (psSphere *) psAlloc(sizeof(psSphere));
+        tmp->r = position->r + tmpR;
+        tmp->r = position->d + tmpD;
+        tmp->rErr = 0.0;
+        tmp->dErr = 0.0;
+
+        // XXX: wrap tmp->r and tmp->d to the allowed range (-PI to PI)
+        // and (0 to 2*PI).
+        return (tmp);
+    }
+    psAbort(__func__, "Unrecognized offset mode\n");
+    return (NULL);
+}