/** @file  psCoord.c
*
*  @brief Contains basic coordinate transformation definitions and operations
*
*  This file defines the basic types for astronomical coordinate
*  transformation
*
*  @ingroup CoordinateTransform
*
*  @author George Gusciora, MHPCC
*
*  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
*  @date $Date: 2004-10-14 01:17:50 $
*
*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
*/
/******************************************************************************/
/*  INCLUDE FILES                                                             */
/******************************************************************************/
#include "psType.h"
#include "psImage.h"
#include "psArray.h"
#include "psList.h"
#include "psCoord.h"
#include "psMemory.h"
#include "psAbort.h"
#include "psTime.h"
#include "psConstants.h"
#include "psError.h"
#include <math.h>
#include <float.h>
/******************************************************************************/
/*  DEFINE STATEMENTS                                                         */
/******************************************************************************/

#define PS_COT(X) (1.0 / atan(X))

/******************************************************************************/
/*  TYPE DEFINITIONS                                                          */
/******************************************************************************/

// None

/*****************************************************************************/
/*  GLOBAL VARIABLES                                                         */
/*****************************************************************************/

// None

/*****************************************************************************/
/*  FILE STATIC VARIABLES                                                    */
/*****************************************************************************/

// None

/*****************************************************************************/
/*  FUNCTION IMPLEMENTATION - LOCAL                                          */
/*****************************************************************************/

static float p_psArg(float x, float y);

/******************************************************************************
XXX: Verify this arc tan function.
XXX: macro
XXX: I don't think this arg() is correct.
 *****************************************************************************/

float p_psArg(float x,
              float y)
{
    if (y!= 0.0) {
        return(atan(y/x));
    } else {
        return(0.5 * M_PI);
    }
    return (0.0);
    /*
        if (x > 0) {
            return (atan(y / x));
        } else if ((x == 0) && (y == 0)) {
            return (0.5 * M_PI);
        } else if ((x == 0) && (y == 0)) {
            return (-0.5 * M_PI);
        } else if ((x == 0) && (y == 0)) {
            return (M_PI + atan(y / x));
        } else if ((x == 0) && (y == 0)) {
            return (-M_PI + atan(y / x));
        }
        psAbort(__func__, "Unacceptable range for (p_psArg(%f, %f).\n", x, y);
    */
}

/*****************************************************************************/
/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
/*****************************************************************************/
psPlane* psPlaneTransformApply(psPlane* out,
                               const psPlaneTransform* transform,
                               const psPlane* coords)
{
    PS_CHECK_NULL_PTR_RETURN_NULL(transform);
    PS_CHECK_NULL_PTR_RETURN_NULL(coords);
    if (out == NULL) {
        out = (psPlane* ) psAlloc(sizeof(psPlane));
    }
    out->x = psDPolynomial2DEval(coords->x, coords->y, transform->x);
    out->y = psDPolynomial2DEval(coords->x, coords->y, transform->y);
    return (out);
}

/******************************************************************************
This transformation takes into account parameters beyond an objects spatial
coordinates: term3 and term4 (magnitude and color).
 *****************************************************************************/
psPlane* psPlaneDistortApply(psPlane* out,
                             const psPlaneDistort* transform,
                             const psPlane* coords,
                             float color,
                             float magnitude)
{
    PS_CHECK_NULL_PTR_RETURN_NULL(transform);
    PS_CHECK_NULL_PTR_RETURN_NULL(coords);
    if (out == NULL) {
        out = (psPlane* ) psAlloc(sizeof(psPlane));
    }
    out->x = psDPolynomial4DEval(coords->x, coords->y, color, magnitude, transform->x);
    out->y = psDPolynomial4DEval(coords->x, coords->y, color, magnitude, transform->y);
    return (out);
}

/******************************************************************************
This function prototype has been modified since the SDRS.
 *****************************************************************************/
psSphereTransform* psSphereTransformAlloc(double NPlat,
        double Xo,
        double xo)
{
    psSphereTransform* tmp = (psSphereTransform* ) psAlloc(sizeof(psSphereTransform));

    tmp->sinPhi = sin(NPlat);
    tmp->cosPhi = cos(NPlat);
    tmp->Xo = Xo;
    tmp->xo = xo;

    return (tmp);
}

