/** @file  psEarthOrientation.c
*
*  @brief Function implementations for earth orientation calculations
*  transformation
*
*  @ingroup EarthOrientation
*
*  @author Robert Daniel DeSonia, MHPCC
*
*  @version $Revision: 1.4.4.1 $ $Name: not supported by cvs2svn $
*  @date $Date: 2005-11-16 20:25:41 $
*
*  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
*/

#include <math.h>
#include <string.h>

#include "psEarthOrientation.h"
#include "psTime.h"
#include "psArray.h"
#include "psPolynomial.h"
#include "psVector.h"
#include "psMetadata.h"
#include "psMetadataConfig.h"
#include "psError.h"
#include "psErrorText.h"
#include "psMemory.h"
#include "psCoord.h"
#include "psConstants.h"

// Sun's Mass (src: Google's Calculator Service)
#define PS_M   1.98892e30 /* kilograms */
// Newton's Gravitational Constant (src:NIST)
#define PS_G   6.6742e-11 /* m^3/kg/s^2 */
// Speed of light in vacuum (src:NIST)
#define PS_C0  299792458.0 /* m/s */
// Average distance from earth to sun
#define PS_AU 149597890000.0 /* meters */
// Modified Julian Day 1/1/2000 00:00:00
#define MJD_2000  51544.0
// Days in Julian century
#define JULIAN_CENTURY 36525.0

static psArray* xTable = NULL;
static psArray* yTable = NULL;
static psArray* sTable = NULL;
static psPolynomial1D* xPoly = NULL;
static psPolynomial1D* yPoly = NULL;
static psPolynomial1D* sPoly = NULL;
static bool eocInitialized = false;

static bool eocInit()
{
    unsigned int nFail = 0;

    // Read config file
    psMetadata* eocMetadata = psMetadataConfigParse(NULL,
                              &nFail,
                              p_psGetConfigFileName(),
                              true);

    if(eocMetadata == NULL) {
        return false;
    } else if(nFail != 0) {
        return false;
    }

    bool success = false;
    // Get table format
    char* tableFormat = psMetadataLookupPtr(&success, eocMetadata, "psLib.eoc.precession.table.format");
    if(! success || tableFormat == NULL) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
                "psLib.eoc.precession.table.format");
        return false;
    }

    char* xTableName = psMetadataLookupPtr(&success, eocMetadata, "psLib.eoc.precession.table.file.x");
    if(! success || xTableName == NULL) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
                "psLib.eoc.precession.x.file");
        return false;
    }

    char* yTableName = psMetadataLookupPtr(&success, eocMetadata, "psLib.eoc.precession.table.file.y");
    if(! success || yTableName == NULL) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
                "psLib.eoc.precession.y.file");
        return false;
    }

    char* sTableName = psMetadataLookupPtr(&success, eocMetadata, "psLib.eoc.precession.table.file.s");
    if(! success || sTableName == NULL) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
                "psLib.eoc.precession.s.file");
        return false;
    }
    xTable = psVectorsReadFromFile(xTableName, tableFormat);
    yTable = psVectorsReadFromFile(yTableName, tableFormat);
    sTable = psVectorsReadFromFile(sTableName, tableFormat);
    if(xTable == NULL || yTable == NULL || sTable == NULL) {
        // XXX: need to move the error message to the error messages file.
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                "Failed to read the precession-nutation model tables.");
        return false;
    }

    psVector* xCoeff = psMetadataLookupPtr(&success, eocMetadata, "psLib.eoc.precession.poly.x");
    if(! success || xCoeff == NULL || xCoeff->type.type != PS_TYPE_F64) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
                "psLib.eoc.precession.poly.x");
        return false;
    }
    psVector* yCoeff = psMetadataLookupPtr(&success, eocMetadata, "psLib.eoc.precession.poly.y");
    if(! success || yCoeff == NULL || yCoeff->type.type != PS_TYPE_F64) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
                "psLib.eoc.precession.poly.y");
        return false;
    }
    psVector* sCoeff = psMetadataLookupPtr(&success, eocMetadata, "psLib.eoc.precession.poly.s");
    if(! success || sCoeff == NULL || sCoeff->type.type != PS_TYPE_F64) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED,
                "psLib.eoc.precession.poly.s");
        return false;
    }
    xPoly = psPolynomial1DAlloc(xCoeff->n, PS_POLYNOMIAL_ORD);
    memcpy(xPoly->coeff, xCoeff->data.F64,PSELEMTYPE_SIZEOF(PS_TYPE_F64)*xCoeff->n);
    yPoly = psPolynomial1DAlloc(yCoeff->n, PS_POLYNOMIAL_ORD);
    memcpy(yPoly->coeff, yCoeff->data.F64,PSELEMTYPE_SIZEOF(PS_TYPE_F64)*yCoeff->n);
    sPoly = psPolynomial1DAlloc(sCoeff->n, PS_POLYNOMIAL_ORD);
    memcpy(sPoly->coeff, sCoeff->data.F64,PSELEMTYPE_SIZEOF(PS_TYPE_F64)*sCoeff->n);

    return true;
}

