Calculating The Layout of a Bifilar Sundial using Vectors

I. INTRODUCTION

Bifilar sundials, by virtue of their offset threads crossing at different heights, present a unique challenge to the dialist who would like to take their 'game' to a higher level. Horizontal dial planes are well-documented, and even vertical, declined dials have readily available solutions. While there are software programs available for calculating arbitrarily oriented bifilar dial planes, they are often either not free, or difficult to use if one is unfamiliar with CAD (Computer Aided Drafting) software. This article attempts to outline a procedure, using 3D vectors, which ultimately leads to a surprisingly simple set of equations for determining the positional data required to layout and construct a bifilar sundial.

Following the derivation does require the reader to be familiar with vectors and the use of dot and cross products. These topics are well-documented and can be easily found on the internet. As such, they are not covered here.

Interested readers may use and download a free, javascript dial calculator using the algorithms outlined in this article. The calculator is intended to be simple and easy to use with any of the major browsers in use at the time of this writing (Firefox, Opera, Chrome, Safari and Internet Explorer). The javascript is documented and contains less than 1000 lines of source code (including the documentation). Individuals familiar with programming are invited to download and tailor it to their own liking.

II. PRELIMINARIES

Basic Problem

The general outline of the task is as follows:

1. Each thread of the bifilar dial represents a unique line in three dimensional space. When the sun is visible to the dial plane there is a third, unique line that starts at the sun and passes through both threads. This results in "crossed" shadows somewhere on the dial plane. The basic goal is to find the equation of the third unique line and where it intersects the dial plane. This point represents the location on the dial plane where the thread shadows cross.

   The basic approach to solving for the point of "crossed shadows" involves three basic steps:

   1. Determine the unique plane in three dimensions formed by the horizontal thread and the solar position direction.
   2. Determine the unique point of intersection between the plane determined in step 1 and the vertical thread.
   3. Determine where the line from the sun, and passing through the point calculated in step 2, intersects the dial plane. This point represents the location of the crossed shadows.

2. Next, the point of intersection, calculated in three dimensions, must be projected onto the 2-dimensional dial plane. This new, 2D coordinate ( DP_x , DP_y ), represents a single data point ultimately used in the construction of the sundial.
3. Steps 1 and 2 are repeated for each relevant solar hour angle (and all relevant solar declinations) until enough data has been calculated to construct the entire sundial.
4. It should be noted that the general technique used in solving this problem, which assumes the threads are oriented at right angles, can easily be extended to threads oriented at just about any orientation that allows the formation of crossing shadows.

Coodinate Systems

• 3-Dimensional
We use the familiar XYZ cartesian coordinate system (right-handed) for all calculations (see Fig. 1). The horizontal plane at the location of the sundial represents the X-Y plane: the positive x-axis is due East, while the positive y-axis is due North. The vertical represents the positive z-axis. All references to "3D" coordinates will use this representation. The azimuth and altitude of the sun are measured as shown in Figure 3.

• 2-Dimensional
Any reference to 2D will refer to the plane of the dial (see Fig. 2). It represents the surface on which the hour lines are ultimately to be drawn. We again use a right-handed cartesian coordinate system. The x-axis is always horizontal, even for arbitrarily oriented dial planes. For a horizontal dial plane (the one with which most people are familiar), the 2D x-axis and y-axis coincide with the 3D x and y axes, respectively. For an arbitrarily oriented dial plane, the x-axis is defined as the intersection between the dial plane and a horizontal plane passing through the origin. The x-axis of the dial plane will, in general, be a horizontal line in the 3D system, confined to the x-y plane and passing through the origin.

• Origin
The origin, whether 3D or 2D, is defined as the point on the dial plane directly beneath the crossed threads. Another way to think of this is to imagine the sun's rays striking the dial plane perpendicularly, and the origin is simply the point where the thread shadows cross. The origin of the 2D system and the origin of the 3D system represent the same physical location.

Fig. 1: XYZ Axes

Fig. 2: 2D Dial Plane

