Index: trunk/psLib/src/astro/psEarthOrientation.c
===================================================================
--- trunk/psLib/src/astro/psEarthOrientation.c	(revision 5793)
+++ trunk/psLib/src/astro/psEarthOrientation.c	(revision 5814)
@@ -8,6 +8,6 @@
  *  @author Robert Daniel DeSonia, MHPCC
  *
- *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-15 02:37:48 $
+ *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-20 05:05:37 $
  *
  *  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
@@ -512,7 +512,5 @@
     S = SEC_TO_RAD(S * 1e-6);
 
-
     // now calculate the non-poly portion from the tables
-
     psF64* cols[17];
     for (int lcv = 0; lcv < 17; lcv++) {
@@ -628,12 +626,12 @@
                 out->x = cols[1][rowNum];
                 out->y = cols[2][rowNum];
-                out->x = SEC_TO_RAD(out->x) * 1e-6;
-                out->y = SEC_TO_RAD(out->y) * 1e-6;
+                out->x = SEC_TO_RAD(out->x) * 1e-3;
+                out->y = SEC_TO_RAD(out->y) * 1e-3;
                 rowNum = numRows;
             } else {
                 out->x = cols[3][rowNum];
                 out->y = cols[4][rowNum];
-                out->x = SEC_TO_RAD(out->x) * 1e-6;
-                out->y = SEC_TO_RAD(out->y) * 1e-6;
+                out->x = SEC_TO_RAD(out->x) * 1e-3;
+                out->y = SEC_TO_RAD(out->y) * 1e-3;
                 rowNum = numRows;
             }
@@ -762,5 +760,5 @@
     T += -2451545.0;
     double theta = 2.0 * M_PI * (0.7790572732640 + 1.00273781191135448 * T);
-    psSphereRot *out = psSphereRotAlloc(theta, 0.0, 0.0);
+    psSphereRot *out = psSphereRotAlloc(-theta, 0.0, 0.0);
     //    psSphereRot *out = psSphereRotInvert(theta, 0.0, 0.0);
 
@@ -773,43 +771,7 @@
 {
     PS_ASSERT_PTR_NON_NULL(time, NULL);
+    PS_ASSERT_INT_WITHIN_RANGE(bulletin, PS_IERS_A, PS_IERS_B, NULL);
 
     psEarthPole *out = psEarthPoleAlloc();
-    //    psSphere *in = NULL;
-    //    double x,y,s;
-    //    psTime *time2 = psTimeAlloc(time->type);
-    //    time2->sec = time->sec;
-    //    time2->nsec = time->nsec;
-    //    time2->leapsecond = time->leapsecond;
-    //XXX:  Time must be converted to TAI before use?
-    /*    if (time->type == PS_TIME_UT1) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "psTime cannot be of type UT1 here.\n");
-            return NULL;
-        } else if (time->type != PS_TIME_TAI) {
-            time2 = psTimeConvert(time2, PS_TIME_TAI);
-        }
-    */
-    //XXX: This may be the wrong idea... ADD says to use 3rd-order polys to interpolate
-    //polar motion coordinates.
-    //    in = p_psTimeGetPoleCoords(time2);
-    //    if (in == NULL) {
-    //        psError(PS_ERR_BAD_PARAMETER_NULL, true,
-    //                "p_psTimeGetPoleCoords return NULL psSphere for non-NULL input time.\n");
-    //        return NULL;
-    //    }
-    //    out->x = in->r;
-    //    out->y = in->d;
-    //    s = psTimeGetUT1Delta(time2, bulletin);
-    //    out->s = s;
-
-    //XXX: Apply polar tide correction here???
-    //    psEarthPole *correction = psEOC_PolarTideCorr(time2);
-    //    out->x += correction->x;
-    //    out->y += correction->y;
-    //    out->s += correction->s;
-
-    //    psFree(time2);
-    //    psFree(in);
-    //    psFree(correction);
-
     out->x = 0.0;
     out->y = 0.0;
@@ -817,12 +779,12 @@
 
     double MJD = psTimeToMJD(time);
-
-    if (MJD < 41684.0 || (MJD > 53434.0 && bulletin == PS_IERS_B) ||
-            (MJD > 53858.0 && bulletin == PS_IERS_A) ) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
-                "Invalid time input.  Date, %lf, is out of range\n", MJD);
-        return out;
-    }
-
+    /*
+        if (MJD < 41684.0 || (MJD > 53434.0 && bulletin == PS_IERS_B) ||
+                (MJD > 53858.0 && bulletin == PS_IERS_A) ) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                    "Invalid time input.  Date, %lf, is out of range\n", MJD);
+            return out;
+        }
+    */
     // Check if EOC data loaded
     if(! eocInitialized) {
@@ -841,26 +803,91 @@
     }
     int numRows = ((psVector*)(finalsTable->data[0]))->n;
