What is Sundial?

A sundial measures time by the position of the sun. The measurements rely on position of a predictable shadow cast by a fixed object on a (usually flat) surface marked with the hours of the day. There are several types of sundials. Traditionally they show local apparent solar time, however presently, in view of commonly used zone time, more practical are designs that measure time referred to a standard meridian (every 1 hour or 15° in geographic longitude). More advanced sundials show the mean solar time. These so called analemmatic sundials correct apparent solar time according to the equation of time (a –14 to +16 minute variation from mean solar time) and are scaled so as to indicate approximate standard or daylight saving time. These usually have hour markings shaped like "figure eights" (analemmas), along which shadow of a particular point on the gnomon (that point is called the nodus) moves in the course of year.

Why This Writing?

It is generally said that to design a certain type of sundial is a relatively simple matter while another requires more involved calculations. To study effects of setting inaccuracies of an existing sundial I was looking for formulae describing a sundial with arbitrarily oriented style or gnomon and arbitrarily oriented dial face. People do offer (for purchase) applications that can output a ready design of certain types of sundials and I got impression that the algorithms they employ behind their computation they keep kind of secret. Therefore I had to devise such a generalized algorithm all by myself. In fact in the end it turned out to be simple, indeed so simple that I decided to give the idea away for anybody interested by means of this page. Here I present a full FORTRAN code with ample comments embedded that should allow for easy adaptation to other programming languages. Hopefully, this shall promote future free access to applications of the above mentioned sort.

Algorithm for an Arbitrary Sundial

Below is a complete subprogram RodShadow coded in Fortran programming language that when called with given location latitude, equatorial position of the Sun, direction of gnomon and orientation of dial face returns the length and direction of the gnomon shadow as cast on that face. The shadow direction is measured on the plane of face and from the direction of largest possible inclination (i.e. from the line of greatest slope, which in case of a circular face is toward its top or bottom), Eastward positively. This reference convention can be easily modified to have the direction measured from North as seen on that plane (see the comment just prior to the end line of the subroutine).

subroutine RodShadow(phi,HA,Dec, af,hf, ag,hg, shadA,shadL) c Procedure to compute shadow length and direction cast onto c arbitrarily oriented plane (face of a dial) by arbitrarily c oriented gnomon illumined by the Sun c Input arguments c phi - geodetic latitude of sundial location [degrees] c HA - hour angle of the Sun [hours] c Dec - declination of the Sun [degrees] c af, hf - azimuth and altitude of the dial/face plane [degrees] c (the azimuth is counted Eastward from North) c ag, hg - azimuth and altitude of the gnomon [degrees] c Output arguments c shadA - shadow direction angle (with respect to face top) [degrees] c shadL - shadow length in units of gnomon length c We shall work in Cartesian system tied to the horizon system of c coordinates with axes directed toward North (x), East (y) and c zenith (z) and origin located at one of the ends of the gnomon. c The length of gnomon will be taken for unit of length. implicit real*8 (a-h,o-z) data d2r/0.0174532925199433d0/ ! degrees to radians factor c First convert equatorial coordinates of the Sun into horizontal ones: c azimuth (AzSun) and altitude (AltSun) call Eq2HorN(phi*d2r,HA*15*d2r,Dec*d2r,AltSun,AzSun) c The Sun direction cosines are thus: cxs = dcos(AltSun)*dcos(AzSun) cys = dcos(AltSun)*dsin(AzSun) czs = dsin(AltSun) c The gnomon cosines cxg = dcos(d2r*hg)*dcos(d2r*ag) cyg = dcos(d2r*hg)*dsin(d2r*ag) czg = dsin(d2r*hg) c Cosines of the normal to the plane of face cxf = dcos(d2r*(hf + 90))*dcos(d2r*af) cyf = dcos(d2r*(hf + 90))*dsin(d2r*af) czf = dsin(d2r*(hf + 90)) c Vector normal to plane containing the Sun and the gnomon and angle c Sun to gnomon or stylus base (or root) to gnomon tip call VecCrossP(cxs,cys,czs, cxg,cyg,czg, px,py,pz,angSG) c Unit vector along the shadow call VecCrossP(px,py,pz, cxf,cyf,czf, sx,sy,sz,thetaSF) c Angle from gnomon to its shadow call VecCrossP(sx,sy,sz, cxg,cyg,czg,qx,qy,qz, thetaGS) c Length of the shadow is the sum of projected gnomon length and c component due to Sun not at right angle to the shadow direction c dcos(thetaGS) + dsin(thetaGS)*dtan(angSG + thetaGS - pi/2), c which can be shown to simplify to shadL = dsin(angSG)/dsin(angSG+thetaGS) c Let a reference vector on the face have these cosines: cxr = dcos(d2r*hf)*dcos(d2r*af) cyr = dcos(d2r*hf)*dsin(d2r*af) czr = dsin(d2r*hf) c i.e. along the greatest slope or towards high/low point of (assumedly) c a circular face, then finally we get the shadow angle through call c VecCrossP(sx,sy,sz, cxr,cyr,czr,x,y,z, shadArad) and scale shadArad c dividing it by d2r. The same obtains with following two lines cosine=sx*cxr + sy*cyr + sz*czr shadA = dacos(dmax1(-1d0,dmin1(1d0,cosine)))/d2r if(sy.lt.cyr) shadA = -shadA c The above makes the angle negative West of face top or face bottom. c An alternative: if(sy.lt.0d0) shadA=-shadA !meaning angle negative c West of the meridian, would require a reference vector (cxr,cyr,czr) c towards North and on the face. Such vector can be obtained through c call VecCrossP(0d0,1d0,0d0,cxf,cyf,czf, cxr,cyr,czr, anglerad) c This alternative with af = 0 was checked to yield, as expected c the same output as the greatest slope reference actually used here. end subroutine VecCrossP(Ax,Ay,Az, Bx,By,Bz, Cxn,Cyn,Czn,theta) c Calculates the cross product of vectors A and B, C = A × B implicit real*8 (a-h,o-z) c Get the product Cx = Ay*Bz - Az*By Cy = Az*Bx - Ax*Bz Cz = Ax*By - Ay*Bx c Normalize it c = dsqrt(Cx*Cx + Cy*Cy + Cz*Cz) Cxn = Cx/c Cyn = Cy/c Czn = Cz/c c Get angle between the two vectors theta = dacos(Ax*Bx + Ay*By + Az*Bz) end subroutine Eq2HorN(fi,t,D,Alt,Az) c Converts equatorial coordinates (D,t - declination, hour angle) into c horizontal coordinates (Alt,Az - altitude, azimuth counted from North) c at geographical latitude fi [all angles expressed in radians] implicit real*8 (a-h,o-z) data TwoPi/6.2831853071795864769d0/ sinfi = dSIN(fi) cosfi = dSQRT(1 - sinfi*sinfi) Az = dATAN2(dSIN(t),sinfi*dCOS(t) - dTAN(D)*cosfi) Alt= dASIN(dSIN(D)*sinfi + dCOS(D)*cosfi*dCOS(t)) Az = dMOD(Az + TwoPi*1.5,TwoPi) end