Fig. 3: Azimuth and Altitude

Dial Face Inclination & Declination

When most people think of a sundial, they envision an object with a horizontal plane and a gnomon of some type, usually slanted, attached at a right angle. But there is nothing to prevent a dial maker from orienting the dial face in a completely arbitrary direction. In general, a dial may be tilted with respect to the vertical (inclination) and rotated azimuthally (declination). This last term should not to be confused with the declination associated with coordinates on the Celestial Sphere. Dial face inclination is defined as the angle measured from the vertical to the dial face normal (see Fig. 4). Dial face declination is the angle measured in the horizontal plane, from due south to the vertical plane which contains the dial's zenith and surface normal (see Fig. 5). Declinations east of south are negative; west of south, positive. The following table should help to illustrate how each measurement is defined.

Fig. 4: Dial Plane Inclination

Fig. 5: Dial Plane Declination

Dial Orientation Inclination Declination horizontal 0° 0° vertical, south-facing 90° 0° vertical, east-facing 90° -90° vertical, west-facing 90° 90° vertical, facing south-west 90° 45°

Spherical Triangle (PZS)

This article does not discuss spherical trigonometry, nor how to solve the PZS triangle. There is already a lot of information available on the web for this, so there is no point in re-inventing the wheel. It is assumed you can already determine the sun's altitude and azimuth at your latitude for the relevant solar hour angles needed to produce a working sundial.

III. CALCULATIONS

In general, calculations may be classified as falling into one of two categories: 1) those which are required only once; and 2) those which must be performed for each position of the sun for which a dial layout data point is desired. In addition, throughout this article no attempt is made to systematically "simplify" the various unit vectors or formulae to expressions involving the basic initial inputs. In many instances, the process of "simplification" leads to longer, more cumbersome equations, although this is not always the case. Purists are free to disagree and rewrite equations as they see fit. In the interest of economy, readers uninterested in the derivation are free to skip directly to the formulas for the dial plane coordinates expressed in terms of low-level inputs.

One-Time Calculations

1 :: Equation and unit normal of the dial plane:
The general equation describing a plane in three dimensions is given by Ax + By + Cz = D. By defining the dial plane as always passing through the origin, this simplifies to:

In addition, a very useful property of the general equation is that the surface normal of the plane is given by: If the coefficients A, B and C are the direction cosines of the vector n, the equation reduces to: If we know the inclination and declination of the dial plane, we can easily calculate the coefficients as follows:

x-direction : A = -sin(Ι) sin(Δ)
y-direction : B = -sin(Ι) cos(Δ)
z-direction : C = cos(Ι)
where Ι = the inclination of the dial plane
and Δ = the declination of the dial plane

2 :: Define the dial plane x-axis:
Earlier it was discussed that the x-axis of the dial plane represents a line which passes through the 3D origin and is confined to the horizontal x-y plane. To determine the equation of this line (in the 3D system), we want to find the intersection between the horizontal plane (i.e. z=0, ) and the plane of the dial (Ax + By + Cz = 0). The direction of the line of intersection between two planes can be determined by taking the cross product of vectors normal to each plane. The unit normal of the dial plane was calculated earlier and the unit normal of the x-y plane is k, so the vector h, representing the dial plane x-axis, can be determine as follows:

Again, h represents a vector in the x-y plane of the 3D coordinate system that also lies completely in the 2D plane of the dial. Normalizing the vector gives the following:

3 :: Define the dial plane y-axis:
Previously, we have defined the dial plane unit normal, n, and x-axis unit vector, h. The cross product of these two vectors results in a third unit vector perpendicular to both. From Fig. 6 it can be seen that this direction, labeled by vector v defines the dial plane y-axis (see Eq. 4).

4 :: Define the horizontal thread:
The height, or offset, of the horizontal thread from the dial plane will be referred to as O_H. This offset allows one to define a unique point on the thread directly above the dial origin by using the dial plane unit normal, n. The position vector of this unique point is given by P_H = O_Hn (Eq. 5a).

