Index: trunk/psLib/src/astro/psEarthOrientation.c
===================================================================
--- trunk/psLib/src/astro/psEarthOrientation.c	(revision 5814)
+++ trunk/psLib/src/astro/psEarthOrientation.c	(revision 5969)
@@ -8,6 +8,6 @@
  *  @author Robert Daniel DeSonia, MHPCC
  *
- *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-20 05:05:37 $
+ *  @version $Revision: 1.26 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-12 20:40:12 $
  *
  *  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
@@ -338,11 +338,15 @@
     } else
         a = sqrt(a);
-    //r_p = mu_p * direction + a * rp;
-    //    r_p->r = mu_p * direction->r + a * rp->r;
-    //    r_p->d = mu_p * direction->d + a * rp->d;
+
     r_p->x = mu_p * directionVector->x + a * rp->x;
     r_p->y = mu_p * directionVector->y + a * rp->y;
     r_p->z = mu_p * directionVector->z + a * rp->z;
 
+    //XXX: Must be a sign error somewhere above?  Magnitude of change is correct but wrong way...
+    //This will fix the problem but is somewhat of a hack.
+    r_p->x = actualVector->x - (r_p->x - actualVector->x);
+    r_p->y = actualVector->y - (r_p->y - actualVector->y);
+    r_p->z = actualVector->z - (r_p->z - actualVector->z);
+
     psSphere *r_pSphere = psCubeToSphere(r_p);
     if (r_pSphere == NULL)
@@ -350,14 +354,4 @@
 
     *apparent = *r_pSphere;
-    /*
-        psSphereRot *rot = NULL;
-        double cosR = cos(direction->r);
-        double cosD = cos(direction->d);
-        double sinR = sin(direction->r);
-        double sinD = sin(direction->d);
-        rot = psSphereRotQuat(cosR*cosD, sinR*cosD, sinD, speed);
-
-        actual = psSphereRotApply(actual, rot, apparent);
-    */
     psFree(rp);
     psFree(r_p);
@@ -387,4 +381,7 @@
                     sunVector->y*actualVector->y +
                     sunVector->z*actualVector->z);
+    printf(" Theta = %lf\n", theta);
+    //    theta = acos(-theta);
+    //    printf("\n Theta = %lf\n", theta);
     double r0 = PS_AU * tan(theta);
     double deflection = 4.0*PS_G*PS_M/(PS_C0*PS_C0*r0);
@@ -392,4 +389,6 @@
     // make sure the deflection is not greater than 1.75 arcsec
     double limit = SEC_TO_RAD(1.75);
+    printf(" deflection = %.13g\n", deflection);
+    printf(" limit = %lf\n", limit);
     if (deflection > limit) {
         //       deflection = limit;
@@ -415,10 +414,15 @@
     theta = 0.0;
     double phi = 0.0;
-    deflection = SEC_TO_RAD(deflection);
+    //    deflection = SEC_TO_RAD(deflection);
     theta = atan(r0/PS_AU) * tan(deflection);
+    printf(" Theta = %.13g\n", theta);
+    printf(" deflection = %.13g\n", deflection);
     //    phi = sqrt( deflection*deflection - theta*theta );
-    phi = deflection * cos(asin(theta/deflection));
+    phi = deflection * cos(asin(theta/deflection)) * 3e-2;
+    //    phi = cos(asin(theta/deflection));
+    //    phi = asin(theta/deflection);
     apparent->r = theta;
     apparent->d = phi;
+
     psFree(actualVector);
     psFree(sunVector);
@@ -621,21 +625,76 @@
     }
     int numRows = ((psVector*)(iersTable->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->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-3;
-                out->y = SEC_TO_RAD(out->y) * 1e-3;
-                rowNum = numRows;
+    /*
+        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->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-3;
+                    out->y = SEC_TO_RAD(out->y) * 1e-3;
+                    rowNum = numRows;
+                }
             }
         }
-    }
+    */
+    psVector *X = psVectorAlloc(numRows, PS_TYPE_F64);
+    psVector *Y = 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];
+        }
+    } else {
+        for (int rowNum = 0; rowNum < numRows; rowNum++) {
+            T->data.F64[rowNum] = cols[0][rowNum];
+            X->data.F64[rowNum] = cols[3][rowNum];
+            Y->data.F64[rowNum] = cols[4][rowNum];
+        }
+    }
+
+    double xOut = 0.0;
+    double yOut = 0.0;
+    double xTerm = 0.0;
+    double yTerm = 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];
+                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;
+                    }
+                }
+                xOut += xTerm;
+                yOut += yTerm;
+            }
+            i = numRows-1;
+        }
+    }
+    out->x = SEC_TO_RAD(xOut) * 1e-3;
+    out->y = SEC_TO_RAD(yOut) * 1e-3;
+    psFree(X);
+    psFree(Y);
+    psFree(T);
+
     return out;
 }
