IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jan 12, 2006, 10:40:13 AM (21 years ago)
Author:
drobbin
Message:

Updated fxns in EOC. Updated/expanded TableParser.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r5824 r5969  
    55*  @author d-Rob, MHPCC
    66*
    7 *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
    8 *  @date $Date: 2005-12-21 06:15:58 $
     7*  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
     8*  @date $Date: 2006-01-12 20:40:13 $
    99*
    1010*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    6464static psSphere *obj = NULL;
    6565static void objSetup(void);
     66static double greatCircle(psSphere *in1, psSphere *in2);
    6667
    6768static void objSetup(void)
    … …  
    8081}
    8182
     83//returns the great circle ANGULAR distance between 2 positions (psSphere's)
     84static double greatCircle(psSphere *in1, psSphere *in2)
     85{
     86    double r1, r2, d1, d2;
     87    r1 = in1->r;
     88    r2 = in2->r;
     89    d1 = in1->d;
     90    d2 = in2->d;
     91    double out = 0.0;
     92    double c1, c2, cd, s1, s2, sum, ac;
     93    c1 = cos(r1);
     94    c2 = cos(r2);
     95    cd = cos(d1-d2);
     96    s1 = sin(r1);
     97    s2 = sin(r2);
     98    sum = c1*c2*cd + s1*s2;
     99    if (sum > 1.0)
     100        sum = 1.0 - (sum - 1.0);
     101
     102    ac = acos(sum);
     103    printf("\n  c1=%lf, c2=%lf, cd=%lf, s1=%lf, s2=%lf, sum=%.19g, ac=%g\n", c1,c2,cd,s1,s2,sum,ac);
     104    //    out = acos(cos(r1)*cos(r2)*cos(d1-d2) + sin(r1)*sin(r2));
     105    //    if (out != 0.0) printf("\n   in greatCircle, out = %lf\n", out);
     106    return out;
     107}
     108
    82109psS32 testAberration(void)
    83110{
    84111
    85112    psSphere *apparent = NULL;
    86     psSphere *actual = psSphereAlloc();
     113    //    psSphere *actual = psSphereAlloc();
    87114    //    psSphere *direction = psSphereAlloc();
    88115    psSphere *empty = NULL;
    89116
    90     //    actual->r = 0.2;
    91     //    actual->d = 0.2;
    92     //    actual->r = DEG_TO_RAD(45.0);
    93     //    actual->d = DEG_TO_RAD(30.0);
    94     actual->r = DEG_TO_RAD(122.9153182445501);
    95     actual->d = DEG_TO_RAD(48.562968978679194);
    96     //    direction->r = 0.2035;
    97     //    direction->d = 0.2035;
    98     //    direction->r = DEG_TO_RAD(48.0);
    99     //    direction->r = DEG_TO_RAD(-157.0);
    100     //    direction->d = DEG_TO_RAD(20.7072);
     117    psCube *actualCube = psCubeAlloc();
     118    //    actual->r = DEG_TO_RAD(122.9153182445501);
     119    //    actual->d = DEG_TO_RAD(48.562968978679194);
     120    //    actualCube->x = -0.3596195125758298;
     121    //    actualCube->y = 0.5555613903455866;
     122    //    actualCube->z = 0.7496834983724809;
     123
     124    //values from After gravity deflection//
     125    actualCube->x = -0.35961949760293604;
     126    actualCube->y = 0.5555613950298085;
     127    actualCube->z = 0.7496835020836093;
     128    psSphere *actual = psCubeToSphere(actualCube);
     129
    101130    psCube *cubeDir = psCubeAlloc();
    102131    cubeDir->x = 5148.713262821658;
    … …  
    106135    cubeDir->y += 248.46429758174693;
    107136    cubeDir->z += 0.09694774143797581;
     137    double length = sqrt(cubeDir->x*cubeDir->x+cubeDir->y*cubeDir->y+cubeDir->z*cubeDir->z);
     138    cubeDir->x = cubeDir->x/length;
     139    cubeDir->y = cubeDir->y/length;
     140    cubeDir->z = cubeDir->z/length;
    108141    psSphere *direction = psCubeToSphere(cubeDir);
     142
    109143    double speed = sqrt(cubeDir->x*cubeDir->x + cubeDir->y*cubeDir->y + cubeDir->z*cubeDir->z);
    110144    // Speed of light in vacuum (src:NIST)
    111145    double c = 299792458.0; /* m/s */
    112     speed = speed / c;
     146    //    speed = speed / c;
     147    speed = length / c;
    113148
    114149    empty = psAberration(empty, apparent, direction, speed);
    … …  
    127162    apparent = psAberration(apparent, actual, direction, speed);
    128163
    129     printf("\nSphere Difference  =  r,d  =  %.13g, %.13g\n",
    130            (actual->r - apparent->r), (actual->d - apparent->d));
     164    //    printf("\nSphere Difference  =  r,d  =  %.13g, %.13g\n",
     165    //           (actual->r - apparent->r), (actual->d - apparent->d));
    131166    psCube *outCube = psSphereToCube(apparent);
    132167    printf(" -- resultCube = x,y,z = %.13g, %.13g, %.13g  -- \n",
    133168           outCube->x, outCube->y, outCube->z);
    134 
     169    //expected cube values
    135170    double x, y, z;
    136171    x = -0.35963388069046304;
    … …  
    138173    z = 0.7497078321908413;
    139174
    140     printf(" -- expectedCube = x,y,z = %.13g, %.13g, %.13g  -- \n\n", x, y, z);
    141 
    142     if ( fabs(x - outCube->x) > DBL_EPSILON || fabs(y - outCube->y) > DBL_EPSILON ||
    143             fabs(z - outCube->z) > DBL_EPSILON ) {
     175    printf(" -- expectedCube = x,y,z = %.13g, %.13g, %.13g  -- \n", x, y, z);
     176    printf("Cube Difference  =  x,y,z  = %.13g, %.13g, %.13g \n\n",
     177           (x - outCube->x), (y - outCube->y), (z - outCube->z) );
     178    psFree(actual);
     179    //        psFree(actualCube);
     180    actualCube->x = x;
     181    actualCube->y = y;
     182    actualCube->z = z;
     183    actual = psCubeToSphere(actualCube);
     184    double xxx = greatCircle(actual, apparent);
     185    printf("   The great circle angular distance between expected & result = %.13g\n",
     186           xxx);
     187
     188    if ( fabs(x - outCube->x) > FLT_EPSILON || fabs(y - outCube->y) > FLT_EPSILON ||
     189            fabs(z - outCube->z) > FLT_EPSILON ) {
    144190        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    145191                "psAberration returned incorrect values.\n");
    146         printf("Cube Difference  =  x,y,z  = %.13g, %.13g, %.13g \n\n",
    147                (x - outCube->x), (y - outCube->y), (z - outCube->z) );
     192        printf("\nMagnitude of expected change = x,y,z = %.13g, %.13g, %.13g \n",
     193               (actualCube->x - x), (actualCube->y - y), (actualCube->z - z) );
     194        printf("Magnitude of actual change = x,y,z = %.13g, %.13g, %.13g \n",
     195               (actualCube->x - outCube->x), (actualCube->y - outCube->y),
     196               (actualCube->z - outCube->z) );
     197        psFree(actual);
     198        //        psFree(actualCube);
     199        actualCube->x = x;
     200        actualCube->y = y;
     201        actualCube->z = z;
     202        actual = psCubeToSphere(actualCube);
     203        x = greatCircle(actual, apparent);
     204        printf("   The great circle angular distance between expected & result = %.13g\n",
     205               x);
    148206        return 3;
    149207    }
    150208
    151209    psFree(outCube);
     210    psFree(actualCube);
    152211
    153212    psFree(cubeDir);
    … …  
    162221{
    163222
    164     psSphere *actual = psSphereAlloc();
     223    //    psSphere *actual = psSphereAlloc();
    165224    psSphere *apparent = NULL;
    166225    //    psSphere *sun = psSphereAlloc();
    167226    psSphere *empty = NULL;
    168227
    169     //    sun->r = 0.2;
    170     //    sun->d = 0.2;
    171     //    actual->r = 0.61001;
    172     //    actual->d = 0.01999;
    173     actual->r = DEG_TO_RAD(122.9153182445501);
    174     actual->d = DEG_TO_RAD(48.562968978679194);
     228    psCube *actualCube = psCubeAlloc();
     229    actualCube->x = -0.3596195125758298;
     230    actualCube->y = 0.5555613903455866;
     231    actualCube->z = 0.7496834983724809;
     232    psSphere *actual = psCubeToSphere(actualCube);
    175233    psCube *sunCube = psCubeAlloc();
    176234    sunCube->x = 1.467797790127511e11;
    … …  
    193251
    194252    apparent = psGravityDeflection(NULL, actual, sun);
     253    psSphere *result2 = psSphereSetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN);
     254    apparent->r *= -1.0;
     255    apparent->d *= -1.0;
    195256    psSphere *result = psSphereSetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN);
    196 
    197     printf("\nActual r,d = %.13g,%.13g    Apparent r,d = %.16g, %.16g \n",
    198            actual->r, actual->d, result->r, result->d);
    199     printf("Sphere Difference  =  r,d  =  %.13g, %.13g\n",
    200            (actual->r - result->r), (actual->d - result->d));
     257    //    psSphere *result = psSphereGetOffset(apparent, actual, PS_SPHERICAL, PS_RADIAN);
     258    printf(" -- actualCube = x,y,z = %.13g, %.13g, %.13g  -- \n",
     259           actualCube->x, actualCube->y, actualCube->z);
    201260    psCube *outCube = psSphereToCube(result);
    202261    printf(" -- resultCube = x,y,z = %.13g, %.13g, %.13g  -- \n",
    203262           outCube->x, outCube->y, outCube->z);
    204     psCube *outCube2 = psSphereToCube(actual);
    205     printf(" -- actualCube = x,y,z = %.13g, %.13g, %.13g  -- \n",
     263    psCube *outCube2 = psSphereToCube(result2);
     264    printf(" -- resultCube2= x,y,z = %.13g, %.13g, %.13g  -- \n",
    206265           outCube2->x, outCube2->y, outCube2->z);
    207 
    208266    double x, y, z;
    209267    x = -0.35961949760293604;
    … …  
    211269    z = 0.7496835020836093;
    212270
    213     printf(" -- expectedCube = x,y,z = %.13g, %.13g, %.13g  -- \n\n", x, y, z);
    214 
    215     if ( fabs(x - outCube->x) > DBL_EPSILON || fabs(y - outCube->y) > DBL_EPSILON ||
    216             fabs(z - outCube->z) > DBL_EPSILON ) {
    217         psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    218                 "psGravityDeflection returned incorrect values.\n");
    219         printf("Cube Difference  =  x,y,z  = %.13g, %.13g, %.13g \n",
    220                (x - outCube->x), (y - outCube->y), (z - outCube->z) );
    221         return 1;
    222     }
     271    printf(" -- expectCube = x,y,z = %.13g, %.13g, %.13g  -- \n\n", x, y, z);
     272
     273    //    if ( fabs(x - outCube->x) > DBL_EPSILON || fabs(y - outCube->y) > DBL_EPSILON ||
     274    //            fabs(z - outCube->z) > DBL_EPSILON ) {
     275    //        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     276    //                "psGravityDeflection returned incorrect values.\n");
     277    printf("expect-actual=  x,y,z  = %.11g, %.11g, %.11g \n",
     278           (x - actualCube->x), (y - actualCube->y), (z - actualCube->z) );
     279    printf("result-actual=  x,y,z  = %.11g, %.11g, %.11g \n",
     280           (outCube->x - actualCube->x), (outCube->y - actualCube->y), (outCube->z - actualCube->z) );
     281    printf("expect-result=  x,y,z  = %.11g, %.11g, %.11g \n",
     282           (x - outCube->x), (y - outCube->y), (z - outCube->z) );
     283    //        return 1;
     284    //    }
     285    psFree(result2);
     286    outCube2->x = x;
     287    outCube2->y = y;
     288    outCube2->z = z;
     289    result2 = psCubeToSphere(outCube2);
     290    psFree(result);
     291    result = psSphereGetOffset(actual, result2, PS_SPHERICAL, PS_RADIAN);
     292    printf("The apparent output sphere = r,d = %.13g, %.13g\n", apparent->r, apparent->d);
     293    printf("The expected output sphere = r,d = %.13g, %.13g\n\n", result->r, result->d);
     294    psFree(result2);
    223295
    224296    psFree(outCube);
    225297    psFree(outCube2);
    226298    psFree(sunCube);
    227 
     299    psFree(actualCube);
    228300    psFree(result);
    229301    psFree(actual);
    … …  
    275347            psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    276348                    "   psEOC_PrecessionModel return incorrect values.\n");
    277             printf("\n  Precession Model output = x,y,s = %.8g, %.8g, %.8g\n",
    278                    pmodel->x, pmodel->y, pmodel->s);
    279             printf("  Expected output = x,y,s = %.13g, %.13g, %.13g\n", x, y, s);
    280             printf("  A difference of:   %.13g, %.13g, %.13g\n\n",
    281                    (pmodel->x - x), (pmodel->y - y), (pmodel->s - s) );
    282349            return 4;
    283350        }
     351        printf("\n  Precession Model output = x,y,s = %.8g, %.8g, %.8g\n",
     352               pmodel->x, pmodel->y, pmodel->s);
     353        printf("  Expected output = x,y,s = %.13g, %.13g, %.13g\n", x, y, s);
     354        printf("  A difference of:   %.13g, %.13g, %.13g\n\n",
     355               (pmodel->x - x), (pmodel->y - y), (pmodel->s - s) );
    284356    }
    285357    psFree(pmodel);
    … …  
    348420        yy = -0.0287618408203125;
    349421        ss = 0.0;
    350         xx = SEC_TO_RAD(xx) * 1e-6;
    351         yy = SEC_TO_RAD(yy) * 1e-6;
     422        xx = SEC_TO_RAD(xx) * 1e-3;
     423        yy = SEC_TO_RAD(yy) * 1e-3;
    352424        //        if ( fabs(pcorr->x - xx) > DBL_EPSILON || fabs(pcorr->y - yy) > DBL_EPSILON
    353425        //                || fabs(pcorr->s - ss) > DBL_EPSILON) {
    354426        //            psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    355427        //                    "   psEOC_PrecessionCorr return incorrect values.\n");
    356         printf("\nPrecession Correction output (IERSB) = x,y,s = %.13g, %.13g, %.13g\n",
     428        printf("PrecessionCorr output (IERSB) = x,y,s = %.13g, %.13g, %.13g\n",
    357429               pcorr->x, pcorr->y, pcorr->s);
    358         printf("Expected output = x,y,s = %.13g, %.13g, %.13g\n\n", xx, yy, ss);
    359         //            printf("  A difference of:   %.13g, %.13g, %.13g\n\n",
    360         //                   (pcorr->x - xx), (pcorr->y - yy), (pcorr->s - ss) );
     430        printf("Expected output               = x,y,s = %.13g, %.13g, %.13g\n\n", xx, yy, ss);
     431        printf("  A difference of:   %.13g, %.13g, %.13g\n\n",
     432               (pcorr->x - xx), (pcorr->y - yy), (pcorr->s - ss) );
    361433        //            return 10;
    362434        //        }
    363435    }
    364436
     437
     438    double xCorr, yCorr;
     439    xCorr = 3.05224300720406e-10;
     440    yCorr = -1.39441339235822e-10;
     441    //pcorr is the *expected* output from PrecessionModel//
    365442    pcorr->x = 2.857175590089105e-4;
    366443    pcorr->y = 2.3968739377734732e-5;
    367444    pcorr->s = -1.3970066457904322e-8;
     445    pcorr->x += xCorr;
     446    pcorr->y += yCorr;
     447    psEarthPole *precess = psEOC_PrecessionModel(time2);
     448    precess->x += xCorr;
     449    precess->y += yCorr;
     450    psSphereRot *precessNutInv = psSphereRot_CEOtoGCRS(precess);
     451    psSphereRot *precessNut = psSphereRotConjugate(NULL, precessNutInv);
    368452    //    pcorr->x += 3.05224300720406e-7;
    369453    //    pcorr->y += -1.39441339235822e-7;
    … …  
    377461        printf("\n Error at CEOtoGCRS, output psSphereRot doesn't match expected.\n");
    378462    }
    379     printf("\n  Output sphere rotation   = %.13g,%.13g,%.13g,%.13g\n",
     463    printf("  Output sphere rotation   = %.13g,%.13g,%.13g,%.13g\n",
    380464           pni->q0, pni->q1, pni->q2, pni->q3);
    381465    printf("  Expected sphere rotation = %.13g,%.13g,%.13g\n", q0,q1,q2);
     466    printf("  MY sphere rotation       = %.13g,%.13g,%.13g\n",
     467           precessNutInv->q0, precessNutInv->q1, precessNutInv->q2);
    382468    psCube *objC = psCubeAlloc();
    383469    //    objC->x = -3.5963388069046304;
    384470    //    objC->y = 0.5555192509816625;
    385471    //    objC->z = 0.7497078321908413;
    386     objSetup();
     472    //    objSetup();
     473
     474    //This is the sphere rotation for the *expected* precession output//
    387475    psSphereRot *pn = psSphereRotConjugate(NULL, pni);
    388     //    psSphere *sphere = psCubeToSphere(objC);
    389     psSphere *sphere = psSphereAlloc();
    390     *sphere = *obj;
    391     psFree(obj);
    392     psSphere *result = psSphereRotApply(NULL, pn, sphere);
    393     objC->x = -0.3598480726985338;
    394     objC->y = 0.5555012823608123;
    395     objC->z = 0.7496183628158023;
    396     psFree(sphere);
    397     sphere = psCubeToSphere(objC);
     476
     477    //    psSphere *sphere = psSphereAlloc();
     478    //    *sphere = *obj;
     479    //    psFree(obj);
     480
     481    //create a psSphere for (from) the start position given in eoc_testing//
     482    objC->x = -0.35963388069046304;
     483    objC->y = 0.5555192509816625;
     484    objC->z = 0.7497078321908413;
     485    psSphere *sphere = psCubeToSphere(objC);
     486
     487    psSphere *expect = psSphereRotApply(NULL, pn, sphere);
     488    //expected results below - stored in:  sphere  //
     489    double x,y,z;
     490    x = -0.3598480726985338;
     491    y = 0.5555012823608123;
     492    z = 0.7496183628158023;
     493    printf("\n<<Expected out      = x,y,z = %.13g, %.13g, %.13g\n", x, y, z);
    398494    psFree(objC);
    399     printf("\n Spheres:  out = %.13g, %.13g,    expect = %.13g, %.13g\n",
    400            result->r, result->d, sphere->r, sphere->d);
    401     double xx = acos(cos(result->r)*cos(sphere->r)*cos(result->d - sphere->d) + sin(result->r)*sin(sphere->r));
    402     printf("GREAT CIRCLE DIFFERENCE = %.13g \n\n", xx);
     495    objC = psSphereToCube(expect);
     496    printf("<<Expected out real = x,y,z = %.13g, %.13g, %.13g\n", objC->x, objC->y, objC->z);
     497    printf("     Difference     =      %.13g, %.13g, %.13g\n", objC->x-x, objC->y-y, objC->z-z);
     498    x = objC->x;
     499    y = objC->y;
     500    z = objC->z;
     501    psSphere *result = psSphereRotApply(NULL, precessNut, sphere);
     502    psFree(objC);
     503    objC = psSphereToCube(result);
     504    printf("<<Resulting out     = x,y,z = %.13g, %.13g, %.13g\n", objC->x, objC->y, objC->z);
     505    printf("     Difference     =      %.13g, %.13g, %.13g\n\n", objC->x-x, objC->y-y, objC->z-z);
     506
     507    double xx = greatCircle(result, expect);
     508    printf("GREAT CIRCLE DIFFERENCE = %.13g \n", xx);
     509    psFree(precess);
     510    psFree(precessNut);
     511    psFree(precessNutInv);
     512    psFree(expect);
     513    psFree(objC);
    403514
    404515    psFree(sphere);
    405516    psFree(result);
    406517    psFree(pn);
    407     //    psFree(obj);
    408518    psFree(pni);
    409 
    410519    psFree(pcorr);
    411 
    412520    psFree(time2);
    413521    if (!p_psEOCFinalize() ) {
    … …  
    443551    }
    444552
     553    //Test for IERS bulletin A.
     554    double x, y, s;
     555    x = -6.454389659777e-07;
     556    y = 2.112606414597e-06;
     557    s = 0.0;
     558    polarMotion = psEOC_GetPolarMotion(in, PS_IERS_A);
     559    if (polarMotion == NULL) {
     560        psError(PS_ERR_BAD_PARAMETER_NULL, false,
     561                "psEOC_GetPolarMotion returned NULL for valid input time.\n");
     562        return 4;
     563    }
     564    if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON
     565            || fabs(polarMotion->s - s) > DBL_EPSILON) {
     566        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     567                "psEOC_GetPolarMotion returned incorrect values.\n");
     568        printf("\n  <>PolarMotion output (IERSA)   = x,y,s = %.13g, %.13g, %.13g\n",
     569               polarMotion->x, polarMotion->y, polarMotion->s);
     570        printf("\n  <>PolarMotion expected (IERSA) = x,y,s = %.13g, %.13g, %.13g\n",
     571               x, y, s);
     572        //        return 5;
     573    }
     574    psFree(polarMotion);
     575
     576
    445577    //Return valid values for correct input time.  Test IERS Bulletin B.
    446578    polarMotion = psEOC_GetPolarMotion(in, PS_IERS_B);
    447     double x, y, s;
    448     x = -6.46014230078e-7;
    449     y = 2.11194535765e-6;
    450     s = -1.66256670849e-6;
     579    x = -6.45381397904e-07;
     580    y = 2.112819726698e-06;
     581    s = 0.0;
    451582
    452583    if (polarMotion == NULL) {
    … …  
    465596        //        return 3;
    466597    }
    467     psFree(polarMotion);
    468 
    469     //Test for IERS bulletin A.
    470     x = -6.46028774489e-7;
    471     y = 2.11172234336e-6;
    472     s = -1.66257785921e-6;
    473     polarMotion = psEOC_GetPolarMotion(in, PS_IERS_A);
    474     if (polarMotion == NULL) {
    475         psError(PS_ERR_BAD_PARAMETER_NULL, false,
    476                 "psEOC_GetPolarMotion returned NULL for valid input time.\n");
    477         return 4;
    478     }
    479     if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON
    480             || fabs(polarMotion->s - s) > DBL_EPSILON) {
    481         psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    482                 "psEOC_GetPolarMotion returned incorrect values.\n");
    483         printf("\n  <>PolarMotion output (IERSA)   = x,y,s = %.13g, %.13g, %.13g\n",
    484                polarMotion->x, polarMotion->y, polarMotion->s);
    485         printf("\n  <>PolarMotion expected (IERSA) = x,y,s = %.13g, %.13g, %.13g\n",
    486                x, y, s);
    487         //        return 5;
    488     }
    489 
    490598
    491599
    … …  
    495603        return 6;
    496604    }
    497     /*
    498         double sum = sqrt(nutationCorr->x*nutationCorr->x + nutationCorr->y*nutationCorr->y
    499                           + nutationCorr->s*nutationCorr->s);
    500         nutationCorr->x = nutationCorr->x / sum;
    501         nutationCorr->y = nutationCorr->y / sum;
    502         nutationCorr->s = nutationCorr->s / sum;
    503     */
    504     printf("  -- NutationCorr = x,y,s = %.13g, %.13g, %.13g\n\n", nutationCorr->x,
    505            nutationCorr->y, nutationCorr->s);
     605
    506606    polarMotion->x += nutationCorr->x;
    507607    polarMotion->y += nutationCorr->y;
    508608    polarMotion->s += nutationCorr->s;
    509     x = -0.13275353774074533;
    510     y = 0.4359436319739848;
    511     s = -4.2376965863576153e-10;
     609    x = -6.43607313124045e-7;
     610    y = 2.11351436973568e-6;
     611    s = -7.39617581324646e-12;
    512612
    513613    if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON
    … …  
    525625    }
    526626
     627    psEarthPole *polarTide = psEOC_PolarTideCorr(in);
     628    polarMotion->x += polarTide->x;
     629    polarMotion->y += polarTide->y;
     630    polarMotion->s += polarTide->s;
     631    psFree(polarTide);
     632    printf("\n  PolarMotion + PolarTideCorr out = x,y,s = %.13g, %.13g, %.13g\n",
     633           polarMotion->x, polarMotion->y, polarMotion->s);
     634
    527635    if (!p_psEOCFinalize() ) {
    528636        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "EOC failed finalization!\n");
    … …  
    602710        return 3;
    603711    } else {
    604         printf("Nutation Correction output = x,y,s = %.8g, %.8g, %.8g\n\n",
     712        printf("Nutation Correction output = x,y,s = %.13g, %.13g, %.13g\n\n",
    605713               nute->x, nute->y, nute->s);
    606714    }
    … …  
    656764    objSetup();
    657765    psSphereRot *earthRot = psSphereRotConjugate(NULL, teoceo);
    658     //    psSphere *result = psSphereRotApply(NULL, earthRot, obj);
    659     psSphere *result = psSphereRotApply(NULL, teoceo, obj);
     766    psSphere *result = psSphereRotApply(NULL, earthRot, obj);
     767    //    psSphere *result = psSphereRotApply(NULL, teoceo, obj);
    660768    psCube *cube = psSphereToCube(result);
    661769
    … …  
    773881    psSphere *expected = psCubeToSphere(expect);
    774882    psFree(expect);
    775     double d = (6.38e6) * acos(cos(result->r)*cos(expected->r)*cos(result->d-expected->d) + sin(result->r)*sin(expected->r));
    776     printf("\n\nGreat circle difference of:  %.13g \n", d);
     883    double d = acos(cos(result->r)*cos(expected->r)*cos(result->d-expected->d) + sin(result->r)*sin(expected->r));
     884    printf("Great circle difference of:  %.13g \n", d);
    777885
    778886    psFree(expected);
    … …  
    781889
    782890
    783     psSphere *test = psSphereAlloc();
    784     test->r = 0.0;
    785     test->d = 0.0;
    786     test = psSphereRotApply(test, rot, test);
    787     printf("\n  Sphere -test- has values r,d = %.8g, %.8g \n", test->r, test->d);
     891    //    psSphere *test = psSphereAlloc();
     892    //    test->r = 0.0;
     893    //    test->d = 0.0;
     894    //    test = psSphereRotApply(test, rot, test);
     895    //    printf("\n  Sphere -test- has values r,d = %.8g, %.8g \n", test->r, test->d);
    788896
    789897    psFree(precessionNutation);
    … …  
    791899    psFree(cube);
    792900
    793     psFree(test);
     901    //    psFree(test);
    794902    psFree(rot);
    795903    psFree(in);
    … …  
    875983    z = 0.7496169753347885;
    876984    printf("\n  Cube -expected- has x,y,z = %.13g,  %.13g, %.13g \n", x, y, z );
     985    temp->x = x;
     986    temp->y = y;
     987    temp->z = z;
     988    psSphere *sphere = psCubeToSphere(temp);
     989    double d = greatCircle(sphere, test);
     990
     991    printf("Great circle difference of:  %.13g \n", d);
     992    psFree(sphere);
    877993
    878994    psFree(newRot);
Note: See TracChangeset for help on using the changeset viewer.