/** @file  tst_psCoord.c
*
*  @brief The code will ...
*
*
*  @author GLG, MHPCC
*
*  @version $Revision: 1.28 $ $Name: not supported by cvs2svn $
*  @date $Date: 2005-04-07 20:27:41 $
*
*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
*/

#include "psTest.h"
#include "pslib_strict.h"
static psS32 testSphereTransformAlloc(void);
static psS32 testPlaneTransformAlloc(void);
static psS32 testPlaneDistortAlloc(void);
static psS32 testPlaneTransformApply(void);
static psS32 testPlaneDistortApply(void);
static psS32 testSphereTransformApply1(void);
static psS32 testSphereTransformApply2(void);
static psS32 testSphereTransformApply3(void);
static psS32 testSphereTransformApply4(void);
static psS32 testSphereTransformApply5(void);
static psS32 testSphereTransformICRSToEcliptic(void);
static psS32 testSphereTransformEclipticToICRS(void);
static psS32 testSphereTransformICRSToGalactic(void);
static psS32 testSphereTransformGalacticToICRS(void);
static psS32 testSphereTransformPrecess(void);
testDescription tests[] = {
                              {testSphereTransformAlloc, 819, "psSphereTransformAlloc()", 0, false},
                              {testPlaneTransformAlloc, 826, "psPlaneTransformAlloc()", 0, false},
                              {testPlaneDistortAlloc, 827, "psPlaneDistortAlloc()", 0, false},
                              {testPlaneTransformApply, 831, "psPlaneTransformApply()", 0, false},
                              {testPlaneDistortApply, 832, "psPlaneDistortApply()", 0, false},
                              {testSphereTransformApply1, 820, "psSphereTransformApply()", 0, false},
                              {testSphereTransformApply2, 820, "psSphereTransformApply()", 0, false},
                              {testSphereTransformApply3, 820, "psSphereTransformApply()", 0, false},
                              {testSphereTransformApply4, 820, "psSphereTransformApply()", 0, false},
                              {testSphereTransformApply5, 820, "psSphereTransformApply()", 0, false},
                              {testSphereTransformICRSToEcliptic, 821, "psSphereTransformICRSToEcliptic()", 0, false},
                              {testSphereTransformEclipticToICRS, 822, "psSphereTransformEclipticToICRS()", 0, false},
                              {testSphereTransformICRSToGalactic, 824, "psSphereTransformICRSToGalactic()", 0, false},
                              {testSphereTransformGalacticToICRS, 823, "psSphereTransformGalacticToICRS()", 0, false},
                              {testSphereTransformPrecess, 825, "psSphereTransformPrecess()", 0, false},
                              {NULL}
                          };

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

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

#define DEG_INC   30.0

#define MJD_1900  15021.0        // Modified Julian Day 1/1/1900 00:00:00
#define MJD_2000  51544.0        // Modified Julian Day 1/1/2000 00:00:00
#define MJD_2100  88069.0        // Modified Julian Day 1/1/2100 00:00:00

#define ERROR_TOL   0.0001

#define ALPHA_P 1.0
#define DELTA_P 2.0
#define PHI_P 3.0

