Make a sundial

equationoftime.R

# The Equation of Time
# based on trigonometric approximation on Wikipedia
equoftime <- function(M){
 -7.657*sin(M) + 9.862*sin(2*M + 3.599)
}

x <- equoftime((1:365)*2*pi/365.242)
y <- ts(x, start = 4, frequency = 365)

# set up the plot
library(zoo)
plot(zoo(y, as.Date("2012-01-04") + 1:365),
     main = "Equation of time",
     xlab = "time", ylab = "minutes",
     type = "n")

# mark off the months
markers <- seq(from = as.Date("2012-01-01"),
               to = as.Date("2013-01-01"),
               by = "month")
abline(v = markers, col = "grey")

# mark off the weeks
markersw <- seq(from = as.Date("2012-01-01"),
               to = as.Date("2013-01-01"),
               by = "week")
abline(v = markersw, lty = 2, lwd = 0.5, col = "grey")

# mark off minutes on y-axis
horizs <- seq(-15, 15, by = 2.5) 
abline(h = horizs, lty = 2, lwd = 0.5, col = "grey")

# x-axis
abline(h = 0)

# solstices (red) and equinoxes (blue)
abline(v = as.Date("2012-03-20"), col = "red", lwd = 2)
abline(v = as.Date("2012-06-20"), col = "blue", lwd = 2)
abline(v = as.Date("2012-09-22"), col = "red", lwd = 2)
abline(v = as.Date("2012-12-21"), col = "blue", lwd = 2)

# finally the equation of time
lines(zoo(y, as.Date("2012-01-04") + 1:365), lwd = 3)

gnomon.R

# Make the gnomon

lat <- 38.2

x <- (0:180)*pi/180
plot(sin(x) ~ cos(x), asp = 1, type = "l", 
     xlab = "Fold in half at dashed line, cut dotted lines to make triangle, angle from circle center rests on plate, pointed South", ylab = "",
     main = bquote(Latitude == .(lat)*degree),
     axes = FALSE, frame.plot = TRUE, 
     xlim = c(-1,1), ylim = c(0,1))
abline(v = 0, lwd = 3, lty = 2)
abline(h = 0, lwd = 3)
abline(a = 0, b = tan(pi/2 - lat*pi/180), lty = 2)
abline(a = 0, b = tan(pi/2 + lat*pi/180), lty = 2)

sundial.R

# Sundial
# enter latitude (degrees N) 

lat <- 38.2
times <- 15*c(1:5)*pi/180
theta <- atan(sin(lat*pi/180) * tan(times))

thlow <- -atan(sin(lat*pi/180) * tan(15*7*pi/180)) - pi/2
thhi  <- atan(sin(lat*pi/180) * tan(15*7*pi/180)) + 3*pi/2


x <- -80:260*pi/180
par(mar=c(5.1,5.1,5.1,5.1))
plot(sin(x) ~ cos(x), asp = 1, type = "l", lty = 2,
     xlab = bquote(Latitude == .(lat)*degree), ylab = "",
     main = "Sundial",
     mar = c(5, 5, 5, 5) + 0.1,
     #sub = bquote(Latitude == .(lat)*degree),
     axes = FALSE, frame.plot = TRUE, 
     xlim = c(-1,1), ylim = c(sin(thlow),1))
abline(v = 0, lwd = 4)
abline(h = 0, lwd = 4)

for (k in 1:5){
  segments(x0 = 0, y0 = 0, x1 = 2*cos(pi/2 - theta[k]), y1 = 2*sin(pi/2 - theta[k]))
  segments(x0 = 0, y0 = 0, x1 = 2*cos(pi/2 + theta[k]), y1 = 2*sin(pi/2 - theta[k]))
}

mtext("XII", at = 0, cex = 2)
mtext("XIII", at = 1/tan(pi/2 - theta[1]), cex = 2)
mtext("XIV", at = 1/tan(pi/2 - theta[2]), cex = 2)
mtext("XV", at = 1/tan(pi/2 - theta[3]), cex = 2)
mtext("XVI",  at = 1/tan(pi/2 - theta[4]), side = 4, las = 1, cex = 2)
mtext("XVII", at = tan(pi/2 - theta[5]), side = 4, las = 1, cex = 2)
mtext("XVIII", at = 0, side = 4, las = 1, cex = 2)

mtext("IX", at = 1/tan(pi/2 + theta[1]), cex = 2)
mtext("X", at = 1/tan(pi/2 + theta[2]), cex = 2)
mtext("IX", at = 1/tan(pi/2 + theta[3]), cex = 2)
mtext("VIII",  at = -1/tan(pi/2 + theta[4]), side = 2, las = 1, cex = 2)
mtext("VII", at = -tan(pi/2 + theta[5]), side = 2, las = 1, cex = 2)
mtext("VI", side = 2, at = 0, las = 2, cex = 2)

segments(x0 = 0, y0 = 0, x1 = 2*cos(thlow), y1 = 2*sin(thlow))
segments(x0 = 0, y0 = 0, x1 = 2*cos(thhi), y1 = 2*sin(thhi))
mtext("XIX", at = cos(thlow), side = 1, cex = 2, padj = 1)
mtext("V", at = cos(thhi), side = 1, cex = 2, padj = 1)