/** @file  tst_psImageManip.c
*
*  @brief The code will ...
*
*
*  @author GLG, MHPCC
*
*  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
*  @date $Date: 2004-10-14 01:53:44 $
*
*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
*/

#include "psTest.h"
#include "pslib.h"

static int test1( void );
static int test2( void );
static int test3( void );
static int test4( void );
static int test5( void );
static int test6( void );
static int test7( void );
static int test8( void );
static int test20( void );
static int test21( void );
static int test40( void );
static int test41( void );

testDescription tests[] = {
                              {test1, 0000, "psSphereTransformAlloc()", 0, false},
                              {test2, 0000, "psPlaneTransformApply()", 0, false},
                              {test3, 0000, "psPlaneDistortApply()", 0, false},
                              {test4, 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},
                              {NULL}
                          };

#define ARCSEC_TO_RAD (M_PI / (180.0 * 60.0 * 60.0))
#define ARCMIN_TO_RAD (M_PI / (180.0 * 60.0))
#define DEG_TO_RAD (M_PI / 180.0)
#define RAD_TO_DEG (180.0 / M_PI)
#define RAD_TO_ARCMIN ((180.0 * 60.0) / M_PI)
#define RAD_TO_ARCSEC ((180.0 * 60.0 * 60.0) / M_PI)


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

    return ( ! runTestSuite( stderr, "psImage", 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\n", #ACTUAL, ACTUAL, EXPECT); \
    testStatus = TESTSTATUS; \
}
#define MY_TINY 0.0001
#define PS_COMPARE_EPSILON_THEN_PRINT_ERROR(ACTUAL, EXPECT, TESTSTATUS) \
if (MY_TINY < fabs(EXPECT - ACTUAL)) { \
    printf("ERROR: %s is %f, should be %f\n", #ACTUAL, ACTUAL, EXPECT); \
    testStatus = TESTSTATUS; \
}

#define PS_COT(X) (1.0 / atan(X))
float p_psArg(float x,
              float y)
{
    if (y!= 0.0) {
        return(atan(y/x));
    } else {
        return(0.5 * M_PI);
    }
    return (0.0);
}

#define NPLAT 1.0
#define X0_BIG 2.0
#define X0_SMALL 2.0

int test1( void )
{
    float tmpF32;
    int testStatus = 0;
    psSphereTransform *myST = psSphereTransformAlloc(NPLAT, X0_BIG, X0_SMALL);

    tmpF32 = sin(NPLAT);
    if (FLT_EPSILON < fabs(tmpF32 - myST->sinPhi)) {
        printf("ERROR: myST->sinPhi is %f, should be %f\n", myST->sinPhi, tmpF32);
        testStatus = 1;
    }
    tmpF32 = cos(NPLAT);
    if (FLT_EPSILON < fabs(tmpF32 - myST->cosPhi)) {
        printf("ERROR: myST->cosPhi is %f, should be %f\n", myST->cosPhi, tmpF32);
        testStatus = 2;
    }

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

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

    psFree(myST);
    return(testStatus);
}

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

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

    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.
int test3( void )
{
    int i;
    int testStatus = 0;
    psPlane in;
    psPlane out;
    psPlaneDistort pt;

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

    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.
int test4( void )
{
    int testStatus = 0;
    psSphere in;
    psSphere out;
    psSphereTransform *myST = psSphereTransformAlloc(0.0, 0.0, 0.0);

    in.r = 45.0 * DEG_TO_RAD;
    in.d = 30.0 * DEG_TO_RAD;
    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 = r * DEG_TO_RAD;
            in.d = d * DEG_TO_RAD;
            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);
        }
    }

    psFree(myST);
    return(testStatus);
}

int test5( void )
{
    int testStatus = 0;
    psTime* now = psTimeGetTime(PS_TIME_UTC);
    struct tm *tm_time = psTimeToTM(now);

    // XXX: This test code is simply a copy of the original source code.
    double year = (double)(1900 + tm_time->tm_year);
    double T = year / 100.0;
    double phi = -23.452294 + 0.013013 * T + 0.000001639 * T * T - 0.000000503 * T * T * T;
    psSphereTransform *myST = psSphereTransformICRStoEcliptic(*now);

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

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

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

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

    psFree(myST);
    psFree(tm_time);
    psFree(now);
    return(testStatus);
}


int test6( void )
{
    int testStatus = 0;
    psTime* now = psTimeGetTime(PS_TIME_UTC);
    struct tm *tm_time = psTimeToTM(now);

    // XXX: This test code is simply a copy of the original source code.
    double year = (double)(1900 + tm_time->tm_year);
    double T = year / 100.0;
    double phi = +23.452294 - 0.013013 * T - 0.000001639 * T * T + 0.000000503 * T * T * T;
    psSphereTransform *myST = psSphereTransformEcliptictoICRS(*now);

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

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

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

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

    psFree(myST);
    psFree(tm_time);
    psFree(now);
    return(testStatus);
}


int test7( void )
{
    int testStatus = 0;
    // XXX: This test code is simply a copy of the original source code.
    double phi = 62.6;
    double Xo = 282.25;
    double xo = 33.0;
    psSphereTransform *myST = psSphereTransformICRStoGalatic();

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

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

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

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

    psFree(myST);
    return(testStatus);
}

int test8( void )
{
    int testStatus = 0;
    // XXX: This test code is simply a copy of the original source code.
    double phi = -62.6;
    double Xo = 33.0;
    double xo = 282.25;
    psSphereTransform *myST = psSphereTransformGalatictoICRS();

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

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

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

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

    psFree(myST);
    return(testStatus);
}

int test20( void )
{
    int testStatus = 0;
    psSphere in;
    double R;
    double expectX;
    double expectY;
    psPlane *out;
    psProjection myProjection;
    myProjection.R = 20.0 * DEG_TO_RAD;
    myProjection.D = 10.0 * DEG_TO_RAD;
    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 = r * DEG_TO_RAD;
            in.d = d * DEG_TO_RAD;
            in.rErr = 0.0;
            in.dErr = 0.0;

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

            R = PS_COT(in.r) * (180.0 / M_PI);
            expectX = R * sin(in.d);
            expectY = R * cos(in.d);
            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);

            R = cos(in.r) * (180.0 / M_PI);
            expectX = R * sin(in.d);
            expectY = R * 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(45.0 + (0.5 * in.r))) * 180.0 / M_PI;

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

            double alpha = 1.0 / ((180.0 / M_PI) * sqrt(1.0 + (cos(in.r) *
                                  cos(0.5 * in.d) * 0.5)));
            expectX = 2.0 * alpha * cos(in.r) * sin(0.5 * in.d);
            expectY = alpha * sin(in.d);

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


        }
    }
    return(testStatus);
}

