Changeset 6425

Timestamp:

Feb 13, 2006, 2:09:27 PM (20 years ago)

Author:

drobbin

Message:

Added psMemThreadSafety. Removed p_lock from psList. Added comments and made minor changes to EOC.

Location:

trunk/psLib/src

Files:

• astro/psEarthOrientation.c (modified) (49 diffs)
• astro/psEarthOrientation.h (modified) (3 diffs)
• astro/psSphereOps.h (modified) (2 diffs)
• sys/psMemory.c (modified) (2 diffs)
• types/psList.c (modified) (11 diffs)
• types/psList.h (modified) (3 diffs)

trunk/psLib/src/astro/psEarthOrientation.c

@@ -8 +8 @@
  *  @author Robert Daniel DeSonia, MHPCC
  *
- *  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-02-06 21:57:04 $
+ *  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-02-14 00:09:26 $
  *
  *  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
@@ -356 +356 @@
                               psSphere *sun)
 {
+    //Check for valid inputs
     PS_ASSERT_PTR_NON_NULL(actual, NULL);
     PS_ASSERT_PTR_NON_NULL(sun, NULL);
 
     psSphere *temp = psSphereAlloc();
-    // calculating the apparent angle from the actual angle and the sun position
-
     psCube* sunVector = psSphereToCube(sun);
     psCube* actualVector = psSphereToCube(actual);
 
     // use dot product to calculate the angle of separation
-    // N.B., assuming the psSphereToCube function returns a unit vector.
-    //    double theta = acos(sunVector->x*actualVector->x +
     double dotProd = (sunVector->x*actualVector->x +
                       sunVector->y*actualVector->y + sunVector->z*actualVector->z);
@@ -378 +375 @@
     theta = acos(dotProd);
 
-    //    theta = acos( cos(sun->d)*cos(actual->d)*cos(sun->r-actual->r) + sin(sun->r)*sin(actual->r) );
-
     printf(" Theta = %.13g\n", theta);
-    //    theta = acos(-theta);
-    //    printf("\n Theta = %lf\n", theta);
     double r0 = PS_AU * tan(theta);
     printf(" r0 = %.19e\n", r0);
@@ -389 +382 @@
     // 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;
         //if deflection is greater than limit, the light rays will hit the sun
         psWarning("Invalid positions.  Light ray will not be seen on earth.\n");
@@ -400 +390 @@
     }
 
-    /*    if (apparent == NULL) {
-            apparent = psSphereAlloc();
-        } else {
-            apparent->r = 0.0;
-            apparent->d = 0.0;
-            apparent->rErr = 0.0;
-            apparent->dErr = 0.0;
-        }
-    */
+
     if (apparent != NULL) {
         psFree(apparent);
     }
     // bend the actual vector away from the sun vector by deflection angle.
-    // XXX: Not sure how to do this.  Dave thinks the formula should be:
-    //      theta = atan(r0/d)*tan(deflection), phi = thete/tan(deflection)
     theta = 0.0;
     double phi = 0.0;
@@ -477 +457 @@
     double t4 = t*t*t*t;
 
-    // N.B., the following formulae are from the ADD
-    // Mean Anomaly of the Moon
+    // The following formulae are from the ADD (& s5.8 of IERS Technical Note 32):
     double F[14];
+    // Mean Anomaly of the Moon (l)
     F[0] = DEG_TO_RAD(134.96340251) +
            SEC_TO_RAD(1717915923.2178)*t +
@@ -485 +465 @@
            SEC_TO_RAD(0.051635)*t3 -
            SEC_TO_RAD(0.00024470)*t4;
-
-    // Mean Anomaly of the Sun
+    // Mean Anomaly of the Sun (l')
     F[1] = DEG_TO_RAD(357.52910918) +
            SEC_TO_RAD(129596581.0481)*t -
@@ -492 +471 @@
            SEC_TO_RAD(0.000136)*t3 -
            SEC_TO_RAD(0.00001149)*t4;
