IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 2635


Ignore:
Timestamp:
Dec 6, 2004, 9:49:36 AM (22 years ago)
Author:
gusciora
Message:

Redoing the coordinate transforms.

Location:
trunk/psLib/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

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

    r2602 r2635  
    1010*  @author GLG, MHPCC
    1111*
    12 *  @version $Revision: 1.38 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2004-12-03 19:43:43 $
     12*  @version $Revision: 1.39 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2004-12-06 19:49:36 $
    1414*
    1515*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    1919/******************************************************************************/
    2020#include "psType.h"
    21 #include "psImage.h"
    22 #include "psArray.h"
    23 #include "psList.h"
    2421#include "psCoord.h"
    2522#include "psMemory.h"
    … …  
    3532
    3633#define PS_COT(X) (1.0 / atan(X))
     34#define DEG_TO_RAD(DEGREES) ((DEGREES) * M_PI / 180.0)
     35#define MIN_TO_RAD(MINUTES) ((MINUTES) * M_PI / (180.0 * 60.0))
     36#define SEC_TO_RAD(SECONDS) ((SECONDS) * M_PI / (180.0 * 60.0 * 60.0))
     37#define RAD_TO_DEG(RADIANS) ((RADIANS) * 180 / M_PI)
     38#define RAD_TO_MIN(RADIANS) ((RADIANS) * 180 * 60.0 / M_PI)
     39#define RAD_TO_SEC(RADIANS) ((RADIANS) * 180 * 60.0 * 60.0 / M_PI)
    3740
    3841/******************************************************************************/
    … …  
    162165
    163166/******************************************************************************
    164 This function prototype has been modified since the SDRS.
    165  *****************************************************************************/
    166 psSphereTransform* psSphereTransformAlloc(double alphaP,
    167         double deltaP,
    168         double phiP)
     167alpha is LONGITUDE
     168delta is LATITUDE
     169 
     170    alphaP: Take the target pole in the source system; calculate its LONGITUDE
     171     in the target system.  That longitude is alphaP.
     172    DeltaP: Take the target pole in the source system; calculate its LATITUDE
     173     in the target system.  That longitude is deltaP.
     174    phiP:   This is the LONGITUDE of the ascending node in the target system.
     175 *****************************************************************************/
     176psSphereTransform* psSphereTransformAlloc(psF64 alphaP,
     177        psF64 deltaP,
     178        psF64 phiP)
    169179{
    170180    psSphereTransform* tmp = (psSphereTransform* ) psAlloc(sizeof(psSphereTransform));
    … …  
    179189
    180190/******************************************************************************
    181 This algorithm comes from an email from Gene.  I assume (x,y) corresponds to
    182 (r,d) in the sphere coordinates.
    183  
    184 XXX: In Gene's email, there are different variables with similar names (y
    185 and Y) and in the code, there's cos_y, cos_Y, and cos(y).  Verify that there
    186 are no typo's.
     191XXX: It's not clear if this code correctly implements what the equations in
     192the ADD imply.  Bug 243 will, hopefully, decide that.
     193 
     194XXX: Once this code is verified, delete Gene's code.
    187195 *****************************************************************************/
    188196psSphere* psSphereTransformApply(psSphere* out,
    … …  
    193201    PS_PTR_CHECK_NULL(coord, NULL);
    194202
    195     double sinY = 0.0;
    196     double cosY = 0.0;
    197     double sinX = 0.0;
    198     double cosX = 0.0;
    199     double x = 0.0;
    200     double y = 0.0;
    201     double dx = 0.0;
    202 
    203203    if (out == NULL) {
    204204        out = (psSphere* ) psAlloc(sizeof(psSphere));
    205205    }
    206206
    207     x = coord->r;
    208     y = coord->d;
    209     dx = x - transform->phiP;
    210     sinY = cos(y) * sin(dx) * transform->sinDeltaP +
    211            sin(y) * transform->cosDeltaP;
    212     cosY = sqrt(1.0 - sinY * sinY);
    213     sinX = (cos(y) * sin(dx) * transform->cosDeltaP -
    214             sin(y) * transform->sinDeltaP) / cos(y);
    215     cosX = cos(y) * cos(dx) / cos(y);
    216 
    217     out->r = atan2(sinX, cosX) + transform->alphaP;
    218     out->d = atan2(sinY, cosY);
    219 
    220     return (out);
    221 }
    222 
     207    psF64 alpha = coord->r;
     208    psF64 delta = coord->d;
     209    psF64 alphaMinusAlphaP = alpha - transform->alphaP;
     210
     211    psF64 eq55 = (sin(delta) * transform->cosDeltaP) -
     212                 (cos(delta) * transform->sinDeltaP * sin(alphaMinusAlphaP));
     213
     214    psF64 eq56 = (cos(delta) * transform->cosDeltaP * sin(alphaMinusAlphaP)) +
     215                 (sin(delta) * transform->sinDeltaP);
     216
     217    psF64 eq57 = cos(delta) * cos(alphaMinusAlphaP);
     218
     219    psF64 theta = asin(eq55);
     220    psF64 phi = atan2(eq56, eq57) + transform->alphaP;
     221
     222    out->r = phi;
     223    out->d = theta;
     224
     225    return(out);
     226
     227    /*
     228    This algorithm comes from an email from Gene.  I assume (x,y) corresponds to
     229    (r,d) in the sphere coordinates.
     230     
     231    XXX: In Gene's email, there are different variables with similar names (y
     232    and Y) and in the code, there's cos_y, cos_Y, and cos(y).  Verify that there
     233    are no typo's.
     234     
     235        psF64 sinY = 0.0;
     236        psF64 cosY = 0.0;
     237        psF64 sinX = 0.0;
     238        psF64 cosX = 0.0;
     239        psF64 x = 0.0;
     240        psF64 y = 0.0;
     241        psF64 dx = 0.0;
     242     
     243        x = coord->r;
     244        y = coord->d;
     245        dx = x - transform->phiP;
     246        sinY = cos(y) * sin(dx) * transform->sinDeltaP +
     247               sin(y) * transform->cosDeltaP;
     248        cosY = sqrt(1.0 - sinY * sinY);
     249        sinX = (cos(y) * sin(dx) * transform->cosDeltaP -
     250                sin(y) * transform->sinDeltaP) / cos(y);
     251        cosX = cos(y) * cos(dx) / cos(y);
     252     
     253        out->r = atan2(sinX, cosX) + transform->alphaP;
     254        out->d = atan2(sinY, cosY);
     255    */
     256}
     257
     258// XXX: the ADD changed.  This code no longer conforms to the ADD.  In fact,
     259// there is no algorithm for this in the ADD.
     260// XXX: This is bug 244
    223261psSphereTransform* psSphereTransformICRSToEcliptic(psTime time)
    224262{
    … …  
    226264
    227265    struct tm *tmTime = psTimeToTM(&time);
    228     double year = (double)(1900 + tmTime->tm_year);
    229     double T = year / 100.0;
    230     double phi = -23.452294 + 0.013013 * T + 0.000001639 * T * T - 0.000000503 * T * T * T;
    231     double alphaP = 0.0;
    232     double phiP = 0.0;
     266    psF64 year = (psF64)(1900 + tmTime->tm_year);
     267    psF64 T = year / 100.0;
     268    psF64 phi = DEG_TO_RAD(-23.452294) +
     269                (DEG_TO_RAD(0.013013) * T) +
     270                (DEG_TO_RAD(0.000001639) * T * T) -
     271                (DEG_TO_RAD(0.000000503) * T * T * T);
     272    psF64 alphaP = 0.0;
     273    psF64 phiP = 0.0;
    233274
    234275    psFree(tmTime);
    … …  
    236277}
    237278
     279// XXX: the ADD changed.  This code no longer conforms to the old ADD.
     280// XXX: Waiting for IfA to sepcify if T should include days/months/etc.
    238281psSphereTransform* psSphereTransformEclipticToICRS(psTime time)
    239282{
    … …  
    241284
    242285    struct tm *tmTime = psTimeToTM(&time);
    243     double year = (double)(1900 + tmTime->tm_year);
    244     double T = year / 100.0;
    245     double phi = +23.452294 - 0.013013 * T - 0.000001639 * T * T + 0.000000503 * T * T * T;
    246     double alphaP = 0.0;
    247     double phiP = 0.0;
     286    psF64 year = (psF64)(1900 + tmTime->tm_year);
     287    psF64 T = year / 100.0;
     288    psF64 alphaP = 0.0;
     289    psF64 deltaP = DEG_TO_RAD(23.0) +
     290                   SEC_TO_RAD(27.0) +
     291                   MIN_TO_RAD(8.0) -
     292                   (SEC_TO_RAD(46.845) * T) -
     293                   (SEC_TO_RAD(0.0059) * T * T) +
     294                   (SEC_TO_RAD(0.00181) * T * T * T);
     295    psF64 phiP = 0.0;
    248296
    249297    psFree(tmTime);
    250     return (psSphereTransformAlloc(phi, alphaP, phiP));
    251 }
    252 
    253 // XXX: SHould these be in radians, not degrees?
     298    return (psSphereTransformAlloc(alphaP, deltaP, phiP));
     299}
     300
     301// XXX: Should these be in radians, not degrees?
     302// XXX: the ADD changed.  This code no longer conforms to the ADD.  In fact,
     303// there is no algorithm for this in the ADD.
     304// XXX: This is bug 245
    254305psSphereTransform* psSphereTransformICRSToGalatic(void)
    255306{
    256     return (psSphereTransformAlloc(62.6, 282.25, 33.0));
    257 }
    258 
    259 // XXX: SHould these be in radians, not degrees?
     307    psF64 alphaP = DEG_TO_RAD(-62.6);
     308    psF64 deltaP = DEG_TO_RAD(33.0);
     309    psF64 phiP = DEG_TO_RAD(282.25);
     310
     311    return (psSphereTransformAlloc(alphaP, deltaP, phiP));
     312}
     313
     314// XXX: Should these be in radians, not degrees?
     315// XXX: the ADD changed.  This code no longer conforms to the old ADD.
    260316psSphereTransform* psSphereTransformGalaticToICRS(void)
    261317{
    262     return (psSphereTransformAlloc(-62.6, 33.0, 282.25));
     318    psF64 alphaP = DEG_TO_RAD(282.85948);
     319    psF64 deltaP = DEG_TO_RAD(62.87175);
     320    psF64 phiP = DEG_TO_RAD(32.93192);
     321
     322    return (psSphereTransformAlloc(alphaP, deltaP, phiP));
    263323}
    264324
    … …  
    273333    PS_PTR_CHECK_NULL(projection, NULL);
    274334
    275     float R = 0.0;
    276     float alpha = 0.0;
     335    //    float R = 0.0;
     336    //    float alpha = 0.0;
    277337    psPlane* tmp = (psPlane* ) psAlloc(sizeof(psPlane));
    278338
    279339    if (projection->type == PS_PROJ_TAN) {
    280         R = PS_COT(coord->r) * (180.0 / M_PI);
    281         tmp->x = R * sin(coord->d);
    282         tmp->y = R * cos(coord->d);
     340        tmp->x = cos(coord->r) * sin(coord->d) / sin(coord->r);
     341        tmp->y = -cos(coord->r) * cos(coord->d) / sin(coord->r);
    283342
    284343    } else if (projection->type == PS_PROJ_SIN) {
    285         R = cos(coord->r) * (180.0 / M_PI);
    286         tmp->x = R * sin(coord->d);
    287         tmp->y = R * cos(coord->d);
     344        tmp->x = cos(coord->d) * sin(coord->d);
     345        tmp->y = -cos(coord->d) * cos(coord->d);
    288346
    289347    } else if (projection->type == PS_PROJ_CAR) {
    … …  
    293351    } else if (projection->type == PS_PROJ_MER) {
    294352        tmp->x = coord->d;
    295         tmp->y = log(tan(45.0 + (0.5 * coord->r))) * 180.0 / M_PI;
    296 
     353        tmp->y = log(tan(DEG_TO_RAD(45.0) + (0.5 * coord->r)));
    297354    } else if (projection->type == PS_PROJ_AIT) {
    298         alpha = 1.0 / ((180.0 / M_PI) * sqrt(1.0 + (cos(coord->r) *
    299                                              cos(0.5 * coord->d) * 0.5)));
    300 
    301         tmp->x = 2.0 * alpha * cos(coord->r) * sin(0.5 * coord->d);
    302         tmp->y = alpha * sin(coord->d);
     355        psF64 tmpF64  = sqrt(0.5 * (1.0 + (cos(coord->r) * cos(0.5 * coord->d))));
     356        tmpF64 = 1.0 / tmpF64;
     357        tmp->x = 2.0 * tmpF64 * cos(coord->r) * sin(0.5 * coord->d);
    303358
    304359    } else if (projection->type == PS_PROJ_PAR) {
    … …  
    337392
    338393    if (projection->type == PS_PROJ_TAN) {
     394        tmp->r = atan2(-coord->y, coord->x);
    339395        R = sqrt((coord->x * coord->x) + (coord->y * coord->y));
    340         tmp->d = atan2(-coord->y, coord->x);
    341         tmp->r = atan(180.0 / (R * M_PI));
     396        tmp->d = atan(1.0 / R);
    342397
    343398    } else if (projection->type == PS_PROJ_SIN) {
    … …  
    434489
    435490        if (unit == PS_ARCSEC) {
    436             tmp->r = (tmp->r * 180.0 * 60.0 * 60.0) / M_PI;
    437             tmp->d = (tmp->d * 180.0 * 60.0 * 60.0) / M_PI;
     491            tmp->r = RAD_TO_SEC(tmp->r);
     492            tmp->d = RAD_TO_SEC(tmp->d);
    438493        } else if (unit == PS_ARCMIN) {
    439             tmp->r = (tmp->r * 180.0 * 60.0) / M_PI;
    440             tmp->d = (tmp->d * 180.0 * 60.0) / M_PI;
     494            tmp->r = RAD_TO_MIN(tmp->r);
     495            tmp->d = RAD_TO_MIN(tmp->d);
    441496        } else if (unit == PS_DEGREE) {
    442             tmp->r = (tmp->r * 180.0) / M_PI;
    443             tmp->d = (tmp->d * 180.0) / M_PI;
     497            tmp->r = RAD_TO_DEG(tmp->r);
     498            tmp->d = RAD_TO_DEG(tmp->d);
    444499        } else if (unit == PS_RADIAN) {}
    445500        else {
    … …  
    480535    psSphere* tmp;
    481536    psProjection proj;
    482     double tmpR = 0.0;
    483     double tmpD = 0.0;
     537    psF64 tmpR = 0.0;
     538    psF64 tmpD = 0.0;
    484539
    485540    if (mode == PS_LINEAR) {
    … …  
    495550        tmp = psDeproject(&lin, &proj);
    496551        return (tmp);
    497 
    498552    } else if (mode == PS_SPHERICAL) {
    499553        if (unit == PS_ARCSEC) {
    500             tmpR = (M_PI * offset->r) / (180.0 * 60.0 * 60.0);
    501             tmpD = (M_PI * offset->d) / (180.0 * 60.0 * 60.0);
     554            tmpR = SEC_TO_RAD(offset->r);
     555            tmpD = SEC_TO_RAD(offset->d);
    502556        } else if (unit == PS_ARCMIN) {
    503             tmpR = (M_PI * offset->r) / (180.0 * 60.0);
    504             tmpD = (M_PI * offset->d) / (180.0 * 60.0);
     557            tmpR = MIN_TO_RAD(offset->r);
     558            tmpD = MIN_TO_RAD(offset->d);
    505559        } else if (unit == PS_DEGREE) {
    506             tmpR = (M_PI * offset->r) / (180.0);
    507             tmpD = (M_PI * offset->d) / (180.0);
     560            tmpR = DEG_TO_RAD(offset->r);
     561            tmpD = DEG_TO_RAD(offset->d);
    508562        } else if (unit == PS_RADIAN) {
    509563            tmpR = offset->r;
    … …  
    533587}
    534588
    535 // XXX: What algorithm should be used?
     589
     590/******************************************************************************
     591psSpherePrecess(coords, fromTime, toTime):
     592 
     593XXX: What algorithm should be used?
     594XXX: Verify that this is correct.
     595XXX: Use static memory for tmpST.
     596 *****************************************************************************/
    536597psSphere *psSpherePrecess(psSphere *coords,
    537598                          const psTime *fromTime,
    538599                          const psTime *toTime)
    539600{
    540     /*
    541         psF64 fromMJD = fromTime->sec/86400.0 + fromTime->usec/86400000000.0 + 40587.0;
    542         psF64 toMJD = toTime->sec/86400.0 + toTime->usec/86400000000.0 + 40587.0;
    543         psF64 = (toMJD - fromMJD) / 36525.0;
    544     */
    545     return(NULL);
    546 }
    547 
     601    psF64 fromMJD = fromTime->sec/86400.0 +
     602                    fromTime->usec/86400000000.0 +
     603                    40587.0;
     604    psF64 toMJD = toTime->sec/86400.0 +
     605                  toTime->usec/86400000000.0 +
     606                  40587.0;
     607    psF64 T = (toMJD - fromMJD) / 36525.0;
     608
     609    psF64 alphaP = DEG_TO_RAD(90.0) - ((DEG_TO_RAD(0.6406161) * T) +
     610                                       (DEG_TO_RAD(0.0000839) * T * T) +
     611                                       (DEG_TO_RAD(0.000005) * T * T * T));
     612
     613    psF64 deltaP = (DEG_TO_RAD(0.5567530) * T) -
     614                   (DEG_TO_RAD(0.0001185) * T * T) -
     615                   (DEG_TO_RAD(0.0000116) * T * T * T);
     616
     617    psF64 phiP = DEG_TO_RAD(90.0) - ((DEG_TO_RAD(0.6406161) * T) +
     618                                     (DEG_TO_RAD(0.0003041) * T * T) +
     619                                     (DEG_TO_RAD(0.0000051) * T * T * T));
     620
     621    psSphereTransform *tmpST = psSphereTransformAlloc(alphaP, deltaP, phiP);
     622
     623    psSphere *out = psSphereTransformApply(NULL, tmpST, coords);
     624
     625    psFree(tmpST);
     626
     627    return(out);
     628}
     629

  • trunk/psLib/src/astro/psCoord.h

    r2600 r2635  
    1010*  @author GLG, MHPCC
    1111*
    12 *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2004-12-02 21:12:51 $
     12*  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2004-12-06 19:49:36 $
    1414*
    1515*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    107107typedef struct
    108108{
     109    double alphaP;                    ///< Longitude of the target system pole in the source system
    109110    double cosDeltaP;                 ///< Cosine of target pole latitude in the source system
    110111    double sinDeltaP;                 ///< Sine of target pole latitude in the source system
    111     double alphaP;                    ///< Longitude of the target system pole in the source system
    112112    double phiP;                      ///< Longitude of the ascending node in the target system
    113113}

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

    r2602 r2635  
    1010*  @author GLG, MHPCC
    1111*
    12 *  @version $Revision: 1.38 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2004-12-03 19:43:43 $
     12*  @version $Revision: 1.39 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2004-12-06 19:49:36 $
    1414*
    1515*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    1919/******************************************************************************/
    2020#include "psType.h"
    21 #include "psImage.h"
    22 #include "psArray.h"
    23 #include "psList.h"
    2421#include "psCoord.h"
    2522#include "psMemory.h"
    … …  
    3532
    3633#define PS_COT(X) (1.0 / atan(X))
     34#define DEG_TO_RAD(DEGREES) ((DEGREES) * M_PI / 180.0)
     35#define MIN_TO_RAD(MINUTES) ((MINUTES) * M_PI / (180.0 * 60.0))
     36#define SEC_TO_RAD(SECONDS) ((SECONDS) * M_PI / (180.0 * 60.0 * 60.0))
     37#define RAD_TO_DEG(RADIANS) ((RADIANS) * 180 / M_PI)
     38#define RAD_TO_MIN(RADIANS) ((RADIANS) * 180 * 60.0 / M_PI)
     39#define RAD_TO_SEC(RADIANS) ((RADIANS) * 180 * 60.0 * 60.0 / M_PI)
    3740
    3841/******************************************************************************/
    … …  
    162165
    163166/******************************************************************************
    164 This function prototype has been modified since the SDRS.
    165  *****************************************************************************/
    166 psSphereTransform* psSphereTransformAlloc(double alphaP,
    167         double deltaP,
    168         double phiP)
     167alpha is LONGITUDE
     168delta is LATITUDE
     169 
     170    alphaP: Take the target pole in the source system; calculate its LONGITUDE
     171     in the target system.  That longitude is alphaP.
     172    DeltaP: Take the target pole in the source system; calculate its LATITUDE
     173     in the target system.  That longitude is deltaP.
     174    phiP:   This is the LONGITUDE of the ascending node in the target system.
     175 *****************************************************************************/
     176psSphereTransform* psSphereTransformAlloc(psF64 alphaP,
     177        psF64 deltaP,
     178        psF64 phiP)
    169179{
    170180    psSphereTransform* tmp = (psSphereTransform* ) psAlloc(sizeof(psSphereTransform));
    … …  
    179189
    180190/******************************************************************************
    181 This algorithm comes from an email from Gene.  I assume (x,y) corresponds to
    182 (r,d) in the sphere coordinates.
    183  
    184 XXX: In Gene's email, there are different variables with similar names (y
    185 and Y) and in the code, there's cos_y, cos_Y, and cos(y).  Verify that there
    186 are no typo's.
     191XXX: It's not clear if this code correctly implements what the equations in
     192the ADD imply.  Bug 243 will, hopefully, decide that.
     193 
     194XXX: Once this code is verified, delete Gene's code.
    187195 *****************************************************************************/
    188196psSphere* psSphereTransformApply(psSphere* out,
    … …  
    193201    PS_PTR_CHECK_NULL(coord, NULL);
    194202
    195     double sinY = 0.0;
    196     double cosY = 0.0;
    197     double sinX = 0.0;
    198     double cosX = 0.0;
    199     double x = 0.0;
    200     double y = 0.0;
    201     double dx = 0.0;
    202 
    203203    if (out == NULL) {
    204204        out = (psSphere* ) psAlloc(sizeof(psSphere));
    205205    }
    206206
    207     x = coord->r;
    208     y = coord->d;
    209     dx = x - transform->phiP;
    210     sinY = cos(y) * sin(dx) * transform->sinDeltaP +
    211            sin(y) * transform->cosDeltaP;
    212     cosY = sqrt(1.0 - sinY * sinY);
    213     sinX = (cos(y) * sin(dx) * transform->cosDeltaP -
    214             sin(y) * transform->sinDeltaP) / cos(y);
    215     cosX = cos(y) * cos(dx) / cos(y);
    216 
    217     out->r = atan2(sinX, cosX) + transform->alphaP;
    218     out->d = atan2(sinY, cosY);
    219 
    220     return (out);
    221 }
    222 
     207    psF64 alpha = coord->r;
     208    psF64 delta = coord->d;
     209    psF64 alphaMinusAlphaP = alpha - transform->alphaP;
     210
     211    psF64 eq55 = (sin(delta) * transform->cosDeltaP) -
     212                 (cos(delta) * transform->sinDeltaP * sin(alphaMinusAlphaP));
     213
     214    psF64 eq56 = (cos(delta) * transform->cosDeltaP * sin(alphaMinusAlphaP)) +
     215                 (sin(delta) * transform->sinDeltaP);
     216
     217    psF64 eq57 = cos(delta) * cos(alphaMinusAlphaP);
     218
     219    psF64 theta = asin(eq55);
     220    psF64 phi = atan2(eq56, eq57) + transform->alphaP;
     221
     222    out->r = phi;
     223    out->d = theta;
     224
     225    return(out);
     226
     227    /*
     228    This algorithm comes from an email from Gene.  I assume (x,y) corresponds to
     229    (r,d) in the sphere coordinates.
     230     
     231    XXX: In Gene's email, there are different variables with similar names (y
     232    and Y) and in the code, there's cos_y, cos_Y, and cos(y).  Verify that there
     233    are no typo's.
     234     
     235        psF64 sinY = 0.0;
     236        psF64 cosY = 0.0;
     237        psF64 sinX = 0.0;
     238        psF64 cosX = 0.0;
     239        psF64 x = 0.0;
     240        psF64 y = 0.0;
     241        psF64 dx = 0.0;
     242     
     243        x = coord->r;
     244        y = coord->d;
     245        dx = x - transform->phiP;
     246        sinY = cos(y) * sin(dx) * transform->sinDeltaP +
     247               sin(y) * transform->cosDeltaP;
     248        cosY = sqrt(1.0 - sinY * sinY);
     249        sinX = (cos(y) * sin(dx) * transform->cosDeltaP -
     250                sin(y) * transform->sinDeltaP) / cos(y);
     251        cosX = cos(y) * cos(dx) / cos(y);
     252     
     253        out->r = atan2(sinX, cosX) + transform->alphaP;
     254        out->d = atan2(sinY, cosY);
     255    */
     256}
     257
     258// XXX: the ADD changed.  This code no longer conforms to the ADD.  In fact,
     259// there is no algorithm for this in the ADD.
     260// XXX: This is bug 244
    223261psSphereTransform* psSphereTransformICRSToEcliptic(psTime time)
    224262{
    … …  
    226264
    227265    struct tm *tmTime = psTimeToTM(&time);
    228     double year = (double)(1900 + tmTime->tm_year);
    229     double T = year / 100.0;
    230     double phi = -23.452294 + 0.013013 * T + 0.000001639 * T * T - 0.000000503 * T * T * T;
    231     double alphaP = 0.0;
    232     double phiP = 0.0;
     266    psF64 year = (psF64)(1900 + tmTime->tm_year);
     267    psF64 T = year / 100.0;
     268    psF64 phi = DEG_TO_RAD(-23.452294) +
     269                (DEG_TO_RAD(0.013013) * T) +
     270                (DEG_TO_RAD(0.000001639) * T * T) -
     271                (DEG_TO_RAD(0.000000503) * T * T * T);
     272    psF64 alphaP = 0.0;
     273    psF64 phiP = 0.0;
    233274
    234275    psFree(tmTime);
    … …  
    236277}
    237278
     279// XXX: the ADD changed.  This code no longer conforms to the old ADD.
     280// XXX: Waiting for IfA to sepcify if T should include days/months/etc.
    238281psSphereTransform* psSphereTransformEclipticToICRS(psTime time)
    239282{
    … …  
    241284
    242285    struct tm *tmTime = psTimeToTM(&time);
    243     double year = (double)(1900 + tmTime->tm_year);
    244     double T = year / 100.0;
    245     double phi = +23.452294 - 0.013013 * T - 0.000001639 * T * T + 0.000000503 * T * T * T;
    246     double alphaP = 0.0;
    247     double phiP = 0.0;
     286    psF64 year = (psF64)(1900 + tmTime->tm_year);
     287    psF64 T = year / 100.0;
     288    psF64 alphaP = 0.0;
     289    psF64 deltaP = DEG_TO_RAD(23.0) +
     290                   SEC_TO_RAD(27.0) +
     291                   MIN_TO_RAD(8.0) -
     292                   (SEC_TO_RAD(46.845) * T) -
     293                   (SEC_TO_RAD(0.0059) * T * T) +
     294                   (SEC_TO_RAD(0.00181) * T * T * T);
     295    psF64 phiP = 0.0;
    248296
    249297    psFree(tmTime);
    250     return (psSphereTransformAlloc(phi, alphaP, phiP));
    251 }
    252 
    253 // XXX: SHould these be in radians, not degrees?
     298    return (psSphereTransformAlloc(alphaP, deltaP, phiP));
     299}
     300
     301// XXX: Should these be in radians, not degrees?
     302// XXX: the ADD changed.  This code no longer conforms to the ADD.  In fact,
     303// there is no algorithm for this in the ADD.
     304// XXX: This is bug 245
    254305psSphereTransform* psSphereTransformICRSToGalatic(void)
    255306{
    256     return (psSphereTransformAlloc(62.6, 282.25, 33.0));
    257 }
    258 
    259 // XXX: SHould these be in radians, not degrees?
     307    psF64 alphaP = DEG_TO_RAD(-62.6);
     308    psF64 deltaP = DEG_TO_RAD(33.0);
     309    psF64 phiP = DEG_TO_RAD(282.25);
     310
     311    return (psSphereTransformAlloc(alphaP, deltaP, phiP));
     312}
     313
     314// XXX: Should these be in radians, not degrees?
     315// XXX: the ADD changed.  This code no longer conforms to the old ADD.
    260316psSphereTransform* psSphereTransformGalaticToICRS(void)
    261317{
    262     return (psSphereTransformAlloc(-62.6, 33.0, 282.25));
     318    psF64 alphaP = DEG_TO_RAD(282.85948);
     319    psF64 deltaP = DEG_TO_RAD(62.87175);
     320    psF64 phiP = DEG_TO_RAD(32.93192);
     321
     322    return (psSphereTransformAlloc(alphaP, deltaP, phiP));
    263323}
    264324
    … …  
    273333    PS_PTR_CHECK_NULL(projection, NULL);
    274334
    275     float R = 0.0;
    276     float alpha = 0.0;
     335    //    float R = 0.0;
     336    //    float alpha = 0.0;
    277337    psPlane* tmp = (psPlane* ) psAlloc(sizeof(psPlane));
    278338
    279339    if (projection->type == PS_PROJ_TAN) {
    280         R = PS_COT(coord->r) * (180.0 / M_PI);
    281         tmp->x = R * sin(coord->d);
    282         tmp->y = R * cos(coord->d);
     340        tmp->x = cos(coord->r) * sin(coord->d) / sin(coord->r);
     341        tmp->y = -cos(coord->r) * cos(coord->d) / sin(coord->r);
    283342
    284343    } else if (projection->type == PS_PROJ_SIN) {
    285         R = cos(coord->r) * (180.0 / M_PI);
    286         tmp->x = R * sin(coord->d);
    287         tmp->y = R * cos(coord->d);
     344        tmp->x = cos(coord->d) * sin(coord->d);
     345        tmp->y = -cos(coord->d) * cos(coord->d);
    288346
    289347    } else if (projection->type == PS_PROJ_CAR) {
    … …  
    293351    } else if (projection->type == PS_PROJ_MER) {
    294352        tmp->x = coord->d;
    295         tmp->y = log(tan(45.0 + (0.5 * coord->r))) * 180.0 / M_PI;
    296 
     353        tmp->y = log(tan(DEG_TO_RAD(45.0) + (0.5 * coord->r)));
    297354    } else if (projection->type == PS_PROJ_AIT) {
    298         alpha = 1.0 / ((180.0 / M_PI) * sqrt(1.0 + (cos(coord->r) *
    299                                              cos(0.5 * coord->d) * 0.5)));
    300 
    301         tmp->x = 2.0 * alpha * cos(coord->r) * sin(0.5 * coord->d);
    302         tmp->y = alpha * sin(coord->d);
     355        psF64 tmpF64  = sqrt(0.5 * (1.0 + (cos(coord->r) * cos(0.5 * coord->d))));
     356        tmpF64 = 1.0 / tmpF64;
     357        tmp->x = 2.0 * tmpF64 * cos(coord->r) * sin(0.5 * coord->d);
    303358
    304359    } else if (projection->type == PS_PROJ_PAR) {
    … …  
    337392
    338393    if (projection->type == PS_PROJ_TAN) {
     394        tmp->r = atan2(-coord->y, coord->x);
    339395        R = sqrt((coord->x * coord->x) + (coord->y * coord->y));
    340         tmp->d = atan2(-coord->y, coord->x);
    341         tmp->r = atan(180.0 / (R * M_PI));
     396        tmp->d = atan(1.0 / R);
    342397
    343398    } else if (projection->type == PS_PROJ_SIN) {
    … …  
    434489
    435490        if (unit == PS_ARCSEC) {
    436             tmp->r = (tmp->r * 180.0 * 60.0 * 60.0) / M_PI;
    437             tmp->d = (tmp->d * 180.0 * 60.0 * 60.0) / M_PI;
     491            tmp->r = RAD_TO_SEC(tmp->r);
     492            tmp->d = RAD_TO_SEC(tmp->d);
    438493        } else if (unit == PS_ARCMIN) {
    439             tmp->r = (tmp->r * 180.0 * 60.0) / M_PI;
    440             tmp->d = (tmp->d * 180.0 * 60.0) / M_PI;
     494            tmp->r = RAD_TO_MIN(tmp->r);
     495            tmp->d = RAD_TO_MIN(tmp->d);
    441496        } else if (unit == PS_DEGREE) {
    442             tmp->r = (tmp->r * 180.0) / M_PI;
    443             tmp->d = (tmp->d * 180.0) / M_PI;
     497            tmp->r = RAD_TO_DEG(tmp->r);
     498            tmp->d = RAD_TO_DEG(tmp->d);
    444499        } else if (unit == PS_RADIAN) {}
    445500        else {
    … …  
    480535    psSphere* tmp;
    481536    psProjection proj;
    482     double tmpR = 0.0;
    483     double tmpD = 0.0;
     537    psF64 tmpR = 0.0;
     538    psF64 tmpD = 0.0;
    484539
    485540    if (mode == PS_LINEAR) {
    … …  
    495550        tmp = psDeproject(&lin, &proj);
    496551        return (tmp);
    497 
    498552    } else if (mode == PS_SPHERICAL) {
    499553        if (unit == PS_ARCSEC) {
    500             tmpR = (M_PI * offset->r) / (180.0 * 60.0 * 60.0);
    501             tmpD = (M_PI * offset->d) / (180.0 * 60.0 * 60.0);
     554            tmpR = SEC_TO_RAD(offset->r);
     555            tmpD = SEC_TO_RAD(offset->d);
    502556        } else if (unit == PS_ARCMIN) {
    503             tmpR = (M_PI * offset->r) / (180.0 * 60.0);
    504             tmpD = (M_PI * offset->d) / (180.0 * 60.0);
     557            tmpR = MIN_TO_RAD(offset->r);
     558            tmpD = MIN_TO_RAD(offset->d);
    505559        } else if (unit == PS_DEGREE) {
    506             tmpR = (M_PI * offset->r) / (180.0);
    507             tmpD = (M_PI * offset->d) / (180.0);
     560            tmpR = DEG_TO_RAD(offset->r);
     561            tmpD = DEG_TO_RAD(offset->d);
    508562        } else if (unit == PS_RADIAN) {
    509563            tmpR = offset->r;
    … …  
    533587}
    534588
    535 // XXX: What algorithm should be used?
     589
     590/******************************************************************************
     591psSpherePrecess(coords, fromTime, toTime):
     592 
     593XXX: What algorithm should be used?
     594XXX: Verify that this is correct.
     595XXX: Use static memory for tmpST.
     596 *****************************************************************************/
    536597psSphere *psSpherePrecess(psSphere *coords,
    537598                          const psTime *fromTime,
    538599                          const psTime *toTime)
    539600{
    540     /*
    541         psF64 fromMJD = fromTime->sec/86400.0 + fromTime->usec/86400000000.0 + 40587.0;
    542         psF64 toMJD = toTime->sec/86400.0 + toTime->usec/86400000000.0 + 40587.0;
    543         psF64 = (toMJD - fromMJD) / 36525.0;
    544     */
    545     return(NULL);
    546 }
    547 
     601    psF64 fromMJD = fromTime->sec/86400.0 +
     602                    fromTime->usec/86400000000.0 +
     603                    40587.0;
     604    psF64 toMJD = toTime->sec/86400.0 +
     605                  toTime->usec/86400000000.0 +
     606                  40587.0;
     607    psF64 T = (toMJD - fromMJD) / 36525.0;
     608
     609    psF64 alphaP = DEG_TO_RAD(90.0) - ((DEG_TO_RAD(0.6406161) * T) +
     610                                       (DEG_TO_RAD(0.0000839) * T * T) +
     611                                       (DEG_TO_RAD(0.000005) * T * T * T));
     612
     613    psF64 deltaP = (DEG_TO_RAD(0.5567530) * T) -
     614                   (DEG_TO_RAD(0.0001185) * T * T) -
     615                   (DEG_TO_RAD(0.0000116) * T * T * T);
     616
     617    psF64 phiP = DEG_TO_RAD(90.0) - ((DEG_TO_RAD(0.6406161) * T) +
     618                                     (DEG_TO_RAD(0.0003041) * T * T) +
     619                                     (DEG_TO_RAD(0.0000051) * T * T * T));
     620
     621    psSphereTransform *tmpST = psSphereTransformAlloc(alphaP, deltaP, phiP);
     622
     623    psSphere *out = psSphereTransformApply(NULL, tmpST, coords);
     624
     625    psFree(tmpST);
     626
     627    return(out);
     628}
     629

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

    r2600 r2635  
    1010*  @author GLG, MHPCC
    1111*
    12 *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
    13 *  @date $Date: 2004-12-02 21:12:51 $
     12*  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
     13*  @date $Date: 2004-12-06 19:49:36 $
    1414*
    1515*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    107107typedef struct
    108108{
     109    double alphaP;                    ///< Longitude of the target system pole in the source system
    109110    double cosDeltaP;                 ///< Cosine of target pole latitude in the source system
    110111    double sinDeltaP;                 ///< Sine of target pole latitude in the source system
    111     double alphaP;                    ///< Longitude of the target system pole in the source system
    112112    double phiP;                      ///< Longitude of the ascending node in the target system
    113113}

  • trunk/psLib/src/dataManip/psConstants.h

    r2583 r2635  
    66 *  @author GLG, MHPCC
    77 *
    8  *  @version $Revision: 1.41 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2004-12-01 19:56:06 $
     8 *  @version $Revision: 1.42 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2004-12-06 19:49:36 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2424 *****************************************************************************/
    2525#define PS_ONE 1.0
    26 
    27 
    2826
    2927/*****************************************************************************

  • trunk/psLib/src/math/psConstants.h

    r2583 r2635  
    66 *  @author GLG, MHPCC
    77 *
    8  *  @version $Revision: 1.41 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2004-12-01 19:56:06 $
     8 *  @version $Revision: 1.42 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2004-12-06 19:49:36 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2424 *****************************************************************************/
    2525#define PS_ONE 1.0
    26 
    27 
    2826
    2927/*****************************************************************************
Note: See TracChangeset for help on using the changeset viewer.