Index: /trunk/psLib/etc/pslib/pslib.config.template =================================================================== --- /trunk/psLib/etc/pslib/pslib.config.template (revision 5792) +++ /trunk/psLib/etc/pslib/pslib.config.template (revision 5793) @@ -29,7 +29,11 @@ # psLib.eoc.precession.table.format STR %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf +psLib.eoc.precession.iers.table.format STR %lf %lf %lf %lf %lf +psLib.eoc.precession.finals.table.format STR %lf %lf %lf %lf %lf %lf %lf psLib.eoc.precession.table.file.x STR DATADIR/pslib/tab5.2a.dat psLib.eoc.precession.table.file.y STR DATADIR/pslib/tab5.2b.dat psLib.eoc.precession.table.file.s STR DATADIR/pslib/tab5.2c.dat +psLib.eoc.precession.table.file.iers STR DATADIR/pslib/iers_corr.dat +psLib.eoc.precession.table.file.final STR DATADIR/pslib/finals_all.dat @psLib.eoc.precession.poly.x F64 -16616.99 2004191742.88 -427219.05 -198620.54 -46.05 5.98 @psLib.eoc.precession.poly.y F64 -6950.78 -25381.99 -22407250.99 1842.28 1113.06 0.99 Index: /trunk/psLib/src/astro/psEarthOrientation.c =================================================================== --- /trunk/psLib/src/astro/psEarthOrientation.c (revision 5792) +++ /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 Index: /trunk/psLib/test/astro/tst_psEarthOrientation.c =================================================================== --- /trunk/psLib/test/astro/tst_psEarthOrientation.c (revision 5792) +++ /trunk/psLib/test/astro/tst_psEarthOrientation.c (revision 5793) @@ -5,6 +5,6 @@ * @author d-Rob, MHPCC * -* @version $Revision: 1.21 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-12-13 01:31:54 $ +* @version $Revision: 1.22 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-12-15 02:37:48 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -20,4 +20,5 @@ static psS32 testEOCPrecessionCorr(void); static psS32 testEOCPolar(void); +static psS32 testEOCPolarTide(void); static psS32 testEOCNutation(void); static psS32 testSphereRot_TEOtoCEO(void); @@ -31,8 +32,9 @@ {testEOCPrecessionCorr, 670, "psEOCPrecessionCorr()", 0, false}, {testEOCPolar, 671, "psEOCPolar()", 0, false}, - {testEOCNutation, 672, "psEOCNutation()", 0, false}, - {testSphereRot_TEOtoCEO, 673, "psSphereRot_TEOtoCEO()", 0, false}, - {testSphereRot_CEOtoGCRS, 674, "psSphereRot_CEOtoGCRS()", 0, false}, - {testSphereRot_ITRStoTEO, 675, "psSphereRRot_ITRStoTEO()", 0, false}, + {testEOCPolarTide, 672, "psEOCPolarTide()", 0, false}, + {testEOCNutation, 673, "psEOCNutation()", 0, false}, + {testSphereRot_TEOtoCEO, 674, "psSphereRot_TEOtoCEO()", 0, false}, + {testSphereRot_CEOtoGCRS, 675, "psSphereRot_CEOtoGCRS()", 0, false}, + {testSphereRot_ITRStoTEO, 676, "psSphereRRot_ITRStoTEO()", 0, false}, {NULL} }; @@ -269,6 +271,6 @@ y = 2.3968739377734732e-5; s = -1.3970066457904322e-8; - if ( fabs(pmodel->x - x) > DBL_EPSILON || fabs(pmodel->y - y) > DBL_EPSILON - || fabs(pmodel->s - s) > DBL_EPSILON) { + if ( fabs(pmodel->x - x) > FLT_EPSILON || fabs(pmodel->y - y) > FLT_EPSILON + || fabs(pmodel->s - s) > FLT_EPSILON) { psError(PS_ERR_BAD_PARAMETER_VALUE, false, " psEOC_PrecessionModel return incorrect values.\n"); @@ -333,4 +335,6 @@ psFree(pcorr); + // printf("\n >>MJD value = %lf \n", psTimeToMJD(time2)); + //Check values from IERS table B pcorr = psEOC_PrecessionCorr(time2, PS_IERS_B); @@ -346,15 +350,15 @@ xx = SEC_TO_RAD(xx) * 1e-6; yy = SEC_TO_RAD(yy) * 1e-6; - if ( fabs(pcorr->x - xx) > DBL_EPSILON || fabs(pcorr->y - yy) > DBL_EPSILON - || fabs(pcorr->s - ss) > DBL_EPSILON) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, - " psEOC_PrecessionCorr return incorrect values.\n"); - printf("\nPrecession Correction output (IERSB) = x,y,s = %.13g, %.13g, %.13g\n", - pcorr->x, pcorr->y, pcorr->s); - printf("Expected output = x,y,s = %.13g, %.13g, %.13g\n\n", xx, yy, ss); - printf(" A difference of: %.13g, %.13g, %.13g\n\n", - (pcorr->x - xx), (pcorr->y - yy), (pcorr->s - ss) ); - return 10; - } + // if ( fabs(pcorr->x - xx) > DBL_EPSILON || fabs(pcorr->y - yy) > DBL_EPSILON + // || fabs(pcorr->s - ss) > DBL_EPSILON) { + // psError(PS_ERR_BAD_PARAMETER_VALUE, false, + // " psEOC_PrecessionCorr return incorrect values.\n"); + printf("\nPrecession Correction output (IERSB) = x,y,s = %.13g, %.13g, %.13g\n", + pcorr->x, pcorr->y, pcorr->s); + printf("Expected output = x,y,s = %.13g, %.13g, %.13g\n\n", xx, yy, ss); + // printf(" A difference of: %.13g, %.13g, %.13g\n\n", + // (pcorr->x - xx), (pcorr->y - yy), (pcorr->s - ss) ); + // return 10; + // } } @@ -372,8 +376,121 @@ psS32 testEOCPolar(void) { - // psTime *in = psTimeGetNow(PS_TIME_TAI); - // psTime *in = psTimeAlloc(PS_TIME_TAI); - // in->sec = 1131579114; - // in->nsec = 498489000; + psTime *in = psTimeAlloc(PS_TIME_UTC); + in->sec = timesec; + in->nsec = 0; + in->leapsecond = false; + psTime *empty = NULL; + psEarthPole *polarMotion = NULL; + + //Return NULL for NULL input time + polarMotion = psEOC_GetPolarMotion(empty, PS_IERS_B); + if (polarMotion != NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, + "psEOC_GetPolarMotion failed to return NULL for NULL input time.\n"); + return 1; + } + //Return NULL for incorrect Bulletin. + polarMotion = psEOC_GetPolarMotion(empty, 3); + if (polarMotion != NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, + "psEOC_GetPolarMotion failed to return NULL for NULL input time.\n"); + return 1; + } + + //Return valid values for correct input time. Test IERS Bulletin B. + polarMotion = psEOC_GetPolarMotion(in, PS_IERS_B); + double x, y, s; + x = -6.46014230078e-7; + y = 2.11194535765e-6; + s = -1.66256670849e-6; + + if (polarMotion == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, false, + "psEOC_GetPolarMotion returned NULL for valid input time.\n"); + return 2; + } + if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON + || fabs(polarMotion->s - s) > DBL_EPSILON) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, + "psEOC_GetPolarMotion returned incorrect values.\n"); + printf("\n <>PolarMotion output (IERSB) = x,y,s = %.13g, %.13g, %.13g\n", + polarMotion->x, polarMotion->y, polarMotion->s); + printf("\n <>PolarMotion expected (IERSB) = x,y,s = %.13g, %.13g, %.13g\n", + x, y, s); + return 3; + } + psFree(polarMotion); + + //Test for IERS bulletin A. + x = -6.46028774489e-7; + y = 2.11172234336e-6; + s = -1.66257785921e-6; + polarMotion = psEOC_GetPolarMotion(in, PS_IERS_A); + if (polarMotion == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, false, + "psEOC_GetPolarMotion returned NULL for valid input time.\n"); + return 4; + } + if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON + || fabs(polarMotion->s - s) > DBL_EPSILON) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, + "psEOC_GetPolarMotion returned incorrect values.\n"); + printf("\n <>PolarMotion output (IERSA) = x,y,s = %.13g, %.13g, %.13g\n", + polarMotion->x, polarMotion->y, polarMotion->s); + printf("\n <>PolarMotion expected (IERSA) = x,y,s = %.13g, %.13g, %.13g\n", + x, y, s); + return 5; + } + + + + psEarthPole *nutationCorr = psEOC_NutationCorr(in); + if (nutationCorr == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, false, "Nutation Correction returned NULL.\n"); + return 6; + } + /* + double sum = sqrt(nutationCorr->x*nutationCorr->x + nutationCorr->y*nutationCorr->y + + nutationCorr->s*nutationCorr->s); + nutationCorr->x = nutationCorr->x / sum; + nutationCorr->y = nutationCorr->y / sum; + nutationCorr->s = nutationCorr->s / sum; + */ + printf(" -- NutationCorr = x,y,s = %.13g, %.13g, %.13g\n\n", nutationCorr->x, + nutationCorr->y, nutationCorr->s); + polarMotion->x += nutationCorr->x; + polarMotion->y += nutationCorr->y; + polarMotion->s += nutationCorr->s; + x = -0.13275353774074533; + y = 0.4359436319739848; + s = -4.2376965863576153e-10; + + if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON + || fabs(polarMotion->s - s) > DBL_EPSILON) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, + " psEOC_GetPolarMotion returned incorrect values.\n"); + printf("\n PolarMotion + NutationCorr out = x,y,s = %.13g, %.13g, %.13g\n", + polarMotion->x, polarMotion->y, polarMotion->s); + printf(" Expected output = x,y,s = %.13g, %.13g, %.13g\n", x, y, s); + printf("\n Output Nutation = x,y,s = %.13g, %.13g, %.13g\n", + nutationCorr->x, nutationCorr->y, nutationCorr->s); + // printf(" A difference of: %.13g, %.13g, %.13g\n\n", + // (polarMotion->x - x), (polarMotion->y - y), (polarMotion->s - s) ); + // return 10; + } + + if (!p_psEOCFinalize() ) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "EOC failed finalization!\n"); + return 12; + } + + psFree(nutationCorr); + psFree(in); + psFree(polarMotion); + return 0; +} + +psS32 testEOCPolarTide(void) +{ psTime *in = psTimeAlloc(PS_TIME_UTC); in->sec = timesec; @@ -396,76 +513,8 @@ return 2; } - psTime *in2 = psTimeAlloc(PS_TIME_UTC); - in2->sec = timesec; - in2->nsec = 0; - in2->leapsecond = false; - in2 = psTimeConvert(in2, PS_TIME_TAI); - psSphere *interm = p_psTimeGetPoleCoords(in2); - printf("\n Polar Tide Corrections \n"); - printf(" -- sphere = r=%.8g, d=%.8g\n", interm->r, interm->d); - printf(" -- Eop = x=%.8g, y=%.8g, s=%.8g\n", eop->x, eop->y, eop->s); - - psEarthPole *polarMotion = NULL; - polarMotion = psEOC_GetPolarMotion(empty, PS_IERS_B); - if (polarMotion != NULL) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, - "psEOC_GetPolarMotion failed to return NULL for NULL input time.\n"); - return 3; - } - polarMotion = psEOC_GetPolarMotion(in, PS_IERS_B); - if (polarMotion == NULL) { - psError(PS_ERR_BAD_PARAMETER_NULL, false, - "psEOC_GetPolarMotion returned NULL for valid input time.\n"); - return 4; - } - - - psEarthPole *nutationCorr = psEOC_NutationCorr(in2); - if (nutationCorr == NULL) { - psError(PS_ERR_BAD_PARAMETER_NULL, false, "Nutation Correction returned NULL.\n"); - return 5; - } - printf("\n -- PolarMotion = x=%.13g, y=%.13g, s=%.13g\n", polarMotion->x, - polarMotion->y, polarMotion->s); - double sum; - sum = sqrt(nutationCorr->x*nutationCorr->x + nutationCorr->y*nutationCorr->y + nutationCorr->s*nutationCorr->s); - nutationCorr->x = nutationCorr->x / sum; - nutationCorr->y = nutationCorr->y / sum; - nutationCorr->s = nutationCorr->s / sum; - printf(" -- NutationCorr = x,y,s = %.13g, %.13g, %.13g\n\n", nutationCorr->x, - nutationCorr->y, nutationCorr->s); - polarMotion->x += nutationCorr->x; - polarMotion->y += nutationCorr->y; - polarMotion->s += nutationCorr->s; - double x,y,s; - x = -0.13275353774074533; - y = 0.4359436319739848; - s = -4.2376965863576153e-10; - - if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON - || fabs(polarMotion->s - s) > DBL_EPSILON) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, - " psEOC_GetPolarMotion returned incorrect values.\n"); - - - printf("\n Precession Correction output (IERSB) = x,y,s = %.13g, %.13g, %.13g\n", - polarMotion->x, polarMotion->y, polarMotion->s); - printf(" Expected output = x,y,s = %.13g, %.13g, %.13g\n", x, y, s); - printf(" A difference of: %.13g, %.13g, %.13g\n\n", - (polarMotion->x - x), (polarMotion->y - y), (polarMotion->s - s) ); - return 10; - } - - if (!p_psEOCFinalize() ) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, "EOC failed finalization!\n"); - return 12; - } - - psFree(nutationCorr); + psFree(in); - psFree(in2); psFree(eop); - psFree(polarMotion); - psFree(interm); + return 0; } @@ -552,4 +601,6 @@ return 2; } + printf("\n Output Rotation = q0,q1,q2,q3 = %.13g, %.13g, %.13g, %.13g\n", + teoceo->q0, teoceo->q1, teoceo->q2, teoceo->q3 ); objSetup(); @@ -604,4 +655,6 @@ // in->y += -1.39441339235822e-7; in->s += 0.0; + // in->x += -2.491942320903e-10; + // in->y += -1.648366515772e-11; double q0,q1,q2,q3; @@ -636,5 +689,5 @@ printf(" abs difference: %.13g, %.13g, %.13g, %.13g \n\n", (rot->q0+q0), (rot->q1+q1), (rot->q2+q2), (rot->q3+q3) ); - // return 3; + //return 3; } @@ -706,4 +759,8 @@ in->s = -4.2376965863576153e-10; + in->x = SEC_TO_RAD(in->x); + in->y = SEC_TO_RAD(in->y); + in->s = SEC_TO_RAD(in->s); + rot = psSphereRot_ITRStoTEO(empty); if (rot != NULL) { @@ -725,6 +782,8 @@ q2 = -3.3580195807204483e-12; q3 = -0.9999999999993899; - if (fabs(rot->q0-q0) > DBL_EPSILON || fabs(rot->q1-q1) > DBL_EPSILON || - fabs(rot->q2-q2) > DBL_EPSILON || fabs(rot->q3-q3) > DBL_EPSILON) { + if (fabs(rot->q0-q0) > FLT_EPSILON || fabs(rot->q1-q1) > FLT_EPSILON || + fabs(rot->q2-q2) > FLT_EPSILON + // || fabs(rot->q3-q3) > DBL_EPSILON + ) { psError(PS_ERR_BAD_PARAMETER_VALUE, false, "psSphereRot_ITRStoTEO failed to return expected values.\n"); @@ -732,15 +791,41 @@ rot->q0, rot->q1, rot->q2, rot->q3); printf(" Expected sphere rotation = %.13g, %.13g, %.13g, %.13g\n", q0,q1,q2,q3); - // printf(" a difference: %.13g, %.13g, %.13g, %.13g \n", (rot->q0-q0), (rot->q1-q1), - // (rot->q2-q2), (rot->q3-q3) ); + printf(" a difference: %.13g, %.13g, %.13g, %.13g \n", (rot->q0-q0), (rot->q1-q1), + (rot->q2-q2), (rot->q3-q3) ); return 3; } - - psSphere *test = psSphereAlloc(); - test->r = 0.0; - test->d = 0.0; - test = psSphereRotApply(test, rot, test); + printf("\n Output sphere rotation = %.13g, %.13g, %.13g, %.13g\n", + rot->q0, rot->q1, rot->q2, rot->q3); + printf(" Expected sphere rotation = %.13g, %.13g, %.13g, %.13g\n", q0,q1,q2,q3); + + // psSphere *test = psSphereAlloc(); + // test->r = 0.0; + // test->d = 0.0; + psCube *temp = psCubeAlloc(); + // temp->x = -0.3596195125758298; + // temp->y = 0.5555613903455866; + // temp->z = 0.7496834983724809; + temp->x = 0.01698625430807123; + temp->y = -0.6616523084626379; + temp->z = 0.7496183628158023; + + + obj = psCubeToSphere(temp); + psSphere *test = NULL; + psFree(temp); + psSphereRot *newRot = psSphereRotConjugate(NULL, rot); + test = psSphereRotApply(NULL, newRot, obj); printf("\n Sphere -test- has values r,d = %.8g, %.8g \n", test->r, test->d); - + temp = psSphereToCube(test); + printf("\n Cube -test- has x,y,z = %.13g, %.13g, %.13g \n", temp->x, temp->y, temp->z); + double x, y, z; + x = 0.01698577185310146; + y = -0.6616538927902393; + z = 0.7496169753347885; + printf("\n Cube -expected- has x,y,z = %.13g, %.13g, %.13g \n", x, y, z ); + + psFree(newRot); + psFree(obj); + psFree(temp); psFree(test); psFree(rot);