Sundialing with Mathematica

Construction of a naive sundial

The shadow's path

If the hand of a sundial is parallel to the axis of the earth then the hour lines are the projections of planes through the bar onto the surface where the shadow is formed. The hand can be either a thin bar or, in the limit, just its end. The planes are separated by 15° for every hour.

The dial

The dial is in the shape of a hyperbolic torus; here are its parts.

\(\(R\_*\) = 2; dψ = π\/8; ppl = 50; ppq = 12; \ 

    ξ = 1\/\@3; 
    η = 1\/2\ 
        \((Cos[\(\((ψ + π\/2)\)\ π\)\/dψ] + 1)\); \)

SetOptions[ParametricPlot3D, DisplayFunction -> Identity];

\(addColor[xyz_, ψ_] := 
      N[Append[xyz, 
          SurfaceColor[Hue[0.35 + 0.223  ψ], Hue[0.35 + 0.223  ψ], 
            1.67]]]; 


    addColor\_l[xyz_] := N[Append[xyz, SurfaceColor[Hue[0], Hue[0], 1.67]]]; 
    

addColor\_r[xyz_] := 
      N[Append[xyz, SurfaceColor[Hue[0.7], Hue[0.7], 1.67]]]\)

(Subscript[parametrization, 1])[R_, φ_, ψ_]
         = {{Cos[ψ], Sin[ψ], 0}, {-Sin[ψ], Cos[ψ], 
            0}, {0, 0, 1}}.
        ({0, R, 0} + {0, Cos[φ], Sin[φ] + ξ});

\(\(𝓅\_1 = 
      ParametricPlot3D[
        Evaluate[
          addColor[parametrization\_1[\(R\_*\), φ, ψ], 
            ψ]], {φ, \(5  π\)\/3, 2\ π}, {ψ, 
          \(-\(π\/2\)\), π\/2}, PlotPoints -> {ppq, ppl}]; \)\)

(Subscript[parametrization, 1  lr])[R_, r_, φ_, ψ_]
         = {{Cos[ψ], Sin[ψ], 0}, {-Sin[ψ], Cos[ψ], 
            0}, {0, 0, 1}}.
        ({0, R, 0} + r {0, 0, ξ} + 
            r {0, Cos[φ], Sin[φ]} + 
            (1 - r) {0, 1, 0});

\(\(𝓅\_\(1  l\) = 
      ParametricPlot3D[
        Evaluate[
          addColor\_l[
            parametrization\_\(1  lr\)[\(R\_*\), η, φ, 
              ψ]]], {φ, \(5  π\)\/3, 2\ π}, {ψ, 
          \(-\(π\/2\)\), \(-\(π\/2\)\) - dψ}, 
        PlotPoints -> {ppq, ppq}]; \)\)

\(\(𝓅\_\(1  r\) = 
      ParametricPlot3D[
        Evaluate[
          addColor\_r[
            parametrization\_\(1  lr\)[\(R\_*\), η, φ, 
              ψ]]], {φ, \(5  π\)\/3, 2\ π}, {ψ, 
          π\/2, π\/2 + dψ}, PlotPoints -> {ppq, ppq}]; 
    \)\)

(Subscript[parametrization, 2])[R_, φ_, ψ_]
         = {{Cos[ψ], Sin[ψ], 0}, {-Sin[ψ], Cos[ψ], 
            0}, {0, 0, 1}}.
        ({0, R, 0} + {0, Cos[φ], Sin[φ] + ξ} + {0, 
              2, 0});

\(\(𝓅\_2 = 
      ParametricPlot3D[
        Evaluate[
          addColor[parametrization\_2[\(R\_*\), φ, ψ], 
            ψ]], {φ, π, \(4  π\)\/3}, {ψ, 
          \(-\(π\/2\)\), π\/2}, PlotPoints -> {ppq, ppl}]; \)\)

(Subscript[parametrization, 2  lr])[R_, r_, φ_, ψ_]
         = {{Cos[ψ], Sin[ψ], 0}, {-Sin[ψ], Cos[ψ], 
            0}, {0, 0, 1}}.
        ({0, R, 0} + r {0, Cos[φ], Sin[φ]} + 
            r {0, 2, ξ} + (1 - r) {0, 1, 0});

\(\(𝓅\_\(2  l\) = 
      ParametricPlot3D[
        Evaluate[
          addColor\_l[
            parametrization\_\(2  lr\)[\(R\_*\), η, φ, 
              ψ]]], {φ, π, \(4  π\)\/3}, {ψ, 
          \(-\(π\/2\)\), \(-\(π\/2\)\) - dψ}, 
        PlotPoints -> {ppq, ppq}]; \)\)

