/** @file  tst_psCoord.c
*
*  @brief The code will ...
*
*
*  @author GLG, MHPCC
*
*  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
*  @date $Date: 2005-10-08 03:51:20 $
*
*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
*/

#include "psTest.h"
#include "pslib_strict.h"
static psS32 testSphereRotAlloc(void);
static psS32 testSphereRotApply1(void);
static psS32 testSphereRotApply2(void);
static psS32 testSphereRotApplyCelestial(void);
static psS32 testSphereRotPrecess(void);

testDescription tests[] = {
                              {testSphereRotAlloc, 819, "psSphereRotAlloc()", 0, false},
                              {testSphereRotApply1, 820, "psSphereRotApply()", 0, false},
                              {testSphereRotApply2, 820, "psSphereRotApply()", 0, false},
                              {testSphereRotApplyCelestial, 820, "psSphereRotApply()", 0, false},
                              //                              {testSphereRotApplyCelestial, 821, "psSphereRotICRSToEcliptic()", true},
                              //                              {testSphereRotApplyCelestial, 822, "psSphereRotEclipticToICRS()", 0, true},
                              //                              {testSphereRotApplyCelestial, 824, "psSphereRotICRSToGalactic()", 0, true},
                              //                              {testSphereRotApplyCelestial, 823, "psSphereRotGalacticToICRS()", 0, true},
                              {testSphereRotPrecess, 825, "psSphereRotPrecess()", 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 0.5 * M_PI
#define DELTA_P 1.0 * M_PI
#define PHI_P 1.5 * M_PI

psS32 testSphereRotAlloc( void )
{
    // Allocate data structure
    psSphereRot* myST = psSphereRotAlloc(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
    double vx = cos(DELTA_P)*cos(ALPHA_P);
    double vy = cos(DELTA_P)*sin(ALPHA_P);
    double vz = sin(DELTA_P);
    double q0 = vx * sin(PHI_P/2.0);
    double q1 = vy * sin(PHI_P/2.0);
    double q2 = vz * sin(PHI_P/2.0);
    double q3 = cos(PHI_P/2.0);
    double len = sqrt(q0*q0 + q1*q1 + q2*q2 + q3*q3);
    printf("\nvz = %lf\n", vz);
    q0 = q0/len;
    q1 = q1/len;
    q2 = q2/len;
    q3 = q3/len;

    if (FLT_EPSILON < fabs(q0 - myST->q0)) {
        psError(PS_ERR_UNKNOWN,true,"myST->q0 is %f, should be %f\n", myST->q0, q0);
        return 2;
    }
    if (FLT_EPSILON < fabs(q1 - myST->q1)) {
        psError(PS_ERR_UNKNOWN,true,"myST->q1 is %f, should be %f\n", myST->q1, q1);
        return 3;
    }
    if (FLT_EPSILON < fabs(q2 - myST->q2)) {
        psError(PS_ERR_UNKNOWN,true,"myST->q2 is %f, should be %f\n", myST->q2, q2);
        return 4;
    }
    if (FLT_EPSILON < fabs(q3 - myST->q3)) {
        psError(PS_ERR_UNKNOWN,true,"myST->q0 is %f, should be %f\n", myST->q3, q3);
        return 5;
    }

    // Free data structure
    psFree(myST);

    return 0;
}

// We do a simple identity transformation on a few RA, DEC pairs.
psS32 testSphereRotApply1( void )
{
    /*
        psSphere *in = psSphereAlloc();
        psSphere *out = psSphereAlloc();
        psSphere *temp = NULL;
        psSphere *rc = NULL;
        psSphereRot *myST = psSphereRotAlloc(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(psSphereRotApply(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 = psSphereRotApply(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 = psSphereRotApply(NULL, NULL, in);
        if (rc != NULL) {
            psError(PS_ERR_UNKNOWN,true,"psSphereRotApply() 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 = psSphereRotApply(NULL, myST, NULL);
        if (rc != NULL) {
            psError(PS_ERR_UNKNOWN,true,"psSphereRotApply() did not return NULL");
            return 6;
        }

        psFree(myST);
        psFree(out);
        psFree(in);
    */
    return 0;
}

/******************************************************************************
testSphereRotApply2(): This test verifies that psSphereRotApply()
works properly.  We create two psSphereRots: 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?
    psSphereRotAlloc(X, Y, 0.0)
    psSphereRotAlloc(-X, -Y, 0.0)
 *****************************************************************************/
#define ERROR_PERCENT 0.01
psS32 testSphereRotApply2( void )
{
    /*
        psS32 testStatus = 0;
        psSphere in;
        psSphere out;
        psSphere out2;
        psSphereRot *mySphereRotForward = NULL;
        psSphereRot *mySphereRotReverse = NULL;


        mySphereRotForward = psSphereRotAlloc(DEG_TO_RAD(22.0),
                                              0.0,
                                              0.0);
        mySphereRotReverse = psSphereRotAlloc(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;

                psSphereRotApply(&out, mySphereRotForward, &in);
                psSphereRotApply(&out2, mySphereRotReverse, &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(mySphereRotForward);
        psFree(mySphereRotReverse);

        mySphereRotForward = psSphereRotAlloc(0.0,
                                              DEG_TO_RAD(33.0),
                                              0.0);
        mySphereRotReverse = psSphereRotAlloc(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;

                psSphereRotApply(&out, mySphereRotForward, &in);
                psSphereRotApply(&out2, mySphereRotReverse, &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(mySphereRotForward);
        psFree(mySphereRotReverse);

        return(testStatus);
        */
    return 0;
}

psS32 testSphereRotApplyCelestial( void)
{
    /*
        int numTestPoints = 3;
        // ICRS coordinates
        double alpha[] = {  0.0,        0.0,     180.0     };
        double delta[] = {  0.0,       90.0,      30.0     };
        //Ecliptic coordinates
        double lambda[] ={  0.0,       90.0,     167.072470};
        double beta[] =  {  0.0,       66.560719, 27.308813};
        // Galactic coordinates
        double l[] =     { 96.337272, 122.93192, 195.639488};
        double b[] =     {-60.188553,  27.12825,  78.353806};

        double t[] =     {  MJD_2000,  MJD_2000,   MJD_2100};

        double TOLERANCE = 0.001;



        for (int x = 0; x < numTestPoints; x++) {

            psTime* time = psTimeFromMJD(t[x]);
            psSphereRot* toEcliptic = psSphereRotICRSToEcliptic(time);
            psSphereRot* fromEcliptic = psSphereRotEclipticToICRS(time);
            psSphereRot* toGalactic = psSphereRotICRSToGalactic();
            psSphereRot* fromGalactic = psSphereRotGalacticToICRS();
            psFree(time);

            // set the ICRS coordinate
            psSphere* icrs = psSphereAlloc();
            icrs->r = DEG_TO_RAD(alpha[x]);
            icrs->d = DEG_TO_RAD(delta[x]);

            // apply/unapply Ecliptic
            psSphere* ecliptic = psSphereRotApply(NULL, toEcliptic, icrs);
            psSphere* icrsFromEcliptic = psSphereRotApply(NULL, fromEcliptic, ecliptic);

            // check ecliptic transforms for correctness
            if (abs(RAD_TO_DEG(ecliptic->r) - lambda[x]) > TOLERANCE ||
                    abs(RAD_TO_DEG(ecliptic->d) - beta[x]) > TOLERANCE) {
                psError(PS_ERR_UNKNOWN, false,
                        "Ecliptic tranformation incorrect.  Result is (%g,%g), expected (%g,%g)",
                        RAD_TO_DEG(ecliptic->r),RAD_TO_DEG(ecliptic->d),
                        lambda[x], beta[x]);
                return 1;
            }
            if (abs(RAD_TO_DEG(icrsFromEcliptic->r) - alpha[x]) > TOLERANCE ||
                    abs(RAD_TO_DEG(icrsFromEcliptic->d) - delta[x]) > TOLERANCE) {
                psError(PS_ERR_UNKNOWN, false,
                        "ICRS for Ecliptic tranformation incorrect.  Result is (%g,%g), expected (%g,%g)",
                        RAD_TO_DEG(icrsFromEcliptic->r),RAD_TO_DEG(icrsFromEcliptic->d),
                        alpha[x], delta[x]);
                return 2;
            }
            psFree(ecliptic);
            psFree(icrsFromEcliptic);

            // apply/unapply Galactic
            psSphere* galactic = psSphereRotApply(NULL, toGalactic, icrs);
            psSphere* icrsFromGalactic = psSphereRotApply(NULL, fromGalactic, galactic);

            // check ecliptic transforms for correctness
            if (abs(RAD_TO_DEG(galactic->r) - l[x]) > TOLERANCE ||
                    abs(RAD_TO_DEG(galactic->d) - b[x]) > TOLERANCE) {
                psError(PS_ERR_UNKNOWN, false,
                        "Galactic tranformation incorrect.  Result is (%g,%g), expected (%g,%g)",
                        RAD_TO_DEG(galactic->r),RAD_TO_DEG(galactic->d),
                        l[x], b[x]);
                return 3;
            }
            if (abs(RAD_TO_DEG(icrsFromGalactic->r) - alpha[x]) > TOLERANCE ||
                    abs(RAD_TO_DEG(icrsFromGalactic->d) - delta[x]) > TOLERANCE) {
                psError(PS_ERR_UNKNOWN, false,
                        "ICRS for Galactic tranformation incorrect.  Result is (%g,%g), expected (%g,%g)",
                        RAD_TO_DEG(icrsFromGalactic->r),RAD_TO_DEG(icrsFromGalactic->d),
                        alpha[x], delta[x]);
                return 4;
            }
            psFree(galactic);
            psFree(icrsFromGalactic);

            psFree(toEcliptic);
            psFree(fromEcliptic);
            psFree(toGalactic);
            psFree(fromGalactic);

        }
    */
    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 testSphereRotPrecess( 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;
}
