diff --git a/CORRECTIONS.md b/CORRECTIONS.md
index 0c545e0..ce676fb 100644
--- a/CORRECTIONS.md
+++ b/CORRECTIONS.md
@@ -175,3 +175,15 @@ Claims made during the exploration that did not survive it, kept so the record s
     not 1.36. Whether the clouds' 1.36 matches the laboratory's 1.3 to 1.4 for a shared reason is not shown. The
     laboratory measures sections of turbulent surfaces (7/3 in 3D for the interface), the clouds a projection of
     supersonic gas, and the two only meet near 4/3.
+
+## 18 October 2026: cosmic homogeneity (issue #17)
+
+21. **"Beyond about 80 to 300 million parsecs (h⁻¹ Mpc, 80 to 300)."** Two errors, and a missing caveat. The 80
+    to 300 h⁻¹ Mpc in Scrimgeour et al. (2012) is the range over which they rule out a fractal with D₂ below
+    2.97, and 300 is the largest scale their survey reached. Their homogeneity scale is 70 to 81 h⁻¹ Mpc: 71 ± 8,
+    70 ± 5, 81 ± 5 and 75 ± 4 at z ≈ 0.2, 0.4, 0.6 and 0.8. And h⁻¹ Mpc is not Mpc: with h = 0.674, 70 to
+    81 h⁻¹ Mpc is 104 to 120 Mpc. The caveat: the scale is a convention, not a fixed length. It depends on the
+    tracer, the epoch and the threshold; a 0.1% threshold instead of 1% puts it at 185 h⁻¹ Mpc for matter today.
+    See [cosmic-homogeneity](cosmic-homogeneity/).
+22. **"The nesting stops, the way a coastline's roughness stops at the size of molecules."** Points the wrong
+    way. A coastline is rough down to a smallest scale; galaxy clustering is clumpy up to a largest one.
diff --git a/README.md b/README.md
index 6ac95ec..e63d631 100644
--- a/README.md
+++ b/README.md
@@ -24,8 +24,11 @@ with Claude (Anthropic), for [Foundation0](https://github.com/Foundation-0/f0).
   that roughness only if they kick it harder than its own momentum: a cold atom in a room-temperature vacuum,
   yes; nothing as warm as, and heavier than, the gas around it ([quantum-paths](quantum-paths/)).
 - **Cosmic structure is fractal only up to a point.** Galaxies cluster in a fractal-like way on small
-  scales, but beyond about 80 to 300 million parsecs (h⁻¹ Mpc, 80 to 300) they are spread almost evenly
-  (Scrimgeour et al., 2012). The nesting stops, the way a coastline's roughness stops at the size of molecules.
+  scales: well inside 5 h⁻¹ Mpc (the surveys' unit, about 1.5 million parsecs), the number of neighbours
+  within a distance r grows as r^1.2 rather than r³. Beyond a radius of about 70 to 80 h⁻¹ Mpc, roughly 100 to
+  120 million parsecs, the galaxies of the WiggleZ survey are spread evenly to within 1% (Scrimgeour et al.,
+  2012). The nesting stops at a largest scale, as a coastline's roughness stops at a smallest one, and where
+  depends on the galaxies and on how even counts as even ([cosmic-homogeneity](cosmic-homogeneity/)).
 - **Turbulence repeats, but not exactly.** Between the size where a flow is stirred and the size where viscosity
   takes over, the same energy flux passes through every scale. In a cascade model it holds to 0.2% over the first
   nine of 21 shells (2.4 decades), and a third-order law that follows from it to 0.1%. The bursts do not scale
diff --git a/cosmic-homogeneity/README.md b/cosmic-homogeneity/README.md
new file mode 100644
index 0000000..7734942
--- /dev/null
+++ b/cosmic-homogeneity/README.md
@@ -0,0 +1,200 @@
+# Cosmic homogeneity: where galaxy clustering stops looking fractal
+
+Status: **exploration, not peer reviewed, no new result claimed.** Issue
+[#17](https://github.com/Foundation-0/fractals/issues/17).
+
+Galaxies cluster: around any galaxy there are more neighbours nearby than an even spread would give. On large
+enough scales the excess fades. This study computes where, from the standard cosmological model, checks it
+against the surveys, and shows how much the answer depends on the definition.
+
+## The measure
+
+Count the galaxies within a distance r of a galaxy, divided by the count an even spread would give:
+N(<r) = 1 + ξ̄(r), with ξ̄ the correlation function averaged over the sphere. For a fractal of dimension D the
+raw count grows as r^D; for an even spread, as r³. The correlation dimension is
+
+```
+D₂(r) = 3 + d ln N / d ln r = 3 + 3 (ξ − ξ̄) / (1 + ξ̄)
+```
+
+Scrimgeour et al. (2012) define the homogeneity scale R_H as the radius beyond which D₂ stays within 1% of 3
+(D₂ = 2.97). A second definition, where N comes within 1% of 1, gives R_H a few h⁻¹ Mpc larger in ΛCDM (in the
+WiggleZ data not always: 78 against 81 at z ≈ 0.6); every number here uses D₂ unless it says N. Lengths are
+comoving and in h⁻¹ Mpc, the surveys' unit (see Units below).
+
+## Small scales: about 1.2, not 2
+
+Galaxy clustering on small scales follows a power law, ξ = (r/r₀)^−γ with r₀ ≈ 5 h⁻¹ Mpc and γ ≈ 1.8
+(Zehavi et al. 2005). Then:
+
+| r (h⁻¹ Mpc) | 0.1 | 1 | 3 | 5 | 10 | 20 |
+|---|---|---|---|---|---|---|
+| D₂ | 1.20 | 1.24 | 1.45 | 1.71 | 2.25 | 2.69 |
+
+Well inside r₀ the neighbour counts grow as r^1.2, a fractal-like dimension of 3 − γ. Beyond about
+1 h⁻¹ Mpc, D₂ never holds still: it climbs steadily toward 3. Between 1 and 30 h⁻¹ Mpc the counts grow on
+average as r^1.9 (endpoint slope; a least-squares fit to log-spaced points gives about 1.86), an average through
+the climb, and one reading of the "D ≈ 2" of Sylos Labini, Pietronero and colleagues, who argue the distribution
+stays fractal to the largest scales. Their 2009 value, D = 2.1 ± 0.1 up to about 30 h⁻¹ Mpc, is about 2σ above
+this average. Their estimator and sample differ, and the power law here is fitted only below about 10 to
+20 h⁻¹ Mpc, so the comparison is rough.
+
+## Large scales: the homogeneity scale
+
+Linear theory: the matter power spectrum of Eisenstein and Hu (1998) with baryon wiggles, normalised to σ₈,
+gives ξ and ξ̄. Galaxies trace matter with a bias b (ξ_g = b²ξ), and redshift-space distortion adds the Kaiser
+factor 1 + 2β/3 + β²/5.
+
+**WiggleZ.** With Scrimgeour et al.'s fiducial cosmology and each slice's measured bias and β, the model
+gives their own ΛCDM homogeneity scales to within about 1 h⁻¹ Mpc. The paper fits the bias separately for each
+statistic (Table 4); the b² here is the fit for D₂:
+
+| z | b² | Model | Their ΛCDM | Measured |
+|---|---|---|---|---|
+| 0.2 | 0.70 | 77.1 | 76 | 71 ± 8 |
+| 0.4 | 0.87 | 78.6 | 78 | 70 ± 5 |
+| 0.6 | 1.25 | 81.5 | 81 | 81 ± 5 |
+| 0.8 | 1.46 | 78.7 | 78 | 75 ± 4 |
+
+The model runs 0.5 to 1.1 above their values: they use a nonlinear spectrum and damp the redshift-space signal
+(see What the model leaves out). R_H barely moves with redshift because the galaxies' bias grows as the
+clustering of matter weakens, as the paper says. The other definition, with the bias fitted to the counts
+(b² = 0.68, 0.87, 1.23, 1.43), gives 80.7, 83.0, 85.3 and 82.5, against their ΛCDM counts of 80, 83, 85 and 82
+(Table 4).
+
+Other surveys use different criteria and tracers. Hogg et al. (2005) found SDSS luminous red galaxy counts
+consistent with r³ beyond about 70 h⁻¹ Mpc, within their errors; for such strongly biased galaxies (b ≈ 2) the
+model's 1% scale is about 130. Yadav et al. (2005) and Sarkar et al. (2009) found 60 to 70 h⁻¹ Mpc from
+multifractal dimensions. Ntelis et al. (2017) used D₂ = 2.97 on BOSS CMASS galaxies, corrected the counts for bias
+and redshift-space distortion, and found 63.3 ± 0.7 h⁻¹ Mpc for the matter distribution, against a ΛCDM 62.0
+(abstract; their Conclusions give 61.9 ± 0.8 for 0.538 < z < 0.592). The 62.0 is a real-space matter value: the
+paper's matter values (Sec. 5.5, Table 5) run from 59 to 67, its galaxy values from 109 to 147. In their
+fiducial cosmology, Planck 2015, the model gives 61.2 for matter at z = 0.57, near the middle of that slice (the
+paper attaches no redshift to 62.0). With Planck 2018 it gives 59.9, mostly because σ₈ is lower (0.811 against 0.831).
+For galaxies of their bias in that slice (b ≈ 1.94) in redshift space the model gives 124; they measure 114 ± 6.
+
+## The scale is a convention
+
+R_H is where clustering falls below a chosen threshold, not a length nature sets.
+
+| How close to even (D₂ = 3(1 − ε)) | ε = 1% | 0.5% | 0.2% | 0.1% |
+|---|---|---|---|---|
+| R_H, matter today | 77 | 95 | 149 | 185 |
+
+| z | Matter | Galaxies, b = 1 | Galaxies, b = 2 |
+|---|---|---|---|
+| 0 | 77 | 87 | 147 |
+| 0.5 | 62 | 76 | 127 |
+| 1 | 50 | 64 | 95 |
+| 2.4 | 29 | 41 | 69 |
+
+(Galaxy columns are in redshift space.) A ten times stricter threshold more than doubles R_H. Matter becomes
+even on smaller scales in the past, because its clustering has had less time to grow. Strongly biased tracers
+look clumpy out to twice the scale. Other definitions give other numbers: Yadav, Bagla and Khandai (2010) put
+an upper limit of about 260 h⁻¹ Mpc on ΛCDM by asking where the deviation from 3 drops below its own
+statistical scatter, and an entropy measure gives about 140 h⁻¹ Mpc for SDSS (Pandey and Sarkar 2015).
+
+The baryon acoustic peak shows up too. The linear ξ has a bump at 99.6 h⁻¹ Mpc (Eisenstein et al. 2005 found
+the peak near 100), and it holds D₂ almost flat between about 100 and 120 h⁻¹ Mpc. In linear theory D₂ even dips
+slightly near 113, and then R_H jumps: for a b = 2 tracer at z = 0.5 in real space (R_H = 99), half a percent
+more clustering amplitude, near ×1.055, moves R_H by 13 h⁻¹ Mpc. Nonlinear motions smooth the peak: fits to BOSS
+galaxies before reconstruction give Σ ≈ 8 h⁻¹ Mpc (Beutler et al. 2017). With Σ = 7 to 10 h⁻¹ Mpc the dip and
+the jump are gone, the peak moves to about 96, and the 1% scale for matter today is unchanged (77).
+
+## Where the README's "80 to 300" came from
+
+From the last sentence of Scrimgeour et al.'s abstract: they exclude a fractal with D₂ below 2.97 "on scales
+from ~80 h⁻¹ Mpc up to the largest scales probed by our measurement, ~300 h⁻¹ Mpc". That is the range they
+could rule a fractal out over, not the scale where evenness starts, and 300 is where their survey ran out. It
+is also not Yadav et al.'s 260.
+
+## Units
+
+h⁻¹ Mpc is a megaparsec divided by h, the Hubble constant in units of 100 km/s/Mpc. With h = 0.674 (Planck
+2018), 70 to 80 h⁻¹ Mpc is 104 to 119 Mpc, about 340 to 390 million light-years. The README read
+"80 to 300 h⁻¹ Mpc" as "80 to 300 million parsecs".
+
+## The coastline
+
+A coastline is rough down to some smallest scale (in the end, grains and molecules); the galaxy distribution
+is clumpy up to some largest scale. Both are fractal only between two limits, but the cosmic limit that
+"stops the nesting" is the upper one.
+
+## Bigger structures
+
+Some structures are claimed to exceed the homogeneity scale: the Huge Large Quasar Group (Clowes et al. 2013,
+about 1.2 Gpc long), the Giant Arc (Lopez et al. 2022, about 1 Gpc) and the Big Ring (Lopez et al. 2024, about
+400 Mpc). Nadathur (2013) showed that the same group finder turns up Gpc "structures" in random points, and
+Sawala et al. (2025) find Giant Arc-like patterns common in ΛCDM simulations; Lopez and Clowes dispute that.
+Homogeneity is a statement about averages, so single large structures do not by themselves contradict it.
+
+## Run it
+
+```
+node cosmic-homogeneity/table.mjs    # D₂(r), R_H against threshold, tracer and redshift, and the WiggleZ comparison
+npm test                             # includes tests/cosmic-homogeneity.mjs
+```
+
+`homogeneity.js` is the whole model. It loads in Node (`require`) and in a browser (`window.HOMOGENEITY`).
+
+## Checks
+
+`tests/cosmic-homogeneity.mjs`:
+
+| Check | Published | Model |
+|---|---|---|
+| WiggleZ ΛCDM R_H, four slices (Scrimgeour et al. 2012) | 76, 78, 81, 78 | 77.1, 78.6, 81.5, 78.7 |
+| Measured WiggleZ R_H within 2σ of the model | 71±8, 70±5, 81±5, 75±4 | yes |
+| WiggleZ ΛCDM R_H by the counts definition (their Table 4) | 80, 83, 85, 82 | 80.7, 83.0, 85.3, 82.5 |
+| Matter R_H at z = 0.57, Planck 2015 (Ntelis et al. 2017, ΛCDM, matter) | 62.0 | 61.2 |
+| Sound horizon at Planck's drag redshift (Planck 2018) | 147.09 Mpc | 147.25 Mpc |
+| Baryon acoustic peak in ξ, sanity check (Eisenstein et al. 2005) | ~100 | 99.6 |
+| D₂ deep inside r₀ (model; γ = 1.8 from Zehavi et al. 2005) | 3 − γ = 1.2 | 1.20 |
+| Endpoint count slope over 1–30 h⁻¹ Mpc (model) | between 1.2 and 3 | 1.90 |
+| Threshold and redshift move R_H (model) | | 77 → 185 for 1% → 0.1%; 77, 62, 50, 29 at z = 0, 0.5, 1, 2.4 |
+| Acoustic plateau, b = 2, z = 0.5: largest R_H change for a 0.5% amplitude step (model) | | 12.8 (linear), 0.8 (smoothed) |
+
+## What the model leaves out
+
+- Linear theory: no nonlinear clustering, which matters below about 20 h⁻¹ Mpc, apart from the optional
+  smoothing of the acoustic peak, and no velocity damping of the redshift-space signal. Scrimgeour et al.'s
+  ΛCDM uses a nonlinear spectrum (CAMB with HALOFIT) and a Lorentzian damping term; with linear Eisenstein–Hu
+  and the Kaiser factor alone the model comes out 0.5 to 1.1 h⁻¹ Mpc above their D₂ values and 0 to 0.7 above
+  their counts values.
+- The Eisenstein–Hu equations, constants and equation numbers were checked against the paper's arXiv LaTeX
+  source (astro-ph/9709112v1) in review on #26: an independent transcription gives the same T(k) to within
+  7 × 10⁻¹⁶. Hu's own code (tf_fit.c) differs slightly from the paper text, moving T by at most 0.06% and R_H by
+  at most 0.02 h⁻¹ Mpc; this model follows the paper. Their fitted drag redshift (1021) is lower than Planck's
+  (1060), by construction (the paper reduces z_d by a factor of 0.96), so their sound horizon is 2.6% larger;
+  at Planck's drag redshift the same formula gives 147.25 Mpc against Planck's 147.09.
+- The WiggleZ bias, β and R_H values and the Ntelis et al. values were checked against the papers' arXiv
+  LaTeX sources in review on #26.
+- The galaxy power law is a fit below about 10 to 20 h⁻¹ Mpc: r₀ and γ change with galaxy luminosity, and γ
+  steepens for the brightest. The 1-to-30 average extrapolates it.
+- Neutrinos are left out: Planck's Ωm includes about 0.06 eV of neutrino mass, which the Eisenstein–Hu fit
+  does not model.
+- Clustering lengths are comoving and in h⁻¹ Mpc. The sound horizon and the large-structure sizes are in Mpc,
+  as their sources give them. The grid resolves ξ out to 300 h⁻¹ Mpc, so R_H is not searched beyond that.
+
+## Sources
+
+- M. I. Scrimgeour et al., *The WiggleZ Dark Energy Survey: the transition to large-scale cosmic homogeneity*, MNRAS 425, 116 (2012). https://arxiv.org/abs/1205.6812
+- D. J. Eisenstein and W. Hu, *Baryonic features in the matter transfer function*, ApJ 496, 605 (1998). https://arxiv.org/abs/astro-ph/9709112
+- Planck Collaboration, *Planck 2018 results. VI. Cosmological parameters*, A&A 641, A6 (2020). https://arxiv.org/abs/1807.06209
+- Planck Collaboration, *Planck 2015 results. XIII. Cosmological parameters*, A&A 594, A13 (2016). https://arxiv.org/abs/1502.01589
+- D. W. Hogg et al., *Cosmic homogeneity demonstrated with luminous red galaxies*, ApJ 624, 54 (2005). https://arxiv.org/abs/astro-ph/0411197
+- J. Yadav, S. Bharadwaj, B. Pandey and T. R. Seshadri, MNRAS 364, 601 (2005). https://arxiv.org/abs/astro-ph/0504315
+- P. Sarkar, J. Yadav, B. Pandey and S. Bharadwaj, MNRAS 399, L128 (2009). https://arxiv.org/abs/0906.3431
+- J. K. Yadav, J. S. Bagla and N. Khandai, *Fractal dimension as a measure of the scale of homogeneity*, MNRAS 405, 2009 (2010). https://arxiv.org/abs/1001.0617
+- P. Ntelis et al., *Exploring cosmic homogeneity with the BOSS DR12 galaxy sample*, JCAP 06 (2017) 019. https://arxiv.org/abs/1702.02159
+- B. Pandey and S. Sarkar, MNRAS 454, 2647 (2015). https://arxiv.org/abs/1507.03124
+- I. Zehavi et al., *The luminosity and color dependence of the galaxy correlation function*, ApJ 630, 1 (2005). https://arxiv.org/abs/astro-ph/0408569
+- D. J. Eisenstein et al., *Detection of the baryon acoustic peak in the large-scale correlation function of SDSS luminous red galaxies*, ApJ 633, 560 (2005). https://arxiv.org/abs/astro-ph/0501171
+- F. Sylos Labini, M. Montuori and L. Pietronero, *Scale-invariance of galaxy clustering*, Phys. Rep. 293, 61 (1998). https://arxiv.org/abs/astro-ph/9711073
+- F. Sylos Labini, N. L. Vasilyev, L. Pietronero and Y. V. Baryshev, EPL 86, 49001 (2009). https://arxiv.org/abs/0805.1132
+- R. G. Clowes et al., MNRAS 429, 2910 (2013). https://arxiv.org/abs/1211.6256
+- S. Nadathur, *Seeing patterns in noise: gigaparsec-scale 'structures' that do not violate homogeneity*, MNRAS 434, 398 (2013). https://arxiv.org/abs/1306.1700
+- F. Beutler et al., *The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: baryon acoustic oscillations in the Fourier space*, MNRAS 464, 3409 (2017). https://arxiv.org/abs/1607.03149
+- A. M. Lopez, R. G. Clowes and G. M. Williger, MNRAS 516, 1557 (2022). https://arxiv.org/abs/2201.06875
+- A. M. Lopez et al., *A Big Ring on the sky*, JCAP (2024). https://arxiv.org/abs/2402.07591
+- T. Sawala et al., MNRAS Letters 541, L22 (2025). https://arxiv.org/abs/2502.03515
diff --git a/cosmic-homogeneity/homogeneity.js b/cosmic-homogeneity/homogeneity.js
new file mode 100644
index 0000000..a55b985
--- /dev/null
+++ b/cosmic-homogeneity/homogeneity.js
@@ -0,0 +1,162 @@
+/* Cosmic homogeneity model: the scale where the clustering of matter and galaxies stops looking fractal. No dependencies.
+ *
+ * Count the galaxies within a distance r of a galaxy and divide by the count for an even spread: N(<r) = 1 + xibar(r),
+ * where xibar is the correlation function averaged over the sphere. A fractal of dimension D gives counts growing as
+ * r^D; an even spread, as r^3. The correlation dimension is
+ *   D2(r) = 3 + dln N / dln r = 3 + 3 (xi - xibar) / (1 + xibar)      (using dxibar/dln r = 3 (xi - xibar)).
+ * The homogeneity scale R_H is where D2 comes within 1% of 3 (D2 = 2.97), or N within 1% of 1 (Scrimgeour et al. 2012).
+ * Linear theory: xi(r) = int k^2 P(k) j0(kr) dk / 2pi^2; xibar is the same with the top-hat window 3 j1(kr)/(kr).
+ * P(k) = A k^ns T(k)^2 with the Eisenstein & Hu (1998) transfer function, normalised to sigma8.
+ * Galaxies: xi_g = b^2 D(z)^2 xi_0 in real space, times the Kaiser factor 1 + 2beta/3 + beta^2/5 in redshift space.
+ * Optionally the acoustic peak is smoothed as nonlinear motions do: P = P_nw + (P - P_nw) exp(-k^2 Sigma^2 / 2), with
+ * P_nw the Eisenstein-Hu no-wiggle spectrum.
+ * Lengths in h^-1 Mpc and wavenumbers in h Mpc^-1, except where a name says Mpc.
+ */
+(function (root) {
+  'use strict';
+  // Planck 2018, TT,TE,EE+lowE+lensing (Planck Collaboration VI, A&A 641, A6 (2020)), with its drag redshift and
+  // sound horizon at the drag epoch (Mpc).
+  const PLANCK18 = { h: 0.6736, Om: 0.3153, Ob: 0.02237 / (0.6736 * 0.6736), ns: 0.9649, sigma8: 0.8111, Tcmb: 2.7255, zdrag: 1059.94, rdrag: 147.09 };
+  // Planck 2015, TT,TE,EE+lowP (Planck Collaboration XIII, A&A 594, A13 (2016), Table 4): the fiducial cosmology of
+  // Ntelis et al. (2017).
+  const PLANCK15 = { h: 0.6727, Om: 0.3156, Ob: 0.02225 / (0.6727 * 0.6727), ns: 0.9645, sigma8: 0.831, Tcmb: 2.7255 };
+  // The fiducial cosmology of Scrimgeour et al. (2012), as in Blake et al. (2011).
+  const SCRIMGEOUR12 = { h: 0.71, Om: 0.27, Ob: 0.04482, ns: 0.96, sigma8: 0.8, Tcmb: 2.725 };
+  // WiggleZ redshift slices (Scrimgeour et al. 2012). Table 4 fits the bias b^2 separately for each statistic: b2D in the
+  // D2 block, b2N in the N(<r) block. beta from Table 3. Measured R_H from D2 = 2.97 (abstract), and the paper's LCDM
+  // values with best-fitting bias by D2 (RHlcdm) and by N (RHlcdmN) from Table 4, h^-1 Mpc.
+  const WIGGLEZ = [
+    { z: 0.2, b2D: 0.70, b2N: 0.68, beta: 0.69, RH: 71, err: 8, RHlcdm: 76, RHlcdmN: 80 },
+    { z: 0.4, b2D: 0.87, b2N: 0.87, beta: 0.73, RH: 70, err: 5, RHlcdm: 78, RHlcdmN: 83 },
+    { z: 0.6, b2D: 1.25, b2N: 1.23, beta: 0.60, RH: 81, err: 5, RHlcdm: 81, RHlcdmN: 85 },
+    { z: 0.8, b2D: 1.46, b2N: 1.43, beta: 0.51, RH: 75, err: 4, RHlcdm: 78, RHlcdmN: 82 },
+  ];
+
+  // Eisenstein & Hu 1998, ApJ 496, 605: matter transfer function with baryon wiggles (eqs. 2-24).
+  function transferEH(par) {
+    const th = par.Tcmb / 2.7, om = par.Om * par.h * par.h, ob = par.Ob * par.h * par.h;
+    const fb = par.Ob / par.Om, fc = 1 - fb;
+    const zeq = 2.5e4 * om * Math.pow(th, -4);                       // eq. 2
+    const keq = 7.46e-2 * om * Math.pow(th, -2);                     // eq. 3, Mpc^-1
+    const b1 = 0.313 * Math.pow(om, -0.419) * (1 + 0.607 * Math.pow(om, 0.674)), b2 = 0.238 * Math.pow(om, 0.223);
+    const zd = 1291 * Math.pow(om, 0.251) / (1 + 0.659 * Math.pow(om, 0.828)) * (1 + b1 * Math.pow(ob, b2)); // eq. 4
+    const R = z => 31.5 * ob * Math.pow(th, -4) * (1e3 / z);         // eq. 5
+    const Req = R(zeq), Rd = R(zd);
+    const soundHorizon = z => { const R1 = R(z); return 2 / (3 * keq) * Math.sqrt(6 / Req) * Math.log((Math.sqrt(1 + R1) + Math.sqrt(R1 + Req)) / (1 + Math.sqrt(Req))); };
+    const s = soundHorizon(zd); // eq. 6, Mpc, at Eisenstein and Hu's fitted drag redshift
+    const ksilk = 1.6 * Math.pow(ob, 0.52) * Math.pow(om, 0.73) * (1 + Math.pow(10.4 * om, -0.95)); // eq. 7, Mpc^-1
+    const a1 = Math.pow(46.9 * om, 0.670) * (1 + Math.pow(32.1 * om, -0.532));
+    const a2 = Math.pow(12.0 * om, 0.424) * (1 + Math.pow(45.0 * om, -0.582));
+    const ac = Math.pow(a1, -fb) * Math.pow(a2, -fb * fb * fb);      // eq. 11
+    const c1 = 0.944 / (1 + Math.pow(458 * om, -0.708)), c2 = Math.pow(0.395 * om, -0.0266);
+    const bc = 1 / (1 + c1 * (Math.pow(fc, c2) - 1));                // eq. 12
+    const y = (1 + zeq) / (1 + zd), sy = Math.sqrt(1 + y);
+    const G = y * (-6 * sy + (2 + 3 * y) * Math.log((sy + 1) / (sy - 1))); // eq. 15
+    const ab = 2.07 * keq * s * Math.pow(1 + Rd, -0.75) * G;         // eq. 14
+    const bnode = 8.41 * Math.pow(om, 0.435);                         // eq. 23
+    const bb = 0.5 + fb + (3 - 2 * fb) * Math.sqrt(Math.pow(17.2 * om, 2) + 1); // eq. 24
+    const T0 = (k, a, b) => { // eqs. 10, 19, 20
+      const q = k / (13.41 * keq), L = Math.log(Math.E + 1.8 * b * q), C = 14.2 / a + 386 / (1 + 69.9 * Math.pow(q, 1.08));
+      return L / (L + C * q * q);
+    };
+    const j0 = x => (x < 1e-4 ? 1 - x * x / 6 : Math.sin(x) / x);
+    function T(kh) {
+      const k = kh * par.h, ks = k * s;
+      const f = 1 / (1 + Math.pow(ks / 5.4, 4));                      // eq. 18
+      const Tc = f * T0(k, 1, bc) + (1 - f) * T0(k, ac, bc);          // eq. 17
+      const st = s / Math.cbrt(1 + Math.pow(bnode / ks, 3));          // eq. 22
+      const Tb = (T0(k, 1, 1) / (1 + Math.pow(ks / 5.2, 2)) + ab / (1 + Math.pow(bb / ks, 3)) * Math.exp(-Math.pow(k / ksilk, 1.4))) * j0(k * st); // eq. 21
+      return fb * Tb + fc * Tc;                                       // eq. 16
+    }
+    // No-wiggle form (eqs. 26, 28-31): the same broad shape without the acoustic oscillations.
+    const sNw = 44.5 * Math.log(9.83 / om) / Math.sqrt(1 + 10 * Math.pow(ob, 0.75)); // eq. 26, Mpc
+    const aG = 1 - 0.328 * Math.log(431 * om) * fb + 0.38 * Math.log(22.3 * om) * fb * fb; // eq. 31
+    function Tnw(kh) {
+      const k = kh * par.h, G = par.Om * par.h * (aG + (1 - aG) / (1 + Math.pow(0.43 * k * sNw, 4))); // eq. 30
+      const q = kh * th * th / G, L = Math.log(2 * Math.E + 1.8 * q), C = 14.2 + 731 / (1 + 62.5 * q); // eqs. 28, 29
+      return L / (L + C * q * q);
+    }
+    return { T, Tnw, s, soundHorizon, zd, zeq };
+  }
+
+  // Linear growth in flat LCDM: D(a) proportional to H(a) int_0^a da' / (a' H(a'))^3, with D = 1 today; f = dln D / dln a.
+  function growth(par) {
+    const E = a => Math.sqrt(par.Om / (a * a * a) + 1 - par.Om);
+    const Dun = a => { const N = 2000; let s = 0; for (let i = 0; i < N; i++) { const x = a * (i + 0.5) / N; s += a / N / Math.pow(x * E(x), 3); } return E(a) * s; };
+    const D0 = Dun(1);
+    const D = z => Dun(1 / (1 + z)) / D0;
+    const f = z => { const a = 1 / (1 + z), e = 1e-4; return Math.log(Dun(a * (1 + e)) / Dun(a * (1 - e))) / Math.log((1 + e) / (1 - e)); };
+    return { D, f };
+  }
+
+  const tophat = x => (x < 1e-3 ? 1 - x * x / 10 : 3 * (Math.sin(x) - x * Math.cos(x)) / (x * x * x));
+  const sinc = x => (x < 1e-4 ? 1 - x * x / 6 : Math.sin(x) / x);
+
+  // A cosmology: today's linear matter power spectrum and the correlation integrals on a log-k grid.
+  // opts.sigmaNL (h^-1 Mpc) smooths the acoustic peak. A Gaussian cutoff exp(-(k kd)^2), kd = 0.25 h^-1 Mpc, lets the
+  // oscillating integrals converge at small r; it changes D2 by under 1e-5 beyond 40 h^-1 Mpc. The log-k grid
+  // resolves xi to 0.3% out to r = 300 h^-1 Mpc and not beyond, so R_H is searched no further than that.
+  function cosmology(par = PLANCK18, opts = {}) {
+    const tf = transferEH(par), gr = growth(par), SNL = opts.sigmaNL || 0;
+    const NK = opts.nk || 8000, k0 = 1e-5, k1 = 100, dl = Math.log(k1 / k0) / NK, kd = opts.kdamp || 0.25, K = [];
+    for (let i = 0; i < NK; i++) K.push(k0 * Math.exp((i + 0.5) * dl));
+    const P0 = k => Math.pow(k, par.ns) * Math.pow(tf.T(k), 2);
+    const Psm = k => { const pnw = Math.pow(k, par.ns) * Math.pow(tf.Tnw(k), 2); return pnw + (P0(k) - pnw) * Math.exp(-k * k * SNL * SNL / 2); };
+    const sum = (w, f) => { let s = 0; for (let i = 0; i < NK; i++) s += w[i] * f(K[i]); return s; };
+    const w0 = K.map(k => k * k * k * P0(k) * dl / (2 * Math.PI * Math.PI));
+    const A = par.sigma8 * par.sigma8 / sum(w0, k => Math.pow(tophat(8 * k), 2));
+    const w = w0.map(x => A * x);
+    const wd = K.map((k, i) => (SNL ? A * k * k * k * Psm(k) * dl / (2 * Math.PI * Math.PI) : w[i]) * Math.exp(-k * k * kd * kd));
+    return {
+      par, tf, growth: gr,
+      P: k => A * P0(k),
+      sigma: R => Math.sqrt(sum(w, k => Math.pow(tophat(k * R), 2))),
+      xi: r => sum(wd, k => sinc(k * r)),
+      xibar: r => sum(wd, k => tophat(k * r)),
+    };
+  }
+
+  // Clustering amplitude relative to today's matter: b^2 D(z)^2, times the Kaiser factor in redshift space.
+  // beta defaults to f(z)/b.
+  function amplitude(cos, { z = 0, b2 = 1, beta, redshiftSpace = false } = {}) {
+    const D = cos.growth.D(z), bt = beta === undefined ? cos.growth.f(z) / Math.sqrt(b2) : beta;
+    return b2 * D * D * (redshiftSpace ? 1 + 2 * bt / 3 + bt * bt / 5 : 1);
+  }
+  function D2(cos, r, amp = 1) {
+    const x = amp * cos.xi(r), xb = amp * cos.xibar(r);
+    return 3 + 3 * (x - xb) / (1 + xb);
+  }
+  const counts = (cos, r, amp = 1) => 1 + amp * cos.xibar(r);
+  // Smallest r beyond which the clustering stays within eps of homogeneity: D2 > 3(1 - eps) (by = 'D2'), or
+  // N < 1 + eps (by = 'N'). Scans outward from 1 to rmax h^-1 Mpc and bisects the last crossing; NaN if the
+  // clustering is already that even at 1 or not yet at rmax.
+  function homogeneityScale(cos, amp = 1, { eps = 0.01, by = 'D2', rmax = 300 } = {}) {
+    const inside = by === 'N' ? r => counts(cos, r, amp) > 1 + eps : r => D2(cos, r, amp) < 3 * (1 - eps);
+    if (!inside(1)) return NaN;
+    let last = null, prev = true;
+    for (let r = 1; r <= rmax; r *= 1.02) { const now = inside(r); if (!now && prev) last = r; prev = now; }
+    if (last === null || prev) return NaN;
+    let lo = last / 1.02, hi = last;
+    for (let i = 0; i < 40; i++) { const mid = Math.sqrt(lo * hi); if (inside(mid)) lo = mid; else hi = mid; }
+    return hi;
+  }
+  // The baryon acoustic peak: the local maximum of xi between 80 and 130 h^-1 Mpc (NaN if xi has none there).
+  function baoPeak(cos) {
+    let best = 80, bv = -Infinity;
+    for (let r = 80; r <= 130; r += 0.1) { const v = cos.xi(r); if (v > bv) { bv = v; best = r; } }
+    return best < 80.05 || best > 129.95 ? NaN : best;
+  }
+  // Galaxies on small scales: xi = (r/r0)^-gamma. Deep inside r0 the counts grow as r^(3 - gamma).
+  function D2powerLaw(r, r0, gamma) {
+    const x = Math.pow(r / r0, -gamma), xb = 3 * x / (3 - gamma);
+    return 3 + 3 * (x - xb) / (1 + xb);
+  }
+  // Endpoint (secant) slope of ln N(<r) (unscaled counts, n r^3 (1 + xibar)) between r1 and r2.
+  function slopePowerLaw(r1, r2, r0, gamma) {
+    const lnN = r => 3 * Math.log(r) + Math.log(1 + 3 * Math.pow(r / r0, -gamma) / (3 - gamma));
+    return (lnN(r2) - lnN(r1)) / Math.log(r2 / r1);
+  }
+
+  const api = { PLANCK18, PLANCK15, SCRIMGEOUR12, WIGGLEZ, transferEH, growth, cosmology, amplitude, D2, counts, homogeneityScale, baoPeak, D2powerLaw, slopePowerLaw };
+  if (typeof module !== 'undefined' && module.exports) module.exports = api; else root.HOMOGENEITY = api;
+})(typeof window !== 'undefined' ? window : globalThis);
diff --git a/cosmic-homogeneity/table.mjs b/cosmic-homogeneity/table.mjs
new file mode 100644
index 0000000..97216dd
--- /dev/null
+++ b/cosmic-homogeneity/table.mjs
@@ -0,0 +1,49 @@
+// Prints the correlation dimension D2(r) and the homogeneity scale for matter and galaxies.
+// Usage: node cosmic-homogeneity/table.mjs
+import { createRequire } from 'node:module';
+const H = createRequire(import.meta.url)('./homogeneity.js');
+
+const pc = H.cosmology(H.PLANCK18);
+console.log('Correlation dimension D2(r), Planck 2018, h^-1 Mpc\n');
+console.log('r'.padStart(6) + 'galaxies, power law'.padStart(22) + 'matter z = 0, linear'.padStart(23));
+for (const r of [0.1, 0.3, 1, 3, 5, 10, 20, 40, 60, 80, 100, 120, 150, 200, 300]) {
+  const lin = r >= 10 ? H.D2(pc, r).toFixed(4) : '-', pl = r <= 20 ? H.D2powerLaw(r, 5, 1.8).toFixed(3) : '-';
+  console.log(String(r).padStart(6) + pl.padStart(22) + lin.padStart(23));
+}
+console.log('\nGalaxy power law xi = (r/5)^-1.8 (Zehavi et al. 2005) up to 20 h^-1 Mpc; linear theory from 10 h^-1 Mpc out.');
+
+console.log('\nHomogeneity scale R_H, h^-1 Mpc: where D2 reaches 3 (1 - eps)\n');
+console.log('eps'.padStart(8) + 'matter z = 0'.padStart(14));
+for (const eps of [0.01, 0.005, 0.002, 0.001]) console.log(String(eps).padStart(8) + H.homogeneityScale(pc, 1, { eps }).toFixed(0).padStart(14));
+
+console.log('\nR_H (eps = 1%) for matter and for galaxies of bias b, h^-1 Mpc\n');
+console.log('z'.padStart(5) + 'matter'.padStart(9) + 'b = 1, z-space'.padStart(16) + 'b = 2, z-space'.padStart(16));
+for (const z of [0, 0.5, 1, 2.4]) {
+  const a = H.amplitude(pc, { z }), a1 = H.amplitude(pc, { z, redshiftSpace: true }), a2 = H.amplitude(pc, { z, b2: 4, redshiftSpace: true });
+  console.log(String(z).padStart(5) + H.homogeneityScale(pc, a).toFixed(0).padStart(9) + H.homogeneityScale(pc, a1).toFixed(0).padStart(16) + H.homogeneityScale(pc, a2).toFixed(0).padStart(16));
+}
+
+console.log('\nWiggleZ (Scrimgeour et al. 2012): their cosmology and measured bias, redshift space, h^-1 Mpc');
+console.log('D2: where D2 = 2.97 (the paper\'s definition); N: where the scaled counts N = 1.01');
+console.log('Each statistic uses its own fitted bias (Table 4): b^2 D2 for D2, b^2 N for N\n');
+console.log('z'.padStart(5) + 'b^2 D2'.padStart(8) + 'model D2'.padStart(10) + 'paper LCDM'.padStart(12) + 'measured'.padStart(11) +
+  'b^2 N'.padStart(7) + 'model N'.padStart(9) + 'paper LCDM'.padStart(12));
+const sc = H.cosmology(H.SCRIMGEOUR12);
+for (const s of H.WIGGLEZ) {
+  const amp = b2 => H.amplitude(sc, { z: s.z, b2, beta: s.beta, redshiftSpace: true });
+  const m = H.homogeneityScale(sc, amp(s.b2D)), n = H.homogeneityScale(sc, amp(s.b2N), { by: 'N' });
+  console.log(String(s.z).padStart(5) + s.b2D.toFixed(2).padStart(8) + m.toFixed(1).padStart(10) + String(s.RHlcdm).padStart(12) + `${s.RH}±${s.err}`.padStart(11) +
+    s.b2N.toFixed(2).padStart(7) + n.toFixed(1).padStart(9) + String(s.RHlcdmN).padStart(12));
+}
+
+console.log('\nBOSS CMASS (Ntelis et al. 2017): matter, real space, z = 0.57, h^-1 Mpc; their LCDM value is 62.0 in Planck 2015\n');
+for (const [name, par] of [['Planck 2015', H.PLANCK15], ['Planck 2018', H.PLANCK18]]) {
+  const c = H.cosmology(par);
+  console.log(name.padStart(12) + H.homogeneityScale(c, H.amplitude(c, { z: 0.57 })).toFixed(1).padStart(8));
+}
+
+console.log('\nThe acoustic peak: D2 for matter today, linear and with the peak smoothed by nonlinear motions (Sigma = 7 h^-1 Mpc)\n');
+const ps = H.cosmology(H.PLANCK18, { sigmaNL: 7 });
+console.log('r'.padStart(6) + 'linear'.padStart(10) + 'smoothed'.padStart(10));
+for (const r of [90, 100, 105, 110, 113, 116, 120, 125]) console.log(String(r).padStart(6) + H.D2(pc, r).toFixed(5).padStart(10) + H.D2(ps, r).toFixed(5).padStart(10));
+console.log(`peak of xi: ${H.baoPeak(pc).toFixed(1)} (linear), ${H.baoPeak(ps).toFixed(1)} (smoothed) h^-1 Mpc`);
diff --git a/tests/cosmic-homogeneity.mjs b/tests/cosmic-homogeneity.mjs
new file mode 100644
index 0000000..7f4f39f
--- /dev/null
+++ b/tests/cosmic-homogeneity.mjs
@@ -0,0 +1,77 @@
+// Checks the cosmic homogeneity model against published numbers and its own closed forms. Run: npm test.
+// Each check states the source, the published or expected value, the model's value and the tolerance.
+import { createRequire } from 'node:module';
+const require = createRequire(import.meta.url);
+const H = require('../cosmic-homogeneity/homogeneity.js');
+
+let fails = 0;
+const fmt = x => (x === 0 ? '0' : x.toExponential(2));
+function within(name, got, want, factor, src) {
+  const ok = got > 0 && got / want <= factor && want / got <= factor;
+  console.log(`${ok ? 'PASS' : 'FAIL'}  ${name}\n      model ${fmt(got)}  expected ${fmt(want)}  (within x${factor})  ${src}`);
+  if (!ok) fails++;
+}
+function check(name, ok, detail, src) {
+  console.log(`${ok ? 'PASS' : 'FAIL'}  ${name}\n      ${detail}  ${src}`);
+  if (!ok) fails++;
+}
+const S12 = 'Scrimgeour et al., MNRAS 425, 116 (2012), arXiv:1205.6812';
+
+// 1. Scrimgeour et al.'s own LCDM prediction for WiggleZ: their fiducial cosmology, and each slice's measured bias b^2
+//    (Table 4, fitted separately for the D2 and the N statistic) and redshift-space beta (Table 3). Linear theory with
+//    the Kaiser factor only.
+const sc = H.cosmology(H.SCRIMGEOUR12);
+const rows = H.WIGGLEZ.map(s => {
+  const amp = b2 => H.amplitude(sc, { z: s.z, b2, beta: s.beta, redshiftSpace: true });
+  return { ...s, model: H.homogeneityScale(sc, amp(s.b2D)), modelN: H.homogeneityScale(sc, amp(s.b2N), { by: 'N' }) };
+});
+check('WiggleZ LCDM homogeneity scale (D2 = 2.97), four slices, within 3 h^-1 Mpc',
+  rows.every(s => Math.abs(s.model - s.RHlcdm) <= 3),
+  rows.map(s => `z ${s.z}: ${s.model.toFixed(1)} vs ${s.RHlcdm}`).join('; '), S12 + ', Table 4 and Conclusions');
+check('Measured WiggleZ homogeneity scales lie within 2 sigma of the model',
+  rows.every(s => Math.abs(s.model - s.RH) <= 2 * s.err),
+  rows.map(s => `z ${s.z}: ${s.RH}±${s.err} vs ${s.model.toFixed(1)}`).join('; '), S12 + ', abstract');
+check('The other definition (N = 1.01) gives a few h^-1 Mpc more, as in their LCDM counts, within 3 h^-1 Mpc',
+  rows.every(s => Math.abs(s.modelN - s.RHlcdmN) <= 3 && s.modelN > s.model),
+  rows.map(s => `z ${s.z}: ${s.modelN.toFixed(1)} vs ${s.RHlcdmN}`).join('; '), S12 + ', Table 4');
+
+// 2. Ntelis et al. (2017), BOSS CMASS galaxies, corrected for bias and redshift-space distortion: their abstract gives
+//    LCDM R_H = 62.0 h^-1 Mpc, a real-space matter value (Sec. 5.5, Table 5), in their fiducial Planck 2015 cosmology.
+//    The paper attaches no redshift to it; z = 0.57 is the CMASS effective redshift.
+const p15 = H.cosmology(H.PLANCK15);
+within('Matter homogeneity scale at z = 0.57, Planck 2015 (h^-1 Mpc)', H.homogeneityScale(p15, H.amplitude(p15, { z: 0.57 })), 62.0, 1.05,
+  'Ntelis et al., JCAP 06 (2017) 019, arXiv:1702.02159, abstract (LCDM value for matter)');
+
+// 3. The linear spectrum itself.
+const pc = H.cosmology(H.PLANCK18);
+within('Sound horizon at Planck\'s drag redshift (z = 1059.94), Eisenstein-Hu eq. 6 (Mpc)', pc.tf.soundHorizon(H.PLANCK18.zdrag), H.PLANCK18.rdrag, 1.01,
+  'Planck 2018 r_drag (TT,TE,EE+lowE+lensing); at their own fitted z_d = 1021 the same formula gives 150.9');
+within('Baryon acoustic peak in the linear correlation function (sanity check)', H.baoPeak(pc), 100, 1.05, 'Eisenstein et al., ApJ 633, 560 (2005), ~100 h^-1 Mpc');
+within('sigma8 normalisation', pc.sigma(8), H.PLANCK18.sigma8, 1.0001, 'model');
+
+// 4. Small scales: galaxies with xi = (r/r0)^-gamma, r0 = 5 h^-1 Mpc and gamma = 1.8 (Zehavi et al., ApJ 630, 1 (2005)).
+//    Deep inside r0 the counts grow as r^(3 - gamma) = r^1.2; between 1 and 30 h^-1 Mpc their endpoint slope is ~1.9.
+const PL = 'model, with the power law of Zehavi et al., ApJ 630, 1 (2005)';
+within('Correlation dimension deep inside r0 tends to 3 - gamma', H.D2powerLaw(0.01, 5, 1.8), 1.2, 1.01, PL);
+check('The counts\' slope over 1-30 h^-1 Mpc is an average through the climb from 1.2 to 3', Math.abs(H.slopePowerLaw(1, 30, 5, 1.8) - 1.9) < 0.05,
+  `slope ${H.slopePowerLaw(1, 30, 5, 1.8).toFixed(2)}; D2 is ${H.D2powerLaw(1, 5, 1.8).toFixed(2)} at 1 and ${H.D2powerLaw(30, 5, 1.8).toFixed(2)} at 30 h^-1 Mpc`, PL);
+
+// 5. The homogeneity scale is a convention, not a length nature sets: it moves with the threshold, the tracer and the epoch.
+const r1 = H.homogeneityScale(pc), r01 = H.homogeneityScale(pc, 1, { eps: 0.001 });
+const rz = [0, 0.5, 1, 2.4].map(z => H.homogeneityScale(pc, H.amplitude(pc, { z })));
+check('A ten times stricter threshold more than doubles R_H; R_H falls with redshift for matter',
+  r01 > 2 * r1 && rz.every((v, i) => i === 0 || v < rz[i - 1]),
+  `matter, z = 0: ${r1.toFixed(0)} (1%) and ${r01.toFixed(0)} (0.1%) h^-1 Mpc; z = 0, 0.5, 1, 2.4: ${rz.map(v => v.toFixed(0)).join(', ')}`, 'model');
+// The acoustic peak makes the linear D2 dip slightly near 113 h^-1 Mpc, so there R_H can jump. Nonlinear motions
+// smooth the peak (Sigma ~ 8 h^-1 Mpc before reconstruction in BOSS, Beutler et al. 2017) and the jump goes away.
+function largestJump(cos) {
+  const a = H.amplitude(cos, { z: 0.5, b2: 4 }); let prev = null, best = 0;
+  for (let f = 1; f <= 1.15001; f += 0.005) { const r = H.homogeneityScale(cos, a * f); if (prev !== null) best = Math.max(best, r - prev); prev = r; }
+  return best;
+}
+const jLin = largestJump(pc), jSm = largestJump(H.cosmology(H.PLANCK18, { sigmaNL: 7 }));
+check('Linear theory: R_H jumps on the acoustic plateau; with the peak smoothed it does not',
+  jLin > 10 && jSm < 2, `b = 2, z = 0.5, real space: largest change for a 0.5% amplitude step ${jLin.toFixed(1)} h^-1 Mpc (linear), ${jSm.toFixed(1)} (Sigma = 7 h^-1 Mpc)`, 'model');
+
+console.log(fails ? `\n${fails} check(s) failed` : '\nall checks pass');
+process.exit(fails ? 1 : 0);