Research: the conserving compact well law (band terms on the orbit tree), for ruling tame-tails-parent-only-rails
Each node pulls exactly as its body inside its outermost child, exactly as one body of its family's mass beyond its own support, and nothing beyond its reach; the root has no reach, so the zone's far field is exactly c(M_zone) however the generator splits it. C4 on every band edge, closed-form potential (grid = K phi), Plummer core. Probe against today, a2, c2, b1 and the tapered draft: ratio, funnel, 12-18 u/s, far field, split invariance, whole-zone depth, and evaluation cost. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Meta Vens
parent
8bc5ef1837
commit
b2f6122dd2
8 changed files with 829 additions and 0 deletions
103
docs/research/gravity-well-law.md
Normal file
103
docs/research/gravity-well-law.md
Normal file
|
|
@ -0,0 +1,103 @@
|
|||
# The conserving compact well law: derivation and probe
|
||||
|
||||
Imagination pass `imagination-ag3-well-law`, campaign aetheria-gravity, session self-2026-10-15-ag3, 2026-10-15, read at master 8bc5ef18. Serves ruling `aetheria-gravity:ruling:tame-tails-parent-only-rails` and revises question `aetheria-gravity:question:tail-taming-law` (now `...:tail-taming-law-r2`). Builds on `gravity-tail-taming.md` and `gravity-smooth-wells.md`: same data, same measures, their scripts reused.
|
||||
|
||||
The operator, verbatim, in order:
|
||||
- "Having all the derivatives shake out nicely isn't just an affectation, the player can feel higher-order discontinuities especially in a game where moving with the grid is such a huge part of the experience."
|
||||
- "this is all a consequence of exaggerating the contributions of small bodies. We're doing that so everything can coexist within one order of magnitude of scale, but that compression has annoying side-effects, like a well getting deeper the more its bodies are split by the zone generator. Two systems, same mass, one of them spent it on planets, that one goes waaaay deeper."
|
||||
- On mass conserved at range: "That would be quite convenient if we can make it work out that way."
|
||||
- "would be even more convenient if we can use this as an optimization and not have to evaluate the little guys at all except when near them"
|
||||
|
||||
Scripts in `gravity-well-law/`:
|
||||
- `well_law.py`: the band term (closed-form pull and potential), and the tapered first draft.
|
||||
- `check_well_law.py` and `check-output.txt`: force = dφ/dr, one-sided derivative jumps at 50 digits, circular speed, float32 error.
|
||||
- `selfpull.py` and `selfpull-output.txt`: the tapered draft's self-pull against the band form.
|
||||
- `probe3.py` and `probe3-output.txt`: the term tree per zone, the measures of `probe2.py`, and split invariance, whole-zone depth, containment and cost. Each law ran in its own process (rng seed 7, 6 draws) on Starfire, 20-95 s.
|
||||
|
||||
## 1. The law
|
||||
|
||||
**One primitive, the band term** U(r; ε, a, b): a Plummer-softened log source whose pull is exact up to radius a and fades to zero over the band [a, b].
|
||||
- Pull toward the centre: F(r) = r/(r²+ε²) · B(s), with s = r².
|
||||
- B(s) = 1 for s ≤ a²; B = 1 − S(t) on the band, t = (s − a²)/(b² − a²); B = 0 for s ≥ b².
|
||||
- S is the order-4 smoothstep, S(t) = 35t⁴ − 84t⁵ + 70t⁶ − 20t⁷.
|
||||
- Potential: 0 beyond b. On the band, ½[R·ln((t+τ)/(1+τ)) − (I(1) − I(t))], where τ = (a²+ε²)/(b²−a²), (1 − S(t)) = (t+τ)·Q(t) + R, and I is the antiderivative of Q. Inside a, ½ln((s+ε²)/(a²+ε²)) plus the band's value at t = 0.
|
||||
- a = b = ∞ is the global term, ½ln(r²+ε²).
|
||||
- The cost is one log, plus a 7-step synthetic division and Horner on the band only.
|
||||
|
||||
**A node's terms.** Every node F of the orbit tree contributes at its own position, F being a body, or a bodyless node with its ring group standing in for the body:
|
||||
|
||||
c(m_F)·U(ε(m_F), L_F, B_F) + c(M_F)·[U(ε(M_F), max(B_F, ½Λ_F), Λ_F) − U(ε(M_F), L_F, B_F)]
|
||||
|
||||
- m_F is the body's mass and M_F the family's. c and ε are the existing `WellStrength` and `RocheLimit` curves, evaluated at m_F or M_F; no new curve.
|
||||
- L_F is the outermost child's distance.
|
||||
- B_F = max(1.5·L_F, max over children of D_C + Λ_C): the family's own support.
|
||||
- Λ_F is the family's reach in its parent: max(nominal, 1.5·B_F). The nominal reach probed is 0.5·D, 1·D or 2·d_Hill, with d_Hill = D·sqrt(c(M_F)/(2c_parent)) (B25's tidal radius, for the family's mass).
|
||||
- A ring member of a bodyless node reaches 4·D, so its partners across the ring feel it whole.
|
||||
- A root's Λ is infinite.
|
||||
- A body with no massive descendants is one term, c(m)·U(ε(m), ½Λ, Λ).
|
||||
|
||||
**Bodyless nodes.** The ring group G is the massive children at the innermost distance. It contributes c(M_G)·[U(ε_G, L_F, B_F) − U(ε_G, D_ring, B_G)], with B_G = max(1.5·D_ring, max over members of D + Λ). L_F is raised to at least B_G.
|
||||
|
||||
**What this guarantees, by construction (not approximated):**
|
||||
- *Inside L_F the node pulls exactly as its body alone*, c(m_F)·d/(d²+ε²). The family terms are zero there, since both of their windows are flat at 1 inside L_F and cancel. So a child's parent-only rail is exactly the body's Plummer pull, and no child feels its own mass through its parent's family term. The tapered draft did feel it: a lone body on a bodyless node was pulled at 38% of its own pull by its own mass (`selfpull.py`; the band form: 0).
|
||||
- *Beyond B_F the family pulls exactly as one body of mass M_F at the node*, c(M_F)·r/(r²+ε(M_F)²), out to max(B_F, ½Λ_F), and fades to nothing at Λ_F. Every descendant term lies inside B_F: the probe counts 0 containment violations of 70,890 (term, ancestor) pairs.
|
||||
- *Split invariance.* Merge a family into one body of its mass in the same slot. Outside B_F the field is identical, term for term, whenever the merge leaves Λ_F unchanged, that is, whenever the nominal reach is at least 1.5·B_F. With the tidal reach this holds for 2382 of the 2427 families whose children sit inside the ruled generator window (half the tidal radius, `moons-within-reach-ruled`), and for 4573 of all 6976 families in today's generator, which places most children outside that window (B25).
|
||||
- *The zone's far field is exactly c(M_zone)*, at 2Z and at 10Z (p10 = p90 = 1.00). The root's Λ is infinite, so beyond B_root the zone is one body of its total mass however the generator split it.
|
||||
- *Smoothness.* φ is a function of r², so it is C^∞ at every centre and wherever B is constant. On each band edge it is C⁴: `check-output.txt` measures φ⁽⁰⁾-φ⁽⁴⁾ jumps of at most 1e-45 (zero at 50 digits) and a φ⁽⁵⁾ jump of 5.5e3 at r = a, the first nonzero order. The soft core is Plummer: no C1 seam at the Roche limit. The order is one constant: order 6 gives C⁶ for 4 more FMAs.
|
||||
- *Grid = K·φ.* The closed form satisfies F = dφ/dr to 1.6e-9 (band) and 1.7e-6 (global term, finite-difference floor). In float32 the potential's error is at most 1.5e-5 per unit strength, about 2e-3 u of drawn height at c = 300 and K = 2.05, against the parity verify's 0.03 u tolerance.
|
||||
- *Time.* Every term rides its node's rail position, and its parameters are fixed when the zone is built, so the field is C^∞ in time. The centre is the node's body, not the family's centre of mass: what is exact is "the merged body in the same slot". A centre of mass would add a moving centroid per family and buy nothing a test can see beyond B_F.
|
||||
|
||||
**The peak slope.** Plummer's peak pull is 0.5·c/ε, at r = ε, which is the Roche limit. Today's harmonic core peaks at 1.0·c/R. `core-at-roche-and-exaggeration-up` holds the peak drawn slope, so K must halve, from 4.0983607 to 2.0491803. Every slope outside the cores then draws 2x steeper than at today's K, while the tamed funnel pulls about 2.6x less.
|
||||
|
||||
## 2. Probe
|
||||
|
||||
The ratio is the free-coasting circular speed over the parent-only rail speed (1 = rails honest). The other measures are as in `gravity-smooth-wells.md` §2. "Split" is the share of in-window families whose field outside B_F is unchanged by merging them, then the error p50/p90 at 1.05·B_F over all 6976 families. "Depth" is the zone's deepest point against one body of the zone's mass, p50/p90.
|
||||
|
||||
| # | Law | Ratio p50/p90: Sun · GasGiant · Planet · Planetoid · Belt | Mid-zone pull Z = 1000 / 2125 / 4182 / 6846 / 10000 | Median-sun planets v p50; share 12-18 | Sun-class rails v p50 | Far field / c(M_zone), 2Z; 10Z | Split | Depth | Pull < 0.5 c/D |
|
||||
|---|---|---|---|---|---|---|---|---|---|
|
||||
| 0 | today | 2.00/2.47 · 1.63/2.23 · 1.24/1.72 · 1.10/1.43 · 1.16/1.55 | 1.59/1.28/0.92/0.68/0.51 | 17.48; 0.92 | 23.62 | 5.54; 5.54 | n/a; 0.53/1.50 | 3.43/4.83 | 46 |
|
||||
| a2 | log-ratio L = 1·D | 1.26/1.55 · 1.12/1.45 · 1.02/1.27 · 1.01/1.20 · 1.02/1.26 | 0.65/0.46/0.29/0.21/0.15 | 16.86; 1.00 | 19.93 | 0.36; 0.02 | n/a; 0.48/1.18 | 1.28/1.45 | 407 |
|
||||
| c2 | Lamb-Oseen + halo L = 1·D | 1.26/1.55 · 1.11/1.45 · 1.02/1.27 · 1.00/1.19 · 1.02/1.26 | 0.66/0.46/0.30/0.21/0.15 | 17.16; 0.99 | 20.04 | 0.36; 0.02 | n/a; 0.48/1.19 | 1.26/1.41 | 146 |
|
||||
| b1 | screened Plummer λ = 0.5·D (smooth k0-half-d) | 1.16/1.49 · 1.04/1.37 · 0.99/1.20 · 0.99/1.16 · 0.99/1.22 | 0.42/0.29/0.19/0.13/0.09 | 16.30; 0.98 | 17.20 | 0.07; 0.00 | n/a; 0.40/0.95 | 1.10/1.18 | 3610 |
|
||||
| t3 | tapered draft, Λ = 2·d_Hill | 1.05/1.27 · 0.99/1.24 · 0.97/1.14 · 0.98/1.13 · 0.98/1.17 | 0.56/0.40/0.25/0.17/0.13 | 17.16; 0.99 | 23.60 | 1.00; 1.00 | 2427/2427; 0.00/0.10 | 0.99/1.31 | 234 |
|
||||
| w1 | band, Λ = 0.5·D | 1.01/1.24 · 1.00/1.24 · 0.99/1.14 · 0.99/1.12 · 0.99/1.15 | 0.55/0.39/0.24/0.17/0.13 | 17.17; 0.99 | 23.60 | 1.00; 1.00 | 2075/2427; 0.52/1.63 | 1.00/1.32 | 234 |
|
||||
| w2 | band, Λ = 1·D | 1.09/1.32 · 1.01/1.30 · 0.97/1.18 · 0.97/1.15 · 0.98/1.19 | 0.61/0.43/0.27/0.19/0.14 | 17.17; 0.99 | 23.60 | 1.00; 1.00 | 2380/2427; 0.00/0.89 | 1.01/1.36 | 234 |
|
||||
| w3 | **band, Λ = 2·d_Hill** | 1.09/1.32 · 1.00/1.29 · 0.97/1.17 · 0.97/1.15 · 0.98/1.19 | 0.60/0.42/0.27/0.19/0.14 | 17.17; 0.99 | 23.60 | 1.00; 1.00 | 2382/2427; 0.00/0.91 | 1.01/1.34 | 234 |
|
||||
|
||||
Readings:
|
||||
- **The depth column is the operator's complaint, measured.** Today a zone's deepest point is 3.43x (p90 4.83x) the depth of one body of the same mass. Under the conserving law it is 1.01x (p90 1.34x). The far field at range is 5.54x today and exactly 1.00x here.
|
||||
- **Rails.** The p50 ratio is 0.97-1.09 for every class (today 1.10-2.00). The p10-p90 spread is 0.77-1.32, from siblings felt inside their reach. The band form is about 0.05 wider at p90 than the tapered draft at the same reach: a sibling is felt whole out to half its reach. Under b1 (screened, not conserving) the suns' children's rails slow from 23.6 to 17.2 u/s and 3610 children fall below half the log pull. Under the conserving law the root has no reach, so its children keep their rails (23.6), and only the 234 children inside their parent's Roche limit (B35's 234) are weak.
|
||||
- **12-18 u/s** holds: planets round a median sun run at 17.2 u/s p50, 99% in band.
|
||||
- **Funnel.** Mid-zone pull falls to 0.38-0.40 of today's in every zone size.
|
||||
- **The 234.** These are the children inside their parent's Roche limit, unchanged from B35. The Plummer core gives them 0.71·√c at ε rather than the harmonic core's solid-body rail. `children-within-reach` moves them out.
|
||||
- **Split.** The band with the tidal reach is exact for 98% of in-window families. The tapered draft is exact for 100%, but at the cost of the self-pull above. Families outside the window (today's generator) differ near B_F (p90 0.91) and never at 2Z.
|
||||
|
||||
## 3. Evaluation cost
|
||||
|
||||
- **Today:** every well everywhere, 17 per sample or texel (p99 58), each one log.
|
||||
- **Conserving law (w3):** 30 terms per zone (p99 114).
|
||||
- *Sim sample* (uniform in the zone disc), walking the tree and skipping any subtree whose reach does not contain the sample: 4 terms evaluated p50 (p99 15, mean 5.6), after 10 reach tests (p99 34).
|
||||
- *Shader texel*, in the 8192 u Surface Displacement window (2048² at 4 u), flat loop with a reach test per term: 3 terms in support p50 (p99 14; per-window mean p50 3.7, p99 12.6). 21 terms meet the window p50 (p99 76), so a CPU cull to the window is optional.
|
||||
- The walk is exact, not Barnes-Hut: the field outside a subtree's reach is zero, so skipping it changes nothing.
|
||||
|
||||
## 4. Consequences the cut carries, or names for later
|
||||
|
||||
- **Rails (orbits-from-field's ring law).** A child's period reads its parent node's own terms at d (`PullAround`). A ring member of a bodyless node adds each ring partner's own terms, projected toward the node at their chord. A binary member at chord 2D therefore feels c(M_partner)·2D/(4D²+ε²), today's half-partner rule in the new law. Gauss's "inside" sum and the recursive `Strength` go. A strengthless orbit on a bodyless ring keeps taking the mean of what its bodied partners feel, until `children-within-reach` deletes that case (`binaries-orbited-from-afar`, follow-up `companion-ring-residuals-r2`).
|
||||
- **K halves** (§1, "The peak slope").
|
||||
- **children-within-reach (named, not done here).** The law owns a reach (Λ) and a family support (B). The generator's "within reach" should read them rather than spell the tidal radius a second time: children inside half d_Hill fit Λ = 2·d_Hill with room (98% exact split). Circumbinary children belong beyond the ring group's B_G = 5·D_ring, which is 2.5x the separation, inside Holman & Wiegert's 2-4x.
|
||||
- **The outer body round a binary** feels c(M_pair+moons), the compressed strength of the summed mass, not today's sum of strengths. That is conservation, and `OrbitsFromFieldTests.AnOuterBodyFeelsTheWholePairAndItsMoons` changes to match.
|
||||
- **gravity-law-one-spelling** stays open: the band arithmetic is still spelled in C# and in the HLSL template. Its coefficients can be generated from one C# array.
|
||||
|
||||
## Prior art
|
||||
|
||||
- **The construction is a partition of unity with charge neutrality.** Ewald 1921 splits each charge into a neutralised short-range part plus a long-range remainder. TreePM codes (Bagla 2002; Springel 2005, GADGET-2) use an erfc split with the long range on a mesh. The fast multipole method (Greengard & Rokhlin 1987) gives exact far fields of clusters by multipole expansion, with a near/far split by well-separatedness. Here the far part of each family is truncated at its monopole by fiat, which is exactly the operator's "mass conserved at range". The near part is exact and compact, so no expansion error arises: the law is defined by the truncation rather than approximated by it.
|
||||
- **Compact smooth windows.** Wendland 1995 (Adv. Comput. Math. 4:389) gives compactly supported C^{2k} radial functions. The smoothstep-in-r² used here is a C^k window of the same family, chosen because ∫ r·B(r²)/(r²+ε²) dr stays elementary (log plus polynomial), so the grid stays exactly K·φ.
|
||||
- **Force switching.** Steinbach & Brooks 1994 prefer switching the force over switching the potential. This law switches the source profile, so the force and the potential agree exactly, and it switches with C⁴ continuity rather than C¹.
|
||||
- **Spheres of influence.** Patched conics (KSP, Laplace's r_SOI = a·(m/M)^(2/5)) let one attractor own each sphere of influence with a kink at the boundary. This law keeps that dominance (inside L_F only the parent pulls on the rail) and replaces the kink with a C⁴ band.
|
||||
|
||||
## Confidence
|
||||
|
||||
- §1: high. The guarantees are construction plus the measured checks: zero containment violations, far field exactly 1.00, C⁴ jumps at 50 digits.
|
||||
- §2: high on internal consistency. Row 0 reproduces the earlier baseline within 0.01 (2.00 against 2.01, from a different rng consumption order). Medium on absolute ratios, for the caveats of `gravity-tail-taming.md`: 6 phase draws, phase-averaged circular speeds, positions from the 2026-10-13 generator.
|
||||
- §3: high on counts. Op costs are by inspection, not compiled HLSL.
|
||||
- The knobs (order 4, band 1.5, fade from ½Λ, ring reach 4) were chosen once, not swept. Each is a single constant.
|
||||
17
docs/research/gravity-well-law/check-output.txt
Normal file
17
docs/research/gravity-well-law/check-output.txt
Normal file
|
|
@ -0,0 +1,17 @@
|
|||
S(t) coefficients (lowest first): [0.0, 0.0, 0.0, 0.0, 35.0, -84.0, 70.0, -20.0]
|
||||
1 - S(t): [1.0, -0.0, -0.0, -0.0, -35.0, 84.0, -70.0, 20.0]
|
||||
eps 79.8 band [ 1000, 1500]: max |F - dphi/dr| / max F = 1.6e-09; phi(b) = 0.0e+00
|
||||
eps 79.8 band [ 0, 1525]: max |F - dphi/dr| / max F = 1.6e-09; phi(b) = 0.0e+00
|
||||
eps 17.0 band [ 40, 60]: max |F - dphi/dr| / max F = 3.7e-09; phi(b) = 0.0e+00
|
||||
eps 32.7 band [ 400, 800]: max |F - dphi/dr| / max F = 2.6e-09; phi(b) = 0.0e+00
|
||||
eps 79.8 band [ 12000, 25000]: max |F - dphi/dr| / max F = 4.3e-07; phi(b) = 0.0e+00
|
||||
global term: max |F - dphi/dr| / max F = 1.7e-06
|
||||
edge r = 1000: jumps of phi^(0..6) x r^n = ['0.0', '6.264e-50', '-4.078e-50', '1.792e-49', '-8.654e-49', '-5470.0', '7.733e+4']
|
||||
edge r = 1500: jumps of phi^(0..6) x r^n = ['0.0', '2.518e-49', '5.419e-48', '7.407e-47', '8.42e-46', '-1.407e+5', '-6.785e+6']
|
||||
odd derivatives at r = 0 (inside branch, must be 0): ['0.0', '0.0', '5.66e-73']
|
||||
median sun: c = 304.9, sqrt(c) = 17.46 u/s, eps = Roche = 79.8
|
||||
v/sqrt(c) at eps, 1.5 eps, 2 eps, 3 eps, 610 u: [0.707, 0.832, 0.894, 0.949, 0.992]
|
||||
peak pull: 0.500 c/eps at r = 1.000 eps (today: 1.0 c/R at R)
|
||||
float32 eps 79.8 band [1000.0, 1500.0]: max abs err 7.28e-06 (unit strength; depth 2.77)
|
||||
float32 eps 17.0 band [40.0, 60.0]: max abs err 1.45e-05 (unit strength; depth 1.15)
|
||||
float32 eps 79.8 band [12000.0, 25000.0]: max abs err 9.06e-07 (unit strength; depth 5.49)
|
||||
69
docs/research/gravity-well-law/check_well_law.py
Normal file
69
docs/research/gravity-well-law/check_well_law.py
Normal file
|
|
@ -0,0 +1,69 @@
|
|||
# Checks of the band term (well_law.py): force = dphi/dr (grid equals field), smoothness class on both band edges
|
||||
# (one-sided derivatives at 50 digits), circular speed round a median sun, and float32 error of the closed form.
|
||||
import os, sys, math
|
||||
import numpy as np, mpmath as mp
|
||||
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
|
||||
sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', 'gravity-tail-taming'))
|
||||
from well_law import *
|
||||
from laws import strength, core_radius
|
||||
|
||||
print('S(t) coefficients (lowest first):', SC.tolist())
|
||||
print('1 - S(t):', FADE.tolist())
|
||||
|
||||
# 1. force = dphi/dr, band terms, the tapered draft (a = 0) and the global term
|
||||
for eps, a, b in [(79.8, 1000.0, 1500.0), (79.8, 0.0, 1525.0), (17.0, 40.0, 60.0), (32.7, 400.0, 800.0), (79.8, 12000.0, 25000.0)]:
|
||||
r = np.linspace(1e-3, 1.2 * b, 400001)
|
||||
p = band_potential(eps, a, b, r); f = band_pull(eps, a, b, r); d = np.gradient(p, r)
|
||||
print(f'eps {eps:6.1f} band [{a:7.0f}, {b:7.0f}]: max |F - dphi/dr| / max F = {np.max(np.abs(d - f)[3:-3]) / f.max():.1e}; '
|
||||
f'phi(b) = {float(band_potential(eps, a, b, np.array(b))):.1e}')
|
||||
r = np.linspace(1e-3, 30000, 200001)
|
||||
p = band_potential(79.8, np.inf, np.inf, r); f = band_pull(79.8, np.inf, np.inf, r)
|
||||
print(f'global term: max |F - dphi/dr| / max F = {np.max(np.abs(np.gradient(p, r) - f)[3:-3]) / f.max():.1e}')
|
||||
|
||||
# 2. smoothness at both band edges: derivatives of each branch at the edge, at 50 digits
|
||||
mp.mp.dps = 50
|
||||
def branches(eps, a, b):
|
||||
eps, a, b = mp.mpf(eps), mp.mpf(a), mp.mpf(b)
|
||||
span = b * b - a * a; tau = (a * a + eps * eps) / span
|
||||
n = len(FADE) - 1; q = [None] * n; q[n - 1] = mp.mpf(FADE[n])
|
||||
for i in range(n - 1, 0, -1): q[i - 1] = mp.mpf(FADE[i]) - tau * q[i]
|
||||
R = mp.mpf(FADE[0]) - tau * q[0]
|
||||
I = lambda t: sum(q[i] * t ** (i + 1) / (i + 1) for i in range(n))
|
||||
band = lambda r: (R * mp.log(((r * r - a * a) / span + tau) / (1 + tau)) - (I(1) - I((r * r - a * a) / span))) / 2
|
||||
inside = lambda r: mp.log((r * r + eps * eps) / (a * a + eps * eps)) / 2 + band(a)
|
||||
return inside, band, (lambda r: mp.mpf(0))
|
||||
eps, a, b = 79.8, 1000.0, 1500.0
|
||||
inside, band, beyond = branches(eps, a, b)
|
||||
for edge, lo, hi in ((a, inside, band), (b, band, beyond)):
|
||||
jumps = [mp.nstr((mp.diff(hi, mp.mpf(edge), n) - mp.diff(lo, mp.mpf(edge), n)) * mp.mpf(edge) ** n, 4) for n in range(7)]
|
||||
print(f'edge r = {edge:.0f}: jumps of phi^(0..6) x r^n = {jumps}')
|
||||
print('odd derivatives at r = 0 (inside branch, must be 0):', [mp.nstr(mp.diff(inside, mp.mpf(0), n), 3) for n in (1, 3, 5)])
|
||||
|
||||
# 3. circular speed round a median sun: the body term is exact Plummer inside a, so v/sqrt(c) = r / sqrt(r^2 + eps^2)
|
||||
m = 14152.4; c = strength(m); e = core_radius(m)
|
||||
print(f'median sun: c = {c:.1f}, sqrt(c) = {math.sqrt(c):.2f} u/s, eps = Roche = {e:.1f}')
|
||||
rr = np.array([e, 1.5 * e, 2 * e, 3 * e, 610.0])
|
||||
print(' v/sqrt(c) at eps, 1.5 eps, 2 eps, 3 eps, 610 u:', np.round(np.sqrt(rr * c * band_pull(e, np.inf, np.inf, rr)) / math.sqrt(c), 3).tolist())
|
||||
rr = np.logspace(0, 4, 40001); F = band_pull(e, np.inf, np.inf, rr)
|
||||
print(f' peak pull: {F.max() * e:.3f} c/eps at r = {rr[np.argmax(F)] / e:.3f} eps (today: 1.0 c/R at R)')
|
||||
|
||||
# 4. float32: the band closed form evaluated in single precision against float64
|
||||
def band32(eps, a, b, r):
|
||||
f32 = np.float32
|
||||
eps, a, b, r = f32(eps), f32(a), f32(b), r.astype(f32)
|
||||
s = r * r; span = b * b - a * a; tau = (a * a + eps * eps) / span
|
||||
n = len(FADE) - 1; q = [None] * n; q[n - 1] = f32(FADE[n])
|
||||
for i in range(n - 1, 0, -1): q[i - 1] = f32(FADE[i]) - tau * q[i]
|
||||
R = f32(FADE[0]) - tau * q[0]
|
||||
def I(t):
|
||||
acc = np.zeros_like(t)
|
||||
for i in range(n - 1, -1, -1): acc = ((acc + q[i] / f32(i + 1)) * t).astype(f32)
|
||||
return acc
|
||||
t = np.clip((s - a * a) / span, 0, 1).astype(f32)
|
||||
bandv = f32(.5) * (R * np.log((t + tau) / (f32(1) + tau)).astype(f32) - (I(np.ones_like(t)) - I(t)))
|
||||
inside = f32(.5) * np.log((s + eps * eps) / (a * a + eps * eps)).astype(f32) + bandv[0] * 0 + f32(.5) * (R * np.log(tau / (f32(1) + tau)) - I(np.ones(1, f32))[0])
|
||||
return np.where(s <= a * a, inside, np.where(s >= b * b, f32(0), bandv))
|
||||
for eps, a, b in [(79.8, 1000.0, 1500.0), (17.0, 40.0, 60.0), (79.8, 12000.0, 25000.0)]:
|
||||
r64 = np.linspace(0, b, 100001)
|
||||
ref = band_potential(eps, a, b, r64)
|
||||
print(f'float32 eps {eps} band [{a}, {b}]: max abs err {np.max(np.abs(band32(eps, a, b, r64) - ref)):.2e} (unit strength; depth {np.abs(ref).max():.2f})')
|
||||
165
docs/research/gravity-well-law/probe3-output.txt
Normal file
165
docs/research/gravity-well-law/probe3-output.txt
Normal file
|
|
@ -0,0 +1,165 @@
|
|||
# python probe3.py bodies.csv 6 <i>, one process per law (rng seed 7 each), Starfire 2026-10-15
|
||||
== [0] today: harmonic core + log (21s, 6 draws)
|
||||
Sun ratio p10/50/90/99 [1.44, 2.0, 2.47, 2.87] parent-only v p10/50/90 [22.06, 23.62, 24.64] n 2211
|
||||
GasGiant ratio p10/50/90/99 [1.1, 1.63, 2.23, 2.71] 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.72, 2.18] parent-only v p10/50/90 [15.56, 17.38, 23.61] n 7445
|
||||
Planetoid ratio p10/50/90/99 [0.95, 1.1, 1.43, 1.8] parent-only v p10/50/90 [11.51, 14.48, 20.27] n 10200
|
||||
Belt ratio p10/50/90/99 [0.96, 1.16, 1.55, 1.98] parent-only v p10/50/90 [11.41, 14.74, 20.29] n 3134
|
||||
planets round a median sun: v p10/50/90 [17.06, 17.48, 17.98]; share in 12-18 0.92; n 1116
|
||||
Sun-class children parent-only v p10/50/90 [22.06, 23.62, 24.64]
|
||||
children with no parent pull: {}; parent pull < 0.5 c/D: {'Planetoid': 42, 'Belt': 3, 'Planet': 1} (total 46)
|
||||
mid-zone inward pull by Z: 1000: 1.59, 2125: 1.28, 4182: 0.92, 6846: 0.68, 10000: 0.51
|
||||
zone far-field c / c(M_zone) p10/50/90: at 2Z [3.69, 5.54, 7.98], at 10Z [3.69, 5.54, 7.98]
|
||||
split invariance (family merged into one body), max |da|/|a| p50/p90/p99: at 1.05 B_f [0.53, 1.5, 5.47], at 2Z [0.07, 0.2, 0.41]; n 6976
|
||||
whole-zone depth over one body of the zone mass, p10/50/90: [2.48, 3.43, 4.83]
|
||||
wells (bodies) per zone p50/p99: [17.0, 58.0]
|
||||
== [1] (a2) log-ratio L=1*D (20s, 6 draws)
|
||||
Sun ratio p10/50/90/99 [0.99, 1.26, 1.55, 1.89] parent-only v p10/50/90 [17.73, 19.93, 22.72] n 2211
|
||||
GasGiant ratio p10/50/90/99 [0.91, 1.12, 1.45, 1.84] parent-only v p10/50/90 [16.19, 19.44, 23.25] n 3566
|
||||
Planet ratio p10/50/90/99 [0.91, 1.02, 1.27, 1.63] parent-only v p10/50/90 [13.89, 16.58, 22.83] n 7445
|
||||
Planetoid ratio p10/50/90/99 [0.9, 1.01, 1.2, 1.65] parent-only v p10/50/90 [9.99, 13.32, 18.04] n 10200
|
||||
Belt ratio p10/50/90/99 [0.9, 1.02, 1.26, 1.65] parent-only v p10/50/90 [9.83, 13.7, 18.99] n 3134
|
||||
planets round a median sun: v p10/50/90 [15.52, 16.86, 17.49]; share in 12-18 1.00; n 1116
|
||||
Sun-class children parent-only v p10/50/90 [17.73, 19.93, 22.72]
|
||||
children with no parent pull: {}; parent pull < 0.5 c/D: {'Planetoid': 329, 'Belt': 47, 'Planet': 31} (total 407)
|
||||
mid-zone inward pull by Z: 1000: 0.65, 2125: 0.46, 4182: 0.29, 6846: 0.21, 10000: 0.15
|
||||
zone far-field c / c(M_zone) p10/50/90: at 2Z [0.25, 0.36, 0.4], at 10Z [0.01, 0.02, 0.02]
|
||||
split invariance (family merged into one body), max |da|/|a| p50/p90/p99: at 1.05 B_f [0.48, 1.18, 4.15], at 2Z [0.0, 0.01, 0.04]; n 6976
|
||||
whole-zone depth over one body of the zone mass, p10/50/90: [1.17, 1.28, 1.45]
|
||||
wells (bodies) per zone p50/p99: [17.0, 58.0]
|
||||
== [2] (c2) Lamb-Oseen + halo L=1*D (26s, 6 draws)
|
||||
Sun ratio p10/50/90/99 [0.99, 1.26, 1.55, 1.9] parent-only v p10/50/90 [17.77, 20.04, 22.76] n 2211
|
||||
GasGiant ratio p10/50/90/99 [0.91, 1.11, 1.45, 1.86] parent-only v p10/50/90 [16.88, 19.7, 23.44] n 3566
|
||||
Planet ratio p10/50/90/99 [0.91, 1.02, 1.27, 1.64] parent-only v p10/50/90 [14.35, 17.0, 23.35] n 7445
|
||||
Planetoid ratio p10/50/90/99 [0.9, 1.0, 1.19, 1.66] parent-only v p10/50/90 [10.37, 13.93, 19.83] n 10200
|
||||
Belt ratio p10/50/90/99 [0.89, 1.02, 1.26, 1.68] parent-only v p10/50/90 [10.18, 14.16, 19.85] n 3134
|
||||
planets round a median sun: v p10/50/90 [16.49, 17.16, 17.75]; share in 12-18 0.99; n 1116
|
||||
Sun-class children parent-only v p10/50/90 [17.77, 20.04, 22.76]
|
||||
children with no parent pull: {}; parent pull < 0.5 c/D: {'Planetoid': 131, 'Belt': 8, 'Planet': 7} (total 146)
|
||||
mid-zone inward pull by Z: 1000: 0.66, 2125: 0.46, 4182: 0.30, 6846: 0.21, 10000: 0.15
|
||||
zone far-field c / c(M_zone) p10/50/90: at 2Z [0.25, 0.36, 0.4], at 10Z [0.01, 0.02, 0.02]
|
||||
split invariance (family merged into one body), max |da|/|a| p50/p90/p99: at 1.05 B_f [0.48, 1.19, 4.16], at 2Z [0.0, 0.01, 0.04]; n 6976
|
||||
whole-zone depth over one body of the zone mass, p10/50/90: [1.15, 1.26, 1.41]
|
||||
wells (bodies) per zone p50/p99: [17.0, 58.0]
|
||||
== [3] (b1) screened Plummer lam=0.5*D (33s, 6 draws)
|
||||
Sun ratio p10/50/90/99 [0.91, 1.16, 1.49, 1.9] parent-only v p10/50/90 [14.3, 17.2, 20.5] n 2211
|
||||
GasGiant ratio p10/50/90/99 [0.85, 1.04, 1.37, 1.81] parent-only v p10/50/90 [14.3, 17.74, 21.89] n 3566
|
||||
Planet ratio p10/50/90/99 [0.89, 0.99, 1.2, 1.62] parent-only v p10/50/90 [12.24, 15.92, 22.05] n 7445
|
||||
Planetoid ratio p10/50/90/99 [0.88, 0.99, 1.16, 1.8] parent-only v p10/50/90 [8.7, 12.75, 17.65] n 10200
|
||||
Belt ratio p10/50/90/99 [0.87, 0.99, 1.22, 1.82] parent-only v p10/50/90 [8.43, 12.9, 18.01] n 3134
|
||||
planets round a median sun: v p10/50/90 [14.71, 16.3, 17.03]; share in 12-18 0.98; n 1116
|
||||
Sun-class children parent-only v p10/50/90 [14.3, 17.2, 20.5]
|
||||
children with no parent pull: {}; parent pull < 0.5 c/D: {'Planetoid': 1326, 'Planet': 530, 'Belt': 420, 'GasGiant': 442, 'Sun': 892} (total 3610)
|
||||
mid-zone inward pull by Z: 1000: 0.42, 2125: 0.29, 4182: 0.19, 6846: 0.13, 10000: 0.09
|
||||
zone far-field c / c(M_zone) p10/50/90: at 2Z [0.02, 0.07, 0.08], at 10Z [0.0, 0.0, 0.0]
|
||||
split invariance (family merged into one body), max |da|/|a| p50/p90/p99: at 1.05 B_f [0.4, 0.95, 3.12], at 2Z [0.0, 0.01, 0.02]; n 6976
|
||||
whole-zone depth over one body of the zone mass, p10/50/90: [1.05, 1.1, 1.18]
|
||||
wells (bodies) per zone p50/p99: [17.0, 58.0]
|
||||
== [4] (t1) tapered draft, Lam = 0.5*D (91s, 6 draws)
|
||||
Sun ratio p10/50/90/99 [0.86, 1.0, 1.19, 1.46] parent-only v p10/50/90 [22.0, 23.65, 24.67] n 2211
|
||||
GasGiant ratio p10/50/90/99 [0.84, 1.0, 1.2, 1.6] parent-only v p10/50/90 [17.58, 22.08, 24.49] n 3566
|
||||
Planet ratio p10/50/90/99 [0.88, 1.0, 1.11, 1.42] parent-only v p10/50/90 [14.91, 16.88, 23.2] n 7445
|
||||
Planetoid ratio p10/50/90/99 [0.9, 1.0, 1.11, 1.42] parent-only v p10/50/90 [10.96, 13.73, 18.14] n 10200
|
||||
Belt ratio p10/50/90/99 [0.88, 0.99, 1.14, 1.43] parent-only v p10/50/90 [10.97, 14.23, 19.45] n 3134
|
||||
planets round a median sun: v p10/50/90 [16.07, 17.13, 17.75]; share in 12-18 0.99; n 1116
|
||||
Sun-class children parent-only v p10/50/90 [22.0, 23.65, 24.67]
|
||||
children with no parent pull: {}; parent pull < 0.5 c/D: {'Planetoid': 214, 'Belt': 11, 'Planet': 9} (total 234)
|
||||
mid-zone inward pull by Z: 1000: 0.55, 2125: 0.38, 4182: 0.23, 6846: 0.17, 10000: 0.12
|
||||
zone far-field c / c(M_zone) p10/50/90: at 2Z [1.0, 1.0, 1.0], at 10Z [1.0, 1.0, 1.0]
|
||||
split invariance (family merged into one body), max |da|/|a| p50/p90/p99: at 1.05 B_f [0.25, 1.13, 3.59], at 2Z [0.0, 0.0, 0.0]; n 6976; families whose reach the merge leaves unchanged 2575/6976, of those inside the ruled half-tidal window 1744/2427
|
||||
whole-zone depth over one body of the zone mass, p10/50/90: [0.99, 0.99, 1.31]
|
||||
wells (bodies) per zone p50/p99: [17.0, 58.0]
|
||||
containment violations (term outside an ancestor support): 0 of 67366
|
||||
terms per zone p50/p99: [29.0, 108.0]
|
||||
sim sample (uniform in the zone disc): terms in support p50/p99/mean [3.0, 12.0]/4.38; tree tests p50/p99 [9.0, 33.0]
|
||||
shader texel (8192 u window): terms in support p50/p99 [3.0, 11.0], mean per texel p50/p99 over windows [2.71, 10.21]; terms meeting the window p50/p99 [19.0, 70.0]
|
||||
== [5] (t2) tapered draft, Lam = 1*D (92s, 6 draws)
|
||||
Sun ratio p10/50/90/99 [0.84, 1.05, 1.27, 1.52] parent-only v p10/50/90 [21.99, 23.6, 24.63] n 2211
|
||||
GasGiant ratio p10/50/90/99 [0.79, 1.0, 1.25, 1.62] parent-only v p10/50/90 [17.58, 22.08, 24.48] n 3566
|
||||
Planet ratio p10/50/90/99 [0.83, 0.97, 1.14, 1.45] parent-only v p10/50/90 [14.9, 16.9, 23.2] n 7445
|
||||
Planetoid ratio p10/50/90/99 [0.85, 0.97, 1.13, 1.43] parent-only v p10/50/90 [10.98, 13.74, 18.14] n 10200
|
||||
Belt ratio p10/50/90/99 [0.84, 0.98, 1.17, 1.47] parent-only v p10/50/90 [11.0, 14.24, 19.45] n 3134
|
||||
planets round a median sun: v p10/50/90 [16.07, 17.16, 17.78]; share in 12-18 0.99; n 1116
|
||||
Sun-class children parent-only v p10/50/90 [21.99, 23.6, 24.63]
|
||||
children with no parent pull: {}; parent pull < 0.5 c/D: {'Planetoid': 214, 'Belt': 11, 'Planet': 9} (total 234)
|
||||
mid-zone inward pull by Z: 1000: 0.58, 2125: 0.40, 4182: 0.25, 6846: 0.18, 10000: 0.13
|
||||
zone far-field c / c(M_zone) p10/50/90: at 2Z [1.0, 1.0, 1.0], at 10Z [1.0, 1.0, 1.0]
|
||||
split invariance (family merged into one body), max |da|/|a| p50/p90/p99: at 1.05 B_f [0.0, 0.1, 0.54], at 2Z [0.0, 0.0, 0.0]; n 6976; families whose reach the merge leaves unchanged 5382/6976, of those inside the ruled half-tidal window 2427/2427
|
||||
whole-zone depth over one body of the zone mass, p10/50/90: [0.99, 0.99, 1.34]
|
||||
wells (bodies) per zone p50/p99: [17.0, 58.0]
|
||||
containment violations (term outside an ancestor support): 0 of 67366
|
||||
terms per zone p50/p99: [29.0, 108.0]
|
||||
sim sample (uniform in the zone disc): terms in support p50/p99/mean [4.0, 13.0]/5.16; tree tests p50/p99 [10.0, 33.0]
|
||||
shader texel (8192 u window): terms in support p50/p99 [3.0, 12.0], mean per texel p50/p99 over windows [3.56, 10.84]; terms meeting the window p50/p99 [21.0, 70.0]
|
||||
== [6] (t3) tapered draft, Lam = 2*d_Hill (91s, 6 draws)
|
||||
Sun ratio p10/50/90/99 [0.83, 1.05, 1.27, 1.52] parent-only v p10/50/90 [21.99, 23.6, 24.63] n 2211
|
||||
GasGiant ratio p10/50/90/99 [0.77, 0.99, 1.24, 1.63] parent-only v p10/50/90 [17.58, 22.08, 24.48] n 3566
|
||||
Planet ratio p10/50/90/99 [0.83, 0.97, 1.14, 1.45] parent-only v p10/50/90 [14.9, 16.9, 23.2] n 7445
|
||||
Planetoid ratio p10/50/90/99 [0.85, 0.98, 1.13, 1.44] parent-only v p10/50/90 [10.96, 13.74, 18.14] n 10200
|
||||
Belt ratio p10/50/90/99 [0.84, 0.98, 1.17, 1.47] parent-only v p10/50/90 [10.99, 14.24, 19.45] n 3134
|
||||
planets round a median sun: v p10/50/90 [16.07, 17.16, 17.78]; share in 12-18 0.99; n 1116
|
||||
Sun-class children parent-only v p10/50/90 [21.99, 23.6, 24.63]
|
||||
children with no parent pull: {}; parent pull < 0.5 c/D: {'Planetoid': 214, 'Belt': 11, 'Planet': 9} (total 234)
|
||||
mid-zone inward pull by Z: 1000: 0.56, 2125: 0.40, 4182: 0.25, 6846: 0.17, 10000: 0.13
|
||||
zone far-field c / c(M_zone) p10/50/90: at 2Z [1.0, 1.0, 1.0], at 10Z [1.0, 1.0, 1.0]
|
||||
split invariance (family merged into one body), max |da|/|a| p50/p90/p99: at 1.05 B_f [0.0, 0.1, 0.55], at 2Z [0.0, 0.0, 0.0]; n 6976; families whose reach the merge leaves unchanged 5295/6976, of those inside the ruled half-tidal window 2427/2427
|
||||
whole-zone depth over one body of the zone mass, p10/50/90: [0.99, 0.99, 1.31]
|
||||
wells (bodies) per zone p50/p99: [17.0, 58.0]
|
||||
containment violations (term outside an ancestor support): 0 of 67366
|
||||
terms per zone p50/p99: [29.0, 108.0]
|
||||
sim sample (uniform in the zone disc): terms in support p50/p99/mean [4.0, 13.0]/5.07; tree tests p50/p99 [10.0, 34.0]
|
||||
shader texel (8192 u window): terms in support p50/p99 [3.0, 12.0], mean per texel p50/p99 over windows [3.56, 10.59]; terms meeting the window p50/p99 [20.0, 70.0]
|
||||
== [7] (w1) conserving compact, Lam = 0.5*D (87s, 6 draws)
|
||||
Sun ratio p10/50/90/99 [0.86, 1.01, 1.24, 1.55] parent-only v p10/50/90 [21.99, 23.6, 24.62] n 2211
|
||||
GasGiant ratio p10/50/90/99 [0.83, 1.0, 1.24, 1.63] parent-only v p10/50/90 [17.58, 22.08, 24.48] n 3566
|
||||
Planet ratio p10/50/90/99 [0.87, 0.99, 1.14, 1.45] parent-only v p10/50/90 [14.9, 16.9, 23.2] n 7445
|
||||
Planetoid ratio p10/50/90/99 [0.89, 0.99, 1.12, 1.42] parent-only v p10/50/90 [11.01, 13.75, 18.14] n 10200
|
||||
Belt ratio p10/50/90/99 [0.87, 0.99, 1.15, 1.45] parent-only v p10/50/90 [11.05, 14.24, 19.45] n 3134
|
||||
planets round a median sun: v p10/50/90 [16.07, 17.17, 17.78]; share in 12-18 0.99; n 1116
|
||||
Sun-class children parent-only v p10/50/90 [21.99, 23.6, 24.62]
|
||||
children with no parent pull: {}; parent pull < 0.5 c/D: {'Planetoid': 214, 'Belt': 11, 'Planet': 9} (total 234)
|
||||
mid-zone inward pull by Z: 1000: 0.55, 2125: 0.39, 4182: 0.24, 6846: 0.17, 10000: 0.13
|
||||
zone far-field c / c(M_zone) p10/50/90: at 2Z [1.0, 1.0, 1.0], at 10Z [1.0, 1.0, 1.0]
|
||||
split invariance (family merged into one body), max |da|/|a| p50/p90/p99: at 1.05 B_f [0.52, 1.63, 5.34], at 2Z [0.0, 0.0, 0.0]; n 6976; families whose reach the merge leaves unchanged 2937/6976, of those inside the ruled half-tidal window 2075/2427
|
||||
whole-zone depth over one body of the zone mass, p10/50/90: [0.99, 1.0, 1.32]
|
||||
wells (bodies) per zone p50/p99: [17.0, 58.0]
|
||||
containment violations (term outside an ancestor support): 0 of 70890
|
||||
terms per zone p50/p99: [30.0, 114.0]
|
||||
sim sample (uniform in the zone disc): terms in support p50/p99/mean [4.0, 14.0]/4.92; tree tests p50/p99 [10.0, 34.0]
|
||||
shader texel (8192 u window): terms in support p50/p99 [3.0, 13.0], mean per texel p50/p99 over windows [2.88, 12.3]; terms meeting the window p50/p99 [21.0, 76.0]
|
||||
== [8] (w2) conserving compact, Lam = 1*D (87s, 6 draws)
|
||||
Sun ratio p10/50/90/99 [0.86, 1.09, 1.32, 1.63] parent-only v p10/50/90 [21.99, 23.6, 24.62] n 2211
|
||||
GasGiant ratio p10/50/90/99 [0.78, 1.01, 1.3, 1.68] parent-only v p10/50/90 [17.58, 22.08, 24.48] n 3566
|
||||
Planet ratio p10/50/90/99 [0.82, 0.97, 1.18, 1.49] parent-only v p10/50/90 [14.9, 16.9, 23.2] n 7445
|
||||
Planetoid ratio p10/50/90/99 [0.83, 0.97, 1.15, 1.45] parent-only v p10/50/90 [11.01, 13.75, 18.14] n 10200
|
||||
Belt ratio p10/50/90/99 [0.83, 0.98, 1.19, 1.49] parent-only v p10/50/90 [11.05, 14.24, 19.45] n 3134
|
||||
planets round a median sun: v p10/50/90 [16.07, 17.17, 17.78]; share in 12-18 0.99; n 1116
|
||||
Sun-class children parent-only v p10/50/90 [21.99, 23.6, 24.62]
|
||||
children with no parent pull: {}; parent pull < 0.5 c/D: {'Planetoid': 214, 'Belt': 11, 'Planet': 9} (total 234)
|
||||
mid-zone inward pull by Z: 1000: 0.61, 2125: 0.43, 4182: 0.27, 6846: 0.19, 10000: 0.14
|
||||
zone far-field c / c(M_zone) p10/50/90: at 2Z [1.0, 1.0, 1.0], at 10Z [1.0, 1.0, 1.0]
|
||||
split invariance (family merged into one body), max |da|/|a| p50/p90/p99: at 1.05 B_f [0.0, 0.89, 3.15], at 2Z [0.0, 0.0, 0.0]; n 6976; families whose reach the merge leaves unchanged 4708/6976, of those inside the ruled half-tidal window 2380/2427
|
||||
whole-zone depth over one body of the zone mass, p10/50/90: [0.99, 1.01, 1.36]
|
||||
wells (bodies) per zone p50/p99: [17.0, 58.0]
|
||||
containment violations (term outside an ancestor support): 0 of 70890
|
||||
terms per zone p50/p99: [30.0, 114.0]
|
||||
sim sample (uniform in the zone disc): terms in support p50/p99/mean [5.0, 16.0]/5.71; tree tests p50/p99 [10.0, 34.0]
|
||||
shader texel (8192 u window): terms in support p50/p99 [3.0, 14.0], mean per texel p50/p99 over windows [3.71, 12.84]; terms meeting the window p50/p99 [21.0, 76.0]
|
||||
== [9] (w3) conserving compact, Lam = 2*d_Hill (88s, 6 draws)
|
||||
Sun ratio p10/50/90/99 [0.84, 1.09, 1.32, 1.62] parent-only v p10/50/90 [21.99, 23.6, 24.62] n 2211
|
||||
GasGiant ratio p10/50/90/99 [0.77, 1.0, 1.29, 1.67] parent-only v p10/50/90 [17.58, 22.08, 24.48] n 3566
|
||||
Planet ratio p10/50/90/99 [0.81, 0.97, 1.17, 1.5] parent-only v p10/50/90 [14.9, 16.9, 23.2] n 7445
|
||||
Planetoid ratio p10/50/90/99 [0.84, 0.97, 1.15, 1.45] parent-only v p10/50/90 [11.01, 13.75, 18.14] n 10200
|
||||
Belt ratio p10/50/90/99 [0.83, 0.98, 1.19, 1.48] parent-only v p10/50/90 [11.05, 14.24, 19.45] n 3134
|
||||
planets round a median sun: v p10/50/90 [16.07, 17.17, 17.78]; share in 12-18 0.99; n 1116
|
||||
Sun-class children parent-only v p10/50/90 [21.99, 23.6, 24.62]
|
||||
children with no parent pull: {}; parent pull < 0.5 c/D: {'Planetoid': 214, 'Belt': 11, 'Planet': 9} (total 234)
|
||||
mid-zone inward pull by Z: 1000: 0.60, 2125: 0.42, 4182: 0.27, 6846: 0.19, 10000: 0.14
|
||||
zone far-field c / c(M_zone) p10/50/90: at 2Z [1.0, 1.0, 1.0], at 10Z [1.0, 1.0, 1.0]
|
||||
split invariance (family merged into one body), max |da|/|a| p50/p90/p99: at 1.05 B_f [0.0, 0.91, 3.16], at 2Z [0.0, 0.0, 0.0]; n 6976; families whose reach the merge leaves unchanged 4573/6976, of those inside the ruled half-tidal window 2382/2427
|
||||
whole-zone depth over one body of the zone mass, p10/50/90: [0.99, 1.01, 1.34]
|
||||
wells (bodies) per zone p50/p99: [17.0, 58.0]
|
||||
containment violations (term outside an ancestor support): 0 of 70890
|
||||
terms per zone p50/p99: [30.0, 114.0]
|
||||
sim sample (uniform in the zone disc): terms in support p50/p99/mean [4.0, 15.0]/5.62; tree tests p50/p99 [10.0, 34.0]
|
||||
shader texel (8192 u window): terms in support p50/p99 [3.0, 14.0], mean per texel p50/p99 over windows [3.69, 12.61]; terms meeting the window p50/p99 [21.0, 76.0]
|
||||
354
docs/research/gravity-well-law/probe3.py
Normal file
354
docs/research/gravity-well-law/probe3.py
Normal file
|
|
@ -0,0 +1,354 @@
|
|||
# Probe of the conserving compact well law (well_law.py) against today's law and the smooth-wells rows a2, c2, b1,
|
||||
# on the measures of ../gravity-smooth-wells/probe2.py (same data, rng seed 7, 6 draws, same ratio/mid-zone/rail code
|
||||
# paths), plus: split invariance (one family merged into one body of its mass), whole-zone merge depth, and
|
||||
# evaluation cost (terms per sample for the sim, per texel for a 2048^2 / 4 u Surface Displacement window).
|
||||
# Rows t1-t3 are the first draft (windows tapered from the centre), kept because they show what the band form buys.
|
||||
# Usage: python probe3.py bodies.csv [draws] [law-index,...]
|
||||
import csv, sys, os, collections, time
|
||||
import numpy as np
|
||||
HERE = os.path.dirname(os.path.abspath(__file__))
|
||||
sys.path.insert(0, HERE)
|
||||
sys.path.insert(0, os.path.join(HERE, '..', 'gravity-smooth-wells'))
|
||||
sys.path.insert(0, os.path.join(HERE, '..', 'gravity-tail-taming'))
|
||||
import well_law as WL
|
||||
from laws2 import (F_today, phi_today, F_logratio, phi_logratio, F_lamb_halo, phi_lamb_halo,
|
||||
F_splummer, phi_splummer, strength, core_radius)
|
||||
|
||||
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 zone_struct(z): # as probe2.zone_struct
|
||||
idx = {int(r['idx']): k for k, r in enumerate(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]
|
||||
return par, D, m, body, cls
|
||||
|
||||
def family(par):
|
||||
n = len(par); fam = np.eye(n, dtype=bool); 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]
|
||||
return fam, depth
|
||||
|
||||
def place(par, D, depth):
|
||||
pos = np.zeros((len(par), 2)); th = rng.uniform(0, 2 * np.pi, len(par))
|
||||
for k in np.argsort(depth, kind='stable'):
|
||||
if par[k] >= 0: pos[k] = pos[par[k]] + D[k] * np.array([np.cos(th[k]), np.sin(th[k])])
|
||||
return pos
|
||||
|
||||
def q(v, ps=(10, 50, 90)):
|
||||
v = np.sort(np.asarray(v, float))
|
||||
return [round(float(v[min(len(v) - 1, int(p / 100 * len(v)))]), 2) for p in ps] if len(v) else []
|
||||
|
||||
# ---------------- the compact law: terms on the orbit tree ----------------
|
||||
# A term is (owner node, strength, eps, a, b): strength x band term U(eps, a, b) at the owner's position.
|
||||
def _tree(par, D, m, body):
|
||||
n = len(par); fam, depth = family(par)
|
||||
kids = [[j for j in range(n) if par[j] == i] for i in range(n)]
|
||||
mb = np.where(body, m, 0.0)
|
||||
MF = np.array([mb[fam[i]].sum() for i in range(n)])
|
||||
cm = np.where(body, strength(np.maximum(m, 1e-12)), 0.0); em = core_radius(np.maximum(m, 1e-12))
|
||||
CF = np.where(MF > 0, strength(np.maximum(MF, 1e-12)), 0.0); eM = core_radius(np.maximum(MF, 1e-12))
|
||||
def ring_of(p):
|
||||
pk = [j for j in kids[p] if MF[j] > 0]
|
||||
dmin = min(D[j] for j in pk)
|
||||
return [j for j in pk if abs(D[j] - dmin) <= 0.1]
|
||||
return n, fam, depth, kids, mb, MF, cm, em, CF, eM, ring_of
|
||||
|
||||
def build_terms_band(par, D, m, body, reach):
|
||||
"""The law: band terms (well_law.py header). reach(node, D, C_F, c_parent_body) is a child's nominal reach."""
|
||||
n, fam, depth, kids, mb, MF, cm, em, CF, eM, ring_of = _tree(par, D, m, body)
|
||||
Lam = np.zeros(n); L = np.zeros(n); B = np.zeros(n); grp = {}
|
||||
for i in sorted(range(n), key=lambda i: -depth[i]): # children before parents
|
||||
L[i] = max([D[j] for j in kids[i]], default=0.0)
|
||||
B[i] = max([WL.BETA * L[i]] + [D[j] + Lam[j] for j in kids[i]])
|
||||
if MF[i] > 0 and not body[i]: # bodyless: its ring group stands in for a body
|
||||
rg = ring_of(i); LG = D[rg[0]]
|
||||
BG = max([WL.BETA * LG] + [D[j] + Lam[j] for j in rg])
|
||||
L[i] = max(L[i], BG); B[i] = max(B[i], WL.BETA * L[i])
|
||||
grp[i] = (MF[rg].sum(), LG, BG)
|
||||
if MF[i] <= 0:
|
||||
Lam[i] = 0.0; continue
|
||||
p = par[i]
|
||||
floor = WL.BETA * B[i] # room for the fade beyond the family's own support
|
||||
if p < 0:
|
||||
Lam[i] = np.inf
|
||||
elif not body[p] and i in ring_of(p) and len(ring_of(p)) >= 2:
|
||||
Lam[i] = max(WL.RING_REACH * D[i], floor)
|
||||
else:
|
||||
Lam[i] = max(reach(i, D[i], CF[i], cm[p] if body[p] else 0.0), floor)
|
||||
terms = []
|
||||
for i in range(n):
|
||||
if MF[i] <= 0: continue
|
||||
fa, fb = (np.inf, np.inf) if np.isinf(Lam[i]) else (max(B[i], WL.GAMMA * Lam[i]), Lam[i])
|
||||
if body[i] and MF[i] <= mb[i] * (1 + 1e-12): # no massive descendants: one term
|
||||
terms.append((i, cm[i], em[i], fa, fb)); continue
|
||||
if body[i]:
|
||||
terms.append((i, cm[i], em[i], L[i], B[i]))
|
||||
else:
|
||||
MG, LG, BG = grp[i]; cg, eg = strength(MG), core_radius(MG)
|
||||
terms += [(i, cg, eg, L[i], B[i]), (i, -cg, eg, LG, BG)]
|
||||
terms += [(i, CF[i], eM[i], fa, fb), (i, -CF[i], eM[i], L[i], B[i])]
|
||||
return terms, Lam, L, B
|
||||
|
||||
def build_terms_tapered(par, D, m, body, reach):
|
||||
"""The first draft: windows tapered from the centre (well_law.py, unit_pull_tapered), KAPPA 2.5, MU 1, RING 5."""
|
||||
KAPPA, MU, RING = 2.5, 1.0, 5.0
|
||||
n, fam, depth, kids, mb, MF, cm, em, CF, eM, ring_of = _tree(par, D, m, body)
|
||||
Lam = np.zeros(n); L = np.zeros(n)
|
||||
for i in sorted(range(n), key=lambda i: -depth[i]):
|
||||
L[i] = max([D[j] + Lam[j] for j in kids[i]], default=0.0)
|
||||
if MF[i] <= 0:
|
||||
Lam[i] = 0.0; continue
|
||||
p = par[i]
|
||||
floor = max(MU * L[i], KAPPA * max([D[j] for j in kids[i]], default=0.0))
|
||||
if p < 0:
|
||||
Lam[i] = np.inf
|
||||
elif not body[p] and i in ring_of(p) and len(ring_of(p)) >= 2:
|
||||
Lam[i] = max(RING * D[i], floor)
|
||||
else:
|
||||
Lam[i] = max(reach(i, D[i], CF[i], cm[p] if body[p] else 0.0), floor)
|
||||
T = []
|
||||
def t(i, c, e, ell): T.append((i, c, e, np.inf if np.isinf(ell) else 0.0, ell))
|
||||
for i in range(n):
|
||||
if MF[i] <= 0: continue
|
||||
if body[i]:
|
||||
if MF[i] <= mb[i] * (1 + 1e-12):
|
||||
t(i, cm[i], em[i], Lam[i])
|
||||
else:
|
||||
t(i, cm[i], em[i], L[i]); t(i, CF[i], eM[i], Lam[i]); t(i, -CF[i], eM[i], L[i])
|
||||
else:
|
||||
rg = ring_of(i); MG = MF[rg].sum(); LG = max(D[j] + Lam[j] for j in rg)
|
||||
if LG < L[i]:
|
||||
t(i, strength(MG), core_radius(MG), L[i]); t(i, -strength(MG), core_radius(MG), LG)
|
||||
t(i, CF[i], eM[i], Lam[i]); t(i, -CF[i], eM[i], L[i])
|
||||
return T, Lam, L, L
|
||||
|
||||
def _arr(T):
|
||||
return (np.array([t[0] for t in T]), np.array([t[1] for t in T]), np.array([t[2] for t in T]),
|
||||
np.array([t[3] for t in T]), np.array([t[4] for t in T]))
|
||||
|
||||
def compact_field(pos, T, x):
|
||||
own, c, e, a, b = _arr(T)
|
||||
d = pos[own][None, :, :] - x[:, None, :]; r = np.linalg.norm(d, axis=2); rs = np.maximum(r, 1e-9)
|
||||
F = c[None, :] * WL.band_pull(e[None, :], a[None, :], b[None, :], rs)
|
||||
return d / rs[:, :, None] * F[:, :, None], own
|
||||
|
||||
def compact_potential(pos, T, x):
|
||||
own, c, e, a, b = _arr(T)
|
||||
r = np.linalg.norm(pos[own][None, :, :] - x[:, None, :], axis=2)
|
||||
return (c[None, :] * WL.band_potential(e[None, :], a[None, :], b[None, :], r)).sum(1)
|
||||
|
||||
def own_pull(T, node, d):
|
||||
return sum(t[1] * float(WL.band_pull(t[2], t[3], t[4], np.array(d))) for t in T if t[0] == node)
|
||||
|
||||
# ---------------- the comparison laws: one well per body ----------------
|
||||
OTHER = {
|
||||
'today': dict(F=F_today, P=phi_today, scale=lambda D, Z: np.ones_like(D)),
|
||||
'a2': dict(F=F_logratio, P=phi_logratio, scale=lambda D, Z: 1.0 * D),
|
||||
'c2': dict(F=F_lamb_halo, P=phi_lamb_halo, scale=lambda D, Z: 1.0 * D),
|
||||
'b1': dict(F=F_splummer, P=phi_splummer, scale=lambda D, Z: 0.5 * D),
|
||||
}
|
||||
def wells_field(law, pos, c, e, s, x):
|
||||
d = pos[None, :, :] - x[:, None, :]; r = np.linalg.norm(d, axis=2); rs = np.maximum(r, 1e-9)
|
||||
F = np.where(r < 1e-9, 0, law['F'](c[None, :], e[None, :], rs, s[None, :]))
|
||||
return d / rs[:, :, None] * F[:, :, None]
|
||||
def wells_potential(law, pos, c, e, s, x):
|
||||
r = np.maximum(np.linalg.norm(pos[None, :, :] - x[:, None, :], axis=2), 1e-9)
|
||||
return np.where(c[None, :] > 0, law['P'](c[None, :], e[None, :], r, s[None, :]), 0).sum(1)
|
||||
|
||||
TIDAL = lambda i, D, CF, cp: (2 * D * np.sqrt(CF / (2 * cp)) if cp > 0 else 0.5 * D)
|
||||
HALF = lambda i, D, CF, cp: 0.5 * D
|
||||
LAWS = [
|
||||
('today: harmonic core + log', 'today', None),
|
||||
('(a2) log-ratio L=1*D', 'a2', None),
|
||||
('(c2) Lamb-Oseen + halo L=1*D', 'c2', None),
|
||||
('(b1) screened Plummer lam=0.5*D', 'b1', None),
|
||||
('(t1) tapered draft, Lam = 0.5*D', 'tapered', HALF),
|
||||
('(t2) tapered draft, Lam = 1*D', 'tapered', lambda i, D, CF, cp: 1.0 * D),
|
||||
('(t3) tapered draft, Lam = 2*d_Hill', 'tapered', TIDAL),
|
||||
('(w1) conserving compact, Lam = 0.5*D', 'band', HALF),
|
||||
('(w2) conserving compact, Lam = 1*D', 'band', lambda i, D, CF, cp: 1.0 * D),
|
||||
('(w3) conserving compact, Lam = 2*d_Hill', 'band', TIDAL),
|
||||
]
|
||||
pick = [int(x) for x in sys.argv[3].split(',')] if len(sys.argv) > 3 else range(len(LAWS))
|
||||
|
||||
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))]
|
||||
pre = {k: zone_struct(z) for k, z in zones.items()}
|
||||
TH8 = np.arange(8) * np.pi / 4; RING8 = np.stack([np.cos(TH8), np.sin(TH8)], 1)
|
||||
|
||||
class ZoneLaw:
|
||||
"""One law on one zone's tree: field, potential, the parent's own pull, and which node owns each contribution."""
|
||||
def __init__(self, kind, reach, par, D, m, body, Z):
|
||||
self.kind, self.par = kind, par
|
||||
self.fam, self.depth = family(par)
|
||||
if kind in ('band', 'tapered'):
|
||||
self.T, self.Lam, self.L, self.B = (build_terms_band if kind == 'band' else build_terms_tapered)(par, D, m, body, reach)
|
||||
else:
|
||||
self.law = OTHER[kind]
|
||||
self.bi = np.where(body)[0]
|
||||
self.c = strength(np.maximum(m, 1e-12))[self.bi]; self.e = core_radius(np.maximum(m, 1e-12))[self.bi]
|
||||
Dz = np.where(par < 0, Z, D)
|
||||
self.s = self.law['scale'](Dz, Z)[self.bi]
|
||||
def field(self, pos, x, exclude_family_of=None):
|
||||
if self.kind in ('band', 'tapered'):
|
||||
A, own = compact_field(pos, self.T, x)
|
||||
else:
|
||||
A = wells_field(self.law, pos[self.bi], self.c, self.e, self.s, x); own = self.bi
|
||||
if exclude_family_of is not None:
|
||||
mask = self.fam[np.asarray(exclude_family_of)][:, own]
|
||||
A = np.where(mask[:, :, None], 0, A)
|
||||
return A.sum(1)
|
||||
def potential(self, pos, x):
|
||||
if self.kind in ('band', 'tapered'): return compact_potential(pos, self.T, x)
|
||||
return wells_potential(self.law, pos[self.bi], self.c, self.e, self.s, x)
|
||||
def own(self, node, d):
|
||||
if self.kind in ('band', 'tapered'): return own_pull(self.T, node, d)
|
||||
j = int(np.where(self.bi == node)[0][0])
|
||||
return float(self.law['F'](self.c[j], self.e[j], np.array(d), self.s[j]))
|
||||
|
||||
def merged(par, D, m, body, f):
|
||||
"""The zone with node f's descendants removed and f's body mass made the family's."""
|
||||
fam, _ = family(par)
|
||||
keep = [j for j in range(len(par)) if not (fam[f, j] and j != f)]
|
||||
remap = {j: k for k, j in enumerate(keep)}
|
||||
mb = np.where(body, m, 0.0)
|
||||
par2 = np.array([remap[par[j]] if par[j] >= 0 else -1 for j in keep])
|
||||
m2 = m[keep].copy(); m2[remap[f]] = mb[fam[f]].sum()
|
||||
body2 = body[keep].copy(); body2[remap[f]] = True
|
||||
return par2, D[keep], m2, body2, keep, remap
|
||||
|
||||
for li in pick:
|
||||
name, kind, reach = LAWS[li]
|
||||
t0 = time.time()
|
||||
ratio = collections.defaultdict(list); vpar = collections.defaultdict(list); vmed = []; vsun = []
|
||||
inward = collections.defaultdict(list); farc = {2: [], 10: []}; weak = collections.Counter(); nopull = collections.Counter()
|
||||
split_err = {'near': [], 'far': []}; split_exact = [0, 0]; split_window = [0, 0]; depth_ratio = []
|
||||
cost_sim = []; cost_sim_nodes = []; cost_tex = []; cost_tex_mean = []; cost_cull = []; nwells = []; nterms = []
|
||||
violations = [0, 0]
|
||||
for key, z in zones.items():
|
||||
par, D, m, body, cls = pre[key]
|
||||
Z = float(z[0]['zoneRadius'])
|
||||
zn = ZoneLaw(kind, reach, par, D, m, body, Z)
|
||||
fam, depth = zn.fam, zn.depth
|
||||
kids = [k for k in range(len(par)) if par[k] >= 0 and body[par[k]]]
|
||||
cbody = np.where(body, strength(np.maximum(m, 1e-12)), 0.0)
|
||||
a_own = {k: zn.own(par[k], D[k]) for k in kids}
|
||||
for k in kids:
|
||||
v = np.sqrt(max(a_own[k], 0) * D[k]); vpar[cls[k]].append(v)
|
||||
if a_own[k] <= 0: nopull[cls[k]] += 1
|
||||
if a_own[k] < 0.5 * cbody[par[k]] / D[k]: weak[cls[k]] += 1
|
||||
if cls[k] == 'Planet' and z[par[k]]['type'] == 'Sun' and SUN_LO <= m[par[k]] <= SUN_HI: vmed.append(v)
|
||||
if cls[k] == 'Sun': vsun.append(v)
|
||||
mz = np.where(body, m, 0).sum()
|
||||
acc_rel = collections.defaultdict(float); mids = []; farc[2].append(0.0); farc[10].append(0.0)
|
||||
for dr in range(DRAWS):
|
||||
pos = place(par, D, depth)
|
||||
Aout = zn.field(pos, pos, exclude_family_of=np.arange(len(par)))
|
||||
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
|
||||
X = 0.5 * Z * RING8
|
||||
AX = zn.field(pos, X)
|
||||
mids.append(np.mean(-(AX * X).sum(1) / np.linalg.norm(X, axis=1)))
|
||||
for fz in (2, 10):
|
||||
XF = fz * Z * RING8; AF = zn.field(pos, XF)
|
||||
farc[fz][-1] += np.mean(-(AF * XF).sum(1)) / DRAWS / strength(mz)
|
||||
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]))
|
||||
|
||||
# ---- split invariance: each non-root body family with massive descendants, merged into one body ----
|
||||
pos = place(par, D, depth)
|
||||
mb = np.where(body, m, 0.0)
|
||||
Lref = zn.B if kind in ('band', 'tapered') else ZoneLaw('band', HALF, par, D, m, body, Z).B
|
||||
for f in range(len(par)):
|
||||
if par[f] < 0 or not body[f] or mb[fam[f]].sum() <= mb[f]: continue
|
||||
par2, D2, m2, body2, keepi, remap = merged(par, D, m, body, f)
|
||||
z2 = ZoneLaw(kind, reach, par2, D2, m2, body2, Z)
|
||||
pos2 = pos[keepi]
|
||||
for tag, X in (('near', pos[f] + 1.05 * Lref[f] * RING8), ('far', 2 * Z * RING8)):
|
||||
a1 = zn.field(pos, X); a2 = z2.field(pos2, X)
|
||||
split_err[tag].append(np.max(np.linalg.norm(a1 - a2, axis=1) / np.maximum(np.linalg.norm(a2, axis=1), 1e-12)))
|
||||
if kind in ('band', 'tapered'):
|
||||
same = int(zn.Lam[f] == z2.Lam[remap[f]])
|
||||
split_exact[1] += 1; split_exact[0] += same
|
||||
# families the ruled generator window would hold: every descendant within half its parent's tidal
|
||||
# radius (B25: d_H = D sqrt(c / 2 c_parent); a parent with no body parent has no window and is skipped)
|
||||
ok = True
|
||||
for j in range(len(par)):
|
||||
if j == f or not fam[f, j]: continue
|
||||
p = par[j]
|
||||
if body[p] and par[p] >= 0 and body[par[p]]:
|
||||
ok &= D[j] <= 0.5 * D[p] * np.sqrt(cbody[p] / (2 * cbody[par[p]]))
|
||||
if ok:
|
||||
split_window[1] += 1; split_window[0] += same
|
||||
# ---- whole zone merged: depth of the deepest point against one body of the zone's mass ----
|
||||
g = np.linspace(-Z, Z, 33); G = np.concatenate([np.stack(np.meshgrid(g, g), -1).reshape(-1, 2), pos[body]])
|
||||
ref = zn.potential(pos, 10 * Z * RING8).mean()
|
||||
z1 = ZoneLaw(kind, reach, np.array([-1]), np.array([0.0]), np.array([mz]), np.array([True]), Z)
|
||||
p1 = np.zeros((1, 2)); ref1 = z1.potential(p1, 10 * Z * RING8).mean()
|
||||
depth_ratio.append((ref - zn.potential(pos, G).min()) / (ref1 - z1.potential(p1, p1).min()))
|
||||
|
||||
# ---- evaluation cost ----
|
||||
nwells.append(int(body.sum()))
|
||||
if kind in ('band', 'tapered'):
|
||||
T = zn.T; own = np.array([t[0] for t in T]); ell = np.array([t[4] for t in T])
|
||||
nterms.append(len(T))
|
||||
# containment: every term lies inside each ancestor's own support L (path-length bound)
|
||||
for t in T:
|
||||
a, path = t[0], 0.0
|
||||
while par[a] >= 0:
|
||||
path += D[a]; a = par[a]
|
||||
violations[1] += 1
|
||||
violations[0] += int(path + t[4] > zn.B[a] * (1 + 1e-9))
|
||||
for dr in range(2):
|
||||
pos = place(par, D, depth)
|
||||
rr = Z * np.sqrt(rng.uniform(0, 1, 64)); th = rng.uniform(0, 2 * np.pi, 64)
|
||||
S = np.stack([rr * np.cos(th), rr * np.sin(th)], 1)
|
||||
dist = np.linalg.norm(pos[own][None, :, :] - S[:, None, :], axis=2)
|
||||
cost_sim += list((dist < ell[None, :]).sum(1))
|
||||
distn = np.linalg.norm(pos[None, :, :] - S[:, None, :], axis=2)
|
||||
inside = distn < zn.Lam[None, :]
|
||||
vis = np.zeros_like(inside)
|
||||
for k in np.argsort(depth, kind='stable'):
|
||||
vis[:, k] = inside[:, k] & (vis[:, par[k]] if par[k] >= 0 else True)
|
||||
ntests = np.zeros(64)
|
||||
for k in range(len(par)):
|
||||
if par[k] >= 0: ntests += vis[:, par[k]]
|
||||
cost_sim_nodes += list(ntests)
|
||||
cx = Z * np.sqrt(rng.uniform()) * np.array([np.cos(th[0]), np.sin(th[0])])
|
||||
gx = cx[0] + np.linspace(-4096, 4096, 64); gy = cx[1] + np.linspace(-4096, 4096, 64)
|
||||
P = np.stack(np.meshgrid(gx, gy), -1).reshape(-1, 2)
|
||||
dist = np.linalg.norm(pos[own][None, :, :] - P[:, None, :], axis=2)
|
||||
cnt = (dist < ell[None, :]).sum(1)
|
||||
cost_tex += list(cnt); cost_tex_mean.append(cnt.mean())
|
||||
dx = np.maximum(np.abs(pos[own][:, 0] - cx[0]) - 4096, 0); dy = np.maximum(np.abs(pos[own][:, 1] - cx[1]) - 4096, 0)
|
||||
cost_cull.append(int((np.hypot(dx, dy) < ell).sum()))
|
||||
|
||||
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])}')
|
||||
vmed = np.array(vmed)
|
||||
print(f' planets round a median sun: v p10/50/90 {q(vmed)}; share in 12-18 {np.mean((vmed >= 12) & (vmed <= 18)):.2f}; n {len(vmed)}')
|
||||
print(f' Sun-class children parent-only v p10/50/90 {q(vsun)}')
|
||||
print(' children with no parent pull: ' + str(dict(nopull)) + '; parent pull < 0.5 c/D: ' + str(dict(weak)) + f' (total {sum(weak.values())})')
|
||||
print(' mid-zone inward pull by Z: ' + ', '.join(f'{Z:.0f}: {np.median(v):.2f}' for Z, v in sorted(inward.items())))
|
||||
print(' zone far-field c / c(M_zone) p10/50/90: at 2Z ' + str(q(farc[2])) + ', at 10Z ' + str(q(farc[10])))
|
||||
print(f' split invariance (family merged into one body), max |da|/|a| p50/p90/p99: at 1.05 B_f {q(split_err["near"], (50, 90, 99))}, at 2Z {q(split_err["far"], (50, 90, 99))}; n {len(split_err["near"])}'
|
||||
+ (f'; families whose reach the merge leaves unchanged {split_exact[0]}/{split_exact[1]}, of those inside the ruled half-tidal window {split_window[0]}/{split_window[1]}' if kind in ('band', 'tapered') else ''))
|
||||
print(f' whole-zone depth over one body of the zone mass, p10/50/90: {q(depth_ratio)}')
|
||||
print(f' wells (bodies) per zone p50/p99: {q(nwells, (50, 99))}')
|
||||
if kind in ('band', 'tapered'):
|
||||
print(f' containment violations (term outside an ancestor support): {violations[0]} of {violations[1]}')
|
||||
print(f' terms per zone p50/p99: {q(nterms, (50, 99))}')
|
||||
print(f' sim sample (uniform in the zone disc): terms in support p50/p99/mean {q(cost_sim, (50, 99))}/{np.mean(cost_sim):.2f}; tree tests p50/p99 {q(cost_sim_nodes, (50, 99))}')
|
||||
print(f' shader texel (8192 u window): terms in support p50/p99 {q(cost_tex, (50, 99))}, mean per texel p50/p99 over windows {q(cost_tex_mean, (50, 99))}; terms meeting the window p50/p99 {q(cost_cull, (50, 99))}')
|
||||
sys.stdout.flush()
|
||||
2
docs/research/gravity-well-law/selfpull-output.txt
Normal file
2
docs/research/gravity-well-law/selfpull-output.txt
Normal file
|
|
@ -0,0 +1,2 @@
|
|||
tapered draft: the node pulls the lone body at 0.380 of its own pull
|
||||
band: the node pulls the lone body at 0.000 of its own pull
|
||||
15
docs/research/gravity-well-law/selfpull.py
Normal file
15
docs/research/gravity-well-law/selfpull.py
Normal file
|
|
@ -0,0 +1,15 @@
|
|||
# Self-pull check: a lone body (mass 100) at 200 u on a bodyless node that orbits a sun at 3000 u. Its rail reads the
|
||||
# node's own terms; a law that counts the body's own mass in the node's family term pulls it toward the node.
|
||||
# Reuses probe3.py's term builders (its definitions only; the probe's main loop is not run).
|
||||
import os, sys
|
||||
import numpy as np
|
||||
HERE = os.path.dirname(os.path.abspath(__file__))
|
||||
src = open(os.path.join(HERE, 'probe3.py'), encoding='utf-8').read()
|
||||
g = {'__file__': os.path.join(HERE, 'probe3.py')}
|
||||
exec(src[:src.index('rng = np.random')], g); g['rng'] = np.random.default_rng(7)
|
||||
exec(src[src.index('def zone_struct'):src.index('pick = [int')], g)
|
||||
par = np.array([-1, 0, 1]); D = np.array([0., 3000., 200.]); m = np.array([5000., 0., 100.]); body = np.array([True, False, True])
|
||||
own = g['strength'](100.0) * 200 / (200 ** 2 + g['core_radius'](100.0) ** 2)
|
||||
for name, build in (('tapered draft', g['build_terms_tapered']), ('band', g['build_terms_band'])):
|
||||
T = build(par, D, m, body, g['TIDAL'])[0]
|
||||
print(f'{name}: the node pulls the lone body at {g["own_pull"](T, 1, 200.0) / own:.3f} of its own pull')
|
||||
104
docs/research/gravity-well-law/well_law.py
Normal file
104
docs/research/gravity-well-law/well_law.py
Normal file
|
|
@ -0,0 +1,104 @@
|
|||
# The conserving compact well law (aetheria-gravity, Imagination pass imagination-ag3-well-law, 2026-10-15).
|
||||
#
|
||||
# One primitive, the band term U(r; eps, a, b): a Plummer-softened log source whose pull is exact up to a and fades
|
||||
# to zero over [a, b] by an order-k smoothstep in s = r^2:
|
||||
# pull (toward the centre) F(r) = r / (r^2 + eps^2) * B(s)
|
||||
# B(s) = 1 for s <= a^2; 1 - S(t), t = (s - a^2) / (b^2 - a^2), on the band; 0 for s >= b^2
|
||||
# S(t) = t^k sum_{j<k} C(k-1+j, j) (1-t)^j (the order-k smoothstep: S - 0 and 1 - S vanish to order k at the ends)
|
||||
# potential: 0 for s >= b^2; on the band 1/2 [R ln((t+tau)/(1+tau)) - (I(1) - I(t))], tau = (a^2+eps^2)/(b^2-a^2),
|
||||
# where (1 - S(t)) = (t + tau) Q(t) + R and I is the antiderivative of Q; inside, 1/2 ln((s+eps^2)/(a^2+eps^2)) plus
|
||||
# the band's value at t = 0. a = b = inf is the global term 1/2 ln(r^2 + eps^2).
|
||||
# The potential is a function of r^2 (C-infinity at the centre and wherever B is constant) and C^k on the circles
|
||||
# r = a and r = b (the k-th derivative of the pull jumps there).
|
||||
#
|
||||
# A zone's potential is a sum of such terms on the orbit tree (see the probe's build_terms): every node F (a body, or
|
||||
# a bodyless node's ring group) at its own position contributes
|
||||
# c(m_F) U(eps(m_F), L_F, B_F) + c(M_F) [ U(eps(M_F), max(B_F, GAMMA Lam_F), Lam_F) - U(eps(M_F), L_F, B_F) ]
|
||||
# m_F is the body's mass, M_F the family's. L_F is the outermost child's distance: inside it the node pulls exactly
|
||||
# as its body alone. B_F >= every member term's reach: beyond it the family pulls exactly as one body of mass M_F,
|
||||
# out to max(B_F, GAMMA Lam_F), and is gone at its reach Lam_F >= BETA B_F. A root's Lam is infinite, so the zone's far field is exactly
|
||||
# c(M_zone). A leaf is one term, c(m) U(eps(m), GAMMA Lam, Lam) (B = 0).
|
||||
#
|
||||
# The first draft, the tapered window W(u) = 1 - S(u) applied from the centre (pull r W(r^2/l^2) / (r^2 + eps^2)), is
|
||||
# kept as unit_pull_tapered / unit_potential_tapered: its family charge leaks inside the family (a child feels a share
|
||||
# of its own mass through its parent's family term), which the band form removes.
|
||||
import math
|
||||
import numpy as np
|
||||
|
||||
K_ORDER = 4 # the smoothstep's order: the potential is C^4 on every band edge
|
||||
BETA = 1.5 # a family's charge is fully on by B_F >= BETA x its outermost child
|
||||
GAMMA = 0.5 # a family is felt whole out to max(B_F, GAMMA x its reach) and fades from there to its reach
|
||||
# (its reach is at least BETA x B_F, so the fade is never narrower than half of B_F)
|
||||
RING_REACH = 4.0 # a ring member of a bodyless node reaches RING_REACH x its distance: whole out to the ring's diameter
|
||||
|
||||
def smoothstep_coeffs(k=K_ORDER):
|
||||
"""S(t) as polynomial coefficients in t, lowest first."""
|
||||
out = np.zeros(2 * k)
|
||||
for j in range(k): # C(k-1+j, j) t^k (1-t)^j
|
||||
cj = math.comb(k - 1 + j, j)
|
||||
for i in range(j + 1):
|
||||
out[k + i] += cj * math.comb(j, i) * (-1) ** i
|
||||
return out
|
||||
|
||||
SC = smoothstep_coeffs()
|
||||
FADE = -SC.copy(); FADE[0] += 1.0 # 1 - S(t)
|
||||
|
||||
def fade(t):
|
||||
t = np.clip(t, 0.0, 1.0)
|
||||
return np.polyval(FADE[::-1], t)
|
||||
|
||||
def band_pull(eps, a, b, r):
|
||||
"""Pull toward the centre per unit strength."""
|
||||
r = np.asarray(r, float); s = r * r
|
||||
a = np.asarray(a, float); b = np.asarray(b, float)
|
||||
inf = np.isinf(b)
|
||||
A2 = np.where(inf, 0.0, a * a); B2 = np.where(inf, 1.0, b * b)
|
||||
t = (s - A2) / np.maximum(B2 - A2, 1e-30)
|
||||
prof = np.where(inf, 1.0, np.where(s <= A2, 1.0, np.where(s >= B2, 0.0, fade(t))))
|
||||
return r * prof / (s + np.asarray(eps) ** 2)
|
||||
|
||||
def _band_coeffs(tau):
|
||||
"""Q (t + tau) + R = 1 - S(t): returns the coefficients q_i (arrays over tau) and R."""
|
||||
n = len(FADE) - 1
|
||||
q = [None] * n
|
||||
q[n - 1] = FADE[n] * np.ones_like(tau)
|
||||
for i in range(n - 1, 0, -1): # synthetic division by (t + tau): q_{i-1} = f_i - tau q_i
|
||||
q[i - 1] = FADE[i] - tau * q[i]
|
||||
R = FADE[0] - tau * q[0]
|
||||
return q, R
|
||||
|
||||
def _qint(q, t): # I(t) = sum q_i t^(i+1) / (i+1), Horner
|
||||
acc = np.zeros_like(t)
|
||||
for i in range(len(q) - 1, -1, -1):
|
||||
acc = (acc + q[i] / (i + 1)) * t
|
||||
return acc
|
||||
|
||||
def band_potential(eps, a, b, r):
|
||||
"""Potential per unit strength: 0 at and beyond b; the global term (b = inf) is 1/2 ln(r^2 + eps^2)."""
|
||||
r = np.asarray(r, float); s = r * r
|
||||
eps = np.asarray(eps, float); a = np.asarray(a, float); b = np.asarray(b, float)
|
||||
shape = np.broadcast(r, eps, a, b).shape
|
||||
s, eps, a, b = (np.broadcast_to(x, shape) for x in (s, eps, a, b))
|
||||
inf = np.isinf(b)
|
||||
A2 = np.where(inf, 0.0, a * a); B2 = np.where(inf, 2.0, b * b)
|
||||
span = B2 - A2
|
||||
tau = (A2 + eps * eps) / span
|
||||
q, R = _band_coeffs(tau)
|
||||
def band(t):
|
||||
return 0.5 * (R * np.log((t + tau) / (1 + tau)) - (_qint(q, np.ones_like(t)) - _qint(q, t)))
|
||||
t = np.clip((s - A2) / span, 0.0, 1.0)
|
||||
inside = 0.5 * np.log((s + eps * eps) / (A2 + eps * eps)) + band(np.zeros_like(t))
|
||||
val = np.where(s <= A2, inside, np.where(s >= B2, 0.0, band(t)))
|
||||
return np.where(inf, 0.5 * np.log(s + eps * eps), val)
|
||||
|
||||
# ---- the first draft: tapered window from the centre (see the header) ----
|
||||
def unit_pull_tapered(eps, ell, r):
|
||||
r = np.asarray(r, float); ell = np.asarray(ell, float)
|
||||
inf = np.isinf(ell)
|
||||
u = np.where(inf, 0.0, r * r / np.where(inf, 1.0, ell) ** 2)
|
||||
win = np.where(inf, 1.0, np.where(u < 1, fade(u), 0.0))
|
||||
return r * win / (r * r + np.asarray(eps) ** 2)
|
||||
|
||||
def unit_potential_tapered(eps, ell, r):
|
||||
"""The tapered term is the band term with a = 0: the same closed form."""
|
||||
return band_potential(eps, np.zeros_like(np.asarray(ell, float)), ell, r)
|
||||
Loading…
Reference in a new issue