Index: trunk/psLib/src/astro/psCoord.c =================================================================== --- trunk/psLib/src/astro/psCoord.c (revision 1325) +++ trunk/psLib/src/astro/psCoord.c (revision 1374) @@ -1,18 +1,18 @@ /** @file psCoord.c - * - * @brief Contains basic coordinate transformation definitions and operations - * - * This file defines the basic types for astronomical coordinate - * transformation - * - * @ingroup AstroImage - * - * @author George Gusciora, MHPCC - * - * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-07-29 02:08:02 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - */ +* +* @brief Contains basic coordinate transformation definitions and operations +* +* This file defines the basic types for astronomical coordinate +* transformation +* +* @ingroup CoordinateTransform +* +* @author George Gusciora, MHPCC +* +* @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ +* @date $Date: 2004-08-04 00:55:17 $ +* +* Copyright 2004 Maui High Performance Computing Center, University of Hawaii +*/ #include "psType.h" @@ -28,20 +28,20 @@ // 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 ); } @@ -49,49 +49,49 @@ // 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, +psSphereTransform *psSphereTransformAlloc( double NPlat, double Xo, - double xo) -{ - psSphereTransform *tmp = (psSphereTransform *) psAlloc(sizeof(psSphereTransform)); - - tmp->sinPhi = sin(NPlat); - tmp->cosPhi = cos(NPlat); + 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 ); } // I understand this one too. -void p_psSphereTransformFree(psSphereTransform *trans) -{ - psFree(trans); +void p_psSphereTransformFree( psSphereTransform *trans ) +{ + psFree( trans ); } @@ -103,7 +103,7 @@ // 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; @@ -114,131 +114,131 @@ double y = 0.0; double dx = 0.0; - - if (out == NULL) { - out = (psSphere *) psAlloc(sizeof(psSphere)); - } - + + if ( out == NULL ) { + out = ( psSphere * ) psAlloc( sizeof( psSphere ) ); + } + x = coord->r; 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); + 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 ); 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 *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; @@ -246,50 +246,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) +psSphere *psSphereGetOffset( const psSphere *restrict position1, + const psSphere *restrict position2, + psSphereOffsetMode mode, + psSphereOffsetUnit unit ) { // psPlane *lin; @@ -298,51 +298,51 @@ double tmpR = 0.0; double tmpD = 0.0; - - 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 - // points on the plane, then deproject? - - // XXX: Do I need to somehow scale this projection? - // project position1? Will it project to (0.0, 0.0)? - proj.R = position1->r; - proj.D = position1->d; - proj.Xs = 1.0; - proj.Ys = 1.0; - proj.type = PS_PROJ_TAN; - - // 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); - } - psAbort(__func__, "Unrecognized offset mode\n"); - return(NULL); + + 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 + // points on the plane, then deproject? + + // XXX: Do I need to somehow scale this projection? + // project position1? Will it project to (0.0, 0.0)? + proj.R = position1->r; + proj.D = position1->d; + proj.Xs = 1.0; + proj.Ys = 1.0; + proj.type = PS_PROJ_TAN; + + // 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 ); + } + psAbort( __func__, "Unrecognized offset mode\n" ); + return ( NULL ); } @@ -351,8 +351,8 @@ // 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; @@ -361,46 +361,46 @@ double tmpR = 0.0; double tmpD = 0.0; - - if (mode == PS_LINEAR) { - proj.R = position->r; - proj.D = position->d; - proj.Xs = 1.0; - proj.Ys = 1.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) { - 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); - } - psAbort(__func__, "Unrecognized offset mode\n"); - return(NULL); -} + + if ( mode == PS_LINEAR ) { + proj.R = position->r; + proj.D = position->d; + proj.Xs = 1.0; + proj.Ys = 1.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 ) { + 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 ); + } + psAbort( __func__, "Unrecognized offset mode\n" ); + return ( NULL ); +}