/******************************************************************************
This algorithm comes from an email from Gene.  I assume (x,y) corresponds to
(r,d) in the sphere coordinates.
 
XXX: In Gene's email, there are different variables with similar names (y
and Y) and in the code, there's cos_y, cos_Y, and cos(y).  Verify that there
are no typo's.
 *****************************************************************************/
psSphere* psSphereTransformApply(psSphere* out,
                                 const psSphereTransform* transform,
                                 const psSphere* coord)
{
    PS_CHECK_NULL_PTR_RETURN_NULL(transform);
    PS_CHECK_NULL_PTR_RETURN_NULL(coord);
    double sinY = 0.0;
    double cosY = 0.0;
    double sinX = 0.0;
    double cosX = 0.0;
    double x = 0.0;
    double y = 0.0;
    double dx = 0.0;

    if (out == NULL) {
        out = (psSphere* ) psAlloc(sizeof(psSphere));
    }

    x = coord->r;
    y = coord->d;
    dx = x - transform->xo;
    sinY = cos(y) * sin(dx) * transform->sinPhi + sin(y) * transform->cosPhi;
    cosY = sqrt(1.0 - sinY * sinY);
    sinX = (cos(y) * sin(dx) * transform->cosPhi - sin(y) * transform->sinPhi) / cos(y);
    cosX = cos(y) * cos(dx) / cos(y);

    out->r = atan2(sinX, cosX) + transform->Xo;
    out->d = atan2(sinY, cosY);

    return (out);
}

psSphereTransform* psSphereTransformICRStoEcliptic(psTime time)
{
    struct tm *tmTime = psTimeToTM(&time);
    double year = (double)(1900 + tmTime->tm_year);
    double T = year / 100.0;
    double phi = -23.452294 + 0.013013 * T + 0.000001639 * T * T - 0.000000503 * T * T * T;
    double Xo = 0.0;
    double xo = 0.0;

    psFree(tmTime);
    return (psSphereTransformAlloc(phi, Xo, xo));
}

psSphereTransform* psSphereTransformEcliptictoICRS(psTime time)
{
    struct tm *tmTime = psTimeToTM(&time);
    double year = (double)(1900 + tmTime->tm_year);
    double T = year / 100.0;
    double phi = +23.452294 - 0.013013 * T - 0.000001639 * T * T + 0.000000503 * T * T * T;
    double Xo = 0.0;
    double xo = 0.0;

    psFree(tmTime);
    return (psSphereTransformAlloc(phi, Xo, xo));
}

psSphereTransform* psSphereTransformICRStoGalatic(void)
{
    return (psSphereTransformAlloc(62.6, 282.25, 33.0));
}

psSphereTransform* psSphereTransformGalatictoICRS(void)
{
    return (psSphereTransformAlloc(-62.6, 33.0, 282.25));
}

/******************************************************************************
XXX: Waiting for the definition of the PS_PROJ_PAR projection.
XXX: Waiting for the definition of the PS_PROJ_GLS projection.
XXX: get rid of the psAborts
 *****************************************************************************/
psPlane* psProject(const psSphere* coord,
                   const psProjection* projection)
{
    PS_CHECK_NULL_PTR_RETURN_NULL(coord);
    PS_CHECK_NULL_PTR_RETURN_NULL(projection);

    float R = 0.0;
    float alpha = 0.0;
    psPlane* tmp = (psPlane* ) psAlloc(sizeof(psPlane));

    if (projection->type == PS_PROJ_TAN) {
        R = PS_COT(coord->r) * (180.0 / M_PI);
        tmp->x = R * sin(coord->d);
        tmp->y = R * cos(coord->d);

    } else if (projection->type == PS_PROJ_SIN) {
        R = cos(coord->r) * (180.0 / M_PI);
        tmp->x = R * sin(coord->d);
        tmp->y = R * cos(coord->d);

    } else if (projection->type == PS_PROJ_CAR) {
        tmp->x = coord->d;
        tmp->y = coord->r;

    } else if (projection->type == PS_PROJ_MER) {
        tmp->x = coord->d;
        tmp->y = log(tan(45.0 + (0.5 * coord->r))) * 180.0 / M_PI;

    } else if (projection->type == PS_PROJ_AIT) {
        alpha = 1.0 / ((180.0 / M_PI) * sqrt(1.0 + (cos(coord->r) *
                                             cos(0.5 * coord->d) * 0.5)));

        tmp->x = 2.0 * alpha * cos(coord->r) * sin(0.5 * coord->d);
        tmp->y = alpha * sin(coord->d);

    } else if (projection->type == PS_PROJ_PAR) {
        psAbort(__func__, "The projection type PS_PROJ_PAR is undefined.\n");

    } else if (projection->type == PS_PROJ_GLS) {
        psAbort(__func__, "The projection type PS_PROJ_GLS is undefined.\n");
    } else {
        psAbort(__func__, "Unknown projection type.\n");
    }

    return (tmp);
}

