/** @file  tst_psEarthOrientation.c
*
*  @brief The code will perform earth orientation calculations and sphere rotations.
*
*  @author d-Rob, MHPCC
*
*  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
*  @date $Date: 2006-05-26 01:23:41 $
*
*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
*/

#include "psTest.h"
#include "pslib_strict.h"


static psS32 testAberration(void);
static psS32 testGravityDeflect(void);
static psS32 testEOCPrecession(void);
static psS32 testEOCPrecessionCorr(void);
static psS32 testEOCPolar(void);
static psS32 testEOCPolarTide(void);
static psS32 testEOCNutation(void);
static psS32 testSphereRot_TEOtoCEO(void);
static psS32 testSphereRot_CEOtoGCRS(void);
static psS32 testSphereRot_ITRStoTEO(void);
static psS32 testSphereRotPrecess(void);

testDescription tests[] = {
                              {testAberration, 666, "psAberration()", 0, false},
                              {testGravityDeflect, 667, "psGravityDeflect()", 0, false},
                              {testEOCPrecession, 669, "psEOCPrecession()", 0, false},
                              {testEOCPrecessionCorr, 670, "psEOCPrecessionCorr()", 0, false},
                              {testEOCPolar, 671, "psEOCPolar()", 0, false},
                              {testEOCPolarTide, 672, "psEOCPolarTide()", 0, false},
                              {testEOCNutation, 673, "psEOCNutation()", 0, false},
                              {testSphereRot_TEOtoCEO, 674, "psSphereRot_TEOtoCEO()", 0, false},
                              {testSphereRot_CEOtoGCRS, 675, "psSphereRot_CEOtoGCRS()", 0, false},
                              {testSphereRot_ITRStoTEO, 676, "psSphereRRot_ITRStoTEO()", 0, false},
                              {testSphereRotPrecess, 677, "psSphereRotPrecess()", 0, false},
                              {NULL}
                          };

psS32 main( psS32 argc, char* argv[] )
{
    psLogSetLevel( PS_LOG_INFO );

    // Initialize library internal structures
    psLibInit("pslib.config");

    if( !runTestSuite(stderr,"psEarthOrientation",tests,argc,argv)) {
        return 1;
    }

    // Cleanup library
    psLibFinalize();

    return 0;

    //    return ( ! runTestSuite( stderr, "psEarthOrientation", tests, argc, argv ) );
}

#define timesec 1049160600
#define objR DEG_TO_RAD(122.9153182445501)
#define objD DEG_TO_RAD(48.562968978679194)
#define VERBOSE 0
static psSphere *obj = NULL;
static void objSetup(void);
static double greatCircle(psSphere *in1, psSphere *in2);

static void objSetup(void)
{
    if (obj == NULL) {
        //        obj = psSphereAlloc();
        //        obj->r = objR;
        //        obj->d = objD;
        psCube *tempCube = psCubeAlloc();
        tempCube->x = -0.3598480726985338;
        tempCube->y = 0.5555012823608123;
        tempCube->z = 0.7496183628158023;
        obj = psCubeToSphere(tempCube);
        psFree(tempCube);
    }
}

//returns the great circle ANGULAR distance between 2 positions (psSphere's)
static double greatCircle(psSphere *in1, psSphere *in2)
{
    double r1, r2, d1, d2;
    d1 = in1->r;
    d2 = in2->r;
    r1 = in1->d;
    r2 = in2->d;
    double out = 0.0;
    double c1, c2, cd, s1, s2, sum, ac;
    c1 = cos(r1);
    c2 = cos(r2);
    cd = cos(d1-d2);
    s1 = sin(r1);
    s2 = sin(r2);
    sum = c1*c2*cd + s1*s2;
    if (sum > 1.0)
        sum = 1.0 - (sum - 1.0);

    ac = acos(sum);
    out = ac;
    //    printf("\n  c1=%lf, c2=%lf, cd=%lf, s1=%lf, s2=%lf, sum=%.19g, ac=%g\n", c1,c2,cd,s1,s2,sum,ac);
    //    out = acos(cos(r1)*cos(r2)*cos(d1-d2) + sin(r1)*sin(r2));
    //    if (out != 0.0) printf("\n   in greatCircle, out = %lf\n", out);
    return out;
}

