IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 2, 2006, 1:19:58 PM (20 years ago)
Author:
drobbin
Message:

Added comments, removed old code & comments, added NULL checks.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r6039 r6309  
    88 *  @author Dave Robbins, MHPCC
    99 *
    10  *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2006-01-18 23:49:06 $
     10 *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2006-02-02 23:19:58 $
    1212 *
    1313 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    4444{
    4545    psSphereRot r,s,t;
    46 
    47     // directly from ADD -- there must be a better way?!
     46    //The following represents: r =a rotation about the z-axis by alphaP, s =a rotation about
     47    // the y-axis by deltaP, and t =a rotation about the z-axis by phiP.
    4848    r.q0=0;
    4949    r.q1=0;
    … …  
    7171bool psMemCheckSphereRot(psPtr ptr)
    7272{
     73    //See if the ptr corresponds to a psSphereRot*
    7374    return ( psMemGetDeallocator(ptr) == (psFreeFunc)sphereRotFree );
    7475}
    7576
    76 
     77//This function is really a second allocate function for psSphereRot's that uses quaternion
     78// component inputs instead of angle inputs as in psSphereRotAlloc.
    7779psSphereRot* psSphereRotQuat(double q0,
    7880                             double q1,
    … …  
    8082                             double q3)
    8183{
     84    //allocate space for a new sphere rotation and set deallocator
    8285    psSphereRot* rot = (psSphereRot*)psAlloc(sizeof(psSphereRot));
    8386    psMemSetDeallocator(rot, (psFreeFunc)sphereRotFree);
    8487
     88    //The magnitude of a rotation quaternion should = 1 so we normalize here in case the
     89    // inputs have not been.
    8590    double len = sqrt(q0*q0 + q1*q1 + q2*q2 + q3*q3);
    8691    rot->q0 = q0 / len;
    … …  
    9297}
    9398
    94 psSphereRot* psSphereRotConjugate(psSphereRot *out, const psSphereRot *in)
    95 {
     99psSphereRot* psSphereRotConjugate(psSphereRot *out,
     100                                  const psSphereRot *in)
     101{
     102    //if input sphere rotation is NULL, return NULL
    96103    if (in == NULL) {
    97104        psError(PS_ERR_BAD_PARAMETER_NULL, true,
    … …  
    99106        return NULL;
    100107    }
     108    //if output sphere rotation is NULL, allocate a new sphere rotation to return
    101109    if (out == NULL) {
    102110        out = (psSphereRot* ) psAlloc(sizeof(psSphereRot));
    103111        psMemSetDeallocator(out, (psFreeFunc)sphereRotFree);
    104112    }
     113    //Since q3 is the magnitude and q0,q1,q2 are direction, the conjugate rotation
     114    // is formed by -q0,-q1,-q2,+q3.  Note that this is the same as q0,q1,q2,-q3.
    105115    out->q0 = -in->q0;
    106116    out->q1 = -in->q1;
    … …  
    111121}
    112122
    113 /******************************************************************************
    114 XXX: We convert Right Ascension angles to the range 0:PI.  Is that acceptable?
    115 XXX: Should we do something for Declination as well?
    116  *****************************************************************************/
    117123psSphere* psSphereRotApply(psSphere* out,
    118124                           const psSphereRot* transform,
    119125                           const psSphere* coord)
    120126{
     127    //Make sure that input coordinates and rotation are not NULL
    121128    PS_ASSERT_PTR_NON_NULL(transform, NULL);
    122129    PS_ASSERT_PTR_NON_NULL(coord, NULL);
    123 
     130    //If output coordinates is NULL allocate a new psSphere
    124131    if (out == NULL) {
    125132        out = psSphereAlloc();
    126133    }
    127     //    double r = -coord->r;
    128     // apply the transform by creating a new psSphereRot from the input coord
     134    //apply the transform by creating a new psSphereRot from the input coord
    129135    // and combining it with the input transform (see ADD)
    130136    double cosD = cos(coord->d);
    … …  
    134140                                 sin(coord->d),
    135141                                 0.0);
    136 
    137     /*    psSphereRot *test = psSphereRotICRSToGalactic();
    138         if (fabs(coord->d-DEG_TO_RAD(27.1283)) < 0.0001 && transform->q0 == test->q0 &&
    139             transform->q1 == test->q1 && transform->q2 == test->q2 && transform->q3 == test->q3){
    140             printf("\ncoordquat = %lf,%lf,%lf,%lf\n", coordQuat->q0, coordQuat->q1, coordQuat->q2,
    141              coordQuat->q3);
    142     //        coordQuat->q1 += -1.0;
    143             test->q3 = - test->q3;
    144             *(psSphereRot**)&transform = test;
    145         }
    146         psFree(test);
    147     */
    148     // Inserted by PAP
     142    //Get the conjugate of the input rotation.  Need to calculate r*p*R to apply rotation
     143    // where R = r-conjugate.
    149144    psSphereRot *conjugate = psSphereRotConjugate(NULL, transform);
    150145    psSphereRot *temp = psSphereRotCombine(NULL, transform, coordQuat);
    151146    psSphereRot *result = psSphereRotCombine(NULL, temp, conjugate);
     147    //From ADD we can find the new sphere coordinates,
     148    // r is calculated as tan^-1 (q1/q0), d is sin^-1 (q2)
    152149    out->r = atan2(result->q1, result->q0);
    153     //    out->r = atan2(result->q1, result->q0);
    154150    out->d = asin(result->q2);
    155151    out->rErr = 0.0;
    156152    out->dErr = 0.0;
    157 
    158     if (out->r < -0.0001) {
     153    //Simply for convention, we make sure all output r-parameters are positive in the range
     154    // of 0 to 2pi
     155    if (out->r < -0.000001) {
    159156        out->r += 2.0 * M_PI;
    160157    }
    … …  
    165162    psFree(coordQuat);
    166163    return out;
    167     // Pau.
    168 
    169     #if 0
    170     //    psSphereRot* coordQuatConjugate = psSphereRotQuat(
    171     //                                       coordQuat->q0, coordQuat->q1, coordQuat->q2, coordQuat->q3);
    172     //    coordQuat = psSphereRotInvert(coordQuat);
    173 
    174     // calculate q=(rp)r'
    175     //    coordQuat = psSphereRotCombine(coordQuat, transform, coordQuat);
    176     //    coordQuat = psSphereRotCombine(coordQuat, coordQuat, coordQuatConjugate);
    177     // N.B., we can recycle coordQuat right away due to the implementation of
    178     // psSphereRotCombine; it puts the input values in a local variable first
    179 
    180     //    out->r = atan2(coordQuat->q1,coordQuat->q0);
    181     //    out->d = asin(coordQuat->q2);
    182 
    183 
    184     //     psSphereRot *inv = psSphereRotInvert(transform);
    185     psSphereRot *inv = (psSphereRot*)psAlloc(sizeof(psSphereRot));
    186     *inv = *transform;
    187     //    psSphereRot *result = psSphereRotCombine(NULL, transform, coordQuat);
    188     //    inv = psSphereRotInvert(inv);
    189     //    psSphereRotCombine(result, result, inv);
    190     //    psSphereRot *result = psSphereRotCombine(NULL, inv, coordQuat);
    191     //    psSphereRotCombine(result, result, transform);
    192     //    psSphereRot *result = psSphereRotCombine(NULL, coordQuat, transform);
    193     //    psSphereRotCombine(result, inv, result);
    194 
    195     psSphereRot *result = psSphereRotCombine(NULL, coordQuat, transform);
    196     inv = psSphereRotInvert(inv);
    197     //    *(psSphereRot**)&transform = psSphereRotInvert( *(psSphereRot**)&transform);
    198     result = psSphereRotCombine(result, inv, result);
    199     //    *(psSphereRot**)&transform = psSphereRotInvert( *(psSphereRot**)&transform);
    200 
    201     //        out->r = atan2(result->q0, result->q1);
    202     out->r = atan2(result->q1, result->q0);
    203     out->d = asin(result->q2);
    204     out->rErr = 0.0;
    205     out->dErr = 0.0;
    206     psFree(inv);
    207     psFree(result);
    208 
    209     psFree(coordQuat);
    210     //    psFree(coordQuatConjugate);
    211 
    212     return out;
    213     #endif
    214164}
    215165
    … …  
    218168                                const psSphereRot* rot2)
    219169{
     170    //Make sure that input rotations are not NULL
    220171    PS_ASSERT_PTR_NON_NULL(rot1, NULL);
    221172    PS_ASSERT_PTR_NON_NULL(rot2, NULL);
    222 
     173    //If output rotation is NULL, allocate a new psSphereRot to return
    223174    if (out == NULL) {
    224175        out = (psSphereRot* ) psAlloc(sizeof(psSphereRot));
    … …  
    235186    double b3 = rot2->q3;
    236187
    237     // following came from ADD
    238     //    out->q0 = b3*a0 + b2*a1 - b1*a2 + b0*a3;
    239     //    out->q1 = b3*a1 - b2*a0 + b1*a3 + b0*a2;
    240     //    out->q2 = b3*a2 + b2*a3 + b1*a0 - b0*a1;
     188    //Combine rot1 & rot2.  Formulas here came from ADD.
    241189    out->q0 = a3*b0 + a0*b3 + a1*b2 - a2*b1;
    242190    out->q1 = a3*b1 - a0*b2 + a1*b3 + a2*b0;
    243191    out->q2 = a3*b2 + a0*b1 - a1*b0 + a2*b3;
    244 
    245192    out->q3 = b3*a3 - b2*a2 - b1*a1 - b0*a0;
    246193
    … …  
    252199                               double phiP)
    253200{
     201    //This function should produce identical results to psSphereRotConjugate in most
     202    // or possibly all situations.  Creates the inverse rotation from the input angles.
    254203    return (psSphereRotAlloc(-phiP, -deltaP, -alphaP));
    255204}
    … …  
    258207{
    259208    psF64 T;
    260 
    261209    // Check for null parameter
    262210    PS_ASSERT_PTR_NON_NULL(time, NULL);
    263 
    264211    // Convert psTime to MJD
    265212    psF64 MJD = psTimeToMJD(time);
    266 
    267213    // Check the specified MJD is greater than 1900
    268214    if ( MJD < MJD_1900 ) {
    … …  
    270216        return NULL;
    271217    }
    272 
    273218    // Calculate number of Julian centuries since 1900
    274219    T = ( MJD - MJD_1900 ) / JULIAN_CENTURY;
    275 
     220    //Formulas for phiP, deltaP, alphaP came from ADD.
    276221    psF64 phiP = - DEG_TO_RAD(270.0);
    277222    psF64 deltaP = - (DEG_TO_RAD(23.0) +
    … …  
    289234{
    290235    psF64 T;
    291 
    292236    // Check for null parameter
    293237    PS_ASSERT_PTR_NON_NULL(time, NULL);
    294 
    295238    // Convert psTime to MJD
    296239    psF64 MJD = psTimeToMJD(time);
    297 
    298240    // Check the specified MJD is greater than 1900
    299241    if ( MJD < MJD_1900 ) {
    … …  
    301243        return NULL;
    302244    }
    303 
    304245    // Calculate number of Julian centuries since 1900
    305246    T = ( MJD - MJD_1900 ) / JULIAN_CENTURY;
    306247
     248    //Formulas for phiP, deltaP, alphaP came from ADD.  Notice that the formulas are the
     249    // same as for EclipticToICRS except that alpha=-phiP, deltaP=-deltaP, phiP=-alphaP to
     250    // produce the inverse rotation as in psSphereRotInverse.
    307251    psF64 alphaP = DEG_TO_RAD(270.0);
    308252    psF64 deltaP = DEG_TO_RAD(23.0) +
    … …  
    317261}
    318262
    319 // XXX: This is bug 245: alphaP swaps with phiP from psSphereTransformGalacticToICRS()
    320263psSphereRot* psSphereRotGalacticToICRS(void)
    321264{
    322     /*    psF64 phiP = - DEG_TO_RAD(180.0-192.85948);
    323         psF64 deltaP = - DEG_TO_RAD(90.0 - 27.12825);
    324         psF64 alphaP = - DEG_TO_RAD(90.0+32.93192);
    325     */
     265    //Formulas for alphaP, deltaP, phiP came from the ADD for ICRSToGalactic.  Notice that
     266    // this is the reason the inverse gets allocated and returned here.
    326267    psF64 alphaP = DEG_TO_RAD(180.0-192.85948);
    327268    psF64 deltaP = DEG_TO_RAD(90.0 - 27.12825);
    … …  
    332273psSphereRot* psSphereRotICRSToGalactic(void)
    333274{
     275    //Formulas for alphaP, deltaP, phiP came from ADD.
    334276    psF64 alphaP = DEG_TO_RAD(180.0-192.85948);
    335277    psF64 deltaP = DEG_TO_RAD(90.0 - 27.12825);
    … …  
    339281}
    340282
    341 /******************************************************************************
    342 The basic idea is to project both positions onto the linear plane, with
    343 position1 at the center, then calculate the linear offset between those
    344 projections.
    345  
    346 XXX: Do I need to check for unacceptable transformation parameters?  Maybe,
    347      if the points are on the North/South Pole, etc?
    348  
    349 XXX: Do I need to somehow scale this projection?
    350  
    351 XXX: Does PS_LINEAR mode make sense?  The result must be returned in psSphere
    352      regardless of the mode.
    353  
    354 XXX: How to compound errors?
    355  *****************************************************************************/
     283//Calculates the difference between coordinates in position1 & 2 and returns this offset.
    356284psSphere* psSphereGetOffset(const psSphere* position1,
    357285                            const psSphere* position2,
    … …  
    359287                            psSphereOffsetUnit unit)
    360288{
     289    //Make sure that input coordinates are not NULL & that mode & unit are valid
    361290    PS_ASSERT_PTR_NON_NULL(position1, NULL);
    362291    PS_ASSERT_PTR_NON_NULL(position2, NULL);
    363 
    364292    // Check positions near 90 degree and issue warnings if necessary
    365293    if (position1->d >= DEG_TO_RAD(90.0)) {
    … …  
    373301        return NULL;
    374302    }
    375 
    376303    // Allocate return structure
    377304    psSphere* tmp = psSphereAlloc();
    … …  
    380307    // onto tangent plane, set point projected into psSphere structure x->r y->d
    381308    if (mode == PS_LINEAR) {
     309        //The basic idea is to project both positions onto the linear plane, with position1
     310        // at the center, then calculate the linear offset between those projections.
    382311        psProjection* proj = psProjectionAlloc(position1->r,
    383312                                               position1->d,
    … …  
    385314                                               1.0,
    386315                                               PS_PROJ_TAN);
    387 
    388316        // Perform projection onto tangent plane
    389317        psPlane* lin = psProject(position2, proj);
    390 
    391318        // Set return values
    392319        tmp->r = lin->x;
    393320        tmp->d = lin->y;
    394 
    395321        // Free data structures allocated
    396322        psFree(proj);
    … …  
    410336        tmp->dErr = 0.0;
    411337
    412         // Convert to desired units
     338        // Convert output to desired units
    413339        if (unit == PS_ARCSEC) {
    414340            tmp->r = RAD_TO_SEC(tmp->r);
    … …  
    428354            return NULL;
    429355        }
     356
    430357        // Invalid mode
    431358    } else {
    432 
    433359        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    434360                PS_ERRORTEXT_psCoord_OFFSET_MODE_UNKNOWN,
    … …  
    442368}
    443369
    444 /******************************************************************************
    445 XXX: Do we need to check for unacceptable transformation parameters?  Maybe,
    446      if the points are on the North/South Pole, etc?
    447  
    448 XXX: Do we need to somehow scale this projection?
    449  
    450 XXX: I copied the algorithm from the ADD exactly.
    451  
    452 XXX: Should we compound errors?
    453  *****************************************************************************/
    454 
     370//Applies the specified offset to the input position coordinates and returns the new position.
    455371psSphere* psSphereSetOffset(const psSphere* position,
    456372                            const psSphere* offset,
    … …  
    458374                            psSphereOffsetUnit unit)
    459375{
     376    //Make sure that input coordinates are not NULL & that mode & unit are valid
    460377    PS_ASSERT_PTR_NON_NULL(position, NULL);
    461378    PS_ASSERT_PTR_NON_NULL(offset, NULL);
    … …  
    469386    // new sphere coordinate
    470387    if (mode == PS_LINEAR) {
    471 
    472388        // Allocate plane coordinate and set coordinate
    473389        psPlane*  lin = psPlaneAlloc();
    474390        lin->x = offset->r;
    475391        lin->y = offset->d;
    476 
    477392        // Allocate and set projection structure
    478393        psProjection* proj = psProjectionAlloc(position->r,
    … …  
    481396                                               1.0,
    482397                                               PS_PROJ_TAN);
    483 
    484398        // Project tangent plane coord to spherical coord
    485399        tmp = psDeproject(lin, proj);
    486 
    487400        // Free data structures used
    488401        psFree(proj);
    … …  
    492405        // to the position and wrap to 0 to 2pi
    493406    } else if (mode == PS_SPHERICAL) {
    494 
    495407        // Convert offset unit to radians
    496408        if (unit == PS_ARCSEC) {
    … …  
    526438        // Invalid mode report error
    527439    } else {
    528 
    529440        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    530441                PS_ERRORTEXT_psCoord_OFFSET_MODE_UNKNOWN,
Note: See TracChangeset for help on using the changeset viewer.