library(acousticTS)

test_that("prolate TMM geometry and retained blocks preserve exact semiaxes", {
  segments <- c(99L, 100L, 101L)
  objects <- lapply(segments, function(count) {
    fls_generate(
      prolate_spheroid(
        semimajor_length = 0.07,
        semiminor_length = 0.01, n_segments = count
      ),
      g_body = 1.03, h_body = 1.02, theta_body = 37 * pi / 180
    )
  })
  initialized <- lapply(objects, function(object) {
    acousticTS:::tmm_initialize(object, 12000,
      boundary = "fixed_rigid",
      density_sw = 1026.8, sound_speed_sw = 1477.3, store_t_matrix = TRUE
    )
  })
  reference <- initialized[[2]]@model_parameters$TMM
  psms <- acousticTS:::psms_initialize(objects[[2]], 12000,
    boundary = "fixed_rigid", density_sw = 1026.8, sound_speed_sw = 1477.3
  )
  psms_inputs <- psms@model_parameters$PSMS
  expect_identical(reference$body$xi, psms_inputs$body$xi)
  expect_identical(reference$body$q, psms_inputs$body$q)
  columns <- c("chi_sw", "m_max", "n_max")
  expect_identical(
    reference$parameters$acoustics[, columns],
    psms_inputs$parameters$acoustics[, columns]
  )
  for (index in c(1L, 3L)) {
    actual <- initialized[[index]]@model_parameters$TMM
    expect_identical(actual$body, reference$body)
    expect_identical(
      actual$parameters$acoustics,
      reference$parameters$acoustics
    )
  }

  # Compare the same finite retained system across profile discretizations =====
  solved <- lapply(initialized, function(object) acousticTS:::TMM(object))
  expected_blocks <- solved[[2]]@model_parameters$TMM$parameters$t_matrix
  expect_length(expected_blocks, 1L)
  expect_gt(length(expected_blocks[[1]]), 0L)
  expect_true(all(is.finite(solved[[2]]@model$TMM$f_bs)))
  expect_gt(Mod(solved[[2]]@model$TMM$f_bs), 0)
  for (index in c(1L, 3L)) {
    expect_identical(
      solved[[index]]@model$TMM$f_bs,
      solved[[2]]@model$TMM$f_bs
    )
    expect_identical(
      solved[[index]]@model_parameters$TMM$parameters$t_matrix,
      expected_blocks
    )
  }

  # Direct and retained branches must initialize the same physical coordinates =
  direct <- acousticTS:::tmm_initialize(objects[[1]], 12000,
    boundary = "fixed_rigid", density_sw = 1026.8, sound_speed_sw = 1477.3
  )
  expect_identical(direct@model_parameters$TMM$body, reference$body)
  expect_identical(
    direct@model_parameters$TMM$parameters$acoustics,
    reference$parameters$acoustics
  )
  invalid <- objects[[1]]
  invalid@shape_parameters$semiminor_length <- 0.07
  expect_error(acousticTS:::tmm_initialize(invalid, 12000,
    boundary = "fixed_rigid"
  ), "finite scalar semiaxes")
})

test_that("exact spherical TMM preserves invalid-angle handling", {
  for (angle in c(NA_real_, NaN, Inf, -Inf)) {
    amplitude <- acousticTS:::tmm_backscatter_cpp(
      38000, angle, "Sphere", .01, "fixed_rigid", 1477.4, 1026.8, 1030, 1490, 8L
    )
    expect_true(all(is.na(amplitude)))
  }
})

test_that("exact spherical TMM preserves retained blocks and modal cutoffs", {
  frequencies <- c(1000, 38000, 180000)
  cutoffs <- c(1L, 8L, 24L)
  for (boundary in c("fixed_rigid", "pressure_release", "liquid_filled", "gas_filled")) {
    for (angle in c(0, 0.7, pi / 2)) {
      object <- if (boundary == "gas_filled") {
        gas_generate(sphere(.01),
          density_fluid = 1.24, sound_speed_fluid = 345,
          theta_body = angle
        )
      } else {
        fls_generate(sphere(.01),
          g_body = 1.0357, h_body = 1.0279,
          theta_body = angle
        )
      }
      direct <- target_strength(object, frequencies, "tmm",
        boundary = boundary,
        n_max = cutoffs, sound_speed_sw = 1477.4, density_sw = 1026.8
      )
      retained <- target_strength(object, frequencies, "tmm",
        boundary = boundary,
        n_max = cutoffs, sound_speed_sw = 1477.4, density_sw = 1026.8,
        store_t_matrix = TRUE
      )
      expect_equal(direct@model$TMM$f_bs, retained@model$TMM$f_bs, tolerance = 2e-10)
      expect_equal(direct@model$TMM$TS, retained@model$TMM$TS, tolerance = 2e-10)
      expect_identical(direct@model$TMM$n_max, cutoffs)
    }
  }
})

test_that("TMM shell spheres agree with modal and retained solutions", {
  grDevices::pdf(NULL)
  on.exit(grDevices::dev.off(), add = TRUE)
  boundaries <- c(
    "shelled_pressure_release",
    "shelled_liquid",
    "shelled_gas",
    "elastic_shelled"
  )
  for (boundary in boundaries) {
    object <- fixture_sphere(boundary)
    reference_model <- if (boundary == "elastic_shelled") {
      "essms"
    } else {
      "sphms"
    }
    reference <- target_strength(
      object,
      c(12000, 38000),
      model = reference_model,
      boundary = boundary,
      density_sw = 1026.8,
      sound_speed_sw = 1477.3
    )
    fit <- target_strength(
      object,
      c(12000, 38000),
      model = "tmm",
      boundary = boundary,
      density_sw = 1026.8,
      sound_speed_sw = 1477.3,
      store_t_matrix = TRUE
    )
    expect_equal(
      fit@model$TMM$TS,
      reference@model[[toupper(reference_model)]]$TS,
      tolerance = 1e-08
    )
    expect_equal(
      fit@model$TMM$sigma_bs,
      Mod(fit@model$TMM$f_bs)^2,
      tolerance = 1e-12
    )
    expect_true(all(is.finite(fit@model$TMM$TS)))
    expect_error(plot(fit, type = "model"), NA)
    scattered <- tmm_scattering(fit)
    expect_equal(scattered$f_scat, fit@model$TMM$f_bs, tolerance = 1e-08)
    if (boundary %in% c("shelled_liquid", "elastic_shelled")) {
      inferred <- target_strength(
        object, c(12e3, 38e3),
        model = "tmm",
        density_sw = 1026.8, sound_speed_sw = 1477.3
      )
      expect_equal(inferred@model$TMM$TS, fit@model$TMM$TS)
      truncated <- target_strength(
        object, c(12e3, 38e3),
        model = "tmm", n_max = 20L,
        density_sw = 1026.8, sound_speed_sw = 1477.3
      )
      expect_equal(truncated@model$TMM$TS, fit@model$TMM$TS, tolerance = 1e-8)
    }
  }
})

test_that("TMM rejects unsupported shell geometries", {
  object <- fixture_ps("elastic_shelled")
  expect_error(
    target_strength(object, 38e3, "tmm"), "Specify 'boundary' explicitly"
  )
  expect_error(
    target_strength(object, 38e3, "tmm", boundary = "elastic_shelled"),
    "Only the following values for 'boundary'"
  )
  expect_error(
    target_strength(object, 38e3, "tmm", boundary = "liquid_filled"),
    "supports spherical fluid shells"
  )
})

test_that("gas-sphere retained scattering can be plotted", {
  object <- target_strength(
    fixture_sphere("gas_filled"), 38e3, "tmm",
    store_t_matrix = TRUE
  )
  grDevices::pdf(NULL)
  on.exit(grDevices::dev.off(), add = TRUE)
  expect_error(plot(object, type = "scattering"), NA)
})

