library(acousticTS)

test_that("spherical derivatives match elementary solutions", {
  # The first three orders have closed forms independent of the compiled
  # recurrence. Symbolic derivatives also check complex arguments.
  j <- list(
    quote(sin(z) / z),
    quote(sin(z) / z^2 - cos(z) / z),
    quote((3 / z^3 - 1 / z) * sin(z) - 3 * cos(z) / z^2)
  )
  y <- list(
    quote(-cos(z) / z),
    quote(-cos(z) / z^2 - sin(z) / z),
    quote(-(3 / z^3 - 1 / z) * cos(z) - 3 * sin(z) / z^2)
  )
  for (z in list(
    c(0.8, 1.7, 3.2, 4.1),
    c(
      0.8 + 0.4i,
      1.7 - 0.3i,
      3.2 + 0.2i,
      4.1 - 0.5i
    )
  )) {
    j_expr <- j
    y_expr <- y
    for (k in 0:4) {
      expected_j <- t(vapply(
        j_expr,
        eval,
        complex(length(z)),
        envir = list(z = z)
      ))
      expected_y <- t(vapply(
        y_expr,
        eval,
        complex(length(z)),
        envir = list(z = z)
      ))
      expect_equal(jsdk(0:2, z, k) + 0i, expected_j, tolerance = 1e-09)
      expect_equal(ysdk(0:2, z, k) + 0i, expected_y, tolerance = 1e-09)
      expect_equal(
        as.complex(jsdk(0:2, z[1:3], k)),
        diag(expected_j[, 1:3]),
        tolerance = 1e-09
      )
      expect_equal(
        hsdk(0:2, z, k),
        expected_j + 1i * expected_y,
        tolerance = 1e-09
      )
      for (l in 0:2) {
        expect_equal(
          as.complex(jsdk(l, z, k)),
          expected_j[l + 1, ],
          tolerance = 1e-09
        )
        expect_equal(
          as.complex(ysdk(l, z, k)),
          expected_y[l + 1, ],
          tolerance = 1e-09
        )
        expect_equal(
          hsdk(l, z, k),
          expected_j[l + 1, ] + 1i * expected_y[l + 1, ],
          tolerance = 1e-09
        )
      }
      j_expr <- lapply(j_expr, D, name = "z")
      y_expr <- lapply(y_expr, D, name = "z")
    }
    for (family in list(
      list(jsd, jsdd, jsdk),
      list(ysd, ysdd, ysdk),
      list(hsd, hsdd, hsdk)
    )) {
      expect_equal(family[[1]](0:2, z), family[[3]](0:2, z, 1))
      expect_equal(family[[2]](0:2, z), family[[3]](0:2, z, 2))
      expect_error(family[[3]](1, z, -1), "non-negative")
    }
  }
  for (z in list(0, 0 + 0i)) {
    expect_equal(as.complex(js(0:2, z)), c(1, 0, 0) + 0i)
    expect_equal(as.complex(jsd(0:2, z)), c(0, 1 / 3, 0) + 0i)
    expect_equal(as.complex(jsdd(0:2, z)), c(-1 / 3, 0, 2 / 15) + 0i)
    expect_true(all(!is.finite(ysd(0:2, z))))
    expect_true(all(is.na(hsdd(0:2, z))))
  }
})

test_that("cylindrical functions preserve analytic identities", {
  x <- c(0.8, 2.1, 3.5, 4.2)
  for (l in c(-2, -1, 0, 1, 2, 0.5)) {
    expect_equal(jc(l, -x), exp(1i * pi * l) * besselJ(x, l), tolerance = 1e-10)
    expect_equal(
      yc(l, -x),
      exp(-1i * pi * l) * besselY(x, l) + 2i * cos(pi * l) * besselJ(x, l),
      tolerance = 1e-10
    )
    expect_equal(
      jc(l, 1i * x),
      exp(1i * pi * l / 2) * besselI(x, l),
      tolerance = 1e-10
    )
    expect_equal(jc(l, -1i * x), Conj(jc(l, 1i * x)), tolerance = 1e-10)
    expect_equal(
      yc(l, 1i * x),
      1i * jc(l, 1i * x) -
        (2 / pi) * exp(-1i * pi * l / 2) * besselK(x, l),
      tolerance = 1e-10
    )
    expect_equal(yc(l, -1i * x), Conj(yc(l, 1i * x)), tolerance = 1e-10)
  }
  for (family in list(
    list(jc, jcd, jcdd, jcdk),
    list(yc, ycd, ycdd, ycdk),
    list(hc, hcd, hcdd, hcdk)
  )) {
    for (z in list(x, 1i * x)) {
      expect_equal(family[[4]](0:2, z, 0), family[[1]](0:2, z))
      expect_equal(family[[4]](0:2, z, 1), family[[2]](0:2, z))
      expect_equal(family[[4]](0:2, z, 2), family[[3]](0:2, z))
      for (l in 0:2) {
        # The Bessel differential equation determines the second
        # derivative.
        expect_equal(
          family[[3]](l, z),
          -family[[2]](l, z) / z - (1 - l^2 / z^2) * family[[1]](l, z),
          tolerance = 1e-09
        )
      }
    }
    expect_error(family[[4]](0, 1, -1), "non-negative")
    expect_error(family[[1]](0 + 1i, 1), "order|orders")
    expect_error(family[[1]](0, 1 + 1i), "complex|Complex")
  }
})

