Index: trunk/psLib/src/astro/psEarthOrientation.c =================================================================== --- trunk/psLib/src/astro/psEarthOrientation.c (revision 5771) +++ trunk/psLib/src/astro/psEarthOrientation.c (revision 5793) @@ -8,6 +8,6 @@ * @author Robert Daniel DeSonia, MHPCC * - * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-12-13 01:31:54 $ + * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-12-15 02:37:48 $ * * Copyright 2005 Maui High Performance Computing Center, University of Hawaii @@ -46,4 +46,6 @@ static psArray* yTable = NULL; static psArray* sTable = NULL; +static psArray* iersTable = NULL; +static psArray* finalsTable = NULL; static psPolynomial1D* xPoly = NULL; static psPolynomial1D* yPoly = NULL; @@ -86,5 +88,6 @@ bool success = false; // Get table format - char* tableFormat = psMetadataLookupStr(&success, eocMetadata, "psLib.eoc.precession.table.format"); + char* tableFormat = psMetadataLookupStr(&success, eocMetadata, + "psLib.eoc.precession.table.format"); if(! success || tableFormat == NULL) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, @@ -93,5 +96,6 @@ } - char* xTableName = psMetadataLookupStr(&success, eocMetadata, "psLib.eoc.precession.table.file.x"); + char* xTableName = psMetadataLookupStr(&success, eocMetadata, + "psLib.eoc.precession.table.file.x"); if(! success || xTableName == NULL) { psFree(tableFormat); @@ -110,5 +114,6 @@ } - char* sTableName = psMetadataLookupStr(&success, eocMetadata, "psLib.eoc.precession.table.file.s"); + char* sTableName = psMetadataLookupStr(&success, eocMetadata, + "psLib.eoc.precession.table.file.s"); if(! success || sTableName == NULL) { psFree(tableFormat); @@ -119,11 +124,72 @@ return false; } + + char* iersTableName = psMetadataLookupStr(&success, eocMetadata, + "psLib.eoc.precession.table.file.iers"); + if(! success || sTableName == NULL) { + psFree(tableFormat); + psFree(xTableName); + psFree(yTableName); + psFree(sTableName); + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, + "psLib.eoc.precession.iers.file"); + return false; + } + + char* finalsTableName = psMetadataLookupStr(&success, eocMetadata, + "psLib.eoc.precession.table.file.final"); + if(! success || sTableName == NULL) { + psFree(tableFormat); + psFree(xTableName); + psFree(yTableName); + psFree(sTableName); + psFree(iersTableName); + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, + "psLib.eoc.precession.finals.file"); + return false; + } + + char* iersTableFormat = psMetadataLookupStr(&success, eocMetadata, + "psLib.eoc.precession.iers.table.format"); + if(! success || iersTableFormat == NULL) { + psFree(tableFormat); + psFree(xTableName); + psFree(yTableName); + psFree(sTableName); + psFree(iersTableName); + psFree(finalsTableName); + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, + "psLib.eoc.precession.iers.table.format"); + return false; + } + + char* finalsTableFormat = psMetadataLookupStr(&success, eocMetadata, + "psLib.eoc.precession.finals.table.format"); + if(! success || iersTableFormat == NULL) { + psFree(tableFormat); + psFree(xTableName); + psFree(yTableName); + psFree(sTableName); + psFree(iersTableName); + psFree(finalsTableName); + psFree(iersTableFormat); + psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psTime_LOOKUP_METADATA_FAILED, + "psLib.eoc.precession.iers.table.format"); + return false; + } + xTable = psVectorsReadFromFile(xTableName, tableFormat); yTable = psVectorsReadFromFile(yTableName, tableFormat); sTable = psVectorsReadFromFile(sTableName, tableFormat); + iersTable = psVectorsReadFromFile(iersTableName, iersTableFormat); + finalsTable = psVectorsReadFromFile(finalsTableName, finalsTableFormat); psFree(tableFormat); + psFree(iersTableFormat); + psFree(finalsTableFormat); psFree(xTableName); psFree(yTableName); psFree(sTableName); + psFree(finalsTableName); + psFree(iersTableName); if(xTable == NULL || yTable == NULL || sTable == NULL) { @@ -131,4 +197,11 @@ psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Failed to read the precession-nutation model tables."); + return false; + } + + if(iersTable == NULL || finalsTable == NULL) { + // XXX: need to move the error message to the error messages file. + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Failed to read the IERS/finals tables."); return false; } @@ -175,9 +248,12 @@ psFree(yTable); psFree(sTable); // There are non dynamic allocated items - + psFree(iersTable); + psFree(finalsTable); xTable = NULL; yTable = NULL; sTable = NULL; + iersTable = NULL; + finalsTable = NULL; psFree(xPoly); @@ -518,130 +594,50 @@ 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 ) { + psEarthPole *out = psEarthPoleAlloc(); + out->x = 0.0; + out->y = 0.0; + out->s = 0.0; + + double MJD = psTimeToMJD(time); + + if (MJD < 41684.0 || MJD > 53334.0) { 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; + "Invalid time input. Date, %lf, is out of range (41684-53334)\n", MJD); + return out; + } + + // Check if EOC data loaded + if(!eocInitialized) { + eocInitialized = p_psEOCInit(); + if(!eocInitialized) { + // XXX: Move error message. + psError(PS_ERR_UNKNOWN, false, + "Could not initialize EOC tables -- check data files."); + return NULL; + } + } + + psF64* cols[5]; + for (int colNum = 0; colNum < 5; colNum++) { + cols[colNum] = ((psVector*)(iersTable->data[colNum]))->data.F64; + } + int numRows = ((psVector*)(iersTable->data[0]))->n; + for (int rowNum = 0; rowNum < numRows; rowNum++) { + if ( (MJD - cols[0][rowNum]) < 1.0 ) { + if (bulletin == PS_IERS_A) { + out->x = cols[1][rowNum]; + out->y = cols[2][rowNum]; + out->x = SEC_TO_RAD(out->x) * 1e-6; + out->y = SEC_TO_RAD(out->y) * 1e-6; + rowNum = numRows; + } else { + out->x = cols[3][rowNum]; + out->y = cols[4][rowNum]; + out->x = SEC_TO_RAD(out->x) * 1e-6; + out->y = SEC_TO_RAD(out->y) * 1e-6; + rowNum = numRows; } - 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; } @@ -779,41 +775,92 @@ psEarthPole *out = psEarthPoleAlloc(); - psSphere *in = NULL; - double s; - psTime *time2 = psTimeAlloc(time->type); - time2->sec = time->sec; - time2->nsec = time->nsec; - time2->leapsecond = time->leapsecond; + // psSphere *in = NULL; + // double x,y,s; + // psTime *time2 = psTimeAlloc(time->type); + // time2->sec = time->sec; + // time2->nsec = time->nsec; + // time2->leapsecond = time->leapsecond; //XXX: Time must be converted to TAI before use? - if (time->type == PS_TIME_UT1) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "psTime cannot be of type UT1 here.\n"); - return NULL; - } else if (time->type != PS_TIME_TAI) { - time2 = psTimeConvert(time2, PS_TIME_TAI); - } - + /* if (time->type == PS_TIME_UT1) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "psTime cannot be of type UT1 here.\n"); + return NULL; + } else if (time->type != PS_TIME_TAI) { + time2 = psTimeConvert(time2, PS_TIME_TAI); + } + */ //XXX: This may be the wrong idea... ADD says to use 3rd-order polys to interpolate //polar motion coordinates. - in = p_psTimeGetPoleCoords(time2); - 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(time2, bulletin); - out->s = s; + // in = p_psTimeGetPoleCoords(time2); + // 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(time2, bulletin); + // out->s = s; //XXX: Apply polar tide correction here??? - psEarthPole *correction = psEOC_PolarTideCorr(time2); - out->x += correction->x; - out->y += correction->y; - out->s += correction->s; - - psFree(time2); - psFree(in); - psFree(correction); + // psEarthPole *correction = psEOC_PolarTideCorr(time2); + // out->x += correction->x; + // out->y += correction->y; + // out->s += correction->s; + + // psFree(time2); + // psFree(in); + // psFree(correction); + + out->x = 0.0; + out->y = 0.0; + out->s = 0.0; + + double MJD = psTimeToMJD(time); + + if (MJD < 41684.0 || (MJD > 53434.0 && bulletin == PS_IERS_B) || + (MJD > 53858.0 && bulletin == PS_IERS_A) ) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Invalid time input. Date, %lf, is out of range\n", MJD); + return out; + } + + // Check if EOC data loaded + if(! eocInitialized) { + eocInitialized = p_psEOCInit(); + if(!eocInitialized) { + // XXX: Move error message. + psError(PS_ERR_UNKNOWN, false, + "Could not initialize EOC tables -- check data files."); + return NULL; + } + } + + psF64* cols[7]; + for (int colNum = 0; colNum < 7; colNum++) { + cols[colNum] = ((psVector*)(finalsTable->data[colNum]))->data.F64; + } + int numRows = ((psVector*)(finalsTable->data[0]))->n; + for (int rowNum = 0; rowNum < numRows; rowNum++) { + if ( (MJD - cols[0][rowNum]) < 1.0 ) { + if (bulletin == PS_IERS_A) { + out->x = cols[1][rowNum]; + out->y = cols[2][rowNum]; + out->s = cols[3][rowNum]; + out->x = SEC_TO_RAD(out->x); + out->y = SEC_TO_RAD(out->y); + out->s = SEC_TO_RAD(out->s); + rowNum = numRows; + } else { + out->x = cols[4][rowNum]; + out->y = cols[5][rowNum]; + out->s = cols[6][rowNum]; + out->x = SEC_TO_RAD(out->x); + out->y = SEC_TO_RAD(out->y); + out->s = SEC_TO_RAD(out->s); + rowNum = numRows; + } + } + } + return out; } @@ -836,5 +883,6 @@ // Calculate number of Julian centuries since 2000 //XXX: NOT SURE IF THIS IS CORRECT FOR THIS SITUATION - double RJD = ( MJD - MJD_2000 ) / JULIAN_CENTURY; + // double RJD = ( MJD - MJD_2000 ) / JULIAN_CENTURY; + double RJD = MJD; //Formula comes from fortran reference @@ -1005,9 +1053,12 @@ X = psPolynomial1DEval(xPoly,t); Y = psPolynomial1DEval(yPoly,t); + X = SEC_TO_RAD(X * 1e-6); + Y = SEC_TO_RAD(Y * 1e-6); + for (int i = 0; i < 10; i++) { tj = pow(t, i); arg = w_l[i]*F[0] + w_l_p[i]*F[1] + w_F[i]*F[2] + w_D[i]*F[3] + w_Omega[i]*F[4]; - X += xp_sin[i] * tj * sin(arg) + xp_cos[i] * tj * cos(arg); - Y += yp_sin[i] * tj * sin(arg) + yp_cos[i] * tj * cos(arg); + X += (xp_sin[i] * 1e-6 * tj * sin(arg) + xp_cos[i] * 1e-6 * cos(arg)) * tj; + Y += (yp_sin[i] * 1e-6 * tj * sin(arg) + yp_cos[i] * 1e-6 * cos(arg)) * tj; } @@ -1079,42 +1130,41 @@ /* - psSphereRot r,s,t; + psSphereRot r,s,t; - // directly from ADD -- there must be a better way?! - r.q0=sin(motion->y/2.0); - r.q1=0; - r.q2=0; - r.q3=cos(motion->y/2.0); + r.q0=sin(motion->y/2.0); + r.q1=0; + r.q2=0; + r.q3=cos(motion->y/2.0); - s.q0=0; - s.q1=sin(motion->x/2.0); - s.q2=0; - s.q3=cos(motion->x/2.0); + s.q0=0; + s.q1=sin(motion->x/2.0); + s.q2=0; + s.q3=cos(motion->x/2.0); - t.q0=0; - t.q1=0; - t.q2=sin(-motion->s/2.0); - t.q3=cos(-motion->s/2.0); + t.q0=0; + t.q1=0; + t.q2=sin(motion->s/2.0); + t.q3=cos(motion->s/2.0); - // calculate t*s*r. - psSphereRot* temp = psSphereRotCombine(NULL,&t,&s); - out = psSphereRotCombine(NULL, temp, &r); - psFree(temp); + // calculate t*s*r. + psSphereRot* temp = psSphereRotCombine(NULL,&t,&s); + out = psSphereRotCombine(NULL, temp, &r); + psFree(temp); - return out; + return out; */ s = -s; //Newest trial - mult. y matrix * x * z - /* A[0][0] = cos(x)*cos(s); - A[1][0] = cos(x)*sin(s); - A[2][0] = -sin(x); - A[0][1] = sin(x)*sin(y)*cos(s) - cos(y)*sin(s); - A[1][1] = sin(x)*sin(y)*sin(s) + cos(y)*cos(s); - A[2][1] = cos(x)*sin(y); - A[0][2] = sin(x)*cos(y)*cos(s) + sin(y)*sin(s); - A[1][2] = sin(x)*cos(y)*sin(s) - sin(y)*cos(s); - A[2][2] = cos(x)*cos(y); - */ + A[0][0] = cos(x)*cos(s); + A[1][0] = cos(x)*sin(s); + A[2][0] = -sin(x); + A[0][1] = sin(x)*sin(y)*cos(s) - cos(y)*sin(s); + A[1][1] = sin(x)*sin(y)*sin(s) + cos(y)*cos(s); + A[2][1] = cos(x)*sin(y); + A[0][2] = sin(x)*cos(y)*cos(s) + sin(y)*sin(s); + A[1][2] = sin(x)*cos(y)*sin(s) - sin(y)*cos(s); + A[2][2] = cos(x)*cos(y); + /* A[0][0] = cos(x)*cos(s); A[0][1] = cos(x)*sin(s); @@ -1138,14 +1188,15 @@ A[2][2] = cos(x)*cos(y); */ - A[0][0] = cos(s)*cos(x); - A[0][1] = sin(s)*cos(y) + cos(s)*sin(x)*sin(y); - A[0][2] = sin(s)*sin(y) - cos(s)*sin(x)*cos(y); - A[1][0] = -sin(s)*cos(x); - A[1][1] = cos(s)*cos(y) - sin(s)*sin(x)*sin(y); - A[1][2] = cos(s)*sin(y) + sin(s)*sin(x)*cos(y); - A[2][0] = sin(x); - A[2][1] = -cos(x)*sin(y); - A[2][2] = cos(x)*cos(y); - + /* + A[0][0] = cos(s)*cos(x); + A[0][1] = sin(s)*cos(y) + cos(s)*sin(x)*sin(y); + A[0][2] = sin(s)*sin(y) - cos(s)*sin(x)*cos(y); + A[1][0] = -sin(s)*cos(x); + A[1][1] = cos(s)*cos(y) - sin(s)*sin(x)*sin(y); + A[1][2] = cos(s)*sin(y) + sin(s)*sin(x)*cos(y); + A[2][0] = sin(x); + A[2][1] = -cos(x)*sin(y); + A[2][2] = cos(x)*cos(y); + */ //Convert rotation matrix to quaternions