/******************************************************************************
XXX: Waiting for the definition of the PS_PROJ_PAR projection.
XXX: Waiting for the definition of the PS_PROJ_GLS projection.
XXX: get rid of the psAborts
 *****************************************************************************/
psSphere* psDeproject(const psPlane* coord,
                      const psProjection* projection)
{
    PS_CHECK_NULL_PTR_RETURN_NULL(coord);
    PS_CHECK_NULL_PTR_RETURN_NULL(projection);

    float R = 0.0;
    float chu = 0.0;
    float chu1 = 0.0;
    float chu2 = 0.0;
    psSphere* tmp = (psSphere* ) psAlloc(sizeof(psSphere));

    if (projection->type == PS_PROJ_TAN) {
        R = sqrt((coord->x * coord->x) + (coord->y * coord->y));
        tmp->d = p_psArg(-coord->y, coord->x);
        tmp->r = atan(180.0 / (R * M_PI));

    } else if (projection->type == PS_PROJ_SIN) {
        R = sqrt((coord->x * coord->x) + (coord->y * coord->y));
        tmp->d = p_psArg(-coord->y, coord->x);
        tmp->r = acos((R * M_PI) / 180.0);

    } else if (projection->type == PS_PROJ_CAR) {
        tmp->d = coord->x;
        tmp->r = coord->y;

    } else if (projection->type == PS_PROJ_MER) {
        tmp->d = coord->x;
        tmp->r = (2.0 * atan(exp((coord->y * M_PI / 180.0)))) - 180.0;

    } else if (projection->type == PS_PROJ_AIT) {
        chu1 = (coord->x * M_PI) / 720.0;
        chu1 *= chu1;
        chu2 = (coord->y * M_PI) / 360.0;
        chu2 *= chu2;
        chu = sqrt(1.0 - chu1 - chu2);
        tmp->d = 2.0 * p_psArg((2.0 * chu * chu) - 1.0, (coord->x * chu * M_PI) / 360.0);
        tmp->r = asin((coord->y * chu * M_PI) / 180.0);

    } else if (projection->type == PS_PROJ_PAR) {
        psAbort(__func__, "The projection type PS_PROJ_PAR is undefined.\n");

    } else if (projection->type == PS_PROJ_GLS) {
        psAbort(__func__, "The projection type PS_PROJ_GLG is undefined.\n");
    } else {
        psAbort(__func__, "Unknown projection type.\n");
    }

    return (tmp);
}

/******************************************************************************
The basic idea is to project both positions onto the linear plane, with
position1 at the center, then calculate the linear offset between those
projections.
 
XXX: Do I need to check for unacceptable transformation parameters?  Maybe,
     if the points are on the North/South Pole, etc?
 
XXX: Do I need to somehow scale this projection?
 
XXX: Does PS_LINEAR mode make sense?  The result must be returned in psSphere
     regardless of the mode.
 
XXX: How to compound errors?
 *****************************************************************************/