-
-    // L - Omega
+    // L - Omega    (L = Mean Longitude of the Moon, Omega = F4)
     F[2] = DEG_TO_RAD(93.27209062) +
            SEC_TO_RAD(1739527262.8478)*t -
@@ -499 +477 @@
            SEC_TO_RAD(0.001037)*t3 +
            SEC_TO_RAD(0.00000417)*t4;
-
-    // Mean Elongation of the Moon from the Sun
+    // Mean Elongation of the Moon from the Sun (D)
     F[3] = DEG_TO_RAD(297.85019547) +
            SEC_TO_RAD(1602961601.2090)*t -
@@ -506 +483 @@
            SEC_TO_RAD(0.006593)*t3 -
            SEC_TO_RAD(0.00003169)*t4;
-
-    // Mean Longitude of the Ascending Node of the Moon
+    // Mean Longitude of the Ascending Node of the Moon (Omega)
     F[4] = DEG_TO_RAD(125.04455501) -
            SEC_TO_RAD(6962890.5431)*t +
@@ -513 +489 @@
            SEC_TO_RAD(0.007702)*t3 -
            SEC_TO_RAD(0.0000593)*t4;
-
+    //F5-F13 are the mean longitudes of the planets
     F[5] = 4.402608842 + 2608.7903141574*t;
     F[6] = 3.176146697 + 1021.3285546211*t;
@@ -528 +504 @@
         eocInitialized = p_psEOCInit();
         if(!eocInitialized) {
-            // XXX: Move error message.
             psError(PS_ERR_UNKNOWN, false,
                     "Could not initialize EOC tables -- check data files.");
@@ -534 +509 @@
         }
     }
-
     // calculate the polynomial portion first
     double X = psPolynomial1DEval(xPoly,t);
     double Y = psPolynomial1DEval(yPoly,t);
     double S = psPolynomial1DEval(sPoly,t);
-
+    //Units from the table & coefficients are in micro-arcseconds so convert to radians.
     X = SEC_TO_RAD(X * 1e-6);
     Y = SEC_TO_RAD(Y * 1e-6);
@@ -545 +519 @@
 
     // now calculate the non-poly portion from the tables
-    psF64* cols[17];
+    psF64* cols[17];  //Used to store all of the table information from tab5.2?.dat
     for (int lcv = 0; lcv < 17; lcv++) {
         cols[lcv] = ((psVector*)(xTable->data[lcv]))->data.F64;
     }
