swear not by the moon, the fickle moon

monthly changes in her circle orb

Tip

If you’re wondering where I’m getting my titles and headings, they come from Shakespeare’s sonnets.

The moon has always fascinated me. Eversince I was a child, I was entranced by it; her constant pursuit everywhere I went and her sneaky tactics during the day. This year, before the solar eclipse, I decided to do a little digging and replicate one algorithm used to determine the phases of the moon. If you’re interested, you can find it here in this very dense book by Jean Meeus.

In case you didn’t know, the times of the new moon, first quarter, full moon, and last quarter are the times at which the moon’s longitude, as viewed from our planet, aligns with the sun to form 0\(^\circ\), 90\(^\circ\), 180\(^\circ\), and 270\(^\circ\), respectively. As the moon orbits the earth and rotates on its axis, it overlaps with the Sun in varying degrees, and thus light is reflected off of it. This means that to calculate the phases of the moon, we would technically need to know the geocentric¹ longitude of the Moon and the Sun…and I’m not doing that🙃 so giving up a sliver of accuracy seems the right solution if we want to avoid an unorthodox amount of math.

¹ as viewed from the Earth

et tu, brute?

Astronomers use Julian Ephemeris Days (JDE), or Dynamical Time (DT), as a measure of time to keep track of celestial bodies, and you probably guessed that its name may be derived from emperor Julius Caesar’s reform to create a calendar that always remained aligned to the sun, but that would be incorrect. It has nothing to do with Brutus’s uncle.

The first year under this calendar is Jan 1, 4,713 BC, which means that 2024 would represent the year 6737 of the Julian era—a pretty nerdy thing to point out, but nobody is complaininng, probably.

To find the JDE for any celestial event, we simply subtract 2,451,550 days from Ephemeris Time², which is the number of days from 4713 BC to 1952…don’t quote me on that😅

² adopted in 1952

The times of the mean phases of the moon are then given by:

\[ JDE = 2,451,550.09766 + 29.530588861 k \\ + 0.00015437 T^{2} \\ - 0.000000150 T^{3} \\ + 0.00000000073 T^{4} \]

where an integer value of \(k = 0\) gives a new moon³, and an increase by .25 gives the next phase, 1st quarter, full moon, and the third or last quarter. Any other value of \(k\) is meaningless.

³ \(k = 0\) corresponds to the new moon of Jan 6, 2000, and negative values of \(k\) are lunar phases before the year 2000.

crux of thought…math

The approximate value of \(k \approx (year - 2000) × 12.3685\) and the \(year\) variable should be expressed in decimals to include the month and day of the year. \(T\) is the time in Julian centuries since 2000, and we can accurately calculate it using this formula \((\frac{k}{1236.85})\), and a term we haven’t seen yet, \(E = 1 - 0.002516T - 0.0000074T^{2}\) that describes the eccentricity of Earth’s orbit around the sun.

The angles for the sun’s and the moon’s mean anomaly⁴, the moon’s latitude and the longitude of the ascending node⁵ of the lunar orbit will also be calculated.

⁴ distance from the center of the moon/sun measured from the side that is closest to earth assuming movement at a uniform speed in a circular orbit.

⁵ the moon has an ascending and descending node where it crosses the ecliptic hemisphere from the south or the north, respectively.

We need to take into account 14 planetary arguments and specific angles and time corrections for each phase of the moon to obtain any past, future and current Julian epochs. First, I’ll create a simple function modf that will split an object of type double into its integer and decimal parts. I’ll use this function to figure out the moon phase later down the road. A second function modr will reduce any angle to the range between 0–360, and a third function deg2rad will convert the degrees to radians, if necessary.

#> Modulus fraction function
modf <- function(k){
     stopifnot(is.numeric(k))
     kint <- trunc(k)
     kfrac <- abs(k - kint)
     kround <- c(kint,kfrac)
     return(kround)
}

#> Modulus reduce angle function
modr <- function(angle){
  stopifnot(is.numeric(angle))
  reduced <- angle%%360
  return(reduced)
}

