begin
    import Pkg
    Pkg.activate(joinpath(@__DIR__, ".."))

    # a writable font cache for CairoMakie (a notebook worker may have no HOME)
    ENV["XDG_CACHE_HOME"] = get(ENV, "XDG_CACHE_HOME", joinpath(homedir(), ".cache"))
end;

using PlutoUI: TableOfContents

using Printf: @sprintf

using Serialization: serialize, deserialize

using Statistics: mean, median, quantile

import Vasculature as V

import CairoMakie as Mke

import GPUSelect

begin
    # licensed: this file and anything made from it stay on lab storage
    const LEG_NRB = get(
        ENV,
        "XCAT_LEG_NRB",
        joinpath(homedir(), "Documents", "dev", "MolloiLab", "xcat-data", "legs", "both_legs_1.nrb"),
    )
    # every cached stage (pass 1, the patches, each pass-2 patch); on a cluster, a shared scratch folder
    const CACHE_DIR = get(ENV, "LIMB_CACHE_DIR", joinpath(homedir(), "xcat-perfused", "limb_pipeline_cache"))
    mkpath(CACHE_DIR)
end;

"Arrays to this machine's GPU (CUDA, Metal, ROCm, oneAPI), or the CPU if there is none."
const DEVICE = GPUSelect.Storage()

"Run `compute` once and keep its result in `CACHE_DIR/name.jls`; afterwards load it."
function cached(compute, name)
    path = joinpath(CACHE_DIR, name * ".jls")
    isfile(path) && return deserialize(path)

    seconds = @elapsed result = compute()
    serialize(path, result)
    @info "$name: $(round(seconds; digits = 1)) s"
    return result
end

"The six face neighbours of a voxel."
const SIX_STEPS = (
    CartesianIndex(1, 0, 0),
    CartesianIndex(-1, 0, 0),
    CartesianIndex(0, 1, 0),
    CartesianIndex(0, -1, 0),
    CartesianIndex(0, 0, 1),
    CartesianIndex(0, 0, -1),
)

begin
    const ENVELOPE = ["leg_right", "foot_right", "foot_right1", "foot_right2", "foot_right3", "foot_right4", "foot_right5"]
    const BONES = ["r_femur", "r_tibia", "r_fibula", "r_patella", "r_talus", "r_calcaneus", "rpelvis1", "rpelvis2", "rpelvis3"]
    const LONG_BONES = ["r_femur", "r_tibia", "r_fibula"] # the fine blockers
    const ARTERIES = "arteries_rleg"
    const VEINS = "veins_rleg"

    "A surface's object name, without the `__occNNN` a repeated name gets."
    group(name) = replace(name, r"__occ\d+$" => "")
    "Foot bones XCAT does not name by side: the right ones are told apart by x."
    side_bone(name) = startswith(name, "metatarsal") || startswith(name, "phalanx")
    "The surfaces this notebook reads."
    wanted(name) = name in ENVELOPE || name in BONES || name in (ARTERIES, VEINS) || (side_bone(name) && !endswith(name, "_inner"))
end;

read_grouped_nrb(path, keep) -> Vector{NurbsSurface}

"""
    read_grouped_nrb(path, keep) -> Vector{NurbsSurface}

XCAT's grouped NRB: per patch a name, an object code line, a patch-scale line, `M :M`, `N :N`, then the knots and control
points of `Vasculature.XCAT.read_nrb`'s format. A name repeated over several patches gets `__occNNN`. Only names `keep`
accepts are parsed.
"""
function read_grouped_nrb(path::AbstractString, keep)
    L = readlines(path)
    n = length(L)
    out = V.XCAT.NurbsSurface[]
    seen = Dict{String, Int}()
    i = 1
    while i <= n - 4
        if isempty(strip(L[i])) || !endswith(rstrip(L[i + 3]), ":M") || !endswith(rstrip(L[i + 4]), ":N")
            i += 1
            continue
        end
        name = String(strip(L[i]))
        M = parse(Int, split(L[i + 3])[1])
        N = parse(Int, split(L[i + 4])[1])
        i += 5
        startswith(strip(L[i]), "U Knot") || error("$name: expected U knots at line $i")
        i += 1
        u = Float64[]
        while !startswith(strip(L[i]), "V Knot")
            push!(u, parse(Float64, strip(L[i])))
            i += 1
        end
        i += 1
        v = Float64[]
        while !startswith(strip(L[i]), "Control")
            push!(v, parse(Float64, strip(L[i])))
            i += 1
        end
        i += 1
        k = get(seen, name, 0) + 1
        seen[name] = k
        if keep(name)
            P = Array{Float64, 3}(undef, 3, N, M)
            for q in 1:(M * N)
                f = split(L[i + q - 1], (',', ' ', '\t'); keepempty = false)
                jj = mod1(q, N)
                ii = cld(q, N)
                for d in 1:3
                    P[d, jj, ii] = parse(Float64, f[d])
                end
            end
            push!(out, V.XCAT.NurbsSurface(k == 1 ? name : @sprintf("%s__occ%03d", name, k), u, v, P))
        end
        i += M * N
    end
    return out
