/** @file  tst_psCoord.c
*
*  @brief The code will ...
*
*
*  @author GLG, MHPCC
*
*  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
*  @date $Date: 2005-01-14 23:25:39 $
*
*  XXX: Must test psSpherePrecess.
*  XXX: THe psProject() and psDeproject() functions do not work fully.  They
*  are not tested here.  The SphereGet/SetOffset() code that rely on them
*  is not tested.
*
*  XXX: Test with a bad offfset unit type in the SphereGet/SetOffset() functions.
*
*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
*/

#include "psTest.h"
#include "pslib.h"
static psS32 test1( void );
static psS32 test1b( void );
static psS32 test1c( void );
static psS32 test2( void );
static psS32 test3( void );
static psS32 test4( void );
static psS32 test4b( void );
static psS32 test5( void );
static psS32 test6( void );
static psS32 test7( void );
static psS32 test8( void );
//static psS32 test20( void );
//static psS32 test21( void );
static psS32 test40( void );
static psS32 test41( void );
//static psS32 test42( void );
//static psS32 test43( void );
testDescription tests[] = {
                              {test1, 0000, "psSphereTransformAlloc()", 0, false},
                              {test1b, 0000, "psPlaneTransformAlloc()", 0, false},
                              {test1c, 0000, "psPlaneDistortAlloc()", 0, false},
                              {test2, 0000, "psPlaneTransformApply()", 0, false},
                              {test3, 0000, "psPlaneDistortApply()", 0, false},
                              {test4, 0000, "psPSphereTransformApply()", 0, false},
                              {test4b, 0000, "psPSphereTransformApply()", 0, false},
                              {test5, 0000, "psSphereTransformICRSToEcliptic()", 0, false},
                              {test6, 0000, "psSphereTransformEclipticToICRS()", 0, false},
                              {test7, 0000, "psSphereTransformICRSToGalatic()", 0, false},
                              {test8, 0000, "psSphereTransformGalaticToICRS()", 0, false},
                              //         {test20, 0000, "psProject()", 0, false},
                              //         {test21, 0000, "psDeProject()", 0, false},
                              {test40, 0000, "psSphereGetOffset()", 0, false},
                              {test41, 0000, "psSphereSetOffset()", 0, false},
                              //         {test42, 0000, "psProject(), psDeproject", 0, false},
                              //         {test43, 0000, "psProject(), psDeproject", 0, false},
                              {NULL}
                          };

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

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