.with_null_pdf_device <- function(expr) {
  grDevices::pdf(NULL)
  on.exit(grDevices::dev.off(), add = TRUE)
  invisible(force(expr))
}

test_that("TMM matches SPHMS across multiple sphere sizes and boundary cases", {
  density_sw <- 1026.8
  sound_speed_sw <- 1477.3
  frequency <- c(12e3, 38e3, 120e3, 180e3)
  radii <- c(0.005, 0.01, 0.018)

  compare_sphere_case <- function(object, boundary, tolerance = 1e-9) {
    tmm_obj <- target_strength(
      object = object,
      frequency = frequency,
      model = "tmm",
      boundary = boundary,
      density_sw = density_sw,
      sound_speed_sw = sound_speed_sw
    )
    sphms_obj <- target_strength(
      object = object,
      frequency = frequency,
      model = "sphms",
      boundary = boundary,
      density_sw = density_sw,
      sound_speed_sw = sound_speed_sw
    )

    expect_equal(
      tmm_obj@model$TMM$f_bs,
      sphms_obj@model$SPHMS$f_bs,
      tolerance = tolerance
    )
    expect_equal(
      tmm_obj@model$TMM$TS,
      sphms_obj@model$SPHMS$TS,
      tolerance = tolerance
    )
  }

  for (radius in radii) {
    compare_sphere_case(
      fls_generate(
        shape = sphere(radius_body = radius, n_segments = 80),
        g_body = 1,
        h_body = 1
      ),
      boundary = "fixed_rigid"
    )

    compare_sphere_case(
      fls_generate(
        shape = sphere(radius_body = radius, n_segments = 80),
        g_body = 1,
        h_body = 1
      ),
      boundary = "pressure_release"
    )

    compare_sphere_case(
      fls_generate(
        shape = sphere(radius_body = radius, n_segments = 80),
        g_body = 1028.9 / density_sw,
        h_body = 1480.3 / sound_speed_sw
      ),
      boundary = "liquid_filled"
    )

    compare_sphere_case(
      gas_generate(
        shape = sphere(radius_body = radius, n_segments = 80),
        density_fluid = 1.24,
        sound_speed_fluid = 345
      ),
      boundary = "gas_filled"
    )
  }
})

test_that("TMM can retain the frequency-specific T-matrix blocks", {
  density_sw <- 1026.8
  sound_speed_sw <- 1477.3
  object <- fls_generate(
    shape = sphere(radius_body = 0.01, n_segments = 80),
    g_body = 1,
    h_body = 1
  )

  object <- target_strength(
    object = object,
    frequency = 38e3,
    model = "tmm",
    boundary = "pressure_release",
    density_sw = density_sw,
    sound_speed_sw = sound_speed_sw,
    store_t_matrix = TRUE
  )

  t_store <- object@model_parameters$TMM$parameters$t_matrix
  expect_length(t_store, 1)
  expect_true(is.list(t_store[[1]]))
  expect_length(t_store[[1]], object@model$TMM$n_max[1] + 1)
  expect_true(all(vapply(t_store[[1]], function(x) !is.null(x$T), logical(1))))
})

test_that(
  "TMM supports oblate spheroids through the spherical-coordinate branch",
  {
    density_sw <- 1026.8
    sound_speed_sw <- 1477.3
    frequency <- c(12e3, 38e3, 70e3, 120e3)

    # Equal polar and equatorial radii recover the exact sphere response =======
    oblate_sphere <- fls_generate(
      shape = oblate_spheroid(
        length_body = 0.02, radius_body = 0.01,
        n_segments = 80
      ),
      density_body = 1028.9,
      sound_speed_body = 1480.3,
      theta_body = pi / 2
    )
    sphere_object <- fls_generate(
      shape = sphere(radius_body = 0.01, n_segments = 80),
      density_body = 1028.9,
      sound_speed_body = 1480.3,
      theta_body = pi / 2
    )

    oblate_tmm <- target_strength(
      object = oblate_sphere,
      frequency = frequency,
      model = "tmm",
      boundary = "liquid_filled",
      density_sw = density_sw,
      sound_speed_sw = sound_speed_sw,
      store_t_matrix = TRUE
    )
    sphere_sphms <- target_strength(
      object = sphere_object,
      frequency = frequency,
      model = "sphms",
      boundary = "liquid_filled",
      density_sw = density_sw,
      sound_speed_sw = sound_speed_sw
    )

    expect_equal(oblate_tmm@model$TMM$TS, sphere_sphms@model$SPHMS$TS,
      tolerance = 1e-9
    )

    t_store <- oblate_tmm@model_parameters$TMM$parameters$t_matrix
    expect_length(t_store, length(frequency))
    expect_true(all(vapply(t_store, is.list, logical(1))))
    has_t_matrix <- vapply(t_store[[1]], function(x) !is.null(x$T), logical(1))
    expect_true(all(has_t_matrix))
  }
)

test_that(
  "Oblate TMM retained blocks support the same post-processing helpers",
  {
    density_sw <- 1026.8
    sound_speed_sw <- 1477.3
    object <- target_strength(
      object = fls_generate(
        shape = oblate_spheroid(
          length_body = 0.012, radius_body = 0.01,
          n_segments = 80
        ),
        density_body = 1028.9,
        sound_speed_body = 1480.3,
        theta_body = pi / 2
      ),
      frequency = c(38e3, 70e3),
      model = "tmm",
      boundary = "liquid_filled",
      density_sw = density_sw,
      sound_speed_sw = sound_speed_sw,
      store_t_matrix = TRUE
    )

    mono <- tmm_scattering(object)
    avg <- tmm_average_orientation(
      object,
      distribution = tmm_orientation_distribution(
        distribution = "uniform",
        lower = 0.45 * pi,
        upper = pi,
        n_theta = 9
      )
    )
    grid <- tmm_scattering_grid(
      object,
      frequency = 70e3,
      theta_body = pi / 2,
      phi_body = pi,
      theta_scatter = pi / 2,
      phi_scatter = 2 * pi
    )
    summary <- tmm_bistatic_summary(
      object,
      frequency = 70e3,
      n_theta = 21,
      n_phi = 41,
      n_psi = 41
    )
    bundle <- tmm_products(
      object,
      frequency = 70e3,
      orientation = tmm_orientation_distribution(
        distribution = "uniform",
        lower = 0.45 * pi,
        upper = pi,
        n_theta = 9
      ),
      bistatic_summary = TRUE,
      n_theta = 21,
      n_phi = 41,
      n_psi = 41
    )

    expect_true(all(is.finite(mono$sigma_scat)))
    expect_true(all(is.finite(avg$sigma_bs)))
    expect_equal(grid$sigma_scat[1, 1], mono$sigma_scat[2], tolerance = 1e-8)
    expect_equal(
      summary$metrics$sigma_bs,
      mono$sigma_scat[2],
      tolerance = 1e-10
    )
    expect_equal(
      bundle$bistatic_summary$metrics$peak_sigma_scat,
      summary$metrics$peak_sigma_scat,
      tolerance = 1e-10
    )
  }
)

