IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 6, 2004, 2:06:06 PM (22 years ago)
Author:
desonia
Message:

another attempt to get astyle to get it right.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/astronomy/psCoord.c

    r1406 r1407  
     1
    12/** @file  psCoord.c
    23*
    … …  
    1011*  @author George Gusciora, MHPCC
    1112*
    12 *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2004-08-06 22:34:05 $
     13*  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
     14*  @date $Date: 2004-08-07 00:06:06 $
    1415*
    1516*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2728#include <float.h>
    2829
    29 static float cot( float x );
    30 static float arg( float x, float y );
     30static float cot(float x);
     31static float arg(float x, float y);
    3132
    3233// This is the only function in this file which I understand.
    33 psPlane *psPlaneTransformApply( psPlane *out,
    34                                 const psPlaneTransform *transform,
    35                                 const psPlane *coords )
    36 {
    37     if ( out == NULL ) {
    38         out = ( psPlane * ) psAlloc( sizeof( psPlane ) );
    39     }
    40     out->x = transform->x->coeff[ 0 ][ 0 ] +
    41              ( transform->x->coeff[ 1 ][ 0 ] * coords->x ) +
    42              ( transform->x->coeff[ 0 ][ 1 ] * coords->y );
    43 
    44     out->y = transform->y->coeff[ 0 ][ 0 ] +
    45              ( transform->y->coeff[ 1 ][ 0 ] * coords->x ) +
    46              ( transform->y->coeff[ 0 ][ 1 ] * coords->y );
    47 
    48     return ( out );
    49 }
    50 
     34psPlane *psPlaneTransformApply(psPlane * out, const psPlaneTransform * transform, const psPlane * coords)
     35{
     36    if (out == NULL) {
     37        out = (psPlane *) psAlloc(sizeof(psPlane));
     38    }
     39    out->x = transform->x->coeff[0][0] +
     40             (transform->x->coeff[1][0] * coords->x) + (transform->x->coeff[0][1] * coords->y);
     41
     42    out->y = transform->y->coeff[0][0] +
     43             (transform->y->coeff[1][0] * coords->x) + (transform->y->coeff[0][1] * coords->y);
     44
     45    return (out);
     46}
    5147
    5248// This transformation takes into account parameters beyond an objects
    5349// spatial coordinates: term3 and term4.
    54 psPlane *psPlaneDistortApply( psPlane *out,
    55                               const psPlaneDistort *transform,
    56                               const psPlane *coords,
    57                               float term3,
    58                               float term4 )
    59 {
    60     if ( out == NULL ) {
    61         out = ( psPlane * ) psAlloc( sizeof( psPlane ) );
    62     }
    63 
    64     out->x = transform->x->coeff[ 0 ][ 0 ][ 0 ][ 0 ] +
    65              ( transform->x->coeff[ 1 ][ 0 ][ 0 ][ 0 ] * coords->x ) +
    66              ( transform->x->coeff[ 0 ][ 1 ][ 0 ][ 0 ] * coords->y ) +
    67              ( transform->x->coeff[ 0 ][ 0 ][ 1 ][ 0 ] * term3 ) +
    68              ( transform->x->coeff[ 0 ][ 0 ][ 0 ][ 1 ] * term4 );
    69 
    70     out->y = transform->y->coeff[ 0 ][ 0 ][ 0 ][ 0 ] +
    71              ( transform->y->coeff[ 1 ][ 0 ][ 0 ][ 0 ] * coords->x ) +
    72              ( transform->y->coeff[ 0 ][ 1 ][ 0 ][ 0 ] * coords->y ) +
    73              ( transform->y->coeff[ 0 ][ 0 ][ 1 ][ 0 ] * term3 ) +
    74              ( transform->y->coeff[ 0 ][ 0 ][ 0 ][ 1 ] * term4 );
    75 
    76     return ( out );
    77 }
    78 
     50psPlane *psPlaneDistortApply(psPlane * out,
     51                             const psPlaneDistort * transform,
     52                             const psPlane * coords, float term3, float term4)
     53{
     54    if (out == NULL) {
     55        out = (psPlane *) psAlloc(sizeof(psPlane));
     56    }
     57
     58    out->x = transform->x->coeff[0][0][0][0] +
     59             (transform->x->coeff[1][0][0][0] * coords->x) +
     60             (transform->x->coeff[0][1][0][0] * coords->y) +
     61             (transform->x->coeff[0][0][1][0] * term3) + (transform->x->coeff[0][0][0][1] * term4);
     62
     63    out->y = transform->y->coeff[0][0][0][0] +
     64             (transform->y->coeff[1][0][0][0] * coords->x) +
     65             (transform->y->coeff[0][1][0][0] * coords->y) +
     66             (transform->y->coeff[0][0][1][0] * term3) + (transform->y->coeff[0][0][0][1] * term4);
     67
     68    return (out);
     69}
    7970
    8071// This function prototype has been modified since the SDRS.
    81 psSphereTransform *psSphereTransformAlloc( double NPlat,
    82         double Xo,
    83         double xo )
    84 {
    85     psSphereTransform * tmp = ( psSphereTransform * ) psAlloc( sizeof( psSphereTransform ) );
    86 
    87     tmp->sinPhi = sin( NPlat );
    88     tmp->cosPhi = cos( NPlat );
     72psSphereTransform *psSphereTransformAlloc(double NPlat, double Xo, double xo)
     73{
     74    psSphereTransform *tmp = (psSphereTransform *) psAlloc(sizeof(psSphereTransform));
     75
     76    tmp->sinPhi = sin(NPlat);
     77    tmp->cosPhi = cos(NPlat);
    8978    tmp->Xo = Xo;
    9079    tmp->xo = xo;
    9180
    92     return ( tmp );
     81    return (tmp);
    9382}
    9483
    … …  
    10089// there are no typo's.
    10190
    102 psSphere *psSphereTransformApply( psSphere *out,
    103                                   const psSphereTransform *transform,
    104                                   const psSphere *coord )
     91psSphere *psSphereTransformApply(psSphere * out, const psSphereTransform * transform, const psSphere * coord)
    10592{
    10693    double sinY = 0.0;
    … …  
    11299    double dx = 0.0;
    113100
    114     if ( out == NULL ) {
    115         out = ( psSphere * ) psAlloc( sizeof( psSphere ) );
     101    if (out == NULL) {
     102        out = (psSphere *) psAlloc(sizeof(psSphere));
    116103    }
    117104
    … …  
    119106    y = coord->d;
    120107    dx = x - transform->xo;
    121     sinY = cos( y ) * sin( dx ) * transform->sinPhi + sin( y ) * transform->cosPhi;
    122     cosY = sqrt( 1.0 - sinY * sinY );
    123     sinX = ( cos( y ) * sin( dx ) * transform->cosPhi - sin( y ) * transform->sinPhi ) /
    124            cos( y );
    125     cosX = cos( y ) * cos( dx ) / cos( y );
    126 
    127     out->r = atan2( sinX, cosX ) + transform->Xo;
    128     out->d = atan2( sinY, cosY );
    129 
    130     return ( out );
    131 }
    132 
    133 psSphereTransform *psSphereTransformICRStoEcliptic( psTime time )
    134 {
    135     struct tm * tmTime = psTimeToTM( time );
    136     double year = ( double ) ( 1900 + tmTime->tm_year );
     108    sinY = cos(y) * sin(dx) * transform->sinPhi + sin(y) * transform->cosPhi;
     109    cosY = sqrt(1.0 - sinY * sinY);
     110    sinX = (cos(y) * sin(dx) * transform->cosPhi - sin(y) * transform->sinPhi) / cos(y);
     111    cosX = cos(y) * cos(dx) / cos(y);
     112
     113    out->r = atan2(sinX, cosX) + transform->Xo;
     114    out->d = atan2(sinY, cosY);
     115
     116    return (out);
     117}
     118
     119psSphereTransform *psSphereTransformICRStoEcliptic(psTime time)
     120{
     121    struct tm *tmTime = psTimeToTM(time);
     122    double year = (double)(1900 + tmTime->tm_year);
    137123    double T = year / 100.0;
    138124    double phi = -23.452294 + 0.013013 * T + 0.000001639 * T * T - 0.000000503 * T * T * T;
    … …  
    140126    double xo = 0.0;
    141127
    142     return ( psSphereTransformAlloc( phi, Xo, xo ) );
    143 }
    144 
    145 psSphereTransform *psSphereTransformEcliptictoICRS( psTime time )
    146 {
    147     struct tm * tmTime = psTimeToTM( time );
    148     double year = ( double ) ( 1900 + tmTime->tm_year );
     128    return (psSphereTransformAlloc(phi, Xo, xo));
     129}
     130
     131psSphereTransform *psSphereTransformEcliptictoICRS(psTime time)
     132{
     133    struct tm *tmTime = psTimeToTM(time);
     134    double year = (double)(1900 + tmTime->tm_year);
    149135    double T = year / 100.0;
    150     double phi = + 23.452294 - 0.013013 * T - 0.000001639 * T * T + 0.000000503 * T * T * T;
     136    double phi = +23.452294 - 0.013013 * T - 0.000001639 * T * T + 0.000000503 * T * T * T;
    151137    double Xo = 0.0;
    152138    double xo = 0.0;
    153139
    154     return ( psSphereTransformAlloc( phi, Xo, xo ) );
    155 }
    156 
    157 psSphereTransform *psSphereTransformICRStoGalatic( void )
    158 {
    159     return ( psSphereTransformAlloc( 62.6, 282.25, 33.0 ) );
    160 }
    161 
    162 psSphereTransform *psSphereTransformGalatictoICRS( void )
    163 {
    164     return ( psSphereTransformAlloc( -62.6, 33.0, 282.25 ) );
     140    return (psSphereTransformAlloc(phi, Xo, xo));
     141}
     142
     143psSphereTransform *psSphereTransformICRStoGalatic(void)
     144{
     145    return (psSphereTransformAlloc(62.6, 282.25, 33.0));
     146}
     147
     148psSphereTransform *psSphereTransformGalatictoICRS(void)
     149{
     150    return (psSphereTransformAlloc(-62.6, 33.0, 282.25));
    165151}
    166152
    167153// XXX: Is this the correct way to calculate this?
    168 float cot( float x )
    169 {
    170     return ( 1.0 / atan( x ) );
     154float cot(float x)
     155{
     156    return (1.0 / atan(x));
    171157}
    172158
    173159// This is some kind of arc tan function.
    174 float arg( float x, float y )
    175 {
    176     if ( x > 0 ) {
    177         return ( atan( y / x ) );
    178     } else
    179         if ( ( x == 0 ) && ( y == 0 ) ) {
    180             return ( 0.5 * M_PI );
    181         } else
    182             if ( ( x == 0 ) && ( y == 0 ) ) {
    183                 return ( -0.5 * M_PI );
    184             } else
    185                 if ( ( x == 0 ) && ( y == 0 ) ) {
    186                     return ( M_PI + atan( y / x ) );
    187                 } else
    188                     if ( ( x == 0 ) && ( y == 0 ) ) {
    189                         return ( -M_PI + atan( y / x ) );
    190                     }
    191 
    192     psAbort( __func__, "Unacceptable range for (arg(%f, %f).\n", x, y );
    193     return ( 0.0 );
     160float arg(float x, float y)
     161{
     162    if (x > 0) {
     163        return (atan(y / x));
     164    } else if ((x == 0) && (y == 0)) {
     165        return (0.5 * M_PI);
     166    } else if ((x == 0) && (y == 0)) {
     167        return (-0.5 * M_PI);
     168    } else if ((x == 0) && (y == 0)) {
     169        return (M_PI + atan(y / x));
     170    } else if ((x == 0) && (y == 0)) {
     171        return (-M_PI + atan(y / x));
     172    }
     173
     174    psAbort(__func__, "Unacceptable range for (arg(%f, %f).\n", x, y);
     175    return (0.0);
    194176}
    195177
    196178// XXX: Waiting for the definition of the PS_PROJ_PAR projection.
    197179// XXX: Waiting for the definition of the PS_PROJ_GLS projection.
    198 psPlane *psProject( const psSphere *coord,
    199                     const psProjection *projection )
     180psPlane *psProject(const psSphere * coord, const psProjection * projection)
    200181{
    201182    float R = 0.0;
    202183    float alpha = 0.0;
    203     psPlane *tmp = ( psPlane * ) psAlloc( sizeof( psPlane ) );
    204 
    205     if ( projection->type == PS_PROJ_TAN ) {
    206         R = cot( coord->r ) * ( 180.0 / M_PI );
    207         tmp->x = R * sin( coord->d );
    208         tmp->y = R * cos( coord->d );
    209 
    210     } else
    211         if ( projection->type == PS_PROJ_SIN ) {
    212             R = cos( coord->r ) * ( 180.0 / M_PI );
    213             tmp->x = R * sin( coord->d );
    214             tmp->y = R * cos( coord->d );
    215 
    216         } else
    217             if ( projection->type == PS_PROJ_CAR ) {
    218                 tmp->x = coord->d;
    219                 tmp->y = coord->r;
    220 
    221             } else
    222                 if ( projection->type == PS_PROJ_MER ) {
    223                     tmp->x = coord->d;
    224                     tmp->y = log( tan( 45.0 + ( 0.5 * coord->r ) ) ) * 180.0 / M_PI;
    225 
    226                 } else
    227                     if ( projection->type == PS_PROJ_AIT ) {
    228                         alpha = 1.0 / ( ( 180.0 / M_PI ) *
    229                                         sqrt( 1.0 + ( cos( coord->r ) * cos( 0.5 * coord->d ) * 0.5 ) ) );
    230 
    231                         tmp->x = 2.0 * alpha * cos( coord->r ) * sin( 0.5 * coord->d );
    232                         tmp->y = alpha * sin( coord->d );
    233 
    234                     } else
    235                         if ( projection->type == PS_PROJ_PAR ) {
    236                             psAbort( __func__, "The projection type PS_PROJ_PAR is undefined.\n" );
    237 
    238                         } else
    239                             if ( projection->type == PS_PROJ_GLS ) {
    240                                 psAbort( __func__, "The projection type PS_PROJ_GLS is undefined.\n" );
    241                             }
    242 
    243     return ( tmp );
     184    psPlane *tmp = (psPlane *) psAlloc(sizeof(psPlane));
     185
     186    if (projection->type == PS_PROJ_TAN) {
     187        R = cot(coord->r) * (180.0 / M_PI);
     188        tmp->x = R * sin(coord->d);
     189        tmp->y = R * cos(coord->d);
     190
     191    } else if (projection->type == PS_PROJ_SIN) {
     192        R = cos(coord->r) * (180.0 / M_PI);
     193        tmp->x = R * sin(coord->d);
     194        tmp->y = R * cos(coord->d);
     195
     196    } else if (projection->type == PS_PROJ_CAR) {
     197        tmp->x = coord->d;
     198        tmp->y = coord->r;
     199
     200    } else if (projection->type == PS_PROJ_MER) {
     201        tmp->x = coord->d;
     202        tmp->y = log(tan(45.0 + (0.5 * coord->r))) * 180.0 / M_PI;
     203
     204    } else if (projection->type == PS_PROJ_AIT) {
     205        alpha = 1.0 / ((180.0 / M_PI) * sqrt(1.0 + (cos(coord->r) * cos(0.5 * coord->d) * 0.5)));
     206
     207        tmp->x = 2.0 * alpha * cos(coord->r) * sin(0.5 * coord->d);
     208        tmp->y = alpha * sin(coord->d);
     209
     210    } else if (projection->type == PS_PROJ_PAR) {
     211        psAbort(__func__, "The projection type PS_PROJ_PAR is undefined.\n");
     212
     213    } else if (projection->type == PS_PROJ_GLS) {
     214        psAbort(__func__, "The projection type PS_PROJ_GLS is undefined.\n");
     215    }
     216
     217    return (tmp);
    244218}
    245219
    246220// XXX: Waiting for the definition of the PS_PROJ_PAR projection.
    247221// XXX: Waiting for the definition of the PS_PROJ_GLS projection.
    248 psSphere *psDeproject( const psPlane *coord,
    249                        const psProjection *projection )
     222psSphere *psDeproject(const psPlane * coord, const psProjection * projection)
    250223{
    251224    float R = 0.0;
    … …  
    253226    float chu1 = 0.0;
    254227    float chu2 = 0.0;
    255     psSphere *tmp = ( psSphere * ) psAlloc( sizeof( psSphere ) );
    256 
    257     if ( projection->type == PS_PROJ_TAN ) {
    258         R = sqrt( ( coord->x * coord->x ) + ( coord->y * coord->y ) );
    259         tmp->d = arg( -coord->y, coord->x );
    260         tmp->r = atan( 180.0 / ( R * M_PI ) );
    261 
    262     } else
    263         if ( projection->type == PS_PROJ_SIN ) {
    264             R = sqrt( ( coord->x * coord->x ) + ( coord->y * coord->y ) );
    265             tmp->d = arg( -coord->y, coord->x );
    266             tmp->r = acos( ( R * M_PI ) / 180.0 );
    267 
    268         } else
    269             if ( projection->type == PS_PROJ_CAR ) {
    270                 tmp->d = coord->x;
    271                 tmp->r = coord->y;
    272 
    273             } else
    274                 if ( projection->type == PS_PROJ_MER ) {
    275                     tmp->d = coord->x;
    276                     tmp->r = ( 2.0 * atan( exp( ( coord->y * M_PI / 180.0 ) ) ) ) - 180.0;
    277 
    278                 } else
    279                     if ( projection->type == PS_PROJ_AIT ) {
    280                         chu1 = ( coord->x * M_PI ) / 720.0;
    281                         chu1 *= chu1;
    282                         chu2 = ( coord->y * M_PI ) / 360.0;
    283                         chu2 *= chu2;
    284                         chu = sqrt( 1.0 - chu1 - chu2 );
    285                         tmp->d = 2.0 * arg( ( 2.0 * chu * chu ) - 1.0,
    286                                             ( coord->x * chu * M_PI ) / 360.0 );
    287                         tmp->r = asin( ( coord->y * chu * M_PI ) / 180.0 );
    288 
    289                     } else
    290                         if ( projection->type == PS_PROJ_PAR ) {
    291                             psAbort( __func__, "The projection type PS_PROJ_PAR is undefined.\n" );
    292 
    293                         } else
    294                             if ( projection->type == PS_PROJ_GLS ) {
    295                                 psAbort( __func__, "The projection type PS_PROJ_GLG is undefined.\n" );
    296                             }
    297 
    298     return ( tmp );
     228    psSphere *tmp = (psSphere *) psAlloc(sizeof(psSphere));
     229
     230    if (projection->type == PS_PROJ_TAN) {
     231        R = sqrt((coord->x * coord->x) + (coord->y * coord->y));
     232        tmp->d = arg(-coord->y, coord->x);
     233        tmp->r = atan(180.0 / (R * M_PI));
     234
     235    } else if (projection->type == PS_PROJ_SIN) {
     236        R = sqrt((coord->x * coord->x) + (coord->y * coord->y));
     237        tmp->d = arg(-coord->y, coord->x);
     238        tmp->r = acos((R * M_PI) / 180.0);
     239
     240    } else if (projection->type == PS_PROJ_CAR) {
     241        tmp->d = coord->x;
     242        tmp->r = coord->y;
     243
     244    } else if (projection->type == PS_PROJ_MER) {
     245        tmp->d = coord->x;
     246        tmp->r = (2.0 * atan(exp((coord->y * M_PI / 180.0)))) - 180.0;
     247
     248    } else if (projection->type == PS_PROJ_AIT) {
     249        chu1 = (coord->x * M_PI) / 720.0;
     250        chu1 *= chu1;
     251        chu2 = (coord->y * M_PI) / 360.0;
     252        chu2 *= chu2;
     253        chu = sqrt(1.0 - chu1 - chu2);
     254        tmp->d = 2.0 * arg((2.0 * chu * chu) - 1.0, (coord->x * chu * M_PI) / 360.0);
     255        tmp->r = asin((coord->y * chu * M_PI) / 180.0);
     256
     257    } else if (projection->type == PS_PROJ_PAR) {
     258        psAbort(__func__, "The projection type PS_PROJ_PAR is undefined.\n");
     259
     260    } else if (projection->type == PS_PROJ_GLS) {
     261        psAbort(__func__, "The projection type PS_PROJ_GLG is undefined.\n");
     262    }
     263
     264    return (tmp);
    299265}
    300266
    301267// XXX: Do I need to check for unacceptable transformation parameters?
    302268// Maybe, if the points are on the North/South Pole, etc?
    303 psSphere *psSphereGetOffset( const psSphere *restrict position1,
    304                              const psSphere *restrict position2,
    305                              psSphereOffsetMode mode,
    306                              psSphereOffsetUnit unit )
    307 {
    308     //    psPlane *lin;
     269psSphere *psSphereGetOffset(const psSphere * restrict position1,
     270                            const psSphere * restrict position2,
     271                            psSphereOffsetMode mode, psSphereOffsetUnit unit)
     272{
     273    // psPlane *lin;
    309274    psProjection proj;
    310275    psSphere *tmp;
    … …  
    312277    double tmpD = 0.0;
    313278
    314     if ( mode == PS_LINEAR ) {
     279    if (mode == PS_LINEAR) {
    315280        // XXX: I have no idea how to construct this.  Maybe project both
    316281        // sperical positions onto the plane, set the origin at one of the
    … …  
    318283
    319284        // XXX: Do I need to somehow scale this projection?
    320         // project position1?  Will it project to (0.0, 0.0)?
     285        // project position1? Will it project to (0.0, 0.0)?
    321286        proj.R = position1->r;
    322287        proj.D = position1->d;
    … …  
    325290        proj.type = PS_PROJ_TAN;
    326291
    327         //        lin = psProject(position2, proj);
    328         //        tmp = psDeproject(lin, proj);
     292        // lin = psProject(position2, proj);
     293        // tmp = psDeproject(lin, proj);
    329294
    330295        // XXX: Do we need to convert units in tmp?
    331         return ( tmp );
    332     } else
    333         if ( mode == PS_SPHERICAL ) {
    334             tmpR = position2->r - position1->r;
    335             tmpD = position2->d - position1->d;
    336 
    337             if ( unit == PS_ARCSEC ) {
    338                 tmpR = ( tmpR * 180.0 * 60.0 * 60.0 ) / M_PI;
    339                 tmpD = ( tmpR * 180.0 * 60.0 * 60.0 ) / M_PI;
    340             } else
    341                 if ( unit == PS_ARCMIN ) {
    342                     tmpR = ( tmpR * 180.0 * 60.0 ) / M_PI;
    343                     tmpD = ( tmpR * 180.0 * 60.0 ) / M_PI;
    344                 } else
    345                     if ( unit == PS_DEGREE ) {
    346                         tmpR = ( tmpR * 180.0 ) / M_PI;
    347                         tmpD = ( tmpR * 180.0 ) / M_PI;
    348                     } else
    349                     if ( unit == PS_RADIAN ) {}
    350                         else {
    351                             psAbort( __func__, "Unknown offset unit: 0x%x\n", unit );
    352                         }
    353 
    354             tmp = ( psSphere * ) psAlloc( sizeof( psSphere ) );
    355             tmp->r = tmpR;
    356             tmp->d = tmpD;
    357             tmp->rErr = 0.0;
    358             tmp->dErr = 0.0;
    359             // XXX: Do we need to wrap these to an acceptable range?
    360             return ( tmp );
     296        return (tmp);
     297    } else if (mode == PS_SPHERICAL) {
     298        tmpR = position2->r - position1->r;
     299        tmpD = position2->d - position1->d;
     300
     301        if (unit == PS_ARCSEC) {
     302            tmpR = (tmpR * 180.0 * 60.0 * 60.0) / M_PI;
     303            tmpD = (tmpR * 180.0 * 60.0 * 60.0) / M_PI;
     304        } else if (unit == PS_ARCMIN) {
     305            tmpR = (tmpR * 180.0 * 60.0) / M_PI;
     306            tmpD = (tmpR * 180.0 * 60.0) / M_PI;
     307        } else if (unit == PS_DEGREE) {
     308            tmpR = (tmpR * 180.0) / M_PI;
     309            tmpD = (tmpR * 180.0) / M_PI;
     310        } else if (unit == PS_RADIAN) {}
     311        else {
     312            psAbort(__func__, "Unknown offset unit: 0x%x\n", unit);
    361313        }
    362     psAbort( __func__, "Unrecognized offset mode\n" );
    363     return ( NULL );
    364 }
    365 
     314
     315        tmp = (psSphere *) psAlloc(sizeof(psSphere));
     316        tmp->r = tmpR;
     317        tmp->d = tmpD;
     318        tmp->rErr = 0.0;
     319        tmp->dErr = 0.0;
     320        // XXX: Do we need to wrap these to an acceptable range?
     321        return (tmp);
     322    }
     323    psAbort(__func__, "Unrecognized offset mode\n");
     324    return (NULL);
     325}
    366326
    367327// XXX: Do I need to check for unacceptable transformation parameters?
    368328// Maybe, if the points are on the North/South Pole, etc?
    369329// XXX: I copied the algorithm from the ADD exactly.
    370 psSphere *psSphereSetOffset( const psSphere *restrict position,
    371                              const psSphere *restrict offset,
    372                              psSphereOffsetMode mode,
    373                              psSphereOffsetUnit unit )
     330psSphere *psSphereSetOffset(const psSphere * restrict position,
     331                            const psSphere * restrict offset,
     332                            psSphereOffsetMode mode, psSphereOffsetUnit unit)
    374333{
    375334    psPlane lin;
    … …  
    379338    double tmpD = 0.0;
    380339
    381     if ( mode == PS_LINEAR ) {
     340    if (mode == PS_LINEAR) {
    382341        proj.R = position->r;
    383342        proj.D = position->d;
    … …  
    389348        lin.y = offset->d;
    390349
    391         tmp = psDeproject( &lin, &proj );
    392         return ( tmp );
    393 
    394     } else
    395         if ( mode == PS_SPHERICAL ) {
    396             if ( unit == PS_ARCSEC ) {
    397                 tmpR = ( M_PI * offset->r ) / ( 180.0 * 60.0 * 60.0 );
    398                 tmpD = ( M_PI * offset->d ) / ( 180.0 * 60.0 * 60.0 );
    399             } else
    400                 if ( unit == PS_ARCMIN ) {
    401                     tmpR = ( M_PI * offset->r ) / ( 180.0 * 60.0 );
    402                     tmpD = ( M_PI * offset->d ) / ( 180.0 * 60.0 );
    403                 } else
    404                     if ( unit == PS_DEGREE ) {
    405                         tmpR = ( M_PI * offset->r ) / ( 180.0 );
    406                         tmpD = ( M_PI * offset->d ) / ( 180.0 );
    407                     } else
    408                         if ( unit == PS_RADIAN ) {
    409                             tmpR = offset->r;
    410                             tmpD = offset->d;
    411                         } else {
    412                             psAbort( __func__, "Unknown offset unit: 0x%x\n", unit );
    413                         }
    414 
    415             tmp = ( psSphere * ) psAlloc( sizeof( psSphere ) );
    416             tmp->r = position->r + tmpR;
    417             tmp->r = position->d + tmpD;
    418             tmp->rErr = 0.0;
    419             tmp->dErr = 0.0;
    420 
    421             // XXX: wrap tmp->r and tmp->d to the allowed range (-PI to PI)
    422             // and (0 to 2*PI).
    423             return ( tmp );
     350        tmp = psDeproject(&lin, &proj);
     351        return (tmp);
     352
     353    } else if (mode == PS_SPHERICAL) {
     354        if (unit == PS_ARCSEC) {
     355            tmpR = (M_PI * offset->r) / (180.0 * 60.0 * 60.0);
     356            tmpD = (M_PI * offset->d) / (180.0 * 60.0 * 60.0);
     357        } else if (unit == PS_ARCMIN) {
     358            tmpR = (M_PI * offset->r) / (180.0 * 60.0);
     359            tmpD = (M_PI * offset->d) / (180.0 * 60.0);
     360        } else if (unit == PS_DEGREE) {
     361            tmpR = (M_PI * offset->r) / (180.0);
     362            tmpD = (M_PI * offset->d) / (180.0);
     363        } else if (unit == PS_RADIAN) {
     364            tmpR = offset->r;
     365            tmpD = offset->d;
     366        } else {
     367            psAbort(__func__, "Unknown offset unit: 0x%x\n", unit);
    424368        }
    425     psAbort( __func__, "Unrecognized offset mode\n" );
    426     return ( NULL );
    427 }
     369
     370        tmp = (psSphere *) psAlloc(sizeof(psSphere));
     371        tmp->r = position->r + tmpR;
     372        tmp->r = position->d + tmpD;
     373        tmp->rErr = 0.0;
     374        tmp->dErr = 0.0;
     375
     376        // XXX: wrap tmp->r and tmp->d to the allowed range (-PI to PI)
     377        // and (0 to 2*PI).
     378        return (tmp);
     379    }
     380    psAbort(__func__, "Unrecognized offset mode\n");
     381    return (NULL);
     382}
Note: See TracChangeset for help on using the changeset viewer.