psS32 testSphereTransformAlloc( void )
{
    // Allocate data structure
    psSphereTransform *myST = psSphereTransformAlloc(ALPHA_P, DELTA_P, PHI_P);
    // Verify null not returned
    if(myST == NULL) {
        psError(PS_ERR_UNKNOWN,true,"Returned NULL with valid parameters");
        return 1;
    }
    // Verify sin member is updated
    psF32 sinDeltaP = sin(DELTA_P);
    if (FLT_EPSILON < fabs(sinDeltaP - myST->sinDeltaP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->sinDeltaP is %f, should be %f\n", myST->sinDeltaP, sinDeltaP);
        return 2;
    }
    // Verify cos member is update
    psF32 cosDeltaP = cos(DELTA_P);
    if (FLT_EPSILON < fabs(cosDeltaP - myST->cosDeltaP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->cosDeltaP is %f, should be %f\n", myST->cosDeltaP, cosDeltaP);
        return 3;
    }
    // Verify alpha is updated
    if (FLT_EPSILON < fabs(ALPHA_P - myST->alphaP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->alphaP is %f, should be %f\n", myST->alphaP, ALPHA_P);
        return 4;
    }
    // Verify phi is updated
    if (FLT_EPSILON < fabs(PHI_P - myST->phiP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->phiP is %f, should be %f\n", myST->phiP, PHI_P);
        return 5;
    }
    // Free data structure
    psFree(myST);

    return 0;
}

#define NW_TERMS 2
#define NX_TERMS 3
#define NY_TERMS 4
#define NZ_TERMS 5

psS32 testPlaneTransformAlloc( void )
{
    // Allocate psPlaneTransform with known x and y terms
    psPlaneTransform *myPT = psPlaneTransformAlloc(NX_TERMS, NY_TERMS);

    // Verify returned value is not NULL
    if(myPT == NULL) {
        psError(PS_ERR_UNKNOWN,true,"Function returned NULL not expected.");
        return 1;
    }

    // Verify returned transform has members set properly
    if (myPT->x->nX != NX_TERMS) {
        psError(PS_ERR_UNKNOWN,true,"myPT->x->nX is %d, should be %d",
                myPT->x->nX, NX_TERMS);
        return 2;
    }
    if (myPT->y->nX != NX_TERMS) {
        psError(PS_ERR_UNKNOWN,true,"myPT->y->nX is %d, should be %d",
                myPT->y->nX, NX_TERMS);
        return 3;
    }
    if (myPT->x->nY != NY_TERMS) {
        psError(PS_ERR_UNKNOWN,true,"myPT->x->nY is %d, should be %d",
                myPT->x->nY, NY_TERMS);
        return 4;
    }
    if (myPT->y->nY != NY_TERMS) {
        psError(PS_ERR_UNKNOWN,true,"myPT->y->nY is %d, should be %d",
                myPT->y->nY, NY_TERMS);
        return 5;
    }

    // Free plane transform
    psFree(myPT);

    // Attempt to specify negative x terms and verify NULL returned and errror message generated
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message: negative x terms");
    myPT = psPlaneTransformAlloc(-1.0, 1.0);
    if (myPT != NULL) {
        psError(PS_ERR_UNKNOWN,true,"psPlaneTransformAlloc() did not return NULL.");
        return 6;
    }

    // Attempt to specify negative y terms and verify NULL returned and error message generated
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message: negative y terms");
    myPT = psPlaneTransformAlloc(1.0, -1.0);
    if (myPT != NULL) {
        psError(PS_ERR_UNKNOWN,true,"psPlaneTransformAlloc() did not return NULL.");
        return 7;
    }

    psFree(myPT);

    return 0;
}


psS32 testPlaneDistortAlloc( void )
{
    // Invoke function with known parameters
    psPlaneDistort *myPD = psPlaneDistortAlloc(NW_TERMS, NX_TERMS, NY_TERMS, NZ_TERMS);

    // Verify NULL is not returned
    if(myPD == NULL) {
        psError(PS_ERR_UNKNOWN,true,"Return of NULL not expected");
        return 1;
    }
    // Verify the terms are properly set after allocation
    if (myPD->x->nW != NW_TERMS) {
        psError(PS_ERR_UNKNOWN,true,"myPD->x->nW is %d, should be %d",
                myPD->x->nW, NW_TERMS);
        return 2;
    }
    if (myPD->y->nW != NW_TERMS) {
        psError(PS_ERR_UNKNOWN,true,"myPD->y->nW is %d, should be %d",
                myPD->y->nW, NW_TERMS);
        return 3;
    }
    if (myPD->x->nX != NX_TERMS) {
        psError(PS_ERR_UNKNOWN,true,"myPD->x->nX is %d, should be %d",
                myPD->x->nX, NX_TERMS);
        return 4;
    }
    if (myPD->y->nX != NX_TERMS) {
        psError(PS_ERR_UNKNOWN,true,"myPD->y->nX is %d, should be %d",
                myPD->y->nX, NX_TERMS);
        return 5;
    }
    if (myPD->x->nY != NY_TERMS) {
        psError(PS_ERR_UNKNOWN,true,"myPD->x->nY is %d, should be %d",
                myPD->x->nY, NY_TERMS);
        return 6;
    }
    if (myPD->y->nY != NY_TERMS) {
        psError(PS_ERR_UNKNOWN,true,"myPD->y->nY is %d, should be %d",
                myPD->y->nY, NY_TERMS);
        return 7;
    }
    if (myPD->x->nZ != NZ_TERMS) {
        psError(PS_ERR_UNKNOWN,true,"myPD->x->nZ is %d, should be %d",
                myPD->x->nZ, NZ_TERMS);
        return 8;
    }
    if (myPD->y->nZ != NZ_TERMS) {
        psError(PS_ERR_UNKNOWN,true,"myPD->y->nZ is %d, should be %d",
                myPD->y->nZ, NZ_TERMS);
        return 9;
    }
    // Free psPlaneTransform
    psFree(myPD);

    // Invoke function with negative x term parameter
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message for negative x term");
    myPD = psPlaneDistortAlloc(-1.0, 1.0, 1.0, 1.0);
    if (myPD != NULL) {
        psError(PS_ERR_UNKNOWN,true,"psPlaneDistortAlloc() did not return NULL.");
        return 10;
    }

    // Invoke function with negative y term parameter
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message for negative y term");
    myPD = psPlaneDistortAlloc(1.0, -1.0, 1.0, 1.0);
    if (myPD != NULL) {
        psError(PS_ERR_UNKNOWN,true,"psPlaneDistortAlloc() did not return NULL.");
        return 11;
    }

    // Invoke function with negative z term parameter
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message for negative z term");
    myPD = psPlaneDistortAlloc(1.0, 1.0, -1.0, 1.0);
    if (myPD != NULL) {
        psError(PS_ERR_UNKNOWN,true,"psPlaneDistortAlloc() did not return NULL.");
        return 12;
    }

    // Invoke function with negative w term parameter
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message for negative w term");
    myPD = psPlaneDistortAlloc(1.0, 1.0, 1.0, -1.0);
    if (myPD != NULL) {
        psError(PS_ERR_UNKNOWN,true,"psPlaneDistortAlloc() did not return NULL.");
    }

    return 0;
}

#define N 10
// We do a simple identity transformation on a few x,y pairs.
psS32 testPlaneTransformApply( void )
{
    psPlane*          out = NULL;
    psDPolynomial2D*  tmp2DPoly = NULL;
    psPlane*          rc;

    // Create input coordinate
    psPlane* in = psPlaneAlloc();
    if(in == NULL) {
        psError(PS_ERR_UNKNOWN,true,"Did not expect NULL return from psPlaneAlloc");
        return 1;
    }

    // Create transform
    psPlaneTransform* pt = psPlaneTransformAlloc(2,2);

    // Set transform coefficients so the x coord input will equal x coord output, same for y
    pt->x->coeff[1][0] = 1.0;
    pt->y->coeff[0][1] = 1.0;

    // Apply transform for several points
    for (psS32 i = 0; i < N; i++) {
        in->x = (psF64) i;
        in->y = (psF64) (i + 5.0);
        in->xErr = 0.0;
        in->yErr = 0.0;

        // Verify function does not return NULL
        out = psPlaneTransformApply(out, pt, in);
        if( out  == NULL) {
            psError(PS_ERR_UNKNOWN,true,"Expected non-NULL return value");
            return 1;
        }
        // Verify the expected values are returned
        if (FLT_EPSILON < fabs(out->x - in->x)) {
            psError(PS_ERR_UNKNOWN,true,"out.x is %lf, should be %lf", out->x, in->x);
            return 3;
        }
        if (FLT_EPSILON < fabs(out->y - in->y)) {
            psError(PS_ERR_UNKNOWN,true,"out.y is %lf, should be %lf", out->y, in->y);
            return 4;
        }
    }
    psFree(out);

    // Verify return null and error message generater for null specified transform
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message for NULL psPlaneTransform");
    rc = psPlaneTransformApply(NULL, NULL, in);
    if (rc != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Function did not return NULL as expected for NULL psPlaneTransform.");
        return 5;
    }

    // Verify return null and error message generated for null psPlaneTransform x coeff
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message for NULL x coeff psPlaneTransform");
    tmp2DPoly = pt->x;
    pt->x = NULL;
    rc = psPlaneTransformApply(NULL, pt, in);
    if (rc != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Function did not return NULL as expected for NULL x coeff psPlaneTransform");
        return 6;
    }
    pt->x = tmp2DPoly;

    // Verify return null and error message generated for null psPlaneTransform y coeff
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message for NULL y coeff psPlaneTransform");
    tmp2DPoly = pt->y;
    pt->y = NULL;
    rc = psPlaneTransformApply(NULL, pt, in);
    if (rc != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Function did not return NULL as expected for NULL y coeff psPlaneTransform");
        return 7;
    }
    pt->y = tmp2DPoly;

    // Verify return null and error message generated for null psPlane coordinate
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message for NULL psPlane");
    rc = psPlaneTransformApply(NULL, pt, NULL);
    if (rc != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Function did not return NULL as expected for NULL psPlane");
        return 8;
    }

    psFree(pt);
    psFree(in);

    return 0;
}

#define COLOR 1.0
#define MAGNITUDE 1.0
// We do a simple identity transformation on a few x,y pairs.  For x and y,
// we add in the COLOR and MAGNITUDE.
psS32 testPlaneDistortApply( void )
{
    psPlane*         out       = NULL;
    psPlane*         rc        = NULL;
    psDPolynomial4D* tmp4DPoly = NULL;

    // Allocate input coordinate
    psPlane* in = psPlaneAlloc();

    // Create psPlaneTransform structure
    psPlaneDistort*  pt = psPlaneDistortAlloc(2, 2, 2, 2);

    pt->x->coeff[1][0][0][0] = 1.0;
    pt->x->coeff[0][0][1][0] = 1.0;
    pt->x->coeff[0][0][0][1] = 1.0;

    pt->y->coeff[0][1][0][0] = 1.0;
    pt->y->coeff[0][0][1][0] = 1.0;
    pt->y->coeff[0][0][0][1] = 1.0;

    for (psS32 i = 0; i < N; i++) {
        in->x = (psF64) i;
        in->y = (psF64) (i + 5.0);
        in->xErr = 0.0;
        in->yErr = 0.0;
        out = psPlaneDistortApply(out, pt, in, COLOR, MAGNITUDE);

        // Verify return value is not NULL
        if(out == NULL) {
            psError(PS_ERR_UNKNOWN,true,"Expected non-NULL return");
            return 1;
        }
        // Verify return value contains expected values
        if (FLT_EPSILON < fabs(out->x - COLOR - MAGNITUDE - in->x)) {
            psError(PS_ERR_UNKNOWN,true,"out->x is %lf, should be %lf",
                    out->x, in->x + COLOR + MAGNITUDE);
            return 2;
        }
        if (FLT_EPSILON < fabs(out->y - COLOR - MAGNITUDE - in->y)) {
            psError(PS_ERR_UNKNOWN,true,"out->y is %lf, should be %lf",
                    out->y, in->y + COLOR + MAGNITUDE);
            return 3;
        }
    }
    psFree(out);

    // Attempt to invoke psPlaneDistortApply with NULL psPlaneDistort
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message for null psPlaneDistort");
    rc = psPlaneDistortApply(NULL, NULL, in, COLOR, MAGNITUDE);
    if (rc != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Did not return NULL.");
        return 4;
    }

    // Attempt to invoke psPlaneDistortApply with NULL psPlaneDistort x member
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message for null x member psPlaneDistort");
    tmp4DPoly = pt->x;
    pt->x = NULL;
    rc = psPlaneDistortApply(NULL, pt, in, COLOR, MAGNITUDE);
    if (rc != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Did not return NULL.");
        return 5;
    }
    pt->x = tmp4DPoly;

    // Attempt to invoke psPlaneDistortApply with NULL psPlaneDistort y member
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message for null y member psPlaneDistort");
    tmp4DPoly = pt->y;
    pt->y = NULL;
    rc = psPlaneDistortApply(NULL, pt, in, COLOR, MAGNITUDE);
    if (rc != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Did not return NULL.");
        return 6;
    }
    pt->y = tmp4DPoly;

    // Attempt to invoke psPlaneDistortApply with NULL input psPlane
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message for null input psPlane");
    rc = psPlaneDistortApply(NULL, pt, NULL, COLOR, MAGNITUDE);
    if (rc != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Did not return NULL");
        return 7;
    }

    psFree(pt);
    psFree(in);

    return 0;
}

// We do a simple identity transformation on a few RA, DEC pairs.
psS32 testSphereTransformApply1( void )
{
    psSphere *in = psSphereAlloc();
    psSphere *out = psSphereAlloc();
    psSphere *temp = NULL;
    psSphere *rc = NULL;
    psSphereTransform *myST = psSphereTransformAlloc(0.0, 0.0, 0.0);

    for (float r=0.0;r<180.0;r+=DEG_INC) {
        for (float d=0.0;d<90.0;d+=DEG_INC) {
            in->r = DEG_TO_RAD(r);
            in->d = DEG_TO_RAD(d);
            in->rErr = 0.0;
            in->dErr = 0.0;

            if(psSphereTransformApply(out, myST, in) != out) {
                psError(PS_ERR_UNKNOWN,true,"Did not return output pointer.");
                return 1;
            }
            if (ERROR_TOL < fabs(out->r - in->r)) {
                psError(PS_ERR_UNKNOWN,true,"out->r is %f, should be %f\n", out->r, in->r);
                return 2;
            }
            if (ERROR_TOL < fabs(out->d - in->d)) {
                psError(PS_ERR_UNKNOWN,true,"out->d is %f, should be %f\n", out->d, in->d);
                return 3;
            }
        }
    }

    // Verify new sphere object is created if out parameter NULL
    temp = psSphereTransformApply(NULL, myST, in);
    if ( temp == NULL) {
        psError(PS_ERR_UNKNOWN,true,"Returned NULL when out parameter was null");
        return 4;
    }
    psFree(temp);

    // Verify NULL returned if transform structure null
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error");
    rc = psSphereTransformApply(NULL, NULL, in);
    if (rc != NULL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply() did not return NULL.");
        return 5;
    }

    // Verify NULL returned when input sphere is NULL
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error");
    rc = psSphereTransformApply(NULL, myST, NULL);
    if (rc != NULL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply() did not return NULL");
        return 6;
    }

    psFree(myST);
    psFree(out);
    psFree(in);

    return 0;
}

/******************************************************************************
testSphereTransformApply2(): This test verifies that psSphereTransformApply()
works properly.  We create two psSphereTransforms: a forward transform and a
reverse transform (which is the mathematical inverse of the forward transform).
We apply both transforms to several spherical coordinates and ensure that the
original input coordinate is obtained after applying both transforms.
 
XXX: We currently test the alpha and delta offsets independently.  Attempts to
test them both concurrently failed.  Determine why this is.  Are the following
spherical transforms not mathematical inverses?
    psSphereTransformAlloc(X, Y, 0.0)
    psSphereTransformAlloc(-X, -Y, 0.0)
 *****************************************************************************/
#define ERROR_PERCENT 0.01
psS32 testSphereTransformApply2( void )
{
    psS32 testStatus = 0;
    psSphere in;
    psSphere out;
    psSphere out2;
    psSphereTransform *mySphereTransformForward = NULL;
    psSphereTransform *mySphereTransformReverse = NULL;


    mySphereTransformForward = psSphereTransformAlloc(DEG_TO_RAD(22.0),
                               0.0,
                               0.0);
    mySphereTransformReverse = psSphereTransformAlloc(DEG_TO_RAD(-22.0),
                               0.0,
                               0.0);

    for (float r=0.1;r<180.0;r+=(DEG_INC/5.0)) {
        for (float d=0.1;d<90.0;d+=(DEG_INC/5.0)) {
            in.r = DEG_TO_RAD(r);
            in.d = DEG_TO_RAD(d);
            in.rErr = 0.0;
            in.dErr = 0.0;

            psSphereTransformApply(&out, mySphereTransformForward, &in);
            psSphereTransformApply(&out2, mySphereTransformReverse, &out);

            if ((fabs((in.r - out2.r) / in.r) > ERROR_PERCENT) ||
                    (fabs((in.d - out2.d) / in.d) > ERROR_PERCENT)) {
                printf("ERROR: \n");
                printf("Input  coords (R, D) are (%f, %f)\n", in.r, in.d);
                printf("Output coords (R, D) are (%f, %f)\n", out2.r, out2.d);
                testStatus = 4;
            }
        }
    }
    psFree(mySphereTransformForward);
    psFree(mySphereTransformReverse);

    mySphereTransformForward = psSphereTransformAlloc(0.0,
                               DEG_TO_RAD(33.0),
                               0.0);
    mySphereTransformReverse = psSphereTransformAlloc(0.0,
                               DEG_TO_RAD(-33.0),
                               0.0);
    for (float r=0.1;r<180.0;r+=(DEG_INC/5.0)) {
        for (float d=0.1;d<90.0;d+=(DEG_INC/5.0)) {
            in.r = DEG_TO_RAD(r);
            in.d = DEG_TO_RAD(d);
            in.rErr = 0.0;
            in.dErr = 0.0;

            psSphereTransformApply(&out, mySphereTransformForward, &in);
            psSphereTransformApply(&out2, mySphereTransformReverse, &out);

            if ((fabs((in.r - out2.r) / in.r) > ERROR_PERCENT) ||
                    (fabs((in.d - out2.d) / in.d) > ERROR_PERCENT)) {
                printf("ERROR: \n");
                printf("Input  coords (R, D) are (%f, %f)\n", in.r, in.d);
                printf("Output coords (R, D) are (%f, %f)\n", out2.r, out2.d);
                testStatus = 4;
            }
        }
    }
    psFree(mySphereTransformForward);
    psFree(mySphereTransformReverse);

    return(testStatus);
}

// Test point 1  ADD  ICRS(0,0) to Galactic (96.337272,-60.188553)
//                              to Ecliptic (0,0) with T=1
#define SPHERE_COORD_TP1_GAL_R           1.681409       //  96.337272 degrees
#define SPHERE_COORD_TP1_GAL_D          -1.050488       // -60.188553 degrees
#define SPHERE_COORD_TP1_ECL_R           0.0
#define SPHERE_COORD_TP1_ECL_D           0.0
#define SPHERE_COORD_TP1_INVERSE_R       0.0
#define SPHERE_COORD_TP1_INVERSE_D       0.0

// Test point 2 ADD  ICRS(0,90) to Galactic (122.93192, 27.12825)
//                              to Ecliptic (90,66.560719) with T=1
#define SPHERE_COORD_TP2_GAL_R           2.145567       // 122.93192 degrees
#define SPHERE_COORD_TP2_GAL_D           0.473477       // 27.12825  degrees
#define SPHERE_COORD_TP2_ECL_R           1.570796       // 90        degrees
#define SPHERE_COORD_TP2_ECL_D           1.161704       // 66.560719 degrees
#define SPHERE_COORD_TP2_INVERSE_R       0.0            // 0         degrees
#define SPHERE_COORD_TP2_INVERSE_D       1.570796       // 90        degrees

// Test point 3 ADD  ICRS(180,30) to Galactic (195.639488, 78.353806)
//                                to Ecliptic (167.072470, 27.308813) with T=2
#define SPHERE_COORD_TP3_GAL_R           3.414553       // 195.639488 degrees
#define SPHERE_COORD_TP3_GAL_D           1.367532       // 78.353806  degrees
#define SPHERE_COORD_TP3_ECL_R           2.915965       // 167.07247  degrees
#define SPHERE_COORD_TP3_ECL_D           0.476629       // 27.308813  degrees
#define SPHERE_COORD_TP3_INVERSE_R       3.141593       // 180.0      degrees
#define SPHERE_COORD_TP3_INVERSE_D       0.523599       // 30.0       degrees

psS32 testSphereTransformApply3( void)
{

    psSphereTransform*  testTransform;
    psSphere*           inputCoord = psSphereAlloc();
    psSphere*           outputCoord = NULL;
    psSphere*           inverseOutputCoord = NULL;

    inputCoord->r = SPHERE_COORD_TP1_INVERSE_R;
    inputCoord->d = SPHERE_COORD_TP1_INVERSE_D;
    inputCoord->rErr = 0.0;
    inputCoord->dErr = 0.0;

    // Set test date time to 1/1/2000 00:00:00
    psTime* testDateTime = psTimeFromMJD(MJD_2000);

    // Create transform ICRS to Galactic
    testTransform = psSphereTransformICRSToGalactic();

    // Perform transform from ICRS to Galactic
    outputCoord = psSphereTransformApply(outputCoord,testTransform,inputCoord);
    // Verify output not NULL and member values set correctly
    if(outputCoord == NULL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply returned NULL");
        return 1;
    }
    if(fabs(outputCoord->r - SPHERE_COORD_TP1_GAL_R) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply return coord r %lg not as expected %lg",
                outputCoord->r, SPHERE_COORD_TP1_GAL_R);
        return 2;
    }
    if(fabs(outputCoord->d - SPHERE_COORD_TP1_GAL_D) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply return coord d %ld not as expected %lg",
                outputCoord->d, SPHERE_COORD_TP1_GAL_D);
        return 3;
    }
    psFree(testTransform);

    // Create inverse transform Galactic to ICRS
    testTransform = psSphereTransformGalacticToICRS();

    // Perform transform from Galactic to ICRS
    inverseOutputCoord = psSphereTransformApply(inverseOutputCoord,testTransform,outputCoord);
    // Verify output not NULL and member values set correctly
    if(inverseOutputCoord == NULL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply returned NULL");
        return 4;
    }

    if(fabs(inverseOutputCoord->r - SPHERE_COORD_TP1_INVERSE_R) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply return coord r %lg not as expected %lg",
                inverseOutputCoord->r, SPHERE_COORD_TP1_INVERSE_R);
        return 5;
    }
    if(fabs(inverseOutputCoord->d - SPHERE_COORD_TP1_INVERSE_D) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply return coord d %ld not as expected %lg",
                inverseOutputCoord->d, SPHERE_COORD_TP1_INVERSE_D);
        return 6;
    }
    psFree(testTransform);

    // Create transform ICRS to Ecliptic
    testTransform = psSphereTransformICRSToEcliptic(testDateTime);

    // Perform transform from ICRS to Ecliptic
    outputCoord = psSphereTransformApply(outputCoord,testTransform,inputCoord);
    // Verify output not NULL and member values set correctly
    if(outputCoord == NULL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply returned NULL");
        return 7;
    }
    if(fabs(outputCoord->r - SPHERE_COORD_TP1_ECL_R) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply return coord r %lg not as expected %lg",
                outputCoord->r, SPHERE_COORD_TP1_ECL_R);
        return 8;
    }
    if(fabs(outputCoord->d - SPHERE_COORD_TP1_ECL_D) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply return coord d %ld not as expected %lg",
                outputCoord->d, SPHERE_COORD_TP1_ECL_D);
        return 9;
    }
    psFree(testTransform);

    // Create inverse transform Ecliptic to ICRS
    testTransform = psSphereTransformEclipticToICRS(testDateTime);

    // Perform transform from Ecliptic to ICRS
    inverseOutputCoord = psSphereTransformApply(inverseOutputCoord,testTransform,outputCoord);
    // Verify output not NULL and member values set correctly
    if(inverseOutputCoord == NULL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply returned NULL");
        return 10;
    }
    if(fabs(inverseOutputCoord->r - SPHERE_COORD_TP1_INVERSE_R) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply return coord r %lg not as expected %lg",
                inverseOutputCoord->r, SPHERE_COORD_TP1_INVERSE_R);
        return 11;
    }
    if(fabs(inverseOutputCoord->d - SPHERE_COORD_TP1_INVERSE_D) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply return coord d %ld not as expected %lg",
                inverseOutputCoord->d, SPHERE_COORD_TP1_INVERSE_D);
        return 12;
    }
    psFree(testTransform);

    psFree(testDateTime);
    psFree(outputCoord);
    psFree(inputCoord);
    psFree(inverseOutputCoord);

    return 0;
}

