/** @file  psCoord.c
*
*  @brief Contains basic coordinate transformation definitions and operations
*
*  This file defines the basic types for astronomical coordinate
*  transformation
*
*  @ingroup CoordinateTransform
*
*  @author GLG, MHPCC
*
*  @version $Revision: 1.46 $ $Name: not supported by cvs2svn $
*  @date $Date: 2005-01-14 23:27:55 $
*
*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
*/
/******************************************************************************/
/*  INCLUDE FILES                                                             */
/******************************************************************************/
#include "psType.h"
#include "psCoord.h"
#include "psMemory.h"
#include "psTime.h"
#include "psConstants.h"
#include "psError.h"
#include "psLogMsg.h"
#include "psAstronomyErrors.h"
#include <math.h>
#include <float.h>
/******************************************************************************/
/*  DEFINE STATEMENTS                                                         */
/******************************************************************************/

// None

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

// None

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

// None

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

// None

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

/*****************************************************************************/
/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
/*****************************************************************************/
static void planeTransformFree(psPlaneTransform *pt)
{
    psFree(pt->x);
    psFree(pt->y);
}

psPlaneTransform* psPlaneTransformAlloc(psS32 n1, psS32 n2)
{
    PS_INT_CHECK_NON_NEGATIVE(n1, NULL);
    PS_INT_CHECK_NON_NEGATIVE(n2, NULL);

    psPlaneTransform *pt = psAlloc(sizeof(psPlaneTransform));
    pt->x = psDPolynomial2DAlloc(n1, n2, PS_POLYNOMIAL_ORD);
    pt->y = psDPolynomial2DAlloc(n1, n2, PS_POLYNOMIAL_ORD);

    p_psMemSetDeallocator(pt, (psFreeFcn) planeTransformFree);
    return(pt);
}

psPlane* psPlaneTransformApply(psPlane* out,
                               const psPlaneTransform* transform,
                               const psPlane* coords)
{
    PS_PTR_CHECK_NULL(transform, NULL);
    PS_PTR_CHECK_NULL(transform->x, NULL);
    PS_PTR_CHECK_NULL(transform->y, NULL);
    PS_PTR_CHECK_NULL(coords, NULL);

    if (out == NULL) {
        out = (psPlane* ) psAlloc(sizeof(psPlane));
    }
    out->x = psDPolynomial2DEval(
                 transform->x,
                 coords->x,
                 coords->y
             );
    out->y = psDPolynomial2DEval(
                 transform->y,
                 coords->x,
                 coords->y
             );
    return (out);
}

static void planeDistortFree(psPlaneDistort *pt)
{
    psFree(pt->x);
    psFree(pt->y);
}

psPlaneDistort* psPlaneDistortAlloc(psS32 n1, psS32 n2, psS32 n3, psS32 n4)
{
    PS_INT_CHECK_NON_NEGATIVE(n1, NULL);
    PS_INT_CHECK_NON_NEGATIVE(n2, NULL);
    PS_INT_CHECK_NON_NEGATIVE(n3, NULL);
    PS_INT_CHECK_NON_NEGATIVE(n4, NULL);

    psPlaneDistort *pt = psAlloc(sizeof(psPlaneDistort));
    pt->x = psDPolynomial4DAlloc(n1, n2, n3, n4, PS_POLYNOMIAL_ORD);
    pt->y = psDPolynomial4DAlloc(n1, n2, n3, n4, PS_POLYNOMIAL_ORD);

    p_psMemSetDeallocator(pt, (psFreeFcn) planeDistortFree);
    return(pt);
}

