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

Timestamp:

Nov 16, 2005, 5:59:05 PM (21 years ago)

Author:

drobbin

Message:

Implemented function NutationCorr and test. Has yet to be verified.

File:

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

trunk/psLib/src/astro/psEarthOrientation.c

@@ -8 +8 @@
  *  @author Robert Daniel DeSonia, MHPCC
  *
- *  @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-11-17 00:07:25 $
+ *  @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-11-17 03:59:05 $
  *
  *  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
@@ -493 +493 @@
     S -= X*Y/2.0;
 
-    psEarthPole* pole = psAlloc(sizeof(psEarthPole));
+    //    psEarthPole* pole = psAlloc(sizeof(psEarthPole));
+    psEarthPole* pole = psEarthPoleAlloc();
     pole->x = X;
     pole->y = Y;
@@ -850 +851 @@
 psEarthPole* psEOC_NutationCorr(psTime *time)
 {
-    return NULL;
+    // 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);
+
+    // Calculate number of Julian centuries since 2000
+    double t = ( MJD - MJD_2000 ) / JULIAN_CENTURY;
+    double t2 = t*t;
+    double t3 = t*t*t;
+    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];
+    F[0] = DEG_TO_RAD(134.96340251) +
+           SEC_TO_RAD(1717915923.2178)*t +
+           SEC_TO_RAD(31.8792)*t2 +
+           SEC_TO_RAD(0.051635)*t3 -
+           SEC_TO_RAD(0.00024470)*t4;
+
+    // Mean Anomaly of the Sun
+    F[1] = DEG_TO_RAD(357.52910918) +
+           SEC_TO_RAD(129596581.0481)*t -
+           SEC_TO_RAD(0.5532)*t2 +
+           SEC_TO_RAD(0.000136)*t3 -
+           SEC_TO_RAD(0.00001149)*t4;
+
+    // L − Omega
+    F[2] = DEG_TO_RAD(93.27209062) +
+           SEC_TO_RAD(1739527262.8478)*t -
+           SEC_TO_RAD(12.7512)*t2 -
+           SEC_TO_RAD(0.001037)*t3 +
+           SEC_TO_RAD(0.00000417)*t4;
+
+    // Mean Elongation of the Moon from the Sun
+    F[3] = DEG_TO_RAD(297.85019547) +
+           SEC_TO_RAD(1602961601.2090)*t -
+           SEC_TO_RAD(6.3706)*t2 +
+           SEC_TO_RAD(0.006593)*t3 -
+           SEC_TO_RAD(0.00003169)*t4;
+
+    // Mean Longitude of the Ascending Node of the Moon
+    F[4] = DEG_TO_RAD(125.04455501) -
+           SEC_TO_RAD(6962890.5431)*t +
+           SEC_TO_RAD(7.4722)*t2 +
+           SEC_TO_RAD(0.007702)*t3 -
+           SEC_TO_RAD(0.0000593)*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 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)
+    // 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);
+    for (int i = 0; i < 10; i++) {
+        //        tj = SEC_TO_RAD(pow(t, 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] * tj * sin(arg) + xp_cos[i] * tj * cos(arg);
+        Y += yp_sin[i] * tj * sin(arg) + yp_cos[i] * tj * cos(arg);
+    }
+
+    psEarthPole *pole = psEarthPoleAlloc();
+    pole->x = X;
+    pole->y = Y;
+    pole->s = SEC_TO_RAD(4.7e-5) * t;   //XXX:  This conv. should include t?
+
+    return pole;
 }