Quirks of Time

(*Authors to supply subtitle with edits.*)

by L. R. King and Todd G. Will

Introduction.

The day with the fewest number of daylight hours, sometimes referred to as the shortest day of the year, is winter solstice, which occurs in the northern hemisphere on December 21st or 22nd.  Given that the number of daylight hours increases after the winter solstice one should naturally expect that the sun would rise earlier on subsequent days.   Surprisingly this is not the case.  Depending on latitude, the sun continues to rise later each day for as much as several weeks past winter solstice despite the fact that the number of daylight hours steadily increases.  We refer to this curious phenomenon as Delayed Early Sunrise (DES) and offer an explanation for it based solely on the fact that the earth's polar axis is tilted, making an approximate angle of 23.5° with its orbital plane.

The explanation comes in two parts.  First, we establish the connection between the time of sunrise and the time of high noon.   By high noon we mean the time when the sun is highest in the sky.  Second, we establish the connection between high noon and twelve o'clock noon.  It is the discrepancy between high noon and twelve o'clock noon (quirk of time!) that accounts for delayed early sunrise in our model.

Choosing coordinates

We choose a coordinate system fixed to the stars with earth at the origin and the z-axis aligned with the north pole.

SetOptions[
        Graphics3D, 
PlotRange -> {{-3, 
              3}, {-3, 3}, {-2, 
              2}}, 
LightSources -> {}, 
Boxed -> 
          False, 
AmbientLight -> 
          White, 
ViewPoint -> {4, 0, 
            1}, 
BoxRatios -> 
          Automatic];

Show[earth, equatorialPlane, xyzAxes, 
        PlotRange -> All];

In this coordinate system the sun orbits the earth counter-clockwise as viewed from the north pole.  To keep the graphics uncomplicated we depict a 12 day year with time t = 0 corresponding to twelve o'clock noon on day 0.  As a further simplification we initially ignore the tilt of the orbital plane and show the sun orbiting in the equatorial plane.