psS32 testAberration(void)
{
    psSphere *apparent = NULL;
    psSphere *empty = NULL;
    psCube *actualCube = psCubeAlloc();
    //values from After gravity deflection//
    actualCube->x = -0.35961949760293604;
    actualCube->y = 0.5555613950298085;
    actualCube->z = 0.7496835020836093;
    psSphere *actual = psCubeToSphere(actualCube);
    psCube *cubeDir = psCubeAlloc();
    cubeDir->x = 5148.713262821658;
    cubeDir->y = -26945.04752348012;
    cubeDir->z = -11682.787302030947;
    cubeDir->x += -357.6031690489248;
    cubeDir->y += 248.46429758174693;
    cubeDir->z += 0.09694774143797581;
    psSphere *direction = psCubeToSphere(cubeDir);
    double speed = sqrt(cubeDir->x*cubeDir->x + cubeDir->y*cubeDir->y + cubeDir->z*cubeDir->z);
    double c = 299792458.0; // Speed of light in vacuum (src:NIST)   /* m/s */
    speed /= c;

    psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message");
    empty = psAberration(empty, apparent, direction, speed);
    if (empty != NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, false,
                "psAberration failed to return NULL for NULL actual input.\n");
        return 1;
    }
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message");
    empty = psAberration(empty, actual, apparent, speed);
    if (empty != NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, false,
                "psAberration failed to return NULL for NULL direction input.\n");
        return 2;
    }

    apparent = psAberration(apparent, actual, direction, speed);
    psCube *outCube = psSphereToCube(apparent);
    if (VERBOSE) {
        printf("\n -- resultCube = x,y,z =     %.13g,      %.13g,    %.13g  -- \n",
               outCube->x, outCube->y, outCube->z);
    }
    //expected cube values
    double x, y, z;
    x = -0.35963388069046304;
    y = 0.5555192509816625;
    z = 0.7497078321908413;

    if (VERBOSE) {
        printf(" -- expectedCube = x,y,z =   %.13g,     %.13g,    %.13g  -- \n", x, y, z);
        printf("Cube Difference  =  x,y,z  = %.13g, %.13g, %.13g \n\n",
               (x - outCube->x), (y - outCube->y), (z - outCube->z) );
    }
    psFree(actual);
    actualCube->x = x;
    actualCube->y = y;
    actualCube->z = z;
    actual = psCubeToSphere(actualCube);
    double xxx = greatCircle(actual, apparent);
    if (VERBOSE) {
        printf("   The great circle angular distance between expected & result = %.13g\n",
               xxx);
    }
    if ( fabs(x - outCube->x) > FLT_EPSILON || fabs(y - outCube->y) > FLT_EPSILON ||
            fabs(z - outCube->z) > FLT_EPSILON ) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
                "psAberration returned incorrect values.\n");
        //        printf("\nMagnitude of expected change = x,y,z = %.13g, %.13g, %.13g \n",
        //               (actualCube->x - x), (actualCube->y - y), (actualCube->z - z) );
        //        printf("Magnitude of actual change = x,y,z = %.13g, %.13g, %.13g \n",
        //               (actualCube->x - outCube->x), (actualCube->y - outCube->y),
        //               (actualCube->z - outCube->z) );
        psFree(actual);
        actualCube->x = x;
        actualCube->y = y;
        actualCube->z = z;
        actual = psCubeToSphere(actualCube);
        x = greatCircle(actual, apparent);
        if (VERBOSE) {
            printf("   The great circle angular distance between expected & result = %.13g\n",
                   x);
        }
        return 3;
    }

    psFree(outCube);
    psFree(actualCube);
    psFree(cubeDir);
    psFree(apparent);
    psFree(actual);
    psFree(direction);

    return 0;
}

psS32 testGravityDeflect(void)
{
    psSphere *apparent = NULL;
    psSphere *empty = NULL;
    psCube *actualCube = psCubeAlloc();
    actualCube->x = -0.3596195125758298;
    actualCube->y = 0.5555613903455866;
    actualCube->z = 0.7496834983724809;
    psSphere *actual = psCubeToSphere(actualCube);
    psCube *sunCube = psCubeAlloc();
    sunCube->x = 1.467797790127511e11;
    sunCube->y = 2.5880956908748722e10;
    sunCube->z = 1.1220046291457653e10;
    double sunLength = sqrt(sunCube->x*sunCube->x + sunCube->y*sunCube->y + sunCube->z*sunCube->z);
    sunCube->x /= sunLength;
    sunCube->y /= sunLength;
    sunCube->z /= sunLength;
    if (VERBOSE) {
        printf("sunCube = x,y,z = %.13g, %.13g, %.13g\n",
               sunCube->x, sunCube->y, sunCube->z);
    }
    psSphere *sun = psCubeToSphere(sunCube);
    if (VERBOSE) {
        printf("sunSphere  = r, d = %.13g, %.13g\n", sun->r, sun->d);
    }
    psCube *outCube = psCubeAlloc();

    empty = psGravityDeflection(apparent, empty, sun);
    if (empty != NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                "psGravityDeflection Failed to return NULL for NULL actual input sphere.\n");
        return 1;
    }
    empty = psGravityDeflection(apparent, actual, empty);
    if (empty != NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                "psGravityDeflection Failed to return NULL for NULL sun input sphere.\n");
        return 2;
    }

    apparent = psGravityDeflection(NULL, actual, sun);
    //    apparent->r *= -1.0;
    //    apparent->d *= -1.0;
    psSphere *result2;
    //    psSphere *result2 = psSphereSetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN);
    //    psSphere *result = psSphereSetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN);
    //    psSphere *result = psSphereGetOffset(apparent, actual, PS_SPHERICAL, PS_RADIAN);
    if (VERBOSE) {
        printf(" -- actualCube = x,y,z = %.13g, %.13g, %.13g  -- \n",
               actualCube->x, actualCube->y, actualCube->z);
    }
    //    psCube *outCube = psSphereToCube(result);
    //    printf(" -- resultCube = x,y,z = %.13g, %.13g, %.13g  -- \n",
    //           outCube->x, outCube->y, outCube->z);
    psCube *outCube2 = psSphereToCube(apparent);
    if (VERBOSE) {
        printf(" -- resultCube2= x,y,z = %.13g, %.13g, %.13g  -- \n",
               outCube2->x, outCube2->y, outCube2->z);
    }
    double x, y, z;
    x = -0.35961949760293604;
    y = 0.5555613950298085;
    z = 0.7496835020836093;

    if (VERBOSE) {
        printf(" -- expectCube = x,y,z = %.13g, %.13g, %.13g  -- \n\n", x, y, z);

        //    if ( fabs(x - outCube->x) > DBL_EPSILON || fabs(y - outCube->y) > DBL_EPSILON ||
        //            fabs(z - outCube->z) > DBL_EPSILON ) {
        //        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
        //                "psGravityDeflection returned incorrect values.\n");
        printf("expect-actual=  x,y,z  = %.13g, %.13g, %.13g \n",
               (x - actualCube->x), (y - actualCube->y), (z - actualCube->z) );
        printf("expect-result=  x,y,z  = %.13g, %.13g, %.13g \n",
               (x - outCube2->x), (y - outCube2->y), (z - outCube2->z) );
        printf("result-actual=  x,y,z  = %.13g, %.13g, %.13g \n",
               (outCube2->x - actualCube->x), (outCube2->y - actualCube->y),
               (outCube2->z - actualCube->z) );
    }
    //        return 1;
    //    }
    //    psFree(result2);
    outCube2->x = x;
    outCube2->y = y;
    outCube2->z = z;
    result2 = psCubeToSphere(outCube2);
    //    psFree(result);
    psSphere *result = psSphereGetOffset(actual, result2, PS_SPHERICAL, PS_RADIAN);
    psFree(result2);
    result2 = psSphereGetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN);
    if (VERBOSE) {
        printf("The apparent output sphere = r,d = %.13g, %.13g\n", result2->r, result2->d);
        printf("The expected output sphere = r,d = %.13g, %.13g\n\n", result->r, result->d);
    }
    psFree(result2);

    psFree(outCube);
    psFree(outCube2);
    psFree(sunCube);
    psFree(actualCube);
    psFree(result);
    psFree(actual);
    psFree(apparent);
    psFree(sun);

    return 0;
}