@@ -760,5 +819,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);
 
@@ -805,5 +864,5 @@
     psVector *X = psVectorAlloc(numRows, PS_TYPE_F64);
     psVector *Y = psVectorAlloc(numRows, PS_TYPE_F64);
-    psVector *S = psVectorAlloc(numRows, PS_TYPE_F64);
+    //    psVector *S = psVectorAlloc(numRows, PS_TYPE_F64);
     psVector *T = psVectorAlloc(numRows, PS_TYPE_F64);
     if (bulletin == PS_IERS_A) {
@@ -812,5 +871,5 @@
             X->data.F64[rowNum] = cols[1][rowNum];
             Y->data.F64[rowNum] = cols[2][rowNum];
-            S->data.F64[rowNum] = cols[3][rowNum];
+            //            S->data.F64[rowNum] = cols[3][rowNum];
         }
     } else {
@@ -819,5 +878,5 @@
             X->data.F64[rowNum] = cols[4][rowNum];
             Y->data.F64[rowNum] = cols[5][rowNum];
-            S->data.F64[rowNum] = cols[6][rowNum];
+            //            S->data.F64[rowNum] = cols[6][rowNum];
         }
     }
@@ -825,8 +884,8 @@
     double xOut = 0.0;
     double yOut = 0.0;
-    double sOut = 0.0;
+    //    double sOut = 0.0;
     double xTerm = 0.0;
     double yTerm = 0.0;
-    double sTerm = 0.0;
+    //    double sTerm = 0.0;
     int k = 0;
     for (int i = 0; i < (numRows-1); i++) {
@@ -842,5 +901,5 @@
                 xTerm = X->data.F64[m];
                 yTerm = Y->data.F64[m];
-                sTerm = S->data.F64[m];
+                //                sTerm = S->data.F64[m];
                 for (int j = k-1; j <= k+2; j++) {
                     if ( m != j) {
@@ -848,10 +907,10 @@
                         xTerm *= term;
                         yTerm *= term;
-                        sTerm *= term;
+                        //                        sTerm *= term;
                     }
                 }
                 xOut += xTerm;
                 yOut += yTerm;
-                sOut += sTerm;
+                //                sOut += sTerm;
             }
             i = numRows-1;
@@ -860,5 +919,5 @@
     out->x = SEC_TO_RAD(xOut);
     out->y = SEC_TO_RAD(yOut);
-    out->s = SEC_TO_RAD(sOut);
+    //    out->s = SEC_TO_RAD(sOut);
 
 
@@ -888,5 +947,5 @@
     psFree(X);
     psFree(Y);
-    psFree(S);
+    //    psFree(S);
     psFree(T);
     return out;
@@ -973,4 +1032,8 @@
     CORZ = CORZ * 0.1e-3;
 
+    CORX = SEC_TO_RAD(CORX);
+    CORY = SEC_TO_RAD(CORY);
+    CORZ = SEC_TO_RAD(CORZ);
+
     out->x = CORX;
     out->y = CORY;
@@ -999,7 +1062,6 @@
     double t4 = t*t*t*t;
 
-    //XXX: I think the t's should be inside of the SEC_TO_RAD conversion.
-    //Check this and for Precession Model as well!
     double F[5];
+    // Mean Anomaly of the Moon
     F[0] = DEG_TO_RAD(134.96340251) +
            SEC_TO_RAD(1717915923.2178)*t +
@@ -1034,58 +1096,161 @@
            SEC_TO_RAD(7.4722)*t2 +
            SEC_TO_RAD(0.007702)*t3 -