rad2deg <- function(rad){
  stopifnot(is.numeric(rad))
  degrees <- (rad*180)/(pi)
  return(degrees)
}

deg2rad <- function(deg){
  stopifnot(is.numeric(deg))
  radians <- (deg*pi)/180
  return(radians)
}

test_it <- modf(-3.35)
test <- modr(400)
print(c(test_it,test))
#> [1] -3.00  0.35 40.00

the code

Note

Almost every function used here except the main function, modf,modr,deg2rad, rad2deg, and the dex.date functions are taken from Jean Meeus’ “Astronomical Algorithms”.

The goal is to calculate the date of any lunar phase any number of years into the future. To do this, we’ll first create functions for the mean anomaly, latitude, and longitude of the moon.

#> Times of lunar phases in JDE
JDE <- function(k,t){
 jde <- 2451550.09766 + (29.530588861*k) + 
   0.00015437*(t^2) - 0.000000150*(t^3) + 
   0.00000000073*(t^4)
 jde <- round(jde,digits = 6)
 return(jde)
}
#> Calculate E
E <- function(k,t){
  e <- 1 - (0.002516*t) - 0.0000074*(t^2)
  e <- round(e,digits = 6)
  return(e)
}

#> Sun's mean anomaly at time JDE
M <- function(k,t){
  m <- 2.5534 + (29.10535670*k) - 0.0000014*(t^2) - 0.00000011*(t^3)
  m <- modr(round(m, digits = 6))
  m <- deg2rad(m)
  return(m)
}
#> Moon’s mean anomaly at time JDE
M.prime <- function(k,t){
  m_prime <- 201.5643 + (385.81693528*k) + 0.0107582*(t^2) + 
    0.00001238*(t^3) - 0.000000058*(t^4)
  m_prime <- modr(round(m_prime,digits = 6))
  m_prime <- deg2rad(m_prime)
  return(m_prime)
}
#> Moon's argument of latitude
F <- function(k,t){
  f <- 160.7108 + (390.67050284*k) - 0.0016118*(t^2) - 
    0.00000227*(t^3) + 0.000000011*(t^4)
  f <- modr(round(f, digits = 6))
  f <- deg2rad(f)
  return(f)
}

omega <- function(k,t){
  o <- 124.7746 - (1.56375588*k) + 0.0020672*(t^2) + 0.00000215*(t^3)
  o <- modr(round(0, digits = 6))
  o <- deg2rad(o)
  return(o)
}

Now it’s time to come up with a function fo each planetary argument (14 total). Most of these argument take \(k\) with one exception \(A_{1}\), which also requires \(T\). The output of the planetary arguments are expressed in degrees and may have to be reduced to 0–360 and converted to radians. Since our main function does use every planetary argument to predict the lunar phases it may be better to provide the arguments as functions instead of vectors. if doing this is too torublesome, we’ll try something different inside the main function.

A1 <- function(k,t){
  a1 <- 299.77 + 0.107408*k - 0.009173*(t^2)
  a1 <- modr(round(a1, digits = 6))
  a1 <- deg2rad(a1)
  return(a1)
}
A2 <- function(k){
  a2 <- 251.88 + 0.016321*k
  a2 <- modr(round(a2,digits = 6))
  a2 <- deg2rad(a2)
  return(a2)
}
A3 <- function(k){
  a3 <- 251.83 + 26.651886*k
  a3 <- modr(round(a3,digits = 6))
  a3 <- deg2rad(a3)
  return(a3)
}

