IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 3, 2004, 2:55:17 PM (22 years ago)
Author:
desonia
Message:

doxygen grouping tweaks.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r1325 r1374  
    11/** @file  psCoord.c
    2  *
    3  *  @brief Contains basic coordinate transformation definitions and operations
    4  *
    5  *  This file defines the basic types for astronomical coordinate
    6  *  transformation
    7  *
    8  *  @ingroup AstroImage
    9  *
    10  *  @author George Gusciora, MHPCC
    11  *
    12  *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2004-07-29 02:08:02 $
    14  *
    15  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    16  */
     2*
     3*  @brief Contains basic coordinate transformation definitions and operations
     4*
     5*  This file defines the basic types for astronomical coordinate
     6*  transformation
     7*
     8*  @ingroup CoordinateTransform
     9*
     10*  @author George Gusciora, MHPCC
     11*
     12*  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2004-08-04 00:55:17 $
     14*
     15*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     16*/
    1717
    1818#include "psType.h"
    … …  
    2828
    2929// This is the only function in this file which I understand.
    30 psPlane *psPlaneTransformApply(psPlane *out,
    31                                const psPlaneTransform *transform,
    32                                const psPlane *coords)
    33 {
    34     if (out == NULL) {
    35         out = (psPlane *) psAlloc(sizeof(psPlane));
    36     }
    37     out->x = transform->x->coeff[0][0] +
    38              (transform->x->coeff[1][0] * coords->x) +
    39              (transform->x->coeff[0][1] * coords->y);
    40 
    41     out->y = transform->y->coeff[0][0] +
    42              (transform->y->coeff[1][0] * coords->x) +
    43              (transform->y->coeff[0][1] * coords->y);
    44 
    45     return(out);
     30psPlane *psPlaneTransformApply( psPlane *out,
     31                                const psPlaneTransform *transform,
     32                                const psPlane *coords )
     33{
     34    if ( out == NULL ) {
     35            out = ( psPlane * ) psAlloc( sizeof( psPlane ) );
     36        }
     37    out->x = transform->x->coeff[ 0 ][ 0 ] +
     38             ( transform->x->coeff[ 1 ][ 0 ] * coords->x ) +
     39             ( transform->x->coeff[ 0 ][ 1 ] * coords->y );
     40             
     41    out->y = transform->y->coeff[ 0 ][ 0 ] +
     42             ( transform->y->coeff[ 1 ][ 0 ] * coords->x ) +
     43             ( transform->y->coeff[ 0 ][ 1 ] * coords->y );
     44             
     45    return ( out );
    4646}
    4747
    … …  
    4949// This transformation takes into account parameters beyond an objects
    5050// spatial coordinates: term3 and term4.
    51 psPlane *psPlaneDistortApply(psPlane *out,
    52                              const psPlaneDistort *transform,
    53                              const psPlane *coords,
    54                              float term3,
    55                              float term4)
    56 {
    57     if (out == NULL) {
    58         out = (psPlane *) psAlloc(sizeof(psPlane));
    59     }
    60 
    61     out->x = transform->x->coeff[0][0][0][0] +
    62              (transform->x->coeff[1][0][0][0] * coords->x) +
    63              (transform->x->coeff[0][1][0][0] * coords->y) +
    64              (transform->x->coeff[0][0][1][0] * term3) +
    65              (transform->x->coeff[0][0][0][1] * term4);
    66 
    67     out->y = transform->y->coeff[0][0][0][0] +
    68              (transform->y->coeff[1][0][0][0] * coords->x) +
    69              (transform->y->coeff[0][1][0][0] * coords->y) +
    70              (transform->y->coeff[0][0][1][0] * term3) +
    71              (transform->y->coeff[0][0][0][1] * term4);
    72 
    73     return(out);
     51psPlane *psPlaneDistortApply( psPlane *out,
     52                              const psPlaneDistort *transform,
     53                              const psPlane *coords,
     54                              float term3,
     55                              float term4 )
     56{
     57    if ( out == NULL ) {
     58            out = ( psPlane * ) psAlloc( sizeof( psPlane ) );
     59        }
     60       
     61    out->x = transform->x->coeff[ 0 ][ 0 ][ 0 ][ 0 ] +
     62             ( transform->x->coeff[ 1 ][ 0 ][ 0 ][ 0 ] * coords->x ) +
     63             ( transform->x->coeff[ 0 ][ 1 ][ 0 ][ 0 ] * coords->y ) +
     64             ( transform->x->coeff[ 0 ][ 0 ][ 1 ][ 0 ] * term3 ) +
     65             ( transform->x->coeff[ 0 ][ 0 ][ 0 ][ 1 ] * term4 );
     66             
     67    out->y = transform->y->coeff[ 0 ][ 0 ][ 0 ][ 0 ] +
     68             ( transform->y->coeff[ 1 ][ 0 ][ 0 ][ 0 ] * coords->x ) +
     69             ( transform->y->coeff[ 0 ][ 1 ][ 0 ][ 0 ] * coords->y ) +
     70             ( transform->y->coeff[ 0 ][ 0 ][ 1 ][ 0 ] * term3 ) +
     71             ( transform->y->coeff[ 0 ][ 0 ][ 0 ][ 1 ] * term4 );
     72             
     73    return ( out );
    7474}
    7575
    7676
    7777// This function prototype has been modified since the SDRS.
    78 psSphereTransform *psSphereTransformAlloc(double NPlat,
     78psSphereTransform *psSphereTransformAlloc( double NPlat,
    7979        double Xo,
    80         double xo)
    81 {
    82     psSphereTransform *tmp = (psSphereTransform *) psAlloc(sizeof(psSphereTransform));
    83 
    84     tmp->sinPhi = sin(NPlat);
    85     tmp->cosPhi = cos(NPlat);
     80        double xo )
     81{
     82    psSphereTransform * tmp = ( psSphereTransform * ) psAlloc( sizeof( psSphereTransform ) );
     83   
     84    tmp->sinPhi = sin( NPlat );
     85    tmp->cosPhi = cos( NPlat );
    8686    tmp->Xo = Xo;
    8787    tmp->xo = xo;
    88 
    89     return(tmp);
     88   
     89    return ( tmp );
    9090}
    9191
    9292// I understand this one too.
    93 void p_psSphereTransformFree(psSphereTransform *trans)
    94 {
    95     psFree(trans);
     93void p_psSphereTransformFree( psSphereTransform *trans )
     94{
     95    psFree( trans );
    9696}
    9797
    … …  
    103103// there are no typo's.
    104104
    105 psSphere *psSphereTransformApply(psSphere *out,
    106                                  const psSphereTransform *transform,
    107                                  const psSphere *coord)
     105psSphere *psSphereTransformApply( psSphere *out,
     106                                  const psSphereTransform *transform,
     107                                  const psSphere *coord )
    108108{
    109109    double sinY = 0.0;
    … …  
    114114    double y = 0.0;
    115115    double dx = 0.0;
    116 
    117     if (out == NULL) {
    118         out = (psSphere *) psAlloc(sizeof(psSphere));
    119     }
    120 
     116   
     117    if ( out == NULL ) {
     118            out = ( psSphere * ) psAlloc( sizeof( psSphere ) );
     119        }
     120       
    121121    x = coord->r;
    122122    y = coord->d;
    123     dx    = x - transform->xo;
    124     sinY = cos(y)*sin(dx)*transform->sinPhi + sin(y)*transform->cosPhi;
    125     cosY = sqrt(1.0 - sinY*sinY);
    126     sinX = (cos(y)*sin(dx)*transform->cosPhi - sin(y)*transform->sinPhi) /
    127            cos(y);
    128     cosX = cos(y)*cos(dx) / cos(y);
    129 
    130     out->r = atan2(sinX, cosX) + transform->Xo;
    131     out->d = atan2(sinY, cosY);
    132 
    133     return(out);
    134 }
    135 
    136 psSphereTransform *psSphereTransformICRStoEcliptic(psTime time)
    137 {
    138     struct tm *tmTime = psTimeToTM(time);
    139     double year = (double) (1900 + tmTime->tm_year);
     123    dx = x - transform->xo;
     124    sinY = cos( y ) * sin( dx ) * transform->sinPhi + sin( y ) * transform->cosPhi;
     125    cosY = sqrt( 1.0 - sinY * sinY );
     126    sinX = ( cos( y ) * sin( dx ) * transform->cosPhi - sin( y ) * transform->sinPhi ) /
     127           cos( y );
     128    cosX = cos( y ) * cos( dx ) / cos( y );
     129   
     130    out->r = atan2( sinX, cosX ) + transform->Xo;
     131    out->d = atan2( sinY, cosY );
     132   
     133    return ( out );
     134}
     135
     136psSphereTransform *psSphereTransformICRStoEcliptic( psTime time )
     137{
     138    struct tm * tmTime = psTimeToTM( time );
     139    double year = ( double ) ( 1900 + tmTime->tm_year );
    140140    double T = year / 100.0;
    141     double phi = -23.452294 + 0.013013*T + 0.000001639*T*T - 0.000000503*T*T*T;
     141    double phi = -23.452294 + 0.013013 * T + 0.000001639 * T * T - 0.000000503 * T * T * T;
    142142    double Xo = 0.0;
    143143    double xo = 0.0;
    144 
    145     return(psSphereTransformAlloc(phi, Xo, xo));
    146 }
    147 
    148 psSphereTransform *psSphereTransformEcliptictoICRS(psTime time)
    149 {
    150     struct tm *tmTime = psTimeToTM(time);
    151     double year = (double) (1900 + tmTime->tm_year);
     144   
     145    return ( psSphereTransformAlloc( phi, Xo, xo ) );
     146}
     147
     148psSphereTransform *psSphereTransformEcliptictoICRS( psTime time )
     149{
     150    struct tm * tmTime = psTimeToTM( time );
     151    double year = ( double ) ( 1900 + tmTime->tm_year );
    152152    double T = year / 100.0;
    153     double phi = +23.452294 - 0.013013*T - 0.000001639*T*T + 0.000000503*T*T*T;
     153    double phi = + 23.452294 - 0.013013 * T - 0.000001639 * T * T + 0.000000503 * T * T * T;
    154154    double Xo = 0.0;
    155155    double xo = 0.0;
    156 
    157     return(psSphereTransformAlloc(phi, Xo, xo));
    158 }
    159 
    160 psSphereTransform *psSphereTransformICRStoGalatic(void)
    161 {
    162     return(psSphereTransformAlloc(62.6, 282.25, 33.0));
    163 }
    164 
    165 psSphereTransform *psSphereTransformGalatictoICRS(void)
    166 {
    167     return(psSphereTransformAlloc(-62.6, 33.0, 282.25));
     156   
     157    return ( psSphereTransformAlloc( phi, Xo, xo ) );
     158}
     159
     160psSphereTransform *psSphereTransformICRStoGalatic( void )
     161{
     162    return ( psSphereTransformAlloc( 62.6, 282.25, 33.0 ) );
     163}
     164
     165psSphereTransform *psSphereTransformGalatictoICRS( void )
     166{
     167    return ( psSphereTransformAlloc( -62.6, 33.0, 282.25 ) );
    168168}
    169169
    170170// XXX: Is this the correct way to calculate this?
    171 float cot(float x)
    172 {
    173     return(1.0 / atan(x));
     171float cot( float x )
     172{
     173    return ( 1.0 / atan( x ) );
    174174}
    175175
    176176// This is some kind of arc tan function.
    177 float arg(float x, float y)
    178 {
    179     if (x > 0) {
    180         return(atan(y/x));
    181     } else if ((x==0) && (y==0)) {
    182         return(0.5 * M_PI);
    183     } else if ((x==0) && (y==0)) {
    184         return(-0.5 * M_PI);
    185     } else if ((x==0) && (y==0)) {
    186         return(M_PI + atan(y/x));
    187     } else if ((x==0) && (y==0)) {
    188         return(-M_PI + atan(y/x));
    189     }
    190 
    191     psAbort(__func__, "Unacceptable range for (arg(%f, %f).\n", x, y);
    192     return(0.0);
     177float arg( float x, float y )
     178{
     179    if ( x > 0 ) {
     180            return ( atan( y / x ) );
     181        } else if ( ( x == 0 ) && ( y == 0 ) ) {
     182            return ( 0.5 * M_PI );
     183        } else if ( ( x == 0 ) && ( y == 0 ) ) {
     184            return ( -0.5 * M_PI );
     185        } else if ( ( x == 0 ) && ( y == 0 ) ) {
     186            return ( M_PI + atan( y / x ) );
     187        } else if ( ( x == 0 ) && ( y == 0 ) ) {
     188            return ( -M_PI + atan( y / x ) );
     189        }
     190       
     191    psAbort( __func__, "Unacceptable range for (arg(%f, %f).\n", x, y );
     192    return ( 0.0 );
    193193}
    194194
    195195// XXX: Waiting for the definition of the PS_PROJ_PAR projection.
    196196// XXX: Waiting for the definition of the PS_PROJ_GLS projection.
    197 psPlane *psProject(const psSphere *coord,
    198                    const psProjection *projection)
     197psPlane *psProject( const psSphere *coord,
     198                    const psProjection *projection )
    199199{
    200200    float R = 0.0;
    201201    float alpha = 0.0;
    202     psPlane *tmp= (psPlane *) psAlloc(sizeof(psPlane));
    203 
    204     if (projection->type == PS_PROJ_TAN) {
    205         R = cot(coord->r) * (180.0 / M_PI);
    206         tmp->x = R * sin(coord->d);
    207         tmp->y = R * cos(coord->d);
    208 
    209     } else if (projection->type == PS_PROJ_SIN) {
    210         R = cos(coord->r) * (180.0 / M_PI);
    211         tmp->x = R * sin(coord->d);
    212         tmp->y = R * cos(coord->d);
    213 
    214     } else if (projection->type == PS_PROJ_CAR) {
    215         tmp->x = coord->d;
    216         tmp->y = coord->r;
    217 
    218     } else if (projection->type == PS_PROJ_MER) {
    219         tmp->x = coord->d;
    220         tmp->y = log(tan(45.0 + (0.5 * coord->r))) * 180.0/M_PI;
    221 
    222     } else if (projection->type == PS_PROJ_AIT) {
    223         alpha = 1.0 / ((180.0 / M_PI) *
    224                        sqrt(1.0 + (cos(coord->r) * cos(0.5 * coord->d) * 0.5)));
    225 
    226         tmp->x = 2.0 * alpha * cos(coord->r) * sin(0.5 * coord->d);
    227         tmp->y = alpha * sin(coord->d);
    228 
    229     } else if (projection->type == PS_PROJ_PAR) {
    230         psAbort(__func__, "The projection type PS_PROJ_PAR is undefined.\n");
    231 
    232     } else if (projection->type == PS_PROJ_GLS) {
    233         psAbort(__func__, "The projection type PS_PROJ_GLS is undefined.\n");
    234     }
    235 
    236     return(tmp);
     202    psPlane *tmp = ( psPlane * ) psAlloc( sizeof( psPlane ) );
     203   
     204    if ( projection->type == PS_PROJ_TAN ) {
     205            R = cot( coord->r ) * ( 180.0 / M_PI );
     206            tmp->x = R * sin( coord->d );
     207            tmp->y = R * cos( coord->d );
     208           
     209        } else if ( projection->type == PS_PROJ_SIN ) {
     210            R = cos( coord->r ) * ( 180.0 / M_PI );
     211            tmp->x = R * sin( coord->d );
     212            tmp->y = R * cos( coord->d );
     213           
     214        } else if ( projection->type == PS_PROJ_CAR ) {
     215            tmp->x = coord->d;
     216            tmp->y = coord->r;
     217           
     218        } else if ( projection->type == PS_PROJ_MER ) {
     219            tmp->x = coord->d;
     220            tmp->y = log( tan( 45.0 + ( 0.5 * coord->r ) ) ) * 180.0 / M_PI;
     221           
     222        } else if ( projection->type == PS_PROJ_AIT ) {
     223            alpha = 1.0 / ( ( 180.0 / M_PI ) *
     224                            sqrt( 1.0 + ( cos( coord->r ) * cos( 0.5 * coord->d ) * 0.5 ) ) );
     225                           
     226            tmp->x = 2.0 * alpha * cos( coord->r ) * sin( 0.5 * coord->d );
     227            tmp->y = alpha * sin( coord->d );
     228           
     229        } else if ( projection->type == PS_PROJ_PAR ) {
     230            psAbort( __func__, "The projection type PS_PROJ_PAR is undefined.\n" );
     231           
     232        } else if ( projection->type == PS_PROJ_GLS ) {
     233            psAbort( __func__, "The projection type PS_PROJ_GLS is undefined.\n" );
     234        }
     235       
     236    return ( tmp );
    237237}
    238238
    239239// XXX: Waiting for the definition of the PS_PROJ_PAR projection.
    240240// XXX: Waiting for the definition of the PS_PROJ_GLS projection.
    241 psSphere *psDeproject(const psPlane *coord,
    242                       const psProjection *projection)
     241psSphere *psDeproject( const psPlane *coord,
     242                       const psProjection *projection )
    243243{
    244244    float R = 0.0;
    … …  
    246246    float chu1 = 0.0;
    247247    float chu2 = 0.0;
    248     psSphere *tmp= (psSphere *) psAlloc(sizeof(psSphere));
    249 
    250     if (projection->type == PS_PROJ_TAN) {
    251         R = sqrt((coord->x * coord->x) + (coord->y * coord->y));
    252         tmp->d = arg(-coord->y, coord->x);
    253         tmp->r = atan(180.0 / (R * M_PI));
    254 
    255     } else if (projection->type == PS_PROJ_SIN) {
    256         R = sqrt((coord->x * coord->x) + (coord->y * coord->y));
    257         tmp->d = arg(-coord->y, coord->x);
    258         tmp->r = acos((R * M_PI) / 180.0);
    259 
    260     } else if (projection->type == PS_PROJ_CAR) {
    261         tmp->d = coord->x;
    262         tmp->r = coord->y;
    263 
    264     } else if (projection->type == PS_PROJ_MER) {
    265         tmp->d = coord->x;
    266         tmp->r = (2.0 * atan(exp((coord->y * M_PI/180.0)))) - 180.0;
    267 
    268     } else if (projection->type == PS_PROJ_AIT) {
    269         chu1 = (coord->x * M_PI) / 720.0;
    270         chu1*= chu1;
    271         chu2 = (coord->y * M_PI) / 360.0;
    272         chu2*= chu2;
    273         chu = sqrt(1.0 - chu1 - chu2);
    274         tmp->d = 2.0 * arg((2.0 * chu * chu) - 1.0,
    275                            (coord->x * chu * M_PI)/360.0);
    276         tmp->r = asin((coord->y * chu * M_PI) / 180.0);
    277 
    278     } else if (projection->type == PS_PROJ_PAR) {
    279         psAbort(__func__, "The projection type PS_PROJ_PAR is undefined.\n");
    280 
    281     } else if (projection->type == PS_PROJ_GLS) {
    282         psAbort(__func__, "The projection type PS_PROJ_GLG is undefined.\n");
    283     }
    284 
    285     return(tmp);
     248    psSphere *tmp = ( psSphere * ) psAlloc( sizeof( psSphere ) );
     249   
     250    if ( projection->type == PS_PROJ_TAN ) {
     251            R = sqrt( ( coord->x * coord->x ) + ( coord->y * coord->y ) );
     252            tmp->d = arg( -coord->y, coord->x );
     253            tmp->r = atan( 180.0 / ( R * M_PI ) );
     254           
     255        } else if ( projection->type == PS_PROJ_SIN ) {
     256            R = sqrt( ( coord->x * coord->x ) + ( coord->y * coord->y ) );
     257            tmp->d = arg( -coord->y, coord->x );
     258            tmp->r = acos( ( R * M_PI ) / 180.0 );
     259           
     260        } else if ( projection->type == PS_PROJ_CAR ) {
     261            tmp->d = coord->x;
     262            tmp->r = coord->y;
     263           
     264        } else if ( projection->type == PS_PROJ_MER ) {
     265            tmp->d = coord->x;
     266            tmp->r = ( 2.0 * atan( exp( ( coord->y * M_PI / 180.0 ) ) ) ) - 180.0;
     267           
     268        } else if ( projection->type == PS_PROJ_AIT ) {
     269            chu1 = ( coord->x * M_PI ) / 720.0;
     270            chu1 *= chu1;
     271            chu2 = ( coord->y * M_PI ) / 360.0;
     272            chu2 *= chu2;
     273            chu = sqrt( 1.0 - chu1 - chu2 );
     274            tmp->d = 2.0 * arg( ( 2.0 * chu * chu ) - 1.0,
     275                                ( coord->x * chu * M_PI ) / 360.0 );
     276            tmp->r = asin( ( coord->y * chu * M_PI ) / 180.0 );
     277           
     278        } else if ( projection->type == PS_PROJ_PAR ) {
     279            psAbort( __func__, "The projection type PS_PROJ_PAR is undefined.\n" );
     280           
     281        } else if ( projection->type == PS_PROJ_GLS ) {
     282            psAbort( __func__, "The projection type PS_PROJ_GLG is undefined.\n" );
     283        }
     284       
     285    return ( tmp );
    286286}
    287287
    288288// XXX: Do I need to check for unacceptable transformation parameters?
    289289// Maybe, if the points are on the North/South Pole, etc?
    290 psSphere *psSphereGetOffset(const psSphere *restrict position1,
    291                             const psSphere *restrict position2,
    292                             psSphereOffsetMode mode,
    293                             psSphereOffsetUnit unit)
     290psSphere *psSphereGetOffset( const psSphere *restrict position1,
     291                             const psSphere *restrict position2,
     292                             psSphereOffsetMode mode,
     293                             psSphereOffsetUnit unit )
    294294{
    295295    //    psPlane *lin;
    … …  
    298298    double tmpR = 0.0;
    299299    double tmpD = 0.0;
    300 
    301     if (mode == PS_LINEAR) {
    302         // XXX: I have no idea how to construct this.  Maybe project both
    303         // sperical positions onto the plane, set the origin at one of the
    304         // points on the plane, then deproject?
    305 
    306         // XXX: Do I need to somehow scale this projection?
    307         // project position1?  Will it project to (0.0, 0.0)?
    308         proj.R = position1->r;
    309         proj.D = position1->d;
    310         proj.Xs = 1.0;
    311         proj.Ys = 1.0;
    312         proj.type = PS_PROJ_TAN;
    313 
    314         //        lin = psProject(position2, proj);
    315         //        tmp = psDeproject(lin, proj);
    316 
    317         // XXX: Do we need to convert units in tmp?
    318         return(tmp);
    319     } else if (mode == PS_SPHERICAL) {
    320         tmpR = position2->r - position1->r;
    321         tmpD = position2->d - position1->d;
    322 
    323         if (unit == PS_ARCSEC) {
    324             tmpR = (tmpR * 180.0 * 60.0 * 60.0) / M_PI;
    325             tmpD = (tmpR * 180.0 * 60.0 * 60.0) / M_PI;
    326         } else if (unit == PS_ARCMIN) {
    327             tmpR = (tmpR * 180.0 * 60.0) / M_PI;
    328             tmpD = (tmpR * 180.0 * 60.0) / M_PI;
    329         } else if (unit == PS_DEGREE) {
    330             tmpR = (tmpR * 180.0) / M_PI;
    331             tmpD = (tmpR * 180.0) / M_PI;
    332         } else if (unit == PS_RADIAN) {}
    333         else {
    334             psAbort(__func__, "Unknown offset unit: 0x%x\n", unit);
    335         }
    336 
    337         tmp = (psSphere *) psAlloc(sizeof(psSphere));
    338         tmp->r = tmpR;
    339         tmp->d = tmpD;
    340         tmp->rErr = 0.0;
    341         tmp->dErr = 0.0;
    342         // XXX: Do we need to wrap these to an acceptable range?
    343         return(tmp);
    344     }
    345     psAbort(__func__, "Unrecognized offset mode\n");
    346     return(NULL);
     300   
     301    if ( mode == PS_LINEAR ) {
     302            // XXX: I have no idea how to construct this.  Maybe project both
     303            // sperical positions onto the plane, set the origin at one of the
     304            // points on the plane, then deproject?
     305           
     306            // XXX: Do I need to somehow scale this projection?
     307            // project position1?  Will it project to (0.0, 0.0)?
     308            proj.R = position1->r;
     309            proj.D = position1->d;
     310            proj.Xs = 1.0;
     311            proj.Ys = 1.0;
     312            proj.type = PS_PROJ_TAN;
     313           
     314            //        lin = psProject(position2, proj);
     315            //        tmp = psDeproject(lin, proj);
     316           
     317            // XXX: Do we need to convert units in tmp?
     318            return ( tmp );
     319        } else if ( mode == PS_SPHERICAL ) {
     320            tmpR = position2->r - position1->r;
     321            tmpD = position2->d - position1->d;
     322           
     323            if ( unit == PS_ARCSEC ) {
     324                    tmpR = ( tmpR * 180.0 * 60.0 * 60.0 ) / M_PI;
     325                    tmpD = ( tmpR * 180.0 * 60.0 * 60.0 ) / M_PI;
     326                } else if ( unit == PS_ARCMIN ) {
     327                    tmpR = ( tmpR * 180.0 * 60.0 ) / M_PI;
     328                    tmpD = ( tmpR * 180.0 * 60.0 ) / M_PI;
     329                } else if ( unit == PS_DEGREE ) {
     330                    tmpR = ( tmpR * 180.0 ) / M_PI;
     331                    tmpD = ( tmpR * 180.0 ) / M_PI;
     332                } else if ( unit == PS_RADIAN ) {}
     333            else {
     334                    psAbort( __func__, "Unknown offset unit: 0x%x\n", unit );
     335                }
     336               
     337            tmp = ( psSphere * ) psAlloc( sizeof( psSphere ) );
     338            tmp->r = tmpR;
     339            tmp->d = tmpD;
     340            tmp->rErr = 0.0;
     341            tmp->dErr = 0.0;
     342            // XXX: Do we need to wrap these to an acceptable range?
     343            return ( tmp );
     344        }
     345    psAbort( __func__, "Unrecognized offset mode\n" );
     346    return ( NULL );
    347347}
    348348
    … …  
    351351// Maybe, if the points are on the North/South Pole, etc?
    352352// XXX: I copied the algorithm from the ADD exactly.
    353 psSphere *psSphereSetOffset(const psSphere *restrict position,
    354                             const psSphere *restrict offset,
    355                             psSphereOffsetMode mode,
    356                             psSphereOffsetUnit unit)
     353psSphere *psSphereSetOffset( const psSphere *restrict position,
     354                             const psSphere *restrict offset,
     355                             psSphereOffsetMode mode,
     356                             psSphereOffsetUnit unit )
    357357{
    358358    psPlane lin;
    … …  
    361361    double tmpR = 0.0;
    362362    double tmpD = 0.0;
    363 
    364     if (mode == PS_LINEAR) {
    365         proj.R = position->r;
    366         proj.D = position->d;
    367         proj.Xs = 1.0;
    368         proj.Ys = 1.0;
    369         proj.type = PS_PROJ_TAN;
    370 
    371         lin.x = offset->r;
    372         lin.y = offset->d;
    373 
    374         tmp = psDeproject(&lin, &proj);
    375         return(tmp);
    376 
    377     } else if (mode == PS_SPHERICAL) {
    378         if (unit == PS_ARCSEC) {
    379             tmpR = (M_PI * offset->r) / (180.0 * 60.0 * 60.0);
    380             tmpD = (M_PI * offset->d) / (180.0 * 60.0 * 60.0);
    381         } else if (unit == PS_ARCMIN) {
    382             tmpR = (M_PI * offset->r) / (180.0 * 60.0);
    383             tmpD = (M_PI * offset->d) / (180.0 * 60.0);
    384         } else if (unit == PS_DEGREE) {
    385             tmpR = (M_PI * offset->r) / (180.0);
    386             tmpD = (M_PI * offset->d) / (180.0);
    387         } else if (unit == PS_RADIAN) {
    388             tmpR = offset->r;
    389             tmpD = offset->d;
    390         } else {
    391             psAbort(__func__, "Unknown offset unit: 0x%x\n", unit);
    392         }
    393 
    394         tmp = (psSphere *) psAlloc(sizeof(psSphere));
    395         tmp->r = position->r + tmpR;
    396         tmp->r = position->d + tmpD;
    397         tmp->rErr = 0.0;
    398         tmp->dErr = 0.0;
    399 
    400         // XXX: wrap tmp->r and tmp->d to the allowed range (-PI to PI)
    401         // and (0 to 2*PI).
    402         return(tmp);
    403     }
    404     psAbort(__func__, "Unrecognized offset mode\n");
    405     return(NULL);
    406 }
     363   
     364    if ( mode == PS_LINEAR ) {
     365            proj.R = position->r;
     366            proj.D = position->d;
     367            proj.Xs = 1.0;
     368            proj.Ys = 1.0;
     369            proj.type = PS_PROJ_TAN;
     370           
     371            lin.x = offset->r;
     372            lin.y = offset->d;
     373           
     374            tmp = psDeproject( &lin, &proj );
     375            return ( tmp );
     376           
     377        } else if ( mode == PS_SPHERICAL ) {
     378            if ( unit == PS_ARCSEC ) {
     379                    tmpR = ( M_PI * offset->r ) / ( 180.0 * 60.0 * 60.0 );
     380                    tmpD = ( M_PI * offset->d ) / ( 180.0 * 60.0 * 60.0 );
     381                } else if ( unit == PS_ARCMIN ) {
     382                    tmpR = ( M_PI * offset->r ) / ( 180.0 * 60.0 );
     383                    tmpD = ( M_PI * offset->d ) / ( 180.0 * 60.0 );
     384                } else if ( unit == PS_DEGREE ) {
     385                    tmpR = ( M_PI * offset->r ) / ( 180.0 );
     386                    tmpD = ( M_PI * offset->d ) / ( 180.0 );
     387                } else if ( unit == PS_RADIAN ) {
     388                    tmpR = offset->r;
     389                    tmpD = offset->d;
     390                } else {
     391                    psAbort( __func__, "Unknown offset unit: 0x%x\n", unit );
     392                }
     393               
     394            tmp = ( psSphere * ) psAlloc( sizeof( psSphere ) );
     395            tmp->r = position->r + tmpR;
     396            tmp->r = position->d + tmpD;
     397            tmp->rErr = 0.0;
     398            tmp->dErr = 0.0;
     399           
     400            // XXX: wrap tmp->r and tmp->d to the allowed range (-PI to PI)
     401            // and (0 to 2*PI).
     402            return ( tmp );
     403        }
     404    psAbort( __func__, "Unrecognized offset mode\n" );
     405    return ( NULL );
     406}
Note: See TracChangeset for help on using the changeset viewer.