#define MY_TINY 0.0001
#define PS_COMPARE_TINY_THEN_PRINT_ERROR(ACTUAL, EXPECT, TESTSTATUS) \
if (MY_TINY < fabs(EXPECT - ACTUAL)) { \
    printf("ERROR: %s is %f, should be %f (testStatus = %d)\n", #ACTUAL, ACTUAL, EXPECT, TESTSTATUS); \
    testStatus = TESTSTATUS; \
}
#define MY_TINY 0.0001
#define PS_COMPARE_EPSILON_THEN_PRINT_ERROR(ACTUAL, EXPECT, TESTSTATUS) \
if (FLT_EPSILON < fabs(EXPECT - ACTUAL)) { \
    printf("ERROR: %s is %f, should be %f (testStatus = %d)\n", #ACTUAL, ACTUAL, EXPECT, TESTSTATUS); \
    testStatus = TESTSTATUS; \
}

#define PS_COT(X) (1.0 / atan(X))
float p_psArg(float x,
              float y)
{
    return((float) atan2((double) x, (double) y));
}

#define ALPHA_P 1.0
#define DELTA_P 2.0
#define PHI_P 3.0
psS32 test1( void )
{
    psS32 testStatus = 0;
    psSphereTransform *myST = psSphereTransformAlloc(ALPHA_P, DELTA_P, PHI_P);

    psF32 sinDeltaP = sin(DELTA_P);
    if (FLT_EPSILON < fabs(sinDeltaP - myST->sinDeltaP)) {
        printf("ERROR: myST->sinDeltaP is %f, should be %f\n", myST->sinDeltaP, sinDeltaP);
        testStatus = 1;
    }
    psF32 cosDeltaP = cos(DELTA_P);
    if (FLT_EPSILON < fabs(cosDeltaP - myST->cosDeltaP)) {
        printf("ERROR: myST->cosDeltaP is %f, should be %f\n", myST->cosDeltaP, cosDeltaP);
        testStatus = 1;
    }

    if (FLT_EPSILON < fabs(ALPHA_P - myST->alphaP)) {
        printf("ERROR: myST->alphaP is %f, should be %f\n", myST->alphaP, ALPHA_P);
        testStatus = 3;
    }

    if (FLT_EPSILON < fabs(PHI_P - myST->phiP)) {
        printf("ERROR: myST->phiP is %f, should be %f\n", myST->phiP, PHI_P);
        testStatus = 4;
    }

    psFree(myST);
    return(testStatus);
}

#define NW_TERMS 2
#define NX_TERMS 3
#define NY_TERMS 4
#define NZ_TERMS 5
psS32 test1b( void )
{
    psS32 testStatus = 0;
    psPlaneTransform *myPT = psPlaneTransformAlloc(NX_TERMS, NY_TERMS);

    if (myPT->x->nX != NX_TERMS) {
        printf("ERROR: myPT->x->nX is %d, should be %d\n", myPT->x->nX, NX_TERMS);
        testStatus = 1;
    }
    if (myPT->y->nX != NX_TERMS) {
        printf("ERROR: myPT->y->nX is %d, should be %d\n", myPT->y->nX, NX_TERMS);
        testStatus = 1;
    }

    if (myPT->x->nY != NY_TERMS) {
        printf("ERROR: myPT->x->nY is %d, should be %d\n", myPT->x->nY, NY_TERMS);
        testStatus = 1;
    }
    if (myPT->y->nY != NY_TERMS) {
        printf("ERROR: myPT->y->nY is %d, should be %d\n", myPT->y->nY, NY_TERMS);
        testStatus = 1;
    }
    psFree(myPT);

    printf("-------------------------------------------------------------------\n");
    printf("Calling psPlaneTransformAlloc() with negative-order polynomial request.  Should generate error, return NULL.\n");
    myPT = psPlaneTransformAlloc(-1.0, 1.0);
    if (myPT != NULL) {
        printf("TEST ERROR: psPlaneTransformAlloc() did not return NULL.\n");
        testStatus = false;
        myPT = NULL;
    }

    printf("-------------------------------------------------------------------\n");
    printf("Calling psPlaneTransformAlloc() with negative-order polynomial request.  Should generate error, return NULL.\n");
    myPT = psPlaneTransformAlloc(1.0, -1.0);
    if (myPT != NULL) {
        printf("TEST ERROR: psPlaneTransformAlloc() did not return NULL.\n");
        testStatus = false;
        myPT = NULL;
    }
    printf("-------------------------------------------------------------------\n");


    psFree(myPT);
    return(testStatus);
}
psS32 test1c( void )
{
    psS32 testStatus = 0;
    psPlaneDistort *myPD = psPlaneDistortAlloc(NW_TERMS, NX_TERMS, NY_TERMS, NZ_TERMS);

    if (myPD->x->nW != NW_TERMS) {
        printf("ERROR: myPD->x->nW is %d, should be %d\n", myPD->x->nW, NW_TERMS);
        testStatus = 1;
    }
    if (myPD->y->nW != NW_TERMS) {
        printf("ERROR: myPD->y->nW is %d, should be %d\n", myPD->y->nW, NW_TERMS);
        testStatus = 1;
    }

    if (myPD->x->nX != NX_TERMS) {
        printf("ERROR: myPD->x->nX is %d, should be %d\n", myPD->x->nX, NX_TERMS);
        testStatus = 1;
    }
    if (myPD->y->nX != NX_TERMS) {
        printf("ERROR: myPD->y->nX is %d, should be %d\n", myPD->y->nX, NX_TERMS);
        testStatus = 1;
    }

    if (myPD->x->nY != NY_TERMS) {
        printf("ERROR: myPD->x->nY is %d, should be %d\n", myPD->x->nY, NY_TERMS);
        testStatus = 1;
    }
    if (myPD->y->nY != NY_TERMS) {
        printf("ERROR: myPD->y->nY is %d, should be %d\n", myPD->y->nY, NY_TERMS);
        testStatus = 1;
    }

    if (myPD->x->nZ != NZ_TERMS) {
        printf("ERROR: myPD->x->nZ is %d, should be %d\n", myPD->x->nZ, NZ_TERMS);
        testStatus = 1;
    }
    if (myPD->y->nZ != NZ_TERMS) {
        printf("ERROR: myPD->y->nZ is %d, should be %d\n", myPD->y->nZ, NZ_TERMS);
        testStatus = 1;
    }

    psFree(myPD);

    printf("-------------------------------------------------------------------\n");
    printf("Calling psPlaneDistortAlloc() with negative-order polynomial request.  Should generate error, return NULL.\n");
    myPD = psPlaneDistortAlloc(-1.0, 1.0, 1.0, 1.0);
    if (myPD != NULL) {
        printf("TEST ERROR: psPlaneDistortAlloc() did not return NULL.\n");
        testStatus = false;
        myPD = NULL;
    }

    printf("-------------------------------------------------------------------\n");
    printf("Calling psPlaneDistortAlloc() with negative-order polynomial request.  Should generate error, return NULL.\n");
    myPD = psPlaneDistortAlloc(1.0, -1.0, 1.0, 1.0);
    if (myPD != NULL) {
        printf("TEST ERROR: psPlaneDistortAlloc() did not return NULL.\n");
        testStatus = false;
        myPD = NULL;
    }

    printf("-------------------------------------------------------------------\n");
    printf("Calling psPlaneDistortAlloc() with negative-order polynomial request.  Should generate error, return NULL.\n");
    myPD = psPlaneDistortAlloc(1.0, 1.0, -1.0, 1.0);
    if (myPD != NULL) {
        printf("TEST ERROR: psPlaneDistortAlloc() did not return NULL.\n");
        testStatus = false;
        myPD = NULL;
    }

    printf("-------------------------------------------------------------------\n");
    printf("Calling psPlaneDistortAlloc() with negative-order polynomial request.  Should generate error, return NULL.\n");
    myPD = psPlaneDistortAlloc(1.0, 1.0, 1.0, -1.0);
    if (myPD != NULL) {
        printf("TEST ERROR: psPlaneDistortAlloc() did not return NULL.\n");
        testStatus = false;
        myPD = NULL;
    }
    printf("-------------------------------------------------------------------\n");

    return(testStatus);
}

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

    pt.x = psDPolynomial2DAlloc(2, 2, PS_POLYNOMIAL_ORD);
    pt.y = psDPolynomial2DAlloc(2, 2, PS_POLYNOMIAL_ORD);
    pt.x->coeff[1][0] = 1.0;
    pt.y->coeff[0][1] = 1.0;
    for (i=0;i<N;i++) {
        in.x = (float) i;
        in.y = (float) (i + 5.0);
        in.xErr = 0.0;
        in.yErr = 0.0;
        psPlaneTransformApply(&out, &pt, &in);

        if (FLT_EPSILON < fabs(out.x - in.x)) {
            printf("ERROR: out.x is %f, should be %f\n", out.x, in.x);
            testStatus = 3;
        }

        if (FLT_EPSILON < fabs(out.y - in.y)) {
            printf("ERROR: out.y is %f, should be %f\n", out.y, in.y);
            testStatus = 4;
        }
        //printf("psPlaneTransformApply(): (%f, %f) -> (%f, %f)\n", in.x, in.y, out.x, out.y);
    }

    printf("-------------------------------------------------------------------\n");
    printf("Calling psPlaneTransformApply() with NULL psPlaneTransform.  Should generate error, return NULL.\n");
    rc = psPlaneTransformApply(NULL, NULL, &in);
    if (rc != NULL) {
        printf("TEST ERROR: psPlaneTransformApply() did not return NULL.\n");
        testStatus = false;
        rc = NULL;
    }

    printf("-------------------------------------------------------------------\n");
    printf("Calling psPlaneTransformApply() with NULL psPlaneTransform->x.  Should generate error, return NULL.\n");
    tmp2DPoly = pt.x;
    pt.x = NULL;
    rc = psPlaneTransformApply(NULL, &pt, &in);
    if (rc != NULL) {
        printf("TEST ERROR: psPlaneTransformApply() did not return NULL.\n");
        testStatus = false;
        rc = NULL;
    }
    pt.x = tmp2DPoly;

    printf("-------------------------------------------------------------------\n");
    printf("Calling psPlaneTransformApply() with NULL psPlaneTransform->y.  Should generate error, return NULL.\n");
    tmp2DPoly = pt.y;
    pt.y = NULL;
    rc = psPlaneTransformApply(NULL, &pt, &in);
    if (rc != NULL) {
        printf("TEST ERROR: psPlaneTransformApply() did not return NULL.\n");
        testStatus = false;
        rc = NULL;
    }
    pt.y = tmp2DPoly;

    printf("-------------------------------------------------------------------\n");
    printf("Calling psPlaneTransformApply() with NULL input coords.  Should generate error, return NULL.\n");
    rc = psPlaneTransformApply(NULL, &pt, NULL);
    if (rc != NULL) {
        printf("TEST ERROR: psPlaneTransformApply() did not return NULL.\n");
        testStatus = false;
        rc = NULL;
    }

    printf("-------------------------------------------------------------------\n");
    psFree(pt.x);
    psFree(pt.y);
    return(testStatus);
}

#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 test3( void )
{
    psS32 i;
    psS32 testStatus = 0;
    psPlane in;
    psPlane out;
    psPlaneDistort pt;
    psPlane *rc = NULL;
    psDPolynomial4D *tmp4DPoly = NULL;

    pt.x = psDPolynomial4DAlloc(2, 2, 2, 2, PS_POLYNOMIAL_ORD);
    pt.y = psDPolynomial4DAlloc(2, 2, 2, 2, PS_POLYNOMIAL_ORD);
    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 (i=0;i<N;i++) {
        in.x = (float) i;
        in.y = (float) (i + 5.0);
        in.xErr = 0.0;
        in.yErr = 0.0;
        psPlaneDistortApply(&out, &pt, &in, COLOR, MAGNITUDE);

        if (FLT_EPSILON < fabs(out.x - COLOR - MAGNITUDE - in.x)) {
            printf("ERROR: out.x is %f, should be %f\n", out.x, in.x + COLOR + MAGNITUDE);
            testStatus = 3;
        }

        if (FLT_EPSILON < fabs(out.y - COLOR - MAGNITUDE - in.y)) {
            printf("ERROR: out.y is %f, should be %f\n", out.y, in.y + COLOR + MAGNITUDE);
            testStatus = 3;
        }
        //printf("psPlaneTransformApply(): (%f, %f) -> (%f, %f)\n", in.x, in.y, out.x, out.y);
    }

    printf("-------------------------------------------------------------------\n");
    printf("Calling psPlaneDistortApply() with NULL psPlaneDistort.  Should generate error, return NULL.\n");
    rc = psPlaneDistortApply(NULL, NULL, &in, COLOR, MAGNITUDE);
    if (rc != NULL) {
        printf("TEST ERROR: psPlaneDistortApply() did not return NULL.\n");
        testStatus = false;
        rc = NULL;
    }

    printf("-------------------------------------------------------------------\n");
    printf("Calling psPlaneDistortApply() with NULL psPlaneDistort->x.  Should generate error, return NULL.\n");
    tmp4DPoly = pt.x;
    pt.x = NULL;
    rc = psPlaneDistortApply(NULL, &pt, &in, COLOR, MAGNITUDE);
    if (rc != NULL) {
        printf("TEST ERROR: psPlaneDistortApply() did not return NULL.\n");
        testStatus = false;
        rc = NULL;
    }
    pt.x = tmp4DPoly;

    printf("-------------------------------------------------------------------\n");
    printf("Calling psPlaneDistortApply() with NULL psPlaneDistort->y.  Should generate error, return NULL.\n");
    tmp4DPoly = pt.y;
    pt.y = NULL;
    rc = psPlaneDistortApply(NULL, &pt, &in, COLOR, MAGNITUDE);
    if (rc != NULL) {
        printf("TEST ERROR: psPlaneDistortApply() did not return NULL.\n");
        testStatus = false;
        rc = NULL;
    }
    pt.y = tmp4DPoly;

    printf("-------------------------------------------------------------------\n");
    printf("Calling psPlaneDistortApply() with NULL input coords.  Should generate error, return NULL.\n");
    rc = psPlaneDistortApply(NULL, &pt, NULL, COLOR, MAGNITUDE);
    if (rc != NULL) {
        printf("TEST ERROR: psPlaneDistortApply() did not return NULL.\n");
        testStatus = false;
        rc = NULL;
    }



    printf("-------------------------------------------------------------------\n");
    psFree(pt.x);
    psFree(pt.y);
    return(testStatus);
}

#define DEG_INC 30.0
// We do a simple identity transformation on a few RA, DEC pairs.
psS32 test4( void )
{
    psS32 testStatus = 0;
    psSphere in;
    psSphere out;
    psSphere *rc = NULL;
    psSphereTransform *myST = psSphereTransformAlloc(0.0, 0.0, 0.0);

    in.r = DEG_TO_RAD(45.0);
    in.d = DEG_TO_RAD(30.0);
    in.rErr = 0.0;
    in.dErr = 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;

            psSphereTransformApply(&out, myST, &in);

            if (FLT_EPSILON < fabs(out.r - in.r)) {
                printf("ERROR: out.r is %f, should be %f\n", out.r, in.r);
                testStatus = 4;
            }

            if (FLT_EPSILON < fabs(out.d - in.d)) {
                printf("ERROR: out.d is %f, should be %f\n", out.d, in.d);
                testStatus = 5;
            }
            // printf("psSphereTransformApply (%f, %f) -> (%f, %f)\n", out.r, out.d, in.r, in.d);
        }
    }

    printf("-------------------------------------------------------------------\n");
    printf("Calling psSphereTransformApply() with NULL psSphereTransform.  Should generate error, return NULL.\n");
    rc = psSphereTransformApply(NULL, NULL, &in);
    if (rc != NULL) {
        printf("TEST ERROR: psSphereTransformApply() did not return NULL.\n");
        testStatus = false;
        rc = NULL;
    }

    printf("-------------------------------------------------------------------\n");
    printf("Calling psSphereTransformApply() with NULL input coords.  Should generate error, return NULL.\n");
    rc = psSphereTransformApply(NULL, myST, NULL);
    if (rc != NULL) {
        printf("TEST ERROR: psSphereTransformApply() did not return NULL.\n");
        testStatus = false;
        rc = NULL;
    }

    printf("-------------------------------------------------------------------\n");
    psFree(myST);
    return(testStatus);
}