end

leg_surfaces = cached("leg_surfaces") do
    S = read_grouped_nrb(LEG_NRB, wanted)
    # the unsided foot bones of the right leg (negative x)
    filter!(s -> !side_bone(group(s.name)) || mean(@view s.P[1, :, :]) < 0, S)
end;

"A NURBS surface's triangle mesh, finer for surfaces with more control points."
mesh(s) = V.XCAT.tessellate(s; nu = clamp(4 * size(s.P, 2), 64, 480), nv = clamp(4 * size(s.P, 3), 64, 480))

"""
The leg on a 1 mm grid: its grid, the skin envelope, the bones, the soft tissue (envelope minus bone), and the fine
blockers (femur, tibia, fibula at 0.4 mm) as `(mask, origin, h)`.
"""
LEG = cached("leg") do
    skin = [mesh(s) for s in leg_surfaces if group(s.name) in ENVELOPE]
    bones = [s for s in leg_surfaces if group(s.name) in BONES || side_bone(group(s.name))]
    around = V.grid_around(reduce(hcat, first.(skin)), 1.0; margin = 3)
    grid = V.VoxelGrid(Float64.(round.(Int, around.origin)), 1.0, around.dims)

    envelope = reduce(.|, (V.voxelize(X, T, grid) for (X, T) in skin))
    bone = reduce(.|, (V.voxelize(mesh(s)..., grid) for s in bones))

    long = [mesh(s) for s in bones if group(s.name) in LONG_BONES]
    fine_grid = V.grid_around(reduce(hcat, first.(long)), 0.4; margin = 2)
    fine = reduce(.|, (V.voxelize(X, T, fine_grid) for (X, T) in long))

    (grid = grid, envelope = envelope, bone = bone, soft = envelope .& .!bone, blockers = (fine, fine_grid.origin, 0.4))
end;

"The soft tissue's volume (mL): the tissue every growth pass and the flow are sized by."
const SOFT_ML = count(LEG.soft) / 1000

"Whether the point `p` (mm) lies in the voxels of `mask` on the leg grid."
function inside(mask, p; grid = LEG.grid)
    k = round.(Int, (p .- grid.origin) ./ grid.h) .+ 1
    return all(1 .<= k .<= size(mask)) && mask[k...]
end

"""
Node spacing of the imaged centrelines (mm): the heart's 2.8 mm × the leg's size scale (its tissue over the 148.1 mL
reference, cube root). A branch leaves its host on the host's centreline, and the collision test exempts only the host,
so nodes closer than the host's final radius would make every new branch hit the next imaged piece.
"""
const STEP = 2.8 * cbrt(SOFT_ML / 148.1)

"A centreline resampled every `step` mm by arc length."
function resample(c::V.XCAT.Centerline, step::Real)
    arc = [0.0; cumsum([sqrt(sum(abs2, c.X[:, k + 1] .- c.X[:, k])) for k in 1:(size(c.X, 2) - 1)])]
    n = max(2, round(Int, arc[end] / step) + 1)
    X = zeros(3, n)
    r = zeros(n)
    for (q, t) in enumerate(range(0, arc[end], length = n))
        k = clamp(searchsortedlast(arc, t), 1, length(arc) - 1)
        w = (t - arc[k]) / max(arc[k + 1] - arc[k], eps())
        X[:, q] .= (1 - w) .* c.X[:, k] .+ w .* c.X[:, k + 1]
        r[q] = (1 - w) * c.r[k] + w * c.r[k + 1]
    end
    return V.XCAT.Centerline(c.name, X, r)