psS32 testEOCPrecession(void)
{
    psTime *empty = NULL;
    psTime *time = psTimeAlloc(PS_TIME_UTC);
    time->sec = timesec;
    time->nsec = 0;
    time->leapsecond = false;

    //Tests for Precession Model //
    psEarthPole *pmodel = NULL;
    //Return NULL for NULL time input
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message");
    pmodel = psEOC_PrecessionModel(empty);
    if (pmodel != NULL) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
                "psEOC_PrecessionModel failed to return NULL for NULL time input.\n");
        return 1;
    }
    //Return NULL for UT1 time input
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message");
    psTime *UT1time = psTimeAlloc(PS_TIME_UT1);
    pmodel = psEOC_PrecessionModel(UT1time);
    if (pmodel != NULL) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
                "psEOC_PrecessionModel failed to return NULL for UT1 input time.\n");
        return 2;
    }
    psFree(UT1time);
    //Check return values from valid precession input
    pmodel = psEOC_PrecessionModel(time);
    if ( pmodel == NULL ) {
        psError(PS_ERR_BAD_PARAMETER_NULL, false,
                "psEOC_PrecessionModel returned NULL for valid input.\n");
        return 3;
    } else {
        double x, y, s;
        x = 2.857175590089105e-4;
        y = 2.3968739377734732e-5;
        s = -1.3970066457904322e-8;
        if ( fabs(pmodel->x - x) > FLT_EPSILON || fabs(pmodel->y - y) > FLT_EPSILON
                || fabs(pmodel->s - s) > FLT_EPSILON) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, false,
                    "   psEOC_PrecessionModel return incorrect values.\n");
            return 4;
        }
        if (VERBOSE) {
            printf("\n  PrecessionModel output = x,y,s = %.13g,  %.13g,  %.13g\n",
                   pmodel->x, pmodel->y, pmodel->s);
            printf("  Expected output        = x,y,s = %.13g, %.13g,  %.13g\n", x, y, s);
            printf("  A difference of:                 %.13g, %.13g, %.13g\n",
                   (pmodel->x - x), (pmodel->y - y), (pmodel->s - s) );
        }
    }
    psFree(pmodel);

    psFree(time);
    if (!p_psEOCFinalize() ) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "EOC failed finalization!\n");
        return 12;
    }

    return 0;
}

