IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jan 17, 2006, 2:41:29 PM (21 years ago)
Author:
drobbin
Message:

Updated Table parser, Fixed misc EOC, worked primarily on GravityDeflection

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r5969 r6030  
    55*  @author d-Rob, MHPCC
    66*
    7 *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
    8 *  @date $Date: 2006-01-12 20:40:13 $
     7*  @version $Revision: 1.26 $ $Name: not supported by cvs2svn $
     8*  @date $Date: 2006-01-18 00:41:29 $
    99*
    1010*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    8585{
    8686    double r1, r2, d1, d2;
    87     r1 = in1->r;
    88     r2 = in2->r;
    89     d1 = in1->d;
    90     d2 = in2->d;
     87    d1 = in1->r;
     88    d2 = in2->r;
     89    r1 = in1->d;
     90    r2 = in2->d;
    9191    double out = 0.0;
    9292    double c1, c2, cd, s1, s2, sum, ac;
    … …  
    101101
    102102    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);
     103    out = ac;
     104    //    printf("\n  c1=%lf, c2=%lf, cd=%lf, s1=%lf, s2=%lf, sum=%.19g, ac=%g\n", c1,c2,cd,s1,s2,sum,ac);
    104105    //    out = acos(cos(r1)*cos(r2)*cos(d1-d2) + sin(r1)*sin(r2));
    105106    //    if (out != 0.0) printf("\n   in greatCircle, out = %lf\n", out);
    … …  
    109110psS32 testAberration(void)
    110111{
    111 
    112112    psSphere *apparent = NULL;
    113     //    psSphere *actual = psSphereAlloc();
    114     //    psSphere *direction = psSphereAlloc();
    115113    psSphere *empty = NULL;
    116 
    117114    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 
    124115    //values from After gravity deflection//
    125116    actualCube->x = -0.35961949760293604;
    … …  
    127118    actualCube->z = 0.7496835020836093;
    128119    psSphere *actual = psCubeToSphere(actualCube);
    129 
    130120    psCube *cubeDir = psCubeAlloc();
    131121    cubeDir->x = 5148.713262821658;
    … …  
    135125    cubeDir->y += 248.46429758174693;
    136126    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;
    141127    psSphere *direction = psCubeToSphere(cubeDir);
    142 
    143128    double speed = sqrt(cubeDir->x*cubeDir->x + cubeDir->y*cubeDir->y + cubeDir->z*cubeDir->z);
    144     // Speed of light in vacuum (src:NIST)
    145     double c = 299792458.0; /* m/s */
    146     //    speed = speed / c;
    147     speed = length / c;
     129    double c = 299792458.0; // Speed of light in vacuum (src:NIST)   /* m/s */
     130    speed /= c;
    148131
    149132    empty = psAberration(empty, apparent, direction, speed);
    … …  
    161144
    162145    apparent = psAberration(apparent, actual, direction, speed);
    163 
    164     //    printf("\nSphere Difference  =  r,d  =  %.13g, %.13g\n",
    165     //           (actual->r - apparent->r), (actual->d - apparent->d));
    166146    psCube *outCube = psSphereToCube(apparent);
    167     printf(" -- resultCube = x,y,z = %.13g, %.13g, %.13g  -- \n",
     147    printf("\n -- resultCube = x,y,z =     %.13g,      %.13g,    %.13g  -- \n",
    168148           outCube->x, outCube->y, outCube->z);
    169149    //expected cube values
    … …  
    173153    z = 0.7497078321908413;
    174154
    175     printf(" -- expectedCube = x,y,z = %.13g, %.13g, %.13g  -- \n", x, y, z);
     155    printf(" -- expectedCube = x,y,z =   %.13g,     %.13g,    %.13g  -- \n", x, y, z);
    176156    printf("Cube Difference  =  x,y,z  = %.13g, %.13g, %.13g \n\n",
    177157           (x - outCube->x), (y - outCube->y), (z - outCube->z) );
    178158    psFree(actual);
    179     //        psFree(actualCube);
    180159    actualCube->x = x;
    181160    actualCube->y = y;
    … …  
    190169        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    191170                "psAberration returned incorrect values.\n");
    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) );
     171        //        printf("\nMagnitude of expected change = x,y,z = %.13g, %.13g, %.13g \n",
     172        //               (actualCube->x - x), (actualCube->y - y), (actualCube->z - z) );
     173        //        printf("Magnitude of actual change = x,y,z = %.13g, %.13g, %.13g \n",
     174        //               (actualCube->x - outCube->x), (actualCube->y - outCube->y),
     175        //               (actualCube->z - outCube->z) );
    197176        psFree(actual);
    198         //        psFree(actualCube);
    199177        actualCube->x = x;
    200178        actualCube->y = y;
    … …  
    209187    psFree(outCube);
    210188    psFree(actualCube);
    211 
    212189    psFree(cubeDir);
    213190    psFree(apparent);
    … …  
    220197psS32 testGravityDeflect(void)
    221198{
    222 
    223     //    psSphere *actual = psSphereAlloc();
    224199    psSphere *apparent = NULL;
    225     //    psSphere *sun = psSphereAlloc();
    226200    psSphere *empty = NULL;
    227 
    228201    psCube *actualCube = psCubeAlloc();
    229202    actualCube->x = -0.3596195125758298;
    … …  
    235208    sunCube->y = 2.5880956908748722e10;
    236209    sunCube->z = 1.1220046291457653e10;
     210    double sunLength = sqrt(sunCube->x*sunCube->x + sunCube->y*sunCube->y + sunCube->z*sunCube->z);
     211    sunCube->x /= sunLength;
     212    sunCube->y /= sunLength;
     213    sunCube->z /= sunLength;
     214    printf("sunCube = x,y,z = %.13g, %.13g, %.13g\n", sunCube->x, sunCube->y, sunCube->z);
    237215    psSphere *sun = psCubeToSphere(sunCube);
     216    printf("sunSphere  = r, d = %.13g, %.13g\n", sun->r, sun->d);
     217    psCube *outCube = psCubeAlloc();
    238218
    239219    empty = psGravityDeflection(apparent, empty, sun);
    … …  
    251231
    252232    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;
    256     psSphere *result = psSphereSetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN);
     233    //    apparent->r *= -1.0;
     234    //    apparent->d *= -1.0;
     235    psSphere *result2;
     236    //    psSphere *result2 = psSphereSetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN);
     237    //    psSphere *result = psSphereSetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN);
    257238    //    psSphere *result = psSphereGetOffset(apparent, actual, PS_SPHERICAL, PS_RADIAN);
    258239    printf(" -- actualCube = x,y,z = %.13g, %.13g, %.13g  -- \n",
    259240           actualCube->x, actualCube->y, actualCube->z);
    260     psCube *outCube = psSphereToCube(result);
    261     printf(" -- resultCube = x,y,z = %.13g, %.13g, %.13g  -- \n",
    262            outCube->x, outCube->y, outCube->z);
    263     psCube *outCube2 = psSphereToCube(result2);
     241    //    psCube *outCube = psSphereToCube(result);
     242    //    printf(" -- resultCube = x,y,z = %.13g, %.13g, %.13g  -- \n",
     243    //           outCube->x, outCube->y, outCube->z);
     244    psCube *outCube2 = psSphereToCube(apparent);
    264245    printf(" -- resultCube2= x,y,z = %.13g, %.13g, %.13g  -- \n",
    265246           outCube2->x, outCube2->y, outCube2->z);
    … …  
    275256    //        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    276257    //                "psGravityDeflection returned incorrect values.\n");
    277     printf("expect-actual=  x,y,z  = %.11g, %.11g, %.11g \n",
     258    printf("expect-actual=  x,y,z  = %.13g, %.13g, %.13g \n",
    278259           (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) );
     260    printf("expect-result=  x,y,z  = %.13g, %.13g, %.13g \n",
     261           (x - outCube2->x), (y - outCube2->y), (z - outCube2->z) );
     262    printf("result-actual=  x,y,z  = %.13g, %.13g, %.13g \n",
     263           (outCube2->x - actualCube->x), (outCube2->y - actualCube->y),
     264           (outCube2->z - actualCube->z) );
    283265    //        return 1;
    284266    //    }
    285     psFree(result2);
     267    //    psFree(result2);
    286268    outCube2->x = x;
    287269    outCube2->y = y;
    288270    outCube2->z = z;
    289271    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);
     272    //    psFree(result);
     273    psSphere *result = psSphereGetOffset(actual, result2, PS_SPHERICAL, PS_RADIAN);
     274    psFree(result2);
     275    result2 = psSphereGetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN);
     276    printf("The apparent output sphere = r,d = %.13g, %.13g\n", result2->r, result2->d);
    293277    printf("The expected output sphere = r,d = %.13g, %.13g\n\n", result->r, result->d);
    294278    psFree(result2);
    … …  
    349333            return 4;
    350334        }
    351         printf("\n  Precession Model output = x,y,s = %.8g, %.8g, %.8g\n",
     335        printf("\n  PrecessionModel output = x,y,s = %.13g,  %.13g,  %.13g\n",
    352336               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",
     337        printf("  Expected output        = x,y,s = %.13g, %.13g, %.13g\n", x, y, s);
     338        printf("  A difference of:                 %.13g, %.13g, %.13g\n",
    355339               (pmodel->x - x), (pmodel->y - y), (pmodel->s - s) );
    356340    }
    … …  
    402386        return 7;
    403387    } else {
    404         printf("\nPrecession Correction output (IERSA) = x,y,s = %.13g, %.13g, %.13g\n\n",
     388        printf("\nPrecessionCorr output (IERSA) = x,y,s = %.13g,  %.13g, %.13g\n",
    405389               pcorr->x, pcorr->y, pcorr->s);
    406390    }
    407391    psFree(pcorr);
    408 
    409     //    printf("\n >>MJD value = %lf \n", psTimeToMJD(time2));
    410392
    411393    //Check values from IERS table B
    … …  
    428410        printf("PrecessionCorr output (IERSB) = x,y,s = %.13g, %.13g, %.13g\n",
    429411               pcorr->x, pcorr->y, pcorr->s);
    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",
     412        printf("Expected output               = x,y,s = %.13g, %.13g, %.13g\n", xx, yy, ss);
     413        printf("          A difference of:              %.13g, %.13g, %.13g\n\n",
    432414               (pcorr->x - xx), (pcorr->y - yy), (pcorr->s - ss) );
    433415        //            return 10;
    434416        //        }
    435417    }
    436 
    437 
     418    //precess is the *actual* output from PrecessionModel + PrecessionCorr
     419    psEarthPole *precess = psEOC_PrecessionModel(time2);
     420    precess->x += pcorr->x;
     421    precess->y += pcorr->y;
    438422    double xCorr, yCorr;
    439423    xCorr = 3.05224300720406e-10;
    … …  
    445429    pcorr->x += xCorr;
    446430    pcorr->y += yCorr;
    447     psEarthPole *precess = psEOC_PrecessionModel(time2);
    448     precess->x += xCorr;
    449     precess->y += yCorr;
    450431    psSphereRot *precessNutInv = psSphereRot_CEOtoGCRS(precess);
    451432    psSphereRot *precessNut = psSphereRotConjugate(NULL, precessNutInv);
    452     //    pcorr->x += 3.05224300720406e-7;
    453     //    pcorr->y += -1.39441339235822e-7;
    454433    double q0, q1, q2;
    455434    q0 = -1.1984522406756289e-5;
    … …  
    461440        printf("\n Error at CEOtoGCRS, output psSphereRot doesn't match expected.\n");
    462441    }
    463     printf("  Output sphere rotation   = %.13g,%.13g,%.13g,%.13g\n",
    464            pni->q0, pni->q1, pni->q2, pni->q3);
    465     printf("  Expected sphere rotation = %.13g,%.13g,%.13g\n", q0,q1,q2);
    466     printf("  MY sphere rotation       = %.13g,%.13g,%.13g\n",
     442    //    printf("  Output from CEOtoGCRS only  = %.13g,%.13g,%.13g,%.13g\n",
     443    //           pni->q0, pni->q1, pni->q2, pni->q3);
     444    printf("  Expected sphere rotation    = %.13g, %.13g, %.13g\n", q0,q1,q2);
     445    printf("  Result sphere rotation      = %.13g, %.13g, %.13g\n",
    467446           precessNutInv->q0, precessNutInv->q1, precessNutInv->q2);
     447    printf("     Difference         =        %.13g, %.13g, %.13g\n\n",
     448           precessNutInv->q0-q0, precessNutInv->q1-q1, precessNutInv->q2-q2);
    468449    psCube *objC = psCubeAlloc();
    469450    //    objC->x = -3.5963388069046304;
    … …  
    491472    y = 0.5555012823608123;
    492473    z = 0.7496183628158023;
    493     printf("\n<<Expected out      = x,y,z = %.13g, %.13g, %.13g\n", x, y, z);
    494     psFree(objC);
    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;
     474    printf("\n<<Expected out       = x,y,z = %.13g, %.13g, %.13g\n", x, y, z);
     475    //    psFree(objC);
     476    //    objC = psSphereToCube(expect);
     477    //printf("<<Expected out (CEO) = x,y,z = %.13g, %.13g, %.13g\n", objC->x, objC->y, objC->z);
     478    //printf("     Difference     =           %.13g, %.13g, %.13g\n", objC->x-x, objC->y-y, objC->z-z);
     479    //    x = objC->x;
     480    //    y = objC->y;
     481    //    z = objC->z;
    501482    psSphere *result = psSphereRotApply(NULL, precessNut, sphere);
    502483    psFree(objC);
    503484    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);
     485    printf("<<Resulting out      = x,y,z = %.13g, %.13g, %.13g\n", objC->x, objC->y, objC->z);
     486    printf("     Difference         =      %.13g, %.13g, %.13g\n\n", objC->x-x, objC->y-y, objC->z-z);
    506487
    507488    double xx = greatCircle(result, expect);
    … …  
    562543        return 4;
    563544    }
    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     }
     545    //    if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON
     546    //            || fabs(polarMotion->s - s) > DBL_EPSILON) {
     547    //        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     548    //                "psEOC_GetPolarMotion returned incorrect values.\n");
     549    printf("  <>PolarMotion output (IERSA)   = x,y,s = %.13g, %.13g, %.13g\n",
     550           polarMotion->x, polarMotion->y, polarMotion->s);
     551    //        printf("  <>PolarMotion expected (IERSA) = x,y,s = %.13g, %.13g, %.13g\n",
     552    //               x, y, s);
     553    //        return 5;
     554    //    }
    574555    psFree(polarMotion);
    575 
    576556
    577557    //Return valid values for correct input time.  Test IERS Bulletin B.
    578558    polarMotion = psEOC_GetPolarMotion(in, PS_IERS_B);
    579     x = -6.45381397904e-07;
    580     y = 2.112819726698e-06;
    581     s = 0.0;
    582 
     559    //    x = -6.45381397904e-07;
     560    //    y = 2.112819726698e-06;
     561    //    s = 0.0;
    583562    if (polarMotion == NULL) {
    584563        psError(PS_ERR_BAD_PARAMETER_NULL, false,
    … …  
    586565        return 2;
    587566    }
    588     if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON
    589             || fabs(polarMotion->s - s) > DBL_EPSILON) {
    590         psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    591                 "psEOC_GetPolarMotion returned incorrect values.\n");
    592         printf("\n  <>PolarMotion output (IERSB)   = x,y,s = %.13g, %.13g, %.13g\n",
    593                polarMotion->x, polarMotion->y, polarMotion->s);
    594         printf("\n  <>PolarMotion expected (IERSB) = x,y,s = %.13g, %.13g, %.13g\n",
    595                x, y, s);
    596         //        return 3;
    597     }
    598 
    599 
     567    //    if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON
     568    //            || fabs(polarMotion->s - s) > DBL_EPSILON) {
     569    //        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     570    //                "psEOC_GetPolarMotion returned incorrect values.\n");
     571    printf("  <>PolarMotion output (IERSB)   = x,y,s = %.13g,  %.13g, %.13g\n",
     572           polarMotion->x, polarMotion->y, polarMotion->s);
     573    //        printf("  <>PolarMotion expected (IERSB) = x,y,s = %.13g, %.13g, %.13g\n",
     574    //               x, y, s);
     575    //        return 3;
     576    //    }
    600577    psEarthPole *nutationCorr = psEOC_NutationCorr(in);
    601578    if (nutationCorr == NULL) {
    … …  
    603580        return 6;
    604581    }
    605 
    606582    polarMotion->x += nutationCorr->x;
    607583    polarMotion->y += nutationCorr->y;
    608584    polarMotion->s += nutationCorr->s;
     585    psEarthPole *polarTide = psEOC_PolarTideCorr(in);
     586    polarMotion->x += polarTide->x;
     587    polarMotion->y += polarTide->y;
     588    psFree(polarTide);
    609589    x = -6.43607313124045e-7;
    610590    y = 2.11351436973568e-6;
    611591    s = -7.39617581324646e-12;
    612 
    613     if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON
    614             || fabs(polarMotion->s - s) > DBL_EPSILON) {
    615         psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    616                 "   psEOC_GetPolarMotion returned incorrect values.\n");
    617         printf("\n  PolarMotion + NutationCorr out = x,y,s = %.13g, %.13g, %.13g\n",
    618                polarMotion->x, polarMotion->y, polarMotion->s);
    619         printf("  Expected output = x,y,s = %.13g, %.13g, %.13g\n", x, y, s);
    620         printf("\n  Output Nutation = x,y,s = %.13g, %.13g, %.13g\n",
    621                nutationCorr->x, nutationCorr->y, nutationCorr->s);
    622         //        printf("  A difference of:   %.13g, %.13g, %.13g\n\n",
    623         //               (polarMotion->x - x), (polarMotion->y - y), (polarMotion->s - s) );
    624         //        return 10;
    625     }
    626 
    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",
     592    //    if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON
     593    //            || fabs(polarMotion->s - s) > DBL_EPSILON) {
     594    //        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     595    //                "   psEOC_GetPolarMotion returned incorrect values.\n");
     596    printf("\n  PolarMotion + NutationCorr out = x,y,s = %.13g, %.13g, %.13g\n",
    633597           polarMotion->x, polarMotion->y, polarMotion->s);
     598    printf("  Expected output                = x,y,s = %.13g,  %.13g, %.13g\n", x, y, s);
     599    printf("                     Difference  = x,y,s = %.13g, %.13g,  %.13g\n",
     600           polarMotion->x - x, polarMotion->y - y, polarMotion->s - s);
     601    //    }
    634602
    635603    if (!p_psEOCFinalize() ) {
    … …  
    679647{
    680648    psTime *in = psTimeAlloc(PS_TIME_UTC);
    681     //    in->sec = 1131579114;
    682     //    in->nsec = 498489000;
    683649    in->sec = timesec;
    684650    in->nsec = 0;
    … …  
    728694    psSphereRot *rot = NULL;
    729695    psTime *empty = NULL;
    730     //    psTime *now = psTimeAlloc(PS_TIME_UT1);
    731     //    now->sec = 1128530000;
    732     //    now->nsec = 931154510;
    733696    psTime *time = psTimeAlloc(PS_TIME_UT1);
    734697    time->sec = timesec-1;
    … …  
    736699    time->leapsecond = false;
    737700
    738     //    psSphereRot *teoceo = psSphereRot_TEOtoCEO(now);
    739 
    740701    //return NULL for NULL input time
    741     rot = psSphereRot_TEOtoCEO(empty);
     702    rot = psSphereRot_TEOtoCEO(empty, NULL);
    742703    if (rot != NULL) {
    743704        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    … …  
    746707    }
    747708
    748     psSphereRot *teoceo = psSphereRot_TEOtoCEO(time);
     709    psEarthPole *polarTideCorr = psEOC_PolarTideCorr(time);
     710    psSphereRot *teoceo = psSphereRot_TEOtoCEO(time, polarTideCorr);
    749711    //Make sure values match for other psTime type
    750712    empty = psTimeAlloc(PS_TIME_UTC);
    … …  
    752714    empty->nsec = 0;
    753715    empty->leapsecond = false;
    754     rot = psSphereRot_TEOtoCEO(empty);
     716
     717    rot = psSphereRot_TEOtoCEO(empty, polarTideCorr);
    755718    if (fabs(rot->q0-teoceo->q0) > DBL_EPSILON || fabs(rot->q1-teoceo->q1) > DBL_EPSILON ||
    756719            fabs(rot->q2-teoceo->q2) > DBL_EPSILON || fabs(rot->q3-teoceo->q3) > DBL_EPSILON) {
    757720        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    758721                "psSphereRot_TEOtoCEO failed to return matching values for different time types.\n");
     722        printf("\n  Output Rotation1 = q0,q1,q2,q3 = %.13g, %.13g, %.13g, %.13g\n",
     723               teoceo->q0, teoceo->q1, teoceo->q2, teoceo->q3 );
     724        printf("\n  Output Rotation2 = q0,q1,q2,q3 = %.13g, %.13g, %.13g, %.13g\n",
     725               rot->q0, rot->q1, rot->q2, rot->q3 );
    759726        return 2;
    760727    }
    … …  
    783750        return 3;
    784751    }
    785     //    psFree(now);
     752
     753    psFree(polarTideCorr);
    786754    psFree(rot);
    787755    psFree(empty);
    … …  
    800768    psEarthPole *empty = NULL;
    801769    psSphereRot *rot = NULL;
    802 
    803     //    in->x = M_PI / 4.0;
    804     //    in->y = M_PI / 6.0;
    805     //    in->s = M_PI / 8.0;
    806 
    807770    in->x = 2.857175590089105e-4;
    808771    in->y = 2.3968739377734732e-5;
    809772    in->s = -1.3970066457904322e-8;
    810     //    in->x += 0.06295703125;
    811     //    in->y += -0.0287618408203125;
    812     //    in->s += 0.0;
    813     //    in->x += 3.05224300720406e-7;
    814     //    in->y += -1.39441339235822e-7;
    815     in->s += 0.0;
    816     //    in->x += -2.491942320903e-10;
    817     //    in->y += -1.648366515772e-11;
    818773
    819774    double q0,q1,q2,q3;
    … …  
    837792    }
    838793
    839     if (fabs(rot->q0-q0) > DBL_EPSILON || fabs(rot->q1-q1) > DBL_EPSILON ||
    840             fabs(rot->q2-q2) > DBL_EPSILON || fabs(rot->q3-q3) > DBL_EPSILON) {
     794    printf("\n  Output sphere rotation   = %.13g, %.13g, %.13g, %.13g\n",
     795           rot->q0, rot->q1, rot->q2, rot->q3);
     796    printf("  Expected sphere rotation = %.13g, %.13g, %.13g, %.13g\n", q0,q1,q2,q3);
     797    printf("  difference:                 %.13g, %.13g, %.13g \n",
     798           (rot->q0-q0), (rot->q1-q1), (rot->q2-q2) );
     799    if (fabs(rot->q0-q0) > FLT_EPSILON || fabs(rot->q1-q1) > FLT_EPSILON ||
     800            fabs(rot->q2-q2) > FLT_EPSILON || fabs(rot->q3+q3) > FLT_EPSILON) {
    841801        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    842802                "psSphereRot_CEOtoGCRS failed to return expected values.\n");
    843         printf("\n  Output sphere rotation   = %.13g,%.13g,%.13g,%.13g\n",
    844                rot->q0, rot->q1, rot->q2, rot->q3);
    845         printf("  Expected sphere rotation = %.13g,%.13g,%.13g,%.13g\n", q0,q1,q2,q3);
    846         printf("  a difference:   %.13g, %.13g, %.13g, %.13g \n", (rot->q0-q0), (rot->q1-q1),
    847                (rot->q2-q2), (rot->q3-q3) );
    848         printf("  abs difference:   %.13g, %.13g, %.13g, %.13g \n\n", (rot->q0+q0), (rot->q1+q1),
    849                (rot->q2+q2), (rot->q3+q3) );
    850         //return 3;
    851     }
    852 
    853 
     803        return 3;
     804    }
    854805    psCube *tempCube = psCubeAlloc();
    855806    tempCube->x = -0.35963388069046304;
    … …  
    858809    obj = psCubeToSphere(tempCube);
    859810    psFree(tempCube);
    860     //    obj = psSphereAlloc();
    861     //    obj->r = objR;
    862     //    obj->d = objD;
    863811    psSphereRot *precessionNutation = psSphereRotConjugate(NULL, rot);
    864812    psSphere *result = psSphereRotApply(NULL, precessionNutation, obj);
    … …  
    868816    y = 0.5555012823608123;
    869817    z = 0.7496183628158023;
    870     printf("\nOutput cube = x,y,z = %.13g, %.13g, %.13g\n", cube->x, cube->y, cube->z);
    871     printf("Expected cube = x,y,z = %.13g, %.13g, %.13g\n", x, y, z);
    872     printf("A difference of:   %.13g, %.13g, %.13g\n\n",
     818    printf("\n  Output cube = x,y,z = %.13g,  %.13g,    %.13g\n", cube->x, cube->y, cube->z);
     819    printf("Expected cube = x,y,z = %.13g,  %.13g,  %.13g\n", x, y, z);
     820    printf("  A difference of:    %.13g, %.13g, %.13g\n\n",
    873821           (x-cube->x), (y-cube->y), (z-cube->z));
    874 
    875     //great sphere difference between result and expected
    876     //Re = 6.38x10^6m
    877822    psCube *expect = psCubeAlloc();
    878823    expect->x = x;
    … …  
    880825    expect->z = z;
    881826    psSphere *expected = psCubeToSphere(expect);
     827    double d = greatCircle(result, expected);
     828    printf("   The great circle angular distance between expected & result = %.13g\n", d);
     829
    882830    psFree(expect);
    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);
    885 
    886831    psFree(expected);
    887832    psFree(obj);
    888 
    889 
    890 
    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);
    896 
    897833    psFree(precessionNutation);
    898834    psFree(result);
    899835    psFree(cube);
    900 
    901     //    psFree(test);
    902836    psFree(rot);
    903837    psFree(in);
     838
    904839    return 0;
    905840}
    … …  
    910845    psEarthPole *empty = NULL;
    911846    psSphereRot *rot = NULL;
    912 
    913     //    in->x = M_PI / 4.0;
    914     //    in->y = M_PI / 6.0;
    915     //    in->s = M_PI / 8.0;
    916847    in->x = -0.13275353774074533;
    917848    in->y = 0.4359436319739848;
    918849    in->s = -4.2376965863576153e-10;
    919 
    920850    in->x = SEC_TO_RAD(in->x);
    921851    in->y = SEC_TO_RAD(in->y);
    … …  
    928858        return 1;
    929859    }
    930 
    931860    rot = psSphereRot_ITRStoTEO(in);
    932861    if (rot == NULL) {
    … …  
    958887    printf("  Expected sphere rotation = %.13g, %.13g, %.13g, %.13g\n", q0,q1,q2,q3);
    959888
    960     //    psSphere *test = psSphereAlloc();
    961     //    test->r = 0.0;
    962     //    test->d = 0.0;
    963889    psCube *temp = psCubeAlloc();
    964890    //    temp->x = -0.3596195125758298;
    … …  
    969895    temp->z = 0.7496183628158023;
    970896
    971 
    972897    obj = psCubeToSphere(temp);
    973898    psSphere *test = NULL;
    … …  
    975900    psSphereRot *newRot = psSphereRotConjugate(NULL, rot);
    976901    test = psSphereRotApply(NULL, newRot, obj);
    977     printf("\n  Sphere -test- has values r,d = %.8g, %.8g \n", test->r, test->d);
    978902    temp = psSphereToCube(test);
    979903    printf("\n  Cube -test- has x,y,z =     %.13g, %.13g, %.13g \n", temp->x, temp->y, temp->z);
    … …  
    988912    psSphere *sphere = psCubeToSphere(temp);
    989913    double d = greatCircle(sphere, test);
    990 
    991914    printf("Great circle difference of:  %.13g \n", d);
     915
    992916    psFree(sphere);
    993 
    994917    psFree(newRot);
    995918    psFree(obj);
    … …  
    998921    psFree(rot);
    999922    psFree(in);
    1000     return 0;
    1001 }
    1002 
     923
     924    return 0;
     925}
     926
Note: See TracChangeset for help on using the changeset viewer.