test_that(
  "TMM matches PSMS across multiple prolate spheroids and boundary cases",
  {
    density_sw <- 1026.8
    sound_speed_sw <- 1477.3
    shape_specs <- list(
      list(
        shape = prolate_spheroid(
          length_body = 0.06,
          radius_body = 0.008,
          n_segments = 80
        ),
        frequency = c(12e3, 18e3, 38e3, 70e3, 100e3, 150e3)
      ),
      list(
        shape = prolate_spheroid(
          length_body = 0.14,
          radius_body = 0.01,
          n_segments = 80
        ),
        frequency = c(12e3, 18e3, 38e3, 70e3, 100e3)
      )
    )

    compare_prolate_case <- function(object, boundary, frequency, max_delta) {
      tmm_obj <- target_strength(
        object = object,
        frequency = frequency,
        model = "tmm",
        boundary = boundary,
        density_sw = density_sw,
        sound_speed_sw = sound_speed_sw
      )
      psms_obj <- target_strength(
        object = object,
        frequency = frequency,
        model = "psms",
        boundary = boundary,
        density_sw = density_sw,
        sound_speed_sw = sound_speed_sw,
        precision = if (
          acousticTS:::.quad_precision_available() &&
            boundary %in% c("liquid_filled", "gas_filled")
        ) {
          "quad"
        } else {
          "double"
        },
        simplify_Amn = FALSE
      )

      expect_true(all(is.finite(tmm_obj@model$TMM$TS)))
      expect_true(all(is.finite(tmm_obj@model$TMM$n_max)))
      expect_true(all(tmm_obj@model$TMM$n_max >= 1))
      expect_lt(
        max(abs(tmm_obj@model$TMM$TS - psms_obj@model$PSMS$TS)),
        max_delta
      )
    }

    for (spec in shape_specs) {
      compare_prolate_case(
        fls_generate(shape = spec$shape, g_body = 1, h_body = 1),
        boundary = "fixed_rigid",
        frequency = spec$frequency,
        max_delta = 1e-10
      )

      compare_prolate_case(
        fls_generate(shape = spec$shape, g_body = 1, h_body = 1),
        boundary = "pressure_release",
        frequency = spec$frequency,
        max_delta = 1e-10
      )

      compare_prolate_case(
        fls_generate(
          shape = spec$shape,
          g_body = 1028.9 / density_sw,
          h_body = 1480.3 / sound_speed_sw
        ),
        boundary = "liquid_filled",
        frequency = spec$frequency,
        max_delta = 1e-10
      )

      compare_prolate_case(
        gas_generate(
          shape = spec$shape,
          density_fluid = 1.24,
          sound_speed_fluid = 345
        ),
        boundary = "gas_filled",
        frequency = spec$frequency,
        max_delta = 1e-10
      )
    }
  }
)

test_that("TMM warns when prolate spheroids request spherical-branch options", {
  density_sw <- 1026.8
  sound_speed_sw <- 1477.3
  object <- fls_generate(
    shape = prolate_spheroid(
      length_body = 0.14, radius_body = 0.01,
      n_segments = 80
    ),
    g_body = 1,
    h_body = 1
  )

  expect_warning(
    target_strength(
      object,
      frequency = 38e3,
      model = "tmm",
      boundary = "fixed_rigid",
      density_sw = density_sw,
      sound_speed_sw = sound_speed_sw,
      n_max = 30,
      store_t_matrix = TRUE
    ),
    "n_max' is ignored"
  )
})

test_that(
  "TMM can retain the frequency-specific T-matrix blocks for prolates",
  {
    density_sw <- 1026.8
    sound_speed_sw <- 1477.3
    object <- fls_generate(
      shape = prolate_spheroid(
        length_body = 0.14, radius_body = 0.01,
        n_segments = 80
      ),
      density_body = 1028.9,
      sound_speed_body = 1480.3,
      theta_body = pi / 2
    )

    object <- target_strength(
      object = object,
      frequency = c(38e3, 100e3),
      model = "tmm",
      boundary = "liquid_filled",
      density_sw = density_sw,
      sound_speed_sw = sound_speed_sw,
      store_t_matrix = TRUE
    )

    t_store <- object@model_parameters$TMM$parameters$t_matrix
    expect_length(t_store, 2)
    expect_true(all(vapply(t_store, is.list, logical(1))))
    has_t_matrix <- vapply(t_store[[1]], function(x) !is.null(x$T), logical(1))
    expect_true(all(has_t_matrix))
  }
)

test_that("Stored prolate TMM keeps the spheroidal modal cutoff", {
  density_sw <- 1026.8
  sound_speed_sw <- 1480
  object <- fls_generate(
    shape = prolate_spheroid(
      length_body = 0.07, radius_body = 0.01,
      n_segments = 80
    ),
    density_body = 1024,
    sound_speed_body = 1480,
    theta_body = pi / 2
  )

  object <- target_strength(
    object = object,
    frequency = 38e3,
    model = "tmm",
    boundary = "pressure_release",
    density_sw = density_sw,
    sound_speed_sw = sound_speed_sw,
    store_t_matrix = TRUE
  )

  acoustics <- object@model_parameters$TMM$parameters$acoustics
  expected_n_max <- acoustics$m_max + ceiling(0.5 * acoustics$chi_sw)

  expect_equal(acoustics$n_max, expected_n_max)
})

test_that(
  paste0(
    "Stored prolate TMM paired-point evaluation matches scalar spheroidal ",
    "calls"
  ),
  {
    density_sw <- 1026.8
    sound_speed_sw <- 1480
    theta_body <- pi / 2
    phi_body <- pi / 3
    object <- target_strength(
      object = fls_generate(
        shape = prolate_spheroid(
          length_body = 0.07, radius_body = 0.01,
          n_segments = 80
        ),
        g_body = 1,
        h_body = 1,
        theta_body = theta_body
      ),
      frequency = 38e3,
      model = "tmm",
      boundary = "pressure_release",
      density_sw = density_sw,
      sound_speed_sw = sound_speed_sw,
      store_t_matrix = TRUE
    )

    model_params <- object@model_parameters$TMM
    acoustics <- model_params$parameters$acoustics[1, , drop = FALSE]
    t_store <- model_params$parameters$t_matrix[[1]]
    theta_scatter <- c(pi / 3, pi / 2, 2 * pi / 3)
    phi_scatter <- c(0, pi / 2, 5 * pi / 4)

    paired <- acousticTS:::.tmm_scattering_points(
      model_params = model_params,
      frequency_idx = 1L,
      shape_parameters = acousticTS::extract(object, "shape_parameters"),
      theta_body = rep(theta_body, length(theta_scatter)),
      phi_body = rep(phi_body, length(theta_scatter)),
      theta_scatter = theta_scatter,
      phi_scatter = phi_scatter
    )
    scalar <- vapply(
      seq_along(theta_scatter),
      function(i) {
        acousticTS:::.tmm_scattering_spheroidal(
          t_store = list(t_store),
          acoustics = acoustics,
          parameters = model_params$parameters,
          theta_body = theta_body,
          phi_body = phi_body,
          theta_scatter = theta_scatter[i],
          phi_scatter = phi_scatter[i]
        )[1]
      },
      complex(1)
    )

    expect_equal(paired, scalar, tolerance = 1e-10)
  }
)

test_that("tmm_scattering reproduces stored monostatic results", {
  density_sw <- 1026.8
  sound_speed_sw <- 1477.3

  sphere_object <- target_strength(
    object = fls_generate(
      shape = sphere(radius_body = 0.01, n_segments = 80),
      density_body = 1028.9,
      sound_speed_body = 1480.3,
      theta_body = pi / 2
    ),
    frequency = c(38e3, 70e3),
    model = "tmm",
    boundary = "liquid_filled",
    density_sw = density_sw,
    sound_speed_sw = sound_speed_sw,
    store_t_matrix = TRUE
  )
  sphere_scat <- tmm_scattering(sphere_object)
  expect_true(all(is.finite(sphere_scat$sigma_scat)))

  prolate_object <- target_strength(
    object = fls_generate(
      shape = prolate_spheroid(
        length_body = 0.14, radius_body = 0.01,
        n_segments = 80
      ),
      density_body = 1028.9,
      sound_speed_body = 1480.3,
      theta_body = pi / 2
    ),
    frequency = c(38e3, 70e3),
    model = "tmm",
    boundary = "liquid_filled",
    density_sw = density_sw,
    sound_speed_sw = sound_speed_sw,
    store_t_matrix = TRUE
  )
  prolate_scat <- tmm_scattering(prolate_object)
  expect_equal(prolate_scat$f_scat, prolate_object@model$TMM$f_bs,
    tolerance =
      1e-10
  )
})