test_that("imaginary cylindrical values match independent real references", {
  x <- c(0.8, 2.1, 3.5, 4.2)
  for (l in c(-2, -1, 0, 1, 2, 0.5)) {
    phase <- complex(real = cospi(l / 2), imaginary = sinpi(l / 2))
    expected_j <- phase * besselI(x, l)
    expected_y <- 1i * expected_j - (2 / pi) * Conj(phase) * besselK(x, l)
    for (direction in c(1, -1)) {
      actual_j <- jc(l, direction * 1i * x)
      actual_y <- yc(l, direction * 1i * x)
      reference_j <- if (direction > 0) expected_j else Conj(expected_j)
      reference_y <- if (direction > 0) expected_y else Conj(expected_y)
      expect_true(all(is.finite(actual_j)))
      expect_true(all(is.finite(actual_y)))
      expect_equal(actual_j, reference_j, tolerance = 1e-12)
      expect_equal(actual_y, reference_y, tolerance = 1e-12)
    }
  }
})

test_that("imaginary cylindrical J values match across orders and arguments", {
  x <- c(0.05, 0.8, 2.1, 3.5, 4.2, 20, 50)
  for (l in c(-5.5, -2, -1, -0.5, 0, 0.5, 1, 2, 5.5, 20)) {
    # Preserve exact zero components when multiplying large Bessel values.
    phase <- complex(real = cospi(l / 2), imaginary = sinpi(l / 2))
    expected <- phase * besselI(x, l)
    for (direction in c(1, -1)) {
      actual <- jc(l, direction * 1i * x)
      reference <- if (direction > 0) expected else Conj(expected)
      expect_true(all(is.finite(actual)))
      expect_equal(actual, reference, tolerance = 1e-12)
    }
  }
})

test_that("spherical sequences preserve zero limits and modal order", {
  sequence <- acousticTS:::spherical_bessel_sequence_matrix_cpp
  z <- c(0, 1e-310, 0.2, pi, 4.5)
  for (name in c("js", "jsd", "hs", "hsd")) {
    fun <- get(name, asNamespace("acousticTS"))
    for (orders in list(0L, c(5L, 0L, 2L, 5L))) {
      expected <- matrix(
        unlist(lapply(z, function(x) {
          as.complex(fun(
            orders,
            x
          ))
        })),
        nrow = length(z),
        byrow = TRUE
      )
      expect_equal(sequence(orders, z, name) + 0i, expected, tolerance = 1e-09)
    }
  }
  expect_error(sequence(integer(), 1, "js"), "non-empty")
  expect_error(sequence(0L, numeric(), "js"), "non-empty")
  expect_error(sequence(-1L, 1, "js"), "non-negative")
  expect_error(sequence(0L, 1, "unknown"), "Unsupported")
  for (family in list(list(js, jsdk), list(ys, ysdk), list(hs, hsdk))) {
    for (z in list(0.8, 0.8 + 0.2i)) {
      expect_equal(family[[2]](0:3, z, 0), family[[1]](0:3, z))
    }
  }
})

