PHP Science Labs - Nutation In Obliquity Of The Ecliptic

Algorithm To Compute The Nutation In Obliquity of the Ecliptic


 T  = (JD - 2451545) / 36525
 T2 = T * T
 T3 = T * T2

 DegToRad = 3.1415926535897932 / 180

 w1 = DegToRad*(297.85036 + 445267.11148*T  - 0.0019142*T2 + (T3/189474))
 w2 = DegToRad*(357.52772 + 35999.05034*T   - 0.0001603*T2 - (T3/300000))
 w3 = DegToRad*(134.96298 + 477198.867398*T + 0.0086972*T2 + (T3/56250))
 w4 = DegToRad*(93.27191  + 483202.017538*T - 0.0036825*T2 + (T3/327270))
 w5 = DegToRad*(125.04452 - 1934.136261*T   + 0.0020708*T2 + (T3/450000))

 w = cos(w5)*(92025 + 8.9*T)
 w = w + cos(2*(w4 - w1 + w5))*(5736 - 3.1*T)
 w = w + cos(2*(w4 + w5))*(977 - 0.5*T)
 w = w + cos(2*w5)*(0.5*T - 895)
 w = w + cos(w2)*(54 - 0.1*T)
 w = w + cos(w2 + 2*(w4 - w1 + w5))*(224 - 0.6*T)
 w = w + cos(w3 + 2*(w4 + w5))*(129 - 0.1*T)
 w = w + cos(2*(w4 - w1 + w5) - w2)*(0.3*T - 95)
 w = w + 200*cos(2*w4 + w5)
 w = w - 70*cos(2*(w4 - w1) + w5)
 w = w - 53*cos(2*(w4 + w5) - w3)
 w = w - 33*cos(w3 + w5)
 w = w + 26*cos(2*(w1 + w4 + w5) - w3)
 w = w + 32*cos(w5 - w3)
 w = w + 27*cos(w3 + 2*w4 + w5)
 w = w - 24*cos(2*(w4 - w3) + w5)
 w = w + 16*cos(2*(w1 + w4 + w5))
 w = w + 13*cos(2*(w3 + w4 + w5))
 w = w - 12*cos(w3 + 2*(w4 - w1 + w5))
 w = w - 10*cos(2*w4 + w5 - w3)
 w = w - 8*cos(2*w1 - w3 + w5)
 w = w + 7*cos(2*(w2 - w1 + w4 + w5))
 w = w - 7*cos(w3)
 w = w + 9*cos(w2 + w5)
 w = w + 7*cos(w3 + w5 - 2*w1)
 w = w + 6*cos(w5 - w2)
 w = w + 5*cos(2*(w1 + w4) - w3 + w5)
 w = w + 3*cos(w3 + 2*(w4 + w1 + w5))
 w = w - 3*cos(w2 + 2*(w4 + w5))
 w = w + 3*cos(2*(w4 + w5) - w2)
 w = w + 3*cos(2*(w1 + w4) + w5)
 w = w - 3*cos(2*(w3 + w4 + w5 - w1))
 w = w - 3*cos(w3 + 2*(w4 - w1) + w5)
 w = w + 3*cos(2*(w1 - w3) + w5)
 w = w + 3*cos(2*w1 + w5)
 w = w + 3*cos(2*(w4 - w1) + w5 - w2)
 w = w + 3*cos(w5 - 2*w1)
 w = w + 3*cos(2*(w3 + w4) + w5)

 dEpsDeg = w / 36000000

PHP function to compute the nutation in obliquity in decimal degrees for given JD argument.