/******************************************************************************
test4b(): 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 test4b( 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);
}

psS32 test5( void )
{
    psS32 testStatus = 0;

    // XXX: This test code is simply a copy of the original source code.
    psTime* now = psTimeGetTime(PS_TIME_UTC);
    struct tm *tm_time = psTimeToTM(now);
    psF64 MJD = psTimeToMJD(now);
    struct tm *tmTime = psTimeToTM(now);
    psF64 year = (psF64)(1900 - MJD + tmTime->tm_year);
    psF64 T = year / 100.0;
    psF64 expectedAlphaP = 0.0;
    psF64 expectedDeltaP = DEG_TO_RAD(23.0) +
                           SEC_TO_RAD(27.0) +
                           MIN_TO_RAD(8.0) -
                           (SEC_TO_RAD(46.845) * T) -
                           (SEC_TO_RAD(0.0059) * T * T) +
                           (SEC_TO_RAD(0.00181) * T * T * T);
    expectedDeltaP*= -1.0;
    psF64 expectedPhiP = 0.0;
    psFree(tmTime);

    psSphereTransform *myST = psSphereTransformICRSToEcliptic(now);

    if (FLT_EPSILON < fabs(sin(expectedDeltaP) - myST->sinDeltaP)) {
        printf("ERROR: myST->sinDeltaP is %f, should be %f\n", myST->sinDeltaP, sin(expectedDeltaP));
        testStatus = 1;
    }

    if (FLT_EPSILON < fabs(cos(expectedDeltaP) - myST->cosDeltaP)) {
        printf("ERROR: myST->cosDeltaP is %f, should be %f\n", myST->cosDeltaP, cos(expectedDeltaP));
        testStatus = 1;
    }

    if (FLT_EPSILON < fabs(expectedAlphaP - myST->alphaP)) {
        printf("ERROR: myST->alphaP is %f, should be %f\n", myST->alphaP, expectedAlphaP);
        testStatus = 3;
    }

    if (FLT_EPSILON < fabs(expectedPhiP - myST->phiP)) {
        printf("ERROR: myST->phiP is %f, should be %f\n", myST->phiP, expectedPhiP);
        testStatus = 4;
    }

    psFree(myST);
    psFree(tm_time);
    psFree(now);

    printf("-------------------------------------------------------------------\n");
    printf("Calling psSphereTransformICRSToEcliptic() with NULL input.  Should generate error, return NULL.\n");
    myST = psSphereTransformICRSToEcliptic(NULL);
    if (myST != NULL) {
        printf("TEST ERROR: psSphereTransformApply() did not return NULL.\n");
        testStatus = false;
        myST = NULL;
    }

    return(testStatus);
}


psS32 test6( void )
{
    psS32 testStatus = 0;
    psTime* now = psTimeGetTime(PS_TIME_UTC);
    struct tm *tm_time = psTimeToTM(now);
    psF64 MJD = psTimeToMJD(now);
    struct tm *tmTime = psTimeToTM(now);
    psF64 year = (psF64)(1900 - MJD + tmTime->tm_year);
    psF64 T = year / 100.0;
    psF64 expectedAlphaP = 0.0;
    psF64 expectedDeltaP = DEG_TO_RAD(23.0) +
                           SEC_TO_RAD(27.0) +
                           MIN_TO_RAD(8.0) -
                           (SEC_TO_RAD(46.845) * T) -
                           (SEC_TO_RAD(0.0059) * T * T) +
                           (SEC_TO_RAD(0.00181) * T * T * T);
    psF64 expectedPhiP = 0.0;
    psSphereTransform *myST = psSphereTransformEclipticToICRS(now);

    if (FLT_EPSILON < fabs(sin(expectedDeltaP) - myST->sinDeltaP)) {
        printf("ERROR: myST->sinDeltaP is %f, should be %f\n", myST->sinDeltaP, sin(expectedDeltaP));
        testStatus = 1;
    }

    if (FLT_EPSILON < fabs(cos(expectedDeltaP) - myST->cosDeltaP)) {
        printf("ERROR: myST->cosDeltaP is %f, should be %f\n", myST->cosDeltaP, cos(expectedDeltaP));
        testStatus = 1;
    }

    if (FLT_EPSILON < fabs(expectedAlphaP - myST->alphaP)) {
        printf("ERROR: myST->alphaP is %f, should be %f\n", myST->alphaP, expectedAlphaP);
        testStatus = 3;
    }

    if (FLT_EPSILON < fabs(expectedPhiP - myST->phiP)) {
        printf("ERROR: myST->phiP is %f, should be %f\n", myST->phiP, expectedPhiP);
        testStatus = 4;
    }

    psFree(myST);
    psFree(tm_time);
    psFree(tmTime);
    psFree(now);

    printf("-------------------------------------------------------------------\n");
    printf("Calling psSphereTransformEclipticToICRS() with NULL input.  Should generate error, return NULL.\n");
    myST = psSphereTransformEclipticToICRS(NULL);
    if (myST != NULL) {
        printf("TEST ERROR: psSphereTransformApply() did not return NULL.\n");
        testStatus = false;
        myST = NULL;
    }

    return(testStatus);
}

psS32 test7( void )
{
    psS32 testStatus = 0;
    // XXX: This test code is simply a copy of the original source code.
    psF64 expectedAlphaP = DEG_TO_RAD(32.93192);
    psF64 expectedDeltaP = DEG_TO_RAD(-62.87175);
    psF64 expectedPhiP = DEG_TO_RAD(282.85948);
    psSphereTransform *myST = psSphereTransformICRSToGalatic();

    if (FLT_EPSILON < fabs(sin(expectedDeltaP) - myST->sinDeltaP)) {
        printf("ERROR: myST->sinDeltaP is %f, should be %f\n", myST->sinDeltaP, sin(expectedDeltaP));
        testStatus = 1;
    }

    if (FLT_EPSILON < fabs(cos(expectedDeltaP) - myST->cosDeltaP)) {
        printf("ERROR: myST->cosDeltaP is %f, should be %f\n", myST->cosDeltaP, cos(expectedDeltaP));
        testStatus = 1;
    }

    if (FLT_EPSILON < fabs(expectedAlphaP - myST->alphaP)) {
        printf("ERROR: myST->alphaP is %f, should be %f\n", myST->alphaP, expectedAlphaP);
        testStatus = 3;
    }

    if (FLT_EPSILON < fabs(expectedPhiP - myST->phiP)) {
        printf("ERROR: myST->phiP is %f, should be %f\n", myST->phiP, expectedPhiP);
        testStatus = 4;
    }

    psFree(myST);
    return(testStatus);
}

psS32 test8( void )
{
    psS32 testStatus = 0;
    // XXX: This test code is simply a copy of the original source code.
    psF64 expectedAlphaP = DEG_TO_RAD(282.85948);
    psF64 expectedDeltaP = DEG_TO_RAD(62.87175);
    psF64 expectedPhiP = DEG_TO_RAD(32.93192);
    psSphereTransform *myST = psSphereTransformGalaticToICRS();

    if (FLT_EPSILON < fabs(sin(expectedDeltaP) - myST->sinDeltaP)) {
        printf("ERROR: myST->sinDeltaP is %f, should be %f\n", myST->sinDeltaP, sin(expectedDeltaP));
        testStatus = 1;
    }

    if (FLT_EPSILON < fabs(cos(expectedDeltaP) - myST->cosDeltaP)) {
        printf("ERROR: myST->cosDeltaP is %f, should be %f\n", myST->cosDeltaP, cos(expectedDeltaP));
        testStatus = 1;
    }

    if (FLT_EPSILON < fabs(expectedAlphaP - myST->alphaP)) {
        printf("ERROR: myST->alphaP is %f, should be %f\n", myST->alphaP, expectedAlphaP);
        testStatus = 3;
    }

    if (FLT_EPSILON < fabs(expectedPhiP - myST->phiP)) {
        printf("ERROR: myST->phiP is %f, should be %f\n", myST->phiP, expectedPhiP);
        testStatus = 4;
    }

    psFree(myST);
    return(testStatus);
}

psS32 test20( void )
{
    psS32 testStatus = 0;
    psSphere in;
    double expectX;
    double expectY;
    psPlane *out;
    psProjection myProjection;
    myProjection.R = DEG_TO_RAD(20.0);
    myProjection.D = DEG_TO_RAD(10.0);
    myProjection.Xs = 1.0;
    myProjection.Ys = 1.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;

            /******************************************************************
             Tangent Plane Projection
             *****************************************************************/
            myProjection.type = PS_PROJ_TAN;
            out = psProject(&in, &myProjection);
            expectX = cos(in.r) * sin(in.d) / sin(in.r);
            expectY = -cos(in.r) * cos(in.d) / sin(in.r);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->x, expectX, 1);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->y, expectY, 2);
            psFree(out);

            /******************************************************************
             PS_PROJ_SIN Projection
             *****************************************************************/
            myProjection.type = PS_PROJ_SIN;
            out = psProject(&in, &myProjection);
            expectX = cos(in.d) * sin(in.d);
            expectY = -cos(in.d) * cos(in.d);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->x, expectX, 3);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->y, expectY, 4);
            psFree(out);

            /******************************************************************
             PS_PROJ_CAR Projection
             *****************************************************************/
            myProjection.type = PS_PROJ_CAR;
            out = psProject(&in, &myProjection);

            expectX = in.d;
            expectY = in.r;
            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->x, expectX, 5);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->y, expectY, 6);
            psFree(out);

            /******************************************************************
             PS_PROJ_MER Projection
             *****************************************************************/
            myProjection.type = PS_PROJ_MER;
            out = psProject(&in, &myProjection);

            expectX = in.d;
            expectY = log(tan(DEG_TO_RAD(45.0) + (0.5 * in.r)));

            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->x, expectX, 7);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->y, expectY, 8);
            psFree(out);

            /******************************************************************
             PS_PROJ_AIT Projection
             *****************************************************************/
            myProjection.type = PS_PROJ_AIT;
            out = psProject(&in, &myProjection);


            psF64 tmpF64  = PS_SQRT_F32(0.5 * (1.0 + (cos(in.r) * cos(0.5 * in.d))));
            tmpF64 = 1.0 / tmpF64;
            expectX = 2.0 * tmpF64 * cos(in.r) * sin(0.5 * in.d);
            expectY = sin(in.r) * tmpF64;

            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->x, expectX, 20);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->y, expectY, 21);
            psFree(out);


        }
    }

    printf("-------------------------------------------------------------------\n");
    printf("Calling psProject() with NULL input coords.  Should generate error, return NULL.\n");
    out = psProject(NULL, &myProjection);
    if (out != NULL) {
        printf("TEST ERROR: psSphereTransformEclipticToICRS() did not return NULL.\n");
        testStatus = false;
        out = NULL;
    }

    printf("-------------------------------------------------------------------\n");
    printf("Calling psProject() with NULL psProjection.  Should generate error, return NULL.\n");
    out = psProject(&in, NULL);
    if (out != NULL) {
        printf("TEST ERROR: psProject() did not return NULL.\n");
        testStatus = false;
        out = NULL;
    }

    printf("-------------------------------------------------------------------\n");

    return(testStatus);
}