test_that("Bessel functions work correctly", {
  # Test jc (cylindrical Bessel function of the first kind)
  expect_equal(Re(jc(0, 1)), besselJ(1, 0))
  expect_equal(Re(jc(1, 2)), besselJ(2, 1))
  expect_equal(Re(jc(0.5, 1.5)), besselJ(1.5, 0.5))

  # Test yc (cylindrical Bessel function of the second kind)
  expect_equal(Re(yc(0, 1)), besselY(1, 0))
  expect_equal(Re(yc(1, 2)), besselY(2, 1))
  expect_equal(Re(yc(0.5, 1.5)), besselY(1.5, 0.5))

  # Test matrix input for `yc`
  expect_equal(
    Re(t(yc(1, matrix(c(1, 2, 3))))),
    t(besselY(c(1, 2, 3), 1))
  )

  # Test jcd (first derivative of jc)
  # Test edge case where n = 0 and l = 1
  expect_equal(Re(jcd(1, 0)), 0.5)
  # Test edge case where n = 0 and l != 1
  expect_equal(Re(jcd(2, 0)), 0.0)
  # Test normal case
  l <- 1
  n <- 2
  expected <- jc(l - 1, n) - (l / n) * jc(l, n)
  expect_equal(jcd(l, n), expected)

  # Test jcd with matrix input
  n_l <- seq(1, 4, 1)
  n_z <- seq(1, 6, 1)
  result <- jcd(n_l, n_z)
  expect_true(is.matrix(result))
  expect_equal(dim(result), c(4, 6))

  # Test ycd (first derivative of yc)
  # Test edge case where n = 0 and l = 1
  expect_true(is.na(Re(ycd(1, 0))))
  # Test edge case where n = 0 and l != 1
  expect_true(is.na(Re(ycd(2, 0))))
  # Test normal case
  l <- 1
  n <- 2
  expected <- yc(l - 1, n) - (l / n) * yc(l, n)
  expect_equal(ycd(l, n), expected)

  # Test matrix input for `ycd`
  expect_equal(
    Re(as.vector(ycd(1, matrix(c(1, 2, 3))))),
    besselY(c(1, 2, 3), 0) - (1 / c(1, 2, 3)) * besselY(c(1, 2, 3), 1)
  )

  # Test jcdd (second derivative of jc)
  l <- 1
  n <- 2
  expected <- 0.25 * (jc(l - 2, n) - 2 * jc(l, n) + jc(l + 2, n))
  expect_equal(jcdd(l, n), expected)
})

test_that("Spherical Bessel functions work correctly", {
  # Test js (spherical Bessel function of the first kind)
  # js should be related to jc by: js(l, n) = sqrt(pi/(2*n)) * jc(l + 0.5, n)
  l <- 1
  n <- 2
  expected <- sqrt(pi / (2 * n)) * jc(l + 0.5, n)
  expect_equal(js(l, n), Re(expected), tolerance = 1e-10)

  # Test edge case
  # ---- Case: n == 0
  expect_equal(js(1, 0), 0)

  # Test error case
  # ---- Case: non-numeric 'l'
  expect_error(js("1", 1), "Inputs must be numeric or complex vectors.")
  # ---- Case: non-numeric 'n'
  expect_error(js(1, "1"), "Inputs must be numeric or complex vectors.")
  # ---- Case: non-numeric 'l' and 'n'
  expect_error(js("1", "1"), "Inputs must be numeric or complex vectors.")

  # Test jsd
  # ---- Case: Single value
  expect_equal(
    jsd(1, 1),
    js(0, 1) - (2) / 1 * js(1, 1)
  )
  # ---- Case: Vector[l]
  expect_equal(
    jsd(c(1, 2, 3), 1),
    js(c(1, 2, 3) - 1, 1) - (c(1, 2, 3) + 1) / 1 * js(c(1, 2, 3), 1)
  )
  # ---- Case: Vector[n]
  expect_equal(
    jsd(1, c(1, 2, 3)),
    js(0, c(1, 2, 3)) - (2) / c(1, 2, 3) * js(1, c(1, 2, 3))
  )

  # Test jsdd
  expect_equal(jsdd(1, 1), -0.177098575)

  # Test ys (spherical Bessel function of the second kind)
  # ys should be related to yc by: ys(l, n) = sqrt(pi/(2*n)) * yc(l + 0.5, n)
  expected_ys <- sqrt(pi / (2 * n)) * yc(l + 0.5, n)
  expect_equal(ys(l, n), Re(expected_ys), tolerance = 1e-10)

  # Test edge-cases
  # ---- Case: n < 0
  expect_equal(ys(1, -1), -sin(1) - cos(1))
  # --- Case: n == 0
  expect_equal(ys(1, 0), -Inf)

  # Test error case
  # ---- Case: non-numeric 'l'
  expect_error(ys("1", 1), "Inputs must be numeric or complex vectors.")
  # ---- Case: non-numeric 'n'
  expect_error(ys(1, "1"), "Inputs must be numeric or complex vectors.")
  # ---- Case: non-numeric 'l' and 'n'
  expect_error(ys("1", "1"), "Inputs must be numeric or complex vectors.")

  # Test ysd
  expect_equal(ysd(1, 1), 2 * sin(1) + cos(1))

  # Test ysdd
  expect_equal(ysdd(1, 1), -5 * sin(1) - 3 * cos(1))

  # Test hs (spherical Hankel function)
  # hs should equal js + 1i * ys
  expected_hs <- js(l, n) + 1i * ys(l, n)
  expect_equal(hs(l, n), expected_hs, tolerance = 1e-10)

  # Test hsd
  expect_equal(hsd(1, 1), (2 + 1i) * exp(1i))

  # Test higher order derivative functions (kth derivatives)
  expect_equal(jsdk(1, 1, 1), jsd(1, 1))
  expect_equal(jsdk(1, 1, 2), jsdd(1, 1))
  expect_equal(jsdk(0, 1, 1), jsd(0, 1))
  expect_equal(jsdk(0, 1, 2), jsdd(0, 1))
})

