Changeset 5771 for trunk/psLib/src/astro/psEarthOrientation.c

Timestamp:

Dec 12, 2005, 3:31:54 PM (21 years ago)

Author:

drobbin

Message:

Updated tests and debugged code. PrecessionModel now working.

File:

• trunk/psLib/src/astro/psEarthOrientation.c (modified) (8 diffs)

trunk/psLib/src/astro/psEarthOrientation.c

@@ -8 +8 @@
  *  @author Robert Daniel DeSonia, MHPCC
  *
- *  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-08 02:49:04 $
+ *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-13 01:31:54 $
  *
  *  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
@@ -432 +432 @@
     double S = psPolynomial1DEval(sPoly,t);
 
+    X = SEC_TO_RAD(X * 1e-6);
+    Y = SEC_TO_RAD(Y * 1e-6);
+    S = SEC_TO_RAD(S * 1e-6);
+
+
     // now calculate the non-poly portion from the tables
 
@@ -448 +453 @@
         double as = cols[1][lcv];
         double ac = cols[2][lcv];
-
-        X += as*tj*sin(arg) + ac*tj*cos(arg);
+        //        printf("as-x, ac-x, = %.13g, %.13g\n", as, ac);
+        as = SEC_TO_RAD(as) * 1e-6;
+        ac = SEC_TO_RAD(ac) * 1e-6;
+        //        X += as*tj*sin(arg) + ac*tj*cos(arg);
+        X += (as*tj*sin(arg) + ac*cos(arg)) * tj;
     }
 
@@ -465 +473 @@
         double as = cols[1][lcv];
         double ac = cols[2][lcv];
-
-        Y += as*tj*sin(arg) + ac*tj*cos(arg);
+        as = SEC_TO_RAD(as) * 1e-6;
+        ac = SEC_TO_RAD(ac) * 1e-6;
+
+        //        Y += as*tj*sin(arg) + ac*tj*cos(arg);
+        Y += (as*tj*sin(arg) + ac*cos(arg)) * tj;
     }
 
@@ -482 +493 @@
         double as = cols[1][lcv];
         double ac = cols[2][lcv];
-
-        S += as*tj*sin(arg) + ac*tj*cos(arg);
+        as = SEC_TO_RAD(as) * 1e-6;
+        ac = SEC_TO_RAD(ac) * 1e-6;
+
+        //        S += as*tj*sin(arg) + ac*tj*cos(arg);
+        S += (as*tj*sin(arg) + ac*cos(arg)) * tj;
     }
 
@@ -700 +714 @@
         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[1][0] = -a*pole->x*pole->y*cos(pole->s) + (1.0 - a*pole->y*pole->y)*sin(pole->s);
@@ -709 +724 @@
     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;
@@ -1010 +1026 @@
     psSphereRot *out = NULL;
 
+    double x,y,s;
+    x = motion->x;
+    y = motion->y;
+    s = motion->s;
+
+
     //Setup Rotation Matrix for transformation (x,y,z rotation)
     //XXX: May need to be (z,y,x as in Mathworld?)