#define SPACE_INC 10.0
psS32 test21( void )
{
    psS32 testStatus = 0;
    psPlane in;
    double R;
    double expectD;
    double expectR;
    psSphere *out;
    psProjection myProjection;
    myProjection.R = DEG_TO_RAD(20.0);
    myProjection.D = DEG_TO_RAD(10.0);
    myProjection.Xs = 1.0;
    myProjection.Ys = 1.0;

    in.xErr = 0.0;
    in.yErr = 0.0;
    for (in.x=0.0; in.x<100.0; in.x+=SPACE_INC) {
        for (in.y=0.0;in.y<100.0;in.y+=SPACE_INC) {

            /******************************************************************
             Tangent Plane Projection
             *****************************************************************/
            myProjection.type = PS_PROJ_TAN;
            out = psDeproject(&in, &myProjection);

            R = PS_SQRT_F32((in.x * in.x) + (in.y * in.y));
            expectD = atan(1.0 / R);
            expectR = atan2(-in.y, in.x);

            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->r, expectR, 1);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->d, expectD, 2);
            psFree(out);

            /******************************************************************
             PS_PROJ_SIN Projection
             *****************************************************************/
            myProjection.type = PS_PROJ_SIN;
            out = psDeproject(&in, &myProjection);

            R = sqrt((in.x * in.x) + (in.y * in.y));
            expectD = p_psArg(-in.y, in.x);
            expectR = acos((R * M_PI) / 180.0);

            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->r, expectR, 3);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->d, expectD, 4);
            psFree(out);

            /******************************************************************
             PS_PROJ_CAR Projection
             *****************************************************************/
            myProjection.type = PS_PROJ_CAR;
            out = psDeproject(&in, &myProjection);

            expectD = in.x;
            expectR = in.y;

            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->r, expectR, 3);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->d, expectD, 4);
            psFree(out);

            /******************************************************************
             PS_PROJ_MER Projection
             *****************************************************************/
            myProjection.type = PS_PROJ_MER;
            out = psDeproject(&in, &myProjection);

            expectD = in.x;
            expectR = (2.0 * atan(exp((in.y * M_PI / 180.0)))) - 180.0;

            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->r, expectR, 3);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->d, expectD, 4);
            psFree(out);


            /******************************************************************
             PS_PROJ_AIT Projection
             *****************************************************************/
            myProjection.type = PS_PROJ_AIT;
            out = psDeproject(&in, &myProjection);

            float chu1 = (in.x * M_PI) / 720.0;
            chu1 *= chu1;
            float chu2 = (in.y * M_PI) / 360.0;
            chu2 *= chu2;
            float chu = sqrt(1.0 - chu1 - chu2);
            expectD = 2.0 * p_psArg((2.0 * chu * chu) - 1.0, (in.x * chu * M_PI) / 360.0);
            expectR = asin((in.y * chu * M_PI) / 180.0);

            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->r, expectR, 3);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(out->d, expectD, 4);
            psFree(out);
        }
    }

    printf("-------------------------------------------------------------------\n");
    printf("Calling psDeproject() with NULL input coords.  Should generate error, return NULL.\n");
    out = psDeproject(NULL, &myProjection);
    if (out != NULL) {
        printf("TEST ERROR: psDeproject() did not return NULL.\n");
        testStatus = false;
        out = NULL;
    }

    printf("-------------------------------------------------------------------\n");
    printf("Calling psDeproject() with NULL psProjection.  Should generate error, return NULL.\n");
    out = psDeproject(&in, NULL);
    if (out != NULL) {
        printf("TEST ERROR: psDeproject() did not return NULL.\n");
        testStatus = false;
        out = NULL;
    }

    printf("-------------------------------------------------------------------\n");

    return(testStatus);
}

