Plotting the Equation of Time: SiderealTime and SunPosition

Question

The equation of time describes the discrepancy between two kinds of solar time. The word equation is used in the medieval sense of reconcile a difference. The two times that differ are the apparent solar time, which directly tracks the motion of the sun, and mean solar time, which tracks a theoretical "mean" sun with noons 24 hours apart.

(from Wikipedia)

In the Mathematica documentation for SiderealTime, there is an "application" sample to "plot the equation of time by finding the difference between the Sun's right ascension at noon and the sidereal time at noon":

sunpos = SunPosition[GeoPosition[{0, 0}], 
   DateRange[DateObject[{2014, 1, 1, 12, 0}, TimeZone -> 0], 
    DateObject[{2014, 12, 31, 12, 0}, TimeZone -> 0], 10], 
   CelestialSystem -> "Equatorial"];

stime = SiderealTime[GeoPosition[{0, 0}], 
   DateRange[DateObject[{2014, 1, 1, 12, 0}, TimeZone -> 0], 
    DateObject[{2014, 12, 31, 12, 0}, TimeZone -> 0], 10]];

equationoftime = 
  TimeSeriesThread[With[{diff = First[#][[1]] - Last[#]},
     UnitConvert[
      Mod[diff, Quantity[24, "HoursOfRightAscension"], 
       Quantity[-12, "HoursOfRightAscension"]], 
      "MinutesOfRightAscension"]] &, {sunpos, stime}];

DateListPlot[equationoftime, FrameLabel -> Automatic]

Mathematica graphics

As you can see, the "amplitude" of the plot is different from the "amplitude" of the plot on Wikipedia. On may sources, graphical or tabular, the first maximum is less than $15^m$ and the last minimum is deeper than $-15^m$.

So, my questions (if someone knows):

What is the reason for these differences
What effects are taken into account by these astronomical functions? (nutation?) For example does SiderealTime return the LAST or the mean ST?
What are the bibliographic sources and the formula used by these astronomical functions?

If you can get hold of a copy of Meeus, you could try implementing his formulae for comparison. (I'm still sad about losing my copy.) — J. M.'s missing motivation, Mar 24 '16 at 14:19
@J.M. - Brilliant book. Three copies on shelf here, replaced with new as copies became beat up from use. That's my animation on the Wikipedia page, built long ago when teaching my then very young daughter Mathematica (and the building of a sundial at the same time). Memories... — ciao, Aug 03 '16 at 06:32
@ciao, nice! Do you still have the code for it lying around somewhere? — J. M.'s missing motivation, Aug 03 '16 at 06:37
@J.M. - Yes, I never delete old code - worst case gets archived to my robotic archive library. But it's on a machine at a home I'm seldom at and have not set up for remote access. Happy to send it when I'm there next, but could be some time... — ciao, Aug 03 '16 at 06:45

score 3 · Answer 1 · edited Apr 13 '17 at 12:55

I dug up my old (version 8) routine for the equation of time, for comparison's sake. The formulae involved are, as always, from Jean Meeus's Astronomical Algorithms.

(* Julian day number *)
Options[jd] = {"Calendar" -> "Gregorian"};

jd[{yr_Integer, mo_Integer, da_?NumericQ, rest___}, opts : OptionsPattern[]] := 
  Module[{y = yr, m = mo, h}, If[m < 3, y--; m += 12];
         h = Switch[OptionValue["Calendar"],
                    "Gregorian", (Quotient[#, 4] - # + 2) &[Quotient[y, 100]],
                    "Julian", 0,
                    _, Return[$Failed]];
         Floor[365.25 y] + Floor[30.6001 (m + 1)] + da +
         FromDMS[PadRight[{rest}, 3]]/24 + 1720994.5 + h]

jd[opts : OptionsPattern[]] := jd[DateList[], opts]

equationOfTime[date_List] := 
 Module[{t = (jd[date] - 2451545)/36525, ε, ℒ0, ℯ, ℳ☉, y},

        (* ε - mean obliquity of the ecliptic *)
        ε = (84381.406 + t (-46.836769 + t (-1.831*^-4 + t (0.0020034 + t (-5.76*^-7 -
             4.34*^-8 t))))) °/3600;

        (* ℒ0 - geometric mean longitude *)
        ℒ0 = Mod[(280.4664567 + t (36000.76982779 + t (3.032028*^-2 + t (1/49931 -
                 t (t/1.53*^4 + t/2*^6))))), 360] °;

        (* ℯ - eccentricity of Earth's orbit *)
        ℯ = 0.0167086342 + t (-0.004203654 + t (-0.00126734 + t (1.444*^-4 +
            t (-2.*^-6 + 3.*^-5 t))));

        (* ℳ☉ - mean solar anomaly *)
        ℳ☉ = Mod[(1.28710479305*^6 + t (1.295965810481*^8 + t (-0.5532 + t (1.36*^-4 -
                  1.149*^-5 t))))/3600, 360] °;

        y = Tan[ε/2]^2;
        4 (-y Sin[2 ℒ0] + 2 ℯ Sin[ℳ☉] - 4 ℯ y Sin[ℳ☉] Cos[2 ℒ0] +
           y^2 Sin[4 ℒ0]/2 + 5 ℯ^2/4 Sin[2 ℳ☉])/°]

(N.B. Altho JulianDate[] is now built-in, for some reason it was unable to reproduce the test dates from Meeus, so I fell back to my old jd[] routine, as previously seen here and here.)

Compare:

dr = DateRange[DateObject[{2014, 1, 1, 12, 0}, TimeZone -> 0],
               DateObject[{2014, 12, 31, 12, 0}, TimeZone -> 0], 10];
sunpos = SunPosition[GeoPosition[{0, 0}], dr, CelestialSystem -> "Equatorial"];
stime = SiderealTime[GeoPosition[{0, 0}], dr];
eotDocs = TimeSeriesThread[With[{diff = First[#][[1]] - Last[#]}, 
          UnitConvert[Mod[diff, Quantity[24, "HoursOfRightAscension"], 
                          Quantity[-12, "HoursOfRightAscension"]], 
                      "MinutesOfRightAscension"]] &, {sunpos, stime}];
eotMeeus = TimeSeries[Transpose[{dr, Composition[equationOfTime, DateList] /@ dr}]];
DateListPlot[{eotDocs, eotMeeus},
             FrameLabel -> Automatic, PlotLegends -> {"Docs", "Meeus"}]

equation of time from different methods

where it can be seen that Meeus's method hews more closely to the figures featured at Wikipedia. (At this juncture, I should note that the equation of time exhibits an ambiguity in sign, so if you are comparing this with other sources, it may well be that equationOfTime[] returns the negative of the convention of your reference.)

As to why there is a difference, without seeing how SunPosition[] and SiderealTime[] were implemented, there is not much to say.

Plotting the Equation of Time: SiderealTime and SunPosition

1 Answers1