# Data from Land Information New Zealand
# See Moon/provenance.txt and Moon/moon.R
lowTides <-
    c("00:55", "01:47", "02:38", "03:29", "04:20", "05:11", "06:05", 
      "07:01", "08:00", "09:01", "10:01", "10:58", "11:51", "00:16", 
      "01:03", "01:44", "02:23", "03:00", "03:36", "04:12", "04:49", 
      "05:29", "06:12", "07:02", "07:59", "09:02", "10:06", "11:08", 
      "12:08", "00:34", "01:29")
phases <- data.frame(date=c("2010-01-01 08:13", "2010-01-07 23:40",
                       "2010-01-15 08:12", "2010-01-23 23:54",
                       "2010-01-30 07:18"),
                     phase=c("Full", "3Q", "New", "1Q", "Full"))
phases$date <- as.POSIXct(as.character(phases$date))
lowTideDate <- as.POSIXct(paste("2010-01-", 
                                sprintf("%02d", 1:31), 
                                " ", lowTides,
                                sep=""))
lowTideHour <- as.POSIXct(paste("2010-01-01 ", 
                                lowTides,
                                sep=""))

mainHours <- ISOdatetime(2010, 1, 1,
                         c(0, 4, 8, 12), 
                         rep(0, 4), 
                         rep(0, 4))



library(jpeg)
moon <- jpeg::readJPEG(system.file("extra", "GPN-2000-000473.jpg",
                             package="RGraphics"))


moonPhase <- function(x, y, phase, size=.05) {
  # size is in inches
  n <- 17
  angle <- seq(0, 2*pi, length=n)
  xx <- x + cos(angle)*xinch(size)
  yy <- y + sin(angle)*yinch(size)
  if (phase == "New")
    fill <- "black"
  else
    fill <- "white"
  polygon(xx, yy, col=fill)
  if (phase == "1Q")
    polygon(xx[(n/4):(n*3/4) + 1],
            yy[(n/4):(n*3/4) + 1],
            col="black")
  if (phase == "3Q")
    polygon(xx[c(1:(n/4 + 1), (n*3/4 + 1):n)],
            yy[c(1:(n/4 + 1), (n*3/4 + 1):n)],
            col="black")
}


# Original image from NASA
# http://grin.hq.nasa.gov/ABSTRACTS/GPN-2000-000473.html
rasterMoon <- jpeg::readJPEG("Moon/GPN-2000-000473.jpg")
par(pin=c(3.5, 1.75), oma=c(0, 3, 0, 0), xaxs="i", yaxs="i", cex=.7)
plot.new()
rect(0, 0, 1, 1, col="black")
rasterImage(rasterMoon, .25, 0, .75, 1*813/703)
par(new=TRUE, xaxs="r", yaxs="r", las=1)
plot(lowTideDate, lowTideHour, type="n",
     ylim=range(mainHours), axes=FALSE, ann=FALSE)
# dashed reference lines
abline(v=phases$date,
       col="white", lty="dashed")
for (subset in list(1:13, 14:29, 30:31)) {
  lines(lowTideDate[subset], lowTideHour[subset],
        lwd=2, col="white")
  points(lowTideDate[subset], lowTideHour[subset],
         pch=16, col="white")
}
box()
axis.POSIXct(1, lowTideDate)
axis.POSIXct(2, at=mainHours, format="%H:%M")
mtext("Time of Low Tide (NZDT)", side=2, line=4, las=0, cex=.7)
mtext("Auckland, New Zealand January 2010", side=1, line=3, cex=.7)
axis(3, at=phases$date, labels=FALSE)
par(xpd=NA)
ymax <- par("usr")[4]
for (i in 1:nrow(phases))
    moonPhase(phases$date[i], ymax + yinch(.2), 
              phases$phase[i])
mtext("Phases of the Moon", side=3, line=3, cex=.7)