Setup and Definitions

We start by establishing the geometry, location and orientation of the gnomon (shadow-casting stick) and face of the sundial.

import sympy as sp
from sympy import sin, cos, tan
from sympy.abc import alpha, psi, sigma, iota, delta
from analemma.algebra import frame, render, util

Fixed Stars and Earth Frames

Define an orthonormal frame (set of basis vectors) \(e_1, e_2, e_3\) that is fixed relative to the "fixed stars", in analemma.algebra.frame.base. Then, let the tilt of the earth's plane (axis) of rotation be \(\alpha\) and measure the earth's rotation by \(\psi\). The Earth frame \(f_1, f_2, f_3\) is then given by analemma.algebra.frame.planet.

.

e1, e2, e3 = frame.base("e")
f1, f2, f3 = frame.planet()
render.frame("f", (f1,f2,f3))

\[\displaystyle \begin{align} f_0 &= \cos{\left (\alpha \right )} \cos{\left (\psi \right )} \boldsymbol{e}_{1} + \sin{\left (\psi \right )} \boldsymbol{e}_{2} + \sin{\left (\alpha \right )} \cos{\left (\psi \right )} \boldsymbol{e}_{3}\nonumber \\ f_1 &= - \sin{\left (\psi \right )} \cos{\left (\alpha \right )} \boldsymbol{e}_{1} + \cos{\left (\psi \right )} \boldsymbol{e}_{2} - \sin{\left (\alpha \right )} \sin{\left (\psi \right )} \boldsymbol{e}_{3}\nonumber \\ f_2 &= - \sin{\left (\alpha \right )} \boldsymbol{e}_{1} + \cos{\left (\alpha \right )} \boldsymbol{e}_{3}\nonumber \end{align}\]

The equatorial plane should only depend on the tilt of the Earth's axis of spin \(\alpha\), not the angle by which it has rotated relative to the fixed stars \(\psi\).

render.expression(r"f_1 \wedge f_2", (f1^f2))

\[\displaystyle \begin{equation} f_1 \wedge f_2 = \cos{\left (\alpha \right )} \boldsymbol{e}_{1}\wedge \boldsymbol{e}_{2} - \sin{\left (\alpha \right )} \boldsymbol{e}_{2}\wedge \boldsymbol{e}_{3} \nonumber \end{equation} \]

Note: the wedge operation "\(\wedge\)" defines a bivector, which encodes the orientation of a plane in space (but also has a magnitude, just as a vector has a direction and a magnitude). This is a concept from Geometric Algebra which is used throughout the analemma package. A good place to learn more is https://bivector.net/.

Surface Frame

Define an orthonormal frame embedded in the Earth's surface, with \(n_1\) pointing South, \(n_2\) pointing East and \(n_3\) pointing up. Note that \(\theta\) is not latitude, but \(90^\circ\) minus latitude.

.

n1, n2, n3 = frame.surface()
render.frame("n", (n1,n2,n3))

\[\displaystyle \begin{align} n_0 &= \left ( \sin{\left (\alpha \right )} \sin{\left (\theta \right )} + \cos{\left (\alpha \right )} \cos{\left (\psi \right )} \cos{\left (\theta \right )}\right ) \boldsymbol{e}_{1} + \sin{\left (\psi \right )} \cos{\left (\theta \right )} \boldsymbol{e}_{2} + \left ( \sin{\left (\alpha \right )} \cos{\left (\psi \right )} \cos{\left (\theta \right )} - \sin{\left (\theta \right )} \cos{\left (\alpha \right )}\right ) \boldsymbol{e}_{3}\nonumber \\ n_1 &= - \sin{\left (\psi \right )} \cos{\left (\alpha \right )} \boldsymbol{e}_{1} + \cos{\left (\psi \right )} \boldsymbol{e}_{2} - \sin{\left (\alpha \right )} \sin{\left (\psi \right )} \boldsymbol{e}_{3}\nonumber \\ n_2 &= \left ( - \sin{\left (\alpha \right )} \cos{\left (\theta \right )} + \sin{\left (\theta \right )} \cos{\left (\alpha \right )} \cos{\left (\psi \right )}\right ) \boldsymbol{e}_{1} + \sin{\left (\psi \right )} \sin{\left (\theta \right )} \boldsymbol{e}_{2} + \left ( \sin{\left (\alpha \right )} \sin{\left (\theta \right )} \cos{\left (\psi \right )} + \cos{\left (\alpha \right )} \cos{\left (\theta \right )}\right ) \boldsymbol{e}_{3}\nonumber \end{align}\]

