Index: trunk/psLib/src/astronomy/psCoord.c =================================================================== --- trunk/psLib/src/astronomy/psCoord.c (revision 4409) +++ trunk/psLib/src/astronomy/psCoord.c (revision 4540) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.78 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-06-27 20:38:11 $ +* @version $Revision: 1.79 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-07-12 19:12:00 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -30,7 +30,4 @@ #include #include -/******************************************************************************/ -/* DEFINE STATEMENTS */ -/******************************************************************************/ // Modified Julian Day 01/01/1900 00:00:00 @@ -40,32 +37,29 @@ #define JULIAN_CENTURY 36525.0 -/******************************************************************************/ -/* TYPE DEFINITIONS */ -/******************************************************************************/ - -// None - -/*****************************************************************************/ -/* GLOBAL VARIABLES */ -/*****************************************************************************/ - -// None - -/*****************************************************************************/ -/* FILE STATIC VARIABLES */ -/*****************************************************************************/ - -// None - -/*****************************************************************************/ -/* FUNCTION IMPLEMENTATION - LOCAL */ -/*****************************************************************************/ - -/*****************************************************************************/ -/* FUNCTION IMPLEMENTATION - PUBLIC */ -/*****************************************************************************/ static void planeFree(psPlane *p) { // There are non dynamic allocated items +} + +static void sphereFree(psSphere *s) +{ + // There are non dynamic allocated items +} + +static void planeTransformFree(psPlaneTransform *pt) +{ + psFree(pt->x); + psFree(pt->y); +} + +static void planeDistortFree(psPlaneDistort *pt) +{ + psFree(pt->x); + psFree(pt->y); +} + +static void projectionFree(psProjection *p) +{ + // There are no dynamically allocated items } @@ -201,10 +195,4 @@ } - -static void sphereFree(psSphere *s) -{ - // There are non dynamic allocated items -} - psSphere* psSphereAlloc(void) { @@ -215,8 +203,41 @@ } -static void planeTransformFree(psPlaneTransform *pt) -{ - psFree(pt->x); - psFree(pt->y); +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; } @@ -260,10 +281,4 @@ return (out); -} - -static void planeDistortFree(psPlaneDistort *pt) -{ - psFree(pt->x); - psFree(pt->y); } @@ -319,51 +334,10 @@ /****************************************************************************** -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:2*PI. - *****************************************************************************/ -psF32 piNormalize(psF32 angle) -{ - while (angle < FLT_EPSILON) { - angle+=M_PI*2; - } - - while (angle >= (M_PI*2)) { - angle-=M_PI*2; - } - 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) +psSphere* psSphereRotApply(psSphere* out, + const psSphereRot* transform, + const psSphere* coord) { PS_ASSERT_PTR_NON_NULL(transform, NULL); @@ -371,26 +345,64 @@ 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->phiP; - out->r = piNormalize(phi); - out->d = theta; - - return(out); -} - -psSphereTransform* psSphereTransformICRSToEcliptic(psTime *time) + 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) +{ + PS_ASSERT_PTR_NON_NULL(rot, NULL); + + double norm = sqrt(rot->q0*rot->q0 + rot->q1*rot->q1 + rot->q2*rot->q2 + rot->q3*rot->q3); + rot->q1 = -rot->q1 / norm; + rot->q2 = -rot->q2 / norm; + rot->q3 = -rot->q3 / norm; + + return rot; +} + +psSphereRot* psSphereRotICRSToEcliptic(const psTime *time) { psF64 T; @@ -421,9 +433,9 @@ // Don't neglect the minus sign on deltaP (bug 244): - return (psSphereTransformAlloc(alphaP, deltaP, phiP)); -} - - -psSphereTransform* psSphereTransformEclipticToICRS(psTime *time) + return (psSphereRotAlloc(alphaP, deltaP, phiP)); +} + + +psSphereRot* psSphereRotEclipticToICRS(const psTime *time) { psF64 T; @@ -453,9 +465,9 @@ psF64 phiP = 0.0; - return (psSphereTransformAlloc(alphaP, -deltaP, phiP)); + return (psSphereRotAlloc(alphaP, -deltaP, phiP)); } // XXX: This is bug 245: alphaP swaps with phiP from psSphereTransformGalacticToICRS() -psSphereTransform* psSphereTransformGalacticToICRS(void) +psSphereRot* psSphereRotGalacticToICRS(void) { psF64 alphaP = DEG_TO_RAD(32.93192); @@ -463,8 +475,8 @@ psF64 phiP = DEG_TO_RAD(282.85948); - return (psSphereTransformAlloc(alphaP, deltaP, phiP)); -} - -psSphereTransform* psSphereTransformICRSToGalactic(void) + return (psSphereRotAlloc(alphaP, deltaP, phiP)); +} + +psSphereRot* psSphereRotICRSToGalactic(void) { psF64 alphaP = DEG_TO_RAD(282.85948); @@ -472,10 +484,5 @@ psF64 phiP = DEG_TO_RAD(32.93192); - return (psSphereTransformAlloc(alphaP, deltaP, phiP)); -} - -void projectionFree(psProjection *p) -{ - // There are no dynamically allocated items + return (psSphereRotAlloc(alphaP, deltaP, phiP)); } @@ -837,8 +844,8 @@ // Create transform with proper constants - psSphereTransform *tmpST = psSphereTransformAlloc(alphaP, deltaP, phiP); + psSphereRot* tmpST = psSphereRotAlloc(alphaP, deltaP, phiP); // Apply transform to coordinates - psSphere *out = psSphereTransformApply(NULL, tmpST, coords); + psSphere *out = psSphereRotApply(NULL, tmpST, coords); psFree(tmpST); Index: trunk/psLib/src/astronomy/psCoord.h =================================================================== --- trunk/psLib/src/astronomy/psCoord.h (revision 4409) +++ trunk/psLib/src/astronomy/psCoord.h (revision 4540) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.37 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-06-27 20:38:11 $ +* @version $Revision: 1.38 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-07-12 19:12:00 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -24,5 +24,5 @@ #include "psList.h" #include "psFunctions.h" -#include "psTime.h" +// N.B. inclusion of psTime.h was done to after the typedefs to handle cross-dependency of typedefs /// @addtogroup CoordinateTransform @@ -38,8 +38,8 @@ typedef struct { - double x; ///< x position - double y; ///< y position - double xErr; ///< Error in x position - double yErr; ///< Error in y position + double x; ///< x position + double y; ///< y position + double xErr; ///< Error in x position + double yErr; ///< Error in y position } psPlane; @@ -52,12 +52,40 @@ * */ -typedef struct psSphere -{ - double r; ///< RA - double d; ///< Dec - double rErr; ///< Error in RA - double dErr; ///< Error in Dec +typedef struct +{ + double r; ///< RA + double d; ///< Dec + double rErr; ///< Error in RA + double dErr; ///< Error in Dec } psSphere; + +/** Cubic Coordinate System + * + */ +typedef struct +{ + double x; ///< cos (DEC) cos (RA) + double y; ///< cos (DEC) sic (RA) + double z; ///< sin (DEC) + double xErr; ///< Error in x + double yErr; ///< Error in y + double zErr; ///< Error in z +} +psCube; + +/** Spherical rotations represent coordinate transformation in 3-D, as well as + * the effects of precession and nutation. The structure contains the + * elements of a quaternion to represent the spherical rotational. + * + */ +typedef struct +{ + double q0; ///< first element of the quaternion + double q1; ///< second element of the quaternion + double q2; ///< third element of the quaternion + double q3; ///< fourth element of the quaternion +} +psSphereRot; /** 2D Polynomial Transform @@ -94,23 +122,4 @@ } psPlaneDistort; - -/** Spherical Transform Definition - * - * We need to be able to convert between ICRS, Galactic and Ecliptic - * coordinates, and potentially between arbitrary spherical coordinate - * systems. All of these basic spherical transformations represent rotations - * of the spherical coordinate reference. We specify a general - * transformation function which takes a structure, psSphereTransform, - * defining the transformation between two spherical coordinate systems - * - */ -typedef struct -{ - double alphaP; ///< Longitude of the target system pole in the source system - double cosDeltaP; ///< Cosine of target pole latitude in the source system - double sinDeltaP; ///< Sine of target pole latitude in the source system - double phiP; ///< Longitude of the ascending node in the target system -} -psSphereTransform; /** Projection type for projection/deprojection @@ -167,4 +176,6 @@ } psSphereOffsetUnit; +#include "psTime.h" + /** Allocates a psPlane * @@ -178,7 +189,28 @@ * @return psSphere* resulting sphere structure. */ - psSphere* psSphereAlloc(void); +/** psSphereRot allocator which defines the rotation in terms of the coordinate + * of the pole and the rotation about that pole. + * + * @return psSphereRot* Newly allocated psSphereRot object + */ +psSphereRot* psSphereRotAlloc( + double alphaP, + double deltaP, + double phiP +); + +/** psSphereRot allocator which defines the rotation from the elements of the + * quaternion. + * + * @return psSphereRot* Newly allocated psSphereRot object + */ +psSphereRot* psSphereRotQuat( + double q0, + double q1, + double q2, + double q3 +); /** Allocates a psPlaneTransform transform. @@ -186,5 +218,4 @@ * @return psPlaneTransform* resulting plane transform */ - psPlaneTransform* psPlaneTransformAlloc( int n1, ///< The order of the x term in the transform. @@ -227,23 +258,32 @@ ); -/** Allocator for psSphereTransform - * - * @return psSphereTransform* newly allocated struct - */ - -psSphereTransform* psSphereTransformAlloc( - double alphaP, ///< north pole latitude - double deltaP, ///< north pole longitude? - double phiP ///< defines the longitude in the input system of the equatorial intersection between the two systems (e.g, the first point of Ares). -); - -/** Applies the psSphereTransform transform for a specified coordinate + +/** Applies the psSphereRot transform for a specified coordinate * * @return psSphere* resulting coordinate based on transform */ -psSphere* psSphereTransformApply( +psSphere* psSphereRotApply( psSphere* out, ///< a psSphere to recycle. If NULL, a new one is generated. - const psSphereTransform* transform,///< the transform to apply + const psSphereRot* transform, ///< the transform to apply const psSphere* coord ///< the coordinate to apply the transform above.x +); + +/** Combines two rotations to produce a single rotation which is equivalent of + * applying the first rotation and then the second. + * + * @return psSphereRot* new psSphereRot transform + */ +psSphereRot* psSphereRotCombine( + psSphereRot* out, + const psSphereRot* rot1, + const psSphereRot* rot2 +); + +/** Inverts a psSphereRot's rotation. + * + * @return psSphereRot* The inverted psSphereRot + */ +psSphereRot* psSphereRotInvert( + psSphereRot* rot ///< the psSphereRot to invert ); @@ -251,8 +291,8 @@ * coordinate systems. * - * @return psSphereTransform* transform for ICRS->Ecliptic coordinate systems - */ -psSphereTransform* psSphereTransformICRSToEcliptic( - psTime *time ///< the time for which the resulting transform will be valid + * @return psSphereRot* transform for ICRS->Ecliptic coordinate systems + */ +psSphereRot* psSphereRotICRSToEcliptic( + const psTime* time ///< the time for which the resulting transform will be valid ); @@ -260,8 +300,8 @@ * coordinate systems. * - * @return psSphereTransform* transform for Ecliptic->ICRS coordinate systems - */ -psSphereTransform* psSphereTransformEclipticToICRS( - psTime *time ///< the time for which the resulting transform will be valid + * @return psSphereRot* transform for Ecliptic->ICRS coordinate systems + */ +psSphereRot* psSphereRotEclipticToICRS( + const psTime* time ///< the time for which the resulting transform will be valid ); @@ -270,5 +310,5 @@ * */ -psSphereTransform* psSphereTransformICRSToGalactic(void); +psSphereRot* psSphereRotICRSToGalactic(void); /** Creates the appropriate transform for converting from Galactic to ICRS @@ -276,5 +316,5 @@ * */ -psSphereTransform* psSphereTransformGalacticToICRS(void); +psSphereRot* psSphereRotGalacticToICRS(void); /** Allocates memory for a psProjection structure Index: trunk/psLib/src/astronomy/psTime.c =================================================================== --- trunk/psLib/src/astronomy/psTime.c (revision 4409) +++ trunk/psLib/src/astronomy/psTime.c (revision 4540) @@ -10,6 +10,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.65 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-06-28 20:17:52 $ + * @version $Revision: 1.66 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-07-12 19:12:00 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -887,5 +887,5 @@ } -struct psSphere* p_psTimeGetPoleCoords(const psTime* time) +psSphere* p_psTimeGetPoleCoords(const psTime* time) { psU32 nTables = 3; @@ -896,5 +896,5 @@ psF64 c = 0.0; psF64 mjdPred = 0.0; - struct psSphere* output = NULL; + psSphere* output = NULL; psLookupStatusType xStatus = PS_LOOKUP_SUCCESS; psLookupStatusType yStatus = PS_LOOKUP_SUCCESS; @@ -908,6 +908,5 @@ PS_ASSERT_INT_WITHIN_RANGE(time->nsec,0,(psU32)((1e9)-1),NULL); - if(time->type != PS_TIME_TAI) - { + if(time->type != PS_TIME_TAI) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_TYPE_INCORRECT, time->type); return NULL; @@ -922,6 +921,5 @@ // Value could not be found through table lookup and interpolation - if(xStatus==PS_LOOKUP_PAST_TOP && yStatus==PS_LOOKUP_PAST_TOP) - { + if(xStatus==PS_LOOKUP_PAST_TOP && yStatus==PS_LOOKUP_PAST_TOP) { // Date too earlier for tables. Get default polar coodinate values from metadata, and issue warning. @@ -942,6 +940,5 @@ y = tableMetadataItem->data.F64; - } else if(xStatus==PS_LOOKUP_PAST_BOTTOM && yStatus==PS_LOOKUP_PAST_BOTTOM) - { + } else if(xStatus==PS_LOOKUP_PAST_BOTTOM && yStatus==PS_LOOKUP_PAST_BOTTOM) { /* Date too late for tables. Issue warning and use following formulae for predicting @@ -1000,6 +997,5 @@ yp->data.F64[4]*sin(c); - } else if(xStatus!=PS_LOOKUP_SUCCESS || yStatus!=PS_LOOKUP_SUCCESS) - { + } else if(xStatus!=PS_LOOKUP_SUCCESS || yStatus!=PS_LOOKUP_SUCCESS) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_INTERPOLATION_FAILED); return NULL; Index: trunk/psLib/src/astronomy/psTime.h =================================================================== --- trunk/psLib/src/astronomy/psTime.h (revision 4409) +++ trunk/psLib/src/astronomy/psTime.h (revision 4540) @@ -11,6 +11,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.33 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-06-28 20:17:52 $ + * @version $Revision: 1.34 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-07-12 19:12:00 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -25,11 +25,8 @@ #include "psType.h" -#include "psImage.h" - -struct psSphere; +// N.B. inclusion of psCoord.h was done to after the typedefs to handle cross-dependency of typedefs /// @addtogroup Time /// @{ - /** Time type. @@ -69,4 +66,6 @@ psTime; +#include "psCoord.h" +#include "psImage.h" /** Initialize time data. @@ -161,7 +160,7 @@ * @return psSphere*: Spherical coordinates of Earth's polar axias. */ -struct psSphere* p_psTimeGetPoleCoords( - const psTime *time ///< psTime determine polar orientation. - ); +psSphere* p_psTimeGetPoleCoords( + const psTime *time ///< psTime determine polar orientation. +); /** Calculate the number of leapseconds between two times.