IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5483


Ignore:
Timestamp:
Nov 7, 2005, 10:52:43 AM (21 years ago)
Author:
drobbin
Message:

update doxygen, implemented sphereRot functions in EOC.

Location:
trunk/psLib
Files:
6 edited

Legend:

Unmodified
Added
Removed

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

    r5466 r5483  
    99*  @author Robert Daniel DeSonia, MHPCC
    1010*
    11 *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
    12 *  @date $Date: 2005-11-03 04:31:22 $
     11*  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
     12*  @date $Date: 2005-11-07 20:52:43 $
    1313*
    1414*  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    5252static bool eocInitialized = false;
    5353
     54static psSphereRot *rotMatrix_To_Quat(double A[3][3]);
     55
     56static void earthPoleFree(psEarthPole *pole)
     57{
     58    // There are non dynamic allocated items
     59}
     60
     61psEarthPole *psEarthPoleAlloc(void)
     62{
     63    psEarthPole* pole = psAlloc(sizeof(psEarthPole));
     64    psMemSetDeallocator(pole, (psFreeFunc) earthPoleFree);
     65    return pole;
     66}
     67
    5468static bool eocInit()
    5569{
    … …  
    139153    psFree(xTable);
    140154    psFree(yTable);
    141     psFree(sTable);
     155    psFree(sTable);    // There are non dynamic allocated items
     156
    142157
    143158    xTable = NULL;
    … …  
    179194        apparent->dErr = 0.0;
    180195    }
    181     psSphere *rp = psSphereAlloc();
    182     psSphere *r_p = psSphereAlloc();
     196    //    psSphere *rp = psSphereAlloc();
     197    psCube *rp = psCubeAlloc();
     198    //    psSphere *r_p = psSphereAlloc();
     199    psCube *r_p = psCubeAlloc();
    183200    double mu = 0.0;
    184201    double mu_p = 0.0;
    … …  
    188205    //    mu = apparent->r * direction->r + apparent->d * direction->d;
    189206    psCube* directionVector = psSphereToCube(direction);
     207    if (directionVector == NULL)
     208        printf("directionVector is null\n");
    190209    psCube* actualVector = psSphereToCube(actual);
     210    if (directionVector == NULL)
     211        printf("actualVector is null\n");
    191212    mu = acos(directionVector->x*actualVector->x +
    192213              directionVector->y*actualVector->y +
    … …  
    194215
    195216    //rp = apparent - mu * direction;
    196     rp->r = actual->r - mu * direction->r;
    197     rp->d = actual->d - mu * direction->d;
     217    //    rp->r = actual->r - mu * direction->r;
     218    //    rp->d = actual->d - mu * direction->d;
     219    rp->x = actualVector->x - mu * directionVector->x;
     220    rp->y = actualVector->y - mu * directionVector->y;
     221    rp->z = actualVector->z - mu * directionVector->z;
    198222
    199223    mu_p = mu + speed * ((mu * mu - 1.0) / (1.0 - speed * mu));
    200224
    201225    //Not sure if this is right.  ADD gets a scalar from division of rp (a vector)
    202     psCube* rpVector = psSphereToCube(rp);
    203     a = sqrt( (1.0 - mu_p * mu_p) /
    204               (sqrt(rpVector->x*rpVector->x
    205                     + rpVector->y*rpVector->y
    206                     + rpVector->z*rpVector->z)) );
    207 
     226    //Paul claims rp is the modulus or square of the vector components
     227    /*    psCube* rpVector = psSphereToCube(rp);
     228        a = sqrt( (1.0 - mu_p * mu_p) / ( rpVector->x*rpVector->x +
     229                rpVector->y*rpVector->y + rpVector->z*rpVector->z)  );
     230    */
     231    a = (1.0 - mu_p * mu_p) / (rp->x * rp->x + rp->y * rp->y + rp->z * rp->z);
     232    if (a < 0.0) {
     233        printf("a is negative\n");
     234        psFree(rp);
     235        psFree(r_p);
     236        psFree(directionVector);
     237        psFree(actualVector);
     238        psFree(apparent)
     239        return NULL;
     240    } else
     241        a = sqrt(a);
    208242    //r_p = mu_p * direction + a * rp;
    209     r_p->r = mu_p * direction->r + a * rp->r;
    210     r_p->d = mu_p * direction->d + a * rp->d;
    211 
    212     *apparent = *r_p;
     243    //    r_p->r = mu_p * direction->r + a * rp->r;
     244    //    r_p->d = mu_p * direction->d + a * rp->d;
     245    r_p->x = mu_p * directionVector->x + a * rp->x;
     246    r_p->y = mu_p * directionVector->y + a * rp->y;
     247    r_p->z = mu_p * directionVector->z + a * rp->z;
     248
     249    psSphere *r_pSphere = psCubeToSphere(r_p);
     250    if (r_pSphere == NULL)
     251        printf("r_pSphere is null\n");
     252
     253    *apparent = *r_pSphere;
    213254    /*
    214255        psSphereRot *rot = NULL;
    … …  
    223264    psFree(rp);
    224265    psFree(r_p);
     266    psFree(r_pSphere)
    225267    psFree(directionVector);
    226268    psFree(actualVector);
    227     psFree(rpVector);
     269    //    psFree(rpVector);
    228270    return apparent;
    229271}
    … …  
    280322    deflection = SEC_TO_RAD(deflection);
    281323    theta = atan(r0/PS_AU) * tan(deflection);
    282     phi = sqrt( deflection*deflection - theta*theta );
     324    //    phi = sqrt( deflection*deflection - theta*theta );
     325    phi = deflection * cos(asin(theta/deflection));
    283326    apparent->r = theta;
    284327    apparent->d = phi;
    … …  
    447490}
    448491
     492static psSphereRot *rotMatrix_To_Quat(double A[3][3])
     493{
     494    int i;
     495    psSphereRot *new = (psSphereRot*)psAlloc(sizeof(psSphereRot));
     496    //    psMemSetDeallocator(new, (psFreeFunc)sphereRotFree);
     497    //Convert rotation matrix to quaternions
     498    double diag_sum[3];
     499    int maxi;
     500    double recip;
     501    diag_sum[0] = 1.0 + A[0][0] - A[1][1] - A[2][2];
     502    diag_sum[1] = 1.0 - A[0][0] + A[1][1] - A[2][2];
     503    diag_sum[2] = 1.0 - A[0][0] - A[1][1] + A[2][2];
     504    diag_sum[3] = 1.0 + A[0][0] + A[1][1] + A[2][2];
     505
     506    maxi = 0;
     507    for (i = 1; i < 4; ++i) {
     508        if (diag_sum[i] > diag_sum[maxi]) {
     509            maxi = i;
     510        }
     511    }
     512
     513    double p = 0.5 * sqrt(diag_sum[maxi]);
     514    recip = 1.0 / (4.0 * p);
     515
     516    if (maxi == 0) {
     517        new->q0 = p;
     518        new->q1 = recip * (A[0][1] + A[1][0]);
     519        new->q2 = recip * (A[2][0] + A[0][2]);
     520        new->q3 = recip * (A[1][2] - A[2][1]);
     521    } else if (maxi == 1) {
     522        new->q0 = recip * (A[0][1] + A[1][0]);
     523        new->q1 = p;
     524        new->q2 = recip * (A[1][2] + A[2][1]);
     525        new->q3 = recip * (A[2][0] - A[0][2]);
     526    } else if (maxi == 2) {
     527        new->q0 = recip * (A[2][0] + A[0][2]);
     528        new->q1 = recip * (A[1][2] + A[2][1]);
     529        new->q2 = p;
     530        new->q3 = recip * (A[0][1] - A[1][0]);
     531    } else if (maxi == 3) {
     532        new->q0 = recip * (A[1][2] - A[2][1]);
     533        new->q1 = recip * (A[2][0] - A[0][2]);
     534        new->q2 = recip * (A[0][1] - A[1][0]);
     535        new->q3 = p;
     536    }
     537    return new;
     538}
    449539
    450540psSphereRot* psSphereRot_CEOtoGCRS(const psEarthPole *pole)
    451541{
    452     return NULL;
     542    PS_ASSERT_PTR_NON_NULL(pole,NULL);
     543    double A[3][3];
     544    psSphereRot *out = NULL;
     545
     546    //Setup the rotation matrix and scalar value, a, as outlined by the ADD
     547    //XXX: Used formula for rotation matrix D from mathworld for z-axis rotation
     548    double a =  1.0 / (1.0 + sqrt(1.0 - (pole->x*pole->x + pole->y*pole->y) ) );
     549    A[0][0] = (1.0 - a*pole->x*pole->x)*cos(pole->s) - a*pole->x*pole->y*sin(pole->s);
     550    A[0][1] = -a*pole->x*pole->y*cos(pole->s) + (1.0 - a*pole->y*pole->y)*sin(pole->s);
     551    A[0][2] = pole->x*cos(pole->s) + pole->y*sin(pole->s);
     552    A[1][0] = -(1.0 - a*pole->x*pole->x)*sin(pole->s) - a*pole->x*pole->y*cos(pole->s);
     553    A[1][1] = a*pole->x*pole->y*sin(pole->s) + (1.0 - a*pole->y*pole->y)*cos(pole->s);
     554    A[1][2] = -pole->x*sin(pole->s) + pole->y*cos(pole->s);
     555    A[2][0] = -pole->x;
     556    A[2][1] = -pole->y;
     557    A[2][2] = 1.0 - a*(pole->x*pole->x + pole->y*pole->y);
     558
     559    out = rotMatrix_To_Quat(A);
     560
     561    return out;
    453562}
    454563
    … …  
    495604psSphereRot* psSphereRot_ITRStoTEO(const psEarthPole* motion)
    496605{
    497     return NULL;
    498 }
     606    PS_ASSERT_PTR_NON_NULL(motion,NULL);
     607    double A[3][3];
     608    psSphereRot *out = NULL;
     609
     610    //Setup Rotation Matrix for transformation (x,y,z rotation)
     611    //XXX: May need to be (z,y,x as in Mathworld?)
     612    A[0][0] = cos(motion->x)*cos(-motion->s);
     613    A[0][1] = sin(motion->y)*sin(motion->x)*cos(-motion->s) + cos(motion->y)*sin(-motion->s);
     614    A[0][2] = -cos(motion->y)*sin(motion->x)*cos(-motion->s) + sin(motion->y)*sin(-motion->s);
     615    A[1][0] = -cos(motion->x)*sin(-motion->s);
     616    A[1][1] = -sin(motion->y)*sin(motion->x)*sin(-motion->s) + cos(motion->y)*cos(-motion->s);
     617    A[1][2] = cos(motion->y)*sin(motion->x)*sin(-motion->s) + sin(motion->y)*cos(-motion->s);
     618    A[2][0] = sin(motion->x);
     619    A[2][1] = -sin(motion->y)*cos(motion->x);
     620    A[2][2] = cos(motion->y)*cos(motion->x);
     621
     622    //Convert rotation matrix to quaternions
     623    double diag_sum[3];
     624    int maxi;
     625    double recip;
     626    diag_sum[0] = 1.0 + A[0][0] - A[1][1] - A[2][2];
     627    diag_sum[1] = 1.0 - A[0][0] + A[1][1] - A[2][2];
     628    diag_sum[2] = 1.0 - A[0][0] - A[1][1] + A[2][2];
     629    diag_sum[3] = 1.0 + A[0][0] + A[1][1] + A[2][2];
     630
     631    maxi = 0;
     632    for (int i = 1; i < 4; ++i) {
     633        if (diag_sum[i] > diag_sum[maxi]) {
     634            maxi = i;
     635        }
     636    }
     637
     638    double p = 0.5 * sqrt(diag_sum[maxi]);
     639    recip = 1.0 / (4.0 * p);
     640
     641    if (maxi == 0) {
     642        out->q0 = p;
     643        out->q1 = recip * (A[0][1] + A[1][0]);
     644        out->q2 = recip * (A[2][0] + A[0][2]);
     645        out->q3 = recip * (A[1][2] - A[2][1]);
     646    } else if (maxi == 1) {
     647        out->q0 = recip * (A[0][1] + A[1][0]);
     648        out->q1 = p;
     649        out->q2 = recip * (A[1][2] + A[2][1]);
     650        out->q3 = recip * (A[2][0] - A[0][2]);
     651    } else if (maxi == 2) {
     652        out->q0 = recip * (A[2][0] + A[0][2]);
     653        out->q1 = recip * (A[1][2] + A[2][1]);
     654        out->q2 = p;
     655        out->q3 = recip * (A[0][1] - A[1][0]);
     656    } else if (maxi == 3) {
     657        out->q0 = recip * (A[1][2] - A[2][1]);
     658        out->q1 = recip * (A[2][0] - A[0][2]);
     659        out->q2 = recip * (A[0][1] - A[1][0]);
     660        out->q3 = p;
     661    }
     662
     663    return out;
     664}

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

    r5466 r5483  
    99*  @author Robert Daniel DeSonia, MHPCC
    1010*
    11 *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
    12 *  @date $Date: 2005-11-03 04:31:22 $
     11*  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
     12*  @date $Date: 2005-11-07 20:52:43 $
    1313*
    1414*  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    4444}
    4545psPrecessMethod;
     46
     47/** Allocates a new psEarthPole structure.  */
     48psEarthPole *psEarthPoleAlloc(void);
    4649
    4750/** Calculates the actual position of a star, given its apparent position and the

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

    r5444 r5483  
    77 *  @author Robert DeSonia, MHPCC
    88 *
    9  *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2005-10-25 00:38:00 $
     9 *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2005-11-07 20:52:43 $
    1111 *
    1212 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    186186 *  coordinate systems.
    187187 *
    188  *  @return psSphereTransform*     transform for ICRS->Ecliptic coordinate systems
     188 *  @return psSphereRot*     transform for ICRS->Ecliptic coordinate systems
    189189 */
    190190psSphereRot* psSphereRotICRSToEcliptic(
    … …  
    195195 *  coordinate systems.
    196196 *
    197  *  @return psSphereTransform*     transform for Ecliptic->ICRS coordinate systems
     197 *  @return psSphereRot*     transform for Ecliptic->ICRS coordinate systems
    198198 */
    199199psSphereRot* psSphereRotEclipticToICRS(
    … …  
    204204 *  coordinate systems.
    205205 *
     206 *  @return psSphereRot*        new sphere rotation for ICRS to Galactic transformations.
    206207 */
    207208psSphereRot* psSphereRotICRSToGalactic(void);
    … …  
    210211 *  coordinate systems.
    211212 *
     213 *  @return psSphereRot*        new sphere rotation for Galactic to ICRS transformations.
    212214 */
    213215psSphereRot* psSphereRotGalacticToICRS(void);

  • trunk/psLib/test/astro/tst_psEarthOrientation.c

    r5466 r5483  
    55*  @author d-Rob, MHPCC
    66*
    7 *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
    8 *  @date $Date: 2005-11-03 04:31:22 $
     7*  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
     8*  @date $Date: 2005-11-07 20:52:43 $
    99*
    1010*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    6363    psSphere *empty = NULL;
    6464
    65     actual->r = 0.2;
    66     actual->d = 0.2;
    67     direction->r = 0.2035;
    68     direction->d = 0.2035;
     65    //    actual->r = 0.2;
     66    //    actual->d = 0.2;
     67    actual->r = DEG_TO_RAD(45.0);
     68    actual->d = DEG_TO_RAD(30.0);
     69    //    direction->r = 0.2035;
     70    //    direction->d = 0.2035;
     71    direction->r = DEG_TO_RAD(48.0);
     72    direction->d = DEG_TO_RAD(20.7072);
    6973
    7074    empty = psAberration(empty, apparent, direction, 0.1);
    … …  
    8286
    8387    apparent = psAberration(apparent, actual, direction, 0.4);
    84     printf("\napparent = r,d  = %.8g, %.8g\n", apparent->r, apparent->d);
     88    if (apparent == NULL)
     89        printf("\nApparent is NULL\n");
     90    else
     91        printf("\napparent = r,d  = %.8g, %.8g\n", apparent->r, apparent->d);
    8592
    8693    psFree(apparent);

  • trunk/psLib/test/astro/verified/tst_psEarthOrientation.stderr

    r5455 r5483  
    2424
    2525---> TESTPOINT PASSED (psEarthOrientation{psGravityDeflect()} | tst_psEarthOrientation.c)
    26 
    27 /***************************** TESTPOINT ******************************************\
    28 *             TestFile: tst_psEarthOrientation.c                                   *
    29 *            TestPoint: psEarthOrientation{psEOCParallax()}                        *
    30 *             TestType: Positive                                                   *
    31 \**********************************************************************************/
    32 
    33 
    34 ---> TESTPOINT PASSED (psEarthOrientation{psEOCParallax()} | tst_psEarthOrientation.c)
    3526
    3627/***************************** TESTPOINT ******************************************\
    … …  
    6354/***************************** TESTPOINT ******************************************\
    6455*             TestFile: tst_psEarthOrientation.c                                   *
    65 *            TestPoint: psEarthOrientation{psSphereRotTransforms()}                *
     56*            TestPoint: psEarthOrientation{psSphereRot_TEOtoCEO()}                 *
    6657*             TestType: Positive                                                   *
    6758\**********************************************************************************/
    6859
    6960
    70 ---> TESTPOINT PASSED (psEarthOrientation{psSphereRotTransforms()} | tst_psEarthOrientation.c)
     61---> TESTPOINT PASSED (psEarthOrientation{psSphereRot_TEOtoCEO()} | tst_psEarthOrientation.c)
    7162
     63/***************************** TESTPOINT ******************************************\
     64*             TestFile: tst_psEarthOrientation.c                                   *
     65*            TestPoint: psEarthOrientation{psSphereRot_CEOtoGCRS()}                *
     66*             TestType: Positive                                                   *
     67\**********************************************************************************/
     68
     69
     70---> TESTPOINT PASSED (psEarthOrientation{psSphereRot_CEOtoGCRS()} | tst_psEarthOrientation.c)
     71
     72/***************************** TESTPOINT ******************************************\
     73*             TestFile: tst_psEarthOrientation.c                                   *
     74*            TestPoint: psEarthOrientation{psSphereRRot_ITRStoTEO()}               *
     75*             TestType: Positive                                                   *
     76\**********************************************************************************/
     77
     78
     79---> TESTPOINT PASSED (psEarthOrientation{psSphereRRot_ITRStoTEO()} | tst_psEarthOrientation.c)
     80

  • trunk/psLib/test/astro/verified/tst_psEarthOrientation.stdout

    r5455 r5483  
    11
    2 apparent = r,d  = 0.10218436, 0.10218436
     2apparent = r,d  = 1.2163241, 1.0535409
    33
    44Actual r,d = 0.2035,0.2035    Apparent r,d = 0.2035000000002, 0.2035000000391
     5
     6The Value of T is = 1128530000.931
     7
     8The Value of T is = 1126078455.931
     9Output sphere rotation = 0.000000,0.000000,0.770554,0.637375
Note: See TracChangeset for help on using the changeset viewer.