\(\(𝓅\_\(2  r\) = 
      ParametricPlot3D[
        Evaluate[
          addColor\_r[
            parametrization\_\(2  lr\)[\(R\_*\), η, φ, 
              ψ]]], {φ, π, \(4  π\)\/3}, {ψ, 
          π\/2, π\/2 + dψ}, PlotPoints -> {ppq, ppq}]; 
    \)\)

(Subscript[parametrization, 3])[R_, φ_, ψ_]
         = {{Cos[ψ], Sin[ψ], 0}, {-Sin[ψ], Cos[ψ], 
            0}, {0, 0, 1}}.
        ({0, R, 0} + {0, Cos[φ], Sin[φ]} + {0, 1, 
              -2 ξ});

\(\(𝓅\_3 = 
      ParametricPlot3D[
        Evaluate[
          addColor[parametrization\_3[\(R\_*\), φ, ψ], 
            ψ]], {φ, π\/3, \(2\ π\)\/3}, {ψ, 
          \(-\(π\/2\)\), π\/2}, PlotPoints -> {ppq, ppl}]; \)\)

(Subscript[parametrization, 3  lr])[R_, r_, φ_, ψ_]
         = {{Cos[ψ], Sin[ψ], 0}, {-Sin[ψ], Cos[ψ], 
            0}, {0, 0, 1}}.
        ({0, R, 0} + r {0, Cos[φ], Sin[φ]} + 
            r {0, 1, -2 ξ} + (1 - r) {0, 1, 0});

\(\(𝓅\_\(3  l\) = 
      ParametricPlot3D[
        Evaluate[
          addColor\_l[
            parametrization\_\(3  lr\)[\(R\_*\), η, φ, 
              ψ]]], {φ, π\/3, \(2\ π\)\/3}, {ψ, 
          \(-\(π\/2\)\), \(-\(π\/2\)\) - dψ}, 
        PlotPoints -> {ppq, ppq}]; \)\)

\(\(𝓅\_\(3  r\) = 
      ParametricPlot3D[
        Evaluate[
          addColor\_r[
            parametrization\_\(3  lr\)[\(R\_*\), η, φ, 
              ψ]]], {φ, π\/3, \(2\ π\)\/3}, {ψ, 
          π\/2, π\/2 + dψ}, PlotPoints -> {ppq, ppq}]; 
    \)\)

wholeDial
       = {EdgeForm[], 
        First/@{(Subscript[𝓅, 1]), (Subscript[𝓅, 2]), (Subscript[𝓅, 3]), 
            (Subscript[𝓅, 1  l]), (Subscript[𝓅, 2  l]), (Subscript[𝓅, 3  l]), 
            (Subscript[𝓅, 1  r]), (Subscript[𝓅, 2  r]), 
            (Subscript[𝓅, 3  r])}};

Show[Graphics3D[wholeDial]];

The hand

The following gives the parametrizations for the sundial's hand.

radius[{l_, h_, a_}, s_] := 
      With[{x = 
            ((3 + 2 a l - (Sqrt[9 + 4 a l + 4 (a^2) (l^2)]))/(4. a))}, 
        ((h (Sqrt[s]) Exp[-a s] (l - s))/((Sqrt[x]) Exp[-a x] 
            (l - x)))]

𝒽 = 3.5; (Subscript[𝒷, max]) = (2/3); a = 1;

(Subscript[parametrization, 1])[r_, h_, φ_, s_] = {r, (r/(Sqrt[3])), h s} + 
        r {Cos[φ], Sin[φ], 0};

(Subscript[𝓅, 1]) = 
      ParametricPlot3D[
        Evaluate[
          (Subscript[parametrization, 1])[
            radius[{𝒽, (Subscript[𝒷, max]), a}, 𝒽 s], 
            𝒽, φ, s]], {φ, π, 
          π + (π/3)}, {s, 0, 1}, PlotPoints -> {ppq, ppl}];

(Subscript[parametrization, 2])[r_, h_, φ_, s_]
         = {-r, (r/(Sqrt[3])), h s} + 
        r {Cos[φ], Sin[φ], 0};

(Subscript[𝓅, 2]) = 
      ParametricPlot3D[
        Evaluate[
          (Subscript[parametrization, 2])[
            radius[{𝒽, (Subscript[𝒷, max]), a}, 𝒽 s], 
            𝒽, φ, s]], {φ, ((5  π)/3), 
          2 π}, {s, 0, 1}, PlotPoints -> {ppq, ppl}];