psS32 testSphereTransformApply4( void)
{

    psSphereTransform*  testTransform;
    psSphere*           inputCoord = psSphereAlloc();
    psSphere*           outputCoord = NULL;
    psSphere*           inverseOutputCoord = NULL;

    inputCoord->r = SPHERE_COORD_TP2_INVERSE_R;
    inputCoord->d = SPHERE_COORD_TP2_INVERSE_D;
    inputCoord->rErr = 0.0;
    inputCoord->dErr = 0.0;

    // Set test date time to 1/1/2000 00:00:00
    psTime* testDateTime = psTimeFromMJD(MJD_2000);

    // Create transform ICRS to Galactic
    testTransform = psSphereTransformICRSToGalactic();

    // Perform transform from ICRS to Galactic
    outputCoord = psSphereTransformApply(outputCoord,testTransform,inputCoord);
    // Verify output not NULL and member values set correctly
    if(outputCoord == NULL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply returned NULL");
        return 1;
    }
    if(fabs(outputCoord->r - SPHERE_COORD_TP2_GAL_R) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply return coord r %lg not as expected %lg",
                outputCoord->r, SPHERE_COORD_TP2_GAL_R);
        return 2;
    }
    if(fabs(outputCoord->d - SPHERE_COORD_TP2_GAL_D) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply return coord d %ld not as expected %lg",
                outputCoord->d, SPHERE_COORD_TP2_GAL_D);
        return 3;
    }
    psFree(testTransform);

    // Create inverse transform Galactic to ICRS
    testTransform = psSphereTransformGalacticToICRS();

    // Perform transform from Galactic to ICRS
    inverseOutputCoord = psSphereTransformApply(inverseOutputCoord,testTransform,outputCoord);
    // Verify output not NULL and member values set correctly
    if(inverseOutputCoord == NULL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply returned NULL");
        return 4;
    }
    // atan2 function introduces some error with y close to zero, increased error tolerance
    if(fabs(inverseOutputCoord->r - SPHERE_COORD_TP2_INVERSE_R) > ERROR_TOL+0.1) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply return coord r %lg not as expected %lg",
                inverseOutputCoord->r, SPHERE_COORD_TP2_INVERSE_R);
        return 5;
    }
    if(fabs(inverseOutputCoord->d - SPHERE_COORD_TP2_INVERSE_D) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply return coord d %ld not as expected %lg",
                inverseOutputCoord->d, SPHERE_COORD_TP2_INVERSE_D);
        return 6;
    }
    psFree(testTransform);

    // Create transform ICRS to Ecliptic
    testTransform = psSphereTransformICRSToEcliptic(testDateTime);

    // Perform transform from ICRS to Ecliptic
    outputCoord = psSphereTransformApply(outputCoord,testTransform,inputCoord);
    // Verify output not NULL and member values set correctly
    if(outputCoord == NULL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply returned NULL");
        return 7;
    }
    if(fabs(outputCoord->r - SPHERE_COORD_TP2_ECL_R) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply return coord r %lg not as expected %lg",
                outputCoord->r, SPHERE_COORD_TP2_ECL_R);
        return 8;
    }
    if(fabs(outputCoord->d - SPHERE_COORD_TP2_ECL_D) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply return coord d %ld not as expected %lg",
                outputCoord->d, SPHERE_COORD_TP2_ECL_D);
        return 9;
    }
    psFree(testTransform);

    // Create inverse transform Ecliptic to ICRS
    testTransform = psSphereTransformEclipticToICRS(testDateTime);

    // Perform transform from Ecliptic to ICRS
    inverseOutputCoord = psSphereTransformApply(inverseOutputCoord,testTransform,outputCoord);
    // Verify output not NULL and member values set correctly
    if(inverseOutputCoord == NULL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply returned NULL");
        return 10;
    }
    if(fabs(inverseOutputCoord->r - SPHERE_COORD_TP2_INVERSE_R) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply return coord r %lg not as expected %lg",
                inverseOutputCoord->r, SPHERE_COORD_TP2_INVERSE_R);
        return 11;
    }
    if(fabs(inverseOutputCoord->d - SPHERE_COORD_TP2_INVERSE_D) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply return coord d %ld not as expected %lg",
                inverseOutputCoord->d, SPHERE_COORD_TP2_INVERSE_D);
        return 12;
    }
    psFree(testTransform);

    psFree(testDateTime);
    psFree(outputCoord);
    psFree(inputCoord);
    psFree(inverseOutputCoord);

    return 0;
}