Orbit Rotor and Meridian Plane

Although we will work with the Earth for concreteness, the derivation applies to any planet. The angle \(\sigma\) measures the progress of the planet around its orbit, with \(\sigma = 0\) when the planet lies along the \(e_1\) axis. Therefore, a vector from the origin to the planet models a ray of light from the star - a sun ray.

s = frame.sunray()
render.expression("s", s)

\[\displaystyle \begin{equation} s = \cos{\left (\sigma \right )} \boldsymbol{e}_{1} + \sin{\left (\sigma \right )} \boldsymbol{e}_{2} \nonumber \end{equation} \]

The meridian plane contains a line of longitude and is defined by \(M = n_1 \wedge n_3\).

M = frame.meridian_plane()
render.expression("M", M)

\[\displaystyle \begin{equation} M = \sin{\left (\alpha \right )} \sin{\left (\psi \right )} \boldsymbol{e}_{1}\wedge \boldsymbol{e}_{2} + \cos{\left (\psi \right )} \boldsymbol{e}_{1}\wedge \boldsymbol{e}_{3} + \sin{\left (\psi \right )} \cos{\left (\alpha \right )} \boldsymbol{e}_{2}\wedge \boldsymbol{e}_{3} \nonumber \end{equation} \]

The noon line is the intersection of the sun ray \(s\) and the meridian plane \(M\), which occurs where \(s \wedge M\) (a 3-d object) vanishes.

render.multivector((s^M).trigsimp())

\[\displaystyle \left ( \sin{\left (\psi \right )} \cos{\left (\alpha \right )} \cos{\left (\sigma \right )} - \sin{\left (\sigma \right )} \cos{\left (\psi \right )}\right ) \boldsymbol{e}_{1}\wedge \boldsymbol{e}_{2}\wedge \boldsymbol{e}_{3}\]

Solving for \(\psi\) gives the angle of Earth's rotation corresponding to noon.

coeff = (s^M).trigsimp().get_coefs(3)[0]
soln = sp.solve(coeff.subs(sin(psi), tan(psi)*cos(psi)), tan(psi))[0]

assert soln.equals(tan(sigma)/cos(alpha))
render.expression( r"\tan(\psi)", soln )

\[\displaystyle \begin{equation} \tan(\psi) = \frac{\tan{\left (\sigma \right )}}{\cos{\left (\alpha \right )}} \nonumber \end{equation} \]

Dial Face and Gnomon

Now we can define the orientation of the sundial. The unit vector \(m_3\) points upward perpendicular to the dial face, while \(m_1\) and \(m_2\) point North and West respectively when the face's inclination \(i\) and declination \(d\) are zero.

.

m1, m2, m3 = frame.dial()
render.frame("m", (m1,m2,m3))

\[\displaystyle \begin{align} m_0 &= \cos{\left (d \right )} \cos{\left (i \right )} \boldsymbol{n}_{1} + \sin{\left (d \right )} \cos{\left (i \right )} \boldsymbol{n}_{2} + \sin{\left (i \right )} \boldsymbol{n}_{3}\nonumber \\ m_1 &= - \sin{\left (d \right )} \boldsymbol{n}_{1} + \cos{\left (d \right )} \boldsymbol{n}_{2}\nonumber \\ m_2 &= - \sin{\left (i \right )} \cos{\left (d \right )} \boldsymbol{n}_{1} - \sin{\left (d \right )} \sin{\left (i \right )} \boldsymbol{n}_{2} + \cos{\left (i \right )} \boldsymbol{n}_{3}\nonumber \end{align}\]