+    //Get the number of rows in the table and loop through the non-poly contributions.
     int numRows = ((psVector*)(xTable->data[0]))->n;
     for (int lcv = 0; lcv < numRows; lcv++) {
-        // arg = sum w_(i,j,k)*F_k
         double arg = 0.0;
+        //Get the argument- from the table and corresponding F-value.  Convert to radians.
         for (int k = 0; k < 14; k++) {
             arg += cols[k+3][lcv]*F[k];
         }
-        double tj = pow(t,cols[0][lcv]);
-        double as = cols[1][lcv];
-        double ac = cols[2][lcv];
-        //        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;
-    }
-
-    for (int lcv = 0; lcv < 17; lcv++) {
-        cols[lcv] = ((psVector*)(yTable->data[lcv]))->data.F64;
-    }
-    numRows = ((psVector*)(yTable->data[0]))->n;
-    for (int lcv = 0; lcv < numRows; lcv++) {
-        // arg = sum w_(i,j,k)*F_k
-        double arg = 0.0;
-        for (int k = 0; k < 14; k++) {
-            arg += cols[k+3][lcv]*F[k];
-        }
+        //The order of t for each part is specified by a column in the tab.dat files.
         double tj = pow(t,cols[0][lcv]);
         double as = cols[1][lcv];
@@ -581 +537 @@
         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;
-    }
-
+        X += (as*tj*sin(arg) + ac*cos(arg)) * tj;
+    }
+    //Do the same procedure from the previous 15 lines for Y.
     for (int lcv = 0; lcv < 17; lcv++) {
-        cols[lcv] = ((psVector*)(sTable->data[lcv]))->data.F64;
-    }
-    numRows = ((psVector*)(sTable->data[0]))->n;
+        cols[lcv] = ((psVector*)(yTable->data[lcv]))->data.F64;
+    }
+    numRows = ((psVector*)(yTable->data[0]))->n;
     for (int lcv = 0; lcv < numRows; lcv++) {
-        // arg = sum w_(i,j,k)*F_k
         double arg = 0.0;
         for (int k = 0; k < 14; k++) {
@@ -601 +554 @@
         as = SEC_TO_RAD(as) * 1e-6;
         ac = SEC_TO_RAD(ac) * 1e-6;
-
-        //        S += as*tj*sin(arg) + ac*tj*cos(arg);
+        Y += (as*tj*sin(arg) + ac*cos(arg)) * tj;
+    }
+    //Again for S.
+    for (int lcv = 0; lcv < 17; lcv++) {
+        cols[lcv] = ((psVector*)(sTable->data[lcv]))->data.F64;
+    }
+    numRows = ((psVector*)(sTable->data[0]))->n;
+    for (int lcv = 0; lcv < numRows; lcv++) {
+        double arg = 0.0;
+        for (int k = 0; k < 14; k++) {
+            arg += cols[k+3][lcv]*F[k];
+        }
+        double tj = pow(t,cols[0][lcv]);
+        double as = cols[1][lcv];
+        double ac = cols[2][lcv];
+        as = SEC_TO_RAD(as) * 1e-6;
+        ac = SEC_TO_RAD(ac) * 1e-6;
         S += (as*tj*sin(arg) + ac*cos(arg)) * tj;
     }
 
-    // now, the tables for S actually gives S + XY/2, so let's get the real S now
+    //the tables for S actually gives S + XY/2, so let's get the real S now. (from ADD)
     S -= X*Y/2.0;
 
-    //    psEarthPole* pole = psAlloc(sizeof(psEarthPole));
+    //Create the output psEarthPole and set the corresponding component values.
     psEarthPole* pole = psEarthPoleAlloc();
     pole->x = X;
@@ -621 +589 @@
                                   psTimeBulletin bulletin)
 {
+    // Check for null parameter or invalid bulletin
     PS_ASSERT_PTR_NON_NULL(time,NULL);
     PS_ASSERT_INT_WITHIN_RANGE(bulletin, PS_IERS_A, PS_IERS_B, NULL);
-
-    psEarthPole *out = psEarthPoleAlloc();
-    out->x = 0.0;
-    out->y = 0.0;
-    out->s = 0.0;
-
+    //Convert the input time to MJD.  If NAN is returned, return NULL for the function.
     double MJD = psTimeToMJD(time);
     if (MJD == NAN) {
@@ -635 +599 @@
         return NULL;
     }
-    /*
-        if (MJD < 41684.0 || MJD > 53334.0) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
-                    "Invalid time input.  Date, %lf, is out of range (41684-53334)\n", MJD);
-            return out;
-        }
-    */
+    //Allocate space for the psEarthPole output.
+    psEarthPole *out = psEarthPoleAlloc();
+
     // 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.");
@@ -653 +612 @@
     }
 
+    //Load the data table and store in the corresponding, newly-allocated vectors.
     psF64* cols[5];
     for (int colNum = 0; colNum < 5; colNum++) {
@@ -658 +618 @@
     }
     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;
-                }
-            }
-        }
-    */
     psVector *X = psVectorAlloc(numRows, PS_TYPE_F64);
     psVector *Y = psVectorAlloc(numRows, PS_TYPE_F64);
@@ -694 +635 @@
     }
 
+    //The following uses lagrange interpolation to calculate the corrections to X & Y.
     double xOut = 0.0;
     double yOut = 0.0;
@@ -724 +666 @@
         }
     }
