Index: trunk/psLib/src/astronomy/psCoord.c =================================================================== --- trunk/psLib/src/astronomy/psCoord.c (revision 4401) +++ 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);