Index: trunk/psLib/src/astro/psEarthOrientation.c =================================================================== --- trunk/psLib/src/astro/psEarthOrientation.c (revision 5749) +++ trunk/psLib/src/astro/psEarthOrientation.c (revision 5771) @@ -8,6 +8,6 @@ * @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,4 +432,9 @@ 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,6 +453,9 @@ 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,6 +473,9 @@ 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,6 +493,9 @@ 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,4 +714,5 @@ 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,4 +724,5 @@ 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,18 +1026,136 @@ 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;