end

imaged_tree(surfaces) -> (tree, r_root, attach_tol)

"""
    imaged_tree(surfaces) -> (tree, r_root, attach_tol)

The imaged patches of one vessel system joined into one tree: each patch's centreline (u and v swapped where the patch
runs the other way), the highest patch as the root, cut where it enters the leg's skin, resampled every `STEP` mm and
joined with the smallest attachment tolerance that joins them all. `r_root` is the root's mean radius over its first
10 mm.
"""
function imaged_tree(surfaces)
    lines = map(surfaces) do s
        c = V.XCAT.centerline(s; step = 1.0)
        # parametrized the other way round (u along the vessel): swap u and v
        median(c.r) > 15 ? V.XCAT.centerline(V.XCAT.NurbsSurface(s.name, s.v, s.u, permutedims(s.P, (1, 3, 2))); step = 1.0) : c
    end
    top = argmax([maximum(@view l.X[3, :]) for l in lines])
    c = lines[top]
    c.X[3, 1] < c.X[3, end] && (c = V.XCAT.Centerline(c.name, c.X[:, end:-1:1], reverse(c.r)))
    k0 = findfirst(k -> inside(LEG.envelope, c.X[:, k]), axes(c.X, 2))
    lines[top] = V.XCAT.Centerline(c.name, c.X[:, k0:end], c.r[k0:end])
    r_root = mean(lines[top].r[1:min(10, end)])
    pieces = [resample(l, STEP) for l in lines]
    for tol in (3.0, 5.0, 8.0, 12.0, 20.0)
        try
            return (tree = only(V.XCAT.base_forest(pieces; roots = [c.name], attach_tol = tol)), r_root, attach_tol = tol)
        catch e
            (e isa ErrorException && occursin("base_forest", e.msg)) || rethrow()
        end
    end
    return error("the imaged $(group(c.name)) patches do not join into one tree within 20 mm")
end

"""
The growth inputs of one vessel system: its imaged tree with penetrating stubs, the stub report, the root radius, and the
lumen voxels outside the tissue that vessels may cross but where no terminal ends.
"""
inlets = cached("inlets") do
    map((arteries = ARTERIES, veins = VEINS)) do name
        base = imaged_tree([s for s in leg_surfaces if group(s.name) == name])
        tree, report = V.XCAT.penetrating_stubs(
            base.tree,
            LEG.grid,
            LEG.soft;
            scale = STEP / 2.8,
            forbid = LEG.bone,
            blockers = (LEG.bone, LEG.grid.origin, LEG.grid.h),
        )
        pass = V.XCAT.vessel_tube(tree, LEG.grid; margin = 1.5) .& .!LEG.soft .& .!LEG.bone
        (tree = tree, report = report, r_root = base.r_root, attach_tol = base.attach_tol, pass = pass)
    end
end;

"The leg's front (x–z) and side (y–z) silhouettes, with segments drawn over them: one `linesegments!` per view."
function leg_figure(trees; title = "", size = (900, 900))
    fig = Mke.Figure(; size)
    g = LEG.grid
    xs = [g.origin[d] .+ (0:(g.dims[d] - 1)) .* g.h for d in 1:3]
    for (column, (view, label)) in enumerate(((1, "front (x–z)"), (2, "side (y–z)")))
        ax = Mke.Axis(fig[1, column]; title = label, aspect = Mke.DataAspect(), xlabel = view == 1 ? "x (mm)" : "y (mm)", ylabel = "z (mm)")
        silhouette = dropdims(any(LEG.soft; dims = 3 - view); dims = 3 - view)
        Mke.heatmap!(ax, xs[view], xs[3], silhouette; colormap = [:white, :gray88])
        for (segments, colour, scale) in trees
            a, b, r = segments
            points = Mke.Point2f[]
            for j in eachindex(r)
                push!(points, Mke.Point2f(a[view, j], a[3, j]), Mke.Point2f(b[view, j], b[3, j]))
            end
            widths = repeat(max.(0.4, scale .* r); inner = 2)
            Mke.linesegments!(ax, points; color = colour, linewidth = widths)
        end
        Mke.hidedecorations!(ax; label = false, ticklabels = false, ticks = false)
    end
    Mke.Label(fig[0, :], title; fontsize = 18, font = :bold)
    return fig
