diff --git a/CLAUDE.md b/CLAUDE.md
new file mode 100644
index 0000000..8274f12
--- /dev/null
+++ b/CLAUDE.md
@@ -0,0 +1,79 @@
+# Horn Simulation Pipeline
+
+Nextflow-orchestrated FEM acoustic horn simulation with Docker-containerised packages.
+
+## Quick reference
+
+```bash
+just                  # list all commands
+just build            # build all Docker images (horn-solver, horn-geometry, horn-analysis)
+just test-local       # run tests locally (horn-core, horn-geometry, horn-analysis)
+just test             # run tests in Docker
+just run              # single mode pipeline (default params)
+just run-auto         # auto mode (all 7 profiles, driver ranking)
+just run-fullauto     # fullauto mode (derives geometry from frequency band)
+just clean            # remove Docker images + Nextflow work dirs
+```
+
+## Running the pipeline
+
+Always use `-profile docker` (the justfile does this automatically). The solver requires `dolfinx` which only exists in the Docker container.
+
+```bash
+# Single mode
+just run --profile conical --throat_radius 0.025 --mouth_radius 0.15 --length 0.3
+
+# Auto mode (mid-horn example)
+just run-auto --target_f_low 250 --target_f_high 6500 --throat_radius 0.025 --mouth_radius 0.2 --length 0.3
+
+# Fullauto mode (geometry derived from frequency band)
+just run-fullauto --target_f_low 500 --target_f_high 4000
+
+# Resume a failed/interrupted run
+just run -resume
+```
+
+## Project structure
+
+```
+main.nf                          # Nextflow pipeline (single/auto/fullauto workflows)
+nextflow.config                  # Docker container mappings per process
+justfile                         # Build, test, run commands
+packages/
+  horn-core/                     # Shared data structures, enums, candidate generation (pure Python, no Docker)
+  horn-geometry/                 # STEP file generation via gmsh (Docker: horn-geometry)
+  horn-solver/                   # FEM solver via dolfinx (Docker: horn-solver)
+  horn-analysis/                 # Plots, reports, scoring, rendering (Docker: horn-analysis)
+```
+
+## Horn profiles
+
+7 flare profiles: `conical`, `exponential`, `hyperbolic`, `tractrix`, `os`, `lecleach`, `cd`
+
+## Testing
+
+```bash
+# Local (fast, no Docker needed — covers horn-core, horn-geometry, horn-analysis)
+just test-local
+
+# Single package in Docker
+just test-package horn-solver
+
+# All packages in Docker
+just test
+```
+
+horn-geometry tests need `gmsh` (installed locally). horn-solver tests need `dolfinx` (Docker only).
+
+## Docker images
+
+Build from repo root (context is `.`):
+```bash
+docker build -t horn-geometry:latest --target production -f ./packages/horn-geometry/Dockerfile .
+```
+
+Note: `gmsh` wheels don't have native `aarch64` Linux builds. On Apple Silicon, build with `--platform linux/amd64`.
+
+## Java for Nextflow
+
+Nextflow requires Java 8-22. If the system Java is too new, source `~/.zshrc` which sets the correct `JAVA_CMD`.
diff --git a/justfile b/justfile
index 53cd1a5..2ae4cea 100644
--- a/justfile
+++ b/justfile
@@ -1,4 +1,5 @@
-packages := "horn-solver horn-geometry horn-analysis"
+docker_packages := "horn-solver horn-geometry horn-analysis"
+local_packages := "horn-core horn-geometry horn-analysis"
 
 # Display help
 default:
@@ -8,30 +9,51 @@ default:
 build:
     #!/usr/bin/env bash
     set -euo pipefail
-    for pkg in {{packages}}; do
-        docker build -t "$pkg:latest" --target production -f "./packages/$pkg/Dockerfile" .
+    for pkg in {{docker_packages}}; do
+        docker build --platform linux/amd64 -t "$pkg:latest" --target production -f "./packages/$pkg/Dockerfile" .
     done
 
-# Run all package tests (build then test)
+# Run all package tests in Docker (build then test)
 test:
     #!/usr/bin/env bash
     set -euo pipefail
-    for pkg in {{packages}}; do
-        docker build -t "$pkg:test" --target test -f "./packages/$pkg/Dockerfile" .
+    for pkg in {{docker_packages}}; do
+        docker build --platform linux/amd64 -t "$pkg:test" --target test -f "./packages/$pkg/Dockerfile" .
     done
-    for pkg in {{packages}}; do
+    for pkg in {{docker_packages}}; do
         echo "Running tests for $pkg..."
         docker run --rm "$pkg:test" pytest "/app/packages/$pkg/tests"
     done
 
-# Build and test a single package: just test-package horn-solver
+# Build and test a single package in Docker: just test-package horn-solver
 test-package pkg:
-    docker build -t "{{pkg}}:test" --target test -f "./packages/{{pkg}}/Dockerfile" .
-    docker run --rm "{{pkg}}:test" pytest "/app/packages/{{pkg}}/tests" -v
+    docker build --platform linux/amd64 -t "{{pkg}}:test" --target test -f "./packages/{{pkg}}/Dockerfile" .
+    docker run --platform linux/amd64 --rm "{{pkg}}:test" pytest "/app/packages/{{pkg}}/tests" -v
 
-# Run the Nextflow pipeline
-run:
-    nextflow run main.nf -profile docker
+# Run tests locally (no Docker — works for horn-core, horn-geometry, horn-analysis)
+test-local:
+    #!/usr/bin/env bash
+    set -euo pipefail
+    echo "Running horn-core tests..."
+    python -m pytest packages/horn-core/tests/ -v
+    echo "Running horn-analysis tests..."
+    python -m pytest packages/horn-analysis/tests/ -v
+    echo "Running horn-geometry tests..."
+    cd packages/horn-geometry && python -m pytest tests/ -v
+
+# Run the Nextflow pipeline (single mode, default params)
+run *ARGS:
+    nextflow run main.nf -profile docker {{ARGS}}
+
+# Run auto mode (unified optimizer): just run-auto --target_f_low 250 --target_f_high 6500
+# Fixed geometry: --mouth_radius 0.15 --length 0.3 (only varies profile)
+# Free geometry: omit mouth_radius/length to derive from frequency band
+run-auto *ARGS:
+    nextflow run main.nf -profile docker --mode auto {{ARGS}}
+
+# Alias for auto mode with all geometry derived (backward compat)
+run-fullauto *ARGS:
+    nextflow run main.nf -profile docker --mode auto {{ARGS}}
 
 # Run Nextflow tests
 test-nextflow:
@@ -42,7 +64,7 @@ clean:
     #!/usr/bin/env bash
     set -euo pipefail
     echo "Cleaning up Docker images..."
-    for pkg in {{packages}}; do
+    for pkg in {{docker_packages}}; do
         docker rmi -f "$pkg:latest" "$pkg:test" || true
     done
     echo "Cleaning up Nextflow files..."
diff --git a/main.nf b/main.nf
index 75878b2..ba14230 100644
--- a/main.nf
+++ b/main.nf
@@ -1,16 +1,17 @@
 #!/usr/bin/env nextflow
 
 // ========================================================================
-// Mode: "single" (default) runs one horn profile; "auto" runs all 3
-// profiles and ranks driver-horn combinations; "fullauto" derives
-// geometry from a target frequency band and explores a grid.
+// Mode: "single" (default) runs one horn profile; "auto" explores a grid
+// of profiles and geometry, ranking driver-horn combinations.
+// "fullauto" is an alias for "auto" with all geometry derived.
 // ========================================================================
 params.mode = "single"
 
-// Horn Geometry
-params.throat_radius = 0.05 // Radius of the horn's throat in meters
-params.mouth_radius = 0.2  // Radius of the horn's mouth in meters
-params.length = 0.5        // Length of the horn in meters
+// Horn Geometry — null means "derive from frequency" in auto mode,
+// or use sensible defaults in single mode
+params.throat_radius = null
+params.mouth_radius = null
+params.length = null
 params.profile = "conical" // Horn flare profile: conical, exponential, hyperbolic, tractrix, os, lecleach, cd
 params.num_sections = 20   // Number of cross-sections for lofting
 
@@ -33,12 +34,13 @@ params.directivity = false
 // Auto mode settings
 params.target_f_low = 500
 params.target_f_high = 4000
-params.drivers_db = "data/drivers.json"
+params.drivers_db = "data/drivers"
 params.top_n = 10
-
-// Fullauto mode settings
-params.num_mouth_radii = 3 // Mouth radius grid points for fullauto
-params.num_lengths = 3     // Length grid points for fullauto
+params.num_mouth_radii = 3 // Mouth radius grid points (when mouth_radius is null)
+params.num_lengths = 3     // Length grid points (when length is null)
+params.lem_top_n = 3       // Number of top candidates to pass from LEM prescreen to FEM
+params.min_diameter = null  // Optional: minimum driver nominal diameter (inches)
+params.max_diameter = null  // Optional: maximum driver nominal diameter (inches)
 
 // ========================================================================
 // Shared processes