A4 <- function(k){
  a4 <- 349.42 + 36.412478*k
  a4 <- modr(round(a4,digits = 6))
  a4 <- deg2rad(a4)
  return(a4)
}
A5 <- function(k){
  a5 <- 84.66 + 18.206239*k
  a5 <- modr(round(a5, digits = 6))
  a5 <- deg2rad(a5)
  return(a5)
}
A6 <- function(k){
  a6 <- 141.74 + 53.303771*k
  a6 <- modr(round(a6,digits = 6))
  a6 <- deg2rad(a6)
  return(a6)
}
A7 <- function(k){
  a7 <- 207.14 + 2.453732*k
  a7 <- modr(round(a7, digits = 6))
  a7 <- deg2rad(a7)
  return(a7)
}
A8 <- function(k){
  a8 <- 154.84 + 7.306860*k
  a8 <- modr(round(a8,digits = 6))
  a8 <- deg2rad(a8)
  return(a8)
}
A9 <- function(k){
  a9 <- 34.52 + 27.261239*k
  a9 <- modr(round(a9,digits = 6))
  a9 <- deg2rad(a9)
  return(a9)
}
A10 <- function(k){
  a10 <- 207.19 + 0.121824*k
  a10 <- modr(round(a10, digits = 6))
  a10 <- deg2rad(a10)
  return(a10)
}
A11 <- function(k){
  a11 <- 291.34 + 1.844379*k
  a11 <- modr(round(a11,digits = 6))
  a11 <- deg2rad(a11)
  return(a11)
}
A12 <- function(k){
  a12 <- 161.72 + 24.198154*k
  a12 <- modr(round(a12,digits = 6))
  a12 <- deg2rad(a12)
  return(a12)
}
A13 <- function(k){
  a13 <- 239.56 + 25.513099*k
  a13 <- modr(round(a13,digits = 6))
  a13 <- deg2rad(a13)
  return(a13)
}
A14 <- function(k){
  a14 <- 331.55 + 3.592518*k
  a14 <- modr(round(a14, digits = 6))
  a14 <- deg2rad(a14)
  return(a14)
}

the main thing

Using the preciously coded functions, we’ll create the main function with and add conditions for each moon phase to accurately implement the time corrections.

Something to keep in mind is that the date variable will have to be expressed as a decimal for the algorithm to work. This means that I’ll have to transform the date outside the main function before doing any computations.

The output will be in date format so we’ll need to look at R’s date manipulation capabilities. The built.in as.Date function won’t work for us because it doesn’t handle time (HH:MM:SS), so we’re left with two choices, the chron package which handles time but doesn’t account for timezones or the built.in POSIX classes that do manipulate timezones. I’m going with POSIXlt so that I can extract time information without pulling from external libraries.

…a detour

I’ll test a simple function to transform a character string in the format "YYYY/MM/DD" to a POSIXlt date and then return a decimal to represent the number of days completed in any given year.

Something to note: Some years will be leap or century years, and this will also have to be dealt with. Additionally, we’re calculating the number of days completed up to the current date, which means we’ll update our function to use day - 1 to find the decimal expression for the specified date.

dex.date <- function(date){
  if (!inherits(date, c("POSIXt", "POSIXct", "POSIXlt", "Date"))) {
    #> Transform to POSIXlt
  POSIX_date <- as.POSIXlt(date)
#> if the format is not found by the function, you can provide                                    
#> your format explicitly.
  }
  year <- POSIX_date$year + 1900 
#> Pull the year attribute and add 1900, default is to subtract                                   
#> 1900 to provided year
  day <- as.numeric(format(POSIX_date,"%j")) # day of the year 0 - 365
  dy_yrs <- ifelse((year %% 4 == 0 & year %% 100 != 0) | (year %% 400 == 0), 366,365.25)
#> Get the time components
  hour <- POSIX_date$hour
  minute <- POSIX_date$min
  second <- POSIX_date$sec
#> Calculate the fraction of the day
  fraction_day <- (hour * 3600 + minute * 60 + second) / (86400)
  dex_date <- year + (day -1 + fraction_day)/dy_yrs
  return(dex_date)
}

#> Let's test it
test <- dex.date("2027/2/19 23:30:45")
print(test)
#> [1] 2027.13684

Alright! The function is accurate up to the fourth significant figure⁶ That’s great, but now I have to create a second function to convert from JDE to Gregorian time. I’ll use October 15, 1582 as the first day in the Gregorian calendar, which means that the number of Julian days from 4,713 BC to Oct. 15, 1582 is \(2299161.05\)