end

"An inlet tree's segments as `(a, b, r)`: each node with its parent."
function base_segments(tree)
    k = 2:length(tree.parent)
    return (tree.X[:, tree.parent[k]], tree.X[:, k], tree.r[k])
end

begin
    "Radius (mm) of a terminal arteriole, the size growth must reach."
    const TERMINAL_RADIUS_MM = 0.01
    "Murray's exponent: the leg's conduits need it (γ = 2.33 starves the beds)."
    const GAMMA = 3.0
    "The terminal count Murray's law gives for the imaged femoral artery and a 10 µm terminal."
    const N_TOTAL = round(Int, (inlets.arteries.r_root / TERMINAL_RADIUS_MM)^GAMMA)
    "Tissue (mL) one pass-1 terminal feeds: one pass-2 patch."
    const PATCH_ML = 10.0
    "Pass-1 terminals: one per patch."
    const N_PASS1 = round(Int, SOFT_ML / PATCH_ML)
end;

"""
The leg's growth preset: `HEART_GROWTH`'s generic values with the leg's own tissue, imaged root and Murray's exponent.
A draft: once the leg is built it belongs in `src/organs.jl` (as `LIMB_GROWTH`), with its sources in
`test/reference/`.
"""
limb_growth(nterm, tissue_ml, root_radius) = V.GrowthSpec(
    V.HEART_GROWTH;
    nterm,
    flow_ml_min = nothing,
    tissue_ml,
    root_caps = [Float64(root_radius)],
    root_fixed = true,
    symmetry_floor = Tuple{Float64, Float64}[],
    gamma = GAMMA,
    device = DEVICE,
)

"""
Grow one tree into `tissue` (a mask on `grid`): `pass` voxels may be crossed but take no terminal. Returns the grown
tree (`V.vtree`) and the seconds it took.
"""
function grow_tree(spec, inlet, tissue, pass, grid, blockers, name)
    dir = mktempdir()
    allow = zeros(UInt16, size(tissue))
    allow[pass] .= V.NOTERM
    domain = V.write_vol(joinpath(dir, "$name.vol"), tissue .| pass, Tuple(round.(Int, grid.origin ./ grid.h)))
    gxl = V.XCAT.write_base_gxl(joinpath(dir, "$name.xml"), inlet)
    seconds = @elapsed engine = V.grow(spec, domain, [gxl]; blockers, territories = allow)
    return (tree = V.vtree(engine), seconds = seconds)
end

pass1 = cached("pass1") do
    map(inlets) do inlet
        spec = limb_growth(N_PASS1, SOFT_ML, inlet.r_root)
        grow_tree(spec, inlet.tree, LEG.soft, inlet.pass, LEG.grid, LEG.blockers, "pass1")
    end
end;

"The terminal segments of a grown tree: the segments without children."
function terminals(tree)
    has_child = falses(length(tree.r))
    foreach(p -> p > 0 && (has_child[p] = true), tree.parent)
    return findall(!, has_child)
end

flood(mask, grid, points) -> (label, voxels, box)

"""
    flood(mask, grid, points) -> (label, voxels, box)

Each voxel of `mask` labelled with the nearest of `points` (mm) through the mask (a breadth-first flood from all of them
at once, six neighbours per step), the voxel count per label, and each label's bounding box (`CartesianIndices`).
"""
function flood(mask, grid, points)
    label = zeros(Int32, size(mask))
    queue = CartesianIndex{3}[]
    for (k, p) in enumerate(eachcol(points))
        c = CartesianIndex(Tuple(clamp.(round.(Int, (p .- grid.origin) ./ grid.h) .+ 1, 1, size(mask))))
        mask[c] && label[c] == 0 || continue
        label[c] = k
        push!(queue, c)
    end
    head = 1
    while head <= length(queue)
        c = queue[head]
        head += 1
        for s in SIX_STEPS
            n = c + s
            checkbounds(Bool, mask, n) && mask[n] && label[n] == 0 || continue
            label[n] = label[c]
            push!(queue, n)
        end
    end
    n = size(points, 2)
    voxels = zeros(Int, n)
    lo = fill(CartesianIndex(size(mask)), n)
    hi = fill(CartesianIndex(1, 1, 1), n)
    for c in CartesianIndices(label)
        k = label[c]
        k == 0 && continue
        voxels[k] += 1
        lo[k] = min(lo[k], c)
        hi[k] = max(hi[k], c)
    end
    return (label = label, voxels = voxels, box = [lo[k]:hi[k] for k in 1:n])