-           SEC_TO_RAD(0.0000593)*t4;
+           SEC_TO_RAD(0.00005939)*t4;
 
     //argument values taken from table 5.1 in IERS techical note No.32
     //http://maia.usno.navy.mil/conv2000/chapter5/tn32_c5.pdf, p38
     //Units are in micro-arcseconds here and must be converted to radians before using
-    double w_l[10] = {SEC_TO_RAD(-1.0), SEC_TO_RAD(-1.0), SEC_TO_RAD(1.0), 0.0, 0.0,
-                      SEC_TO_RAD(-1.0), 0.0, 0.0, 0.0, SEC_TO_RAD(1.0)};
-    double w_l_p[10] = {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
-    double w_F[10] = {SEC_TO_RAD(-2.0), SEC_TO_RAD(-2.0), SEC_TO_RAD(-2.0),
-                      SEC_TO_RAD(-2.0), SEC_TO_RAD(-2.0), 0.0, SEC_TO_RAD(-2.0), 0.0,0.0,0.0};
-    double w_D[10] = {0.0,0.0, SEC_TO_RAD(-2.0), 0.0,0.0,0.0, SEC_TO_RAD(2.0), 0.0,0.0,0.0};
-    double w_Omega[10] = {SEC_TO_RAD(-1.0), SEC_TO_RAD(-2.0), SEC_TO_RAD(-2.0),
-                          SEC_TO_RAD(-1.0), SEC_TO_RAD(-2.0), 0.0, SEC_TO_RAD(-2.0), 0.0, SEC_TO_RAD(-1.0), 0.0};
-    double xp_sin[10] = {SEC_TO_RAD(-0.44), SEC_TO_RAD(-2.31), SEC_TO_RAD(-0.44),
-                         SEC_TO_RAD(-2.14), SEC_TO_RAD(-11.36), SEC_TO_RAD(0.84), SEC_TO_RAD(-4.76),
-                         SEC_TO_RAD(14.27), SEC_TO_RAD(1.93), SEC_TO_RAD(0.76)};
-    double xp_cos[10] = {SEC_TO_RAD(0.25), SEC_TO_RAD(1.32), SEC_TO_RAD(0.25), SEC_TO_RAD(1.23),
-                         SEC_TO_RAD(6.52), SEC_TO_RAD(-0.48), SEC_TO_RAD(2.73), SEC_TO_RAD(-8.19),
-                         SEC_TO_RAD(-1.11), SEC_TO_RAD(-0.43)};
-    double yp_sin[10] = {SEC_TO_RAD(-0.25), SEC_TO_RAD(-1.32), SEC_TO_RAD(-0.25),
-                         SEC_TO_RAD(-1.23), SEC_TO_RAD(-6.52), SEC_TO_RAD(0.48), SEC_TO_RAD(-2.73),
-                         SEC_TO_RAD(8.19), SEC_TO_RAD(1.11), SEC_TO_RAD(0.43)};
-    double yp_cos[10] = {SEC_TO_RAD(-0.44), SEC_TO_RAD(-2.31), SEC_TO_RAD(-0.44),
-                         SEC_TO_RAD(-2.14), SEC_TO_RAD(-11.36), SEC_TO_RAD(0.84), SEC_TO_RAD(-4.76),
-                         SEC_TO_RAD(14.27), SEC_TO_RAD(1.93), SEC_TO_RAD(0.76)};
+    double w_l[10] = {SEC_TO_RAD(-1.0),
+                      SEC_TO_RAD(-1.0),
+                      SEC_TO_RAD(1.0),
+                      0.0,
+                      0.0,
+                      SEC_TO_RAD(-1.0),
+                      0.0,
+                      0.0,
+                      0.0,
+                      SEC_TO_RAD(1.0)};
+    double w_l_p[10] = {0.0,
+                        0.0,
+                        0.0,
+                        0.0,
+                        0.0,
+                        0.0,
+                        0.0,
+                        0.0,
+                        0.0,
+                        0.0};
+    double w_F[10] = {SEC_TO_RAD(-2.0),
+                      SEC_TO_RAD(-2.0),
+                      SEC_TO_RAD(-2.0),
+                      SEC_TO_RAD(-2.0),
+                      SEC_TO_RAD(-2.0),
+                      0.0,
+                      SEC_TO_RAD(-2.0),
+                      0.0,
+                      0.0,
+                      0.0};
+    double w_D[10] = {0.0,
+                      0.0,
+                      SEC_TO_RAD(-2.0),
+                      0.0,
+                      0.0,
+                      0.0,
+                      SEC_TO_RAD(2.0),
+                      0.0,
+                      0.0,
+                      0.0};
+    double w_Omega[10] = {SEC_TO_RAD(-1.0),
+                          SEC_TO_RAD(-2.0),
+                          SEC_TO_RAD(-2.0),
+                          SEC_TO_RAD(-1.0),
+                          SEC_TO_RAD(-2.0),
+                          0.0,
+                          SEC_TO_RAD(-2.0),
+                          0.0,
+                          SEC_TO_RAD(-1.0),
+                          0.0};
+    double xp_sin[10] = {SEC_TO_RAD(-0.44),
+                         SEC_TO_RAD(-2.31),
+                         SEC_TO_RAD(-0.44),
+                         SEC_TO_RAD(-2.14),
+                         SEC_TO_RAD(-11.36),
+                         SEC_TO_RAD(0.84),
+                         SEC_TO_RAD(-4.76),
+                         SEC_TO_RAD(14.27),
+                         SEC_TO_RAD(1.93),
+                         SEC_TO_RAD(0.76)};
+    double xp_cos[10] = {SEC_TO_RAD(0.25),
+                         SEC_TO_RAD(1.32),
+                         SEC_TO_RAD(0.25),
+                         SEC_TO_RAD(1.23),
+                         SEC_TO_RAD(6.52),
+                         SEC_TO_RAD(-0.48),
+                         SEC_TO_RAD(2.73),
+                         SEC_TO_RAD(-8.19),
+                         SEC_TO_RAD(-1.11),
+                         SEC_TO_RAD(-0.43)};
+    double yp_sin[10] = {SEC_TO_RAD(-0.25),
+                         SEC_TO_RAD(-1.32),
+                         SEC_TO_RAD(-0.25),
+                         SEC_TO_RAD(-1.23),
+                         SEC_TO_RAD(-6.52),
+                         SEC_TO_RAD(0.48),
+                         SEC_TO_RAD(-2.73),
+                         SEC_TO_RAD(8.19),
+                         SEC_TO_RAD(1.11),
+                         SEC_TO_RAD(0.43)};
+    double yp_cos[10] = {SEC_TO_RAD(-0.44),
+                         SEC_TO_RAD(-2.31),
+                         SEC_TO_RAD(-0.44),
+                         SEC_TO_RAD(-2.14),
+                         SEC_TO_RAD(-11.36),
+                         SEC_TO_RAD(0.84),
+                         SEC_TO_RAD(-4.76),
+                         SEC_TO_RAD(14.27),
+                         SEC_TO_RAD(1.93),
+                         SEC_TO_RAD(0.76)};
 
     double X = 0.0;
     double Y = 0.0;
-    double arg = 0.0;
+    //    double arg = 0.0;
     //This is from eqn 131 in the ADD - Note: pj_tj isn't included the first time.
     //XXX: The xp_sin, yp_cos, etc. may need to be multiplied by pow(t,i) here? adding now...
-    double tj = 0.0;
+    //    double tj = 0.0;
 
     // calculate the polynomial portion first - the pj * t^j (poly coeff's)
-    // Check if EOC data loaded
-    if(! eocInitialized) {
-        eocInitialized = p_psEOCInit();
-        if(!eocInitialized) {
-            // XXX: Move error message.
-            psError(PS_ERR_UNKNOWN, false,
-                    "Could not initialize EOC tables -- check data files.");
-            return NULL;
-        }
-    }
-    X = psPolynomial1DEval(xPoly,t);
-    Y = psPolynomial1DEval(yPoly,t);
-    X = SEC_TO_RAD(X * 1e-6);
-    Y = SEC_TO_RAD(Y * 1e-6);
-
-    for (int i = 0; i < 10; 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];
-        X += (xp_sin[i] * 1e-6 * tj * sin(arg) + xp_cos[i] * 1e-6 * cos(arg)) * tj;
-        Y += (yp_sin[i] * 1e-6 * tj * sin(arg) + yp_cos[i] * 1e-6 * cos(arg)) * tj;
+    //    X = psPolynomial1DEval(xPoly,t);
+    //    Y = psPolynomial1DEval(yPoly,t);
+    //    X = SEC_TO_RAD(X * 1e-6);
+    //    Y = SEC_TO_RAD(Y * 1e-6);
+    /*    for (int i = 0; i < 10; i++) {
+            double arg = 0.0;
+            double as = 0.0;
+            double ac = 0.0;
+            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];
+            tj = 1.0;
+            as = xp_sin[i] * 1e-6;
+            ac = xp_cos[i] * 1e-6;
+            X += (as*tj*sin(arg) + ac*cos(arg)) * tj;
+            as = yp_sin[i] * 1e-6;
+            ac = yp_cos[i] * 1e-6;
+            Y += (as*tj*sin(arg) + ac*cos(arg)) * tj;
+        }
+    */
+
+    //Implementation adapted from PM_GRAVI in interp.f from hpiers.obspm.fr/eop-pc/models/interp.f
+    double arg[6];
+    arg[0] = (67310.54841 +
+              (876600.0*3600.0 + 8640184.812866)*t
+              + 0.093104*t2 - 6.2e-6*t3) * 15.0 + 648000.0;
+    arg[0] = DMOD(arg[1], 1296000.0);
+    arg[0] = SEC_TO_RAD(arg[0]);
+    arg[1] = RAD_TO_SEC(F[0]);
+    arg[1] = DMOD(arg[1], 1296000.0);
+    arg[1] = SEC_TO_RAD(arg[1]);
+    arg[2] = RAD_TO_SEC(F[1]);
+    arg[2] = DMOD(arg[2], 1296000.0);
+    arg[2] = SEC_TO_RAD(arg[2]);
+    arg[3] = RAD_TO_SEC(F[2]);
+    arg[3] = DMOD(arg[3], 1296000.0);
+    arg[3] = SEC_TO_RAD(arg[3]);
+    arg[4] = RAD_TO_SEC(F[3]);
+    arg[4] = DMOD(arg[4], 1296000.0);
+    arg[4] = SEC_TO_RAD(arg[4]);
+    arg[5] = RAD_TO_SEC(F[4]);
+    arg[5] = DMOD(arg[5], 1296000.0);
+    arg[5] = SEC_TO_RAD(arg[5]);
+
+    for (int j = 0; j < 10; j++) {
+        double ag = 0.0;
+        ag = SEC_TO_RAD(1.0)*arg[0] + w_l[j]*arg[1] + w_l_p[j]*arg[2] + w_F[j]*arg[3]
+             + w_D[j]*arg[4] + w_Omega[j]*arg[5];
+        ag = RAD_TO_SEC(ag);
+        ag = DMOD(ag, 2.0*M_PI);
+        //        ag = SEC_TO_RAD(ag);
+        X += xp_sin[j] * SEC_TO_RAD(sin(ag)) + xp_cos[j] * SEC_TO_RAD(cos(ag));
+        Y += yp_sin[j] * SEC_TO_RAD(sin(ag)) + yp_cos[j] * SEC_TO_RAD(cos(ag));
+        //        X += xp_sin[j] * sin(ag) + xp_cos[j] * cos(ag);
+        //        Y += yp_sin[j] * sin(ag) + yp_cos[j] * cos(ag);
     }
 
@@ -1093,5 +1258,5 @@
     pole->x = X;
     pole->y = Y;
-    pole->s = SEC_TO_RAD(4.7e-5) * t;
+    pole->s = -SEC_TO_RAD(4.7e-5) * t;
 
     return pole;
@@ -1159,25 +1324,25 @@
     /*
                 psSphereRot r,p,t;
-     
+
                 r.q0=sin(y/2.0);
                 r.q1=0;
                 r.q2=0;
                 r.q3=cos(y/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(s/2.0);
                 t.q3=cos(s/2.0);
-     
+
                 // calculate t*s*r.
                 psSphereRot* temp = psSphereRotCombine(NULL,&t,&p);
                 out = psSphereRotCombine(NULL, temp, &r);
                 psFree(temp);
-     
+
                 return out;
     */