psS32 testSphereTransformApply5( void)
{
    psSphereTransform*  testTransform;
    psSphere*           inputCoord = psSphereAlloc();
    psSphere*           outputCoord = NULL;
    psSphere*           inverseOutputCoord = NULL;

    inputCoord->r = SPHERE_COORD_TP3_INVERSE_R;
    inputCoord->d = SPHERE_COORD_TP3_INVERSE_D;
    inputCoord->rErr = 0.0;
    inputCoord->dErr = 0.0;

    // Set test date time to 1/1/2100 00:00:00
    psTime* testDateTime = psTimeFromMJD(MJD_2100);

    // Create transform ICRS to Galactic
    testTransform = psSphereTransformICRSToGalactic();

    // Perform transform from ICRS to Galactic
    outputCoord = psSphereTransformApply(outputCoord,testTransform,inputCoord);
    // Verify output not NULL and member values set correctly
    if(outputCoord == NULL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply returned NULL");
        return 1;
    }
    if(fabs(outputCoord->r - SPHERE_COORD_TP3_GAL_R) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply return coord r %lg not as expected %lg",
                outputCoord->r, SPHERE_COORD_TP3_GAL_R);
        return 2;
    }
    if(fabs(outputCoord->d - SPHERE_COORD_TP3_GAL_D) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply return coord d %ld not as expected %lg",
                outputCoord->d, SPHERE_COORD_TP3_GAL_D);
        return 3;
    }
    psFree(testTransform);

    // Create inverse transform Galactic to ICRS
    testTransform = psSphereTransformGalacticToICRS();

    // Perform transform from Galactic to ICRS
    inverseOutputCoord = psSphereTransformApply(inverseOutputCoord,testTransform,outputCoord);
    // Verify output not NULL and member values set correctly
    if(inverseOutputCoord == NULL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply returned NULL");
        return 4;
    }
    if(fabs(inverseOutputCoord->r - SPHERE_COORD_TP3_INVERSE_R) > ERROR_TOL+0.1) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply return coord r %lg not as expected %lg",
                inverseOutputCoord->r, SPHERE_COORD_TP3_INVERSE_R);
        return 5;
    }
    if(fabs(inverseOutputCoord->d - SPHERE_COORD_TP3_INVERSE_D) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply return coord d %ld not as expected %lg",
                inverseOutputCoord->d, SPHERE_COORD_TP3_INVERSE_D);
        return 6;
    }
    psFree(testTransform);

    // Create transform ICRS to Ecliptic
    testTransform = psSphereTransformICRSToEcliptic(testDateTime);

    // Perform transform from ICRS to Ecliptic
    outputCoord = psSphereTransformApply(outputCoord,testTransform,inputCoord);
    // Verify output not NULL and member values set correctly
    if(outputCoord == NULL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply returned NULL");
        return 7;
    }
    if(fabs(outputCoord->r - SPHERE_COORD_TP3_ECL_R) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply return coord r %lg not as expected %lg",
                outputCoord->r, SPHERE_COORD_TP3_ECL_R);
        return 8;
    }
    if(fabs(outputCoord->d - SPHERE_COORD_TP3_ECL_D) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply return coord d %ld not as expected %lg",
                outputCoord->d, SPHERE_COORD_TP3_ECL_D);
        return 9;
    }
    psFree(testTransform);

    // Create inverse transform Ecliptic to ICRS
    testTransform = psSphereTransformEclipticToICRS(testDateTime);

    // Perform transform from Ecliptic to ICRS
    inverseOutputCoord = psSphereTransformApply(inverseOutputCoord,testTransform,outputCoord);
    // Verify output not NULL and member values set correctly
    if(inverseOutputCoord == NULL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply returned NULL");
        return 10;
    }
    if(fabs(inverseOutputCoord->r - SPHERE_COORD_TP3_INVERSE_R) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply return coord r %lg not as expected %lg",
                inverseOutputCoord->r, SPHERE_COORD_TP3_INVERSE_R);
        return 11;
    }
    if(fabs(inverseOutputCoord->d - SPHERE_COORD_TP3_INVERSE_D) > ERROR_TOL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformApply return coord d %ld not as expected %lg",
                inverseOutputCoord->d, SPHERE_COORD_TP3_INVERSE_D);
        return 12;
    }
    psFree(testTransform);

    psFree(testDateTime);
    psFree(outputCoord);
    psFree(inputCoord);
    psFree(inverseOutputCoord);

    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