test_that(
  paste0(
    "Stored and unstored prolate TMM runs agree and stay tied to the PSMS ",
    "benchmark"
  ),
  {
    density_sw <- 1026.8
    sound_speed_sw <- 1477.3
    frequency <- c(38e3, 70e3, 100e3)
    object <- fls_generate(
      shape = prolate_spheroid(
        length_body = 0.14, radius_body = 0.01,
        n_segments = 80
      ),
      density_body = 1028.9,
      sound_speed_body = 1480.3,
      theta_body = pi / 2
    )

    tmm_direct <- target_strength(
      object = object,
      frequency = frequency,
      model = "tmm",
      boundary = "liquid_filled",
      density_sw = density_sw,
      sound_speed_sw = sound_speed_sw,
      store_t_matrix = FALSE
    )
    tmm_stored <- target_strength(
      object = object,
      frequency = frequency,
      model = "tmm",
      boundary = "liquid_filled",
      density_sw = density_sw,
      sound_speed_sw = sound_speed_sw,
      store_t_matrix = TRUE
    )
    psms_obj <- target_strength(
      object = object,
      frequency = frequency,
      model = "psms",
      boundary = "liquid_filled",
      density_sw = density_sw,
      sound_speed_sw = sound_speed_sw,
      precision = if (acousticTS:::.quad_precision_available()) {
        "quad"
      } else {
        "double"
      },
      simplify_Amn = FALSE
    )
    expected_f_bs <- (-2i / wavenumber(frequency, sound_speed_sw)) *
      psms_obj@model$PSMS$f_bs

    expect_lt(max(Mod(tmm_direct@model$TMM$f_bs - expected_f_bs)), 1e-6)
    expect_lt(max(Mod(tmm_stored@model$TMM$f_bs - expected_f_bs)), 1e-6)
    expect_lt(
      max(Mod(tmm_stored@model$TMM$f_bs - tmm_direct@model$TMM$f_bs)),
      1e-6
    )
    expect_equal(tmm_stored@model$TMM$sigma_bs, tmm_direct@model$TMM$sigma_bs,
      tolerance = 1e-6
    )
    expect_equal(tmm_stored@model$TMM$TS, tmm_direct@model$TMM$TS,
      tolerance =
        5e-4
    )
  }
)

tmm_finite_rigid_release_field <- function(acoustics, body, boundary) {
  vapply(seq_len(nrow(acoustics)), function(index) {
    acoustic <- acoustics[index, , drop = FALSE]
    orders <- seq.int(0L, acoustic$m_max)
    m <- rep(orders, times = acoustic$n_max - orders + 1L)
    n <- unlist(lapply(orders, function(order) seq.int(order, acoustic$n_max)))
    # Evaluate the finite separated series independently through public APIs ===
    incident <- Smn(m, n, acoustic$chi_sw, cos(body$theta_body),
      normalize = TRUE
    )$value
    receiving <- Smn(m, n, acoustic$chi_sw, cos(body$theta_scatter),
      normalize = TRUE
    )$value
    first <- Rmn(m, n, acoustic$chi_sw, body$xi, kind = 1)
    outgoing <- Rmn(m, n, acoustic$chi_sw, body$xi, kind = 3)
    component <- if (boundary == "fixed_rigid") "derivative" else "value"
    coefficient <- -diag(first[[component]]) / diag(outgoing[[component]])
    multiplicity <- ifelse(m == 0L, 1, 2)
    (-2i / acoustic$k_sw) * sum(multiplicity * coefficient * incident *
      receiving * cos(m * (body$phi_scatter - body$phi_body)))
  }, complex(1))
}

tmm_rebuild_prolate_support <- function(object, extra_orders, extra_degrees) {
  # Rebuild test objects internally since public n_max is ignored for prolates =
  acoustics <- object@model_parameters$TMM$parameters$acoustics
  acoustics$m_max <- acoustics$m_max + extra_orders
  acoustics$n_max <- acoustics$n_max + extra_degrees
  object@model_parameters$TMM$parameters$acoustics <- acoustics
  acousticTS:::TMM(object)
}

test_that(
  "Stored prolate pressure release preserves finite support and converges",
  {
    density_sw <- 1026.8
    sound_speed_sw <- 1480
    theta_body <- pi / 2
    phi_body <- pi / 2
    scatter_angles <- data.frame(
      theta_scatter = c(pi / 2, pi / 2, pi / 2, pi / 3, 2 * pi / 3),
      phi_scatter = c(0, pi / 2, pi, pi / 2, 3 * pi / 2)
    )

    object <- target_strength(
      object = fls_generate(
        shape = prolate_spheroid(
          length_body = 0.07, radius_body = 0.01,
          n_segments = 80
        ),
        g_body = 1,
        h_body = 1,
        theta_body = theta_body
      ),
      frequency = c(38e3, 70e3),
      model = "tmm",
      boundary = "pressure_release",
      density_sw = density_sw,
      sound_speed_sw = sound_speed_sw,
      store_t_matrix = TRUE
    )

    acoustics <- object@model_parameters$TMM$parameters$acoustics
    body <- object@model_parameters$TMM$body
    medium <- object@model_parameters$TMM$medium
    quad <- gauss_legendre(n = 96, a = -1, b = 1)
    enlarged <- list(
      tmm_rebuild_prolate_support(object, 8L, 8L),
      tmm_rebuild_prolate_support(object, 8L, 16L),
      tmm_rebuild_prolate_support(object, 16L, 16L)
    )

    reference_acoustics <-
      enlarged[[length(enlarged)]]@model_parameters$TMM$parameters$acoustics

    for (i in seq_len(nrow(scatter_angles))) {
      body$theta_body <- theta_body
      body$phi_body <- phi_body
      body$theta_scatter <- scatter_angles$theta_scatter[i]
      body$phi_scatter <- scatter_angles$phi_scatter[i]

      exact_raw <- prolate_spheroid_fbs(
        acoustics = reference_acoustics,
        body = body,
        medium = medium,
        integration_pts = quad,
        precision = "double",
        Amn_method = "Amn_pressure_release",
        adaptive = FALSE
      )
      exact_f <- (-2i / acoustics$k_sw) * exact_raw
      stored_f <- tmm_scattering(
        object,
        theta_body = theta_body,
        phi_body = phi_body,
        theta_scatter = body$theta_scatter,
        phi_scatter = body$phi_scatter
      )$f_scat

      finite_f <- tmm_finite_rigid_release_field(acoustics, body, "pressure_release")
      expect_equal(stored_f, finite_f, tolerance = 1e-10)
      # Compare enlarged fields against the same fixed-support PSMS reference ==
      previous <- NULL
      for (fitted in enlarged) {
        enlarged_f <- tmm_scattering(fitted,
          theta_body = theta_body, phi_body = phi_body,
          theta_scatter = body$theta_scatter,
          phi_scatter = body$phi_scatter
        )$f_scat
        if (!is.null(previous)) {
          expect_equal(enlarged_f, previous, tolerance = 1e-10)
        }
        previous <- enlarged_f
      }
      expect_equal(previous, exact_f, tolerance = 1e-10)
    }
  }
)

