tboudreaux/MeanField_SERiF
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

feat(utils/laneemden): 1D lane-emden solver

Browse Source 
simple rk4 lane emden solver. reimplemented from work in Dartmouth A115
This commit is contained in:
Emily Boudreaux
2026-11-09 13:39:32 -05:00
parent ba457e7d10
commit 7134f99d7c
5 changed files with 515 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
src/include/serif/utils/laneemden/laneemden.hpp Normal file
View File

@@ -0,0 +1,4 @@
#pragma once
#include "serif/utils/laneemden/solution.hpp"
#include "serif/utils/laneemden/scaling.hpp"

63
src/include/serif/utils/laneemden/scaling.hpp Normal file
View File

@@ -0,0 +1,63 @@
#pragma once
#include <expected>
#include <optional>
#include <string>
#include "serif/dimensions/dimensions.hpp"
#include "serif/utils/laneemden/solution.hpp"
namespace serif::utils::laneemden {
/** @brief A physical conversion of one finite-radius dimensionless polytrope, in the caller's consistent units. */
struct PhysicalScales final {
double n;
double radial_scale;
dimensions::QuantityValue<dimensions::Mass> mass;
dimensions::QuantityValue<dimensions::Length> radius;
dimensions::QuantityValue<dimensions::MassDensity> central_density;
dimensions::QuantityValue<dimensions::SpecificEnthalpy> central_enthalpy;
dimensions::QuantityValue<dimensions::GravitationalConstant> gravitational_constant;
// K has n-dependent dimensions; absent for n = 0 or when this optional derived value is not representable.
std::optional<double> polytropic_constant;
};
/** @brief Physical scalar fields; gravity is the outward potential derivative, not the inward force component. */
struct PhysicalSample final {
dimensions::QuantityValue<dimensions::MassDensity> density;
dimensions::QuantityValue<dimensions::Pressure> pressure;
dimensions::QuantityValue<dimensions::SpecificEnthalpy> enthalpy;
dimensions::QuantityValue<dimensions::GravitationalPotential> potential;
dimensions::QuantityValue<dimensions::Mass> enclosed_mass;
dimensions::QuantityValue<dimensions::Acceleration> gravity;
};
/** @brief Fix a reference star by mass and central density; its implied K is an output, never a change to an EOS. */
[[nodiscard]] std::expected<PhysicalScales, std::string> scale_from_mass(
const Solution &solution,
dimensions::QuantityValue<dimensions::Mass> mass,
dimensions::QuantityValue<dimensions::MassDensity> central_density,
dimensions::QuantityValue<dimensions::GravitationalConstant> gravitational_constant);
/** @brief Fix a physical star by central enthalpy and density, including the incompressible n = 0 limit. */
[[nodiscard]] std::expected<PhysicalScales, std::string> scale_from_enthalpy(
const Solution &solution,
dimensions::QuantityValue<dimensions::SpecificEnthalpy> central_enthalpy,
dimensions::QuantityValue<dimensions::MassDensity> central_density,
dimensions::QuantityValue<dimensions::GravitationalConstant> gravitational_constant);
/** @brief Fix a polytropic EOS through its K and central density; n must be strictly positive. */
[[nodiscard]] std::expected<PhysicalScales, std::string> scale_from_polytropic_constant(
const Solution &solution,
double polytropic_constant,
dimensions::QuantityValue<dimensions::MassDensity> central_density,
dimensions::QuantityValue<dimensions::GravitationalConstant> gravitational_constant);
/**
* @brief Evaluate physical fields, continuing into vacuum beyond the stellar surface.
*
* The potential vanishes at infinity. Density, pressure, and enthalpy are zero at and
* beyond the surface, including the exterior side of the incompressible density jump.
*/
[[nodiscard]] std::expected<PhysicalSample, std::string> evaluate_physical(
const Solution &solution, const PhysicalScales &scales, dimensions::QuantityValue<dimensions::Length> radius);
}

141
src/include/serif/utils/laneemden/solution.hpp Normal file
View File

