When I was working for a small wind turbine manufacturer I was asked to calculate the number of hours per year the houses in the neighborhood of our prototype wind farm would be exposed to shadow flicker from our wind turbine. Now, there are software that calculates this sort of thing with a few clicks of the mouse (the dominant code beeing WindPro). However, these programs cost money, and they are no fun. Instead I solved the problem myself. And to do that I first had to immerse myself in the theory of sundials.

First consider a flagpole casting a shadow on a flat field. There are only a few parameters that influence where the tip of the shadow will be at a particular time. There are some astronomical parameters, namely: the earth axial tilt, the earth rotational speed, the length of the solar year, the eccentricity of the orbital ellipse, and the angle between the major axis of the orbital ellipse and the line of solstices.

Then there are some parameters of the site: height of pole, latitude and longitude. And there is a convention where the year is said to begin relating to the winter solstice.

The rest is a matter of geometry really. Basically one must keep track of exactly how earth is facing the sun at any particular instant in time.

Coordinate tranformations

This is how I think about it. I want to find a transformation from geographic coordinates (which rotate with earth) to a coordinate system with origo in the center of the earth that always has the sun in the direction of the y-axis and the z-axis normal to the orbital plane. In this system it should be easy to find the tip of the shadow by finding the intersection between the line passing from the sun, through the tip of the pole and the sphere that is earth. Once this is done we go back to geographic coordinates by applying the inverse tranformation.

To implement this we first want a coordinate transformation from the geographic coordinate system (longitude $\phi$, latitude $\delta$, height $h$) rotating with earth to a cartesian system rotating with earth. This is a simple adaption of the spherical coordinate transformation. The only difference being that the latitude replaces the normal angle from the north pole.

$$ \begin{array}{rcl} x & = & (R+h) \sin(\pi/2-\delta)\cos \phi \ y & = & (R+h) \sin(\pi/2-\delta)\sin \phi \ z & = & (R+h) \cos (\pi/2 - \delta) \ \end{array} $$

What we want to do next is to transform a point expressed this way to a nonrotating coordinate system. The transformation matrix for rotation about the z-axis is

$$ R_z(\Omega t) = \left(\begin{array}{ccc} \cos \Omega t & -\sin \Omega t & 0 \ \sin \Omega t & \cos \Omega t & 0 \ 0 & 0 & 1 \end{array}\right) $$

where $\Omega$ is the rotational speed of earth $2 \pi / 86164.09890369 $ rad/s.

Next we want to express this in a system aliged with the earth orbital plane. Thus we want to rotate about the universal y-axis by the axial tilt angle $\alpha = 23.4{}^\circ$. The transformation matrix for this rotation is

$$ R_x(\alpha) = \left(\begin{array}{ccc} 1 & 0 & 0 \ 0 & \cos \alpha & - \sin \alpha \ 0 & \sin \alpha & \cos \alpha \end{array} \right) $$

Finally we want to rotate this system so we face the sun.

$$ R_z(-\theta(t)) = \left( \begin{array}{ccc} \cos \theta(t) & \sin \theta(t) & 0 \ -\sin \theta(t) & \cos \theta(t) & 0 \ 0 & 0 & 1 \end{array} \right) $$

where $\theta(t)$ is the orbital angle. This is a function that must be computed based on orbital mechanics.

The total transformation from ($\phi$,$\delta$) is thus

$$ \left(\begin{array}{c} X \ Y \ Z \end{array}\right) = T(\phi,\delta) = R_z(-\theta(t)) R_x(\alpha) R_z(\Omega t) \left(\begin{array}{c} (R+h) \sin(\pi/2-\delta)\cos \phi \ (R+h) \sin(\pi/2-\delta)\sin \phi \ (R+h) \cos (\pi/2 - \delta) \end{array}\right) $$

The inverse of this tranformation is obtained by

$$ \begin{array}{rcl} h & = & \sqrt{x^2+y^2+z^2}-R \ \phi & = & \mbox{arctan}\left(\frac{y}{x}\right) \ \delta & = & \pi/2 - \mbox{arccos}\left(\frac{z}{r}\right) \end{array} $$

where $$(x,y,z)$$ is obtained by inversing the rotations

$$ \left(\begin{array}{c} x \ y \ z \end{array}\right) = \left(R_z(-\theta(t)) R_x(\alpha) R_z(\Omega t)\right)^{-1} \left(\begin{array}{c} X \ Y \ Z \end{array}\right) $$

Intersection between ray and earth surface

Assume the earth surface is a perfect sphere centered in origo. Also assume all light rays from the sun are coming in the -Y-direction. Now consider a point $\vec{p}$ representing the location of the flagpole tip. The tip may or may not cast a shadow on earth.

To find the location of the shadow simply apply the formula for line sphere intersection

$$ \ d = -\hat{Y} \cdot \vec{p} - \sqrt{\left(-\hat{Y} \cdot \vec{p}\right)^2-\vec{p} \cdot \vec{p}+R^2} $$

where $R$ is the earth radius and $d$ is the distance from $\vec{p}$ to the intersection. $\hat{Y}$ is a unit vector in the Y-direction. If $d$ is real and greater than zero the location of the shadow is

$$ \vec{s} = \vec{p} - d \hat{Y}. $$

I have implemented this theory as a python module. The source is over at gitHub

A sundial for my back yard

As an application of this theory I have decided to construct a sundial custom made for my backyard. This is a bit of a project of its own and I will delay the details for a later blog. The following picture is what I have come up with so far. This sundial shows approximate date and also takes daylight saving into account.