⁶ We can test this by using the decimal_date function in the lubridate package. We get 2027.1369.

convert.JDE <- function(jde){
  jde <- jde + 0.5
  Z <- trunc(jde)
  F <- jde - Z
  if (Z < 2299161) {
    A <- Z
  } else {
    alpha <- trunc((Z - 1867216.25) / 36524.25)
    A <- Z + 1 + alpha - trunc(alpha / 4)
  }
  B <- A + 1524
  C <- trunc((B - 122.1) / 365.25)
  D <- trunc(365.25 * C)
  E <- trunc((B - D) / 30.6001)

  day <- B - D - trunc(30.6001 * E) + F
  month <- ifelse(E < 14, E - 1, E - 13)
  year <- ifelse(month > 2, C - 4716, C - 4715)

  # Extract the fractional day part for time
  hour <- (day - trunc(day)) * 24
  minute <- (hour - trunc(hour)) * 60
  second <- (minute - trunc(minute)) * 60
  
  # Create a POSIXlt object for the Gregorian date and time
  gregorian_date <- as.POSIXlt(sprintf("%04d-%02d-%02d %02d:%02d:%06.3f",
                                       year, month, trunc(day),
                                       trunc(hour), trunc(minute), second),
                               format = "%Y-%m-%d %H:%M:%OS", tz = "UTC")
  return(gregorian_date)
}

After checking the accuracy of the convert.JDE function by using the following examples, we can create a third function to convert from Gregorian time to JDE. This is extra work, but it assures us that our JDE calculation is accurate and error-free🤞

Dates for testing JDE conversion function

Calendar	JDE
2000/01/01 12:00:00	2451545.0
1999/01/01	2451179.5
1600/01/01	2305447.5
837/04/10 7:12:00	2026871.8
1987/06/19 12:00:00	2446966.0
1900/01/01	2415020.5

The JDE function we created cannot be used here because the variable \(k\) is not needed for calendar transformation. If we were to use it, we would always get a Julian date that’s off by 1 - 7 days…I’ve already tried it.

greg.convert<- function(date){
  if (!inherits(date, c("POSIXt", "POSIXct", "POSIXlt", "Date"))) {
    #> Transform to POSIXlt
  POSIX_date <- as.POSIXlt(date)
  }
  Y <- POSIX_date$year + 1900
  M <- POSIX_date$mon + 1
  D <- POSIX_date$mday
  H <- POSIX_date$hour
  mn <- POSIX_date$min
  S <- POSIX_date$sec
  Y_up <- ifelse(M > 2, Y,Y - 1)
  M_up <- ifelse(M > 2,M, M + 12)
  A <- floor(Y_up/100)
  if(Y_up < 1582){ # If looking at Julian calendar
    B = 0
  } else{
    B = 2 - A + floor(A/4)
  }
  D_dec <- (H + mn/60 + S/3600)/24
  jde <- floor(365.25*(Y_up + 4716)) + floor(30.6001*(M_up +1)) + D + D_dec + B - 1524.5
  return(jde)
}
testy <- greg.convert("2027/02/19 23:30:45, tz = UTC")
print(testy)
#> [1] 2461456.48

We can use the Dates table to test the accuracy of this function. One example that caught my attention is the date 3837-04-10 07:12:00. This particular date returned \(3122594.8\) and the convert.JDE function was 3837-04-10 07:11:60 EDT, which made me think that it was a day with a leap second, but those are unpredictable and not known so far in advance. We’ll have to deal with this error as well.

the main function

With all the parts in place, we can get back to writing the main function. Ideally, a date will be entered, we will calculate \(k\), then \(T\) to obtain the JDE date which, by the way, is not accurately representing the date of the chosen lunar phase yet. Our algorithm will them find \(E, M, M^\prime, F\), and \(\Omega\). The fourteen planetary corrections will be applied next, and additional correction for the selected lunar phase will be calculated. We’ll then use these corrections to adjust the JDE date and finally convert it to the Gregorian date using the convert.JDE function. We can check the accuracy of our results by using the greg.convert function. Here we go:

