IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5642


Ignore:
Timestamp:
Nov 30, 2005, 1:50:40 PM (21 years ago)
Author:
drobbin
Message:

Revised/Updated tests to use provided test values. NONE of the values match. fabulous.

Location:
trunk/psLib
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r5626 r5642  
    88 *  @author Robert Daniel DeSonia, MHPCC
    99 *
    10  *  @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2005-11-30 02:17:17 $
     10 *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2005-11-30 23:50:40 $
    1212 *
    1313 *  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    727727    }
    728728    double T = (double)(in->sec) + (double)(in->nsec / 1e9);
    729     printf("\nThe Value of T is = %.13g\n", T);
     729    printf("\nThe Value of T is = %.19g\n", T);
    730730    T += -2451545.0;
    731     printf("\nThe Value of T is = %.13g\n", T);
    732731    double theta = 2.0 * M_PI * (0.7790572732640 + 1.00273781191135448 * T);
     732    printf("\nThe Value of theta is = %.13g\n", theta);
     733    //    psSphereRot *out = psSphereRotInvert(theta, 0.0, 0.0);
    733734    psSphereRot *out = psSphereRotInvert(theta, 0.0, 0.0);
    734735
    … …  
    745746    psSphere *in = NULL;
    746747    double s;
     748    psTime *time2 = psTimeAlloc(time->type);
     749    time2->sec = time->sec;
     750    time2->nsec = time->nsec;
     751    time2->leapsecond = time->leapsecond;
     752    //XXX:  Time must be converted to TAI before use?
     753    if (time->type == PS_TIME_UT1) {
     754        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "psTime cannot be of type UT1 here.\n");
     755        return NULL;
     756    } else if (time->type != PS_TIME_TAI) {
     757        time2 = psTimeConvert(time2, PS_TIME_TAI);
     758    }
    747759
    748760    //XXX: This may be the wrong idea... ADD says to use 3rd-order polys to interpolate
    749761    //polar motion coordinates.
    750     in = p_psTimeGetPoleCoords(time);
     762    in = p_psTimeGetPoleCoords(time2);
    751763    if (in == NULL) {
    752764        psError(PS_ERR_BAD_PARAMETER_NULL, true,
    … …  
    757769    out->y = in->d;
    758770
    759     s = psTimeGetUT1Delta(time, bulletin);
     771    s = psTimeGetUT1Delta(time2, bulletin);
    760772    out->s = s;
    761773
    762774    //XXX: Apply polar tide correction here???
    763     psEarthPole *correction = psEOC_PolarTideCorr(time);
     775    psEarthPole *correction = psEOC_PolarTideCorr(time2);
    764776    out->x += correction->x;
    765777    out->y += correction->y;
    766778    out->s += correction->s;
    767779
    768 
     780    psFree(time2);
    769781    psFree(in);
    770782    psFree(correction);

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

    r5626 r5642  
    55*  @author d-Rob, MHPCC
    66*
    7 *  @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
    8 *  @date $Date: 2005-11-30 02:17:17 $
     7*  @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
     8*  @date $Date: 2005-11-30 23:50:40 $
    99*
    1010*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    1313#include "psTest.h"
    1414#include "pslib_strict.h"
     15
    1516
    1617static psS32 testAberration(void);
    … …  
    5455}
    5556
     57#define timesec 1049160600
     58
    5659psS32 testAberration(void)
    5760{
     61
    5862    psSphere *apparent = NULL;
    5963    psSphere *actual = psSphereAlloc();
    60     psSphere *direction = psSphereAlloc();
     64    //    psSphere *direction = psSphereAlloc();
    6165    psSphere *empty = NULL;
    6266
    6367    //    actual->r = 0.2;
    6468    //    actual->d = 0.2;
    65     actual->r = DEG_TO_RAD(45.0);
    66     actual->d = DEG_TO_RAD(30.0);
     69    //    actual->r = DEG_TO_RAD(45.0);
     70    //    actual->d = DEG_TO_RAD(30.0);
     71    actual->r = DEG_TO_RAD(122.9153182445501);
     72    actual->d = DEG_TO_RAD(48.562968978679194);
    6773    //    direction->r = 0.2035;
    6874    //    direction->d = 0.2035;
    6975    //    direction->r = DEG_TO_RAD(48.0);
    70     direction->r = DEG_TO_RAD(-157.0);
    71     direction->d = DEG_TO_RAD(20.7072);
    72 
    73     empty = psAberration(empty, apparent, direction, 0.1);
     76    //    direction->r = DEG_TO_RAD(-157.0);
     77    //    direction->d = DEG_TO_RAD(20.7072);
     78    psCube *cubeDir = psCubeAlloc();
     79    cubeDir->x = 5148.713262821658;
     80    cubeDir->y = -26945.04752348012;
     81    cubeDir->z = -11682.787302030947;
     82    cubeDir->x += -357.6031690489248;
     83    cubeDir->y += 248.46429758174693;
     84    cubeDir->z += 0.09694774143797581;
     85    psSphere *direction = psCubeToSphere(cubeDir);
     86    double speed = sqrt(cubeDir->x*cubeDir->x + cubeDir->y*cubeDir->y + cubeDir->z*cubeDir->z);
     87    // Speed of light in vacuum (src:NIST)
     88    double c;
     89    c = 299792458.0; /* m/s */
     90    speed = speed / c;
     91
     92    empty = psAberration(empty, apparent, direction, speed);
    7493    if (empty != NULL) {
    7594        psError(PS_ERR_BAD_PARAMETER_NULL, false,
    … …  
    7796        return 1;
    7897    }
    79     empty = psAberration(empty, actual, apparent, 0.1);
     98    empty = psAberration(empty, actual, apparent, speed);
    8099    if (empty != NULL) {
    81100        psError(PS_ERR_BAD_PARAMETER_NULL, false,
    … …  
    84103    }
    85104
    86     apparent = psAberration(apparent, actual, direction, 0.4);
     105    apparent = psAberration(apparent, actual, direction, speed);
    87106    if (apparent == NULL)
    88107        printf("\nApparent is NULL\n");
    89     else
     108    else {
    90109        printf("\napparent = r,d  = %.8g, %.8g\n", apparent->r, apparent->d);
    91 
     110        psFree(cubeDir);
     111        cubeDir = psSphereToCube(apparent);
     112        printf("  ><><Apparent = x,y,z = %.13g, %.13g, %.13g \n", cubeDir->x, cubeDir->y, cubeDir->z);
     113    }
     114    psFree(cubeDir);
    92115    psFree(apparent);
    93116    psFree(actual);
    … …  
    159182{
    160183    psTime *empty = NULL;
    161     psTime *time = psTimeAlloc(PS_TIME_TAI);
    162     time->sec = 1131579114;
    163     time->nsec = 498489000;
     184    psTime *time = psTimeAlloc(PS_TIME_UTC);
     185    time->sec = timesec;
     186    time->nsec = 0;
    164187    time->leapsecond = false;
    165188
    … …  
    189212        return 3;
    190213    } else {
    191         printf("Precession Model output = x,y,s = %.8g, %.8g, %.8g\n\n",
     214        printf("Precession Model output = x,y,s = %.8g, %.8g, %.8g\n",
    192215               pmodel->x, pmodel->y, pmodel->s);
     216        double x, y, s;
     217        x = 2.857175590089105e-4;
     218        y = 2.3968739377734732e-5;
     219        s = -1.3970066457904322e-8;
     220        printf("Expected output = x,y,s = %.13g, %.13g, %.13g\n\n", x, y, s);
    193221    }
    194222    psFree(pmodel);
    195223
    196 
     224    //XXX:  MAY Be an Error here... Time for Precession Correction has to be TAI or seg faults
     225    psTime *time2 = psTimeAlloc(PS_TIME_UTC);
     226    *time2 = *time;
     227    time2 = psTimeConvert(time2, PS_TIME_TAI);
    197228    //Tests for Precession Correction function//
    198229    //Return NULL for NULL time input
    … …  
    206237
    207238    //Return NULL for Invalid IERS table
    208     pcorr = psEOC_PrecessionCorr(time, 3);
     239    pcorr = psEOC_PrecessionCorr(time2, 3);
    209240    if (pcorr != NULL) {
    210241        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    … …  
    215246
    216247    //Check values from IERS table A
    217     pcorr = psEOC_PrecessionCorr(time, PS_IERS_A);
     248    pcorr = psEOC_PrecessionCorr(time2, PS_IERS_A);
    218249    if ( pcorr == NULL ) {
    219250        psError(PS_ERR_BAD_PARAMETER_NULL, false,
    … …  
    227258
    228259    //Check values from IERS table B
    229     pcorr = psEOC_PrecessionCorr(time, PS_IERS_B);
     260    pcorr = psEOC_PrecessionCorr(time2, PS_IERS_B);
    230261    if ( pcorr == NULL ) {
    231262        psError(PS_ERR_BAD_PARAMETER_NULL, false,
    … …  
    235266        printf("Precession Correction output (IERSB) = x,y,s = %.8g, %.8g, %.8g\n",
    236267               pcorr->x, pcorr->y, pcorr->s);
     268        double xx, yy, ss;
     269        xx = 0.06295703125;
     270        yy = -0.0287618408203125;
     271        ss = 0.0;
     272        printf("Expected output = x,y,s = %.13g, %.13g, %.13g\n\n", xx, yy, ss);
    237273    }
    238274    psFree(pcorr);
    239275
    240276    psFree(time);
     277    psFree(time2);
    241278    if (!p_psEOCFinalize() ) {
    242279        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "EOC failed finalization!\n");
    … …  
    250287{
    251288    //    psTime *in = psTimeGetNow(PS_TIME_TAI);
    252     psTime *in = psTimeAlloc(PS_TIME_TAI);
    253     in->sec = 1131579114;
    254     in->nsec = 498489000;
     289    //    psTime *in = psTimeAlloc(PS_TIME_TAI);
     290    //    in->sec = 1131579114;
     291    //    in->nsec = 498489000;
     292    psTime *in = psTimeAlloc(PS_TIME_UTC);
     293    in->sec = timesec;
     294    in->nsec = 0;
    255295    in->leapsecond = false;
    256296    psTime *empty = NULL;
    … …  
    284324        return 4;
    285325    }
    286 
    287     psSphere *interm = p_psTimeGetPoleCoords(in);
     326    psTime *in2 = psTimeAlloc(PS_TIME_UTC);
     327    in2->sec = timesec;
     328    in2->nsec = 0;
     329    in2->leapsecond = false;
     330    in2 = psTimeConvert(in2, PS_TIME_TAI);
     331    psSphere *interm = p_psTimeGetPoleCoords(in2);
    288332    printf(" -- sphere = r=%.8g, d=%.8g\n", interm->r, interm->d);
    289333    printf(" -- Eop = x=%.8g, y=%.8g, s=%.8g\n", eop->x, eop->y, eop->s);
    290     printf(" -- PolarMotion = x=%.8g, y=%.8g, s=%.8g\n", polarMotion->x,
     334
     335
     336    psEarthPole *nutationCorr = psEOC_NutationCorr(in2);
     337    if (nutationCorr == NULL) {
     338        psError(PS_ERR_BAD_PARAMETER_NULL, false, "Nutation Correction returned NULL.\n");
     339        return 5;
     340    }
     341    printf("\n -- PolarMotion = x=%.13g, y=%.13g, s=%.13g\n", polarMotion->x,
    291342           polarMotion->y, polarMotion->s);
     343    printf("  -- NutationCorr = x,y,s = %.13g, %.13g, %.13g\n", nutationCorr->x,
     344           nutationCorr->y, nutationCorr->s);
     345    polarMotion->x += nutationCorr->x;
     346    polarMotion->y += nutationCorr->y;
     347    polarMotion->s += nutationCorr->s;
     348    printf("\n  -- Output Earthpole = x,y,s = %.13g, %.13g, %.13g\n", polarMotion->x,
     349           polarMotion->y, polarMotion->s);
     350    double x,y,s;
     351    x = -0.13275353774074533;
     352    y = 0.4359436319739848;
     353    s = -4.2376965863576153e-10;
     354    printf("  -- Expected Earthpole = x,y,s = %.13g, %.13g, %.13g\n", x, y, s);
     355
     356    if (!p_psEOCFinalize() ) {
     357        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "EOC failed finalization!\n");
     358        return 12;
     359    }
     360
     361    psFree(nutationCorr);
    292362    psFree(in);
     363    psFree(in2);
    293364    psFree(eop);
    294365    psFree(polarMotion);
    … …  
    345416psS32 testSphereRot_TEOtoCEO(void)
    346417{
    347     psTime *now = psTimeAlloc(PS_TIME_UT1);
    348     now->sec = 1128530000;
    349     now->nsec = 931154510;
    350     psSphereRot *teoceo = psSphereRot_TEOtoCEO(now);
     418    //    psTime *now = psTimeAlloc(PS_TIME_UT1);
     419    //    now->sec = 1128530000;
     420    //    now->nsec = 931154510;
     421    psTime *time = psTimeAlloc(PS_TIME_UTC);
     422    time->sec = timesec;
     423    time->nsec = 0;
     424    time->leapsecond = false;
     425    //    psSphereRot *teoceo = psSphereRot_TEOtoCEO(now);
     426    psSphereRot *teoceo = psSphereRot_TEOtoCEO(time);
    351427    printf("Output sphere rotation = %lf,%lf,%lf,%lf\n",
    352428           teoceo->q0, teoceo->q1, teoceo->q2, teoceo->q3);
    353429
    354     psFree(now);
     430    //    psFree(now);
     431    psFree(time);
    355432    psFree(teoceo);
    356433    return 0;
    … …  
    363440    psSphereRot *rot = NULL;
    364441
    365     in->x = M_PI / 4.0;
    366     in->y = M_PI / 6.0;
    367     in->s = M_PI / 8.0;
     442    //    in->x = M_PI / 4.0;
     443    //    in->y = M_PI / 6.0;
     444    //    in->s = M_PI / 8.0;
     445
     446    in->x = 2.857175590089105e-4;
     447    in->y = 2.3968739377734732e-5;
     448    in->s = -1.3970066457904322e-8;
     449    //    in->x += 0.06295703125;
     450    //    in->y += -0.0287618408203125;
     451    //    in->s += 0.0;
     452    double q0,q1,q2,q3;
     453    q0 = -1.1984522406756289e-5;
     454    q1 = 1.4285893358610674e-4;
     455    q2 = 1.2191193518914336e-10;
     456    q3 = -0.9999999897238481;
    368457
    369458    rot = psSphereRot_CEOtoGCRS(empty);
    … …  
    380469        return 2;
    381470    }
     471    printf("Output sphere rotation = %.13g,%.13g,%.13g,%.13g\n",
     472           rot->q0, rot->q1, rot->q2, rot->q3);
     473    printf("Expected sphere rotation = %.13g,%.13g,%.13g,%.13g\n", q0,q1,q2,q3);
    382474
    383475    psSphere *test = psSphereAlloc();
    … …  
    399491    psSphereRot *rot = NULL;
    400492
    401     in->x = M_PI / 4.0;
    402     in->y = M_PI / 6.0;
    403     in->s = M_PI / 8.0;
     493    //    in->x = M_PI / 4.0;
     494    //    in->y = M_PI / 6.0;
     495    //    in->s = M_PI / 8.0;
     496    in->x = -0.13275353774074533;
     497    in->y = 0.4359436319739848;
     498    in->s = -4.2376965863576153e-10;
    404499
    405500    rot = psSphereRot_ITRStoTEO(empty);
    … …  
    416511        return 2;
    417512    }
     513    printf("Output sphere rotation = %.13g,%.13g,%.13g,%.13g\n",
     514           rot->q0, rot->q1, rot->q2, rot->q3);
     515    double q0,q1,q2,q3;
     516    q0 = -1.0567571848664005e-6;
     517    q1 = 3.218036562931509e-7;
     518    q2 = -3.3580195807204483e-12;
     519    q3 = -0.9999999999993899;
     520    printf("Expected sphere rotation = %.13g,%.13g,%.13g,%.13g\n", q0,q1,q2,q3);
    418521
    419522    psSphere *test = psSphereAlloc();

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

    r5626 r5642  
    11
    2 apparent = r,d  = 0.65340071, 0.2463343
     2apparent = r,d  = 2.1452242, 0.84754694
     3  ><><Apparent = x,y,z = -0.3596051297757, 0.5556035294057, 0.7496591683899
    34
    45 <<<<SUN = r,d = 0.1745310859793, 0.0751383461293
    … …  
    1617 -- outCube = x,y,z = -0.3596195125763, 0.5555613903461, 0.7496834983719  --
    1718 -- outCube = x,y,z = 1, 1.888856707723e-13, -9.220762225745e-13  --
    18 Precession Model output = x,y,s = 1.1557751e+08, 8865375.4, -5.1231901e+14
     19Precession Model output = x,y,s = 58931469, 4942375.1, -1.4563071e+14
     20Expected output = x,y,s = 0.0002857175590089, 2.396873937773e-05, -1.39700664579e-08
    1921
    20 Precession Correction output (IERSA) = x,y,s = 0.26026728, 0.26026728, -0.61282361
    21 Precession Correction output (IERSB) = x,y,s = 0.26026728, 0.26026728, 0
    22  -- sphere = r=0.26026728, d=0.26026728
    23  -- Eop = x=-0.00045805944, y=0.00013368871, s=-2.986086e-05
    24  -- PolarMotion = x=0.25980922, y=0.26040097, s=-2.986086e-05
     22Precession Correction output (IERSA) = x,y,s = 0.25564256, 0.25564256, -0.34299142
     23Precession Correction output (IERSB) = x,y,s = 0.25564256, 0.25564256, -0.34298841
     24Expected output = x,y,s = 0.06295703125, -0.02876184082031, 0
     25
     26 -- sphere = r=0.25564256, d=0.25564256
     27 -- Eop = x=-0.00052201938, y=0.00015531693, s=-2.8841699e-05
     28
     29 -- PolarMotion = x=0.2551205394888, y=0.2557978757764, s=-0.3430172541992
     30  -- NutationCorr = x,y,s = 65064305.90552, -31402.7500932, 7.399292756179e-12
     31
     32  -- Output Earthpole = x,y,s = 65064306.16064, -31402.49429532, -0.3430172541918
     33  -- Expected Earthpole = x,y,s = -0.1327535377407, 0.435943631974, -4.237696586358e-10
    2534Nutation Correction output = x,y,s = 1.1740641e+08, -85355.474, 1.335034e-11
    2635
    2736
    28 The Value of T is = 1128530000.931
     37The Value of T is = 1049160599.657008886
    2938
    30 The Value of T is = 1126078455.931
    31 Output sphere rotation = 0.000000,0.000000,0.770554,0.637375
     39The Value of theta is = 6594672635.129
     40Output sphere rotation = 0.000000,0.000000,0.145165,0.989407
     41Output sphere rotation = 1.198437080989e-05,-0.0001428587808888,-6.985033157173e-09,0.9999999897239
     42Expected sphere rotation = -1.198452240676e-05,0.0001428589335861,1.219119351891e-10,-0.9999999897238
    3243
    33   Sphere -test- has values r,d = 0.048163027, 1.1836384
     44  Sphere -test- has values r,d = -1.7394212e-08, 0.00028571756
     45Output sphere rotation = 0.2157736658448,-0.06475858774092,-0.01434342985264,0.9741880293156
     46Expected sphere rotation = -1.056757184866e-06,3.218036562932e-07,-3.35801958072e-12,-0.9999999999994
    3447
    35   Sphere -test- has values r,d = 6.2758795, -0.85886273
     48  Sphere -test- has values r,d = 6.226856, 0.12027398
Note: See TracChangeset for help on using the changeset viewer.