psS32 testSphereTransformPrecess( void )
{
    psSphere*     inputCoord  = psSphereAlloc();
    psSphere*     outputCoord = NULL;
    psTime*       fromTime    = psTimeFromMJD(MJD_2100);
    psTime*       toTime      = psTimeFromMJD(MJD_1900);

    // 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
    outputCoord = psSpherePrecess(inputCoord, fromTime, toTime);
    // 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);

    // 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
    outputCoord = psSpherePrecess(inputCoord, fromTime, toTime);
    // 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);

    // 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
    outputCoord = psSpherePrecess(inputCoord, fromTime, toTime);
    // 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);

    // Invoke precess with invalid parameter
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message");
    outputCoord = psSpherePrecess(inputCoord, fromTime, NULL);
    if(outputCoord != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Did not return NULL with invalid input");
        return 10;
    }

    // Invoke precess with invalid parameter
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message");
    outputCoord = psSpherePrecess(inputCoord, NULL, toTime);
    if(outputCoord != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Did not return NULL with invalid input");
        return 11;
    }

    // Invoke precess with invalid parameter
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message");
    outputCoord = psSpherePrecess(NULL, fromTime, toTime);
    if(outputCoord != NULL) {
        psError(PS_ERR_UNKNOWN,true,"Did not return NULL with invalid input");
        return 12;
    }

    // Free objects
    psFree(fromTime);
    psFree(toTime);
    psFree(inputCoord);

    return 0;
}