(Subscript[parametrization, 3])[r_, h_, φ_, s_]
         = { 0, -((2  r)/(Sqrt[3])), h s} + 
        r {Cos[φ], Sin[φ], 0};

(Subscript[𝓅, 3]) = 
      ParametricPlot3D[
        Evaluate[
          (Subscript[parametrization, 3])[
            radius[{𝒽, (Subscript[𝒷, max]), a}, 𝒽 s], 
            𝒽, φ, s]], {φ, (π/3), 
          ((2  π)/3)}, {s, 0, 1}, PlotPoints -> {ppq, ppl}];

handPolys = 
      Cases[{(Subscript[𝓅, 1]), (Subscript[𝓅, 2]), (Subscript[𝓅, 3])}, _Polygon, 
        ∞];

(Subscript[ℓ, 1]) = 
      ParametricPlot3D[
        Evaluate[
          (Subscript[parametrization, 1])[
            radius[{𝒽, (Subscript[𝒷, max]), a}, 𝒽 s], 
            𝒽, π, s]], {s, 0, 1}, PlotPoints -> 100];

(Subscript[ℓ, 2]) = 
      ParametricPlot3D[
        Evaluate[
          (Subscript[parametrization, 2])[
            radius[{𝒽, (Subscript[𝒷, max]), a}, 𝒽 s], 
            𝒽, ((5  π)/3), s]], {s, 0, 1}, 
        PlotPoints -> 100];

(Subscript[ℓ, 3]) = 
      ParametricPlot3D[
        Evaluate[
          (Subscript[parametrization, 3])[
            radius[{𝒽, (Subscript[𝒷, max]), a}, 𝒽 s], 
            𝒽, (π/3), s]], {s, 0, 1}, PlotPoints -> 100];

handEdges = 
      Map[{1.01, 1.01, 1} #&, 
        Cases[{(Subscript[ℓ, 1]), (Subscript[ℓ, 2]), (Subscript[ℓ, 3])}, _Line, 
          ∞], {-2}];

