Changeset 5749

Timestamp:

Dec 7, 2005, 4:49:04 PM (21 years ago)

Author:

drobbin

Message:

updated tests, got TEOtoCEO working and verified. Currently working on CEOtoGCRS.

Location:

trunk/psLib

Files:

• src/astro/psEarthOrientation.c (modified) (2 diffs)
• test/astro/tst_psEarthOrientation.c (modified) (7 diffs)

trunk/psLib/src/astro/psEarthOrientation.c

@@ -8 +8 @@
  *  @author Robert Daniel DeSonia, MHPCC
  *
- *  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-02 04:40:00 $
+ *  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-08 02:49:04 $
  *
  *  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
@@ -690 +690 @@
     //XXX: Used formula for rotation matrix D from mathworld for z-axis rotation
     double a =  1.0 / (1.0 + sqrt(1.0 - (pole->x*pole->x + pole->y*pole->y) ) );
+    /*    A[0][0] = (1.0 - a*pole->x*pole->x)*cos(pole->s) - a*pole->x*pole->y*sin(pole->s);
+        A[0][1] = -a*pole->x*pole->y*cos(pole->s) + (1.0 - a*pole->y*pole->y)*sin(pole->s);
+        A[0][2] = pole->x*cos(pole->s) + pole->y*sin(pole->s);
+        A[1][0] = -(1.0 - a*pole->x*pole->x)*sin(pole->s) - a*pole->x*pole->y*cos(pole->s);
+        A[1][1] = a*pole->x*pole->y*sin(pole->s) + (1.0 - a*pole->y*pole->y)*cos(pole->s);
+        A[1][2] = -pole->x*sin(pole->s) + pole->y*cos(pole->s);
+        A[2][0] = -pole->x;
+        A[2][1] = -pole->y;
+        A[2][2] = 1.0 - a*(pole->x*pole->x + pole->y*pole->y);
+    */
     A[0][0] = (1.0 - a*pole->x*pole->x)*cos(pole->s) - a*pole->x*pole->y*sin(pole->s);
-    A[0][1] = -a*pole->x*pole->y*cos(pole->s) + (1.0 - a*pole->y*pole->y)*sin(pole->s);
-    A[0][2] = pole->x*cos(pole->s) + pole->y*sin(pole->s);
-    A[1][0] = -(1.0 - a*pole->x*pole->x)*sin(pole->s) - a*pole->x*pole->y*cos(pole->s);
+    A[1][0] = -a*pole->x*pole->y*cos(pole->s) + (1.0 - a*pole->y*pole->y)*sin(pole->s);
+    A[2][0] = pole->x*cos(pole->s) + pole->y*sin(pole->s);
+    A[0][1] = -(1.0 - a*pole->x*pole->x)*sin(pole->s) - a*pole->x*pole->y*cos(pole->s);
     A[1][1] = a*pole->x*pole->y*sin(pole->s) + (1.0 - a*pole->y*pole->y)*cos(pole->s);
-    A[1][2] = -pole->x*sin(pole->s) + pole->y*cos(pole->s);
-    A[2][0] = -pole->x;
-    A[2][1] = -pole->y;
+    A[2][1] = -pole->x*sin(pole->s) + pole->y*cos(pole->s);
+    A[0][2] = -pole->x;
+    A[1][2] = -pole->y;
     A[2][2] = 1.0 - a*(pole->x*pole->x + pole->y*pole->y);
-
+    double x, y, s;
+    x = pole->x;
+    y = pole->y;
+    s = pole->s;
+    /*    A[0][0] = (1.0 - a*x*x)*cos(s) + a*x*y*sin(s);
+        A[0][1] = (1.0 - a*x*x)*sin(s) - a*x*y*cos(s);
+        A[0][2] = x;
+        A[1][0] = -a*x*y*cos(s) - (1.0 - a*y*y)*sin(s);
+        A[1][1] = -a*x*y*sin(s) + (1.0 - a*y*y)*cos(s);
+        A[1][2] = y;
+        A[2][0] = -x*cos(s) + y*sin(s);
+        A[2][1] = -x*sin(s) - y*cos(s);
+        A[2][2] = 1.0 - a * (x*x + y*y);
+    */
+    /*    A[0][0] = (1.0 - a*x*x)*cos(s) + a*x*y*sin(s);
+        A[1][0] = (1.0 - a*x*x)*sin(s) - a*x*y*cos(s);
+        A[2][0] = x;
+        A[0][1] = -a*x*y*cos(s) - (1.0 - a*y*y)*sin(s);
+        A[1][1] = -a*x*y*sin(s) + (1.0 - a*y*y)*cos(s);
+        A[2][1] = y;
+        A[0][2] = -x*cos(s) + y*sin(s);
+        A[1][2] = -x*sin(s) - y*cos(s);
+        A[2][2] = 1.0 - a * (x*x + y*y);
+    */
     out = rotMatrix_To_Quat(A);
 

trunk/psLib/test/astro/tst_psEarthOrientation.c

@@ -5 +5 @@
 *  @author d-Rob, MHPCC
 *
-*  @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-12-07 02:13:40 $
+*  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-12-08 02:49:04 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -66 +66 @@
 {
     if (obj == NULL) {
-        obj = psSphereAlloc();
-        obj->r = objR;
-        obj->d = objD;
+        //        obj = psSphereAlloc();
+        //        obj->r = objR;
+        //        obj->d = objD;
+        psCube *tempCube = psCubeAlloc();
+        tempCube->x = -0.3598480726985338;
+        tempCube->y = 0.5555012823608123;
+        tempCube->z = 0.7496183628158023;
+        obj = psCubeToSphere(tempCube);
+        psFree(tempCube);
     }
 }
@@ -544 +550 @@
         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;
@@ -558 +560 @@
     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) {
+    if ( fabs(x-cube->x) > FLT_EPSILON  || fabs(y-cube->y) > FLT_EPSILON ||
+            fabs(z-cube->z) > FLT_EPSILON) {
         psError(PS_ERR_BAD_PARAMETER_VALUE, false,
                 "psSphereRot_TEOtoCEO returned incorrect values.\n");
@@ -570 +572 @@
     }
     //    psFree(now);
+    psFree(rot);
+    psFree(empty);
     psFree(result);
     psFree(obj);
@@ -615 +619 @@
     }
 
+    //    psCube *tempCube = psCubeAlloc();
+    //    tempCube->x = -0.35963388069046304;
+    //    tempCube->y = 0.5555192509816625;
+    //    tempCube->z = 0.749078321908413;
+    //    obj = psCubeToSphere(tempCube);
+    //    psFree(tempCube);
+    obj = psSphereAlloc();
+    obj->r = objR;
+    obj->d = objD;
+    psSphereRot *precessionNutation = psSphereRotConjugate(NULL, rot);
+    psSphere *result = psSphereRotApply(NULL, precessionNutation, obj);
+    psCube *cube = psSphereToCube(result);
+    double x, y, z;
+    x = -0.3598480726985338;
+    y = 0.5555012823608123;
+    z = 0.7496183628158023;
+    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));
+    psFree(obj);
+
+
     if (fabs(rot->q0-q0) > DBL_EPSILON || fabs(rot->q1-q1) > DBL_EPSILON ||
             fabs(rot->q2-q2) > DBL_EPSILON || fabs(rot->q3-q3) > DBL_EPSILON) {
@@ -635 +662 @@
     printf("\n  Sphere -test- has values r,d = %.8g, %.8g \n", test->r, test->d);
 
+    psFree(precessionNutation);
+    psFree(result);
+    psFree(cube);
+
     psFree(test);
     psFree(rot);