Given this frame, we can form the unit bivector \(G_n\) that encodes the dial face directly.

Gn = frame.dialface()
render.expression( "G_n", Gn)

\[\displaystyle \begin{equation} G_n = \cos{\left (i \right )} \boldsymbol{n}_{1}\wedge \boldsymbol{n}_{2} + \sin{\left (d \right )} \sin{\left (i \right )} \boldsymbol{n}_{1}\wedge \boldsymbol{n}_{3} - \sin{\left (i \right )} \cos{\left (d \right )} \boldsymbol{n}_{2}\wedge \boldsymbol{n}_{3} \nonumber \end{equation} \]

The gnomon \(g_n\) expressed relative to the planet's surface frame is given by:

.

gn = frame.gnomon("n", zero_decl=False)

# extra manipulation to display exactly as in paper
render.expression("g_n", sp.collect(sp.trigsimp(gn.obj), -sin(iota)))

\[\displaystyle \begin{equation} g_n = - \sin{\left (\iota \right )} \left(\sin{\left (\delta \right )} \boldsymbol{n}_{2} + \cos{\left (\delta \right )} \boldsymbol{n}_{1}\right) + \cos{\left (\iota \right )} \boldsymbol{n}_{3} \nonumber \end{equation} \]

Projected onto the fixed-stars basis, the gnomon is

nn1, nn2, nn3 = frame.base("n")
g = util.project_vector(gn, target_frame=(nn1, nn2, nn3), render_frame=(n1, n2, n3))
g = g.trigsimp()

render.align("g", g)

\[\displaystyle \begin{align} g & = - \sin{\left (\alpha \right )} \sin{\left (\iota \right )} \sin{\left (\theta \right )} \cos{\left (\delta \right )} - \sin{\left (\alpha \right )} \cos{\left (\iota \right )} \cos{\left (\theta \right )} + \sin{\left (\delta \right )} \sin{\left (\iota \right )} \sin{\left (\psi \right )} \cos{\left (\alpha \right )} - \sin{\left (\iota \right )} \cos{\left (\alpha \right )} \cos{\left (\delta \right )} \cos{\left (\psi \right )} \cos{\left (\theta \right )} + \sin{\left (\theta \right )} \cos{\left (\alpha \right )} \cos{\left (\iota \right )} \cos{\left (\psi \right )} e_0 \nonumber \\ &- \sin{\left (\delta \right )} \sin{\left (\iota \right )} \cos{\left (\psi \right )} - \sin{\left (\iota \right )} \sin{\left (\psi \right )} \cos{\left (\delta \right )} \cos{\left (\theta \right )} + \sin{\left (\psi \right )} \sin{\left (\theta \right )} \cos{\left (\iota \right )} e_1 \nonumber \\ &\sin{\left (\alpha \right )} \sin{\left (\delta \right )} \sin{\left (\iota \right )} \sin{\left (\psi \right )} - \sin{\left (\alpha \right )} \sin{\left (\iota \right )} \cos{\left (\delta \right )} \cos{\left (\psi \right )} \cos{\left (\theta \right )} + \sin{\left (\alpha \right )} \sin{\left (\theta \right )} \cos{\left (\iota \right )} \cos{\left (\psi \right )} + \sin{\left (\iota \right )} \sin{\left (\theta \right )} \cos{\left (\alpha \right )} \cos{\left (\delta \right )} + \cos{\left (\alpha \right )} \cos{\left (\iota \right )} \cos{\left (\theta \right )} e_2& \nonumber \end{align}\]

The gnomon lies in the meridian plane when the following trivector vanishes:

M_wedge_g = M^g
assert M_wedge_g.obj.trigsimp().equals((sin(delta)*sin(iota)*e1^e2^e3).obj)
render.multivector( M_wedge_g )

\[\displaystyle \sin{\left (\delta \right )} \sin{\left (\iota \right )} \boldsymbol{e}_{1}\wedge \boldsymbol{e}_{2}\wedge \boldsymbol{e}_{3}\]

Next, we calculate the hour angle, measured between the face of the dial and the plane containing the sun ray and the gnomon