psS32 test40( void )
{
    psS32 testStatus = 0;
    psSphere position1;
    psSphere position2;
    psSphere *offset = NULL;

    position1.r = DEG_TO_RAD(90.0);
    position1.d = DEG_TO_RAD(45.0);
    position1.rErr = 0.0;
    position1.dErr = 0.0;

    for (float r=0.0;r<180.0;r+=DEG_INC) {
        for (float d=0.0;d<90.0;d+=DEG_INC) {
            position2.r = DEG_TO_RAD(r);
            position2.d = DEG_TO_RAD(d);
            position2.rErr = 0.0;
            position2.dErr = 0.0;

            offset = psSphereGetOffset(&position1, &position2,
                                       PS_SPHERICAL, PS_RADIAN);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->r, (position2.r - position1.r), 1);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->d, (position2.d - position1.d), 1);
            psFree(offset);

            offset = psSphereGetOffset(&position1, &position2,
                                       PS_SPHERICAL, PS_DEGREE);
            offset->r = DEG_TO_RAD(offset->r);
            offset->d = DEG_TO_RAD(offset->d);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->r, (position2.r - position1.r), 2);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->d, (position2.d - position1.d), 3);
            psFree(offset);

            offset = psSphereGetOffset(&position1, &position2,
                                       PS_SPHERICAL, PS_ARCMIN);
            offset->r = MIN_TO_RAD(offset->r);
            offset->d = MIN_TO_RAD(offset->d);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->r, (position2.r - position1.r), 2);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->d, (position2.d - position1.d), 3);
            psFree(offset);

            offset = psSphereGetOffset(&position1, &position2,
                                       PS_SPHERICAL, PS_ARCSEC);
            offset->r = SEC_TO_RAD(offset->r);
            offset->d = SEC_TO_RAD(offset->d);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->r, (position2.r - position1.r), 2);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->d, (position2.d - position1.d), 3);
            psFree(offset);

            //HEY            /* XXX: This code does not work correctly.
            //            printf("--------------- (%f, %f) to (%f, %f) is (%f, %f) ---------------\n",
            //                   position1.r, position1.d, position2.r, position2.d,
            //                   offset->r, offset->d);
            //            offset = psSphereGetOffset(&position1, &position2, PS_LINEAR, 0);
            //            PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->r, (position2.r - position1.r), 1);
            //            PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->d, (position2.d - position1.d), 1);
            //            psFree(offset);
        }
    }

    printf("-------------------------------------------------------------------\n");
    printf("Calling psSphereGetOffset() with NULL position1.  Should generate error, return NULL.\n");
    offset = psSphereGetOffset(NULL, &position2, PS_LINEAR, 0);
    if (offset != NULL) {
        printf("TEST ERROR: psSphereGetOffset() did not return NULL.\n");
        testStatus = false;
        offset = NULL;
    }

    printf("-------------------------------------------------------------------\n");
    printf("Calling psSphereGetOffset() with NULL position2.  Should generate error, return NULL.\n");
    offset = psSphereGetOffset(&position1, NULL, PS_LINEAR, 0);
    if (offset != NULL) {
        printf("TEST ERROR: psSphereGetOffset() did not return NULL.\n");
        testStatus = false;
        offset = NULL;
    }

    printf("-------------------------------------------------------------------\n");
    printf("Calling psSphereGetOffset() with bogus offset mode.  Should generate error, return NULL.\n");
    offset = psSphereGetOffset(&position1, &position2, 0x54321, 0);
    if (offset != NULL) {
        printf("TEST ERROR: psSphereGetOffset() did not return NULL.\n");
        testStatus = false;
        offset = NULL;
    }

    printf("-------------------------------------------------------------------\n");
    printf("Calling psSphereGetOffset() with bogus offset units.  Should generate error, return NULL.\n");
    offset = psSphereGetOffset(&position1, &position2, PS_SPHERICAL, 0x54321);
    if (offset != NULL) {
        printf("TEST ERROR: psSphereGetOffset() did not return NULL.\n");
        testStatus = false;
        offset = NULL;
    }

    printf("-------------------------------------------------------------------\n");
    return(testStatus);
}