/******************************************************************************
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_PTR_CHECK_NULL(transform, NULL);
    PS_PTR_CHECK_NULL(transform->x, NULL);
    PS_PTR_CHECK_NULL(transform->y, NULL);
    PS_PTR_CHECK_NULL(coords, NULL);

    if (out == NULL) {
        out = (psPlane* ) psAlloc(sizeof(psPlane));
    }
    out->x = psDPolynomial4DEval(
                 transform->x,
                 coords->x,
                 coords->y,
                 color,
                 magnitude
             );
    out->y = psDPolynomial4DEval(
                 transform->y,
                 coords->x,
                 coords->y,
                 color,
                 magnitude
             );
    return (out);
}

/******************************************************************************
alpha is LONGITUDE
delta is LATITUDE
 
    alphaP: Take the target pole in the source system; calculate its LONGITUDE
     in the target system.  That longitude is alphaP.
    DeltaP: Take the target pole in the source system; calculate its LATITUDE
     in the target system.  That longitude is deltaP.
    phiP:   This is the LONGITUDE of the ascending node in the target system.
 *****************************************************************************/
psSphereTransform* psSphereTransformAlloc(psF64 alphaP,
        psF64 deltaP,
        psF64 phiP)
{
    psSphereTransform* tmp = (psSphereTransform* ) psAlloc(sizeof(psSphereTransform));

    tmp->cosDeltaP = cos(deltaP);
    tmp->sinDeltaP = sin(deltaP);
    tmp->alphaP = alphaP;
    tmp->phiP = phiP;

    return (tmp);
}

/******************************************************************************
XXX: Private Function.
 
piNormalize(): take an input angle in radians and convert it to the range 0:PI.
 *****************************************************************************/
psF32 piNormalize(psF32 angle)
{
    while (angle < FLT_EPSILON) {
        angle+=M_PI;
    }

    while (angle >= M_PI) {
        angle-=M_PI;
    }
    return(angle);
}

/******************************************************************************
XXX: We convert Right Ascension angles to the range 0:PI.  Is that acceptable?
XXX: Should we do something for Declination as well?
 *****************************************************************************/
psSphere* psSphereTransformApply(psSphere* out,
                                 const psSphereTransform* transform,
                                 const psSphere* coord)
{
    PS_PTR_CHECK_NULL(transform, NULL);
    PS_PTR_CHECK_NULL(coord, NULL);

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

    psF64 alpha = coord->r;
    psF64 delta = coord->d;
    psF64 alphaMinusAlphaP = alpha - transform->alphaP;

    psF64 eq55 = (sin(delta) * transform->cosDeltaP) -
                 (cos(delta) * transform->sinDeltaP * sin(alphaMinusAlphaP));

    psF64 eq56 = (cos(delta) * transform->cosDeltaP * sin(alphaMinusAlphaP)) +
                 (sin(delta) * transform->sinDeltaP);

    psF64 eq57 = cos(delta) * cos(alphaMinusAlphaP);

    psF64 theta = asin(eq55);
    psF64 phi = atan2(eq56, eq57) + transform->alphaP;


    out->r = piNormalize(phi);
    out->d = theta;

    return(out);
}

psSphereTransform* psSphereTransformICRSToEcliptic(psTime *time)
{
    PS_PTR_CHECK_NULL(time, NULL);
    psF64 MJD = psTimeToMJD(time);
    struct tm *tmTime = psTimeToTM(time);
    psF64 year = (psF64)(1900 - MJD + tmTime->tm_year);
    psF64 T = year / 100.0;
    psF64 alphaP = 0.0;
    psF64 deltaP = DEG_TO_RAD(23.0) +
                   SEC_TO_RAD(27.0) +
                   MIN_TO_RAD(8.0) -
                   (SEC_TO_RAD(46.845) * T) -
                   (SEC_TO_RAD(0.0059) * T * T) +
                   (SEC_TO_RAD(0.00181) * T * T * T);
    psF64 phiP = 0.0;

    psFree(tmTime);

    // Don't neglect the minus sign on deltaP (bug 244):
    return (psSphereTransformAlloc(alphaP, -deltaP, phiP));
}