-    A[0][0] = cos(motion->x)*cos(-motion->s);
-    A[0][1] = sin(motion->y)*sin(motion->x)*cos(-motion->s) + cos(motion->y)*sin(-motion->s);
-    A[0][2] = -cos(motion->y)*sin(motion->x)*cos(-motion->s) + sin(motion->y)*sin(-motion->s);
-    A[1][0] = -cos(motion->x)*sin(-motion->s);
-    A[1][1] = -sin(motion->y)*sin(motion->x)*sin(-motion->s) + cos(motion->y)*cos(-motion->s);
-    A[1][2] = cos(motion->y)*sin(motion->x)*sin(-motion->s) + sin(motion->y)*cos(-motion->s);
-    A[2][0] = sin(motion->x);
-    A[2][1] = -sin(motion->y)*cos(motion->x);
-    A[2][2] = cos(motion->y)*cos(motion->x);
+    /*   A[0][0] = cos(motion->x)*cos(-motion->s);
+        A[0][1] = sin(motion->y)*sin(motion->x)*cos(-motion->s) + cos(motion->y)*sin(-motion->s);
+        A[0][2] = -cos(motion->y)*sin(motion->x)*cos(-motion->s) + sin(motion->y)*sin(-motion->s);
+        A[1][0] = -cos(motion->x)*sin(-motion->s);
+        A[1][1] = -sin(motion->y)*sin(motion->x)*sin(-motion->s) + cos(motion->y)*cos(-motion->s);
+        A[1][2] = cos(motion->y)*sin(motion->x)*sin(-motion->s) + sin(motion->y)*cos(-motion->s);
+        A[2][0] = sin(motion->x);
+        A[2][1] = -sin(motion->y)*cos(motion->x);
+        A[2][2] = cos(motion->y)*cos(motion->x);
+    */
+    /*
+        A[0][0] = cos(motion->x)*cos(-motion->s);
+        A[1][0] = sin(motion->y)*sin(motion->x)*cos(-motion->s) + cos(motion->y)*sin(-motion->s);
+        A[2][0] = -cos(motion->y)*sin(motion->x)*cos(-motion->s) + sin(motion->y)*sin(-motion->s);
+        A[0][1] = -cos(motion->x)*sin(-motion->s);
+        A[1][1] = -sin(motion->y)*sin(motion->x)*sin(-motion->s) + cos(motion->y)*cos(-motion->s);
+        A[2][1] = cos(motion->y)*sin(motion->x)*sin(-motion->s) + sin(motion->y)*cos(-motion->s);
+        A[0][2] = sin(motion->x);
+        A[1][2] = -sin(motion->y)*cos(motion->x);
+        A[2][2] = cos(motion->y)*cos(motion->x);
+    */
+    /*
+        A[0][0] = cos(s)*cos(x);
+        A[0][1] = sin(s)*sin(y) + cos(s)*sin(x)*sin(y);
+        A[0][2] = sin(s)*sin(y) - cos(s)*sin(x)*cos(y);
+        A[1][0] = -sin(s)*cos(x);
+        A[1][1] = cos(s)*cos(y) - sin(s)*sin(x)*sin(y);
+        A[1][2] = cos(s)*sin(y) + sin(s)*sin(x)*cos(y);
+        A[2][0] = sin(x);
+        A[2][1] = -cos(x)*sin(y);
+        A[2][2] = cos(x)*cos(y);
+    */
+    /*
+        A[0][0] = cos(s)*cos(x);
+        A[1][0] = sin(s)*sin(y) + cos(s)*sin(x)*sin(y);
+        A[2][0] = sin(s)*sin(y) - cos(s)*sin(x)*cos(y);
+        A[0][1] = -sin(s)*cos(x);
+        A[1][1] = cos(s)*cos(y) - sin(s)*sin(x)*sin(y);
+        A[2][1] = cos(s)*sin(y) + sin(s)*sin(x)*cos(y);
+        A[0][2] = sin(x);
+        A[1][2] = -cos(x)*sin(y);
+        A[2][2] = cos(x)*cos(y);
+    */
+
+    /*
+        psSphereRot r,s,t;
+     
+        // directly from ADD -- there must be a better way?!
+        r.q0=sin(motion->y/2.0);
+        r.q1=0;
+        r.q2=0;
+        r.q3=cos(motion->y/2.0);
+     
+        s.q0=0;
+        s.q1=sin(motion->x/2.0);
+        s.q2=0;
+        s.q3=cos(motion->x/2.0);
+     
+        t.q0=0;
+        t.q1=0;
+        t.q2=sin(-motion->s/2.0);
+        t.q3=cos(-motion->s/2.0);
+     
+        // calculate t*s*r.
+        psSphereRot* temp = psSphereRotCombine(NULL,&t,&s);
+        out = psSphereRotCombine(NULL, temp, &r);
+        psFree(temp);
+     
+        return out;
+    */
+
+    s = -s;
+    //Newest trial - mult. y matrix * x * z
+    /*    A[0][0] = cos(x)*cos(s);
+        A[1][0] = cos(x)*sin(s);
+        A[2][0] = -sin(x);
+        A[0][1] = sin(x)*sin(y)*cos(s) - cos(y)*sin(s);
+        A[1][1] = sin(x)*sin(y)*sin(s) + cos(y)*cos(s);
+        A[2][1] = cos(x)*sin(y);
+        A[0][2] = sin(x)*cos(y)*cos(s) + sin(y)*sin(s);
+        A[1][2] = sin(x)*cos(y)*sin(s) - sin(y)*cos(s);
+        A[2][2] = cos(x)*cos(y);
+    */
+    /*    A[0][0] = cos(x)*cos(s);
+        A[0][1] = cos(x)*sin(s);
+        A[0][2] = -sin(x);
+        A[1][0] = sin(x)*sin(y)*cos(s) - cos(y)*sin(s);
+        A[1][1] = sin(x)*sin(y)*sin(s) + cos(y)*cos(s);
+        A[1][2] = cos(x)*sin(y);
+        A[2][0] = sin(x)*cos(y)*cos(s) + sin(y)*sin(s);
+        A[2][1] = sin(x)*cos(y)*sin(s) - sin(y)*cos(s);
+        A[2][2] = cos(x)*cos(y);
+    */
+    //New trial - mult z * y * x
+    /*    A[0][0] = cos(s)*cos(x);
+        A[1][0] = sin(s)*cos(y) + cos(s)*sin(x)*sin(y);
+        A[2][0] = sin(s)*sin(y) - cos(s)*sin(x)*cos(y);
+        A[0][1] = -sin(s)*cos(x);
+        A[1][1] = cos(s)*cos(y) - sin(s)*sin(x)*sin(y);
+        A[2][1] = cos(s)*sin(y) + sin(s)*sin(x)*cos(y);
+        A[0][2] = sin(x);
+        A[1][2] = -cos(x)*sin(y);
+        A[2][2] = cos(x)*cos(y);
+    */
+    A[0][0] = cos(s)*cos(x);
+    A[0][1] = sin(s)*cos(y) + cos(s)*sin(x)*sin(y);
+    A[0][2] = sin(s)*sin(y) - cos(s)*sin(x)*cos(y);
+    A[1][0] = -sin(s)*cos(x);
+    A[1][1] = cos(s)*cos(y) - sin(s)*sin(x)*sin(y);
+    A[1][2] = cos(s)*sin(y) + sin(s)*sin(x)*cos(y);
+    A[2][0] = sin(x);
+    A[2][1] = -cos(x)*sin(y);
+    A[2][2] = cos(x)*cos(y);
+
 
     //Convert rotation matrix to quaternions
     out = rotMatrix_To_Quat(A);
+
+    //    out = psSphereRotAlloc(motion->y, motion->x, motion->s);
+
+    //    out->q0 = SEC_TO_RAD(out->q0);
+    //    out->q1 = SEC_TO_RAD(out->q1);
+    //    out->q2 = SEC_TO_RAD(out->q2) + (out->q0+out->q1/2.0);
+    //    out->q3 = sqrt(1.0 - (out->q0*out->q0 + out->q1*out->q1 + out->q2*out->q2) );
     /*    double diag_sum[3];
         int maxi;