Index: /trunk/psLib/src/astro/psSphereOps.c =================================================================== --- /trunk/psLib/src/astro/psSphereOps.c (revision 4708) +++ /trunk/psLib/src/astro/psSphereOps.c (revision 4708) @@ -0,0 +1,457 @@ +/** @file psSphereOps.c + * + * @brief Contains spherical rotation and offset operations + * + * @ingroup CoordinateTransform + * + * @author Robert DeSonia, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-08-05 01:40:45 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ + +#include +#include + +#include "psSphereOps.h" +#include "psType.h" +#include "psCoord.h" +#include "psMemory.h" +#include "psTime.h" +#include "psConstants.h" +#include "psError.h" +#include "psLogMsg.h" +#include "psErrorText.h" + + +// Modified Julian Day 01/01/1900 00:00:00 +#define MJD_1900 15021.0 + +// Days in Julian century +#define JULIAN_CENTURY 36525.0 + + +psSphereRot* psSphereRotAlloc(double alphaP, + double deltaP, + double phiP) +{ + psSphereRot r,s,t; + + // directly from ADD -- there must be a better way?! + r.q0=0; + r.q1=0; + r.q2=sin(alphaP/2.0); + r.q3=cos(alphaP/2.0); + + s.q0=0; + s.q1=sin(deltaP/2.0); + s.q2=0; + s.q3=cos(deltaP/2.0); + + t.q0=0; + t.q1=0; + t.q2=sin(phiP/2.0); + t.q3=cos(phiP/2.0); + + // calculate t*s*r. + psSphereRot* result = psSphereRotCombine(NULL,&t,&s); + psSphereRotCombine(result,result,&r); + + return result; +} + +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) + double cosD = cos(coord->d); + psSphereRot* coordQuat = psSphereRotQuat( + cosD*cos(coord->r), + cosD*sin(coord->r), + sin(coord->d), + 0.0); + psSphereRot* coordQuatConjugate = psSphereRotQuat( + coordQuat->q0, coordQuat->q1, coordQuat->q2, coordQuat->q3); + coordQuat = psSphereRotInvert(coordQuat); + + // calculate q=(rp)r' + coordQuat = psSphereRotCombine(coordQuat, transform, coordQuat); + coordQuat = psSphereRotCombine(coordQuat, coordQuat, coordQuatConjugate); + // N.B., we can recycle coordQuat right away due to the implementation of + // psSphereRotCombine; it puts the input values in a local variable first + + out->r = atan2(coordQuat->q1,coordQuat->q0); + out->d = asin(coordQuat->q2); + + psFree(coordQuat); + psFree(coordQuatConjugate); + + 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) +{ + if (rot == NULL) { + return NULL; + // XXX: Error? + } + + rot->q0 = -rot->q0; + rot->q1 = -rot->q1; + rot->q2 = -rot->q2; + + return rot; +} + + +psSphereRot* psSphereRotEclipticToICRS(const 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 = DEG_TO_RAD(270.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 = DEG_TO_RAD(90.0); + + return (psSphereRotAlloc(alphaP, deltaP, phiP)); +} + +psSphereRot* psSphereRotICRSToEcliptic(const psTime *time) +{ + return psSphereRotInvert(psSphereRotEclipticToICRS(time)); +} + +// XXX: This is bug 245: alphaP swaps with phiP from psSphereTransformGalacticToICRS() +psSphereRot* psSphereRotGalacticToICRS(void) +{ + psF64 alphaP = DEG_TO_RAD(180.0-192.85948); + psF64 deltaP = DEG_TO_RAD(90.0-62.87175); + psF64 phiP = DEG_TO_RAD(90.0+32.93192); + + return (psSphereRotAlloc(alphaP, deltaP, phiP)); +} + +psSphereRot* psSphereRotICRSToGalactic(void) +{ + return psSphereRotInvert(psSphereRotGalacticToICRS()); +} + +/****************************************************************************** +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* position1, + const psSphere* position2, + psSphereOffsetMode mode, + psSphereOffsetUnit unit) +{ + PS_ASSERT_PTR_NON_NULL(position1, NULL); + PS_ASSERT_PTR_NON_NULL(position2, NULL); + + // Check positions near 90 degree and issue warnings if necessary + if (position1->d >= DEG_TO_RAD(90.0)) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: psDeproject(): position1->d is larger than 90 degrees. Returning NULL."); + 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."); + return NULL; + } + + // Allocate return structure + psSphere* tmp = psSphereAlloc(); + + // Mode is LINEAR - Use first position as projection center and project second point + // onto tangent plane, set point projected into psSphere structure x->r y->d + if (mode == PS_LINEAR) { + psProjection* proj = psProjectionAlloc(position1->r, + position1->d, + 1.0, + 1.0, + PS_PROJ_TAN); + + // Perform projection onto tangent plane + psPlane* lin = psProject(position2, proj); + + // Set return values + tmp->r = lin->x; + tmp->d = lin->y; + + // Free data structures allocated + psFree(proj); + psFree(lin); + + // Mode is SPHERICAL - Get difference between positiion 1 and position 2 and convert + // offset value from radians to desired units and return + } else if (mode == PS_SPHERICAL) { + 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; + + // Convert to desired units + 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; + } + // Invalid mode + } else { + + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psCoord_OFFSET_MODE_UNKNOWN, + mode); + psFree(tmp); + return NULL; + } + + // Return value + return tmp; +} + +/****************************************************************************** +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* position, + const psSphere* offset, + psSphereOffsetMode mode, + psSphereOffsetUnit unit) +{ + PS_ASSERT_PTR_NON_NULL(position, NULL); + PS_ASSERT_PTR_NON_NULL(offset, NULL); + + psSphere* tmp; + psF64 tmpR = 0.0; + psF64 tmpD = 0.0; + + // If mode is linear then set position to projection center + // and offset to linear coordinate then deproject to obtain + // new sphere coordinate + if (mode == PS_LINEAR) { + + // Allocate plane coordinate and set coordinate + psPlane* lin = psPlaneAlloc(); + lin->x = offset->r; + lin->y = offset->d; + + // Allocate and set projection structure + psProjection* proj = psProjectionAlloc(position->r, + position->d, + 1.0, + 1.0, + PS_PROJ_TAN); + + // Project tangent plane coord to spherical coord + tmp = psDeproject(lin, proj); + + // Free data structures used + psFree(proj); + psFree(lin); + + // If mode is spherical then convert offset to radians, add the offset + // to the position and wrap to 0 to 2pi + } else if (mode == PS_SPHERICAL) { + + // Convert offset unit to radians + 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; + } + + // Allocate sphere structure to return + tmp = psSphereAlloc(); + + // Add offset and wrap to 0 to 2PI if necessary + 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; + + // Invalid mode report error + } else { + + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psCoord_OFFSET_MODE_UNKNOWN, + mode); + return NULL; + } + + return tmp; +} + +/****************************************************************************** +psSpherePrecess(coords, fromTime, toTime): + +XXX: Use static memory for tmpST. + *****************************************************************************/ +psSphere *psSpherePrecess(psSphere *coords, + const psTime *fromTime, + const psTime *toTime) +{ + // Check input for NULL pointers + PS_ASSERT_PTR_NON_NULL(coords, NULL); + PS_ASSERT_PTR_NON_NULL(fromTime, NULL); + PS_ASSERT_PTR_NON_NULL(toTime, NULL); + + // Calculate Julian centuries + psF64 fromMJD = psTimeToMJD(fromTime); + psF64 toMJD = psTimeToMJD(toTime); + psF64 T = (toMJD - fromMJD) / JULIAN_CENTURY; + + // Calculate conversion constants + 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)); + + // Create transform with proper constants + psSphereRot* tmpST = psSphereRotAlloc(alphaP, deltaP, phiP); + + // Apply transform to coordinates + psSphere *out = psSphereRotApply(NULL, tmpST, coords); + + psFree(tmpST); + + return(out); +} + Index: /trunk/psLib/src/astro/psSphereOps.h =================================================================== --- /trunk/psLib/src/astro/psSphereOps.h (revision 4708) +++ /trunk/psLib/src/astro/psSphereOps.h (revision 4708) @@ -0,0 +1,217 @@ +/** @file psSphereOps.h + * + * @brief Contains spherical rotation and offset operations + * + * @ingroup CoordinateTransform + * + * @author Robert DeSonia, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-08-05 01:40:45 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ + +#ifndef PS_SPHERE_H +#define PS_SPHERE_H + +/// @addtogroup CoordinateTransform +/// @{ + +#include "psCoord.h" + +/** Spherical Rotation 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. + * + */ +typedef struct +{ + double q0; + double q1; + double q2; + double q3; +} +psSphereRot; + +/** Mode for Offset calculation between two sky positions + * + * @see psSphereGetOffset, psSphereSetOffset + * + */ +typedef enum { + PS_SPHERICAL, ///< offset corresponds to an angular offset + PS_LINEAR ///< offset corresponds to a linear offset +} psSphereOffsetMode; + +/** The units of the offset + * + * @see psSphereGetOffset, psSphereSetOffset + * + */ +typedef enum { + PS_ARCSEC, ///< Arcseconds + PS_ARCMIN, ///< Arcminutes + PS_DEGREE, ///< Degrees + PS_RADIAN ///< Radians +} psSphereOffsetUnit; + + +/** Allocator for psSphereRot + * + * @return psSphereRot* newly allocated psSphereRot + */ + +psSphereRot* psSphereRotAlloc( + 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). +); + +/** Allocator for psSphereRot given quaternions. + * + * Normalizes a given quaternion and returns the cooresponding newly allocated + * psSphereRot object. + * + * @return psSphereRot* newly allocated psSphereRot + */ +psSphereRot* psSphereRotQuat( + double q0, ///< q0 + double q1, ///< q1 + double q2, ///< q2 + double q3 ///< q3 +); + +/** Applies the psSphereTransform transform for a specified coordinate + * + * @return psSphere* resulting coordinate based on transform + */ +psSphere* psSphereRotApply( + psSphere* out, ///< a psSphere to recycle. If NULL, a new one is generated. + const psSphereRot* transform, ///< the transform to apply + const psSphere* coord ///< the coordinate to apply the transform above.x +); + +/** Combines two rotations. + * + * Combines two rotations to produce a single rotation which is the + * equivalent of applying the first rotation and then the second. The output + * rotation may be supplied, or will be allocated if NULL. + * + * @return psSphereRot* New psSphereRot that is the combination of the input psSphereRots + */ +psSphereRot* psSphereRotCombine( + psSphereRot* out, ///< a psSphereRot to recycle or NULL + const psSphereRot* rot1, ///< first rotation to combine + const psSphereRot* rot2 ///< second rotation to combine +); + +/** Inverts the rotation. + * + * @return psSphereRot* Inverted input psSphereRot + */ +psSphereRot* psSphereRotInvert( + psSphereRot *rot +); + +/** Converts a psCube to a psSphere + * + * @return psSphere* New psSphere that is equivalent to the input psCube + * + */ +psSphere* psCubeToSphere( + const psCube* cube ///< a cubic coordinate to convert +); + +/** Converts a psSphere to a psCube + * + * @return psCube* New psCube that is equivalent to the input psSphere + * + */ +psCube* psSphereToCube( + const psSphere* sphere ///< a spherical coordinate to convert +); + +/** Determines the offset (RA,Dec) on the sky between two positions. + * + * Both an offset mode and an offset unit may be defined. The mode may be + * either PS_SPHERICAL, in which case the specified offset corresponds to an + * offset in angles, or it may be PS_LINEAR, in which case the offset + * corresponds to a linear offset in a local projection. The offset unit may + * be in one of PS_ARCSEC, PS_ARCMIN, PS_DEGREE, and PS_RADIAN, which + * specifies the units of the offset only. + * + * @return psSphere* the offset between position1 and position2 + */ +psSphere* psSphereGetOffset( + const psSphere* position1, ///< first position for calculating offset + const psSphere* position2, ///< second position for calculating offset + psSphereOffsetMode mode, ///< type of offset can be PS_SPHERICAL or PS_LINEAR + psSphereOffsetUnit unit ///< specifies the units of offset only +); + +/** Applies the given offset to a coordinate. + * + * Both an offset mode and an offset unit may be defined. The mode may be + * either PS_SPHERICAL, in which case the specified offset corresponds to an + * offset in angles, or it may be PS_LINEAR, in which case the offset + * corresponds to a linear offset in a local projection. The offset unit may + * be in one of PS_ARCSEC, PS_ARCMIN, PS_DEGREE, and PS_RADIAN, which + * specifies the units of the offset only. + * + * @return psSphere* the original position with the given offset applied. + */ +psSphere* psSphereSetOffset( + const psSphere* position, ///< coordinate of origin + const psSphere* offset, ///< coordinate of offset to apply + psSphereOffsetMode mode, ///< corresponds to an offset in angles or local projection + psSphereOffsetUnit unit ///< specifies the units of offset only +); + +/** Creates the appropriate transform for converting from ICRS to Ecliptic + * coordinate systems. + * + * @return psSphereTransform* transform for ICRS->Ecliptic coordinate systems + */ +psSphereRot* psSphereRotICRSToEcliptic( + const psTime *time ///< the time for which the resulting transform will be valid +); + +/** Creates the appropriate transform for converting from Ecliptic to ICRS + * coordinate systems. + * + * @return psSphereTransform* transform for Ecliptic->ICRS coordinate systems + */ +psSphereRot* psSphereRotEclipticToICRS( + const psTime *time ///< the time for which the resulting transform will be valid +); + +/** Creates the appropriate transform for converting from ICRS to Galactic + * coordinate systems. + * + */ +psSphereRot* psSphereRotICRSToGalactic(void); + +/** Creates the appropriate transform for converting from Galactic to ICRS + * coordinate systems. + * + */ +psSphereRot* psSphereRotGalacticToICRS(void); + +/** Generates the complete spherical rotation to account for precession + * between two times. The equinoxes shall be Julian equinoxes. + * + * @return psSphere* the resulting spherical rotation + */ +psSphere* psSpherePrecess( + psSphere *coords, ///< coordinates (modified in-place) + const psTime *fromTime, ///< equinox of coords input + const psTime *toTime ///< equinox of coords output +); + +/// @} + +#endif // #ifndef PS_SPHERE_H