Index: trunk/psLib/src/astro/psEarthOrientation.c =================================================================== --- trunk/psLib/src/astro/psEarthOrientation.c (revision 5493) +++ trunk/psLib/src/astro/psEarthOrientation.c (revision 5507) @@ -8,6 +8,6 @@ * @author Robert Daniel DeSonia, MHPCC * -* @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-11-10 00:13:50 $ +* @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-11-12 03:37:34 $ * * Copyright 2005 Maui High Performance Computing Center, University of Hawaii @@ -487,9 +487,137 @@ } - psEarthPole *psEOC_PrecessionCorr(const psTime *time, psTimeBulletin bulletin) { - return NULL; + PS_ASSERT_PTR_NON_NULL(time,NULL); + PS_ASSERT_INT_WITHIN_RANGE(bulletin, PS_IERS_A, PS_IERS_B, NULL); + + int tableColumn; + double mjd; + psLookupStatusType status; + char* tableName[1] = {"dailyTable"}; + double x, y, s; + // psMetadataItem *tableMetadataItem = NULL; + // psVector *temp = NULL; + // double calc = 0.0; + // double t = 0.0; + // psMetadata *timeMetadata = p_psTimeGetMetadata(); + + mjd = psTimeToMJD(time); + + if (bulletin == PS_IERS_A) { + tableColumn = 1; + } else { + tableColumn = 4; + } + + x = p_psTimeSearchTables(mjd, tableColumn, tableName, 1, &status); + if ( status == PS_LOOKUP_ERROR ) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Time table search returned unsuccessful status.\n"); + return NULL; + } else if (status != PS_LOOKUP_SUCCESS) { + psSphere *sphere = p_psTimeGetPoleCoords(time); + x = sphere->r; + psFree(sphere); + } else { + x = SEC_TO_RAD(x); + } + + + /* }else if(status == PS_LOOKUP_PAST_TOP) { + // Date too early for tables. Get default time delta value from metadata, and issue warning. + psLogMsg(__func__,PS_LOG_WARN,PS_ERRORTEXT_psTime_TIME_PREDATES_TABLES,mjd,"xp"); + + // Lookup value from time metadata loaded from pslib.config + tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.before.xp"); + if(tableMetadataItem == NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, + "psLib.time.before.xp"); + return NULL; + } + x = tableMetadataItem->data.F64; + + } else if(status == PS_LOOKUP_PAST_BOTTOM) { + tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.predict.mjd"); + if(tableMetadataItem == NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, + "psLib.time.predict.mjd"); + return NULL; + } + mjdPred = tableMetadataItem->data.F64; + + // Get xp + tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.predict.xp"); + if(tableMetadataItem == NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.predict.xp"); + return NULL; + } + xp = (psVector*)tableMetadataItem->data.V; + PS_ASSERT_PTR_NON_NULL(xp,NULL); + + // Get yp + tableMetadataItem = psMetadataLookup(timeMetadata, "psLib.time.predict.yp"); + if(tableMetadataItem == NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, "psLib.time.predict.yp"); + return NULL; + } + yp = (psVector*)tableMetadataItem->data.V; + PS_ASSERT_PTR_NON_NULL(yp,NULL); + + // Calculate "a" and "c" constants + a = TWOPI*(mjd - mjdPred)/365.25; + c = TWOPI*(mjd - mjdPred)/435.0; + + // Calculate x and y polar coordinates + x = xp->data.F64[0] + + xp->data.F64[1]*cos(a) + + xp->data.F64[2]*sin(a) + + xp->data.F64[3]*cos(c) + + xp->data.F64[4]*sin(c); + + + + + + } + */ + + + tableColumn++; + y = p_psTimeSearchTables(mjd, tableColumn, tableName, 1, &status); + if ( status == PS_LOOKUP_ERROR ) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Time table search returned unsuccessful status.\n"); + return NULL; + } else if (status != PS_LOOKUP_SUCCESS) { + psSphere *sphere = p_psTimeGetPoleCoords(time); + y = sphere->d; + psFree(sphere); + } else { + y = SEC_TO_RAD(y); + } + + tableColumn++; + s = p_psTimeSearchTables(mjd, tableColumn, tableName, 1, &status); + if ( status == PS_LOOKUP_ERROR ) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Time table search returned unsuccessful status.\n"); + return NULL; + } else if (status != PS_LOOKUP_SUCCESS) { + s = psTimeGetUT1Delta(time, bulletin); + } else { + s = SEC_TO_RAD(s); + } + + //XXX: Do I need to do something for nutation terms w/periods of less than 2 days??? + psEarthPole *out = psEarthPoleAlloc(); + out->x = x; + out->y = y; + out->s = s; + + // psFree(tableMetadataItem); + // psFree(timeMetadata); + return out; } @@ -566,5 +694,4 @@ } - psSphereRot* psSphereRot_TEOtoCEO(const psTime *time) { @@ -589,6 +716,34 @@ psTimeBulletin bulletin) { - - return NULL; + PS_ASSERT_PTR_NON_NULL(time, NULL); + + psEarthPole *out = psEarthPoleAlloc(); + psSphere *in = NULL; + double s; + + //XXX: This may be the wrong idea... ADD says to use 3rd-order polys to interpolate + //polar motion coordinates. + in = p_psTimeGetPoleCoords(time); + if (in == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + "p_psTimeGetPoleCoords return NULL psSphere for non-NULL input time.\n"); + return NULL; + } + out->x = in->r; + out->y = in->d; + + s = psTimeGetUT1Delta(time, bulletin); + out->s = s; + + //XXX: Apply polar tide correction here??? + psEarthPole *correction = psEOC_PolarTideCorr(time); + out->x += correction->x; + out->y += correction->y; + out->s += correction->s; + + + psFree(in); + psFree(correction); + return out; } @@ -711,5 +866,5 @@ diag_sum[2] = 1.0 - A[0][0] - A[1][1] + A[2][2]; diag_sum[3] = 1.0 + A[0][0] + A[1][1] + A[2][2]; - + maxi = 0; for (int i = 1; i < 4; ++i) { @@ -718,8 +873,8 @@ } } - + double p = 0.5 * sqrt(diag_sum[maxi]); recip = 1.0 / (4.0 * p); - + if (maxi == 0) { out->q0 = p;