Index: /trunk/psLib/src/astro/psEarthOrientation.c =================================================================== --- /trunk/psLib/src/astro/psEarthOrientation.c (revision 5656) +++ /trunk/psLib/src/astro/psEarthOrientation.c (revision 5657) @@ -8,6 +8,6 @@ * @author Robert Daniel DeSonia, MHPCC * - * @version $Revision: 1.20 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-30 23:50:40 $ + * @version $Revision: 1.21 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-12-02 04:40:00 $ * * Copyright 2005 Maui High Performance Computing Center, University of Hawaii @@ -302,16 +302,7 @@ // calculating the apparent angle from the actual angle and the sun position - // first, calculate the angle between the sun vector and the actual vector - - // Moving to cartesian first: XXX -- is this required? psCube* sunVector = psSphereToCube(sun); psCube* actualVector = psSphereToCube(actual); - printf("\n Actual sphere = %.13g, %.13g", actual->r, actual->d); - printf("\n Actual cube = %.13g, %.13g, %.13g \n", actualVector->x, actualVector->y, actualVector->z); - printf("\n SUN sphere = %.13g, %.13g", sun->r, sun->d); - printf("\n SUN cube = %.13g, %.13g, %.13g \n", sunVector->x, sunVector->y, sunVector->z); - psSphere *new = psCubeToSphere(actualVector); - printf(" New sphere = %.13g, %.13g", new->r, new->d); - psFree(new); + // use dot product to calculate the angle of separation // N.B., assuming the psSphereToCube function returns a unit vector. @@ -320,9 +311,7 @@ sunVector->y*actualVector->y + sunVector->z*actualVector->z); - printf("\n Theta = %.13g ", theta); double r0 = PS_AU * tan(theta); - printf(" r0 = %.13g", r0); double deflection = 4.0*PS_G*PS_M/(PS_C0*PS_C0*r0); - printf(" deflection = %.13g \n", deflection); + // make sure the deflection is not greater than 1.75 arcsec double limit = SEC_TO_RAD(1.75); @@ -352,8 +341,6 @@ deflection = SEC_TO_RAD(deflection); theta = atan(r0/PS_AU) * tan(deflection); - printf(" deflection = %.13g Theta=%.13g ", deflection, theta); // phi = sqrt( deflection*deflection - theta*theta ); phi = deflection * cos(asin(theta/deflection)); - printf(" Phi = %.13g \n", phi); apparent->r = theta; apparent->d = phi; @@ -726,11 +713,10 @@ in = psTimeConvert(in, PS_TIME_UT1); } - double T = (double)(in->sec) + (double)(in->nsec / 1e9); - printf("\nThe Value of T is = %.19g\n", T); + // double T = (double)(in->sec) + (double)(in->nsec / 1e9); + double T = psTimeToJD(in); T += -2451545.0; double theta = 2.0 * M_PI * (0.7790572732640 + 1.00273781191135448 * T); - printf("\nThe Value of theta is = %.13g\n", theta); + psSphereRot *out = psSphereRotAlloc(theta, 0.0, 0.0); // psSphereRot *out = psSphereRotInvert(theta, 0.0, 0.0); - psSphereRot *out = psSphereRotInvert(theta, 0.0, 0.0); psFree(in); Index: /trunk/psLib/test/astro/tst_psEarthOrientation.c =================================================================== --- /trunk/psLib/test/astro/tst_psEarthOrientation.c (revision 5656) +++ /trunk/psLib/test/astro/tst_psEarthOrientation.c (revision 5657) @@ -5,6 +5,6 @@ * @author d-Rob, MHPCC * -* @version $Revision: 1.17 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-11-30 23:50:40 $ +* @version $Revision: 1.18 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-12-02 04:40:00 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -56,4 +56,17 @@ #define timesec 1049160600 +#define objR DEG_TO_RAD(122.9153182445501) +#define objD DEG_TO_RAD(48.562968978679194) +static psSphere *obj = NULL; +static void objSetup(void); + +static void objSetup(void) +{ + if (obj == NULL) { + obj = psSphereAlloc(); + obj->r = objR; + obj->d = objD; + } +} psS32 testAberration(void) @@ -86,6 +99,5 @@ 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 */ + double c = 299792458.0; /* m/s */ speed = speed / c; @@ -104,12 +116,29 @@ apparent = psAberration(apparent, actual, direction, speed); - if (apparent == NULL) - printf("\nApparent is NULL\n"); - else { - printf("\napparent = r,d = %.8g, %.8g\n", apparent->r, apparent->d); - psFree(cubeDir); - cubeDir = psSphereToCube(apparent); - printf(" ><>x, cubeDir->y, cubeDir->z); - } + + 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); + + double x, y, z; + x = -0.35963388069046304; + y = 0.5555192509816625; + 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 ) { + 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) ); + return 3; + } + + psFree(outCube); + psFree(cubeDir); psFree(apparent); @@ -138,11 +167,5 @@ sunCube->y = 2.5880956908748722e10; sunCube->z = 1.1220046291457653e10; - // double sum = sqrt(sunCube->x*sunCube->x + sunCube->y*sunCube->y + sunCube->z*sunCube->z); - // sunCube->x = sunCube->x/sum; - // sunCube->y = sunCube->y/sum; - // sunCube->z = sunCube->z/sum; - // printf("\n <<<x, sunCube->y, sunCube->z); psSphere *sun = psCubeToSphere(sunCube); - printf("\n <<<r, sun->d); empty = psGravityDeflection(apparent, empty, sun); @@ -161,11 +184,33 @@ apparent = psGravityDeflection(NULL, actual, sun); 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)); psCube *outCube = psSphereToCube(result); - printf(" -- outCube = x,y,z = %.13g, %.13g, %.13g -- \n", outCube->x, outCube->y, outCube->z); + 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", + outCube2->x, outCube2->y, outCube2->z); + + double x, y, z; + x = -0.35961949760293604; + y = 0.5555613950298085; + 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; + } + psFree(outCube); - psCube *outCube2 = psSphereToCube(apparent); - printf(" -- outCube = x,y,z = %.13g, %.13g, %.13g -- \n", outCube2->x, outCube2->y, outCube2->z); psFree(outCube2); psFree(sunCube); @@ -416,17 +461,67 @@ psS32 testSphereRot_TEOtoCEO(void) { + psSphereRot *rot = NULL; + psTime *empty = NULL; // psTime *now = psTimeAlloc(PS_TIME_UT1); // now->sec = 1128530000; // now->nsec = 931154510; - psTime *time = psTimeAlloc(PS_TIME_UTC); - time->sec = timesec; - time->nsec = 0; + psTime *time = psTimeAlloc(PS_TIME_UT1); + time->sec = timesec-1; + time->nsec = 657017200; time->leapsecond = false; + // psSphereRot *teoceo = psSphereRot_TEOtoCEO(now); + + //return NULL for NULL input time + rot = psSphereRot_TEOtoCEO(empty); + if (rot != NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, + "psSphereRot_TEOtoCEO failed to return NULL for NULL input time.\n"); + return 1; + } + psSphereRot *teoceo = psSphereRot_TEOtoCEO(time); - printf("Output sphere rotation = %lf,%lf,%lf,%lf\n", + //Make sure values match for other psTime type + empty = psTimeAlloc(PS_TIME_UTC); + empty->sec = timesec; + empty->nsec = 0; + empty->leapsecond = false; + rot = psSphereRot_TEOtoCEO(empty); + if (fabs(rot->q0-teoceo->q0) > DBL_EPSILON || fabs(rot->q1-teoceo->q1) > DBL_EPSILON || + fabs(rot->q2-teoceo->q2) > DBL_EPSILON || fabs(rot->q3-teoceo->q3) > DBL_EPSILON) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, + "psSphereRot_TEOtoCEO failed to return matching values for different time types.\n"); + return 2; + } + printf("\nOutput sphere rotation = %.13g, %.13g, %.13g, %.13g\n\n", teoceo->q0, teoceo->q1, teoceo->q2, teoceo->q3); + objSetup(); + psSphereRot *earthRot = psSphereRotConjugate(NULL, teoceo); + psSphere *result = psSphereRotApply(NULL, earthRot, obj); + // psCube *cube = psSphereToCube(obj); + psCube *cube = psSphereToCube(result); + + + double x, y, z; + x = 0.01698625430807123; + y = -0.6616523084626379; + z = 0.7496183628158023; + if ( fabs(x-cube->x) > DBL_EPSILON || fabs(y-cube->y) > DBL_EPSILON || + fabs(z-cube->z) > DBL_EPSILON) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, + "psSphereRot_TEOtoCEO returned incorrect values.\n"); + printf("\nOutput cube = x,y,z = %.13g, %.13g, %.13g\n", + cube->x, cube->y, cube->z); + printf("Expected cube = x,y,z = %.13g, %.13g, %.13g\n", x, y, z); + printf("A difference of: %.13g, %.13g, %.13g\n\n", + (x-cube->x), (y-cube->y), (z-cube->z)); + return 3; + } // psFree(now); + psFree(result); + psFree(obj); + psFree(earthRot); + psFree(cube); psFree(time); psFree(teoceo); Index: /trunk/psLib/test/astro/verified/tst_psEarthOrientation.stderr =================================================================== --- /trunk/psLib/test/astro/verified/tst_psEarthOrientation.stderr (revision 5656) +++ /trunk/psLib/test/astro/verified/tst_psEarthOrientation.stderr (revision 5657) @@ -9,6 +9,8 @@