// Thes values calculated from ADD-09 equations
#define EXPECT_COS_DELTAP_2000    0.917482
#define EXPECT_SIN_DELTAP_2000    0.397777
#define EXPECT_COS_DELTAP_2100    0.917572
#define EXPECT_SIN_DELTAP_2100    0.397567

psS32 testSphereTransformICRSToEcliptic( void )
{
    psF64  expectedPhiP = 0.0;
    psF64  expectedAlphaP = 0.0;
    psF64  expectedCosDeltaP = EXPECT_COS_DELTAP_2000;
    psF64  expectedSinDeltaP = EXPECT_SIN_DELTAP_2000;

    // Set test date time to 1/1/2000 00:00:00
    psTime* testDateTime = psTimeFromMJD(MJD_2000);

    // Invoke function to set psSphereTransforma objec
    psSphereTransform *myST = psSphereTransformICRSToEcliptic(testDateTime);

    // Verify expected values for the specified time
    if (ERROR_TOL < fabs(expectedSinDeltaP - myST->sinDeltaP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->sinDeltaP is %f, should be %f",
                myST->sinDeltaP, expectedSinDeltaP);
        return 1;
    }
    if (ERROR_TOL < fabs(expectedCosDeltaP - myST->cosDeltaP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->cosDeltaP is %f, should be %f",
                myST->cosDeltaP, expectedCosDeltaP);
        return 2;
    }
    if (ERROR_TOL < fabs(expectedAlphaP - myST->alphaP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->alphaP is %f, should be %f",
                myST->alphaP, expectedAlphaP);
        return 3;
    }
    if (ERROR_TOL < fabs(expectedPhiP - myST->phiP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->phiP is %f, should be %f",
                myST->phiP, expectedPhiP);
        return 4;
    }
    psFree(myST);
    psFree(testDateTime);

    expectedCosDeltaP = EXPECT_COS_DELTAP_2100;
    expectedSinDeltaP = EXPECT_SIN_DELTAP_2100;

    // Set test date time to 1/1/2100 00:00:00
    testDateTime = psTimeFromMJD(MJD_2100);

    // Invoke function to set psSphereTransforma objec
    myST = psSphereTransformICRSToEcliptic(testDateTime);

    // Verify expected values for the specified time
    if (ERROR_TOL < fabs(expectedSinDeltaP - myST->sinDeltaP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->sinDeltaP is %f, should be %f",
                myST->sinDeltaP, expectedSinDeltaP);
        return 5;
    }
    if (ERROR_TOL < fabs(expectedCosDeltaP - myST->cosDeltaP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->cosDeltaP is %f, should be %f",
                myST->cosDeltaP, expectedCosDeltaP);
        return 6;
    }
    if (ERROR_TOL < fabs(expectedAlphaP - myST->alphaP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->alphaP is %f, should be %f",
                myST->alphaP, expectedAlphaP);
        return 7;
    }
    if (ERROR_TOL < fabs(expectedPhiP - myST->phiP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->phiP is %f, should be %f",
                myST->phiP, expectedPhiP);
        return 8;
    }
    psFree(myST);
    psFree(testDateTime);

    // Verify if argument psTime is NULL, function returns NULL and error message generated
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message");
    myST = psSphereTransformICRSToEcliptic(NULL);
    if (myST != NULL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformICRSToEcliptic() did not return NULL.");
        return 9;
    }

    // Verify if argument psTime is less than 1900, function return NULL and error msg generated
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message");
    testDateTime = psTimeFromMJD(MJD_1900 - 1.0);
    myST = psSphereTransformICRSToEcliptic(testDateTime);
    if (myST != NULL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformICRSToEcliptic() did not return NULL.");
        return 10;
    }
    psFree(testDateTime);

    return 0;
}


