Sundial calculations

December 28, 2013
By

(This article was first published on Dan Kelley Blog/R, and kindly contributed to R-bloggers)

After experimenting with calculations for what I eventually came to realize were analemma-based sundials (with shadow cast by a vertical pole), I remembered that the common sundial has a wedge as the shadow-maker. A bit of research told me that the wedge is called a gnomon. It is a right triangle with one vertex (the centre vertex, shall we say) having angle equal to the local latitude. If this wedge is placed upright on a horizontal plane with the centre vertex aligned south and the 90deg vertex aligned north, then the shadow will produce a line that indicates the hour of the day. This works throughout the year, with approximately quarter-hour adjustments being required through the seasons.

The R code given below the diagram creates an outline of the expected edge of the shadow of the gnomon. To illustrate the variation in angle through the year (which relates to the so-called equation of time), colours are used to indicate four significant times during the year.

Printed at full scale, the diagram might be suitable for laying out the horizontal scale for a sundial. Naturally, readers must alter the latitude and longitude if they do not live in Halifax, Nova Scotia.

A few notes:

  1. The gnomon hypotenuse will point to the pole star (Polaris) when the apparatus is aligned properly towards the north.
  2. Calling the function with debug=1 will show dots along the radial lines. These are the shadows of virtual points lying along the hypotenuse of the gnomon, and provide a check against errors in the calculation (since they should lie along a line if the gnomon angle matches the latitude).
  3. Noon is not aligned with North because the longitude is not an even multiple of 15 degrees.
  4. The length of the shadow provides extra information, but here this information is not shown (the lengths are normalized in lines 33 to 35 of the code.)

graph

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
## gnonom style sundial
if (!interactive())
    png("sundial_with_gnomon.png", width=7, height=6, unit="in", 
        res=200, pointsize=13)
library(oce)
sundial <- function(lon=-63.60, lat=44.65,
                    days=c("2014-03-20", "2014-06-20", "2014-09-23", "2014-12-21"),
                    keys=c("Spring equinox", "Summer solstice",
                           "Autumn equinox", "Winter solstice"),
                    debug=0)
{
    col <- 1:4
    glwd <- 8
    timezone <- floor(0.5 + lon / 15)
    L <- 1                           # horiz gnomon length (arbitrary)
    nhours <- 24
    first <- TRUE
    for (season in 1:4) {
        hours <- seq.POSIXt(as.POSIXct(days[season], tz="UTC"),
                            by="1 hour", length.out=nhours)
        for (hour in seq_along(hours)) {
            t <- hours[hour]
            tlocal <- t + 3600 * timezone
            sa <- sunAngle(t, lon=lon, lat=lat)
            gy <- seq(0, L, length.out=10)
            gx <- rep(0, length(gy))
            gz <- gy * tan(lat * pi / 180)
            R <- gz / tan(-sa$altitude * pi / 180) # radius of shadow
            theta <- (90 - sa$azimuth) * pi / 180
            par(mar=rep(2, 4))
            x <- gx + R * cos(theta)
            y <- gy + R * sin(theta)
            len <- max(sqrt(x^2 + y^2))
            x <- x / len * L
            y <- y / len * L
            if (sa$altitude > 0) {
                if (first) {
                    first <- FALSE
                    D <- L * 5
                    plot(x, y, type='n', pch=20, asp=1,
                         xlim=c(-L, L), ylim=c(-L/5, L),
                         axes=FALSE, col=col[season], xlab="", ylab="")
                    ## Draw gnomon as a gray bar
                    segments(0, 0, 0, L, lwd=glwd, col='gray')
                    grid()
                    abline(h=0, lwd=1.5*par('lwd'), lty='dotted')
                    abline(v=0, lwd=1.5*par('lwd'), lty='dotted')
                    mtext("South", side=1, at=0)
                    mtext("West", side=2, at=0)
                    mtext("North", side=3, at=0)
                    mtext("East", side=4, at=0)
                    legend("topright", lwd=glwd, col="gray",
                           legend="Gnomon")
                    legend("bottomright",
                           legend=sprintf("%.3fE %.3fN", lon, lat))
                    legend("topleft", lwd=1, col=1:4, cex=3/4,
                           legend=keys)
                    box()
                    points(0, 0, pch=20, cex=2)
                    segments(0, 0, x, y, col=col[season])
                } else {
                    segments(0, 0, x, y, col=col[season])
                    if (debug)
                        points(x, y)
                }
                if (season==2) {
                    xend <- tail(x, 1)
                    yend <- tail(y, 1)
                    text(1.05*xend, 1.05*yend, format(tlocal, "%H"))
                }

            }
        }
    }
}
if (!interactive())
    dev.off()

To leave a comment for the author, please follow the link and comment on their blog: Dan Kelley Blog/R.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Sponsors

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)