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