/*

   *************************************************
   Compute the nutation in Epsilon, Delta Epsilon or
   nutation in obliquity of the ecliptic, in degrees.

   This is the correction to apply to the MEAN value
   of the obliquity to obtain the TRUE obliquity.

   The general precision level is  ±0.001 arc sec.

   This function is based on the 1980 IAU theory.

   $JD = JD number corresponding to any date and time

   $t  = Julian time factor as reckoned from J2000
   $t2 = $t to the power of 2
   $t3 = $t to the power of 3

   ==============================================
   Nutation correction elements (1980 IAU Theory)

   $w1 = Mean elongation of the moon from the sun
   $w2 = Mean anomaly of the sun
   $w3 = Mean anomaly of the moon
   $w4 = Argument of latitude of the moon
   $w5 = Longitude of lunar ascending node

   $w6 = Nutation series accumulator

*/

   Function Delta_Epsilon ($JD)
{

// Compute time factor in Julian centuries as
// reckoned from J2000 and its powers.
   $T1 = ($JD - 2451545) / 36525;
   $T2 = $T1 * $T1;
   $T3 = $T2 * $T1;

// -----------------------------------------
// Compute mean lunar elongation in radians.

   $w1 = 297.85036 + 445267.11148*$T1 - 0.0019142*$T2 + ($T3 / 189474);
   $w1 = Deg2Rad($w1);

// --------------------------------------
// Compute mean solar anomaly in radians.

   $w2 = 357.52772 + 35999.05034*$T1 - 0.0001603*$T2 - ($T3 / 300000);
   $w2 = Deg2Rad($w2);

// --------------------------------------
// Compute mean lunar anomaly in radians.

   $w3 = 134.96298 + 477198.867398*$T1 + 0.0086972*$T2 + ($T3 / 56250);
   $w3 = Deg2Rad($w3);

// ----------------------------------------------
// Compute lunar argument of latitude in radians.

   $w4 = 93.27191 + 483202.017538*$T1 - 0.0036825*$T2 + ($T3 / 327270);
   $w4 = Deg2Rad($w4);

// -----------------------------------------------------
// Compute longitude of lunar ascending node in radians.

   $w5 = 125.04452 - 1934.136261*$T1 + 0.0020708*$T2 + ($T3 / 450000);
   $w5 = Deg2Rad($w5);

// ---------------------------------------------------
// Compute the nutation in obliquity (ArcSec * 10000).

   $w6  = cos($w5)*(92025 + 8.9*$T1);
   $w6 += cos(2*($w4 - $w1 + $w5))*(5736 - 3.1*$T1);
   $w6 += cos(2*($w4 + $w5))*(977 - 0.5*$T1);
   $w6 += cos(2*$w5)*(0.5*$T1 - 895);
   $w6 += cos($w2)*(54 - 0.1*$T1);
   $w6 += cos($w2 + 2*($w4 - $w1 + $w5))*(224 - 0.6*$T1);
   $w6 += cos($w3 + 2*($w4 + $w5))*(129 - 0.1*$T1);
   $w6 += cos(2*($w4 - $w1 + $w5) - $w2)*(0.3*$T1 - 95);
   $w6 += 200*cos(2*$w4 + $w5);
   $w6 -=  70*cos(2*($w4 - $w1) + $w5);
   $w6 -=  53*cos(2*($w4 + $w5) - $w3);
   $w6 -=  33*cos($w3 + $w5);
   $w6 +=  26*cos(2*($w1 + $w4 + $w5) - $w3);
   $w6 +=  32*cos($w5 - $w3);
   $w6 +=  27*cos($w3 + 2*$w4 + $w5);
   $w6 -=  24*cos(2*($w4 - $w3) + $w5);
   $w6 +=  16*cos(2*($w1 + $w4 + $w5));
   $w6 +=  13*cos(2*($w3 + $w4 + $w5));
   $w6 -=  12*cos($w3 + 2*($w4 - $w1 + $w5));
   $w6 -=  10*cos(2*$w4 + $w5 - $w3);
   $w6 -=   8*cos(2*$w1 - $w3 + $w5);
   $w6 +=   7*cos(2*($w2 - $w1 + $w4 + $w5));
   $w6 -=   7*cos($w3);
   $w6 +=   9*cos($w2 + $w5);
   $w6 +=   7*cos($w3 + $w5 - 2*$w1);
   $w6 +=   6*cos($w5 - $w2);
   $w6 +=   5*cos(2*($w1 + $w4) - $w3 + $w5);
   $w6 +=   3*cos($w3 + 2*($w4 + $w1 + $w5));
   $w6 -=   3*cos($w2 + 2*($w4 + $w5));
   $w6 +=   3*cos(2*($w4 + $w5) - $w2);
   $w6 +=   3*cos(2*($w1 + $w4) + $w5);
   $w6 -=   3*cos(2*($w3 + $w4 + $w5 - $w1));
   $w6 -=   3*cos($w3 + 2*($w4 - $w1) + $w5);
   $w6 +=   3*cos(2*($w1 - $w3) + $w5);
   $w6 +=   3*cos(2*$w1 + $w5);
   $w6 +=   3*cos(2*($w4 - $w1) + $w5 - $w2);
   $w6 +=   3*cos($w5 - 2*$w1);
   $w6 +=   3*cos(2*($w3 + $w4) + $w5);

// -----------------------------------
// Done.  Return nutation in obliquity
// expressed in degrees.

   return $w6 / 36000000;

}