psSphereTransform* psSphereTransformEclipticToICRS(psTime *time)
{
    PS_PTR_CHECK_NULL(time, NULL);
    psF64 MJD = psTimeToMJD(time);
    struct tm *tmTime = psTimeToTM(time);
    psF64 year = (psF64)(1900 - MJD + tmTime->tm_year);
    psF64 T = year / 100.0;
    psF64 alphaP = 0.0;
    psF64 deltaP = DEG_TO_RAD(23.0) +
                   SEC_TO_RAD(27.0) +
                   MIN_TO_RAD(8.0) -
                   (SEC_TO_RAD(46.845) * T) -
                   (SEC_TO_RAD(0.0059) * T * T) +
                   (SEC_TO_RAD(0.00181) * T * T * T);
    psF64 phiP = 0.0;

    psFree(tmTime);
    return (psSphereTransformAlloc(alphaP, deltaP, phiP));
}

// XXX: This is bug 245: alphaP swaps with phiP from psSphereTransformGalaticToICRS()
psSphereTransform* psSphereTransformICRSToGalatic(void)
{
    psF64 alphaP = DEG_TO_RAD(32.93192);
    psF64 deltaP = DEG_TO_RAD(-62.87175);
    psF64 phiP = DEG_TO_RAD(282.85948);

    return (psSphereTransformAlloc(alphaP, deltaP, phiP));
}

psSphereTransform* psSphereTransformGalaticToICRS(void)
{
    psF64 alphaP = DEG_TO_RAD(282.85948);
    psF64 deltaP = DEG_TO_RAD(62.87175);
    psF64 phiP = DEG_TO_RAD(32.93192);

    return (psSphereTransformAlloc(alphaP, deltaP, phiP));
}

/******************************************************************************
XXX: Waiting for the definition of the PS_PROJ_PAR projection.
XXX: Waiting for the definition of the PS_PROJ_GLS projection.
XXX: Must apply scaling at the end.
 *****************************************************************************/
psPlane* psProject(const psSphere* coord,
                   const psProjection* projection)
{
    PS_PTR_CHECK_NULL(coord, NULL);
    PS_PTR_CHECK_NULL(projection, NULL);

    psPlane* tmp = (psPlane* ) psAlloc(sizeof(psPlane));
    if (projection->type == PS_PROJ_TAN) {

        /*
                // ********************************************
                // From the ADD
                // ********************************************
                // psLogMsg(__func__, PS_LOG_WARN, "WARNING: psProject(): These projections don't work.");
                // return(NULL);
                psF32 sinTheta;
                psF32 cosThetaCosPhi;
                psF32 cosThetaSinPhi;
         
                sinTheta = (sin(coord->d) * sin(projection->D)) + 
                           (cos(coord->d) * cos(projection->D) * (cos(coord->r - projection->R)));
                cosThetaCosPhi = (sin(coord->d) * cos(projection->D)) - 
                                 (cos(coord->d) * sin(projection->D) * (cos(coord->r - projection->R)));
                cosThetaSinPhi = - cos(coord->d) * sin(coord->r - projection->R);
         
                tmp->x =  -cosThetaSinPhi / sinTheta;
                tmp->y = cosThetaCosPhi / sinTheta;
         
        */
        // ********************************************
        // From the mathworks: http://mathworld.wolfram.com/GnomonicProjection.html
        // delta_0 and phi_1 are the projection centers, not the point being projected.
        // (delta_0, phi_1) == (projection->R, projection->D)
        // ********************************************
        psF32 cosC = (sin(projection->D) * sin(coord->d)) +
                     (cos(projection->D) * cos(coord->d) * cos(coord->r - projection->R));
        tmp->x = (cos(coord->d) * sin(coord->r - projection->R)) / cosC;
        tmp->y = ((cos(projection->D) * sin(coord->d)) -
                  (sin(projection->D) * cos(coord->d) * cos(coord->r - projection->R))) /
                 cosC;

    } else if (projection->type == PS_PROJ_SIN) {
        // psLogMsg(__func__, PS_LOG_WARN, "WARNING: psProject(): These projections don't work.");
        // return(NULL);

        tmp->x = cos(coord->d) * sin(coord->d);
        tmp->y = -cos(coord->d) * 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(DEG_TO_RAD(45.0) + (0.5 * coord->r)));
    } else if (projection->type == PS_PROJ_AIT) {
        psF64 tmpF64  = PS_SQRT_F32(0.5 * (1.0 + (cos(coord->r) * cos(0.5 * coord->d))));
        tmpF64 = 1.0 / tmpF64;
        tmp->x = 2.0 * tmpF64 * cos(coord->r) * sin(0.5 * coord->d);
        tmp->y = sin(coord->r) * tmpF64;

    } else if (projection->type == PS_PROJ_PAR) {
        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNDEFINED,
                "PS_PROJ_PAR");

    } else if (projection->type == PS_PROJ_GLS) {
        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNDEFINED,
                "PS_PROJ_GLS");
    } else {
        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNKNOWN,
                projection->type);
    }

    return (tmp);
}