end

"The arterial patches: each pass-1 arterial terminal's segment and its tissue."
PATCHES = cached("patches") do
    t = pass1.arteries.tree
    ends = terminals(t)
    f = flood(LEG.soft, LEG.grid, t.b[:, ends])
    (segment = ends, f...)
end;

"The segments of a tree (its `parent` vector) in breadth-first order from the root, and each segment's children."
function walk(parent)
    children = [Int[] for _ in eachindex(parent)]
    order = Int[]
    for j in eachindex(parent)
        parent[j] == 0 ? push!(order, j) : push!(children[parent[j]], j)
    end
    for j in order # `order` grows as it is read
        append!(order, children[j])
    end
    return order, children
end

downstream(parent, weight) -> (order, children, n)

"""
    downstream(parent, weight) -> (order, children, n)

The terminals below each segment: a terminal `j` counts `weight(j)` (the final terminals it stands for), every other
segment the sum over its children.
"""
function downstream(parent, weight)
    order, children = walk(parent)
    n = zeros(length(parent))
    for j in reverse(order)
        n[j] = isempty(children[j]) ? weight(j) : sum(n[c] for c in children[j])
    end
    return order, children, n
end

"Each patch's final terminal count: its share of `N_TOTAL` by tissue volume, so every terminal carries the same flow."
const PATCH_TERMINALS = round.(Int, N_TOTAL .* PATCHES.voxels ./ sum(PATCHES.voxels))

"""
The pass-1 arteries at their Murray radii: each segment's final terminal count below it (`below`) and its radius
`TERMINAL_RADIUS_MM · below^(1/γ)` (mm). A pass-1 terminal counts its patch's terminals (at least one).
"""
PASS1_MURRAY = let t = pass1.arteries.tree, patch_of = Dict(j => k for (k, j) in enumerate(PATCHES.segment))
    _, _, below = downstream(t.parent, j -> max(PATCH_TERMINALS[patch_of[j]], 1))
    (below = below, radius = TERMINAL_RADIUS_MM .* below .^ (1 / GAMMA))
end;

grow_patch(k) -> (tree, seconds, nterm, ml, root_radius, segment)

"""
    grow_patch(k) -> (tree, seconds, nterm, ml, root_radius, segment)

Pass 2 for patch `k`: the patch cropped to its box (and the pass-1 segment that feeds it), grown from that segment at
its Murray radius with the patch's share of `N_TOTAL` terminals. The neighbouring patches' tissue within
`PATCH_MARGIN_MM` and the segment's lumen may be crossed but take no terminal, so the patch's terminals reach its border
and neighbouring patches meet without a seam.
"""
function grow_patch(k; tree = pass1.arteries.tree)
    j = PATCHES.segment[k]
    a = Float64.(tree.a[:, j])
    b = Float64.(tree.b[:, j])
    r = Float64(PASS1_MURRAY.radius[j])
    g = LEG.grid
    # the patch's box, grown to hold its feeding segment and the margin
    pad = ceil(Int, r / g.h) + ceil(Int, PATCH_MARGIN_MM / g.h) + 2
    ka = round.(Int, (a .- g.origin) ./ g.h) .+ 1
    lo = max.(Tuple(min(first(PATCHES.box[k]), CartesianIndex(Tuple(ka)))) .- pad, 1)
    hi = min.(Tuple(max(last(PATCHES.box[k]), CartesianIndex(Tuple(ka)))) .+ pad, g.dims)
    box = CartesianIndices(Tuple(UnitRange.(lo, hi)))
    grid = V.VoxelGrid(g.origin .+ (Tuple(first(box)) .- 1) .* g.h, g.h, size(box))

    tissue = PATCHES.label[box] .== k
    soft = LEG.soft[box]
    margin = dilate(tissue, round(Int, PATCH_MARGIN_MM / g.h)) .& soft .& .!tissue
    inlet = V.XCAT.BaseTree("patch_$k", [0, 1], hcat(a, b), [r, r], ["pass1", "pass1"])
    pass = (V.XCAT.vessel_tube(inlet, grid; margin = 1.5) .& .!LEG.bone[box] .| margin) .& .!tissue
    ml = count(tissue) / 1000
    nterm = PATCH_TERMINALS[k]
    grown = grow_tree(limb_growth(nterm, ml, r), inlet, tissue, pass, grid, LEG.blockers, "patch_$k")
    return (tree = grown.tree, seconds = grown.seconds, nterm = nterm, ml = ml, root_radius = r, segment = j)