#define SPACE_INC 10.0
int test21( void )
{
    int testStatus = 0;
    psPlane in;
    double R;
    double expectD;
    double expectR;
    psSphere *out;
    psProjection myProjection;
    myProjection.R = 20.0 * DEG_TO_RAD;
    myProjection.D = 10.0 * DEG_TO_RAD;
    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 = sqrt((in.x * in.x) + (in.y * in.y));
            expectD = p_psArg(-in.y, in.x);
            expectR = atan(180.0 / (R * M_PI));

            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);
        }
    }
    return(testStatus);
}

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

    position1.r = 90.0 * DEG_TO_RAD;
    position1.d = 45.0 * DEG_TO_RAD;
    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 = r * DEG_TO_RAD;
            position2.d = d * DEG_TO_RAD;
            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 = offset->r * DEG_TO_RAD;
            offset->d = offset->d * DEG_TO_RAD;
            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 = offset->r * ARCMIN_TO_RAD;
            offset->d = offset->d * ARCMIN_TO_RAD;
            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 = offset->r * ARCSEC_TO_RAD;
            offset->d = offset->d * ARCSEC_TO_RAD;
            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);

            /* XXX: This code does not work correctly.
                        offset = psSphereGetOffset(&position1, &position2,
                                                   PS_LINEAR, 0);
                        printf("--------------- (%f, %f) to (%f, %f) is (%f, %f) ---------------\n",
                                position1.r, position1.d, position2.r, position2.d,
                                offset->r, offset->d);
                        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);
            */
        }
    }

    return(testStatus);
}

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

    position1.r = 90.0 * DEG_TO_RAD;
    position1.d = 45.0 * DEG_TO_RAD;
    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 = r * DEG_TO_RAD;
            offset.d = d * DEG_TO_RAD;
            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), 1);
            psFree(position2);

            tmpOffset.r = offset.r * RAD_TO_DEG;
            tmpOffset.d = offset.d * RAD_TO_DEG;
            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), 1);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->d, (position1.d + offset.d), 1);
            psFree(position2);

            tmpOffset.r = offset.r * RAD_TO_ARCMIN;
            tmpOffset.d = offset.d * RAD_TO_ARCMIN;
            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), 1);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->d, (position1.d + offset.d), 1);
            psFree(position2);

            tmpOffset.r = offset.r * RAD_TO_ARCSEC;
            tmpOffset.d = offset.d * RAD_TO_ARCSEC;
            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), 1);
            PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->d, (position1.d + offset.d), 1);
            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);
            */
        }
    }

    return(testStatus);
}