rotate[{x_, y_, z_}, α_] := 
      {{Cos[#], Sin[#], 0}, {-Sin[#], Cos[#], 0}, {0, 0, 1}}&[
        N[((z 2 π α)/𝒽)]];

𝓇 = N[{{(Sqrt[2]), 0, 0}, {0, 1, 1}, {0, -1, 1}}/(Sqrt[2])];

𝒻[p_] := 𝓇.(rotate[p, 1].p) + {0, -0.3, 0};

wholeHand
       = {{EdgeForm[], SurfaceColor[Hue[0.11], Hue[0.8], 1.78], 
          Map[𝒻, handPolys, {-2}]}, {Thickness[0.003], 
          Map[𝒻, handEdges, {-2}]}};

Show[Graphics3D[wholeHand], PlotRange -> All, 
      Axes -> True];

Projecting the hour lines onto the dial

We project the shadow of the hand onto the inside face of the dial.

ppl = 20; ppq = 8;

\(\(zBottom = 
      \(Simplify[parametrization\_1[\(R\_*\), \(5  π\)\/3, 0]]
          \)[[\(-1\)]]; \)\)

\(\(bottomRadius = 
      N[\(\@\(#.#\)&\)[{0, 0, zBottom} - 
            N[parametrization\_1[\(R\_*\), \(5  π\)\/3, 0]]]]; \)\)

\(\(𝓅\_1 = 
      ParametricPlot3D[
        Evaluate[parametrization\_1[\(R\_*\), φ, ψ]], {
          φ, 2 \( 5  π\)\/3, 2\ π}, {ψ, 
          \(-\(π\/2\)\), π\/2}, PlotPoints -> {ppq, ppl}]; \)\)

\(𝒸\_1[α_, R_, φ_] = 
      ArcSin[\((
            \(-\((Cos[α]\ 
                    \((4\ \@3\ R + 4\ \@3\ Cos[φ] + 
                        12\ R\ Sin[φ] + 
                        12\ Cos[φ]\ Sin[φ])\)\ 
                    Sin[α])\)\) - 
              √\((
                  Cos[α]\^2\ 
                      \((4\ \@3\ R + 4\ \@3\ Cos[φ] + 
                            12\ R\ Sin[φ] + 
                            12\ Cos[φ]\ Sin[φ])\)^2\ 
                      Sin[α]\^2 + 
                    4\ \((\(-1\) + 3\ R\^2 + 6\ R\ Cos[φ] + 
                          3\ Cos[φ]\^2 - 2\ \@3\ Sin[φ] - 
                          3\ Sin[φ]\^2)\)\ Sin[α]\^2\ 
                      \((12\ R\^2\ Cos[α]\^2 + 
                          24\ R\ Cos[φ]\ Cos[α]\^2 + 
                          12\ Cos[φ]\^2\ Cos[α]\^2 + 
                          3\ R\^2\ Sin[α]\^2 + 
                          6\ R\ Cos[φ]\ Sin[α]\^2 + 
                          3\ Cos[φ]\^2\ Sin[α]\^2)\))\))\)/
          \((2\ \((
                12\ R\^2\ Cos[α]\^2 + 
                  24\ R\ Cos[φ]\ Cos[α]\^2 + 
                  12\ Cos[φ]\^2\ Cos[α]\^2 + 
                  3\ R\^2\ Sin[α]\^2 + 
                  6\ R\ Cos[φ]\ Sin[α]\^2 + 
                  3\ Cos[φ]\^2\ Sin[α]\^2)\))\)]; 


    𝒸\_2[α_, R_, φ_] = 
      ArcSin[\((
            \(-\((Cos[α]\ 
                    \((4\ \@3\ R + 4\ \@3\ Cos[φ] + 
                        12\ R\ Sin[φ] + 
                        12\ Cos[φ]\ Sin[φ])\)\ 
                    Sin[α])\)\) + 
              √\((
                  Cos[α]\^2\ 
                      \((4\ \@3\ R + 4\ \@3\ Cos[φ] + 
                            12\ R\ Sin[φ] + 
                            12\ Cos[φ]\ Sin[φ])\)^2\ 
                      Sin[α]\^2 + 
                    4\ \((\(-1\) + 3\ R\^2 + 6\ R\ Cos[φ] + 
                          3\ Cos[φ]\^2 - 2\ \@3\ Sin[φ] - 
                          3\ Sin[φ]\^2)\)\ Sin[α]\^2\ 
                      \((12\ R\^2\ Cos[α]\^2 + 
                          24\ R\ Cos[φ]\ Cos[α]\^2 + 
                          12\ Cos[φ]\^2\ Cos[α]\^2 + 
                          3\ R\^2\ Sin[α]\^2 + 
                          6\ R\ Cos[φ]\ Sin[α]\^2 + 
                          3\ Cos[φ]\^2\ Sin[α]\^2)\))\))\)/
          \((2\ \((
                12\ R\^2\ Cos[α]\^2 + 
                  24\ R\ Cos[φ]\ Cos[α]\^2 + 
                  12\ Cos[φ]\^2\ Cos[α]\^2 + 
                  3\ R\^2\ Sin[α]\^2 + 
                  6\ R\ Cos[φ]\ Sin[α]\^2 + 
                  3\ Cos[φ]\^2\ Sin[α]\^2)\))\)]; \)

hourLine[α_, R_, φ_] = 
      (Subscript[parametrization, 1])[R, φ, 
        If[α < 0, Evaluate[(Subscript[𝒸, 1])[α, R, φ]], 
          Evaluate[(Subscript[𝒸, 2])[α, R, φ]]]];

\(\(hourLines
       = {GrayLevel[0], Thickness[0.005], 
        Cases[Show[
            Table[ParametricPlot3D[
                Evaluate[
                  hourLine[N[α], 0.995\ \(R\_*\), φ]], {
                  φ, \(5  π\)\/3, 2\ π}, 
                PlotPoints -> 100, PlotRange -> All, 
                AxesLabel -> {"x", "y", "z"}, 
                DisplayFunction -> Identity], {α, 
                \(-75\)\ °, 75\ °, 15\ °}], 
            DisplayFunction -> $DisplayFunction], _Line, ∞]}; 
    \)\)

The bottom

At the bottom we use a triangle subdivsion due Radin.

bottomPolys = 
      MapIndexed[{
            Hue[Mod[#2[[1]], 5] 0.8], 
            Polygon[Map[{-1, -1} (#1 - {0, 2})&, #1, {
                  -2}]]}&, 
        Nest[Flatten[
              Apply[({{#5, #2, #4}, {#4, #7, #5}, {#7, #1, #3}, {#6, #4, #1}, {#6, #4, #7}}&)@@{#1, #2, #3, (#1 + #2)/2,  
                      (3 #2 + 2 #3)/5,  (5 #1 + #2 + 4 #3)/10,  
                      (#2 + 4 #3)/5}&, #1, {1}], 
              1]&, {{{0., 0.}, {0., 2.}, {-1., 0.}}}, 3]];

We generate a three-dimensional version of this pattern by reflection. We fit the boundaries into the circular boundary of the  dial.

to3DCurved[{ψ1_, ψ2_}, p_] := 
      Module[{φ = ArcTan@@p, r = (Sqrt[p.p]), ψMin, ψMax, 
          newφ, newr}, {ψMin, ψMax} = 
          Sort[N[{ψ1, ψ2}]]; 
        newφ = 
          N[ψ1 + 
              (((ψ2 - ψ1) 
                  (φ - ArcTan[2]))/((π/2) - ArcTan[2]))]; 
        newr = N[((bottomRadius r)/(2 Csc[φ]))]; {
          newr Cos[newφ], newr Sin[newφ], zBottom}]

to3DCurved[{ψ1_, ψ2_}, {0., 0.}] := {0., 0., zBottom}

to3DCurved[{ψ1_, ψ2_}] := {EdgeForm[], 
        bottomPolys
           /. {Polygon[l_] :> 
              Polygon[to3DCurved[{ψ1, ψ2}, #1]&/@l], 
            h_Hue :> SurfaceColor[h, h, 2.7]}}

\(\(hourLineTab = 
      N[Join[{{\(-bottomRadius\), 0, zBottom}}, 
          N[Table[hourLine[N[α], 0.99\ \(R\_*\), \(5  π\)\/3], {
                α, \(-75\)\ °, 75\ °, 
                15\ °}]], {{bottomRadius, 0, zBottom}}]]; \)\)

bottomPolys = 
      to3DCurved/@
        Sort/@Partition[(N[ArcTan@@Drop[#1, -1]]&)/@hourLineTab, 2, 
            1];

Here are the hour line shadows.

bottomHourLines
       = {Thickness[0.005], (Line[{{0, 0, zBottom}, #1}]&)/@hourLineTab};

Show[Graphics3D[{bottomPolys, bottomHourLines}], 
      PlotRange -> All];

The Roman numbers

ℛ[n] represents the Roman number n as a set of line primitives.

\(ℛ[1] = {Line[{{0, \(-1\)}, {0, 1}}], 
        Line[{{\(-0.2\), \(-1\)}, {0.2, \(-1\)}}], 
        Line[{{\(-0.2\), 1}, {0.2, 1}}]}; 


    ℛ[2] = 
      Flatten[{Map[{0.3, 0} + #&, ℛ[1], {\(-2\)}], 
          Map[{\(-0.3\), 0} + #&, ℛ[1], {\(-2\)}]}]; 


    ℛ[3] = 
      Flatten[{ℛ[1], 
          Map[{0.6, 0} + #&, ℛ[1], {\(-2\)}], 
          Map[{\(-0.6\), 0} + #&, ℛ[1], {\(-2\)}]}]; 


    ℛ[5] = {Line[{{\(-0.5\), 1}, {0, \(-1\)}}], 
        Line[{{0, \(-1\)}, {0.5, 1}}], Line[{{\(-0.7\), 1}, {\(-0.3\), 1}}], 
        Line[{{0.7, 1}, {0.3, 1}}]}; 


    ℛ[4] = 
      Flatten[{Map[{0.4, 0} + #&, ℛ[5], {\(-2\)}], 
          Map[{\(-0.7\), 0} + #&, ℛ[1], {\(-2\)}]}]; 


    ℛ[6] = 
      Flatten[{Map[{0.4, 0} + #&, ℛ[1], {\(-2\)}], 
          Map[{\(-0.7\), 0} + #&, ℛ[5], {\(-2\)}]}]; 


    ℛ[7] = 
      Flatten[{Map[{0.4, 0} + #&, ℛ[1], {\(-2\)}], 
          Map[{1.0, 0} + #&, ℛ[1], {\(-2\)}], 
          Map[{\(-0.7\), 0} + #&, ℛ[5], {\(-2\)}]}]; 


    ℛ[8] = 
      Flatten[{Map[{0.4, 0} + #&, ℛ[1], {\(-2\)}], 
          Map[{1.0, 0} + #&, ℛ[1], {\(-2\)}], 
          Map[{1.6, 0} + #&, ℛ[1], {\(-2\)}], 
          Map[{\(-0.7\), 0} + #&, ℛ[5], {\(-2\)}]}]; 


    ℛ[10] = {Line[{{\(-0.5\), 1}, {0.5, \(-1\)}}], 
        Line[{{0.5, 1}, {\(-0.5\), \(-1\)}}], 
        Line[{{\(-0.7\), 1}, {\(-0.3\), 1}}], Line[{{0.7, 1}, {0.3, 1}}], 
        Line[{{\(-0.7\), \(-1\)}, {\(-0.3\), \(-1\)}}], 
        Line[{{0.7, \(-1\)}, {0.3, \(-1\)}}]}; 


    ℛ[9] = 
      Flatten[{Map[{0.4, 0} + #&, ℛ[10], {\(-2\)}], 
          Map[{\(-0.7\), 0} + #&, ℛ[1], {\(-2\)}]}]; 


    ℛ[11] = 
      Flatten[{Map[{0.4, 0} + #&, ℛ[1], {\(-2\)}], 
          Map[{\(-0.7\), 0} + #&, ℛ[10], {\(-2\)}]}]; 


    ℛ[12] = 
      Flatten[{Map[{0.4, 0} + #&, ℛ[1], {\(-2\)}], 
          Map[{1.0, 0} + #&, ℛ[1], {\(-2\)}], 
          Map[{\(-0.7\), 0} + #&, ℛ[10], {\(-2\)}]}]; \)

Do[
      mp =  (Min[#] + Max[#])/2&[
          First/@Level[ℛ[i], {-2}]]; 
      ℛ[i] = 
        Map[# - {mp, 0}&, ℛ[i], {-2}], {i, 1, 12}];

Show[
      Graphics[Table[
          Map[# + {4  i, 0}&, ℛ[i], {-2}], {i, 1, 12}]], 
      AspectRatio -> 0.14];

Projecting the Roman numbers

We project the two-dimensional Roman numbers onto the curved three-dimensional dial surface.

refine[l_, n_] := 
      l /. Line[{a_, b_}] :> 
          Line[Table[a + ((i )/n) (b - a), {i, 0, n}]]

Do[
      ℛRefined[i] = refine[ℛ[i], 30], {i, 1, 
        12}];

\(mℛ[α_] := 
      N[N[{\(17  π\)\/9, 
            If[α < 0, 
              Evaluate[
                𝒸\_1[α, 0.99\ \(R\_*\), \(17  π\)\/9]], 
              Evaluate[
                𝒸\_2[α, 0.99\ \(R\_*\), \(17  π\)\/9]]]}]]
        \)

\(latitudeToHour[α_] := Mod[12 + α\/\(15\ °\), 12]; 
    
latitudeToHour[0] = 12; \)

\(\(romanHours
       = {Thickness[0.003], GrayLevel[1], 
        Table[Apply[parametrization\_1[0.985\ \(R\_*\), ##]&, 
            Map[mℛ[α] + {0.12, 0.04}\ Reverse[#1]&, 
              ℛRefined[latitudeToHour[α]], {\(-2\)}], {
              \(-2\)}], {α, \(-75\)\ °, 75\ °, 
            15\ °}]}; \)\)

Show[
      Graphics3D[romanHours[[-1]]], 
      PlotRange -> All, Boxed -> False];

Putting all the pieces together

Putting all of the pieces together we get the following sundial.

Show[
      Graphics3D[{wholeDial, wholeHand, hourLines, bottomPolys, 
          bottomHourLines, romanHours}], PlotRange -> All, 
      Boxed -> False];

Calculating an advanced sundial—correct to the minute

Background from celestial mechanics

The eccentricity of the earth's orbit around the sun results in a time-dependent angular velocity. That, together with the tilt of the earth's axis, means that length of a sun-day depends on the time of the year. A sun-day may be either a few seconds longer or shorter than 24 hours. So the position of the shadow of the hand's end depends on the time of the year, too. After a few weeks the sundial constructed above might be wrong up to 16 minutes.

Calculating the orbit of the sun

Find the position of the sun relative to the earth by solving the Kepler problem in terms of the polar coordinates {ψ(t),  r(t)} with eccentricity κ.

\(TraditionalForm
      \`\(ψ\^′\)(t) = 
        \((κ\ \(cos(ψ(t))\) + 1)\)\^2
            \/\@\((1 + κ\^2)\)\^3\)        \(TraditionalForm
      \`r(t) = \(1 - κ\^2\)\/\(1 + κ\ \(cos(ψ(t))\)\)\)

We can't get a closed form solution for ψ(t), so we calculate a series solution for it.

(Subscript[ψ, ser])=t+(∑_(j = 1)^2) (Subscript[ψ, j])[t] (κ^j)+O[κ]^3

(
        ψ'[t] == 
          ((((1 + κ Cos[ψ[t]]))^2)/((Sqrt[(1 - (κ^2))])^3)))
         /. {ψ[t]  -> (Subscript[ψ, ser]), 
        ψ'[t]  -> (∂_t) (Subscript[ψ, ser])}

LogicalExpand[%]

DSolve[
      Apply[List, % ∧ (Subscript[ψ, 1])[0] == 0 ∧ 
          (Subscript[ψ, 2])[0] == 0], {(Subscript[ψ, 1])[t], (Subscript[ψ, 2])[t]}, t]

So we have the following series solution (Subscript[ψ, α]) to the second order in κ.

(Subscript[ψ, α]) = Normal[(Subscript[ψ, ser]) /. First[%]]

We compare this with the numerical solution (Subscript[ψ, 𝒩]) using the actual value of \(κ\&~\) for the earth's orbit:

\(κ\&~ = 0.0016722; 

    ψ\_𝒩 = 
      ψ[t] /. 
        \(NDSolve[{
              \(ψ'\)[t] == 
                \((1 + κ\&~\ Cos[ψ[t]])\)\^2
                    \/\@\((1 - κ\&~\^2)\)\^3, ψ[0] == 0}, 

            ψ[t], {t, 0, 2  π}, \ 
				
            PrecisionGoal -> 14, WorkingPrecision -> 26]
            \)[[1]]\)

\(\(Plot[
      Evaluate[\((ψ\_α - ψ\_𝒩)\) /. 
          κ -> κ\&~], {t, 0, 2  π}]; \)\)

Here is the maximal deviation of the series solution.

\(FindMinimum[
      Evaluate[\((ψ\_α - ψ\_𝒩)\) /. 
          κ -> κ\&~], {t, \(3\/2\) π}]\)

The equation of time

A sun-day deviates from 24 hours by only a few second per day , but over weeks this adds up to a quarter of an hour.

\(ψ[t_] = 
      t + 2  κ\ \ Sin[t] + \(5\/4\) κ\^2\ \ Sin[2  t]; 

    δ[t_] = ψ[t] - t; 

    h[t_] := Sin[ψ[t] - α]\ Sin[ϵ]; 

    r[t_] := \@\(1 - h[t]\^2\); 

    μ[t_] := 
      ArcTan[\(\(-Sin[δ[t]]\) + 
            2 \( Sin[ϵ\/\(2\ \)]\^2\) Sin[ψ[t] - α] 
              Cos[t - α]\)\/\(\(+Cos[δ[t]]\) + 
            2 \( Sin[ϵ\/\(2\ \)]\^2\) Sin[ψ[t] - α] 
              Sin[t - α]\)]; \)

We will work with this numerical data.

κ = 0.016722;    \(ϵ = 23.45  °; \ \ \ \ 
     (*\ ecliptic - earth\ axis\ angle\ *) 
α = 78.5  °; 
    \ \ \ \ \ \  (*\ spring\ point\ *) 
Θ = 40  °; 
    \ \ \ \ \ \ \ \ \ \ \  (*\ altitude\ *) 
ϕ = 15  °; 
    \ \ \ \ \ \ \ \ \ \ \ \  (*\ longitude\ *) \)

The equation of time describes the deviation of "naive sundial time" from real time.

Plot[Evaluate[((1440 μ[t])/(2 π))], {t, 0, 2 π}, 
      Axes -> {True, False}, Frame -> True, 
      FrameLabel -> {None, "min"}, 
      FrameTicks 
        -> {
          Join[(Append[#1, {0, 0}]&)/@
              Transpose[{
                  Table[i (( 2 π)/12) - ((2 π)/(2 12 )), {i, 
                      12}], {"Jan", "Feb", "Mar", "Apr", 
                    "May", "Jun", "Jul", "Aug", "Sep", 
                    "Oct", "Nov", "Dec"}}], 
            {#, ""}&/@Table[((i 2 π)/12), {i, 12}]], 
          Automatic, None, None}, 
      PlotStyle -> {Hue[0], Thickness[0.01]}];

The sun's path through the year

This describes the path of the sun through the year and on any given day.

\(x\_1[s_, t_] := \(+r[t]\)\ Sin[s + μ[t]]; 

    x\_2[s_, t_] := 
      \(+r[t]\)\ Cos[s + μ[t]] Sin[Θ] + 
        h[t]\ Cos[Θ]; 

    x\_3[s_, t_] := 
      \(-r[t]\)\ Cos[s + μ[t]] Cos[Θ] + 
        h[t]\ Sin[Θ]; \)

X[s_, t_] = {(Subscript[x, 1])[s, t], (Subscript[x, 2])[s, t], (Subscript[x, 3])[s, t]};

Needs["Graphics`Shapes`"]

Needs["Graphics`Colors`"]

blueSphere
       = {EdgeForm[], 
        SurfaceColor[RGBColor[0, 0, 1], RGBColor[0, 0, 1], 1.78], 
        Sphere[0.99, 20, 20]};

Show[
      GraphicsArray[{
          Show[{Graphics3D[{EdgeForm[], 
                  SurfaceColor[RGBColor[0, 0, 1], RGBColor[0, 0, 1], 1.78], 
                  Sphere[0.99, 20, 20]}], 
              Table[ParametricPlot3D[
                  Evaluate[Append[X[s, t], Yellow]], {t, 0, 2 π}], {s, 
                  0, 2 π, ((2 π)/24)}]}, Boxed -> False, 
            DisplayFunction -> Identity], 
          Show[{Graphics3D[{EdgeForm[], 
                  SurfaceColor[RGBColor[0, 0, 1], RGBColor[0, 0, 1], 1.78], 
                  Sphere[0.99, 20, 20]}], 
              Table[ParametricPlot3D[
                  Evaluate[Append[X[s, t], Yellow]], {s, 0, 2 π}], {t, 
                  0, 2 π, ((2 π)/12)}]}, Boxed -> False, 
            DisplayFunction -> Identity]}]];

This is the path of the sun as seen form a fixed point on the earth.

ParametricPlot[{μ[t], ArcSin[h[t]]}/°, {t, 0, 2  π}
         ];

Analemmas of the dial—from sunrise to sunset

This is the position of the shadow on the earth.

localhour[S_] := 
      (π/12) (S + ((ϕ - 15)/15))

shadowPosition[S_, t_] := {
        -(Subscript[x, 1])[localhour[S], t]/(Subscript[x, 3])[localhour[S], t], 
        -(Subscript[x, 2])[localhour[S], t]/(Subscript[x, 3])[localhour[S], t]}

These are the times of the sunrise and sunset over the year.

\(sunriseAndSunset[S_] := 
      \(sunriseAndSunset[S] = 
        \((If[#1 === {}, {0, 2\ π}, {1.02, 0.98}\ #1]&)\)[
          \((t /. FindRoot[Evaluate[x\_3[localhour[S], t] == 0], {t, #1}]&)
                \)/@\((2\ π\ 
                  Flatten[
                    Position[
                      Apply[Times, 
                        Partition[
                          Table[Evaluate[x\_3[localhour[S], t]], {t, 0, 
                              2\ π, 2\ π/1000}], 2, 1], {1}], _?
                        Negative]])\)/1000]\); 

    Do[sunriseAndSunset[S], {S, 8, 18}]; \)

This is the shadow of the sun at the hours.

\(\(hourLines = 
      Table[Table[
          If[
				
            \(sunriseAndSunset[S]\)[[1]] 
              ≤ t ≤ 
              \(sunriseAndSunset[S]\)[[2
                ]], shadowPosition[S, t], {}], {S, 8, 
            18}], {t, \(2\ π\)\/\(2\ 12\), 2\ π, \(2\ π\)\/12}]; 
    \)\)

These are the analemmas and the hourlines.

\(\(Show[
      Table[\(({Graphics[{Thickness[0.001], Hue[0.78], 
                    Line/@DeleteCases[hourLines, {}, {2}]}], #1, 
                Graphics[
                  Text[FontForm[S, {"Helvetica", 10}], 
                    1.06\ #1[[1, 1, 1, 1, 
                        \(\((Position[#1, Max[#1]]&)\)[
                            \((#1.#1&)\)/@#1[[1, 1, 1, 1
                                ]]]
                            \)[[1, 1
                          ]]]]]]}&)\)[
          ParametricPlot[Evaluate[shadowPosition[S, t]], 
            Evaluate[Prepend[sunriseAndSunset[S], t]], 
				
            DisplayFunction -> Identity, PlotPoints -> 100]], {S, 
          7, 19}], DisplayFunction -> $DisplayFunction, 
      PlotRange -> {{\(-4\), 4}, {\(-1\), 4}}, Axes -> False, 
      Frame -> True, FrameTicks -> False]; \)\)

Web hyperlinks to sundials

*Sundial1

*Sundial2

References

C. Blatter.  Elem. Math. 49, 155 (1994).
A. W. Fuller.  Math. Gazette XLI, 9 (1957).
R. R. J. Rohr. Die Sonnenuhr, Callwey, München, 1982.
E. Salkowski.  Sitzungsber. Berlin. Math. Ges. 10, 23 (1910).
W. Schaub.  Die Sterne 203 (1957).
G. Scheffers.  Sitzungsber. Berlin. Math. Ges. 8, 122 (1909).
H. Schumacher.  Sonnenuhren, Callwey, München, 1973.
A. Zenkert.  Faszination Sonnenuhr, Verlag Technik, Berlin, 1984.

Initialization

Off[General::"spell1"]; Off[General::"spell"];