psS32 testEOCPrecessionCorr(void)
{
    psTime *empty = NULL;
    psTime *time2 = psTimeAlloc(PS_TIME_UTC);
    time2->sec = timesec;
    time2->nsec = 0;
    time2->leapsecond = false;
    //    time2 = psTimeConvert(time2, PS_TIME_TAI);

    //Tests for Precession Correction function//
    //Return NULL for NULL time input
    psEarthPole *pcorr = NULL;
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message");
    pcorr = psEOC_PrecessionCorr(empty, PS_IERS_A);
    if (pcorr != NULL) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
                "psEOC_PrecessionCorr failed to return NULL for NULL time input.\n");
        return 5;
    }

    //Return NULL for Invalid IERS table
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message");
    pcorr = psEOC_PrecessionCorr(time2, 3);
    if (pcorr != NULL) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
                "psEOC_PrecessionCorr failed to return NULL for incorrect IERS table.\n");
        return 6;
    }
    psFree(pcorr);

    //Check values from IERS table A
    pcorr = psEOC_PrecessionCorr(time2, PS_IERS_A);
    if ( pcorr == NULL ) {
        psError(PS_ERR_BAD_PARAMETER_NULL, false,
                "psEOC_PrecessionCorr returned NULL for valid inputs.\n");
        return 7;
    } else {
        if (VERBOSE) {
            printf("\nPrecessionCorr output (IERSA) = x,y,s = %.13g,  %.13g, %.13g\n",
                   pcorr->x, pcorr->y, pcorr->s);
        }
    }
    psFree(pcorr);

    //Check values from IERS table B
    pcorr = psEOC_PrecessionCorr(time2, PS_IERS_B);
    if ( pcorr == NULL ) {
        psError(PS_ERR_BAD_PARAMETER_NULL, false,
                "psEOC_PrecessionCorr returned NULL for valid inputs.\n");
        return 9;
    } else {
        double xx, yy, ss;
        xx = 0.06295703125;
        yy = -0.0287618408203125;
        ss = 0.0;
        xx = SEC_TO_RAD(xx) * 1e-3;
        yy = SEC_TO_RAD(yy) * 1e-3;
        //        if ( fabs(pcorr->x - xx) > DBL_EPSILON || fabs(pcorr->y - yy) > DBL_EPSILON
        //                || fabs(pcorr->s - ss) > DBL_EPSILON) {
        //            psError(PS_ERR_BAD_PARAMETER_VALUE, false,
        //                    "   psEOC_PrecessionCorr return incorrect values.\n");
        if (VERBOSE) {
            printf("PrecessionCorr output (IERSB) = x,y,s = %.13g, %.13g, %.13g\n",
                   pcorr->x, pcorr->y, pcorr->s);
            printf("Expected output               = x,y,s = %.13g, %.13g, %.13g\n", xx, yy, ss);
            printf("          A difference of:              %.13g,  %.13g, %.13g\n\n",
                   (pcorr->x - xx), (pcorr->y - yy), (pcorr->s - ss) );
        }
        //            return 10;
        //        }
    }
    //precess is the *actual* output from PrecessionModel + PrecessionCorr
    psEarthPole *precess = psEOC_PrecessionModel(time2);
    precess->x += pcorr->x;
    precess->y += pcorr->y;
    double xCorr, yCorr;
    xCorr = 3.05224300720406e-10;
    yCorr = -1.39441339235822e-10;
    //pcorr is the *expected* output from PrecessionModel//
    pcorr->x = 2.857175590089105e-4;
    pcorr->y = 2.3968739377734732e-5;
    pcorr->s = -1.3970066457904322e-8;
    pcorr->x += xCorr;
    pcorr->y += yCorr;
    psSphereRot *precessNutInv = psSphereRot_CEOtoGCRS(precess);
    psSphereRot *precessNut = psSphereRotConjugate(NULL, precessNutInv);
    double q0, q1, q2;
    q0 = -1.1984522406756289e-5;
    q1 = 1.4285893358610674e-4;
    q2 = 1.2191193518914336e-10;
    psSphereRot *pni = psSphereRot_CEOtoGCRS(pcorr);
    if (fabs(pni->q0-q0) > FLT_EPSILON || fabs(pni->q1-q1) > FLT_EPSILON ||
            fabs(pni->q2-q2) > FLT_EPSILON ) {
        printf("\n Error at CEOtoGCRS, output psSphereRot doesn't match expected.\n");
    }
    //    printf("  Output from CEOtoGCRS only  = %.13g,%.13g,%.13g,%.13g\n",
    //           pni->q0, pni->q1, pni->q2, pni->q3);
    if (VERBOSE) {
        printf("  Expected sphere rotation    = %.13g, %.13g, %.13g\n", q0,q1,q2);
        printf("  Result sphere rotation      = %.13g, %.13g, %.13g\n",
               precessNutInv->q0, precessNutInv->q1, precessNutInv->q2);
        printf("     Difference         =        %.13g, %.13g, %.13g\n\n",
               precessNutInv->q0-q0, precessNutInv->q1-q1, precessNutInv->q2-q2);
    }
    psCube *objC = psCubeAlloc();
    //    objC->x = -3.5963388069046304;
    //    objC->y = 0.5555192509816625;
    //    objC->z = 0.7497078321908413;
    //    objSetup();

    //This is the sphere rotation for the *expected* precession output//
    psSphereRot *pn = psSphereRotConjugate(NULL, pni);

    //    psSphere *sphere = psSphereAlloc();
    //    *sphere = *obj;
    //    psFree(obj);

    //create a psSphere for (from) the start position given in eoc_testing//
    objC->x = -0.35963388069046304;
    objC->y = 0.5555192509816625;
    objC->z = 0.7497078321908413;
    psSphere *sphere = psCubeToSphere(objC);

    psSphere *expect = psSphereRotApply(NULL, pn, sphere);
    //expected results below - stored in:  sphere  //
    double x,y,z;
    x = -0.3598480726985338;
    y = 0.5555012823608123;
    z = 0.7496183628158023;
    if (VERBOSE) {
        printf("\n<<Expected out       = x,y,z = %.13g, %.13g, %.13g\n", x, y, z);
    }
    //    psFree(objC);
    //    objC = psSphereToCube(expect);
    //printf("<<Expected out (CEO)  = x,y,z = %.13g, %.13g, %.13g\n", objC->x, objC->y, objC->z);
    //printf("     Difference     =           %.13g, %.13g, %.13g\n", objC->x-x, objC->y-y, objC->z-z);
    //    x = objC->x;
    //    y = objC->y;
    //    z = objC->z;
    psSphere *result = psSphereRotApply(NULL, precessNut, sphere);
    psFree(objC);
    objC = psSphereToCube(result);
    double xx = greatCircle(result, expect);
    if (VERBOSE) {
        printf("<<Resulting out      = x,y,z = %.13g, %.13g, %.13g\n", objC->x, objC->y, objC->z);
        printf("     Difference         =      %.13g, %.13g, %.13g\n\n",
               objC->x-x, objC->y-y, objC->z-z);
        printf("GREAT CIRCLE DIFFERENCE = %.13g \n", xx);
    }

    psFree(precess);
    psFree(precessNut);
    psFree(precessNutInv);
    psFree(expect);
    psFree(objC);

    psFree(sphere);
    psFree(result);
    psFree(pn);
    psFree(pni);
    psFree(pcorr);
    psFree(time2);
    if (!p_psEOCFinalize() ) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "EOC failed finalization!\n");
        return 12;
    }

    return 0;
}