psS32 test41( void )
{
    psS32 testStatus = 0;
    psSphere position1;
    psSphere *position2;
    psSphere offset;
    psSphere tmpOffset;

    position1.r = DEG_TO_RAD(90.0);
    position1.d = DEG_TO_RAD(45.0);
    position1.rErr = 0.0;
    position1.dErr = 0.0;

    for (float r=0.0;r<180.0;r+=DEG_INC) {
        for (float d=0.0;d<90.0;d+=DEG_INC) {
            offset.r = DEG_TO_RAD(r);
            offset.d = DEG_TO_RAD(d);
            offset.rErr = 0.0;
            offset.dErr = 0.0;

            position2 = psSphereSetOffset(&position1, &offset,
                                          PS_SPHERICAL, PS_RADIAN);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->r, (position1.r + offset.r), 1);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->d, (position1.d + offset.d), 2);
            psFree(position2);

            tmpOffset.r = RAD_TO_DEG(offset.r);
            tmpOffset.d = RAD_TO_DEG(offset.d);
            tmpOffset.rErr = 0.0;
            tmpOffset.dErr = 0.0;
            position2 = psSphereSetOffset(&position1, &tmpOffset,
                                          PS_SPHERICAL, PS_DEGREE);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->r, (position1.r + offset.r), 3);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->d, (position1.d + offset.d), 4);
            psFree(position2);

            tmpOffset.r = RAD_TO_MIN(offset.r);
            tmpOffset.d = RAD_TO_MIN(offset.d);
            tmpOffset.rErr = 0.0;
            tmpOffset.dErr = 0.0;
            position2 = psSphereSetOffset(&position1, &tmpOffset,
                                          PS_SPHERICAL, PS_ARCMIN);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->r, (position1.r + offset.r), 5);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->d, (position1.d + offset.d), 6);
            psFree(position2);

            tmpOffset.r = RAD_TO_SEC(offset.r);
            tmpOffset.d = RAD_TO_SEC(offset.d);
            tmpOffset.rErr = 0.0;
            tmpOffset.dErr = 0.0;
            position2 = psSphereSetOffset(&position1, &tmpOffset,
                                          PS_SPHERICAL, PS_ARCSEC);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->r, (position1.r + offset.r), 7);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->d, (position1.d + offset.d), 8);
            psFree(position2);

            /* XXX: This code does not work correctly.
                        position2 = psSphereSetOffset(&position1, &offset,
                                                       PS_LINEAR, 0);
                        printf("--------------- (%f, %f) and (%f, %f) is (%f, %f) ---------------\n",
                                position1.r, position1.d, offset.r, offset.d, 
                                position2->r, position2->d);
                        PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->r, (position1.r + offset.r), 1);
                        PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->d, (position1.d + offset.d), 1);
                        psFree(position2);
            */
        }
    }

    printf("-------------------------------------------------------------------\n");
    printf("Calling psSphereSetOffset() with NULL position1.  Should generate error, return NULL.\n");
    position2 = psSphereSetOffset(NULL, &offset, PS_LINEAR, 0);
    if (position2 != NULL) {
        printf("TEST ERROR: psSphereSetOffset() did not return NULL.\n");
        testStatus = false;
        position2 = NULL;
    }

    printf("-------------------------------------------------------------------\n");
    printf("Calling psSphereSetOffset() with NULL offset.  Should generate error, return NULL.\n");
    position2 = psSphereSetOffset(&position1, NULL, PS_LINEAR, 0);
    if (position2 != NULL) {
        printf("TEST ERROR: psSphereSetOffset() did not return NULL.\n");
        testStatus = false;
        position2 = NULL;
    }

    printf("-------------------------------------------------------------------\n");
    printf("Calling psSphereSetOffset() with bogus offset mode.  Should generate error, return NULL.\n");
    position2 = psSphereSetOffset(&position1, &offset, 0x54321, 0);
    if (position2 != NULL) {
        printf("TEST ERROR: psSphereSetOffset() did not return NULL.\n");
        testStatus = false;
        position2 = NULL;
    }

    printf("-------------------------------------------------------------------\n");
    printf("Calling psSphereSetOffset() with bogus offset unit.  Should generate error, return NULL.\n");
    position2 = psSphereSetOffset(&position1, &offset, PS_SPHERICAL, 0x54321);
    if (position2 != NULL) {
        printf("TEST ERROR: psSphereSetOffset() did not return NULL.\n");
        testStatus = false;
        position2 = NULL;
    }

    printf("-------------------------------------------------------------------\n");
    return(testStatus);
}