psS32 testSphereTransformEclipticToICRS( void )
{
    psF64  expectedPhiP = 0.0;
    psF64  expectedAlphaP = 0.0;
    psF64  expectedCosDeltaP = EXPECT_COS_DELTAP_2000;
    psF64  expectedSinDeltaP = -EXPECT_SIN_DELTAP_2000;

    // Set test date time to 1/1/2000 00:00:00
    psTime* testDateTime = psTimeFromMJD(MJD_2000);

    // Invoke function to set psSphereTransforma object
    psSphereTransform *myST = psSphereTransformEclipticToICRS(testDateTime);

    // Verify expected values for the specified time
    if (ERROR_TOL < fabs(expectedSinDeltaP - myST->sinDeltaP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->sinDeltaP is %f, should be %f",
                myST->sinDeltaP, expectedSinDeltaP);
        return 1;
    }
    if (ERROR_TOL < fabs(expectedCosDeltaP - myST->cosDeltaP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->cosDeltaP is %f, should be %f",
                myST->cosDeltaP, expectedCosDeltaP);
        return 2;
    }
    if (ERROR_TOL < fabs(expectedAlphaP - myST->alphaP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->alphaP is %f, should be %f",
                myST->alphaP, expectedAlphaP);
        return 3;
    }
    if (ERROR_TOL < fabs(expectedPhiP - myST->phiP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->phiP is %f, should be %f",
                myST->phiP, expectedPhiP);
        return 4;
    }
    psFree(myST);
    psFree(testDateTime);

    expectedCosDeltaP = EXPECT_COS_DELTAP_2100;
    expectedSinDeltaP = -EXPECT_SIN_DELTAP_2100;

    // Set test date time to 1/1/2100 00:00:00
    testDateTime = psTimeFromMJD(MJD_2100);

    // Invoke function to set psSphereTransforma objec
    myST = psSphereTransformEclipticToICRS(testDateTime);

    // Verify expected values for the specified time
    if (ERROR_TOL < fabs(expectedSinDeltaP - myST->sinDeltaP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->sinDeltaP is %f, should be %f",
                myST->sinDeltaP, expectedSinDeltaP);
        return 5;
    }
    if (ERROR_TOL < fabs(expectedCosDeltaP - myST->cosDeltaP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->cosDeltaP is %f, should be %f",
                myST->cosDeltaP, expectedCosDeltaP);
        return 6;
    }
    if (ERROR_TOL < fabs(expectedAlphaP - myST->alphaP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->alphaP is %f, should be %f",
                myST->alphaP, expectedAlphaP);
        return 7;
    }
    if (ERROR_TOL < fabs(expectedPhiP - myST->phiP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->phiP is %f, should be %f",
                myST->phiP, expectedPhiP);
        return 8;
    }
    psFree(myST);
    psFree(testDateTime);

    // Verify if argument psTime is NULL, function returns NULL and error message generated
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message");
    myST = psSphereTransformEclipticToICRS(NULL);
    if (myST != NULL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformEclipticICRS() did not return NULL.");
        return 9;
    }

    // Verify if argument psTime is less than 1900, function return NULL and error msg generated
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message");
    testDateTime = psTimeFromMJD(MJD_1900 - 1.0);
    myST = psSphereTransformEclipticToICRS(testDateTime);
    if (myST != NULL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformEclipticToICRS() did not return NULL.");
        return 10;
    }
    psFree(testDateTime);

    return 0;
}