+    //Convert the values from milli-arcsecond to radian.
     out->x = SEC_TO_RAD(xOut) * 1e-3;
     out->y = SEC_TO_RAD(yOut) * 1e-3;
@@ -736 +679 @@
 {
     int i;
-    psSphereRot *new = (psSphereRot*)psAlloc(sizeof(psSphereRot));
-    //    psMemSetDeallocator(new, (psFreeFunc)sphereRotFree);
-    //Convert rotation matrix to quaternions
+    psSphereRot *new = psSphereRotAlloc(0.0, 0.0, 0.0);
+    new->q3 = 0.0;
+
+    //Convert rotation matrix to quaternions.  Formula directly from ADD.
     double diag_sum[3];
     int maxi;
@@ -783 +727 @@
 psSphereRot* psSphereRot_CEOtoGCRS(const psEarthPole *pole)
 {
+    // Check for null parameter
     PS_ASSERT_PTR_NON_NULL(pole,NULL);
     double A[3][3];
@@ -788 +733 @@
 
     //Setup the rotation matrix and scalar value, a, as outlined by the ADD
-    //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[1][0] = -a*pole->x*pole->y*cos(pole->s) + (1.0 - a*pole->y*pole->y)*sin(pole->s);
@@ -811 +744 @@
     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);
-    */
+    //Convert the rotation matrix to a psSphereRot (quaternions)
     out = rotMatrix_To_Quat(A);
 
@@ -843 +753 @@
                                   psEarthPole *tidalCorr)
 {
+    // Check for null parameter
     PS_ASSERT_PTR_NON_NULL(time,NULL);
-    psTime *in = psTimeAlloc(time->type);
-    *in = *time;
+    //Create a copy of the input time that can be manipulated/changed (input time is const).
+    psTime *in = p_psTimeCopy(time);
     if (in->type != PS_TIME_UT1) {
         in = psTimeConvert(in, PS_TIME_UT1);
     }
+    //Check if tidal corrections should be included.
+    //If so, make sure values are positive and in the correct range.
     if (tidalCorr != NULL && tidalCorr->s != 0.0) {
         int nsec = in->nsec + (int)(tidalCorr->s * 1e9);
+        if (nsec > 1e9) {
+            nsec += -1e9;
+            in->sec += 1;
+        }
         if (nsec < 0) {
             in->sec += -1;
@@ -858 +775 @@
         }
     }
+    //Calculate the Julian Date from the input time in UT1 format.
     double T = psTimeToJD(in);
     if (T == NAN) {
@@ -865 +783 @@
     }
     T += -2451545.0;
+    //Formula for theta comes directly from the ADD.  Create output psSphereRot from theta.
     double theta = 2.0 * M_PI * (0.7790572732640 + 1.00273781191135448 * T);
     psSphereRot *out = psSphereRotAlloc(theta, 0.0, 0.0);
@@ -1003 +922 @@
 static double DMOD(double x, double y)
 {
+    //Internal function for calculating the remainder of a double quotient.
+    //used often in the algorithm for psEOC_PolarTideCorr which is the Ray model of
+    //Simon et. al. from its fortran implementation.
     double value = x - y * trunc(x/y);
     return value;
@@ -1011 +933 @@
     // Check for null parameter
     PS_ASSERT_PTR_NON_NULL(time, NULL);
-    psEarthPole *out = psEarthPoleAlloc();
-
     // Convert psTime to MJD
     double MJD = psTimeToMJD(time);
@@ -1020 +940 @@
         return NULL;
     }
-
-    // Calculate number of Julian centuries since 2000
-    double RJD = MJD;
-
-    //Formula comes from fortran reference
-    //DMOD in fortran ref. = double remainder -> x - y * trunc(x/y)
+    psEarthPole *out = psEarthPoleAlloc();
+
+    //Formula comes from fortran reference of the Ray model of Simon et. al.
     double T, L, LPRIME, CAPF, CAPD, OMEGA, THETA, CORX, CORY, CORZ;
     double ARG1, ARG2, ARG3, ARG4, ARG5, ARG6, ARG7, ARG8;
     double T2, T3, T4;
-    T = (RJD - 51544.5) / 36525.0;
+    // Calculate number of Julian centuries since 2000
+    T = (MJD - 51544.5) / 36525.0;
     T2 = T*T;
     T3 = T*T*T;
@@ -1086 +1004 @@
     CORX = SEC_TO_RAD(CORX);
     CORY = SEC_TO_RAD(CORY);
-    //    CORZ = SEC_TO_RAD(CORZ);
 
     out->x = CORX;
@@ -1099 +1016 @@
     // Check for null parameter
     PS_ASSERT_PTR_NON_NULL(time, NULL);
-    /*    if (time->type == PS_TIME_UT1) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
-                    "Invalid time input.  Time cannot be of type UT1.\n");
-            return NULL;
-        }
-    */
     // Convert psTime to MJD
     double MJD = psTimeToMJD(time);
@@ -1120 +1031 @@
     double t4 = t*t*t*t;
 
+    // The following formulae are from the ADD (& s5.8 of IERS Technical Note 32):
     double F[5];
     // Mean Anomaly of the Moon
@@ -1135 +1047 @@
            SEC_TO_RAD(0.00001149)*t4;
 
-    // L − Omega
+    // L − Omega    (L = Mean Longitude of the Moon, Omega = F4)
     F[2] = DEG_TO_RAD(93.27209062) +
            SEC_TO_RAD(1739527262.8478)*t -
@@ -1252 +1164 @@
     double X = 0.0;
     double Y = 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;
-
-    // calculate the polynomial portion first - the pj * t^j (poly coeff's)
-    //    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
@@ -1306 +1194 @@
         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);
     }
 
