Index: trunk/psLib/test/astro/tst_psEarthOrientation.c =================================================================== --- trunk/psLib/test/astro/tst_psEarthOrientation.c (revision 5824) +++ trunk/psLib/test/astro/tst_psEarthOrientation.c (revision 5969) @@ -5,6 +5,6 @@ * @author d-Rob, MHPCC * -* @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-12-21 06:15:58 $ +* @version $Revision: 1.25 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-01-12 20:40:13 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -64,4 +64,5 @@ static psSphere *obj = NULL; static void objSetup(void); +static double greatCircle(psSphere *in1, psSphere *in2); static void objSetup(void) @@ -80,23 +81,51 @@ } +//returns the great circle ANGULAR distance between 2 positions (psSphere's) +static double greatCircle(psSphere *in1, psSphere *in2) +{ + double r1, r2, d1, d2; + r1 = in1->r; + r2 = in2->r; + d1 = in1->d; + d2 = in2->d; + double out = 0.0; + double c1, c2, cd, s1, s2, sum, ac; + c1 = cos(r1); + c2 = cos(r2); + cd = cos(d1-d2); + s1 = sin(r1); + s2 = sin(r2); + sum = c1*c2*cd + s1*s2; + if (sum > 1.0) + sum = 1.0 - (sum - 1.0); + + ac = acos(sum); + printf("\n c1=%lf, c2=%lf, cd=%lf, s1=%lf, s2=%lf, sum=%.19g, ac=%g\n", c1,c2,cd,s1,s2,sum,ac); + // out = acos(cos(r1)*cos(r2)*cos(d1-d2) + sin(r1)*sin(r2)); + // if (out != 0.0) printf("\n in greatCircle, out = %lf\n", out); + return out; +} + psS32 testAberration(void) { psSphere *apparent = NULL; - psSphere *actual = psSphereAlloc(); + // psSphere *actual = 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(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); + psCube *actualCube = psCubeAlloc(); + // actual->r = DEG_TO_RAD(122.9153182445501); + // actual->d = DEG_TO_RAD(48.562968978679194); + // actualCube->x = -0.3596195125758298; + // actualCube->y = 0.5555613903455866; + // actualCube->z = 0.7496834983724809; + + //values from After gravity deflection// + actualCube->x = -0.35961949760293604; + actualCube->y = 0.5555613950298085; + actualCube->z = 0.7496835020836093; + psSphere *actual = psCubeToSphere(actualCube); + psCube *cubeDir = psCubeAlloc(); cubeDir->x = 5148.713262821658; @@ -106,9 +135,15 @@ cubeDir->y += 248.46429758174693; cubeDir->z += 0.09694774143797581; + double length = sqrt(cubeDir->x*cubeDir->x+cubeDir->y*cubeDir->y+cubeDir->z*cubeDir->z); + cubeDir->x = cubeDir->x/length; + cubeDir->y = cubeDir->y/length; + cubeDir->z = cubeDir->z/length; 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 = 299792458.0; /* m/s */ - speed = speed / c; + // speed = speed / c; + speed = length / c; empty = psAberration(empty, apparent, direction, speed); @@ -127,10 +162,10 @@ apparent = psAberration(apparent, actual, direction, speed); - printf("\nSphere Difference = r,d = %.13g, %.13g\n", - (actual->r - apparent->r), (actual->d - apparent->d)); + // printf("\nSphere Difference = r,d = %.13g, %.13g\n", + // (actual->r - apparent->r), (actual->d - apparent->d)); psCube *outCube = psSphereToCube(apparent); printf(" -- resultCube = x,y,z = %.13g, %.13g, %.13g -- \n", outCube->x, outCube->y, outCube->z); - + //expected cube values double x, y, z; x = -0.35963388069046304; @@ -138,16 +173,40 @@ z = 0.7497078321908413; - printf(" -- expectedCube = x,y,z = %.13g, %.13g, %.13g -- \n\n", x, y, z); - - if ( fabs(x - outCube->x) > DBL_EPSILON || fabs(y - outCube->y) > DBL_EPSILON || - fabs(z - outCube->z) > DBL_EPSILON ) { + printf(" -- expectedCube = x,y,z = %.13g, %.13g, %.13g -- \n", x, y, z); + printf("Cube Difference = x,y,z = %.13g, %.13g, %.13g \n\n", + (x - outCube->x), (y - outCube->y), (z - outCube->z) ); + psFree(actual); + // psFree(actualCube); + actualCube->x = x; + actualCube->y = y; + actualCube->z = z; + actual = psCubeToSphere(actualCube); + double xxx = greatCircle(actual, apparent); + printf(" The great circle angular distance between expected & result = %.13g\n", + xxx); + + if ( fabs(x - outCube->x) > FLT_EPSILON || fabs(y - outCube->y) > FLT_EPSILON || + fabs(z - outCube->z) > FLT_EPSILON ) { psError(PS_ERR_BAD_PARAMETER_VALUE, false, "psAberration returned incorrect values.\n"); - printf("Cube Difference = x,y,z = %.13g, %.13g, %.13g \n\n", - (x - outCube->x), (y - outCube->y), (z - outCube->z) ); + printf("\nMagnitude of expected change = x,y,z = %.13g, %.13g, %.13g \n", + (actualCube->x - x), (actualCube->y - y), (actualCube->z - z) ); + printf("Magnitude of actual change = x,y,z = %.13g, %.13g, %.13g \n", + (actualCube->x - outCube->x), (actualCube->y - outCube->y), + (actualCube->z - outCube->z) ); + psFree(actual); + // psFree(actualCube); + actualCube->x = x; + actualCube->y = y; + actualCube->z = z; + actual = psCubeToSphere(actualCube); + x = greatCircle(actual, apparent); + printf(" The great circle angular distance between expected & result = %.13g\n", + x); return 3; } psFree(outCube); + psFree(actualCube); psFree(cubeDir); @@ -162,15 +221,14 @@ { - psSphere *actual = psSphereAlloc(); + // psSphere *actual = psSphereAlloc(); psSphere *apparent = NULL; // psSphere *sun = psSphereAlloc(); psSphere *empty = NULL; - // sun->r = 0.2; - // sun->d = 0.2; - // actual->r = 0.61001; - // actual->d = 0.01999; - actual->r = DEG_TO_RAD(122.9153182445501); - actual->d = DEG_TO_RAD(48.562968978679194); + psCube *actualCube = psCubeAlloc(); + actualCube->x = -0.3596195125758298; + actualCube->y = 0.5555613903455866; + actualCube->z = 0.7496834983724809; + psSphere *actual = psCubeToSphere(actualCube); psCube *sunCube = psCubeAlloc(); sunCube->x = 1.467797790127511e11; @@ -193,17 +251,17 @@ apparent = psGravityDeflection(NULL, actual, sun); + psSphere *result2 = psSphereSetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN); + apparent->r *= -1.0; + apparent->d *= -1.0; psSphere *result = psSphereSetOffset(actual, apparent, PS_SPHERICAL, PS_RADIAN); - - printf("\nActual r,d = %.13g,%.13g Apparent r,d = %.16g, %.16g \n", - actual->r, actual->d, result->r, result->d); - printf("Sphere Difference = r,d = %.13g, %.13g\n", - (actual->r - result->r), (actual->d - result->d)); + // psSphere *result = psSphereGetOffset(apparent, actual, PS_SPHERICAL, PS_RADIAN); + printf(" -- actualCube = x,y,z = %.13g, %.13g, %.13g -- \n", + actualCube->x, actualCube->y, actualCube->z); psCube *outCube = psSphereToCube(result); printf(" -- resultCube = x,y,z = %.13g, %.13g, %.13g -- \n", outCube->x, outCube->y, outCube->z); - psCube *outCube2 = psSphereToCube(actual); - printf(" -- actualCube = x,y,z = %.13g, %.13g, %.13g -- \n", + psCube *outCube2 = psSphereToCube(result2); + printf(" -- resultCube2= x,y,z = %.13g, %.13g, %.13g -- \n", outCube2->x, outCube2->y, outCube2->z); - double x, y, z; x = -0.35961949760293604; @@ -211,19 +269,33 @@ z = 0.7496835020836093; - printf(" -- expectedCube = x,y,z = %.13g, %.13g, %.13g -- \n\n", x, y, z); - - if ( fabs(x - outCube->x) > DBL_EPSILON || fabs(y - outCube->y) > DBL_EPSILON || - fabs(z - outCube->z) > DBL_EPSILON ) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, - "psGravityDeflection returned incorrect values.\n"); - printf("Cube Difference = x,y,z = %.13g, %.13g, %.13g \n", - (x - outCube->x), (y - outCube->y), (z - outCube->z) ); - return 1; - } + printf(" -- expectCube = x,y,z = %.13g, %.13g, %.13g -- \n\n", x, y, z); + + // if ( fabs(x - outCube->x) > DBL_EPSILON || fabs(y - outCube->y) > DBL_EPSILON || + // fabs(z - outCube->z) > DBL_EPSILON ) { + // psError(PS_ERR_BAD_PARAMETER_VALUE, false, + // "psGravityDeflection returned incorrect values.\n"); + printf("expect-actual= x,y,z = %.11g, %.11g, %.11g \n", + (x - actualCube->x), (y - actualCube->y), (z - actualCube->z) ); + printf("result-actual= x,y,z = %.11g, %.11g, %.11g \n", + (outCube->x - actualCube->x), (outCube->y - actualCube->y), (outCube->z - actualCube->z) ); + printf("expect-result= x,y,z = %.11g, %.11g, %.11g \n", + (x - outCube->x), (y - outCube->y), (z - outCube->z) ); + // return 1; + // } + psFree(result2); + outCube2->x = x; + outCube2->y = y; + outCube2->z = z; + result2 = psCubeToSphere(outCube2); + psFree(result); + result = psSphereGetOffset(actual, result2, PS_SPHERICAL, PS_RADIAN); + printf("The apparent output sphere = r,d = %.13g, %.13g\n", apparent->r, apparent->d); + printf("The expected output sphere = r,d = %.13g, %.13g\n\n", result->r, result->d); + psFree(result2); psFree(outCube); psFree(outCube2); psFree(sunCube); - + psFree(actualCube); psFree(result); psFree(actual); @@ -275,11 +347,11 @@ psError(PS_ERR_BAD_PARAMETER_VALUE, false, " psEOC_PrecessionModel return incorrect values.\n"); - printf("\n Precession Model output = x,y,s = %.8g, %.8g, %.8g\n", - pmodel->x, pmodel->y, pmodel->s); - printf(" Expected output = x,y,s = %.13g, %.13g, %.13g\n", x, y, s); - printf(" A difference of: %.13g, %.13g, %.13g\n\n", - (pmodel->x - x), (pmodel->y - y), (pmodel->s - s) ); return 4; } + printf("\n Precession Model output = x,y,s = %.8g, %.8g, %.8g\n", + pmodel->x, pmodel->y, pmodel->s); + printf(" Expected output = x,y,s = %.13g, %.13g, %.13g\n", x, y, s); + printf(" A difference of: %.13g, %.13g, %.13g\n\n", + (pmodel->x - x), (pmodel->y - y), (pmodel->s - s) ); } psFree(pmodel); @@ -348,22 +420,34 @@ yy = -0.0287618408203125; ss = 0.0; - xx = SEC_TO_RAD(xx) * 1e-6; - yy = SEC_TO_RAD(yy) * 1e-6; + xx = SEC_TO_RAD(xx) * 1e-3; + yy = SEC_TO_RAD(yy) * 1e-3; // 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", + printf("PrecessionCorr 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) ); + 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; // } } + + double xCorr, yCorr; + xCorr = 3.05224300720406e-10; + yCorr = -1.39441339235822e-10; + //pcorr is the *expected* output from PrecessionModel// pcorr->x = 2.857175590089105e-4; pcorr->y = 2.3968739377734732e-5; pcorr->s = -1.3970066457904322e-8; + pcorr->x += xCorr; + pcorr->y += yCorr; + psEarthPole *precess = psEOC_PrecessionModel(time2); + precess->x += xCorr; + precess->y += yCorr; + psSphereRot *precessNutInv = psSphereRot_CEOtoGCRS(precess); + psSphereRot *precessNut = psSphereRotConjugate(NULL, precessNutInv); // pcorr->x += 3.05224300720406e-7; // pcorr->y += -1.39441339235822e-7; @@ -377,37 +461,61 @@ printf("\n Error at CEOtoGCRS, output psSphereRot doesn't match expected.\n"); } - printf("\n Output sphere rotation = %.13g,%.13g,%.13g,%.13g\n", + printf(" Output sphere rotation = %.13g,%.13g,%.13g,%.13g\n", pni->q0, pni->q1, pni->q2, pni->q3); printf(" Expected sphere rotation = %.13g,%.13g,%.13g\n", q0,q1,q2); + printf(" MY sphere rotation = %.13g,%.13g,%.13g\n", + precessNutInv->q0, precessNutInv->q1, precessNutInv->q2); psCube *objC = psCubeAlloc(); // objC->x = -3.5963388069046304; // objC->y = 0.5555192509816625; // objC->z = 0.7497078321908413; - objSetup(); + // objSetup(); + + //This is the sphere rotation for the *expected* precession output// psSphereRot *pn = psSphereRotConjugate(NULL, pni); - // psSphere *sphere = psCubeToSphere(objC); - psSphere *sphere = psSphereAlloc(); - *sphere = *obj; - psFree(obj); - psSphere *result = psSphereRotApply(NULL, pn, sphere); - objC->x = -0.3598480726985338; - objC->y = 0.5555012823608123; - objC->z = 0.7496183628158023; - psFree(sphere); - sphere = psCubeToSphere(objC); + + // psSphere *sphere = psSphereAlloc(); + // *sphere = *obj; + // psFree(obj); + + //create a psSphere for (from) the start position given in eoc_testing// + objC->x = -0.35963388069046304; + objC->y = 0.5555192509816625; + objC->z = 0.7497078321908413; + psSphere *sphere = psCubeToSphere(objC); + + psSphere *expect = psSphereRotApply(NULL, pn, sphere); + //expected results below - stored in: sphere // + double x,y,z; + x = -0.3598480726985338; + y = 0.5555012823608123; + z = 0.7496183628158023; + printf("\n<r, result->d, sphere->r, sphere->d); - double xx = acos(cos(result->r)*cos(sphere->r)*cos(result->d - sphere->d) + sin(result->r)*sin(sphere->r)); - printf("GREAT CIRCLE DIFFERENCE = %.13g \n\n", xx); + objC = psSphereToCube(expect); + printf("<x, objC->y, objC->z); + printf(" Difference = %.13g, %.13g, %.13g\n", objC->x-x, objC->y-y, objC->z-z); + x = objC->x; + y = objC->y; + z = objC->z; + psSphere *result = psSphereRotApply(NULL, precessNut, sphere); + psFree(objC); + objC = psSphereToCube(result); + printf("<x, objC->y, objC->z); + printf(" Difference = %.13g, %.13g, %.13g\n\n", objC->x-x, objC->y-y, objC->z-z); + + double xx = greatCircle(result, expect); + printf("GREAT CIRCLE DIFFERENCE = %.13g \n", xx); + psFree(precess); + psFree(precessNut); + psFree(precessNutInv); + psFree(expect); + psFree(objC); psFree(sphere); psFree(result); psFree(pn); - // psFree(obj); psFree(pni); - psFree(pcorr); - psFree(time2); if (!p_psEOCFinalize() ) { @@ -443,10 +551,33 @@ } + //Test for IERS bulletin A. + double x, y, s; + x = -6.454389659777e-07; + y = 2.112606414597e-06; + s = 0.0; + 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; + } + psFree(polarMotion); + + //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; + x = -6.45381397904e-07; + y = 2.112819726698e-06; + s = 0.0; if (polarMotion == NULL) { @@ -465,27 +596,4 @@ // 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; - } - @@ -495,19 +603,11 @@ 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; + x = -6.43607313124045e-7; + y = 2.11351436973568e-6; + s = -7.39617581324646e-12; if ( fabs(polarMotion->x - x) > DBL_EPSILON || fabs(polarMotion->y - y) > DBL_EPSILON @@ -525,4 +625,12 @@ } + psEarthPole *polarTide = psEOC_PolarTideCorr(in); + polarMotion->x += polarTide->x; + polarMotion->y += polarTide->y; + polarMotion->s += polarTide->s; + psFree(polarTide); + printf("\n PolarMotion + PolarTideCorr out = x,y,s = %.13g, %.13g, %.13g\n", + polarMotion->x, polarMotion->y, polarMotion->s); + if (!p_psEOCFinalize() ) { psError(PS_ERR_BAD_PARAMETER_VALUE, false, "EOC failed finalization!\n"); @@ -602,5 +710,5 @@ return 3; } else { - printf("Nutation Correction output = x,y,s = %.8g, %.8g, %.8g\n\n", + printf("Nutation Correction output = x,y,s = %.13g, %.13g, %.13g\n\n", nute->x, nute->y, nute->s); } @@ -656,6 +764,6 @@ objSetup(); psSphereRot *earthRot = psSphereRotConjugate(NULL, teoceo); - // psSphere *result = psSphereRotApply(NULL, earthRot, obj); - psSphere *result = psSphereRotApply(NULL, teoceo, obj); + psSphere *result = psSphereRotApply(NULL, earthRot, obj); + // psSphere *result = psSphereRotApply(NULL, teoceo, obj); psCube *cube = psSphereToCube(result); @@ -773,6 +881,6 @@ psSphere *expected = psCubeToSphere(expect); psFree(expect); - double d = (6.38e6) * acos(cos(result->r)*cos(expected->r)*cos(result->d-expected->d) + sin(result->r)*sin(expected->r)); - printf("\n\nGreat circle difference of: %.13g \n", d); + double d = acos(cos(result->r)*cos(expected->r)*cos(result->d-expected->d) + sin(result->r)*sin(expected->r)); + printf("Great circle difference of: %.13g \n", d); psFree(expected); @@ -781,9 +889,9 @@ - psSphere *test = psSphereAlloc(); - test->r = 0.0; - test->d = 0.0; - test = psSphereRotApply(test, rot, test); - printf("\n Sphere -test- has values r,d = %.8g, %.8g \n", test->r, test->d); + // psSphere *test = psSphereAlloc(); + // test->r = 0.0; + // test->d = 0.0; + // test = psSphereRotApply(test, rot, test); + // printf("\n Sphere -test- has values r,d = %.8g, %.8g \n", test->r, test->d); psFree(precessionNutation); @@ -791,5 +899,5 @@ psFree(cube); - psFree(test); + // psFree(test); psFree(rot); psFree(in); @@ -875,4 +983,12 @@ z = 0.7496169753347885; printf("\n Cube -expected- has x,y,z = %.13g, %.13g, %.13g \n", x, y, z ); + temp->x = x; + temp->y = y; + temp->z = z; + psSphere *sphere = psCubeToSphere(temp); + double d = greatCircle(sphere, test); + + printf("Great circle difference of: %.13g \n", d); + psFree(sphere); psFree(newRot);