The winter solstice has been on many minds lately. The days are about to start getting longer … but just how fast will they do that?

This post provides R that calculates and graphs day length and its variation, using `uniroot()`

to find sunrises and sunsets as indicated by solar altitude that is calculated with `sunAngle()`

in the oce package.

The day of the solstice is indicated with vertical lines. All times are in UTC, which is the conventional system for scientific work and the one required by `sunAngle()`

.

The first step in making the graph shown above is to load the `oce`

library and create a function that measures daylength by finding sunrise and sunset times. Note that `uniroot()`

, which is used to find times of zero solar altitude, needs lower and upper limits on `t`

, and these are calculated by looking back and then forward a half-day. This works well for application to Halifax, but in other timezones other offsets would be needed. Interested readers might want to devised a method based on the longitude, which can be transformed into a timezone.

1
2
3
4
5
6
7
8
9
10

library(oce)
daylength <- function(t, lon=-63.60, lat=44.65)
{
t <- as.numeric(t)
alt <- function(t)
sunAngle(t, longitude=lon, latitude=lat)$altitude
rise <- uniroot(alt, lower=t-86400/2, upper=t)$root
set <- uniroot(alt, lower=t, upper=t+86400/2)$root
set - rise
}

Next, use `lappy()`

to find the daylength for December, 2013.

1
2
3

t0 <- as.POSIXct("2013-12-01 12:00:00", tz="UTC")
t <- seq.POSIXt(t0, by="1 day", length.out=1*31)
dayLength <- unlist(lapply(t, daylength))

Set up to plot two panels, with narrowed margins.

1

par(mfrow=c(2,1), mar=c(3, 3, 1, 1), mgp=c(2, 0.7, 0))

Plot daylength in the top panel, indicating the day of the solstice with vertical lines.

1
2
3
4
5

plot(t, dayLength/3600, type='o', pch=20,
xlab="", ylab="Day length (hours)")
grid()
solstice <- as.POSIXct("2013-12-21", tz="UTC")
abline(v=solstice+c(0, 86400))

Plot daylength difference in the bottom panel, again indicating the solstice.

1
2
3
4

plot(t[-1], diff(dayLength), type='o', pch=20,
xlab="Day in 2013", ylab="Seconds gained per day")
grid()
abline(v=solstice+c(0, 86400))