@@ -77,6 +79,7 @@ process run_simulation {
     path "results_${band_index}.csv"
 
     script:
+    def sim_length = params.length ?: 0.5
     def band_width = (params.max_freq - params.min_freq) / (params.num_bands as double)
     def min_f = params.min_freq + band_width * band_index
     def max_f = params.min_freq + band_width * (band_index + 1)
@@ -89,7 +92,7 @@ process run_simulation {
         --min-freq ${min_f} \
         --max-freq ${max_f} \
         --num-intervals ${num_intervals_per_band} \
-        --length ${params.length} \
+        --length ${sim_length} \
         --mesh-size ${params.mesh_size} \
         --radiation-model ${params.radiation_model}
     """
@@ -208,26 +211,6 @@ process render_horn_3d {
     """
 }
 
-process render_auto_horn_3d {
-    publishDir "${params.outdir}/auto", mode: 'copy'
-
-    input:
-    val profile
-
-    output:
-    path "horn_3d_${profile}.png"
-
-    script:
-    """
-    python3 -m horn_analysis.horn_render \
-        horn_3d_${profile}.png \
-        --throat-radius ${params.throat_radius} \
-        --mouth-radius ${params.mouth_radius} \
-        --length ${params.length} \
-        --profile ${profile}
-    """
-}
-
 process run_simulation_directivity {
     input:
     tuple path(horn_step), val(band_index)
@@ -236,6 +219,7 @@ process run_simulation_directivity {
     path "directivity_${band_index}.csv"
 
     script:
+    def sim_length = params.length ?: 0.5
     def band_width = (params.max_freq - params.min_freq) / (params.num_bands as double)
     def min_f = params.min_freq + band_width * band_index
     def max_f = params.min_freq + band_width * (band_index + 1)
@@ -248,7 +232,7 @@ process run_simulation_directivity {
         --min-freq ${min_f} \
         --max-freq ${max_f} \
         --num-intervals ${num_intervals_per_band} \
-        --length ${params.length} \
+        --length ${sim_length} \
         --mesh-size ${params.mesh_size} \
         --radiation-model bem \
         --compute-directivity \
@@ -305,13 +289,16 @@ process generate_single_report {
     path "single_report.html"
 
     script:
+    def throat_r = params.throat_radius ?: 0.05
+    def mouth_r = params.mouth_radius ?: 0.2
+    def horn_len = params.length ?: 0.5
     """
     horn-single-report \
         --kpis ${kpis_json} \
         --final-csv ${final_csv} \
-        --throat-radius ${params.throat_radius} \
-        --mouth-radius ${params.mouth_radius} \
-        --length ${params.length} \
+        --throat-radius ${throat_r} \
+        --mouth-radius ${mouth_r} \
+        --length ${horn_len} \
         --profile ${params.profile} \
         --spl-png ${spl_png} \
         --impedance-png ${impedance_png} \
@@ -342,13 +329,16 @@ process generate_single_report_with_directivity {
     path "single_report.html"
 
     script:
+    def throat_r = params.throat_radius ?: 0.05
+    def mouth_r = params.mouth_radius ?: 0.2
+    def horn_len = params.length ?: 0.5
     """
     horn-single-report \
         --kpis ${kpis_json} \
         --final-csv ${final_csv} \
-        --throat-radius ${params.throat_radius} \
-        --mouth-radius ${params.mouth_radius} \
-        --length ${params.length} \
+        --throat-radius ${throat_r} \
+        --mouth-radius ${mouth_r} \
+        --length ${horn_len} \
         --profile ${params.profile} \
         --spl-png ${spl_png} \
         --impedance-png ${impedance_png} \
@@ -364,7 +354,7 @@ process generate_single_report_with_directivity {
 }
 
 // ========================================================================
-// Auto mode processes
+// Unified auto mode processes
 // ========================================================================
 
 process prescreen_drivers {
@@ -373,225 +363,83 @@ process prescreen_drivers {
     input:
     val target_f_low
     val target_f_high
-    val mouth_radius
-    val length
     path drivers_db
 
     output:
     path "prescreen_result.json"
 
     script:
+    def mouth_flag = params.mouth_radius != null ? "--mouth-radius ${params.mouth_radius}" : ""
+    def length_flag = params.length != null ? "--length ${params.length}" : ""
+    def min_dia_flag = params.min_diameter != null ? "--min-diameter ${params.min_diameter}" : ""
+    def max_dia_flag = params.max_diameter != null ? "--max-diameter ${params.max_diameter}" : ""
     """
     python3 -m horn_analysis.prescreen \
         --drivers-db ${drivers_db} \
         --target-f-low ${target_f_low} \
         --target-f-high ${target_f_high} \
-        --mouth-radius ${mouth_radius} \
-        --length ${length} \
+        ${mouth_flag} \
+        ${length_flag} \
+        ${min_dia_flag} \
+        ${max_dia_flag} \
         --output prescreen_result.json
     """
 }
 
-process generate_auto_geometry {
-    input:
-    tuple val(profile), val(throat_radius)
-
-    output:
-    tuple val(profile), path("horn_${profile}.step")
-
-    script:
-    """
-    python3 -m horn_geometry.generator \
-        --throat-radius ${throat_radius} \
-        --mouth-radius ${params.mouth_radius} \
-        --length ${params.length} \
-        --profile ${profile} \
-        --num-sections ${params.num_sections} \
-        --output-file horn_${profile}.step
-    """
-}
-
-process run_auto_simulation {
-    input:
-    tuple val(profile), path(horn_step), val(band_index)
-
-    output:
-    tuple val(profile), path("results_${profile}_${band_index}.csv")
-
-    script:
-    def band_width = (params.max_freq - params.min_freq) / (params.num_bands as double)
-    def min_f = params.min_freq + band_width * band_index
-    def max_f = params.min_freq + band_width * (band_index + 1)
-    def num_intervals_per_band = Math.ceil(params.num_intervals / (params.num_bands as double)) as int
-    """
-    echo "Running ${profile} band ${band_index}: ${min_f} Hz to ${max_f} Hz"
-    python3 -m horn_solver.solver \
-        --step-file ${horn_step} \
-        --output-file results_${profile}_${band_index}.csv \
-        --min-freq ${min_f} \
-        --max-freq ${max_f} \
-        --num-intervals ${num_intervals_per_band} \
-        --length ${params.length} \
-        --mesh-size ${params.mesh_size} \
-        --radiation-model ${params.radiation_model}
-    """
-}
-
-process merge_auto_results {
-    publishDir "${params.outdir}/auto", mode: 'copy'
-
-    input:
-    tuple val(profile), path(csv_files)
-
-    output:
-    tuple val(profile), path("${profile}_results.csv")
-
-    script:
-    """
-    python3 -c "
-import pandas as pd, glob
-files = glob.glob('results_${profile}_*.csv')
-df = pd.concat((pd.read_csv(f) for f in files), ignore_index=True)
-df.sort_values(by='frequency').to_csv('${profile}_results.csv', index=False)
-"
-    """
-}
-
-process score_and_rank {
+process derive_auto_geometry {
     publishDir "${params.outdir}/auto", mode: 'copy'
 
     input:
-    path solver_csvs
     path prescreen_json
-    path drivers_db
 
     output:
-    path "ranked_results.json"
+    tuple path("candidates.csv"), path("design.json")
 
     script:
+    def mouth_flag = params.mouth_radius != null ? "--mouth-radius ${params.mouth_radius}" : ""
+    def length_flag = params.length != null ? "--length ${params.length}" : ""
     """
-    python3 -c "
-import json, glob
-from pathlib import Path
-from horn_drivers.loader import load_drivers
-from horn_analysis.rank_pipeline import rank_horn_drivers
-from horn_analysis.scoring import TargetSpec
-
-prescreen = json.loads(Path('${prescreen_json}').read_text())
-throat_radius = prescreen['throat_radius_m']
-target = TargetSpec(f_low_hz=${params.target_f_low}, f_high_hz=${params.target_f_high})
-
-# Load only pre-screened drivers
-all_drivers = load_drivers('${drivers_db}')
-driver_ids = set(prescreen['drivers'])
-drivers = [d for d in all_drivers if d.driver_id in driver_ids]
-
-all_results = []
-for csv_path in sorted(glob.glob('*_results.csv')):
-    profile = Path(csv_path).stem.replace('_results', '')
-    results = rank_horn_drivers(
-        solver_csv=csv_path,
-        horn_label=profile,
-        throat_radius=throat_radius,
-        drivers=drivers,
-        target=target,
-        top_n=${params.top_n},
-    )
-    all_results.extend(results)
-
-# Sort all by composite score and take overall top N
-all_results.sort(key=lambda r: r['composite_score'], reverse=True)
-all_results = all_results[:${params.top_n}]
-
-Path('ranked_results.json').write_text(json.dumps(all_results, indent=2))
-print(f'Ranked {len(all_results)} driver-horn combinations')
-"
+    python3 -m horn_core.geometry_designer \
+        --target-f-low ${params.target_f_low} \
+        --target-f-high ${params.target_f_high} \
+        --prescreen-json ${prescreen_json} \
+        --num-mouth-radii ${params.num_mouth_radii} \
+        --num-lengths ${params.num_lengths} \
+        ${mouth_flag} \
+        ${length_flag} \
+        --output candidates.csv \
+        --design-json design.json
     """
 }
 
-process generate_auto_report {
+process lem_prescreen {
     publishDir "${params.outdir}/auto", mode: 'copy'
 
     input:
-    path ranked_json
-    path solver_csvs
-    path drivers_db
-    path prescreen_json
-
-    output:
-    path "report/auto_ranking.json"
-    path "report/auto_comparison.png"
-    path "report/auto_summary.txt"
-    path "report/auto_report.html"
-
-    script:
-    """
-    python3 -c "
-import json
-from pathlib import Path
-from horn_drivers.loader import load_drivers
-from horn_analysis.scoring import TargetSpec
-from horn_analysis.auto_report import generate_auto_report
-
-prescreen = json.loads(Path('${prescreen_json}').read_text())
-throat_radius = prescreen['throat_radius_m']
-
-all_ranked = json.loads(Path('${ranked_json}').read_text())
-
-import glob
-solver_csvs = {}
-for csv_path in sorted(glob.glob('*_results.csv')):
-    profile = Path(csv_path).stem.replace('_results', '')
-    solver_csvs[profile] = csv_path
-
-driver_list = load_drivers('${drivers_db}')
-drivers = {d.driver_id: d for d in driver_list}
-
-target = TargetSpec(f_low_hz=${params.target_f_low}, f_high_hz=${params.target_f_high})
-
-generate_auto_report(
-    all_ranked=all_ranked,
-    solver_csvs=solver_csvs,
-    drivers=drivers,
-    throat_radius=throat_radius,
-    target=target,
-    output_dir='report',
-    top_n=5,
-    mouth_radius=${params.mouth_radius},
-    horn_length=${params.length},
-)
-"
-    """
-}
-
-// ========================================================================
-// Fullauto mode processes
-// ========================================================================
-
-process derive_fullauto_geometry {
-    publishDir "${params.outdir}/fullauto", mode: 'copy'
-
-    input:
+    path candidates_csv
     path prescreen_json
+    path drivers_db
+    path design_json
 
     output:
-    path "candidates.csv"
-    path "design.json"
+    tuple path("lem_results.json"), path("lem_filtered_candidates.csv")
 
     script:
     """
-    python3 -m horn_core.geometry_designer \
+    python3 -m horn_analysis.lem_prescreen \
+        --candidates-csv ${candidates_csv} \
+        --prescreen-json ${prescreen_json} \
+        --drivers-db ${drivers_db} \
+        --design-json ${design_json} \
         --target-f-low ${params.target_f_low} \
         --target-f-high ${params.target_f_high} \
-        --prescreen-json ${prescreen_json} \
-        --num-mouth-radii ${params.num_mouth_radii} \
-        --num-lengths ${params.num_lengths} \
-        --output candidates.csv \
-        --design-json design.json
+        --top-n ${params.lem_top_n} \
+        --output lem_results.json \
+        --filtered-csv lem_filtered_candidates.csv
     """
 }
 
-process generate_fullauto_geometry {
+process generate_candidate_geometry {
     input:
     tuple val(candidate_id), val(profile), val(throat_radius), val(mouth_radius), val(length)
 
@@ -610,7 +458,9 @@ process generate_fullauto_geometry {
     """
 }
 
-process run_fullauto_simulation {
+process run_candidate_simulation {
+    errorStrategy 'ignore'
+
     input:
     tuple val(candidate_id), val(profile), val(mouth_radius), val(length), path(horn_step), val(band_index), val(sim_min_freq), val(sim_max_freq)
 
@@ -635,8 +485,8 @@ process run_fullauto_simulation {
     """
 }
 
-process merge_fullauto_results {
-    publishDir "${params.outdir}/fullauto", mode: 'copy'
+process merge_candidate_results {
+    publishDir "${params.outdir}/auto", mode: 'copy'
 
     input:
     tuple val(candidate_id), path(csv_files)
@@ -655,8 +505,8 @@ df.sort_values(by='frequency').to_csv('${candidate_id}_results.csv', index=False
     """
 }
 
-process score_and_rank_fullauto {
-    publishDir "${params.outdir}/fullauto", mode: 'copy'
+process score_and_rank {
+    publishDir "${params.outdir}/auto", mode: 'copy'
 
     input:
     path solver_csvs
@@ -691,14 +541,18 @@ with open('${candidates_csv}') as f:
     for row in csv.DictReader(f):
         candidates_lookup[row['candidate_id']] = row
 
+csv_files = sorted(glob.glob('*_results.csv'))
+total_candidates = len(csv_files)
+
 all_results = []
-for csv_path in sorted(glob.glob('*_results.csv')):
+for csv_path in csv_files:
     candidate_id = Path(csv_path).stem.replace('_results', '')
     cand = candidates_lookup.get(candidate_id, {})
+    cand_throat = float(cand.get('throat_radius', throat_radius))
     results = rank_horn_drivers(
         solver_csv=csv_path,
         horn_label=candidate_id,
-        throat_radius=throat_radius,
+        throat_radius=cand_throat,
         drivers=drivers,
         target=target,
         top_n=${params.top_n},
@@ -707,21 +561,28 @@ for csv_path in sorted(glob.glob('*_results.csv')):
     for r in results:
         r['mouth_radius'] = float(cand.get('mouth_radius', 0))
         r['length'] = float(cand.get('length', 0))
+        r['throat_radius'] = float(cand.get('throat_radius', 0))
         r['profile'] = cand.get('profile', '')
     all_results.extend(results)
 
 # Sort all by composite score and take overall top N
+total_scored = len(all_results)
 all_results.sort(key=lambda r: r['composite_score'], reverse=True)
 all_results = all_results[:${params.top_n}]
 
-Path('ranked_results.json').write_text(json.dumps(all_results, indent=2))
-print(f'Ranked {len(all_results)} driver-horn combinations')
+output = {
+    'total_candidates': total_candidates,
+    'total_scored': total_scored,
+    'results': all_results,
+}
+Path('ranked_results.json').write_text(json.dumps(output, indent=2))
+print(f'Ranked {total_scored} driver-horn combinations ({total_candidates} geometries)')
 "
     """
 }
 
-process generate_fullauto_report {
-    publishDir "${params.outdir}/fullauto", mode: 'copy'
+process generate_auto_report {
+    publishDir "${params.outdir}/auto", mode: 'copy'
 
     input:
     path ranked_json
@@ -729,6 +590,7 @@ process generate_fullauto_report {
     path drivers_db
     path prescreen_json
     path design_json
+    path lem_results_json
 
     output:
     path "report/auto_ranking.json"
@@ -749,7 +611,13 @@ prescreen = json.loads(Path('${prescreen_json}').read_text())
 throat_radius = prescreen['throat_radius_m']
 design = json.loads(Path('${design_json}').read_text())
 
-all_ranked = json.loads(Path('${ranked_json}').read_text())
+ranked_data = json.loads(Path('${ranked_json}').read_text())
+all_ranked = ranked_data['results']
+total_candidates = ranked_data.get('total_candidates', len(all_ranked))
+total_scored = ranked_data.get('total_scored', len(all_ranked))
+
+# Load LEM prescreen results
+lem_results = json.loads(Path('${lem_results_json}').read_text())
 
 solver_csvs = {}
 for csv_path in sorted(glob.glob('*_results.csv')):
@@ -770,6 +638,9 @@ generate_auto_report(
     output_dir='report',
     top_n=5,
     derived_geometry=design,
+    total_candidates=total_candidates,
+    total_scored=total_scored,
+    lem_results=lem_results,
 )
 "
     """
@@ -780,11 +651,16 @@ generate_auto_report(
 // ========================================================================
 
 workflow single {
+    // Apply defaults for single mode when params are null
+    def throat_r = params.throat_radius ?: 0.05
+    def mouth_r = params.mouth_radius ?: 0.2
+    def horn_len = params.length ?: 0.5
+
     // 1. Generate geometry once
     ch_step_file = generate_geometry(
-        params.throat_radius,
-        params.mouth_radius,
-        params.length,
+        throat_r,
+        mouth_r,
+        horn_len,
         params.profile,
         params.num_sections
     )
@@ -816,13 +692,13 @@ workflow single {
 
     // 10. 3D horn geometry render (runs in parallel with simulation)
     render_horn_3d(
-        params.throat_radius,
-        params.mouth_radius,
-        params.length,
+        throat_r,
+        mouth_r,
+        horn_len,
         params.profile
     )
 
-    // 11. Directivity (opt-in, requires BEM) — parallelized across frequency bands
+    // 11. Directivity (opt-in, requires BEM) -- parallelized across frequency bands
     if (params.directivity) {
         ch_dir_band_indices = Channel.from(0..<params.num_bands)
         ch_dir_inputs = ch_step_file.combine(ch_dir_band_indices)
@@ -852,90 +728,42 @@ workflow single {
 
 workflow auto {
     // 1. Pre-screen drivers
-    ch_drivers_db = Channel.fromPath(params.drivers_db)
+    ch_drivers_db = Channel.value(file(params.drivers_db))
     ch_prescreen = prescreen_drivers(
         params.target_f_low,
         params.target_f_high,
-        params.mouth_radius,
-        params.length,
         ch_drivers_db,
     )
 
-    // 2. Generate geometry for each profile using throat radius from prescreen
-    ch_profiles = Channel.from("conical", "exponential", "hyperbolic", "tractrix", "os", "lecleach", "cd")
-
-    // Extract throat radius from prescreen result
-    ch_throat_radius = ch_prescreen.map { json_file ->
-        def data = new groovy.json.JsonSlurper().parse(json_file)
-        data.throat_radius_m
-    }
-
-    // Combine profiles with throat radius
-    ch_geom_inputs = ch_profiles.combine(ch_throat_radius)
-
-    // 3. Generate geometries (3 parallel jobs)
-    ch_geometries = generate_auto_geometry(ch_geom_inputs)
-
-    // 4. Create band indices and combine with geometries
-    ch_band_indices = Channel.from(0..<params.num_bands)
-    ch_sim_inputs = ch_geometries.combine(ch_band_indices)
-
-    // 5. Run simulations (3 profiles x num_bands parallel jobs)
-    ch_band_results = run_auto_simulation(ch_sim_inputs)
-
-    // 6. Group by profile and merge
-    ch_grouped = ch_band_results.groupTuple()
-    ch_merged = merge_auto_results(ch_grouped)
-
-    // 7. Score and rank all combinations
-    ch_all_csvs = ch_merged.map { profile, csv -> csv }.collect()
-    ch_ranked = score_and_rank(
-        ch_all_csvs,
-        ch_prescreen,
-        ch_drivers_db,
-    )
+    // 2. Derive geometry grid from frequency band + prescreen throat radii
+    //    Fixed params (mouth_radius, length) are passed via CLI flags in the process
+    ch_geom_derived = derive_auto_geometry(ch_prescreen)
+    ch_candidates_csv = ch_geom_derived.map { csv, json -> csv }
+    ch_design_json = ch_geom_derived.map { csv, json -> json }
 
-    // 8. Generate report
-    generate_auto_report(
-        ch_ranked,
-        ch_all_csvs,
-        ch_drivers_db,
+    // 3. LEM/Webster prescreening — score all candidates analytically,
+    //    pass only the top N to expensive STEP + FEM stages
+    ch_lem = lem_prescreen(
+        ch_candidates_csv,
         ch_prescreen,
-    )
-
-    // 9. 3D horn geometry renders (one per profile, parallel)
-    ch_render_profiles = Channel.from("conical", "exponential", "hyperbolic", "tractrix", "os", "lecleach", "cd")
-    render_auto_horn_3d(ch_render_profiles)
-}
-
-workflow fullauto {
-    // 1. Pre-screen drivers (no mouth/length needed)
-    ch_drivers_db = Channel.fromPath(params.drivers_db)
-    ch_prescreen = prescreen_drivers(
-        params.target_f_low,
-        params.target_f_high,
-        0,  // placeholder — not used by prescreen filtering
-        0,  // placeholder — not used by prescreen filtering
         ch_drivers_db,
+        ch_design_json,
     )
+    ch_lem_results = ch_lem.map { json, csv -> json }
+    ch_filtered_csv = ch_lem.map { json, csv -> csv }
 
-    // 2. Derive geometry from frequency band + prescreen throat radius
-    ch_geom_derived = derive_fullauto_geometry(ch_prescreen)
-    ch_candidates_csv = ch_geom_derived.map { csv, json -> csv }
-    ch_design_json = ch_geom_derived.map { csv, json -> json }
-
-    // 3. Parse candidates CSV into channel of tuples
-    ch_candidates = ch_candidates_csv
+    // 4. Parse filtered candidates CSV into channel of tuples
+    ch_candidates = ch_filtered_csv
         .splitCsv(header: true)
         .map { row ->
             tuple(row.candidate_id, row.profile, row.throat_radius as double,
                   row.mouth_radius as double, row.length as double)
         }
 
-    // 4. Generate STEP geometry for each candidate
-    ch_geometries = generate_fullauto_geometry(ch_candidates)
+    // 5. Generate STEP geometry for filtered candidates only
+    ch_geometries = generate_candidate_geometry(ch_candidates)
 
-    // 5. Read sim freq range from design.json and combine with band indices
+    // 6. Read sim freq range from design.json and combine with band indices
     ch_sim_range = ch_design_json.map { json_file ->
         def data = new groovy.json.JsonSlurper().parse(json_file)
         tuple(data.sim_freq_range[0] as double, data.sim_freq_range[1] as double)
@@ -946,36 +774,35 @@ workflow fullauto {
         .combine(ch_band_indices)
         .combine(ch_sim_range)
 
-    // 6. Run FEM simulations (candidates x bands)
-    ch_band_results = run_fullauto_simulation(ch_sim_inputs)
+    // 7. Run FEM simulations (filtered candidates x bands)
+    ch_band_results = run_candidate_simulation(ch_sim_inputs)
 
-    // 7. Group by candidate_id and merge bands
+    // 8. Group by candidate_id and merge bands
     ch_grouped = ch_band_results.groupTuple()
-    ch_merged = merge_fullauto_results(ch_grouped)
+    ch_merged = merge_candidate_results(ch_grouped)
 
-    // 8. Score and rank all driver-horn combinations
+    // 9. Score and rank all driver-horn combinations
     ch_all_csvs = ch_merged.map { candidate_id, csv -> csv }.collect()
-    ch_ranked = score_and_rank_fullauto(
+    ch_ranked = score_and_rank(
         ch_all_csvs,
         ch_prescreen,
         ch_drivers_db,
-        ch_candidates_csv,
+        ch_filtered_csv,
     )
 
-    // 9. Generate report with design summary
-    generate_fullauto_report(
+    // 10. Generate report with design summary + LEM stats
+    generate_auto_report(
         ch_ranked,
         ch_all_csvs,
         ch_drivers_db,
         ch_prescreen,
         ch_design_json,
+        ch_lem_results,
     )
 }
 
 workflow {
-    if (params.mode == "fullauto") {
-        fullauto()
-    } else if (params.mode == "auto") {
+    if (params.mode == "fullauto" || params.mode == "auto") {
         auto()
     } else {
         single()
diff --git a/nextflow.config b/nextflow.config
index 6740133..d6d06df 100644
--- a/nextflow.config
+++ b/nextflow.config
@@ -36,18 +36,24 @@ profiles {
             withName: generate_phase_plot {
                 container = 'horn-analysis:latest'
             }
-            // Auto mode processes
+            // Unified auto mode processes
             withName: prescreen_drivers {
                 container = 'horn-analysis:latest'
             }
-            withName: generate_auto_geometry {
+            withName: derive_auto_geometry {
+                container = 'horn-analysis:latest'
+            }
+            withName: lem_prescreen {
+                container = 'horn-analysis:latest'
+            }
+            withName: generate_candidate_geometry {
                 container = 'horn-geometry:latest'
             }
-            withName: run_auto_simulation {
+            withName: run_candidate_simulation {
                 container = 'horn-solver:latest'
                 cpus = 2
             }
-            withName: merge_auto_results {
+            withName: merge_candidate_results {
                 container = 'horn-analysis:latest'
             }
             withName: score_and_rank {
@@ -63,9 +69,6 @@ profiles {
             withName: render_horn_3d {
                 container = 'horn-analysis:latest'
             }
-            withName: render_auto_horn_3d {
-                container = 'horn-analysis:latest'
-            }
             withName: run_simulation_directivity {
                 container = 'horn-bem-solver:latest'
                 cpus = 2
@@ -82,26 +85,6 @@ profiles {
             withName: generate_single_report_with_directivity {
                 container = 'horn-analysis:latest'
             }
-            // Fullauto mode processes
-            withName: derive_fullauto_geometry {
-                container = 'horn-analysis:latest'
-            }
-            withName: generate_fullauto_geometry {
-                container = 'horn-geometry:latest'
-            }
-            withName: run_fullauto_simulation {
-                container = 'horn-solver:latest'
-                cpus = 2
-            }
-            withName: merge_fullauto_results {
-                container = 'horn-analysis:latest'
-            }
-            withName: score_and_rank_fullauto {
-                container = 'horn-analysis:latest'
-            }
-            withName: generate_fullauto_report {
-                container = 'horn-analysis:latest'
-            }
         }
     }
     standard {
diff --git a/packages/horn-analysis/pyproject.toml b/packages/horn-analysis/pyproject.toml
index 4ca4b41..891a9fb 100644
--- a/packages/horn-analysis/pyproject.toml
+++ b/packages/horn-analysis/pyproject.toml
@@ -26,6 +26,7 @@ horn-dashboard = "horn_analysis.dashboard:main"
 horn-render = "horn_analysis.horn_render:main"
 horn-directivity-plot = "horn_analysis.directivity_plot:main"
 horn-single-report = "horn_analysis.single_report:main"
+horn-lem-prescreen = "horn_analysis.lem_prescreen:main"
 
 [project.optional-dependencies]
 test = [
diff --git a/packages/horn-analysis/src/horn_analysis/auto_report.py b/packages/horn-analysis/src/horn_analysis/auto_report.py
index ac93d02..990d630 100644
--- a/packages/horn-analysis/src/horn_analysis/auto_report.py
+++ b/packages/horn-analysis/src/horn_analysis/auto_report.py
@@ -30,6 +30,9 @@ def generate_auto_report(
     mouth_radius: float | None = None,
     horn_length: float | None = None,
     derived_geometry: Optional[dict] = None,
+    total_candidates: int | None = None,
+    total_scored: int | None = None,
+    lem_results: Optional[dict] = None,
 ) -> Path:
     """Generate the auto-select report with rankings, plots, and CSVs.
 
@@ -44,6 +47,8 @@ def generate_auto_report(
         mouth_radius: Horn mouth radius in metres (for report display).
         horn_length: Horn length in metres (for report display).
         derived_geometry: Optional dict from geometry_designer (fullauto mode).
+        total_candidates: Total number of geometry candidates simulated.
+        total_scored: Total number of driver-horn combinations scored.
 
     Returns:
         Path to the output directory.
@@ -74,7 +79,8 @@ def generate_auto_report(
         solver_spl = df["spl"].values
         z_real = df["z_real"].values
         z_imag = df["z_imag"].values
-        throat_area = np.pi * throat_radius ** 2
+        cand_throat = result.get("throat_radius", throat_radius)
+        throat_area = np.pi * cand_throat ** 2
 
         p_throat = compute_driver_response(drv, freq, z_real, z_imag, throat_area)
         coupled_spl = scale_solver_spl(solver_spl, p_throat)
@@ -98,17 +104,30 @@ def generate_auto_report(
         )
 
     # 4. Human-readable summary
+    scored_display = total_scored if total_scored is not None else len(all_ranked)
     lines = [
         "Horn Driver Auto-Select Results",
         "=" * 40,
         f"Target: {target.f_low_hz:.0f} Hz - {target.f_high_hz:.0f} Hz",
         f"Throat radius: {throat_radius:.4f} m",
         f"Profiles evaluated: {', '.join(solver_csvs.keys())}",
-        f"Total candidates scored: {len(all_ranked)}",
+        f"Total candidates scored: {scored_display}",
+    ]
+
+    if lem_results:
+        lines.extend([
+            "",
+            "LEM Prescreening:",
+            f"  Candidates evaluated by LEM: {lem_results.get('total_evaluated', 'N/A')}",
+            f"  Driver-horn pairs scored: {lem_results.get('total_pairs', 'N/A')}",
+            f"  Passed to FEM: {len(lem_results.get('filtered_candidate_ids', []))}",
+        ])
+
+    lines.extend([
         "",
         f"Top {len(top_results)} Results:",
         "-" * 40,
-    ]
+    ])
 
     for rank, result in enumerate(top_results, 1):
         lines.append(
@@ -141,6 +160,9 @@ def generate_auto_report(
         mouth_radius=mouth_radius,
         length=horn_length,
         derived_geometry=derived_geometry,
+        total_candidates=total_candidates,
+        total_scored=total_scored,
+        lem_results=lem_results,
     )
     (out / "auto_report.html").write_text(html_report)
 
diff --git a/packages/horn-analysis/src/horn_analysis/compare.py b/packages/horn-analysis/src/horn_analysis/compare.py
index d7e6d60..f75b39e 100644
--- a/packages/horn-analysis/src/horn_analysis/compare.py
+++ b/packages/horn-analysis/src/horn_analysis/compare.py
@@ -39,12 +39,16 @@ def plot_multi_comparison(
     else:
         fig, ax_plot = plot_theme.create_figure(figsize=(12, 8))
 
+    _linestyles = ["-", "--", "-.", ":", (0, (3, 1, 1, 1))]
+
     all_freq = []
     all_spl = []
     for i, (csv_path, label) in enumerate(file_label_pairs):
         df = pd.read_csv(csv_path)
         color = plot_theme.MULTI_COLORS[i % len(plot_theme.MULTI_COLORS)]
-        ax_plot.plot(df["frequency"], df["spl"], label=label, color=color, linewidth=1.4)
+        ls = _linestyles[i % len(_linestyles)]
+        ax_plot.plot(df["frequency"], df["spl"], label=label,
+                     color=color, linestyle=ls, linewidth=1.8)
         all_freq.extend(df["frequency"].values)
         all_spl.extend(df["spl"].values)
 
diff --git a/packages/horn-analysis/src/horn_analysis/directivity_plot.py b/packages/horn-analysis/src/horn_analysis/directivity_plot.py
index af59281..a4f192d 100644
--- a/packages/horn-analysis/src/horn_analysis/directivity_plot.py
+++ b/packages/horn-analysis/src/horn_analysis/directivity_plot.py
@@ -87,15 +87,30 @@ def plot_polar_directivity(
     ax.set_thetamin(0)
     ax.set_thetamax(180)
 
+    # Determine if data is already relative (on-axis near 0 dB) or absolute
+    db_range = 40  # show 40 dB of dynamic range
+    all_on_axis = df[df["theta_deg"] == df["theta_deg"].min()]["spl_db"]
+    ref_level = all_on_axis.max()
+    # Shift so peak on-axis = 0 dB, then offset so plot radius is non-negative
+    # Display: radius = spl_db - (ref_level - db_range), clipped at 0
+    r_floor = ref_level - db_range
+
     colors = plot_theme.MULTI_COLORS
     for i, freq in enumerate(frequencies):
         sub = df[df["frequency"] == freq].sort_values("theta_deg")
         theta_rad = np.radians(sub["theta_deg"].values)
         spl = sub["spl_db"].values
+        r = np.clip(spl - r_floor, 0, None)
         label = f"{freq:.0f} Hz" if freq < 1000 else f"{freq / 1000:.1f} kHz"
-        ax.plot(theta_rad, spl, color=colors[i % len(colors)], linewidth=1.3, label=label)
+        ax.plot(theta_rad, r, color=colors[i % len(colors)], linewidth=1.3, label=label)
+
+    # Custom radial tick labels showing actual dB values
+    r_ticks = np.linspace(0, db_range, 5)
+    ax.set_rticks(r_ticks)
+    ax.set_yticklabels([f"{v + r_floor:.0f}" for v in r_ticks], fontsize=7)
+    ax.set_rlim(0, db_range + 2)
 
-    ax.set_title("Polar Directivity", pad=20)
+    ax.set_title("Polar Directivity (dB)", pad=20)
     ax.legend(loc="lower left", fontsize=7, bbox_to_anchor=(1.05, 0))
 
     plot_theme.save_figure(fig, output_file)
@@ -148,7 +163,7 @@ def plot_directivity_contour(
     ax.set_ylabel("Angle (degrees)")
     ax.set_title("Directivity Contour")
     cbar = fig.colorbar(mesh, ax=ax, pad=0.02)
-    cbar.set_label("SPL (dB)")
+    cbar.set_label("Relative SPL (dB)")
     plot_theme.setup_grid(ax)
 
     plot_theme.save_figure(fig, output_file)
@@ -265,8 +280,10 @@ def compute_directivity_index(
         theta_deg = sub["theta_deg"].values
         spl = sub["spl_db"].values
 
-        # Convert SPL to linear pressure squared (relative)
-        p_sq = 10 ** (spl / 10)
+        # Convert relative dB to linear pressure squared ratio
+        # spl_db is relative to on-axis (0 dB = on-axis), so
+        # p_sq_ratio = 10^(spl_db / 10) gives p²/p²_on_axis
+        p_sq = 10 ** (spl / 10.0)
 
         theta_rad = np.radians(theta_deg)
 
@@ -275,6 +292,7 @@ def compute_directivity_index(
         p_sq_on = p_sq[idx_on]
 
         # Numerical integration using trapezoidal rule
+        # DI = p²_on / <p²> where <p²> = ∫ p² sin(θ) dθ / ∫ sin(θ) dθ
         integrand = p_sq * np.sin(theta_rad)
         numerator = _trapezoid(integrand, theta_rad)
         denominator = _trapezoid(np.sin(theta_rad), theta_rad)
diff --git a/packages/horn-analysis/src/horn_analysis/horn_render.py b/packages/horn-analysis/src/horn_analysis/horn_render.py
index 0da4d91..526fdab 100644
--- a/packages/horn-analysis/src/horn_analysis/horn_render.py
+++ b/packages/horn-analysis/src/horn_analysis/horn_render.py
@@ -220,6 +220,9 @@ def fig_to_b64_3d(
     n_theta = kwargs.pop("n_theta", 60)
     show_profile = kwargs.pop("show_profile", True)
     figsize = kwargs.pop("figsize", (14, 6))
+    title = kwargs.pop("title", None)
+
+    base_title = title or f"{profile.capitalize()} Horn"
 
     z_vals = np.linspace(0, length, n_z)
     r_vals = _radius_profile(z_vals, throat_radius, mouth_radius, length, profile)
@@ -245,7 +248,7 @@ def fig_to_b64_3d(
     ax3d.set_xlabel("Axial position (m)")
     ax3d.set_ylabel("X (m)")
     ax3d.set_zlabel("Y (m)")
-    ax3d.set_title(f"{profile.capitalize()} Horn — 3D View")
+    ax3d.set_title("3D View", fontsize=9)
     ax3d.view_init(elev=20, azim=-60)
 
     max_range = max(length, 2 * mouth_radius) / 2
@@ -263,10 +266,12 @@ def fig_to_b64_3d(
         ax2d.plot(z_mm, -r_mm, color=plot_theme.COLORS["primary"], linewidth=1.4)
         ax2d.set_xlabel("Axial position (mm)")
         ax2d.set_ylabel("Radius (mm)")
-        ax2d.set_title(f"{profile.capitalize()} Horn — Wall Profile")
+        ax2d.set_title("Wall Profile", fontsize=9)
         ax2d.set_aspect("equal", adjustable="datalim")
         plot_theme.setup_grid(ax2d)
 
+    fig.suptitle(base_title, fontsize=11, fontweight="bold")
+
     fig.tight_layout()
     fig.savefig(buf, format="png", dpi=150, bbox_inches="tight")
     plt.close(fig)
diff --git a/packages/horn-analysis/src/horn_analysis/html_report.py b/packages/horn-analysis/src/horn_analysis/html_report.py
index e72d8e2..7733912 100644
--- a/packages/horn-analysis/src/horn_analysis/html_report.py
+++ b/packages/horn-analysis/src/horn_analysis/html_report.py
@@ -53,6 +53,9 @@ def _fmt(value, fmt: str = ".1f", fallback: str = "\u2014") -> str:
 
 # -- Plot generators --------------------------------------------------------
 
+_COMPARISON_LINESTYLES = ["-", "--", "-.", ":", (0, (3, 1, 1, 1))]
+
+
 def _plot_coupled_spl_comparison(
     csv_pairs: List[Tuple[str, str]],
     target: TargetSpec,
@@ -61,13 +64,12 @@ def _plot_coupled_spl_comparison(
     fig, ax = plot_theme.create_figure(figsize=(11, 5.5))
 
     all_freq = []
-    for csv_path, label in csv_pairs:
+    for i, (csv_path, label) in enumerate(csv_pairs):
         df = pd.read_csv(csv_path)
-        # Extract profile name from label (last word in parens)
-        profile = label.rsplit("(", 1)[-1].rstrip(")") if "(" in label else ""
-        style = plot_theme.profile_style(profile)
+        color = plot_theme.MULTI_COLORS[i % len(plot_theme.MULTI_COLORS)]
+        ls = _COMPARISON_LINESTYLES[i % len(_COMPARISON_LINESTYLES)]
         ax.plot(df["frequency"], df["spl"], label=label,
-                color=style["color"], linestyle=style["linestyle"], linewidth=1.4)
+                color=color, linestyle=ls, linewidth=1.8)
         all_freq.extend(df["frequency"].values)
 
     plot_theme.target_band_span(ax, target)
@@ -82,6 +84,37 @@ def _plot_coupled_spl_comparison(
     return plot_theme.fig_to_b64(fig)
 
 
+_SHORT_TO_PROFILE = {
+    "con": "conical",
+    "exp": "exponential",
+    "hyp": "hyperbolic",
+    "tra": "tractrix",
+    "os": "os",
+    "lec": "lecleach",
+    "cd": "cd",
+}
+
+
+def _pick_representative_per_profile(solver_csvs: Dict[str, str]) -> Dict[str, str]:
+    """Select one representative candidate per horn profile.
+
+    Groups candidates by profile (extracted from IDs like ``auto_hyp_0012``)
+    and picks the median candidate from each group.
+    """
+    groups: Dict[str, list] = {}
+    for cid, path in solver_csvs.items():
+        parts = cid.split("_")
+        short = parts[1] if len(parts) >= 2 else cid
+        profile = _SHORT_TO_PROFILE.get(short, short)
+        groups.setdefault(profile, []).append((cid, path))
+
+    result: Dict[str, str] = {}
+    for profile, items in sorted(groups.items()):
+        items.sort()
+        result[profile] = items[len(items) // 2][1]
+    return result
+
+
 def _plot_raw_profile_spl(
     solver_csvs: Dict[str, str],
     target: TargetSpec,
@@ -302,9 +335,10 @@ def _render_drivers_rows(drivers: Dict[str, DriverParameters]) -> str:
   .card .label {{ font-size: 0.75em; text-transform: uppercase; color: #64748b; letter-spacing: 0.05em; }}
   .card .value {{ font-size: 1.5em; font-weight: 700; color: #0f172a; margin-top: 4px; }}
   /* Tables */
-  table {{ width: 100%; border-collapse: collapse; font-size: 0.85em; background: #fff; border-radius: 8px; overflow: hidden; }}
+  .table-wrap {{ overflow-x: auto; -webkit-overflow-scrolling: touch; margin: 0 -4px; padding: 0 4px; }}
+  table {{ width: 100%; border-collapse: collapse; font-size: 0.85em; background: #fff; border-radius: 8px; overflow: hidden; min-width: 800px; }}
   th {{ background: #f1f5f9; text-align: left; padding: 10px 12px; font-weight: 600; white-space: nowrap; }}
-  td {{ padding: 8px 12px; border-top: 1px solid #e2e8f0; }}
+  td {{ padding: 8px 12px; border-top: 1px solid #e2e8f0; white-space: nowrap; }}
   tr:hover {{ background: #f8fafc; }}
   /* Plots */
   .plot {{ text-align: center; margin: 16px 0; }}
@@ -344,6 +378,7 @@ def _render_drivers_rows(drivers: Dict[str, DriverParameters]) -> str:
 {geometry_section}
 
 <h2>Rankings</h2>
+<div class="table-wrap">
 <table>
 <thead>
 <tr>
@@ -357,6 +392,7 @@ def _render_drivers_rows(drivers: Dict[str, DriverParameters]) -> str:
 {rankings_rows}
 </tbody>
 </table>
+</div>
 
 <h2>Coupled SPL — Top Candidates</h2>
 <div class="plot"><img src="{plot_coupled_spl}" alt="Coupled SPL comparison"></div>
@@ -371,6 +407,7 @@ def _render_drivers_rows(drivers: Dict[str, DriverParameters]) -> str:
 </div>
 
 <h2>Pre-Screened Drivers</h2>
+<div class="table-wrap">
 <table>
 <thead>
 <tr>
@@ -382,6 +419,7 @@ def _render_drivers_rows(drivers: Dict[str, DriverParameters]) -> str:
 {drivers_rows}
 </tbody>
 </table>
+</div>
 
 <div class="footer">Horn Auto-Select Report &mdash; generated {timestamp}</div>
 
@@ -394,23 +432,46 @@ def _render_drivers_rows(drivers: Dict[str, DriverParameters]) -> str:
 # -- Public entry point -----------------------------------------------------
 
 def _render_design_summary(derived_geometry: dict) -> str:
-    """Render the Design Summary section for fullauto mode."""
+    """Render the Design Summary section with optimization parameters."""
     mr = derived_geometry.get("mouth_radius_range", [])
     lr = derived_geometry.get("length_range", [])
     sr = derived_geometry.get("sim_freq_range", [])
+
+    # Determine fixed vs optimized for mouth radius and length
+    mr_fixed = len(mr) == 2 and abs(mr[1] - mr[0]) < 1e-6
+    lr_fixed = len(lr) == 2 and abs(lr[1] - lr[0]) < 1e-6
+
+    if mr_fixed:
+        mouth_desc = f'{_fmt(mr[0], ".3f")} m <span style="color:#059669;font-weight:600">(fixed)</span>'
+    else:
+        mouth_desc = (
+            f'{_fmt(mr[0] if mr else None, ".4f")} — '
+            f'{_fmt(mr[1] if len(mr) > 1 else None, ".4f")} m '
+            f'<span style="color:#2563eb;font-weight:600">(optimized)</span>'
+        )
+
+    if lr_fixed:
+        length_desc = f'{_fmt(lr[0], ".3f")} m <span style="color:#059669;font-weight:600">(fixed)</span>'
+    else:
+        length_desc = (
+            f'{_fmt(lr[0] if lr else None, ".4f")} — '
+            f'{_fmt(lr[1] if len(lr) > 1 else None, ".4f")} m '
+            f'<span style="color:#2563eb;font-weight:600">(optimized, &lambda;/4 to &lambda;/2)</span>'
+        )
+
     return (
         '<h2>Design Summary</h2>\n'
         '<div class="design-summary"><dl>'
+        f'<dt>Target frequency band</dt><dd>{_fmt(derived_geometry.get("target_f_low"), ".0f")} — '
+        f'{_fmt(derived_geometry.get("target_f_high"), ".0f")} Hz</dd>'
+        f'<dt>Mouth radius</dt><dd>{mouth_desc}</dd>'
+        f'<dt>Length</dt><dd>{length_desc}</dd>'
         f'<dt>Ideal mouth radius</dt><dd>{_fmt(derived_geometry.get("ideal_mouth_radius"), ".4f")} m '
         f'(from c&#8320;/2&pi;f<sub>low</sub>)</dd>'
-        f'<dt>Mouth radius range</dt><dd>{_fmt(mr[0] if mr else None, ".4f")} — '
-        f'{_fmt(mr[1] if len(mr) > 1 else None, ".4f")} m (&plusmn;30%)</dd>'
-        f'<dt>Length range</dt><dd>{_fmt(lr[0] if lr else None, ".4f")} — '
-        f'{_fmt(lr[1] if len(lr) > 1 else None, ".4f")} m '
-        f'(&lambda;/4 to &lambda;/2)</dd>'
         f'<dt>Simulation freq range</dt><dd>{_fmt(sr[0] if sr else None, ".0f")} — '
         f'{_fmt(sr[1] if len(sr) > 1 else None, ".0f")} Hz (&plusmn;0.5 octave)</dd>'
-        f'<dt>Candidate count</dt><dd>{derived_geometry.get("candidate_count", "—")}</dd>'
+        f'<dt>Profiles</dt><dd>7 (conical, exponential, hyperbolic, tractrix, OS, Le Cl&eacute;ac&#8217;h, CD)</dd>'
+        f'<dt>Geometry candidates</dt><dd>{derived_geometry.get("candidate_count", "—")}</dd>'
         '</dl></div>'
     )
 
@@ -426,6 +487,9 @@ def generate_html_report(
     mouth_radius: Optional[float] = None,
     length: Optional[float] = None,
     derived_geometry: Optional[dict] = None,
+    total_candidates: Optional[int] = None,
+    total_scored: Optional[int] = None,
+    lem_results: Optional[dict] = None,
 ) -> str:
     """Generate a self-contained HTML report string.
 
@@ -440,6 +504,8 @@ def generate_html_report(
         mouth_radius: Horn mouth radius in metres (enables 3D geometry renders).
         length: Horn length in metres (enables 3D geometry renders).
         derived_geometry: Optional dict from geometry_designer (fullauto mode).
+        total_candidates: Total geometry candidates simulated.
+        total_scored: Total driver-horn combinations scored (before top-N).
 
     Returns:
         Complete HTML document as a string.
@@ -459,9 +525,10 @@ def generate_html_report(
         min((r.get("passband_ripple_db", 99) for r in top_results), default=None), ".1f"
     ) if top_results else "\u2014"
 
-    # Generate 3D horn geometry renders (when dimensions are provided)
+    # Generate 3D horn geometry renders from top-ranked candidates
     geometry_html = ""
     if mouth_radius is not None and length is not None:
+        # Fixed geometry: render each profile
         geom_imgs = []
         for profile in sorted(solver_csvs.keys()):
             b64 = fig_to_b64_3d(
@@ -476,12 +543,46 @@ def generate_html_report(
                 f'alt="{html.escape(profile.capitalize())} horn geometry"></div>'
             )
         geometry_html = "\n".join(geom_imgs)
+    elif top_results:
+        # Variable geometry (auto/fullauto): render unique geometries from top results
+        seen = set()
+        geom_imgs = []
+        for r in top_results:
+            profile = r.get("profile", "") or r.get("horn_label", "")
+            horn_label = r.get("horn_label", "")
+            mr = r.get("mouth_radius")
+            ln = r.get("length")
+            tr = r.get("throat_radius", throat_radius)
+            if profile and mr and ln and mr > tr:
+                key = (profile, round(tr, 6), round(mr, 6), round(ln, 6))
+                if key not in seen:
+                    seen.add(key)
+                    title = (
+                        f"{profile.capitalize()} ({horn_label})"
+                        f" — throat {tr*1000:.1f}, mouth {mr*1000:.1f}, L {ln*1000:.0f} mm"
+                    )
+                    b64 = fig_to_b64_3d(
+                        throat_radius=tr,
+                        mouth_radius=mr,
+                        length=ln,
+                        profile=profile,
+                        title=title,
+                        figsize=(12, 5),
+                    )
+                    geom_imgs.append(
+                        f'<div class="plot"><img src="{b64}" '
+                        f'alt="{html.escape(title)}"></div>'
+                    )
+        geometry_html = "\n".join(geom_imgs)
+
+    # Reduce to one representative per profile for raw/impedance/phase plots
+    representative_csvs = _pick_representative_per_profile(solver_csvs)
 
     # Generate plots
     plot_coupled_spl = _plot_coupled_spl_comparison(csv_pairs, target)
-    plot_raw_spl = _plot_raw_profile_spl(solver_csvs, target)
-    plot_impedance = _plot_profile_impedance(solver_csvs)
-    plot_phase = _plot_profile_phase(solver_csvs)
+    plot_raw_spl = _plot_raw_profile_spl(representative_csvs, target)
+    plot_impedance = _plot_profile_impedance(representative_csvs)
+    plot_phase = _plot_profile_phase(representative_csvs)
 
     # Render tables
     rankings_rows = _render_rankings_rows(top_results, drivers, show_geometry=show_geometry)
@@ -489,6 +590,20 @@ def generate_html_report(
 
     # Conditional sections for fullauto
     design_summary_section = _render_design_summary(derived_geometry) if show_geometry else ""
+
+    # LEM prescreen stats
+    if lem_results:
+        lem_total = lem_results.get("total_evaluated", "—")
+        lem_pairs = lem_results.get("total_pairs", "—")
+        lem_passed = len(lem_results.get("filtered_candidate_ids", []))
+        design_summary_section += (
+            '<h2>LEM Prescreening</h2>\n'
+            '<div class="design-summary"><dl>'
+            f'<dt>Candidates evaluated (LEM/Webster)</dt><dd>{lem_total}</dd>'
+            f'<dt>Driver-horn pairs scored</dt><dd>{lem_pairs}</dd>'
+            f'<dt>Candidates passed to FEM</dt><dd>{lem_passed}</dd>'
+            '</dl></div>'
+        )
     geometry_header_cols = '<th>Mouth R (m)</th><th>Length (m)</th>' if show_geometry else ""
 
     # Mouth/Length display: for fullauto show "varies", for auto show fixed value
@@ -512,6 +627,8 @@ def generate_html_report(
     else:
         geometry_section = ""
 
+    n_scored = total_scored if total_scored is not None else len(all_ranked)
+
     return _HTML_TEMPLATE.format_map({
         "target_low": target.f_low_hz,
         "target_high": target.f_high_hz,
@@ -520,7 +637,7 @@ def generate_html_report(
         "horn_length": length_display,
         "profiles": ", ".join(sorted(solver_csvs.keys())),
         "timestamp": timestamp,
-        "n_scored": len(all_ranked),
+        "n_scored": n_scored,
         "n_top": len(top_results),
         "best_score": best_score,
         "best_bw": best_bw,
diff --git a/packages/horn-analysis/src/horn_analysis/kpi.py b/packages/horn-analysis/src/horn_analysis/kpi.py
index f14d913..7e999b8 100644
--- a/packages/horn-analysis/src/horn_analysis/kpi.py
+++ b/packages/horn-analysis/src/horn_analysis/kpi.py
@@ -18,6 +18,7 @@
 import pandas as pd
 from scipy.interpolate import interp1d
 from scipy.optimize import brentq
+from scipy.signal import savgol_filter
 
 
 @dataclass
@@ -62,9 +63,15 @@ def spl_minus_threshold(f):
 
     # Find f3_low: search from lowest freq up to peak
     f3_low = _find_crossing(spl_minus_threshold, freq[0], freq[peak_idx], direction="rising")
+    if f3_low is None and spl_minus_threshold(freq[0]) >= 0:
+        # Response is above -3 dB at the low measurement boundary
+        f3_low = float(freq[0])
 
     # Find f3_high: search from peak to highest freq
     f3_high = _find_crossing(spl_minus_threshold, freq[peak_idx], freq[-1], direction="falling")
+    if f3_high is None and spl_minus_threshold(freq[-1]) >= 0:
+        # Response is above -3 dB at the high measurement boundary
+        f3_high = float(freq[-1])
 
     # Derived KPIs
     bandwidth_hz = None
@@ -76,12 +83,20 @@ def spl_minus_threshold(f):
         bandwidth_hz = f3_high - f3_low
         bandwidth_octaves = np.log2(f3_high / f3_low) if f3_low > 0 else None
 
-        # Passband: frequencies within [f3_low, f3_high]
-        mask = (freq >= f3_low) & (freq <= f3_high)
-        if np.any(mask):
-            passband_spl = spl[mask]
-            passband_ripple = float(np.max(passband_spl) - np.min(passband_spl))
-            avg_sensitivity = float(np.mean(passband_spl))
+        # Passband: resample onto uniform log grid and smooth to remove
+        # band-stitching artifacts before computing ripple/sensitivity
+        n_passband = max(200, len(freq) * 2)
+        freq_uniform = np.geomspace(f3_low, f3_high, n_passband)
+        spl_passband = np.array([float(spl_interp(f)) for f in freq_uniform])
+
+        # Light Savitzky-Golay smoothing to suppress band-boundary glitches
+        # Window must be odd and < n_passband; 11 points is ~5% of 200
+        sg_window = min(11, n_passband if n_passband % 2 == 1 else n_passband - 1)
+        if sg_window >= 5:
+            spl_passband = savgol_filter(spl_passband, sg_window, polyorder=3)
+
+        passband_ripple = float(np.max(spl_passband) - np.min(spl_passband))
+        avg_sensitivity = float(np.mean(spl_passband))
 
     return HornKPI(
         peak_spl_db=peak_spl,
diff --git a/packages/horn-analysis/src/horn_analysis/lem_prescreen.py b/packages/horn-analysis/src/horn_analysis/lem_prescreen.py
new file mode 100644
index 0000000..c794fcf
--- /dev/null
+++ b/packages/horn-analysis/src/horn_analysis/lem_prescreen.py
@@ -0,0 +1,235 @@
+"""LEM (Lumped Element Model) prescreening for auto pipeline.
+
+Combines Webster/TMM throat impedance with driver coupling and scoring
+to rapidly rank all geometry candidates analytically. Only the top-N
+candidates proceed to expensive FEM simulation.
+
+Lives in horn-analysis (not horn-core) because it imports transfer_function,
+kpi, and scoring from horn-analysis.
+"""
+
+import argparse
+import csv
+import json
+from pathlib import Path
+from typing import Dict, List, Optional
+
+import numpy as np
+
+from horn_core.candidates import CandidateGeometry
+from horn_core.parameters import DriverParameters
+from horn_core.profiles import get_radius_func
+from horn_core.webster import compute_throat_impedance_tmm
+from horn_analysis.kpi import extract_kpis_from_arrays
+from horn_analysis.scoring import TargetSpec, compute_selection_score
+from horn_analysis.transfer_function import compute_driver_response
+
+
+P_REF = 20e-6  # reference pressure for dB SPL (Pa)
+
+
+def _spl_from_pressure(p_throat: np.ndarray) -> np.ndarray:
+    """Convert complex throat pressure to SPL (dB re 20 µPa)."""
+    p_mag = np.abs(p_throat)
+    p_mag = np.maximum(p_mag, 1e-30)
+    return 20.0 * np.log10(p_mag / P_REF)
+
+
+def lem_prescreen_candidates(
+    candidates: List[CandidateGeometry],
+    drivers: List[DriverParameters],
+    target_f_low: float,
+    target_f_high: float,
+    sim_freq_range: tuple,
+    num_frequencies: int = 100,
+    top_n: int = 10,
+    n_segments: int = 200,
+) -> dict:
+    """Score all (candidate, driver) pairs using TMM + driver coupling.
+
+    For each pair:
+    1. get_radius_func() → compute_throat_impedance_tmm() → z_real, z_imag
+    2. compute_driver_response(driver, freq, z_real, z_imag, throat_area) → p_throat
+    3. SPL = 20·log₁₀(|p_throat| / 20µPa)
+    4. extract_kpis_from_arrays() → compute_selection_score() → composite score
+
+    Ranks all scores, collects unique top-N candidate_ids.
+
+    Args:
+        candidates: List of CandidateGeometry from the candidates CSV.
+        drivers: Pre-screened drivers.
+        target_f_low: Target low frequency (Hz).
+        target_f_high: Target high frequency (Hz).
+        sim_freq_range: (min_freq, max_freq) for simulation.
+        num_frequencies: Number of frequency points for TMM evaluation.
+        top_n: Number of top candidates to pass to FEM.
+        n_segments: Number of TMM segments per horn.
+
+    Returns:
+        Dict with keys: total_evaluated, total_pairs, top_n, rankings (list),
+        filtered_candidate_ids (list of unique candidate IDs).
+    """
+    frequencies = np.linspace(sim_freq_range[0], sim_freq_range[1], num_frequencies)
+    target = TargetSpec(f_low_hz=target_f_low, f_high_hz=target_f_high)
+
+    all_scores = []
+
+    # Cache TMM results per candidate (independent of driver)
+    tmm_cache: Dict[str, tuple] = {}
+
+    for cand in candidates:
+        cache_key = cand.candidate_id
+        if cache_key not in tmm_cache:
+            radius_func = get_radius_func(
+                cand.profile, cand.throat_radius, cand.mouth_radius, cand.length
+            )
+            z_real, z_imag = compute_throat_impedance_tmm(
+                frequencies=frequencies,
+                radius_func=radius_func,
+                length=cand.length,
+                throat_radius=cand.throat_radius,
+                mouth_radius=cand.mouth_radius,
+                n_segments=n_segments,
+            )
+            tmm_cache[cache_key] = (z_real, z_imag)
+
+        z_real, z_imag = tmm_cache[cache_key]
+        throat_area = np.pi * cand.throat_radius**2
+
+        for drv in drivers:
+            p_throat = compute_driver_response(
+                drv, frequencies, z_real, z_imag, throat_area
+            )
+            coupled_spl = _spl_from_pressure(p_throat)
+            kpi = extract_kpis_from_arrays(frequencies, coupled_spl)
+            score = compute_selection_score(
+                kpi, target,
+                driver_id=drv.driver_id,
+                horn_label=cand.candidate_id,
+            )
+
+            all_scores.append({
+                "candidate_id": cand.candidate_id,
+                "profile": cand.profile,
+                "throat_radius": cand.throat_radius,
+                "mouth_radius": cand.mouth_radius,
+                "length": cand.length,
+                "driver_id": drv.driver_id,
+                "composite_score": score.composite_score,
+                "bandwidth_coverage": score.bandwidth_coverage,
+                "passband_ripple_db": score.passband_ripple_db,
+                "avg_sensitivity_db": score.avg_sensitivity_db,
+            })
+
+    # Sort by composite score descending
+    all_scores.sort(key=lambda x: x["composite_score"], reverse=True)
+
+    # Collect unique top-N candidate IDs (by best score per candidate)
+    seen_ids = set()
+    filtered_ids = []
+    for entry in all_scores:
+        cid = entry["candidate_id"]
+        if cid not in seen_ids:
+            seen_ids.add(cid)
+            filtered_ids.append(cid)
+            if len(filtered_ids) >= top_n:
+                break
+
+    return {
+        "total_evaluated": len(candidates),
+        "total_pairs": len(all_scores),
+        "top_n": top_n,
+        "filtered_candidate_ids": filtered_ids,
+        "rankings": all_scores,
+    }
+
+
+def _load_candidates_csv(csv_path: str) -> List[CandidateGeometry]:
+    """Load candidates from a CSV file."""
+    candidates = []
+    with open(csv_path) as f:
+        for row in csv.DictReader(f):
+            candidates.append(CandidateGeometry(
+                candidate_id=row["candidate_id"],
+                profile=row["profile"],
+                throat_radius=float(row["throat_radius"]),
+                mouth_radius=float(row["mouth_radius"]),
+                length=float(row["length"]),
+            ))
+    return candidates
+
+
+def _write_filtered_csv(
+    candidates: List[CandidateGeometry],
+    filtered_ids: List[str],
+    output_path: str,
+) -> None:
+    """Write filtered candidates CSV containing only the top-N."""
+    id_set = set(filtered_ids)
+    filtered = [c for c in candidates if c.candidate_id in id_set]
+
+    with open(output_path, "w", newline="") as f:
+        writer = csv.writer(f)
+        writer.writerow(["candidate_id", "profile", "throat_radius", "mouth_radius", "length"])
+        for c in filtered:
+            writer.writerow([c.candidate_id, c.profile, c.throat_radius, c.mouth_radius, c.length])
+
+
+def main():
+    """CLI entry point for LEM prescreening."""
+    parser = argparse.ArgumentParser(
+        description="LEM/Webster prescreening for auto pipeline candidates.",
+    )
+    parser.add_argument("--candidates-csv", required=True, help="Input candidates CSV.")
+    parser.add_argument("--prescreen-json", required=True, help="Driver prescreen result JSON.")
+    parser.add_argument("--drivers-db", required=True, help="Driver database path.")
+    parser.add_argument("--design-json", required=True, help="Design JSON with sim_freq_range.")
+    parser.add_argument("--target-f-low", type=float, required=True, help="Target low freq (Hz).")
+    parser.add_argument("--target-f-high", type=float, required=True, help="Target high freq (Hz).")
+    parser.add_argument("--top-n", type=int, default=10, help="Number of top candidates for FEM.")
+    parser.add_argument("--num-frequencies", type=int, default=100, help="Frequency points for TMM.")
+    parser.add_argument("--output", required=True, help="Output LEM results JSON.")
+    parser.add_argument("--filtered-csv", required=True, help="Output filtered candidates CSV.")
+    args = parser.parse_args()
+
+    from horn_drivers.loader import load_drivers
+
+    # Load candidates
+    candidates = _load_candidates_csv(args.candidates_csv)
+    print(f"Loaded {len(candidates)} candidates from {args.candidates_csv}")
+
+    # Load pre-screened driver IDs and full driver data
+    prescreen = json.loads(Path(args.prescreen_json).read_text())
+    driver_ids = set(prescreen["drivers"])
+    all_drivers = load_drivers(args.drivers_db)
+    drivers = [d for d in all_drivers if d.driver_id in driver_ids]
+    print(f"Using {len(drivers)} pre-screened drivers")
+
+    # Load simulation frequency range from design JSON
+    design = json.loads(Path(args.design_json).read_text())
+    sim_freq_range = tuple(design["sim_freq_range"])
+
+    # Run LEM prescreening
+    results = lem_prescreen_candidates(
+        candidates=candidates,
+        drivers=drivers,
+        target_f_low=args.target_f_low,
+        target_f_high=args.target_f_high,
+        sim_freq_range=sim_freq_range,
+        num_frequencies=args.num_frequencies,
+        top_n=args.top_n,
+    )
+
+    # Write outputs
+    Path(args.output).write_text(json.dumps(results, indent=2))
+    _write_filtered_csv(candidates, results["filtered_candidate_ids"], args.filtered_csv)
+
+    print(f"LEM prescreening complete:")
+    print(f"  {results['total_evaluated']} candidates × {len(drivers)} drivers = {results['total_pairs']} pairs")
+    print(f"  Top {len(results['filtered_candidate_ids'])} candidates passed to FEM")
+    print(f"  Results: {args.output}")
+    print(f"  Filtered CSV: {args.filtered_csv}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/packages/horn-analysis/src/horn_analysis/plot_theme.py b/packages/horn-analysis/src/horn_analysis/plot_theme.py
index c8c77cf..768edcf 100644
--- a/packages/horn-analysis/src/horn_analysis/plot_theme.py
+++ b/packages/horn-analysis/src/horn_analysis/plot_theme.py
@@ -28,14 +28,14 @@
 }
 
 MULTI_COLORS = [
-    "#1f4e79",
-    "#c0392b",
-    "#27864e",
-    "#7b4ea3",
-    "#d4831a",
-    "#2a9d8f",
-    "#6c5b7b",
-    "#c97b3d",
+    "#2563eb",  # vivid blue
+    "#dc2626",  # vivid red
+    "#16a34a",  # vivid green
+    "#9333ea",  # vivid purple
+    "#ea580c",  # vivid orange
+    "#0891b2",  # cyan
+    "#c026d3",  # magenta
+    "#ca8a04",  # amber
 ]
 
 PROFILE_STYLES = {
diff --git a/packages/horn-analysis/src/horn_analysis/prescreen.py b/packages/horn-analysis/src/horn_analysis/prescreen.py
index 7147f7c..c306e01 100644
--- a/packages/horn-analysis/src/horn_analysis/prescreen.py
+++ b/packages/horn-analysis/src/horn_analysis/prescreen.py
@@ -7,7 +7,8 @@
 import argparse
 import json
 import math
-from dataclasses import dataclass, asdict
+import re
+from dataclasses import dataclass, field
 from pathlib import Path
 from typing import List, Tuple, Optional
 
@@ -26,8 +27,9 @@ class PrescreenConfig:
     mouth_radius_m: Optional[float] = None
     length_m: Optional[float] = None
     min_ebp: float = 50.0
-    horn_load_factor: float = 10.0
-    sd_ratio_range: Tuple[float, float] = (0.3, 3.0)
+    ka_max: float = 2 * math.pi  # ~6.28, absolute cap on throat ka at f_high
+    min_nominal_diameter_in: Optional[float] = None
+    max_nominal_diameter_in: Optional[float] = None
 
 
 @dataclass
@@ -35,28 +37,43 @@ class PrescreenResult:
     """Result of driver pre-screening."""
     drivers: List[DriverParameters]
     throat_radius_m: float
+    throat_radii_m: List[float]
     count: int
 
     def to_dict(self) -> dict:
         return {
             "drivers": [d.driver_id for d in self.drivers],
             "throat_radius_m": self.throat_radius_m,
+            "throat_radii_m": self.throat_radii_m,
             "count": self.count,
         }
 
 
+def _parse_diameter_inches(d: Optional[str]) -> Optional[float]:
+    """Parse a nominal diameter string like '4in' or '6.5' to inches."""
+    if not d:
+        return None
+    m = re.match(r"(\d+(?:\.\d+)?)", d)
+    return float(m.group(1)) if m else None
+
+
 def prescreen_drivers(
     drivers: List[DriverParameters],
     config: PrescreenConfig,
 ) -> PrescreenResult:
     """Filter drivers to candidates suitable for the target horn.
 
+    Throat radius is decoupled from driver size — the throat is a property
+    of the horn (constrained by acoustics at f_high), not the driver.
+
     Filtering criteria:
     1. fs_hz < target_f_low_hz * 1.5 -- driver resonance in or near target band
     2. fs_hz / qes > min_ebp -- horn suitability (Efficiency Bandwidth Product)
-    3. Upper freq capability: f_piston * horn_load_factor >= target_f_high
-    4. Representative throat radius = median(sqrt(Sd/pi)) of passing drivers
-    5. Filter drivers whose effective radius is outside sd_ratio_range of representative
+    3. Optional nominal diameter filter (user-specified min/max)
+    4. Driver must physically fit in the horn mouth
+
+    Throat radii are derived from acoustic constraints (ka ≤ ka_max at f_high),
+    not from driver Sd.
 
     Args:
         drivers: Full list of drivers from the database.
@@ -67,6 +84,13 @@ def prescreen_drivers(
     """
     candidates = []
 
+    # Max driver radius: driver must fit inside the horn mouth
+    if config.mouth_radius_m is not None and config.mouth_radius_m > 0:
+        max_driver_radius = config.mouth_radius_m
+    else:
+        ideal_mouth = _SPEED_OF_SOUND / (2 * math.pi * config.target_f_low_hz)
+        max_driver_radius = ideal_mouth
+
     for drv in drivers:
         # 1. Resonance frequency check
         if drv.fs_hz >= config.target_f_low_hz * 1.5:
@@ -78,50 +102,59 @@ def prescreen_drivers(
             if ebp < config.min_ebp:
                 continue
         else:
-            # Cannot compute EBP, skip
             continue
 
-        # 3. Upper frequency capability from Sd
-        #    f_piston = c / (2π·a_eff) where a_eff = √(Sd/π)
-        #    Horn loading extends usable range by ~10× above piston breakup
-        if drv.sd_m2 > 0:
-            a_eff = math.sqrt(drv.sd_m2 / math.pi)
-            f_piston = _SPEED_OF_SOUND / (2 * math.pi * a_eff)
-            if f_piston * config.horn_load_factor < config.target_f_high_hz:
-                continue
+        # 3. Optional nominal diameter filter
+        if config.min_nominal_diameter_in is not None or config.max_nominal_diameter_in is not None:
+            dia_in = _parse_diameter_inches(drv.nominal_diameter)
+            if dia_in is not None:
+                if config.min_nominal_diameter_in is not None and dia_in < config.min_nominal_diameter_in:
+                    continue
+                if config.max_nominal_diameter_in is not None and dia_in > config.max_nominal_diameter_in:
+                    continue
+
+        # 4. Driver must fit in the horn mouth
+        drv_radius = math.sqrt(drv.effective_throat_area / math.pi)
+        if drv_radius > max_driver_radius:
+            continue
 
         candidates.append(drv)
 
     if not candidates:
-        return PrescreenResult(drivers=[], throat_radius_m=0.0, count=0)
+        return PrescreenResult(drivers=[], throat_radius_m=0.0, throat_radii_m=[], count=0)
 
-    # 4. Compute representative throat radius
-    radii = [math.sqrt(d.sd_m2 / math.pi) for d in candidates if d.sd_m2 > 0]
-    if not radii:
-        return PrescreenResult(drivers=[], throat_radius_m=0.0, count=0)
+    # Derive throat radii from acoustic constraints (ka ≤ ka_max at f_high)
+    a_acoustic_max = _SPEED_OF_SOUND * config.ka_max / (
+        2 * math.pi * config.target_f_high_hz
+    )
 
-    representative_radius = float(np.median(radii))
+    # Acoustic range: fractions of the maximum acoustic throat radius
+    acoustic_radii = [a_acoustic_max * f for f in [0.3, 0.65, 1.0]]
 
-    # 5. Filter by Sd ratio relative to representative
-    lo, hi = config.sd_ratio_range
-    filtered = []
-    for drv in candidates:
-        if drv.sd_m2 <= 0:
-            continue
-        drv_radius = math.sqrt(drv.sd_m2 / math.pi)
-        ratio = drv_radius / representative_radius
-        if lo <= ratio <= hi:
-            filtered.append(drv)
+    # Also include driver-matched radii for direct-coupling scenarios
+    driver_radii = [
+        math.sqrt(d.effective_throat_area / math.pi)
+        for d in candidates if d.effective_throat_area > 0
+    ]
+    direct_radii = [r for r in driver_radii if r <= a_acoustic_max]
+
+    all_radii = sorted(set(
+        round(r, 6) for r in acoustic_radii + direct_radii
+    ))
+
+    # Keep at most 5 to limit combinatorial explosion
+    if len(all_radii) > 5:
+        indices = np.linspace(0, len(all_radii) - 1, 5, dtype=int)
+        all_radii = [all_radii[i] for i in indices]
 
-    # Recompute representative from final set
-    if filtered:
-        final_radii = [math.sqrt(d.sd_m2 / math.pi) for d in filtered]
-        representative_radius = float(np.median(final_radii))
+    representative_radius = all_radii[len(all_radii) // 2]
+    throat_radii_m = all_radii
 
     return PrescreenResult(
-        drivers=filtered,
+        drivers=candidates,
         throat_radius_m=representative_radius,
-        count=len(filtered),
+        throat_radii_m=throat_radii_m,
+        count=len(candidates),
     )
 
 
@@ -136,8 +169,12 @@ def main():
     parser.add_argument("--mouth-radius", type=float, default=None, help="Horn mouth radius (m). Optional for fullauto mode.")
     parser.add_argument("--length", type=float, default=None, help="Horn length (m). Optional for fullauto mode.")
     parser.add_argument("--min-ebp", type=float, default=50.0, help="Minimum EBP threshold.")
-    parser.add_argument("--horn-load-factor", type=float, default=10.0,
-                        help="Multiplier on piston frequency for horn-loaded upper limit estimate.")
+    parser.add_argument("--ka-max", type=float, default=2 * math.pi,
+                        help="Absolute cap on throat ka at f_high (default 2π ≈ 6.28).")
+    parser.add_argument("--min-diameter", type=float, default=None,
+                        help="Minimum driver nominal diameter (inches).")
+    parser.add_argument("--max-diameter", type=float, default=None,
+                        help="Maximum driver nominal diameter (inches).")
     parser.add_argument("--output", type=str, default="prescreen_result.json", help="Output JSON file.")
     args = parser.parse_args()
 
@@ -152,7 +189,9 @@ def main():
         mouth_radius_m=args.mouth_radius,
         length_m=args.length,
         min_ebp=args.min_ebp,
-        horn_load_factor=args.horn_load_factor,
+        ka_max=args.ka_max,
+        min_nominal_diameter_in=args.min_diameter,
+        max_nominal_diameter_in=args.max_diameter,
     )
 
     result = prescreen_drivers(drivers, config)
diff --git a/packages/horn-analysis/src/horn_analysis/scoring.py b/packages/horn-analysis/src/horn_analysis/scoring.py
index a9a0e49..62c5ee4 100644
--- a/packages/horn-analysis/src/horn_analysis/scoring.py
+++ b/packages/horn-analysis/src/horn_analysis/scoring.py
@@ -90,6 +90,10 @@ def compute_selection_score(
         + weights["sensitivity"] * float(sensitivity_score)
     )
 
+    # --- Bandwidth floor: non-functional combos score zero ---
+    if bandwidth_coverage < 0.10:
+        composite = 0.0
+
     return SelectionScore(
         driver_id=driver_id,
         horn_label=horn_label,
diff --git a/packages/horn-analysis/tests/test_analysis.py b/packages/horn-analysis/tests/test_analysis.py
index bd19ec3..0653738 100644
--- a/packages/horn-analysis/tests/test_analysis.py
+++ b/packages/horn-analysis/tests/test_analysis.py
@@ -118,8 +118,8 @@ def test_to_dict_and_json(self, bandpass_csv):
         json_str = json.dumps(d)
         assert len(json_str) > 0
 
-    def test_flat_response_no_f3(self, tmp_path):
-        """A flat response at the peak should have no -3dB crossing below the peak."""
+    def test_flat_response_covers_full_range(self, tmp_path):
+        """A flat response covers the full measurement range (never drops 3 dB)."""
         from horn_analysis.kpi import extract_kpis
 
         csv_path = tmp_path / "flat.csv"
@@ -129,9 +129,9 @@ def test_flat_response_no_f3(self, tmp_path):
 
         kpis = extract_kpis(str(csv_path))
         assert kpis.peak_spl_db == pytest.approx(90.0)
-        # With perfectly flat response, there's no -3dB crossing
-        assert kpis.f3_low_hz is None
-        assert kpis.f3_high_hz is None
+        # Flat response is above -3 dB at both edges → full range
+        assert kpis.f3_low_hz == pytest.approx(100.0)
+        assert kpis.f3_high_hz == pytest.approx(10000.0)
 
 
 class TestMultiComparison:
diff --git a/packages/horn-analysis/tests/test_kpi_arrays.py b/packages/horn-analysis/tests/test_kpi_arrays.py
index f7bd7ec..d455ef7 100644
--- a/packages/horn-analysis/tests/test_kpi_arrays.py
+++ b/packages/horn-analysis/tests/test_kpi_arrays.py
@@ -33,7 +33,11 @@ def test_returns_horn_kpi(self, bandpass_data):
         assert isinstance(result, HornKPI)
 
     def test_matches_csv_version(self, bandpass_data, bandpass_csv):
-        """Array version should produce identical results to CSV version."""
+        """Array version should produce similar results to CSV version.
+
+        Note: passband_ripple and avg_sensitivity may differ slightly due to
+        resampling onto a uniform log grid + Savgol smoothing.
+        """
         freq, spl = bandpass_data
         from_arrays = extract_kpis_from_arrays(freq, spl)
         from_csv = extract_kpis(bandpass_csv)
@@ -57,11 +61,12 @@ def test_matches_csv_version(self, bandpass_data, bandpass_csv):
         if from_csv.bandwidth_octaves is not None:
             assert from_arrays.bandwidth_octaves == pytest.approx(from_csv.bandwidth_octaves, rel=1e-6)
 
+        # Ripple and sensitivity use resampled+smoothed data, so allow wider tolerance
         if from_csv.passband_ripple_db is not None:
-            assert from_arrays.passband_ripple_db == pytest.approx(from_csv.passband_ripple_db, rel=1e-6)
+            assert from_arrays.passband_ripple_db == pytest.approx(from_csv.passband_ripple_db, abs=0.5)
 
         if from_csv.average_sensitivity_db is not None:
-            assert from_arrays.average_sensitivity_db == pytest.approx(from_csv.average_sensitivity_db, rel=1e-6)
+            assert from_arrays.average_sensitivity_db == pytest.approx(from_csv.average_sensitivity_db, abs=0.3)
 
     def test_peak_detection(self, bandpass_data):
         """Should detect the peak correctly."""
@@ -71,10 +76,45 @@ def test_peak_detection(self, bandpass_data):
         assert result.peak_frequency_hz > 0
 
     def test_flat_response(self):
-        """Flat response should have no -3dB crossings."""
+        """Flat response covers the full measurement range (never drops 3 dB)."""
         freq = np.geomspace(100, 10000, 50)
         spl = np.full_like(freq, 90.0)
         result = extract_kpis_from_arrays(freq, spl)
         assert result.peak_spl_db == pytest.approx(90.0)
-        assert result.f3_low_hz is None
-        assert result.f3_high_hz is None
+        assert result.f3_low_hz == pytest.approx(100.0)
+        assert result.f3_high_hz == pytest.approx(10000.0)
+
+    def test_band_boundary_artifacts_suppressed(self):
+        """Band-stitching glitches should not inflate ripple measurement.
+
+        Simulates 8-band FEM output with 1 dB discontinuities at each
+        band boundary. The true response is flat at 90 dB, so the
+        measured ripple should be well below the raw artifact magnitude.
+        """
+        num_bands = 8
+        points_per_band = 13
+        f_min, f_max = 500.0, 8000.0
+        bands = np.geomspace(f_min, f_max, num_bands + 1)
+
+        freq_all = []
+        spl_all = []
+        for i in range(num_bands):
+            band_freq = np.geomspace(bands[i], bands[i + 1], points_per_band)
+            band_spl = np.full_like(band_freq, 90.0)
+            # Inject a 1 dB glitch at the start of each band (except first)
+            if i > 0:
+                band_spl[0] -= 1.0
+            freq_all.append(band_freq)
+            spl_all.append(band_spl)
+
+        freq = np.concatenate(freq_all)
+        spl = np.concatenate(spl_all)
+        # Sort (like the merge process does)
+        order = np.argsort(freq)
+        freq, spl = freq[order], spl[order]
+
+        result = extract_kpis_from_arrays(freq, spl)
+        # Without smoothing, ripple would be ~1.0 dB (the injected glitch).
+        # With smoothing, it should be well under 0.5 dB.
+        assert result.passband_ripple_db < 0.5
+        assert result.average_sensitivity_db == pytest.approx(90.0, abs=0.2)
diff --git a/packages/horn-analysis/tests/test_lem_prescreen.py b/packages/horn-analysis/tests/test_lem_prescreen.py
new file mode 100644
index 0000000..3f47e86
--- /dev/null
+++ b/packages/horn-analysis/tests/test_lem_prescreen.py
@@ -0,0 +1,245 @@
+"""Tests for LEM/Webster prescreening orchestrator."""
+
+import json
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+from horn_core.candidates import CandidateGeometry
+from horn_core.parameters import DriverParameters
+from horn_analysis.lem_prescreen import (
+    lem_prescreen_candidates,
+    _spl_from_pressure,
+    _load_candidates_csv,
+    _write_filtered_csv,
+)
+
+
+def _make_driver(
+    driver_id="drv1",
+    fs_hz=500.0,
+    qes=0.4,
+    qms=5.0,
+    sd_m2=0.0008,
+    re_ohm=6.0,
+    bl_tm=8.0,
+    mms_kg=0.003,
+    le_h=0.0005,
+):
+    return DriverParameters(
+        driver_id=driver_id,
+        manufacturer="Test",
+        model_name=driver_id,
+        fs_hz=fs_hz,
+        re_ohm=re_ohm,
+        bl_tm=bl_tm,
+        sd_m2=sd_m2,
+        mms_kg=mms_kg,
+        le_h=le_h,
+        qms=qms,
+        qes=qes,
+    )
+
+
+def _make_candidate(cid, profile="conical", throat=0.025, mouth=0.15, length=0.3):
+    return CandidateGeometry(
+        candidate_id=cid,
+        profile=profile,
+        throat_radius=throat,
+        mouth_radius=mouth,
+        length=length,
+    )
+
+
+class TestSplFromPressure:
+    def test_reference_pressure(self):
+        """20 µPa should give 0 dB SPL."""
+        p = np.array([20e-6])
+        spl = _spl_from_pressure(p)
+        assert spl[0] == pytest.approx(0.0, abs=0.01)
+
+    def test_1_pa(self):
+        """1 Pa should give 94 dB SPL."""
+        p = np.array([1.0])
+        spl = _spl_from_pressure(p)
+        assert spl[0] == pytest.approx(94.0, abs=0.1)
+
+
+class TestLemPrescreenCandidates:
+    @pytest.fixture
+    def driver(self):
+        return _make_driver()
+
+    @pytest.fixture
+    def candidates(self):
+        return [
+            _make_candidate("c1", "conical", 0.025, 0.15, 0.3),
+            _make_candidate("c2", "exponential", 0.025, 0.15, 0.3),
+            _make_candidate("c3", "hyperbolic", 0.025, 0.15, 0.3),
+        ]
+
+    def test_output_keys(self, candidates, driver):
+        result = lem_prescreen_candidates(
+            candidates=candidates,
+            drivers=[driver],
+            target_f_low=500,
+            target_f_high=4000,
+            sim_freq_range=(354, 5657),
+            num_frequencies=30,
+            top_n=2,
+        )
+        assert "total_evaluated" in result
+        assert "total_pairs" in result
+        assert "top_n" in result
+        assert "filtered_candidate_ids" in result
+        assert "rankings" in result
+
+    def test_top_n_respected(self, candidates, driver):
+        result = lem_prescreen_candidates(
+            candidates=candidates,
+            drivers=[driver],
+            target_f_low=500,
+            target_f_high=4000,
+            sim_freq_range=(354, 5657),
+            num_frequencies=30,
+            top_n=2,
+        )
+        assert len(result["filtered_candidate_ids"]) <= 2
+
+    def test_total_evaluated_matches_input(self, candidates, driver):
+        result = lem_prescreen_candidates(
+            candidates=candidates,
+            drivers=[driver],
+            target_f_low=500,
+            target_f_high=4000,
+            sim_freq_range=(354, 5657),
+            num_frequencies=30,
+            top_n=10,
+        )
+        assert result["total_evaluated"] == 3
+
+    def test_total_pairs(self, candidates, driver):
+        """3 candidates × 1 driver = 3 pairs."""
+        result = lem_prescreen_candidates(
+            candidates=candidates,
+            drivers=[driver],
+            target_f_low=500,
+            target_f_high=4000,
+            sim_freq_range=(354, 5657),
+            num_frequencies=30,
+            top_n=10,
+        )
+        assert result["total_pairs"] == 3
+
+    def test_multiple_drivers(self, candidates):
+        drivers = [_make_driver("d1"), _make_driver("d2")]
+        result = lem_prescreen_candidates(
+            candidates=candidates,
+            drivers=drivers,
+            target_f_low=500,
+            target_f_high=4000,
+            sim_freq_range=(354, 5657),
+            num_frequencies=30,
+            top_n=10,
+        )
+        assert result["total_pairs"] == 6  # 3 candidates × 2 drivers
+
+    def test_filtered_ids_are_valid(self, candidates, driver):
+        result = lem_prescreen_candidates(
+            candidates=candidates,
+            drivers=[driver],
+            target_f_low=500,
+            target_f_high=4000,
+            sim_freq_range=(354, 5657),
+            num_frequencies=30,
+            top_n=10,
+        )
+        valid_ids = {c.candidate_id for c in candidates}
+        for cid in result["filtered_candidate_ids"]:
+            assert cid in valid_ids
+
+    def test_rankings_sorted_by_score(self, candidates, driver):
+        result = lem_prescreen_candidates(
+            candidates=candidates,
+            drivers=[driver],
+            target_f_low=500,
+            target_f_high=4000,
+            sim_freq_range=(354, 5657),
+            num_frequencies=30,
+            top_n=10,
+        )
+        scores = [r["composite_score"] for r in result["rankings"]]
+        assert scores == sorted(scores, reverse=True)
+
+    def test_ranking_entries_have_expected_keys(self, candidates, driver):
+        result = lem_prescreen_candidates(
+            candidates=candidates,
+            drivers=[driver],
+            target_f_low=500,
+            target_f_high=4000,
+            sim_freq_range=(354, 5657),
+            num_frequencies=30,
+            top_n=10,
+        )
+        for entry in result["rankings"]:
+            assert "candidate_id" in entry
+            assert "driver_id" in entry
+            assert "composite_score" in entry
+            assert "profile" in entry
+
+
+class TestKnownGoodGeometry:
+    """A well-matched horn should rank higher than a poorly-matched one."""
+
+    def test_good_beats_bad(self):
+        driver = _make_driver(fs_hz=300, sd_m2=0.0008)
+        good = _make_candidate("good", "exponential", 0.025, 0.15, 0.3)
+        bad = _make_candidate("bad", "conical", 0.025, 0.03, 0.05)  # tiny horn
+
+        result = lem_prescreen_candidates(
+            candidates=[good, bad],
+            drivers=[driver],
+            target_f_low=500,
+            target_f_high=4000,
+            sim_freq_range=(354, 5657),
+            num_frequencies=50,
+            top_n=2,
+        )
+        # The well-matched horn should appear first in filtered IDs
+        assert result["filtered_candidate_ids"][0] == "good"
+
+
+class TestCsvIO:
+    def test_roundtrip(self, tmp_path):
+        candidates = [
+            _make_candidate("c1", "conical"),
+            _make_candidate("c2", "exponential"),
+            _make_candidate("c3", "hyperbolic"),
+        ]
+        csv_path = str(tmp_path / "candidates.csv")
+        # Write
+        import csv
+        with open(csv_path, "w", newline="") as f:
+            writer = csv.writer(f)
+            writer.writerow(["candidate_id", "profile", "throat_radius", "mouth_radius", "length"])
+            for c in candidates:
+                writer.writerow([c.candidate_id, c.profile, c.throat_radius, c.mouth_radius, c.length])
+
+        loaded = _load_candidates_csv(csv_path)
+        assert len(loaded) == 3
+        assert loaded[0].candidate_id == "c1"
+
+    def test_filtered_csv(self, tmp_path):
+        candidates = [
+            _make_candidate("c1", "conical"),
+            _make_candidate("c2", "exponential"),
+            _make_candidate("c3", "hyperbolic"),
+        ]
+        out_path = str(tmp_path / "filtered.csv")
+        _write_filtered_csv(candidates, ["c1", "c3"], out_path)
+
+        loaded = _load_candidates_csv(out_path)
+        assert len(loaded) == 2
+        ids = {c.candidate_id for c in loaded}
+        assert ids == {"c1", "c3"}
diff --git a/packages/horn-analysis/tests/test_prescreen.py b/packages/horn-analysis/tests/test_prescreen.py
index 1314539..d8aa540 100644
--- a/packages/horn-analysis/tests/test_prescreen.py
+++ b/packages/horn-analysis/tests/test_prescreen.py
@@ -6,7 +6,14 @@
 import numpy as np
 
 from horn_core.parameters import DriverParameters
-from horn_analysis.prescreen import PrescreenConfig, PrescreenResult, prescreen_drivers
+from horn_analysis.prescreen import (
+    PrescreenConfig,
+    PrescreenResult,
+    prescreen_drivers,
+    _parse_diameter_inches,
+)
+
+_SPEED_OF_SOUND = 343.0
 
 
 def _make_driver(
@@ -16,6 +23,7 @@ def _make_driver(
     qms=5.0,
     sd_m2=0.0008,
     driver_type="compression",
+    nominal_diameter=None,
     **kwargs,
 ):
     """Helper to create a DriverParameters with sensible defaults."""
@@ -32,6 +40,7 @@ def _make_driver(
         qms=qms,
         qes=qes,
         driver_type=driver_type,
+        nominal_diameter=nominal_diameter,
     )
 
 
@@ -47,52 +56,87 @@ def config():
 
 class TestPrescreenDrivers:
     def test_basic_filtering(self, config):
-        """Drivers with suitable fs, EBP, and type should pass."""
+        """Drivers with suitable fs and EBP should pass."""
         good = _make_driver("good", fs_hz=400, qes=0.4, driver_type="compression")
-        # EBP = 400/0.4 = 1000 > 50 -- passes
         result = prescreen_drivers([good], config)
         assert result.count == 1
         assert result.drivers[0].driver_id == "good"
 
     def test_fs_too_high_filtered(self, config):
         """Driver with fs above target_f_low * 1.5 should be filtered."""
-        # target_f_low = 500, threshold = 750
         bad = _make_driver("bad_fs", fs_hz=800, qes=0.4)
         result = prescreen_drivers([bad], config)
         assert result.count == 0
 
     def test_ebp_too_low_filtered(self, config):
         """Driver with EBP below threshold should be filtered."""
-        # EBP = 400 / 10.0 = 40 < 50
         bad = _make_driver("bad_ebp", fs_hz=400, qes=10.0)
         result = prescreen_drivers([bad], config)
         assert result.count == 0
 
-    def test_large_cone_filtered_for_high_freq(self, config):
-        """Large cone (e.g. 15") should be filtered when f_piston*factor < target_f_high."""
-        # 15" driver: Sd ~ 855 cm² = 0.0855 m²
-        # f_piston = 343/(2π·√(0.0855/π)) = ~331 Hz, ×10 = 3310 < 4000 → filtered
+    def test_driver_larger_than_mouth_filtered(self, config):
+        """Driver with effective radius > mouth_radius should be filtered.
+
+        config has mouth_radius=0.2m.
+        18" cone: sqrt(0.1/π) = 0.178m < 0.2m → passes.
+        Giant driver: sqrt(0.15/π) = 0.219m > 0.2m → filtered.
+        """
+        giant = _make_driver("giant", fs_hz=400, qes=0.4, sd_m2=0.15, driver_type="cone")
+        result = prescreen_drivers([giant], config)
+        assert result.count == 0
+
+    def test_large_cone_now_passes_with_decoupled_throat(self, config):
+        """15" cone now passes since throat is decoupled from driver.
+
+        Old behavior: max_throat = 0.2/2.0 = 0.1m, 15" filtered.
+        New behavior: max_driver = 0.2m (mouth), 15" cone radius=0.165m < 0.2m → passes.
+        """
         big_cone = _make_driver("big_cone", fs_hz=400, qes=0.4, sd_m2=0.0855, driver_type="cone")
         result = prescreen_drivers([big_cone], config)
-        assert result.count == 0
+        assert result.count == 1
 
-    def test_small_cone_passes_high_freq(self, config):
-        """Small cone (e.g. 6") should pass when f_piston*factor >= target_f_high."""
-        # 6" driver: Sd ~ 130 cm² = 0.0130 m²
-        # f_piston = 343/(2π·√(0.0130/π)) = ~849 Hz, ×10 = 8490 > 4000 → passes
+    def test_small_cone_passes(self, config):
+        """6" cone should pass: fits in mouth easily."""
         small_cone = _make_driver("small_cone", fs_hz=400, qes=0.4, sd_m2=0.0130, driver_type="cone")
         result = prescreen_drivers([small_cone], config)
         assert result.count == 1
 
+    def test_8inch_cone_passes_low_freq_horn(self):
+        """8" cone should pass for 200-2kHz horn with mouth_radius=0.25m."""
+        config = PrescreenConfig(
+            target_f_low_hz=200,
+            target_f_high_hz=2000,
+            mouth_radius_m=0.25,
+            length_m=0.4,
+        )
+        cone_8 = _make_driver("cone_8", fs_hz=100, qes=0.4, sd_m2=0.0214, driver_type="cone")
+        result = prescreen_drivers([cone_8], config)
+        assert result.count == 1
+
+    def test_8inch_cone_passes_small_mouth(self):
+        """8" cone now passes even with smaller mouth (decoupled throat).
+
+        mouth_radius=0.15m → max_driver = 0.15m.
+        8" cone: sqrt(0.0214/π) = 0.0825m < 0.15m → passes.
+        """
+        config = PrescreenConfig(
+            target_f_low_hz=300,
+            target_f_high_hz=7000,
+            mouth_radius_m=0.15,
+            length_m=0.3,
+        )
+        cone_8 = _make_driver("cone_8", fs_hz=200, qes=0.4, sd_m2=0.0214, driver_type="cone")
+        result = prescreen_drivers([cone_8], config)
+        assert result.count == 1
+
     def test_compression_passes_low_freq(self):
-        """Compression drivers should pass low-freq targets (high f_piston)."""
+        """Compression drivers should pass low-freq targets."""
         config = PrescreenConfig(
             target_f_low_hz=100,
             target_f_high_hz=1000,
             mouth_radius_m=0.3,
             length_m=0.8,
         )
-        # Sd=20cm²=0.0020 m², f_piston ~2186 Hz, ×10 = 21860 > 1000
         comp = _make_driver("comp", fs_hz=80, qes=0.4, sd_m2=0.0020, driver_type="compression")
         result = prescreen_drivers([comp], config)
         assert result.count == 1
@@ -110,29 +154,6 @@ def test_both_types_pass_mid_band(self):
         result = prescreen_drivers([comp, cone], config)
         assert result.count == 2
 
-    def test_representative_throat_radius(self, config):
-        """Throat radius should be median of sqrt(Sd/pi)."""
-        d1 = _make_driver("d1", fs_hz=400, qes=0.4, sd_m2=0.0008)
-        d2 = _make_driver("d2", fs_hz=400, qes=0.4, sd_m2=0.0012)
-        d3 = _make_driver("d3", fs_hz=400, qes=0.4, sd_m2=0.0010)
-
-        result = prescreen_drivers([d1, d2, d3], config)
-        assert result.count == 3
-
-        expected_radii = sorted([math.sqrt(s / math.pi) for s in [0.0008, 0.0012, 0.0010]])
-        expected_median = expected_radii[1]  # median of 3
-        assert result.throat_radius_m == pytest.approx(expected_median, rel=1e-3)
-
-    def test_sd_ratio_filtering(self, config):
-        """Drivers with extreme Sd ratios should be filtered."""
-        # Create 3 normal drivers and 1 with very different Sd
-        normal = [_make_driver(f"n{i}", fs_hz=400, qes=0.4, sd_m2=0.001) for i in range(3)]
-        outlier = _make_driver("outlier", fs_hz=400, qes=0.4, sd_m2=0.05)  # 50x larger
-
-        result = prescreen_drivers(normal + [outlier], config)
-        ids = [d.driver_id for d in result.drivers]
-        assert "outlier" not in ids
-
     def test_empty_driver_list(self, config):
         """Empty input should return empty result."""
         result = prescreen_drivers([], config)
@@ -147,12 +168,266 @@ def test_all_filtered_out(self, config):
         assert result.count == 0
 
     def test_to_dict(self, config):
-        """PrescreenResult.to_dict should contain driver IDs."""
+        """PrescreenResult.to_dict should contain driver IDs and throat_radii_m."""
         d = _make_driver("test1", fs_hz=400, qes=0.4)
         result = prescreen_drivers([d], config)
         d_dict = result.to_dict()
         assert "drivers" in d_dict
         assert "throat_radius_m" in d_dict
+        assert "throat_radii_m" in d_dict
         assert "count" in d_dict
         assert d_dict["count"] == 1
         assert "test1" in d_dict["drivers"]
+
+    def test_10inch_cone_passes_with_adequate_mouth(self):
+        """10" cone passes when mouth is large enough (decoupled throat).
+
+        mouth_radius=0.2m → max_driver = 0.2m.
+        10" cone: sqrt(0.0346/π) = 0.105m < 0.2m → passes.
+        """
+        config = PrescreenConfig(
+            target_f_low_hz=300,
+            target_f_high_hz=7000,
+            mouth_radius_m=0.2,
+            length_m=0.3,
+        )
+        cone_10 = _make_driver("cone_10", fs_hz=50, qes=0.35, sd_m2=0.0346, driver_type="cone")
+        result = prescreen_drivers([cone_10], config)
+        assert result.count == 1
+
+    def test_compression_driver_small_exit_passes(self, config):
+        """Compression driver with small exit area passes."""
+        comp = _make_driver("comp_small", fs_hz=400, qes=0.4, sd_m2=0.0020,
+                            driver_type="compression")
+        comp.exit_area_m2 = 0.001
+        result = prescreen_drivers([comp], config)
+        assert result.count == 1
+
+    def test_6inch_cone_passes_mid_freq(self):
+        """6" cone should pass for a mid-frequency target with adequate mouth."""
+        config = PrescreenConfig(
+            target_f_low_hz=300,
+            target_f_high_hz=3000,
+            mouth_radius_m=0.2,
+            length_m=0.4,
+        )
+        cone_6 = _make_driver("cone_6", fs_hz=80, qes=0.4, sd_m2=0.0130, driver_type="cone")
+        result = prescreen_drivers([cone_6], config)
+        assert result.count == 1
+
+    def test_8inch_cone_passes_low_freq_large_mouth(self):
+        """8" cone passes for low-frequency target with large mouth."""
+        config = PrescreenConfig(
+            target_f_low_hz=200,
+            target_f_high_hz=2000,
+            mouth_radius_m=0.30,
+            length_m=0.5,
+        )
+        cone_8 = _make_driver("cone_8", fs_hz=80, qes=0.4, sd_m2=0.0214, driver_type="cone")
+        result = prescreen_drivers([cone_8], config)
+        assert result.count == 1
+
+    def test_max_driver_derived_from_target_freq(self):
+        """When mouth_radius is None, max_driver is derived from ideal mouth at f_low.
+
+        f_low=300 Hz → ideal_mouth = 343/(2π·300) = 0.182m
+        6" cone: radius=0.0643m < 0.182m → passes.
+        Giant: radius=0.219m > 0.182m → filtered.
+        """
+        config = PrescreenConfig(
+            target_f_low_hz=300,
+            target_f_high_hz=3000,
+            mouth_radius_m=None,
+        )
+        cone_6 = _make_driver("cone_6", fs_hz=80, qes=0.4, sd_m2=0.0130, driver_type="cone")
+        giant = _make_driver("giant", fs_hz=80, qes=0.4, sd_m2=0.15, driver_type="cone")
+        result = prescreen_drivers([cone_6, giant], config)
+        ids = [d.driver_id for d in result.drivers]
+        assert "cone_6" in ids
+        assert "giant" not in ids
+
+    def test_effective_throat_area_with_exit_area(self, config):
+        """Driver with exit_area_m2 should use it for driver radius calculation."""
+        d = _make_driver("phase_plug", fs_hz=400, qes=0.4, sd_m2=0.0214)
+        d.exit_area_m2 = 0.001
+        result = prescreen_drivers([d], config)
+        assert result.count == 1
+
+
+class TestAcousticThroatRadii:
+    """Test that throat radii are derived from acoustic constraints."""
+
+    def test_throat_radii_from_acoustics(self):
+        """Throat radii should be based on ka_max at f_high, not driver Sd."""
+        config = PrescreenConfig(
+            target_f_low_hz=300,
+            target_f_high_hz=10000,
+            mouth_radius_m=0.2,
+        )
+        # ka_max = 2π, f_high = 10kHz
+        # a_acoustic_max = 343 * 2π / (2π * 10000) = 343/10000 = 0.0343m
+        a_max = _SPEED_OF_SOUND * config.ka_max / (2 * math.pi * config.target_f_high_hz)
+
+        d = _make_driver("comp", fs_hz=400, qes=0.4, sd_m2=0.0008)
+        result = prescreen_drivers([d], config)
+
+        # Acoustic radii at 0.3, 0.65, 1.0 of a_max should be present
+        assert len(result.throat_radii_m) >= 3
+        assert result.throat_radii_m[-1] == pytest.approx(round(a_max, 6), rel=1e-3)
+
+    def test_throat_radii_include_acoustic_fractions(self, ):
+        """throat_radii_m should include 0.3, 0.65, 1.0 fractions of a_acoustic_max."""
+        config = PrescreenConfig(
+            target_f_low_hz=500,
+            target_f_high_hz=4000,
+            mouth_radius_m=0.2,
+        )
+        a_max = _SPEED_OF_SOUND * config.ka_max / (2 * math.pi * config.target_f_high_hz)
+        expected_acoustic = [round(a_max * f, 6) for f in [0.3, 0.65, 1.0]]
+
+        d = _make_driver("comp", fs_hz=400, qes=0.4, sd_m2=0.0008)
+        result = prescreen_drivers([d], config)
+
+        for r in expected_acoustic:
+            assert r in result.throat_radii_m
+
+    def test_throat_radii_include_direct_coupling(self):
+        """Small drivers that fit within a_acoustic_max should appear in throat_radii."""
+        config = PrescreenConfig(
+            target_f_low_hz=500,
+            target_f_high_hz=4000,
+            mouth_radius_m=0.2,
+        )
+        a_max = _SPEED_OF_SOUND * config.ka_max / (2 * math.pi * config.target_f_high_hz)
+
+        # Small compression driver: radius well below a_max
+        d = _make_driver("small_comp", fs_hz=400, qes=0.4, sd_m2=0.0008)
+        drv_radius = round(math.sqrt(0.0008 / math.pi), 6)
+        assert drv_radius < a_max  # sanity check
+
+        result = prescreen_drivers([d], config)
+        assert drv_radius in result.throat_radii_m
+
+    def test_large_driver_radius_excluded_from_throat_radii(self):
+        """Driver radii exceeding a_acoustic_max should not appear in throat_radii."""
+        config = PrescreenConfig(
+            target_f_low_hz=300,
+            target_f_high_hz=10000,
+            mouth_radius_m=0.2,
+        )
+        a_max = _SPEED_OF_SOUND * config.ka_max / (2 * math.pi * config.target_f_high_hz)
+        # 0.0343m max. 6" cone radius = 0.0643m > a_max
+        cone_6 = _make_driver("cone_6", fs_hz=200, qes=0.4, sd_m2=0.0130, driver_type="cone")
+        drv_radius = round(math.sqrt(0.0130 / math.pi), 6)
+        assert drv_radius > a_max
+
+        result = prescreen_drivers([cone_6], config)
+        assert drv_radius not in result.throat_radii_m
+
+    def test_representative_is_middle_of_radii(self):
+        """Representative radius should be the middle element of throat_radii_m."""
+        config = PrescreenConfig(
+            target_f_low_hz=500,
+            target_f_high_hz=4000,
+            mouth_radius_m=0.2,
+        )
+        d = _make_driver("comp", fs_hz=400, qes=0.4, sd_m2=0.0008)
+        result = prescreen_drivers([d], config)
+        mid_idx = len(result.throat_radii_m) // 2
+        assert result.throat_radius_m == result.throat_radii_m[mid_idx]
+
+    def test_throat_radii_capped_at_five(self):
+        """At most 5 throat radii to limit combinatorial explosion."""
+        config = PrescreenConfig(
+            target_f_low_hz=500,
+            target_f_high_hz=4000,
+            mouth_radius_m=0.2,
+        )
+        # Create many drivers with varied Sd to generate many direct_radii
+        drivers = [
+            _make_driver(f"d{i}", fs_hz=400, qes=0.4, sd_m2=0.0001 * (i + 1))
+            for i in range(20)
+        ]
+        result = prescreen_drivers(drivers, config)
+        assert len(result.throat_radii_m) <= 5
+
+
+class TestDiameterFilter:
+    """Test optional nominal diameter filtering."""
+
+    def test_min_diameter_filter(self):
+        """Drivers below min_nominal_diameter_in should be filtered."""
+        config = PrescreenConfig(
+            target_f_low_hz=300,
+            target_f_high_hz=3000,
+            mouth_radius_m=0.2,
+            min_nominal_diameter_in=4.0,
+        )
+        small = _make_driver("small", fs_hz=200, qes=0.4, sd_m2=0.0020,
+                             nominal_diameter="2in")
+        big = _make_driver("big", fs_hz=200, qes=0.4, sd_m2=0.0130,
+                           nominal_diameter="6in")
+        result = prescreen_drivers([small, big], config)
+        ids = [d.driver_id for d in result.drivers]
+        assert "small" not in ids
+        assert "big" in ids
+
+    def test_max_diameter_filter(self):
+        """Drivers above max_nominal_diameter_in should be filtered."""
+        config = PrescreenConfig(
+            target_f_low_hz=300,
+            target_f_high_hz=3000,
+            mouth_radius_m=0.3,
+            max_nominal_diameter_in=6.0,
+        )
+        small = _make_driver("small", fs_hz=200, qes=0.4, sd_m2=0.0020,
+                             nominal_diameter="2in")
+        big = _make_driver("big", fs_hz=200, qes=0.4, sd_m2=0.0346,
+                           nominal_diameter="10in")
+        result = prescreen_drivers([small, big], config)
+        ids = [d.driver_id for d in result.drivers]
+        assert "small" in ids
+        assert "big" not in ids
+
+    def test_diameter_filter_skips_missing_diameter(self):
+        """Drivers without nominal_diameter should not be filtered by diameter."""
+        config = PrescreenConfig(
+            target_f_low_hz=300,
+            target_f_high_hz=3000,
+            mouth_radius_m=0.2,
+            min_nominal_diameter_in=4.0,
+        )
+        no_dia = _make_driver("no_dia", fs_hz=200, qes=0.4, sd_m2=0.0020)
+        result = prescreen_drivers([no_dia], config)
+        assert result.count == 1
+
+    def test_no_diameter_filter_by_default(self):
+        """Without diameter config, all drivers pass diameter check."""
+        config = PrescreenConfig(
+            target_f_low_hz=300,
+            target_f_high_hz=3000,
+            mouth_radius_m=0.2,
+        )
+        d = _make_driver("d", fs_hz=200, qes=0.4, sd_m2=0.0020, nominal_diameter="1in")
+        result = prescreen_drivers([d], config)
+        assert result.count == 1
+
+
+class TestParseDiameterInches:
+    def test_simple_integer(self):
+        assert _parse_diameter_inches("4in") == 4.0
+
+    def test_float_value(self):
+        assert _parse_diameter_inches("6.5in") == 6.5
+
+    def test_bare_number(self):
+        assert _parse_diameter_inches("8") == 8.0
+
+    def test_none_input(self):
+        assert _parse_diameter_inches(None) is None
+
+    def test_empty_string(self):
+        assert _parse_diameter_inches("") is None
+
+    def test_no_match(self):
+        assert _parse_diameter_inches("abc") is None
diff --git a/packages/horn-analysis/tests/test_scoring.py b/packages/horn-analysis/tests/test_scoring.py
index 9893852..84032be 100644
--- a/packages/horn-analysis/tests/test_scoring.py
+++ b/packages/horn-analysis/tests/test_scoring.py
@@ -89,6 +89,36 @@ def test_low_sensitivity(self):
         assert score.avg_sensitivity_db == 80.0
 
 
+class TestBandwidthFloor:
+    def test_low_coverage_zeroes_composite(self):
+        """Bandwidth coverage < 10% should zero out the composite score."""
+        # Horn covers only 3.6% of target: (530-770) / (300-7000) = 240/6700 ≈ 3.6%
+        kpi = _make_kpi(
+            f3_low_hz=530.0,
+            f3_high_hz=770.0,
+            passband_ripple_db=1.0,
+            average_sensitivity_db=100.0,
+        )
+        target = TargetSpec(f_low_hz=300.0, f_high_hz=7000.0)
+        score = compute_selection_score(kpi, target)
+        assert score.bandwidth_coverage < 0.10
+        assert score.composite_score == 0.0
+
+    def test_above_floor_scores_normally(self):
+        """Bandwidth coverage >= 10% should score normally (not zeroed)."""
+        # Horn covers ~14%: (500-1450) / (500-7000) = 950/6500 ≈ 14.6%
+        kpi = _make_kpi(
+            f3_low_hz=500.0,
+            f3_high_hz=1450.0,
+            passband_ripple_db=2.0,
+            average_sensitivity_db=92.0,
+        )
+        target = TargetSpec(f_low_hz=500.0, f_high_hz=7000.0)
+        score = compute_selection_score(kpi, target)
+        assert score.bandwidth_coverage >= 0.10
+        assert score.composite_score > 0.0
+
+
 class TestCompositeScore:
     def test_perfect_score(self):
         """Full coverage, zero ripple, max sensitivity → composite ≈ 1.0."""
diff --git a/packages/horn-core/pyproject.toml b/packages/horn-core/pyproject.toml
index dfcbe84..c670028 100644
--- a/packages/horn-core/pyproject.toml
+++ b/packages/horn-core/pyproject.toml
@@ -3,7 +3,8 @@ name = "horn-core"
 version = "0.1.0"
 description = "Core data structures for the horn simulation project."
 dependencies = [
-    "numpy"
+    "numpy",
+    "scipy",
 ]
 
 [project.scripts]
diff --git a/packages/horn-core/src/horn_core/geometry_designer.py b/packages/horn-core/src/horn_core/geometry_designer.py
index c736f6a..fe1b9d4 100644
--- a/packages/horn-core/src/horn_core/geometry_designer.py
+++ b/packages/horn-core/src/horn_core/geometry_designer.py
@@ -87,6 +87,89 @@ def derive_simulation_freq_range(
     return (sim_min, sim_max)
 
 
+def generate_auto_candidates(
+    target_f_low: float,
+    target_f_high: float,
+    throat_radii: List[float],
+    mouth_radius: Optional[float] = None,
+    length: Optional[float] = None,
+    num_mouth_radii: int = 3,
+    num_lengths: int = 3,
+    profiles: Optional[List[str]] = None,
+) -> tuple:
+    """Generate geometry candidates for unified auto mode.
+
+    When mouth_radius or length is provided, it is fixed (1 value).
+    When None, the range is derived from target_f_low.
+
+    Args:
+        target_f_low: Target low-frequency cutoff (Hz).
+        target_f_high: Target high-frequency cutoff (Hz).
+        throat_radii: Throat radii (m) from prescreen (1-3 values).
+        mouth_radius: Fixed mouth radius (m), or None to derive range.
+        length: Fixed horn length (m), or None to derive range.
+        num_mouth_radii: Grid points for mouth radius (used when mouth_radius is None).
+        num_lengths: Grid points for length (used when length is None).
+        profiles: Horn profile types. Defaults to all seven.
+
+    Returns:
+        Tuple of (candidates, derived_geometry).
+    """
+    if profiles is None:
+        profiles = DEFAULT_PROFILES
+
+    # Mouth radius: fixed or derived
+    if mouth_radius is not None:
+        mouth_range = (mouth_radius, mouth_radius)
+        mouth_radii = [mouth_radius]
+    else:
+        mouth_range = derive_mouth_radius_range(target_f_low)
+        mouth_radii = np.linspace(mouth_range[0], mouth_range[1], num_mouth_radii).tolist()
+
+    # Length: fixed or derived
+    if length is not None:
+        length_range = (length, length)
+        lengths = [length]
+    else:
+        length_range = derive_length_range(target_f_low)
+        lengths = np.linspace(length_range[0], length_range[1], num_lengths).tolist()
+
+    sim_range = derive_simulation_freq_range(target_f_low, target_f_high)
+
+    candidates = []
+    idx = 0
+    for profile in profiles:
+        for r_throat in throat_radii:
+            for r_mouth in mouth_radii:
+                for horn_length in lengths:
+                    if r_mouth <= r_throat:
+                        continue
+                    candidate_id = f"auto_{profile[:3]}_{idx:04d}"
+                    candidates.append(
+                        CandidateGeometry(
+                            candidate_id=candidate_id,
+                            profile=profile,
+                            throat_radius=r_throat,
+                            mouth_radius=round(r_mouth, 6),
+                            length=round(horn_length, 6),
+                        )
+                    )
+                    idx += 1
+
+    ideal_mouth = derive_mouth_radius(target_f_low)
+    derived = DerivedGeometry(
+        target_f_low=target_f_low,
+        target_f_high=target_f_high,
+        ideal_mouth_radius=ideal_mouth,
+        mouth_radius_range=mouth_range,
+        length_range=length_range,
+        sim_freq_range=sim_range,
+        candidate_count=len(candidates),
+    )
+
+    return candidates, derived
+
+
 def generate_fullauto_candidates(
     target_f_low: float,
     target_f_high: float,
@@ -173,6 +256,14 @@ def main():
         required=True,
         help="Prescreen result JSON (provides throat_radius_m).",
     )
+    parser.add_argument(
+        "--mouth-radius", type=float, default=None,
+        help="Fixed mouth radius (m). If omitted, derive range from f_low.",
+    )
+    parser.add_argument(
+        "--length", type=float, default=None,
+        help="Fixed horn length (m). If omitted, derive range from f_low.",
+    )
     parser.add_argument(
         "--num-mouth-radii", type=int, default=3, help="Mouth radius grid points."
     )
@@ -191,12 +282,14 @@ def main():
     args = parser.parse_args()
 
     prescreen = json.loads(Path(args.prescreen_json).read_text())
-    throat_radius = prescreen["throat_radius_m"]
+    throat_radii = prescreen.get("throat_radii_m") or [prescreen["throat_radius_m"]]
 
-    candidates, derived = generate_fullauto_candidates(
+    candidates, derived = generate_auto_candidates(
         target_f_low=args.target_f_low,
         target_f_high=args.target_f_high,
-        throat_radii=[throat_radius],
+        throat_radii=throat_radii,
+        mouth_radius=args.mouth_radius,
+        length=args.length,
         num_mouth_radii=args.num_mouth_radii,
         num_lengths=args.num_lengths,
     )
diff --git a/packages/horn-core/src/horn_core/parameters.py b/packages/horn-core/src/horn_core/parameters.py
index c4acf17..b34785e 100644
--- a/packages/horn-core/src/horn_core/parameters.py
+++ b/packages/horn-core/src/horn_core/parameters.py
@@ -56,11 +56,22 @@ class DriverParameters:
     cms_m_per_n: Optional[float] = field(default=None, repr=False)
     rms_kg_per_s: Optional[float] = field(default=None, repr=False)
 
+    # Phase plug exit area (m²) — for drivers like 8" cones with phase plugs
+    # where the effective throat area is much smaller than Sd
+    exit_area_m2: Optional[float] = None
+
     # Optional metadata
     driver_type: Optional[str] = None
     nominal_diameter: Optional[str] = None  # e.g. "18in", "15in", "1in"
     xmax_m: Optional[float] = None
     nominal_impedance_ohm: Optional[float] = None
+    power_w: Optional[float] = None           # RMS / AES / continuous power (W)
+    peak_power_w: Optional[float] = None      # Peak / program power (W)
+
+    @property
+    def effective_throat_area(self) -> float:
+        """Return the effective throat area: exit_area_m2 if set, else sd_m2."""
+        return self.exit_area_m2 if self.exit_area_m2 is not None else self.sd_m2
 
     def __post_init__(self):
         omega_s = 2.0 * np.pi * self.fs_hz
diff --git a/packages/horn-core/src/horn_core/profiles.py b/packages/horn-core/src/horn_core/profiles.py
new file mode 100644
index 0000000..1fffd7a
--- /dev/null
+++ b/packages/horn-core/src/horn_core/profiles.py
@@ -0,0 +1,133 @@
+"""Pure-Python radius functions for all 7 horn flare profiles.
+
+Extracted from horn-geometry/generator.py without the gmsh dependency.
+Only requires numpy. Each function maps axial position z ∈ [0, L] to
+the local horn radius r(z).
+"""
+
+import math
+from typing import Callable
+
+import numpy as np
+
+
+def conical_radius(z: float, r_t: float, r_m: float, L: float) -> float:
+    """Linear flare: r(z) = r_t + (r_m - r_t) * z / L."""
+    return r_t + (r_m - r_t) * z / L
+
+
+def exponential_radius(z: float, r_t: float, r_m: float, L: float) -> float:
+    """Exponential flare: r(z) = r_t * exp(m * z / L), m = ln(r_m / r_t)."""
+    m = math.log(r_m / r_t)
+    return r_t * math.exp(m * z / L)
+
+
+def hyperbolic_radius(z: float, r_t: float, r_m: float, L: float) -> float:
+    """Hyperbolic (hypex) flare: r(z) = r_t * cosh(m * z / L), m = acosh(r_m / r_t)."""
+    m = np.arccosh(r_m / r_t)
+    return float(r_t * np.cosh(m * z / L))
+
+
+def _build_tractrix_interp(r_t: float, r_m: float, L: float):
+    """Pre-compute the tractrix interpolation arrays."""
+    t = np.linspace(np.pi - 1e-6, np.pi / 2, 500)
+    y, x = np.sin(t), np.log(np.tan(t / 2)) + np.cos(t)
+    x -= x[0]
+    x_n, y_n = x / x[-1], y / y[-1]
+    return x_n * L, y_n
+
+
+def tractrix_radius(z: float, r_t: float, r_m: float, L: float) -> float:
+    """Tractrix flare: rapid initial expansion, decelerating toward mouth."""
+    x_scaled, y_n = _build_tractrix_interp(r_t, r_m, L)
+    return r_t + (r_m - r_t) * float(np.interp(z, x_scaled, y_n))
+
+
+def os_radius(z: float, r_t: float, r_m: float, L: float) -> float:
+    """Oblate spheroidal (OS / Geddes) flare: r(z) = sqrt(r_t² + (z·tanθ)²)."""
+    theta = math.atan2(math.sqrt(r_m**2 - r_t**2), L)
+    return math.sqrt(r_t**2 + (z * math.tan(theta)) ** 2)
+
+
+def _build_lecleach_interp(r_t: float, r_m: float, L: float):
+    """Pre-compute the Le Cléac'h interpolation arrays."""
+    t = np.linspace(np.pi - 1e-6, np.pi / 2, 500)
+    y, x = np.sin(t), np.log(np.tan(t / 2)) + np.cos(t)
+    x -= x[0]
+    idx = np.searchsorted(y, r_t / r_m)
+    x_c, y_c = x[idx:] - x[idx], y[idx:]
+    return x_c / x_c[-1] * L, y_c / y_c[-1] * r_m
+
+
+def lecleach_radius(z: float, r_t: float, r_m: float, L: float) -> float:
+    """Le Cléac'h flare: tractrix with a = r_m, clipped at r >= r_t."""
+    x_scaled, y_scaled = _build_lecleach_interp(r_t, r_m, L)
+    return float(np.interp(z, x_scaled, y_scaled))
+
+
+def cd_radius(z: float, r_t: float, r_m: float, L: float) -> float:
+    """Constant Directivity (CD) horn: exponential throat (30%) + conical body."""
+    frac = 0.3
+    z_t = frac * L
+    r_trans = r_t * (r_m / r_t) ** frac
+
+    if z <= z_t:
+        return r_t * math.exp(math.log(r_trans / r_t) * z / z_t)
+    else:
+        return r_trans + (r_m - r_trans) * (z - z_t) / (L - z_t)
+
+
+_PROFILE_DISPATCH = {
+    "conical": conical_radius,
+    "exponential": exponential_radius,
+    "hyperbolic": hyperbolic_radius,
+    "tractrix": tractrix_radius,
+    "os": os_radius,
+    "lecleach": lecleach_radius,
+    "cd": cd_radius,
+}
+
+
+def get_radius_func(
+    profile: str, r_t: float, r_m: float, L: float
+) -> Callable[[float], float]:
+    """Return a closure r(z) for the given profile and geometry.
+
+    Args:
+        profile: One of "conical", "exponential", "hyperbolic", "tractrix",
+                 "os", "lecleach", "cd".
+        r_t: Throat radius (m).
+        r_m: Mouth radius (m).
+        L: Horn length (m).
+
+    Returns:
+        Callable that maps z ∈ [0, L] → radius (m).
+    """
+    if profile not in _PROFILE_DISPATCH:
+        raise ValueError(
+            f"Unknown profile '{profile}'. Choose from: {list(_PROFILE_DISPATCH.keys())}"
+        )
+
+    raw_func = _PROFILE_DISPATCH[profile]
+
+    # For profiles with expensive interpolation setup, pre-compute once
+    if profile == "tractrix":
+        x_scaled, y_n = _build_tractrix_interp(r_t, r_m, L)
+
+        def _tractrix(z: float) -> float:
+            return r_t + (r_m - r_t) * float(np.interp(z, x_scaled, y_n))
+
+        return _tractrix
+
+    if profile == "lecleach":
+        x_scaled, y_scaled = _build_lecleach_interp(r_t, r_m, L)
+
+        def _lecleach(z: float) -> float:
+            return float(np.interp(z, x_scaled, y_scaled))
+
+        return _lecleach
+
+    def _radius(z: float) -> float:
+        return raw_func(z, r_t, r_m, L)
+
+    return _radius
diff --git a/packages/horn-core/src/horn_core/webster.py b/packages/horn-core/src/horn_core/webster.py
new file mode 100644
index 0000000..0de650f
--- /dev/null
+++ b/packages/horn-core/src/horn_core/webster.py
@@ -0,0 +1,144 @@
+"""Transfer Matrix Method (TMM) solver for horn throat impedance.
+
+Slices the horn into N thin cylindrical segments, builds a 2×2 transfer
+matrix per segment using plane-wave propagation at the local cross-section
+area, and cascades from throat → mouth. A radiation impedance boundary
+condition is applied at the mouth.
+
+With fine segmentation (N ≥ 200) this converges to the Webster horn
+equation solution, giving millisecond-speed approximations of the FEM
+solver output suitable for prescreening large candidate grids.
+"""
+
+import numpy as np
+from scipy.special import j1, struve
+from typing import Callable, Tuple
+
+# Speed of sound and air density at ~20 °C
+C0 = 343.0
+RHO0 = 1.225
+
+
+def piston_radiation_impedance(k: float, a: float) -> complex:
+    """Flanged circular piston radiation impedance normalised by ρc.
+
+    Z_rad / (ρc) = R1(2ka) + j·X1(2ka)
+    where R1(x) = 1 - 2·J1(x)/x,  X1(x) = 2·H1(x)/x.
+
+    Duplicated from horn_solver/radiation.py to avoid dolfinx dependency.
+    """
+    ka = k * a
+    x = 2.0 * ka
+
+    if x < 1e-12:
+        r1 = x**2 / 8.0
+        x1 = 4.0 * x / (3.0 * np.pi)
+    else:
+        r1 = 1.0 - 2.0 * j1(x) / x
+        x1 = 2.0 * struve(1, x) / x
+
+    return complex(r1, x1)
+
+
+def unflanged_radiation_impedance(k: float, a: float) -> complex:
+    """Unflanged pipe radiation impedance (Levine-Schwinger), normalised by ρc."""
+    ka = k * a
+    return complex(0.25 * ka**2, 0.6133 * ka)
+
+
+def _segment_matrix(k: float, S: float, dz: float) -> np.ndarray:
+    """2×2 transfer matrix for a cylindrical segment (plane-wave).
+
+    Convention: [p_in, U_in] = T · [p_out, U_out]
+    where p is pressure and U is volume velocity.
+
+    Args:
+        k: Wavenumber (rad/m).
+        S: Cross-section area (m²).
+        dz: Segment length (m).
+    """
+    cos_kl = np.cos(k * dz)
+    sin_kl = np.sin(k * dz)
+    Zc = RHO0 * C0 / S  # characteristic acoustic impedance
+
+    return np.array([
+        [cos_kl, 1j * Zc * sin_kl],
+        [1j * sin_kl / Zc, cos_kl],
+    ], dtype=complex)
+
+
+def compute_throat_impedance_tmm(
+    frequencies: np.ndarray,
+    radius_func: Callable[[float], float],
+    length: float,
+    throat_radius: float,
+    mouth_radius: float,
+    n_segments: int = 200,
+    radiation_model: str = "flanged_piston",
+) -> Tuple[np.ndarray, np.ndarray]:
+    """Compute specific acoustic impedance at the horn throat using TMM.
+
+    Slices the horn into cylindrical segments at their midpoint radius,
+    builds the cascaded transfer matrix from throat to mouth, and applies
+    a radiation impedance boundary condition at the mouth.
+
+    Args:
+        frequencies: Array of frequencies in Hz.
+        radius_func: Callable mapping z ∈ [0, length] to radius (m).
+        length: Horn length (m).
+        throat_radius: Throat radius (m).
+        mouth_radius: Mouth radius (m).
+        n_segments: Number of segments (≥50 recommended).
+        radiation_model: "flanged_piston" or "unflanged".
+
+    Returns:
+        Tuple of (z_real, z_imag) arrays — specific acoustic impedance
+        at the throat in Pa·s/m, matching the FEM solver CSV format.
+    """
+    # Segment boundaries and midpoint radii
+    z_edges = np.linspace(0, length, n_segments + 1)
+    dz = length / n_segments
+    z_mid = 0.5 * (z_edges[:-1] + z_edges[1:])
+    r_mid = np.array([radius_func(z) for z in z_mid])
+    S_mid = np.pi * r_mid**2
+
+    z_real_out = np.zeros(len(frequencies))
+    z_imag_out = np.zeros(len(frequencies))
+
+    a_mouth = radius_func(length)
+    S_throat = np.pi * throat_radius**2
+
+    for i_f, freq in enumerate(frequencies):
+        k = 2.0 * np.pi * freq / C0
+
+        # Radiation impedance at the mouth (normalised → acoustic)
+        if radiation_model == "flanged_piston":
+            z_rad_norm = piston_radiation_impedance(k, a_mouth)
+        else:
+            z_rad_norm = unflanged_radiation_impedance(k, a_mouth)
+
+        S_mouth = np.pi * a_mouth**2
+        Z_load = z_rad_norm * RHO0 * C0 / S_mouth  # acoustic impedance
+
+        # Cascade transfer matrices: T_total = T_0 · T_1 · ... · T_{N-1}
+        # Maps (p, U) at throat to (p, U) at mouth
+        T = np.eye(2, dtype=complex)
+        for seg in range(n_segments):
+            T_seg = _segment_matrix(k, S_mid[seg], dz)
+            T = T @ T_seg
+
+        # Throat acoustic impedance from transfer matrix:
+        # p_throat = T11*p_mouth + T12*U_mouth
+        # U_throat = T21*p_mouth + T22*U_mouth
+        # At mouth: p_mouth = Z_load * U_mouth
+        # So: Z_throat = p_throat/U_throat
+        #   = (T11*Z_load + T12) / (T21*Z_load + T22)
+        Z_throat_acoustic = (T[0, 0] * Z_load + T[0, 1]) / (T[1, 0] * Z_load + T[1, 1])
+
+        # Convert acoustic → specific acoustic impedance
+        Z_throat_specific = Z_throat_acoustic * S_throat
+
+        z_real_out[i_f] = np.real(Z_throat_specific)
+        z_imag_out[i_f] = np.imag(Z_throat_specific)
+
+    return z_real_out, z_imag_out
diff --git a/packages/horn-core/tests/test_geometry_designer.py b/packages/horn-core/tests/test_geometry_designer.py
index c3013ca..2b29869 100644
--- a/packages/horn-core/tests/test_geometry_designer.py
+++ b/packages/horn-core/tests/test_geometry_designer.py
@@ -12,6 +12,7 @@
     derive_mouth_radius,
     derive_mouth_radius_range,
     derive_simulation_freq_range,
+    generate_auto_candidates,
     generate_fullauto_candidates,
 )
 
@@ -236,3 +237,91 @@ def test_high_freq_produces_small_horn(self):
         )
         assert derived.ideal_mouth_radius == pytest.approx(0.01365, rel=1e-2)
         assert derived.length_range[1] < 0.05  # half-wave at 4000 Hz ~0.043m
+
+
+class TestGenerateAutoCandidates:
+    def test_all_free(self):
+        """With no fixed params, behaves like fullauto: 7×1×3×3 = 63."""
+        candidates, derived = generate_auto_candidates(
+            target_f_low=500,
+            target_f_high=4000,
+            throat_radii=[0.025],
+        )
+        assert len(candidates) == 63
+        assert derived.candidate_count == 63
+
+    def test_fixed_mouth_radius(self):
+        """Fixed mouth radius → 7×1×1×3 = 21."""
+        candidates, derived = generate_auto_candidates(
+            target_f_low=500,
+            target_f_high=4000,
+            throat_radii=[0.025],
+            mouth_radius=0.15,
+        )
+        assert len(candidates) == 21
+        # All candidates should have the fixed mouth radius
+        for c in candidates:
+            assert c.mouth_radius == pytest.approx(0.15)
+        # Derived range should be fixed
+        assert derived.mouth_radius_range[0] == derived.mouth_radius_range[1] == 0.15
+
+    def test_fixed_length(self):
+        """Fixed length → 7×1×3×1 = 21."""
+        candidates, derived = generate_auto_candidates(
+            target_f_low=500,
+            target_f_high=4000,
+            throat_radii=[0.025],
+            length=0.3,
+        )
+        assert len(candidates) == 21
+        for c in candidates:
+            assert c.length == pytest.approx(0.3)
+        assert derived.length_range[0] == derived.length_range[1] == 0.3
+
+    def test_both_fixed(self):
+        """Fixed mouth + length → 7×1×1×1 = 7 (like old auto)."""
+        candidates, derived = generate_auto_candidates(
+            target_f_low=500,
+            target_f_high=4000,
+            throat_radii=[0.025],
+            mouth_radius=0.15,
+            length=0.3,
+        )
+        assert len(candidates) == 7
+        profiles = {c.profile for c in candidates}
+        assert profiles == {"conical", "exponential", "hyperbolic", "tractrix", "os", "lecleach", "cd"}
+
+    def test_multiple_throat_radii_fixed_geom(self):
+        """3 throat radii with fixed mouth+length → 7×3×1×1 = 21."""
+        candidates, _ = generate_auto_candidates(
+            target_f_low=500,
+            target_f_high=4000,
+            throat_radii=[0.015, 0.020, 0.025],
+            mouth_radius=0.15,
+            length=0.3,
+        )
+        assert len(candidates) == 21
+
+    def test_candidate_ids_prefixed_auto(self):
+        """Candidate IDs should start with 'auto_'."""
+        candidates, _ = generate_auto_candidates(
+            target_f_low=500,
+            target_f_high=4000,
+            throat_radii=[0.025],
+            mouth_radius=0.15,
+            length=0.3,
+        )
+        for c in candidates:
+            assert c.candidate_id.startswith("auto_")
+
+    def test_sim_freq_range_always_set(self):
+        """Simulation freq range should be set regardless of fixed params."""
+        _, derived = generate_auto_candidates(
+            target_f_low=500,
+            target_f_high=4000,
+            throat_radii=[0.025],
+            mouth_radius=0.15,
+            length=0.3,
+        )
+        assert derived.sim_freq_range[0] < 500
+        assert derived.sim_freq_range[1] > 4000
diff --git a/packages/horn-core/tests/test_profiles.py b/packages/horn-core/tests/test_profiles.py
new file mode 100644
index 0000000..bd3fed6
--- /dev/null
+++ b/packages/horn-core/tests/test_profiles.py
@@ -0,0 +1,149 @@
+"""Tests for pure-Python horn profile radius functions."""
+
+import math
+
+import numpy as np
+import pytest
+
+from horn_core.profiles import (
+    conical_radius,
+    exponential_radius,
+    hyperbolic_radius,
+    tractrix_radius,
+    os_radius,
+    lecleach_radius,
+    cd_radius,
+    get_radius_func,
+)
+
+
+THROAT_R = 0.025
+MOUTH_R = 0.15
+LENGTH = 0.3
+ALL_PROFILES = ["conical", "exponential", "hyperbolic", "tractrix", "os", "lecleach", "cd"]
+
+
+class TestEndpoints:
+    """Each profile must satisfy r(0) == throat_radius and r(L) == mouth_radius."""
+
+    @pytest.fixture(params=ALL_PROFILES)
+    def profile(self, request):
+        return request.param
+
+    def test_throat_radius(self, profile):
+        func = get_radius_func(profile, THROAT_R, MOUTH_R, LENGTH)
+        # Le Cléac'h clips the tractrix curve, so throat radius is approximate
+        tol = 0.02 if profile == "lecleach" else 1e-3
+        assert func(0.0) == pytest.approx(THROAT_R, rel=tol)
+
+    def test_mouth_radius(self, profile):
+        func = get_radius_func(profile, THROAT_R, MOUTH_R, LENGTH)
+        tol = 0.02 if profile == "lecleach" else 1e-3
+        assert func(LENGTH) == pytest.approx(MOUTH_R, rel=tol)
+
+
+class TestMonotonicity:
+    """Horn profiles should expand monotonically (r(z2) >= r(z1) for z2 > z1)."""
+
+    @pytest.fixture(params=ALL_PROFILES)
+    def profile(self, request):
+        return request.param
+
+    def test_monotonic_expansion(self, profile):
+        func = get_radius_func(profile, THROAT_R, MOUTH_R, LENGTH)
+        z_vals = np.linspace(0, LENGTH, 200)
+        r_vals = [func(z) for z in z_vals]
+        for i in range(1, len(r_vals)):
+            assert r_vals[i] >= r_vals[i - 1] - 1e-10, (
+                f"{profile}: non-monotonic at z={z_vals[i]:.4f}, "
+                f"r[{i-1}]={r_vals[i-1]:.6f} > r[{i}]={r_vals[i]:.6f}"
+            )
+
+
+class TestConicalAnalytical:
+    """Conical profile has exact analytical form."""
+
+    def test_midpoint(self):
+        z_mid = LENGTH / 2
+        expected = THROAT_R + (MOUTH_R - THROAT_R) / 2
+        assert conical_radius(z_mid, THROAT_R, MOUTH_R, LENGTH) == pytest.approx(expected)
+
+    def test_quarter_point(self):
+        z = LENGTH / 4
+        expected = THROAT_R + (MOUTH_R - THROAT_R) / 4
+        assert conical_radius(z, THROAT_R, MOUTH_R, LENGTH) == pytest.approx(expected)
+
+
+class TestExponentialAnalytical:
+    """Exponential profile: r(z) = r_t * exp(m*z/L)."""
+
+    def test_midpoint(self):
+        m = math.log(MOUTH_R / THROAT_R)
+        z_mid = LENGTH / 2
+        expected = THROAT_R * math.exp(m * 0.5)
+        assert exponential_radius(z_mid, THROAT_R, MOUTH_R, LENGTH) == pytest.approx(expected)
+
+
+class TestHyperbolicAnalytical:
+    """Hyperbolic profile: r(z) = r_t * cosh(m*z/L)."""
+
+    def test_midpoint(self):
+        m = np.arccosh(MOUTH_R / THROAT_R)
+        z_mid = LENGTH / 2
+        expected = float(THROAT_R * np.cosh(m * 0.5))
+        assert hyperbolic_radius(z_mid, THROAT_R, MOUTH_R, LENGTH) == pytest.approx(expected, rel=1e-6)
+
+
+class TestGetRadiusFunc:
+    def test_invalid_profile_raises(self):
+        with pytest.raises(ValueError, match="Unknown profile"):
+            get_radius_func("parabolic", THROAT_R, MOUTH_R, LENGTH)
+
+    def test_closure_returns_float(self):
+        for profile in ALL_PROFILES:
+            func = get_radius_func(profile, THROAT_R, MOUTH_R, LENGTH)
+            result = func(LENGTH / 2)
+            assert isinstance(result, float)
+
+    def test_different_geometries(self):
+        """Different throat/mouth should give different mid-section radii."""
+        f1 = get_radius_func("exponential", 0.01, 0.1, 0.3)
+        f2 = get_radius_func("exponential", 0.01, 0.2, 0.3)
+        assert f1(0.15) < f2(0.15)
+
+
+class TestCrossValidation:
+    """Verify profiles.py matches horn-geometry/generator.py radius functions."""
+
+    def test_conical_matches_generator(self):
+        """Conical is trivial — exact match expected."""
+        from horn_core.profiles import conical_radius
+        for z in np.linspace(0, LENGTH, 50):
+            expected = THROAT_R + (MOUTH_R - THROAT_R) * z / LENGTH
+            assert conical_radius(z, THROAT_R, MOUTH_R, LENGTH) == pytest.approx(expected, abs=1e-12)
+
+    def test_exponential_matches_generator(self):
+        from horn_core.profiles import exponential_radius
+        m = np.log(MOUTH_R / THROAT_R)
+        for z in np.linspace(0, LENGTH, 50):
+            expected = THROAT_R * np.exp(m * z / LENGTH)
+            assert exponential_radius(z, THROAT_R, MOUTH_R, LENGTH) == pytest.approx(float(expected), rel=1e-10)
+
+    def test_os_matches_generator(self):
+        from horn_core.profiles import os_radius
+        theta = math.atan2(math.sqrt(MOUTH_R**2 - THROAT_R**2), LENGTH)
+        for z in np.linspace(0, LENGTH, 50):
+            expected = math.sqrt(THROAT_R**2 + (z * math.tan(theta))**2)
+            assert os_radius(z, THROAT_R, MOUTH_R, LENGTH) == pytest.approx(expected, rel=1e-10)
+
+    def test_cd_matches_generator(self):
+        from horn_core.profiles import cd_radius
+        frac = 0.3
+        z_t = frac * LENGTH
+        r_trans = THROAT_R * (MOUTH_R / THROAT_R) ** frac
+        for z in np.linspace(0, LENGTH, 50):
+            if z <= z_t:
+                expected = THROAT_R * math.exp(math.log(r_trans / THROAT_R) * z / z_t)
+            else:
+                expected = r_trans + (MOUTH_R - r_trans) * (z - z_t) / (LENGTH - z_t)
+            assert cd_radius(z, THROAT_R, MOUTH_R, LENGTH) == pytest.approx(expected, rel=1e-10)
diff --git a/packages/horn-core/tests/test_webster.py b/packages/horn-core/tests/test_webster.py
new file mode 100644
index 0000000..2579a8e
--- /dev/null
+++ b/packages/horn-core/tests/test_webster.py
@@ -0,0 +1,153 @@
+"""Tests for the TMM/Webster horn throat impedance solver."""
+
+import math
+
+import numpy as np
+import pytest
+
+from horn_core.profiles import get_radius_func
+from horn_core.webster import (
+    C0,
+    RHO0,
+    compute_throat_impedance_tmm,
+    piston_radiation_impedance,
+    unflanged_radiation_impedance,
+)
+
+
+class TestRadiationImpedance:
+    """Verify radiation impedance matches horn_solver/radiation.py."""
+
+    def test_flanged_piston_small_ka(self):
+        """For small ka, R1 ≈ (ka)²/2, X1 ≈ 8ka/(3π)."""
+        k = 2 * math.pi * 100 / C0  # 100 Hz
+        a = 0.01  # small radius
+        z = piston_radiation_impedance(k, a)
+        ka = k * a
+        # Taylor: R1(2ka) ≈ (2ka)²/8 = ka²/2
+        assert z.real == pytest.approx(ka**2 / 2, rel=0.1)
+        assert z.imag > 0  # reactance is positive
+
+    def test_flanged_piston_large_ka(self):
+        """For large ka, R1 → 1."""
+        k = 2 * math.pi * 10000 / C0
+        a = 0.1  # large radius
+        z = piston_radiation_impedance(k, a)
+        assert z.real == pytest.approx(1.0, abs=0.1)
+
+    def test_unflanged_small_ka(self):
+        k = 2 * math.pi * 100 / C0
+        a = 0.01
+        z = unflanged_radiation_impedance(k, a)
+        ka = k * a
+        assert z.real == pytest.approx(0.25 * ka**2, rel=1e-6)
+        assert z.imag == pytest.approx(0.6133 * ka, rel=1e-6)
+
+
+class TestCylindricalPipe:
+    """Cylindrical horn (r_t == r_m) should behave like a simple pipe."""
+
+    def test_real_impedance_positive(self):
+        """Re(Z) > 0 at all frequencies (energy conservation)."""
+        r = 0.025
+        L = 0.3
+        freq = np.linspace(200, 8000, 50)
+        func = get_radius_func("conical", r, r + 1e-10, L)  # near-cylindrical
+        z_real, z_imag = compute_throat_impedance_tmm(
+            freq, func, L, r, r + 1e-10, n_segments=100
+        )
+        # Real part should be non-negative (radiation resistance at mouth)
+        assert np.all(z_real >= -1e-6), f"Negative Re(Z) found: min={z_real.min()}"
+
+
+class TestConicalHorn:
+    """TMM for a conical horn — compare against known properties."""
+
+    THROAT_R = 0.025
+    MOUTH_R = 0.15
+    LENGTH = 0.3
+
+    def test_real_impedance_positive(self):
+        """Re(Z_throat) > 0 at all frequencies."""
+        freq = np.linspace(300, 8000, 80)
+        func = get_radius_func("conical", self.THROAT_R, self.MOUTH_R, self.LENGTH)
+        z_real, z_imag = compute_throat_impedance_tmm(
+            freq, func, self.LENGTH, self.THROAT_R, self.MOUTH_R
+        )
+        assert np.all(z_real > -1e-6)
+
+    def test_impedance_magnitude_bounded(self):
+        """Impedance magnitude should be bounded by ρc (= 420 Pa·s/m)."""
+        freq = np.linspace(500, 8000, 80)
+        func = get_radius_func("conical", self.THROAT_R, self.MOUTH_R, self.LENGTH)
+        z_real, z_imag = compute_throat_impedance_tmm(
+            freq, func, self.LENGTH, self.THROAT_R, self.MOUTH_R
+        )
+        z_mag = np.sqrt(z_real**2 + z_imag**2)
+        # Specific impedance should be in a reasonable range
+        # (not zero, not wildly larger than rho*c)
+        rho_c = RHO0 * C0
+        assert np.all(z_mag < 10 * rho_c), f"Impedance too large: max={z_mag.max()}"
+
+    def test_high_frequency_approaches_rho_c(self):
+        """At high frequencies where ka >> 1, throat impedance → ρc."""
+        freq = np.array([15000.0, 20000.0])
+        func = get_radius_func("conical", self.THROAT_R, self.MOUTH_R, self.LENGTH)
+        z_real, z_imag = compute_throat_impedance_tmm(
+            freq, func, self.LENGTH, self.THROAT_R, self.MOUTH_R, n_segments=400
+        )
+        rho_c = RHO0 * C0
+        # At high freq the horn looks like a matched load
+        # Real part should approach ρc within a factor of ~2
+        for i in range(len(freq)):
+            assert 0.1 * rho_c < z_real[i] < 5 * rho_c
+
+
+class TestSegmentConvergence:
+    """TMM should converge as n_segments increases."""
+
+    def test_convergence(self):
+        freq = np.linspace(500, 4000, 30)
+        func = get_radius_func("exponential", 0.025, 0.15, 0.3)
+
+        z_50 = compute_throat_impedance_tmm(freq, func, 0.3, 0.025, 0.15, n_segments=50)
+        z_200 = compute_throat_impedance_tmm(freq, func, 0.3, 0.025, 0.15, n_segments=200)
+        z_500 = compute_throat_impedance_tmm(freq, func, 0.3, 0.025, 0.15, n_segments=500)
+
+        # Error between 200 and 500 should be smaller than between 50 and 200
+        err_50_200 = np.mean(np.abs(z_50[0] - z_200[0]))
+        err_200_500 = np.mean(np.abs(z_200[0] - z_500[0]))
+        assert err_200_500 < err_50_200, (
+            f"Not converging: err(50,200)={err_50_200:.4f}, err(200,500)={err_200_500:.4f}"
+        )
+
+
+class TestAllProfiles:
+    """Basic sanity checks across all 7 profiles."""
+
+    PROFILES = ["conical", "exponential", "hyperbolic", "tractrix", "os", "lecleach", "cd"]
+
+    @pytest.fixture(params=PROFILES)
+    def profile(self, request):
+        return request.param
+
+    def test_produces_valid_output(self, profile):
+        """TMM returns finite, correctly shaped arrays for all profiles."""
+        freq = np.linspace(500, 4000, 20)
+        func = get_radius_func(profile, 0.025, 0.15, 0.3)
+        z_real, z_imag = compute_throat_impedance_tmm(
+            freq, func, 0.3, 0.025, 0.15, n_segments=100
+        )
+        assert z_real.shape == freq.shape
+        assert z_imag.shape == freq.shape
+        assert np.all(np.isfinite(z_real))
+        assert np.all(np.isfinite(z_imag))
+
+    def test_real_part_positive(self, profile):
+        """Re(Z) should be positive (passive system)."""
+        freq = np.linspace(500, 4000, 30)
+        func = get_radius_func(profile, 0.025, 0.15, 0.3)
+        z_real, _ = compute_throat_impedance_tmm(
+            freq, func, 0.3, 0.025, 0.15, n_segments=200
+        )
+        assert np.all(z_real > -1e-3), f"{profile}: Re(Z) negative, min={z_real.min()}"
diff --git a/packages/horn-drivers/src/horn_drivers/loader.py b/packages/horn-drivers/src/horn_drivers/loader.py
index 2ee819e..cd1b0a0 100644
--- a/packages/horn-drivers/src/horn_drivers/loader.py
+++ b/packages/horn-drivers/src/horn_drivers/loader.py
@@ -34,6 +34,7 @@ def _driver_from_dict(d: dict) -> DriverParameters:
     le_h = params.get("le_h") or _mh_to_h(params.get("le_mh"))
     sd_m2 = params.get("sd_m2") or params.get("sd_sq_meters")
     xmax_m = params.get("xmax_m") or _mm_to_m(params.get("xmax_mm"))
+    exit_area_m2 = params.get("exit_area_m2") or _cm2_to_m2(params.get("exit_area_cm2"))
 
     return DriverParameters(
         driver_id=d.get("driver_id", "unknown"),
@@ -48,10 +49,13 @@ def _driver_from_dict(d: dict) -> DriverParameters:
         qms=params.get("qms"),
         qes=params.get("qes"),
         qts=params.get("qts"),
+        exit_area_m2=exit_area_m2,
         driver_type=d.get("driver_type"),
         nominal_diameter=d.get("nominal_diameter"),
         xmax_m=xmax_m,
         nominal_impedance_ohm=params.get("nominal_impedance_ohm"),
+        power_w=params.get("power_w"),
+        peak_power_w=params.get("peak_power_w"),
     )
 
 
@@ -125,3 +129,7 @@ def _mh_to_h(val: Optional[float]) -> Optional[float]:
 
 def _mm_to_m(val: Optional[float]) -> Optional[float]:
     return val * 1e-3 if val is not None else None
+
+
+def _cm2_to_m2(val: Optional[float]) -> Optional[float]:
+    return val * 1e-4 if val is not None else None
diff --git a/packages/horn-drivers/src/horn_drivers/scraper.py b/packages/horn-drivers/src/horn_drivers/scraper.py
index f39445b..72c5977 100644
--- a/packages/horn-drivers/src/horn_drivers/scraper.py
+++ b/packages/horn-drivers/src/horn_drivers/scraper.py
@@ -127,6 +127,10 @@ def _parse_data_woofer(json_str: str) -> Optional[dict]:
         si["nominal_impedance_ohm"] = float(raw["z"])
     if raw.get("xmax") is not None and raw["xmax"] > 0:
         si["xmax_m"] = float(raw["xmax"]) * 1e-3  # mm -> m
+    if raw.get("pmax") is not None and raw["pmax"] > 0:
+        pmax = float(raw["pmax"])
+        si["peak_power_w"] = pmax             # pmax is program power
+        si["power_w"] = pmax / 2.0            # RMS ≈ program / 2
 
     return si
 
diff --git a/packages/horn-geometry/src/horn_geometry/generator.py b/packages/horn-geometry/src/horn_geometry/generator.py
index 5a56a69..5373e18 100644
--- a/packages/horn-geometry/src/horn_geometry/generator.py
+++ b/packages/horn-geometry/src/horn_geometry/generator.py
@@ -7,6 +7,45 @@
 import argparse
 
 
+def _adaptive_z_positions(
+    radius_func: Callable[[float], float],
+    length: float,
+    num_sections: int,
+) -> np.ndarray:
+    """Compute z-positions that concentrate sections where the radius changes fast.
+
+    Uses a blended metric: 50 % uniform arc-length + 50 % radius-change.  This
+    ensures adequate section density everywhere (no huge gaps in flat regions)
+    while still packing extra sections into steep regions (e.g. the mouth of a
+    tractrix).  Degenerates to uniform spacing for constant-gradient profiles.
+    """
+    n_samples = 5000
+    z_fine = np.linspace(0, length, n_samples)
+    r_fine = np.array([radius_func(z) for z in z_fine])
+
+    dr = np.abs(np.diff(r_fine))
+    dz = np.diff(z_fine)
+
+    # Blend: uniform component (dz) + curvature component (dr)
+    total_dr = dr.sum()
+    total_dz = dz.sum()
+
+    if total_dr == 0:
+        return np.linspace(0, length, num_sections)
+
+    # Normalise both to the same scale, then blend 50/50
+    metric = 0.5 * (dz / total_dz) + 0.5 * (dr / total_dr)
+    cum_metric = np.concatenate(([0.0], np.cumsum(metric)))
+
+    targets = np.linspace(0, cum_metric[-1], num_sections)
+    z_positions = np.interp(targets, cum_metric, z_fine)
+
+    # Guarantee exact endpoints
+    z_positions[0] = 0.0
+    z_positions[-1] = length
+    return z_positions
+
+
 def _loft_horn_profile(
     radius_func: Callable[[float], float],
     length: float,
@@ -30,7 +69,7 @@ def _loft_horn_profile(
     gmsh.model.add(model_name)
     gmsh.option.setNumber("General.Terminal", 1)
 
-    z_positions = np.linspace(0, length, num_sections)
+    z_positions = _adaptive_z_positions(radius_func, length, num_sections)
     curve_loops = []
 
     for z in z_positions:
diff --git a/packages/horn-geometry/tests/test_profiles.py b/packages/horn-geometry/tests/test_profiles.py
index 4beb0fb..d1c3c2f 100644
--- a/packages/horn-geometry/tests/test_profiles.py
+++ b/packages/horn-geometry/tests/test_profiles.py
@@ -244,6 +244,41 @@ def test_throat_mouth_radii_match(self, tmp_path):
             )
 
 
+@pytest.mark.skipif(gmsh is None, reason="gmsh not available")
+class TestProfileMeshable:
+    """Verify that every profile produces a geometry that gmsh can actually mesh."""
+
+    THROAT_R = 0.025
+    MOUTH_R = 0.15
+    LENGTH = 0.3
+
+    @pytest.fixture(params=["conical", "exponential", "hyperbolic", "tractrix", "os", "lecleach", "cd"])
+    def profile_name(self, request):
+        return request.param
+
+    def test_all_profiles_meshable(self, profile_name, tmp_path):
+        step = create_horn(
+            profile=profile_name,
+            throat_radius=self.THROAT_R,
+            mouth_radius=self.MOUTH_R,
+            length=self.LENGTH,
+            output_file=tmp_path / f"{profile_name}.step",
+            num_sections=20,
+        )
+        gmsh.initialize()
+        gmsh.option.setNumber("General.Terminal", 0)
+        gmsh.model.add("mesh_check")
+        gmsh.model.occ.importShapes(str(step))
+        gmsh.model.occ.synchronize()
+        gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.01)
+        try:
+            gmsh.model.mesh.generate(3)
+            nodes = gmsh.model.mesh.getNodes()
+            assert len(nodes[0]) > 0, f"{profile_name}: mesh produced no nodes"
+        finally:
+            gmsh.finalize()
+
+
 @pytest.mark.skipif(gmsh is None, reason="gmsh not available")
 class TestCreateHornDispatch:
     """Test the create_horn dispatch function."""
diff --git a/packages/horn-solver/Dockerfile b/packages/horn-solver/Dockerfile
index 4072a10..c7cf4e8 100644
--- a/packages/horn-solver/Dockerfile
+++ b/packages/horn-solver/Dockerfile
@@ -18,7 +18,7 @@ COPY ./packages /app/packages
 RUN pip install --no-deps /app/packages/horn-core \
     /app/packages/horn-drivers \
     && pip install /app/packages/horn-solver \
-    && pip install bempp-cl
+    && pip install 'bempp-cl==0.3.1'
 
 # --- Test Stage ---
 FROM production as test
diff --git a/packages/horn-solver/src/horn_solver/bem_coupling.py b/packages/horn-solver/src/horn_solver/bem_coupling.py
index 5cab566..7123734 100644
--- a/packages/horn-solver/src/horn_solver/bem_coupling.py
+++ b/packages/horn-solver/src/horn_solver/bem_coupling.py
@@ -23,6 +23,27 @@
     from bempp.api.external import fenicsx as bempp_fenicsx
     from bempp.api.assembly.blocked_operator import BlockedDiscreteOperator
 
+    # Patch bempp's fenics_space_info for dolfinx 0.8 compatibility.
+    # bempp-cl 0.3.x calls element.family() which doesn't exist on
+    # dolfinx 0.8's _BasixElement; the family lives on basix_element instead.
+    _orig_space_info = bempp_fenicsx.fenics_space_info
+
+    def _patched_space_info(fenics_space):
+        element = fenics_space.ufl_element()
+        if hasattr(element, "family"):
+            return _orig_space_info(fenics_space)
+        # dolfinx 0.8: pull family name from basix_element
+        import basix
+        _family_map = {
+            basix.ElementFamily.P: "Lagrange",
+        }
+        be = element.basix_element
+        family = _family_map.get(be.family, str(be.family))
+        degree = be.degree
+        return (family, degree)
+
+    bempp_fenicsx.fenics_space_info = _patched_space_info
+
     BEMPP_AVAILABLE = True
 except ImportError:
     BEMPP_AVAILABLE = False
diff --git a/packages/horn-solver/src/horn_solver/solver.py b/packages/horn-solver/src/horn_solver/solver.py
index 231e666..7692db9 100644
--- a/packages/horn-solver/src/horn_solver/solver.py
+++ b/packages/horn-solver/src/horn_solver/solver.py
@@ -124,7 +124,19 @@ def create_mesh_from_step(step_file: str, mesh_size: float, horn_length: float)
     gmsh.option.setNumber("Mesh.MeshSizeMin", mesh_size)
     gmsh.option.setNumber("Mesh.MeshSizeMax", mesh_size)
     gmsh.model.mesh.generate(3)
-    
+
+    # Check mesh size before converting to dolfinx (avoid OOM)
+    node_tags, _, _ = gmsh.model.mesh.getNodes()
+    n_nodes = len(node_tags)
+    max_nodes = int(os.environ.get("HORN_MAX_MESH_NODES", "100000"))
+    if n_nodes > max_nodes:
+        gmsh.finalize()
+        raise RuntimeError(
+            f"Mesh too large: {n_nodes} nodes exceeds limit of {max_nodes}. "
+            f"Skipping this candidate to avoid OOM."
+        )
+    print(f"Mesh nodes: {n_nodes} (limit: {max_nodes})")
+
     # Note the change here: we are now capturing all three return values.
     domain, cell_tags, facet_tags = gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=3)
     gmsh.finalize()
@@ -357,8 +369,12 @@ def run_simulation(
                     k,
                     ff_directions,
                 )
-                p_ref = 20e-6
-                spl_far = 20 * np.log10(np.abs(p_far) / p_ref + 1e-12)
+                # Far-field operators return pattern values with arbitrary
+                # scaling.  Normalise so that the on-axis (theta=0) magnitude
+                # is 0 dB, giving a relative radiation pattern in dB.
+                p_far_abs = np.abs(p_far)
+                p_on_axis = p_far_abs[0] if p_far_abs[0] > 0 else 1e-30
+                spl_far = 20 * np.log10(p_far_abs / p_on_axis + 1e-30)
                 for angle_deg, spl_val in zip(directivity_angles, spl_far):
                     directivity_results.append({
                         "frequency": frequency,