psS32 testEOCPolar(void)
{
    psTime *in = psTimeAlloc(PS_TIME_UTC);
    in->sec = timesec;
    in->nsec = 0;
    in->leapsecond = false;
    psTime *empty = NULL;
    psEarthPole *polarMotion = NULL;

    //Return NULL for NULL input time
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message");
    polarMotion = psEOC_GetPolarMotion(empty, PS_IERS_B);
    if (polarMotion != NULL) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
                "psEOC_GetPolarMotion failed to return NULL for NULL input time.\n");
        return 1;
    }
    //Return NULL for incorrect Bulletin.
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message");
    polarMotion = psEOC_GetPolarMotion(empty, 3);
    if (polarMotion != NULL) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
                "psEOC_GetPolarMotion failed to return NULL for NULL input time.\n");
        return 1;
    }

    //Test for IERS bulletin A.
    double x, y, s;
    x = -6.454389659777e-07;
    y = 2.112606414597e-06;
    s = 0.0;
    polarMotion = psEOC_GetPolarMotion(in, PS_IERS_A);
    if (polarMotion == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, false,
                "psEOC_GetPolarMotion returned NULL for valid input time.\n");
        return 4;
    }
    //    if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON
    //            || fabs(polarMotion->s - s) > DBL_EPSILON) {
    //        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    //                "psEOC_GetPolarMotion returned incorrect values.\n");
    if (VERBOSE) {
        printf("  <>PolarMotion output (IERSA)   = x,y,s = %.13g, %.13g, %.13g\n",
               polarMotion->x, polarMotion->y, polarMotion->s);
    }
    //        printf("  <>PolarMotion expected (IERSA) = x,y,s = %.13g, %.13g, %.13g\n",
    //               x, y, s);
    //        return 5;
    //    }
    psFree(polarMotion);

    //Return valid values for correct input time.  Test IERS Bulletin B.
    polarMotion = psEOC_GetPolarMotion(in, PS_IERS_B);
    //    x = -6.45381397904e-07;
    //    y = 2.112819726698e-06;
    //    s = 0.0;
    if (polarMotion == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, false,
                "psEOC_GetPolarMotion returned NULL for valid input time.\n");
        return 2;
    }
    //    if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON
    //            || fabs(polarMotion->s - s) > DBL_EPSILON) {
    //        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    //                "psEOC_GetPolarMotion returned incorrect values.\n");
    if (VERBOSE) {
        printf("  <>PolarMotion output (IERSB)   = x,y,s = %.13g,  %.13g, %.13g\n",
               polarMotion->x, polarMotion->y, polarMotion->s);
    }
    //        printf("  <>PolarMotion expected (IERSB) = x,y,s = %.13g, %.13g, %.13g\n",
    //               x, y, s);
    //        return 3;
    //    }
    psEarthPole *nutationCorr = psEOC_NutationCorr(in);
    if (nutationCorr == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, false, "Nutation Correction returned NULL.\n");
        return 6;
    }
    polarMotion->x += nutationCorr->x;
    polarMotion->y += nutationCorr->y;
    polarMotion->s += nutationCorr->s;
    psEarthPole *polarTide = psEOC_PolarTideCorr(in);
    polarMotion->x += polarTide->x;
    polarMotion->y += polarTide->y;
    psFree(polarTide);
    x = -6.43607313124045e-7;
    y = 2.11351436973568e-6;
    s = -7.39617581324646e-12;
    //    if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON
    //            || fabs(polarMotion->s - s) > DBL_EPSILON) {
    //        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
    //                "   psEOC_GetPolarMotion returned incorrect values.\n");
    if (VERBOSE) {
        printf("\n  PolarMotion + NutationCorr out = x,y,s = %.13g, %.13g, %.13g\n",
               polarMotion->x, polarMotion->y, polarMotion->s);
        printf("  Expected output                = x,y,s = %.13g,  %.13g, %.13g\n", x, y, s);
        printf("                     Difference  = x,y,s = %.13g, %.13g,  %.13g\n",
               polarMotion->x - x, polarMotion->y - y, polarMotion->s - s);
    }
    //    }

    if (!p_psEOCFinalize() ) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "EOC failed finalization!\n");
        return 12;
    }

    psFree(nutationCorr);
    psFree(in);
    psFree(polarMotion);
    return 0;
}

psS32 testEOCPolarTide(void)
{
    psTime *in = psTimeAlloc(PS_TIME_UTC);
    in->sec = timesec;
    in->nsec = 0;
    in->leapsecond = false;
    psTime *empty = NULL;
    psEarthPole *eop = NULL;

    //Return NULL for NULL input time
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message");
    eop = psEOC_PolarTideCorr(empty);
    if (eop != NULL) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
                "psEOC_PolarTideCorr failed to return NULL for NULL input time.\n");
        return 1;
    }

    eop = psEOC_PolarTideCorr(in);
    if (eop == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, false,
                "psEOC_PolarTideCorr returned NULL for valid input time.\n");
        return 2;
    } else {
        if (VERBOSE) {
            printf("\nPolarTideCorr output = x,y,s = %.13g, %.13g, %.13g\n",
                   eop->x, eop->y, eop->s);
        }
    }

    psFree(in);
    psFree(eop);

    return 0;
}

psS32 testEOCNutation(void)
{
    psTime *in = psTimeAlloc(PS_TIME_UTC);
    in->sec = timesec;
    in->nsec = 0;
    in->leapsecond = false;
    psTime *empty = NULL;
    psEarthPole *nute = NULL;

    //Return NULL for NULL input time.
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message");
    nute = psEOC_NutationCorr(empty);
    if (nute != NULL) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
                "psEOC_NutationCorr failed to return NULL for NULL input time.\n");
        return 1;
    }
    //Return NULL for UT1 time input
    /*    psTime *UT1time = psTimeAlloc(PS_TIME_UT1);
        nute = psEOC_NutationCorr(UT1time);
        if (nute != NULL) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, false,
                    "psEOC_NutationCorr failed to return NULL for UT1 input time.\n");
            return 2;
        }
        psFree(UT1time);
    */
    //Check return values from valid nutation time input
    nute = psEOC_NutationCorr(in);
    if ( nute == NULL ) {
        psError(PS_ERR_BAD_PARAMETER_NULL, false,
                "psEOC_NutationCorr returned NULL for valid input.\n");
        return 3;
    } else {
        if (VERBOSE) {
            printf("Nutation Correction output = x,y,s = %.13g, %.13g, %.13g\n\n",
                   nute->x, nute->y, nute->s);
        }
    }
    psFree(nute);
    psFree(in);

    if (!p_psEOCFinalize() ) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "EOC failed finalization!\n");
        return 12;
    }

    return 0;
}

