R/love.R

In earthtide: Parallel Implementation of 'ETERNA 3.40' for Prediction and Analysis of Earth Tides

# love
#
# returns the love numbers - direct port of eterna 3.40 code
# http://igets.u-strasbg.fr/soft_and_tool.php
#
# @param latitude latitude the station latitude
# @param elevation elevation the station elevation
#
# @return list of love numbers
#
# @keywords internal
#
love <- function(latitude, elevation) {
  dg0 <- c(
    1.1576, 1.1542, 1.1600, 1.0728, 1.0728, 1.0728,
    1.0728, 1.0363, 1.0363, 1.0363, 1.0363, 1.0363
  )
  dgp <- c(
    -0.0016, -0.0018, -0.0010, 0.0, 0.0, 0.0, -0.0010,
    0.0, 0.0, 0.0, 0.0, -0.000315
  )
  dgm <- c(0.0054, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
  dh0 <- c(
    0.6165, 0.6069, 0.6133, 0.2946, 0.2946, 0.2946,
    0.2946, 0.1807, 0.1807, 0.1807, 0.1807, 0.1807
  )
  dhp <- c(
    0.0007, 0.0007, 0.0005, 0.0, 0.0, 0.0, 0.0003,
    0.0, 0.0, 0.0, 0.0, 0.00015
  )
  dhm <- c(
    0.0018, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
    0.0, 0.0
  )
  dk0 <- c(
    0.3068, 0.3009, 0.3034, 0.0942, 0.0942, 0.0942,
    0.0942, 0.0427, 0.0427, 0.0427, 0.0427, 0.0427
  )
  dkp <- c(
    0.0015, 0.0014, 0.0009, 0.0, 0.0, 0.0, 0.0007,
    0.0, 0.0, 0.0, 0.0, 0.00066
  )
  dkm <- c(-0.0004, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)

  # Shida-numbers
  dl0 <- c(
    0.08400, 0.08410, 0.08520, 0.01490, 0.01490, 0.01490,
    0.01490, 0.01000, 0.01000, 0.01000, 0.01000, 0.01000
  )
  dlp <- c(
    -0.0020, -0.0020, -0.0010, 0.00000, 0.00000, 0.00000,
    0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000
  )
  dlm <- c(
    0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
    0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000
  )
  latitude_p <- c(
    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
    0.0, 0.0
  )
  latitude_m <- c(
    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
    0.0, 0.0
  )

  dea <- 6378136.3
  dee <- 6.69439795140e-3

  domr <- 15.073729
  dom0 <- 13.943036
  dgr <- -0.000625
  dhr <- -0.002505
  dkr <- -0.001261
  dlr <-  0.0000781

  dclat <- cos(latitude * pi / 180)
  dslat <- sin(latitude * pi / 180)
  dn <- dea / sqrt(1.0 - dee * dslat^2) # curvature

  dpsi <- 180 / pi * atan(((dn * (1.0 - dee) + elevation) * dslat) /
    ((dn + elevation) * dclat)) # geo_latitude
  dthet <- 90.0 - dpsi # theta
  dct <- cos(dthet * pi / 180)
  dct2 <- dct * dct

  # 0.335410 = 3/(4*sqrt(5))
  latitude_p[1] <- 0.335410 * (35.0 * dct2 * dct2 - 30.0 * dct2 + 3.0) /
    (3.0 * dct2 - 1.0)  # 2,0 +

  latitude_m[1]  <- 0.894427 / (3.0 * dct2 - 1.0)   # 2,0 -
  latitude_p[2]  <- 0.612372 * (7.0 * dct2 - 3.0)   # 2,1 +
  latitude_p[3]  <- 0.866025 * (7.0 * dct2 - 1.0)   # 2,2 +
  latitude_p[7]  <- 0.829156 * (9.0 * dct2 - 1.0)   # 3,3 +
  latitude_p[12] <- 0.806226 * (11.0 * dct2 - 1.0)  # 4,4 +

  # shida numbers are not reliable - disregard latitude dependence?
  # ET34-X-V71 manual (2018) which suggests 10 as the cutoff
  # For horizontal displacements this number may need to be lower
  # if(any(abs(latitude_p) > 10)) {
  #   latitude_p[] <- 0
  #   latitude_m[] <- 0
  # }

  dglat <- dg0 + dgp * latitude_p + dgm * latitude_m
  dhlat <- dh0 + dhp * latitude_p + dhm * latitude_m
  dklat <- dk0 + dkp * latitude_p + dkm * latitude_m
  dllat <- dl0 + dlp * latitude_p + dlm * latitude_m

  dtlat <- 1.0 + dk0 - dh0 + latitude_p * (dkp - dhp) + latitude_m * (dkm - dhm)
  dtr <- dkr - dhr

  return(list(
    dom0 = dom0, domr = domr, dgr = dgr, dhr = dhr,
    dkr = dkr, dlr = dlr, dtr = dtr, dglat = dglat,
    dhlat = dhlat, dklat = dklat, dllat = dllat,
    dtlat = dtlat
  ))
}