test_that("Stored prolate rigid scattering preserves finite support and converges", {
  density_sw <- 1026.8
  sound_speed_sw <- 1480
  theta_body <- pi / 2
  phi_body <- pi / 2
  scatter_angles <- data.frame(
    theta_scatter = c(pi / 2, pi / 2, pi / 2, pi / 3, 2 * pi / 3),
    phi_scatter = c(0, pi / 2, pi, pi / 2, 3 * pi / 2)
  )

  object <- target_strength(
    object = fls_generate(
      shape = prolate_spheroid(
        length_body = 0.07, radius_body = 0.01,
        n_segments = 80
      ),
      g_body = 1,
      h_body = 1,
      theta_body = theta_body
    ),
    frequency = c(38e3, 70e3),
    model = "tmm",
    boundary = "fixed_rigid",
    density_sw = density_sw,
    sound_speed_sw = sound_speed_sw,
    store_t_matrix = TRUE
  )

  acoustics <- object@model_parameters$TMM$parameters$acoustics
  body <- object@model_parameters$TMM$body
  medium <- object@model_parameters$TMM$medium
  quad <- gauss_legendre(n = 96, a = -1, b = 1)
  enlarged <- list(
    tmm_rebuild_prolate_support(object, 8L, 8L),
    tmm_rebuild_prolate_support(object, 8L, 16L),
    tmm_rebuild_prolate_support(object, 16L, 16L)
  )

  reference_acoustics <-
    enlarged[[length(enlarged)]]@model_parameters$TMM$parameters$acoustics

  for (i in seq_len(nrow(scatter_angles))) {
    body$theta_body <- theta_body
    body$phi_body <- phi_body
    body$theta_scatter <- scatter_angles$theta_scatter[i]
    body$phi_scatter <- scatter_angles$phi_scatter[i]

    exact_raw <- prolate_spheroid_fbs(
      acoustics = reference_acoustics,
      body = body,
      medium = medium,
      integration_pts = quad,
      precision = "double",
      Amn_method = "Amn_fixed_rigid",
      adaptive = FALSE
    )
    exact_f <- (-2i / acoustics$k_sw) * exact_raw
    stored_f <- tmm_scattering(
      object,
      theta_body = theta_body,
      phi_body = phi_body,
      theta_scatter = body$theta_scatter,
      phi_scatter = body$phi_scatter
    )$f_scat

    finite_f <- tmm_finite_rigid_release_field(acoustics, body, "fixed_rigid")
    expect_equal(stored_f, finite_f, tolerance = 1e-10)
    # Compare enlarged fields against the same fixed-support PSMS reference ====
    previous <- NULL
    for (fitted in enlarged) {
      enlarged_f <- tmm_scattering(fitted,
        theta_body = theta_body, phi_body = phi_body,
        theta_scatter = body$theta_scatter,
        phi_scatter = body$phi_scatter
      )$f_scat
      if (!is.null(previous)) {
        expect_equal(enlarged_f, previous, tolerance = 1e-10)
      }
      previous <- enlarged_f
    }
    expect_equal(previous, exact_f, tolerance = 1e-10)
  }
})

test_that("tmm_average_orientation reduces to the supplied single angle", {
  density_sw <- 1026.8
  sound_speed_sw <- 1477.3
  object <- target_strength(
    object = fls_generate(
      shape = prolate_spheroid(
        length_body = 0.14, radius_body = 0.01,
        n_segments = 80
      ),
      density_body = 1028.9,
      sound_speed_body = 1480.3,
      theta_body = pi / 2
    ),
    frequency = c(38e3, 70e3),
    model = "tmm",
    boundary = "liquid_filled",
    density_sw = density_sw,
    sound_speed_sw = sound_speed_sw,
    store_t_matrix = TRUE
  )

  one_angle <- tmm_scattering(object, theta_body = pi / 2)
  avg_angle <- tmm_average_orientation(object, theta_body = pi / 2)
  expect_equal(avg_angle$sigma_bs, one_angle$sigma_scat, tolerance = 1e-10)
  expect_equal(avg_angle$TS, db(avg_angle$sigma_bs), tolerance = 1e-12)
})

test_that(
  "TMM orientation distributions normalize and agree with explicit averaging",
  {
    density_sw <- 1026.8
    sound_speed_sw <- 1477.3
    object <- target_strength(
      object = fls_generate(
        shape = prolate_spheroid(
          length_body = 0.14, radius_body = 0.01,
          n_segments = 80
        ),
        density_body = 1028.9,
        sound_speed_body = 1480.3,
        theta_body = pi / 2
      ),
      frequency = c(38e3, 70e3),
      model = "tmm",
      boundary = "liquid_filled",
      density_sw = density_sw,
      sound_speed_sw = sound_speed_sw,
      store_t_matrix = TRUE
    )

    quadrature_dist <- tmm_orientation_distribution(
      distribution = "quadrature",
      theta_body = c(0.6 * pi, 0.75 * pi, pi),
      weights = c(1, 2, 1),
      phi_body = pi
    )
    expect_equal(sum(quadrature_dist$weights), 1, tolerance = 1e-12)

    avg_from_dist <- tmm_average_orientation(object,
      distribution =
        quadrature_dist
    )
    avg_explicit <- tmm_average_orientation(
      object,
      theta_body = quadrature_dist$theta_body,
      weights = quadrature_dist$weights,
      phi_body = quadrature_dist$phi_body
    )
    expect_equal(
      avg_from_dist$sigma_bs,
      avg_explicit$sigma_bs,
      tolerance = 1e-10
    )
    expect_equal(avg_from_dist$TS, avg_explicit$TS, tolerance = 1e-10)

    pdf_dist <- tmm_orientation_distribution(
      distribution = "pdf",
      theta_body = seq(0.5 * pi, pi, length.out = 9),
      pdf = function(theta) sin(theta)^2
    )
    expect_equal(sum(pdf_dist$weights), 1, tolerance = 1e-12)
    expect_true(all(pdf_dist$weights >= 0))

    trunc_dist <- tmm_orientation_distribution(
      distribution = "truncated_normal",
      mean_theta = 0.75 * pi,
      sd_theta = 0.08 * pi,
      lower = 0.5 * pi,
      upper = pi,
      n_theta = 11
    )
    expect_equal(sum(trunc_dist$weights), 1, tolerance = 1e-12)
    expect_true(all(trunc_dist$theta_body >= 0.5 * pi))
    expect_true(all(trunc_dist$theta_body <= pi))
  }
)

test_that("TMM scattering plots work from stored blocks", {
  density_sw <- 1026.8
  sound_speed_sw <- 1477.3

  sphere_object <- target_strength(
    object = fls_generate(
      shape = sphere(radius_body = 0.01, n_segments = 80),
      density_body = 1028.9,
      sound_speed_body = 1480.3,
      theta_body = pi / 2
    ),
    frequency = c(38e3, 70e3),
    model = "tmm",
    boundary = "liquid_filled",
    density_sw = density_sw,
    sound_speed_sw = sound_speed_sw,
    store_t_matrix = TRUE
  )
  expect_error(
    .with_null_pdf_device(
      plot(
        sphere_object,
        type = "scattering",
        frequency = 38e3,
        vary = "theta_scatter"
      )
    ),
    NA
  )
  expect_error(
    .with_null_pdf_device(
      plot(
        sphere_object,
        type = "scattering",
        frequency = 38e3,
        heatmap = TRUE,
        n_theta = 31,
        n_phi = 61
      )
    ),
    NA
  )

  prolate_object <- target_strength(
    object = fls_generate(
      shape = prolate_spheroid(
        length_body = 0.14, radius_body = 0.01,
        n_segments = 80
      ),
      density_body = 1028.9,
      sound_speed_body = 1480.3,
      theta_body = pi / 2
    ),
    frequency = c(38e3, 70e3),
    model = "tmm",
    boundary = "liquid_filled",
    density_sw = density_sw,
    sound_speed_sw = sound_speed_sw,
    store_t_matrix = TRUE
  )
  expect_error(
    .with_null_pdf_device(
      plot(
        prolate_object,
        type = "scattering",
        frequency = 70e3,
        vary = "phi_scatter",
        quantity = "sigma_scat"
      )
    ),
    NA
  )
  expect_error(
    .with_null_pdf_device(
      plot(
        prolate_object,
        type = "scattering",
        frequency = 70e3,
        polar = TRUE,
        n_theta = 31,
        n_phi = 61
      )
    ),
    NA
  )
})