Taking advantage of the fact that the dial plane x-axis and the horizontal thread are parallel lines, we can write the equation defining the horizontal thread as shown in Eq. 5b. No term for the z-direction is given because the direction cosine of the dial plane's x-axis is undefined (i.e. h lies somewhere in the horizontal plane of the 3D system).

5 :: Define the vertical thread:
The vertical thread can be defined in a manner similar to that used for the horizontal thread. The height, or offset, of the vertical thread from the dial plane will be referred to as O_V. The position vector of a point on the vertical thread directly above the dial plane origin, P_V = O_Vn, is shown in Eq. 6a. The set of equations which uniquely define the vertical thread as a line in three dimensions is given in Eq. 6b.

Repetitive Calculations

This section outlines the calculations required for each solar position (i.e. each combination of hour angle and solar declination) used in laying out the sundial.

1 :: Calculate the solar position unit vector:
This operation is straight-forward and is merely a determination of the sun's direction cosines from its (previously determined) altitude and azimuth:

Direction cosines:
x-dir : α = cos(alt) sin(azm)
y-dir : β = cos(alt) cos(azm)
z-dir : γ = sin(alt)
where, alt = solar altitude
and azm = solar azimuth

The unit position vector becomes: It should be noted that unit vector s points from the dial toward the sun. This representation is an important consideration in the next section.

2 :: Determine visibility of the sun:
Given the dial plane may have any orientation (i.e. inclined and/or declined), one cannot be sure the sun will shine on the dial plane just because it is above the horizon (e.g. a vertical, west-facing dial plane does not "see" the sun until it has transited the local meridian). For a particular solar position to cast a shadow on the dial plane, the sun must be both: 1) above the horizon; and 2) in front of the dial plane. We can quickly and easily check for solar visibility by calculating a single dot product for each condition:

• If the dot product between the solar position vector and the vertical is positive (Eq. 8a), the first condition is satified, and the sun is above the horizon.
• Similarly, if the dot product between the solar position vector and the dial plane normal is positive (Eq. 8b), the second condition is satisfied, and the sun is in front of the dial plane.