bool p_psEOCFinalize(void)
{
    psFree(xTable);
    psFree(yTable);
    psFree(sTable);

    xTable = NULL;
    yTable = NULL;
    sTable = NULL;

    psFree(xPoly);
    psFree(yPoly);
    psFree(sPoly);

    xPoly = NULL;
    yPoly = NULL;
    sPoly = NULL;

    eocInitialized = false;

    return true;
}

psSphere *psAberration(psSphere *apparent, const psSphere *actual, const psSphere *direction, double speed)
{


    return NULL;
}

psSphere *psGravityDeflection(psSphere *apparent, psSphere *actual, psSphere *sun)
{
    // calculating the apparent angle from the actual angle and the sun position

    // first, calculate the angle between the sun vector and the actual vector

    // Moving to cartesian first:  XXX -- is this required?
    psCube* sunVector = psSphereToCube(sun);
    psCube* actualVector = psSphereToCube(actual);

    // use dot product to calculate the angle of separation
    // N.B., assuming the psSphereToCube function returns a unit vector.
    double theta = acos(sunVector->x*actualVector->x +
                        sunVector->y*actualVector->y +
                        sunVector->z*actualVector->z);

    double r0 = PS_AU * tan(theta);

    double deflection = 4.0*PS_G*PS_M/(PS_C0*PS_C0*r0);

    // make sure the deflection is not greater than 1.75 arcsec
    double limit = SEC_TO_RAD(1.75);
    if (deflection > limit) {
        deflection = limit;
    }

    // bend the actual vector away from the sun vector by deflection angle.
    // XXX: Not sure how to do this.  Dave thinks the formula should be:
    //      theta = atan(r0/d)*tan(deflection), phi = thete/tan(deflection)

    return NULL;
}


double psEOC_ParallaxFactor(const psSphere *coords, const psTime *time)
{


    return NAN;
}