#define ITG_EXPECT_ALPHAP      4.93683
#define ITG_EXPECT_SINDELTAP   0.88999
#define ITG_EXPECT_COSDELTAP   0.45598
#define ITG_EXPECT_PHIP        0.57477

psS32 testSphereTransformICRSToGalactic( void )
{
    // Invoke function
    psSphereTransform *myST = psSphereTransformICRSToGalactic();

    // Verify did not return NULL
    if(myST == NULL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformICRSToGalactic returned NULL");
        return 1;
    }
    // Verify returned structure has members set properly
    if (ERROR_TOL < fabs(ITG_EXPECT_SINDELTAP - myST->sinDeltaP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->sinDeltaP is %f, should be %f",
                myST->sinDeltaP, ITG_EXPECT_SINDELTAP);
        return 2;
    }
    if (ERROR_TOL < fabs(ITG_EXPECT_COSDELTAP - myST->cosDeltaP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->cosDeltaP is %f, should be %f",
                myST->cosDeltaP, ITG_EXPECT_COSDELTAP);
        return 3;
    }
    if (ERROR_TOL < fabs(ITG_EXPECT_ALPHAP - myST->alphaP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->alphaP is %f, should be %f",
                myST->alphaP, ITG_EXPECT_ALPHAP);
        return 4;
    }
    if (ERROR_TOL < fabs(ITG_EXPECT_PHIP - myST->phiP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->phiP is %f, should be %f",
                myST->phiP, ITG_EXPECT_PHIP);
        return 5;
    }
    // Free data structure created
    psFree(myST);

    return 0;
}

#define GTI_EXPECT_PHIP        4.93683
#define GTI_EXPECT_SINDELTAP  -0.88999
#define GTI_EXPECT_COSDELTAP   0.45598
#define GTI_EXPECT_ALPHAP      0.57477

psS32 testSphereTransformGalacticToICRS( void )
{
    // Invoke function
    psSphereTransform *myST = psSphereTransformGalacticToICRS();

    // Verify did not return NULL
    if(myST == NULL) {
        psError(PS_ERR_UNKNOWN,true,"psSphereTransformGalacticToICRS returned NULL");
        return 1;
    }
    // Verify returned structure has members set properly
    if (ERROR_TOL < fabs(GTI_EXPECT_SINDELTAP - myST->sinDeltaP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->sinDeltaP is %f, should be %f",
                myST->sinDeltaP, GTI_EXPECT_SINDELTAP);
        return 2;
    }
    if (ERROR_TOL < fabs(GTI_EXPECT_COSDELTAP - myST->cosDeltaP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->cosDeltaP is %f, should be %f",
                myST->cosDeltaP, GTI_EXPECT_COSDELTAP);
        return 3;
    }
    if (ERROR_TOL < fabs(GTI_EXPECT_ALPHAP - myST->alphaP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->alphaP is %f, should be %f",
                myST->alphaP, GTI_EXPECT_ALPHAP);
        return 4;
    }
    if (ERROR_TOL < fabs(GTI_EXPECT_PHIP - myST->phiP)) {
        psError(PS_ERR_UNKNOWN,true,"myST->phiP is %f, should be %f",
                myST->phiP, GTI_EXPECT_PHIP);
        return 5;
    }
    // Free data structure created
    psFree(myST);

    return 0;
}