@@ -0,0 +1,141 @@
#pragma once
#include <cstddef>
#include <expected>
#include <optional>
#include <string>
#include <string_view>
#include <vector>
#include <format>
#include "serif/utils/error/errors.hpp"
namespace serif::utils::laneemden {
/** @brief Dimensionless integration controls; the surface tolerance does not bound the global integration error. */
struct IntegrationOptions final {
double step = 0.01;
double maximum_radius = 100.0;
std::size_t maximum_steps = 100000;
double surface_tolerance = 1e-10;
std::size_t central_series_terms = 12;
/** @brief Compare the numerical prescription when tracking reference dependencies. */
friend bool operator==(const IntegrationOptions &, const IntegrationOptions &) = default;
};
/** @brief The regular Lane–Emden solution and its first derivative at one dimensionless radius. */
struct Sample final {
double xi;
double theta;
double derivative;
};
/** @brief A regular profile, ending at its first zero or the requested integration radius. */
struct Solution final {
double n;
std::vector<Sample> samples;
std::optional<Sample> surface;
};
/**
* @brief Integrate theta'' + 2 theta'/xi + theta^n = 0 from the regular center.
*
* The central series follows Boudreaux's original power-series recurrence. Both dependent
* variables now advance through the same four Runge–Kutta stages. A missing surface is a
* valid finite-interval result, including the n >= 5 solutions. Exhausting the step budget
* or producing an unrepresentable numerical state returns a failure instead.
*/
[[nodiscard]] std::expected<Solution, std::string> solve(double n, IntegrationOptions options = {});
/** @brief Evaluate a solved dimensionless profile with cubic Hermite interpolation, without extrapolation. */
[[nodiscard]] std::expected<Sample, std::string> evaluate(const Solution &solution, double xi);
}
template <>
struct std::formatter<serif::utils::laneemden::IntegrationOptions> {
template <typename FormatContext>
auto format(const serif::utils::laneemden::IntegrationOptions &options, FormatContext &ctx) const {
return std::format_to(ctx.out(), "step={}, maximum_radius={}, maximum_steps={}, surface_tolerance={}, central_series_terms={}",
options.step, options.maximum_radius, options.maximum_steps, options.surface_tolerance, options.central_series_terms);
}
};
/** @brief Select all sample values, or one named component with :xi, :theta, or :derivative. */
template <>
struct std::formatter<serif::utils::laneemden::Sample> {
enum class Component { all, xi, theta, derivative };
Component component{Component::all};
/** @brief Validate the selector during literal format checking or runtime format parsing. */
constexpr auto parse(std::format_parse_context &ctx) {
using namespace serif::utils::errors;
component = Component::all;
auto end = ctx.begin();
while (end != ctx.end() && *end != '}') {
++end;
}
if (end == ctx.begin()) {
return end;
}
const std::string_view selection(ctx.begin(), end);
if (selection == "all") {
component = Component::all;
} else if (selection == "xi") {
component = Component::xi;
} else if (selection == "theta") {
component = Component::theta;
} else if (selection == "derivative") {
component = Component::derivative;
} else {
serif_error<std::format_error, SERiFErrorCode::INVALID_ARGUMENT>(
"Unknown Lane–Emden sample selector '{}'. Use an empty specifier, 'all', 'xi', 'theta', or 'derivative'.", selection);
}
return end;
}
/** @brief Keep the labeled tuple for complete samples and emit a bare number for a selected component. */
template <typename FormatContext>
auto format(const serif::utils::laneemden::Sample &sample, FormatContext &ctx) const {
switch (component) {
case Component::xi:
return std::format_to(ctx.out(), "{}", sample.xi);
case Component::theta:
return std::format_to(ctx.out(), "{}", sample.theta);
case Component::derivative:
return std::format_to(ctx.out(), "{}", sample.derivative);
default:
return std::format_to(ctx.out(), "(ξ={}, θ={}, dθ/dξ={})", sample.xi, sample.theta, sample.derivative);
}
}
};
/** @brief Describe the solution while applying one shared selection to its profile and optional surface sample. */
template <>
struct std::formatter<serif::utils::laneemden::Solution> {
std::formatter<serif::utils::laneemden::Sample> sample_formatter;
/** @brief Delegate selection to Sample so the two formatters cannot develop different grammars. */
constexpr auto parse(std::format_parse_context &ctx) {
return sample_formatter.parse(ctx);
}
/** @brief Stream the profile through its parsed sample formatter without building intermediate strings. */
template <typename FormatContext>
auto format(const serif::utils::laneemden::Solution &solution, FormatContext &ctx) const {
ctx.advance_to(std::format_to(ctx.out(), "Solution(n={}, samples=[", solution.n));
for (std::size_t index = 0; index < solution.samples.size(); ++index) {
if (index != 0) {
ctx.advance_to(std::format_to(ctx.out(), ", "));
}
ctx.advance_to(sample_formatter.format(solution.samples[index], ctx));
}
ctx.advance_to(std::format_to(ctx.out(), "], surface="));
if (solution.surface) {
ctx.advance_to(sample_formatter.format(*solution.surface, ctx));
} else {
ctx.advance_to(std::format_to(ctx.out(), "none"));
}
return std::format_to(ctx.out(), ")");
}
};