psSphere* psSphereGetOffset(const psSphere* restrict position1,
                            const psSphere* restrict position2,
                            psSphereOffsetMode mode,
                            psSphereOffsetUnit unit)
{
    PS_CHECK_NULL_PTR_RETURN_NULL(position1);
    PS_CHECK_NULL_PTR_RETURN_NULL(position2);

    psPlane* lin;
    psProjection proj;
    psSphere* tmp;

    if (mode == PS_LINEAR) {
        proj.R = position1->r;
        proj.D = position1->d;
        proj.Xs = 1.0;
        proj.Ys = 1.0;
        proj.type = PS_PROJ_TAN;

        lin = psProject(position2, &proj);
        tmp = psDeproject(lin, &proj);
        psFree(lin);

        // XXX: Do we need to convert units in tmp?
        return (tmp);
    } else if (mode == PS_SPHERICAL) {
        tmp = (psSphere* ) psAlloc(sizeof(psSphere));
        tmp->r = position2->r - position1->r;
        tmp->d = position2->d - position1->d;

        // Wrap these to an acceptable range.  This assumes that all
        // angles are in radians.
        tmp->r = fmod(tmp->r, 2*M_PI);
        tmp->d = fmod(tmp->d, 2*M_PI);
        tmp->rErr = 0.0;
        tmp->dErr = 0.0;

        if (unit == PS_ARCSEC) {
            tmp->r = (tmp->r * 180.0 * 60.0 * 60.0) / M_PI;
            tmp->d = (tmp->d * 180.0 * 60.0 * 60.0) / M_PI;
        } else if (unit == PS_ARCMIN) {
            tmp->r = (tmp->r * 180.0 * 60.0) / M_PI;
            tmp->d = (tmp->d * 180.0 * 60.0) / M_PI;
        } else if (unit == PS_DEGREE) {
            tmp->r = (tmp->r * 180.0) / M_PI;
            tmp->d = (tmp->d * 180.0) / M_PI;
        } else if (unit == PS_RADIAN) {}
        else {
            psAbort(__func__, "Unknown offset unit: 0x%x\n", unit);
        }

        return (tmp);
    }
    psAbort(__func__, "Unrecognized offset mode\n");
    return (NULL);
}

/******************************************************************************
XXX: Do we need to check for unacceptable transformation parameters?  Maybe,
     if the points are on the North/South Pole, etc?
 
XXX: Do we need to somehow scale this projection?
 
XXX: I copied the algorithm from the ADD exactly.
 
XXX: Should we compound errors?
 *****************************************************************************/

psSphere* psSphereSetOffset(const psSphere* restrict position,
                            const psSphere* restrict offset,
                            psSphereOffsetMode mode,
                            psSphereOffsetUnit unit)
{
    PS_CHECK_NULL_PTR_RETURN_NULL(position);
    PS_CHECK_NULL_PTR_RETURN_NULL(offset);

    psPlane lin;
    psSphere* tmp;
    psProjection proj;
    double tmpR = 0.0;
    double tmpD = 0.0;

    if (mode == PS_LINEAR) {
        proj.R = position->r;
        proj.D = position->d;
        proj.Xs = 1.0;
        proj.Ys = 1.0;
        proj.type = PS_PROJ_TAN;

        lin.x = offset->r;
        lin.y = offset->d;

        tmp = psDeproject(&lin, &proj);
        return (tmp);

    } else if (mode == PS_SPHERICAL) {
        if (unit == PS_ARCSEC) {
            tmpR = (M_PI * offset->r) / (180.0 * 60.0 * 60.0);
            tmpD = (M_PI * offset->d) / (180.0 * 60.0 * 60.0);
        } else if (unit == PS_ARCMIN) {
            tmpR = (M_PI * offset->r) / (180.0 * 60.0);
            tmpD = (M_PI * offset->d) / (180.0 * 60.0);
        } else if (unit == PS_DEGREE) {
            tmpR = (M_PI * offset->r) / (180.0);
            tmpD = (M_PI * offset->d) / (180.0);
        } else if (unit == PS_RADIAN) {
            tmpR = offset->r;
            tmpD = offset->d;
        } else {
            psAbort(__func__, "Unknown offset unit: 0x%x\n", unit);
        }

        tmp = (psSphere* ) psAlloc(sizeof(psSphere));
        tmp->r = position->r + tmpR;
        tmp->r = fmod(tmp->r, 2.0*M_PI);
        tmp->d = position->d + tmpD;
        tmp->d = fmod(tmp->d, 2.0*M_PI);
        tmp->rErr = 0.0;
        tmp->dErr = 0.0;

        return (tmp);
    }
    psAbort(__func__, "Unrecognized offset mode\n");
    return (NULL);
}