CPP function to compute the nutation in obliquity in decimal degrees for given JD argument.


/*

   ========================================================
   This CPP/C++ function computes the nutation in obliquity
   of the ecliptic (delta epsilon) in decimal degrees.

   The theoretical precision is to about ±0.001"

   Source code provided under the provisions of the
   GNU General Public License (GPL), version 3.
   http://www.fsf.org/licensing/licenses/gpl.html

   -----
   REFS.

   Astronomical Algorithms - Second Edition
   Jean Meeus - 1998
   Chapter 22
   Willmann-Bell, Inc.
   ISBN: 0-943396-61-1



   Explanatory Supplement to the Astronomical Almanac
   U.S. Naval Observatory
   ISBN 0-935702-68-7

   The 1980 IAU Theory of Nutation and Reduction
   p111-p114


   --------
   ARGUMENT

   JD = JD number for a given date and time


   ----------
   VARIABLES

   T = Julian centuries reckoned from J2000.0
       and corresponding to the JD value.

   T2 = T to the power of 2
   T3 = T to the power of 3

   w  = Nutation series accumulator
   w1 = Moon's mean elongation from the sun
   w2 = Sun's mean anomaly
   w3 = Moon's mean anomaly
   w4 = Moon's argument of latitude
   w5 = Moon's longitude of ascending node


   ------
   ERRORS

   This function does NOT check for errors.

*/

   double delta_epsilon (double JD)
{

// -------------------------
// Define working variables.

   double T, T2, T3, DegToRad;
   double w, w1, w2, w3, w4, w5;

// ---------------------------------
// Compute (T) in Julian centuries
// from J2000.0 and powers of (T).

   T = (JD - 2451545.0) / 36525.0;
   T2 = T * T;
   T3 = T * T2;

// ------------------------------------------------
// Define conversion factor for degrees to radians.
// Multiply degrees by this factor to get radians.

   DegToRad = 3.1415926535897932 / 180;

// -----------------------------------------------
// Compute mean elongation of the moon in radians.

   w1 = DegToRad*(297.85036 + 445267.11148*T
      - 0.0019142*T2 + (T3 / 189474));

// -------------------------------------------
// Compute mean anomaly of the sun in radians.

   w2 = DegToRad*(357.52772 + 35999.05034*T
      - 0.0001603*T2 - (T3 / 300000));

// --------------------------------------------
// Compute mean anomaly of the moon in radians.

   w3 = DegToRad*(134.96298 + 477198.867398*T
      + 0.0086972*T2 + (T3 / 56250));

// -----------------------------------------------
// Compute moon's argument of latitude in radians.

   w4 = DegToRad*(93.27191 + 483202.017538*T
      - 0.0036825*T2 + (T3 / 327270));

// ------------------------------------------------------
// Compute longitude of moon's ascending node in radians.

   w5 = DegToRad*(125.04452 - 1934.136261*T
      + 0.0020708*T2 + (T3 / 450000));

// ---------------------------------------------
// Compute the nutation in obliquity in terms of
// arc seconds * 10000

   w  = cos(w5)*(92025 + 8.9*T);

   w += cos(2*(w4 - w1 + w5))*(5736 - 3.1*T);
   w += cos(2*(w4 + w5))*(977 - 0.5*T);
   w += cos(2*w5)*(0.5*T - 895);
   w += cos(w2)*(54 - 0.1*T);
   w -= 7*cos(w3);
   w += cos(w2 + 2*(w4 - w1 + w5))*(224 - 0.6*T);
   w += 200*cos(2*w4 + w5);
   w += cos(w3 + 2*(w4 + w5))*(129 - 0.1*T);
   w += cos(2*(w4 - w1 + w5) - w2)*(0.3*T - 95);
   w -= 70*cos(2*(w4 - w1) + w5);
   w -= 53*cos(2*(w4 + w5) - w3);
   w -= 33*cos(w3 + w5);
   w += 26*cos(2*(w1 + w4 + w5) - w3);
   w += 32*cos(w5 - w3);
   w += 27*cos(w3 + 2*w4 + w5);
   w -= 24*cos(2*(w4 - w3) + w5);
   w += 16*cos(2*(w1 + w4 + w5));
   w += 13*cos(2*(w3 + w4 + w5));
   w -= 12*cos(w3 + 2*(w4 - w1 + w5));
   w -= 10*cos(2*w4 + w5 - w3);
   w -= 8*cos(2*w1 - w3 + w5);
   w += 7*cos(2*(w2 - w1 + w4 + w5));
   w += 9*cos(w2 + w5);
   w += 7*cos(w3 + w5 - 2*w1);
   w += 6*cos(w5 - w2);
   w += 5*cos(2*(w1 + w4) - w3 + w5);
   w += 3*cos(w3 + 2*(w4 + w1 + w5));
   w -= 3*cos(w2 + 2*(w4 + w5));
   w += 3*cos(2*(w4 + w5) - w2);
   w += 3*cos(2*(w1 + w4) + w5);
   w -= 3*cos(2*(w3 + w4 + w5 - w1));
   w -= 3*cos(w3 + 2*(w4 - w1) + w5);
   w += 3*cos(2*(w1 - w3) + w5);
   w += 3*cos(2*w1 + w5);
   w += 3*cos(2*(w4 - w1) + w5 - w2);
   w += 3*cos(w5 - 2*w1);
   w += 3*cos(2*(w3 + w4) + w5);

// -----------------------------------
// Done.  Return nutation in obliquity
// expressed in decimal degrees.

   return w / 36000000.0;

}

