#!/usr/bin/env python3 """PERITONEAL CAVITY GEOMETRY — ovoid shape, fluid height, per-layer slabs, lumen pressure. The cavity is modelled as an ELLIPSOID (ovoid) with semi-axes a (width), b (depth), c (half height), lying with its long axes horizontal. Fluid fills it FROM THE BOTTOM to a height h, so: * the fluid HEIGHT is a function of instilled VOLUME and is NOT linear in it — a litre added to an empty cavity raises the level far more than a litre added to a half-full one, because the horizontal cross-section is widest at mid-height; * horizontal SLABS of equal height therefore hold UNEQUAL volumes (thin at the bottom, thick in the middle, thin at the top). That is what makes a real inter-layer concentration gradient possible: the layer at the inflow lumen is small and receives all the fresh dialysate, so it runs the highest glucose, and concentration falls with distance from that lumen; * each slab has its own membrane CONTACT AREA (the ellipsoid wall it touches), so transport per layer is geometric, not an equal split; * the HYDROSTATIC PRESSURE at any lumen follows from how much fluid stands above it. A pressure transducer in the inflow or outflow lumen therefore reads a proxy for cavity VOLUME — with no blood-side measurement of any kind, so it does not disturb the concentration-blind contract. Nothing here reads a patient concentration. Pure geometry and hydrostatics. Run: python3 physics_geometry.py (self-tests) """ import math RHO_G_CMH2O = 1.0 # SOURCE: textbook-standard (water density × g) # 1 cm of water column = 1 cmH2O by definition CMH2O_PER_MMHG = 1.35951 # SOURCE: textbook-standard (1 mmHg = 1.35951 cmH2O) def cavity_axes(weight_kg, fill_capacity_L=None): """Semi-axes (a, b, c) in metres of the ellipsoid peritoneal cavity for a given body weight. Anthropometric anchor: an adult peritoneal cavity tolerates roughly 6 L before intra-abdominal pressure becomes limiting; clinical PD fills are 1.5-3 L. A 70 kg adult is taken as a = 0.150 m (width), b = 0.100 m (depth), c = 0.100 m (half height), giving (4/3)*pi*a*b*c = 6.28 L of total capacity. Linear dimensions scale as weight^(1/3) so capacity scales linearly with body mass. """ s = (weight_kg / 70.0) ** (1.0 / 3.0) a, b, c = 0.150 * s, 0.100 * s, 0.100 * s if fill_capacity_L is not None: # rescale isotropically to a stated capacity cur = 4.0 / 3.0 * math.pi * a * b * c * 1000.0 k = (fill_capacity_L / cur) ** (1.0 / 3.0) a, b, c = a * k, b * k, c * k return a, b, c def capacity_L(a, b, c): """Total cavity volume (L) if filled to the top.""" return 4.0 / 3.0 * math.pi * a * b * c * 1000.0 def volume_at_height(h, a, b, c): """Volume (L) of fluid standing to height h (m) above the cavity floor. Cross-section at height z is an ellipse of semi-axes a*s, b*s with s = sqrt(1 - u^2), u = (z - c)/c, so A(z) = pi*a*b*(1 - u^2) and V(h) = pi*a*b*c * [ (t - t^3/3) + 2/3 ], t = (h - c)/c which is exact, not a quadrature. """ h = min(max(h, 0.0), 2.0 * c) t = (h - c) / c return math.pi * a * b * c * ((t - t ** 3 / 3.0) + 2.0 / 3.0) * 1000.0 def height_at_volume(V_L, a, b, c): """Invert volume_at_height: fluid height (m) for a given volume (L). Monotone -> bisection.""" if V_L <= 0.0: return 0.0 cap = capacity_L(a, b, c) if V_L >= cap: return 2.0 * c lo, hi = 0.0, 2.0 * c for _ in range(60): # 60 halvings: exact to ~1e-18 m mid = 0.5 * (lo + hi) if volume_at_height(mid, a, b, c) < V_L: lo = mid else: hi = mid return 0.5 * (lo + hi) def _radius(z, a, b, c): """Equivalent (geometric-mean) horizontal radius of the wall at height z.""" u = (z - c) / c s2 = max(0.0, 1.0 - u * u) return math.sqrt(a * b) * math.sqrt(s2) def _dS_dz(z, a, b, c): """Wall area per unit height at height z — EXACT analytic integrand for the surface-of-revolution approximation (no numerical derivative). Rewritten 2026-08-07 after an adversarial review measured the FLOOR slab's area 2.13% low. The old code formed dr/dz by a FORWARD DIFFERENCE and evaluated it at the pole, where dr/dz -> infinity and r -> 0; the product is finite but the finite difference is not, so the error landed precisely on the floor slab — the smallest slab, the one the drain sits in, and the one the whole area/volume argument leans on. There is no need to approximate anything. With r(z) = R*sqrt(1-u^2), R = sqrt(a*b), u = (z-c)/c, we have dr/dz = -R*u / (c*sqrt(1-u^2)), so dS/dz = 2*pi*r*sqrt(1 + (dr/dz)^2) = 2*pi*R*sqrt( (1 - u^2) + (R/c)^2 * u^2 ) in which the sqrt(1-u^2) factors cancel exactly. The result is smooth and finite for every z including both poles, so Simpson converges immediately instead of fighting a singularity. """ u = (z - c) / c u = min(1.0, max(-1.0, u)) R = math.sqrt(a * b) return 2.0 * math.pi * R * math.sqrt(max(0.0, (1.0 - u * u) + (R / c) ** 2 * u * u)) def slab_geometry(h, a, b, c, n_layers): """Split the fluid column into `n_layers` slabs of EQUAL HEIGHT and return, per slab: (z_mid, volume_L, wall_area_m2) listed BOTTOM FIRST (index 0 = floor of the cavity). Volume comes from the exact ellipsoid integral. Wall area is the surface of revolution of the equivalent radius, dS = 2*pi*r*sqrt(1 + (dr/dz)^2) dz, integrated numerically per slab — this is the membrane the slab is actually in contact with, so transport scales geometrically rather than by an equal split. """ out = [] if h <= 0.0 or n_layers < 1: return out dz = h / n_layers for i in range(n_layers): z0, z1 = i * dz, (i + 1) * dz v = volume_at_height(z1, a, b, c) - volume_at_height(z0, a, b, c) # wall area: composite Simpson over the slab of the EXACT integrand (see _dS_dz). m = 8 area = 0.0 for k in range(m): za = z0 + (z1 - z0) * k / m zb = z0 + (z1 - z0) * (k + 1) / m zm = 0.5 * (za + zb) area += (zb - za) / 6.0 * (_dS_dz(za, a, b, c) + 4.0 * _dS_dz(zm, a, b, c) + _dS_dz(zb, a, b, c)) out.append((0.5 * (z0 + z1), max(v, 0.0), max(area, 0.0))) return out def lumen_pressure_cmH2O(h, z_lumen, iap_baseline_cmH2O=3.0): # SOURCE: textbook-standard (resting intra-abdominal pressure ~3 cmH2O) """Gauge pressure (cmH2O) a transducer in a lumen at height z_lumen reads. If the lumen sits BELOW the fluid surface it reads the standing column above it plus the resting intra-abdominal pressure. If it sits above the surface it reads the baseline only — which is itself the signal that the cavity has drained past that lumen. """ if z_lumen >= h: return iap_baseline_cmH2O return iap_baseline_cmH2O + (h - z_lumen) * 100.0 * RHO_G_CMH2O def volume_from_pressure_L(p_cmH2O, z_lumen, a, b, c, iap_baseline_cmH2O=3.0, strict=False): """The device's INVERSE reading: infer cavity volume (L) from a lumen pressure transducer. This is the point of the geometry — a pressure tap in the drain line is a volume sensor that touches no blood and reads no concentration. RESOLUTION LIMIT (documented 2026-08-07 after an adversarial review). A lumen reads only the column standing ABOVE it. Once the fluid line falls below the lumen the transducer reads the baseline and the volume is genuinely UNRESOLVED — the reading is identical for every fill from empty up to z_lumen. The old code returned volume_at_height(z_lumen) in that case, which is a confident-looking point estimate that can be wildly wrong: with the lumen at 10 cm a true 0.8 L fill infers 3.14 L (+293%), and at 12 cm a 0.3 L fill infers 4.07 L (+1257%). The shipped configuration puts the drain lumen on the floor (z=0), where this cannot arise, but this is a public API. `strict=True` returns None when the reading is unresolved rather than an upper bound. The default keeps the old return type for callers that expect a float, but the value is now correctly understood as "at most this much", not "this much". """ col = (p_cmH2O - iap_baseline_cmH2O) / (100.0 * RHO_G_CMH2O) if col <= 0.0: if strict: return None # unresolved: fluid line is below the lumen return volume_at_height(z_lumen, a, b, c) # UPPER BOUND only, not an estimate return volume_at_height(min(2.0 * c, z_lumen + col), a, b, c) if __name__ == '__main__': ok = True def chk(cond, msg): global ok print((' PASS ' if cond else ' FAIL ') + msg) ok = ok and cond a, b, c = cavity_axes(70.0) cap = capacity_L(a, b, c) print(f"70 kg cavity: a={a*100:.1f} cm b={b*100:.1f} cm c={c*100:.1f} cm " f"capacity {cap:.2f} L full height {2*c*100:.1f} cm\n") print("1. volume <-> height inversion") for V in (0.5, 1.0, 2.0, 3.4, 5.0): h = height_at_volume(V, a, b, c) back = volume_at_height(h, a, b, c) chk(abs(back - V) < 1e-9, f"V={V:.2f} L -> h={h*100:5.2f} cm -> V={back:.6f} L") print("\n2. height is NON-LINEAR in volume (the reason slabs differ)") h1 = height_at_volume(1.0, a, b, c) h2 = height_at_volume(2.0, a, b, c) h3 = height_at_volume(3.0, a, b, c) chk((h2 - h1) < h1, f"first litre raises {h1*100:.2f} cm, second only {(h2-h1)*100:.2f} cm") print(f" third litre raises {(h3-h2)*100:.2f} cm") print("\n3. slab volumes sum to the total, and are UNEQUAL") for V in (1.5, 3.4): h = height_at_volume(V, a, b, c) sl = slab_geometry(h, a, b, c, 10) tot = sum(s[1] for s in sl) vols = [s[1] for s in sl] chk(abs(tot - V) < 1e-9, f"V={V} L: slabs sum to {tot:.9f} L") chk(max(vols) / min(vols) > 1.5, f"V={V} L: slab volumes span {min(vols)*1000:.0f}-{max(vols)*1000:.0f} mL " f"(ratio {max(vols)/min(vols):.2f})") print("\n4. wall area per slab, and the area/volume gradient that drives per-layer transport") h = height_at_volume(3.4, a, b, c) sl = slab_geometry(h, a, b, c, 10) ar = [s[2] for s in sl] vo = [s[1] for s in sl] chk(all(x > 0 for x in ar), f"all slab areas > 0 (total {sum(ar)*1e4:.0f} cm2)") av = [(a_ * 1e4) / (v * 1000.0) for a_, v in zip(ar, vo)] # m^2 -> cm^2, L -> mL # The FLOOR slab is the curved bottom cap: little volume, a lot of wall, so the highest # area-to-volume ratio in the cavity. Upper slabs are wide and shallow-walled. That ratio is # what makes each layer transport differently and is the point of resolving them. chk(av[0] == max(av), f"floor slab has the highest area/volume ({av[0]:.4f} cm2/mL)") chk(av[0] / av[-1] > 5.0, f"area/volume spans {av[-1]:.4f}-{av[0]:.4f} cm2/mL (ratio {av[0]/av[-1]:.1f}x)") print("\n5. lumen pressure and its inverse") for V in (0.8, 1.5, 2.5, 3.4): h = height_at_volume(V, a, b, c) p = lumen_pressure_cmH2O(h, 0.0) Vb = volume_from_pressure_L(p, 0.0, a, b, c) chk(abs(Vb - V) < 1e-6, f"V={V:.2f} L -> h={h*100:5.2f} cm -> floor lumen reads {p:6.2f} cmH2O -> V={Vb:.4f} L") print("\n6. a lumen ABOVE the fluid line reads baseline only (drained-past signal)") h = height_at_volume(0.5, a, b, c) chk(abs(lumen_pressure_cmH2O(h, 0.15) - 3.0) < 1e-12, f"lumen at 15 cm with fluid at {h*100:.1f} cm reads baseline 3.00 cmH2O") print("\n" + ("ALL PASS" if ok else "FAILURES")) raise SystemExit(0 if ok else 1)