Index: /trunk/psLib/src/astro/psEarthOrientation.c =================================================================== --- /trunk/psLib/src/astro/psEarthOrientation.c (revision 5641) +++ /trunk/psLib/src/astro/psEarthOrientation.c (revision 5642) @@ -8,6 +8,6 @@ * @author Robert Daniel DeSonia, MHPCC * - * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-30 02:17:17 $ + * @version $Revision: 1.20 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-11-30 23:50:40 $ * * Copyright 2005 Maui High Performance Computing Center, University of Hawaii @@ -727,8 +727,9 @@ } double T = (double)(in->sec) + (double)(in->nsec / 1e9); - printf("\nThe Value of T is = %.13g\n", T); + printf("\nThe Value of T is = %.19g\n", T); T += -2451545.0; - printf("\nThe Value of T is = %.13g\n", T); double theta = 2.0 * M_PI * (0.7790572732640 + 1.00273781191135448 * T); + printf("\nThe Value of theta is = %.13g\n", theta); + // psSphereRot *out = psSphereRotInvert(theta, 0.0, 0.0); psSphereRot *out = psSphereRotInvert(theta, 0.0, 0.0); @@ -745,8 +746,19 @@ psSphere *in = NULL; double 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); + } //XXX: This may be the wrong idea... ADD says to use 3rd-order polys to interpolate //polar motion coordinates. - in = p_psTimeGetPoleCoords(time); + in = p_psTimeGetPoleCoords(time2); if (in == NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, true, @@ -757,14 +769,14 @@ out->y = in->d; - s = psTimeGetUT1Delta(time, bulletin); + s = psTimeGetUT1Delta(time2, bulletin); out->s = s; //XXX: Apply polar tide correction here??? - psEarthPole *correction = psEOC_PolarTideCorr(time); + psEarthPole *correction = psEOC_PolarTideCorr(time2); out->x += correction->x; out->y += correction->y; out->s += correction->s; - + psFree(time2); psFree(in); psFree(correction); Index: /trunk/psLib/test/astro/tst_psEarthOrientation.c =================================================================== --- /trunk/psLib/test/astro/tst_psEarthOrientation.c (revision 5641) +++ /trunk/psLib/test/astro/tst_psEarthOrientation.c (revision 5642) @@ -5,6 +5,6 @@ * @author d-Rob, MHPCC * -* @version $Revision: 1.16 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-11-30 02:17:17 $ +* @version $Revision: 1.17 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-11-30 23:50:40 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -13,4 +13,5 @@ #include "psTest.h" #include "pslib_strict.h" + static psS32 testAberration(void); @@ -54,22 +55,40 @@ } +#define timesec 1049160600 + psS32 testAberration(void) { + psSphere *apparent = NULL; psSphere *actual = psSphereAlloc(); - psSphere *direction = psSphereAlloc(); + // psSphere *direction = psSphereAlloc(); psSphere *empty = NULL; // actual->r = 0.2; // actual->d = 0.2; - actual->r = DEG_TO_RAD(45.0); - actual->d = DEG_TO_RAD(30.0); + // actual->r = DEG_TO_RAD(45.0); + // actual->d = DEG_TO_RAD(30.0); + actual->r = DEG_TO_RAD(122.9153182445501); + actual->d = DEG_TO_RAD(48.562968978679194); // direction->r = 0.2035; // direction->d = 0.2035; // direction->r = DEG_TO_RAD(48.0); - direction->r = DEG_TO_RAD(-157.0); - direction->d = DEG_TO_RAD(20.7072); - - empty = psAberration(empty, apparent, direction, 0.1); + // direction->r = DEG_TO_RAD(-157.0); + // direction->d = DEG_TO_RAD(20.7072); + psCube *cubeDir = psCubeAlloc(); + cubeDir->x = 5148.713262821658; + cubeDir->y = -26945.04752348012; + cubeDir->z = -11682.787302030947; + cubeDir->x += -357.6031690489248; + cubeDir->y += 248.46429758174693; + cubeDir->z += 0.09694774143797581; + psSphere *direction = psCubeToSphere(cubeDir); + double speed = sqrt(cubeDir->x*cubeDir->x + cubeDir->y*cubeDir->y + cubeDir->z*cubeDir->z); + // Speed of light in vacuum (src:NIST) + double c; + c = 299792458.0; /* m/s */ + speed = speed / c; + + empty = psAberration(empty, apparent, direction, speed); if (empty != NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, false, @@ -77,5 +96,5 @@ return 1; } - empty = psAberration(empty, actual, apparent, 0.1); + empty = psAberration(empty, actual, apparent, speed); if (empty != NULL) { psError(PS_ERR_BAD_PARAMETER_NULL, false, @@ -84,10 +103,14 @@ } - apparent = psAberration(apparent, actual, direction, 0.4); + apparent = psAberration(apparent, actual, direction, speed); if (apparent == NULL) printf("\nApparent is NULL\n"); - else + else { printf("\napparent = r,d = %.8g, %.8g\n", apparent->r, apparent->d); - + psFree(cubeDir); + cubeDir = psSphereToCube(apparent); + printf(" ><>x, cubeDir->y, cubeDir->z); + } + psFree(cubeDir); psFree(apparent); psFree(actual); @@ -159,7 +182,7 @@ { psTime *empty = NULL; - psTime *time = psTimeAlloc(PS_TIME_TAI); - time->sec = 1131579114; - time->nsec = 498489000; + psTime *time = psTimeAlloc(PS_TIME_UTC); + time->sec = timesec; + time->nsec = 0; time->leapsecond = false; @@ -189,10 +212,18 @@ return 3; } else { - printf("Precession Model output = x,y,s = %.8g, %.8g, %.8g\n\n", + printf("Precession Model output = x,y,s = %.8g, %.8g, %.8g\n", pmodel->x, pmodel->y, pmodel->s); + double x, y, s; + x = 2.857175590089105e-4; + y = 2.3968739377734732e-5; + s = -1.3970066457904322e-8; + printf("Expected output = x,y,s = %.13g, %.13g, %.13g\n\n", x, y, s); } psFree(pmodel); - + //XXX: MAY Be an Error here... Time for Precession Correction has to be TAI or seg faults + psTime *time2 = psTimeAlloc(PS_TIME_UTC); + *time2 = *time; + time2 = psTimeConvert(time2, PS_TIME_TAI); //Tests for Precession Correction function// //Return NULL for NULL time input @@ -206,5 +237,5 @@ //Return NULL for Invalid IERS table - pcorr = psEOC_PrecessionCorr(time, 3); + pcorr = psEOC_PrecessionCorr(time2, 3); if (pcorr != NULL) { psError(PS_ERR_BAD_PARAMETER_VALUE, false, @@ -215,5 +246,5 @@ //Check values from IERS table A - pcorr = psEOC_PrecessionCorr(time, PS_IERS_A); + pcorr = psEOC_PrecessionCorr(time2, PS_IERS_A); if ( pcorr == NULL ) { psError(PS_ERR_BAD_PARAMETER_NULL, false, @@ -227,5 +258,5 @@ //Check values from IERS table B - pcorr = psEOC_PrecessionCorr(time, PS_IERS_B); + pcorr = psEOC_PrecessionCorr(time2, PS_IERS_B); if ( pcorr == NULL ) { psError(PS_ERR_BAD_PARAMETER_NULL, false, @@ -235,8 +266,14 @@ printf("Precession Correction output (IERSB) = x,y,s = %.8g, %.8g, %.8g\n", pcorr->x, pcorr->y, pcorr->s); + double xx, yy, ss; + xx = 0.06295703125; + yy = -0.0287618408203125; + ss = 0.0; + printf("Expected output = x,y,s = %.13g, %.13g, %.13g\n\n", xx, yy, ss); } psFree(pcorr); psFree(time); + psFree(time2); if (!p_psEOCFinalize() ) { psError(PS_ERR_BAD_PARAMETER_VALUE, false, "EOC failed finalization!\n"); @@ -250,7 +287,10 @@ { // psTime *in = psTimeGetNow(PS_TIME_TAI); - psTime *in = psTimeAlloc(PS_TIME_TAI); - in->sec = 1131579114; - in->nsec = 498489000; + // 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; @@ -284,11 +324,42 @@ return 4; } - - psSphere *interm = p_psTimeGetPoleCoords(in); + 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(" -- 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); - printf(" -- PolarMotion = x=%.8g, y=%.8g, s=%.8g\n", polarMotion->x, + + + 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); + printf(" -- NutationCorr = x,y,s = %.13g, %.13g, %.13g\n", nutationCorr->x, + nutationCorr->y, nutationCorr->s); + polarMotion->x += nutationCorr->x; + polarMotion->y += nutationCorr->y; + polarMotion->s += nutationCorr->s; + printf("\n -- Output Earthpole = x,y,s = %.13g, %.13g, %.13g\n", polarMotion->x, + polarMotion->y, polarMotion->s); + double x,y,s; + x = -0.13275353774074533; + y = 0.4359436319739848; + s = -4.2376965863576153e-10; + printf(" -- Expected Earthpole = x,y,s = %.13g, %.13g, %.13g\n", x, y, s); + + 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); @@ -345,12 +416,18 @@ psS32 testSphereRot_TEOtoCEO(void) { - psTime *now = psTimeAlloc(PS_TIME_UT1); - now->sec = 1128530000; - now->nsec = 931154510; - psSphereRot *teoceo = psSphereRot_TEOtoCEO(now); + // psTime *now = psTimeAlloc(PS_TIME_UT1); + // now->sec = 1128530000; + // now->nsec = 931154510; + psTime *time = psTimeAlloc(PS_TIME_UTC); + time->sec = timesec; + time->nsec = 0; + time->leapsecond = false; + // psSphereRot *teoceo = psSphereRot_TEOtoCEO(now); + psSphereRot *teoceo = psSphereRot_TEOtoCEO(time); printf("Output sphere rotation = %lf,%lf,%lf,%lf\n", teoceo->q0, teoceo->q1, teoceo->q2, teoceo->q3); - psFree(now); + // psFree(now); + psFree(time); psFree(teoceo); return 0; @@ -363,7 +440,19 @@ psSphereRot *rot = NULL; - in->x = M_PI / 4.0; - in->y = M_PI / 6.0; - in->s = M_PI / 8.0; + // in->x = M_PI / 4.0; + // in->y = M_PI / 6.0; + // in->s = M_PI / 8.0; + + in->x = 2.857175590089105e-4; + in->y = 2.3968739377734732e-5; + in->s = -1.3970066457904322e-8; + // in->x += 0.06295703125; + // in->y += -0.0287618408203125; + // in->s += 0.0; + double q0,q1,q2,q3; + q0 = -1.1984522406756289e-5; + q1 = 1.4285893358610674e-4; + q2 = 1.2191193518914336e-10; + q3 = -0.9999999897238481; rot = psSphereRot_CEOtoGCRS(empty); @@ -380,4 +469,7 @@ return 2; } + printf("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(); @@ -399,7 +491,10 @@ psSphereRot *rot = NULL; - in->x = M_PI / 4.0; - in->y = M_PI / 6.0; - in->s = M_PI / 8.0; + // in->x = M_PI / 4.0; + // in->y = M_PI / 6.0; + // in->s = M_PI / 8.0; + in->x = -0.13275353774074533; + in->y = 0.4359436319739848; + in->s = -4.2376965863576153e-10; rot = psSphereRot_ITRStoTEO(empty); @@ -416,4 +511,12 @@ return 2; } + printf("Output sphere rotation = %.13g,%.13g,%.13g,%.13g\n", + rot->q0, rot->q1, rot->q2, rot->q3); + double q0,q1,q2,q3; + q0 = -1.0567571848664005e-6; + q1 = 3.218036562931509e-7; + q2 = -3.3580195807204483e-12; + q3 = -0.9999999999993899; + printf("Expected sphere rotation = %.13g,%.13g,%.13g,%.13g\n", q0,q1,q2,q3); psSphere *test = psSphereAlloc(); Index: /trunk/psLib/test/astro/verified/tst_psEarthOrientation.stdout =================================================================== --- /trunk/psLib/test/astro/verified/tst_psEarthOrientation.stdout (revision 5641) +++ /trunk/psLib/test/astro/verified/tst_psEarthOrientation.stdout (revision 5642) @@ -1,4 +1,5 @@ -apparent = r,d = 0.65340071, 0.2463343 +apparent = r,d = 2.1452242, 0.84754694 + ><>