/******************************************************************************
XXX: Waiting for algorithms.
XXX: Waiting for the definition of the PS_PROJ_PAR projection.
XXX: Waiting for the definition of the PS_PROJ_GLS projection.
XXX: Must apply scaling at the beginning.
 *****************************************************************************/
psSphere* psDeproject(const psPlane* coord,
                      const psProjection* projection)
{
    PS_PTR_CHECK_NULL(coord, NULL);
    PS_PTR_CHECK_NULL(projection, NULL);
    //    psLogMsg(__func__, PS_LOG_WARN, "WARNING: psDeproject(): These projections don't work.");
    //    return(NULL);
    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) {
        /*
                // ********************************************
                // From the ADD
                // ********************************************
                psF32 phi;
                psF32 theta;
                psF32 sinDelta;
                psF32 cosDeltaCosAlphaMinusAlphaP;
                psF32 cosDeltaSinAlphaMinusAlphaP;
         
                phi = atan2(-coord->y, coord->x);
                R = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y));
                theta = atan(1.0 / R);
         
                sinDelta = (sin(theta) * sin(projection->D)) + 
                           (cos(theta) * cos(projection->D) * cos(phi));
                cosDeltaCosAlphaMinusAlphaP = 
                           (sin(theta) * cos(projection->D)) - 
                           (cos(theta) * sin(projection->D) * cos(phi));
                cosDeltaSinAlphaMinusAlphaP = -cos(theta) * sin(phi);
                tmp->d = asin(sinDelta);
                tmp->r = atan2(cosDeltaSinAlphaMinusAlphaP, cosDeltaCosAlphaMinusAlphaP) + projection->R;
        */
        // ********************************************
        // From the mathworks: http://mathworld.wolfram.com/GnomonicProjection.html
        // delta_0 and phi_1 are the projection centers, not the point being projected.
        // (delta_0, phi_1) == (projection->R, projection->D)
        // XXX: figure out how to use the two-argument atan2() here.
        // ********************************************
        psF32 row = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y));
        psF32 C = atan(row);
        tmp->d = asin((cos(C) * sin(projection->D)) +
                      ((coord->y * sin(C) * cos(projection->D)) / row));
        psF32 tmpAtan = atan((coord->x * sin(C)) /
                             ((row * cos(projection->D) * cos(C)) -
                              (coord->y * sin(projection->D) * sin(C))));
        tmp->r = projection->R + tmpAtan;


    } else if (projection->type == PS_PROJ_SIN) {
        R = PS_SQRT_F32((coord->x * coord->x) + (coord->y * coord->y));
        tmp->d = atan2(-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 = PS_SQRT_F32(1.0 - chu1 - chu2);
        tmp->d = 2.0 * atan2((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) {
        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNDEFINED,
                "PS_PROJ_PAR");

    } else if (projection->type == PS_PROJ_GLS) {
        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNDEFINED,
                "PS_PROJ_GLS");
    } else {
        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNKNOWN,
                projection->type);
    }

    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_PTR_CHECK_NULL(position1, NULL);
    PS_PTR_CHECK_NULL(position2, NULL);

    if (position1->d >= DEG_TO_RAD(90.0)) {
        psLogMsg(__func__, PS_LOG_WARN, "WARNING: psDeproject(): position1->d is larger than 90 degrees.  Returning NULL.\n");
        return(NULL);
    }
    if (position2->d >= DEG_TO_RAD(90.0)) {
        psLogMsg(__func__, PS_LOG_WARN, "WARNING: psDeproject(): position2->d is larger than 90 degrees.  Returning NULL.\n");
        return(NULL);
    }

    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(position1, &proj);
        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 = RAD_TO_SEC(tmp->r);
            tmp->d = RAD_TO_SEC(tmp->d);
        } else if (unit == PS_ARCMIN) {
            tmp->r = RAD_TO_MIN(tmp->r);
            tmp->d = RAD_TO_MIN(tmp->d);
        } else if (unit == PS_DEGREE) {
            tmp->r = RAD_TO_DEG(tmp->r);
            tmp->d = RAD_TO_DEG(tmp->d);
        } else if (unit == PS_RADIAN) {}
        else {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                    PS_ERRORTEXT_psCoord_UNITS_UNKNOWN,
                    unit);
            psFree(tmp);
            return(NULL);
        }

        return (tmp);
    }

    psError(PS_ERR_BAD_PARAMETER_VALUE, true,
            PS_ERRORTEXT_psCoord_OFFSET_MODE_UNKNOWN,
            mode);
    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_PTR_CHECK_NULL(position, NULL);
    PS_PTR_CHECK_NULL(offset, NULL);

    psPlane lin;
    psSphere* tmp;
    psProjection proj;
    psF64 tmpR = 0.0;
    psF64 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 = SEC_TO_RAD(offset->r);
            tmpD = SEC_TO_RAD(offset->d);
        } else if (unit == PS_ARCMIN) {
            tmpR = MIN_TO_RAD(offset->r);
            tmpD = MIN_TO_RAD(offset->d);
        } else if (unit == PS_DEGREE) {
            tmpR = DEG_TO_RAD(offset->r);
            tmpD = DEG_TO_RAD(offset->d);
        } else if (unit == PS_RADIAN) {
            tmpR = offset->r;
            tmpD = offset->d;
        } else {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                    PS_ERRORTEXT_psCoord_UNITS_UNKNOWN,
                    unit);
            return(NULL);
        }

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

    psError(PS_ERR_BAD_PARAMETER_VALUE, true,
            PS_ERRORTEXT_psCoord_OFFSET_MODE_UNKNOWN,
            mode);

    return (NULL);
}



