diff --git a/README.md b/README.md
index 3312cd1..170dd53 100644
--- a/README.md
+++ b/README.md
@@ -22,7 +22,7 @@ HASA-PAM is a Julia project for GPU-accelerated transcranial passive acoustic ma
 
 The repository implements a reconstruction pipeline for localising acoustic sources inside the skull. The motivating application is non-invasive targeted drug delivery in the brain: if ultrasound-controllable drug carriers can be aggregated and uncaged at a chosen target, a future closed-loop system would also need feedback about whether drug-carrier or bubble clusters formed and where they are located.
 
-This project addresses that feedback problem with transcranial PAM. Receiver recordings are time-reversed through an acoustic model, and the reconstructed intensity field is used to estimate source locations. The implementation compares a fast homogeneous angular-spectrum baseline with a skull-aware Heterogeneous Angular Spectrum Approach (HASA) reconstruction that uses a CT-derived speed-of-sound model. The current benchmark is simulation-based and uses synthetic source emissions.
+This project addresses that feedback problem with transcranial PAM. While the primary focus is PAM, transcranial **focusing** (geometric and HASA delay estimation with optional k-Wave verification) is also fully implemented and accessible via CLI for both 2D and 3D domains. Receiver recordings are time-reversed through an acoustic model, and the reconstructed intensity field is used to estimate source locations. The implementation compares a fast homogeneous angular-spectrum baseline with a skull-aware Heterogeneous Angular Spectrum Approach (HASA) reconstruction that uses a CT-derived speed-of-sound model. The current benchmark is simulation-based and uses synthetic source emissions.
 
 In the reported synthetic 3D benchmark, skull correction reduced mean localisation error from **1.77 mm** to **0.49 mm**, while the GPU implementation reduced corrected reconstruction march time from **13.28 s** to **0.23 s**.
 
@@ -213,7 +213,10 @@ For more detailed setup instructions, workflow notes, and CLI parameter referenc
 
 ### Command line entry points
 
-The repository is meant to be used primarily from the command line. The maintained PAM entry point is `scripts/run_pam.jl`; focusing scripts are kept for the earlier transcranial focusing workflow.
+The repository exposes two primary command line entry points:
+
+- **`scripts/run_pam.jl`** — the main PAM workflow: simulate emissions, reconstruct source distributions, evaluate accuracy.
+- **`scripts/run_focus.jl`** — transcranial focusing: compute geometric or HASA transmit delays, optionally verify with k-Wave, for 2D and 3D domains.
 
 ### Setup
 
@@ -256,10 +259,28 @@ julia --project=. scripts/run_pam.jl \
   --recon-use-gpu=true
 ```
 
-### Quick focusing example
+### Quick focusing examples
+
+2D skull-backed HASA focus:
 
 ```bash
-julia --project=. scripts/run_focus.jl --estimator=hasa --medium=skull_in_water --slice-index=250
+julia --project=. scripts/run_focus.jl \
+  --estimator=hasa \
+  --aberrator=skull \
+  --slice-index=250
+```
+
+3D HASA focus in water (no k-Wave, fast smoke test):
+
+```bash
+julia --project=. scripts/run_focus.jl \
+  --dimension=3 \
+  --estimator=hasa \
+  --aberrator=water \
+  --focus-mm=60:0:0 \
+  --receiver-aperture-y-mm=60 \
+  --receiver-aperture-z-mm=60 \
+  --run-kwave=false
 ```
 
 ### Tests
diff --git a/docs/generate_cli_reference.jl b/docs/generate_cli_reference.jl
index 047cccb..a5e1da8 100644
--- a/docs/generate_cli_reference.jl
+++ b/docs/generate_cli_reference.jl
@@ -68,6 +68,11 @@ function generate_cli_reference(path=joinpath(@__DIR__, "src", "cli", "parameter
             println(io)
             _write_option_table(io, category_options)
             println(io)
+            if category == "Medium"
+                println(io, "!!! warning \"Scan-specific defaults\"")
+                println(io, "    The defaults for `--slice-index`, `--skull-transducer-distance-mm`, and `--hu-bone-thr` are calibrated to the specific CT scan used during development. Using a different scan without adjusting these values will likely produce an incorrect skull medium.")
+                println(io)
+            end
         end
     end
 
diff --git a/docs/src/cli/overview.md b/docs/src/cli/overview.md
index 4216a54..98967ce 100644
--- a/docs/src/cli/overview.md
+++ b/docs/src/cli/overview.md
@@ -1,12 +1,12 @@
 # CLI Overview
 
-The maintained entry point, and the main objective of this repository, is PAM:
+The primary entry point for passive acoustic mapping:
 
 ```bash
 julia --project=. scripts/run_pam.jl --option=value
 ```
 
-The focusing scripts are also available for reference to the earlier focusing workflow:
+Transcranial focusing (geometric and HASA delay estimation, 2D and 3D):
 
 ```bash
 julia --project=. scripts/run_focus.jl --option=value
diff --git a/docs/src/cli/parameters.md b/docs/src/cli/parameters.md
index ec10884..7a10ded 100644
--- a/docs/src/cli/parameters.md
+++ b/docs/src/cli/parameters.md
@@ -135,6 +135,9 @@ For practical guidance, start with [Running PAM](@ref) before tuning individual
 | `--skull-transducer-distance-mm` | `30` | mm | skull |   | Distance from the receiver/transducer plane to the outer skull surface. |
 | `--hu-bone-thr` | `200` | HU | skull |   | Hounsfield-unit threshold used to identify bone in CT data. |
 
+!!! warning "Scan-specific defaults"
+    The defaults for `--slice-index`, `--skull-transducer-distance-mm`, and `--hu-bone-thr` are calibrated to the specific CT scan used during development. Using a different scan without adjusting these values will likely produce an incorrect skull medium.
+
 ## Simulation
 
 | Option | Default | Value | Applies to | Choices | Description |
@@ -156,3 +159,4 @@ For practical guidance, start with [Running PAM](@ref) before tuning individual
 | `--recon-min-window-energy-ratio` | `0.001` | ratio | windowed |   | Skips windows whose energy is below this fraction of the maximum window energy. |
 | `--recon-progress` | `false` | bool | PAM |   | Prints reconstruction progress updates. |
 | `--window-batch` | `1` | integer | windowed GPU |   | Number of reconstruction windows batched together on the GPU. |
+
diff --git a/docs/src/cli/run-focus.md b/docs/src/cli/run-focus.md
index 5f78203..bbe1e64 100644
--- a/docs/src/cli/run-focus.md
+++ b/docs/src/cli/run-focus.md
@@ -1,6 +1,6 @@
 # Running Focus
 
-The focusing scripts are retained only as reference material from the earlier transcranial focusing workflow. The main objective of this repository is PAM; these scripts are useful for comparing geometric and HASA focusing corrections through water or a CT-backed skull medium.
+`scripts/run_focus.jl` computes transcranial transmit delays (geometric or HASA) for 2D and 3D domains, optionally verifying focus quality with k-Wave. It supports both water and CT-backed skull media. `scripts/compare_focus_estimators.jl` runs both estimators side-by-side and writes a comparison figure.
 
 ## Single Focusing Case
 
@@ -9,7 +9,7 @@ Default skull-backed HASA case:
 ```bash
 julia --project=. scripts/run_focus.jl \
   --estimator=hasa \
-  --medium=skull_in_water \
+  --aberrator=skull \
   --slice-index=250
 ```
 
@@ -17,23 +17,53 @@ Centered target at 60 mm below the transducer:
 
 ```bash
 julia --project=. scripts/run_focus.jl \
+  --dimension=2 \
   --estimator=hasa \
-  --medium=skull_in_water \
+  --aberrator=skull \
   --placement=fixed_transducer \
   --slice-index=250 \
-  --lateral-cm=0.0 \
-  --focal-cm=6.0
+  --focus-mm=60:0
+```
+
+3D water-backed HASA setup dry enough for checking the focusing delays without k-Wave:
+
+```bash
+julia --project=. scripts/run_focus.jl \
+  --dimension=3 \
+  --estimator=hasa \
+  --aberrator=water \
+  --focus-mm=60:0:0 \
+  --receiver-aperture-y-mm=60 \
+  --receiver-aperture-z-mm=60 \
+  --run-kwave=false
+```
+
+3D skull-backed k-Wave outward propagation without inward verification:
+
+```bash
+julia --project=. scripts/run_focus.jl \
+  --dimension=3 \
+  --estimator=hasa \
+  --aberrator=skull \
+  --outward-propagation=kwave \
+  --verify-inward-kwave=false \
+  --focus-mm=30:0:0 \
+  --slice-index=250 \
+  --skull-transducer-distance-mm=5 \
+  --receiver-aperture-y-mm=20 \
+  --receiver-aperture-z-mm=20
 ```
 
 Target 30 mm below the inner skull:
 
 ```bash
 julia --project=. scripts/run_focus.jl \
+  --dimension=2 \
   --estimator=hasa \
-  --medium=skull_in_water \
+  --aberrator=skull \
   --placement=fixed_focus_depth \
   --slice-index=250 \
-  --focal-cm=6.0 \
+  --focus-mm=60:0 \
   --focus-depth-from-inner-skull-mm=30
 ```
 
@@ -50,13 +80,23 @@ julia --project=. scripts/compare_focus_estimators.jl \
 - `--ct-path`: DICOM folder for CT-backed skull runs.
 - `--slice-index`: CT slice used for the 2D focusing medium.
 - `--frequency-mhz`: transmit frequency.
-- `--focal-cm`: transducer-to-focus distance for fixed-transducer placement, or transducer distance to the resolved target for fixed-depth placement.
-- `--lateral-cm`: lateral target offset.
-- `--aperture-cm`: transducer aperture width.
+- `--dimension`: `2` or `3`.
+- `--focus-mm`: target coordinates. Use `depth:lateral` in 2D and `depth:y:z` in 3D.
+- `--receiver-aperture-mm`: transducer aperture width. In 3D, `--receiver-aperture-y-mm` and `--receiver-aperture-z-mm` override each axis independently.
 - `--estimator`: `geometric` or `hasa`.
-- `--medium`: `water` or `skull_in_water`.
-- `--placement`: `auto`, `fixed_transducer`, or `fixed_focus_depth`.
-- `--focus-depth-from-inner-skull-mm`: target depth below the inner skull for fixed-depth placement.
+- `--aberrator`: `none`, `water`, or `skull`.
+- `--placement`: `auto`, `fixed_transducer`, or `fixed_focus_depth` (2D only).
+- `--focus-depth-from-inner-skull-mm`: target depth below the inner skull for fixed-depth placement (2D only).
+- `--skull-transducer-distance-mm`: gap between the outer skull surface and the transducer plane.
+- `--outward-propagation`: `angular_spectrum` or `kwave`.
+- `--verify-inward-kwave`: set `false` to skip the return pressure simulation after computing transmit delays.
+- `--run-kwave`: legacy shorthand; setting it to `false` also skips inward verification.
+- `--kwave-use-gpu`: set `true` to run k-Wave on GPU.
+- `--transverse-mm`: transverse grid extent (mm). In 3D, `--transverse-y-mm` and `--transverse-z-mm` override each axis independently.
+- `--dx-mm`, `--dz-mm`: lateral and axial grid spacing. `--dy-mm` overrides lateral spacing for the y axis in 3D.
+- `--axial-mm`: axial domain extent; defaults to `auto` (derived from focus depth and padding).
+- `--axial-padding`: multiplicative padding applied to the focus depth when computing the axial domain size (default `1.5`).
 - `--out-dir`: override the automatically generated output directory.
 
 `scripts/run_focus.jl` writes `summary.json`, `result.jld2`, and `pressure.png`. `scripts/compare_focus_estimators.jl` writes `summary.json` and `comparison.png`.
+Focus runs also write `sensor_timings.png`, which visualizes the computed transmit delay and normalized amplitude at the active sensor/aperture elements.
diff --git a/docs/src/getting-started.md b/docs/src/getting-started.md
index c0cb960..1d52d72 100644
--- a/docs/src/getting-started.md
+++ b/docs/src/getting-started.md
@@ -22,6 +22,9 @@ The private CT data used during development is not distributed with this reposit
 
 Override this with `--ct-path=/path/to/dicom-folder` for skull-backed runs. Homogeneous runs with `--aberrator=none` do not need CT data.
 
