diff --git a/README.md b/README.md
index ce962eb..7b7175c 100644
--- a/README.md
+++ b/README.md
@@ -16,8 +16,8 @@ and of the methods of its three references:
 The roadmap, assumptions, critical assessment and full reference list are in
 [`docs/IMPLEMENTATION_PLAN.md`](docs/IMPLEMENTATION_PLAN.md).
 
-**Status:** work packages 1 (foundations), 2 (forward modelling) and 3 (ZO-VSP processing, Q
-and the sonic calibration) are implemented. WP4–WP8 (sonic and walkaway anisotropy,
+**Status:** work packages 1 (foundations), 2 (forward modelling), 3 (ZO-VSP processing, Q and
+the sonic calibration) and 4 (sonic anisotropy) are implemented. WP5–WP8 (walkaway anisotropy,
 Q-compensated imaging, the integrated loop, real-data I/O) are still to do.
 
 ![Synthetic well](docs/figures/wp1_well.png)
@@ -95,6 +95,17 @@ f, S = attenuation.first_arrival_spectra(vsp.t, down, pk)
 dts = attenuation.centroid_shift(f, S, (10.0, 100.0))           # t* relative to 450 m
 q = attenuation.interval_q(depths, t30, dts, boundaries=[600.0, 740.0, 800.0])
 drift = calibration.drift(logs.depth, logs.vp, depths, t30)     # VSP minus integrated sonic
+
+# Thomsen gamma from the sonic logs (WP4)
+from iwt import sonic
+from iwt.synthetic import caliper_log, stoneley_log
+
+radius, caliper = caliper_log(model, logs.depth)                # washouts in the shales
+st = stoneley_log(model, logs.depth, sonic.Mud(1540.0, 1250.0), radius)   # 500 Hz, 46 mm tool
+tool = sonic.Borehole(radius=0.108, tool_radius=0.046)
+res = sonic.sonic_anisotropy(logs.vp, logs.vs, logs.rho, st.slowness, sonic.Mud(1500.0, 1250.0),
+                             tool, frequency=500.0, radius=caliper)   # exact TI Stoneley mode
+res.gamma                                                       # (C66 - C44) / (2 C44)
 ```
 
 ## What WP1 provides
@@ -146,6 +157,19 @@ about 1400 slownesses) with `workers=4`.
 | drift and calibration | `iwt.calibration` | sonic time, drift curve, classic knee calibration (slowness shift or scaling), Backus vertical velocity, physical calibration (Kjartansson dispersion + Backus + Thomsen anisotropy term) |
 | figures | `iwt.plotting.depth_profiles`, `trace_columns`; `examples/wp3_zovsp.py` | `docs/figures/wp3_*.png` |
 
+## What WP4 provides
+
+| Plan item | Code | Notes |
+|---|---|---|
+| C33, C44, C66 and gamma (vertical well) | `iwt.sonic.sonic_anisotropy` | Walsh et al. (2007); C66 from the Stoneley slowness, by the tube-wave relation or with the dispersion removed |
+| tube wave with a tool | `tube_wave_slowness`, `c66_from_tube_wave` | the annulus compliance with a concentric tool of given area modulus (Norris, 1990); C66 is the effective modulus (Norris & Sinha, 1993) |
+| exact Stoneley mode | `stoneley_slowness`, `ti_stiffness` | fluid-filled borehole along the axis of a TI formation, with a tool; leaky modes on the outgoing sheet (complex slowness); validated against the tube wave and the Scholte wave |
+| uncertainty | `gamma_jacobian`, `gamma_uncertainty` | exact derivatives for mud slowness and density, radius, tool radius and compliance, Stoneley and dipole slowness, density |
+| mud slowness | `calibrate_mud_slowness` | fit in an interval of assumed gamma |
+| resolution and frequency matching | `running_average`; `shear_frequency`, `qs` | dipole C44 averaged to the Stoneley aperture and moved to the Stoneley frequency |
+| synthetic logs | `iwt.synthetic.caliper_log`, `stoneley_log` | washouts; Stoneley slowness from the exact mode in the Backus medium of the array aperture |
+| figures | `examples/wp4_sonic_anisotropy.py` | `docs/figures/wp4_*.png` |
+
 ## The synthetic well
 
 The Hanssen data are proprietary, and the article gives no velocities, densities, anisotropy or
@@ -425,6 +449,93 @@ Residual drift after calibration, against the 30 Hz phase delays:
 
 ![Drift](docs/figures/wp3_drift.png)
 
+## WP4 results (synthetic well; `python examples/wp4_sonic_anisotropy.py`, about 50 s)
+
+**The data.**
+
+- **Stoneley log:** "measured" at 500 Hz with a rigid 46 mm-radius tool, in mud of 1540 m/s and
+  1250 kg/m³ (assumption 15). It is computed with the exact TI borehole mode in the Backus
+  average of the model over a 1.8 m array, with 0.3 % noise.
+- **Borehole:** the bit radius is 108 mm. Washouts enlarge it by up to about 10 % in the
+  overburden shales, and the caliper reads it with 1 mm noise.
+- **Processing:** assumes a mud velocity of 1500 m/s unless the mud is calibrated.
+- **Truth:** γ of the Backus medium over the same 1.8 m.
+
+**The Stoneley wave is leaky in the overburden shales.** In the Lower Cretaceous shale and
+the Hekkingen, the tube wave is faster than the shear wave travelling along the borehole: it is
+controlled by C66 (γ ≈ 0.3), the shear wave by C44 (Vs0 ≈ 950–1050 m/s). The Stoneley wave
+therefore radiates shear energy into the formation.
+
+Slowness relative to the tube wave, unit averages:
+
+| Unit | 500 Hz | 1 kHz | radiation 1/Q at 500 Hz | at 500 Hz, ε = δ = 0 | at 500 Hz, isotropic with μ = C66 |
+|---|---|---|---|---|---|
+| Lower Cretaceous shale (leaky) | 1.0188 | 1.0596 | 0.0099 | 1.0205 | 1.0077 |
+| Hekkingen (leaky) | 1.0138 | 1.0398 | 0.0034 | 1.0157 | 1.0049 |
+| Stø | 1.0007 | 0.9983 | 0 | 1.0007 | 1.0004 |
+| Fruholmen | 1.0000 | 0.9972 | 0 | 1.0000 | 0.9995 |
+
+- **The leaky mode turns into a trapped one at about 1.2–2 kHz.** Near that transition the measured
+  slowness alone is ambiguous (see the notes), so the Stoneley wave is best analysed at
+  ≤ 500 Hz.
+- **The measured C44 matters more than the unknown C11 and C13.** At 500 Hz, the TI mode with
+  ε = δ = 0 is 0.2 % off. The usual isotropic substitute, which ignores the measured C44, is
+  1.1 % off. In slowness that error is amplified about 3.5 times in C66.
+
+![Stoneley dispersion](docs/figures/wp4_dispersion.png)
+
+**The processing ladder.** Each step adds one correction. Errors are the γ estimate minus the
+true sonic-scale γ, as unit means:
+
+| Step | L. Cret. | Hekkingen | Fuglen | Stø | Nordmela | Fruholmen | Snadd | RMS |
+|---|---|---|---|---|---|---|---|---|
+| true γ | 0.267 | 0.303 | 0.159 | 0.023 | 0.095 | 0.063 | 0.122 | |
+| A tube wave, no tool, nominal mud, bit size | −0.167 | −0.144 | −0.096 | −0.059 | −0.078 | −0.039 | −0.045 | 0.088 |
+| B + tool term | −0.034 | +0.003 | +0.028 | +0.044 | +0.036 | +0.078 | +0.083 | 0.092 |
+| C + caliper | −0.050 | −0.019 | +0.020 | +0.044 | +0.036 | +0.077 | +0.083 | 0.092 |
+| D + dispersion (exact TI mode, ε = δ = 0) | +0.006 | +0.027 | +0.034 | +0.045 | +0.044 | +0.077 | +0.084 | 0.093 |
+| D′ as D, isotropic substitute | −0.030 | −0.007 | +0.025 | +0.043 | +0.038 | +0.073 | +0.078 | 0.088 |
+| E + dipole C44 moved to 500 Hz, matched resolution | +0.053 | +0.065 | +0.068 | +0.058 | +0.061 | +0.095 | +0.104 | 0.104 |
+| F + mud calibrated in the Stø, assumed isotropic (1563 m/s) | −0.012 | −0.010 | −0.027 | −0.023 | −0.024 | −0.035 | −0.041 | 0.045 |
+| **G as E with the true mud (1540 m/s)** | **+0.009** | **+0.014** | **+0.003** | **+0.003** | **+0.003** | **+0.004** | **+0.003** | **0.029** |
+
+- **The tool term is the largest single correction.** It raises γ by 0.10–0.15.
+- **With all corrections and the right mud slowness, every unit is within 0.014** (G). The
+  remaining RMS of 0.029 is noise and thin-bed spikes.
+- **The mud slowness decides the result.** A 2.7 % error (1500 against 1540 m/s) shifts γ by
+  +0.05 to +0.10 (E), more than every other correction. It also hides whether the other
+  corrections are right: D and E look worse than C only because the mud error dominates.
+- **Calibrating the mud in the Stø under the assumption that it is isotropic** gives 1563 m/s:
+  its true γ is 0.023. The result is then γ errors of −0.01 to −0.04. In a real well the
+  Stø is permeable, which would also raise its Stoneley slowness (Winkler et al., 1989) and
+  bias such a calibration further.
+
+![Sonic gamma](docs/figures/wp4_gamma.png)
+
+**The uncertainty budget.** Linear propagation for the final processing, with 1σ input
+errors:
+
+- mud slowness 2 %, mud density 1 %;
+- Stoneley slowness 0.3 %, dipole slowness 1 %;
+- density 15 kg/m³;
+- radius and tool radius 1 mm each.
+
+It is checked by Monte Carlo (4000 draws), which agrees to 0.001–0.003.
+
+| Unit (median) | L. Cret. | Hekkingen | Fuglen | Stø | Nordmela | Fruholmen | Snadd |
+|---|---|---|---|---|---|---|---|
+| σ(γ), total | 0.033 | 0.041 | 0.039 | 0.035 | 0.038 | 0.050 | 0.055 |
+| from the mud slowness | 0.026 | 0.035 | 0.034 | 0.032 | 0.034 | 0.047 | 0.052 |
+| from the dipole | 0.014 | 0.015 | 0.012 | 0.010 | 0.011 | 0.010 | 0.011 |
+| from the Stoneley | 0.008 | 0.010 | 0.009 | 0.008 | 0.008 | 0.010 | 0.011 |
+
+- **The mud slowness dominates, as the plan expected.** Its effect grows in fast rock, where
+  S_T² − S_f² is small.
+- **In the fast sands and the Snadd, σ(γ) ≈ 0.04–0.055.** That is as large as γ itself, so
+  sonic γ there is only an upper-bound indicator.
+- **In the overburden shales (γ ≈ 0.3) it is well resolved.** Those are the units where the
+  article needed the anisotropy.
+
 ## Assumptions
 
 1. **"LFP article" means Lewis et al. (2016)**, the only article in [`References/`](References/).
@@ -474,6 +585,20 @@ Residual drift after calibration, against the 30 Hz phase delays:
     with no Q, which is enough at vertical incidence.
 14. **The well is vertical,** so the anisotropy term of the physical calibration is zero. It is
     implemented with Thomsen's (1986) weak-anisotropy relation for deviated wells.
+15. **WP4 tool and borehole parameters (assumed):**
+    - the Stoneley slowness is measured at 500 Hz, a single-frequency phase slowness. Real
+      array processing gives a spectrum-weighted average of the dispersion curve (Tang, 2003);
+    - the tool is concentric and rigid, with a 46 mm radius;
+    - the mud is 1540 m/s and 1250 kg/m³; processing assumes 1500 m/s and the correct
+      density, from the mud weight;
+    - the dipole reads the axial shear velocity at 2 kHz; flexural dispersion is assumed
+      already removed.
+16. **The Stoneley wave sees the Backus medium of the 1.8 m array aperture** (arithmetic C66).
+    This is exact for layers much thinner than the borehole radius. For thicker layers the
+    tube wave averages the compliance 1/C66 instead, which matters only next to the stringers.
+17. **Qs = Qp / 2** moves the dipole C44 to the Stoneley frequency, because the VSP gives only
+    the P-wave Q.
+18. **No permeability or mudcake effect** on the Stoneley wave, and no eccentricity.
 
 ## Notes found while implementing
 
@@ -502,6 +627,19 @@ Residual drift after calibration, against the 30 Hz phase delays:
 - **The plan's formula 1/Q_int = 1/Q_app − 1/Q_strat holds only for consistent measurements.**
   Q_app and Q_strat must use the same band and estimator. The second-order interaction (the
   multiples are themselves attenuated) is at or below 1 × 10⁻³ in 1/Q here.
+- **The Stoneley wave is leaky in slow VTI shales** (WP4). The tube wave, set by C66, outruns
+  the axial shear wave, set by C44, so the mode radiates; its slowness is complex. The
+  isotropic substitute with shear modulus C66 misses this, and with it more than half of
+  the dispersion at 500 Hz in the shallow shales.
+- **Close to the leaky–trapped transition, the real Stoneley slowness does not determine C66
+  uniquely.** One leaky and one trapped C66 can give the same slowness; the Stoneley attenuation
+  would tell them apart. In tests with random formations this affected 1 sample in 2400, and
+  about 4 % of the samples had no followable root at 500 Hz.
+- **Gamma combines two tools with different frequencies and apertures.**
+  - The dipole's higher frequency (2 kHz against 500 Hz, Qs ≈ 25) biases γ low by 0.01–0.03.
+  - Its shorter aperture, next to stringers, creates spurious γ spikes.
+  - Both are fixed by moving and averaging C44 to the Stoneley frequency and aperture before
+    γ is formed.
 - **The plan's WP1 check "Backus averaging of isotropic layers gives ε ≥ 0" was wrong.** The test
   suite found counter-examples, in about 10 % of random stacks. The correct statements are
   γ ≥ 0 and ε ≥ δ (η ≥ 0); ε can be negative when Vp contrasts are small and Vs contrasts large
@@ -515,8 +653,9 @@ iwt/
   model.py        ElasticModel1D; Thomsen <-> Cij; stability; phase velocities; Backus; time-depth
   qmodel.py       constant-Q (Kjartansson) dispersion and attenuation
   rockphysics.py  Gassmann fluid substitution, Reuss fluid mixing
-  synthetic.py    the Wisting-like true model and its wireline logs
+  synthetic.py    the Wisting-like true model, its wireline logs, caliper and Stoneley logs
   calibration.py  drift curve; knee and physical (dispersion + Backus + anisotropy) sonic calibration
+  sonic.py        C33, C44, C66 and gamma; tube wave with tool; exact TI Stoneley mode; uncertainty
   plotting.py     well, anisotropy, AVO, gather, depth-profile and trace displays (matplotlib)
   modeling/
     elastic.py       VTI plane waves with constant Q: slownesses, eigenvectors, interface scattering
@@ -529,7 +668,8 @@ iwt/
 examples/wp1_synthetic_well.py      prints the WP1 numbers and writes docs/figures/wp1_*.png
 examples/wp2_forward_modelling.py   AVO, normal-incidence check, multiples, ZO-VSP; docs/figures/wp2_*.png
 examples/wp3_zovsp.py               ZO-VSP processing, Q, drift calibration; docs/figures/wp3_*.png
-tests/            pytest suite (93 tests)
+examples/wp4_sonic_anisotropy.py    Stoneley dispersion, gamma processing ladder, uncertainty; wp4_*.png
+tests/            pytest suite (111 tests)
 docs/             implementation plan and figures
 References/       the article (PDF)
 ```
@@ -542,7 +682,9 @@ Run the tests with `pytest`. The plotting tests are skipped when matplotlib is m
 * Backus, G. E. (1962). Long-wave elastic anisotropy produced by horizontal layering. *Journal of Geophysical Research*, 67(11), 4427–4440. https://doi.org/10.1029/JZ067i011p04427
 * Banik, N. C., Lerche, I., & Shuey, R. T. (1985). Stratigraphic filtering, Part I: Derivation of the O'Doherty–Anstey formula. *Geophysics*, 50(12), 2768–2774. https://doi.org/10.1190/1.1441897
 * Berryman, J. G., Grechka, V., & Berge, P. A. (1999). Analysis of Thomsen parameters for finely layered VTI media. *Geophysical Prospecting*, 47(6), 959–978. https://doi.org/10.1046/j.1365-2478.1999.00163.x
+* Biot, M. A. (1952). Propagation of elastic waves in a cylindrical bore containing a fluid. *Journal of Applied Physics*, 23(9), 997–1005. https://doi.org/10.1063/1.1702365
 * Cavalca, M., Fletcher, R., & Riedel, M. (2013). Q-compensation in complex media—Ray-based and wavefield extrapolation approaches. *SEG Technical Program Expanded Abstracts 2013*, 3831–3835. https://doi.org/10.1190/segam2013-1130.1
+* Cheng, C. H., & Toksöz, M. N. (1981). Elastic wave propagation in a fluid-filled borehole and synthetic acoustic logs. *Geophysics*, 46(7), 1042–1053. https://doi.org/10.1190/1.1441242
 * Fryer, G. J., & Frazer, L. N. (1984). Seismic waves in stratified anisotropic media. *Geophysical Journal International*, 78(3), 691–710. https://doi.org/10.1111/j.1365-246X.1984.tb05065.x
 * Gassmann, F. (1951). Über die Elastizität poröser Medien. *Vierteljahrsschrift der Naturforschenden Gesellschaft in Zürich*, 96, 1–23.
 * Goetz, J. F., Dupal, L., & Bowler, J. (1979). An investigation into discrepancies between sonic log and seismic check shot velocities. *APEA Journal*, 19(1), 131–141. https://doi.org/10.1071/AJ78014
@@ -558,6 +700,8 @@ Run the tests with `pytest`. The plotting tests are skipped when matplotlib is m
 * Marion, D., Mukerji, T., & Mavko, G. (1994). Scale effects on velocity dispersion: From ray to effective medium theories in stratified media. *Geophysics*, 59(10), 1613–1619. https://doi.org/10.1190/1.1443550
 * Mavko, G., Mukerji, T., & Dvorkin, J. (2009). *The Rock Physics Handbook* (2nd ed.). Cambridge University Press. https://doi.org/10.1017/CBO9780511626753
 * Molyneux, J. B., & Schmitt, D. R. (1999). First-break timing: Arrival onset times by direct correlation. *Geophysics*, 64(5), 1492–1501. https://doi.org/10.1190/1.1444653
+* Norris, A. N. (1990). The speed of a tube wave. *Journal of the Acoustical Society of America*, 87(1), 414–417. https://doi.org/10.1121/1.399262
+* Norris, A. N., & Sinha, B. K. (1993). Weak elastic anisotropy and the tube wave. *Geophysics*, 58(8), 1091–1098. https://doi.org/10.1190/1.1443493
 * O'Doherty, R. F., & Anstey, N. A. (1971). Reflections on amplitudes. *Geophysical Prospecting*, 19(3), 430–458. https://doi.org/10.1111/j.1365-2478.1971.tb00610.x
 * Quan, Y., & Harris, J. M. (1997). Seismic attenuation tomography using the frequency shift method. *Geophysics*, 62(3), 895–905. https://doi.org/10.1190/1.1444197
 * Richards, P. G., & Menke, W. (1983). The apparent attenuation of a scattering medium. *Bulletin of the Seismological Society of America*, 73(4), 1005–1021. https://doi.org/10.1785/BSSA0730041005
@@ -566,10 +710,15 @@ Run the tests with `pytest`. The plotting tests are skipped when matplotlib is m
 * Seeman, B., & Horowicz, L. (1983). Vertical seismic profiling: Separation of upgoing and downgoing acoustic waves in a stratified medium. *Geophysics*, 48(5), 555–568. https://doi.org/10.1190/1.1441486
 * Shuey, R. T. (1985). A simplification of the Zoeppritz equations. *Geophysics*, 50(4), 609–614. https://doi.org/10.1190/1.1441936
 * Stewart, R. R., Huddleston, P. D., & Kan, T. K. (1984). Seismic versus sonic velocities: A vertical seismic profiling study. *Geophysics*, 49(8), 1153–1168. https://doi.org/10.1190/1.1441745
+* Tang, X. M. (2003). Determining formation shear-wave transverse isotropy from borehole Stoneley-wave measurements. *Geophysics*, 68(1), 118–126. https://doi.org/10.1190/1.1543199
 * Thomsen, L. (1986). Weak elastic anisotropy. *Geophysics*, 51(10), 1954–1966. https://doi.org/10.1190/1.1442051
+* Valero, H.-P., Djikpesse, H. A., & Sinha, B. K. (2009). Estimation of borehole fluid slowness using sonic array waveforms. *SEG Technical Program Expanded Abstracts 2009*, 361–365. https://doi.org/10.1190/1.3255617
 * Walsh, J., Sinha, B., Plona, T., & Ammerman, M. (2007). Derivation of anisotropy parameters in a shale using borehole sonic data. *SEG Technical Program Expanded Abstracts 2007*, 323–327. https://doi.org/10.1190/1.2792435
+* White, J. E. (1983). *Underground Sound: Application of Seismic Waves*. Elsevier.
+* White, J. E., & Tongtaow, C. (1981). Cylindrical waves in transversely isotropic media. *Journal of the Acoustical Society of America*, 70(4), 1147–1155. https://doi.org/10.1121/1.386946
 * White, R. E., & Simm, R. (2003). Tutorial: Good practice in well ties. *First Break*, 21(10). https://doi.org/10.3997/1365-2397.21.10.25640
+* Winkler, K. W., Liu, H.-L., & Johnson, D. L. (1989). Permeability and borehole Stoneley waves: Comparison between experiment and theory. *Geophysics*, 54(1), 66–75. https://doi.org/10.1190/1.1442578
 * OMV (2014). OMV finds hydrocarbons at Hanssen well in Norway (press release). https://www.omv.com/en/media/press-releases/2014/omv-finds-hydrocarbons-at-hanssen-well-in-norway
 * Sokkeldirektoratet. Wellbore 7324/7-2 factpage. https://factpages.sodir.no/en/wellbore/pageview/exploration/all/7450
 
-DOIs were checked against the Scite index; the WP3 references were also checked against publisher records (the White & Simm page range could not be confirmed and is omitted). The plan lists further references for WP4–WP8.
+DOIs were checked against the Scite index; the WP3 and WP4 references were also checked against publisher records (the White & Simm page range could not be confirmed and is omitted). The exact tool-term notation of Norris (1990) could not be checked against its full text; the form used here was derived from the annulus compliance. The plan lists further references for WP5–WP8.
diff --git a/docs/IMPLEMENTATION_PLAN.md b/docs/IMPLEMENTATION_PLAN.md
index 94b0831..f75cda9 100644
--- a/docs/IMPLEMENTATION_PLAN.md
+++ b/docs/IMPLEMENTATION_PLAN.md
@@ -1,6 +1,6 @@
 # Implementation plan: integrated borehole and surface seismic well tie (Lewis et al., 2016)
 
-Status: **WP1, WP2 and WP3 implemented** (see the README); WP4–WP8 to do. This document says what will be built, in what
+Status: **WP1–WP4 implemented** (see the README); WP5–WP8 to do. This document says what will be built, in what
 order, how each part will be checked, and where the method can be improved. It follows the
 conventions of the earlier repositories in this series (for example `ES_PHYS_P010`): a Python
 package, a `paper.py` that records where each published number comes from, figures redrawn on
@@ -237,7 +237,24 @@ effect is the Backus term, and using Q_app for the dispersion as well counts it
 - recover a known drift from a known Q and frequency pair;
 - the corridor stack correlates with the true reflectivity.
 
-### WP4: Sonic anisotropy (Walsh et al., 2007)
+### WP4: Sonic anisotropy (Walsh et al., 2007) (done)
+
+*As built:* `iwt.sonic`, plus synthetic Stoneley and caliper logs in `iwt.synthetic`. Changes
+from the plan:
+
+- **Tool term:** built as planned. **Eccentricity is not modelled.** Norris (1990) reports
+  that it lowers the tube-wave speed, but the size of the effect was not quantified here.
+- **C66 citation:** that the tube wave senses C66 is now cited to Norris & Sinha (1993).
+  Norris (1990) is cited for the tool term.
+- **Added: the exact Stoneley mode.** The borehole mode is solved exactly in the TI formation,
+  with a tool (White & Tongtaow, 1981, for the TI fields). This removes the dispersion bias
+  of the tube-wave formula.
+  - The Stoneley wave in the slow Wisting-like shales turned out to be *leaky*. The tube
+    wave, set by C66, outruns the axial shear wave, set by C44.
+  - The usual isotropic substitute, with C66 as its shear modulus, misses this.
+- **Added: the dipole's frequency and resolution.** The dipole C44 is moved to the Stoneley
+  frequency and averaged to the Stoneley aperture before gamma is formed.
+- **Added: mud-slowness calibration** in a reference interval.
 
 **Builds** the vertical-well case:
 
@@ -479,6 +496,7 @@ All DOIs were checked against the Scite index on 2026-10-13, except where noted.
 * Kjartansson, E. (1979). Constant Q-wave propagation and attenuation. *Journal of Geophysical Research*, 84(B9), 4737–4748. https://doi.org/10.1029/JB084iB09p04737
 * Miller, D. E., & Spencer, C. (1994). An exact inversion for anisotropic moduli from phase slowness data. *Journal of Geophysical Research*, 99(B11), 21651–21657. https://doi.org/10.1029/94JB01848
 * Norris, A. N. (1990). The speed of a tube wave. *Journal of the Acoustical Society of America*, 87(1), 414–417. https://doi.org/10.1121/1.399262
+* Norris, A. N., & Sinha, B. K. (1993). Weak elastic anisotropy and the tube wave. *Geophysics*, 58(8), 1091–1098. https://doi.org/10.1190/1.1443493 (added in WP4)
 * O'Doherty, R. F., & Anstey, N. A. (1971). Reflections on amplitudes. *Geophysical Prospecting*, 19(3), 430–458. https://doi.org/10.1111/j.1365-2478.1971.tb00610.x
 * Quan, Y., & Harris, J. M. (1997). Seismic attenuation tomography using the frequency shift method. *Geophysics*, 62(3), 895–905. https://doi.org/10.1190/1.1444197
 * Rüger, A. (1997). P-wave reflection coefficients for transversely isotropic models with vertical and horizontal axis of symmetry. *Geophysics*, 62(3), 713–722. https://doi.org/10.1190/1.1444181
@@ -490,6 +508,7 @@ All DOIs were checked against the Scite index on 2026-10-13, except where noted.
 * Walden, A. T., & White, R. E. (1998). Seismic wavelet estimation: A frequency domain solution to a geophysical noisy input-output problem. *IEEE Transactions on Geoscience and Remote Sensing*, 36(1), 287–297. https://doi.org/10.1109/36.655337
 * Wang, Y. (2002). A stable and efficient approach of inverse Q filtering. *Geophysics*, 67(2), 657–663. https://doi.org/10.1190/1.1468627
 * Wang, Y. (2006). Inverse Q-filter for seismic resolution enhancement. *Geophysics*, 71(3), V51–V60. https://doi.org/10.1190/1.2192912
+* White, J. E., & Tongtaow, C. (1981). Cylindrical waves in transversely isotropic media. *Journal of the Acoustical Society of America*, 70(4), 1147–1155. https://doi.org/10.1121/1.386946 (added in WP4)
 * White, R. E., & Simm, R. (2003). Tutorial: Good practice in well ties. *First Break*, 21(10). https://doi.org/10.3997/1365-2397.21.10.25640
 
 **Facts about the Hanssen well and data access (web sources)**
diff --git a/docs/figures/wp4_dispersion.png b/docs/figures/wp4_dispersion.png
new file mode 100644
index 0000000..3e62491
Binary files /dev/null and b/docs/figures/wp4_dispersion.png differ
diff --git a/docs/figures/wp4_gamma.png b/docs/figures/wp4_gamma.png
new file mode 100644
index 0000000..04ef4b9
Binary files /dev/null and b/docs/figures/wp4_gamma.png differ
diff --git a/examples/wp4_sonic_anisotropy.py b/examples/wp4_sonic_anisotropy.py
new file mode 100644
index 0000000..3c2d3f0
--- /dev/null
+++ b/examples/wp4_sonic_anisotropy.py
@@ -0,0 +1,230 @@
+"""Work package 4: Thomsen's gamma from the sonic logs of the vertical synthetic well.
+
+C44 comes from the dipole shear log, C66 from the Stoneley log through the tube-wave
+relation of Norris (1990), and gamma = (C66 - C44) / (2 C44) (Walsh et al., 2007).
+
+1. Stoneley dispersion in the units of the well: the exact borehole mode in the TI
+   formation, against the two shortcuts (eps = delta = 0, and an isotropic formation
+   with shear modulus C66).  The slow overburden shales make the Stoneley wave leaky.
+2. A processing ladder, each step adding one correction, against the true
+   sonic-scale gamma: tool term, caliper, dispersion, the dipole's frequency and
+   resolution, and the mud slowness (nominal, calibrated in the Stø, or true).
+3. The gamma uncertainty budget by linear propagation, checked by Monte Carlo.
+
+The Stoneley log is "measured" at 500 Hz with a rigid 46 mm tool in mud of
+1540 m/s and 1250 kg/m^3; processing assumes 1500 m/s unless calibrated.
+
+Writes ``docs/figures/wp4_*.png`` and prints the numbers quoted in the README.
+
+    python examples/wp4_sonic_anisotropy.py
+"""
+
+from pathlib import Path
+
+import numpy as np
+
+from iwt import sonic as S
+from iwt.model import thomsen_to_cij
+from iwt.plotting import unit_tops
+from iwt.synthetic import caliper_log, stoneley_log, wisting_like_well
+
+FIG = Path(__file__).resolve().parents[1] / "docs" / "figures"
+MUD_TRUE = S.Mud(1540.0, 1250.0)
+MUD_NOMINAL = S.Mud(1500.0, 1250.0)
+BIT = 0.108
+TOOL = S.Borehole(BIT, 0.046)
+F_ST = 500.0
+APERTURE = 1.8
+UNITS = [(470.0, 600.0, "Lower Cretaceous shale"), (600.0, 631.0, "Hekkingen"),
+         (631.0, 672.0, "Fuglen"), (672.0, 702.0, "Stø"), (702.0, 740.0, "Nordmela"),
+         (740.0, 800.0, "Fruholmen"), (800.0, 1670.0, "Snadd")]
+SIGMA = {  # 1-sigma input errors for the budget
+    "s_stoneley": 0.003, "s_shear": 0.01, "rho": 15.0, "mud_slowness": 0.02,
+    "mud_density": 0.01, "radius": 0.001, "tool_radius": 0.001,
+}
+
+
+def per_unit(z, x, fn=np.nanmean):
+    return np.array([fn(x[(z >= a) & (z < b)]) for a, b, _ in UNITS])
+
+
+def dispersion_curves(model):
+    """Stoneley dispersion for unit averages: exact TI, eps = delta = 0, equivalent isotropic."""
+    f = np.geomspace(50.0, 3000.0, 40)
+    blk = model.block([z for z, _ in unit_tops(model)[1:]] + [694.5])
+    out = {}
+    print(f"1. Stoneley dispersion (unit Backus averages, rigid tool, mud {MUD_TRUE.velocity:.0f} m/s):"
+          " slowness relative to the tube wave at 500 Hz / 1 kHz, and the radiation 1/Q at 500 Hz")
+    for name in ("Lower Cretaceous shale", "Hekkingen", "Stø", "Fruholmen"):
+        i = list(blk.label).index(name if name != "Stø" else "Stø (oil)")
+        vp0, vs0, rho = blk.vp0[i], blk.vs0[i], blk.rho[i]
+        c = thomsen_to_cij(vp0, vs0, rho, blk.eps[i], blk.delta[i], blk.gamma[i])
+        cij = (c.c11, c.c13, c.c33, c.c44, c.c66)
+        st = S.tube_wave_slowness(c.c66, MUD_TRUE.slowness, MUD_TRUE.density, BIT, TOOL.tool_radius)
+        args = (MUD_TRUE.slowness, MUD_TRUE.density, BIT, TOOL.tool_radius)
+        exact = S.stoneley_slowness(f, cij, rho, *args, full=True)
+        ti0 = S.stoneley_slowness(f, S.ti_stiffness(vp0, vs0, rho, c.c66), rho, *args)
+        iso = S.stoneley_slowness(f, S.ti_stiffness(vp0, np.sqrt(c.c66 / rho), rho, c.c66), rho, *args)
+        out[name] = dict(f=f, exact=exact.real / st, ti0=ti0 / st, iso=iso / st,
+                         inv_q=2 * exact.imag / exact.real)
+        k500, k1k = np.argmin(abs(f - 500)), np.argmin(abs(f - 1000))
+        leaky = "leaky" if st < 1.0 / vs0 else "trapped"
+        print(f"  {name:24s} gamma {blk.gamma[i]:.2f}, {leaky:7s}: exact {exact.real[k500] / st:.4f} / "
+              f"{exact.real[k1k] / st:.4f}; eps=delta=0 {ti0[k500] / st:.4f}; equivalent isotropic "
+              f"{iso[k500] / st:.4f}; 1/Q {2 * exact.imag[k500] / exact.real[k500]:.4f}")
+    return out
+
+
+def main():
+    model, logs = wisting_like_well()
+    z = logs.depth
+    radius, caliper = caliper_log(model, z)
+    st = stoneley_log(model, z, MUD_TRUE, radius, TOOL.tool_radius, frequency=F_ST,
+                      aperture=APERTURE)
+    truth = model.backus(APERTURE, respect=[model.interfaces[1]]).at(z)["gamma"]
+    curves = dispersion_curves(model)
+
+    # resolution matching and the dipole's frequency (Qs = Qp / 2 assumed)
+    c44_m = S.running_average(z, logs.rho * logs.vs ** 2, APERTURE, harmonic=True)
+    rho_m = S.running_average(z, logs.rho, APERTURE)
+    vs_m = np.sqrt(c44_m / rho_m)
+    qs = 0.5 * model.qp[model.layer_index(z)]
+    f_dip = logs.meta["f_sonic_s"]
+    base = (logs.vp, logs.vs, logs.rho, st.slowness)
+    matched = (logs.vp, vs_m, rho_m, st.slowness)
+    ref = (z >= 672.0) & (z < 702.0)                     # the Stø, assumed isotropic
+    sf_cal = S.calibrate_mud_slowness(logs.vp[ref], vs_m[ref], rho_m[ref], st.slowness[ref],
+                                      MUD_NOMINAL.density, TOOL, 0.0, F_ST, caliper[ref],
+                                      shear_frequency=f_dip, qs=qs[ref])
+    mud_cal = S.Mud(1.0 / sf_cal, MUD_NOMINAL.density)
+    ladder = {
+        "A tube wave, no tool, nominal mud, bit size":
+            S.sonic_anisotropy(*base, MUD_NOMINAL, S.Borehole(BIT)),
+        "B + tool term (Norris, 1990)": S.sonic_anisotropy(*base, MUD_NOMINAL, TOOL),
+        "C + caliper": S.sonic_anisotropy(*base, MUD_NOMINAL, TOOL, radius=caliper),
+        "D + dispersion (exact TI mode, eps = delta = 0)":
+            S.sonic_anisotropy(*base, MUD_NOMINAL, TOOL, F_ST, caliper),
+        "D' (as D, equivalent isotropic formation)":
+            S.sonic_anisotropy(*base, MUD_NOMINAL, TOOL, F_ST, caliper, equivalent_isotropic=True),
+        "E + dipole C44 at 500 Hz, matched resolution":
+            S.sonic_anisotropy(*matched, MUD_NOMINAL, TOOL, F_ST, caliper,
+                               shear_frequency=f_dip, qs=qs),
+        f"F + mud calibrated in the Stø ({mud_cal.velocity:.0f} m/s)":
+            S.sonic_anisotropy(*matched, mud_cal, TOOL, F_ST, caliper, shear_frequency=f_dip, qs=qs),
+        "G as E with the true mud (1540 m/s)":
+            S.sonic_anisotropy(*matched, MUD_TRUE, TOOL, F_ST, caliper, shear_frequency=f_dip, qs=qs),
+    }
+    print("\n2. Gamma error (estimate - truth) per unit, mean; RMS over the well; samples without a solution")
+    print("  " + " " * 48 + " ".join(f"{n[:9]:>9s}" for _, _, n in UNITS))
+    print("  " + f"{'true gamma':48s}" + " ".join(f"{x:9.3f}" for x in per_unit(z, truth)))
+    for name, r in ladder.items():
+        e = r.gamma - truth
+        print(f"  {name:48s}" + " ".join(f"{x:+9.3f}" for x in per_unit(z, e))
+              + f"   RMS {np.sqrt(np.nanmean(e ** 2)):.3f}, NaN {int(np.isnan(e).sum())}")
+    final = list(ladder.values())[6]
+
+    # 3. uncertainty budget of the final processing (tube-wave relation at the corrected slowness)
+    s_tube = st.slowness / final.dispersion
+    J = S.gamma_jacobian(s_tube, 1.0 / vs_m, rho_m, mud_cal.slowness, mud_cal.density, caliper,
+                         TOOL.tool_radius, 0.0)
+    sig = dict(SIGMA)
+    sig["s_stoneley"] = SIGMA["s_stoneley"] * s_tube
+    sig["s_shear"] = SIGMA["s_shear"] / vs_m
+    sig["mud_slowness"] = SIGMA["mud_slowness"] * mud_cal.slowness
+    sig["mud_density"] = SIGMA["mud_density"] * mud_cal.density
+    total, contrib = S.gamma_uncertainty(J, sig)
+    print("\n3. Gamma standard deviation per unit (linear propagation) and its largest terms")
+    print("  " + " " * 30 + " ".join(f"{n[:9]:>9s}" for _, _, n in UNITS))
+    print("  " + f"{'total':30s}" + " ".join(f"{x:9.3f}" for x in per_unit(z, total, np.nanmedian)))
+    for k in ("mud_slowness", "s_stoneley", "s_shear", "radius", "tool_radius", "mud_density", "rho"):
+        print("  " + f"{k:30s}" + " ".join(f"{x:9.3f}" for x in per_unit(z, contrib[k], np.nanmedian)))
+    rng = np.random.default_rng(1)
+    print("  Monte Carlo check (4000 draws) at the unit midpoints, linear / sampled std:")
+    for a, b, name in UNITS:
+        i = np.searchsorted(z, 0.5 * (a + b))
+        x = dict(s=s_tube[i], ss=1 / vs_m[i], rho=rho_m[i], sf=mud_cal.slowness,
+                 rf=mud_cal.density, b=caliper[i], a=TOOL.tool_radius)
+        n = 4000
+        d = {k: x[k] + sd * rng.standard_normal(n) for k, sd in
+             (("s", sig["s_stoneley"][i]), ("ss", sig["s_shear"][i]), ("rho", sig["rho"]),
+              ("sf", sig["mud_slowness"]), ("rf", sig["mud_density"]), ("b", sig["radius"]),
+              ("a", sig["tool_radius"]))}
+        g = S.gamma_from_moduli(S.c66_from_tube_wave(d["s"], d["sf"], d["rf"], d["b"], d["a"]),
+                                d["rho"] / d["ss"] ** 2)
+        print(f"    {name:24s} {total[i]:.3f} / {np.nanstd(g):.3f}")
+
+    try:
+        import matplotlib
+
+        matplotlib.use("Agg")
+        from iwt.plotting import PALETTE, depth_profiles
+    except ImportError:
+        print("\nmatplotlib not installed: no figures")
+        return
+    import matplotlib.pyplot as plt
+
+    FIG.mkdir(parents=True, exist_ok=True)
+    tops = [(t, n) for t, n in unit_tops(model)[1:] if t > 440.0]
+    us = 1e6
+    order = list(ladder)
+
+    def smooth(x):
+        x = np.where(np.isfinite(x), x, np.nan)
+        ok = np.isfinite(x)
+        return np.interp(z, z[ok], S.running_average(z[ok], x[ok], 5.0))
+    fig = depth_profiles({
+        "logs": ("slowness (us/m)", {
+            "monopole P (12 kHz)": (us / logs.vp, z),
+            "dipole shear (2 kHz)": (us / logs.vs, z),
+            "Stoneley (500 Hz)": (us * st.slowness, z),
+            "mud (nominal)": (np.full(z.size, us * MUD_NOMINAL.slowness), z)}),
+        "gamma (5 m running mean)": ("gamma", {
+            "true (sonic scale)": (smooth(truth), z),
+            "A: tube wave, nominal mud, no tool": (smooth(ladder[order[0]].gamma), z),
+            "F: all corrections, calibrated mud": (smooth(final.gamma), z),
+            "G: all corrections, true mud": (smooth(ladder[order[7]].gamma), z)}),
+        "gamma error budget (1 sigma)": ("standard deviation", {
+            "total": (total, z),
+            "mud slowness (2 %)": (contrib["mud_slowness"], z),
+            "Stoneley (0.3 %) + dipole (1 %)":
+                (np.hypot(contrib["s_stoneley"], contrib["s_shear"]), z),
+            "radius + tool (1 mm each)": (np.hypot(contrib["radius"], contrib["tool_radius"]), z)}),
+    }, top=460.0, bottom=1670.0, tops=tops, title="Sonic gamma in the vertical synthetic well "
+       "(Stoneley at 500 Hz, 46 mm tool; thin-bed spikes clipped)", width=3.4)
+    fig.axes[1].set_xlim(-0.2, 0.6)
+    fig.axes[2].set_xlim(0.0, 0.12)
+    fig.savefig(FIG / "wp4_gamma.png", dpi=150)
+
+    fig, axes = plt.subplots(1, 2, figsize=(9.6, 4.2))
+    fig.patch.set_facecolor(PALETTE["surface"])
+    styles = ("-", "--", "-.", ":")
+    for k, (name, c) in enumerate(curves.items()):
+        col = PALETTE["series"][k]
+        axes[0].plot(c["f"], c["exact"], color=col, linestyle="-", linewidth=1.4, label=f"{name}",
+                     marker="osD^"[k], markevery=6, markersize=3.5)
+        axes[0].plot(c["f"], c["iso"], color=col, linestyle=":", linewidth=1.2)
+        axes[1].plot(c["f"], c["inv_q"], color=col, linestyle=styles[k], linewidth=1.4, label=name)
+    axes[0].plot([], [], color=PALETTE["ink2"], linestyle=":", label="dotted: equivalent isotropic")
+    axes[0].set_ylabel("Stoneley slowness / tube-wave slowness", color=PALETTE["ink2"], fontsize=8)
+    axes[1].set_ylabel("radiation attenuation 1/Q", color=PALETTE["ink2"], fontsize=8)
+    for ax, t in zip(axes, ("Dispersion: exact TI mode (solid)", "Leakage into the formation")):
+        ax.set_xscale("log")
+        ax.set_xlabel("frequency (Hz)", color=PALETTE["ink2"], fontsize=8)
+        ax.set_title(t, fontsize=8.5, color=PALETTE["ink"], loc="left")
+        ax.tick_params(colors=PALETTE["muted"], labelsize=7, length=2)
+        ax.grid(True, color=PALETTE["grid"], linewidth=0.5)
+        ax.set_facecolor(PALETTE["surface"])
+        for spine in ax.spines.values():
+            spine.set_color(PALETTE["axis"])
+        leg = ax.legend(fontsize=6.5, frameon=False, loc="upper left")
+        for txt in leg.get_texts():
+            txt.set_color(PALETTE["ink2"])
+    fig.suptitle("Stoneley wave in the units of the synthetic well (unit averages, 46 mm rigid tool)",
+                 fontsize=9.5, color=PALETTE["ink"], x=0.02, ha="left")
+    fig.tight_layout()
+    fig.savefig(FIG / "wp4_dispersion.png", dpi=150)
+    print(f"\nfigures written to {FIG}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/iwt/__init__.py b/iwt/__init__.py
index bdc4153..2729305 100644
--- a/iwt/__init__.py
+++ b/iwt/__init__.py
@@ -12,11 +12,12 @@
 substitution (:mod:`iwt.rockphysics`) and a synthetic Wisting-like well
 (:mod:`iwt.synthetic`).  Work package 2 adds forward modelling
 (:mod:`iwt.modeling`), work package 3 ZO-VSP processing and Q
-(:mod:`iwt.vsp`) and the sonic-to-seismic calibration (:mod:`iwt.calibration`).
+(:mod:`iwt.vsp`) and the sonic-to-seismic calibration (:mod:`iwt.calibration`),
+and work package 4 Thomsen's gamma from the sonic logs (:mod:`iwt.sonic`).
 Plotting (:mod:`iwt.plotting`) needs matplotlib.
 """
 
-from . import calibration, paper, qmodel, rockphysics
+from . import calibration, paper, qmodel, rockphysics, sonic
 from .model import (
     ElasticModel1D,
     Stiffness,
@@ -29,13 +30,14 @@
 )
 from .synthetic import WellLogs, synthetic_logs, wisting_like_model, wisting_like_well
 
-__version__ = "0.3.0"
+__version__ = "0.4.0"
 
 __all__ = [
     "calibration",
     "paper",
     "qmodel",
     "rockphysics",
+    "sonic",
     "ElasticModel1D",
     "Stiffness",
     "backus_average",
diff --git a/iwt/sonic.py b/iwt/sonic.py
new file mode 100644
index 0000000..09da7be
--- /dev/null
+++ b/iwt/sonic.py
@@ -0,0 +1,517 @@
+"""Sonic anisotropy in a vertical well: C33, C44, C66 and Thomsen's gamma.
+
+In a vertical well through a VTI formation (Walsh et al., 2007):
+
+* the monopole P wave travels along the symmetry axis, ``C33 = rho Vp^2``;
+* the dipole shear waves travel along the axis, ``C44 = rho Vs^2``;
+* the low-frequency Stoneley (tube) wave senses an effective shear modulus
+  (Norris & Sinha, 1993), which for a borehole along the axis of a VTI
+  formation is ``C66``: the borehole wall expands in plane strain in the
+  isotropic plane;
+* ``gamma = (C66 - C44) / (2 C44)``.
+
+**The tube-wave relation.**  At zero frequency the Stoneley wave is a tube wave
+whose slowness follows from the compliance of the fluid annulus.  For a
+borehole of radius b with a concentric tool of radius a and area modulus K_t
+(``lambda + mu`` of a solid rod in plane strain; inf for a rigid tool),
+
+    S_T^2 = S_f^2 + rho_f (b^2 / C66 + a^2 / K_t) / (b^2 - a^2),
+
+with S_f and rho_f the mud slowness and density.  Without a tool this is the
+classic ``S_f^2 + rho_f / mu`` (White, 1983).  A rigid tool of 46 mm radius in
+a 108 mm borehole raises the C66 term by 22 %, so the tool term is not small.
+Norris (1990) treats the tool, casing and eccentricity; an eccentred tool or
+borehole lowers the tube-wave speed further, which is not modelled here.
+
+**Dispersion.**  Sonic tools measure the Stoneley wave at 0.5-3 kHz, where it is
+no longer a tube wave.  :func:`stoneley_slowness` solves the exact modal
+equation of a fluid-filled borehole, with an optional concentric tool, along the
+axis of a TI formation (Biot, 1952, and Cheng & Toksoz, 1981, for the isotropic
+case; White & Tongtaow, 1981, for the TI fields).  The dispersion depends on
+C11 and C13 as well, which a vertical well does not measure, so
+:func:`sonic_anisotropy` assumes Thomsen eps and delta (default 0).  In the slow
+Wisting-like shales the 1 kHz slowness is about 3 % above the tube-wave value,
+which biases C66 by about -10 % if ignored.  Where the tube wave (set by C66)
+is faster than the shear wave along the borehole (set by C44), the Stoneley
+wave radiates shear energy into the formation and is leaky; this happens in the
+slow, strongly anisotropic Wisting-like overburden shales.  An "equivalent
+isotropic" formation with shear modulus C66, a common shortcut (cf. Tang,
+2003), misses this and gets the dispersion wrong there.
+
+**Uncertainty.**  C66 depends on the *difference* ``S_T^2 - S_f^2``, so errors
+in the mud slowness are amplified, most strongly in fast formations.  The mud
+slowness is not measured in the hole at sonic frequencies (Valero et al.,
+2009); :func:`calibrate_mud_slowness` fits it in an interval of assumed gamma.
+Permeability raises the Stoneley slowness in sands (Winkler et al., 1989), which
+is not modelled here and would bias both C66 and such a calibration.
+:func:`gamma_jacobian` gives the exact derivatives of gamma with respect to
+every input, for linear error propagation.
+
+References
+----------
+Biot, M. A., 1952, Propagation of elastic waves in a cylindrical bore containing
+a fluid: Journal of Applied Physics, 23(9), 997-1005.
+https://doi.org/10.1063/1.1702365
+Cheng, C. H., & Toksoz, M. N., 1981, Elastic wave propagation in a fluid-filled
+borehole and synthetic acoustic logs: Geophysics, 46(7), 1042-1053.
+https://doi.org/10.1190/1.1441242
+Norris, A. N., 1990, The speed of a tube wave: JASA, 87(1), 414-417.
+https://doi.org/10.1121/1.399262
+Norris, A. N., & Sinha, B. K., 1993, Weak elastic anisotropy and the tube wave:
+Geophysics, 58(8), 1091-1098. https://doi.org/10.1190/1.1443493
+Tang, X. M., 2003, Determining formation shear-wave transverse isotropy from
+borehole Stoneley-wave measurements: Geophysics, 68(1), 118-126.
+https://doi.org/10.1190/1.1543199
+Valero, H.-P., Djikpesse, H. A., & Sinha, B. K., 2009, Estimation of borehole
+fluid slowness using sonic array waveforms: SEG Expanded Abstracts, 361-365.
+https://doi.org/10.1190/1.3255617
+Walsh, J., Sinha, B., Plona, T., & Ammerman, M., 2007, Derivation of anisotropy
+parameters in a shale using borehole sonic data: SEG Expanded Abstracts,
+323-327. https://doi.org/10.1190/1.2792435
+White, J. E., 1983, Underground Sound: Application of Seismic Waves: Elsevier.
+White, J. E., & Tongtaow, C., 1981, Cylindrical waves in transversely isotropic
+media: JASA, 70(4), 1147-1155. https://doi.org/10.1121/1.386946
+Winkler, K. W., Liu, H.-L., & Johnson, D. L., 1989, Permeability and borehole
+Stoneley waves: Comparison between experiment and theory: Geophysics, 54(1),
+66-75. https://doi.org/10.1190/1.1442578
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+from typing import NamedTuple, Optional
+
+import numpy as np
+
+__all__ = [
+    "Mud",
+    "Borehole",
+    "tube_wave_slowness",
+    "c66_from_tube_wave",
+    "stoneley_slowness",
+    "ti_stiffness",
+    "gamma_from_moduli",
+    "running_average",
+    "SonicAnisotropy",
+    "sonic_anisotropy",
+    "GAMMA_INPUTS",
+    "gamma_jacobian",
+    "gamma_uncertainty",
+    "calibrate_mud_slowness",
+]
+
+
+@dataclass(frozen=True)
+class Mud:
+    """Borehole fluid: velocity (m/s) and density (kg/m^3)."""
+
+    velocity: float = 1500.0
+    density: float = 1200.0
+
+    @property
+    def slowness(self) -> float:
+        return 1.0 / self.velocity
+
+
+@dataclass(frozen=True)
+class Borehole:
+    """Borehole radius and a concentric tool (radius 0: no tool; modulus inf: rigid)."""
+
+    radius: float = 0.108
+    tool_radius: float = 0.0
+    tool_modulus: float = np.inf
+
+
+# -- the tube wave (zero frequency) -----------------------------------------------
+
+
+def _tool_term(tool_radius, tool_modulus):
+    a = np.asarray(tool_radius, dtype=float)
+    with np.errstate(divide="ignore"):
+        kt_inv = np.where(np.isinf(tool_modulus), 0.0, 1.0 / np.asarray(tool_modulus, float))
+    return a * a * kt_inv
+
+
+def tube_wave_slowness(c66, mud_slowness, mud_density, radius, tool_radius=0.0,
+                       tool_modulus=np.inf):
+    """Zero-frequency Stoneley (tube-wave) slowness (s/m), Norris (1990)."""
+    b2 = np.asarray(radius, float) ** 2
+    a2 = np.asarray(tool_radius, float) ** 2
+    return np.sqrt(np.asarray(mud_slowness, float) ** 2
+                   + mud_density * (b2 / np.asarray(c66, float)
+                                    + _tool_term(tool_radius, tool_modulus)) / (b2 - a2))
+
+
+def c66_from_tube_wave(s_tube, mud_slowness, mud_density, radius, tool_radius=0.0,
+                       tool_modulus=np.inf):
+    """C66 (Pa) from the tube-wave slowness: the inverse of :func:`tube_wave_slowness`.
+
+    NaN where the slowness is not larger than the tool and mud terms allow.
+    """
+    b2 = np.asarray(radius, float) ** 2
+    a2 = np.asarray(tool_radius, float) ** 2
+    den = ((np.asarray(s_tube, float) ** 2 - np.asarray(mud_slowness, float) ** 2)
+           * (b2 - a2) / mud_density - _tool_term(tool_radius, tool_modulus))
+    with np.errstate(divide="ignore", invalid="ignore"):
+        return np.where(den > 0, b2 / den, np.nan)
+
+
+def gamma_from_moduli(c66, c44):
+    """Thomsen's gamma = (C66 - C44) / (2 C44)."""
+    return 0.5 * (np.asarray(c66, float) / np.asarray(c44, float) - 1.0)
+
+
+# -- the exact Stoneley mode --------------------------------------------------------
+
+
+def _radial(S, c11, c13, c33, c44, rho):
+    """Squared radial slownesses Y (two roots) and polarizations w of the formation waves.
+
+    For (complex) axial slowness S the axisymmetric fields of a TI formation
+    (axis along the borehole) are ``u_r = U K1(kappa r)``,
+    ``u_z = i w U K0(kappa r)``, with ``kappa^2 = omega^2 Y`` and Y a root of the
+    Christoffel equation with an imaginary horizontal slowness.  Returns
+    ``(Y, w)``, each of shape (2, len(S)).
+    """
+    a2 = c11 * c44
+    a1 = c11 * (rho - c33 * S ** 2) + c44 * (rho - c44 * S ** 2) + S ** 2 * (c13 + c44) ** 2
+    a0 = (rho - c44 * S ** 2) * (rho - c33 * S ** 2)
+    sq = np.sqrt(a1 * a1 - 4.0 * a2 * a0 + 0j)
+    Y = np.array([(-a1 - sq) / (2.0 * a2), (-a1 + sq) / (2.0 * a2)])
+    kap = np.sqrt(Y)
+    c = c13 + c44
+    # two equivalent forms of the polarization; use the better-conditioned one
+    d1 = rho + c44 * Y - c33 * S ** 2
+    d2 = S * kap * c
+    with np.errstate(divide="ignore", invalid="ignore"):
+        w1 = S * kap * c / d1
+        w2 = (rho + c11 * Y - c44 * S ** 2) / d2
+    w = np.where(np.abs(d1) >= 1e-6 * np.abs(c11 * Y), w1, w2)
+    return Y, w
+
+
+def _kappa(Y):
+    """Radial wavenumber / omega: evanescent (Re > 0) or outgoing (Hankel-1) radiation."""
+    ev = np.sqrt(Y)
+    out = -1j * np.sqrt(-Y)
+    return np.where(Y.real >= 0.0, ev, out)
+
+
+def _period(S, w, cij, rho, sf, rhof, b, a, kt):
+    """Scaled period determinant of the axisymmetric borehole modes at (complex) slownesses S.
+
+    A formation wave with Re(Y) < 0 radiates; it takes the outgoing branch, so
+    the root is a leaky mode with complex slowness (Im S > 0: decaying along the
+    borehole).
+    """
+    from scipy.special import ive, kve
+
+    c11, c13, c33, c44, c66 = cij
+    S = np.atleast_1d(np.asarray(S, dtype=complex))
+    k = w * S
+    Y, w_y = _radial(S, c11, c13, c33, c44, rho)
+    q = _kappa(Y)
+    pol = np.where(Y.real >= 0.0, w_y, w_y * (q / np.sqrt(Y)))
+    kap = w * q
+    f = w * np.sqrt(S ** 2 - sf ** 2)
+    K0, K1 = kve(0, kap * b), kve(1, kap * b)
+    ur = -K1                                                     # (fluid - solid) u_r
+    srr = (c11 * kap + c13 * k * pol) * K0 + 2.0 * c66 * K1 / b  # -(solid sigma_rr)
+    srz = c44 * (k - kap * pol) * K1                             # solid sigma_rz / i
+    pf = rhof * w ** 2
+    sc = c44 * w                                                 # row scale of stresses
+    zero = np.zeros_like(S)
+    col = [[ur[j], srr[j] / sc, srz[j] / sc] for j in range(2)]
+    if a <= 0.0:
+        M = np.array([
+            [f * ive(1, f * b), col[0][0], col[1][0]],
+            [-pf * ive(0, f * b) / sc, col[0][1], col[1][1]],
+            [zero, col[0][2], col[1][2]],
+        ])
+    else:
+        ea = np.exp(-f * (b - a))
+        c = 0.0 if np.isinf(kt) else pf * a / (2.0 * kt)
+        M = np.array([
+            [f * ive(1, f * b), -f * kve(1, f * b) * ea, col[0][0], col[1][0]],
+            [-pf * ive(0, f * b) / sc, -pf * kve(0, f * b) * ea / sc, col[0][1], col[1][1]],
+            [zero, zero, col[0][2], col[1][2]],
+            [(f * ive(1, f * a) + c * ive(0, f * a)) * ea,
+             -f * kve(1, f * a) + c * kve(0, f * a), zero, zero],
+        ])
+    return np.linalg.det(np.moveaxis(M, -1, 0))
+
+
+def _track(frequency: float, cij, rho, sf, rhof, b, a, kt, start: float = 20.0,
+           ratio: float = 1.15, max_jump: float = 0.03):
+    """Stoneley root at one frequency for many formations at once (complex, NaN on failure).
+
+    The root is followed from the tube-wave slowness at ``frequency / start`` in
+    geometric frequency steps of ``ratio``, with a few secant iterations per
+    step.  A formation is dropped (NaN) when its root jumps by more than
+    ``max_jump`` in one step, or when it leaves the physical region.
+    """
+    c11, c13, c33, c44, c66, rho, b = np.broadcast_arrays(
+        *(np.atleast_1d(np.asarray(c, float)) for c in (*cij, rho, b)))
+    cij = (c11, c13, c33, c44, c66)
+    n_steps = max(1, int(np.ceil(np.log(start) / np.log(ratio))))
+    freqs = frequency / start * start ** (np.arange(n_steps + 1) / n_steps)
+    S = tube_wave_slowness(c66, sf, rhof, b, a, kt).astype(complex)
+    S = np.where(np.isfinite(S), S, np.nan)
+    for fr in freqs:
+        w = 2.0 * np.pi * fr
+        live = np.isfinite(S)
+        if not live.any():
+            break
+        idx = np.flatnonzero(live)
+        sub = tuple(c[idx] for c in cij)
+
+        def g(x):
+            with np.errstate(all="ignore"):
+                return _period(x, w, sub, rho[idx], sf, rhof, b[idx], a, kt)
+
+        x0 = S[idx]
+        x1 = x0 * (1.0 + 1e-6)
+        f0, f1 = g(x0), g(x1)
+        done = np.zeros(idx.size, dtype=bool)
+        for _ in range(15):
+            act = ~done
+            with np.errstate(all="ignore"):
+                x2 = x1[act] - f1[act] * (x1[act] - x0[act]) / (f1[act] - f0[act])
+            x2 = np.where(np.isfinite(x2), x2, x1[act])
+            x0[act], f0[act] = x1[act], f1[act]
+            x1[act] = x2
+            sub_act = tuple(c[act] for c in sub)
+            with np.errstate(all="ignore"):
+                f1[act] = _period(x2, w, sub_act, rho[idx][act], sf, rhof, b[idx][act], a, kt)
+            # converged samples are frozen: further secant steps on a zero residual are noise
+            done |= np.abs(x1 - x0) <= 1e-12 * np.abs(x1)
+            if done.all():
+                break
+        ok = (done & np.isfinite(x1) & (x1.real > sf) & (x1.imag > -1e-9 * np.abs(x1))
+              & (np.abs(x1 - S[idx]) < max_jump * np.abs(S[idx])))
+        S[idx] = np.where(ok, x1, np.nan)
+    return S
+
+
+def stoneley_slowness(frequency, cij, rho, mud_slowness, mud_density, radius,
+                      tool_radius=0.0, tool_modulus=np.inf, full=False):
+    """Stoneley phase slowness (s/m) of a fluid-filled borehole along the axis of a TI formation.
+
+    Root of the exact modal (period) equation of a borehole of ``radius`` in a
+    transversely isotropic formation with stiffnesses
+    ``cij = (C11, C13, C33, C44, C66)`` (Pa) and density ``rho``, with an
+    optional concentric tool of ``tool_radius`` whose wall follows the pressure
+    with area modulus ``tool_modulus`` (inf: rigid).  For an isotropic formation
+    this is the equation of Biot (1952) and Cheng & Toksoz (1981); the TI fields
+    follow White & Tongtaow (1981).
+
+    Either ``frequency`` (Hz) is an array and the formation is one set of
+    scalars (a dispersion curve), or ``frequency`` is a scalar and the
+    stiffnesses, ``rho`` and ``radius`` are arrays (a log; solved for all
+    depths at once).
+
+    The root is followed in frequency from the tube-wave limit (Norris, 1990).
+    Where the Stoneley wave is faster than a formation wave travelling along the
+    borehole it radiates into the formation and is leaky, with a complex
+    slowness.  ``full=True`` returns the complex slowness (Im S > 0; the
+    Stoneley attenuation from radiation is ``1/Q = 2 Im S / Re S``), otherwise
+    its real part.  NaN where the root cannot be followed, which happens where a
+    leaky mode turns into a trapped one.  With no tool and an isotropic
+    formation the slowness tends to the Scholte wave of a flat fluid-solid
+    interface at high frequency (Biot, 1952).
+    """
+    args = (rho, mud_slowness, mud_density, radius, tool_radius, tool_modulus)
+    if np.ndim(frequency):
+        freq = np.asarray(frequency, dtype=float)
+        out = np.array([_track(fr, cij, *args)[0] for fr in freq.ravel()]).reshape(freq.shape)
+    else:
+        out = _track(float(frequency), cij, *args)
+        if not any(np.ndim(x) for x in (*cij, rho, radius)):
+            out = out.reshape(())[()]
+    return out if full else np.real(out)
+
+
+def ti_stiffness(vp0, vs0, rho, c66, eps=0.0, delta=0.0):
+    """``(C11, C13, C33, C44, C66)`` from vertical velocities, C66 and Thomsen eps, delta."""
+    c33, c44 = rho * vp0 ** 2, rho * vs0 ** 2
+    c11 = c33 * (1.0 + 2.0 * eps)
+    c13 = np.sqrt(2.0 * c33 * (c33 - c44) * delta + (c33 - c44) ** 2) - c44
+    return c11, c13, c33, c44, c66
+
+
+# -- the anisotropy -----------------------------------------------------------------
+
+
+def running_average(depth, values, length: float, harmonic: bool = False) -> np.ndarray:
+    """Centred running (arithmetic or harmonic) average over ``length`` (m) on a regular grid.
+
+    Use it to bring two logs to a common resolution before combining them: the
+    dipole C44 (harmonic, as Backus) and the density to the Stoneley array
+    aperture, for example.  The window is shortened at the ends.
+    """
+    d = np.asarray(depth, dtype=float)
+    x = np.asarray(values, dtype=float)
+    dz = float(np.median(np.diff(d)))
+    n = max(1, int(round(length / dz)) | 1)
+    if n == 1:
+        return x.copy()
+    v = 1.0 / x if harmonic else x
+    k = np.ones(n)
+    m = np.convolve(v, k, mode="same") / np.convolve(np.ones_like(v), k, mode="same")
+    return 1.0 / m if harmonic else m
+
+
+
+class SonicAnisotropy(NamedTuple):
+    """Stiffnesses (Pa), gamma, and the dispersion ratio S(f) / S(0) used for C66."""
+
+    c33: np.ndarray
+    c44: np.ndarray
+    c66: np.ndarray
+    gamma: np.ndarray
+    dispersion: np.ndarray
+
+
+def sonic_anisotropy(vp, vs, rho, s_stoneley, mud: Mud, borehole: Borehole,
+                     frequency: Optional[float] = None, radius=None, eps=0.0, delta=0.0,
+                     equivalent_isotropic: bool = False, shear_frequency: Optional[float] = None,
+                     qs=None, max_iter: int = 30, tol: float = 1e-9) -> SonicAnisotropy:
+    """C33, C44, C66 and gamma from monopole P, dipole shear and Stoneley logs.
+
+    ``vp`` and ``vs`` (m/s) and ``rho`` (kg/m^3) are the logs, ``s_stoneley``
+    (s/m) the Stoneley slowness measured at ``frequency`` (Hz).  ``radius``
+    (array, e.g. the caliper) overrides ``borehole.radius``.
+
+    With ``frequency=None`` the Stoneley slowness is taken as the tube-wave
+    slowness.  Otherwise C66 is found by a secant iteration (in log C66, for all
+    depths at once) such that the exact borehole mode (:func:`stoneley_slowness`)
+    reproduces the measured slowness.  The formation is TI with C33 and C44 from
+    the logs and C11, C13 from the assumed Thomsen ``eps`` and ``delta``; with
+    ``equivalent_isotropic`` it is isotropic with shear modulus C66 instead (the
+    usual shortcut, which ignores the measured C44).  ``dispersion`` is
+    ``S(f) / S(0)`` at the solution; NaN where no solution was found.
+
+    The dipole measures C44 at its own frequency (about 2 kHz), above the
+    Stoneley frequency.  With ``shear_frequency`` and ``qs``, C44 is moved to
+    the Stoneley ``frequency`` with the constant-Q law before gamma is formed;
+    otherwise gamma is biased low by about ``(1 + 2 gamma) g ln(f_dipole / f)``
+    with ``g = arctan(1/Qs) / pi`` (0.02-0.03 for Qs = 25 and a factor 4).
+    """
+    vp, vs, rho = (np.asarray(x, dtype=float) for x in (vp, vs, rho))
+    s_st = np.asarray(s_stoneley, dtype=float)
+    shape = np.broadcast(vp, vs, rho, s_st).shape
+    b = np.broadcast_to(borehole.radius if radius is None else np.asarray(radius, float), shape)
+    vp, vs, rho, s_st = (np.broadcast_to(x, shape).astype(float) for x in (vp, vs, rho, s_st))
+    eps, delta = (np.broadcast_to(np.asarray(x, float), shape) for x in (eps, delta))
+    sf, rf = mud.slowness, mud.density
+    a, kt = borehole.tool_radius, borehole.tool_modulus
+    c33, c44 = rho * vp ** 2, rho * vs ** 2
+    if frequency is not None and shear_frequency is not None and qs is not None:
+        from .qmodel import dispersion_exponent
+
+        c44 = c44 * (frequency / shear_frequency) ** (2.0 * dispersion_exponent(qs))
+        vs = np.sqrt(c44 / rho)
+    c66 = c66_from_tube_wave(s_st, sf, rf, b, a, kt)
+    ratio = np.ones(shape)
+    if frequency is not None:
+
+        def model(c):
+            if equivalent_isotropic:
+                cij = ti_stiffness(vp, np.sqrt(c / rho), rho, c)
+            else:
+                cij = ti_stiffness(vp, vs, rho, c, eps, delta)
+            return stoneley_slowness(frequency, cij, rho, sf, rf, b, a, kt)
+
+        x0 = np.log(c66)
+        f0 = model(c66) - s_st
+        # where the model is undefined (too soft: past the leaky-trapped transition), stiffen
+        for _ in range(10):
+            bad = ~np.isfinite(f0) & np.isfinite(x0)
+            if not bad.any():
+                break
+            x0 = np.where(bad, x0 + 0.1, x0)
+            f0 = np.where(bad, model(np.exp(x0)) - s_st, f0)
+        x1 = x0 + 0.02
+        f1 = model(np.exp(x1)) - s_st
+        for _ in range(max_iter):
+            with np.errstate(all="ignore"):
+                x2 = x1 - f1 * (x1 - x0) / (f1 - f0)
+            x2 = np.where(np.isfinite(x2), np.clip(x2, x1 - 0.3, x1 + 0.3), x1)
+            x0, f0, x1 = x1, f1, x2
+            f1 = model(np.exp(x1)) - s_st
+            if np.all(~np.isfinite(f1) | (np.abs(x1 - x0) < tol)):
+                break
+        c66 = np.where(np.isfinite(f1) & (np.abs(f1) < 1e-6 * s_st), np.exp(x1), np.nan)
+        ratio = s_st / tube_wave_slowness(c66, sf, rf, b, a, kt)
+    return SonicAnisotropy(c33, c44, c66, gamma_from_moduli(c66, c44), ratio)
+
+
+#: Inputs of :func:`gamma_jacobian`, in order.
+GAMMA_INPUTS = ("s_stoneley", "s_shear", "rho", "mud_slowness", "mud_density", "radius",
+                "tool_radius", "tool_compliance")
+
+
+def gamma_jacobian(s_stoneley, s_shear, rho, mud_slowness, mud_density, radius,
+                   tool_radius=0.0, tool_compliance=0.0):
+    """Exact derivatives of gamma (tube-wave relation) with respect to :data:`GAMMA_INPUTS`.
+
+    ``s_stoneley`` is the (dispersion-corrected) tube-wave slowness and
+    ``tool_compliance`` = 1 / K_t.  Returns an array ``[..., 8]`` of partial
+    derivatives (per s/m, per kg/m^3, per m, per 1/Pa).
+    """
+    S, Ss, rho, Sf, rf, b, a, ct = np.broadcast_arrays(
+        *(np.asarray(x, dtype=float) for x in (s_stoneley, s_shear, rho, mud_slowness,
+                                                mud_density, radius, tool_radius,
+                                                tool_compliance)))
+    e = S ** 2 - Sf ** 2
+    ann = b ** 2 - a ** 2
+    den = e * ann / rf - a ** 2 * ct
+    c66 = b ** 2 / den
+    c44 = rho / Ss ** 2
+    dg_dc66 = 0.5 / c44
+    dg_dc44 = -0.5 * c66 / c44 ** 2
+    dc66 = {
+        "s_stoneley": -b ** 2 * (2.0 * S * ann / rf) / den ** 2,
+        "mud_slowness": b ** 2 * (2.0 * Sf * ann / rf) / den ** 2,
+        "mud_density": b ** 2 * (e * ann / rf ** 2) / den ** 2,
+        "radius": 2.0 * b / den - b ** 2 * (2.0 * b * e / rf) / den ** 2,
+        "tool_radius": -b ** 2 * (-2.0 * a * e / rf - 2.0 * a * ct) / den ** 2,
+        "tool_compliance": b ** 2 * a ** 2 / den ** 2,
+    }
+    dc44 = {"s_shear": -2.0 * rho / Ss ** 3, "rho": 1.0 / Ss ** 2}
+    return np.stack([dg_dc66 * dc66.get(k, 0.0) + dg_dc44 * dc44.get(k, 0.0)
+                     for k in GAMMA_INPUTS], axis=-1)
+
+
+def gamma_uncertainty(jacobian, sigmas: dict):
+    """Linear propagation: ``(std, contributions)`` of gamma for independent input errors.
+
+    ``sigmas`` maps names in :data:`GAMMA_INPUTS` to standard deviations (scalars
+    or arrays); ``contributions`` maps each name to its share ``|dg/dx| sigma_x``.
+    """
+    contrib = {k: np.abs(jacobian[..., GAMMA_INPUTS.index(k)]) * np.asarray(v, float)
+               for k, v in sigmas.items()}
+    total = np.sqrt(sum(c ** 2 for c in contrib.values()))
+    return total, contrib
+
+
+def calibrate_mud_slowness(vp, vs, rho, s_stoneley, mud_density: float, borehole: Borehole,
+                           gamma_ref: float = 0.0, frequency: Optional[float] = None,
+                           radius=None, bounds=(1.0 / 1800.0, 1.0 / 1300.0), **kwargs) -> float:
+    """Mud slowness (s/m) for which the mean gamma over the given samples equals ``gamma_ref``.
+
+    The samples should come from an interval believed to have a known gamma
+    (usually assumed isotropic).  Any error in that belief maps one to one into
+    every gamma of the well, and permeability raises the Stoneley slowness in
+    sands, which this calibration would wrongly absorb into the mud slowness.
+    Further keyword arguments (``eps``, ``delta``, ``shear_frequency``, ``qs``,
+    ...) go to :func:`sonic_anisotropy`.
+    """
+    from scipy.optimize import brentq
+
+    def misfit(sf):
+        g = sonic_anisotropy(vp, vs, rho, s_stoneley, Mud(1.0 / sf, mud_density), borehole,
+                             frequency, radius, **kwargs).gamma
+        return float(np.nanmean(g)) - gamma_ref
+
+    return brentq(misfit, bounds[0], bounds[1], xtol=1e-12)
diff --git a/iwt/synthetic.py b/iwt/synthetic.py
index 8d17820..f028a54 100644
--- a/iwt/synthetic.py
+++ b/iwt/synthetic.py
@@ -34,6 +34,11 @@
 
 In a vertical VTI well the monopole sonic measures Vp0 and both dipole modes
 measure Vs0.
+
+For WP4, :func:`caliper_log` adds washouts in the overburden shales and
+:func:`stoneley_log` a Stoneley slowness log, computed with the exact
+borehole mode of :mod:`iwt.sonic` in the running Backus average of the model
+over the array aperture, at the Stoneley frequency.
 """
 
 from __future__ import annotations
@@ -60,6 +65,10 @@
     "wisting_like_model",
     "synthetic_logs",
     "wisting_like_well",
+    "WASHOUT",
+    "caliper_log",
+    "StoneleyLog",
+    "stoneley_log",
 ]
 
 
@@ -321,3 +330,77 @@ def wisting_like_well(seed: int = 2016, log_seed: Optional[int] = 7):
     """``(model, logs)`` for the default synthetic well."""
     model = wisting_like_model(seed=seed)
     return model, synthetic_logs(model, seed=log_seed)
+
+
+#: Assumed washout (fractional enlargement of the bit radius) per unit: the soft,
+#: shallow shales wash out, the sands and the Snadd hardly.
+WASHOUT = {"Lower Cretaceous shale": 0.06, "Hekkingen": 0.08, "Fuglen": 0.03}
+
+
+def caliper_log(model: ElasticModel1D, depth, bit_radius: float = 0.108,
+                washout: Optional[dict] = None, noise: float = 0.001,
+                seed: Optional[int] = 5):
+    """``(radius, caliper)``: true borehole radius and its measurement (m) at ``depth``.
+
+    The enlargement in a unit is ``washout[unit]`` times a positive, smooth
+    random factor (mean about 1, correlation length 2 m); ``noise`` (m) is the
+    caliper's standard deviation.
+    """
+    washout = WASHOUT if washout is None else washout
+    depth = np.asarray(depth, dtype=float)
+    rng = np.random.default_rng(seed)
+    dz = float(np.median(np.diff(depth))) if depth.size > 1 else 1.0
+    factor = np.clip(1.0 + 0.5 * _ar1(depth.size, dz, 2.0, rng), 0.0, None)
+    names = np.array([str(lab).replace(" stringer", "").replace(" (oil)", "")
+                      for lab in model.at(depth)["label"]])
+    frac = np.array([washout.get(n, 0.0) for n in names])
+    radius = bit_radius * (1.0 + frac * factor)
+    caliper = radius + noise * rng.standard_normal(depth.size)
+    return radius, caliper
+
+
+@dataclass
+class StoneleyLog:
+    """Stoneley phase slowness (s/m) and radiation attenuation 1/Q at ``frequency`` (Hz)."""
+
+    depth: np.ndarray
+    slowness: np.ndarray
+    inv_q: np.ndarray
+    frequency: float
+    meta: dict = field(default_factory=dict)
+
+
+def stoneley_log(model: ElasticModel1D, depth, mud, radius, tool_radius: float = 0.046,
+                 tool_modulus: float = np.inf, frequency: float = 500.0,
+                 aperture: float = 1.8, noise: float = 0.003,
+                 seed: Optional[int] = 11) -> StoneleyLog:
+    """Stoneley slowness log "measured" in a vertical well through ``model``.
+
+    The model is dispersed from ``model.f_ref`` to ``frequency`` (constant Q),
+    Backus averaged over the array ``aperture`` (m), and the exact Stoneley mode
+    of the TI formation is computed with :func:`iwt.sonic.stoneley_slowness` for
+    the ``mud`` (:class:`iwt.sonic.Mud`), borehole ``radius`` (m, scalar or per
+    depth) and a concentric tool.  ``noise`` is the relative standard deviation
+    of the slowness.  Where the mode cannot be followed the log is NaN.
+    """
+    from .model import thomsen_to_cij
+    from .sonic import stoneley_slowness
+
+    depth = np.asarray(depth, dtype=float)
+    m = model.copy()
+    m.vp0 = phase_velocity(model.vp0, model.qp, frequency, model.f_ref)
+    m.vs0 = np.where(model.is_fluid, 0.0,
+                     phase_velocity(np.where(model.is_fluid, 1.0, model.vs0), model.qs, frequency,
+                                    model.f_ref))
+    eff = m.backus(aperture, respect=[model.interfaces[np.argmax(~model.is_fluid)]]).at(depth)
+    c = thomsen_to_cij(eff["vp0"], eff["vs0"], eff["rho"], eff["eps"], eff["delta"], eff["gamma"])
+    rad = np.broadcast_to(np.asarray(radius, dtype=float), depth.shape)
+    s = stoneley_slowness(frequency, (c.c11, c.c13, c.c33, c.c44, c.c66), eff["rho"],
+                          mud.slowness, mud.density, rad, tool_radius, tool_modulus, full=True)
+    inv_q = 2.0 * s.imag / s.real
+    slow = s.real
+    if seed is not None:
+        slow = slow * (1.0 + noise * np.random.default_rng(seed).standard_normal(depth.size))
+    meta = dict(mud_velocity=mud.velocity, mud_density=mud.density, tool_radius=tool_radius,
+                tool_modulus=tool_modulus, aperture=aperture, noise=noise, seed=seed)
+    return StoneleyLog(depth, slow, inv_q, frequency, meta)
diff --git a/pyproject.toml b/pyproject.toml
index 3f56519..58a8e35 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "iwt"
-version = "0.3.0"
+version = "0.4.0"
 description = "Integrated borehole and surface seismic well tie (Lewis et al., 2016, 78th EAGE) in Python"
 readme = "README.md"
 requires-python = ">=3.9"
diff --git a/tests/test_sonic.py b/tests/test_sonic.py
new file mode 100644
index 0000000..2b49dfd
--- /dev/null
+++ b/tests/test_sonic.py
@@ -0,0 +1,221 @@
+"""Tests of the sonic anisotropy (iwt.sonic) and the synthetic Stoneley and caliper logs."""
+
+import numpy as np
+import pytest
+
+from iwt import sonic as S
+from iwt.model import thomsen_to_cij
+from iwt.qmodel import phase_velocity
+
+MUD = S.Mud(1540.0, 1250.0)
+B = 0.108
+A = 0.046
+SHALE = (2400.0, 950.0, 2280.0, 0.20, 0.08, 0.24)       # leaky: slow vertical shear
+SAND = (3350.0, 1800.0, 2420.0, 0.05, 0.02, 0.06)       # fast
+
+
+def _cij(vp0, vs0, rho, eps, delta, gamma):
+    c = thomsen_to_cij(vp0, vs0, rho, eps, delta, gamma)
+    return (c.c11, c.c13, c.c33, c.c44, c.c66)
+
+
+def _iso(vp, vs, rho):
+    return S.ti_stiffness(vp, vs, rho, rho * vs ** 2)
+
+
+# -- tube wave ------------------------------------------------------------------------
+
+
+@pytest.mark.parametrize("a, kt", [(0.0, np.inf), (A, np.inf), (A, 8e9)])
+def test_tube_wave_round_trip(a, kt):
+    c66 = np.array([2.0e9, 5.0e9, 9.0e9])
+    s = S.tube_wave_slowness(c66, MUD.slowness, MUD.density, B, a, kt)
+    np.testing.assert_allclose(S.c66_from_tube_wave(s, MUD.slowness, MUD.density, B, a, kt),
+                               c66, rtol=1e-12)
+    if a == 0.0:   # White's formula
+        np.testing.assert_allclose(s ** 2, MUD.slowness ** 2 + MUD.density / c66, rtol=1e-12)
+    else:          # the tool always slows the tube wave
+        assert np.all(s > S.tube_wave_slowness(c66, MUD.slowness, MUD.density, B))
+
+
+def test_gamma_is_zero_in_an_isotropic_formation():
+    vp, vs, rho = 2800.0, 1400.0, 2350.0
+    s0 = S.tube_wave_slowness(rho * vs ** 2, MUD.slowness, MUD.density, B, A)
+    r = S.sonic_anisotropy(vp, vs, rho, s0, MUD, S.Borehole(B, A))
+    assert abs(float(r.gamma)) < 1e-12
+    s1 = S.stoneley_slowness(500.0, _iso(vp, vs, rho), rho, MUD.slowness, MUD.density, B, A)
+    r = S.sonic_anisotropy(vp, vs, rho, s1, MUD, S.Borehole(B, A), frequency=500.0)
+    assert abs(float(r.gamma)) < 1e-8
+    assert float(r.dispersion) != 1.0
+
+
+# -- the exact borehole mode --------------------------------------------------------
+
+
+@pytest.mark.parametrize("a, kt", [(0.0, np.inf), (A, np.inf), (A, 5e9)])
+def test_stoneley_low_frequency_limit_is_the_tube_wave(a, kt):
+    for p in (SHALE, SAND):
+        cij = _cij(*p)
+        s = S.stoneley_slowness(0.5, cij, p[2], MUD.slowness, MUD.density, B, a, kt)
+        st = S.tube_wave_slowness(cij[4], MUD.slowness, MUD.density, B, a, kt)
+        assert s == pytest.approx(st, rel=1e-6)
+
+
+def test_stoneley_high_frequency_limit_is_the_scholte_wave():
+    from scipy.optimize import brentq
+
+    vp, vs, rho = 3000.0, 1600.0, 2400.0
+    sp2, ss2, sf2 = 1 / vp ** 2, 1 / vs ** 2, MUD.slowness ** 2
+
+    def scholte(x):
+        a, b, f = np.sqrt(x ** 2 - sp2), np.sqrt(x ** 2 - ss2), np.sqrt(x ** 2 - sf2)
+        return (2 * x ** 2 - ss2) ** 2 - 4 * x ** 2 * a * b + MUD.density / rho * ss2 ** 2 * a / f
+
+    lo = max(np.sqrt(ss2), np.sqrt(sf2)) * (1 + 1e-9)
+    grid = np.linspace(lo, 2 * lo, 4000)
+    d = scholte(grid)
+    j = np.flatnonzero(np.sign(d[1:]) != np.sign(d[:-1]))[0]
+    s_sch = brentq(scholte, grid[j], grid[j + 1])
+    s = S.stoneley_slowness(1e6, _iso(vp, vs, rho), rho, MUD.slowness, MUD.density, B)
+    assert s == pytest.approx(s_sch, rel=1e-4)
+
+
+def test_leaky_stoneley_in_slow_anisotropic_shale():
+    """The tube wave (C66) outruns the axial shear wave (C44): the mode radiates."""
+    f = np.array([200.0, 500.0, 1000.0])
+    s = S.stoneley_slowness(f, _cij(*SHALE), SHALE[2], MUD.slowness, MUD.density, B, A, full=True)
+    assert np.all(s.imag > 0) and np.all(np.diff(s.imag) > 0)
+    s_tube = S.tube_wave_slowness(_cij(*SHALE)[4], MUD.slowness, MUD.density, B, A)
+    assert s_tube < 1.0 / SHALE[1]
+    sand = S.stoneley_slowness(f, _cij(*SAND), SAND[2], MUD.slowness, MUD.density, B, A, full=True)
+    np.testing.assert_allclose(sand.imag, 0.0, atol=1e-12)
+    # the measured C44 matters more than the unknown eps and delta at low frequency
+    vp0, vs0, rho, eps, delta, gamma = SHALE
+    c66 = _cij(*SHALE)[4]
+    ti0 = S.stoneley_slowness(500.0, S.ti_stiffness(vp0, vs0, rho, c66), rho, MUD.slowness,
+                              MUD.density, B, A)
+    iso = S.stoneley_slowness(500.0, S.ti_stiffness(vp0, np.sqrt(c66 / rho), rho, c66), rho,
+                              MUD.slowness, MUD.density, B, A)
+    err_ti0, err_iso = abs(ti0 / s[1].real - 1), abs(iso / s[1].real - 1)
+    assert err_ti0 < 4e-3 and err_iso > 8e-3 and err_iso > 3 * err_ti0
+
+
+@pytest.mark.parametrize("p", [SHALE, SAND])
+def test_dispersion_corrected_inversion_recovers_c66(p):
+    cij = _cij(*p)
+    vp0, vs0, rho, eps, delta, gamma = p
+    s = S.stoneley_slowness(500.0, cij, rho, MUD.slowness, MUD.density, B, A)
+    bh = S.Borehole(B, A)
+    r = S.sonic_anisotropy(vp0, vs0, rho, s, MUD, bh, frequency=500.0, eps=eps, delta=delta)
+    assert float(r.c66) == pytest.approx(cij[4], rel=1e-7)
+    assert float(r.gamma) == pytest.approx(gamma, abs=1e-7)
+    naive = S.sonic_anisotropy(vp0, vs0, rho, s, MUD, bh)
+    if p is SHALE:
+        assert float(naive.gamma) < gamma - 0.05     # dispersion ignored: gamma far too low
+
+
+def test_inversion_of_a_log_at_once():
+    rng = np.random.default_rng(3)
+    n = 40
+    vp = 2800 + 400 * rng.random(n)
+    vs = 1300 + 400 * rng.random(n)
+    rho = 2300 + 100 * rng.random(n)
+    gam = 0.15 * rng.random(n)
+    c66 = rho * vs ** 2 * (1 + 2 * gam)
+    s = S.stoneley_slowness(500.0, S.ti_stiffness(vp, vs, rho, c66), rho, MUD.slowness,
+                            MUD.density, B, A)
+    assert s.shape == (n,) and np.all(np.isfinite(s))
+    r = S.sonic_anisotropy(vp, vs, rho, s, MUD, S.Borehole(B, A), frequency=500.0)
+    np.testing.assert_allclose(r.gamma, gam, atol=1e-7)
+
+
+def test_shear_frequency_correction():
+    vp0, vs0, rho, eps, delta, gamma = SAND
+    qs = 30.0
+    vs_dipole = phase_velocity(vs0, qs, 2000.0, 500.0)     # dipole sees a stiffer C44
+    s = S.stoneley_slowness(500.0, _cij(*SAND), rho, MUD.slowness, MUD.density, B, A)
+    bh = S.Borehole(B, A)
+    biased = S.sonic_anisotropy(vp0, vs_dipole, rho, s, MUD, bh, 500.0, eps=eps, delta=delta)
+    fixed = S.sonic_anisotropy(vp0, vs_dipole, rho, s, MUD, bh, 500.0, eps=eps, delta=delta,
+                               shear_frequency=2000.0, qs=qs)
+    assert float(biased.gamma) < gamma - 0.01
+    assert float(fixed.gamma) == pytest.approx(gamma, abs=1e-6)
+
+
+# -- uncertainty and calibration ------------------------------------------------------
+
+
+def test_gamma_jacobian_against_finite_differences():
+    x0 = dict(s_stoneley=9.3e-4, s_shear=1 / 1100.0, rho=2300.0, mud_slowness=1 / 1500.0,
+              mud_density=1200.0, radius=0.108, tool_radius=0.046, tool_compliance=1 / 20e9)
+
+    def gamma(**kw):
+        c66 = S.c66_from_tube_wave(kw["s_stoneley"], kw["mud_slowness"], kw["mud_density"],
+                                   kw["radius"], kw["tool_radius"], 1 / kw["tool_compliance"])
+        return S.gamma_from_moduli(c66, kw["rho"] / kw["s_shear"] ** 2)
+
+    J = S.gamma_jacobian(**x0)
+    for i, k in enumerate(S.GAMMA_INPUTS):
+        h = abs(x0[k]) * 1e-6
+        xp, xm = dict(x0), dict(x0)
+        xp[k] += h
+        xm[k] -= h
+        assert J[i] == pytest.approx((gamma(**xp) - gamma(**xm)) / (2 * h), rel=1e-6), k
+
+
+def test_mud_slowness_dominates_the_gamma_error_in_fast_rock():
+    s_sand = S.tube_wave_slowness(_cij(*SAND)[4], MUD.slowness, MUD.density, B, A)
+    J = S.gamma_jacobian(s_sand, 1 / SAND[1], SAND[2], MUD.slowness, MUD.density, B, A, 0.0)
+    sig = {"s_stoneley": 0.003 * s_sand, "s_shear": 0.01 / SAND[1], "rho": 15.0,
+           "mud_slowness": 0.02 * MUD.slowness, "mud_density": 0.01 * MUD.density,
+           "radius": 0.001, "tool_radius": 0.001}
+    total, contrib = S.gamma_uncertainty(J, sig)
+    assert np.isclose(total, np.sqrt(sum(c ** 2 for c in contrib.values())))
+    assert max(contrib, key=contrib.get) == "mud_slowness"
+    assert contrib["mud_slowness"] > 2 * sorted(contrib.values())[-2]
+
+
+def test_calibrate_mud_slowness_with_a_known_reference():
+    rng = np.random.default_rng(4)
+    vp = 3000 + 200 * rng.random(20)
+    vs = 1600 + 100 * rng.random(20)
+    rho = np.full(20, 2300.0)
+    gam = 0.03
+    s = S.tube_wave_slowness(rho * vs ** 2 * (1 + 2 * gam), MUD.slowness, MUD.density, B, A)
+    sf = S.calibrate_mud_slowness(vp, vs, rho, s, MUD.density, S.Borehole(B, A), gamma_ref=gam)
+    assert sf == pytest.approx(MUD.slowness, rel=1e-9)
+    # assuming the reference is isotropic when it is not makes the mud too fast
+    sf0 = S.calibrate_mud_slowness(vp, vs, rho, s, MUD.density, S.Borehole(B, A), gamma_ref=0.0)
+    assert sf0 < MUD.slowness * 0.995
+
+
+def test_running_average_harmonic_matches_backus_c44():
+    z = np.arange(0.0, 99.0, 0.5)
+    c44 = np.where(np.arange(z.size) % 3 == 0, 2e9, 6e9)
+    h = S.running_average(z, c44, 16.5, harmonic=True)        # 33 samples = 11 periods
+    assert np.allclose(h[50:150], 1.0 / np.mean(1.0 / c44), rtol=1e-12)
+    assert np.allclose(S.running_average(z, c44, 0.0), c44)
+
+
+# -- synthetic logs -----------------------------------------------------------------
+
+
+def test_synthetic_caliper_and_stoneley_logs():
+    from iwt.synthetic import caliper_log, stoneley_log, wisting_like_model
+
+    model = wisting_like_model().subset(0.0, 800.0)
+    z = np.arange(470.0, 780.0, 2.0)
+    rad, cal = caliper_log(model, z)
+    labels = model.at(z)["label"]
+    sand = np.array([str(lab).startswith(("Stø", "Nordmela", "Fruholmen")) for lab in labels])
+    assert np.all(rad[sand] == 0.108) and np.all(rad[~sand] >= 0.108)
+    assert np.max(rad) > 0.112
+    st = stoneley_log(model, z, MUD, rad, seed=None)
+    assert np.isfinite(st.slowness).mean() > 0.95
+    shale = np.array([str(lab).startswith("Lower Cretaceous") for lab in labels])
+    assert np.nanmedian(st.inv_q[shale]) > 1e-3            # leaky in the slow shale
+    assert np.nanmax(np.abs(st.inv_q[sand])) < 1e-9
+    st2 = stoneley_log(model, z, MUD, rad, seed=11)
+    ok = np.isfinite(st.slowness)
+    assert np.array_equal(ok, np.isfinite(st2.slowness))
+    assert 0.002 < np.std(st2.slowness[ok] / st.slowness[ok] - 1) < 0.004