IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5493


Ignore:
Timestamp:
Nov 9, 2005, 2:13:51 PM (21 years ago)
Author:
drobbin
Message:

Implemented functions PolarTideCorr, SphereRot CEOtoGCRS & ITRStoTEO (& tests) - still need verification!

Location:
trunk/psLib
Files:
5 edited

Legend:

Unmodified
Added
Removed

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

    r5483 r5493  
    22*
    33*  @brief Function implementations for earth orientation calculations
    4 *  transformation
    54*
    65*  @ingroup EarthOrientation
    … …  
    98*  @author Robert Daniel DeSonia, MHPCC
    109*
    11 *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
    12 *  @date $Date: 2005-11-07 20:52:43 $
     10*  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
     11*  @date $Date: 2005-11-10 00:13:50 $
    1312*
    1413*  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    6362    psEarthPole* pole = psAlloc(sizeof(psEarthPole));
    6463    psMemSetDeallocator(pole, (psFreeFunc) earthPoleFree);
     64    pole->x = 0.0;
     65    pole->y = 0.0;
     66    pole->s = 0.0;
    6567    return pole;
    6668}
    … …  
    210212    if (directionVector == NULL)
    211213        printf("actualVector is null\n");
    212     mu = acos(directionVector->x*actualVector->x +
    213               directionVector->y*actualVector->y +
    214               directionVector->z*actualVector->z);
     214    //    mu = acos(directionVector->x*actualVector->x +
     215    mu = (directionVector->x*actualVector->x +
     216          directionVector->y*actualVector->y +
     217          directionVector->z*actualVector->z);
    215218
    216219    //rp = apparent - mu * direction;
    … …  
    288291    // use dot product to calculate the angle of separation
    289292    // N.B., assuming the psSphereToCube function returns a unit vector.
    290     double theta = acos(sunVector->x*actualVector->x +
    291                         sunVector->y*actualVector->y +
    292                         sunVector->z*actualVector->z);
     293    //    double theta = acos(sunVector->x*actualVector->x +
     294    double theta = (sunVector->x*actualVector->x +
     295                    sunVector->y*actualVector->y +
     296                    sunVector->z*actualVector->z);
    293297
    294298    double r0 = PS_AU * tan(theta);
    … …  
    322326    deflection = SEC_TO_RAD(deflection);
    323327    theta = atan(r0/PS_AU) * tan(deflection);
    324     //    phi = sqrt( deflection*deflection - theta*theta );
    325     phi = deflection * cos(asin(theta/deflection));
     328    phi = sqrt( deflection*deflection - theta*theta );
     329    //    phi = deflection * cos(asin(theta/deflection));
    326330    apparent->r = theta;
    327331    apparent->d = phi;
    … …  
    582586}
    583587
    584 
    585588psEarthPole* psEOC_GetPolarMotion(const psTime *time,
    586589                                  psTimeBulletin bulletin)
    587590{
     591
    588592    return NULL;
    589593}
    590594
     595static double DMOD(double x, double y)
     596{
     597    double value = x - y * trunc(x/y);
     598    return value;
     599}
    591600
    592601psEarthPole* psEOC_PolarTideCorr(const psTime *time)
    593602{
     603    // Check for null parameter
     604    PS_ASSERT_PTR_NON_NULL(time, NULL);
     605    psEarthPole *out = psEarthPoleAlloc();
     606
     607    // Convert psTime to MJD
     608    double MJD = psTimeToMJD(time);
     609
     610    // Calculate number of Julian centuries since 2000
     611    //XXX: NOT SURE IF THIS IS CORRECT FOR THIS SITUATION
     612    double RJD = ( MJD - MJD_2000 ) / JULIAN_CENTURY;
     613
     614    //Formula comes from fortran reference
     615    //DMOD in fortran ref. = double remainder -> x - y * trunc(x/y)
     616    double T, L, LPRIME, CAPF, CAPD, OMEGA, THETA, CORX, CORY, CORZ;
     617    double ARG1, ARG2, ARG3, ARG4, ARG5, ARG6, ARG7, ARG8;
     618    double T2, T3, T4;
     619    T = (RJD - 51544.5) / 36525.0;
     620    T2 = T*T;
     621    T3 = T*T*T;
     622    T4 = T*T*T*T;
     623    L = -0.0002447 * T4 + 0.051635 * T3 + 31.8792 * T2 + 1717915923.2178 * T + 485868.249036;
     624    L = DMOD(L, 1296000.0);
     625    LPRIME = -0.00001149 * T4 - 0.000136 * T3 - 0.5532 * T2 + 129596581.0481 * T + 1287104.79305;
     626    LPRIME = DMOD(LPRIME, 1296000.0);
     627    CAPF = 0.00000417 * T4 - 0.001037 * T3 - 12.7512 * T2 + 1739527262.8478 * T + 335779.526232;
     628    CAPF = DMOD(CAPF, 1296000.0);
     629    CAPD = -0.00003169 * T4 + 0.006593 * T3 - 6.3706 * T2 + 1602961601.209 * T + 1072260.70369;
     630    CAPD = DMOD(CAPD, 1296000.0);
     631    OMEGA = -0.00005939 * T4 + 0.007702 * T3 + 7.4722 * T2 - 6962890.2665 * T + 450160.398036;
     632    OMEGA = DMOD(OMEGA, 1296000.0);
     633    THETA = (67310.54841 + (876600.0 * 3600.0 + 8640184.812866) * T + 0.093104 * T2 -
     634             6.2e-6 * T3) * 15.0 + 648000.0;
     635    ARG7 = DMOD((-L - 2.0 * CAPF - 2.0 * OMEGA + THETA) * M_PI / 648000.0, 2.0 * M_PI)
     636           - M_PI / 2.0;
     637    ARG1 = DMOD((-2.0 * CAPF - 2.0 * OMEGA + THETA) * M_PI / 648000.0, 2.0 * M_PI) - M_PI / 2.0;
     638    ARG2 = DMOD((-2.0 * CAPF + 2.0 * CAPD - 2.0 * OMEGA + THETA) * M_PI / 648000.0, 2.0 * M_PI)
     639           - M_PI / 2.0;
     640    ARG3 = DMOD(THETA * M_PI / 648000.0, 2.0 * M_PI) - M_PI / 2.0;
     641    ARG4 = DMOD((-L - 2.0 * CAPF - 2.0 * OMEGA + 2.0 * THETA) * M_PI / 648000.0, 2.0 * M_PI);
     642    ARG5 = DMOD((-2.0 * CAPF - 2.0 * OMEGA + 2.0 * THETA) * M_PI / 648000.0, 2.0 * M_PI);
     643    ARG6 = DMOD((-2.0 * CAPF + 2.0 * CAPD - 2.0 * OMEGA + 2.0 * THETA) * M_PI / 648000.0,
     644                2.0 * M_PI);
     645    ARG8 = DMOD((2.0 * THETA) * M_PI / 648000.0, 2.0 * M_PI);
     646    CORX = -0.026 * sin(ARG7) + 0.006 * cos(ARG7)
     647           -0.133 * sin(ARG1) + 0.049 * cos(ARG1)
     648           -0.050 * sin(ARG2) + 0.025 * cos(ARG2)
     649           -0.152 * sin(ARG3) + 0.078 * cos(ARG3)
     650           -0.057 * sin(ARG4) - 0.013 * cos(ARG4)
     651           -0.330 * sin(ARG5) - 0.028 * cos(ARG5)
     652           -0.145 * sin(ARG6) + 0.064 * cos(ARG6)
     653           -0.036 * sin(ARG8) + 0.017 * cos(ARG8);
     654    CORY = -0.006 * sin(ARG7) - 0.026 * cos(ARG7)
     655           -0.049 * sin(ARG1) - 0.133 * cos(ARG1)
     656           -0.025 * sin(ARG2) - 0.050 * cos(ARG2)
     657           -0.078 * sin(ARG3) - 0.152 * cos(ARG3)
     658           +0.011 * sin(ARG4) + 0.033 * cos(ARG4)
     659           +0.037 * sin(ARG5) + 0.196 * cos(ARG5)
     660           +0.059 * sin(ARG6) + 0.087 * cos(ARG6)
     661           +0.018 * sin(ARG8) + 0.022 * cos(ARG8);
     662    CORZ =  0.0245 * sin(ARG7) + 0.0503 * cos(ARG7)
     663            +0.1210 * sin(ARG1) + 0.1605 * cos(ARG1)
     664            +0.0286 * sin(ARG2) + 0.0516 * cos(ARG2)
     665            +0.0864 * sin(ARG3) + 0.1771 * cos(ARG3)
     666            -0.0380 * sin(ARG4) - 0.0154 * cos(ARG4)
     667            -0.1617 * sin(ARG5) - 0.0720 * cos(ARG5)
     668            -0.0759 * sin(ARG6) - 0.0004 * cos(ARG6)
     669            -0.0196 * sin(ARG8) - 0.0038 * cos(ARG8);
     670    CORX = CORX * 1.0e-3;
     671    CORY = CORY * 1.0e-3;
     672    CORZ = CORZ * 0.1e-3;
     673
     674    out->x = CORX;
     675    out->y = CORY;
     676    out->s = CORZ;
     677
     678    return out;
     679}
     680
     681psEarthPole* psEOC_NutationCorr(psTime *time)
     682{
    594683    return NULL;
    595684}
    596 
    597 
    598 psEarthPole* psEOC_NutationCorr(psTime *time)
    599 {
    600     return NULL;
    601 }
    602 
    603685
    604686psSphereRot* psSphereRot_ITRStoTEO(const psEarthPole* motion)
    … …  
    621703
    622704    //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;
     705    out = rotMatrix_To_Quat(A);
     706    /*    double diag_sum[3];
     707        int maxi;
     708        double recip;
     709        diag_sum[0] = 1.0 + A[0][0] - A[1][1] - A[2][2];
     710        diag_sum[1] = 1.0 - A[0][0] + A[1][1] - A[2][2];
     711        diag_sum[2] = 1.0 - A[0][0] - A[1][1] + A[2][2];
     712        diag_sum[3] = 1.0 + A[0][0] + A[1][1] + A[2][2];
     713     
     714        maxi = 0;
     715        for (int i = 1; i < 4; ++i) {
     716            if (diag_sum[i] > diag_sum[maxi]) {
     717                maxi = i;
     718            }
    635719        }
    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 
     720     
     721        double p = 0.5 * sqrt(diag_sum[maxi]);
     722        recip = 1.0 / (4.0 * p);
     723     
     724        if (maxi == 0) {
     725            out->q0 = p;
     726            out->q1 = recip * (A[0][1] + A[1][0]);
     727            out->q2 = recip * (A[2][0] + A[0][2]);
     728            out->q3 = recip * (A[1][2] - A[2][1]);
     729        } else if (maxi == 1) {
     730            out->q0 = recip * (A[0][1] + A[1][0]);
     731            out->q1 = p;
     732            out->q2 = recip * (A[1][2] + A[2][1]);
     733            out->q3 = recip * (A[2][0] - A[0][2]);
     734        } else if (maxi == 2) {
     735            out->q0 = recip * (A[2][0] + A[0][2]);
     736            out->q1 = recip * (A[1][2] + A[2][1]);
     737            out->q2 = p;
     738            out->q3 = recip * (A[0][1] - A[1][0]);
     739        } else if (maxi == 3) {
     740            out->q0 = recip * (A[1][2] - A[2][1]);
     741            out->q1 = recip * (A[2][0] - A[0][2]);
     742            out->q2 = recip * (A[0][1] - A[1][0]);
     743            out->q3 = p;
     744        }
     745    */
    663746    return out;
    664747}

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

    r5444 r5493  
    88 *  @author Dave Robbins, MHPCC
    99 *
    10  *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2005-10-25 00:38:00 $
     10 *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2005-11-10 00:13:50 $
    1212 *
    1313 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii

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

    r5483 r5493  
    55*  @author d-Rob, MHPCC
    66*
    7 *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
    8 *  @date $Date: 2005-11-07 20:52:43 $
     7*  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
     8*  @date $Date: 2005-11-10 00:13:51 $
    99*
    1010*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    6969    //    direction->r = 0.2035;
    7070    //    direction->d = 0.2035;
    71     direction->r = DEG_TO_RAD(48.0);
     71    //    direction->r = DEG_TO_RAD(48.0);
     72    direction->r = DEG_TO_RAD(-157.0);
    7273    direction->d = DEG_TO_RAD(20.7072);
    7374
    … …  
    108109    sun->r = 0.2;
    109110    sun->d = 0.2;
    110     actual->r = 0.2035;
    111     actual->d = 0.2035;
     111    actual->r = 0.61001;
     112    actual->d = 0.01999;
    112113
    113114    empty = psGravityDeflection(apparent, empty, sun);
    … …  
    126127    apparent = psGravityDeflection(apparent, actual, sun);
    127128    psSphere *result = psSphereSetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN);
    128     printf("\nActual r,d = %.13g,%.13g    Apparent r,d = %.13g, %.13g \n",
     129    printf("\nActual r,d = %.13g,%.13g    Apparent r,d = %.16g, %.16g \n",
    129130           actual->r, actual->d, result->r, result->d);
    130131    psFree(result);
    … …  
    144145psS32 testEOCPolar(void)
    145146{
    146 
     147    //    psTime *in = psTimeGetNow(PS_TIME_TAI);
     148    psTime *in = psTimeAlloc(PS_TIME_TAI);
     149    in->sec = 1131579114;
     150    in->nsec = 498489000;
     151    in->leapsecond = false;
     152    psTime *empty = NULL;
     153    psEarthPole *eop = NULL;
     154
     155    eop = psEOC_PolarTideCorr(empty);
     156    if (eop != NULL) {
     157        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     158                "psEOC_PolarTideCorr failed to return NULL for NULL input time.\n");
     159        return 1;
     160    }
     161
     162    eop = psEOC_PolarTideCorr(in);
     163    if (eop == NULL) {
     164        psError(PS_ERR_BAD_PARAMETER_NULL, false,
     165                "psEOC_PolarTideCorr returned NULL for valid input time.\n");
     166        return 2;
     167    }
     168
     169    printf(" -- Eop = x=%.8g, y=%.8g, s=%.8g\n", eop->x, eop->y, eop->s);
     170    psFree(in);
     171    psFree(eop);
    147172    return 0;
    148173}
    … …  
    170195psS32 testSphereRot_CEOtoGCRS(void)
    171196{
     197    psEarthPole *in = psEarthPoleAlloc();
     198    psEarthPole *empty = NULL;
     199    psSphereRot *rot = NULL;
     200
     201    in->x = M_PI / 4.0;
     202    in->y = M_PI / 6.0;
     203    in->s = M_PI / 8.0;
     204
     205    rot = psSphereRot_CEOtoGCRS(empty);
     206    if (rot != NULL) {
     207        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     208                "psSphereRot_CEOtoGCRS failed to return NULL for NULL input psEarthPole.\n");
     209        return 1;
     210    }
     211
     212    rot = psSphereRot_CEOtoGCRS(in);
     213    if (rot == NULL) {
     214        psError(PS_ERR_BAD_PARAMETER_NULL, false,
     215                "psSphereRot_CEOtoGCRS return NULL for valid psEarthPole input.\n");
     216        return 2;
     217    }
     218
     219    psSphere *test = psSphereAlloc();
     220    test->r = 0.0;
     221    test->d = 0.0;
     222    test = psSphereRotApply(test, rot, test);
     223    printf("\n  Sphere -test- has values r,d = %.8g, %.8g \n", test->r, test->d);
     224
     225    psFree(test);
     226    psFree(rot);
     227    psFree(in);
    172228    return 0;
    173229}
    … …  
    175231psS32 testSphereRot_ITRStoTEO(void)
    176232{
    177     return 0;
    178 }
    179 
     233    psEarthPole *in = psEarthPoleAlloc();
     234    psEarthPole *empty = NULL;
     235    psSphereRot *rot = NULL;
     236
     237    in->x = M_PI / 4.0;
     238    in->y = M_PI / 6.0;
     239    in->s = M_PI / 8.0;
     240
     241    rot = psSphereRot_ITRStoTEO(empty);
     242    if (rot != NULL) {
     243        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     244                "psSphereRot_ITRStoTEO failed to return NULL for NULL input psEarthPole.\n");
     245        return 1;
     246    }
     247
     248    rot = psSphereRot_ITRStoTEO(in);
     249    if (rot == NULL) {
     250        psError(PS_ERR_BAD_PARAMETER_NULL, false,
     251                "psSphereRot_ITRStoTEO return NULL for valid psEarthPole input.\n");
     252        return 2;
     253    }
     254
     255    psSphere *test = psSphereAlloc();
     256    test->r = 0.0;
     257    test->d = 0.0;
     258    test = psSphereRotApply(test, rot, test);
     259    printf("\n  Sphere -test- has values r,d = %.8g, %.8g \n", test->r, test->d);
     260
     261    psFree(test);
     262    psFree(rot);
     263    psFree(in);
     264    return 0;
     265}
     266

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

    r5483 r5493  
    4040\**********************************************************************************/
    4141
     42<DATE><TIME>|<HOST>|E|psEOC_PolarTideCorr (FILE:LINENO)
     43    Unallowable operation: time is NULL.
    4244
    4345---> TESTPOINT PASSED (psEarthOrientation{psEOCPolar()} | tst_psEarthOrientation.c)
    … …  
    6769\**********************************************************************************/
    6870
     71<DATE><TIME>|<HOST>|E|psSphereRot_CEOtoGCRS (FILE:LINENO)
     72    Unallowable operation: pole is NULL.
    6973
    7074---> TESTPOINT PASSED (psEarthOrientation{psSphereRot_CEOtoGCRS()} | tst_psEarthOrientation.c)
    … …  
    7680\**********************************************************************************/
    7781
     82<DATE><TIME>|<HOST>|E|psSphereRot_ITRStoTEO (FILE:LINENO)
     83    Unallowable operation: motion is NULL.
    7884
    7985---> TESTPOINT PASSED (psEarthOrientation{psSphereRRot_ITRStoTEO()} | tst_psEarthOrientation.c)

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

    r5483 r5493  
    11
    2 apparent = r,d  = 1.2163241, 1.0535409
     2apparent = r,d  = 1.0900423, 0.75502353
    33
    4 Actual r,d = 0.2035,0.2035    Apparent r,d = 0.2035000000002, 0.2035000000391
     4Actual r,d = 0.61001,0.01999    Apparent r,d = 0.6100100000001365, 0.01999000000006505
     5 -- Eop = x=-0.00045805944, y=0.00013368871, s=-2.986086e-05
    56
    67The Value of T is = 1128530000.931
    … …  
    89The Value of T is = 1126078455.931
    910Output sphere rotation = 0.000000,0.000000,0.770554,0.637375
     11
     12  Sphere -test- has values r,d = 0.048163027, 1.1836384
     13
     14  Sphere -test- has values r,d = 6.2758795, -0.85886273
Note: See TracChangeset for help on using the changeset viewer.