GameCult/Aetheria
Watch
9
0
Fork
You've already forked Aetheria
 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Merge gravity smooth-wells research

Browse source 
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Meta Vens 2026-10-15 21:44:31 +02:00
commit 8bc5ef1837
8 changed files with 635 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

129
docs/research/gravity-smooth-wells.md Normal file
View file

@ -0,0 +1,129 @@
# Smooth wells: C-infinity laws for the grid and the rails
Eyes pass `eyes-ag3-smooth-wells`, campaign aetheria-gravity, session self-2026-10-15-ag3, 2026-10-15, read at master 39ebf7f6. Facts with evidence pointers; no recommendation. Extends `gravity-tail-taming.md` (same data, metrics, phases and table format) under question `aetheria-gravity:question:tail-taming-law` and ruling `aetheria-gravity:ruling:tame-tails-parent-only-rails`.
Operator, verbatim: "Since both orbits and well shape are negotiable, is there a solution available that's more mathematically elegant? 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."
Scripts in `gravity-smooth-wells/`:
- `laws2.py`: the laws (force and potential). It imports c, the core radius, today's law and the A&S K0/K1 from `../gravity-tail-taming/laws.py`.
- `check_laws2.py`: force = dφ/dr for each law, the shader form of Ein against mpmath, and circular-speed shapes. Output: `check-output.txt` (first lines) and the text below.
- `seams.py`: jumps in φ…φ''' at each piecewise seam, with each branch differentiated analytically in mpmath at 50 digits. Output: `seams-output.txt`.
- `probe2.py`: `../gravity-tail-taming/probe.py` with whole-law pulls. Output: `probe2-output.txt`.
Every law has c = 45.08-calibrated strength (c_sun = 304.9 at the median sun) and ε = today's Roche core radius 7.32·m^0.25. Law 0 of `probe2.py` reproduces the tail-taming baseline exactly (Sun ratio p50 2.01, mid-zone 1.59 … 0.51). Each law ran in 9-16 s on Starfire, one process per law.
## 1. Smoothness class
A radial potential φ(r) is C^∞ on the plane exactly when it is a smooth function of r². Each law below is a function of r² with no seam, so it is real-analytic everywhere, the core included. The A&S polynomial forms are piecewise, so "the HLSL form" is a separate column.
| Law | φ on the plane (exact) | HLSL/C# form | Seam location |
|---|---|---|---|
| Today: harmonic core + log | C^1 | C^1 | r = R. φ'' jumps by 2c/R² (inside +c/R², outside −c/R²); φ''' jumps by 2c/R³ |
| Tail-taming K0 rows (harmonic core + K0) | C^1 (core seam) | C^1 at R, plus the A&S seam | r = R, and r = 2λ |
| Tail-taming cap | C^1 at R and at L (φ'' jumps by c/L² at L) | same | |
| Tail-taming force switch | C^1 at R; C^2 at L and 2L | same | |
| (a) log-ratio | C^ω | C^∞: elementary (log, div), no branch | none |
| (b) screened Plummer | C^ω | C^0 in the formula: A&S piecewise | √(r²+ε²) = 2λ |
| (c) Lamb-Oseen, with or without halo | C^ω (Ein is entire) | φ: A&S E1 is piecewise. Force c(1−e^(−r²/ε²))/r is elementary and C^∞ | φ seam at r = ε (u = r²/ε² = 1) |
| (d) Plummer log | C^ω | C^∞: one log, no branch | none |
| (e) conserving halo | C^ω | C^∞: two logs, no branch | none |
**Measured seam jumps** (`seams-output.txt`; jump = high branch − low branch; "rel" = |jump| / |true derivative|):
| Seam | φ | φ' | φ'' | φ''' |
|---|---|---|---|---|
| A&S K0 at x = 2, in K0^(n) (multiply by −c/λ^n) | 7.6e-9 (rel 6.7e-8) | 2.5e-7 (rel 1.8e-6) | 8.4e-7 (rel 4.6e-6) | 7.7e-5 (rel 2.9e-4) |
| A&S K1 at x = 2 (force path), in K1^(n) (multiply by c/λ^(n+1)) | 1.5e-9 (rel 1.1e-8) | 1.9e-7 (rel 1.0e-6) | 5.1e-6 (rel 1.9e-5) | 4.6e-5 (rel 1.0e-4) |
| A&S E1 → Lamb-Oseen φ at r = ε, units c/ε^n | 2.2e-9 (rel 5.5e-9) | 2.2e-6 (rel 3.5e-6) | 1.3e-4 (rel 1.2e-3) | 5.1e-3 (rel 5.4e-3) |
| Today's core at r = R, units c/R^n | 0 | 0 | 2.0 (rel 2.0) | 2.0 (rel 1.0) |
On the same units (c/ε², with ε = R), the Lamb-Oseen shader seam's φ'' jump is 6.4e-5 of today's. The K0 seam's φ'' jump is 8.4e-7·c/λ². With λ ≥ 0.5·D ≫ R, that is below 1e-6 of today's 2c/R². The K0 value jump at the seam (6.7e-8 relative) is about one float32 ulp of K0(2) = 0.114 (float32 ε = 6e-8).
## 2. Probe (same measures as `gravity-tail-taming.md` Question 4; 6 draws, rng seed 7)
**Family (e), conserving halo** (added on the coordinator's brief, after the operator's observation, verbatim: "this is all a consequence of exaggerating the contributions of small bodies... 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."):
- φ_i = (c_near/2)·ln(r²+ε²) − ((c_near − c_far)/2)·ln(r²+L²).
- c_near is today's compressed 45.08-calibrated m_i^0.2 strength.
- c_far,i = c(M_F)·m_i/M_F, where F is the body's parent family (the parent and all its descendants; for the root, the whole zone) and only non-belt, non-empty bodies carry mass.
- L = k·ext, where ext is the body's own family extent: the longest summed-distance path from the body to any descendant. A body with no descendants falls back to D (Z for a root).
- c_far/c_near = (m_i/M_F)^0.8 ≤ 1, so the halo strength is never negative.
**Far-field measure** (every law): the inward pull at 8 points on a circle of radius 2Z (and 10Z) round the zone centre, times the radius, averaged over the 6 draws. It is divided by c(M_zone), the strength of the zone's body mass merged into one body under the same calibration.
"outer" = the body's outermost direct-child distance (any child type, belts and empties included); a body with no children falls back to D (Z for a root). The Sun class in "Ratio" and in the rail-speed line is **Sun-class children, i.e. suns orbiting a sun** (`probe.py` keys by the child's class). `gravity-tail-taming.md` labels the same numbers "children of root and rosette suns", and that label does not match its code.
| # | 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 | Sun-class children's parent-only rail v p10/p50/p90 | Children with no parent pull | Children with parent pull < 0.5·c/D | Zone far-field c / c(M_zone), p50/p90 at 2Z; at 10Z |
|---|---|---|---|---|---|---|---|---|
| 0 | today (harmonic core + log) | 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 | 22.1/23.6/24.6 | 0 | 46 | 5.54/7.98; 5.54/7.98 |
| d | Plummer log, no taming | 2.01/2.46 · 1.63/2.22 · 1.25/1.72 · 1.13/1.44 · 1.19/1.55 | 1.57 / 1.26 / 0.91 / 0.67 / 0.51 | 16.1/17.2/17.8; 0.99 | 22.0/23.6/24.6 | 0 | 234 | 5.54/7.98; 5.54/7.98 |
| a1 | log-ratio L = 0.5·D | 1.17/1.57 · 1.04/1.37 · 0.99/1.19 · 0.99/1.16 · 0.99/1.24 | 0.34 / 0.24 / 0.15 / 0.11 / 0.08 | 14.4/16.3/17.1; 0.96 | 12.3/15.4/19.5 | 18 (L ≤ ε for 285 bodies: repulsive) | 5838 | 0.10/0.11; 0.00/0.00 |
| a2 | log-ratio L = 1·D | 1.26/1.55 · 1.12/1.43 · 1.02/1.26 · 1.01/1.20 · 1.02/1.26 | 0.64 / 0.45 / 0.29 / 0.20 / 0.15 | 15.5/16.9/17.5; 1.00 | 17.7/19.9/22.7 | 0 (L ≤ ε for 34 bodies) | 407 | 0.36/0.40; 0.02/0.02 |
| a3 | log-ratio L = 1.5·outer | 1.17/1.41 · 1.12/1.39 · 1.05/1.31 · 1.03/1.27 · 1.05/1.32 | 0.74 / 0.46 / 0.27 / 0.18 / 0.14 | 14.1/15.2/16.7; 1.00 | 18.4/20.6/23.1 | 0 (7 bodies) | 637 | 0.34/0.48; 0.02/0.02 |
| a4 | log-ratio L = 2·outer | 1.22/1.46 · 1.16/1.44 · 1.07/1.35 · 1.04/1.28 · 1.08/1.35 | 0.85 / 0.54 / 0.32 / 0.22 / 0.16 | 15.1/15.9/17.0; 1.00 | 19.8/22.0/23.7 | 0 (7 bodies) | 326 | 0.52/0.72; 0.03/0.04 |
| b1 | screened Plummer λ = 0.5·D | 1.16/1.49 · 1.04/1.34 · 0.99/1.19 · 0.99/1.16 · 1.00/1.22 | 0.41 / 0.29 / 0.18 / 0.13 / 0.09 | 14.7/16.3/17.0; 0.98 | 14.3/17.2/20.5 | 0 | 3610 | 0.07/0.08; 0.00/0.00 |
| b2 | screened Plummer λ = 1·D | 1.28/1.56 · 1.15/1.46 · 1.04/1.29 · 1.02/1.22 · 1.04/1.28 | 0.70 / 0.50 / 0.32 / 0.23 / 0.16 | 15.6/16.8/17.5; 1.00 | 18.7/20.8/22.9 | 0 | 297 | 0.45/0.50; 0.00/0.00 |
| c0 | Lamb-Oseen core, log tail | 2.01/2.46 · 1.62/2.21 · 1.24/1.71 · 1.11/1.43 · 1.18/1.54 | 1.59 / 1.27 / 0.91 / 0.67 / 0.51 | 17.0/17.4/18.0; 0.93 | 22.1/23.6/24.6 | 0 | 121 | 5.54/7.98; 5.54/7.98 |
| c1 | Lamb-Oseen + halo L = 0.5·D | 1.17/1.57 · 1.04/1.37 · 0.99/1.19 · 0.99/1.15 · 0.99/1.22 | 0.35 / 0.24 / 0.16 / 0.11 / 0.08 | 14.8/16.6/17.5; 0.97 | 12.3/15.4/19.6 | 1 | 4944 | 0.10/0.12; 0.00/0.00 |
| c2 | Lamb-Oseen + halo L = 1·D | 1.26/1.55 · 1.12/1.44 · 1.02/1.26 · 1.00/1.19 · 1.02/1.26 | 0.66 / 0.46 / 0.29 / 0.21 / 0.15 | 16.5/17.2/17.8; 0.99 | 17.8/20.0/22.8 | 0 | 146 | 0.36/0.40; 0.02/0.02 |
| c3 | Lamb-Oseen + halo L = 1.5·outer | 1.17/1.41 · 1.11/1.39 · 1.04/1.31 · 1.02/1.26 · 1.05/1.32 | 0.75 / 0.46 / 0.28 / 0.18 / 0.14 | 14.3/15.8/17.2; 1.00 | 18.4/20.6/23.1 | 0 | 136 | 0.34/0.48; 0.02/0.02 |
| c4 | Lamb-Oseen + halo L = 2·outer | 1.22/1.46 · 1.15/1.44 · 1.06/1.34 · 1.04/1.27 · 1.07/1.34 | 0.87 / 0.54 / 0.32 / 0.22 / 0.16 | 15.4/16.4/17.5; 1.00 | 19.8/22.0/23.7 | 0 | 122 | 0.52/0.73; 0.03/0.04 |
| e1 | conserving halo L = 1.5·ext | 1.25/1.50 · 1.17/1.44 · 1.06/1.33 · 1.04/1.28 · 1.07/1.34 | 0.86 / 0.58 / 0.37 / 0.25 / 0.19 | 14.5/15.9/16.9; 1.00 | 21.7/23.4/24.5 | 0 (7 bodies) | 431 | 1.42/1.99; 1.28/1.84 |
| e2 | conserving halo L = 2·ext | 1.30/1.55 · 1.20/1.49 · 1.09/1.37 · 1.05/1.30 · 1.09/1.36 | 0.95 / 0.64 / 0.40 / 0.27 / 0.21 | 15.3/16.3/17.1; 1.00 | 21.8/23.5/24.5 | 0 (7 bodies) | 270 | 1.52/2.07; 1.28/1.85 |
Notes from `probe2-output.txt`:
- *Core alone (rows d, c0 against 0)*: replacing the harmonic core with a C^∞ core moves no ratio or mid-zone figure by more than 0.02. It lowers rail speeds for children near their parent's core. Planetoids' p10 drops 11.5 → 11.0 (Plummer) and 11.4 (Lamb-Oseen). The count of children with a parent pull below half the log value rises from 46 to 234 (Plummer) and 121 (Lamb-Oseen).
- *Screening with the smooth core (b) against the tail-taming K0 rows (harmonic core)*: λ = 0.5·D ratio p50s 1.16/1.04/0.99/0.99/1.00, against 1.15/1.04/0.99/0.99/0.99. λ = 1·D gives 1.28/1.15/1.04/1.02/1.04, against 1.28/1.14/1.03/1.01/1.03. The mid-zone figures differ by ≤ 0.03. The median-sun planets go from 16.7 to 16.3 (0.5·D) and from 17.2 to 16.8 (1·D).
- *L = 0.5·D*: L ≤ ε for 285 of 27,165 bodies, where the log-ratio force c·r·(L²−ε²)/(…) is repulsive. 18 children then have no inward parent pull, and the ratio n drops by 2-12 per class.
- *Split-deepening (far-field column)*:
- *Today, Plummer and Lamb-Oseen*: the zone's far field is Σc_i = 5.54× (p50) and 7.98× (p90) the merged body's c(M_zone), the same at 2Z and 10Z. This is the n^(1−p) effect the operator names.
- *Screened and halo laws (a, b, c')*: at 2Z the zone's far field is 0.07-0.52× c(M_zone), and at 10Z ≤ 0.04×. Their far field is neutral by construction.
- *(e)*: 1.42×/1.99× at 2Z and 1.28×/1.84× at 10Z (p50/p90), not 1.0. The static sum Σc_far,i / c(M_zone) over each zone gives the same 1.02/1.27/1.84 (p10/50/90; one root in each of 1500 zones). The excess comes from nesting. A moon's c_far is its share of its planet family's compressed c(M_F), and the compressed strength per unit mass is higher for a smaller family. So the shares sum to c(M_F) only within one family level, not across the zone. If every body instead took its share of the zone total, c(M_zone)·m_i/M_zone, the shares would sum to exactly c(M_zone).
- *(e) rails*:
- It keeps today's Sun-class rails (21.7/23.4/24.5 against 22.1/23.6/24.6), because c_near governs inside L.
- It cuts the Sun-class ratio p50 from 2.01 to 1.25 (1.5·ext) or 1.30 (2·ext), and the mid-zone pull from 1.59 to 0.86-0.95 at Z = 1000.
- Planets round a median sun: 15.9 and 16.3 p50, all in 12-18.
- *Planets round all suns*, share in 12-18 u/s: 0.51 today, 0.53-0.64 for the smooth laws (highest: log-ratio 1.5·outer at 0.64).
## 3. Closed forms, circular speed and shader cost
Pull F is toward the body; v² = r·F; s = √(r²+ε²); u = r²/ε². Cost is per well per texel. The force vector is F·Δ/r, with Δ the offset to the well, so "no sqrt" below means the vector form needs no square root.
| Law | φ | F | v²(r) | Plateau v ≈ √c | φ cost (grid) | Force-vector cost |
|---|---|---|---|---|---|---|
| Today | c r²/2R² (r<R); c(ln(r/R)+½) | c r/R² ; c/r | c r²/R² ; c | exactly √c for all r ≥ R | 1 log (of r²), 1 select | 1 div, 1 select, no sqrt |
| (d) Plummer | (c/2)·ln(r²+ε²) | c r/(r²+ε²) | c r²/(r²+ε²) | v ≥ 0.9√c from 2.07ε, rising to √c at infinity; 0.71√c at ε | 1 log | 1 div, no sqrt |
| (a) log-ratio | (c/2)·ln[(r²+ε²)/(r²+L²)] | c r(L²−ε²)/[(r²+ε²)(r²+L²)] | c r²(L²−ε²)/[(r²+ε²)(r²+L²)] | peak v² = c(L−ε)/(L+ε) at r = √(εL); v ≥ 0.9√c requires L ≥ 9.5ε. Median sun (ε = 79.8): L = 1000 gives a peak of 0.923√c, ≥0.9 on 185-431; L = 4000 gives a peak of 0.980√c, ≥0.9 on 166-1926 | 1 log + 1 div (or 2 logs) | 1 div (2 FMA products), no sqrt |
| (b) screened Plummer | −c·K0(s/λ) | (c/λ)·K1(s/λ)·r/s | (c r²/(λ s))·K1(s/λ) | peak 0.923√c at r = 300 for λ = 1000 (median sun); ≥0.9 on 190-485 | 1 sqrt + A&S K0: low branch 1 log + ~13 FMA, high branch 1 exp + 1 rsqrt + 1 div + 7 FMA; branchless does both | 1 sqrt + A&S K1 (same shape) + 1 div |
| (c) Lamb-Oseen | (c/2)·Ein(u), Ein(u) = E1(u)+ln u+γ | c(1−e^(−u))/r | c(1−e^(−u)) | v ≥ 0.9√c from 1.29ε; 0.80√c at ε, 0.99√c at 2ε | A&S: low branch a degree-5 polynomial in u (the ln cancels); high branch 1 exp + 2 div + 1 log; branchless does both | 1 exp + 1 div, no sqrt, no branch |
| (c') Lamb-Oseen + halo | (c/2)·[Ein(u) − ln(r²+L²)] | c(1−e^(−u))/r − c r/(r²+L²) | c[(1−e^(−u)) − r²/(r²+L²)] | median sun, L = 1000: peak 0.981√c at r = 181; ≥0.9 on 105-484 | as (c) + 1 log | 1 exp + 2 div, no sqrt |
| (e) conserving halo | (c_n/2)·ln(r²+ε²) − ((c_n−c_f)/2)·ln(r²+L²) | c_n r/(r²+ε²) − (c_n−c_f) r/(r²+L²) | c_n r²/(r²+ε²) − (c_n−c_f) r²/(r²+L²) | two plateaus: v → √c_n for ε ≪ r ≪ L and v → √c_f for r ≫ L; the halo term removes (c_n−c_f)/2 of v² at r = L | 2 logs | 1 div (single fraction over (r²+ε²)(r²+L²)), no sqrt; two constants per well (c_f, L) beyond today's |
Further facts:
- *Peak pull* (units c/ε): today 1.0 (at R), Lamb-Oseen 0.638 (at r ≈ 1.12ε), Plummer 0.5 (at r = ε).
- *Far fields*: (e) keeps a log tail of strength c_far, F → c_f/r. (a) and (c') fall as φ ≈ −(c/2)(L²−ε²)/r², F ≈ c(L²−ε²)/r³, which is algebraic. (b) falls as e^(−r/λ)/√r. (c) and (d) keep the log tail.
- *Sources (2D Laplacian of φ)*: (d) 2cε²/(r²+ε²)², the surface density of a projected Plummer sphere. (c) a Gaussian, (2c/ε²)·e^(−u). (a) the same Plummer core minus a Plummer halo of width L with equal strength, so the net source is zero ("neutral"). (c') a Gaussian core minus the same halo.
- *Gauss's law*: (a), (b), (c') and (e) are not harmonic outside the core: the halo (a, c') or screening (b) is spread over all r. A ring of siblings therefore no longer acts as its centre. (c) and (d) are harmonic to within e^(−u) and ε²/r² respectively, so the enclosed-strength funnel stays: rows c0 and d equal row 0.
- *Check*: `check_laws2.py` shows max |F − dφ/dr| / max F of 1.9e-7 (Plummer, log-ratio, conserving halo), 1.2e-5 (screened Plummer, through the A&S pair), 5e-6 (Lamb-Oseen, through A&S E1), and 1.1e-4 for today (the finite difference straddles the kink at R). The shader form of Ein agrees with mpmath to 2.3e-7 absolute over u from 1e-4 to 100.
## 4. Literature
- **Plummer softening in 2D.** The potential (c/2)·ln(r²+ε²) is the log potential of a projected Plummer sphere, Σ ∝ ε²/(r²+ε²)² (Plummer 1911, MNRAS 71:460; Binney & Tremaine, *Galactic Dynamics*, §2.2). Krasny 1986 (J. Comput. Phys. 65:292) regularised 2D point vortices by adding δ² to the Biot-Savart denominator, which gives the kernel r/(r²+δ²) of row d (the "vortex blob" kernel).
- **Kuzmin disc.** Kuzmin 1956 (Astron. Zh. 33:27) and Toomre 1963 (ApJ 138:385): a razor-thin disc with Σ ∝ a/(R²+a²)^(3/2), whose in-plane potential under 3D gravity is −GM/√(R²+a²). It is the 3D-gravity counterpart of Plummer softening, not a 2D log law. It shares the r²+a² form but not the kernel.
- **Lamb-Oseen vortex.** Oseen 1912 and Lamb, *Hydrodynamics* (6th ed. 1932, §334a): a Gaussian vorticity core with swirl velocity (Γ/2πr)(1−e^(−r²/r_c²)). This is row c's force law, and Gaussian blobs are a standard choice in vortex methods. Beale & Majda 1985 (J. Comput. Phys. 58:188, DOI 10.1016/0021-9991(85)90176-7) build smooth explicit 2D/3D velocity kernels of arbitrary order. I did not confirm from the abstract whether they compare against Gaussian blobs.
- **Charge-neutral smoothing.** Ewald 1921 (Ann. Phys. 369:253) adds to each point charge an opposite Gaussian cloud of equal charge, so the real-space sum (erfc kernel) is short-ranged. In 2D, a point log minus an equal Gaussian is −(c/2)·E1(r²/σ²), which is row c' with a Gaussian halo in place of the Plummer halo. Wolf et al. 1999 (J. Chem. Phys. 110:8254) neutralise each cutoff sphere for a pairwise 1/r sum. Debye & Hückel 1923 (Phys. Z. 24:185): a mobile neutralising cloud screens each charge. In 2D the screened Green's function is K0, which is row b.
- **Dehnen 2001** (MNRAS 324:273, astro-ph/0011568): Plummer softening gives significantly larger force errors than kernels of finite extent. "Compensating" kernels (e.g. his K1) reduce the error further by exceeding the Newtonian force near the softening length, which offsets the too-small force inside. None of the laws here overshoot: (d) and (c) approach c/r from below, and (a), (b), (c') and (e) stay below c_near/r at all r. Dehnen's error measure is a 3D N-body sampling bias and is not this probe's metric.
## Confidence
- §1: high. The seam jumps are branch derivatives at 50 digits, and the smoothness classes follow from analyticity in r².
- §2: high on internal consistency (law 0 reproduces the tail-taming baseline digit for digit). Medium on absolute ratios, for the same caveats as `gravity-tail-taming.md` (6 phase draws, phase-averaged circular speeds, generator positions). "outer" is my reading of "outermost child distance of the body's family" (direct children only).
- §3: high on formulas and checks. Op counts are by inspection, not compiled HLSL.
- §4: high for Plummer, Lamb-Oseen, Ewald, Debye-Hückel, Krasny 1986 and the Dehnen 2001 abstract claims. Medium on the Kuzmin 1956 and Lamb section citations and on the Wolf 1999 page. Dehnen's kernel definitions were not read beyond the abstract.
## Sources
- Abramowitz & Stegun §5.1.53-5.1.54 (E1), §9.8.1-9.8.8 (I0, I1, K0, K1).
- Dehnen 2001: https://arxiv.org/abs/astro-ph/0011568 ; its K1 kernel's use downstream: https://arxiv.org/pdf/1310.4502
- Beale & Majda 1985: https://scholars.duke.edu/publication/759540 ; Beale & Majda 1982, Math. Comp. 39:29: https://www.ams.org/mcom/1982-39-159/S0025-5718-1982-0658213-7/
- Krasny 1986: https://dept.math.lsa.umich.edu/~krasny/paper_jcp_1986.pdf
- Plummer 1911; Kuzmin 1956; Toomre 1963; Lamb 1932; Ewald 1921; Wolf et al. 1999; Debye & Hückel 1923, as cited above.

8
docs/research/gravity-smooth-wells/check-output.txt Normal file
View file

@ -0,0 +1,8 @@
today max |F - dphi/dr| / max F = 1.1e-04
plummer max |F - dphi/dr| / max F = 1.9e-07
logratio L=800 max |F - dphi/dr| / max F = 1.9e-07
splummer lam=800 max |F - dphi/dr| / max F = 1.2e-05
lamb max |F - dphi/dr| / max F = 5.0e-06
lamb+halo L=800 max |F - dphi/dr| / max F = 5.1e-06
conserve L=800 max |F - dphi/dr| / max F = 1.9e-07
Ein shader form vs mpmath: max abs err 2.3e-07

32
docs/research/gravity-smooth-wells/check_laws2.py Normal file
View file

@ -0,0 +1,32 @@
# Checks: force = dphi/dr for every smooth law (grid equals field); Ein shader form against mpmath; circular-speed shape.
import numpy as np, mpmath as mp, os, sys
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
from laws2 import *
c, eps = 304.9, 7.32 * 14152.4 ** 0.25 # median sun
r = np.linspace(0.01, 20000, 400001)
for name, F, P, s in [('today', F_today, phi_today, 0), ('plummer', F_plummer, phi_plummer, 0),
('logratio L=800', F_logratio, phi_logratio, 800.0), ('splummer lam=800', F_splummer, phi_splummer, 800.0),
('lamb', F_lamb, phi_lamb, 0), ('lamb+halo L=800', F_lamb_halo, phi_lamb_halo, 800.0)]:
ph = P(c, eps, r, s); d = np.gradient(ph, r); f = F(c, eps, r, s)
print(f'{name:18s} max |F - dphi/dr| / max F = {np.max(np.abs(f - d)[5:-5]) / np.abs(f).max():.1e}')
us = np.concatenate([np.logspace(-4, 0, 50), np.logspace(0, 2, 50)])
ref = np.array([float(mp.ein(u)) if hasattr(mp, 'ein') else float(mp.e1(u) + mp.log(u) + mp.euler) for u in us])
ph = phi_conserve(c, eps, r, 800.0, 60.0); d = np.gradient(ph, r); f = F_conserve(c, eps, r, 800.0, 60.0)
print(f'{"conserve L=800":18s} max |F - dphi/dr| / max F = {np.max(np.abs(f - d)[5:-5]) / np.abs(f).max():.1e}')
print('Ein shader form vs mpmath: max abs err', f'{np.max(np.abs(ein(us) - ref)):.1e}')
# circular speed v^2 = r F: peak and plateau
print(f'median sun: c = {c}, sqrt(c) = {np.sqrt(c):.2f} u/s, eps = {eps:.1f}')
rr = np.logspace(0, 5, 20001)
for name, F, s in [('today', F_today, 0), ('plummer', F_plummer, 0), ('lamb', F_lamb, 0),
('logratio L=1000', F_logratio, 1000.0), ('logratio L=4000', F_logratio, 4000.0),
('splummer lam=1000', F_splummer, 1000.0), ('lamb+halo L=1000', F_lamb_halo, 1000.0)]:
v = np.sqrt(np.maximum(rr * F(c, eps, rr, s), 0)) / np.sqrt(c)
i = np.argmax(v); band = rr[v >= 0.9]
print(f'{name:18s} v/sqrt(c): peak {v[i]:.3f} at r = {rr[i]:.0f} ({rr[i]/eps:.2f} eps); v >= 0.9 sqrt(c) on r in '
f'[{band.min():.0f}, {band.max():.0f}]' if len(band) else f'{name:18s} peak {v[i]:.3f} at r = {rr[i]:.0f}; never 0.9')
for x in (eps, 2 * eps, 5 * eps):
pass
print('v/sqrt(c) at r = eps, 2eps, 3eps:')
for name, F in [('today', F_today), ('plummer', F_plummer), ('lamb', F_lamb)]:
print(f' {name:8s}', [round(float(np.sqrt(x * F(c, eps, x, 0) / c)), 3) for x in (eps, 2 * eps, 3 * eps)])
print('max pull (units c/eps):', {n: round(float(np.max(F(c, eps, rr, 0)) * eps / c), 3) for n, F in [('today', F_today), ('plummer', F_plummer), ('lamb', F_lamb)]})

66
docs/research/gravity-smooth-wells/laws2.py Normal file
View file

@ -0,0 +1,66 @@
# Smooth (C-infinity) one-well laws (aetheria-gravity, Eyes pass eyes-ag3-smooth-wells, 2026-10-15).
# Each law: pull F(c, eps, r, s) toward the body (u/s^2) for ALL r >= 0 (the core is part of the law, no seam),
# and potential phi(c, eps, r, s) with dphi/dr = F (phi up to a constant). eps = today's Roche core radius 7.32*m^0.25.
# Reuses ../gravity-tail-taming/laws.py for c, eps, the A&S K0/K1 and today's piecewise law.
import os, sys
import numpy as np
sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', 'gravity-tail-taming'))
from laws import k0, k1, strength, core_radius, pull, tail_log, P0 # noqa: F401
EULER = 0.5772156649015329
# ---- exponential integral E1, A&S 5.1.53 (0 < u <= 1) and 5.1.54 (u >= 1): the shader form, piecewise at u = 1 ----
def e1_lo(u): # E1(u) = -ln u + poly(u), |err| < 2e-7
return -np.log(u) + (-0.57721566 + u*(0.99999193 + u*(-0.24991055 + u*(0.05519968 + u*(-0.00976004 + u*0.00107857)))))
def e1_hi(u): # u e^u E1(u) = P4/Q4, |err| < 2e-8
p = 0.2677737343 + u*(8.6347608925 + u*(18.0590169730 + u*(8.5733287401 + u)))
q = 3.9584969228 + u*(21.0996530827 + u*(25.6329561486 + u*(9.5733223454 + u)))
return np.exp(-u) / u * p / q
def ein(u):
"""Ein(u) = integral_0^u (1 - e^-t)/t dt = E1(u) + ln u + gamma; entire. Shader form: A&S pieces, seam at u = 1."""
u = np.asarray(u, float)
ul = np.clip(u, 1e-300, 1.0); uh = np.maximum(u, 1.0)
lo = e1_lo(ul) + np.log(ul) + EULER # the -ln u cancels: lo is a pure polynomial in u
hi = e1_hi(uh) + np.log(uh) + EULER
return np.where(u <= 1, lo, hi)
# ---- laws: F(c, eps, r, s) and phi(c, eps, r, s); s is the per-body scale (L or lambda), unused where noted ----
def F_today(c, eps, r, s): # harmonic core to eps, c/r beyond (C1 potential at eps)
return pull(tail_log, c, eps, r, s)
def phi_today(c, eps, r, s):
r = np.asarray(r, float)
return np.where(r < eps, c * r * r / (2 * eps * eps), c * (np.log(np.maximum(r, 1e-300) / eps) + 0.5))
def F_plummer(c, eps, r, s): # (d) Plummer-softened log, no taming
return c * r / (r * r + eps * eps)
def phi_plummer(c, eps, r, s):
return 0.5 * c * np.log(r * r + eps * eps)
def F_logratio(c, eps, r, L): # (a) neutralised Plummer: Plummer core minus Plummer halo of width L, same strength
r2 = r * r
return c * r * (L * L - eps * eps) / ((r2 + eps * eps) * (r2 + L * L))
def phi_logratio(c, eps, r, L):
return 0.5 * c * np.log((r * r + eps * eps) / (r * r + L * L))
def F_splummer(c, eps, r, lam): # (b) screened Plummer: phi = -c K0(sqrt(r^2+eps^2)/lam)
q = np.sqrt(r * r + eps * eps)
return c / lam * k1(q / lam) * r / q
def phi_splummer(c, eps, r, lam):
return -c * k0(np.sqrt(r * r + eps * eps) / lam)
def F_lamb(c, eps, r, s): # (c) Lamb-Oseen core (Gaussian source), log tail
r = np.asarray(r, float); rs = np.maximum(r, 1e-12)
return c * (-np.expm1(-(rs / eps) ** 2)) / rs
def phi_lamb(c, eps, r, s):
return 0.5 * c * ein((np.asarray(r, float) / eps) ** 2)
def F_lamb_halo(c, eps, r, L): # (c') Lamb-Oseen core minus the log-ratio's Plummer halo at L
return F_lamb(c, eps, r, L) - c * r / (r * r + L * L)
def phi_lamb_halo(c, eps, r, L):
return phi_lamb(c, eps, r, L) - 0.5 * c * np.log(r * r + L * L)
def F_conserve(c, eps, r, L, cf): # (e) conserving halo: compressed c near, mass-additive c_far beyond L
r2 = r * r
return c * r / (r2 + eps * eps) - (c - cf) * r / (r2 + L * L)
def phi_conserve(c, eps, r, L, cf):
return 0.5 * c * np.log(r * r + eps * eps) - 0.5 * (c - cf) * np.log(r * r + L * L)

195
docs/research/gravity-smooth-wells/probe2-output.txt Normal file
View file

@ -0,0 +1,195 @@
== [0] today: harmonic core + log (10s, 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: {}
bodies with scale s <= eps (law today: harmonic core + log): 27165 of 27165
children whose parent pull < 0.5 c/D (v < 0.71 of the log rail): {'Planetoid': 42, 'Belt': 3, 'Planet': 1}
zone far-field c over merged-body c(M_zone), p10/50/90: at 2Z [3.69, 5.54, 7.98], at 10Z [3.69, 5.54, 7.98]
mid-zone inward pull by Z: 1000: 1.59, 2125: 1.27, 4182: 0.91, 6846: 0.67, 10000: 0.51
== [1] (d) Plummer log, no taming (10s, 6 draws)
Sun ratio p10/50/90/99 [1.44, 2.01, 2.46, 2.81] parent-only v p10/50/90 [21.99, 23.6, 24.62] n 2211
GasGiant ratio p10/50/90/99 [1.11, 1.63, 2.22, 2.67] parent-only v p10/50/90 [17.58, 22.08, 24.48] n 3566
Planet ratio p10/50/90/99 [0.97, 1.25, 1.72, 2.15] parent-only v p10/50/90 [14.9, 16.9, 23.2] n 7445
Planetoid ratio p10/50/90/99 [0.95, 1.13, 1.44, 1.77] parent-only v p10/50/90 [11.01, 13.75, 18.14] n 10200
Belt ratio p10/50/90/99 [0.97, 1.19, 1.55, 2.01] parent-only v p10/50/90 [11.05, 14.24, 19.45] n 3134
planets round suns: parent-only v p10/50/90 [15.76, 17.36, 23.33]; share in 12-18 u/s 0.58
planets round a median sun (mass p40-p60, 10607-20254): v p10/50/90 [16.07, 17.17, 17.78]; share in 12-18 0.99; n 1116
children with no parent pull at their distance: {}
bodies with scale s <= eps (law (d) Plummer log, no taming): 27165 of 27165
children whose parent pull < 0.5 c/D (v < 0.71 of the log rail): {'Planetoid': 214, 'Belt': 11, 'Planet': 9}
zone far-field c over merged-body c(M_zone), p10/50/90: at 2Z [3.69, 5.54, 7.98], at 10Z [3.69, 5.54, 7.98]
mid-zone inward pull by Z: 1000: 1.57, 2125: 1.26, 4182: 0.91, 6846: 0.67, 10000: 0.51
== [2] (a) log-ratio L=0.5*D (9s, 6 draws)
Sun ratio p10/50/90/99 [0.89, 1.17, 1.57, 2.01] parent-only v p10/50/90 [12.26, 15.37, 19.45] n 2211
GasGiant ratio p10/50/90/99 [0.86, 1.04, 1.37, 1.87] parent-only v p10/50/90 [12.55, 17.03, 21.52] n 3566
Planet ratio p10/50/90/99 [0.89, 0.99, 1.19, 1.58] parent-only v p10/50/90 [11.11, 15.78, 21.98] n 7443
Planetoid ratio p10/50/90/99 [0.87, 0.99, 1.16, 1.98] parent-only v p10/50/90 [7.68, 12.64, 17.63] n 10188
Belt ratio p10/50/90/99 [0.87, 0.99, 1.24, 2.04] parent-only v p10/50/90 [7.29, 12.73, 17.86] n 3130
planets round suns: parent-only v p10/50/90 [11.39, 16.23, 22.23]; share in 12-18 u/s 0.60
planets round a median sun (mass p40-p60, 10607-20254): v p10/50/90 [14.35, 16.25, 17.08]; share in 12-18 0.96; n 1116
children with no parent pull at their distance: {'Planetoid': 12, 'Belt': 4, 'Planet': 2}
bodies with scale s <= eps (law (a) log-ratio L=0.5*D): 285 of 27165
children whose parent pull < 0.5 c/D (v < 0.71 of the log rail): {'Planetoid': 1902, 'Planet': 955, 'Belt': 643, 'GasGiant': 905, 'Sun': 1433}
zone far-field c over merged-body c(M_zone), p10/50/90: at 2Z [0.07, 0.1, 0.11], at 10Z [0.0, 0.0, 0.0]
mid-zone inward pull by Z: 1000: 0.34, 2125: 0.24, 4182: 0.15, 6846: 0.11, 10000: 0.08
== [3] (a) log-ratio L=1*D (9s, 6 draws)
Sun ratio p10/50/90/99 [0.98, 1.26, 1.55, 1.82] parent-only v p10/50/90 [17.73, 19.93, 22.72] n 2211
GasGiant ratio p10/50/90/99 [0.92, 1.12, 1.43, 1.79] 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.26, 1.58] 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.61] 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.75] parent-only v p10/50/90 [9.83, 13.7, 18.99] n 3134
planets round suns: parent-only v p10/50/90 [14.31, 17.05, 22.93]; share in 12-18 u/s 0.61
planets round a median sun (mass p40-p60, 10607-20254): v p10/50/90 [15.52, 16.86, 17.49]; share in 12-18 1.00; n 1116
children with no parent pull at their distance: {}
bodies with scale s <= eps (law (a) log-ratio L=1*D): 34 of 27165
children whose parent pull < 0.5 c/D (v < 0.71 of the log rail): {'Planetoid': 329, 'Belt': 47, 'Planet': 31}
zone far-field c over merged-body c(M_zone), p10/50/90: at 2Z [0.25, 0.36, 0.4], at 10Z [0.01, 0.02, 0.02]
mid-zone inward pull by Z: 1000: 0.64, 2125: 0.45, 4182: 0.29, 6846: 0.20, 10000: 0.15
== [4] (a) log-ratio L=1.5*outer (9s, 6 draws)
Sun ratio p10/50/90/99 [0.98, 1.17, 1.41, 1.7] parent-only v p10/50/90 [18.39, 20.55, 23.07] n 2211
GasGiant ratio p10/50/90/99 [0.93, 1.12, 1.39, 1.83] parent-only v p10/50/90 [15.02, 19.76, 23.51] n 3566
Planet ratio p10/50/90/99 [0.92, 1.05, 1.31, 1.62] parent-only v p10/50/90 [12.84, 15.44, 22.92] n 7445
Planetoid ratio p10/50/90/99 [0.89, 1.03, 1.27, 1.58] parent-only v p10/50/90 [9.19, 12.35, 18.07] n 10200
Belt ratio p10/50/90/99 [0.89, 1.05, 1.32, 1.72] parent-only v p10/50/90 [9.1, 12.48, 19.08] n 3134
planets round suns: parent-only v p10/50/90 [13.69, 15.9, 23.0]; share in 12-18 u/s 0.64
planets round a median sun (mass p40-p60, 10607-20254): v p10/50/90 [14.13, 15.16, 16.65]; share in 12-18 1.00; n 1116
children with no parent pull at their distance: {}
bodies with scale s <= eps (law (a) log-ratio L=1.5*outer): 7 of 27165
children whose parent pull < 0.5 c/D (v < 0.71 of the log rail): {'Planetoid': 485, 'Belt': 98, 'GasGiant': 9, 'Planet': 45}
zone far-field c over merged-body c(M_zone), p10/50/90: at 2Z [0.26, 0.34, 0.48], at 10Z [0.01, 0.02, 0.02]
mid-zone inward pull by Z: 1000: 0.74, 2125: 0.46, 4182: 0.27, 6846: 0.18, 10000: 0.14
== [5] (a) log-ratio L=2*outer (9s, 6 draws)
Sun ratio p10/50/90/99 [1.01, 1.22, 1.46, 1.74] parent-only v p10/50/90 [19.75, 21.97, 23.67] n 2211
GasGiant ratio p10/50/90/99 [0.94, 1.16, 1.44, 1.85] parent-only v p10/50/90 [16.0, 20.82, 23.9] n 3566
Planet ratio p10/50/90/99 [0.93, 1.07, 1.35, 1.64] parent-only v p10/50/90 [13.71, 15.89, 23.05] n 7445
Planetoid ratio p10/50/90/99 [0.9, 1.04, 1.28, 1.57] parent-only v p10/50/90 [9.91, 12.89, 18.11] n 10200
Belt ratio p10/50/90/99 [0.91, 1.08, 1.35, 1.71] parent-only v p10/50/90 [9.84, 13.14, 19.26] n 3134
planets round suns: parent-only v p10/50/90 [14.61, 16.38, 23.13]; share in 12-18 u/s 0.63
planets round a median sun (mass p40-p60, 10607-20254): v p10/50/90 [15.13, 15.94, 16.98]; share in 12-18 1.00; n 1116
children with no parent pull at their distance: {}
bodies with scale s <= eps (law (a) log-ratio L=2*outer): 7 of 27165
children whose parent pull < 0.5 c/D (v < 0.71 of the log rail): {'Planetoid': 277, 'Belt': 36, 'Planet': 13}
zone far-field c over merged-body c(M_zone), p10/50/90: at 2Z [0.4, 0.52, 0.72], at 10Z [0.02, 0.03, 0.04]
mid-zone inward pull by Z: 1000: 0.85, 2125: 0.54, 4182: 0.32, 6846: 0.22, 10000: 0.16
== [6] (b) screened Plummer lam=0.5*D (16s, 6 draws)
Sun ratio p10/50/90/99 [0.9, 1.16, 1.49, 1.85] parent-only v p10/50/90 [14.3, 17.2, 20.5] n 2211
GasGiant ratio p10/50/90/99 [0.86, 1.04, 1.34, 1.78] parent-only v p10/50/90 [14.3, 17.74, 21.89] n 3566
Planet ratio p10/50/90/99 [0.88, 0.99, 1.19, 1.54] 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.75] parent-only v p10/50/90 [8.7, 12.75, 17.65] n 10200
Belt ratio p10/50/90/99 [0.87, 1.0, 1.22, 1.85] parent-only v p10/50/90 [8.43, 12.9, 18.01] n 3134
planets round suns: parent-only v p10/50/90 [12.61, 16.37, 22.22]; share in 12-18 u/s 0.63
planets round a median sun (mass p40-p60, 10607-20254): v p10/50/90 [14.71, 16.3, 17.03]; share in 12-18 0.98; n 1116
children with no parent pull at their distance: {}
bodies with scale s <= eps (law (b) screened Plummer lam=0.5*D): 285 of 27165
children whose parent pull < 0.5 c/D (v < 0.71 of the log rail): {'Planetoid': 1326, 'Planet': 530, 'Belt': 420, 'GasGiant': 442, 'Sun': 892}
zone far-field c over merged-body c(M_zone), p10/50/90: at 2Z [0.02, 0.07, 0.08], at 10Z [0.0, 0.0, 0.0]
mid-zone inward pull by Z: 1000: 0.41, 2125: 0.29, 4182: 0.18, 6846: 0.13, 10000: 0.09
== [7] (b) screened Plummer lam=1*D (15s, 6 draws)
Sun ratio p10/50/90/99 [1.0, 1.28, 1.56, 1.81] parent-only v p10/50/90 [18.68, 20.76, 22.88] n 2211
GasGiant ratio p10/50/90/99 [0.93, 1.15, 1.46, 1.8] parent-only v p10/50/90 [16.95, 19.89, 23.31] n 3566
Planet ratio p10/50/90/99 [0.92, 1.04, 1.29, 1.6] parent-only v p10/50/90 [14.02, 16.57, 22.77] n 7445
Planetoid ratio p10/50/90/99 [0.91, 1.02, 1.22, 1.57] parent-only v p10/50/90 [10.15, 13.32, 18.0] n 10200
Belt ratio p10/50/90/99 [0.91, 1.04, 1.28, 1.7] parent-only v p10/50/90 [10.08, 13.75, 18.94] n 3134
planets round suns: parent-only v p10/50/90 [14.58, 17.04, 22.86]; share in 12-18 u/s 0.61
planets round a median sun (mass p40-p60, 10607-20254): v p10/50/90 [15.57, 16.81, 17.45]; share in 12-18 1.00; n 1116
children with no parent pull at their distance: {}
bodies with scale s <= eps (law (b) screened Plummer lam=1*D): 34 of 27165
children whose parent pull < 0.5 c/D (v < 0.71 of the log rail): {'Planetoid': 265, 'Belt': 21, 'Planet': 11}
zone far-field c over merged-body c(M_zone), p10/50/90: at 2Z [0.22, 0.45, 0.5], at 10Z [0.0, 0.0, 0.0]
mid-zone inward pull by Z: 1000: 0.70, 2125: 0.50, 4182: 0.32, 6846: 0.23, 10000: 0.16
== [8] (c) Lamb-Oseen core, log tail (10s, 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.66] parent-only v p10/50/90 [18.2, 22.2, 24.57] n 3566
Planet ratio p10/50/90/99 [0.97, 1.24, 1.71, 2.15] parent-only v p10/50/90 [15.5, 17.34, 23.59] n 7445
Planetoid ratio p10/50/90/99 [0.95, 1.11, 1.43, 1.76] parent-only v p10/50/90 [11.42, 14.34, 20.11] n 10200
Belt ratio p10/50/90/99 [0.96, 1.18, 1.54, 2.0] parent-only v p10/50/90 [11.39, 14.69, 20.25] n 3134
planets round suns: parent-only v p10/50/90 [16.26, 17.82, 23.64]; share in 12-18 u/s 0.53
planets round a median sun (mass p40-p60, 10607-20254): v p10/50/90 [17.03, 17.44, 17.96]; share in 12-18 0.93; n 1116
children with no parent pull at their distance: {}
bodies with scale s <= eps (law (c) Lamb-Oseen core, log tail): 27165 of 27165
children whose parent pull < 0.5 c/D (v < 0.71 of the log rail): {'Planetoid': 116, 'Belt': 3, 'Planet': 2}
zone far-field c over merged-body c(M_zone), p10/50/90: at 2Z [3.69, 5.54, 7.98], at 10Z [3.69, 5.54, 7.98]
mid-zone inward pull by Z: 1000: 1.59, 2125: 1.27, 4182: 0.91, 6846: 0.67, 10000: 0.51
== [9] (c) Lamb-Oseen + halo L=0.5*D (10s, 6 draws)
Sun ratio p10/50/90/99 [0.89, 1.17, 1.57, 2.02] parent-only v p10/50/90 [12.31, 15.4, 19.64] n 2211
GasGiant ratio p10/50/90/99 [0.85, 1.04, 1.37, 1.89] parent-only v p10/50/90 [12.69, 17.53, 21.7] n 3566
Planet ratio p10/50/90/99 [0.89, 0.99, 1.19, 1.6] parent-only v p10/50/90 [11.57, 16.25, 22.57] n 7445
Planetoid ratio p10/50/90/99 [0.86, 0.99, 1.15, 1.97] parent-only v p10/50/90 [8.38, 13.27, 19.14] n 10199
Belt ratio p10/50/90/99 [0.86, 0.99, 1.22, 2.06] parent-only v p10/50/90 [7.85, 13.19, 19.18] n 3134
planets round suns: parent-only v p10/50/90 [11.79, 16.67, 22.86]; share in 12-18 u/s 0.55
planets round a median sun (mass p40-p60, 10607-20254): v p10/50/90 [14.76, 16.63, 17.49]; share in 12-18 0.97; n 1116
children with no parent pull at their distance: {'Planetoid': 1}
bodies with scale s <= eps (law (c) Lamb-Oseen + halo L=0.5*D): 285 of 27165
children whose parent pull < 0.5 c/D (v < 0.71 of the log rail): {'Planetoid': 1321, 'Planet': 807, 'Belt': 512, 'GasGiant': 881, 'Sun': 1423}
zone far-field c over merged-body c(M_zone), p10/50/90: at 2Z [0.07, 0.1, 0.12], at 10Z [0.0, 0.0, 0.0]
mid-zone inward pull by Z: 1000: 0.35, 2125: 0.24, 4182: 0.16, 6846: 0.11, 10000: 0.08
== [10] (c) Lamb-Oseen + halo L=1*D (11s, 6 draws)
Sun ratio p10/50/90/99 [0.98, 1.26, 1.55, 1.83] parent-only v p10/50/90 [17.77, 20.04, 22.76] n 2211
GasGiant ratio p10/50/90/99 [0.91, 1.12, 1.44, 1.8] 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.26, 1.59] parent-only v p10/50/90 [14.35, 17.0, 23.35] n 7445
Planetoid ratio p10/50/90/99 [0.89, 1.0, 1.19, 1.63] parent-only v p10/50/90 [10.37, 13.93, 19.83] n 10200
Belt ratio p10/50/90/99 [0.9, 1.02, 1.26, 1.77] parent-only v p10/50/90 [10.18, 14.16, 19.85] n 3134
planets round suns: parent-only v p10/50/90 [14.83, 17.51, 23.44]; share in 12-18 u/s 0.56
planets round a median sun (mass p40-p60, 10607-20254): v p10/50/90 [16.49, 17.16, 17.75]; share in 12-18 0.99; n 1116
children with no parent pull at their distance: {}
bodies with scale s <= eps (law (c) Lamb-Oseen + halo L=1*D): 34 of 27165
children whose parent pull < 0.5 c/D (v < 0.71 of the log rail): {'Planetoid': 131, 'Belt': 8, 'Planet': 7}
zone far-field c over merged-body c(M_zone), p10/50/90: at 2Z [0.25, 0.36, 0.4], at 10Z [0.01, 0.02, 0.02]
mid-zone inward pull by Z: 1000: 0.66, 2125: 0.46, 4182: 0.29, 6846: 0.21, 10000: 0.15
== [11] (c) Lamb-Oseen + halo L=1.5*outer (10s, 6 draws)
Sun ratio p10/50/90/99 [0.98, 1.17, 1.41, 1.71] parent-only v p10/50/90 [18.44, 20.57, 23.1] n 2211
GasGiant ratio p10/50/90/99 [0.93, 1.11, 1.39, 1.86] parent-only v p10/50/90 [15.63, 19.92, 23.6] n 3566
Planet ratio p10/50/90/99 [0.92, 1.04, 1.31, 1.63] parent-only v p10/50/90 [13.32, 15.95, 23.4] n 7445
Planetoid ratio p10/50/90/99 [0.88, 1.02, 1.26, 1.6] parent-only v p10/50/90 [9.56, 13.01, 19.9] n 10200
Belt ratio p10/50/90/99 [0.88, 1.05, 1.32, 1.72] parent-only v p10/50/90 [9.47, 12.88, 19.93] n 3134
planets round suns: parent-only v p10/50/90 [13.86, 16.58, 23.48]; share in 12-18 u/s 0.60
planets round a median sun (mass p40-p60, 10607-20254): v p10/50/90 [14.29, 15.83, 17.22]; share in 12-18 1.00; n 1116
children with no parent pull at their distance: {}
bodies with scale s <= eps (law (c) Lamb-Oseen + halo L=1.5*outer): 7 of 27165
children whose parent pull < 0.5 c/D (v < 0.71 of the log rail): {'Planetoid': 129, 'Belt': 5, 'Planet': 2}
zone far-field c over merged-body c(M_zone), p10/50/90: at 2Z [0.26, 0.34, 0.48], at 10Z [0.01, 0.02, 0.02]
mid-zone inward pull by Z: 1000: 0.75, 2125: 0.46, 4182: 0.28, 6846: 0.18, 10000: 0.14
== [12] (c) Lamb-Oseen + halo L=2*outer (10s, 6 draws)
Sun ratio p10/50/90/99 [1.01, 1.22, 1.46, 1.74] parent-only v p10/50/90 [19.8, 21.98, 23.7] n 2211
GasGiant ratio p10/50/90/99 [0.94, 1.15, 1.44, 1.9] parent-only v p10/50/90 [16.76, 20.91, 23.99] n 3566
Planet ratio p10/50/90/99 [0.92, 1.06, 1.34, 1.65] parent-only v p10/50/90 [14.25, 16.36, 23.46] n 7445
Planetoid ratio p10/50/90/99 [0.9, 1.04, 1.27, 1.58] parent-only v p10/50/90 [10.27, 13.51, 19.99] n 10200
Belt ratio p10/50/90/99 [0.9, 1.07, 1.34, 1.73] parent-only v p10/50/90 [10.18, 13.59, 20.05] n 3134
planets round suns: parent-only v p10/50/90 [14.88, 16.96, 23.53]; share in 12-18 u/s 0.59
planets round a median sun (mass p40-p60, 10607-20254): v p10/50/90 [15.36, 16.4, 17.45]; share in 12-18 1.00; n 1116
children with no parent pull at their distance: {}
bodies with scale s <= eps (law (c) Lamb-Oseen + halo L=2*outer): 7 of 27165
children whose parent pull < 0.5 c/D (v < 0.71 of the log rail): {'Planetoid': 117, 'Belt': 3, 'Planet': 2}
zone far-field c over merged-body c(M_zone), p10/50/90: at 2Z [0.4, 0.52, 0.73], at 10Z [0.02, 0.03, 0.04]
mid-zone inward pull by Z: 1000: 0.87, 2125: 0.54, 4182: 0.32, 6846: 0.22, 10000: 0.16
== [13] (e) conserving halo L=1.5*ext (9s, 6 draws)
Sun ratio p10/50/90/99 [1.03, 1.25, 1.5, 1.73] parent-only v p10/50/90 [21.66, 23.43, 24.48] n 2211
GasGiant ratio p10/50/90/99 [0.96, 1.17, 1.44, 1.88] parent-only v p10/50/90 [16.23, 21.96, 24.4] n 3566
Planet ratio p10/50/90/99 [0.93, 1.06, 1.33, 1.63] parent-only v p10/50/90 [13.42, 16.32, 23.18] n 7445
Planetoid ratio p10/50/90/99 [0.9, 1.04, 1.28, 1.6] parent-only v p10/50/90 [9.36, 12.54, 18.14] n 10200
Belt ratio p10/50/90/99 [0.91, 1.07, 1.34, 1.75] parent-only v p10/50/90 [9.3, 12.64, 19.43] n 3134
planets round suns: parent-only v p10/50/90 [14.29, 16.45, 23.32]; share in 12-18 u/s 0.63
planets round a median sun (mass p40-p60, 10607-20254): v p10/50/90 [14.52, 15.86, 16.92]; share in 12-18 1.00; n 1116
children with no parent pull at their distance: {}
bodies with scale s <= eps (law (e) conserving halo L=1.5*ext): 7 of 27165
children whose parent pull < 0.5 c/D (v < 0.71 of the log rail): {'Planetoid': 356, 'Belt': 54, 'Planet': 20, 'GasGiant': 1}
zone far-field c over merged-body c(M_zone), p10/50/90: at 2Z [1.18, 1.42, 1.99], at 10Z [1.03, 1.28, 1.84]
mid-zone inward pull by Z: 1000: 0.86, 2125: 0.58, 4182: 0.37, 6846: 0.25, 10000: 0.19
== [14] (e) conserving halo L=2*ext (9s, 6 draws)
Sun ratio p10/50/90/99 [1.06, 1.3, 1.55, 1.79] parent-only v p10/50/90 [21.82, 23.49, 24.53] n 2211
GasGiant ratio p10/50/90/99 [0.97, 1.2, 1.49, 1.92] parent-only v p10/50/90 [16.79, 22.01, 24.44] n 3566
Planet ratio p10/50/90/99 [0.94, 1.09, 1.37, 1.66] parent-only v p10/50/90 [14.07, 16.43, 23.2] n 7445
Planetoid ratio p10/50/90/99 [0.91, 1.05, 1.3, 1.58] parent-only v p10/50/90 [10.01, 13.0, 18.14] n 10200
Belt ratio p10/50/90/99 [0.92, 1.09, 1.36, 1.75] parent-only v p10/50/90 [9.97, 13.26, 19.44] n 3134
planets round suns: parent-only v p10/50/90 [14.91, 16.67, 23.32]; share in 12-18 u/s 0.62
planets round a median sun (mass p40-p60, 10607-20254): v p10/50/90 [15.32, 16.31, 17.14]; share in 12-18 1.00; n 1116
children with no parent pull at their distance: {}
bodies with scale s <= eps (law (e) conserving halo L=2*ext): 7 of 27165
children whose parent pull < 0.5 c/D (v < 0.71 of the log rail): {'Planetoid': 239, 'Belt': 22, 'Planet': 9}
zone far-field c over merged-body c(M_zone), p10/50/90: at 2Z [1.24, 1.52, 2.07], at 10Z [1.03, 1.28, 1.85]
mid-zone inward pull by Z: 1000: 0.95, 2125: 0.64, 4182: 0.40, 6846: 0.27, 10000: 0.21

150
docs/research/gravity-smooth-wells/probe2.py Normal file
View file

@ -0,0 +1,150 @@
# Smooth-well probe: ../gravity-tail-taming/probe.py with whole-law pulls F(c, eps, r, s) from laws2.py (core inside the
# law, eps = Roche core radius) and a per-body 'outer' scale (outermost direct child distance; falls back to D).
# Same data, metrics, phases (rng seed 7, 6 draws) and output format; law 0 must reproduce the tail-taming baseline.
# Usage: python probe2.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 laws2 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)
W = lambda F: (lambda c, eps, r, s, cf: F(c, eps, r, s)) # laws without a far strength ignore cf
ONE = lambda m, D, Z, c, cp, o, e: np.ones_like(m)
# (name, F(c, eps, r, s), p, scale(m, D, Z, c, c_parent, outer, own family extent) -> per-body s, subtract disc)
LAWS = [
('today: harmonic core + log', W(F_today), 0.2, ONE, False),
('(d) Plummer log, no taming', W(F_plummer), 0.2, ONE, False),
('(a) log-ratio L=0.5*D', W(F_logratio), 0.2, lambda m, D, Z, c, cp, o, e: 0.5 * D, False),
('(a) log-ratio L=1*D', W(F_logratio), 0.2, lambda m, D, Z, c, cp, o, e: 1.0 * D, False),
('(a) log-ratio L=1.5*outer', W(F_logratio), 0.2, lambda m, D, Z, c, cp, o, e: 1.5 * o, False),
('(a) log-ratio L=2*outer', W(F_logratio), 0.2, lambda m, D, Z, c, cp, o, e: 2.0 * o, False),
('(b) screened Plummer lam=0.5*D', W(F_splummer), 0.2, lambda m, D, Z, c, cp, o, e: 0.5 * D, False),
('(b) screened Plummer lam=1*D', W(F_splummer), 0.2, lambda m, D, Z, c, cp, o, e: 1.0 * D, False),
('(c) Lamb-Oseen core, log tail', W(F_lamb), 0.2, ONE, False),
('(c) Lamb-Oseen + halo L=0.5*D', W(F_lamb_halo), 0.2, lambda m, D, Z, c, cp, o, e: 0.5 * D, False),
('(c) Lamb-Oseen + halo L=1*D', W(F_lamb_halo), 0.2, lambda m, D, Z, c, cp, o, e: 1.0 * D, False),
('(c) Lamb-Oseen + halo L=1.5*outer', W(F_lamb_halo), 0.2, lambda m, D, Z, c, cp, o, e: 1.5 * o, False),
('(c) Lamb-Oseen + halo L=2*outer', W(F_lamb_halo), 0.2, lambda m, D, Z, c, cp, o, e: 2.0 * o, False),
('(e) conserving halo L=1.5*ext', F_conserve, 0.2, lambda m, D, Z, c, cp, o, e: 1.5 * e, False),
('(e) conserving halo L=2*ext', F_conserve, 0.2, lambda m, D, Z, c, cp, o, e: 2.0 * e, 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, cf):
"""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, tail(c[None, :], R[None, :], rs, s[None, :], cf[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); farc = {2: [], 10: []}; nLeR = [0, 0]; weak = collections.Counter()
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)
outer = np.array([max([D[j] for j in range(len(par)) if par[j] == k], default=Dz[k]) for k in range(len(par))])
# own family extent: max path length from the body to any descendant; no descendants -> D (Z for a root)
n_ = len(par); plen = np.zeros((n_, n_))
for j in range(n_):
a, acc = j, 0.0
while par[a] >= 0:
acc += D[a]; a = par[a]; plen[a, j] = acc
ext = np.where(plen.max(axis=1) > 0, plen.max(axis=1), Dz)
# c_far: the body's linear mass share of its parent family's compressed strength (root: the whole zone)
mb = np.where(body, m, 0.0)
MF = np.array([mb[fam[par[k]]].sum() if par[k] >= 0 else mb.sum() for k in range(n_)])
cf = np.where(body, strength(np.maximum(MF, 1e-9), p) * mb / np.maximum(MF, 1e-9), 0.0)
s = scale(m, Dz, Z, c, cpar, outer, ext)
nLeR[0] += int(np.sum(body & (s <= R))); nLeR[1] += int(np.sum(body))
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(tail(c[par[k]], R[par[k]], D[k], s[par[k]], cf[par[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.5 * c[par[k]] / D[k]: weak[cls[k]] += 1
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 = []; farc[2].append(0.0); farc[10].append(0.0)
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, cf[bi])
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, cf[bi]).sum(axis=1)
for fz in (2, 10): # far-field strength: mean inward pull x radius at 2Z and 10Z, over the merged body's c
XF = fz * Z * np.stack([np.cos(th), np.sin(th)], 1)
AF = field(pos[bi], c[bi], R[bi], s[bi], tail, XF, cf[bi]).sum(axis=1)
farc[fz][-1] += np.mean(-(AF * XF).sum(1)) / DRAWS / strength(mb.sum(), p)
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(f' bodies with scale s <= eps (law {name}): {nLeR[0]} of {nLeR[1]}')
print(' children whose parent pull < 0.5 c/D (v < 0.71 of the log rail): ' + str(dict(weak)))
print(' zone far-field c over merged-body c(M_zone), p10/50/90: at 2Z ' + str(q(farc[2])) + ', at 10Z ' + str(q(farc[10])))
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()

20
docs/research/gravity-smooth-wells/seams-output.txt Normal file
View file

@ -0,0 +1,20 @@
A&S K0 at x = 2 (phi = -c K0(r/lam); jumps in K0^(n), multiply by -c/lam^n):
d^0: jump -7.58e-9 1/lam^0 true value 0.113894 |jump|/|true| 6.65e-8
d^1: jump -2.48e-7 1/lam^1 true value -0.139866 |jump|/|true| 1.77e-6
d^2: jump 8.43e-7 1/lam^2 true value 0.183827 |jump|/|true| 4.59e-6
d^3: jump -7.68e-5 1/lam^3 true value -0.266746 |jump|/|true| 0.000288
A&S K1 at x = 2 (force = (c/lam) K1(r/lam)):
d^0: jump 1.48e-9 1/lam^0 true value 0.139866 |jump|/|true| 1.06e-8
d^1: jump 1.92e-7 1/lam^1 true value -0.183827 |jump|/|true| 1.04e-6
d^2: jump -5.14e-6 1/lam^2 true value 0.266746 |jump|/|true| 1.93e-5
d^3: jump 4.64e-5 1/lam^3 true value -0.44408 |jump|/|true| 0.000105
A&S E1 at u = 1 (phi = (c/2) Ein(r^2/eps^2), seam at r = eps; jumps in units c/eps^n):
d^0: jump 2.2e-9 1/eps^0 true value 0.3983 |jump|/|true| 5.52e-9
d^1: jump -2.19e-6 1/eps^1 true value 0.632121 |jump|/|true| 3.46e-6
d^2: jump -0.000128 1/eps^2 true value 0.103638 |jump|/|true| 0.00124
d^3: jump -0.00512 1/eps^3 true value -0.943036 |jump|/|true| 0.00543
Today's core at r = R (phi: r^2/2R^2 inside, ln(r/R)+1/2 beyond; units c/R^n):
d^0: jump 0.0 1/R^0 true value 0.5 |jump|/|true| 0.0
d^1: jump 0.0 1/R^1 true value 1.0 |jump|/|true| 0.0
d^2: jump -2.0 1/R^2 true value -1.0 |jump|/|true| 2.0
d^3: jump 2.0 1/R^3 true value 2.0 |jump|/|true| 1.0

35
docs/research/gravity-smooth-wells/seams.py Normal file
View file

@ -0,0 +1,35 @@
# Seam jumps of the piecewise shader forms (A&S K0/K1 at x = 2; A&S E1 at u = 1) and of today's core at R.
# Each branch is evaluated as its own analytic function in mpmath (50 digits) and differentiated at the seam,
# so the jump is the branch mismatch itself, not finite-difference noise.
import mpmath as mp
mp.mp.dps = 50
def i0(x): t = (x / 3.75) ** 2; return 1 + t*(3.5156229 + t*(3.0899424 + t*(1.2067492 + t*(0.2659732 + t*(0.0360768 + t*0.0045813)))))
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_lo(x): y = (x / 2) ** 2; return -mp.log(x / 2) * i0(x) + (-0.57721566 + y*(0.42278420 + y*(0.23069756 + y*(0.03488590 + y*(0.00262698 + y*(0.00010750 + y*0.00000740))))))
def k0_hi(x): z = 2 / x; return mp.exp(-x) / mp.sqrt(x) * (1.25331414 + z*(-0.07832358 + z*(0.02189568 + z*(-0.01062446 + z*(0.00587872 + z*(-0.00251540 + z*0.00053208))))))
def k1_lo(x): y = (x / 2) ** 2; return mp.log(x / 2) * i1(x) + (1 / x) * (1 + y*(0.15443144 + y*(-0.67278579 + y*(-0.18156897 + y*(-0.01919402 + y*(-0.00110404 + y*(-0.00004686)))))))
def k1_hi(x): z = 2 / x; return mp.exp(-x) / mp.sqrt(x) * (1.25331414 + z*(0.23498619 + z*(-0.03655620 + z*(0.01504268 + z*(-0.00780353 + z*(0.00325614 + z*(-0.00068245)))))))
def e1_lo(u): return -mp.log(u) + (-0.57721566 + u*(0.99999193 + u*(-0.24991055 + u*(0.05519968 + u*(-0.00976004 + u*0.00107857)))))
def e1_hi(u):
p = 0.2677737343 + u*(8.6347608925 + u*(18.0590169730 + u*(8.5733287401 + u)))
q = 3.9584969228 + u*(21.0996530827 + u*(25.6329561486 + u*(9.5733223454 + u)))
return mp.exp(-u) / u * p / q
def report(title, lo, hi, true, x0, unit):
print(title)
for n in range(4):
a, b, t = mp.diff(lo, x0, n), mp.diff(hi, x0, n), mp.diff(true, x0, n)
j = b - a
print(f' d^{n}: jump {mp.nstr(j, 3):>10s} {unit}^{n} true value {mp.nstr(t, 6):>10s} |jump|/|true| {mp.nstr(abs(j / t), 3)}')
# phi = -c K0(r/lam): phi^(n) jump = -c/lam^n * (K0_hi^(n) - K0_lo^(n)) at x = 2, i.e. at r = 2 lam
report('A&S K0 at x = 2 (phi = -c K0(r/lam); jumps in K0^(n), multiply by -c/lam^n):', k0_lo, k0_hi, lambda x: mp.besselk(0, x), 2, '1/lam')
report('A&S K1 at x = 2 (force = (c/lam) K1(r/lam)):', k1_lo, k1_hi, lambda x: mp.besselk(1, x), 2, '1/lam')
# Lamb-Oseen phi = (c/2) Ein(r^2/eps^2): seam at r = eps. Differentiate in r (eps = 1, c = 1) -> units c/eps^n
ein_lo = lambda r: (e1_lo(r * r) + mp.log(r * r) + mp.euler) / 2
ein_hi = lambda r: (e1_hi(r * r) + mp.log(r * r) + mp.euler) / 2
ein_true = lambda r: (mp.e1(r * r) + mp.log(r * r) + mp.euler) / 2
report('A&S E1 at u = 1 (phi = (c/2) Ein(r^2/eps^2), seam at r = eps; jumps in units c/eps^n):', ein_lo, ein_hi, ein_true, 1, '1/eps')
# today's core at R (c = 1, R = 1): inside r^2/2, outside ln r + 1/2
report("Today's core at r = R (phi: r^2/2R^2 inside, ln(r/R)+1/2 beyond; units c/R^n):",
lambda r: r * r / 2, lambda r: mp.log(r) + 0.5, lambda r: mp.log(r) + 0.5, 1, '1/R')