+!!! warning "Scan-specific defaults"
+    Parameters such as `--slice-index`, `--skull-transducer-distance-mm`, and `--hu-bone-thr` are set to values that work with the specific CT scan used during development of this project. They will likely need to be adjusted for a different scan to ensure the skull medium is constructed correctly.
+
 ## First PAM Run
 
 A small homogeneous point-source run is the fastest way to confirm the CLI works:
diff --git a/docs/src/index.md b/docs/src/index.md
index 6869321..d36ea21 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -1,7 +1,7 @@
 # TranscranialFUS Documentation
 
-`TranscranialFUS` is a Julia project whose main objective is passive acoustic mapping (PAM) for transcranial ultrasound. The maintained workflow is `scripts/run_pam.jl`; the older focusing scripts are documented only as reference material.
+`TranscranialFUS` is a Julia project whose primary objective is passive acoustic mapping (PAM) for transcranial ultrasound. The main workflow is `scripts/run_pam.jl`. Transcranial **focusing** (geometric and HASA delay estimation, with optional k-Wave verification) is also fully implemented and available via `scripts/run_focus.jl` for both 2D and 3D domains.
 
-Start with [Getting Started](@ref), then use [Running PAM](@ref) for examples and [PAM CLI Parameters](@ref) for the generated option reference. Supporting pages cover CLI conventions, workflow, outputs, benchmark, troubleshooting, and focusing scripts kept from the earlier workflow.
+Start with [Getting Started](@ref), then use [Running PAM](@ref) for examples and [PAM CLI Parameters](@ref) for the generated option reference. See [Running Focus](@ref) for the focusing CLI. Supporting pages cover CLI conventions, workflow, outputs, benchmark, and troubleshooting.
 
 The code is research-oriented and assumes access to local CT data for skull-backed examples. Homogeneous water examples do not require CT data.
diff --git a/scripts/README.md b/scripts/README.md
new file mode 100644
index 0000000..e8a6041
--- /dev/null
+++ b/scripts/README.md
@@ -0,0 +1,119 @@
+# Script Examples
+
+## 3D Focusing
+
+Fast delay/HASA smoke test without running k-Wave:
+
+```bash
+julia --project=. scripts/run_focus.jl \
+  --dimension=3 \
+  --aberrator=water \
+  --estimator=hasa \
+  --focus-mm=12:0:0 \
+  --dx-mm=1 \
+  --dy-mm=1 \
+  --dz-mm=1 \
+  --transverse-mm=21 \
+  --receiver-aperture-mm=11 \
+  --run-kwave=false \
+  --out-dir=/private/tmp/focus3d_smoke
+```
+
+Use k-Wave for the outward focus-to-aperture propagation, then skip the inward
+verification simulation:
+
+```bash
+julia --project=. scripts/run_focus.jl \
+  --dimension=3 \
+  --aberrator=skull \
+  --estimator=hasa \
+  --outward-propagation=kwave \
+  --verify-inward-kwave=false \
+  --focus-mm=30:0:0 \
+  --slice-index=250 \
+  --skull-transducer-distance-mm=5 \
+  --receiver-aperture-y-mm=20 \
+  --receiver-aperture-z-mm=20
+```
+
+3D water focusing with k-Wave enabled:
+
+```bash
+julia --project=. scripts/run_focus.jl \
+  --dimension=3 \
+  --aberrator=water \
+  --estimator=geometric \
+  --focus-mm=60:0:0 \
+  --transverse-y-mm=80 \
+  --transverse-z-mm=80 \
+  --receiver-aperture-y-mm=60 \
+  --receiver-aperture-z-mm=60
+```
+
+Kwave outward, no verification with mac skull path
+```bash
+julia --project=. scripts/run_focus.jl \
+  --ct-path="$HOME/INI_code/Ultrasound/DIRU_20240404_human_skull_phase_correction_1_2_(skull_Normal)/DICOM/PAT_0000/STD_0000/SER_0002/OBJ_0001" \
+  --dimension=3 \
+  --aberrator=skull \
+  --focus-mm=30:0:0 \
+  --outward-propagation=kwave \
+  --verify-inward-kwave=false \
+  --slice-index=250 \
+  --skull-transducer-distance-mm=5 \
+  --receiver-aperture-y-mm=20 \
+  --receiver-aperture-z-mm=20 \
+  --transverse-y-mm=30 \
+  --transverse-z-mm=30 \
+  --axial-padding=1.2
+```
+
+3D CT-backed skull focusing:
+
+```bash
+julia --project=. scripts/run_focus.jl \
+  --ct-path="$HOME/INI_code/Ultrasound/DIRU_20240404_human_skull_phase_correction_1_2_(skull_Normal)/DICOM/PAT_0000/STD_0000/SER_0002/OBJ_0001" \
+  --dimension=3 \
+  --aberrator=skull \
+  --estimator=hasa \
+  --focus-mm=30:0:0 \
+  --slice-index=250 \
+  --skull-transducer-distance-mm=5 \
+  --receiver-aperture-y-mm=20 \
+  --receiver-aperture-z-mm=20
+```
+
+## 2D Focusing
+
+Fast 2D smoke test without running k-Wave:
+
+```bash
+julia --project=. scripts/run_focus.jl \
+  --dimension=2 \
+  --aberrator=water \
+  --estimator=hasa \
+  --focus-mm=12:0 \
+  --dx-mm=1 \
+  --dz-mm=1 \
+  --transverse-mm=21 \
+  --receiver-aperture-mm=11 \
+  --run-kwave=false \
+  --out-dir=/private/tmp/focus2d_smoke
+```
+
+2D CT-backed skull focusing at fixed depth below the inner skull:
+
+```bash
+julia --project=. scripts/run_focus.jl \
+  --dimension=2 \
+  --aberrator=skull \
+  --estimator=hasa \
+  --placement=fixed_focus_depth \
+  --focus-mm=60:0 \
+  --focus-depth-from-inner-skull-mm=30 \
+  --slice-index=250
+```
+
+`--focus-mm` is `depth:lateral` in 2D and `depth:y:z` in 3D. Outputs are written under `outputs/` unless `--out-dir` is provided.
+`--outward-propagation=kwave` estimates transmit delays by simulating from the focus to the aperture; the default `angular_spectrum` uses the geometric/HASA path. `--verify-inward-kwave=false` skips the return pressure-field verification step.
+Every focus run writes `sensor_timings.png`, showing the per-element transmit delay and normalized amplitude at the active sensor/aperture.
diff --git a/scripts/run_focus.jl b/scripts/run_focus.jl
index 158750a..c211c86 100644
--- a/scripts/run_focus.jl
+++ b/scripts/run_focus.jl
@@ -1,54 +1,136 @@
 #!/usr/bin/env julia
 
-using Pkg
-Pkg.activate(joinpath(@__DIR__, ".."))
+if abspath(PROGRAM_FILE) == @__FILE__
+    import Pkg
+    Pkg.activate(joinpath(@__DIR__, ".."))
+end
 
 using Dates
 using CairoMakie
 using JLD2
 using JSON3
 using TranscranialFUS
+import TranscranialFUS: parse_bool, parse_dimension, parse_receiver_aperture_mm
 
-function parse_cli(args)
-    opts = Dict{String, String}(
-        "ct-path" => DEFAULT_CT_PATH,
-        "slice-index" => "250",
-        "frequency-mhz" => "0.5",
-        "focal-cm" => "6.0",
-        "lateral-cm" => "0.0",
-        "aperture-cm" => "10.0",
-        "estimator" => "hasa",
-        "medium" => "skull_in_water",
-        "placement" => "auto",
-    )
+const FOCUS_DEFAULTS = Dict{String, String}(
+    "dimension" => "2",
+    "ct-path" => DEFAULT_CT_PATH,
+    "slice-index" => "250",
+    "frequency-mhz" => "0.5",
+    "focus-mm" => "",
+    "estimator" => "hasa",
+    "aberrator" => "skull",
+    "placement" => "auto",
+    "focus-depth-from-inner-skull-mm" => "",
+    "skull-transducer-distance-mm" => "5",
+    "hu-bone-thr" => "200",
+    "axial-mm" => "auto",
+    "transverse-mm" => "120",
+    "transverse-y-mm" => "",
+    "transverse-z-mm" => "",
+    "dx-mm" => "0.2",
+    "dy-mm" => "",
+    "dz-mm" => "0.2",
+    "receiver-aperture-mm" => "100",
+    "receiver-aperture-y-mm" => "",
+    "receiver-aperture-z-mm" => "",
+    "axial-padding" => "1.5",
+    "t-max-us" => "220",
+    "dt-ns" => "40",
+    "record" => "p_rms",
+    "kwave-use-gpu" => "false",
+    "outward-propagation" => "angular_spectrum",
+    "verify-inward-kwave" => "true",
+    "run-kwave" => "true",
+    "out-dir" => "",
+)
 
+function parse_focus_cli(args)
+    opts = copy(FOCUS_DEFAULTS)
+    provided = Set{String}()
     for arg in args
         startswith(arg, "--") || error("Unsupported argument format: $arg")
         parts = split(arg[3:end], "="; limit=2)
         length(parts) == 2 || error("Arguments must use --name=value, got: $arg")
+        haskey(opts, parts[1]) || error("Unknown focus CLI option: --$(parts[1])")
         opts[parts[1]] = parts[2]
+        push!(provided, parts[1])
+    end
+    return opts, provided
+end
+
+parse_estimator(s) = lowercase(strip(s)) == "geometric" ? GEOMETRIC :
+    lowercase(strip(s)) == "hasa" ? HASA :
+    error("--estimator must be geometric or hasa, got: $s")
+
+function parse_focus_aberrator(s)
+    value = lowercase(strip(s))
+    value in ("none", "water") && return WATER, :water
+    value == "skull" && return SKULL_IN_WATER, :skull
+    error("--aberrator must be none, water, or skull, got: $s")
+end
+
+function parse_record_mode(s)
+    value = Symbol(strip(s))
+    value in (:p_rms, :p) || error("--record must be p_rms or p, got: $s")
+    return value
+end
+
+function parse_outward_propagation(s)
+    value = replace(lowercase(strip(s)), "-" => "_")
+    value in ("angular_spectrum", "kwave") || error("--outward-propagation must be angular_spectrum or kwave, got: $s")
+    return Symbol(value)
+end
+
+function parse_focus_mm(spec::AbstractString, dimension::Int)
+    text = strip(spec)
+    if isempty(text)
+        text = dimension == 3 ? "60:0:0" : "60:0"
+    end
+    values = [parse(Float64, strip(item)) * 1e-3 for item in split(text, ":")]
+    expected = dimension == 3 ? 3 : 2
+    length(values) == expected || error("--focus-mm must have $expected colon-separated values for $(dimension)D.")
+    return values
+end
+
+function resolve_focus_placement(opts, dimension::Int, medium_type::MediumType)
+    requested_focus_depth = parse_optional_mm(opts["focus-depth-from-inner-skull-mm"])
+    placement = parse_placement_mode(opts["placement"])
+    if dimension == 3
+        isnothing(requested_focus_depth) || error("3D focus currently supports fixed transducer placement only; omit --focus-depth-from-inner-skull-mm.")
+        placement == :fixed_focus_depth && error("3D focus currently supports fixed transducer placement only; use --placement=fixed_transducer or --placement=auto.")
+        return :fixed_transducer, nothing
     end
-    return opts
+    return resolve_placement_mode(
+        placement,
+        medium_type;
+        focus_depth_from_inner_skull=requested_focus_depth,
+    )
+end
+
+function parse_optional_mm(spec::AbstractString)
+    isempty(strip(spec)) && return nothing
+    return parse(Float64, spec) * 1e-3
 end
 
-parse_est(s) = lowercase(s) == "geometric" ? GEOMETRIC : lowercase(s) == "hasa" ? HASA : error("Unknown estimator: $s")
-parse_medium(s) = lowercase(s) == "water" ? WATER : lowercase(s) == "skull_in_water" ? SKULL_IN_WATER : error("Unknown medium: $s")
+function parse_axial_dim(spec::AbstractString)
+    lowercase(strip(spec)) == "auto" && return nothing
+    return parse(Float64, spec) * 1e-3
+end
 
 slug_value(x; digits::Int=1) = replace(string(round(Float64(x); digits=digits)), "-" => "m", "." => "p")
 
