Variation of the Equation of Time with the Precession of the Equinoxes

(*
The original Mathematica code for the Clock of the Long Now's Equation of Time Cam was by W. Daniel Hillis after
Calendrical Calculations By Nachum Dershowitz and Edward M.Reingold.Cambridge University Press,1997.
*)

(*
A ratioan appoximation to the length of the tropical year (JD2000)
*)

tropicalyear = 365.242189 ;

(*
a rational approximation for the ratio of synodic month to tropical year
*)

synodicmonth = 29.5305888531 - 0.000000216 d / 36525;

start = 36525;

solarday = 
  (86400 - 0.0015 T)/86400 /. T -> d / 36525;

(*
How far off is the month calcation? What is the best constant?
*)

montherror[jdays_] := 
  Integrate[ 1/ (synodicmonth /. d -> start ), 
    {d, start, (start + jdays)}] - 
  Integrate[ (solarday/synodicmonth) ,
    {d, start, (start + jdays)}]

jcen= 3652425;

newmontherror[jdays_] := .007 + 
  Integrate[ 0.033863174755042884,
    {d, start, (start + jdays)} ]- Integrate[ (solarday/synodicmonth), 
    {d, start, (start + jdays)}]

montherror[jcs_] := 
  Integrate[ 36525/(29.5305888531 - 0.000000216) ,
        {T, 1, (1 + jcs)} ] - 
  Integrate[ ( 36525 ( 86400 - 0.0015 T) / 86400) /
        (29.5305888531 - 0.000000216 T),
        {T, 1, (1 + jcs )}]

newmontherror[jcs_] := .000325521 + 
   Integrate[ 1236.8529296875, {T, 1, (1 + jcs)}] - 
   Integrate[ (36525 ( 86400 - 0.0015 T)/86400) /
         (29.5305888531 - 0.000000216 T),
         {T, 1, (1 + jcs)}]

(*
utility functions
*)

Poly[x_, a_] :=
	a[[1]] +Sum[a[[i]] x^(i - 1), {i, 2, Length[a]}] 

DegreesMod360[theta_] :=
	Mod[theta, 360]

RadiansToDegrees[theta_] :=
	DegreesMod360[theta / Pi 180]

DegreesToRadians[theta_] :=
	DegreesMod360[theta] Pi / 180

SinDegrees[theta_] :=
	Sin[DegreesToRadians[theta]]

CosineDegrees[theta_] :=
	Cos[DegreesToRadians[theta]]

TangentDegrees[theta_] :=
	Tan[DegreesToRadians[theta]]

ArcTanDegrees[x_, quad_] :=
	Module[{deg},
		deg = RadiansToDegrees[ArcTan[x]];
		If[quad == 1 || quad == 4, deg, deg + 180] 	]

ArcSinDegrees[x_] :=
	RadiansToDegrees[ArcSin[x]]

ArcCosDegrees[x_] :=
	RadiansToDegrees[ArcCos[x]]

(*
This is the julian day number of Noon, Jan. 1, 2000
*)
J2000 = 2451545.0;

(*
equation of time 


(in Julian days) (adapted from Dershowitz & Reingold)
*)

EquationOfTime[jd_] :=
	Module[{c, longitude, anomaly, inclination, eccentricity, y}, 		
  c = (jd - J2000) / 36525;
		longitude = Poly[c, {280.46645, 36000.76983, 0.0003032}]; 		
  anomaly =Poly[c, {357.52910, 35999.05030, -0.0001559,-0.00000048}]; 		
  inclination =
   Poly[c, {23.43929111, -0.013004167, -0.00000016389,0.0000005036}]; 		
  eccentricity =Poly[c, {0.016708617, -0.000042037, -0.0000001236}]; 		
  y = TangentDegrees[inclination / 2]^2;
  N[(y SinDegrees[2 longitude] + -2 eccentricity SinDegrees[anomaly] +
	4 eccentricity y SinDegrees[anomaly]CosineDegrees[2 longitude] +
    -0.5 y^2 SinDegrees[4 longitude] +
	-1.25 eccentricity^2 SinDegrees[2 anomaly])/ (2 Pi)]
	]

cam[y_,c_,f_] := {(2 +f) Sin[2 Pi y],(2 + f) Cos[2 Pi y], c/10}
rotrate = N[(171 173)/(9 9)];
MakePolygons[vl_List] :=
  Block[{l = vl,
      l1 = Map[RotateLeft, vl],
      mesh},
    mesh = {l, RotateLeft[l], RotateLeft[l1], l1};
    mesh = Map[Drop[#, -1]&, mesh, {1}];
    mesh = Map[Drop[#, -1]&, mesh, {2}];
    Polygon [ Transpose[ Map[Flatten[#, 1]&, mesh] ]]
  ];

facePolys = 
  MakePolygons[
   Table[
    N@Evaluate[
     cam[y, c/2, 
      (30 24)/30 EquationOfTime[((c 100 + y)) rotrate ]]]
      , {y, 0, 1, 1/52}, {c, (-5), 105, 2}]];

faceLine[c_] := Line[Table[
			Evaluate[
     		cam[y, c/2,
      		(30 24)/30 EquationOfTime[((c 100 + y)) rotrate ]]
           ],
      		{y, 0, 1, 1/365}]];

Manipulate[
Module[{c = (year-2000)/100}, 
Grid[{{Plot[60 24  EquationOfTime[t],{t,J2000+year-2000,J2000 + year - 2000 + 365}, ImageSize->{295,250},Axes->{True,False},Frame->True,FrameLabel->{None,"minutes","year "<>ToString[ year],None},Exclusions->None,FrameTicks->{None,{-15, -10, -5, 0, 5, 10, 15},None,{-15, -10, -5, 0, 5, 10, 15}}],Graphics3D[{EdgeForm[], facePolys, faceLine[c]}, ImageSize->{295,250}, Axes->{False, False, True}, AxesLabel->{None, None, "year"},Ticks->None]}}]], 
{{year, 02000, "year"}, 02000, 12000,1,
Appearance->"Labeled"},SaveDefinitions->True]

CAPTION

The Clock of the Long Now is a mechanical computer designed to keep accurate time for 10,000 years. The Clock relies on the constancy of the length of a day to keep it in sync via an optical trigger that fires at local solar noon. The equation of time converts between local solar time and the Earth's rotational period. The version of the equation of time implemented here accounts for the variation of sidereal time with the precession of the equinoxes over a period of 26,000 years. The cylindrical plot of the equation of time shown here in 3D is implemented as a mechanical cam that adjusts the timing of the optical sync trigger.

DETAILS

For more information, visit The Long Now Foundation website.

THIS NOTEBOOK IS THE SOURCE CODE FROM

"Variation of the Equation of Time with the Precession of the Equinoxes" from the Wolfram Demonstrations Project
 http://demonstrations.wolfram.com/VariationOfTheEquationOfTimeWithThePrecessionOfTheEquinoxes/

Contributed by: Stewart Dickson

A full-function Wolfram Mathematica system (Version 6 or higher) is required to edit this notebook.
GET WOLFRAM MATHEMATICA »

© Wolfram Demonstrations Project & Contributors   |   Terms of Use   |   Make a new version of this Demonstration »