Changeset 5626

Timestamp:

Nov 29, 2005, 4:17:17 PM (21 years ago)

Author:

drobbin

Message:

Edited psCubeToSphere, tests. EOC - gravity still fails. Others still unchecked.

Location:

trunk/psLib

Files:

• src/astro/psCoord.c (modified) (2 diffs)
• src/astro/psEarthOrientation.c (modified) (6 diffs)
• test/astro/tst_psEarthOrientation.c (modified) (3 diffs)
• test/astro/tst_psSphereOps.c (modified) (2 diffs)
• test/astro/verified/tst_psEarthOrientation.stdout (modified) (1 diff)

trunk/psLib/src/astro/psCoord.c

@@ -10 +10 @@
 *  @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 +1137 @@
     //    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);

trunk/psLib/src/astro/psEarthOrientation.c

@@ -8 +8 @@
  *  @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 +307 @@
     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 +320 @@
                     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 +352 @@
     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 +959 @@
     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 +968 @@
     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;

trunk/psLib/test/astro/tst_psEarthOrientation.c

@@ -5 +5 @@
 *  @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 +102 @@
     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 <<<<SUNvec = x,y,z = %.13g, %.13g %.13g", sunCube->x, sunCube->y, sunCube->z);
+    psSphere *sun = psCubeToSphere(sunCube);
+    printf("\n <<<<SUN = r,d = %.13g, %.13g\n", sun->r, sun->d);
 
     empty = psGravityDeflection(apparent, empty, sun);
@@ -123 +136 @@
     }
 
-    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);

trunk/psLib/test/astro/tst_psSphereOps.c

@@ -5 +5 @@
 *  @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 +494 @@
         }
     */
-    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");

trunk/psLib/test/astro/verified/tst_psEarthOrientation.stdout

@@ -1 +1 @@
 
-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 
+ <<<<SUN = r,d = 0.1745310859793, 0.0751383461293
+
+ Actual sphere = 2.145277004504, 0.8475837032218
+ Actual cube = -0.3596195125758, 0.5555613903456, 0.7496834983725 
+
+ SUN sphere = 0.1745310859793, 0.0751383461293
+ SUN cube = 0.9820293796219, 0.1731564131522, 0.0750676637684 
+ New sphere = 2.145277004504, 0.8475837032218
+ Theta = -0.2006809204081    r0 = -30431058247.35   deflection = -1.941413573621e-07 
+        deflection = -9.412238611831e-13  Theta=1.888856707723e-13    Phi = -9.220762225745e-13 
+
+Actual r,d = 2.145277004504,0.8475837032218    Apparent r,d = 2.145277004504245, 0.8475837032208979 
+ -- outCube = x,y,z = -0.3596195125763, 0.5555613903461, 0.7496834983719  -- 
+ -- outCube = x,y,z = 1, 1.888856707723e-13, -9.220762225745e-13  -- 
 Precession Model output = x,y,s = 1.1557751e+08, 8865375.4, -5.1231901e+14