@@ -1316 +1201 @@
     pole->x = X;
     pole->y = Y;
+    //The value of s is simply: s = -4.7e-5 * t as specified by the ADD and IERS 32.
     pole->s = -SEC_TO_RAD(4.7e-5) * t;
 
@@ -1323 +1209 @@
 psSphereRot* psSphereRot_ITRStoTEO(const psEarthPole* motion)
 {
+    // Check for null parameter
     PS_ASSERT_PTR_NON_NULL(motion,NULL);
     double A[3][3];
@@ -1333 +1220 @@
     s = -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[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,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;
-    */
-
-    //Newest trial - mult. y matrix * x * z
-
+    //Setup Rotation matrix.  Rotation is constructed by rotation about the X-axis by y,
+    // about the Y-axis by x, and about the Z-axis by s' (where s' = -s).
     A[0][0] = cos(x)*cos(s);
     A[1][0] = cos(x)*sin(s);
@@ -1418 +1232 @@
     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);
@@ -1484 +1264 @@
     psTime *from = NULL;
     psTime *to = NULL;
+    if (fromMJD > toMJD) {
+        psWarning("From time is later than to time in psSpherePrecess.\n");
+    }
 
     if (mode == PS_PRECESS_ROUGH) {
@@ -1525 +1308 @@
             to = p_psTimeCopy(toTime);
         }
-
+        //Calculate the earthpoles and quaternions corresponding to each time (from, to).
+        //Combine the quaternions to produce the output psSphereRot.
         psEarthPole *fromEP = psEOC_PrecessionModel(from);
         psSphereRot *fromRot = psSphereRot_CEOtoGCRS(fromEP);
@@ -1557 +1341 @@
             to = p_psTimeCopy(toTime);
         }
-
+        //Calculate the earthpoles and quaternions corresponding to each time (from, to).
+        //Add in the precession corrections from IERS bulletin A.
+        //Combine the quaternions to produce the output psSphereRot.
         psEarthPole *fromEP = psEOC_PrecessionModel(from);
         psEarthPole *fromCorr = psEOC_PrecessionCorr(from, PS_IERS_A);
         fromEP->x += fromCorr->x;
         fromEP->y += fromCorr->y;
