/** @file  tst_psCoord.c
*
*  @brief The code will ...
*
*
*  @author GLG, MHPCC
*
*  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
*  @date $Date: 2005-07-23 02:44:09 $
*
*  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 1.0
#define DELTA_P 2.0
#define PHI_P 3.0

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

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

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