-function default_output_dir(opts, placement_mode, focus_depth)
+function default_focus_output_dir(opts, dimension::Int, placement_mode, focus_depth)
     timestamp = Dates.format(Dates.now(), "yyyymmdd_HHMMSS")
     parts = String[
         timestamp,
         "run_focus",
+        "$(dimension)d",
         lowercase(opts["estimator"]),
-        lowercase(opts["medium"]),
+        lowercase(opts["aberrator"]),
         String(placement_mode),
         "f$(slug_value(parse(Float64, opts["frequency-mhz"]); digits=2))mhz",
-        "z$(slug_value(parse(Float64, opts["focal-cm"]); digits=1))cm",
-        "x$(slug_value(parse(Float64, opts["lateral-cm"]); digits=1))cm",
-        "ap$(slug_value(parse(Float64, opts["aperture-cm"]); digits=1))cm",
-        "slice$(parse(Int, opts["slice-index"]))",
+        "focus$(replace(strip(isempty(strip(opts["focus-mm"])) ? (dimension == 3 ? "60:0:0" : "60:0") : opts["focus-mm"]), ":" => "_"))mm",
     ]
     if !isnothing(focus_depth)
         push!(parts, "fd$(slug_value(focus_depth * 1e3; digits=1))mm")
@@ -56,17 +138,16 @@ function default_output_dir(opts, placement_mode, focus_depth)
     return joinpath(pwd(), "outputs", join(parts, "_"))
 end
 
-function target_mm(kgrid, cfg)
+function target_mm(kgrid::KGrid2D, cfg::SimulationConfig)
     x_tgt_mm = (kgrid.x_vec[cfg.trans_index] - cfg.z_focus) * 1e3
     y_tgt_mm = (kgrid.y_vec[fld(kgrid.Ny, 2) + 1] + cfg.x_focus) * 1e3
     return x_tgt_mm, y_tgt_mm
 end
 
-function overlay_skull!(ax, c_map, kgrid)
+function overlay_skull_2d!(ax, c_map, kgrid)
     c_map === nothing && return nothing
     skull_mask = skull_mask_from_c_columnwise(c_map; mask_outside=false)
     any(skull_mask) || return nothing
-
     overlay = fill(Float64(NaN), size(c_map))
     overlay[skull_mask] .= Float64.(c_map[skull_mask]) ./ maximum(Float64.(c_map[skull_mask]))
     return heatmap!(
@@ -81,7 +162,7 @@ function overlay_skull!(ax, c_map, kgrid)
     )
 end
 