calc.MoonPhase <- function(k){
  stopifnot(is.numeric(k))
  
  t <- k / 1236.85
  
  # Mean phase JDE
  jde <- JDE(k, t)
  
  # Anomaly, latitude, and longitude of Sun/Moon
  e <- E(k, t)
  m <- M(k, t)
  mpr <- M.prime(k, t)
  f <- F(k, t)
  om <- omega(k, t)
  
  kfrac <- modf(k)[2]
  kint <- modf(k)[1]
  
  planet_args <- (
    0.000325 * sin(A1(k, t)) + 
    0.000165 * sin(A2(k)) + 
    0.000164 * sin(A3(k)) +
    0.000126 * sin(A4(k)) + 
    0.000110 * sin(A5(k)) + 
    0.000062 * sin(A6(k)) +
    0.000060 * sin(A7(k)) + 
    0.000056 * sin(A8(k)) + 
    0.000047 * sin(A9(k)) +
    0.000042 * sin(A10(k)) + 
    0.000040 * sin(A11(k)) + 
    0.000037 * sin(A12(k)) +
    0.000035 * sin(A13(k)) + 
    0.000023 * sin(A14(k))
  )
  
  if (kfrac == 0) {
    adj <- (
      -0.40720 * sin(mpr) + 
      0.17241 * e * sin(m) + 
      0.01608 * sin(2 * mpr) + 
      0.01039 * sin(2 * f) + 
      0.00739 * e * sin(mpr - m) - 
      0.00514 * e * sin(mpr + m) +
      0.00208 * (e^2) * sin(2 * m) - 
      0.00111 * sin(mpr - 2 * f) - 
      0.00057 * sin(mpr + 2 * f) +
      0.00056 * e * sin(2 * mpr + m) - 
      0.00042 * sin(3 * mpr) + 
      0.00042 * e * sin(m + 2 * f) +
      0.00038 * e * sin(m - 2 * f) - 
      0.00024 * e * sin(2 * mpr - m) - 
      0.00017 * sin(om) -
      0.00007 * sin(mpr + 2 * m) + 
      0.00004 * sin(2 * mpr - 2 * f) + 
      0.00004 * sin(3 * m) +
      0.00003 * sin(mpr + m - 2 * f) + 
      0.00003 * sin(2 * mpr + 2 * f) - 
      0.00003 * sin(mpr + m + 2 * f) +
      0.00003 * sin(mpr - m + 2 * f) - 
      0.00002 * sin(mpr - m - 2 * f) - 
      0.00002 * sin(3 * mpr + m) +
      0.00002 * sin(4 * mpr)
      )
    
    true_jde <- (jde + adj) + planet_args
     true_date <- convert.JDE(true_jde)
  return(true_date)
    
  } else if (kfrac == .25 || kfrac == .75) {
    adj2 <- (
      -0.62801 * sin(mpr) + 
      0.17172 * e * sin(m) - 
      0.01183 * e * sin(mpr + m) + 
      0.00862 * sin(2 * mpr) + 
      0.00804 * sin(2 * f) + 
      0.00454 * e * sin(mpr - m) + 
      0.00204 * (e^2) * sin(2 * m) - 
      0.00180 * sin(mpr - 2 * f) - 
      0.00070 * sin(mpr + 2 * f) -
      0.00040 * sin(3 * mpr) - 
      0.00034 * e * sin(2 * mpr - m) + 
      0.00032 * e * sin(m + 2 * f) + 
      0.00032 * e * sin(m - 2 * f) - 
      0.00028 * (e^2) * sin(mpr + 2 * m) + 
      0.00027 * e * sin(2 * mpr + m) - 
      0.00017 * sin(om) - 
      0.00005 * sin(mpr - m - 2 * f) + 
      0.00004 * sin(2 * mpr + 2 * f) - 
      0.00004 * sin(mpr + m + 2 * f) + 
      0.00004 * sin(mpr - 2 * m) + 
      0.00003 * sin(mpr + m - 2 * f) + 
      0.00003 * sin(3 * m) + 
      0.00002 * sin(2 * mpr - 2 * f) + 
      0.00002 * sin(mpr - m + 2 * f) - 
      0.00002 * sin(3 * mpr + m)
      )
    
    w <- (
      0.00306 - 0.00038 * e * cos(m) + 
      0.00026 * cos(mpr) - 
      0.00002 * cos(mpr - m) +
         0.00002 * cos(mpr + m) +
      0.00002 * cos(2 * f)
    )
    
    if (kfrac == .25) {
      if(kint > 0){
        true_jde <- (jde + adj2) + w + planet_args
      } else {  # First Quarter before 2000 AD
        true_jde <- (jde + adj2) - w + planet_args
      }
    } else {  # kfrac == 0.75
    if (kint > 0) {  # Last Quarter after 2000 AD
        true_jde <- (jde + adj2) - w + planet_args
      } else {  # Last Quarter before 2000 AD
        true_jde <- (jde + adj2) + w + planet_args
      }
    }
    true_date <- convert.JDE(true_jde)
    return(true_date)
    
  } else if (kfrac == .5) {
    adj3 <- (
      -0.40614 * sin(mpr) + 
      0.17302 * e * sin(m) + 
      0.01614 * sin(2 * mpr) + 
            0.01043 * sin(2 * f) + 
      0.00734 * e * sin(mpr - m) - 
      0.00514 * e * sin(mpr + m) + 
            0.00209 * (e^2) * sin(2 * m) - 
      0.00111 * sin(mpr - 2 * f) - 
      0.00057 * sin(mpr + 2 * f) + 
            0.00056 * e * sin(2 * mpr + m) - 
      0.00042 * sin(3 * mpr) + 
      0.00042 * e * sin(m + 2 * f) + 
            0.00038 * e * sin(m - 2 * f) - 
      0.00024 * e * sin(2 * mpr - m) - 
      0.00017 * sin(om) - 
            0.00007 * sin(mpr + 2 * m) + 
      0.00004 * sin(2 * mpr - 2 * f) + 
      0.00004 * sin(3 * m) + 
            0.00003 * sin(mpr + m - 2 * f) + 
      0.00003 * sin(2 * mpr + 2 * f) - 
      0.00003 * sin(mpr + m + 2 * f) + 
            0.00003 * sin(mpr - m + 2 * f) -
      0.00002 * sin(mpr - m - 2 * f) - 
            0.00002 * sin(3 * mpr + m) + 
      0.00002 * sin(4 * mpr)
      )
    
    true_jde <- (jde + adj3) + planet_args
    true_date <- convert.JDE(true_jde)
  return(true_date)
    
  } else if (kfrac > 0 && kfrac < .25) {
    return(cat("Moon is in Waxing Crescent Phase. No valid date returned"))
  } else if (kfrac > .25 && kfrac < .5) {
    return(cat("Moon is in Waxing Gibbous Phase. No valid date returned"))
  } else if (kfrac > .5 && kfrac < .75) {
    return(cat("Moon is in Waning Gibbous Phase. No valid date returned"))
  } else {
    return(cat("Moon is in Waning Crescent Phase. No valid date returned"))
  }
}

