Index: trunk/psLib/src/astro/psEarthOrientation.c =================================================================== --- trunk/psLib/src/astro/psEarthOrientation.c (revision 4541) +++ trunk/psLib/src/astro/psEarthOrientation.c (revision 4541) @@ -0,0 +1,79 @@ +/** @file psEarthOrientation.c +* +* @brief Function implementations for earth orientation calculations +* transformation +* +* @ingroup EarthOrientation +* +* @author Robert Daniel DeSonia, MHPCC +* +* @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-07-12 19:27:27 $ +* +* Copyright 2005 Maui High Performance Computing Center, University of Hawaii +*/ + +psSphere *psAberration(psSphere *apparent, const psSphere *actual, const psSphere *direction, double speed) +{ + return NULL; +} + +psSphere *psGravityDeflection(psSphere *apparent, psSphere *actual, psSphere *sun) +{ + return NULL; +} + + +double psEOC_ParallaxFactor(const psSphere *coords, const psTime *time) +{ + return NAN; +} + +psEarthPole *psEOC_PrecessionModel(const psTime *time) +{ + return NULL; +} + + +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; +} + + Index: trunk/psLib/src/astro/psEarthOrientation.h =================================================================== --- trunk/psLib/src/astro/psEarthOrientation.h (revision 4541) +++ trunk/psLib/src/astro/psEarthOrientation.h (revision 4541) @@ -0,0 +1,88 @@ +/** @file psEarthOrientation.h +* +* @brief Function prototypes for earth orientation calculations +* transformation +* +* @ingroup EarthOrientation +* +* @author Robert Daniel DeSonia, MHPCC +* +* @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-07-12 19:27:27 $ +* +* Copyright 2005 Maui High Performance Computing Center, University of Hawaii +*/ + +#ifndef PS_EARTH_ORIENTATION +#define PS_EARTH_ORIENTATION + +#include "psCoord.h" +#include "psTime.h" + +typedef struct +{ + double x; + double y; + double s; +} +psEarthPole; + +typedef enum { + PS_PRECESS_ROUGH, + PS_PRECESS_COMPLETE, + PS_PRECESS_IAU2000A, +} psPrecessMethod; + +psSphere *psAberration( + psSphere *apparent, + const psSphere *actual, + const psSphere *direction, + double speed +); + +psSphere *psGravityDeflection( + psSphere *apparent, + psSphere *actual, + psSphere *sun +); + +double psEOC_ParallaxFactor( + const psSphere *coords, + const psTime *time +); + +psEarthPole *psEOC_PrecessionModel( + const psTime *time +); + +psEarthPole *psEOC_PrecessionCorr( + const psTime *time, + psTimeBulletin bulletin +); + +psSphereRot *psSphereRot_CEOtoGCRS( + const psEarthPole *pole +); + +psSphereRot *psSphereRot_TEOtoCEO( + const psTime *time +); + +psEarthPole *psEOC_GetPolarMotion( + const psTime *time, + psTimeBulletin bulletin +); + +psEarthPole *psEOC_PolarTideCorr( + const psTime *time +); + +psEarthPole *psEOC_NutationCorr( + psTime *time +); + +psSphereRot *psSphereRot_ITRStoTEO( + const psEarthPole *motion +); + +#endif // #ifndef PS_EARTH_ORIENTATION Index: trunk/psLib/src/astro/psSphere.c =================================================================== --- trunk/psLib/src/astro/psSphere.c (revision 4541) +++ trunk/psLib/src/astro/psSphere.c (revision 4541) @@ -0,0 +1,230 @@ +/** @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.1 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-07-12 19:27:27 $ +* +* Copyright 2004-2005 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 "psAstrometry.h" +#include "psMatrix.h" +#include +#include + +static void sphereFree(psSphere *s) +{ + // There are non dynamic allocated items +} + +psSphere* psSphereAlloc(void) +{ + psSphere *s = psAlloc(sizeof(psSphere)); + + psMemSetDeallocator(s, (psFreeFcn) sphereFree); + return(s); +} + +psSphereRot* psSphereRotAlloc(double alphaP, + double deltaP, + double phiP) +{ + psSphereRot* rot = psAlloc(sizeof(psSphereRot)); + + double cosDelta = cos(deltaP); + double halfPhi = phiP / 2.0; + double sinHalfPhi = sin(halfPhi); + + // equations are directly from ADD + double vx = cosDelta*cos(alphaP); + double vy = cosDelta*sin(alphaP); + double vz = sin(deltaP); + + rot->q0 = vx*sinHalfPhi; + rot->q1 = vy*sinHalfPhi; + rot->q2 = vz*sinHalfPhi; + rot->q3 = cos(halfPhi); + + return rot; +} + +psSphereRot* psSphereRotQuat(double q0, + double q1, + double q2, + double q3) +{ + psSphereRot* rot = psAlloc(sizeof(psSphereRot)); + + double len = sqrt(q0*q0 + q1*q1 + q2*q2 + q3*q3); + rot->q0 = q0 / len; + rot->q1 = q1 / len; + rot->q2 = q2 / len; + rot->q3 = q3 / len; + + return rot; +} + +/****************************************************************************** +XXX: We convert Right Ascension angles to the range 0:PI. Is that acceptable? +XXX: Should we do something for Declination as well? + *****************************************************************************/ +psSphere* psSphereRotApply(psSphere* out, + const psSphereRot* transform, + const psSphere* coord) +{ + PS_ASSERT_PTR_NON_NULL(transform, NULL); + PS_ASSERT_PTR_NON_NULL(coord, NULL); + + if (out == NULL) { + out = psSphereAlloc(); + } + + + // apply the transform by creating a new psSphereRot from the input coord + // and combining it with the input transform (see ADD) + psSphereRot* coordRot = psSphereRotAlloc(coord->r, coord->d, 0); + coordRot->q3 = 0.0; + coordRot = psSphereRotCombine(coordRot, transform, coordRot); + // N.B., we can recycle coordRot right away due to the implementation of + // psSphereRotCombine puts the values of coordRot in a local variable first + + out->r = atan2(coordRot->q1,coordRot->q0); + out->d = atan2(coordRot->q2,sqrt(coordRot->q1*coordRot->q1+coordRot->q0*coordRot->q0)); + + return out; +} + +psSphereRot* psSphereRotCombine(psSphereRot* out, + const psSphereRot* rot1, + const psSphereRot* rot2) +{ + PS_ASSERT_PTR_NON_NULL(rot1, NULL); + PS_ASSERT_PTR_NON_NULL(rot2, NULL); + + if (out == NULL) { + out = (psSphereRot* ) psAlloc(sizeof(psSphereRot)); + } + + double a0 = rot1->q0; + double a1 = rot1->q1; + double a2 = rot1->q2; + double a3 = rot1->q3; + double b0 = rot2->q0; + double b1 = rot2->q1; + double b2 = rot2->q2; + double b3 = rot2->q3; + + // following came from ADD + out->q0 = b3*a0 + b2*a1 - b1*a2 + b0*a3; + out->q1 = b3*a1 - b2*a0 + b1*a3 + b0*a2; + out->q2 = b3*a2 + b2*a3 + b1*a0 - b0*a1; + out->q3 = b3*a3 - b3*a2 - b1*a1 - b0*a0; + + return out; +} + +psSphereRot *psSphereRotInvert(psSphereRot *rot) +{ +} + +psSphereTransform* psSphereTransformICRSToEcliptic(psTime *time) +{ + psF64 T; + + // Check for null parameter + PS_ASSERT_PTR_NON_NULL(time, NULL); + + // Convert psTime to MJD + psF64 MJD = psTimeToMJD(time); + + // Check the specified MJD is greater than 1900 + if ( MJD < MJD_1900 ) { + psError(PS_ERR_BAD_PARAMETER_TYPE,true,PS_ERRORTEXT_psCoord_INVALID_MJD); + return NULL; + } + + // Calculate number of Julian centuries since 1900 + T = ( MJD - MJD_1900 ) / JULIAN_CENTURY; + + psF64 alphaP = 0.0; + psF64 deltaP = DEG_TO_RAD(23.0) + + MIN_TO_RAD(27.0) + + SEC_TO_RAD(8.26) - + (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; + + // Don't neglect the minus sign on deltaP (bug 244): + return (psSphereTransformAlloc(alphaP, deltaP, phiP)); +} + + +psSphereTransform* psSphereTransformEclipticToICRS(psTime *time) +{ + psF64 T; + + // Check for null parameter + PS_ASSERT_PTR_NON_NULL(time, NULL); + + // Convert psTime to MJD + psF64 MJD = psTimeToMJD(time); + + // Check the specified MJD is greater than 1900 + if ( MJD < MJD_1900 ) { + psError(PS_ERR_BAD_PARAMETER_TYPE,true,PS_ERRORTEXT_psCoord_INVALID_MJD); + return NULL; + } + + // Calculate number of Julian centuries since 1900 + T = ( MJD - MJD_1900 ) / JULIAN_CENTURY; + + psF64 alphaP = 0.0; + psF64 deltaP = DEG_TO_RAD(23.0) + + MIN_TO_RAD(27.0) + + SEC_TO_RAD(8.26) - + (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; + + return (psSphereTransformAlloc(alphaP, -deltaP, phiP)); +} + +// XXX: This is bug 245: alphaP swaps with phiP from psSphereTransformGalacticToICRS() +psSphereTransform* psSphereTransformGalacticToICRS(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* psSphereTransformICRSToGalactic(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)); +} +