/** @file  tst_psEarthOrientation.c
*
*  @brief The code will perform earth orientation calculations and sphere rotations.
*
*  @author d-Rob, MHPCC
*
*  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
*  @date $Date: 2005-12-20 05:05:37 $
*
*  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);

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},
                              {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)
static psSphere *obj = NULL;
static void objSetup(void);

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);
    }
}

psS32 testAberration(void)
{

    psSphere *apparent = NULL;
    psSphere *actual = psSphereAlloc();
    //    psSphere *direction = psSphereAlloc();
    psSphere *empty = NULL;

    //    actual->r = 0.2;
    //    actual->d = 0.2;
    //    actual->r = DEG_TO_RAD(45.0);
    //    actual->d = DEG_TO_RAD(30.0);
    actual->r = DEG_TO_RAD(122.9153182445501);
    actual->d = DEG_TO_RAD(48.562968978679194);
    //    direction->r = 0.2035;
    //    direction->d = 0.2035;
    //    direction->r = DEG_TO_RAD(48.0);
    //    direction->r = DEG_TO_RAD(-157.0);
    //    direction->d = DEG_TO_RAD(20.7072);
    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);
    // Speed of light in vacuum (src:NIST)
    double c = 299792458.0; /* m/s */
    speed = speed / c;

    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;
    }
    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);

    printf("\nSphere Difference  =  r,d  =  %.13g, %.13g\n",
           (actual->r - apparent->r), (actual->d - apparent->d));
    psCube *outCube = psSphereToCube(apparent);
    printf(" -- resultCube = x,y,z = %.13g, %.13g, %.13g  -- \n",
           outCube->x, outCube->y, outCube->z);

    double x, y, z;
    x = -0.35963388069046304;
    y = 0.5555192509816625;
    z = 0.7497078321908413;

    printf(" -- expectedCube = 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,
                "psAberration returned incorrect values.\n");
        printf("Cube Difference  =  x,y,z  = %.13g, %.13g, %.13g \n\n",
               (x - outCube->x), (y - outCube->y), (z - outCube->z) );
        return 3;
    }

    psFree(outCube);

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

    return 0;
}