-    for (int rowNum = 0; rowNum < numRows; rowNum++) {
-        if ( (MJD - cols[0][rowNum]) < 1.0 ) {
-            if (bulletin == PS_IERS_A) {
-                out->x = cols[1][rowNum];
-                out->y = cols[2][rowNum];
-                out->s = cols[3][rowNum];
-                out->x = SEC_TO_RAD(out->x);
-                out->y = SEC_TO_RAD(out->y);
-                out->s = SEC_TO_RAD(out->s);
-                rowNum = numRows;
-            } else {
-                out->x = cols[4][rowNum];
-                out->y = cols[5][rowNum];
-                out->s = cols[6][rowNum];
-                out->x = SEC_TO_RAD(out->x);
-                out->y = SEC_TO_RAD(out->y);
-                out->s = SEC_TO_RAD(out->s);
-                rowNum = numRows;
+    psVector *X = psVectorAlloc(numRows, PS_TYPE_F64);
+    psVector *Y = psVectorAlloc(numRows, PS_TYPE_F64);
+    psVector *S = psVectorAlloc(numRows, PS_TYPE_F64);
+    psVector *T = psVectorAlloc(numRows, PS_TYPE_F64);
+    if (bulletin == PS_IERS_A) {
+        for (int rowNum = 0; rowNum < numRows; rowNum++) {
+            T->data.F64[rowNum] = cols[0][rowNum];
+            X->data.F64[rowNum] = cols[1][rowNum];
+            Y->data.F64[rowNum] = cols[2][rowNum];
+            S->data.F64[rowNum] = cols[3][rowNum];
+        }
+    } else {
+        for (int rowNum = 0; rowNum < numRows; rowNum++) {
+            T->data.F64[rowNum] = cols[0][rowNum];
+            X->data.F64[rowNum] = cols[4][rowNum];
+            Y->data.F64[rowNum] = cols[5][rowNum];
+            S->data.F64[rowNum] = cols[6][rowNum];
+        }
+    }
+
+    double xOut = 0.0;
+    double yOut = 0.0;
+    double sOut = 0.0;
+    double xTerm = 0.0;
+    double yTerm = 0.0;
+    double sTerm = 0.0;
+    int k = 0;
+    for (int i = 0; i < (numRows-1); i++) {
+        if (MJD >= T->data.F64[i] && MJD < T->data.F64[i+1]) {
+            k = i;
+            if (k < 2) {
+                k = 2;
             }
-        }
-    }
-
+            if (k > (numRows-2)) {
+                k = numRows-2;
+            }
+            for (int m = k-1; m <= k+2; m++) {
+                xTerm = X->data.F64[m];
+                yTerm = Y->data.F64[m];
+                sTerm = S->data.F64[m];
+                for (int j = k-1; j <= k+2; j++) {
+                    if ( m != j) {
+                        double term = (MJD - T->data.F64[j])/(T->data.F64[m] - T->data.F64[j]);
+                        xTerm *= term;
+                        yTerm *= term;
+                        sTerm *= term;
+                    }
+                }
+                xOut += xTerm;
+                yOut += yTerm;
+                sOut += sTerm;
+            }
+            i = numRows-1;
+        }
+    }
+    out->x = SEC_TO_RAD(xOut);
+    out->y = SEC_TO_RAD(yOut);
+    out->s = SEC_TO_RAD(sOut);
+
+
+    /*
+        for (int rowNum = 0; rowNum < numRows; rowNum++) {
+            if ( (MJD - cols[0][rowNum]) < 1.0 ) {
+                if (bulletin == PS_IERS_A) {
+                    out->x = cols[1][rowNum];
+                    out->y = cols[2][rowNum];
+                    out->s = cols[3][rowNum];
+                    out->x = SEC_TO_RAD(out->x);
+                    out->y = SEC_TO_RAD(out->y);
+                    out->s = SEC_TO_RAD(out->s);
+                    rowNum = numRows;
+                } else {
+                    out->x = cols[4][rowNum];
+                    out->y = cols[5][rowNum];
+                    out->s = cols[6][rowNum];
+                    out->x = SEC_TO_RAD(out->x);
+                    out->y = SEC_TO_RAD(out->y);
+                    out->s = SEC_TO_RAD(out->s);
+                    rowNum = numRows;
+                }
+            }
+        }
+    */
+    psFree(X);
+    psFree(Y);
+    psFree(S);
+    psFree(T);
     return out;
 }
