Index: /trunk/psLib/src/astro/Makefile.am =================================================================== --- /trunk/psLib/src/astro/Makefile.am (revision 4619) +++ /trunk/psLib/src/astro/Makefile.am (revision 4620) @@ -6,5 +6,6 @@ libpslibastro_la_SOURCES = \ psTime.c \ - psCoord.c + psCoord.c \ + psSphereOps.c EXTRA_DIST = astro.i @@ -13,3 +14,4 @@ pslibinclude_HEADERS = \ psTime.h \ - psCoord.h + psCoord.h \ + psSphereOps.h Index: /trunk/psLib/src/astro/psCoord.c =================================================================== --- /trunk/psLib/src/astro/psCoord.c (revision 4619) +++ /trunk/psLib/src/astro/psCoord.c (revision 4620) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.81 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-07-27 00:36:01 $ +* @version $Revision: 1.82 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-07-27 19:55:15 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -29,42 +29,25 @@ #include #include -/******************************************************************************/ -/* DEFINE STATEMENTS */ -/******************************************************************************/ - -// Modified Julian Day 01/01/1900 00:00:00 -#define MJD_1900 15021.0 - -// Days in Julian century -#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 cubeFree(psCube* c) +{ + // There are non dynamic allocated items +} + +static void planeTransformFree(psPlaneTransform *pt) +{ + psFree(pt->x); + psFree(pt->y); } @@ -195,26 +178,16 @@ { psPlane *p = psAlloc(sizeof(psPlane)); - psMemSetDeallocator(p, (psFreeFunc) planeFree); + return(p); } -static void sphereFree(psSphere *s) -{ - // There are non dynamic allocated items -} - psSphere* psSphereAlloc(void) { psSphere *s = psAlloc(sizeof(psSphere)); - psMemSetDeallocator(s, (psFreeFunc) sphereFree); + return(s); -} - -static void cubeFree(psCube *c) -{ - // There are non dynamic allocated items } @@ -225,10 +198,4 @@ psMemSetDeallocator(c, (psFreeFunc) cubeFree); return(c); -} - -static void planeTransformFree(psPlaneTransform *pt) -{ - psFree(pt->x); - psFree(pt->y); } @@ -330,27 +297,4 @@ } -/****************************************************************************** -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); -} /****************************************************************************** @@ -369,120 +313,4 @@ } 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_ASSERT_PTR_NON_NULL(transform, NULL); - PS_ASSERT_PTR_NON_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->phiP; - out->r = piNormalize(phi); - out->d = theta; - - return(out); -} - -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)); } @@ -617,246 +445,4 @@ } -/****************************************************************************** -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 - psSphereTransform *tmpST = psSphereTransformAlloc(alphaP, deltaP, phiP); - - // Apply transform to coordinates - psSphere *out = psSphereTransformApply(NULL, tmpST, coords); - - psFree(tmpST); - - return(out); -} - /***************************************************************************** multiplyDPoly2D(trans1, trans2): Takes two 2-D polynomials as input and Index: /trunk/psLib/src/astro/psCoord.h =================================================================== --- /trunk/psLib/src/astro/psCoord.h (revision 4619) +++ /trunk/psLib/src/astro/psCoord.h (revision 4620) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.39 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-07-27 00:36:01 $ +* @version $Revision: 1.40 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-07-27 19:55:15 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -112,23 +112,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 @@ -163,31 +144,8 @@ psProjection; -/** 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; - /** Allocates a psPlane * * @return psPlane* resulting plane structure. */ - psPlane* psPlaneAlloc(void); @@ -196,6 +154,11 @@ * @return psSphere* resulting sphere structure. */ - psSphere* psSphereAlloc(void); + +/** Allocates a psCube + * + * @return psCube* resulting cubic structure. + */ +psCube* psCubeAlloc(void); /** Allocates a psCube @@ -250,55 +213,4 @@ ); -/** 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 - * - * @return psSphere* resulting coordinate based on transform - */ -psSphere* psSphereTransformApply( - psSphere* out, ///< a psSphere to recycle. If NULL, a new one is generated. - const psSphereTransform* transform,///< the transform to apply - const psSphere* coord ///< the coordinate to apply the transform above.x -); - -/** Creates the appropriate transform for converting from ICRS to Ecliptic - * coordinate systems. - * - * @return psSphereTransform* transform for ICRS->Ecliptic coordinate systems - */ -psSphereTransform* psSphereTransformICRSToEcliptic( - 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 - */ -psSphereTransform* psSphereTransformEclipticToICRS( - psTime *time ///< the time for which the resulting transform will be valid -); - -/** Creates the appropriate transform for converting from ICRS to Galactic - * coordinate systems. - * - */ -psSphereTransform* psSphereTransformICRSToGalactic(void); - -/** Creates the appropriate transform for converting from Galactic to ICRS - * coordinate systems. - * - */ -psSphereTransform* psSphereTransformGalacticToICRS(void); - /** Allocates memory for a psProjection structure * @@ -329,51 +241,4 @@ const psPlane* coord, ///< coordinate to project const psProjection* projection ///< parameters of the projection -); - -/** 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 -); - -/** 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 ); Index: unk/psLib/src/astro/psSphere.c =================================================================== --- /trunk/psLib/src/astro/psSphere.c (revision 4619) +++ (revision ) @@ -1,230 +1,0 @@ -/** @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)); -} - Index: /trunk/psLib/src/math/psStats.c =================================================================== --- /trunk/psLib/src/math/psStats.c (revision 4619) +++ /trunk/psLib/src/math/psStats.c (revision 4620) @@ -14,6 +14,6 @@ * stats->binsize * - * @version $Revision: 1.137 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-07-12 19:12:01 $ + * @version $Revision: 1.138 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-07-27 19:55:16 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii Index: /trunk/psLib/src/pslib_strict.h =================================================================== --- /trunk/psLib/src/pslib_strict.h (revision 4619) +++ /trunk/psLib/src/pslib_strict.h (revision 4620) @@ -9,6 +9,6 @@ * @author Eric Van Alst, MHPCC * -* @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-07-12 19:12:00 $ +* @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-07-27 19:55:15 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -30,4 +30,5 @@ #include "psTime.h" #include "psCoord.h" +#include "psSphereOps.h" /// @}