139
src/lib/serif/utils/laneemden/scaling.cpp Normal file
View File

@@ -0,0 +1,139 @@
#include "serif/utils/laneemden/scaling.hpp"
#include <algorithm>
#include <cmath>
#include <numbers>
namespace serif::utils::laneemden {
namespace {
constexpr long double four_pi = 4.0L * std::numbers::pi_v<long double>;
/** @brief Validate a positive physical scale without assuming a particular unit system. */
bool positive_finite(double value) {
return std::isfinite(value) && value > 0.0;
}
/** @brief Require the finite surface needed to assign a stellar radius and total mass. */
bool finite_surface(const Solution &solution) {
return std::isfinite(solution.n) && solution.n >= 0.0 && solution.surface && positive_finite(solution.surface->xi) &&
std::isfinite(solution.surface->derivative) && solution.surface->derivative < 0.0;
}
/** @brief Collect the shared scaling identities, keeping intermediate products representable in physical units. */
std::expected<PhysicalScales, std::string> make_scales(
const Solution &solution, long double radial_scale, double central_density, double gravitational_constant) {
const long double surface_xi = solution.surface->xi;
const long double dimensionless_mass = -surface_xi * surface_xi * solution.surface->derivative;
const double radius = static_cast<double>(radial_scale * surface_xi);
const double mass = static_cast<double>(four_pi * central_density * radial_scale * radial_scale * radial_scale * dimensionless_mass);
const double enthalpy = static_cast<double>(four_pi * gravitational_constant * central_density * radial_scale * radial_scale);
const double alpha = static_cast<double>(radial_scale);
if (!positive_finite(alpha) || !positive_finite(radius) || !positive_finite(mass) || !positive_finite(enthalpy)) {
return std::unexpected("Lane–Emden physical scales overflow or underflow; choose representable consistent units.");
}
std::optional<double> polytropic_constant;
if (solution.n > 0.0) {
const double constant = static_cast<double>(enthalpy / ((solution.n + 1.0L) * std::pow(static_cast<long double>(central_density), 1.0L / solution.n)));
if (positive_finite(constant)) {
polytropic_constant = constant;
}
}
return PhysicalScales{
.n = solution.n,
.radial_scale = alpha,
.mass = dimensions::QuantityValue<dimensions::Mass>{mass},
.radius = dimensions::QuantityValue<dimensions::Length>{radius},
.central_density = dimensions::QuantityValue<dimensions::MassDensity>{central_density},
.central_enthalpy = dimensions::QuantityValue<dimensions::SpecificEnthalpy>{enthalpy},
.gravitational_constant = dimensions::QuantityValue<dimensions::GravitationalConstant>{gravitational_constant},
.polytropic_constant = polytropic_constant
};
}
}
std::expected<PhysicalScales, std::string> scale_from_mass(
const Solution &solution,
dimensions::QuantityValue<dimensions::Mass> mass,
dimensions::QuantityValue<dimensions::MassDensity> central_density,
dimensions::QuantityValue<dimensions::GravitationalConstant> gravitational_constant) {
if (!finite_surface(solution) || !positive_finite(mass.value()) || !positive_finite(central_density.value()) || !positive_finite(gravitational_constant.value())) {
return std::unexpected("Lane–Emden mass scaling needs a finite surface and positive finite mass, central density, and G.");
}
const long double xi = solution.surface->xi;
const long double dimensionless_mass = -xi * xi * solution.surface->derivative;
const long double alpha = std::cbrt(mass.value() / (four_pi * central_density.value() * dimensionless_mass));
return make_scales(solution, alpha, central_density.value(), gravitational_constant.value());
}
std::expected<PhysicalScales, std::string> scale_from_enthalpy(
const Solution &solution,
dimensions::QuantityValue<dimensions::SpecificEnthalpy> central_enthalpy,
dimensions::QuantityValue<dimensions::MassDensity> central_density,
dimensions::QuantityValue<dimensions::GravitationalConstant> gravitational_constant) {
if (!finite_surface(solution) || !positive_finite(central_enthalpy.value()) || !positive_finite(central_density.value()) || !positive_finite(gravitational_constant.value())) {
return std::unexpected("Lane–Emden enthalpy scaling needs a finite surface and positive finite enthalpy, central density, and G.");
}
const long double alpha = std::sqrt(central_enthalpy.value() / (four_pi * gravitational_constant.value() * central_density.value()));
return make_scales(solution, alpha, central_density.value(), gravitational_constant.value());
}
std::expected<PhysicalScales, std::string> scale_from_polytropic_constant(
const Solution &solution, double polytropic_constant,
dimensions::QuantityValue<dimensions::MassDensity> central_density,
dimensions::QuantityValue<dimensions::GravitationalConstant> gravitational_constant) {
if (!finite_surface(solution) || solution.n <= 0.0 || !positive_finite(polytropic_constant) || !positive_finite(central_density.value()) ||
!positive_finite(gravitational_constant.value())) {
return std::unexpected("Lane–Emden EOS scaling needs n > 0, a finite surface, and positive finite K, central density, and G.");
}
const double enthalpy = static_cast<double>((solution.n + 1.0L) * polytropic_constant * std::pow(static_cast<long double>(central_density.value()), 1.0L / solution.n));
auto scales = scale_from_enthalpy(solution, dimensions::QuantityValue<dimensions::SpecificEnthalpy>{enthalpy}, central_density, gravitational_constant);
if (scales) {
scales->polytropic_constant = polytropic_constant;
}
return scales;
}
std::expected<PhysicalSample, std::string> evaluate_physical(
const Solution &solution, const PhysicalScales &scales, dimensions::QuantityValue<dimensions::Length> radius) {
if (!finite_surface(solution) || solution.n != scales.n || !std::isfinite(radius.value()) || radius.value() < 0.0 ||
!positive_finite(scales.radial_scale) || !positive_finite(scales.radius.value()) || !positive_finite(scales.mass.value()) ||
!positive_finite(scales.central_density.value()) || !positive_finite(scales.central_enthalpy.value()) || !positive_finite(scales.gravitational_constant.value())) {
return std::unexpected("Lane–Emden physical evaluation requires matching finite-surface scales and a finite nonnegative radius.");
}
const long double gravitational_mass = static_cast<long double>(scales.gravitational_constant.value()) * scales.mass.value();
double density = 0.0;
double pressure = 0.0;
double enthalpy = 0.0;
double mass = scales.mass.value();
double potential;
double gravity;
if (radius.value() >= scales.radius.value()) {
potential = static_cast<double>(-gravitational_mass / radius.value());
gravity = static_cast<double>(gravitational_mass / radius.value() / radius.value());
} else {
const double xi = std::min(radius.value() / scales.radial_scale, solution.surface->xi);
const auto sample = evaluate(solution, xi);
if (!sample) {
return std::unexpected(sample.error());
}
const double theta = std::max(0.0, sample->theta);
density = scales.central_density.value() * std::pow(theta, solution.n);
enthalpy = scales.central_enthalpy.value() * theta;
pressure = static_cast<double>(static_cast<long double>(density) * enthalpy / (solution.n + 1.0L));
potential = static_cast<double>(-gravitational_mass / scales.radius.value() - enthalpy);
gravity = -scales.central_enthalpy.value() / scales.radial_scale * sample->derivative;
mass = static_cast<double>(-four_pi * scales.central_density.value() * std::pow(static_cast<long double>(scales.radial_scale), 3) * xi * xi * sample->derivative);
}
if (!std::isfinite(density) || !std::isfinite(pressure) || !std::isfinite(enthalpy) || !std::isfinite(potential) || !std::isfinite(gravity) || !std::isfinite(mass)) {
return std::unexpected("Lane–Emden physical fields overflow in the supplied units.");
}
return PhysicalSample{
.density = dimensions::QuantityValue<dimensions::MassDensity>{density},
.pressure = dimensions::QuantityValue<dimensions::Pressure>{pressure},
.enthalpy = dimensions::QuantityValue<dimensions::SpecificEnthalpy>{enthalpy},
.potential = dimensions::QuantityValue<dimensions::GravitationalPotential>{potential},
.enclosed_mass = dimensions::QuantityValue<dimensions::Mass>{mass},
.gravity = dimensions::QuantityValue<dimensions::Acceleration>{gravity}
};
}
}