-        fromEP->s += fromCorr->s;
         psSphereRot *fromRot = psSphereRot_CEOtoGCRS(fromEP);
         psEarthPole *toEP = psEOC_PrecessionModel(to);
@@ -1568 +1353 @@
         toEP->x += toCorr->x;
         toEP->y += toCorr->y;
-        toEP->s += toCorr->s;
         psSphereRot *toRot = psSphereRot_CEOtoGCRS(toEP);
         psSphereRot *fromConj = psSphereRotConjugate(NULL, fromRot);
@@ -1600 +1384 @@
         }
 
+        //Calculate the earthpoles and quaternions corresponding to each time (from, to).
+        //Add in the precession corrections from IERS bulletin B.
+        //Combine the quaternions to produce the output psSphereRot.
         psEarthPole *fromEP = psEOC_PrecessionModel(from);
         psEarthPole *fromCorr = psEOC_PrecessionCorr(from, PS_IERS_B);
         fromEP->x += fromCorr->x;
         fromEP->y += fromCorr->y;
-        fromEP->s += fromCorr->s;
         psSphereRot *fromRot = psSphereRot_CEOtoGCRS(fromEP);
         psEarthPole *toEP = psEOC_PrecessionModel(to);
@@ -1610 +1396 @@
         toEP->x += toCorr->x;
         toEP->y += toCorr->y;
-        toEP->s += toCorr->s;
         psSphereRot *toRot = psSphereRot_CEOtoGCRS(toEP);
         psSphereRot *fromConj = psSphereRotConjugate(NULL, fromRot);
@@ -1625 +1410 @@
         return out;
     }
-    // Apply transform to coordinates
-    //    psSphere *out = psSphereRotApply(NULL, tmpST, coords);
-    //    if (out->r < -0.0001) {
-    //        out->r += 2.0 * M_PI;
-    //    }
-}
+}

trunk/psLib/src/astro/psEarthOrientation.h

@@ -9 +9 @@
 *  @author Robert Daniel DeSonia, MHPCC
 *
-*  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-01-23 20:04:31 $
+*  @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-02-14 00:09:27 $
 *
 *  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
@@ -131 +131 @@
 /** Calculates the rotation of the Earth about the CIP.
  *
- *  @return psSphereRot*:       sphereical rotation
+ *  If tidalCorr is non-NULL, use the S-component to provide tidal corrections to the
+ *  UT1 time.  If tidalCorr is NULL, no corrections are made & this step is skipped.
+ *
+ *  @return psSphereRot*:       spherical rotation of the Earth about the CIP.
  */
 psSphereRot *psSphereRot_TEOtoCEO(
@@ -187 +190 @@
  */
 psSphereRot* psSpherePrecess(
-    //    psSphere *coords,                  ///< coordinates (modified in-place)
     const psTime *fromTime,            ///< equinox of coords input
     const psTime *toTime,              ///< equinox of coords output

trunk/psLib/src/astro/psSphereOps.h

@@ -8 +8 @@
  *  @author David Robbins, MHPCC
  *
- *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-01-23 23:52:15 $
+ *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-02-14 00:09:27 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -20 +20 @@
 /// @{
 
-//#include "psCoord.h"
 #include "psTime.h"
 

trunk/psLib/src/sys/psMemory.c

@@ -8 +8 @@
 *  @author Robert Lupton, Princeton University
 *
-*  @version $Revision: 1.68 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-02-10 02:44:46 $
+*  @version $Revision: 1.69 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-02-14 00:09:27 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -1076 +1076 @@
     bool out = safeThreads;
     safeThreads = safe;
-    if ( out != p_psListThreadSafety(safe) ) {
-        psWarning("Thread safety setting in psMemory didn't match that of psList.\n");
-    }
     return out;
 }

trunk/psLib/src/types/psList.c

@@ -6 +6 @@
  *  @author Robert Daniel DeSonia, MHPCC
  *
- *  @version $Revision: 1.46 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-02-10 02:44:46 $
+ *  @version $Revision: 1.47 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-02-14 00:09:27 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -33 +33 @@
 static psBool listIteratorRemove(psListIterator* iterator);
 
-//private boolean for enabling/disabling thread safety.  Default = enabled.
-static bool safeThreads = true;
-
 static void listFree(psList* list)
 {
     if (list == NULL) {
         return;
-    }
-
-    if (safeThreads) {
-        pthread_mutex_lock(&list->p_lock);
-        //        pthread_mutex_lock(&p_lock);
     }
 
@@ -69 +61 @@
         ptr = next;
     }
-
-    if (safeThreads) {
-        pthread_mutex_unlock(&list->p_lock);
-        pthread_mutex_destroy(&list->p_lock);
-        //        pthread_mutex_unlock(&p_lock);
-        //        pthread_mutex_destroy(&p_lock);
-    } else {
-        safeThreads = true;
-    }
 }
 