-function save_pressure_plot(path::AbstractString, pressure, c_map, kgrid, cfg, stats, title::AbstractString)
+function save_pressure_plot_2d(path::AbstractString, pressure, c_map, kgrid, cfg, stats, title::AbstractString)
     fig = Figure(size=(900, 600))
     ax = Axis(
         fig[1, 1];
@@ -91,7 +172,7 @@ function save_pressure_plot(path::AbstractString, pressure, c_map, kgrid, cfg, s
         aspect=DataAspect(),
     )
     hm = heatmap!(ax, kgrid.y_vec .* 1e3, kgrid.x_vec .* 1e3, pressure'; colormap=:viridis)
-    overlay_skull!(ax, c_map, kgrid)
+    overlay_skull_2d!(ax, c_map, kgrid)
     Colorbar(fig[1, 2], hm; label="Pressure")
 
     x_tgt_mm, y_tgt_mm = target_mm(kgrid, cfg)
@@ -103,79 +184,285 @@ function save_pressure_plot(path::AbstractString, pressure, c_map, kgrid, cfg, s
     save(path, fig)
 end
 
-opts = parse_cli(ARGS)
-medium_type = parse_medium(opts["medium"])
-requested_focus_depth = haskey(opts, "focus-depth-from-inner-skull-mm") ?
-    parse(Float64, opts["focus-depth-from-inner-skull-mm"]) * 1e-3 :
-    nothing
-placement_mode, focus_depth = resolve_placement_mode(
-    opts["placement"],
-    medium_type;
-    focus_depth_from_inner_skull=requested_focus_depth,
-)
-out_dir = if haskey(opts, "out-dir") && !isempty(strip(opts["out-dir"]))
-    opts["out-dir"]
-else
-    default_output_dir(opts, placement_mode, focus_depth)
-end
-mkpath(out_dir)
-
-hu_vol, spacing_m = load_default_ct(ct_path=opts["ct-path"])
-ct_info = CTInfo(hu_vol, spacing_m)
-@info "Loaded CT volume" size=size(hu_vol) spacing_m spacing_info=ct_info
-
-cfg = SimulationConfig(
-    fc=parse(Float64, opts["frequency-mhz"]) * 1e6,
-    z_focus=parse(Float64, opts["focal-cm"]) * 1e-2,
-    x_focus=parse(Float64, opts["lateral-cm"]) * 1e-2,
-    trans_aperture=parse(Float64, opts["aperture-cm"]) * 1e-2,
-    focus_depth_from_inner_skull=focus_depth,
-)
+function save_pressure_plot_3d(path::AbstractString, pressure, c_map, kgrid, cfg, stats, title::AbstractString)
+    max_y = dropdims(maximum(pressure; dims=3); dims=3)
+    max_z = dropdims(maximum(pressure; dims=2); dims=2)
+    fig = Figure(size=(1200, 520))
+    ax1 = Axis(fig[1, 1]; title="$title depth-y", xlabel="Y [mm]", ylabel="Axial [mm]", aspect=DataAspect())
+    hm1 = heatmap!(ax1, kgrid.y .* 1e3, kgrid.x .* 1e3, max_y'; colormap=:viridis)
+    Colorbar(fig[1, 2], hm1; label="Pressure")
+    ax2 = Axis(fig[1, 3]; title="depth-z", xlabel="Z [mm]", ylabel="Axial [mm]", aspect=DataAspect())
+    hm2 = heatmap!(ax2, kgrid.z .* 1e3, kgrid.x .* 1e3, max_z'; colormap=:viridis)
+    Colorbar(fig[1, 4], hm2; label="Pressure")
 
-stats, pressure, c_map, kgrid, cfg, hasa_info = run_focus_case(
-    hu_vol,
-    cfg,
-    medium_type,
-    parse_est(opts["estimator"]),
-    SweepSettings();
-    slice_index=parse(Int, opts["slice-index"]),
-    return_c=true,
-)
+    focus_x = kgrid.x[target_index(cfg)] * 1e3
+    scatter!(ax1, [cfg.lateral_y_focus * 1e3], [focus_x]; color=:red, marker=:x, markersize=16, strokewidth=2)
+    scatter!(ax2, [cfg.lateral_z_focus * 1e3], [focus_x]; color=:red, marker=:x, markersize=16, strokewidth=2)
+    if stats !== nothing
+        peak_x, peak_y, peak_z = stats[:peak_mm]
+        scatter!(ax1, [peak_y], [peak_x]; color=:white, marker=:circle, markersize=12, strokecolor=:white, strokewidth=2)
+        scatter!(ax2, [peak_z], [peak_x]; color=:white, marker=:circle, markersize=12, strokecolor=:white, strokewidth=2)
+    end
+    save(path, fig)
+end
 
-summary = Dict(
-    "out_dir" => out_dir,
-    "ct_path" => opts["ct-path"],
-    "slice_index" => parse(Int, opts["slice-index"]),
-    "spacing_m" => spacing_m,
-    "placement_mode" => String(placement_mode),
-    "stats" => stats,
-    "config" => Dict(
-        "fc" => cfg.fc,
-        "x_focus" => cfg.x_focus,
-        "z_focus" => cfg.z_focus,
-        "trans_aperture" => cfg.trans_aperture,
-        "focus_depth_from_inner_skull" => cfg.focus_depth_from_inner_skull,
-    ),
-)
+function save_sensor_timing_plot_2d(path::AbstractString, kgrid::KGrid2D, hasa_info, title::AbstractString)
+    tau_us = Float64.(hasa_info[:tau]) .* 1e6
+    amp = Float64.(hasa_info[:amplitudes])
+    cols = hasa_info[:transducer_cols]
+    y_mm = kgrid.y_vec[cols] .* 1e3
+
+    fig = Figure(size=(900, 520))
+    ax1 = Axis(fig[1, 1]; title="$title delay", xlabel="Sensor position [mm]", ylabel="Delay [us]")
+    lines!(ax1, y_mm, tau_us; color=:black, linewidth=2.0)
+    scatter!(ax1, y_mm, tau_us; color=:black, markersize=4)
+    ax2 = Axis(fig[2, 1]; title="Normalized amplitude", xlabel="Sensor position [mm]", ylabel="Amplitude")
+    lines!(ax2, y_mm, amp; color=:steelblue, linewidth=2.0)
+    scatter!(ax2, y_mm, amp; color=:steelblue, markersize=4)
+    save(path, fig)
+    return path
+end
 
-open(joinpath(out_dir, "summary.json"), "w") do io
-    JSON3.pretty(io, summary)
+function save_sensor_timing_plot_3d(path::AbstractString, kgrid, hasa_info, title::AbstractString)
+    tau_us = Float64.(hasa_info[:tau]) .* 1e6
+    amp = Float64.(hasa_info[:amplitudes])
+    range_y = hasa_info[:transducer_cols_y]
+    range_z = hasa_info[:transducer_cols_z]
+    y_mm = collect(kgrid.y)[range_y] .* 1e3
+    z_mm = collect(kgrid.z)[range_z] .* 1e3
+
+    fig = Figure(size=(1200, 520))
+    ax1 = Axis(fig[1, 1]; title="$title delay", xlabel="Z [mm]", ylabel="Y [mm]", aspect=DataAspect())
+    hm1 = heatmap!(ax1, z_mm, y_mm, tau_us'; colormap=:viridis)
+    Colorbar(fig[1, 2], hm1; label="Delay [us]")
+    ax2 = Axis(fig[1, 3]; title="Normalized amplitude", xlabel="Z [mm]", ylabel="Y [mm]", aspect=DataAspect())
+    hm2 = heatmap!(ax2, z_mm, y_mm, amp'; colormap=:magma, colorrange=(0, max(maximum(amp), eps(Float64))))
+    Colorbar(fig[1, 4], hm2; label="Amplitude")
+    save(path, fig)
+    return path
 end
 
-hasa_save = copy(hasa_info)
-hasa_save[:fig] = nothing
-@save joinpath(out_dir, "result.jld2") pressure c_map kgrid cfg stats hasa_save
+function load_focus_ct_if_needed(opts, aberrator::Symbol)
+    aberrator == :skull || return nothing, nothing, nothing
+    hu_vol, spacing_m = load_default_ct(ct_path=opts["ct-path"])
+    return hu_vol, spacing_m, CTInfo(hu_vol, spacing_m)
+end
 
-if pressure !== nothing && ndims(pressure) == 2
-    save_pressure_plot(
-        joinpath(out_dir, "pressure.png"),
-        pressure,
-        c_map,
-        kgrid,
+function focus_dry_plan(args::AbstractVector{<:AbstractString})
+    opts, _ = parse_focus_cli(String.(args))
+    dimension = parse_dimension(opts["dimension"])
+    medium_type, aberrator = parse_focus_aberrator(opts["aberrator"])
+    placement_mode, focus_depth = resolve_focus_placement(opts, dimension, medium_type)
+    focus = parse_focus_mm(opts["focus-mm"], dimension)
+    out_dir = isempty(strip(opts["out-dir"])) ?
+        default_focus_output_dir(opts, dimension, placement_mode, focus_depth) :
+        opts["out-dir"]
+    return Dict(
+        :branch => dimension == 3 ? :focus3d : :focus2d,
+        :out_dir => out_dir,
+        :aberrator => aberrator,
+        :estimator => parse_estimator(opts["estimator"]),
+        :placement_mode => placement_mode,
+        :focus_m => focus,
+        :outward_propagation => parse_outward_propagation(opts["outward-propagation"]),
+        :verify_inward_kwave => parse_bool(opts["verify-inward-kwave"]) && parse_bool(opts["run-kwave"]),
+    )
+end
+
+function main(args::AbstractVector{<:AbstractString}=ARGS; dry_run::Bool=false)
+    dry_run && return focus_dry_plan(args)
+
+    opts, _ = parse_focus_cli(String.(args))
+    dimension = parse_dimension(opts["dimension"])
+    medium_type, aberrator = parse_focus_aberrator(opts["aberrator"])
+    estimator = parse_estimator(opts["estimator"])
+    record = parse_record_mode(opts["record"])
+    placement_mode, focus_depth = resolve_focus_placement(opts, dimension, medium_type)
+
+    out_dir = isempty(strip(opts["out-dir"])) ?
+        default_focus_output_dir(opts, dimension, placement_mode, focus_depth) :
+        opts["out-dir"]
+    mkpath(out_dir)
+    hu_vol, spacing_m, ct_info = load_focus_ct_if_needed(opts, aberrator)
+    if ct_info !== nothing
+        @info "Loaded CT volume" size=size(hu_vol) spacing_m spacing_info=ct_info
+    end
+
+    verify_inward_kwave = parse_bool(opts["verify-inward-kwave"]) && parse_bool(opts["run-kwave"])
+    animation_settings = verify_inward_kwave ? nothing : AnimationSettings(run_kwave=false, Nt=round(Int, parse(Float64, opts["t-max-us"]) * 1e-6 / (parse(Float64, opts["dt-ns"]) * 1e-9)))
+    sweep = SweepSettings(
+        use_gpu=parse_bool(opts["kwave-use-gpu"]),
+        record=record,
+        outward_solver=parse_outward_propagation(opts["outward-propagation"]),
+        verify_inward_kwave=verify_inward_kwave,
+    )
+
+    if dimension == 3
+        focus = parse_focus_mm(opts["focus-mm"], 3)
+        dy_mm = isempty(strip(opts["dy-mm"])) ? parse(Float64, opts["dz-mm"]) : parse(Float64, opts["dy-mm"])
+        transverse_y_mm = isempty(strip(opts["transverse-y-mm"])) ? parse(Float64, opts["transverse-mm"]) : parse(Float64, opts["transverse-y-mm"])
+        transverse_z_mm = isempty(strip(opts["transverse-z-mm"])) ? parse(Float64, opts["transverse-mm"]) : parse(Float64, opts["transverse-z-mm"])
+        aperture_y_spec = isempty(strip(opts["receiver-aperture-y-mm"])) ? opts["receiver-aperture-mm"] : opts["receiver-aperture-y-mm"]
+        aperture_z_spec = isempty(strip(opts["receiver-aperture-z-mm"])) ? opts["receiver-aperture-mm"] : opts["receiver-aperture-z-mm"]
+
+        cfg = SimulationConfig3D(
+            fc=parse(Float64, opts["frequency-mhz"]) * 1e6,
+            axial_focus=focus[1],
+            lateral_y_focus=focus[2],
+            lateral_z_focus=focus[3],
+            trans_aperture_y=parse_receiver_aperture_mm(aperture_y_spec),
+            trans_aperture_z=parse_receiver_aperture_mm(aperture_z_spec),
+            trans_skull_dist=parse(Float64, opts["skull-transducer-distance-mm"]) * 1e-3,
+            dx=parse(Float64, opts["dx-mm"]) * 1e-3,
+            dy=dy_mm * 1e-3,
+            dz=parse(Float64, opts["dz-mm"]) * 1e-3,
+            axial_dim=parse_axial_dim(opts["axial-mm"]),
+            transverse_dim_y=transverse_y_mm * 1e-3,
+            transverse_dim_z=transverse_z_mm * 1e-3,
+            axial_padding=parse(Float64, opts["axial-padding"]),
+            t_max=parse(Float64, opts["t-max-us"]) * 1e-6,
+            dt=parse(Float64, opts["dt-ns"]) * 1e-9,
+        )
+
+        stats, pressure, c_map, kgrid, cfg, hasa_info, medium_info = run_focus_case_3d(
+            hu_vol,
+            spacing_m,
+            cfg,
+            medium_type,
+            estimator,
+            sweep;
+            ct_path=opts["ct-path"],
+            slice_index_z=parse(Int, opts["slice-index"]),
+            hu_bone_thr=parse(Int, opts["hu-bone-thr"]),
+            animation_settings=animation_settings,
+            return_c=true,
+        )
+        timing_path = save_sensor_timing_plot_3d(
+            joinpath(out_dir, "sensor_timings.png"),
+            kgrid,
+            hasa_info,
+            "3D sensor timing",
+        )
+
+        summary = Dict(
+            "out_dir" => out_dir,
+            "dimension" => 3,
+            "ct_path" => opts["ct-path"],
+            "slice_index" => parse(Int, opts["slice-index"]),
+            "spacing_m" => spacing_m,
+            "placement_mode" => String(placement_mode),
+            "outward_propagation" => opts["outward-propagation"],
+            "verify_inward_kwave" => verify_inward_kwave,
+            "sensor_timing_figure" => timing_path,
+            "stats" => stats,
+            "medium" => TranscranialFUS.run_pam_medium_summary(medium_info),
+            "config" => Dict(
+                "fc" => cfg.fc,
+                "axial_focus" => cfg.axial_focus,
+                "lateral_y_focus" => cfg.lateral_y_focus,
+                "lateral_z_focus" => cfg.lateral_z_focus,
+                "trans_aperture_y" => cfg.trans_aperture_y,
+                "trans_aperture_z" => cfg.trans_aperture_z,
+                "trans_skull_dist" => cfg.trans_skull_dist,
+                "dx" => cfg.dx,
+                "dy" => cfg.dy,
+                "dz" => cfg.dz,
+                "axial_dim" => cfg.axial_dim,
+                "transverse_dim_y" => cfg.transverse_dim_y,
+                "transverse_dim_z" => cfg.transverse_dim_z,
+                "t_max" => cfg.t_max,
+                "dt" => cfg.dt,
+            ),
+        )
+        open(joinpath(out_dir, "summary.json"), "w") do io
+            JSON3.pretty(io, summary)
+        end
+        @save joinpath(out_dir, "result.jld2") pressure c_map kgrid cfg stats hasa_info medium_info
+        if pressure !== nothing && ndims(pressure) == 3
+            save_pressure_plot_3d(joinpath(out_dir, "pressure.png"), pressure, c_map, kgrid, cfg, stats, "3D Focus Run")
+        end
+        println("Saved 3D focus outputs to $out_dir")
+        return 0
+    end
+
+    focus = parse_focus_mm(opts["focus-mm"], 2)
+    cfg = SimulationConfig(
+        fc=parse(Float64, opts["frequency-mhz"]) * 1e6,
+        z_focus=focus[1],
+        x_focus=focus[2],
+        trans_aperture=parse_receiver_aperture_mm(opts["receiver-aperture-mm"]),
+        trans_skull_dist=parse(Float64, opts["skull-transducer-distance-mm"]) * 1e-3,
+        focus_depth_from_inner_skull=focus_depth,
+        dx=parse(Float64, opts["dx-mm"]) * 1e-3,
+        dz=parse(Float64, opts["dz-mm"]) * 1e-3,
+        transverse_dim=parse(Float64, opts["transverse-mm"]) * 1e-3,
+        axial_padding=parse(Float64, opts["axial-padding"]),
+        t_max=parse(Float64, opts["t-max-us"]) * 1e-6,
+        dt=parse(Float64, opts["dt-ns"]) * 1e-9,
+    )
+    if lowercase(strip(opts["axial-mm"])) == "auto"
+        cfg.axial_dim = max(cfg.axial_dim, cfg.z_focus + (cfg.PML_GUARD + 1) * cfg.dx)
+    else
+        cfg.axial_dim = parse(Float64, opts["axial-mm"]) * 1e-3
+    end
+
+    stats, pressure, c_map, kgrid, cfg, hasa_info = run_focus_case(
+        something(hu_vol, zeros(Float32, 1, 1, 1)),
         cfg,
-        stats,
-        "Focus Run: $(opts["estimator"]) in $(opts["medium"])",
+        medium_type,
+        estimator,
+        sweep;
+        slice_index=parse(Int, opts["slice-index"]),
+        hu_bone_thr=parse(Int, opts["hu-bone-thr"]),
+        animation_settings=animation_settings,
+        return_c=true,
     )
+    timing_path = save_sensor_timing_plot_2d(
+        joinpath(out_dir, "sensor_timings.png"),
+        kgrid,
+        hasa_info,
+        "2D sensor timing",
+    )
+
+    summary = Dict(
+        "out_dir" => out_dir,
+        "dimension" => 2,
+        "ct_path" => opts["ct-path"],
+        "slice_index" => parse(Int, opts["slice-index"]),
+        "spacing_m" => spacing_m,
+        "placement_mode" => String(placement_mode),
+        "outward_propagation" => opts["outward-propagation"],
+        "verify_inward_kwave" => verify_inward_kwave,
+        "sensor_timing_figure" => timing_path,
+        "stats" => stats,
+        "config" => Dict(
+            "fc" => cfg.fc,
+            "x_focus" => cfg.x_focus,
+            "z_focus" => cfg.z_focus,
+            "trans_aperture" => cfg.trans_aperture,
+            "trans_skull_dist" => cfg.trans_skull_dist,
+            "focus_depth_from_inner_skull" => cfg.focus_depth_from_inner_skull,
+            "dx" => cfg.dx,
+            "dz" => cfg.dz,
+            "axial_dim" => cfg.axial_dim,
+            "transverse_dim" => cfg.transverse_dim,
+            "t_max" => cfg.t_max,
+            "dt" => cfg.dt,
+        ),
+    )
+    open(joinpath(out_dir, "summary.json"), "w") do io
+        JSON3.pretty(io, summary)
+    end
+    hasa_save = copy(hasa_info)
+    hasa_save[:fig] = nothing
+    @save joinpath(out_dir, "result.jld2") pressure c_map kgrid cfg stats hasa_save
+    if pressure !== nothing && ndims(pressure) == 2
+        save_pressure_plot_2d(joinpath(out_dir, "pressure.png"), pressure, c_map, kgrid, cfg, stats, "2D Focus Run")
+    end
+    println("Saved 2D focus outputs to $out_dir")
+    return 0
 end
 
-println("Saved outputs to $out_dir")
+if abspath(PROGRAM_FILE) == @__FILE__
+    exit(main())
+end
diff --git a/src/TranscranialFUS.jl b/src/TranscranialFUS.jl
index 8fdf547..15c0b6b 100644
--- a/src/TranscranialFUS.jl
+++ b/src/TranscranialFUS.jl
@@ -22,14 +22,15 @@ const DEFAULT_ROI_SIZE_XYZ = (705, 360, 450)
 
 export DEFAULT_CT_PATH, DEFAULT_ROI_INDEX_XYZ, DEFAULT_ROI_SIZE_XYZ
 export source_summary
-export CTInfo, KGrid2D, SimulationConfig, SweepSettings, AnimationSettings, MediumType, Est
+export CTInfo, KGrid2D, SimulationConfig, SimulationConfig3D, SweepSettings, AnimationSettings, MediumType, Est
 export WATER, SKULL_IN_WATER, GEOMETRIC, HASA
 export parse_placement_mode, resolve_placement_mode
 export omega, Nx, Nz, Nt, Nx_hasa, target_index, Nz_active, active_col_range, set_z_focus!
 export load_roi_resample_xy, load_default_ct
 export hu_to_rho_c, find_skull_boundaries, skull_mask_from_c_columnwise
 export make_medium_fixed_distance_from_skull, make_medium_fixed_transducer, make_medium
-export plot_hasa_results, focus, analyse_focus_2d, run_focus_case, kwave_available
+export plot_hasa_results, focus, analyse_focus_2d, analyse_focus_3d, run_focus_case, run_focus_case_3d, kwave_available
+export focus_pam_config_3d, focus_grid_3d, active_col_range_y, active_col_range_z, make_focus_medium_3d, skull_mask_from_c_columnwise_3d
 export EmissionSource2D, PointSource2D, BubbleCluster2D, PAMConfig, PAMWindowConfig, SourceVariabilityConfig, fit_pam_config, pam_Nx, pam_Ny, pam_Nt, pam_grid, receiver_row, receiver_col_range, depth_coordinates
 export emission_frequencies
 export make_squiggle_bubble_sources, make_pam_medium, source_grid_index, simulate_point_sources, simulate_point_sources_3d, reconstruct_pam, reconstruct_pam_windowed, find_pam_peaks
@@ -42,5 +43,6 @@ export find_pam_peaks_3d, pam_truth_mask_3d, source_detection_stats_3d, threshol
 
 include("focus.jl")
 include("pam.jl")
+include("focus/3d.jl")
 
 end
diff --git a/src/focus/3d.jl b/src/focus/3d.jl
new file mode 100644
index 0000000..4c193fa
--- /dev/null
+++ b/src/focus/3d.jl
@@ -0,0 +1,699 @@
+mutable struct SimulationConfig3D
+    fc::Float64
+    axial_focus::Float64
+    lateral_y_focus::Float64
+    lateral_z_focus::Float64
+    trans_aperture_y::Union{Nothing, Float64}
+    trans_aperture_z::Union{Nothing, Float64}
+    trans_index::Int
+    trans_skull_dist::Float64
+    focus_depth_from_inner_skull::Union{Nothing, Float64}
+    dx::Float64
+    dy::Float64
+    dz::Float64
+    axial_dim::Float64
+    transverse_dim_y::Float64
+    transverse_dim_z::Float64
+    axial_padding::Float64
+    t_max::Float64
+    dt::Float64
+    c0::Float64
+    rho0::Float64
+    PML_GUARD::Int
+end
+
+function SimulationConfig3D(;
+    fc::Real=5e5,
+    axial_focus::Real=50e-3,
+    lateral_y_focus::Real=0.0,
+    lateral_z_focus::Real=0.0,
+    trans_aperture_y::Union{Nothing, Real}=60e-3,
+    trans_aperture_z::Union{Nothing, Real}=60e-3,
+    trans_index::Union{Nothing, Integer}=nothing,
+    trans_skull_dist::Real=5e-3,
+    focus_depth_from_inner_skull::Union{Nothing, Real}=nothing,
+    dx::Real=0.2e-3,
+    dy::Real=0.2e-3,
+    dz::Real=0.2e-3,
+    axial_dim::Union{Nothing, Real}=nothing,
+    transverse_dim_y::Real=80e-3,
+    transverse_dim_z::Real=80e-3,
+    axial_padding::Real=1.5,
+    t_max::Real=220e-6,
+    dt::Real=40e-9,
+    c0::Real=1500.0,
+    rho0::Real=1000.0,
+    PML_GUARD::Integer=20,
+)
+    dx_f = Float64(dx)
+    dy_f = Float64(dy)
+    dz_f = Float64(dz)
+    guard = Int(PML_GUARD)
+    tx = isnothing(trans_index) ? guard + 1 : Int(trans_index)
+    tx >= 1 || error("trans_index must be >= 1.")
+
+    axial_focus_f = Float64(axial_focus)
+    axial_dim_f = isnothing(axial_dim) ?
+        (tx + round(Int, Float64(axial_padding) * axial_focus_f / dx_f)) * dx_f :
+        Float64(axial_dim)
+    transverse_y_f = Float64(transverse_dim_y)
+    transverse_z_f = Float64(transverse_dim_z)
+    ap_y = isnothing(trans_aperture_y) ? nothing : Float64(trans_aperture_y)
+    ap_z = isnothing(trans_aperture_z) ? nothing : Float64(trans_aperture_z)
+    !isnothing(ap_y) && ap_y > transverse_y_f && error("Y transducer aperture must be <= transverse Y grid size.")
+    !isnothing(ap_z) && ap_z > transverse_z_f && error("Z transducer aperture must be <= transverse Z grid size.")
+    focus_depth = isnothing(focus_depth_from_inner_skull) ? nothing : Float64(focus_depth_from_inner_skull)
+
+    return SimulationConfig3D(
+        Float64(fc),
+        axial_focus_f,
+        Float64(lateral_y_focus),
+        Float64(lateral_z_focus),
+        ap_y,
+        ap_z,
+        tx,
+        Float64(trans_skull_dist),
+        focus_depth,
+        dx_f,
+        dy_f,
+        dz_f,
+        axial_dim_f,
+        transverse_y_f,
+        transverse_z_f,
+        Float64(axial_padding),
+        Float64(t_max),
+        Float64(dt),
+        Float64(c0),
+        Float64(rho0),
+        guard,
+    )
+end
+
+omega(cfg::SimulationConfig3D) = 2π * cfg.fc
+Nx(cfg::SimulationConfig3D) = round(Int, cfg.axial_dim / cfg.dx)
+Ny(cfg::SimulationConfig3D) = round(Int, cfg.transverse_dim_y / cfg.dy)
+Nz(cfg::SimulationConfig3D) = round(Int, cfg.transverse_dim_z / cfg.dz)
+Nt(cfg::SimulationConfig3D) = round(Int, cfg.t_max / cfg.dt)
+Nx_hasa(cfg::SimulationConfig3D) = round(Int, cfg.axial_focus / cfg.dx)
+target_index(cfg::SimulationConfig3D) = cfg.trans_index + Nx_hasa(cfg)
+
+function focus_pam_config_3d(cfg::SimulationConfig3D)
+    return PAMConfig3D(
+        dx=cfg.dx,
+        dy=cfg.dy,
+        dz=cfg.dz,
+        axial_dim=cfg.axial_dim,
+        transverse_dim_y=cfg.transverse_dim_y,
+        transverse_dim_z=cfg.transverse_dim_z,
+        dt=cfg.dt,
+        t_max=cfg.t_max,
+        c0=cfg.c0,
+        rho0=cfg.rho0,
+        receiver_row=cfg.trans_index,
+        receiver_aperture_y=cfg.trans_aperture_y,
+        receiver_aperture_z=cfg.trans_aperture_z,
+        PML_GUARD=cfg.PML_GUARD,
+    )
+end
+
+function focus_grid_3d(cfg::SimulationConfig3D; nt::Integer=Nt(cfg))
+    return (
+        x = range((1 - cfg.trans_index) * cfg.dx; step=cfg.dx, length=Nx(cfg)),
+        y = _kwave_axis(Ny(cfg), cfg.dy),
+        z = _kwave_axis(Nz(cfg), cfg.dz),
+        t = range(0.0; step=cfg.dt, length=Int(nt)),
+        Nt = Int(nt),
+    )
+end
+
+function active_col_range_y(cfg::SimulationConfig3D)
+    ny = Ny(cfg)
+    isnothing(cfg.trans_aperture_y) && return 1:ny
+    n_active = clamp(round(Int, cfg.trans_aperture_y / cfg.dy), 1, ny)
+    mid = fld(ny, 2) + 1
+    start_col = mid - fld(n_active, 2)
+    return start_col:(start_col + n_active - 1)
+end
+
+function active_col_range_z(cfg::SimulationConfig3D)
+    nz = Nz(cfg)
+    isnothing(cfg.trans_aperture_z) && return 1:nz
+    n_active = clamp(round(Int, cfg.trans_aperture_z / cfg.dz), 1, nz)
+    mid = fld(nz, 2) + 1
+    start_col = mid - fld(n_active, 2)
+    return start_col:(start_col + n_active - 1)
+end
+
+_fftshift2(a::AbstractMatrix) = circshift(a, (fld(size(a, 1), 2), fld(size(a, 2), 2)))
+_ifftshift2(a::AbstractMatrix) = circshift(a, (-fld(size(a, 1), 2), -fld(size(a, 2), 2)))
+
+function _unwrap_phase_matrix(a::AbstractMatrix{<:Real}; discont::Real=π)
+    out = _unwrap_phase_rows(a; discont=discont)
+    for col in axes(out, 2)
+        out[:, col] = _unwrap_phase(view(out, :, col); discont=discont)
+    end
+    return out
+end
+
+function _focus_wavenumbers(n::Int, d::Float64)
+    start_val = -fld(n, 2)
+    end_val = ceil(Int, n / 2) - 1
+    return collect(start_val:end_val) .* (2π / (n * d))
+end
+
+function _focus_seed_3d(y::AbstractVector{<:Real}, z::AbstractVector{<:Real}, cfg::SimulationConfig3D)
+    sigma = max(2 * max(cfg.dy, cfg.dz), 1.0e-3)
+    p0 = Matrix{Float64}(undef, length(y), length(z))
+    @inbounds for iz in eachindex(z), iy in eachindex(y)
+        r2 = (Float64(y[iy]) - cfg.lateral_y_focus)^2 + (Float64(z[iz]) - cfg.lateral_z_focus)^2
+        p0[iy, iz] = 1e3 * exp(-r2 / (2 * sigma^2))
+    end
+    return p0
+end
+
+function _simulate_kwave_3d(
+    c::AbstractArray{<:Real, 3},
+    rho::AbstractArray{<:Real, 3},
+    tau::AbstractMatrix{<:Real},
+    amplitudes::AbstractMatrix{<:Real},
+    cfg::SimulationConfig3D,
+    kgrid,
+    col_range_y::UnitRange{Int},
+    col_range_z::UnitRange{Int};
+    record::Union{Symbol, AbstractString}=:p_rms,
+    use_gpu::Bool=false,
+    animation_settings::Union{Nothing, AnimationSettings}=nothing,
+)
+    mods = _kwave_modules()
+    pc = mods.PythonCall
+    np = mods.np
+    deepcopy_py = mods.copy.deepcopy
+
+    rec_symbol = _normalize_record(record)
+    nx, ny, nz = size(c)
+    nt = isnothing(animation_settings) ? Nt(cfg) : animation_settings.Nt
+    num_cycles = isnothing(animation_settings) ? 40 : animation_settings.num_cycles
+    use_envelope = isnothing(animation_settings)
+    rise_fall = max(1, floor(Int, 0.1 * num_cycles))
+    pml_guard = max(4, cfg.PML_GUARD)
+
+    py_kgrid = mods.kgrid.kWaveGrid(pc.Py([nx, ny, nz]), pc.Py([cfg.dx, cfg.dy, cfg.dz]))
+    py_kgrid.dt = cfg.dt
+    py_kgrid.Nt = nt
+
+    medium = mods.kmedium.kWaveMedium(
+        sound_speed=np.array(pc.Py(Float32.(c)), dtype=np.float32),
+        density=np.array(pc.Py(Float32.(rho)), dtype=np.float32),
+    )
+
+    source = mods.ksource.kSource()
+    src_mask_jl = falses(nx, ny, nz)
+    src_mask_jl[cfg.trans_index, col_range_y, col_range_z] .= true
+    source.p_mask = np.array(pc.Py(Array(src_mask_jl)), dtype=pc.pybuiltins.bool)
+
+    tau_flat = vec(Float64.(tau))
+    amp_flat = vec(Float64.(amplitudes))
+    int_delay = floor.(Int, tau_flat ./ cfg.dt)
+    frac_delay = tau_flat .- int_delay .* cfg.dt
+    int_delay_py = np.array(pc.Py(Int64.(int_delay)), dtype=np.int64)
+    envelope_py = np.array(pc.Py(Int64[rise_fall, rise_fall]), dtype=np.int64)
+
+    if use_envelope
+        burst_py = mods.signals.tone_burst(
+            sample_freq=1 / cfg.dt,
+            signal_freq=cfg.fc,
+            num_cycles=num_cycles,
+            envelope=envelope_py,
+            signal_offset=int_delay_py,
+        )
+    else
+        burst_py = mods.signals.tone_burst(
+            sample_freq=1 / cfg.dt,
+            signal_freq=cfg.fc,
+            num_cycles=num_cycles,
+            signal_offset=int_delay_py,
+        )
+    end
+    burst = pc.pyconvert(Matrix{Float64}, burst_py)
+    drive_weights = 1e5 .* amp_flat .* exp.(-1im .* omega(cfg) .* frac_delay)
+    source.p = np.array(pc.Py(real.(reshape(drive_weights, :, 1) .* burst)), dtype=np.float64)
+    source.p_frequency_ref = cfg.fc
+    source.medium = medium
+
+    sensor = mods.ksensor.kSensor()
+    sensor.mask = np.ones((nx, ny, nz), dtype=pc.pybuiltins.bool)
+    sensor.record = pc.pybuiltins.list((String(rec_symbol),))
+
+    sim_opts = mods.simopts.SimulationOptions(
+        pml_inside=true,
+        pml_size=pml_guard,
+        data_recast=false,
+        save_to_disk=true,
+    )
+    exec_opts = mods.execopts.SimulationExecutionOptions(
+        is_gpu_simulation=use_gpu,
+        delete_data=true,
+    )
+
+    data = mods.kspace3d.kspaceFirstOrder3D(
+        kgrid=deepcopy_py(py_kgrid),
+        medium=deepcopy_py(medium),
+        source=deepcopy_py(source),
+        sensor=deepcopy_py(sensor),
+        simulation_options=sim_opts,
+        execution_options=exec_opts,
+    )
+
+    if rec_symbol == :p_rms
+        reshaped = np.reshape(data["p_rms"], (nx, ny, nz), order="F")
+        return pc.pyconvert(Array{Float64, 3}, reshaped)
+    end
+
+    reshaped = np.reshape(data["p"], (nx, ny, nz, nt), order="F")
+    return pc.pyconvert(Array{Float64, 4}, reshaped)
+end
+
+function _relative_delays_from_phase_matrix(phases::AbstractMatrix{<:Real}, fc::Real)
+    unwrapped = _unwrap_phase_matrix(phases)
+    travel_like = .-unwrapped ./ (2π * Float64(fc))
+    return maximum(travel_like) .- travel_like
+end
+
+function _outward_kwave_nt_3d(
+    c::AbstractArray{<:Real, 3},
+    cfg::SimulationConfig3D,
+    kgrid,
+    col_range_y::UnitRange{Int},
+    col_range_z::UnitRange{Int},
+)
+    y_aperture = collect(kgrid.y)[col_range_y]
+    z_aperture = collect(kgrid.z)[col_range_z]
+    max_dy = maximum(abs.(y_aperture .- cfg.lateral_y_focus))
+    max_dz = maximum(abs.(z_aperture .- cfg.lateral_z_focus))
+    max_distance = sqrt(cfg.axial_focus^2 + max_dy^2 + max_dz^2)
+    c_min = max(minimum(Float64.(c)), eps(Float64))
+    required_t = max_distance / c_min + _OUTWARD_KWAVE_NUM_CYCLES / cfg.fc + _OUTWARD_KWAVE_TIME_MARGIN
+    return clamp(ceil(Int, required_t / cfg.dt), 1, kgrid.Nt)
+end
+
+function _outward_kwave_focus_3d(
+    c::AbstractArray{<:Real, 3},
+    rho::AbstractArray{<:Real, 3},
+    cfg::SimulationConfig3D,
+    kgrid,
+    col_range_y::UnitRange{Int},
+    col_range_z::UnitRange{Int};
+    use_gpu::Bool=false,
+)
+    mods = _kwave_modules()
+    pc = mods.PythonCall
+    np = mods.np
+    deepcopy_py = mods.copy.deepcopy
+
+    nx, ny, nz = size(c)
+    nt = _outward_kwave_nt_3d(c, cfg, kgrid, col_range_y, col_range_z)
+    pml_guard = max(4, cfg.PML_GUARD)
+
+    py_kgrid = mods.kgrid.kWaveGrid(pc.Py([nx, ny, nz]), pc.Py([cfg.dx, cfg.dy, cfg.dz]))
+    py_kgrid.dt = cfg.dt
+    py_kgrid.Nt = nt
+
+    medium = mods.kmedium.kWaveMedium(
+        sound_speed=np.array(pc.Py(Float32.(c)), dtype=np.float32),
+        density=np.array(pc.Py(Float32.(rho)), dtype=np.float32),
+    )
+
+    source = mods.ksource.kSource()
+    src_mask_jl = falses(nx, ny, nz)
+    src_row = target_index(cfg)
+    src_col_y = argmin(abs.(collect(kgrid.y) .- cfg.lateral_y_focus))
+    src_col_z = argmin(abs.(collect(kgrid.z) .- cfg.lateral_z_focus))
+    src_mask_jl[src_row, src_col_y, src_col_z] = true
+    source.p_mask = np.array(pc.Py(Array(src_mask_jl)), dtype=pc.pybuiltins.bool)
+    zero_offset_py = np.array(pc.Py(Int64[0]), dtype=np.int64)
+    burst_py = mods.signals.tone_burst(
+        sample_freq=1 / cfg.dt,
+        signal_freq=cfg.fc,
+        num_cycles=_OUTWARD_KWAVE_NUM_CYCLES,
+        signal_offset=zero_offset_py,
+    )
+    burst = pc.pyconvert(Matrix{Float64}, burst_py)
+    source.p = np.array(pc.Py(reshape(vec(burst), 1, :)), dtype=np.float64)
+    source.p_frequency_ref = cfg.fc
+    source.medium = medium
+
+    sensor = mods.ksensor.kSensor()
+    sensor_mask_jl = falses(nx, ny, nz)
+    sensor_mask_jl[cfg.trans_index, col_range_y, col_range_z] .= true
+    sensor.mask = np.array(pc.Py(Array(sensor_mask_jl)), dtype=pc.pybuiltins.bool)
+    sensor.record = pc.pybuiltins.list(("p",))
+
+    sim_opts = mods.simopts.SimulationOptions(
+        pml_inside=true,
+        pml_size=pml_guard,
+        data_recast=false,
+        save_to_disk=true,
+    )
+    exec_opts = mods.execopts.SimulationExecutionOptions(
+        is_gpu_simulation=use_gpu,
+        delete_data=true,
+    )
+
+    data = mods.kspace3d.kspaceFirstOrder3D(
+        kgrid=deepcopy_py(py_kgrid),
+        medium=deepcopy_py(medium),
+        source=deepcopy_py(source),
+        sensor=deepcopy_py(sensor),
+        simulation_options=sim_opts,
+        execution_options=exec_opts,
+    )
+    raw = pc.pyconvert(Matrix{Float64}, np.array(data["p"], dtype=np.float64))
+    n_sensors = length(col_range_y) * length(col_range_z)
+    sensor_data = if size(raw, 1) == n_sensors && size(raw, 2) == nt
+        raw
+    elseif size(raw, 1) == nt && size(raw, 2) == n_sensors
+        permutedims(raw)
+    else
+        error("Unexpected 3D outward sensor data shape $(size(raw)); expected ($n_sensors, $nt) or ($nt, $n_sensors).")
+    end
+    phasors = reshape(_sensor_phasors(sensor_data, cfg.dt, cfg.fc), length(col_range_y), length(col_range_z))
+    tau = _relative_delays_from_phase_matrix(angle.(phasors), cfg.fc)
+    amplitudes = abs.(phasors)
+    amplitudes ./= max(maximum(amplitudes), eps(Float64))
+    return tau, amplitudes, Dict{Symbol, Any}(
+        :solver => :kwave,
+        :source_index => (src_row, src_col_y, src_col_z),
+        :nt => nt,
+        :t_max => nt * cfg.dt,
+        :phasors => phasors,
+    )
+end
+
+function _focus_impl_3d(
+    c::AbstractArray{<:Real, 3},
+    rho::AbstractArray{<:Real, 3},
+    est_type::Est,
+    cfg::SimulationConfig3D,
+    sweep_settings::SweepSettings,
+    animation_settings::Union{Nothing, AnimationSettings},
+)
+    size(c) == size(rho) || error("Density volume must have the same size as the sound-speed volume.")
+    size(c) == (Nx(cfg), Ny(cfg), Nz(cfg)) || error("Sound-speed volume size $(size(c)) does not match SimulationConfig3D ($(Nx(cfg)), $(Ny(cfg)), $(Nz(cfg))).")
+    focus_idx = target_index(cfg)
+    focus_idx <= Nx(cfg) || error("Focus lies outside the computational grid.")
+    cfg.trans_index < focus_idx || error("Focus must be deeper than the transducer plane.")
+
+    nt = isnothing(animation_settings) ? Nt(cfg) : animation_settings.Nt
+    kgrid = focus_grid_3d(cfg; nt=nt)
+
+    row_range = cfg.trans_index:focus_idx
+    c_hasa = Float64.(reverse(c[row_range, :, :]; dims=1))
+    nprop = size(c_hasa, 1)
+    y = collect(kgrid.y)
+    z = collect(kgrid.z)
+
+    ky = reshape(_focus_wavenumbers(length(y), cfg.dy), :, 1)
+    kz = reshape(_focus_wavenumbers(length(z), cfg.dz), 1, :)
+    k0 = omega(cfg) / cfg.c0
+    kx = sqrt.(complex.(k0^2 .- ky .^ 2 .- kz .^ 2, 0.0))
+    real_mask = real.(kx ./ k0) .> 0.0
+    weighting = Float64.(real_mask)
+    inv_2ikx = zeros(ComplexF64, size(kx))
+    inv_2ikx[real_mask] .= 1.0 ./ (2im .* kx[real_mask])
+
+    p0 = _focus_seed_3d(y, z, cfg)
+    p0_hat = _fftshift2(fft(p0))
+    p_hat = Array{ComplexF64}(undef, nprop, length(y), length(z))
+    p_hat[1, :, :] .= p0_hat
+    propagator = exp.(1im .* kx .* cfg.dx)
+
+    c_prime = c_hasa .- cfg.c0
+    mu = (1 .+ c_prime ./ cfg.c0) .^ -2
+    lam = (k0^2) .* (1 .- mu)
+
+    current = copy(p0_hat)
+    @views for xi in 2:nprop
+        pz = ifft(_ifftshift2(current))
+        pz_hat = _fftshift2(fft(lam[xi, :, :] .* pz))
+        p1 = current .* propagator .+ propagator .* inv_2ikx .* pz_hat .* cfg.dx
+        p1[.!real_mask] .= 0.0
+        p1 .*= weighting
+        p_hat[xi, :, :] .= p1
+        current = p1
+    end
+
+    p_field = Array{ComplexF64}(undef, nprop, length(y), length(z))
+    @views for xi in 1:nprop
+        p_field[xi, :, :] .= ifft(_ifftshift2(p_hat[xi, :, :]))
+    end
+    trans_plane = p_field[end, :, :]
+    phase_plane = _unwrap_phase_matrix(angle.(trans_plane); discont=π)
+    amp_plane = abs.(trans_plane)
+    amplitudes = amp_plane ./ max(maximum(amp_plane), eps(Float64))
+
+    delays = phase_plane ./ omega(cfg)
+    tau_kw_h = maximum(delays) .- delays
+    delays_geo = Matrix{Float64}(undef, length(y), length(z))
+    @inbounds for iz in eachindex(z), iy in eachindex(y)
+        delays_geo[iy, iz] = sqrt(
+            cfg.axial_focus^2 +
+            (y[iy] - cfg.lateral_y_focus)^2 +
+            (z[iz] - cfg.lateral_z_focus)^2,
+        ) / cfg.c0
+    end
+    tau_kw_geo = maximum(delays_geo) .- delays_geo
+
+    outward_solver = _validate_outward_solver(sweep_settings.outward_solver)
+    outward_info = Dict{Symbol, Any}(:solver => outward_solver)
+
+    range_y = active_col_range_y(cfg)
+    range_z = active_col_range_z(cfg)
+
+    if outward_solver == :kwave && est_type != GEOMETRIC
+        tau_active, amp_active, outward_info = _outward_kwave_focus_3d(
+            c,
+            rho,
+            cfg,
+            kgrid,
+            range_y,
+            range_z;
+            use_gpu=sweep_settings.use_gpu,
+        )
+    elseif est_type == GEOMETRIC
+        tau = copy(tau_kw_geo)
+        amplitudes .= 1.0
+        tau_active = tau[range_y, range_z]
+        amp_active = amplitudes[range_y, range_z]
+    else
+        tau = sweep_settings.use_hasa_phase ? copy(tau_kw_h) : copy(tau_kw_geo)
+        sweep_settings.use_hasa_amplitude || (amplitudes .= 1.0)
+        tau_active = tau[range_y, range_z]
+        amp_active = amplitudes[range_y, range_z]
+    end
+
+    all(isfinite, tau_active) || error("3D focus produced non-finite transmit delays.")
+    all(isfinite, amp_active) || error("3D focus produced non-finite transmit amplitudes.")
+
+    pressure = nothing
+    run_kwave = sweep_settings.verify_inward_kwave && (isnothing(animation_settings) || animation_settings.run_kwave)
+    if run_kwave
+        pressure = _simulate_kwave_3d(
+            c,
+            rho,
+            tau_active,
+            amp_active,
+            cfg,
+            kgrid,
+            range_y,
+            range_z;
+            record=sweep_settings.record,
+            use_gpu=sweep_settings.use_gpu,
+            animation_settings=animation_settings,
+        )
+    end
+
+    hasa_info = Dict{Symbol, Any}(
+        :p_hasa => p_field,
+        :c_hasa => c_hasa,
+        :tau => tau_active,
+        :tau_geo => tau_kw_geo[range_y, range_z],
+        :amplitudes => amp_active,
+        :transducer_cols_y => range_y,
+        :transducer_cols_z => range_z,
+        :outward => outward_info,
+    )
+    return pressure, hasa_info, kgrid
+end
+
+function focus(
+    c::AbstractArray{<:Real, 3},
+    rho::AbstractArray{<:Real, 3},
+    est_type::Est,
+    cfg::SimulationConfig3D,
+    sweep_settings::SweepSettings,
+    animation_settings::Union{Nothing, AnimationSettings},
+)
+    return _focus_impl_3d(c, rho, est_type, cfg, sweep_settings, animation_settings)
+end
+
+function focus(
+    c::AbstractArray{<:Real, 3},
+    rho::AbstractArray{<:Real, 3},
+    est_type::Est,
+    cfg::SimulationConfig3D,
+    sweep_settings::SweepSettings;
+    animation_settings::Union{Nothing, AnimationSettings}=nothing,
+)
+    return _focus_impl_3d(c, rho, est_type, cfg, sweep_settings, animation_settings)
+end
+
+function skull_mask_from_c_columnwise_3d(
+    c::AbstractArray{<:Real, 3};
+    c_water::Real=1500.0,
+    tol::Real=5.0,
+    min_thick_rows::Integer=2,
+    mask_outside::Bool=true,
+)
+    nx, ny, nz = size(c)
+    diff_mask = abs.(Float64.(c) .- Float64(c_water)) .> Float64(tol)
+    skull_mask = falses(nx, ny, nz)
+    @inbounds for iz in 1:nz, iy in 1:ny
+        rows = findall(diff_mask[:, iy, iz])
+        isempty(rows) && continue
+        first_row = first(rows)
+        last_row = last(rows)
+        if last_row - first_row + 1 >= min_thick_rows
+            skull_mask[first_row:last_row, iy, iz] .= true
+            if mask_outside
+                skull_mask[last_row:end, iy, iz] .= true
+            end
+        end
+    end
+    return skull_mask
+end
+
+function analyse_focus_3d(
+    p::AbstractArray{<:Real, 3},
+    kgrid,
+    cfg::SimulationConfig3D;
+    thr_db::Real=3.0,
+    exclude_mask::Union{Nothing, AbstractArray{Bool, 3}}=nothing,
+    exclude_skull::Bool=false,
+)
+    focus_idx = target_index(cfg)
+    x_focus_global = kgrid.x[focus_idx]
+    y_focus_global = cfg.lateral_y_focus
+    z_focus_global = cfg.lateral_z_focus
+    rows = (cfg.trans_index + 1):size(p, 1)
+    p_roi = Float64.(p[rows, :, :])
+    if exclude_skull && !isnothing(exclude_mask)
+        p_roi = ifelse.(exclude_mask[rows, :, :], -Inf, p_roi)
+    end
+
+    idx_local = Tuple(argmax(p_roi))
+    idx_peak = (first(rows) + idx_local[1] - 1, idx_local[2], idx_local[3])
+    p_peak = Float64(p[idx_peak...])
+    x_pk = kgrid.x[idx_peak[1]]
+    y_pk = kgrid.y[idx_peak[2]]
+    z_pk = kgrid.z[idx_peak[3]]
+    idx_peak_global = Tuple(argmax(p))
+
+    d_axial_mm = (x_focus_global - x_pk) * 1e3
+    d_y_mm = (y_focus_global - y_pk) * 1e3
+    d_z_mm = (z_focus_global - z_pk) * 1e3
+    error_mm = sqrt(d_axial_mm^2 + d_y_mm^2 + d_z_mm^2)
+
+    thr = p_peak / (10^(thr_db / 20))
+    mask = p_roi .>= thr
+    volume_mm3 = count(mask) * cfg.dx * cfg.dy * cfg.dz * 1e9
+    rows_masked = findall(vec(any(mask; dims=(2, 3))))
+    axial_len_mm = isempty(rows_masked) ? 0.0 :
+        (maximum(kgrid.x[first(rows) .+ rows_masked .- 1]) - minimum(kgrid.x[first(rows) .+ rows_masked .- 1])) * 1e3
+    y_cols = findall(mask[idx_local[1], :, idx_local[3]])
+    z_cols = findall(mask[idx_local[1], idx_local[2], :])
+    y_diam_mm = isempty(y_cols) ? 0.0 : (maximum(kgrid.y[y_cols]) - minimum(kgrid.y[y_cols])) * 1e3
+    z_diam_mm = isempty(z_cols) ? 0.0 : (maximum(kgrid.z[z_cols]) - minimum(kgrid.z[z_cols])) * 1e3
+
+    return Dict{Symbol, Any}(
+        :p_peak => p_peak,
+        :error_mm => error_mm,
+        :d_axial_mm => d_axial_mm,
+        :d_y_mm => d_y_mm,
+        :d_z_mm => d_z_mm,
+        :focal_volume_mm3 => volume_mm3,
+        :axial_len_mm => axial_len_mm,
+        :y_diam_mm => y_diam_mm,
+        :z_diam_mm => z_diam_mm,
+        :peak_mm => (x_pk * 1e3, y_pk * 1e3, z_pk * 1e3),
+        :is_global_peak => idx_peak == idx_peak_global,
+    )
+end
+
+function make_focus_medium_3d(
+    hu_vol::Union{Nothing, AbstractArray{<:Real, 3}},
+    spacing_m::Union{Nothing, NTuple{3, <:Real}},
+    cfg::SimulationConfig3D,
+    medium_type::MediumType;
+    ct_path::AbstractString=DEFAULT_CT_PATH,
+    slice_index_z::Integer=250,
+    hu_bone_thr::Integer=200,
+)
+    aberrator = medium_type == WATER ? :water : medium_type == SKULL_IN_WATER ? :skull : error("Unsupported medium type: $medium_type")
+    c, rho, info = make_pam_medium_3d(
+        focus_pam_config_3d(cfg);
+        aberrator=aberrator,
+        hu_vol=hu_vol,
+        spacing_m=spacing_m,
+        ct_path=ct_path,
+        slice_index_z=slice_index_z,
+        skull_to_transducer=cfg.trans_skull_dist,
+        hu_bone_thr=hu_bone_thr,
+    )
+    return c, rho, info
+end
+
+function run_focus_case_3d(
+    hu_vol::Union{Nothing, AbstractArray{<:Real, 3}},
+    spacing_m::Union{Nothing, NTuple{3, <:Real}},
+    cfg::SimulationConfig3D,
+    medium_type::MediumType,
+    est_type::Est,
+    sweep_settings::SweepSettings;
+    ct_path::AbstractString=DEFAULT_CT_PATH,
+    slice_index_z::Integer=250,
+    hu_bone_thr::Integer=200,
+    animation_settings::Union{Nothing, AnimationSettings}=nothing,
+    return_c::Bool=false,
+)
+    c, rho, info = make_focus_medium_3d(
+        hu_vol,
+        spacing_m,
+        cfg,
+        medium_type;
+        ct_path=ct_path,
+        slice_index_z=slice_index_z,
+        hu_bone_thr=hu_bone_thr,
+    )
+
+    pressure, hasa_info, kgrid = focus(c, rho, est_type, cfg, sweep_settings, animation_settings)
+    stats = nothing
+    if pressure !== nothing && ndims(pressure) == 3
+        exclude_mask = nothing
+        if sweep_settings.exclude_skull
+            exclude_mask = skull_mask_from_c_columnwise_3d(
+                c;
+                c_water=cfg.c0,
+                mask_outside=sweep_settings.mask_outside,
+            )
+        end
+        stats = analyse_focus_3d(
+            pressure,
+            kgrid,
+            cfg;
+            exclude_mask=exclude_mask,
+            exclude_skull=sweep_settings.exclude_skull,
+        )
+    end
+
+    return stats, pressure, return_c ? c : nothing, kgrid, cfg, hasa_info, info
+end
diff --git a/src/focus/focus.jl b/src/focus/focus.jl
index aaf7f80..5a7cdd2 100644
--- a/src/focus/focus.jl
+++ b/src/focus/focus.jl
@@ -18,6 +18,9 @@ end
 @enum MediumType WATER SKULL_IN_WATER
 @enum Est GEOMETRIC HASA
 
+const _OUTWARD_KWAVE_NUM_CYCLES = 5
+const _OUTWARD_KWAVE_TIME_MARGIN = 5e-6
+
 function parse_placement_mode(s::AbstractString)
     norm = replace(lowercase(s), "-" => "_")
     norm in ("auto", "fixed_transducer", "fixed_focus_depth") || error("Unknown placement mode: $s")
@@ -154,6 +157,8 @@ Base.@kwdef struct SweepSettings
     exclude_skull::Bool = true
     mask_outside::Bool = true
     record::Symbol = :p_rms
+    outward_solver::Symbol = :angular_spectrum
+    verify_inward_kwave::Bool = true
 end
 
 omega(cfg::SimulationConfig) = 2π * cfg.fc
@@ -256,6 +261,119 @@ function plot_hasa_results(
     return fig
 end
 
+function _validate_outward_solver(solver::Symbol)
+    solver in (:angular_spectrum, :kwave) || error("Unsupported outward propagation solver: $solver")
+    return solver
+end
+
+function _kwave_phase_reference(nt::Int, dt::Real, fc::Real)
+    t = collect(0:(nt - 1)) .* Float64(dt)
+    return exp.(-1im .* (2π * Float64(fc)) .* t)
+end
+
+function _sensor_phasors(sensor_data::AbstractMatrix{<:Real}, dt::Real, fc::Real)
+    ref = _kwave_phase_reference(size(sensor_data, 2), dt, fc)
+    return vec(Matrix{Float64}(sensor_data) * ref)
+end
+
+function _outward_kwave_nt_2d(c::AbstractMatrix{<:Real}, cfg::SimulationConfig, kgrid::KGrid2D, col_range::UnitRange{Int})
+    y_aperture = kgrid.y_vec[col_range]
+    max_distance = sqrt(cfg.z_focus^2 + maximum(abs.(y_aperture .- cfg.x_focus))^2)
+    c_min = max(minimum(Float64.(c)), eps(Float64))
+    required_t = max_distance / c_min + _OUTWARD_KWAVE_NUM_CYCLES / cfg.fc + _OUTWARD_KWAVE_TIME_MARGIN
+    return clamp(ceil(Int, required_t / cfg.dt), 1, kgrid.Nt)
+end
+
+function _relative_delays_from_phase(phases::AbstractVector{<:Real}, fc::Real)
+    unwrapped = _unwrap_phase(phases)
+    travel_like = .-unwrapped ./ (2π * Float64(fc))
+    return maximum(travel_like) .- travel_like
+end
+
+function _outward_kwave_focus_2d(
+    c::AbstractMatrix{<:Real},
+    rho::AbstractMatrix{<:Real},
+    cfg::SimulationConfig,
+    kgrid::KGrid2D,
+    col_range::UnitRange{Int};
+    use_gpu::Bool=false,
+)
+    mods = _kwave_modules()
+    pc = mods.PythonCall
+    np = mods.np
+    deepcopy_py = mods.copy.deepcopy
+
+    nx, ny = size(c)
+    nt = _outward_kwave_nt_2d(c, cfg, kgrid, col_range)
+    pml_guard = max(4, cfg.PML_GUARD)
+
+    py_kgrid = mods.kgrid.kWaveGrid(pc.Py([nx, ny]), pc.Py([cfg.dx, cfg.dz]))
+    py_kgrid.dt = cfg.dt
+    py_kgrid.Nt = nt
+
+    medium = mods.kmedium.kWaveMedium(
+        sound_speed=np.array(pc.Py(Float32.(c)), dtype=np.float32),
+        density=np.array(pc.Py(Float32.(rho)), dtype=np.float32),
+    )
+
+    source = mods.ksource.kSource()
+    src_mask = _py_bool_matrix(np, nx, ny)
+    src_row = target_index(cfg)
+    src_col = argmin(abs.(kgrid.y_vec .- cfg.x_focus))
+    src_mask[src_row - 1, src_col - 1] = true
+    source.p_mask = src_mask
+    zero_offset_py = np.array(pc.Py(Int64[0]), dtype=np.int64)
+    burst_py = mods.signals.tone_burst(
+        sample_freq=1 / cfg.dt,
+        signal_freq=cfg.fc,
+        num_cycles=_OUTWARD_KWAVE_NUM_CYCLES,
+        signal_offset=zero_offset_py,
+    )
+    burst = pc.pyconvert(Matrix{Float64}, burst_py)
+    source.p = np.array(pc.Py(reshape(vec(burst), 1, :)), dtype=np.float64)
+    source.p_frequency_ref = cfg.fc
+    source.medium = medium
+
+    sensor = mods.ksensor.kSensor()
+    sensor_mask = _py_bool_matrix(np, nx, ny)
+    sensor_mask[cfg.trans_index - 1, (first(col_range) - 1):(last(col_range) - 1)] = true
+    sensor.mask = sensor_mask
+    sensor.record = pc.pybuiltins.list(("p",))
+
+    sim_opts = mods.simopts.SimulationOptions(
+        pml_inside=true,
+        pml_size=pml_guard,
+        data_recast=false,
+        save_to_disk=true,
+    )
+    exec_opts = mods.execopts.SimulationExecutionOptions(
+        is_gpu_simulation=use_gpu,
+        delete_data=true,
+    )
+
+    data = mods.kspace.kspaceFirstOrder2D(
+        kgrid=deepcopy_py(py_kgrid),
+        medium=deepcopy_py(medium),
+        source=deepcopy_py(source),
+        sensor=deepcopy_py(sensor),
+        simulation_options=sim_opts,
+        execution_options=exec_opts,
+    )
+    raw = pc.pyconvert(Array, np.array(data["p"], dtype=np.float64))
+    sensor_data = _as_sensor_matrix(raw, length(col_range), nt)
+    phasors = _sensor_phasors(sensor_data, cfg.dt, cfg.fc)
+    tau = _relative_delays_from_phase(angle.(phasors), cfg.fc)
+    amplitudes = abs.(phasors)
+    amplitudes ./= max(maximum(amplitudes), eps(Float64))
+    return tau, amplitudes, Dict{Symbol, Any}(
+        :solver => :kwave,
+        :source_index => (src_row, src_col),
+        :nt => nt,
+        :t_max => nt * cfg.dt,
+        :phasors => phasors,
+    )
+end
+
 function _focus_impl(
     c::AbstractMatrix{<:Real},
     rho::AbstractMatrix{<:Real},
@@ -360,22 +478,40 @@ function _focus_impl(
 
     fig = sweep_settings.generate_figure ? plot_hasa_results(z, x, cfg, asamap, delays, delays_h, amplitudes) : nothing
 
-    if est_type == GEOMETRIC
+    outward_solver = _validate_outward_solver(sweep_settings.outward_solver)
+    outward_info = Dict{Symbol, Any}(:solver => outward_solver)
+
+    if outward_solver == :kwave && est_type != GEOMETRIC
+        col_range = active_col_range(cfg)
+        tau, amplitudes, outward_info = _outward_kwave_focus_2d(
+            c,
+            rho,
+            cfg,
+            kgrid,
+            col_range;
+            use_gpu=sweep_settings.use_gpu,
+        )
+    elseif est_type == GEOMETRIC
         tau = copy(tau_kw_geo)
         amplitudes = ones(length(amplitudes))
+        col_range = active_col_range(cfg)
+        tau = tau[col_range]
+        amplitudes = amplitudes[col_range]
     else
         tau = sweep_settings.use_hasa_phase ? copy(tau_kw_h) : copy(tau_kw_geo)
         if !sweep_settings.use_hasa_amplitude
             amplitudes = ones(length(amplitudes))
         end
+        col_range = active_col_range(cfg)
+        tau = tau[col_range]
+        amplitudes = amplitudes[col_range]
     end
 
-    col_range = active_col_range(cfg)
-    tau = tau[col_range]
-    amplitudes = amplitudes[col_range]
+    all(isfinite, tau) || error("2D focus produced non-finite transmit delays.")
+    all(isfinite, amplitudes) || error("2D focus produced non-finite transmit amplitudes.")
 
     pressure = nothing
-    run_kwave = isnothing(animation_settings) || animation_settings.run_kwave
+    run_kwave = sweep_settings.verify_inward_kwave && (isnothing(animation_settings) || animation_settings.run_kwave)
     if run_kwave
         pressure = _simulate_kwave(
             c,
@@ -398,6 +534,8 @@ function _focus_impl(
         :tau => tau,
         :tau_geo => tau_kw_geo[col_range],
         :amplitudes => amplitudes,
+        :transducer_cols => col_range,
+        :outward => outward_info,
     )
     return pressure, hasa_info, kgrid
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index adafa6b..17d046e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -74,6 +74,59 @@ end
     tau = hasa_info[:tau]
     @test tau ≈ reverse(tau) atol=1e-9
     @test all(hasa_info[:amplitudes] .≈ 1.0)
+    kgrid_focus = KGrid2D(Nx(cfg), Nz(cfg), cfg.dx, cfg.dz; dt=cfg.dt, Nt=Nt(cfg))
+    @test TranscranialFUS._outward_kwave_nt_2d(c, cfg, kgrid_focus, active_col_range(cfg)) < Nt(cfg)
+end
+
+@testset "3D geometric delay symmetry" begin
+    cfg = SimulationConfig3D(
+        axial_focus=0.012,
+        dx=1e-3,
+        dy=1e-3,
+        dz=1e-3,
+        transverse_dim_y=0.021,
+        transverse_dim_z=0.021,
+        trans_aperture_y=0.011,
+        trans_aperture_z=0.011,
+        axial_padding=1.5,
+        PML_GUARD=3,
+    )
+    c = fill(Float32(cfg.c0), Nx(cfg), TranscranialFUS.Ny(cfg), Nz(cfg))
+    rho = fill(Float32(cfg.rho0), size(c))
+    _, hasa_info, _ = focus(
+        c,
+        rho,
+        GEOMETRIC,
+        cfg,
+        SweepSettings();
+        animation_settings=AnimationSettings(run_kwave=false, Nt=80),
+    )
+    tau = hasa_info[:tau]
+    @test tau ≈ reverse(tau; dims=1) atol=1e-9
+    @test tau ≈ reverse(tau; dims=2) atol=1e-9
+    @test all(hasa_info[:amplitudes] .≈ 1.0)
+    @test TranscranialFUS._outward_kwave_nt_3d(c, cfg, focus_grid_3d(cfg), active_col_range_y(cfg), active_col_range_z(cfg)) < Nt(cfg)
+
+    c[5:8, :, :] .= 2200.0f0
+    _, hasa_info_hasa, _ = focus(
+        c,
+        rho,
+        HASA,
+        cfg,
+        SweepSettings();
+        animation_settings=AnimationSettings(run_kwave=false, Nt=80),
+    )
+    @test all(isfinite, hasa_info_hasa[:tau])
+    @test all(isfinite, hasa_info_hasa[:amplitudes])
+
+    pressure_skip, _, _ = focus(
+        c,
+        rho,
+        GEOMETRIC,
+        cfg,
+        SweepSettings(verify_inward_kwave=false),
+    )
+    @test pressure_skip === nothing
 end
 
 @testset "Placement resolution" begin
@@ -113,6 +166,7 @@ end
 end
 
 include("common/kwave_wrapper.jl")
+include("scripts/run_focus.jl")
 include("scripts/run_pam.jl")
 
 @testset "Focus analysis" begin
@@ -134,6 +188,54 @@ include("scripts/run_pam.jl")
     @test stats[:focal_area_mm2] > 0
 end
 
+@testset "3D focus analysis and medium" begin
+    cfg = SimulationConfig3D(
+        axial_focus=0.012,
+        dx=1e-3,
+        dy=1e-3,
+        dz=1e-3,
+        transverse_dim_y=0.021,
+        transverse_dim_z=0.011,
+        trans_aperture_y=0.011,
+        trans_aperture_z=0.007,
+        trans_skull_dist=0.008,
+        axial_padding=2.0,
+        PML_GUARD=3,
+    )
+    pam_cfg = focus_pam_config_3d(cfg)
+    @test receiver_row(pam_cfg) == cfg.trans_index
+
+    hu_vol = synthetic_hu_volume()
+    c, rho, info = make_focus_medium_3d(
+        hu_vol,
+        (1e-3, 1e-3, 1e-3),
+        cfg,
+        SKULL_IN_WATER;
+        slice_index_z=2,
+    )
+    @test size(c) == (Nx(cfg), TranscranialFUS.Ny(cfg), Nz(cfg))
+    @test size(rho) == size(c)
+    @test info[:outer_row] == cfg.trans_index + 8
+
+    kgrid = focus_grid_3d(cfg)
+    row_tgt = target_index(cfg)
+    col_y_tgt = argmin(abs.(collect(kgrid.y) .- cfg.lateral_y_focus))
+    col_z_tgt = argmin(abs.(collect(kgrid.z) .- cfg.lateral_z_focus))
+    p = Array{Float64}(undef, size(c))
+    σ = 1.5e-3
+    for i in 1:size(p, 1), j in 1:size(p, 2), k in 1:size(p, 3)
+        p[i, j, k] = exp(-(
+            (kgrid.x[i] - kgrid.x[row_tgt])^2 +
+            (kgrid.y[j] - kgrid.y[col_y_tgt])^2 +
+            (kgrid.z[k] - kgrid.z[col_z_tgt])^2
+        ) / (2σ^2))
+    end
+    stats = analyse_focus_3d(p, kgrid, cfg)
+    @test stats[:error_mm] < 0.8
+    @test stats[:p_peak] ≈ maximum(p)
+    @test stats[:focal_volume_mm3] > 0
+end
+
 include("pam/runtests.jl")
 
 @testset "k-Wave smoke tests" begin
diff --git a/test/scripts/run_focus.jl b/test/scripts/run_focus.jl
new file mode 100644
index 0000000..6b8fa39
--- /dev/null
+++ b/test/scripts/run_focus.jl
@@ -0,0 +1,74 @@
+@testset "run_focus script entrypoint" begin
+    include("../../scripts/run_focus.jl")
+
+    mktempdir() do dir
+        out2d = joinpath(dir, "focus2d")
+        args2d = [
+            "--dimension=2",
+            "--aberrator=water",
+            "--focus-mm=12:0",
+            "--out-dir=$out2d",
+            "--run-kwave=false",
+        ]
+        dry2d = focus_dry_plan(args2d)
+        @test dry2d[:branch] == :focus2d
+        @test dry2d[:out_dir] == out2d
+        @test dry2d[:aberrator] == :water
+        @test dry2d[:focus_m] == [0.012, 0.0]
+        @test dry2d[:verify_inward_kwave] == false
+        @test dry2d[:outward_propagation] == :angular_spectrum
+        @test main(args2d; dry_run=true) == dry2d
+        smoke2d = [
+            "--dimension=2",
+            "--aberrator=water",
+            "--focus-mm=12:0",
+            "--dx-mm=1",
+            "--dz-mm=1",
+            "--transverse-mm=21",
+            "--receiver-aperture-mm=11",
+            "--out-dir=$(joinpath(dir, "focus2d_smoke"))",
+            "--run-kwave=false",
+        ]
+        @test main(smoke2d) == 0
+        @test isfile(joinpath(dir, "focus2d_smoke", "sensor_timings.png"))
+
+        out3d = joinpath(dir, "focus3d")
+        args3d = [
+            "--dimension=3",
+            "--aberrator=water",
+            "--focus-mm=12:1:-2",
+            "--receiver-aperture-y-mm=11",
+            "--receiver-aperture-z-mm=9",
+            "--outward-propagation=kwave",
+            "--verify-inward-kwave=false",
+            "--out-dir=$out3d",
+        ]
+        dry3d = focus_dry_plan(args3d)
+        @test dry3d[:branch] == :focus3d
+        @test dry3d[:out_dir] == out3d
+        @test dry3d[:focus_m] == [0.012, 0.001, -0.002]
+        @test dry3d[:placement_mode] == :fixed_transducer
+        @test dry3d[:outward_propagation] == :kwave
+        @test dry3d[:verify_inward_kwave] == false
+
+        skull3d = focus_dry_plan([
+            "--dimension=3",
+            "--aberrator=skull",
+            "--estimator=hasa",
+            "--focus-mm=45:0:0",
+            "--slice-index=250",
+            "--skull-transducer-distance-mm=5",
+            "--receiver-aperture-y-mm=20",
+            "--receiver-aperture-z-mm=20",
+            "--run-kwave=false",
+        ])
+        @test skull3d[:branch] == :focus3d
+        @test skull3d[:placement_mode] == :fixed_transducer
+
+        @test_throws ErrorException main(["--focal-cm=6.0"]; dry_run=true)
+        @test_throws ErrorException main(["--medium=water"]; dry_run=true)
+        @test_throws ErrorException main(["--dimension=3", "--placement=fixed_focus_depth"]; dry_run=true)
+        @test_throws ErrorException main(["--dimension=3", "--focus-depth-from-inner-skull-mm=30"]; dry_run=true)
+        @test_throws ErrorException main(["--outward-propagation=bad"]; dry_run=true)
+    end
+end