tests/testthat/test_astronomical.R

In dankelley/oce: Analysis of Oceanographic Data

# vim:textwidth=80:expandtab:shiftwidth=4:softtabstop=4

# References used in this file:
##
# 1. Meeus, Jean, 1982. Astronomical formulae for calculators. Willmann-Bell. Richmond VA, USA. 201 pages.
# 2. Meeus, Jean, 1991. Astronomical algorithms.  Willmann-Bell, Richmond VA, USA. 429 pages.

library(oce)

RPD <- atan2(1, 1) / 45 # radians per degree

test_that("Angles", {
    # A randomly-chosen example on page 99 of Meeus (1991).
    # Hundreds of other cases could be found in that book, or his
    # earlier one, but the fomulae is self-evident so a single
    # value ought to be enough to ensure against editing errors that
    # might accidentally alter degree2hms().
    hms <- angle2hms(177.74208)
    expect_equal(hms$hourDecimal, 177.74208 * 24 / 360)
    expect_equal(hms$hour, 11)
    expect_equal(hms$minute, 50)
    expect_equal(hms$second, 58.0992)
    expect_equal(hms$string, "11h50m58s.10")
})


test_that("Times", {
    # [1] chapter 3 page 24-25
    # FIXME: previously this had the unintelligble tz="ET" but it is *exact* as is
    t <- ISOdatetime(1957, 10, 4, hour = 0, min = 0, sec = 0, tz = "UTC") + 0.81 * 86400
    expect_equal(julianDay(t), 2436116.31, tolerance = 0.01)
    # [1] example 15.a
    t <- ISOdatetime(1978, 11, 13, 4, 35, 0, tz = "UTC")
    jd <- julianDay(t)
    jca <- julianCenturyAnomaly(jd)
    expect_equal(jd, 2443825.69, tolerance = 0.01)
    expect_equal(jca, 0.788656810, tolerance = 1e-7) # fractional error 3e-8
    # [1] page 40
    t <- ISOdatetime(1978, 11, 13, 0, 0, 0, tz = "UTC")
    expect_equal(siderealTime(t), 3.4503696, tolerance = 0.0000001)
    t <- ISOdatetime(1978, 11, 13, 4, 34, 0, tz = "UTC")
    expect_equal(siderealTime(t), 8.0295394, tolerance = 0.0000001)
})

test_that("Moon", {
    # [2] example 45.a (pages 312-313)
    # Do not check too many digits, because the code does not have all terms
    # in formulae.
    # NB. this block also tests eclipticalToEquatorial(), julianDay(),
    # and julianCenturyAnomaly().
    t <- ISOdatetime(1992, 04, 12, 0, 0, 0, tz = "UTC")
    m <- moonAngle(t, 0, 0) # lat and lon arbitrary
    expect_equal(m$lambda, 133.162659, tolerance = 0.0002)
    expect_equal(m$beta, -3.229127, tolerance = 0.0002)
    ## expect_equal(abs(m$obliquity - 23.440636) < 0.001)
    expect_equal(m$rightAscension, 134.688473, tolerance = 0.02)
    expect_equal(m$declination, 13.768366, tolerance = 0.001)
    expect_equal(m$diameter, 0.991990, tolerance = 0.0001)
    expect_equal(m$distance, 368405.6, tolerance = 0.1)
    # moon illuminated fraction [1] ex 31.b page 156
    illfrac <- (1 + cos(RPD * 105.8493)) / 2
    expect_equal(moonAngle(ISOdatetime(1979, 12, 25, 0, 0, 0, tz = "UTC"), 0, 0)$illuminatedFraction, illfrac, tolerance = 0.001)
    # Local time
    tlocal <- t
    attributes(tlocal)$tzone <- ""
    mlocal <- moonAngle(tlocal, 0, 0)
    expect_identical(m, mlocal)
    # Numerical time
    expect_identical(m, moonAngle(as.numeric(t), 0, 0))
    # time as a Sys.Date value
    moonAngle(Sys.Date(), -63.5744, 44.6479)
    # time as a character value
    moonAngle("2019-01-20", -63.5744, 44.6479)
    # How close to "total full" was the time of the lunar eclipse for the
    # "super Blood full moon" of 2019, in Halifax, NA?
    expect_gt(moonAngle("2019-01-21 00:12:00", -63.5744, 44.6479)$illuminatedFraction, 0.999)
})