test_that("TMM scattering plots require stored blocks", {
  density_sw <- 1026.8
  sound_speed_sw <- 1477.3

  object <- target_strength(
    object = fls_generate(
      shape = sphere(radius_body = 0.01, n_segments = 80),
      density_body = 1028.9,
      sound_speed_body = 1480.3,
      theta_body = pi / 2
    ),
    frequency = 38e3,
    model = "tmm",
    boundary = "liquid_filled",
    density_sw = density_sw,
    sound_speed_sw = sound_speed_sw
  )

  expect_error(
    .with_null_pdf_device(
      plot(object, type = "scattering", frequency = 38e3)
    ),
    "Stored T-matrix blocks are required"
  )
})

test_that("prolate TMM angular batches preserve pointwise complex fields", {
  precisions <- c("double", if (acousticTS:::.quad_precision_available()) "quad")
  theta <- c(0, 0.3 * pi, 0.7 * pi, pi)
  phi <- c(0.2, 1.4, 2.8)
  for (precision in precisions) {
    object <- target_strength(
      fls_generate(
        prolate_spheroid(length_body = 0.14, radius_body = 0.01),
        density_body = 1028.9, sound_speed_body = 1480.3,
        theta_body = 1.1
      ),
      38e3, "tmm",
      boundary = "liquid_filled", precision = precision,
      density_sw = 1026.8, sound_speed_sw = 1477.3, store_t_matrix = TRUE
    )
    grid <- tmm_scattering_grid(
      object,
      theta_body = 1.1, phi_body = 0.4,
      theta_scatter = theta, phi_scatter = phi
    )
    expected <- matrix(0i, length(theta), length(phi))
    for (i in seq_along(theta)) {
      for (j in seq_along(phi)) {
        point <- tmm_scattering(
          object,
          theta_body = 1.1, phi_body = 0.4,
          theta_scatter = theta[i], phi_scatter = phi[j]
        )
        expected[i, j] <- point$f_scat
      }
    }
    expect_equal(grid$f_scat, expected, tolerance = 1e-12)
  }
})

test_that("tmm_scattering_grid reproduces stored monostatic scattering", {
  density_sw <- 1026.8
  sound_speed_sw <- 1477.3

  object <- target_strength(
    object = fls_generate(
      shape = prolate_spheroid(
        length_body = 0.14, radius_body = 0.01,
        n_segments = 80
      ),
      density_body = 1028.9,
      sound_speed_body = 1480.3,
      theta_body = pi / 2
    ),
    frequency = c(38e3, 70e3),
    model = "tmm",
    boundary = "liquid_filled",
    density_sw = density_sw,
    sound_speed_sw = sound_speed_sw,
    store_t_matrix = TRUE
  )

  mono <- tmm_scattering(
    object,
    theta_body = pi / 2,
    phi_body = pi,
    theta_scatter = pi / 2,
    phi_scatter = 2 * pi
  )
  grid <- tmm_scattering_grid(
    object,
    frequency = 70e3,
    theta_body = pi / 2,
    phi_body = pi,
    theta_scatter = pi / 2,
    phi_scatter = 2 * pi
  )

  expect_equal(grid$frequency, 70e3)
  expect_equal(grid$sigma_scat[1, 1], mono$sigma_scat[2], tolerance = 1e-10)
  expect_equal(grid$sigma_scat_dB[1, 1], mono$sigma_scat_dB[2],
    tolerance =
      1e-10
  )
})

test_that(
  paste0(
    "stored spherical TMM point scattering stays on the monostatic amplitude ",
    "scale"
  ),
  {
    density_sw <- 1026.8
    sound_speed_sw <- 1477.3

    object <- target_strength(
      object = fls_generate(
        shape = sphere(radius_body = 0.01, n_segments = 80),
        g_body = 1,
        h_body = 1,
        theta_body = pi / 2
      ),
      frequency = c(38e3, 70e3),
      model = "tmm",
      boundary = "pressure_release",
      density_sw = density_sw,
      sound_speed_sw = sound_speed_sw,
      store_t_matrix = TRUE
    )

    mono <- tmm_scattering(
      object,
      theta_body = pi / 2,
      phi_body = pi / 2,
      theta_scatter = pi / 2,
      phi_scatter = 3 * pi / 2
    )
    grid <- tmm_scattering_grid(
      object,
      frequency = 70e3,
      theta_body = pi / 2,
      phi_body = pi / 2,
      theta_scatter = pi / 2,
      phi_scatter = 3 * pi / 2
    )

    expect_equal(
      mono$f_scat[2],
      object@model$TMM$f_bs[2],
      tolerance = 1e-10
    )
    expect_equal(
      mono$sigma_scat[2],
      object@model$TMM$sigma_bs[2],
      tolerance = 1e-10
    )
    expect_equal(
      grid$f_scat[1, 1],
      mono$f_scat[2],
      tolerance = 1e-10
    )
  }
)

test_that(
  "TMM bistatic summaries and product bundles are internally consistent",
  {
    density_sw <- 1026.8
    sound_speed_sw <- 1477.3

    object <- target_strength(
      object = fls_generate(
        shape = prolate_spheroid(
          length_body = 0.14, radius_body = 0.01,
          n_segments = 80
        ),
        density_body = 1028.9,
        sound_speed_body = 1480.3,
        theta_body = pi / 2
      ),
      frequency = c(38e3, 70e3, 100e3),
      model = "tmm",
      boundary = "liquid_filled",
      density_sw = density_sw,
      sound_speed_sw = sound_speed_sw,
      store_t_matrix = TRUE
    )

    summary_70 <- tmm_bistatic_summary(
      object,
      frequency = 70e3,
      n_theta = 21,
      n_phi = 41,
      n_psi = 41,
      include_grid = TRUE
    )
    monostatic_70 <- tmm_scattering(
      object,
      theta_body = pi / 2,
      phi_body = pi,
      theta_scatter = pi / 2,
      phi_scatter = 2 * pi
    )

    expect_equal(summary_70$metrics$sigma_bs, monostatic_70$sigma_scat[2],
      tolerance = 1e-10
    )
    expect_equal(summary_70$metrics$TS, db(summary_70$metrics$sigma_bs),
      tolerance = 1e-12
    )
    expect_equal(
      summary_70$metrics$peak_sigma_scat,
      max(summary_70$grid$sigma_scat),
      tolerance = 1e-10
    )
    expect_true(is.finite(summary_70$metrics$backscatter_lobe_width))
    expect_true(all(summary_70$sector_integrals$integrated_sigma_scat >= 0))

    orientation <- tmm_orientation_distribution(
      distribution = "uniform",
      lower = 0.5 * pi,
      upper = pi,
      n_theta = 9
    )
    products <- tmm_products(
      object,
      frequency = 70e3,
      orientation = orientation,
      bistatic_summary = TRUE,
      n_theta = 21,
      n_phi = 41,
      n_psi = 41
    )

    expect_equal(products$monostatic$f_scat, tmm_scattering(object)$f_scat,
      tolerance = 1e-10
    )
    expect_equal(
      products$orientation_average$sigma_bs,
      tmm_average_orientation(object, distribution = orientation)$sigma_bs,
      tolerance = 1e-10
    )
    expect_equal(
      products$bistatic_summary$metrics$peak_sigma_scat,
      summary_70$metrics$peak_sigma_scat,
      tolerance = 1e-10
    )
  }
)