We can test the function using a few examples of past, present and future lunar phases. We’ll use dex.date to convert the year to decimal format before finding \(k\). Once we find the value of \(k\), we’ll round up or down depending on the moon phase⁷ we’re trying to estimate. There should be a small discrepancy between actual and predicted lunar phase dates because:

⁷ \(k = 0\) for the new moon, \(k = .25\) for first quarter, \(k = .50\) for full moon and \(k = .75\) for the last quarter.

• We’ve ignored several arguments and corrections.
• The U.S Naval Observatory doesn’t include seconds in their date and time estimates.
• Future dates estimates are as accurate as the amazing,but fallible algorithms that predicted them.
• There are no high precision records of the orbital positions of the sun and the moon prior to the invention of the telescope in 1609 A.D.

the very error of the moon 🌓

Important

In the interest of space, I’ve hidden the data generation code. You can download a .txt version of it here if you’d like to follow along. - USN Dates

USN_dates <- data.frame(Dates = c(
"1930/01/08 03:11:00",
"1930/01/14 22:21:00",
"1930/01/21 16:07:00",
"1930/01/29 19:07:00",
"1930/11/28 06:18:00",
"1930/12/06 00:40:00",
"1930/12/12 20:07:00",
"1930/12/20 01:24:00",
"1971/03/04 02:01:00",
"1971/03/12 02:34:00",
"1971/03/20 02:30:00",
"1971/05/24 12:32:00",
"1971/09/27 17:17:00",
"1971/10/04 12:19:00",
"1971/10/11 05:29:00",
"1971/12/17 19:03:00",
"1701/01/17 09:42:00",
"1701/01/24 12:52:00",
"1701/01/31 04:21:00",
"1701/02/07 22:57:00",
"1701/11/07 12:03:00",
"1701/11/15 16:53:00",
"1701/11/23 06:09:00",
"1701/11/29 21:53:00",
"2010/09/15 05:50:00",
"2010/09/23 09:17:00",
"2010/10/01 03:52:00",
"2010/10/07 18:44:00",
"2010/11/13 16:39:00",
"2010/11/21 17:27:00",
"2010/11/28 20:36:00",
"2010/12/05 17:36:00",
"2023/04/27 21:20:00",
"2023/05/05 17:34:00",
"2023/05/12 14:28:00",
"2023/05/19 15:53:00",
"2023/08/24 09:57:00",
"2023/08/31 01:35:00",
"2023/09/06 22:21:00",
"2023/09/15 01:40:00",
"2072/01/28 10:36:00",
"2072/02/04 04:56:00",
"2072/02/10 23:27:00",
"2072/02/19 02:04:00",
"2100/01/18 12:35:00",
"2100/01/26 02:51:00",
"2100/02/01 21:17:00",
"2100/02/09 04:56:00",
"1833/03/28 22:50:00",
"1833/04/04 14:43:00",
"1833/04/12 00:09:00",
"1833/05/19 13:36:00"),
                        Phases = rep(c(
"First Quarter", "Full Moon", "Last Quarter", "New Moon"),13)
)
#> file_path <- here::here("06","moon","USN_dates.txt")
#> write.table(USN_dates, file = file_path, sep = "\t", col.names = TRUE)