/******************************************************************************
psSpherePrecess(coords, fromTime, toTime):
 
XXX: Use static memory for tmpST.
 *****************************************************************************/
psSphere *psSpherePrecess(psSphere *coords,
                          const psTime *fromTime,
                          const psTime *toTime)
{
    psF64 fromMJD = psTimeToMJD(fromTime);
    psF64 toMJD = psTimeToMJD(toTime);
    psF64 T = (toMJD - fromMJD) / 36525.0;

    psF64 alphaP = DEG_TO_RAD(90.0) - ((DEG_TO_RAD(0.6406161) * T) +
                                       (DEG_TO_RAD(0.0000839) * T * T) +
                                       (DEG_TO_RAD(0.000005) * T * T * T));

    psF64 deltaP = (DEG_TO_RAD(0.5567530) * T) -
                   (DEG_TO_RAD(0.0001185) * T * T) -
                   (DEG_TO_RAD(0.0000116) * T * T * T);

    psF64 phiP = DEG_TO_RAD(90.0) - ((DEG_TO_RAD(0.6406161) * T) +
                                     (DEG_TO_RAD(0.0003041) * T * T) +
                                     (DEG_TO_RAD(0.0000051) * T * T * T));

    psSphereTransform *tmpST = psSphereTransformAlloc(alphaP, deltaP, phiP);

    psSphere *out = psSphereTransformApply(NULL, tmpST, coords);

    psFree(tmpST);

    return(out);
}