end

begin
    "Neighbouring tissue (mm) a patch may grow through but not end in, so patches meet without a seam."
    const PATCH_MARGIN_MM = 3.0

    "`mask` grown `steps` voxels through its six face neighbours."
    function dilate(mask, steps)
        out = copy(mask)
        for _ in 1:steps
            grown = copy(out)
            for c in CartesianIndices(out)
                out[c] || continue
                for s in SIX_STEPS
                    n = c + s
                    checkbounds(Bool, out, n) && (grown[n] = true)
                end
            end
            out = grown
        end
        return out
    end
end

"""
The region grown here: the pass-1 branch with 8 to 16 patches below it whose end lies nearest 30% of the way up the leg
(the calf), and the patches below it.
"""
REGION = let t = pass1.arteries.tree, z = t.b[3, PATCHES.segment]
    _, _, patches = downstream(t.parent, j -> 1)
    height = minimum(z) + 0.3 * (maximum(z) - minimum(z))
    candidates = findall(n -> 8 <= n <= 16, patches)
    branch = candidates[argmin(abs.(t.b[3, candidates] .- height))]
    function under(j)
        while j != 0 && j != branch
            j = t.parent[j]
        end
        return j == branch
    end
    (branch = branch, patches = [k for (k, j) in enumerate(PATCHES.segment) if under(j) && PATCHES.voxels[k] > 0])
end

"Pass 2 for every patch of the region, each cached on its own."
pass2 = [
    cached("patch_$k") do
        grow_patch(k)
    end for k in REGION.patches
];

stitch(tree, radius, grown) -> (parent, a, b, r, pass, from)

"""
    stitch(tree, radius, grown) -> (parent, a, b, r, pass, from)

One tree from the pass-1 `tree` (at its Murray `radius`) and the pass-2 trees `grown`: each grown tree replaces the
pass-1 terminal segment it grew from and hangs from that segment's parent. `pass` is each segment's pass (1 or 2) and
`from` a pass-1 segment's index in `tree` (0 for a grown one).
"""
function stitch(tree, radius, grown)
    replaced = Set(p.segment for p in grown)
    keep = [j for j in eachindex(tree.r) if !(j in replaced)]
    index = zeros(Int, length(tree.r))
    index[keep] = eachindex(keep)
    parent = Int[p == 0 ? 0 : index[p] for p in tree.parent[keep]]
    a = [tree.a[:, keep]]
    b = [tree.b[:, keep]]
    r = [Float32.(radius[keep])]
    pass = [fill(1, length(keep))]
    from = [keep]
    n = length(keep)
    for p in grown
        t = p.tree
        up = index[tree.parent[p.segment]]
        append!(parent, [q == 0 ? up : q + n for q in t.parent])
        push!(a, t.a)
        push!(b, t.b)
        push!(r, t.r)
        push!(pass, fill(2, length(t.r)))
        push!(from, zeros(Int, length(t.r)))
        n += length(t.r)
    end
    return (
        parent = parent, a = reduce(hcat, a), b = reduce(hcat, b), r = reduce(vcat, r),
        pass = reduce(vcat, pass), from = reduce(vcat, from),
    )
end

"The leg's arterial tree with the region complete to 10 µm: pass 1 at its Murray radii, the region's patches grafted in."
stitched = stitch(pass1.arteries.tree, PASS1_MURRAY.radius, pass2);