psEarthPole *psEOC_PrecessionModel(const psTime *time)
{
    // Check for null parameter
    PS_ASSERT_PTR_NON_NULL(time, NULL);

    // Convert psTime to MJD
    double MJD = psTimeToMJD(time);

    // Calculate number of Julian centuries since 2000
    double t = ( MJD - MJD_2000 ) / JULIAN_CENTURY;
    double t2 = t*t;
    double t3 = t*t*t;
    double t4 = t*t*t*t;

    // N.B., the following formulae are from the ADD
    // Mean Anomaly of the Moon
    double F[14];
    F[0] = DEG_TO_RAD(134.96340251) +
           SEC_TO_RAD(1717915923.2178)*t +
           SEC_TO_RAD(31.8792)*t2 +
           SEC_TO_RAD(0.051635)*t3 -
           SEC_TO_RAD(0.00024470)*t4;

    // Mean Anomaly of the Sun
    F[1] = DEG_TO_RAD(357.52910918) +
           SEC_TO_RAD(129596581.0481)*t -
           SEC_TO_RAD(0.5532)*t2 +
           SEC_TO_RAD(0.000136)*t3 -
           SEC_TO_RAD(0.00001149)*t4;

    // L Omega
    F[2] = DEG_TO_RAD(93.27209062) +
           SEC_TO_RAD(1739527262.8478)*t -
           SEC_TO_RAD(12.7512)*t2 -
           SEC_TO_RAD(0.001037)*t3 +
           SEC_TO_RAD(0.00000417)*t4;

    // Mean Elongation of the Moon from the Sun
    F[3] = DEG_TO_RAD(297.85019547) +
           SEC_TO_RAD(1602961601.2090)*t -
           SEC_TO_RAD(6.3706)*t2 +
           SEC_TO_RAD(0.006593)*t3 -
           SEC_TO_RAD(0.00003169)*t4;

    // Mean Longitude of the Ascending Node of the Moon
    F[4] = DEG_TO_RAD(125.04455501) -
           SEC_TO_RAD(6962890.5431)*t +
           SEC_TO_RAD(7.4722)*t2 +
           SEC_TO_RAD(0.007702)*t3 -
           SEC_TO_RAD(0.0000593)*t4;

    F[5] = 4.402608842 + 2608.7903141574*t;
    F[6] = 3.176146697 + 1021.3285546211*t;
    F[7] = 1.753470314 + 628.3075849991*t;
    F[8] = 6.203480913 + 334.0612426700*t;
    F[9] = 0.599546497 + 52.9690962641*t;
    F[10] = 0.874016757 + 21.3299104960*t;
    F[11] = 5.481293872 + 7.4781598567*t;
    F[12] = 5.311886287 + 3.8133035638*t;
    F[13] = 0.024381750*t + 0.00000538691*t2;

    // Check if EOC data loaded
    if(! eocInitialized) {
        eocInitialized = eocInit();
        if(!eocInitialized) {
            // XXX: Move error message.
            psError(PS_ERR_UNKNOWN, false,
                    "Could not initialize EOC tables -- check data files.");
            return NULL;
        }
    }

    // calculate the polynomial portion first
    double X = psPolynomial1DEval(xPoly,t);
    double Y = psPolynomial1DEval(yPoly,t);
    double S = psPolynomial1DEval(sPoly,t);

    // now calculate the non-poly portion from the tables

    psF64* cols[17];
    for (int lcv = 0; lcv < 17; lcv++) {
        cols[lcv] = ((psVector*)(xTable->data[lcv]))->data.F64;
    }
    int numRows = ((psVector*)(xTable->data[0]))->n;
    for (int lcv = 0; lcv < numRows; lcv++) {
        // arg = sum w_(i,j,k)*F_k
        double arg = 0.0;
        for (int k = 0; k < 14; k++) {
            arg += cols[k+3][lcv]*F[k];
        }
        double tj = pow(t,cols[0][lcv]);
        double as = cols[1][lcv];
        double ac = cols[2][lcv];

        X += as*tj*sin(arg) + ac*tj*cos(arg);
    }

    for (int lcv = 0; lcv < 17; lcv++) {
        cols[lcv] = ((psVector*)(yTable->data[lcv]))->data.F64;
    }
    numRows = ((psVector*)(yTable->data[0]))->n;
    for (int lcv = 0; lcv < numRows; lcv++) {
        // arg = sum w_(i,j,k)*F_k
        double arg = 0.0;
        for (int k = 0; k < 14; k++) {
            arg += cols[k+3][lcv]*F[k];
        }
        double tj = pow(t,cols[0][lcv]);
        double as = cols[1][lcv];
        double ac = cols[2][lcv];

        Y += as*tj*sin(arg) + ac*tj*cos(arg);
    }

    for (int lcv = 0; lcv < 17; lcv++) {
        cols[lcv] = ((psVector*)(sTable->data[lcv]))->data.F64;
    }
    numRows = ((psVector*)(sTable->data[0]))->n;
    for (int lcv = 0; lcv < numRows; lcv++) {
        // arg = sum w_(i,j,k)*F_k
        double arg = 0.0;
        for (int k = 0; k < 14; k++) {
            arg += cols[k+3][lcv]*F[k];
        }
        double tj = pow(t,cols[0][lcv]);
        double as = cols[1][lcv];
        double ac = cols[2][lcv];

        S += as*tj*sin(arg) + ac*tj*cos(arg);
    }

    // now, the tables for S actually gives S + XY/2, so let's get the real S now
    S -= X*Y/2.0;

    psEarthPole* pole = psAlloc(sizeof(psEarthPole));
    pole->x = X;
    pole->y = Y;
    pole->s = S;

    return pole;
}


psEarthPole *psEOC_PrecessionCorr(const psTime *time, psTimeBulletin bulletin)
{
    return NULL;
}


psSphereRot* psSphereRot_CEOtoGCRS(const psEarthPole *pole)
{
    return NULL;
}


psSphereRot* psSphereRot_TEOtoCEO(const psTime *time)
{
    return NULL;
}


psEarthPole* psEOC_GetPolarMotion(const psTime *time, psTimeBulletin bulletin)
{
    return NULL;
}


psEarthPole* psEOC_PolarTideCorr(const psTime *time)
{
    return NULL;
}


psEarthPole* psEOC_NutationCorr(psTime *time)
{
    return NULL;
}


psSphereRot* psSphereRot_ITRStoTEO(const psEarthPole* motion)
{
    return NULL;
}
