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