/******************************************************************************
test43(): Not really a test.  We simply project many points and print the
results.
 *****************************************************************************/
psS32 test43( void )
{
    psS32 testStatus = 0;
    psSphere projCenter;
    psSphere startCoords;
    psPlane *planeCoords;

    projCenter.r = DEG_TO_RAD(33.0);
    projCenter.d = DEG_TO_RAD(22.0);
    projCenter.r = 0.0;
    projCenter.d = 0.0;
    // Create the psProjection plane:
    psProjection myProj;
    myProj.R = projCenter.r;
    myProj.D = projCenter.d;
    myProj.Xs = 1.0;
    myProj.Ys = 1.0;
    myProj.type = PS_PROJ_TAN;

    printf("    projCenter is (%f, %f)\n", projCenter.r, projCenter.d);
    for (float r=-90.0;r<90.0;r+=DEG_INC/5.0) {
        for (float d=-90.0;d<90.0;d+=DEG_INC/5.0) {
            printf("----------------------------------------------------------\n");
            printf("Angles (R, D) are (%.2f, %.2f)\n", r, d);
            startCoords.r = DEG_TO_RAD(r);
            startCoords.d = DEG_TO_RAD(d);

            planeCoords = psProject(&startCoords, &myProj);
            printf("    startCoords (R, D) is (%f, %f)\n", startCoords.r, startCoords.d);
            printf("        plane Coords (X, Y) is (%f, %f)\n", planeCoords->x, planeCoords->y);
            psFree(planeCoords);
        }
    }
    printf("----------------------------------------------------------\n");
    return(testStatus);
}

