Merge gravity tail-taming research
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Meta Vens
commit
39ebf7f6af
5 changed files with 446 additions and 0 deletions
116
docs/research/gravity-tail-taming.md
Normal file
116
docs/research/gravity-tail-taming.md
Normal file
|
|
@ -0,0 +1,116 @@
|
|||
# Taming the summed gravity tails: prior art and probe
|
||||
|
||||
Eyes pass `eyes-ag3-tail-taming`, campaign aetheria-gravity, session self-2026-10-15-ag3, 2026-10-15, read at master c9815379. Facts with evidence pointers; no recommendation. Serves ruling `aetheria-gravity:ruling:tame-tails-parent-only-rails` (operator: "tame the tails then, keep parent-only rails"), under `field-strength-ruled` (planets 12-18 u/s round a median sun) and `core-at-roche-and-exaggeration-up` (core radius = Roche limit 7.32·m^0.25).
|
||||
|
||||
Scripts: `gravity-tail-taming/laws.py` (the candidate laws, force and potential), `check_laws.py` (K0/K1 accuracy, force = -dφ/dr for every law), `probe.py` (the measures), output in `probe-output.txt`. Data: `docs/research/2026-10-13-orbit-probe/data/bodies.csv` (1501 zones, 31,640 nodes), as B35/B38 used. Each law ran in 8-18 s on Starfire.
|
||||
|
||||
## The law today and the problem
|
||||
|
||||
`GravityField` (Assets/Scripts/ServerShared/GravityField.cs at c9815379): pull c/r beyond the core radius R, c·r/R² inside, potential c·(ln(r/R) + 1/2) beyond and c·r²/(2R²) inside. Map B19: c = 45.08·m^0.2. `Gravity.cginc` is generated from it and evaluates `gravityPotential` once per well per texel (one `log`).
|
||||
|
||||
The log potential is harmonic in the plane, so 2D Gauss's law holds: the phase-averaged pull at distance d from a parent is the strength enclosed within d, and siblings outside cancel (map B37, B38). Every sibling inside the child's distance counts in full at any range, which is the funnel (B21) and the 1.2-2x ratio (B38).
|
||||
|
||||
## Question 1: screened and truncated potentials
|
||||
|
||||
**2D screened Poisson (Yukawa, modified Helmholtz).** The Green's function of ∇²φ − φ/λ² = source in the plane is (1/2π)·K0(r/λ). K0(x) ≈ −ln(x/2) − γ for small x and ≈ sqrt(π/2x)·e^(−x) for large x (Abramowitz & Stegun 9.6.8, 9.7.2). The pull, (c/λ)·K1(r/λ), is c/r near and exponentially small beyond λ. Two physical systems obey exactly this law. Abrikosov vortices in a London superconductor interact through K0(r/λ_L) (the London penetration depth). Quasi-geostrophic point vortices in the equivalent-barotropic model interact through K0(r/L_d) (the Rossby deformation radius). Gauss's law fails, so a ring no longer acts as its centre. Graf's addition theorem for modified Bessel functions (DLMF §10.44) gives the ring average of K0 over a circle of radius s, seen from distance d. It is K0(d/λ)·I0(s/λ) outside the ring (d > s). Inside the ring (d < s) it is I0(d/λ)·K0(s/λ), which pulls outward: weak, but not zero as under the log law.
|
||||
|
||||
**Softening (changes the core, not the tail).** Plummer softening φ = −Gm/sqrt(r² + ε²) and spline kernels (Monaghan & Lattanzio 1985, used in GADGET-2: exactly Newtonian beyond 2.8ε) only reshape the inside. Aetheria's harmonic core at the Roche limit is of this family. None of them changes the tail, so none changes the enclosed/parent ratio.
|
||||
|
||||
**Split, not screening.** TreePM codes split the Newtonian potential: a short-range part, φ·erfc(r/2r_s), is summed pairwise and truncated at about 4.5 r_s, and the long-range part is solved on a mesh (Bagla 2002; Springel 2005, GADGET-2, MNRAS 364:1105). The total force stays Newtonian. The split exists for speed and does not change the physics, so it offers no precedent for weakening the far field.
|
||||
|
||||
**Molecular-dynamics cutoffs** (Allen & Tildesley, *Computer Simulation of Liquids*; Steinbach & Brooks 1994, J. Comput. Chem. 15:667, DOI 10.1002/jcc.540150702):
|
||||
- *Hard truncation*: the force jumps at r_c, which gives impulsive kicks and energy drift under any integrator.
|
||||
- *Shifted potential*: the energy is continuous but the force still jumps.
|
||||
- *Shifted force*: the force goes linearly to zero at r_c. The force is continuous, but the potential is altered at every r < r_c.
|
||||
- *Potential switching* (polynomial S(r) on the energy over [r_on, r_c]): leaves a force "bump" larger than the unswitched force inside the window.
|
||||
- *Force switching*: S(r) applies to the force, and a constant keeps the potential continuous. Steinbach & Brooks rate force switching and force shifting the best atom-based cutoffs, and add that a longer cutoff matters more than the choice of method.
|
||||
- Damped shifted force (Fennell & Gezelter 2006, J. Chem. Phys. 124:234104) is the electrostatics descendant: an erfc damping plus a force shift.
|
||||
|
||||
All of these are conservative pair potentials, so a potential exists and a drawn grid can show it exactly. What they trade away is accuracy beyond r_c, which here is the point.
|
||||
|
||||
**Orbits.** Only the Kepler and harmonic potentials have all bound orbits closed (Bertrand's theorem). The log law already precesses (rosettes, docs/research/rosette-relative-equilibria.md), and every candidate below does too. In the probe, the cap and switch laws, being C0 and C1 in force, put a slope break on the drawn grid at L (and at 2L for the switch). The switch also leaves the force exactly zero beyond 2L, so the grid there is flat and a child out there has no rail pull (measured below).
|
||||
|
||||
**λ per body.** The laws were probed with a per-body screening length λ = k·D (D = the body's distance from its parent, the zone radius Z for a root) and λ = 2·d_Hill (B25's tidal radius D·sqrt(c/2c_parent); falls back to D without a body parent). Since d_Hill ≈ 0.5·D for a planet round a median sun, λ = 2·d_Hill and λ = 1·D come out within 0.03 of each other on every measure (table). A body's own children (inside half its Hill radius, B25) sit at r/λ ≤ ~0.25 for λ = 1·D, where K1 still gives the full c/r (x·K1(x) = 0.94 at x = 0.25, 0.60 at x = 1, 0.28 at x = 2). Its siblings sit at r/λ ≈ 1-2, where their pull is cut to 0.60-0.28 of c/r and falls off as e^(−x) beyond. Because λ differs per body, the summed field is not the solution of one screened Poisson equation. It is a sum of per-body kernels, still a conservative potential.
|
||||
|
||||
## Question 2: the hyperbolic-plane Green's function
|
||||
|
||||
Confirmed, with a caveat. On the hyperbolic plane (curvature −1) the Green's function (1/2π)·ln coth(r/2) is harmonic away from the source, and the divergence theorem holds on the manifold. The flux through a geodesic circle is the enclosed strength, and the circumference is 2π·sinh r, so the field is c_enc/sinh r: a ring acts as its centre outside and pulls nothing inside. The 1/sinh factor is common to parent and enclosed siblings, so the enclosed/parent ratio is the flat one; it decays every body equally. Caveat: Aetheria's positions are Euclidean. On flat distances the kernel is not harmonic: its flat Laplacian is (r·cosh r − sinh r)/(r·sinh² r) > 0, a distributed screening halo. So used that way it behaves as a screened law with one global length, asymptotically e^(−r/ℓ). Probe [10] (ℓ = Z/2): the suns' children ratio p50 moves only from 2.01 to 1.98, because ℓ is much larger than sibling separations. The mid-zone inward pull falls by about 20%.
|
||||
|
||||
## Question 3: games that sum (or refuse to sum) many wells
|
||||
|
||||
| Game | How local dominance is kept | Artefact accepted | Source |
|
||||
|---|---|---|---|
|
||||
| Kerbal Space Program | Patched conics: exactly one attractor inside its sphere of influence, r_SOI = a·(m/M)^(2/5) (Laplace) | No Lagrange points, no perturbations; trajectories kink at SOI boundaries; Keplerian "on rails" outside physics range | KSP forum "on rails" threads; patched-conic literature (Wikibooks Astrodynamics) |
|
||||
| Outer Wilds | Each body is a gravity volume; planets use **linear** falloff, the sun and some moons inverse-square; planets feel only the sun; a moon feels only its planet plus an added acceleration matching the planet's orbit; volume priorities override rather than sum | Planets don't perturb each other (developer Logan Ver Hoef: real mutual pull would make them "escape orbit"); non-physical falloff chosen for feel | New Horizons mod schema `GravityVolumeInfo` (fields `fallOff`, `upperRadius`, `priority`, `layer`); PC Gamer "Recreating gravity in games is harder than you'd think" |
|
||||
| Super Mario Galaxy 1/2 | Gravity areas (sphere, plane, cylinder...) with priorities; the highest-priority area wins; equal priorities are combined | Gravity direction switches discontinuously at area boundaries; no superposition across priorities | Luma's Workshop wiki "Gravity"; Super Mario Wiki "Gravity field" |
|
||||
| Universe Sandbox | Full N-body integration, selectable integrators | No local dominance: true superposition, with integrator drift | Universe Sandbox wiki "N-Body Simulation" |
|
||||
| Elite Dangerous | Orbits set by the Stellar Forge generator; no source found for runtime propagation (forum claims only) | Unverified | Frontier forum thread; PC Gamer Stellar Forge article |
|
||||
| Children of a Dead Earth | No source found for n-body versus patched conics | Unverified | none found |
|
||||
| Spore (space stage) | No source found for any summed gravity | Unverified | none found |
|
||||
| Geometry Wars warp grid | Spring-mass grid displaced by short-range forces from black holes, bullets and respawns; not a potential; border anchored, damped | The grid is a presentation, not the field ships feel; effects are local by construction | Tuts+ "Make a Neon Vector Shooter: The Warping Grid"; gamedev.net "Fluid Grid" thread |
|
||||
|
||||
## Question 4: the probe
|
||||
|
||||
**Measures** (`probe.py`), over 6 random-phase placements per zone:
|
||||
- *Ratio*: the radial pull toward the parent at the child's position, phase-averaged. It is computed from every body outside the child's own family, minus the field at the parent from every body outside the parent's family (the parent rides its own rail), so tides and outer siblings count. It is reported as a circular-speed ratio sqrt(a_rel / a_parent_own), by the child's class, p10/p50/p90/p99. Unlike B38's Gauss-only figure it can fall below 1.
|
||||
- *Mid-zone*: the inward radial pull at 8 points at 0.5·Z (B21's funnel measure), median over zones, by Z.
|
||||
- *Planets round a median sun*: the parent-only rail speed sqrt(a·d) for planets whose sun has a mass between the sun p40 and p60 (10607-20254, n = 1116).
|
||||
|
||||
Every law shares c (calibrated so that c_sun = 304.9 at the median sun mass 14152) and the harmonic core at R = 7.32·m^0.25. The baseline (law 0) reproduces B38 (suns' children p50 2.01 against B38's 1.98).
|
||||
|
||||
Comparability with B21: B21 measured total |a| (not inward) with core 44.9·m^0.25, and its alternatives at p = 0.175. Here every row is inward pull at p = 0.2 with the Roche core.
|
||||
|
||||
| # | Law | Ratio p50 / p90: Sun · GasGiant · Planet · Planetoid · Belt | Mid-zone inward pull, Z = 1000 / 2125 / 4182 / 6846 / 10000 (u/s²) | Planets round a median sun, v p10/p50/p90 (u/s); share in 12-18 | Children left with no parent pull | HLSL cost per well per texel |
|
||||
|---|---|---|---|---|---|---|
|
||||
| 0 | log (today) | 2.01/2.46 · 1.62/2.21 · 1.24/1.71 · 1.11/1.43 · 1.18/1.55 | 1.59 / 1.27 / 0.91 / 0.67 / 0.51 | 17.1/17.5/18.0; 0.92 | 0 | 1 log |
|
||||
| 1 | K0, λ = 0.5·D | 1.15/1.50 · 1.04/1.37 · 0.99/1.19 · 0.99/1.14 · 0.99/1.20 | 0.42 / 0.29 / 0.18 / 0.13 / 0.09 | 15.3/16.7/17.5; 0.99 | 0 | K0 approx (below) |
|
||||
| 2 | K0, λ = 1·D | 1.28/1.56 · 1.14/1.46 · 1.03/1.29 · 1.01/1.21 · 1.03/1.27 | 0.73 / 0.52 / 0.32 / 0.22 / 0.16 | 16.6/17.2/17.7; 0.99 | 0 | K0 approx |
|
||||
| 3 | K0, λ = 2·d_Hill | 1.25/1.54 · 1.12/1.46 · 1.02/1.27 · 1.00/1.21 · 1.02/1.25 | 0.68 / 0.48 / 0.31 / 0.22 / 0.16 | 16.6/17.2/17.7; 0.99 | 0 | K0 approx |
|
||||
| 4 | cap c·L/r² beyond L = 100·m^0.25 | 1.37/1.74 · 1.21/1.64 · 1.08/1.42 · 1.04/1.32 · 1.08/1.40 | 1.06 / 0.70 / 0.39 / 0.19 / 0.11 | 16.4/17.4/18.0; 0.93 | 0 | 1 log or 1 div, branch |
|
||||
| 5 | cap c·L/r², L = 300·m^0.25 | 1.50/1.84 · 1.41/1.92 · 1.18/1.64 · 1.08/1.41 · 1.13/1.51 | 1.52 / 1.09 / 0.63 / 0.38 / 0.25 | 17.1/17.5/18.0; 0.92 | 0 | same |
|
||||
| 6 | force switch, exact to L, zero at 2L, L = 100·m^0.25 | 1.02/1.58 · 1.04/1.49 · 1.01/1.30 · 1.01/1.27 · 1.02/1.33 (Sun row over 1443 of 2211) | 0.78 / 0.49 / 0.27 / 0.10 / 0.01 | 16.9/17.4/18.0; 0.90 | 962 (Sun 768, GasGiant 75, Planetoid 78, Belt 33, Planet 8) | 1 log + cubic, branch |
|
||||
| 7 | force switch, L = 300·m^0.25 | 1.23/1.61 · 1.28/1.80 · 1.15/1.61 · 1.08/1.40 · 1.12/1.48 | 1.50 / 1.01 / 0.48 / 0.26 / 0.16 | 17.1/17.5/18.0; 0.92 | 0 | same |
|
||||
| 8 | uniform-disc subtraction (zone total c over Z) | 1.33/1.77 · 1.31/1.80 · 1.14/1.48 · 1.07/1.38 · 1.13/1.47 | 1.00 / 0.77 / 0.55 / 0.42 / 0.32 | 17.1/17.5/18.0; 0.92 | 0 | 1 log; plus one disc term per texel |
|
||||
| 9 | exponent p = 0.5 (recalibrated at the median sun) | 1.24/1.40 · 1.14/1.37 · 1.05/1.24 · 1.03/1.26 · 1.05/1.30 | 0.67 / 0.87 / 0.73 / 0.61 / 0.51 | 16.5/17.5/18.8; 0.68 | 0 | 1 log |
|
||||
| 10 | hyperbolic kernel on flat distances, ℓ = Z/2 | 1.98/2.50 · 1.59/2.22 · 1.21/1.69 · 1.09/1.41 · 1.15/1.51 | 1.28 / 1.03 / 0.72 / 0.54 / 0.41 | 17.0/17.4/17.9; 0.96 | 0 | 1 log + 1 exp |
|
||||
|
||||
**Other numbers from `probe-output.txt`:**
|
||||
- *Parent-only rail speed of the suns' children* (children of root and rosette suns), p10/p50/p90 u/s: today 22.1/23.6/24.6. K0 λ = 0.5·D gives 14.3/17.3/20.5, K0 λ = 1·D gives 18.7/20.8/22.9, and cap L = 100 gives 14.5/17.8/22.2. Screening lowers the parent's own pull at the outer children of a root, so their rails slow.
|
||||
- *Planetoids' parent-only speed* under K0 λ = 0.5·D: p10 9.2 against 11.5 today.
|
||||
- *p = 0.5*: suns' children run 31-41 u/s, planetoids p10 6.2 u/s, and planets round all suns reach p90 37.3 u/s.
|
||||
- *Force switch at L = 100·m^0.25*: 768 of the 2211 children of suns (35%) sit beyond their parent's 2L, where the parent pulls nothing, so a parent-only rail has no period (B35's "no field, no motion").
|
||||
- *Planets round all suns*, v p50: 16.9-17.9 across every law except p = 0.5 (18.5). The share in 12-18 is 0.47-0.57 for every law, today included, because heavy root suns put p90 near 23.6 u/s.
|
||||
|
||||
**Grid equals field.** Every candidate has a closed-form potential, and `check_laws.py` confirms force = −dφ/dr to the finite-difference floor (4.3e-5 relative) for each one. The potentials (tail part, r ≥ R, with the harmonic core matched in value and slope at R as today):
|
||||
- K0: −c·K0(r/λ).
|
||||
- Cap: c·ln r inside L, and c·(ln L + 1 − L/r) beyond.
|
||||
- Force switch: c·[s0·ln r + s1·r + s2·r²/2 + s3·r³/3] on [L, 2L], with the s_i the switch polynomial's coefficients (`laws.py: phi_switch`), and a constant beyond 2L.
|
||||
- Disc: the zone term C·r²/(2Z²) inside Z.
|
||||
- Hyperbolic: −c·ln coth(r/2ℓ).
|
||||
|
||||
Each fits the per-well reference subtraction `Gravity.cginc` already does (B22). Matching the core needs F_tail(R) per well, which is one more per-well constant in the `GravityWell` buffer. The same holds for λ or L (a `float` per well).
|
||||
|
||||
**K0 in HLSL.** Abramowitz & Stegun 9.8.5-9.8.6 give K0, and 9.8.7-9.8.8 give K1, as two branches:
|
||||
- x ≤ 2: a 7-term polynomial in (x/2)², plus −ln(x/2)·I0(x), with I0 from A&S 9.8.1 (6 terms).
|
||||
- x > 2: e^(−x)/sqrt(x) times a 7-term polynomial in 2/x.
|
||||
|
||||
A&S state |ε| < 1e-8 (9.8.5) and < 1.9e-7 (9.8.6). Against direct quadrature, `check_laws.py` measures 5.8e-8 relative for K0 and 2.6e-8 for K1, in float64. In float32 the shader's precision floor dominates. Per well per texel the cost is about 14 FMAs plus `log` or `exp`+`rsqrt` (the same routine is Numerical Recipes' `bessk0`/`bessk1`). Written branchless (both halves evaluated, one selected) it is about double. A 1D lookup texture in x = r/λ is the common alternative (one sample per well). The grid only needs the potential, K0; the force path needs K1.
|
||||
|
||||
## Confidence
|
||||
|
||||
- Q1: high on the math and the MD/N-body facts; medium on the exact DLMF section for Graf's theorem.
|
||||
- Q2: high.
|
||||
- Q3: high for KSP, Outer Wilds (schema plus developer quotes) and Universe Sandbox; medium for Mario Galaxy (community wiki marked in progress); low or unverified for Elite, CoaDE and Spore.
|
||||
- Q4: high that the probe is internally consistent (the baseline reproduces B38). Medium on absolute ratios: 6 phase draws, phase-averaged circular speeds rather than integrated free orbits, and positions from the 2026-10-13 orbit probe's generator (B35: unchanged since c8d40eaf).
|
||||
|
||||
## Sources
|
||||
|
||||
- Abramowitz & Stegun, *Handbook of Mathematical Functions*, §9.6-9.8; NIST DLMF ch. 10 (§10.44 addition theorems).
|
||||
- Steinbach & Brooks 1994, J. Comput. Chem. 15(7):667-683, https://api.crossref.org/works/10.1002%2FJCC.540150702 ; force switch and bump discussion: https://arxiv.org/html/2410.14187v2
|
||||
- Fennell & Gezelter 2006, J. Chem. Phys. 124:234104. Allen & Tildesley, *Computer Simulation of Liquids*. Springel 2005, MNRAS 364:1105 (GADGET-2). Bagla 2002, J. Astrophys. Astron. 23:185 (TreePM). Monaghan & Lattanzio 1985, A&A 149:135.
|
||||
- Outer Wilds: https://nh.outerwildsmods.com/schemas/body-schema/defs/gravityvolumeinfo ; https://pcgamer.com/how-to-build-a-black-hole ; https://thephysicsmill.com/2024/09/20/an-astrophysicist-attempts-to-measure-the-physics-of-outer-wilds/
|
||||
- Super Mario Galaxy: https://lumasworkshop.com/wiki/Gravity ; https://mariowiki.com/Gravity_field
|
||||
- KSP: https://forum.kerbalspaceprogram.com/topic/123592-a-clarification-of-quoton-railsquot ; https://en.wikibooks.org/wiki/Astrodynamics/Interplanetary_Trajectories
|
||||
- Universe Sandbox: https://universesandbox.gamepedia.com/N-Body_Simulation
|
||||
- Elite Dangerous: https://forums.frontier.co.uk/threads/are-the-orbits-of-the-planets-physically-correct.394502/post-6203566 ; https://www.pcgamer.com/uk/the-mind-bending-science-behind-the-planets-of-elite-dangerous/
|
||||
- Geometry Wars grid: https://code.tutsplus.com/make-a-neon-vector-shooter-in-xna-the-warping-grid--gamedev-9904t ; https://www.gamedev.net/forums/topic/555317-fluid-grid-aka-geometry-wars-background/4569783/
|
||||
14
docs/research/gravity-tail-taming/check_laws.py
Normal file
14
docs/research/gravity-tail-taming/check_laws.py
Normal file
|
|
@ -0,0 +1,14 @@
|
|||
# Checks: A&S K0/K1 against a direct integral; each law's force equals -dphi/dr (grid equals field).
|
||||
import numpy as np, sys, os
|
||||
sys.path.insert(0, os.path.dirname(__file__))
|
||||
from laws import *
|
||||
x = np.array([0.05, 0.3, 1.0, 1.99, 2.01, 4.0, 10.0])
|
||||
t = np.linspace(0, 12, 200001)
|
||||
ref0 = np.array([np.trapezoid(np.exp(-xi * np.cosh(t)), t) for xi in x])
|
||||
ref1 = np.array([np.trapezoid(np.exp(-xi * np.cosh(t)) * np.cosh(t), t) for xi in x])
|
||||
print('K0 rel err', np.max(np.abs(k0(x) / ref0 - 1)), ' K1 rel err', np.max(np.abs(k1(x) / ref1 - 1)))
|
||||
r = np.linspace(20, 5000, 20000); c = 300.0
|
||||
for name, tf, pf, s in [('log', tail_log, phi_log, 0), ('k0', tail_k0, phi_k0, 800.0), ('cap', tail_cap, phi_cap, 1000.0),
|
||||
('switch', tail_switch, phi_switch, 1000.0), ('hyp', tail_hyp, phi_hyp, 1500.0)]:
|
||||
ph = pf(c, r, s); dphi = np.gradient(ph, r); F = tf(c, r, s)
|
||||
print(f'{name:7s} max |F - dphi/dr| / max F = {np.max(np.abs(F - dphi)[5:-5]) / F.max():.2e}')
|
||||
79
docs/research/gravity-tail-taming/laws.py
Normal file
79
docs/research/gravity-tail-taming/laws.py
Normal file
|
|
@ -0,0 +1,79 @@
|
|||
# Candidate one-well laws for tail taming (aetheria-gravity, Eyes pass eyes-ag3-tail-taming, 2026-10-15).
|
||||
# Each law gives the pull F(r) toward the body (u/s^2) and a potential phi(r) with -dphi/dr = -F (phi up to a constant).
|
||||
# Shared: c = G*m^p, core radius R = Roche = 7.32*m^0.25 with a harmonic core F = F_tail(R)*r/R inside R (GravityField at 09ec04a8+).
|
||||
# K0/K1 use Abramowitz & Stegun 9.8.5-9.8.8 polynomial approximations, the same ones an HLSL port would use.
|
||||
import numpy as np
|
||||
|
||||
G0, P0, RO = 45.08, 0.2, 7.32
|
||||
CSUN, MSUN = 0.5 * 609.8, 14152.4 # field-strength-ruled calibration (B19)
|
||||
|
||||
def strength(m, p=P0):
|
||||
G = CSUN / MSUN ** p
|
||||
return G * m ** p
|
||||
|
||||
def core_radius(m):
|
||||
return RO * m ** 0.25
|
||||
|
||||
# ---- modified Bessel functions, A&S 9.8.1-9.8.8 ----
|
||||
def _i0(x):
|
||||
t = (x / 3.75) ** 2
|
||||
small = 1 + t*(3.5156229 + t*(3.0899424 + t*(1.2067492 + t*(0.2659732 + t*(0.0360768 + t*0.0045813)))))
|
||||
return small # valid |x| <= 3.75; only used for x <= 2
|
||||
|
||||
def _i1(x):
|
||||
t = (x / 3.75) ** 2
|
||||
return x * (0.5 + t*(0.87890594 + t*(0.51498869 + t*(0.15084934 + t*(0.02658733 + t*(0.00301532 + t*0.00032411))))))
|
||||
|
||||
def k0(x):
|
||||
x = np.asarray(x, float)
|
||||
xs = np.minimum(x, 2.0); y = (xs / 2) ** 2
|
||||
lo = -np.log(xs / 2) * _i0(xs) + (-0.57721566 + y*(0.42278420 + y*(0.23069756 + y*(0.03488590 + y*(0.00262698 + y*(0.00010750 + y*0.00000740))))))
|
||||
xl = np.maximum(x, 2.0); z = 2 / xl
|
||||
hi = np.exp(-xl) / np.sqrt(xl) * (1.25331414 + z*(-0.07832358 + z*(0.02189568 + z*(-0.01062446 + z*(0.00587872 + z*(-0.00251540 + z*0.00053208))))))
|
||||
return np.where(x <= 2, lo, hi)
|
||||
|
||||
def k1(x):
|
||||
x = np.asarray(x, float)
|
||||
xs = np.minimum(x, 2.0); y = (xs / 2) ** 2
|
||||
lo = np.log(xs / 2) * _i1(xs) + (1 / xs) * (1 + y*(0.15443144 + y*(-0.67278579 + y*(-0.18156897 + y*(-0.01919402 + y*(-0.00110404 + y*(-0.00004686)))))))
|
||||
xl = np.maximum(x, 2.0); z = 2 / xl
|
||||
hi = np.exp(-xl) / np.sqrt(xl) * (1.25331414 + z*(0.23498619 + z*(-0.03655620 + z*(0.01504268 + z*(-0.00780353 + z*(0.00325614 + z*(-0.00068245)))))))
|
||||
return np.where(x <= 2, lo, hi)
|
||||
|
||||
# ---- tail laws: F_tail(c, r, s) with s a per-body scale parameter (array broadcastable) ----
|
||||
def tail_log(c, r, s): # today: F = c/r
|
||||
return c / r
|
||||
def tail_k0(c, r, lam): # 2D screened Poisson: phi = -c K0(r/lam), F = (c/lam) K1(r/lam) -> c/r near
|
||||
return c / lam * k1(r / lam)
|
||||
def tail_cap(c, r, L): # B21's cap: c/r inside L, c L / r^2 beyond (force continuous at L)
|
||||
return np.where(r < L, c / r, c * L / r ** 2)
|
||||
def _switch(r, a, b): # force switch 1 -> 0 over [a, b], cubic smoothstep (C1)
|
||||
u = np.clip((r - a) / (b - a), 0, 1)
|
||||
return 1 - u * u * (3 - 2 * u)
|
||||
def tail_switch(c, r, L): # force-switched cutoff: exact c/r to L, smooth to zero at 2L
|
||||
return c / r * _switch(r, L, 2 * L)
|
||||
def tail_hyp(c, r, ell): # hyperbolic-plane Green's function used on flat distances: F = c/(ell sinh(r/ell))
|
||||
return c / (ell * np.sinh(np.minimum(r / ell, 80)))
|
||||
|
||||
def pull(tail, c, R, r, s):
|
||||
"""Pull toward the body at distance r: tail outside the core, linear inside, continuous at R."""
|
||||
r = np.asarray(r, float)
|
||||
rr = np.maximum(r, R)
|
||||
ft = tail(c, rr, s)
|
||||
return np.where(r < R, ft * r / R, ft)
|
||||
|
||||
# Potentials (closed forms, for the HLSL statement and the grid-equals-field check), tail part only, r >= R.
|
||||
def phi_log(c, r, s): return c * np.log(r)
|
||||
def phi_k0(c, r, lam): return -c * k0(r / lam)
|
||||
def phi_cap(c, r, L): return np.where(r < L, c * np.log(r), c * (np.log(L) + 1 - L / r))
|
||||
def phi_hyp(c, r, ell): return -c * np.log(1 / np.tanh(np.minimum(r / (2 * ell), 40)))
|
||||
def phi_switch(c, r, L):
|
||||
# integral of (c/r) * (1 - 3u^2 + 2u^3), u = (r-L)/L on [L, 2L]: expand in r, integrate term by term
|
||||
a, b = L, 2 * L
|
||||
def inner(x): # antiderivative of S(x)/x on [a, b]
|
||||
u = (x - a) / (b - a); d = b - a
|
||||
# S = 1 - 3u^2 + 2u^3 with u = (x-a)/d; expand to polynomial in x: S = s0 + s1 x + s2 x^2 + s3 x^3
|
||||
s3 = 2 / d**3; s2 = -3 / d**2 - 6 * a / d**3; s1 = 6 * a / d**2 + 6 * a * a / d**3; s0 = 1 - 3 * a * a / d**2 - 2 * a**3 / d**3
|
||||
return s0 * np.log(x) + s1 * x + s2 * x * x / 2 + s3 * x**3 / 3
|
||||
rc = np.clip(r, a, b)
|
||||
return c * np.where(r < a, np.log(np.maximum(r, 1e-9)), np.log(a) + inner(rc) - inner(a))
|
||||
110
docs/research/gravity-tail-taming/probe-output.txt
Normal file
110
docs/research/gravity-tail-taming/probe-output.txt
Normal file
|
|
@ -0,0 +1,110 @@
|
|||
== [0] log (today) (11s, 6 draws)
|
||||
Sun ratio p10/50/90/99 [1.44, 2.01, 2.46, 2.81] parent-only v p10/50/90 [22.06, 23.62, 24.64] n 2211
|
||||
GasGiant ratio p10/50/90/99 [1.1, 1.62, 2.21, 2.67] parent-only v p10/50/90 [18.27, 22.2, 24.57] n 3566
|
||||
Planet ratio p10/50/90/99 [0.97, 1.24, 1.71, 2.16] parent-only v p10/50/90 [15.56, 17.38, 23.61] n 7445
|
||||
Planetoid ratio p10/50/90/99 [0.95, 1.11, 1.43, 1.76] parent-only v p10/50/90 [11.51, 14.48, 20.27] n 10200
|
||||
Belt ratio p10/50/90/99 [0.95, 1.18, 1.55, 2.01] parent-only v p10/50/90 [11.41, 14.74, 20.29] n 3134
|
||||
planets round suns: parent-only v p10/50/90 [16.3, 17.88, 23.65]; share in 12-18 u/s 0.51
|
||||
planets round a median sun (mass p40-p60, 10607-20254): v p10/50/90 [17.06, 17.48, 17.98]; share in 12-18 0.92; n 1116
|
||||
children with no parent pull at their distance: {}
|
||||
mid-zone inward pull by Z: 1000: 1.59, 2125: 1.27, 4182: 0.91, 6846: 0.67, 10000: 0.51
|
||||
== [1] K0 lam=0.5*D (18s, 6 draws)
|
||||
Sun ratio p10/50/90/99 [0.88, 1.15, 1.5, 1.81] parent-only v p10/50/90 [14.32, 17.26, 20.53] n 2211
|
||||
GasGiant ratio p10/50/90/99 [0.83, 1.04, 1.37, 1.84] parent-only v p10/50/90 [14.37, 18.21, 22.01] n 3566
|
||||
Planet ratio p10/50/90/99 [0.89, 0.99, 1.19, 1.6] parent-only v p10/50/90 [12.7, 16.4, 22.68] n 7445
|
||||
Planetoid ratio p10/50/90/99 [0.86, 0.99, 1.14, 1.81] parent-only v p10/50/90 [9.21, 13.52, 19.8] n 10200
|
||||
Belt ratio p10/50/90/99 [0.84, 0.99, 1.2, 1.87] parent-only v p10/50/90 [8.79, 13.42, 19.6] n 3134
|
||||
planets round suns: parent-only v p10/50/90 [12.99, 16.87, 22.94]; share in 12-18 u/s 0.57
|
||||
planets round a median sun (mass p40-p60, 10607-20254): v p10/50/90 [15.31, 16.68, 17.47]; share in 12-18 0.99; n 1116
|
||||
children with no parent pull at their distance: {}
|
||||
mid-zone inward pull by Z: 1000: 0.42, 2125: 0.29, 4182: 0.18, 6846: 0.13, 10000: 0.09
|
||||
== [2] K0 lam=1*D (12s, 6 draws)
|
||||
Sun ratio p10/50/90/99 [1.0, 1.28, 1.56, 1.8] parent-only v p10/50/90 [18.72, 20.77, 22.92] n 2211
|
||||
GasGiant ratio p10/50/90/99 [0.93, 1.14, 1.46, 1.87] parent-only v p10/50/90 [17.12, 20.03, 23.43] n 3566
|
||||
Planet ratio p10/50/90/99 [0.92, 1.03, 1.29, 1.65] parent-only v p10/50/90 [14.59, 17.03, 23.34] n 7445
|
||||
Planetoid ratio p10/50/90/99 [0.89, 1.01, 1.21, 1.62] parent-only v p10/50/90 [10.61, 14.09, 20.1] n 10200
|
||||
Belt ratio p10/50/90/99 [0.89, 1.03, 1.27, 1.66] parent-only v p10/50/90 [10.42, 14.27, 20.05] n 3134
|
||||
planets round suns: parent-only v p10/50/90 [15.13, 17.56, 23.45]; share in 12-18 u/s 0.55
|
||||
planets round a median sun (mass p40-p60, 10607-20254): v p10/50/90 [16.56, 17.16, 17.72]; share in 12-18 0.99; n 1116
|
||||
children with no parent pull at their distance: {}
|
||||
mid-zone inward pull by Z: 1000: 0.73, 2125: 0.52, 4182: 0.32, 6846: 0.22, 10000: 0.16
|
||||
== [3] K0 lam=2*d_Hill (12s, 6 draws)
|
||||
Sun ratio p10/50/90/99 [0.98, 1.25, 1.54, 1.83] parent-only v p10/50/90 [18.72, 20.77, 22.92] n 2211
|
||||
GasGiant ratio p10/50/90/99 [0.91, 1.12, 1.46, 1.85] parent-only v p10/50/90 [17.12, 20.03, 23.43] n 3566
|
||||
Planet ratio p10/50/90/99 [0.9, 1.02, 1.27, 1.63] parent-only v p10/50/90 [14.59, 17.06, 23.34] n 7445
|
||||
Planetoid ratio p10/50/90/99 [0.89, 1.0, 1.21, 1.63] parent-only v p10/50/90 [10.48, 14.06, 20.1] n 10200
|
||||
Belt ratio p10/50/90/99 [0.89, 1.02, 1.25, 1.7] parent-only v p10/50/90 [10.25, 14.24, 20.05] n 3134
|
||||
planets round suns: parent-only v p10/50/90 [15.16, 17.57, 23.45]; share in 12-18 u/s 0.55
|
||||
planets round a median sun (mass p40-p60, 10607-20254): v p10/50/90 [16.58, 17.18, 17.74]; share in 12-18 0.99; n 1116
|
||||
children with no parent pull at their distance: {}
|
||||
mid-zone inward pull by Z: 1000: 0.68, 2125: 0.48, 4182: 0.31, 6846: 0.22, 10000: 0.16
|
||||
== [4] cap c*L/r^2, L=100m^.25 (9s, 6 draws)
|
||||
Sun ratio p10/50/90/99 [1.02, 1.37, 1.74, 2.18] parent-only v p10/50/90 [14.47, 17.81, 22.17] n 2211
|
||||
GasGiant ratio p10/50/90/99 [0.94, 1.21, 1.64, 2.15] parent-only v p10/50/90 [17.78, 20.24, 23.64] n 3566
|
||||
Planet ratio p10/50/90/99 [0.92, 1.08, 1.42, 1.84] parent-only v p10/50/90 [15.37, 17.21, 23.59] n 7445
|
||||
Planetoid ratio p10/50/90/99 [0.91, 1.04, 1.32, 1.66] parent-only v p10/50/90 [11.28, 14.43, 20.27] n 10200
|
||||
Belt ratio p10/50/90/99 [0.91, 1.08, 1.4, 1.74] parent-only v p10/50/90 [11.15, 14.64, 20.29] n 3134
|
||||
planets round suns: parent-only v p10/50/90 [16.12, 17.77, 23.64]; share in 12-18 u/s 0.52
|
||||
planets round a median sun (mass p40-p60, 10607-20254): v p10/50/90 [16.35, 17.39, 17.95]; share in 12-18 0.93; n 1116
|
||||
children with no parent pull at their distance: {}
|
||||
mid-zone inward pull by Z: 1000: 1.06, 2125: 0.70, 4182: 0.39, 6846: 0.19, 10000: 0.11
|
||||
== [5] cap c*L/r^2, L=300m^.25 (9s, 6 draws)
|
||||
Sun ratio p10/50/90/99 [1.14, 1.5, 1.84, 2.18] parent-only v p10/50/90 [22.06, 23.62, 24.64] n 2211
|
||||
GasGiant ratio p10/50/90/99 [1.0, 1.41, 1.92, 2.4] parent-only v p10/50/90 [18.27, 22.2, 24.57] n 3566
|
||||
Planet ratio p10/50/90/99 [0.95, 1.18, 1.64, 2.06] parent-only v p10/50/90 [15.56, 17.38, 23.61] n 7445
|
||||
Planetoid ratio p10/50/90/99 [0.93, 1.08, 1.41, 1.77] parent-only v p10/50/90 [11.51, 14.48, 20.27] n 10200
|
||||
Belt ratio p10/50/90/99 [0.94, 1.13, 1.51, 1.91] parent-only v p10/50/90 [11.41, 14.74, 20.29] n 3134
|
||||
planets round suns: parent-only v p10/50/90 [16.3, 17.88, 23.65]; share in 12-18 u/s 0.51
|
||||
planets round a median sun (mass p40-p60, 10607-20254): v p10/50/90 [17.06, 17.48, 17.98]; share in 12-18 0.92; n 1116
|
||||
children with no parent pull at their distance: {}
|
||||
mid-zone inward pull by Z: 1000: 1.52, 2125: 1.09, 4182: 0.63, 6846: 0.38, 10000: 0.25
|
||||
== [6] switch L..2L, L=100m^.25 (10s, 6 draws)
|
||||
Sun ratio p10/50/90/99 [0.66, 1.02, 1.58, 4.99] parent-only v p10/50/90 [0.0, 13.03, 22.28] n 1443
|
||||
GasGiant ratio p10/50/90/99 [0.83, 1.04, 1.49, 2.1] parent-only v p10/50/90 [16.47, 20.26, 23.68] n 3491
|
||||
Planet ratio p10/50/90/99 [0.87, 1.01, 1.3, 1.77] parent-only v p10/50/90 [15.39, 17.25, 23.6] n 7437
|
||||
Planetoid ratio p10/50/90/99 [0.87, 1.01, 1.27, 1.62] parent-only v p10/50/90 [11.22, 14.43, 20.27] n 10122
|
||||
Belt ratio p10/50/90/99 [0.87, 1.02, 1.33, 1.63] parent-only v p10/50/90 [11.09, 14.66, 20.29] n 3101
|
||||
planets round suns: parent-only v p10/50/90 [16.15, 17.78, 23.65]; share in 12-18 u/s 0.51
|
||||
planets round a median sun (mass p40-p60, 10607-20254): v p10/50/90 [16.92, 17.4, 17.95]; share in 12-18 0.90; n 1116
|
||||
children with no parent pull at their distance: {'Planetoid': 78, 'Belt': 33, 'Planet': 8, 'GasGiant': 75, 'Sun': 768}
|
||||
mid-zone inward pull by Z: 1000: 0.78, 2125: 0.49, 4182: 0.27, 6846: 0.10, 10000: 0.01
|
||||
== [7] switch L..2L, L=300m^.25 (9s, 6 draws)
|
||||
Sun ratio p10/50/90/99 [0.93, 1.23, 1.61, 2.11] parent-only v p10/50/90 [22.06, 23.62, 24.64] n 2211
|
||||
GasGiant ratio p10/50/90/99 [0.92, 1.28, 1.8, 2.37] parent-only v p10/50/90 [18.27, 22.2, 24.57] n 3566
|
||||
Planet ratio p10/50/90/99 [0.93, 1.15, 1.61, 1.99] parent-only v p10/50/90 [15.56, 17.38, 23.61] n 7445
|
||||
Planetoid ratio p10/50/90/99 [0.93, 1.08, 1.4, 1.74] parent-only v p10/50/90 [11.51, 14.48, 20.27] n 10200
|
||||
Belt ratio p10/50/90/99 [0.92, 1.12, 1.48, 1.89] parent-only v p10/50/90 [11.41, 14.74, 20.29] n 3134
|
||||
planets round suns: parent-only v p10/50/90 [16.3, 17.88, 23.65]; share in 12-18 u/s 0.51
|
||||
planets round a median sun (mass p40-p60, 10607-20254): v p10/50/90 [17.06, 17.48, 17.98]; share in 12-18 0.92; n 1116
|
||||
children with no parent pull at their distance: {}
|
||||
mid-zone inward pull by Z: 1000: 1.50, 2125: 1.01, 4182: 0.48, 6846: 0.26, 10000: 0.16
|
||||
== [8] uniform-disc subtraction (9s, 6 draws)
|
||||
Sun ratio p10/50/90/99 [0.87, 1.33, 1.77, 2.27] parent-only v p10/50/90 [22.06, 23.62, 24.64] n 2211
|
||||
GasGiant ratio p10/50/90/99 [0.93, 1.31, 1.8, 2.48] parent-only v p10/50/90 [18.27, 22.2, 24.57] n 3566
|
||||
Planet ratio p10/50/90/99 [0.91, 1.14, 1.48, 1.93] parent-only v p10/50/90 [15.56, 17.38, 23.61] n 7445
|
||||
Planetoid ratio p10/50/90/99 [0.92, 1.07, 1.38, 1.68] parent-only v p10/50/90 [11.51, 14.48, 20.27] n 10200
|
||||
Belt ratio p10/50/90/99 [0.91, 1.13, 1.47, 1.75] parent-only v p10/50/90 [11.41, 14.74, 20.29] n 3134
|
||||
planets round suns: parent-only v p10/50/90 [16.3, 17.88, 23.65]; share in 12-18 u/s 0.51
|
||||
planets round a median sun (mass p40-p60, 10607-20254): v p10/50/90 [17.06, 17.48, 17.98]; share in 12-18 0.92; n 1116
|
||||
children with no parent pull at their distance: {}
|
||||
mid-zone inward pull by Z: 1000: 1.00, 2125: 0.77, 4182: 0.55, 6846: 0.42, 10000: 0.32
|
||||
== [9] exponent p=0.5 (8s, 6 draws)
|
||||
Sun ratio p10/50/90/99 [1.07, 1.24, 1.4, 1.65] parent-only v p10/50/90 [31.31, 37.15, 41.3] n 2211
|
||||
GasGiant ratio p10/50/90/99 [1.0, 1.14, 1.37, 1.87] parent-only v p10/50/90 [19.55, 31.81, 41.01] n 3566
|
||||
Planet ratio p10/50/90/99 [0.98, 1.05, 1.24, 1.58] parent-only v p10/50/90 [13.09, 17.28, 37.12] n 7445
|
||||
Planetoid ratio p10/50/90/99 [0.94, 1.03, 1.26, 2.32] parent-only v p10/50/90 [6.19, 11.0, 25.38] n 10200
|
||||
Belt ratio p10/50/90/99 [0.93, 1.05, 1.3, 2.33] parent-only v p10/50/90 [6.04, 11.43, 25.43] n 3134
|
||||
planets round suns: parent-only v p10/50/90 [14.72, 18.53, 37.29]; share in 12-18 u/s 0.47
|
||||
planets round a median sun (mass p40-p60, 10607-20254): v p10/50/90 [16.48, 17.5, 18.77]; share in 12-18 0.68; n 1116
|
||||
children with no parent pull at their distance: {}
|
||||
mid-zone inward pull by Z: 1000: 0.67, 2125: 0.87, 4182: 0.73, 6846: 0.61, 10000: 0.51
|
||||
== [10] hyperbolic kernel ell=Z/2 (9s, 6 draws)
|
||||
Sun ratio p10/50/90/99 [1.39, 1.98, 2.5, 2.93] parent-only v p10/50/90 [18.65, 20.75, 23.26] n 2211
|
||||
GasGiant ratio p10/50/90/99 [1.05, 1.59, 2.22, 2.7] parent-only v p10/50/90 [17.47, 20.27, 23.71] n 3566
|
||||
Planet ratio p10/50/90/99 [0.95, 1.21, 1.69, 2.17] parent-only v p10/50/90 [15.06, 17.31, 23.46] n 7445
|
||||
Planetoid ratio p10/50/90/99 [0.94, 1.09, 1.41, 1.76] parent-only v p10/50/90 [11.49, 14.45, 20.17] n 10200
|
||||
Belt ratio p10/50/90/99 [0.95, 1.15, 1.51, 1.91] parent-only v p10/50/90 [11.4, 14.66, 20.15] n 3134
|
||||
planets round suns: parent-only v p10/50/90 [15.7, 17.78, 23.53]; share in 12-18 u/s 0.53
|
||||
planets round a median sun (mass p40-p60, 10607-20254): v p10/50/90 [17.0, 17.42, 17.91]; share in 12-18 0.96; n 1116
|
||||
children with no parent pull at their distance: {}
|
||||
mid-zone inward pull by Z: 1000: 1.28, 2125: 1.03, 4182: 0.72, 6846: 0.54, 10000: 0.41
|
||||
127
docs/research/gravity-tail-taming/probe.py
Normal file
127
docs/research/gravity-tail-taming/probe.py
Normal file
|
|
@ -0,0 +1,127 @@
|
|||
# Tail-taming probe over docs/research/2026-10-13-orbit-probe/data/bodies.csv (1501 zones, 31,640 nodes).
|
||||
# For each candidate law: (1) relative radial pull on a child, phase-averaged over random phases, over its parent's own
|
||||
# pull -> circular-speed ratio sqrt(a_rel / a_parent); (2) B21's mid-zone inward pull; (3) parent-only rail speeds.
|
||||
# a_rel = radial component, toward the parent, of [field at the child from every body outside the child's family]
|
||||
# minus [field at the parent from every body outside the parent's family] (the parent rides its own rail),
|
||||
# with the disc term included where a law subtracts one.
|
||||
# Usage: python probe.py bodies.csv [draws] [law-index,...]
|
||||
import csv, sys, os, collections, time
|
||||
import numpy as np
|
||||
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
|
||||
from laws import *
|
||||
|
||||
rng = np.random.default_rng(7)
|
||||
DRAWS = int(sys.argv[2]) if len(sys.argv) > 2 else 6
|
||||
rows = [r for r in csv.DictReader(open(sys.argv[1])) if r.get('type')]
|
||||
zones = collections.defaultdict(list)
|
||||
for r in rows: zones[(r['density'], r['seed'])].append(r)
|
||||
|
||||
def hill2(D, c, cpar): # 2 x B25's tidal radius d_H = D sqrt(c / (2 c_parent)); a root or bodyless-parent child: D
|
||||
ok = cpar > 0
|
||||
return np.where(ok, 2 * D * np.sqrt(c / (2 * np.where(ok, cpar, 1.0))), D)
|
||||
|
||||
ONE = lambda m, D, Z, c, cp: np.ones_like(m)
|
||||
# (name, tail, p, scale(m, D, Z, c, c_parent) -> per-body parameter, subtract disc)
|
||||
LAWS = [
|
||||
('log (today)', tail_log, 0.2, ONE, False),
|
||||
('K0 lam=0.5*D', tail_k0, 0.2, lambda m, D, Z, c, cp: 0.5 * D, False),
|
||||
('K0 lam=1*D', tail_k0, 0.2, lambda m, D, Z, c, cp: 1.0 * D, False),
|
||||
('K0 lam=2*d_Hill', tail_k0, 0.2, lambda m, D, Z, c, cp: hill2(D, c, cp), False),
|
||||
('cap c*L/r^2, L=100m^.25', tail_cap, 0.2, lambda m, D, Z, c, cp: 100 * m ** 0.25, False),
|
||||
('cap c*L/r^2, L=300m^.25', tail_cap, 0.2, lambda m, D, Z, c, cp: 300 * m ** 0.25, False),
|
||||
('switch L..2L, L=100m^.25', tail_switch, 0.2, lambda m, D, Z, c, cp: 100 * m ** 0.25, False),
|
||||
('switch L..2L, L=300m^.25', tail_switch, 0.2, lambda m, D, Z, c, cp: 300 * m ** 0.25, False),
|
||||
('uniform-disc subtraction', tail_log, 0.2, ONE, True),
|
||||
('exponent p=0.5', tail_log, 0.5, ONE, False),
|
||||
('hyperbolic kernel ell=Z/2', tail_hyp, 0.2, lambda m, D, Z, c, cp: 0.5 * Z * np.ones_like(m), False),
|
||||
]
|
||||
pick = [int(x) for x in sys.argv[3].split(',')] if len(sys.argv) > 3 else range(len(LAWS))
|
||||
|
||||
def zone_struct(z):
|
||||
idx = {int(r['idx']): k for k, r in enumerate(z)}
|
||||
n = len(z)
|
||||
par = np.array([idx.get(int(r['parentIdx']), -1) for r in z])
|
||||
D = np.array([float(r['distance']) for r in z]); m = np.array([float(r['mass']) for r in z])
|
||||
body = np.array([r['empty'] == '0' and r['belt'] == '0' for r in z])
|
||||
cls = ['Belt' if r['belt'] == '1' else r['type'] for r in z]
|
||||
fam = np.eye(n, dtype=bool) # fam[i, j]: j is i or a descendant of i
|
||||
depth = np.zeros(n, int)
|
||||
for j in range(n):
|
||||
a = par[j]
|
||||
while a >= 0:
|
||||
fam[a, j] = True; depth[j] += 1; a = par[a]
|
||||
order = list(np.argsort(depth, kind='stable'))
|
||||
return par, D, m, body, cls, fam, order
|
||||
|
||||
def place(par, D, order):
|
||||
pos = np.zeros((len(par), 2)); th = rng.uniform(0, 2 * np.pi, len(par))
|
||||
for k in order:
|
||||
if par[k] >= 0: pos[k] = pos[par[k]] + D[k] * np.array([np.cos(th[k]), np.sin(th[k])])
|
||||
return pos
|
||||
|
||||
def field(pos_b, c, R, s, tail, x):
|
||||
"""Per-body acceleration contributions at points x: (n, nb, 2), each toward its body."""
|
||||
d = pos_b[None, :, :] - x[:, None, :]; r = np.linalg.norm(d, axis=2)
|
||||
rs = np.maximum(r, 1e-6)
|
||||
F = np.where(r < 1e-6, 0, pull(tail, c[None, :], R[None, :], rs, s[None, :]))
|
||||
return d / rs[:, :, None] * F[:, :, None]
|
||||
|
||||
def disc_acc(x, C, Z): # field of a uniform disc of total strength C over radius Z centred at 0 (inward)
|
||||
r = np.linalg.norm(x, axis=1, keepdims=True); rs = np.maximum(r, 1e-9)
|
||||
return -x / rs * np.where(r < Z, C * r / Z ** 2, C / rs)
|
||||
|
||||
def q(v, ps=(10, 50, 90)):
|
||||
v = np.sort(np.asarray(v))
|
||||
return [round(float(v[min(len(v) - 1, int(p / 100 * len(v)))]), 2) for p in ps] if len(v) else []
|
||||
|
||||
pre = {k: zone_struct(z) for k, z in zones.items()}
|
||||
sunm = np.sort([float(r['mass']) for r in rows if r['type'] == 'Sun' and r['empty'] == '0'])
|
||||
SUN_LO, SUN_HI = sunm[int(.4 * len(sunm))], sunm[int(.6 * len(sunm))] # 'median sun': mass between p40 and p60
|
||||
for li in pick:
|
||||
name, tail, p, scale, disc = LAWS[li]
|
||||
t0 = time.time()
|
||||
ratio = collections.defaultdict(list); vpar = collections.defaultdict(list); vsp = []; vmed = []
|
||||
inward = collections.defaultdict(list)
|
||||
for key, z in zones.items():
|
||||
par, D, m, body, cls, fam, order = pre[key]
|
||||
Z = float(z[0]['zoneRadius'])
|
||||
c = np.where(body, strength(m, p), 0.0); R = core_radius(m)
|
||||
cpar = np.array([c[par[k]] if par[k] >= 0 else 0.0 for k in range(len(par))])
|
||||
Dz = np.where(par < 0, Z, D)
|
||||
s = scale(m, Dz, Z, c, cpar)
|
||||
bi = np.where(body)[0]; bpos = {j: t for t, j in enumerate(bi)}
|
||||
kids = [k for k in range(len(par)) if par[k] >= 0 and body[par[k]]]
|
||||
a_own = {k: float(pull(tail, c[par[k]], R[par[k]], D[k], s[par[k]])) for k in kids}
|
||||
for k in kids:
|
||||
v = np.sqrt(a_own[k] * D[k]); vpar[cls[k]].append(v)
|
||||
if cls[k] == 'Planet' and z[par[k]]['type'] == 'Sun':
|
||||
vsp.append(v)
|
||||
if SUN_LO <= m[par[k]] <= SUN_HI: vmed.append(v)
|
||||
acc_rel = collections.defaultdict(float); mids = []
|
||||
fb = fam[:, bi][:, :, None]
|
||||
for dr in range(DRAWS):
|
||||
pos = place(par, D, order)
|
||||
A = field(pos[bi], c[bi], R[bi], s[bi], tail, pos)
|
||||
Aout = np.where(fb, 0, A).sum(axis=1) # at node k, from bodies outside k's family
|
||||
if disc: Aout = Aout - disc_acc(pos, c.sum(), Z)
|
||||
for k in kids:
|
||||
pk = par[k]; u = pos[pk] - pos[k]; u /= np.linalg.norm(u)
|
||||
acc_rel[k] += np.dot(Aout[k] - Aout[pk], u) / DRAWS # Aout[k] includes the parent's own pull
|
||||
th = np.arange(8) * np.pi / 4; X = 0.5 * Z * np.stack([np.cos(th), np.sin(th)], 1)
|
||||
AX = field(pos[bi], c[bi], R[bi], s[bi], tail, X).sum(axis=1)
|
||||
if disc: AX = AX - disc_acc(X, c.sum(), Z)
|
||||
mids.append(np.mean(-(AX * X).sum(1) / np.linalg.norm(X, axis=1)))
|
||||
inward[Z].append(np.mean(mids))
|
||||
for k in kids:
|
||||
if a_own[k] > 0: ratio[cls[k]].append(np.sqrt(max(acc_rel[k], 0) / a_own[k]))
|
||||
vsp = np.array(vsp)
|
||||
print(f'== [{li}] {name} ({time.time() - t0:.0f}s, {DRAWS} draws)')
|
||||
for k in ['Sun', 'GasGiant', 'Planet', 'Planetoid', 'Belt']:
|
||||
print(f' {k:9s} ratio p10/50/90/99 {q(ratio[k], (10, 50, 90, 99))} parent-only v p10/50/90 {q(vpar[k])} n {len(ratio[k])}')
|
||||
print(f' planets round suns: parent-only v p10/50/90 {q(vsp)}; share in 12-18 u/s {np.mean((vsp >= 12) & (vsp <= 18)):.2f}')
|
||||
vmed = np.array(vmed)
|
||||
print(f' planets round a median sun (mass p40-p60, {SUN_LO:.0f}-{SUN_HI:.0f}): v p10/50/90 {q(vmed)}; share in 12-18 {np.mean((vmed >= 12) & (vmed <= 18)):.2f}; n {len(vmed)}')
|
||||
nz = {k: sum(1 for x in vpar[k] if x == 0) for k in vpar}
|
||||
print(' children with no parent pull at their distance: ' + str({k: v for k, v in nz.items() if v}))
|
||||
print(' mid-zone inward pull by Z: ' + ', '.join(f'{Z:.0f}: {np.median(v):.2f}' for Z, v in sorted(inward.items())))
|
||||
sys.stdout.flush()
|
||||
Loading…
Reference in a new issue