psS32 testGravityDeflect(void)
{

    psSphere *actual = psSphereAlloc();
    psSphere *apparent = NULL;
    //    psSphere *sun = psSphereAlloc();
    psSphere *empty = NULL;

    //    sun->r = 0.2;
    //    sun->d = 0.2;
    //    actual->r = 0.61001;
    //    actual->d = 0.01999;
    actual->r = DEG_TO_RAD(122.9153182445501);
    actual->d = DEG_TO_RAD(48.562968978679194);
    psCube *sunCube = psCubeAlloc();
    sunCube->x = 1.467797790127511e11;
    sunCube->y = 2.5880956908748722e10;
    sunCube->z = 1.1220046291457653e10;
    psSphere *sun = psCubeToSphere(sunCube);

    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);
    psSphere *result = psSphereSetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN);

    printf("\nActual r,d = %.13g,%.13g    Apparent r,d = %.16g, %.16g \n",
           actual->r, actual->d, result->r, result->d);
    printf("Sphere Difference  =  r,d  =  %.13g, %.13g\n",
           (actual->r - result->r), (actual->d - result->d));
    psCube *outCube = psSphereToCube(result);
    printf(" -- resultCube = x,y,z = %.13g, %.13g, %.13g  -- \n",
           outCube->x, outCube->y, outCube->z);
    psCube *outCube2 = psSphereToCube(actual);
    printf(" -- actualCube = 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;

    printf(" -- expectedCube = 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("Cube Difference  =  x,y,z  = %.13g, %.13g, %.13g \n",
               (x - outCube->x), (y - outCube->y), (z - outCube->z) );
        return 1;
    }

    psFree(outCube);
    psFree(outCube2);
    psFree(sunCube);

    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
    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
    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");
            printf("\n  Precession Model output = x,y,s = %.8g, %.8g, %.8g\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\n",
                   (pmodel->x - x), (pmodel->y - y), (pmodel->s - s) );
            return 4;
        }
    }
    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;

    //XXX:  MAY Be an Error here... Time for Precession Correction has to be TAI or seg faults
    time2 = psTimeConvert(time2, PS_TIME_TAI);
    //Tests for Precession Correction function//
    //Return NULL for NULL time input
    psEarthPole *pcorr = NULL;
    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
    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 {
        printf("\nPrecession Correction output (IERSA) = x,y,s = %.13g, %.13g, %.13g\n\n",
               pcorr->x, pcorr->y, pcorr->s);
    }
    psFree(pcorr);

    //    printf("\n >>MJD value = %lf \n", psTimeToMJD(time2));

    //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-6;
        yy = SEC_TO_RAD(yy) * 1e-6;
        //        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");
        printf("\nPrecession Correction 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\n", xx, yy, ss);
        //            printf("  A difference of:   %.13g, %.13g, %.13g\n\n",
        //                   (pcorr->x - xx), (pcorr->y - yy), (pcorr->s - ss) );
        //            return 10;
        //        }
    }

    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
    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.
    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;
    }

    //Return valid values for correct input time.  Test IERS Bulletin B.
    polarMotion = psEOC_GetPolarMotion(in, PS_IERS_B);
    double x, y, s;
    x = -6.46014230078e-7;
    y = 2.11194535765e-6;
    s = -1.66256670849e-6;

    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");
        printf("\n  <>PolarMotion output (IERSB)   = x,y,s = %.13g, %.13g, %.13g\n",
               polarMotion->x, polarMotion->y, polarMotion->s);
        printf("\n  <>PolarMotion expected (IERSB) = x,y,s = %.13g, %.13g, %.13g\n",
               x, y, s);
        //        return 3;
    }
    psFree(polarMotion);

    //Test for IERS bulletin A.
    x = -6.46028774489e-7;
    y = 2.11172234336e-6;
    s = -1.66257785921e-6;
    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");
        printf("\n  <>PolarMotion output (IERSA)   = x,y,s = %.13g, %.13g, %.13g\n",
               polarMotion->x, polarMotion->y, polarMotion->s);
        printf("\n  <>PolarMotion expected (IERSA) = x,y,s = %.13g, %.13g, %.13g\n",
               x, y, s);
        //        return 5;
    }



    psEarthPole *nutationCorr = psEOC_NutationCorr(in);
    if (nutationCorr == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, false, "Nutation Correction returned NULL.\n");
        return 6;
    }
    /*
        double sum = sqrt(nutationCorr->x*nutationCorr->x + nutationCorr->y*nutationCorr->y
                          + nutationCorr->s*nutationCorr->s);
        nutationCorr->x = nutationCorr->x / sum;
        nutationCorr->y = nutationCorr->y / sum;
        nutationCorr->s = nutationCorr->s / sum;
    */
    printf("  -- NutationCorr = x,y,s = %.13g, %.13g, %.13g\n\n", nutationCorr->x,
           nutationCorr->y, nutationCorr->s);
    polarMotion->x += nutationCorr->x;
    polarMotion->y += nutationCorr->y;
    polarMotion->s += nutationCorr->s;
    x = -0.13275353774074533;
    y = 0.4359436319739848;
    s = -4.2376965863576153e-10;

    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");
        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("\n  Output Nutation = x,y,s = %.13g, %.13g, %.13g\n",
               nutationCorr->x, nutationCorr->y, nutationCorr->s);
        //        printf("  A difference of:   %.13g, %.13g, %.13g\n\n",
        //               (polarMotion->x - x), (polarMotion->y - y), (polarMotion->s - s) );
        //        return 10;
    }

    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;

    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 {
        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 = 1131579114;
    //    in->nsec = 498489000;
    in->sec = timesec;
    in->nsec = 0;
    in->leapsecond = false;
    psTime *empty = NULL;
    psEarthPole *nute = NULL;

    //Return NULL for NULL input time.
    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 {
        printf("Nutation Correction output = x,y,s = %.8g, %.8g, %.8g\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 *now = psTimeAlloc(PS_TIME_UT1);
    //    now->sec = 1128530000;
    //    now->nsec = 931154510;
    psTime *time = psTimeAlloc(PS_TIME_UT1);
    time->sec = timesec-1;
    time->nsec = 657017200;
    time->leapsecond = false;

    //    psSphereRot *teoceo = psSphereRot_TEOtoCEO(now);

    //return NULL for NULL input time
    rot = psSphereRot_TEOtoCEO(empty);
    if (rot != NULL) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
                "psSphereRot_TEOtoCEO failed to return NULL for NULL input time.\n");
        return 1;
    }

    psSphereRot *teoceo = psSphereRot_TEOtoCEO(time);
    //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);
    if (fabs(rot->q0-teoceo->q0) > DBL_EPSILON || fabs(rot->q1-teoceo->q1) > DBL_EPSILON ||
            fabs(rot->q2-teoceo->q2) > DBL_EPSILON || fabs(rot->q3-teoceo->q3) > DBL_EPSILON) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
                "psSphereRot_TEOtoCEO failed to return matching values for different time types.\n");
        return 2;
    }
    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");
        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(now);
    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 = M_PI / 4.0;
    //    in->y = M_PI / 6.0;
    //    in->s = M_PI / 8.0;

    in->x = 2.857175590089105e-4;
    in->y = 2.3968739377734732e-5;
    in->s = -1.3970066457904322e-8;
    //    in->x += 0.06295703125;
    //    in->y += -0.0287618408203125;
    //    in->s += 0.0;
    //    in->x += 3.05224300720406e-7;
    //    in->y += -1.39441339235822e-7;
    in->s += 0.0;
    //    in->x += -2.491942320903e-10;
    //    in->y += -1.648366515772e-11;

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

    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 return NULL for valid psEarthPole input.\n");
        return 2;
    }

    if (fabs(rot->q0-q0) > DBL_EPSILON || fabs(rot->q1-q1) > DBL_EPSILON ||
            fabs(rot->q2-q2) > DBL_EPSILON || fabs(rot->q3-q3) > DBL_EPSILON) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
                "psSphereRot_CEOtoGCRS 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("  abs difference:   %.13g, %.13g, %.13g, %.13g \n\n", (rot->q0+q0), (rot->q1+q1),
               (rot->q2+q2), (rot->q3+q3) );
        //return 3;
    }


    psCube *tempCube = psCubeAlloc();
    tempCube->x = -0.35963388069046304;
    tempCube->y = 0.5555192509816625;
    tempCube->z = 0.7497078321908413;
    obj = psCubeToSphere(tempCube);
    psFree(tempCube);
    //    obj = psSphereAlloc();
    //    obj->r = objR;
    //    obj->d = objD;
    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;
    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));

    //great sphere difference between result and expected
    //Re = 6.38x10^6m
    psCube *expect = psCubeAlloc();
    expect->x = x;
    expect->y = y;
    expect->z = z;
    psSphere *expected = psCubeToSphere(expect);
    psFree(expect);
    double d = (6.38e6) * acos(cos(result->r)*cos(expected->r)*cos(result->d-expected->d) + sin(result->r)*sin(expected->r));
    printf("\n\nGreat circle difference of:  %.13g \n", d);

    psFree(expected);
    psFree(obj);



    psSphere *test = psSphereAlloc();
    test->r = 0.0;
    test->d = 0.0;
    test = psSphereRotApply(test, rot, test);
    printf("\n  Sphere -test- has values r,d = %.8g, %.8g \n", test->r, test->d);

    psFree(precessionNutation);
    psFree(result);
    psFree(cube);

    psFree(test);
    psFree(rot);
    psFree(in);
    return 0;
}

psS32 testSphereRot_ITRStoTEO(void)
{
    psEarthPole *in = psEarthPoleAlloc();
    psEarthPole *empty = NULL;
    psSphereRot *rot = NULL;

    //    in->x = M_PI / 4.0;
    //    in->y = M_PI / 6.0;
    //    in->s = M_PI / 8.0;
    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);

    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
       ) {
        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) );
        return 3;
    }
    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);

    //    psSphere *test = psSphereAlloc();
    //    test->r = 0.0;
    //    test->d = 0.0;
    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);
    printf("\n  Sphere -test- has values r,d = %.8g, %.8g \n", test->r, test->d);
    temp = psSphereToCube(test);
    printf("\n  Cube -test- has x,y,z =     %.13g, %.13g, %.13g \n", temp->x, temp->y, temp->z);
    double x, y, z;
    x = 0.01698577185310146;
    y = -0.6616538927902393;
    z = 0.7496169753347885;
    printf("\n  Cube -expected- has x,y,z = %.13g,  %.13g, %.13g \n", x, y, z );

    psFree(newRot);
    psFree(obj);
    psFree(temp);
    psFree(test);
    psFree(rot);
    psFree(in);
    return 0;
}