test_that("Sun", {
    # Testing against values that worked on 2016-12-06;
    # FIXME: replace by numbers from [2] if they can be found.
    t <- ISOdatetime(1992, 04, 12, 0, 0, 0, tz = "UTC")
    s <- sunAngle(t, 0, 0) # lat and lon arbitrary
    expect_equal(s$azimuth, 358.624168865654)
    expect_equal(s$altitude, -81.3030909018573)
    expect_equal(s$diameter, 0.531871759139675)
    expect_equal(s$distance, 1.00249729533012)
    # Local time
    tlocal <- t
    attributes(tlocal)$tzone <- ""
    slocal <- sunAngle(tlocal, 0, 0)
    expect_identical(s, slocal)
    # Numerical time
    expect_identical(s, sunAngle(as.numeric(t), 0, 0))
    # Character time
    expect_identical(s, sunAngle("1992-04-12 00:00:00", 0, 0))
})

test_that("Sun Declination and Right Ascension", {
    # Example 24.a in Meeus (1991) (page 158 PDF, 153 print)
    # This is *apparent* declination and right ascension
    time <- as.POSIXct("1992-10-13 00:00:00", tz = "UTC")
    a <- sunDeclinationRightAscension(time, apparent = TRUE)
    expect_equal(a$declination, -7.78507, tolerance = 0.00004)
    expect_equal(a$rightAscension, -161.61919, tolerance = 0.00003)
    b <- sunDeclinationRightAscension(time)
    # check against previous results, to protect aginst code-drift errors
    # This is *actual* declination and right ascension
    expect_equal(b$declination, -7.785464443, tolerance = 0.000000001)
    expect_equal(b$rightAscension, -161.6183305, tolerance = 0.0000001)
})

test_that("Handle multiples correctly (issue 2178)", {
    d <- data.frame(
        t = c("2023-01-01 00:00", "2024-01-01 01:00", "2023-01-01 02:00"),
        lat = c(45, 46, 47),
        lon = c(-65, -66, -67),
        ref = rep(TRUE, 3)
    )
    a <- data.frame(with(d, sunAngle(t, lon, lat, ref)))
    for (i in seq_len(nrow(d))) {
        ai <- data.frame(with(d[i, ], sunAngle(t, lon, lat, ref)))
        for (name in names(a)) {
            expect_equal(a[i, name], ai[[name]])
        }
    }
})

test_that("Meeuse Ex 45.a", {
    # Test against example 45a (page 312) in Meeuse (1991). Note
    # that those values are for an updated set of formulae, and
    # this is not what we have in oce (which uses Meeuse 1982
    # formulae.
    #
    # I set up these tests to help me to decide whether to recode
    # oce to the new formulae.  I think the differences are too
    # small to justify that effort, which is quite substantial,
    # involving typing a lot of formulae and numbers in from a book.
    #
    # I am not using expect_equal() because I want to display
    # the mismatch concretely.
    #
    # Summary: the worst differences are under 0.01 degrees,
    # which I judge to be tolerable for simple tasks.
    t <- as.POSIXct("1992-04-12 00:00:00", tz = "UTC")
    ma <- moonAngle(t, longitude = 0, latitude = 0)
    expect_true(abs(ma$rightAscension - 134.688473) < 0.017) # alpha='apparent right ascension'
    expect_true(abs(ma$declination - 13.768366) < 0.0046) # delta='declinaiton'
    expect_true(abs(ma$lambda - 133.162659) < 0.013) # longitude
    expect_true(abs(ma$beta - (-3.229127)) < 0.00011)
    expect_true(abs(ma$distance - 368409.7) < 7.5) # delta='distance [km]'
    expect_true(abs(ma$diameter - 0.991990) < 0.000021) # pi
})