156 lines
7.0 KiB
C++
156 lines
7.0 KiB
C++
|
#include <cmath>
|
|||
|
|
#include <limits>
|
|||
|
|
#include <numbers>
|
|||
|
|
|
|||
|
|
#include <catch2/catch_approx.hpp>
|
|||
|
|
#include <catch2/catch_test_macros.hpp>
|
|||
|
|
|
|||
|
|
#include "serif/utils/laneemden/solution.hpp"
|
|||
|
|
|
|||
|
|
namespace {
|
|||
|
|
namespace laneemden = serif::utils::laneemden;
|
|||
|
|
|
|||
|
|
/** @brief The three analytic regular-center Lane–Emden profiles and their derivatives. */
|
|||
|
|
laneemden::Sample analytic(double n, double xi) {
|
|||
|
|
if (n == 0.0) {
|
|||
|
|
return {.xi = xi, .theta = 1.0 - xi * xi / 6.0, .derivative = -xi / 3.0};
|
|||
|
|
}
|
|||
|
|
if (n == 1.0) {
|
|||
|
|
if (xi == 0.0) {
|
|||
|
|
return {.xi = 0.0, .theta = 1.0, .derivative = 0.0};
|
|||
|
|
}
|
|||
|
|
return {.xi = xi, .theta = std::sin(xi) / xi, .derivative = (xi * std::cos(xi) - std::sin(xi)) / (xi * xi)};
|
|||
|
|
}
|
|||
|
|
const double denominator = 1.0 + xi * xi / 3.0;
|
|||
|
|
return {.xi = xi, .theta = 1.0 / std::sqrt(denominator), .derivative = -xi / (3.0 * std::pow(denominator, 1.5))};
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
TEST_CASE("Lane–Emden matches all regular-center analytic solutions", "[laneemden][analytic]") {
|
|||
|
|
for (double n : {0.0, 1.0, 5.0}) {
|
|||
|
|
CAPTURE(n);
|
|||
|
|
const auto solution = laneemden::solve(n, {.maximum_radius = 20.0});
|
|||
|
|
REQUIRE(solution);
|
|||
|
|
REQUIRE(solution->samples.front().xi == 0.0);
|
|||
|
|
REQUIRE(solution->samples.front().theta == 1.0);
|
|||
|
|
REQUIRE(solution->samples.front().derivative == 0.0);
|
|||
|
|
const double endpoint = solution->samples.back().xi;
|
|||
|
|
for (int index = 0; index <= 400; ++index) {
|
|||
|
|
const double xi = endpoint * static_cast<double>(index) / 400.0;
|
|||
|
|
const auto actual = laneemden::evaluate(*solution, xi);
|
|||
|
|
REQUIRE(actual);
|
|||
|
|
const laneemden::Sample exact = analytic(n, xi);
|
|||
|
|
CHECK(actual->theta == Catch::Approx(exact.theta).margin(3e-10));
|
|||
|
|
CHECK(actual->derivative == Catch::Approx(exact.derivative).margin(3e-9));
|
|||
|
|
}
|
|||
|
|
if (n < 5.0) {
|
|||
|
|
REQUIRE(solution->surface);
|
|||
|
|
const double exact_radius = n == 0.0 ? std::sqrt(6.0) : std::numbers::pi;
|
|||
|
|
CHECK(solution->surface->xi == Catch::Approx(exact_radius).margin(2e-9));
|
|||
|
|
CHECK(solution->surface->derivative == Catch::Approx(analytic(n, exact_radius).derivative).margin(2e-9));
|
|||
|
|
CHECK(solution->samples.back().theta == 0.0);
|
|||
|
|
} else {
|
|||
|
|
CHECK_FALSE(solution->surface);
|
|||
|
|
CHECK(endpoint == 20.0);
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
TEST_CASE("Lane–Emden advances the coupled system at fourth order", "[laneemden][convergence]") {
|
|||
|
|
double previous_error = 0.0;
|
|||
|
|
for (double step : {0.08, 0.04, 0.02, 0.01}) {
|
|||
|
|
const auto solution = laneemden::solve(1.0, {.step = step, .maximum_radius = 2.0});
|
|||
|
|
REQUIRE(solution);
|
|||
|
|
const double error = std::abs(solution->samples.back().theta - analytic(1.0, 2.0).theta);
|
|||
|
|
if (previous_error > 0.0) {
|
|||
|
|
CHECK(previous_error / error > 14.0);
|
|||
|
|
CHECK(previous_error / error < 18.0);
|
|||
|
|
}
|
|||
|
|
previous_error = error;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
TEST_CASE("Lane–Emden finds stellar surfaces for fractional and centrally concentrated profiles", "[laneemden][surface]") {
|
|||
|
|
struct Reference final {
|
|||
|
|
double n;
|
|||
|
|
double radius;
|
|||
|
|
double mass;
|
|||
|
|
double mass_tolerance;
|
|||
|
|
};
|
|||
|
|
// Standard first-zero values provide an independent endpoint check beyond the analytic indices.
|
|||
|
|
for (const Reference reference : {
|
|||
|
|
Reference{.n = 0.5, .radius = 2.752698054065, .mass = 3.78865118488, .mass_tolerance = 2e-8},
|
|||
|
|
Reference{.n = 1.5, .radius = 3.653753736219, .mass = 2.714055120109, .mass_tolerance = 2e-9},
|
|||
|
|
Reference{.n = 3.0, .radius = 6.896848619377, .mass = 2.018235950966, .mass_tolerance = 2e-9}}) {
|
|||
|
|
CAPTURE(reference.n);
|
|||
|
|
const auto solution = laneemden::solve(reference.n, {.step = 0.0005});
|
|||
|
|
REQUIRE(solution);
|
|||
|
|
REQUIRE(solution->surface);
|
|||
|
|
const laneemden::Sample surface = *solution->surface;
|
|||
|
|
CHECK(surface.xi == Catch::Approx(reference.radius).margin(2e-9));
|
|||
|
|
CHECK(-surface.xi * surface.xi * surface.derivative == Catch::Approx(reference.mass).margin(reference.mass_tolerance));
|
|||
|
|
for (std::size_t index = 1; index < solution->samples.size(); ++index) {
|
|||
|
|
const laneemden::Sample &previous = solution->samples[index - 1];
|
|||
|
|
const laneemden::Sample ¤t = solution->samples[index];
|
|||
|
|
REQUIRE(current.xi > previous.xi);
|
|||
|
|
REQUIRE(current.theta >= 0.0);
|
|||
|
|
REQUIRE(current.theta < previous.theta);
|
|||
|
|
REQUIRE(std::isfinite(current.derivative));
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
TEST_CASE("Fractional-index surface values converge without taking nonreal powers", "[laneemden][convergence]") {
|
|||
|
|
double previous_error = 0.0;
|
|||
|
|
for (double step : {0.04, 0.02, 0.01, 0.005}) {
|
|||
|
|
const auto solution = laneemden::solve(0.5, {.step = step});
|
|||
|
|
REQUIRE(solution);
|
|||
|
|
REQUIRE(solution->surface);
|
|||
|
|
const laneemden::Sample surface = *solution->surface;
|
|||
|
|
const double error = std::abs(-surface.xi * surface.xi * surface.derivative - 3.78865118488);
|
|||
|
|
if (previous_error > 0.0) {
|
|||
|
|
CHECK(error < 0.45 * previous_error);
|
|||
|
|
}
|
|||
|
|
previous_error = error;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
TEST_CASE("Lane–Emden distinguishes finite interval profiles from failures", "[laneemden]") {
|
|||
|
|
for (double n : {1.0, 5.0, 6.0}) {
|
|||
|
|
const auto solution = laneemden::solve(n, {.maximum_radius = 1.0});
|
|||
|
|
REQUIRE(solution);
|
|||
|
|
CHECK_FALSE(solution->surface);
|
|||
|
|
CHECK(solution->samples.back().xi == 1.0);
|
|||
|
|
CHECK(solution->samples.back().theta > 0.0);
|
|||
|
|
}
|
|||
|
|
const auto exhausted = laneemden::solve(3.0, {.maximum_steps = 2});
|
|||
|
|
CHECK_FALSE(exhausted);
|
|||
|
|
const auto tiny_interval = laneemden::solve(1.0, {.maximum_radius = 1e-6});
|
|||
|
|
REQUIRE(tiny_interval);
|
|||
|
|
CHECK(tiny_interval->samples.size() == 2);
|
|||
|
|
CHECK(tiny_interval->samples.back().derivative == Catch::Approx(-1e-6 / 3.0).margin(1e-18));
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
TEST_CASE("Lane–Emden rejects invalid controls and evaluation radii", "[laneemden][rejection]") {
|
|||
|
|
const double infinity = std::numeric_limits<double>::infinity();
|
|||
|
|
const double nan = std::numeric_limits<double>::quiet_NaN();
|
|||
|
|
for (double invalid : {-1.0, infinity, nan}) {
|
|||
|
|
CHECK_FALSE(laneemden::solve(invalid));
|
|||
|
|
}
|
|||
|
|
for (double invalid : {0.0, -1.0, infinity, nan}) {
|
|||
|
|
CHECK_FALSE(laneemden::solve(1.0, {.step = invalid}));
|
|||
|
|
CHECK_FALSE(laneemden::solve(1.0, {.maximum_radius = invalid}));
|
|||
|
|
CHECK_FALSE(laneemden::solve(1.0, {.surface_tolerance = invalid}));
|
|||
|
|
}
|
|||
|
|
CHECK_FALSE(laneemden::solve(1.0, {.maximum_steps = 0}));
|
|||
|
|
CHECK_FALSE(laneemden::solve(1.0, {.surface_tolerance = 0.01}));
|
|||
|
|
CHECK_FALSE(laneemden::solve(1.0, {.central_series_terms = 5}));
|
|||
|
|
CHECK_FALSE(laneemden::solve(1.0, {.central_series_terms = 65}));
|
|||
|
|
const auto solution = laneemden::solve(1.0);
|
|||
|
|
REQUIRE(solution);
|
|||
|
|
for (double invalid : {-1.0, infinity, nan, 4.0}) {
|
|||
|
|
CHECK_FALSE(laneemden::evaluate(*solution, invalid));
|
|||
|
|
}
|
|||
|
|
CHECK_FALSE(laneemden::evaluate(laneemden::Solution{.n = 1.0, .samples = {}, .surface = {}}, 0.0));
|
|||
|
|
}
|