Index: trunk/psLib/test/astro/tst_psEarthOrientation.c =================================================================== --- trunk/psLib/test/astro/tst_psEarthOrientation.c (revision 5771) +++ 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);