VB.NET function to compute the nutation in obliquity in decimal degrees for given JD argument.


 Public Function Delta_Epsilon(ByVal JD As Double) As Double

     Dim T, T2, T3, w, w1, w2, w3, w4, w5, DegToRad As Double

     DegToRad = Math.PI / 180

     ' -------------------------------------------
     ' COMPUTE TIME FACTOR (T) IN JULIAN CENTURIES
     ' AS RECKONED FROM J2000 AND POWERS OF T.

     T = (JD - 2451545.0) / 36525.0
     T2 = T * T
     T3 = T2 * T

     ' -----------------------------------------
     ' COMPUTE MEAN LUNAR ELONGATION IN RADIANS.

     w1 = DegToRad * (297.85036 + 445267.11148 _
                   * T - 0.0019142 * T2 + (T3 / 189474))

     ' --------------------------------------
     ' COMPUTE MEAN SOLAR ANOMALY IN RADIANS.

     w2 = DegToRad * (357.52772 + 35999.05034 _
                   * T - 0.0001603 * T2 - (T3 / 300000))

     ' --------------------------------------
     ' COMPUTE MEAN LUNAR ANOMALY IN RADIANS.

     w3 = DegToRad * (134.96298 + 477198.867398 _
                   * T + 0.0086972 * T2 + (T3 / 56250))

     ' ----------------------------------------------
     ' COMPUTE LUNAR ARGUMENT OF LATITUDE IN RADIANS.

     w4 = DegToRad * (93.27191 + 483202.017538 _
                   * T - 0.0036825 * T2 + (T3 / 327270))

     ' -----------------------------------------------------
     ' COMPUTE LONGITUDE OF LUNAR ASCENDING NODE IN RADIANS.

     w5 = DegToRad * (125.04452 - 1934.136261 _
                   * T + 0.0020708 * T2 + (T3 / 450000))

     ' --------------------------------------------------
     ' COMPUTE THE NUTATION IN OBLIQUITY (ARCSEC * 10000)

     w = Math.Cos(w5)*(92025 + 8.9*T)
     w += Math.Cos(2*(w4 - w1 + w5))*(5736 - 3.1*T)
     w += Math.Cos(2*(w4 + w5))*(977 - 0.5*T)
     w += Math.Cos(2*w5)*(0.5*T - 895)
     w += Math.Cos(w2)*(54 - 0.1*T)
     w += Math.Cos(w2 + 2*(w4 - w1 + w5))*(224 - 0.6*T)
     w += Math.Cos(w3 + 2*(w4 + w5))*(129 - 0.1*T)
     w += Math.Cos(2*(w4 - w1 + w5) - w2)*(0.3*T - 95)
     w += 200*Math.Cos(2*w4 + w5)
     w -= 70*Math.Cos(2*(w4 - w1) + w5)
     w -= 53*Math.Cos(2*(w4 + w5) - w3)
     w -= 33*Math.Cos(w3 + w5)
     w += 26*Math.Cos(2*(w1 + w4 + w5) - w3)
     w += 32*Math.Cos(w5 - w3)
     w += 27*Math.Cos(w3 + 2*w4 + w5)
     w -= 24*Math.Cos(2*(w4 - w3) + w5)
     w += 16*Math.Cos(2*(w1 + w4 + w5))
     w += 13*Math.Cos(2*(w3 + w4 + w5))
     w -= 12*Math.Cos(w3 + 2*(w4 - w1 + w5))
     w -= 10*Math.Cos(2*w4 + w5 - w3)
     w -= 8*Math.Cos(2*w1 - w3 + w5)
     w += 7*Math.Cos(2*(w2 - w1 + w4 + w5))
     w -= 7*Math.Cos(w3)
     w += 9*Math.Cos(w2 + w5)
     w += 7*Math.Cos(w3 + w5 - 2*w1)
     w += 6*Math.Cos(w5 - w2)
     w += 5*Math.Cos(2*(w1 + w4) - w3 + w5)
     w += 3*Math.Cos(w3 + 2*(w4 + w1 + w5))
     w -= 3*Math.Cos(w2 + 2*(w4 + w5))
     w += 3*Math.Cos(2*(w4 + w5) - w2)
     w += 3*Math.Cos(2*(w1 + w4) + w5)
     w -= 3*Math.Cos(2*(w3 + w4 + w5 - w1))
     w -= 3*Math.Cos(w3 + 2*(w4 - w1) + w5)
     w += 3*Math.Cos(2*(w1 - w3) + w5)
     w += 3*Math.Cos(2*w1 + w5)
     w += 3*Math.Cos(2*(w4 - w1) + w5 - w2)
     w += 3*Math.Cos(w5 - 2*w1)
     w += 3*Math.Cos(2*(w3 + w4) + w5)

     ' -----------------------------------
     ' DONE.  RETURN NUTATION IN OBLIQUITY
     ' EXPRESSED IN DECIMAL DEGREES.

     Return w / 36000000.0

 End Function

Nutation In Obliquity Of The Ecliptic - IAU 1980 Series
Date: ±Y M D	Time UT: H M S	ΔT: ± H M S

Calendar Mode: Julian Gregorian