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; + // } }