library(oce)
par(mfrow = c(3, 2))
y <- 2012
m <- 4
days <- 1:3
# Halifax sunrise/sunset (see e.g. https://www.timeanddate.com/worldclock)
rises <- ISOdatetime(y, m, days, c(13, 15, 16), c(55, 04, 16), 0, tz = "UTC") + 3 * 3600 # ADT
sets <- ISOdatetime(y, m, days, c(3, 4, 4), c(42, 15, 45), 0, tz = "UTC") + 3 * 3600
azrises <- c(69, 75, 82)
azsets <- c(293, 288, 281)
latitude <- 44.65
longitude <- -63.6
for (i in 1:3) {
t <- ISOdatetime(y, m, days[i], 0, 0, 0, tz = "UTC") + seq(0, 24 * 3600, 3600 / 4)
ma <- moonAngle(t, longitude, latitude)
oce.plot.ts(t, ma$altitude, type = "l", mar = c(2, 3, 1, 1), cex = 1 / 2, ylab = "Altitude")
abline(h = 0)
points(rises[i], 0, col = "red", pch = 3, lwd = 2, cex = 1.5)
points(sets[i], 0, col = "blue", pch = 3, lwd = 2, cex = 1.5)
oce.plot.ts(t, ma$azimuth, type = "l", mar = c(2, 3, 1, 1), cex = 1 / 2, ylab = "Azimuth")
points(rises[i], -180 + azrises[i], col = "red", pch = 3, lwd = 2, cex = 1.5)
points(sets[i], -180 + azsets[i], col = "blue", pch = 3, lwd = 2, cex = 1.5)
}
Run the code above in your browser using DataLab