test_that("Default cylinder TMM matches FCMS through the cylindrical backend", {
  old_opt <- options(acousticTS.warn_tmm_cylinder = FALSE)
  on.exit(options(old_opt), add = TRUE)
  density_sw <- 1026.8
  sound_speed_sw <- 1477.3
  frequency <- c(12e3, 18e3, 38e3, 70e3, 100e3, 150e3, 200e3)

  compare_cylinder_case <- function(object, boundary, tolerance = 1e-12) {
    tmm_obj <- target_strength(
      object = object,
      frequency = frequency,
      model = "tmm",
      boundary = boundary,
      density_sw = density_sw,
      sound_speed_sw = sound_speed_sw
    )
    fcms_obj <- target_strength(
      object = object,
      frequency = frequency,
      model = "fcms",
      boundary = boundary,
      density_sw = density_sw,
      sound_speed_sw = sound_speed_sw
    )

    expect_true(all(is.finite(tmm_obj@model$TMM$TS)))
    expect_true(
      all(
        tmm_obj@model_parameters$TMM$parameters$coordinate_system ==
          "cylindrical"
      )
    )
    expect_equal(tmm_obj@model$TMM$f_bs, fcms_obj@model$FCMS$f_bs,
      tolerance =
        tolerance
    )
    expect_equal(tmm_obj@model$TMM$TS, fcms_obj@model$FCMS$TS,
      tolerance =
        tolerance
    )
  }

  cylinder_fls <- fls_generate(
    shape = cylinder(length_body = 0.07, radius_body = 0.01, n_segments = 80),
    density_body = 1028.9,
    sound_speed_body = 1480.3,
    theta_body = pi / 2
  )
  cylinder_gas <- gas_generate(
    shape = cylinder(length_body = 0.07, radius_body = 0.01, n_segments = 80),
    density_fluid = 1.24,
    sound_speed_fluid = 345,
    theta_body = pi / 2
  )

  compare_cylinder_case(cylinder_fls, "fixed_rigid")
  compare_cylinder_case(cylinder_fls, "pressure_release")
  compare_cylinder_case(cylinder_fls, "liquid_filled")
  compare_cylinder_case(cylinder_gas, "gas_filled")
})

test_that("Cylinder TMM warns that the branch is still experimental", {
  density_sw <- 1026.8
  sound_speed_sw <- 1477.3

  expect_warning(
    target_strength(
      object = fls_generate(
        shape = cylinder(
          length_body = 0.07, radius_body = 0.01, n_segments =
            80
        ),
        density_body = 1028.9,
        sound_speed_body = 1480.3,
        theta_body = pi / 2
      ),
      frequency = 38e3,
      model = "tmm",
      boundary = "pressure_release",
      density_sw = density_sw,
      sound_speed_sw = sound_speed_sw
    ),
    "Cylinder 'TMM' support remains experimental"
  )
})

test_that(
  "Cylinder TMM stores cylindrical-family retained state when requested",
  {
    old_opt <- options(acousticTS.warn_tmm_cylinder = FALSE)
    on.exit(options(old_opt), add = TRUE)
    density_sw <- 1026.8
    sound_speed_sw <- 1477.3
    object <- target_strength(
      object = fls_generate(
        shape = cylinder(
          length_body = 0.07,
          radius_body = 0.01,
          n_segments = 80
        ),
        density_body = 1028.9,
        sound_speed_body = 1480.3,
        theta_body = pi / 2
      ),
      frequency = c(38e3, 70e3),
      model = "tmm",
      boundary = "fixed_rigid",
      density_sw = density_sw,
      sound_speed_sw = sound_speed_sw,
      store_t_matrix = TRUE
    )

    t_store <- object@model_parameters$TMM$parameters$t_matrix

    expect_identical(
      object@model_parameters$TMM$parameters$coordinate_system,
      "cylindrical"
    )
    expect_true(is.list(t_store))
    expect_length(t_store, 2)
    expect_true(all(vapply(t_store, is.list, logical(1))))
    expect_true(all(vapply(t_store, function(x) {
      identical(
        x$family,
        "cylindrical_mono"
      )
    }, logical(1))))
  }
)

test_that("Default cylinder TMM matches FCMS across incident-angle sweeps", {
  old_opt <- options(acousticTS.warn_tmm_cylinder = FALSE)
  on.exit(options(old_opt), add = TRUE)
  density_sw <- 1026.8
  sound_speed_sw <- 1477.3
  frequency <- c(38e3, 70e3, 120e3, 200e3)
  theta_vals <- c(0.55 * pi, 0.7 * pi, 0.85 * pi, pi - 0.05)

  compare_angle_case <- function(boundary, gas = FALSE) {
    for (theta_i in theta_vals) {
      object <- if (gas) {
        gas_generate(
          shape = cylinder(
            length_body = 0.07, radius_body = 0.01,
            n_segments = 80
          ),
          density_fluid = 1.24,
          sound_speed_fluid = 345,
          theta_body = theta_i
        )
      } else {
        fls_generate(
          shape = cylinder(
            length_body = 0.07, radius_body = 0.01,
            n_segments = 80
          ),
          density_body = 1028.9,
          sound_speed_body = 1480.3,
          theta_body = theta_i
        )
      }

      tmm_obj <- target_strength(
        object = object,
        frequency = frequency,
        model = "tmm",
        boundary = boundary,
        density_sw = density_sw,
        sound_speed_sw = sound_speed_sw
      )
      fcms_obj <- target_strength(
        object = object,
        frequency = frequency,
        model = "fcms",
        boundary = boundary,
        density_sw = density_sw,
        sound_speed_sw = sound_speed_sw
      )

      expect_equal(tmm_obj@model$TMM$f_bs, fcms_obj@model$FCMS$f_bs,
        tolerance = 1e-12
      )
      expect_equal(tmm_obj@model$TMM$TS, fcms_obj@model$FCMS$TS,
        tolerance =
          1e-12
      )
    }
  }

  compare_angle_case("fixed_rigid")
  compare_angle_case("pressure_release")
  compare_angle_case("liquid_filled")
  compare_angle_case("gas_filled", gas = TRUE)
})

test_that("Cylinder TMM retained state supports exact monostatic reuse only", {
  old_opt <- options(acousticTS.warn_tmm_cylinder = FALSE)
  on.exit(options(old_opt), add = TRUE)
  density_sw <- 1026.8
  sound_speed_sw <- 1477.3
  object <- target_strength(
    object = fls_generate(
      shape = cylinder(length_body = 0.07, radius_body = 0.01, n_segments = 80),
      density_body = 1028.9,
      sound_speed_body = 1480.3,
      theta_body = pi / 2
    ),
    frequency = c(38e3, 70e3),
    model = "tmm",
    boundary = "fixed_rigid",
    density_sw = density_sw,
    sound_speed_sw = sound_speed_sw,
    store_t_matrix = TRUE
  )

  mono <- tmm_scattering(object)
  avg <- tmm_average_orientation(
    object,
    distribution = tmm_orientation_distribution(
      distribution = "uniform",
      lower = 0.45 * pi,
      upper = pi,
      n_theta = 9
    )
  )
  bundle <- tmm_products(
    object,
    frequency = 70e3,
    orientation = tmm_orientation_distribution(
      distribution = "uniform",
      lower = 0.45 * pi,
      upper = pi,
      n_theta = 9
    ),
    bistatic_summary = FALSE
  )

  expect_true(all(is.finite(mono$sigma_scat)))
  expect_true(all(is.finite(avg$sigma_bs)))
  expect_equal(mono$f_scat, object@model$TMM$f_bs, tolerance = 1e-12)
  expect_equal(bundle$monostatic$f_scat, mono$f_scat, tolerance = 1e-12)
  expect_equal(
    bundle$orientation_average$sigma_bs,
    avg$sigma_bs,
    tolerance = 1e-12
  )
  expect_error(
    tmm_scattering_grid(
      object,
      frequency = 70e3,
      theta_body = pi / 2,
      phi_body = pi,
      theta_scatter = pi / 2,
      phi_scatter = 2 * pi
    ),
    "Stored cylindrical TMM grid evaluations are not available yet"
  )
  expect_error(
    tmm_bistatic_summary(
      object,
      frequency = 70e3,
      n_theta = 21,
      n_phi = 41,
      n_psi = 41
    ),
    "Stored cylindrical TMM bistatic summaries are not available yet"
  )
  expect_error(
    tmm_products(
      object,
      frequency = 70e3,
      orientation = tmm_orientation_distribution(
        distribution = "uniform",
        lower = 0.45 * pi,
        upper = pi,
        n_theta = 9
      ),
      bistatic_summary = TRUE,
      n_theta = 21,
      n_phi = 41,
      n_psi = 41
    ),
    "Stored cylindrical TMM bistatic summaries are not available yet"
  )
})

