Time AND date sundial

Photo taken by me at 18:00 BST May 20, 2018

I was out and about in Cambridge and I saw a modern sundial on the side of a building in Tennis Court road.  It occurred to me then that I should be able to build a sundial which tells you both the time - in GMT - and the date.  (Although you do need to know whether it is before or after the summer solstice!).  In fact all you need is to know one out of compass bearing, time, date and you should be able to work out the other two!

To run it you need to do the following

1. copy the text into sundial.py and chmod +x sundial.py
2. install pre-requisite packages: sudo apt install python-numpy python-matplotlib
3. run it: ./sundial.py

The result is a printout like the one shown.  If you want to adapt the picture for your locale just edit the parameters passed to plot_fixed_lat_long().

When I printed out my first sundial I was surprised to discover that the trajectory of the shadow is straight on the equinoxes.  After a day of pondering I realized why: these are the only two days of the year when the Earth's axis makes a right angle with the light coming from the Sun.  So as the Earth rotates the shadow scores out a line of constant latitude!  Anyway,

Here's my code

#!/usr/bin/python -i
"""
 Coordinate system 1:

   zero:   centre of earth on winter solstice
   z-axis: towards sun on winter solstice
   x-axis: direction of motion of earth on winter solstice
   y-axis: defined by right hand rule

 Coordinate system 2:

   zero:    position on earth's surface where unit rod is placed vertically
   x-axis:  east direction
   y-axis:  north direction
   z-axis:  radial direction

 Variables:

   h:       hours since midday GMT
   d:       days since winter solstice
   n:       latitude in radians north of equator
   e:       longitude in radians east of Greenwich 
   
 Let i,j,k denote the unit row vectors of CS2 expressed in CS1.
 Let A = (i,j,k) be a 3x3 matrix, and let B = inv(A).
 Then if a and b represent the same translation from the PoV of CS1 and CS2:

       a = b*A
       b = a*B

 Now in CS1 a unit vector in the direction of sunlight is

       u = matrix([sin(theta),0,-cos(theta)])       (where theta = 2*pi*d/365)

 So in CS2 the unit vector in the direction of sunlight is

       v = u*B

 Scale up by k = 1/v[2] and project into the i,j plane and you get the position of
 the tip of the unit length rod on the paper, i.e. in CS2

  x-coord of shadow = v[0]/v[2]
  y-coord of shadow = v[1]/v[2]
"""

import matplotlib.pyplot as plt
from numpy import pi, matrix, sin, cos, linalg
import numpy as np
import datetime

def rotate_about_x(phi):
    '''
    Return matrix R with property that for any row matrix w, w*R is w 
    rotated about the x axis by phi
    '''
    return matrix([[1,0,0],
                   [0,cos(phi),sin(phi)],
                   [0,-sin(phi),cos(phi)]])

def rotate_about_y(phi):
    '''
    Return matrix R with property that for any row matrix w, w*R is w 
    rotated about the y axis by phi
    '''
    return matrix([[cos(phi),0,-sin(phi)],
                   [0,1,0],
                   [sin(phi),0,cos(phi)]])
    

def rotate_about_z(phi):
    '''
    Return matrix R with property that for any row matrix w, w*R is w 
    rotated about the y axis by phi
    '''
    return matrix([[cos(phi),sin(phi),0],
                   [-sin(phi),cos(phi),0],
                   [0,0,1]])

def get_shadow_xy(h,d,n,e,w=None):
    '''
    Return x,y,z where
    
    x = x-coordinate of shadow of tip of rod in CS2
    y = y-coordinate of shadow of tip of rod in CS2
    z > 0 if shadow is cast in x-y plane, z < 0 otherwise

    h:       hours since midday GMT
    d:       days since winter solstice
    n:       latitude in radians north of equator
    e:       longitude in radians east of Greenwich 
    w:       if w == None then CS2 not on wall, otherwise CS2 on wall facing w radians anticlockwise from north
    '''
    
    # how far round sun is earth
    theta = 2*pi*d/365

    # Calculate matrix A whose rows are the unit vectors i,j,k of CS2 represented in CS1 coordinates
    A = matrix([[1.,0.,0.],[0.,1.,0.],[0.,0.,1.]])
    if w != None:
        A = A*rotate_about_x(pi/2)
        A = A*rotate_about_z(pi + w)
    A = A*rotate_about_x(-n)
    A = A*rotate_about_y(theta + e + h*2*pi/24)
    A = A*rotate_about_x(-pi*23.5/180)          # axial tilt

    # unit vector in direction of light from sun
    u = matrix([-sin(theta), 0, -cos(theta)])

    # convert u into CS2 coordinates
    v = linalg.solve(A.transpose(),u.transpose()).flatten().tolist()[0]

    # return position of shadow on paper
    return -v[0]/v[2], -v[1]/v[2], -v[2]