Clear[sun, days]
\(\(sun[t_, days_:  365.25, tilt_:  23.5  Degree] = \ 
        RotationMatrix3D[\(-Pi\)/2, tilt, 
            Pi/2] . {Cos[\(\(2  Pi\)\/days\) t], Sin[\(\(2  Pi\)\/days\) t], 
            0};\)
\)
Clear[sunOrbit]
sunOrbit[t_, days_:  365.25, tilt_:  23.5  Degree, 
        col_:  0.3] := 
{orbitalPlane[tilt, 
          col],  (*Defined in initialization *) 
Graphics3D[{Hue[col, 1, 
              0.6], 
dateLabel[t, days, 
              tilt],   (*Defined in initialization *) 
{PointSize[0.01], Point[#], PointSize[0.02], Point[2.7  #], 
                  Line[{#, 2.7  #}]} & /@   Table[
                sun[s, days, tilt], {s, 0, t, 
                  days/12}]}]
}

noonMeridian[t_, 
        days_] := 
Graphics3D[{Red, Thickness[0.003], 
          Line[Table[
              1.01  RotationMatrix3D[((-2 π t)/days), 0, 
                    0] . {Cos[Pi s], 0, Sin[Pi s]}, {s, -1/2, 1/2, 
                0.1}]]
}]

<< Graphics`Animation`
SetOptions[Animate, Closed -> True, 
        Frames -> 12];
Animate[{  earth, sunOrbit[t, 12, 0, 0], noonMeridian[t, 12]}, {t, 
          0, 12}];

In addition to the sun's orbit about the earth, the animation gives a potentially misleading picture of the rotation of the earth on its axis.  Any confusion about this motion can be cleared up by focusing on the earth's line of longitude shown in red.  We refer to this meridian as the twelve o'clock meridian since the animation shows its location at twelve o'clock noon on each day.

Note that at twelve o'clock noon on day 0 when t = 0, the sun appears highest in the sky to folks living on the twelve o'clock meridian.  During the first day, as the sun proceeds through 1/12 of its orbit, the animation seems to show the earth making 1/12 of a revolution on its axis.  In fact, during the first day the earth rotates counter-clockwise on its axis 13/12 of a revolution causing the sun to set in the west, rise in the east and again appear directly overhead to viewers on the twelve o'clock meridian at twelve o'clock noon on day 1.  In the next 24 hours the earth rotates another 13/12 = 1 + (1/12) of revolution, the extra 1/12 of a revolution again allowing the twelve o'clock meridian to catch up to the new location of the sun at twelve o'clock noon.

Tracking the twelve o'clock meridian gives us a connection between clock time and astronomical reality.  In our currently simplified model, in which the sun orbits in the equatorial plane, clock time also appears to be connected to the orbit of the sun.  In particular we note in this untilted model that throughout the year viewers on the twelve o'clock meridian experience high noon at twelve o'clock noon.  For this reason we refer to the sun in this simplified model as the "clock-sun".  In fact, the clock sun for a normal 365.25 day year gives the correct time if we choose the 12 o'clock meridian to be the prime meridian (0 degrees longitude).  We will see shortly that orbital tilt complicates the relationship between high noon and twelve o'clock noon, introducing what we call the quirk of time.

Now look at the more realistic orbit of the sun incorporating the 23.5° tilt of the orbital plane.

Animate[{ earth, equatorialPlane,  sunOrbit[t, 12]}, {t, 0, 
          12}];

The orbital plane has been tipped along the axis defined by the equinoxes so that winter solstice occurs at twelve o'clock noon on day 0. This puts the summer solstice on day 6 and the spring and fall equinoxes on days 3 and 9, respectively.  Note how the declination of the sun, that is the angle the sun makes with the equatorial plane, varies from minimum of -23.5° at winter solstice to a maximum of 23.5  ° at summer solstice.  It is this seasonal variation that makes sunrise times latitude dependent.

Here is a look at how the sun illuminates the earth throughout the year.

litEarth[t_, days_:  365.25, tilt_:  23.5  Degree] := 
        Module[{sunDir = sun[t, days, tilt], n1, 
            mat}, 
n1 = 
            Cross[sunDir, {0, 0, 1.0}]; 
n1 = 
            n1/Sqrt[n1 . n1]; 
mat = 
            Transpose[{n1, Cross[n1, sunDir], 
                sunDir}]; 
{halfLight /. {x__Real} :> mat . {x},    (* 
              halfLight defined in initialization *) 
WireFrame[
              Sphere[1, 12, 16]]
}];

Animate[{ equatorialPlane, litEarth[t, 12],  
        sunOrbit[t, 12]}, 
{t, 0, 12}]

Here's the same animation viewed looking down from the north pole.

Animate[{ equatorialPlane, litEarth[t, 12],  
        sunOrbit[t, 12]}, 
{t, 0, 
        12}, 
ViewPoint -> {0, 0, 
          5}, 
ViewVertical -> {-1, 0, 
          0}]

The animation above shows seasonal variations in the lighting of the earth, NOT the daily variations of night and day.  Night and day are caused by the undepicted rotation of the earth on its axis.

One can see that at northern latitudes the fewest hours of daylight occur on the day of winter solstice and the greatest number on the day of summer solstice.  Consequently the angle from the sun to the location of sunrise (the boundary of blue and yellow) is smallest at winter solstice and greatest at summer solstice.  Since the earth rotates on its axis at a constant rate, focusing on this angle will provide a connection between sunrise time and the time of high noon.  To compute this angle we first need to compute the declination of the sun.

To find the declination at any time t first note that our parametrization for the sun can be viewed as a unit vector.  Thus the cosine of the angle between the sun and the north pole can be computed as the dot product of the sun vector and <0,0,1>.  Subtracting this angle from (π/2) gives the declination for any time t, declination[t].

declination[t_, days_:  365.25, tilt_:  23.5  Degree] = 
      FullSimplify[ 
        Pi/2 - ArcCos[sun[t, days, tilt] . {0, 0, 1}]]

From sunrise to high noon.

As mentioned above, the connection between high noon and sunrise will be made by focusing on the angle that the earth rotates through from sunrise to high noon.  We call this angle ρ and note that it depends on both the latitude λ and the declination of the sun which is determined by time.  The following illustrates the angle ρ at winter solstice for latitude λ≈25° N.

\(\(Show[
litEarth[0, 
          12], 
Graphics3D[{Red, 
            Thickness[
              0.005], 
Text["ρ", {0.2  Cos[0.5], 
                0.2  Sin[0.5], 1}], 
Line[
              Join[{{0, 0, 1}}, 
                Table[{0.7  Cos[t], 0.7  Sin[t], 1}, {t, 0, 1.16, 
                    1.16/6}], {{0, 0, 
                    1}}]]
}], 
ViewPoint -> {0, 0, 5}, 
ViewVertical -> {\(-1\), 0, 
            0}, 
PlotRange -> All];\)\)

In the previous section we developed a function, declination[t], to compute the declination of the sun at any time t.  For simplicity we fix the sun in the xz-plane and use the declination function to create a unit vector that points to the sun with the correct declination at time t:  
	{Cos[declination[t]], 0, Sin[declination[t]]}.
	
In our coordinate system the points on earth at latitude λ are swept out by the vector 
 	{Cos[θ]Cos[λ],Sin[θ]Cos[λ],Sin[λ]} as θ  ranges from 0 to 2π.

Sunrise and sunset occur at this latitude at the values of θ for which these two vectors are perpendicular.  Mathematica finds both of these values, the second of which corresponds to sunrise.

Clear[days, tilt]
ρ[t_, λ_, days_:  365.25, tilt_:  23.5  Degree] = 
      FullSimplify[Solve[
            0 == {Cos[θ] Cos[λ], Sin[θ] Cos[λ], 
                  Sin[λ]} . {Cos[declination[t, days, tilt]], 0, 
                  Sin[declination[t, days, 
                      tilt]]}, θ]
[[2, 1, 
            2]] ]

Given that the earth rotates on its axis roughly once every 24 hours we have that sunrise on day t should occur roughly (24/(2  π)) ρ[t, λ] hours before high noon on day t.

There are two small errors in our reasoning.  First, given our coordinate system, the earth rotates 2π ((days + 1)/days) radians each 24 hours (days = number of days in a year).  Second, the declination of the sun at high noon, which we use to predict sunrise time, differs slightly from the declination of the sun at sunrise, which occurred roughly six hours prior to high noon.  Although these errors would be significant for the 12 day year we display in our graphics, they are both negligible in the calculations we perform for the real year with 365 days.

From high noon to twelve o'clock noon.

We saw above that sunrise at latitude λ on day t occurs roughly (12/π) ρ[t, λ] hours before high noon.  Here's a plot for our 365 day year showing sunrise times at 31° N  latitude assuming that high noon occurs at twelve o'clock noon.

Plot[ 
        12 - (12/π) ρ[t, 31  Degree], {t, -200, 
          200}];

Recalling that day t = 0 corresponds to the winter solstice we see that there is no Delayed Early Sunrise.  That is, the latest sunrise occurs on the winter solstice at approximately 7 a.m. and earlier sunrises start immediately on the following days.  The reality of DES forces the conclusion that high noon and twelve o'clock noon do not always coincide.

As stated in the beginning our explanation for DES rests on the tilt of the orbital plane.  To see the effect of the tilt look at the following plot which depicts the orbit of the sun with and without the orbital tilt.

Show[
        earth, 
sunOrbit[11, 12, 0, 
          0], 
sunOrbit[11, 12, 40  Degree, 
          0.3]
];

The red graphics show a 12 day orbit with no orbital tilt.  The green graphics show a 12 day orbit with an exaggerated orbital tilt of 40°.  Look at the same plot from above the north pole.

Show[
        earth, 
sunOrbit[11, 12, 0, 
          0], 
sunOrbit[11, 12, 40  Degree, 
          0.3], 
ViewPoint -> {0, 0, 
            5}, 
ViewVertical -> {-1, 0, 
            0}];

Recall that the labels 0 through 11 mark the location of the sun at twelve o'clock noon on each day.  The red graphics depict the fictitious clock-sun which orbits the earth in the equatorial plane.  The red lines depict the positions of the twelve o'clock meridian. Recall that at twelve o'clock noon on each day the location of the clock-sun coincides with the location of twelve o'clock meridian.

The green graphics show an effect of an exaggerated orbital tilt of 40°. The green lines depict the positions of lines of longitude along which viewers experience high noon.  The mis-match of the two sets of lines is evident in the graphic.  At time t = 0 both the red sun and green sun are at high noon for viewers on the twelve o'clock meridian.  In 24 hours at twelve o'clock noon on day 1 the twelve o'clock meridian has made 13/12 of a revolution to once again align with the fictitious clock-sun.  But  at this time the green sun is not overhead.  The earth must continue its counter-clockwise rotation for several more minutes before the green sun is overhead.  Thus on the days following winter solstice high noon occurs after twelve o'clock noon.  The graphic shows that the same behavior continues until the spring equinox when once again high noon and twelve o'clock noon coincide.

At twelve o'clock noon on day 4 the folks on the twelve o'clock meridian have rotated past the green sun and thus past the position of high noon.  Thus on days 4 and 5 high noon occurs before twelve o'clock noon.  Days 7 and 8 following summer solstice are similar to those following winter solstice in that high noon again occurs after twelve o'clock noon.  On days 10 and 11 following the spring equinox high noon again occurs before twelve o'clock noon.

To measure the time between high noon and twelve o'clock noon we will measure the angle the earth rotates to reach one from the other.  We refer to this angle as quirk[t] and to measure it we begin by projecting the location of the tilted sun onto the equatorial plane and then normalizing.

<< LinearAlgebra`Orthogonalization`
projectedSun[t_, days_, tilt_]=Normalize[  (({{1, 0, 0}, {0, 1, 0}, {0, 0, 0}})).sun[t, days, tilt]]

The direction of the twelve o'clock meridian on day t is given by {Cos[ ((2 π t)/days)], 
        Sin[((2 π t)/days) ], 0}.  We could get the unsigned angle between high noon and the twelve o'clock meridian by computing the ArcCosine of these two vectors.  However the following trick of adding and subtracting (π/2) will give us the signed angle that we need.

quirk[t_, days_:  365.25, tilt_:  23.5  Degree] = 
      FullSimplify[
(π/2) - 
          ArcCos[{Cos[((2 π t)/days) + (π/2)], 
                Sin[((2 π t)/days) + (π/2)], 0} .  
              projectedSun[t, days, tilt]]
]

Again relying on the constant rate of the earth's rotation we plot the minutes from twelve o'clock noon to high noon on the days following the winter solstice.

Plot[((60 24)/(2  π)) quirk[t], {t, 0, 365.25}, 
        PlotLabel -> "Minutes from 12:00 to high noon"];

Finally we see the root cause for DES.  The tilt of the orbital plane causes the sun to be increasingly tardy relative to the clock for nearly 50 days following the winter solstice.  To see the effect on sunrise we update our previous formula,  sunrise[t, λ] = 
        highnoon[t]  - (12/π) ρ[t, λ] , replacing highnoon with 12 + (12/π) quirk[t]

Clear[sunrise]
sunrise[t_, λ_, days_:  365.25, 
          tilt_:  23.5  Degree] = 
12 + (12/π) 
            quirk[t, days, tilt]  - (12/π) ρ[t, λ, days, 
              tilt];

Here is a plot showing sunrise times for various latitude on the days following winter solstice.

\(\(Show[
GraphicsArray[
Partition[
Table[
Plot[
                sunrise[t, λ\ Degree], {t, \(-10\), 
                  50}, 
DisplayFunction -> 
                  Identity, 
PlotLabel -> "Latitude " <> ToString[λ]], 
{λ, 0, 60, 
                12}], 
3]
], 
        DisplayFunction -> $DisplayFunction
];\)\)

We have seen that after winter solstice both quirk[t] and ρ[t,λ] increase.  Thus the change in  
		sunrise[t] = 
        12 + (12/π) 
            quirk[t]   - (12/π) ρ[
              t, λ]  
is determined by the relative change in quirk and ρ.  
As long as the increase in quirk is larger than the increase in ρ, sunrise time continues to be later after winter solstice.  The delay in earlier sunrise is greatest at the equator where ρ doesn't change at all.  Because the increase in ρ becomes more pronounced with greater latitude, the delayed early sunrise is of shorter duration for higher latitudes.

About Our Model.

Our analysis of delayed early sunrise, concentrating as is does on the tilt of the earth , omits consideration of the earth's non-uniform motion along its elliptical orbit.   The earth is moving at its fastest in early January when it is closest to the sun (perihelion), gradually slowing until early July when it is most distant from the sun (aphelion) at wich time the process is reversed. At perihelion the earth's rotation needs a little more time to catch up to the sun, tending to make the sun late; at apehelion just the opposite happens. Because aphelion and perihelion are off-set from the solstices, the two major influences on quirk - non-uniform motion and  tilt - are accordingly off-set. The result is the nice symmetry of our plot of quirk[t] is lost. Moreover, the "real" plot of quirk[t] differs from ours in that the amplitude is greater (about 15 versus 10  minutes). Other (finer) points: we set our clock at 12:00 noon precisely at winter solstice for computational convenience and played the role of an observer at 0 degrees longitude (12 o'clock meridian), avoiding corrections due to civil time zones.

Graphics Primitives

<< Graphics`Shapes`
\(<< Graphics`Colors`
\)
ring = ParametricPlot3D[
            r {Cos[ t], Sin[ t], 0}, {r, 2.3, 2.8}, {t, 0, 
              2  Pi}, 
PlotPoints -> {2, 25}, 
            DisplayFunction -> 
              Identity][[1]];

equatorialPlane = 
        Graphics3D[{SurfaceColor[Hue[0, 0.2, 1]], 
            ring}];

equator = 
        Graphics3D[{
Thickness[0.007], 
            Blue, 
Line[ 
              Table[{Cos[Pi t], Sin[Pi t], 0}, {t, 0, 2, 
                  1/12}]]
}];

northPole = 
        Graphics3D[{
Thickness[0.007], 
            Blue, 
Line[{{0, 0, 1}, {0, 0, 
                  1.3}}]
}];

earth = {equator, 
          northPole, 
Graphics3D[{SurfaceColor[SkyBlue], 
              Sphere[1, 12, 
                16]}]
};

xyzAxes = 
        Graphics3D[{
Text["x", {4.4, 0, 
                0}], 
Text["y", {0, 3.2, 
                0}], 
Text["z", {0, 0, 
                2.2}], 
Gold, 
            Thickness[
              0.007], 
Line[{{1, 0, 0}, {4, 0, 
                  0}}], 
Line[{{0, 1, 0}, {0, 3, 
                  0}}], 
Line[{{0, 0, 1}, {0, 0, 
                  2}}]}];


dateLabel[t_, days_, tilt_] := 
    Table[Text[s, (3  sun[s, days, tilt])], {s, 0, t, 
        days/12}]



orbitalPlane[tilt_, col_] := 
RotateShape[
          Graphics3D[{SurfaceColor[Hue[col, 0.2, 1]], 
              ring}], 
0, tilt, 
          Pi/2];


halfLight = 
        Graphics3D[
          Join[
{EdgeForm[], 
              SurfaceColor[Yellow]}, 
ParametricPlot3D[
                Evaluate[0.97 {Cos[s] Cos[u], Sin[s] Cos[u], Sin[u]}], {s, 
                  0, 2  Pi}, {u, 0, Pi/2}, PlotPoints -> {17, 10}, 
                DisplayFunction -> 
                  Identity][[1]], 
{EdgeForm[], 
              SurfaceColor[Blue]}, 
ParametricPlot3D[
                Evaluate[0.97 {Cos[s] Cos[u], Sin[s] Cos[u], Sin[u]}], {s, 
                  0, 2  Pi}, {u, -Pi/2, 0}, PlotPoints -> {17, 10}, 
                DisplayFunction -> 
                  Identity][[1]]
]];

References

Article
King, L. R. 1996. Clock Time vs. Sun Time, The UMAP Journal, 17.2: 123-144.

Websites
U. S. Naval Observatory Data Services: 	http://aa.usno.navy.mil/AA/data/
The Analemma:				http://www.analemma.com/

About the authors

L. R. King and  Todd G. Will
Department of Mathematics
Davidson College
Box 1719
Davidson NC 28036

ELECTRONIC SUBSCRIPTIONS

Included in the distribution for each electronic subscription is the file qirks.nb containing Mathematica code for the material described in this article.