Index: trunk/psLib/src/astro/psEarthOrientation.c =================================================================== --- trunk/psLib/src/astro/psEarthOrientation.c (revision 6039) +++ trunk/psLib/src/astro/psEarthOrientation.c (revision 6184) @@ -8,6 +8,6 @@ * @author Robert Daniel DeSonia, MHPCC * - * @version $Revision: 1.29 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-01-18 23:49:06 $ + * @version $Revision: 1.30 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-01-23 20:04:31 $ * * Copyright 2005 Maui High Performance Computing Center, University of Hawaii @@ -1462,49 +1462,175 @@ } -/****************************************************************************** -psSpherePrecess(coords, fromTime, toTime): - -XXX: Use static memory for tmpST. - *****************************************************************************/ -psSphere *psSpherePrecess(psSphere *coords, - const psTime *fromTime, - const psTime *toTime) +psSphereRot *psSpherePrecess(const psTime *fromTime, + const psTime *toTime, + psPrecessMethod mode) { // 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); - + if (fromTime == NULL && toTime == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + "Invalid time inputs. fromTime & toTime cannot both be NULL.\n"); + return NULL; + } + PS_ASSERT_INT_WITHIN_RANGE(mode, PS_PRECESS_ROUGH, PS_PRECESS_IAU2000A, NULL); + psF64 fromMJD, toMJD; // 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) + - psF64 alphaP = DEG_TO_RAD(180.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) + - psF64 phiP = DEG_TO_RAD(180.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); - + //If either input time is NULL, assume it to be J2000 -> from SDRS as of rev 18 + if (fromTime == NULL) { + fromMJD = MJD_2000; + } else { + fromMJD = psTimeToMJD(fromTime); + } + if (toTime == NULL) { + toMJD = MJD_2000; + } else { + toMJD = psTimeToMJD(toTime); + } + psTime *from = NULL; + psTime *to = NULL; + + if (mode == PS_PRECESS_ROUGH) { + //For PS_PRECESS_ROUGH, no time/earthpole corrections are used. This is the + //lowest level of detail mode. + psF64 T = (toMJD - fromMJD) / JULIAN_CENTURY; + + // Calculate conversion constants + // psF64 alphaP = DEG_TO_RAD(90.0) - ((DEG_TO_RAD(0.6406161) * T) + + psF64 alphaP = DEG_TO_RAD(180.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) + + psF64 phiP = DEG_TO_RAD(180.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); + return tmpST; + } else if (mode == PS_PRECESS_IAU2000A) { + //For IAU2000A mode, run psEOC_PrecessionModel and then psSphereRot_CEOtoGCRS for + //each time. Then difference the resulting rotations by adding the inverse + //rotation corresponding to fromTime to the toTime rotation. + + //Since the time inputs are allowed to be NULL, either convert the MJD time + //or copy to non-NULL time + if (fromTime == NULL) { + from = psTimeFromMJD(fromMJD); + } else { + from = p_psTimeCopy(fromTime); + } + if (toTime == NULL) { + to = psTimeFromMJD(toMJD); + } else { + to = p_psTimeCopy(toTime); + } + + psEarthPole *fromEP = psEOC_PrecessionModel(from); + psSphereRot *fromRot = psSphereRot_CEOtoGCRS(fromEP); + psEarthPole *toEP = psEOC_PrecessionModel(to); + psSphereRot *toRot = psSphereRot_CEOtoGCRS(toEP); + psSphereRot *fromConj = psSphereRotConjugate(NULL, fromRot); + psSphereRot *out = psSphereRotCombine(NULL, toRot, fromConj); + psFree(from); + psFree(to); + psFree(fromEP); + psFree(fromRot); + psFree(toEP); + psFree(toRot); + psFree(fromConj); + return out; + } else if (mode == PS_PRECESS_COMPLETE_A) { + //For PS_PRECESS_COMPLETE_A the same procedure as IAU2000A is used but with + //additional earthpole corrections from psEOC_PrecessionCorr. The corrections + //for COMPLETE_A come from the IERS Bulletin A. + + //Since the time inputs are allowed to be NULL, either convert the MJD time + //or copy to non-NULL time + if (fromTime == NULL) { + from = psTimeFromMJD(fromMJD); + } else { + from = p_psTimeCopy(fromTime); + } + if (toTime == NULL) { + to = psTimeFromMJD(toMJD); + } else { + to = p_psTimeCopy(toTime); + } + + psEarthPole *fromEP = psEOC_PrecessionModel(from); + psEarthPole *fromCorr = psEOC_PrecessionCorr(from, PS_IERS_A); + fromEP->x += fromCorr->x; + fromEP->y += fromCorr->y; + fromEP->s += fromCorr->s; + psSphereRot *fromRot = psSphereRot_CEOtoGCRS(fromEP); + psEarthPole *toEP = psEOC_PrecessionModel(to); + psEarthPole *toCorr = psEOC_PrecessionCorr(to, PS_IERS_A); + toEP->x += toCorr->x; + toEP->y += toCorr->y; + toEP->s += toCorr->s; + psSphereRot *toRot = psSphereRot_CEOtoGCRS(toEP); + psSphereRot *fromConj = psSphereRotConjugate(NULL, fromRot); + psSphereRot *out = psSphereRotCombine(NULL, toRot, fromConj); + psFree(from); + psFree(to); + psFree(fromEP); + psFree(fromCorr); + psFree(fromRot); + psFree(toEP); + psFree(toCorr); + psFree(toRot); + psFree(fromConj); + return out; + } else { //mode == PS_PRECESS_COMPLETE_B + //For PS_PRECESS_COMPLETE_B the same procedure as IAU2000A is used but with + //additional earthpole corrections from psEOC_PrecessionCorr. The corrections + //for COMPLETE_B come from the IERS Bulletin B. + + //Since the time inputs are allowed to be NULL, either convert the MJD time + //or copy to non-NULL time + if (fromTime == NULL) { + from = psTimeFromMJD(fromMJD); + } else { + from = p_psTimeCopy(fromTime); + } + if (toTime == NULL) { + to = psTimeFromMJD(toMJD); + } else { + to = p_psTimeCopy(toTime); + } + + psEarthPole *fromEP = psEOC_PrecessionModel(from); + psEarthPole *fromCorr = psEOC_PrecessionCorr(from, PS_IERS_B); + fromEP->x += fromCorr->x; + fromEP->y += fromCorr->y; + fromEP->s += fromCorr->s; + psSphereRot *fromRot = psSphereRot_CEOtoGCRS(fromEP); + psEarthPole *toEP = psEOC_PrecessionModel(to); + psEarthPole *toCorr = psEOC_PrecessionCorr(to, PS_IERS_B); + toEP->x += toCorr->x; + toEP->y += toCorr->y; + toEP->s += toCorr->s; + psSphereRot *toRot = psSphereRot_CEOtoGCRS(toEP); + psSphereRot *fromConj = psSphereRotConjugate(NULL, fromRot); + psSphereRot *out = psSphereRotCombine(NULL, toRot, fromConj); + psFree(from); + psFree(to); + psFree(fromEP); + psFree(fromCorr); + psFree(fromRot); + psFree(toEP); + psFree(toCorr); + psFree(toRot); + psFree(fromConj); + return out; + } // Apply transform to coordinates - psSphere *out = psSphereRotApply(NULL, tmpST, coords); - if (out->r < -0.0001) { - out->r += 2.0 * M_PI; - } - - psFree(tmpST); - - return(out); + // psSphere *out = psSphereRotApply(NULL, tmpST, coords); + // if (out->r < -0.0001) { + // out->r += 2.0 * M_PI; + // } } Index: trunk/psLib/src/astro/psEarthOrientation.h =================================================================== --- trunk/psLib/src/astro/psEarthOrientation.h (revision 6039) +++ trunk/psLib/src/astro/psEarthOrientation.h (revision 6184) @@ -9,6 +9,6 @@ * @author Robert Daniel DeSonia, MHPCC * -* @version $Revision: 1.13 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-01-18 23:49:06 $ +* @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-01-23 20:04:31 $ * * Copyright 2005 Maui High Performance Computing Center, University of Hawaii @@ -40,5 +40,6 @@ typedef enum { PS_PRECESS_ROUGH, ///< roughest, lowest level of detail - PS_PRECESS_COMPLETE, ///< complete level of detail + PS_PRECESS_COMPLETE_A, ///< complete level of detail using IERS A + PS_PRECESS_COMPLETE_B, ///< complete level of detail using IERS B PS_PRECESS_IAU2000A, ///< highest level of detail } @@ -179,10 +180,15 @@ * between two times. The equinoxes shall be Julian equinoxes. * - * @return psSphere* the resulting spherical rotation + * If NULL is provided for either time, it is assumed to have the reference + * equinox value of J2000. The mode argument is used to specify the level of + * detail used in the calculation. + * + * @return psSphere*: the resulting spherical rotation */ -psSphere* psSpherePrecess( - psSphere *coords, ///< coordinates (modified in-place) +psSphereRot* psSpherePrecess( + // psSphere *coords, ///< coordinates (modified in-place) const psTime *fromTime, ///< equinox of coords input - const psTime *toTime ///< equinox of coords output + const psTime *toTime, ///< equinox of coords output + psPrecessMethod mode ///< level of detail to use ); Index: trunk/psLib/src/astro/psTime.c =================================================================== --- trunk/psLib/src/astro/psTime.c (revision 6039) +++ trunk/psLib/src/astro/psTime.c (revision 6184) @@ -10,6 +10,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.77 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-01-18 20:59:31 $ + * @version $Revision: 1.78 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-01-23 20:04:31 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -920,4 +920,5 @@ PS_ASSERT_PTR_NON_NULL(time, NULL); psTime *out = NULL; + //Also see psEOC_PolarTideCorr for more info. // Convert psTime to MJD @@ -982,10 +983,9 @@ out = psTimeAlloc(time->type); *out = *time; - // out->sec = time->sec; - // out->nsec = time->nsec; - // out->leapsecond = time->leapsecond; if (out->type != PS_TIME_UT1) { out = psTimeConvert(out, PS_TIME_UT1); } + //see if corrections include seconds or just nano-seconds + //nano-seconds must be converted (scaled) to integer values & total cannot be negative if (fabs(CORZ) > 1.0) { int sec = (int)CORZ; @@ -1806,2 +1806,10 @@ } +psTime *p_psTimeCopy(const psTime *inTime) +{ + PS_ASSERT_PTR_NON_NULL(inTime, NULL); + psTime *outTime = psTimeAlloc(inTime->type); + *outTime = *inTime; + return outTime; +} + Index: trunk/psLib/src/astro/psTime.h =================================================================== --- trunk/psLib/src/astro/psTime.h (revision 6039) +++ trunk/psLib/src/astro/psTime.h (revision 6184) @@ -11,6 +11,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.44 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-01-18 20:59:31 $ + * @version $Revision: 1.45 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-01-23 20:04:31 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -432,4 +432,12 @@ psF64 psTimerStop(void); +/** Private function for copying a psTime. + * + * @return psTime*: New copy of existing psTime. + */ +psTime *p_psTimeCopy( + const psTime *inTime ///< input time to copy. +); + /// @}