TCfA Sundial

The Torun Centre for Astronomy has an equatorial sundial (gnomon directed to the North celestial pole) with a circular face that is not perpendicular to the gnomon but is inclined at certain angle with respect to it as well as to the horizon plane. It shows the zonal apparent solar time, rather than the local one. It means in particular that at the time of local noon it indicates 14 minutes less than 12 hours, since at nearest zonal meridian (the one 15° East of Greenwich and 14 minutes of time West of the sundial location) the noon will arrive 14 minutes later.

The design below has been calculated by calling the RodShadow procedure for every half-hour between 5:30 and 18:00. The hour angle (HA) corresponding to these time marks range from –6.5 to 6.0 but here were increased by 14 minutes of the location longitude correction (e.g. for 5:30 mark on the dial face the HA was set to –6.5 + 14/60 or –6.2666...).

Fig. 1. Design of the face of the TCfA sundial. Note the different units (cm) of top and right scale labels (these units were calculated assuming 67.2 cm for the gnomon length, which itself is the unit of labels at the other two scales). Each encircled dot (in red) marks exact direction of gnomon shadow at the hour (of CET when the Equation of Time is zero) rendered in blue. The shadow extends from the red cross (gnomon footprint) and in this case usually goes out far beyond the face (see Fig. 3). The continuum line is a spline approximating circular edge of the face (as drawn by OriginPro v. 7.5 application of OriginLab).

Fig. 2. Effects of inaccurate setting of the TCfA sundial. The gray color shows the design of Fig. 1 above while the red and blue marks indicate corresponding shadow directions calculated for incorrectly mounted dial (inclination of the gnomon with respect to face is assumed correct and fixed to 32.7°). The red dots correspond to the dial inclination/altitude in error of 5° (RodShadow arguments hf and hg both increased by 5), and the blue stars correspond to the entire dial rotated in azimuth by 5° East from North (RodShadow arguments af and ag both increased by 5).

The results presented in Fig. 2 demonstrate that inaccuracies in adjusting of a sundial orientation of the order of a few degrees are serious source of errors (in this example they cause errors of some 2 minutes of time per one degree in orientation) and cannot be compensated by rotation of its face around its normal (i.e. while keeping it in the same plane). Note that inclination error has opposite effect before and after noon (i.e. it makes the sundial both morning and evening fast or slow depending on sense of the rotation) while rotation in azimuth — the same.

Computing for an Analemmatic Sundial

A sundial that shows mean solar time (in practice the time our clocks and watches show) obtains the same way as for apparent time except that a –14 to 16 minute correction must be taken into account. The correction, called the Equation of Time, depends on date of year (and also, although very slightly, on the year itself). The apparent hour of the day is the hour angle of the Sun plus 12 hours. The apparent Sun advences ahead of the mean Sun by the Equation of Time. Therefore the required modification of an algorithm for the common sundial, such as considered above in this report, cosists only in increasing of the hour angle by the Equation of Time. Of course, if we desire the sundial to indicate a zone time, we must also add to this angle the diffrence in geographic longitude between our location and the zone meridian (assuming the difference to be positive East of the meridian and expressed in units of time, i.e. each degree converted into 4 minutes).

When using the above procedure RodShadow for chosen hour of the day HD, one first finds the declination of the Sun Dec and the Equation of Time EqT for a particular date and calculates the hour angle:

Fig. 3. Track of shadow of the gnomon tip on the plane of face of the TCfA sundial at integer hours of Central European Time (CET) computed approximately every 10 days of 2006.

Fig. 4. One of analemmas of Fig. 3, that of 12 hours, expanded for details. Note that the horizontal axis here is magnified by a factor of about 5 with respect to vertical axis. The red labels along the curve mark the date, month (1 to 12 stand for January to December) and day (1st, 11th and 21st day of month) of the year, for which the shadow end position (red circle) has been determined. The horizontal extent of this analemma corresponds to 30 minutes of time. The origin on the horizontal scales (labels 0,00) coincides with the direction of North (local apparent noon), so that during most of February the TCfA sundial, though designed for the mean zone time (CET), indicates quite correctly also the apparent solar time.

    (Posted 2006.09.27; last revised 2006.12.29.)