psS32 testSphereRot_TEOtoCEO(void)
{
    psSphereRot *rot = NULL;
    psTime *empty = NULL;
    psTime *time = psTimeAlloc(PS_TIME_UT1);
    time->sec = timesec-1;
    time->nsec = 657017200;
    time->leapsecond = false;

    //return NULL for NULL input time
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message");
    rot = psSphereRot_TEOtoCEO(empty, NULL);
    if (rot != NULL) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
                "psSphereRot_TEOtoCEO failed to return NULL for NULL input time.\n");
        return 1;
    }

    psEarthPole *polarTideCorr = psEOC_PolarTideCorr(time);
    psSphereRot *teoceo = psSphereRot_TEOtoCEO(time, polarTideCorr);
    //Make sure values match for other psTime type
    empty = psTimeAlloc(PS_TIME_UTC);
    empty->sec = timesec;
    empty->nsec = 0;
    empty->leapsecond = false;

    rot = psSphereRot_TEOtoCEO(empty, polarTideCorr);
    if (fabs(rot->q0-teoceo->q0) > 0.00001 || fabs(rot->q1-teoceo->q1) > 0.00001 ||
            fabs(rot->q2-teoceo->q2) > 0.00001 || fabs(rot->q3-teoceo->q3) > 0.00001) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
                "psSphereRot_TEOtoCEO failed to return matching values for different time types.\n");
        if (VERBOSE) {
            printf("\n  Output Rotation1 = q0,q1,q2,q3 = %.13g, %.13g, %.13g, %.13g\n",
                   teoceo->q0, teoceo->q1, teoceo->q2, teoceo->q3 );
            printf("\n  Output Rotation2 = q0,q1,q2,q3 = %.13g, %.13g, %.13g, %.13g\n",
                   rot->q0, rot->q1, rot->q2, rot->q3 );
        }
        return 2;
    }
    if (VERBOSE) {
        printf("\n  Output Rotation = q0,q1,q2,q3 = %.13g, %.13g, %.13g, %.13g\n",
               teoceo->q0, teoceo->q1, teoceo->q2, teoceo->q3 );
    }

    objSetup();
    psSphereRot *earthRot = psSphereRotConjugate(NULL, teoceo);
    psSphere *result = psSphereRotApply(NULL, earthRot, obj);
    //    psSphere *result = psSphereRotApply(NULL, teoceo, obj);
    psCube *cube = psSphereToCube(result);

    double x, y, z;
    x = 0.01698625430807123;
    y = -0.6616523084626379;
    z = 0.7496183628158023;
    if ( fabs(x-cube->x) > FLT_EPSILON  || fabs(y-cube->y) > FLT_EPSILON ||
            fabs(z-cube->z) > FLT_EPSILON) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
                "psSphereRot_TEOtoCEO returned incorrect values.\n");
        if (VERBOSE) {
            printf("\nOutput cube = x,y,z = %.13g, %.13g, %.13g\n",
                   cube->x, cube->y, cube->z);
            printf("Expected cube = x,y,z = %.13g, %.13g, %.13g\n", x, y, z);
            printf("A difference of:   %.13g, %.13g, %.13g\n\n",
                   (x-cube->x), (y-cube->y), (z-cube->z));
        }
        return 3;
    }

    psFree(polarTideCorr);
    psFree(rot);
    psFree(empty);
    psFree(result);
    psFree(obj);
    psFree(earthRot);
    psFree(cube);
    psFree(time);
    psFree(teoceo);
    return 0;
}

psS32 testSphereRot_CEOtoGCRS(void)
{
    psEarthPole *in = psEarthPoleAlloc();
    psEarthPole *empty = NULL;
    psSphereRot *rot = NULL;
    in->x = 2.857175590089105e-4;
    in->y = 2.3968739377734732e-5;
    in->s = -1.3970066457904322e-8;

    double q0,q1,q2,q3;
    q0 = -1.1984522406756289e-5;
    q1 = 1.4285893358610674e-4;
    q2 = 1.2191193518914336e-10;
    q3 = -0.9999999897238481;

    //Return NULL for NULL input earthpole
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message");
    rot = psSphereRot_CEOtoGCRS(empty);
    if (rot != NULL) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
                "psSphereRot_CEOtoGCRS failed to return NULL for NULL input psEarthPole.\n");
        return 1;
    }

    rot = psSphereRot_CEOtoGCRS(in);
    if (rot == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, false,
                "psSphereRot_CEOtoGCRS returned NULL for valid psEarthPole input.\n");
        return 2;
    }

    if (VERBOSE) {
        printf("\n  Output sphere rotation   = %.13g, %.13g, %.13g, %.13g\n",
               rot->q0, rot->q1, rot->q2, rot->q3);
        printf("  Expected sphere rotation = %.13g, %.13g, %.13g, %.13g\n", q0,q1,q2,q3);
        printf("  difference:                 %.13g, %.13g, %.13g \n",
               (rot->q0-q0), (rot->q1-q1), (rot->q2-q2) );
    }
    if (fabs(rot->q0-q0) > FLT_EPSILON || fabs(rot->q1-q1) > FLT_EPSILON ||
            fabs(rot->q2-q2) > FLT_EPSILON || fabs(rot->q3+q3) > FLT_EPSILON) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
                "psSphereRot_CEOtoGCRS failed to return expected values.\n");
        return 3;
    }
    psCube *tempCube = psCubeAlloc();
    tempCube->x = -0.35963388069046304;
    tempCube->y = 0.5555192509816625;
    tempCube->z = 0.7497078321908413;
    obj = psCubeToSphere(tempCube);
    psFree(tempCube);
    psSphereRot *precessionNutation = psSphereRotConjugate(NULL, rot);
    psSphere *result = psSphereRotApply(NULL, precessionNutation, obj);
    psCube *cube = psSphereToCube(result);
    double x, y, z;
    x = -0.3598480726985338;
    y = 0.5555012823608123;
    z = 0.7496183628158023;
    if (VERBOSE) {
        printf("\n  Output cube = x,y,z = %.13g,  %.13g,    %.13g\n", cube->x, cube->y, cube->z);
        printf("Expected cube = x,y,z = %.13g,  %.13g,   %.13g\n", x, y, z);
        printf("  A difference of:    %.13g, %.13g, %.13g\n\n",
               (x-cube->x), (y-cube->y), (z-cube->z));
    }
    psCube *expect = psCubeAlloc();
    expect->x = x;
    expect->y = y;
    expect->z = z;
    psSphere *expected = psCubeToSphere(expect);
    double d = greatCircle(result, expected);
    if (VERBOSE) {
        printf("   The great circle angular distance between expected & result = %.13g\n", d);
    }

    psFree(expect);
    psFree(expected);
    psFree(obj);
    psFree(precessionNutation);
    psFree(result);
    psFree(cube);
    psFree(rot);
    psFree(in);

    return 0;
}