test_that(
  "Cylinder TMM rejects unsupported general-angle stored scattering requests",
  {
    old_opt <- options(acousticTS.warn_tmm_cylinder = FALSE)
    on.exit(options(old_opt), add = TRUE)
    density_sw <- 1026.8
    sound_speed_sw <- 1477.3

    object <- target_strength(
      object = fls_generate(
        shape = cylinder(
          length_body = 0.07,
          radius_body = 0.01,
          n_segments = 80
        ),
        density_body = 1028.9,
        sound_speed_body = 1480.3,
        theta_body = pi / 2
      ),
      frequency = c(38e3, 70e3),
      model = "tmm",
      boundary = "fixed_rigid",
      density_sw = density_sw,
      sound_speed_sw = sound_speed_sw,
      store_t_matrix = TRUE
    )

    expect_error(
      tmm_scattering(
        object,
        theta_body = pi / 2,
        phi_body = pi,
        theta_scatter = pi / 2,
        phi_scatter = pi / 2
      ),
      "currently available only for the exact monostatic direction"
    )

    expect_error(
      .with_null_pdf_device(
        suppressWarnings(plot(
          object,
          type = "scattering",
          frequency = 70e3,
          polar = TRUE,
          n_theta = 21,
          n_phi = 41
        ))
      ),
      "Stored cylindrical TMM grid evaluations are not available yet"
    )
  }
)

test_that(
  "tmm_diagnostics provides physics-based checks for stored TMM objects",
  {
    old_opt <- options(acousticTS.warn_tmm_cylinder = FALSE)
    on.exit(options(old_opt), add = TRUE)
    density_sw <- 1026.8
    sound_speed_sw <- 1477.3

    sphere_obj <- target_strength(
      object = fls_generate(
        shape = sphere(radius_body = 0.01, n_segments = 80),
        g_body = 1,
        h_body = 1,
        theta_body = pi / 2
      ),
      frequency = c(38e3, 70e3),
      model = "tmm",
      boundary = "fixed_rigid",
      density_sw = density_sw,
      sound_speed_sw = sound_speed_sw,
      store_t_matrix = TRUE
    )

    prolate_obj <- target_strength(
      object = fls_generate(
        shape = prolate_spheroid(
          length_body = 0.14,
          radius_body = 0.01,
          n_segments = 80
        ),
        density_body = 1028.9,
        sound_speed_body = 1480.3,
        theta_body = pi / 2
      ),
      frequency = c(38e3, 70e3),
      model = "tmm",
      boundary = "liquid_filled",
      density_sw = density_sw,
      sound_speed_sw = sound_speed_sw,
      store_t_matrix = TRUE
    )

    cylinder_obj <- target_strength(
      object = fls_generate(
        shape = cylinder(
          length_body = 0.07,
          radius_body = 0.01,
          n_segments = 80
        ),
        density_body = 1028.9,
        sound_speed_body = 1480.3,
        theta_body = pi / 2
      ),
      frequency = c(38e3, 70e3),
      model = "tmm",
      boundary = "fixed_rigid",
      density_sw = density_sw,
      sound_speed_sw = sound_speed_sw,
      store_t_matrix = TRUE
    )

    sphere_diag <- tmm_diagnostics(sphere_obj, n_theta = 31, n_phi = 61)
    prolate_diag <- tmm_diagnostics(prolate_obj, n_theta = 31, n_phi = 61)
    cylinder_diag <- tmm_diagnostics(cylinder_obj, n_theta = 31, n_phi = 61)

    expect_true(all(c(
      "shape", "coordinate_system", "boundary",
      "frequency"
    ) %in% names(sphere_diag$summary)))
    expect_equal(length(sphere_diag$block_metrics), 2L)
    expect_true(all(is.finite(sphere_diag$summary$monostatic_rel_residual)))
    expect_true(all(sphere_diag$summary$reciprocity_rel_residual < 1e-10))
    expect_true(all(is.finite(
      sphere_diag$summary$optical_theorem_rel_residual
    )))
    expect_true(
      all(is.na(sphere_diag$summary$continuation_max_abs_second_diff_TS))
    )

    expect_true(all(prolate_diag$summary$monostatic_rel_residual < 1e-10))
    expect_true(all(prolate_diag$summary$reciprocity_rel_residual < 1e-5))
    expect_true(all(is.na(prolate_diag$summary$optical_theorem_rel_residual)))
    expect_true(all(is.finite(
      prolate_diag$summary$continuation_max_abs_step_TS
    )))
    expect_true(
      all(is.finite(prolate_diag$summary$continuation_max_abs_second_diff_TS))
    )
    expect_false(any(prolate_diag$summary$continuation_any_nonfinite))
    expect_true(is.data.frame(prolate_diag$continuation))
    expect_equal(min(prolate_diag$continuation$aspect_ratio), 1,
      tolerance =
        1e-12
    )
    expect_equal(
      max(prolate_diag$continuation$aspect_ratio),
      prolate_diag$summary$continuation_target_aspect_ratio[1],
      tolerance = 1e-12
    )

    expect_true(all(is.finite(cylinder_diag$summary$monostatic_rel_residual)))
    expect_true(all(is.na(cylinder_diag$summary$reciprocity_rel_residual)))
    expect_true(all(is.na(cylinder_diag$summary$optical_theorem_rel_residual)))
    expect_true(all(is.na(cylinder_diag$summary$min_block_rcond)))
    expect_true(all(is.na(cylinder_diag$summary$max_block_transpose_residual)))
    block_rows <- vapply(cylinder_diag$block_metrics, nrow, integer(1))
    expect_true(all(block_rows == 0L))
    expect_true(
      all(is.na(cylinder_diag$summary$continuation_max_abs_second_diff_TS))
    )
  }
)

test_that(
  paste0(
    "TMM sphere-to-spheroid continuation starts from the exact equal-volume ",
    "sphere"
  ),
  {
    density_sw <- 1026.8
    sound_speed_sw <- 1480

    prolate_obj <- target_strength(
      object = fls_generate(
        shape = prolate_spheroid(
          length_body = 0.07, radius_body = 0.01,
          n_segments = 80
        ),
        g_body = 1,
        h_body = 1,
        theta_body = pi / 2
      ),
      frequency = c(38e3, 70e3),
      model = "tmm",
      boundary = "pressure_release",
      density_sw = density_sw,
      sound_speed_sw = sound_speed_sw,
      store_t_matrix = TRUE
    )
    diag <- tmm_diagnostics(prolate_obj,
      continuation_steps = 5L, n_theta = 21,
      n_phi = 41
    )
    first_step <- diag$continuation[diag$continuation$step == 1, ]

    r_eq <- ((0.07 / 2) * 0.01^2)^(1 / 3)
    sphere_obj <- target_strength(
      object = fls_generate(
        shape = sphere(radius_body = r_eq, n_segments = 80),
        g_body = 1,
        h_body = 1,
        theta_body = pi / 2
      ),
      frequency = c(38e3, 70e3),
      model = "tmm",
      boundary = "pressure_release",
      density_sw = density_sw,
      sound_speed_sw = sound_speed_sw
    )

    expect_true(all(first_step$shape == "Sphere"))
    expect_equal(first_step$TS, sphere_obj@model$TMM$TS, tolerance = 1e-10)
  }
)

test_that("TMM rejects unsupported shapes", {
  object <- fls_generate(
    shape = arbitrary(
      x_body = seq(0, 0.04, length.out = 50),
      zU_body = 0.005 * sin(seq(0, pi, length.out = 50)),
      zL_body = -0.005 * sin(seq(0, pi, length.out = 50)),
      n_segments = 50
    ),
    density_body = 1028.9,
    sound_speed_body = 1480.3
  )

  expect_error(
    target_strength(object,
      frequency = 38e3, model = "tmm", boundary =
        "liquid_filled"
    ),
    "supports the following shape-types"
  )
})