However, the lunar phase estimates are not too far off from the observed dates. Let’s find out how accurate our function really is by calculating the mean prediction error for a bigger set of dates taken from the U.S Naval Observatory

Mean error of time difference (in mins) between predicted and actual dates

Year	First Quarter	Full Moon	Last Quarter	New Moon
1701	0.920	0.930	0.330	0.200
1833	0.380	0.560	0.010	0.290
1930	0.220	0.310	0.265	0.610
1971	0.700	0.550	0.590	0.275
2010	0.465	1.150	1.150	0.955
2023	1.350	1.565	1.475	1.305
2072	2.090	1.640	1.950	1.720
2100	1.920	1.830	2.120	2.070

Visualizing it may provide a different perspective on the time discrepancy between the actual and the predicted dates of the lunar phases. We can see how the algorithm is more accurate as we reach the present and loses accuracy as we move away from it.

take thy flight

Finally, we can accurately predict when the werewolves will turn 🐺, just kidding!

In all seriousness, the calc.MoonPhase function is accurate to within 2 minutes for a period of 400 years, but deteriorates as the time interval increases⁸. So, if an error of a few minutes is not important, we may use the algorithm to accurately predict the date of any lunar phase. However, we must also keep in mind that there will always be a small prediction error in every algorithm because the time interval between consecutive lunations will vary as the Sun perturbs the lunar path, and we cannot account for every force or action exerted on our “little” satellite.

⁸ the time difference in the expected date and the official date for the full moon of 0033/04/03 is 2 hrs, 50 mins and 24 seconds.