bool testNonEqualF32(psF32 num1, psF32 num2)
{
    if ((fabs(num1) <= FLT_EPSILON) && (fabs(num2) <= FLT_EPSILON)) {
        return(false);
    }

    if (fabs((num2 - num1) / num1) <= ERROR_PERCENT) {
        return(false);
    }
    return(true);
}

#define NUM_DEGREES 5
/******************************************************************************
test42(): This test attempts to verify the psProject() and psDeproject()
functions.  We create an arbitrary psProjection.  We then project several
spherical coordinates onto the plane and then immediately deproject them; we
ensure that the resulting spherical coordinates are identical to the start
coordinates.
 
Currently, this test works only for angles near the first quadrant.
 *****************************************************************************/
psS32 test42( void )
{
    psS32 testStatus = 0;
    psSphere projCenter;
    psSphere startCoords;
    psSphere *endCoords;
    psPlane *planeCoords;

    #define PROG_CENTER_R_DEG 15.0
    #define PROG_CENTER_D_DEG 20.0

    projCenter.r = DEG_TO_RAD(PROG_CENTER_R_DEG);
    projCenter.d = DEG_TO_RAD(PROG_CENTER_D_DEG);
    // Create the psProjection plane:
    psProjection myProj;
    myProj.R = projCenter.r;
    myProj.D = projCenter.d;
    myProj.Xs = 1.0;
    myProj.Ys = 1.0;
    myProj.type = PS_PROJ_TAN;

    #define RADIUS_DEG 30.0
    #define NUM_ANGLES 30
    #define MY_DEG_INC ((2.0 * RADIUS_DEG) / NUM_ANGLES)

    for (float r= (PROG_CENTER_R_DEG - RADIUS_DEG);
            r<=(PROG_CENTER_R_DEG + RADIUS_DEG);
            r+=MY_DEG_INC) {
        printf("------------------------ r is %.1f ------------------------\n", r);
        for (float d= (PROG_CENTER_D_DEG - RADIUS_DEG);
                d<=(PROG_CENTER_D_DEG + RADIUS_DEG);
                d+=MY_DEG_INC) {
            startCoords.r = DEG_TO_RAD(r);
            startCoords.d = DEG_TO_RAD(d);
            printf("----- startCoords (R, D) is (%f, %f) is deg (%.1f, %.1f)  -----\n", startCoords.r, startCoords.d, r, d);

            planeCoords = psProject(&startCoords, &myProj);
            endCoords = psDeproject(planeCoords, &myProj);

            if (testNonEqualF32(endCoords->r, startCoords.r) ||
                    testNonEqualF32(endCoords->d, startCoords.d)) {
                printf("ERROR: \n");
                printf("    startCoords (R, D) is (%f, %f)\n", startCoords.r, startCoords.d);
                //                printf("        plane Coords is (%f, %f)\n", planeCoords->x, planeCoords->y);
                printf("    endCoords (R, D) is (%f, %f)\n", endCoords->r, endCoords->d);
            }
            psFree(planeCoords);
            psFree(endCoords);
        }
    }
    return(testStatus);
}

// This code will ...

