86 lines
3.7 KiB
Python
86 lines
3.7 KiB
Python
#!/usr/bin/env python3
|
||
"""Independent n=3 Lane–Emden integration; no FEM or production code changes."""
|
||
from fractions import Fraction
|
||
from pathlib import Path
|
||
import json
|
||
import math
|
||
import sys
|
||
|
||
import scipy
|
||
from scipy.integrate import solve_ivp
|
||
|
||
|
||
def central_series(x, terms=10):
|
||
# Recurrence from theta'' + 2 theta'/x + theta^3 = 0 avoids
|
||
# a singular initial-value evaluation and a hand-transcribed series.
|
||
coeff = [Fraction(1)]
|
||
for k in range(terms - 1):
|
||
cube = sum(
|
||
coeff[i] * coeff[j] * coeff[k-i-j]
|
||
for i in range(k+1) for j in range(k-i+1)
|
||
)
|
||
coeff.append(-cube / ((2*k+2)*(2*k+3)))
|
||
theta = math.fsum(float(a)*x**(2*k) for k, a in enumerate(coeff))
|
||
mu = -math.fsum(2*k*float(a)*x**(2*k+1) for k, a in enumerate(coeff) if k)
|
||
return theta, mu
|
||
|
||
|
||
def solve(rtol, atol, maximum_step, start):
|
||
def rhs(x, y):
|
||
theta, mu = y
|
||
return -mu / x**2, x**2 * theta**3
|
||
|
||
def zero(x, y):
|
||
return y[0]
|
||
zero.terminal = True
|
||
zero.direction = -1
|
||
|
||
result = solve_ivp(rhs, (start, 10.0), central_series(start),
|
||
events=zero, method="DOP853", rtol=rtol, atol=atol,
|
||
max_step=maximum_step)
|
||
if not result.success or len(result.t_events[0]) != 1:
|
||
raise RuntimeError(result.message)
|
||
xi = float(result.t_events[0][0])
|
||
mu = float(result.y_events[0][0, 1])
|
||
kcrit = math.pi / (4.0*math.pi*mu)**(2.0/3.0)
|
||
rho_for_unit_radius = xi**3/(4.0*math.pi*mu)
|
||
mass = 4.0*math.pi*(driver_k/math.pi)**1.5*mu
|
||
radius = xi*math.sqrt(driver_k/math.pi)*driver_rho**(-1.0/3.0)
|
||
return dict(rtol=rtol, atol=atol, maximum_step=maximum_step, start=start,
|
||
steps=len(result.t)-1, function_evaluations=result.nfev,
|
||
xi1=xi, mu1=mu, K_for_mass_one=kcrit,
|
||
rho_c_for_mass_and_radius_one=rho_for_unit_radius,
|
||
M_LE_at_driver_K=mass, M_target_minus_M_LE=1.0-mass,
|
||
R_LE_at_driver_K_and_rho_c=radius,
|
||
K_driver_relative_to_critical_minus_one=driver_k/kcrit-1.0,
|
||
virial_bound_T_over_absW_max=0.5*(1.0-mass**(2.0/3.0)))
|
||
|
||
|
||
metadata = Path(sys.argv[1]).resolve()
|
||
values = dict(line.split("=", 1) for line in metadata.read_text().splitlines() if "=" in line)
|
||
if float(values["n"]) != 3.0 or float(values["M"]) != 1.0:
|
||
raise ValueError("This diagnostic uses the campaign's n=3, M=1, G=1 convention")
|
||
driver_k = float(values["K"])
|
||
driver_rho = float(values["rho_c"])
|
||
checks = [solve(3e-12, 1e-14, 0.2, 1e-3),
|
||
solve(3e-13, 1e-15, 0.1, 1e-3),
|
||
solve(3e-14, 1e-16, 0.05, 1e-3),
|
||
solve(3e-14, 1e-16, 0.025, 1e-2),
|
||
solve(3e-14, 1e-16, 0.025, 1e-4)]
|
||
mu = checks[-1]["mu1"]
|
||
def mass_for_k(k):
|
||
return 4.0*math.pi*(k/math.pi)**1.5*mu
|
||
|
||
out = dict(metadata=str(metadata), python=sys.version, scipy=scipy.__version__,
|
||
convention="G=1; theta'= -mu/xi^2; mu'=xi^2 theta^3; first theta zero",
|
||
driver={name: values[name] for name in ("n", "K", "rho_c", "M", "J", "xi1", "mu1")},
|
||
checks=checks,
|
||
mass_check_spread=max(x["M_LE_at_driver_K"] for x in checks)-min(x["M_LE_at_driver_K"] for x in checks),
|
||
full_metadata_K_binary64_ulp=math.ulp(driver_k),
|
||
M_LE_for_K_adjacent_binary64=[mass_for_k(math.nextafter(driver_k, -math.inf)),
|
||
mass_for_k(math.nextafter(driver_k, math.inf))],
|
||
K_quoted_12_decimal_rounding_interval=[0.3639390599515, 0.3639390599525],
|
||
M_LE_for_K_quoted_rounding_interval=[mass_for_k(0.3639390599515), mass_for_k(0.3639390599525)],
|
||
limitation="Tolerance/step/start sensitivity check, not an interval-arithmetic error certificate.")
|
||
print(json.dumps(out, indent=2))
|