When using the outline of calculations presented in this article as part of an executable program, these two very simple checks will quickly discriminate between solar positions which do (or don't) cast a shadow on the dial plane. Even though the amount of time saved is probably small when using today's average computer, it is still a good idea to use them. For spreadsheet applications, however, this step may not be applicable.

3 :: Calculate the plane containing the solar direction & the horizontal thread:
We know that two intersecting lines define a plane. Mathematically, the cross product of the unit vectors defining the direction of each line produce a third vector which represents the normal to the plane defined by the intersecting lines. In this step, we want to determine the plane formed by the dial's horizontal thread and the solar direction unit vectors. The variable φ will be used to refer to this plane. The general equation describing this plane is given in Eq. 9a. The coefficients A_φ, B_φ and C_φ can be determined by taking the cross product between the solar direction unit vector, s, and the unit vector defining the direction of the horizontal thread, h. This operation is shown in Eq. 9b.

While the result of Eq. 9b determines the normal to φ, it does not completely define the plane. The coefficient D_φ has not yet been determined. To do so, we make use of the fact that point P_H (Eq. 5a) lies within φ, and so must satisfy the general plane equation in 9a. Using the results of Eqs. 5a and 9b, D_φ can be calculated as shown in Eq. 9c. By defining ψ = αA + βB + γC we can simplify Eq. 9c a bit further to Eq. 9d. It is worth noting that ψ is nothing more than s • n, the dot product between the solar position unit vector and the dial plane unit normal. The plane φ is now completely determined.

4 :: Calculate the point of intersection between plane φ and the vertical thread:
We now wish to calculate the point on the vertical thread where it intersects plane φ. To determine this point, which will be referred to as g, it is necessary to solve equations 6b and 9a simultaneously. Doing so yields the following equations for determining the coordinates of g.

All of the quantities on the right-hand side of Eq. 10a are known, so it is now possible to calculate the z-axis coordinate, g_z, of the point where φ intersects the vertical thread. Substituting back into equations 10b and 10c yields g_x and g_y, respectively.

5 :: Calculate the 3D position where the shadows cross on the dial plane:
We are now ready to calculate the 3D position where the shadows of the horizontal and vertical threads cross on the dial plane. In the previous step, we calculated the point on the vertical thread intersected by the the plane formed by the solar direction and the horizontal thread. So, the 3D line (Eq. 11a) connecting this point, g, on the vertical thread and the sun also passes through the horizontal thread. The intersection of this line with the dial plane represents the position on the dial plane where the shadows of the threads cross.

To determine the coordinates of this point, ( c_x , c_y , c_z ), we must solve Eqs. 1 and 11a simultaneously, resulting in Eqs. 11b → 11d.

6 :: Calculate dial plane coordinate → DP_x:

Before proceeding, it is worthwhile to be reminded that the position vector of the point of intersection calculated in the previous step, c = c_xi + c_yj + c_zk, is a vector which lies in the plane of the dial. Realizing this, it becomes a relatively simple matter to determine DP_x by taking the dot product of the 3D coordinate position vector c with the unit vector, h, which represents the x-axis in the plane of the dial (See Fig. 7). The operation yields the component of c along the dial plane's x-axis, which is the dial plane x-coordinate DP_x (see Eq. 12).

7 :: Calculate dial plane coordinate → DP_y:
Computation of the 2D y-coordinate, DP_y, proceeds in a manner similar to that for DP_x: by having previously defined the unit vector v parallel to the dial plane y-axis, one can now determine the component of c along the 2D y-axis by calculating the dot product c • v, as shown in Eq. 13.

8 :: Simplified formulas for DP_x & DP_y: TOP
Performing all the necessary algebraic substitutions yields expressions for each dial plane coordinate in terms of the solar direction cosines (α, β, γ), the dial plane unit normal direction cosines (A, B, C), and the horizontal and vertical thread offsets, O_H and O_V, respectively. These formulas are given in Eqs. 14a & 14b. Note that for a horizontal dial both A & B are zero, so an inclination of 0° results in division by zero; however, these formulas can still be used. Simply use a very small value for the inclination (e.g. inclination = 0.001°). (A close inspection of the formulas using the original inputs of solar azimuth & altitude, and dial inclination & declination reveals why this operation is possible. It will not be demonstrated here. Interested or skeptical readers are invited to see for themselves.)

IV. CONCLUSION

For difficult to calculate inclined/declined bifilar sundials we find that a technique using primarily vector operations can yield a surprisingly simple set of final equations for the dial plane coordinates (Eqs. 14a & 14b). The steps necessary to arrive at these expressions are summarized as follows:

Step	Description	Calculation
1†	Determination of dial plane coefficients using the values of inclination and declination.	A = -sin(Ι)sin(Δ) B = -sin(Ι)cos(Δ) C = cos(Ι) where Ι = the inclination of the dial plane and Δ = the declination of the dial plane
2†	Determine the dial plane x-axis and the 3D line equations defining the horizontal thread.	Dial Plane X-Axis: Horizontal Thread:
3†	Determine the dial plane y-axis and the 3D line equations defining the vertical thread.	Dial Plane Y-Axis: Vertical Thread:
4‡	Determine the solar direction cosines.	α = cos(alt) × sin(azm) β = cos(alt) × cos(azm) γ = sin(alt) where, alt = solar altitude and azm = solar azimuth
5‡	Determine the equation of the plane (φ) formed by the horizontal thread and solar direction unit vector.	, ,     where,
6‡	Determine the point of intersection, g, between plane φ and the vertical thread.
7‡	Determine the point of intersection, ( cx , cy , cz ), between the dial plane and the line formed by point g and the solar direction unit vector, s.
8‡	Determine the dial plane coordinates ( DPx , DPy ).
† : "one time" calculation. ‡ : calculation required for each plottable solar position.