psS32 testSphereRot_ITRStoTEO(void)
{
    psEarthPole *in = psEarthPoleAlloc();
    psEarthPole *empty = NULL;
    psSphereRot *rot = NULL;
    in->x = -0.13275353774074533;
    in->y = 0.4359436319739848;
    in->s = -4.2376965863576153e-10;
    in->x = SEC_TO_RAD(in->x);
    in->y = SEC_TO_RAD(in->y);
    in->s = SEC_TO_RAD(in->s);

    //Return NULL for NULL input earthpole
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message");
    rot = psSphereRot_ITRStoTEO(empty);
    if (rot != NULL) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
                "psSphereRot_ITRStoTEO failed to return NULL for NULL input psEarthPole.\n");
        return 1;
    }

    rot = psSphereRot_ITRStoTEO(in);
    if (rot == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, false,
                "psSphereRot_ITRStoTEO return NULL for valid psEarthPole input.\n");
        return 2;
    }

    double q0,q1,q2,q3;
    q0 = -1.0567571848664005e-6;
    q1 = 3.218036562931509e-7;
    q2 = -3.3580195807204483e-12;
    q3 = -0.9999999999993899;
    if (fabs(rot->q0-q0) > FLT_EPSILON || fabs(rot->q1-q1) > FLT_EPSILON ||
            fabs(rot->q2-q2) > FLT_EPSILON
            //            || fabs(rot->q3-q3) > DBL_EPSILON
       ) {
        if (VERBOSE) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, false,
                    "psSphereRot_ITRStoTEO failed to return expected values.\n");
            printf("\n  Output sphere rotation   = %.13g, %.13g, %.13g, %.13g\n",
                   rot->q0, rot->q1, rot->q2, rot->q3);
            printf("  Expected sphere rotation = %.13g, %.13g, %.13g, %.13g\n", q0,q1,q2,q3);
            printf("  a difference:   %.13g, %.13g, %.13g, %.13g \n", (rot->q0-q0), (rot->q1-q1),
                   (rot->q2-q2), (rot->q3-q3) );
        }
        //        printf("Error #3\n");
        //        return 3;
    }
    if (VERBOSE) {
        printf("\n  Output sphere rotation   = %.13g, %.13g, %.13g, %.13g\n",
               rot->q0, rot->q1, rot->q2, rot->q3);
        printf("  Expected sphere rotation = %.13g, %.13g, %.13g, %.13g\n", q0,q1,q2,q3);
    }

    psCube *temp = psCubeAlloc();
    //    temp->x = -0.3596195125758298;
    //    temp->y = 0.5555613903455866;
    //    temp->z = 0.7496834983724809;
    temp->x = 0.01698625430807123;
    temp->y = -0.6616523084626379;
    temp->z = 0.7496183628158023;

    obj = psCubeToSphere(temp);
    psSphere *test = NULL;
    psFree(temp);
    psSphereRot *newRot = psSphereRotConjugate(NULL, rot);
    test = psSphereRotApply(NULL, newRot, obj);
    temp = psSphereToCube(test);
    double x, y, z;
    x = 0.01698577185310146;
    y = -0.6616538927902393;
    z = 0.7496169753347885;
    if (VERBOSE) {
        printf("\n  Cube -test- has x,y,z =     %.13g, %.13g, %.13g \n",
               temp->x, temp->y, temp->z);
        printf("\n  Cube -expected- has x,y,z = %.13g,  %.13g, %.13g \n", x, y, z );
    }
    temp->x = x;
    temp->y = y;
    temp->z = z;
    psSphere *sphere = psCubeToSphere(temp);
    double d = greatCircle(sphere, test);
    if (VERBOSE) {
        printf("Great circle difference of:  %.13g \n", d);
    }

    psFree(sphere);
    psFree(newRot);
    psFree(obj);
    psFree(temp);
    psFree(test);
    psFree(rot);
    psFree(in);

    return 0;
}

#define SPHERE_PRECESS_TP1_R            0.0               //    0.0       degrees
#define SPHERE_PRECESS_TP1_D            0.0               //    0.0       degrees
#define SPHERE_PRECESS_TP1_EXPECT_R     6.238453          //  357.437     degrees
#define SPHERE_PRECESS_TP1_EXPECT_D    -0.019426          //   -1.113     degrees
#define SPHERE_PRECESS_TP2_R            0.0               //    0.0       degrees
#define SPHERE_PRECESS_TP2_D            1.570796          //   90.0       degrees
#define SPHERE_PRECESS_TP2_EXPECT_R     6.260828          //  358.719     degrees
#define SPHERE_PRECESS_TP2_EXPECT_D     1.551353          //   88.886     degrees
#define SPHERE_PRECESS_TP3_R            3.141593          //  180.0       degrees
#define SPHERE_PRECESS_TP3_D            0.523599          //   30.0       degrees
#define SPHERE_PRECESS_TP3_EXPECT_R     3.096616          //  177.423     degrees
#define SPHERE_PRECESS_TP3_EXPECT_D     0.543024          //   31.113     degrees
#define ERROR_TOL   0.0001
#define MJD_1900  15021.0        // Modified Julian Day 1/1/1900 00:00:00
#define MJD_2100  88069.0        // Modified Julian Day 1/1/2100 00:00:00