test_that("Spherical Bessel sequence matrix matches scalar API", {
  orders <- 0:28
  z_j <- c(0.15, 1.2, 9)
  seq_j <- acousticTS:::spherical_bessel_sequence_matrix_cpp(orders, z_j, "js")
  ref_j <- t(vapply(z_j, function(zz) js(orders, zz), numeric(length(orders))))
  expect_equal(seq_j, ref_j, tolerance = 1e-10)

  seq_jd <- acousticTS:::spherical_bessel_sequence_matrix_cpp(
    orders, z_j,
    "jsd"
  )
  ref_jd <- t(vapply(
    z_j, function(zz) jsd(orders, zz),
    numeric(length(orders))
  ))
  expect_equal(seq_jd, ref_jd, tolerance = 1e-9)

  z_h <- c(1.2, 9)
  seq_h <- acousticTS:::spherical_bessel_sequence_matrix_cpp(orders, z_h, "hs")
  ref_h <- t(vapply(z_h, function(zz) hs(orders, zz), complex(length(orders))))
  expect_equal(seq_h, ref_h, tolerance = 1e-9)

  seq_hd <- acousticTS:::spherical_bessel_sequence_matrix_cpp(
    orders, z_h,
    "hsd"
  )
  ref_hd <- t(vapply(
    z_h, function(zz) hsd(orders, zz),
    complex(length(orders))
  ))
  expect_equal(seq_hd, ref_hd, tolerance = 1e-8)
})

test_that("Spherical Bessel functions support complex arguments", {
  z <- 1 + 0.5i

  expect_equal(
    js(0, z),
    sin(z) / z,
    tolerance = 1e-12
  )
  expect_equal(
    js(1, z),
    sin(z) / z^2 - cos(z) / z,
    tolerance = 1e-12
  )
  expect_equal(
    ys(0, z),
    -cos(z) / z,
    tolerance = 1e-12
  )
  expect_equal(
    ys(1, z),
    -cos(z) / z^2 - sin(z) / z,
    tolerance = 1e-12
  )
  expect_equal(
    hs(1, z),
    js(1, z) + 1i * ys(1, z),
    tolerance = 1e-12
  )
  expect_equal(
    jsd(0, z),
    -js(1, z),
    tolerance = 1e-12
  )
  expect_equal(
    ysd(0, z),
    -ys(1, z),
    tolerance = 1e-12
  )
  expect_equal(
    hsd(0, z),
    -hs(1, z),
    tolerance = 1e-12
  )
})

test_that("Hankel functions work correctly", {
  # Test hc (cylindrical Hankel function)
  # hc should equal jc + 1i * yc
  l <- 1
  n <- 2
  expected <- jc(l, n) + 1i * yc(l, n)
  expect_equal(hc(l, n), expected, tolerance = 1e-10)

  # Test hcd (first derivative of hc)
  # hcd should equal jcd + 1i * ycd
  expected_hcd <- jcd(l, n) + 1i * ycd(l, n)
  expect_equal(hcd(l, n), expected_hcd, tolerance = 1e-10)

  # Test hcdd
  expect_equal(hcdd(1, 1), hc(1, 1) - hc(0, 1))
  # ---- Case: vector[n]
  expect_equal(
    hcdd(1, c(1, 2, 3)),
    c(
      hc(1, 1) - hc(0, 1),
      0.5 * (-hc(0, 2) - hc(1, 2)),
      1 / 9 * (-3 * hc(0, 3) - 7 * hc(1, 3))
    )
  )
  # ---- Case: matrix[n]
  expect_equal(
    as.vector(hcdd(1, matrix(c(1, 2, 3)))),
    c(
      hc(1, 1) - hc(0, 1),
      0.5 * (-hc(0, 2) - hc(1, 2)),
      1 / 9 * (-3 * hc(0, 3) - 7 * hc(1, 3))
    )
  )
  # ---- Case: n == 0
  expect_true(
    is.na(Re(hcdd(1, 0))) & is.na(Im(hcdd(1, 0)))
  )


  # Test hcdk
  # ---- Case: first derivative
  expect_equal(hcdk(1, 1, 1), hcd(1, 1))
  # ---- Case: second derivative
  expect_equal(hcdk(1, 1, 2), hcdd(1, 1))
  # ---- Case: vector[n]
  expect_equal(hcdk(1, c(1, 2, 3), 1), hcd(1, c(1, 2, 3)))
  # ---- Case: matrix[n]
  expect_equal(
    hcdk(1, matrix(c(1, 2, 3)), 1),
    hcd(1, matrix(c(1, 2, 3)))
  )
})