def days_since_dec21_to_label(days_since_dec21):
    '''
    convert days_since_dec21 to label, e.g. "Dec 23/Dec 19"
    '''
    dec21 = datetime.date(2000,12,21)
    date1 = dec21 + datetime.timedelta(days=days_since_dec21)
    date2 = dec21 + datetime.timedelta(days=365-days_since_dec21)
    if (date1.day,date1.month) == (date2.day,date2.month):
        return date1.strftime('%b %d');
    return date1.strftime('%b %d')+'/'+date2.strftime('%b %d')

def equation_of_time(days_since_dec21):
    '''
    Return correction in minutes for eccentricity of orbit.
    See https://www.pveducation.org/pvcdrom/2-properties-sunlight/solar-time.
    '''
    d = days_since_dec21 - 10
    B = 2*pi*(d-81)/365.
    return 9.87*sin(2*B) - 7.53*cos(B) -1.5*sin(B)

def plot_fixed_lat_long(n,e,w=None):
    '''
    Plot transit of shadow for a fixed latitude and longitude

    n:       latitude in radians north of equator
    e:       longitude in radians east of Greenwich 
    w:       if w == None then CS2 not on wall, otherwise CS2 on wall facing w radians anticlockwise from north
    '''

    # Set up canvas
    plt.ion()
    plt.figure()
    ax = plt.gca()
    ax.set_aspect('equal')
    plt.scatter(0,0, marker='x',zorder=3, color='black')

    x_min = -4.5
    x_max = +4.5
    y_min = -4.5
    y_max = +4.5
    if w:
        y_max = 0.0

    # lines labelled by date
    days_since_dec21 = map(lambda n: 365*n/12., range(7))
    for d in days_since_dec21:
        x_coords = []
        y_coords = []
        hours_since_midday_GMT = np.arange(-12,12,0.1)
        for h in hours_since_midday_GMT:
            (x, y, z) = get_shadow_xy(h,d,n,e,w)
            if z <= 0 or x < x_min or x > x_max or y < y_min or y > y_max:
                continue
            x_coords.append(x)
            y_coords.append(y)
        N = len(x_coords)
        if N < 2:
            continue        
        plt.plot(x_coords,y_coords, color='black')
        
        # add label
        label = days_since_dec21_to_label(d);
        x = x_coords[0]
        y = y_coords[0]
        plt.text(x,y, label, size=9, color='black', ha="center", va="center", bbox=dict(ec='1',fc='1'))

    # lines labelled by hour
    for h in range(-12,12):

        x_coords = []
        y_coords = []
        for d in range(365):

            # Correct for eccentricity of orbit
            h_new = h+equation_of_time(d)/60.
            
            (x, y, z) = get_shadow_xy(h_new,d,n,e,w)
            if z <= 0 or x < x_min or x > x_max or y < y_min or y > y_max:
                continue
            x_coords.append(x)
            y_coords.append(y)
        N = len(x_coords)
        if N < 2:
            continue
        plt.plot(x_coords,y_coords, color='black')

        # add label
        label="%02d00" % (12 + h)
        if h == 0:
            label="1200\nGMT"
        x = x_coords[0]
        y = y_coords[0]
        plt.text(x,y, label, size=9, color='black', ha="center", va="center", bbox=dict(ec='1',fc='1'))
        
if __name__ == '__main__':
    n = 0.91129839
    e = 0.0019329521
    w = None
#    w = 1.93 # back of home, cambridge
#    w = pi*150/180. # tennis court road
    plot_fixed_lat_long(n,e,w)

    # Optionally plot out triangle to cut out and stand up as dial
    plt.plot([0,1,0],[0,0,-1], color='black')