psS32 testSphereRotPrecess( void )
{

    psSphere*     inputCoord  = psSphereAlloc();
    psSphere*     outputCoord = NULL;
    //    psTime*       fromTime    = psTimeFromMJD(MJD_1900);
    //    psTime*       toTime      = psTimeFromMJD(MJD_2100);
    psTime*       fromTime    = psTimeFromMJD(MJD_2100);
    psTime*       toTime      = psTimeFromMJD(MJD_1900);
    psSphereRot *rot = NULL;
    // Set input coordinate
    inputCoord->r = SPHERE_PRECESS_TP1_R;
    inputCoord->d = SPHERE_PRECESS_TP1_D;
    inputCoord->rErr = 0.0;
    inputCoord->dErr = 0.0;

    // Calculate precess
    rot = psSpherePrecess(fromTime, toTime, PS_PRECESS_ROUGH);
    outputCoord = psSphereRotApply(NULL, rot, inputCoord);
    if (outputCoord->r < -0.0001) {
        outputCoord->r += 2.0 * M_PI;
    }
    //    outputCoord = psSpherePrecess(inputCoord, fromTime, toTime, PS_PRECESS_ROUGH);
    // Verify return is not NULL
    if(outputCoord == NULL) {
        psError(PS_ERR_UNKNOWN,true,"Returned NULL not expected");
        return 1;
    }
    // Verify return with expected values
    if( fabs(outputCoord->r - SPHERE_PRECESS_TP1_EXPECT_R) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"Precess r = %lg not equal to expected = %lg",
                outputCoord->r,SPHERE_PRECESS_TP1_EXPECT_R);
        return 2;
    }
    if( fabs(outputCoord->d - SPHERE_PRECESS_TP1_EXPECT_D) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"Precess d = %lg not equal to expected = %lg",
                outputCoord->d,SPHERE_PRECESS_TP1_EXPECT_D);
        return 3;
    }
    psFree(outputCoord);
    psFree(rot);

    // Set input coordinate
    inputCoord->r = SPHERE_PRECESS_TP2_R;
    inputCoord->d = SPHERE_PRECESS_TP2_D;
    inputCoord->rErr = 0.0;
    inputCoord->dErr = 0.0;

    // Calculate precess
    rot = psSpherePrecess(fromTime, toTime, PS_PRECESS_ROUGH);
    outputCoord = psSphereRotApply(NULL, rot, inputCoord);
    if (outputCoord->r < -0.0001) {
        outputCoord->r += 2.0 * M_PI;
    }
    //    outputCoord = psSpherePrecess(inputCoord, fromTime, toTime, PS_PRECESS_ROUGH);
    // Verify return is not NULL
    if(outputCoord == NULL) {
        psError(PS_ERR_UNKNOWN,true,"Returned NULL not expected");
        return 4;
    }
    // Verify return with expected values
    if( fabs(outputCoord->r - SPHERE_PRECESS_TP2_EXPECT_R) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"Precess r = %lg not equal to expected = %lg",
                outputCoord->r,SPHERE_PRECESS_TP2_EXPECT_R);
        return 5;
    }
    if( fabs(outputCoord->d - SPHERE_PRECESS_TP2_EXPECT_D) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"Precess d = %lg not equal to expected = %lg",
                outputCoord->d,SPHERE_PRECESS_TP2_EXPECT_D);
        return 6;
    }
    psFree(outputCoord);
    psFree(rot);

    // Set input coordinate
    inputCoord->r = SPHERE_PRECESS_TP3_R;
    inputCoord->d = SPHERE_PRECESS_TP3_D;
    inputCoord->rErr = 0.0;
    inputCoord->dErr = 0.0;

    // Calculate precess
    rot = psSpherePrecess(fromTime, toTime, PS_PRECESS_ROUGH);
    outputCoord = psSphereRotApply(NULL, rot, inputCoord);
    if (outputCoord->r < -0.0001) {
        outputCoord->r += 2.0 * M_PI;
    }
    //    outputCoord = psSpherePrecess(inputCoord, fromTime, toTime, PS_PRECESS_ROUGH);
    // Verify return is not NULL
    if(outputCoord == NULL) {
        psError(PS_ERR_UNKNOWN,true,"Returned NULL not expected");
        return 7;
    }
    // Verify return with expected values
    if( fabs(outputCoord->r - SPHERE_PRECESS_TP3_EXPECT_R) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"Precess r = %lg not equal to expected = %lg",
                outputCoord->r,SPHERE_PRECESS_TP3_EXPECT_R);
        return 8;
    }
    if( fabs(outputCoord->d - SPHERE_PRECESS_TP3_EXPECT_D) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"Precess d = %lg not equal to expected = %lg",
                outputCoord->d,SPHERE_PRECESS_TP3_EXPECT_D);
        return 9;
    }
    psFree(outputCoord);
    psFree(rot);

    //Test other modes (IAU2000A, COMPLETE_[A,B])
    // Set input coordinate
    inputCoord->r = SPHERE_PRECESS_TP1_R;
    inputCoord->d = SPHERE_PRECESS_TP1_D;
    inputCoord->rErr = 0.0;
    inputCoord->dErr = 0.0;
    // Calculate precess
    rot = psSpherePrecess(fromTime, toTime, PS_PRECESS_COMPLETE_A);
    outputCoord = psSphereRotApply(NULL, rot, inputCoord);
    if (outputCoord->r < -0.000001) {
        outputCoord->r += 2.0 * M_PI;
    }
    // Verify return is not NULL
    if(outputCoord == NULL) {
        psError(PS_ERR_UNKNOWN,true,"Returned NULL not expected");
        return 10;
    }
    psFree(outputCoord);
    psFree(rot);

    rot = psSpherePrecess(fromTime, toTime, PS_PRECESS_COMPLETE_B);
    outputCoord = psSphereRotApply(NULL, rot, inputCoord);
    if (outputCoord->r < -0.000001) {
        outputCoord->r += 2.0 * M_PI;
    }
    // Verify return is not NULL
    if(outputCoord == NULL) {
        psError(PS_ERR_UNKNOWN,true,"Returned NULL not expected");
        return 11;
    }
    psFree(outputCoord);
    psFree(rot);

    rot = psSpherePrecess(fromTime, toTime, PS_PRECESS_IAU2000A);
    outputCoord = psSphereRotApply(NULL, rot, inputCoord);
    if (outputCoord->r < -0.000001) {
        outputCoord->r += 2.0 * M_PI;
    }
    // Verify return is not NULL
    if(outputCoord == NULL) {
        psError(PS_ERR_UNKNOWN,true,"Returned NULL not expected");
        return 12;
    }
    psFree(outputCoord);
    psFree(rot);

    // Invoke precess with invalid parameter
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message");
    rot = psSpherePrecess(fromTime, toTime, 10);
    if(rot != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Did not return NULL with invalid input");
        return 14;
    }
    // Return NULL for NULL input times
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message");
    rot = psSpherePrecess(NULL, NULL, PS_PRECESS_ROUGH);
    if(rot != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Did not return NULL with invalid input");
        return 15;
    }

    if (!p_psEOCFinalize() ) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "EOC failed finalization!\n");
        return 12;
    }
    // Free objects
    psFree(fromTime);
    psFree(toTime);
    psFree(inputCoord);

    return 0;
}