"The final terminals a terminal of `stitched` stands for: one if it was grown in pass 2, its patch's count if not."
terminal_weight = let patch_of = Dict(j => k for (k, j) in enumerate(PATCHES.segment))
    j -> stitched.pass[j] == 2 ? 1.0 : Float64(max(PATCH_TERMINALS[patch_of[stitched.from[j]]], 1))
end

pressure_drop(tree, q, weight) -> (terminals, drop)

"""
    pressure_drop(tree, q, weight) -> (terminals, drop)

Pressure lost (mmHg) from the root to each terminal when the tree carries `q` mL/min, each terminal taking a share in
proportion to `weight(j)`, the final terminals it stands for. Every vessel at its radius (fully dilated), Poiseuille
with the Pries–Secomb in vivo viscosity at the discharge hematocrit.
"""
function pressure_drop(tree, q, weight; phys = LIMB_FLOW)
    order, children, n = downstream(tree.parent, weight)
    unit = q / 60 / n[order[1]] # mL/s per final terminal
    P = zeros(length(tree.r))
    for j in order
        L = Float32(sqrt(sum(abs2, tree.b[:, j] .- tree.a[:, j])))
        r = tree.r[j]
        drop = V.Flow.poiseuille(L, r, phys.plasma_viscosity_mPas * V.Flow.mu_rel(2000 * r, phys.hd)) * unit * n[j]
        P[j] = (tree.parent[j] == 0 ? 0.0 : P[tree.parent[j]]) + drop
    end
    ends = findall(isempty, children)
    return (terminals = ends, drop = P[ends])
end

"""
The leg's flow physiology: the fields that differ from the heart's, each with its source. It moves into `src/organs.jl`
as a preset once its sources are in `test/reference/`, as every preset's are. The heart-only fields take the leg's
no-change values: no measured blood-volume rise under adenosine (`fb_dilated = fb`) and a constant PS.
"""
const LIMB_FLOW = V.FlowSpec(
    V.HEART_FLOW;
    p_out = 11.0f0, # mmHg. Common femoral vein, supine rest, catheter (Pearson 2011)
    fb = 0.037f0, # mL/mL. Muscle 3.3 mL/100 g by PET (Raitakari 1995), fat and skin (Leggett 1991), + conduits
    fb_dilated = 0.037f0, # no measured rise under adenosine
    share_art = 0.15f0, # systemic arteries : veins = 13 : 64 (Guyton & Hall), applied to 1 - share_cap
    share_cap = 0.11f0, # muscle capillary volume 0.48 of muscle blood (Gidlöf 1988; Raitakari 1995)
    share_ven = 0.74f0, # the remainder
    cap_length_mm = 0.44f0, # functional capillary path 438 µm (Sarelius 1986)
    cap_diameter_um = 5.3f0, # human quadriceps capillary lumen (Gidlöf 1988)
    ve = 0.11f0, # muscle 0.07, fat 0.14, skin 0.46 by composition (Goodall 2020; Thompson 2005; Levitt 2003)
    ps = 0.014f0, # forearm EDTA PS (Keiding 1988) scaled to iodinated contrast (Levitt 2002)
    ps_slope = 0.0f0,
    rho = 1.02f0, # ICRU-44 muscle, adipose, skin by composition (Hubbell & Seltzer 1995)
    r_total_dilated = 300.0f0, # 90 mmHg / PET flow under adenosine (Heinonen 2010, 2014)
    d_active = 300.0f0, # muscle feed arteries regulate (Williams & Segal 1993; Ives 2012)
    oef = 0.315f0, # supine femoral venous SO2 67.1% vs arterial 97.9% (Pearson 2011)
    mvo2 = 0.00146f0, # 2.6 mL/min/100 mL (Koch Esteves 2023) × O2 content × oef / rho
)

"""
The injection, the same as the heart's and lungs' (Hadjiabdolhamid et al. 2024, Shojazadeh et al. 2026): 0.5 mL/kg of
contrast at 5 mL/s, then 0.5 mL/kg of a 30:70 contrast:saline chaser at the same rate, for the body model's 70 kg.
"""
const INJECTION = let kg = 70, conc = 350
    V.Body.Injection(; phases = [(0.5 * kg, 5, conc), (0.5 * kg, 5, 0.3 * conc)])
end