168
src/lib/serif/utils/laneemden/solution.cpp Normal file
View File

@@ -0,0 +1,168 @@
#include "serif/utils/laneemden/solution.hpp"
#include <algorithm>
#include <array>
#include <cmath>
#include <limits>
namespace serif::utils::laneemden {
namespace {
/** @brief Evaluate Boudreaux's coupled a_k and c_k series iteratively, avoiding exponential recursive work. */
Sample central_sample(double n, double xi, std::size_t terms) {
std::array<double, 64> coefficients{};
std::array<double, 64> powers{};
coefficients[0] = 1.0;
powers[0] = 1.0;
for (std::size_t order = 1; order < terms; ++order) {
if (order >= 2) {
coefficients[order] = -powers[order - 2] / static_cast<double>(order * (order + 1));
}
double sum = 0.0;
for (std::size_t k = 1; k <= order; ++k) {
sum += (static_cast<double>(k) * (n + 1.0) - static_cast<double>(order)) * coefficients[k] * powers[order - k];
}
powers[order] = sum / static_cast<double>(order);
}
double theta = coefficients[terms - 1];
double derivative = 0.0;
for (std::size_t order = terms - 1; order > 0; --order) {
derivative = derivative * xi + theta;
theta = theta * xi + coefficients[order - 1];
}
return {.xi = xi, .theta = theta, .derivative = derivative};
}
/** @brief Evaluate the regular positive-density equation, using vacuum continuation only while locating the first zero. */
Sample slope(double n, const Sample &sample) {
const double source = std::pow(std::max(0.0, sample.theta), n);
return {.xi = 1.0, .theta = sample.derivative, .derivative = -2.0 * sample.derivative / sample.xi - source};
}
/** @brief Construct a Runge–Kutta stage with both coupled variables evaluated at the same radius. */
Sample stage(const Sample &initial, const Sample &rate, double step) {
return {
.xi = initial.xi + step,
.theta = initial.theta + step * rate.theta,
.derivative = initial.derivative + step * rate.derivative
};
}
/** @brief Advance the coupled Lane–Emden system with classical fourth-order Runge–Kutta. */
Sample advance(double n, const Sample &initial, double step) {
const Sample first = slope(n, initial);
const Sample second = slope(n, stage(initial, first, 0.5 * step));
const Sample third = slope(n, stage(initial, second, 0.5 * step));
const Sample fourth = slope(n, stage(initial, third, step));
return {
.xi = initial.xi + step,
.theta = initial.theta + step * (first.theta + 2.0 * second.theta + 2.0 * third.theta + fourth.theta) / 6.0,
.derivative = initial.derivative + step * (first.derivative + 2.0 * second.derivative + 2.0 * third.derivative + fourth.derivative) / 6.0
};
}
/** @brief Reject numerical states that have overflowed or left the regular positive branch. */
bool valid_interior(const Sample &sample) {
return std::isfinite(sample.xi) && std::isfinite(sample.theta) && std::isfinite(sample.derivative) &&
sample.theta > 0.0 && sample.derivative < 0.0;
}
/** @brief Resolve the first zero from a short bracket without raising negative theta to a fractional power. */
std::expected<Sample, std::string> locate_surface(double n, const Sample &initial, double width, double tolerance) {
double lower = 0.0;
double upper = width;
Sample result = advance(n, initial, upper);
if (!std::isfinite(result.theta) || result.theta > 0.0) {
return std::unexpected("Lane–Emden surface could not be bracketed; reduce the integration step.");
}
for (std::size_t iteration = 0; iteration < 80 && upper - lower > tolerance; ++iteration) {
const double middle = lower + 0.5 * (upper - lower);
result = advance(n, initial, middle);
if (result.theta > 0.0) {
lower = middle;
} else {
upper = middle;
}
}
result = advance(n, initial, lower + 0.5 * (upper - lower));
result.theta = 0.0;
if (!std::isfinite(result.derivative) || result.derivative >= 0.0) {
return std::unexpected("Lane–Emden surface derivative is not finite and negative; reduce the integration step.");
}
return result;
}
}
std::expected<Solution, std::string> solve(const double n, IntegrationOptions options) {
if (!std::isfinite(n) || n < 0.0) {
return std::unexpected("Lane–Emden index must be finite and nonnegative.");
}
if (!std::isfinite(options.step) || options.step <= 0.0 || !std::isfinite(options.maximum_radius) || options.maximum_radius <= 0.0 ||
!std::isfinite(options.surface_tolerance) || options.surface_tolerance <= 0.0 || options.surface_tolerance >= 0.01 ||
options.maximum_steps == 0 || options.central_series_terms < 6 || options.central_series_terms > 64) {
return std::unexpected("Lane–Emden integration needs positive finite step and maximum radius, a surface tolerance in (0, 0.01), a nonzero step budget, and 6 to 64 series terms.");
}
Solution solution{.n = n, .samples = {{.xi = 0.0, .theta = 1.0, .derivative = 0.0}}, .surface = std::nullopt};
const double first_radius = std::min({options.step, options.maximum_radius, 0.01 / std::sqrt(std::max(1.0, n))});
Sample current = central_sample(n, first_radius, options.central_series_terms);
if (!valid_interior(current)) {
return std::unexpected("Lane–Emden central series is not representable for the supplied index and step.");
}
solution.samples.push_back(current);
for (std::size_t iteration = 0; current.xi < options.maximum_radius; ++iteration) {
if (iteration >= options.maximum_steps) {
return std::unexpected("Lane–Emden integration exhausted maximum_steps before reaching the surface or maximum_radius.");
}
const double distance = -current.theta / current.derivative;
const double tolerance = options.surface_tolerance * std::max(1.0, current.xi);
if (distance <= tolerance && current.xi + 2.0 * distance <= options.maximum_radius) {
const double root_tolerance = std::max(4.0 * std::numeric_limits<double>::epsilon() * current.xi, tolerance * 0.01);
const auto surface = locate_surface(n, current, 2.0 * distance, root_tolerance);
if (!surface) {
return std::unexpected(surface.error());
}
solution.surface = *surface;
solution.samples.push_back(*surface);
return solution;
}
// Shorten only the final approach; no accepted sample extends past the first zero.
const double step = std::min({options.step, 0.5 * distance, options.maximum_radius - current.xi});
if (current.xi + step == current.xi) {
return std::unexpected("Lane–Emden step cannot advance the radius at machine precision; relax surface_tolerance.");
}
current = advance(n, current, step);
if (!valid_interior(current)) {
return std::unexpected("Lane–Emden integration left the finite positive branch; reduce the integration step.");
}
solution.samples.push_back(current);
}
return solution;
}
std::expected<Sample, std::string> evaluate(const Solution &solution, double xi) {
if (!std::isfinite(xi) || xi < 0.0 || solution.samples.empty() || xi > solution.samples.back().xi) {
return std::unexpected("Lane–Emden evaluation radius must lie within the stored finite dimensionless profile.");
}
const auto upper = std::lower_bound(solution.samples.begin(), solution.samples.end(), xi,
[](const Sample &sample, double radius) { return sample.xi < radius; });
if (upper == solution.samples.begin() || upper->xi == xi) {
return *upper;
}
const Sample &left = *(upper - 1);
const Sample &right = *upper;
const double width = right.xi - left.xi;
if (!std::isfinite(width) || width <= 0.0) {
return std::unexpected("Lane–Emden evaluation requires a profile with increasing finite radii produced by solve.");
}
const double position = (xi - left.xi) / width;
const double quadratic = 3.0 * (right.theta - left.theta) - width * (2.0 * left.derivative + right.derivative);
const double cubic = 2.0 * (left.theta - right.theta) + width * (left.derivative + right.derivative);
return Sample{
.xi = xi,
.theta = left.theta + position * (width * left.derivative + position * (quadratic + position * cubic)),
.derivative = left.derivative + position * (2.0 * quadratic + 3.0 * position * cubic) / width
};
}
}