@@ -1081,4 +1108,5 @@
     y = motion->y;
     s = motion->s;
+    s = -s;
 
 
@@ -1130,31 +1158,31 @@
 
     /*
-            psSphereRot r,s,t;
+                psSphereRot r,p,t;
      
-            r.q0=sin(motion->y/2.0);
-            r.q1=0;
-            r.q2=0;
-            r.q3=cos(motion->y/2.0);
+                r.q0=sin(y/2.0);
+                r.q1=0;
+                r.q2=0;
+                r.q3=cos(y/2.0);
      
-            s.q0=0;
-            s.q1=sin(motion->x/2.0);
-            s.q2=0;
-            s.q3=cos(motion->x/2.0);
+                p.q0=0;
+                p.q1=sin(x/2.0);
+                p.q2=0;
+                p.q3=cos(x/2.0);
      
-            t.q0=0;
-            t.q1=0;
-            t.q2=sin(motion->s/2.0);
-            t.q3=cos(motion->s/2.0);
+                t.q0=0;
+                t.q1=0;
+                t.q2=sin(s/2.0);
+                t.q3=cos(s/2.0);
      
-            // calculate t*s*r.
-            psSphereRot* temp = psSphereRotCombine(NULL,&t,&s);
-            out = psSphereRotCombine(NULL, temp, &r);
-            psFree(temp);
+                // calculate t*s*r.
+                psSphereRot* temp = psSphereRotCombine(NULL,&t,&p);
+                out = psSphereRotCombine(NULL, temp, &r);
+                psFree(temp);
      
-            return out;
-    */
-
-    s = -s;
+                return out;
+    */
+
     //Newest trial - mult. y matrix * x * z
+
     A[0][0] = cos(x)*cos(s);
     A[1][0] = cos(x)*sin(s);
@@ -1167,13 +1195,14 @@
     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);
+    /*
+            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
@@ -1203,50 +1232,4 @@
     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;
-        double recip;
-        diag_sum[0] = 1.0 + A[0][0] - A[1][1] - A[2][2];
-        diag_sum[1] = 1.0 - A[0][0] + A[1][1] - A[2][2];
-        diag_sum[2] = 1.0 - A[0][0] - A[1][1] + A[2][2];
-        diag_sum[3] = 1.0 + A[0][0] + A[1][1] + A[2][2];
-
-        maxi = 0;
-        for (int i = 1; i < 4; ++i) {
-            if (diag_sum[i] > diag_sum[maxi]) {
-                maxi = i;
-            }
-        }
-
-        double p = 0.5 * sqrt(diag_sum[maxi]);
-        recip = 1.0 / (4.0 * p);
-
-        if (maxi == 0) {
-            out->q0 = p;
-            out->q1 = recip * (A[0][1] + A[1][0]);
-            out->q2 = recip * (A[2][0] + A[0][2]);
-            out->q3 = recip * (A[1][2] - A[2][1]);
-        } else if (maxi == 1) {
-            out->q0 = recip * (A[0][1] + A[1][0]);
-            out->q1 = p;
-            out->q2 = recip * (A[1][2] + A[2][1]);
-            out->q3 = recip * (A[2][0] - A[0][2]);
-        } else if (maxi == 2) {
-            out->q0 = recip * (A[2][0] + A[0][2]);
-            out->q1 = recip * (A[1][2] + A[2][1]);
-            out->q2 = p;
-            out->q3 = recip * (A[0][1] - A[1][0]);
-        } else if (maxi == 3) {
-            out->q0 = recip * (A[1][2] - A[2][1]);
-            out->q1 = recip * (A[2][0] - A[0][2]);
-            out->q2 = recip * (A[0][1] - A[1][0]);
-            out->q3 = p;
-        }
-    */
     return out;
 }