Index: /trunk/psLib/src/astro/psCoord.c =================================================================== --- /trunk/psLib/src/astro/psCoord.c (revision 5625) +++ /trunk/psLib/src/astro/psCoord.c (revision 5626) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.94 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-11-30 02:00:00 $ +* @version $Revision: 1.95 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-11-30 02:17:17 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -1137,8 +1137,23 @@ // sphere->rErr = arctan(cube->xErr/cube->yErr); // sphere->dErr = arctan(sqrt(cube->xErr*cube->xErr + cube->yErr*cube->yErr)/cube->zErr); - sphere->r = 1 / (atan(cube->x/cube->y)); - sphere->d = 1 / (atan(sqrt(cube->x*cube->x + cube->y*cube->y)/cube->z)); - sphere->rErr = 1 / (atan(cube->xErr/cube->yErr)); - sphere->dErr = 1 / (atan(sqrt(cube->xErr*cube->xErr + cube->yErr*cube->yErr)/cube->zErr)); + // sphere->r = 1 / (atan(cube->x/cube->y)); + // sphere->d = 1 / (atan(sqrt(cube->x*cube->x + cube->y*cube->y)/cube->z)); + // sphere->rErr = 1 / (atan(cube->xErr/cube->yErr)); + // sphere->dErr = 1 / (atan(sqrt(cube->xErr*cube->xErr + cube->yErr*cube->yErr)/cube->zErr)); + psCube *cube2 = psCubeAlloc(); + *cube2 = *cube; + double mag = sqrt(cube->x*cube->x + cube->y*cube->y + cube->z*cube->z); + if (mag > 1.0) { + cube2->x = cube2->x/mag; + cube2->y = cube2->y/mag; + cube2->z = cube2->z/mag; + } + + sphere->r = atan2(cube2->y, cube2->x); + sphere->d = asin(cube2->z); + // sphere->d = atan2((cube->x*cube->x + cube->y*cube->y), cube->z); + sphere->rErr = atan2(cube2->yErr, cube2->xErr); + sphere->dErr = asin(cube2->zErr); + psFree(cube2); return(sphere); Index: /trunk/psLib/src/astro/psEarthOrientation.c =================================================================== --- /trunk/psLib/src/astro/psEarthOrientation.c (revision 5625) +++ /trunk/psLib/src/astro/psEarthOrientation.c (revision 5626) @@ -8,6 +8,6 @@ * @author Robert Daniel DeSonia, MHPCC * - * @version $Revision: 1.18 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-17 03:59:05 $ + * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-11-30 02:17:17 $ * * Copyright 2005 Maui High Performance Computing Center, University of Hawaii @@ -307,5 +307,11 @@ 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. @@ -314,8 +320,9 @@ 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); @@ -345,6 +352,8 @@ deflection = SEC_TO_RAD(deflection); theta = atan(r0/PS_AU) * tan(deflection); - phi = sqrt( deflection*deflection - theta*theta ); - // phi = deflection * cos(asin(theta/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; @@ -950,5 +959,4 @@ Y = psPolynomial1DEval(yPoly,t); for (int i = 0; i < 10; i++) { - // tj = SEC_TO_RAD(pow(t, 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]; @@ -960,5 +968,5 @@ pole->x = X; pole->y = Y; - pole->s = SEC_TO_RAD(4.7e-5) * t; //XXX: This conv. should include t? + pole->s = SEC_TO_RAD(4.7e-5) * t; return pole; Index: /trunk/psLib/test/astro/tst_psEarthOrientation.c =================================================================== --- /trunk/psLib/test/astro/tst_psEarthOrientation.c (revision 5625) +++ /trunk/psLib/test/astro/tst_psEarthOrientation.c (revision 5626) @@ -5,6 +5,6 @@ * @author d-Rob, MHPCC * -* @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-11-17 04:04:30 $ +* @version $Revision: 1.16 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-11-30 02:17:17 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -102,11 +102,24 @@ psSphere *actual = psSphereAlloc(); psSphere *apparent = NULL; - psSphere *sun = psSphereAlloc(); + // psSphere *sun = psSphereAlloc(); psSphere *empty = NULL; - sun->r = 0.2; - sun->d = 0.2; - actual->r = 0.61001; - actual->d = 0.01999; + // 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 *sunCube = psCubeAlloc(); + sunCube->x = 1.467797790127511e11; + 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); @@ -123,8 +136,16 @@ } - apparent = psGravityDeflection(apparent, actual, sun); + 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); + psCube *outCube = psSphereToCube(result); + printf(" -- outCube = x,y,z = %.13g, %.13g, %.13g -- \n", outCube->x, outCube->y, outCube->z); + 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); + psFree(result); psFree(actual); Index: /trunk/psLib/test/astro/tst_psSphereOps.c =================================================================== --- /trunk/psLib/test/astro/tst_psSphereOps.c (revision 5625) +++ /trunk/psLib/test/astro/tst_psSphereOps.c (revision 5626) @@ -5,6 +5,6 @@ * @author d-Rob, MHPCC * -* @version $Revision: 1.13 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-11-15 00:04:30 $ +* @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-11-30 02:17:17 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -494,8 +494,8 @@ } */ - output->r *= -1.0; - output->d *= -1.0; + // output->r *= -1.0; + // output->d *= -1.0; empty = psSphereGetOffset(origin, output, PS_LINEAR, PS_RADIAN); - if ( fabs(offset->r - empty->r) > 0.0001 || abs(offset->d - empty->r) > 0.0001 ) { + if ( fabs(offset->r - empty->r) > 0.0001 || fabs(offset->d - empty->d) > 0.0001 ) { psError(PS_ERR_BAD_PARAMETER_VALUE, false, "psSphereGetOffset failed to return correct linear offset values.\n"); Index: /trunk/psLib/test/astro/verified/tst_psEarthOrientation.stdout =================================================================== --- /trunk/psLib/test/astro/verified/tst_psEarthOrientation.stdout (revision 5625) +++ /trunk/psLib/test/astro/verified/tst_psEarthOrientation.stdout (revision 5626) @@ -1,6 +1,19 @@ -apparent = r,d = 1.0900423, 0.75502353 +apparent = r,d = 0.65340071, 0.2463343 -Actual r,d = 0.61001,0.01999 Apparent r,d = 0.6100100000001365, 0.01999000000006505 + <<<