@@ -110 +93 @@
     int index = iterator->index;
 
-    if (safeThreads) {
-        pthread_mutex_lock(&list->p_lock);
-        //        pthread_mutex_lock(&p_lock);
-    }
     if (elem == list->head) {        // head of list?
         list->head = elem->next;
@@ -138 +117 @@
     list->n--;
 
-    if (safeThreads) {
-        pthread_mutex_unlock(&list->p_lock);
-        //        pthread_mutex_unlock(&p_lock);
-    }
-
     // OK, delete orphaned list element and its data
     psFree(elem->data);
@@ -164 +138 @@
     psListIteratorAlloc(list,PS_LIST_HEAD,true);
 
-    if (safeThreads) {
-        pthread_mutex_init(&(list->p_lock), NULL);
-        //        pthread_mutex_init(&p_lock, NULL);
-    }
-
     if (data != NULL) {
         psListAdd(list, PS_LIST_TAIL, data);
@@ -179 +148 @@
 {
     return ( psMemGetDeallocator(ptr) == (psFreeFunc)listFree );
-}
-
-bool p_psListThreadSafety(bool safe)
-{
-    bool out = safeThreads;
-    safeThreads = safe;
-    return out;
 }
 
@@ -347 +309 @@
 
     psListElem* elem = psAlloc(sizeof(psListElem));
-
-    if (safeThreads) {
-        pthread_mutex_lock(&list->p_lock);
-        //        pthread_mutex_lock(&p_lock);
-    }
 
     // set the new list element's attributes
@@ -387 +344 @@
     }
 
-    if (safeThreads) {
-        pthread_mutex_unlock(&list->p_lock);
-        //        pthread_mutex_unlock(&p_lock);
-    }
-
     return true;
 }
@@ -426 +378 @@
 
     psListElem* elem = psAlloc(sizeof(psListElem));
-
-    if (safeThreads) {
-        pthread_mutex_lock(&list->p_lock);
-        //        pthread_mutex_lock(&p_lock);
-    }
 
     // set the new list element's attributes
@@ -466 +413 @@
     }
 
-    if (safeThreads) {
-        pthread_mutex_unlock(&list->p_lock);
-        //        pthread_mutex_unlock(&p_lock);
-    }
-
     return true;
 }

trunk/psLib/src/types/psList.h

@@ -7 +7 @@
  *  @ingroup LinkedList
  *
- *  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-02-10 02:44:46 $
+ *  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-02-14 00:09:27 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -59 +59 @@
     ///< used internally to improve performance when using indexed access, all
     ///< others are user-level iterators created by psListIteratorAlloc.
-    pthread_mutex_t p_lock;            ///< mutex to lock a node during changes
     void *lock;                        ///< Optional lock for thread safety
 }
@@ -94 +93 @@
 ;
 
-/** Private function used by psMemThreadSafety to activate/deactivate thread safety.
- *
- *  @return bool:       The previous value of the thread safety.
- */
-bool p_psListThreadSafety(
-    bool safe                          ///< current value of thread safety to set
-);
-
 /** Creates a psList linked list object.
  *