from __future__ import annotations

import os
import time
from dataclasses import dataclass

import numpy as np

import cfsem

GRID_N = 124 if os.getenv("CFSEM_TESTING") else 228
DEFAULT_SOURCE_COUNT = 63 if os.getenv("CFSEM_TESTING") else 159
MIN_SOURCE_COUNT = 8
MAX_SOURCE_COUNT = 255 if os.getenv("CFSEM_TESTING") else 100_000
MAX_DIRECT_SELF_DIPOLE_SOURCES = 100_000
MAX_DIRECT_SELF_FILAMENT_SOURCES = 100_000
MAX_DIRECT_SELF_BOUNDARY_SOURCES = 10_000
CURRENT = 1.0
DEFAULT_SOURCE_RADIUS = 0.02
HELICAL_WIRE_RADIUS = 0.055
TRIANGLE_STRIP_WIDTH = 0.08
SOURCE_SPAN = 1.4
DISTRIBUTED_SOURCE_RADIUS = 0.5 * SOURCE_SPAN
DISTRIBUTED_SOURCE_CENTER_OFFSET = 2.5 * SOURCE_SPAN
DEFAULT_TWIST_PITCH = 0.36
DEFAULT_HELIX_WIDTH = HELICAL_WIRE_RADIUS
DEFAULT_BEND_CURVATURE = 2.0 / SOURCE_SPAN
DEFAULT_LOOP_FRACTION = 0.5
DEFAULT_THETA = 0.05
LOG10_MIN_SOURCE_COUNT = float(np.log10(MIN_SOURCE_COUNT))
LOG10_DEFAULT_SOURCE_COUNT = float(np.log10(DEFAULT_SOURCE_COUNT))
LOG10_MAX_SOURCE_COUNT = float(np.log10(MAX_SOURCE_COUNT))
MAX_PLOTTED_PATH_POINTS = 400
MAX_PLOTTED_AABBS = 500
DISTRIBUTED_SOURCE_SEED = 1729


def contiguous_triple(
    x: np.ndarray,
    y: np.ndarray,
    z: np.ndarray,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Return x/y/z component arrays that satisfy the Rust bindings' slice requirements."""

    return (
        np.ascontiguousarray(x),
        np.ascontiguousarray(y),
        np.ascontiguousarray(z),
    )


@dataclass(frozen=True)
class Geometry:
    centerline: np.ndarray
    helix: np.ndarray
    xyzfil: tuple[np.ndarray, np.ndarray, np.ndarray]
    dlxyzfil: tuple[np.ndarray, np.ndarray, np.ndarray]
    current: np.ndarray
    wire_radius: np.ndarray
    volume_xyzfil: tuple[np.ndarray, np.ndarray, np.ndarray]
    volume_dlxyzfil: tuple[np.ndarray, np.ndarray, np.ndarray]
    volume_current: np.ndarray
    volume_wire_radius: np.ndarray
    dipole_loc: tuple[np.ndarray, np.ndarray, np.ndarray]
    dipole_moment: tuple[np.ndarray, np.ndarray, np.ndarray]
    dipole_outer_radius: np.ndarray
    volume_dipole_loc: tuple[np.ndarray, np.ndarray, np.ndarray]
    volume_dipole_moment: tuple[np.ndarray, np.ndarray, np.ndarray]
    volume_dipole_outer_radius: np.ndarray
    strip_nodes: np.ndarray
    strip_triangles: np.ndarray
    strip_stream_function: np.ndarray
    volume_strip_nodes: np.ndarray
    volume_strip_triangles: np.ndarray
    volume_strip_stream_function: np.ndarray
    obs: tuple[np.ndarray, np.ndarray, np.ndarray]
    obs_grid: tuple[np.ndarray, np.ndarray]
    extent: float


def fixed_span_arc_centerline(span: float, curvature: float, loop_fraction: float, n: int) -> np.ndarray:
    curvature = max(0.0, float(curvature))
    x = np.linspace(-0.5 * span, 0.5 * span, n, endpoint=True)
    loop_fraction = max(0.0, min(1.0, float(loop_fraction)))
    if curvature <= 1.0e-12 or loop_fraction <= 1.0e-12:
        return np.vstack((x, np.zeros_like(x), np.zeros_like(x)))

    max_curvature = 2.0 / span
    curvature = min(curvature, max_curvature)
    radius = 1.0 / curvature
    half_angle = np.arcsin(0.5 * span / radius)
    half_angle = min(np.pi, half_angle * loop_fraction / DEFAULT_LOOP_FRACTION)
    if half_angle <= 1.0e-12:
        return np.vstack((x, np.zeros_like(x), np.zeros_like(x)))
    theta = np.linspace(-half_angle, half_angle, n, endpoint=True)
    x = radius * np.sin(theta)
    z = radius * np.cos(theta)
    return np.vstack((x, np.zeros_like(x), z))


def section_observation_plane(
    extent: float,
    n: int,
) -> tuple[tuple[np.ndarray, np.ndarray, np.ndarray], tuple[np.ndarray, np.ndarray]]:
    x = np.linspace(-extent, extent, n)
    z = np.linspace(-extent, extent, n)
    xg, zg = np.meshgrid(x, z, indexing="xy")
    yg = np.zeros_like(xg)
    return contiguous_triple(xg.ravel(), yg.ravel(), zg.ravel()), (xg, zg)


def build_geometry(
    geometry_layout: str,
    n_centerline: int,
    grid_n: int,
    twist_pitch: float,
    helix_width: float,
    bend_curvature: float,
    loop_fraction: float,
    source_radius: float,
) -> Geometry:
    source_radius = max(0.0, float(source_radius))
    centerline = fixed_span_arc_centerline(SOURCE_SPAN, bend_curvature, loop_fraction, n_centerline)
    helix = np.asarray(
        cfsem.filament_helix_path(
            path=centerline,
            helix_start_offset=(0.0, helix_width, 0.0),
            twist_pitch=float(twist_pitch),
            angle_offset=0.0,
        )
    )
    segment_centers = 0.5 * (helix[:, :-1] + helix[:, 1:])
    source_centroid = np.mean(segment_centers, axis=1, keepdims=True)
    centerline = centerline - source_centroid
    helix = helix - source_centroid
    starts = helix[:, :-1].T
    ends = helix[:, 1:].T
    dl = ends - starts
    xyzfil = contiguous_triple(starts[:, 0], starts[:, 1], starts[:, 2])
    dlxyzfil = contiguous_triple(dl[:, 0], dl[:, 1], dl[:, 2])
    current = np.full(starts.shape[0], CURRENT)
    wire_radius = np.full(starts.shape[0], source_radius)
    dipole_loc, dipole_moment, dipole_outer_radius = build_segment_dipoles(
        starts,
        dl,
        current,
        source_radius,
    )
    (
        volume_dipole_loc,
        volume_dipole_moment,
        volume_dipole_outer_radius,
        volume_xyzfil,
        volume_dlxyzfil,
        volume_current,
        volume_wire_radius,
        volume_strip_nodes,
        volume_strip_triangles,
        volume_strip_stream_function,
    ) = build_volume_sources(
        max(1, starts.shape[0]),
        DISTRIBUTED_SOURCE_RADIUS,
        source_radius,
    )
    strip_nodes, strip_triangles, strip_stream_function = build_triangle_strip(
        helix,
        TRIANGLE_STRIP_WIDTH,
        CURRENT,
    )
    extent = 8.0 * (SOURCE_SPAN + 3.0 * helix_width)
    obs, obs_grid = section_observation_plane(extent, grid_n)
    return Geometry(
        centerline,
        helix,
        xyzfil,
        dlxyzfil,
        current,
        wire_radius,
        volume_xyzfil,
        volume_dlxyzfil,
        volume_current,
        volume_wire_radius,
        dipole_loc,
        dipole_moment,
        dipole_outer_radius,
        volume_dipole_loc,
        volume_dipole_moment,
        volume_dipole_outer_radius,
        strip_nodes,
        strip_triangles,
        strip_stream_function,
        volume_strip_nodes,
        volume_strip_triangles,
        volume_strip_stream_function,
        obs,
        obs_grid,
        extent,
    )


def build_volume_sources(
    n: int,
    radius: float,
    source_radius: float,
) -> tuple[
    tuple[np.ndarray, np.ndarray, np.ndarray],
    tuple[np.ndarray, np.ndarray, np.ndarray],
    np.ndarray,
    tuple[np.ndarray, np.ndarray, np.ndarray],
    tuple[np.ndarray, np.ndarray, np.ndarray],
    np.ndarray,
    np.ndarray,
    np.ndarray,
    np.ndarray,
    np.ndarray,
]:
    rng = np.random.default_rng(DISTRIBUTED_SOURCE_SEED)
    angles = 2.0 * np.pi * rng.random(n)
    radii = radius * np.sqrt(rng.random(n))
    disk_offsets = np.array(
        [
            [-DISTRIBUTED_SOURCE_CENTER_OFFSET, 0.0, -DISTRIBUTED_SOURCE_CENTER_OFFSET],
            [DISTRIBUTED_SOURCE_CENTER_OFFSET, 0.0, -DISTRIBUTED_SOURCE_CENTER_OFFSET],
            [-DISTRIBUTED_SOURCE_CENTER_OFFSET, 0.0, DISTRIBUTED_SOURCE_CENTER_OFFSET],
            [DISTRIBUTED_SOURCE_CENTER_OFFSET, 0.0, DISTRIBUTED_SOURCE_CENTER_OFFSET],
        ]
    )
    disk_ids = np.arange(n) % disk_offsets.shape[0]
    loc = np.column_stack((radii * np.cos(angles), np.zeros(n), radii * np.sin(angles)))
    loc = loc + disk_offsets[disk_ids]

    central_moment = np.array([0.0, 0.0, 1.0])
    r_norm = np.linalg.norm(loc, axis=1)
    r_hat = loc / np.maximum(r_norm[:, None], 1.0e-30)
    m_dot_r = r_hat @ central_moment
    direction = 3.0 * r_hat * m_dot_r[:, None] - central_moment[None, :]
    direction_norm = np.linalg.norm(direction, axis=1)
    direction = direction / np.maximum(direction_norm[:, None], 1.0e-30)
    moment_scale = CURRENT * radius**3 / max(1, n)
    moments = moment_scale * direction
    segment_length = 2.0 * radius / max(4.0, n ** (1.0 / 3.0))
    starts = loc - 0.5 * segment_length * direction
    dl = segment_length * direction
    strip_nodes, strip_triangles, strip_stream_function = build_distributed_triangle_strips(
        loc,
        direction,
        segment_length,
        TRIANGLE_STRIP_WIDTH,
        CURRENT,
    )
    return (
        contiguous_triple(loc[:, 0], loc[:, 1], loc[:, 2]),
        contiguous_triple(moments[:, 0], moments[:, 1], moments[:, 2]),
        np.full(n, source_radius),
        contiguous_triple(starts[:, 0], starts[:, 1], starts[:, 2]),
        contiguous_triple(dl[:, 0], dl[:, 1], dl[:, 2]),
        np.full(n, CURRENT),
        np.full(n, source_radius),
        strip_nodes,
        strip_triangles,
        strip_stream_function,
    )


def build_distributed_triangle_strips(
    centers: np.ndarray,
    directions: np.ndarray,
    length: float,
    strip_width: float,
    stream_function_jump: float,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    n = centers.shape[0]
    reference = np.tile(np.array([0.0, 0.0, 1.0]), (n, 1))
    parallel = np.abs(np.sum(reference * directions, axis=1)) > 0.9
    reference[parallel] = np.array([0.0, 1.0, 0.0])
    width_dir = np.cross(directions, reference)
    width_norm = np.linalg.norm(width_dir, axis=1)
    width_dir = width_dir / np.maximum(width_norm[:, None], 1.0e-30)

    half_length = 0.5 * length * directions
    half_width = 0.5 * strip_width * width_dir
    nodes = np.empty((4 * n, 3))
    nodes[0::4] = centers - half_length - half_width
    nodes[1::4] = centers - half_length + half_width
    nodes[2::4] = centers + half_length - half_width
    nodes[3::4] = centers + half_length + half_width

    triangles = np.empty((2 * n, 3), dtype=np.int64)
    for i in range(n):
        base = 4 * i
        triangles[2 * i] = [base, base + 2, base + 1]
        triangles[2 * i + 1] = [base + 1, base + 2, base + 3]

    stream_function = np.empty(4 * n)
    stream_function[0::4] = 0.0
    stream_function[2::4] = 0.0
    stream_function[1::4] = stream_function_jump
    stream_function[3::4] = stream_function_jump
    return nodes, triangles, stream_function


def build_segment_dipoles(
    starts: np.ndarray,
    dl: np.ndarray,
    current: np.ndarray,
    source_radius: float,
) -> tuple[tuple[np.ndarray, np.ndarray, np.ndarray], tuple[np.ndarray, np.ndarray, np.ndarray], np.ndarray]:
    centers = starts + 0.5 * dl
    # This is a compact source distribution for comparing direct and hierarchical dipole kernels.
    moments = current[:, None] * dl
    return (
        contiguous_triple(centers[:, 0], centers[:, 1], centers[:, 2]),
        contiguous_triple(moments[:, 0], moments[:, 1], moments[:, 2]),
        np.full(starts.shape[0], source_radius),
    )


def build_triangle_strip(
    helix: np.ndarray,
    strip_width: float,
    stream_function_jump: float,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    offset = np.array([[0.0], [0.5 * strip_width], [0.0]])
    lower = (helix - offset).T
    upper = (helix + offset).T
    n_path = helix.shape[1]
    nodes = np.empty((2 * n_path, 3))
    nodes[0::2] = lower
    nodes[1::2] = upper

    triangles = np.empty((2 * (n_path - 1), 3), dtype=np.int64)
    for i in range(n_path - 1):
        lower0 = 2 * i
        upper0 = lower0 + 1
        lower1 = lower0 + 2
        upper1 = lower0 + 3
        triangles[2 * i] = [lower0, lower1, upper0]
        triangles[2 * i + 1] = [upper0, lower1, upper1]

    stream_function = np.empty(2 * n_path)
    stream_function[0::2] = 0.0
    stream_function[1::2] = stream_function_jump
    return nodes, triangles, stream_function


def source_count_from_log10(log10_source_count: float) -> int:
    count = int(round(10.0 ** float(log10_source_count)))
    return max(MIN_SOURCE_COUNT, min(MAX_SOURCE_COUNT, count))


def field_magnitude(field: tuple[np.ndarray, np.ndarray, np.ndarray]) -> np.ndarray:
    return np.sqrt(field[0] * field[0] + field[1] * field[1] + field[2] * field[2])


def relative_error(
    hierarchical: tuple[np.ndarray, np.ndarray, np.ndarray],
    direct: tuple[np.ndarray, np.ndarray, np.ndarray],
) -> np.ndarray:
    err = field_magnitude(
        (
            hierarchical[0] - direct[0],
            hierarchical[1] - direct[1],
            hierarchical[2] - direct[2],
        )
    )
    ref = np.maximum(field_magnitude(direct), 1e-30)
    return err / ref


def linear_filament_centers(
    geometry: Geometry,
    geometry_layout: str,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    xyzfil, dlxyzfil, _current, _wire_radius = filament_source_arrays(geometry, geometry_layout)
    return (
        np.ascontiguousarray(xyzfil[0] + 0.5 * dlxyzfil[0]),
        np.ascontiguousarray(xyzfil[1] + 0.5 * dlxyzfil[1]),
        np.ascontiguousarray(xyzfil[2] + 0.5 * dlxyzfil[2]),
    )


def triangle_centroids(
    geometry: Geometry,
    geometry_layout: str,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    nodes, triangles, _stream_function = boundary_source_arrays(geometry, geometry_layout)
    tri_nodes = nodes[triangles]
    centroids = np.mean(tri_nodes, axis=1)
    return contiguous_triple(centroids[:, 0], centroids[:, 1], centroids[:, 2])


def dipole_source_arrays(
    geometry: Geometry,
    geometry_layout: str,
) -> tuple[tuple[np.ndarray, np.ndarray, np.ndarray], tuple[np.ndarray, np.ndarray, np.ndarray], np.ndarray]:
    if geometry_layout == "distributed":
        return (
            geometry.volume_dipole_loc,
            geometry.volume_dipole_moment,
            geometry.volume_dipole_outer_radius,
        )
    return geometry.dipole_loc, geometry.dipole_moment, geometry.dipole_outer_radius


def filament_source_arrays(
    geometry: Geometry,
    geometry_layout: str,
) -> tuple[
    tuple[np.ndarray, np.ndarray, np.ndarray],
    tuple[np.ndarray, np.ndarray, np.ndarray],
    np.ndarray,
    np.ndarray,
]:
    if geometry_layout == "distributed":
        return (
            geometry.volume_xyzfil,
            geometry.volume_dlxyzfil,
            geometry.volume_current,
            geometry.volume_wire_radius,
        )
    return geometry.xyzfil, geometry.dlxyzfil, geometry.current, geometry.wire_radius


def boundary_source_arrays(
    geometry: Geometry,
    geometry_layout: str,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    if geometry_layout == "distributed":
        return (
            np.ascontiguousarray(geometry.volume_strip_nodes),
            np.ascontiguousarray(geometry.volume_strip_triangles),
            np.ascontiguousarray(geometry.volume_strip_stream_function),
        )
    return (
        np.ascontiguousarray(geometry.strip_nodes),
        np.ascontiguousarray(geometry.strip_triangles),
        np.ascontiguousarray(geometry.strip_stream_function),
    )


def max_direct_self_sources(source_geometry: str) -> int:
    if source_geometry == "dipole":
        return MAX_DIRECT_SELF_DIPOLE_SOURCES
    if source_geometry == "boundary":
        return MAX_DIRECT_SELF_BOUNDARY_SOURCES
    return MAX_DIRECT_SELF_FILAMENT_SOURCES


def solve_self_fields(
    geometry: Geometry,
    geometry_layout: str,
    source_geometry: str,
    construction_method: str,
    theta: float,
    par: bool,
) -> dict[str, object]:
    if source_geometry == "dipole":
        dipole_loc, dipole_moment, dipole_outer_radius = dipole_source_arrays(geometry, geometry_layout)
        self_obs = dipole_loc
        source_count = dipole_outer_radius.size
    elif source_geometry == "boundary":
        strip_nodes, strip_triangles, strip_stream_function = boundary_source_arrays(
            geometry, geometry_layout
        )
        self_obs = triangle_centroids(geometry, geometry_layout)
        source_count = strip_triangles.shape[0]
    else:
        xyzfil, dlxyzfil, current, wire_radius = filament_source_arrays(geometry, geometry_layout)
        self_obs = linear_filament_centers(geometry, geometry_layout)
        source_count = current.size

    interactions = source_count * source_count
    direct_time: float | None = None
    direct_b = None
    direct_a = None
    if source_count <= max_direct_self_sources(source_geometry):
        t0 = time.perf_counter()
        if source_geometry == "dipole":
            dipole_loc, dipole_moment, dipole_outer_radius = dipole_source_arrays(geometry, geometry_layout)
            direct_b = cfsem.flux_density_dipole(
                dipole_loc,
                dipole_moment,
                self_obs,
                par=par,
                outer_radius=dipole_outer_radius,
            )
            direct_a = cfsem.vector_potential_dipole(
                dipole_loc,
                dipole_moment,
                self_obs,
                par=par,
                outer_radius=dipole_outer_radius,
            )
        elif source_geometry == "boundary":
            strip_nodes, strip_triangles, strip_stream_function = boundary_source_arrays(
                geometry, geometry_layout
            )
            self_obs_array = np.column_stack(self_obs)
            direct_b = cfsem.flux_density_triangle_mesh(
                self_obs_array,
                strip_nodes,
                strip_triangles,
                strip_stream_function,
                par=par,
            )
            direct_a = cfsem.vector_potential_triangle_mesh(
                self_obs_array,
                strip_nodes,
                strip_triangles,
                strip_stream_function,
                par=par,
            )
        else:
            xyzfil, dlxyzfil, current, wire_radius = filament_source_arrays(geometry, geometry_layout)
            direct_b = cfsem.flux_density_linear_filament(
                self_obs,
                xyzfil,
                dlxyzfil,
                current,
                wire_radius,
                par=par,
            )
            direct_a = cfsem.vector_potential_linear_filament(
                self_obs,
                xyzfil,
                dlxyzfil,
                current,
                wire_radius,
                par=par,
            )
        direct_time = time.perf_counter() - t0

    if source_geometry == "dipole":
        dipole_loc, dipole_moment, dipole_outer_radius = dipole_source_arrays(geometry, geometry_layout)
        result_b = cfsem.flux_density_dipole_hierarchical(
            dipole_loc,
            dipole_moment,
            self_obs,
            dipole_outer_radius,
            theta=theta,
            construction_method=construction_method,
            par=par,
        )
        result_a = cfsem.vector_potential_dipole_hierarchical(
            dipole_loc,
            dipole_moment,
            self_obs,
            dipole_outer_radius,
            theta=theta,
            construction_method=construction_method,
            par=par,
        )
    elif source_geometry == "boundary":
        strip_nodes, strip_triangles, strip_stream_function = boundary_source_arrays(
            geometry, geometry_layout
        )
        self_obs_array = np.column_stack(self_obs)
        result_b = cfsem.flux_density_triangle_mesh_hierarchical(
            self_obs_array,
            strip_nodes,
            strip_triangles,
            strip_stream_function,
            theta=theta,
            construction_method=construction_method,
            par=par,
        )
        result_a = cfsem.vector_potential_triangle_mesh_hierarchical(
            self_obs_array,
            strip_nodes,
            strip_triangles,
            strip_stream_function,
            theta=theta,
            construction_method=construction_method,
            par=par,
        )
    else:
        xyzfil, dlxyzfil, current, wire_radius = filament_source_arrays(geometry, geometry_layout)
        result_b = cfsem.flux_density_linear_filament_hierarchical(
            self_obs,
            xyzfil,
            dlxyzfil,
            current,
            wire_radius,
            theta=theta,
            construction_method=construction_method,
            par=par,
        )
        result_a = cfsem.vector_potential_linear_filament_hierarchical(
            self_obs,
            xyzfil,
            dlxyzfil,
            current,
            wire_radius,
            theta=theta,
            construction_method=construction_method,
            par=par,
        )
    hierarchical_b = result_b.field
    hierarchical_a = result_a.field
    build_time = result_b.diagnostics.construction_time + result_a.diagnostics.construction_time
    eval_time = result_b.diagnostics.evaluation_time + result_a.diagnostics.evaluation_time

    return {
        "direct_b": direct_b,
        "direct_a": direct_a,
        "hierarchical_b": hierarchical_b,
        "hierarchical_a": hierarchical_a,
        "interactions": interactions,
        "direct_time": direct_time,
        "direct_skipped": direct_time is None,
        "hierarchical_build_time": build_time,
        "hierarchical_eval_time": eval_time,
    }


def solve_fields(
    geometry: Geometry,
    geometry_layout: str,
    source_geometry: str,
    construction_method: str,
    theta: float,
    par: bool,
    calc_self_field: bool,
    overlay_aabbs: bool,
    overlay_tree_links: bool,
) -> dict[str, object]:
    direct_build_time = 0.0

    t0 = time.perf_counter()
    if source_geometry == "dipole":
        dipole_loc, dipole_moment, dipole_outer_radius = dipole_source_arrays(geometry, geometry_layout)
        direct_b = cfsem.flux_density_dipole(
            dipole_loc,
            dipole_moment,
            geometry.obs,
            par=par,
            outer_radius=dipole_outer_radius,
        )
        direct_a = cfsem.vector_potential_dipole(
            dipole_loc,
            dipole_moment,
            geometry.obs,
            par=par,
            outer_radius=dipole_outer_radius,
        )
    elif source_geometry == "boundary":
        strip_nodes, strip_triangles, strip_stream_function = boundary_source_arrays(
            geometry, geometry_layout
        )
        obs_array = np.column_stack(geometry.obs)
        direct_b = cfsem.flux_density_triangle_mesh(
            obs_array,
            strip_nodes,
            strip_triangles,
            strip_stream_function,
            par=par,
        )
        direct_a = cfsem.vector_potential_triangle_mesh(
            obs_array,
            strip_nodes,
            strip_triangles,
            strip_stream_function,
            par=par,
        )
    else:
        xyzfil, dlxyzfil, current, wire_radius = filament_source_arrays(geometry, geometry_layout)
        direct_b = cfsem.flux_density_linear_filament(
            geometry.obs,
            xyzfil,
            dlxyzfil,
            current,
            wire_radius,
            par=par,
        )
        direct_a = cfsem.vector_potential_linear_filament(
            geometry.obs,
            xyzfil,
            dlxyzfil,
            current,
            wire_radius,
            par=par,
        )
    direct_time = time.perf_counter() - t0

    if source_geometry == "dipole":
        dipole_loc, dipole_moment, dipole_outer_radius = dipole_source_arrays(geometry, geometry_layout)
        source_count = dipole_outer_radius.size
        result_b = cfsem.flux_density_dipole_hierarchical(
            dipole_loc,
            dipole_moment,
            geometry.obs,
            dipole_outer_radius,
            theta=theta,
            construction_method=construction_method,
            par=par,
            extra_diagnostics=True,
        )
        result_a = cfsem.vector_potential_dipole_hierarchical(
            dipole_loc,
            dipole_moment,
            geometry.obs,
            dipole_outer_radius,
            theta=theta,
            construction_method=construction_method,
            par=par,
            extra_diagnostics=True,
        )
    elif source_geometry == "boundary":
        strip_nodes, strip_triangles, strip_stream_function = boundary_source_arrays(
            geometry, geometry_layout
        )
        source_count = strip_triangles.shape[0]
        obs_array = np.column_stack(geometry.obs)
        result_b = cfsem.flux_density_triangle_mesh_hierarchical(
            obs_array,
            strip_nodes,
            strip_triangles,
            strip_stream_function,
            theta=theta,
            construction_method=construction_method,
            par=par,
            extra_diagnostics=True,
        )
        result_a = cfsem.vector_potential_triangle_mesh_hierarchical(
            obs_array,
            strip_nodes,
            strip_triangles,
            strip_stream_function,
            theta=theta,
            construction_method=construction_method,
            par=par,
            extra_diagnostics=True,
        )
    else:
        xyzfil, dlxyzfil, current, wire_radius = filament_source_arrays(geometry, geometry_layout)
        source_count = current.size
        result_b = cfsem.flux_density_linear_filament_hierarchical(
            geometry.obs,
            xyzfil,
            dlxyzfil,
            current,
            wire_radius,
            theta=theta,
            construction_method=construction_method,
            par=par,
            extra_diagnostics=True,
        )
        result_a = cfsem.vector_potential_linear_filament_hierarchical(
            geometry.obs,
            xyzfil,
            dlxyzfil,
            current,
            wire_radius,
            theta=theta,
            construction_method=construction_method,
            par=par,
            extra_diagnostics=True,
        )
    hierarchical_b = result_b.field
    hierarchical_a = result_a.field
    accepted_levels_b = result_b.diagnostics.accepted_levels
    accepted_levels_a = result_a.diagnostics.accepted_levels
    build_time = result_b.diagnostics.construction_time + result_a.diagnostics.construction_time
    eval_time = result_b.diagnostics.evaluation_time + result_a.diagnostics.evaluation_time
    source_tree_aabbs = result_b.diagnostics.source_tree if overlay_aabbs or overlay_tree_links else None

    results: dict[str, object] = {
        "direct_b": direct_b,
        "direct_a": direct_a,
        "direct_build_time": direct_build_time,
        "hierarchical_b": hierarchical_b,
        "hierarchical_a": hierarchical_a,
        "direct_time": direct_time,
        "build_time": build_time,
        "eval_time": eval_time,
        "source_count": source_count,
        "source_target_interactions": source_count * geometry.obs[0].size,
        "geometry_layout": geometry_layout,
    }
    results["accepted_levels_b"] = accepted_levels_b
    results["accepted_levels_a"] = accepted_levels_a
    if source_tree_aabbs is not None:
        results["source_tree_aabbs"] = source_tree_aabbs
    if calc_self_field:
        results["self_field"] = solve_self_fields(
            geometry,
            geometry_layout=geometry_layout,
            source_geometry=source_geometry,
            construction_method=construction_method,
            theta=theta,
            par=par,
        )
    return results


def heatmap_values(values: np.ndarray, geometry: Geometry) -> np.ndarray:
    return values.reshape(geometry.obs_grid[0].shape)


def decimate_path_for_plot(path: np.ndarray, max_points: int = MAX_PLOTTED_PATH_POINTS) -> np.ndarray:
    if path.shape[1] <= max_points:
        return path
    step = max(1, int(np.ceil(path.shape[1] / max_points)))
    decimated = path[:, ::step]
    if decimated.shape[1] == 0 or not np.array_equal(decimated[:, -1], path[:, -1]):
        decimated = np.column_stack((decimated, path[:, -1]))
    return decimated


def selected_aabb_indices(
    levels: np.ndarray,
    max_boxes: int = MAX_PLOTTED_AABBS,
) -> np.ndarray:
    """Return a level-coherent prefix of source-tree node indices for plotting."""

    if levels.size == 0:
        return np.array([], dtype=np.int64)
    indices_by_level: list[np.ndarray] = []
    count = 0
    for level in np.unique(levels.astype(np.int64)):
        level_indices = np.flatnonzero(levels == float(level))
        if count + level_indices.size > max_boxes:
            break
        indices_by_level.append(level_indices)
        count += level_indices.size
    return np.concatenate(indices_by_level) if indices_by_level else np.array([0], dtype=np.int64)


def aabb_overlay_path(
    aabbs: tuple[np.ndarray, ...],
    max_boxes: int = MAX_PLOTTED_AABBS,
) -> tuple[list[float | None], list[float | None]]:
    min_x, _min_y, min_z, max_x, _max_y, max_z, levels = aabbs[:7]
    indices = selected_aabb_indices(levels, max_boxes)

    xs: list[float | None] = []
    zs: list[float | None] = []
    for idx in indices:
        i = int(idx)
        xs.extend([min_x[i], max_x[i], max_x[i], min_x[i], min_x[i], None])
        zs.extend([min_z[i], min_z[i], max_z[i], max_z[i], min_z[i], None])
    return xs, zs


def aabb_child_link_path(
    aabbs: tuple[np.ndarray, ...],
    max_boxes: int = MAX_PLOTTED_AABBS,
) -> tuple[list[float | None], list[float | None]]:
    """Return parent-to-child AABB center links for the plotted source-tree nodes."""

    if len(aabbs) < 9:
        return [], []
    min_x, _min_y, min_z, max_x, _max_y, max_z, levels, left_child, right_child = aabbs[:9]
    indices = selected_aabb_indices(levels, max_boxes)
    if indices.size == 0:
        return [], []
    plotted = set(int(index) for index in indices)
    center_x = 0.5 * (min_x + max_x)
    center_z = 0.5 * (min_z + max_z)

    xs: list[float | None] = []
    zs: list[float | None] = []
    for parent in indices:
        parent_index = int(parent)
        for child in (int(left_child[parent_index]), int(right_child[parent_index])):
            if child not in plotted or child >= min_x.size:
                continue
            xs.extend([float(center_x[parent_index]), float(center_x[child]), None])
            zs.extend([float(center_z[parent_index]), float(center_z[child]), None])
    return xs, zs


def aabb_center_marker_points(
    aabbs: tuple[np.ndarray, ...],
    max_boxes: int = MAX_PLOTTED_AABBS,
) -> tuple[list[float], list[float]]:
    """Return center points for the plotted source-tree AABBs."""

    min_x, _min_y, min_z, max_x, _max_y, max_z, levels = aabbs[:7]
    indices = selected_aabb_indices(levels, max_boxes)
    center_x = 0.5 * (min_x + max_x)
    center_z = 0.5 * (min_z + max_z)
    return (
        [float(center_x[int(index)]) for index in indices],
        [float(center_z[int(index)]) for index in indices],
    )


def source_geometry_overlay_path(
    geometry: Geometry,
    geometry_layout: str,
) -> tuple[list[float | None], list[float | None], list[float | None], list[float | None]]:
    if geometry_layout == "distributed":
        xyzfil, dlxyzfil, _current, _wire_radius = filament_source_arrays(geometry, geometry_layout)
        count = xyzfil[0].size
        if count <= MAX_PLOTTED_PATH_POINTS:
            indices = np.arange(count)
        else:
            rng = np.random.default_rng(DISTRIBUTED_SOURCE_SEED)
            indices = np.sort(rng.choice(count, size=MAX_PLOTTED_PATH_POINTS, replace=False))

        segment_x: list[float | None] = []
        segment_z: list[float | None] = []
        marker_x: list[float | None] = []
        marker_z: list[float | None] = []
        for idx in indices:
            i = int(idx)
            x0 = float(xyzfil[0][i])
            z0 = float(xyzfil[2][i])
            x1 = float(xyzfil[0][i] + dlxyzfil[0][i])
            z1 = float(xyzfil[2][i] + dlxyzfil[2][i])
            segment_x.extend([x0, x1, None])
            segment_z.extend([z0, z1, None])
            marker_x.append(0.5 * (x0 + x1))
            marker_z.append(0.5 * (z0 + z1))
        return segment_x, segment_z, marker_x, marker_z

    centerline_plot = decimate_path_for_plot(geometry.centerline)
    helix_plot = decimate_path_for_plot(geometry.helix)
    return (
        list(centerline_plot[0]),
        list(centerline_plot[2]),
        list(helix_plot[0]),
        list(helix_plot[2]),
    )


def make_figure(
    geometry: Geometry,
    results: dict[str, object],
    field: str,
    show_error: bool,
    overlay_aabbs: bool,
    overlay_tree_links: bool,
):
    import plotly.graph_objects as go
    from plotly.subplots import make_subplots

    direct = results[f"direct_{field}"]
    hierarchical = results[f"hierarchical_{field}"]
    assert isinstance(direct, tuple)
    assert isinstance(hierarchical, tuple)

    left = np.log10(np.maximum(field_magnitude(direct), 1e-30))
    middle = np.log10(np.maximum(field_magnitude(hierarchical), 1e-30))
    right = np.log10(np.maximum(relative_error(hierarchical, direct), 1e-16)) if show_error else middle - left
    right_title = "log10 relative error" if show_error else "log10 magnitude difference"
    accepted_levels = results.get(f"accepted_levels_{field}")
    show_level_diagnostic = isinstance(accepted_levels, np.ndarray)
    subplot_titles = (
        ("Direct", "Hierarchical", right_title, "Accepted Source Level")
        if show_level_diagnostic
        else ("Direct", "Hierarchical", right_title, "")
    )

    fig = make_subplots(
        rows=2,
        cols=2,
        subplot_titles=subplot_titles,
        horizontal_spacing=0.08,
        vertical_spacing=0.13,
    )
    xg, zg = geometry.obs_grid
    geometry_layout = str(results.get("geometry_layout", "helical"))
    geom_x0, geom_z0, geom_x1, geom_z1 = source_geometry_overlay_path(geometry, geometry_layout)
    source_tree_overlay = results.get("source_tree_aabbs")
    aabb_overlay = source_tree_overlay if overlay_aabbs else None
    link_overlay = source_tree_overlay if overlay_tree_links else None
    aabb_x: list[float | None] = []
    aabb_z: list[float | None] = []
    link_x: list[float | None] = []
    link_z: list[float | None] = []
    if isinstance(aabb_overlay, tuple):
        aabb_x, aabb_z = aabb_overlay_path(aabb_overlay)
    if isinstance(link_overlay, tuple):
        link_x, link_z = aabb_child_link_path(link_overlay)
    traces = [
        (left, "log10 |direct|"),
        (middle, "log10 |hierarchical|"),
        (right, right_title),
    ]
    if show_level_diagnostic:
        traces.append((accepted_levels, "accepted levels"))
    for trace_index, (values, title) in enumerate(traces):
        row = trace_index // 2 + 1
        col = trace_index % 2 + 1
        colorscale = "Viridis"
        if trace_index == 2:
            colorscale = [[0.0, "#2c7bb6"], [0.5, "#ffffbf"], [1.0, "#d7191c"]] if show_error else "RdBu"
        if trace_index == 3:
            colorscale = "Cividis"
        fig.add_trace(
            go.Heatmap(
                x=xg[0, :],
                y=zg[:, 0],
                z=heatmap_values(values, geometry),
                showscale=trace_index == 2,
                colorbar={"title": title} if trace_index == 2 else None,
                colorscale=colorscale,
            ),
            row=row,
            col=col,
        )
        fig.add_trace(
            go.Scatter(
                x=geom_x0,
                y=geom_z0,
                mode="lines",
                line={"color": "white", "width": 2, "dash": "dash"},
                showlegend=False,
                hoverinfo="skip",
            ),
            row=row,
            col=col,
        )
        fig.add_trace(
            go.Scatter(
                x=geom_x1,
                y=geom_z1,
                mode="markers" if geometry_layout == "distributed" else "lines",
                line={"color": "black", "width": 1},
                marker={"color": "black", "size": 3} if geometry_layout == "distributed" else None,
                showlegend=False,
                hoverinfo="skip",
            ),
            row=row,
            col=col,
        )
        if link_x:
            fig.add_trace(
                go.Scatter(
                    x=link_x,
                    y=link_z,
                    mode="lines",
                    line={"color": "rgba(10, 70, 160, 0.45)", "width": 1},
                    showlegend=False,
                    hoverinfo="skip",
                ),
                row=row,
                col=col,
            )
        if aabb_x:
            fig.add_trace(
                go.Scatter(
                    x=aabb_x,
                    y=aabb_z,
                    mode="lines",
                    line={"color": "rgba(20, 20, 20, 0.32)", "width": 1},
                    showlegend=False,
                    hoverinfo="skip",
                ),
                row=row,
                col=col,
            )

    for row in range(1, 3):
        for col in range(1, 3):
            axis_index = (row - 1) * 2 + col
            x_axis_name = "x" if axis_index == 1 else f"x{axis_index}"
            fig.update_xaxes(
                title="x [m]" if row == 2 else "",
                showgrid=False,
                row=row,
                col=col,
            )
            fig.update_yaxes(
                title="z [m]" if row == 1 and col == 1 else "",
                scaleanchor=x_axis_name,
                scaleratio=1.0,
                showgrid=False,
                row=row,
                col=col,
            )
    fig.update_layout(
        template="plotly_white",
        height=620,
        margin={"l": 50, "r": 50, "t": 68, "b": 42},
        title=f"{'B-field' if field == 'b' else 'A-field'} comparison on the centerline plane",
    )
    return fig


def make_hierarchical_detail_figure(
    geometry: Geometry,
    results: dict[str, object],
    field: str,
    overlay_aabbs: bool,
    overlay_tree_links: bool,
):
    import plotly.graph_objects as go

    hierarchical = results[f"hierarchical_{field}"]
    assert isinstance(hierarchical, tuple)
    values = np.log10(np.maximum(field_magnitude(hierarchical), 1e-30))
    xg, zg = geometry.obs_grid
    geometry_layout = str(results.get("geometry_layout", "helical"))
    geom_x0, geom_z0, geom_x1, geom_z1 = source_geometry_overlay_path(geometry, geometry_layout)
    source_tree_overlay = results.get("source_tree_aabbs")
    aabb_x: list[float | None] = []
    aabb_z: list[float | None] = []
    link_x: list[float | None] = []
    link_z: list[float | None] = []
    node_x: list[float] = []
    node_z: list[float] = []
    if overlay_aabbs and isinstance(source_tree_overlay, tuple):
        aabb_x, aabb_z = aabb_overlay_path(source_tree_overlay)
    if overlay_tree_links and isinstance(source_tree_overlay, tuple):
        link_x, link_z = aabb_child_link_path(source_tree_overlay)
        node_x, node_z = aabb_center_marker_points(source_tree_overlay)

    fig = go.Figure()
    fig.add_trace(
        go.Heatmap(
            x=xg[0, :],
            y=zg[:, 0],
            z=heatmap_values(values, geometry),
            colorscale="Viridis",
            colorbar={"title": "log10 |hierarchical|"},
        )
    )
    fig.add_trace(
        go.Scatter(
            x=geom_x0,
            y=geom_z0,
            mode="lines",
            line={"color": "white", "width": 2, "dash": "dash"},
            showlegend=False,
            hoverinfo="skip",
        )
    )
    fig.add_trace(
        go.Scatter(
            x=geom_x1,
            y=geom_z1,
            mode="markers" if geometry_layout == "distributed" else "lines",
            line={"color": "black", "width": 2},
            marker={"color": "black", "size": 3} if geometry_layout == "distributed" else None,
            showlegend=False,
            hoverinfo="skip",
        )
    )
    if link_x:
        fig.add_trace(
            go.Scatter(
                x=link_x,
                y=link_z,
                mode="lines",
                line={"color": "white", "width": 2},
                showlegend=False,
                hoverinfo="skip",
            )
        )
    if aabb_x:
        fig.add_trace(
            go.Scatter(
                x=aabb_x,
                y=aabb_z,
                mode="lines",
                line={"color": "black", "width": 2},
                showlegend=False,
                hoverinfo="skip",
            )
        )
    if node_x:
        fig.add_trace(
            go.Scatter(
                x=node_x,
                y=node_z,
                mode="markers",
                marker={
                    "symbol": "circle",
                    "size": 5,
                    "color": "white",
                    "line": {"color": "black", "width": 1},
                },
                showlegend=False,
                hoverinfo="skip",
            )
        )

    fig.update_xaxes(title="x [m]")
    fig.update_yaxes(title="z [m]", scaleanchor="x", scaleratio=1.0)
    fig.update_layout(
        template="plotly_white",
        height=620,
        margin={"l": 50, "r": 55, "t": 54, "b": 45},
        title=f"Standalone hierarchical {'B-field' if field == 'b' else 'A-field'} on the centerline plane",
    )
    return fig


def make_self_field_figure(results: dict[str, object], field: str):
    import plotly.graph_objects as go
    from plotly.subplots import make_subplots

    self_field = results.get("self_field")
    fig = make_subplots(
        rows=1,
        cols=2,
        subplot_titles=("Self Magnitude", "Self Relative Error"),
        horizontal_spacing=0.08,
    )
    fig.update_layout(
        template="plotly_white",
        height=190,
        margin={"l": 40, "r": 40, "t": 48, "b": 24},
        title=f"Self-field {'B' if field == 'b' else 'A'} by source index",
    )
    for axis in fig.select_xaxes():
        axis.update(title="source index")
    for axis in fig.select_yaxes():
        axis.update(title="log10 magnitude")

    if not isinstance(self_field, dict):
        fig.add_annotation(
            text="Enable Calculate self-field to show source-index diagnostics",
            xref="paper",
            yref="paper",
            x=0.5,
            y=0.5,
            showarrow=False,
        )
        return fig

    direct = self_field[f"direct_{field}"]
    hierarchical = self_field[f"hierarchical_{field}"]
    assert isinstance(hierarchical, tuple)
    source_index = np.arange(hierarchical[0].size)
    hierarchical_log = np.log10(np.maximum(field_magnitude(hierarchical), 1e-30))
    fig.add_trace(
        go.Scatter(
            x=source_index,
            y=hierarchical_log,
            mode="lines",
            name="Hierarchical",
            line={"color": "#3b6fb6"},
        ),
        row=1,
        col=1,
    )

    if direct is None:
        fig.add_annotation(
            text="Direct self-field skipped",
            xref="x2 domain",
            yref="y2 domain",
            x=0.5,
            y=0.5,
            showarrow=False,
        )
        return fig

    assert isinstance(direct, tuple)
    direct_log = np.log10(np.maximum(field_magnitude(direct), 1e-30))
    error_log = np.log10(np.maximum(relative_error(hierarchical, direct), 1e-16))
    fig.add_trace(
        go.Scatter(
            x=source_index,
            y=direct_log,
            mode="lines",
            name="Direct",
            line={"color": "#555"},
        ),
        row=1,
        col=1,
    )
    fig.add_trace(
        go.Scatter(
            x=source_index,
            y=error_log,
            mode="lines",
            name="Error",
            line={"color": "#b63b4a"},
        ),
        row=1,
        col=2,
    )
    fig.update_yaxes(title="log10 relative error", row=1, col=2)
    return fig


def speedup_text(direct_time: float, hierarchical_time: float) -> str:
    if hierarchical_time <= 0.0:
        return "n/a"
    return f"{direct_time / hierarchical_time:.2f}x"


def result_float(results: dict[str, object], key: str) -> float:
    value = results[key]
    if isinstance(value, int | float | np.floating):
        return float(value)
    raise TypeError(f"{key} must be numeric")


def make_app():
    from dash import Dash, Input, Output, dcc, html

    app = Dash(__name__)
    sidebar_style = {
        "width": "320px",
        "minWidth": "320px",
        "height": "100vh",
        "overflowY": "auto",
        "padding": "16px",
        "borderRight": "1px solid #d9dde3",
        "boxSizing": "border-box",
        "background": "#f7f8fa",
    }
    content_style = {
        "flex": "1 1 auto",
        "minWidth": "0",
    }
    app.layout = html.Div(
        [
            html.Div(
                [
                    html.Label("Source kernel"),
                    dcc.Dropdown(
                        id="source-geometry",
                        value="linear",
                        clearable=False,
                        options=[
                            {"label": "Linear filament", "value": "linear"},
                            {"label": "Dipole", "value": "dipole"},
                            {"label": "Boundary-element triangle strip", "value": "boundary"},
                        ],
                    ),
                    html.Label("Geometry"),
                    dcc.Dropdown(
                        id="geometry-layout",
                        value="helical",
                        clearable=False,
                        options=[
                            {"label": "Helical loop", "value": "helical"},
                            {"label": "Distributed disk", "value": "distributed"},
                        ],
                    ),
                    html.Label("Twist pitch"),
                    dcc.Slider(
                        id="twist-pitch",
                        min=0.12,
                        max=0.8,
                        step=0.01,
                        value=DEFAULT_TWIST_PITCH,
                        marks={0.12: "0.12", 0.36: "0.36", 0.8: "0.8"},
                    ),
                    html.Label("Helix width"),
                    dcc.Slider(
                        id="helix-width",
                        min=0.0,
                        max=0.16,
                        step=0.005,
                        value=DEFAULT_HELIX_WIDTH,
                        marks={0.0: "0", 0.055: "0.055", 0.16: "0.16"},
                    ),
                    html.Label("Bend curvature"),
                    dcc.Slider(
                        id="bend-curvature",
                        min=0.0,
                        max=DEFAULT_BEND_CURVATURE,
                        step=0.05,
                        value=DEFAULT_BEND_CURVATURE,
                        marks={
                            0.0: "0",
                            0.5: "0.5",
                            1.0: "1.0",
                            round(DEFAULT_BEND_CURVATURE, 2): f"{DEFAULT_BEND_CURVATURE:.2f}",
                        },
                    ),
                    html.Label("Loop fraction"),
                    dcc.Slider(
                        id="loop-fraction",
                        min=0.0,
                        max=1.0,
                        step=0.01,
                        value=DEFAULT_LOOP_FRACTION,
                        marks={0.0: "0", 0.25: "0.25", 0.5: "0.5", 0.75: "0.75", 1.0: "1"},
                    ),
                    html.Label("Source radius"),
                    dcc.Slider(
                        id="source-radius",
                        min=0.0,
                        max=0.08,
                        step=0.001,
                        value=DEFAULT_SOURCE_RADIUS,
                        marks={0.0: "0", 0.02: "0.02", 0.05: "0.05", 0.08: "0.08"},
                    ),
                    html.Label("Construction"),
                    dcc.Dropdown(
                        id="construction",
                        value="longest_axis",
                        clearable=False,
                        options=[
                            {"label": "Longest-axis source tree", "value": "longest_axis"},
                            {"label": "Morton/LBVH source tree", "value": "morton_lbvh"},
                        ],
                    ),
                    html.Label("Field"),
                    dcc.RadioItems(
                        id="field",
                        value="b",
                        inline=True,
                        options=[
                            {"label": "B", "value": "b"},
                            {"label": "A", "value": "a"},
                        ],
                    ),
                    html.Label("Theta"),
                    dcc.Slider(
                        id="theta",
                        min=0.0,
                        max=0.6,
                        step=0.01,
                        value=DEFAULT_THETA,
                        marks={round(i * 0.1, 1): f"{i * 0.1:.1f}" for i in range(7)},
                    ),
                    html.Label("Sources"),
                    dcc.Slider(
                        id="source-count",
                        min=LOG10_MIN_SOURCE_COUNT,
                        max=LOG10_MAX_SOURCE_COUNT,
                        step=0.01,
                        value=LOG10_DEFAULT_SOURCE_COUNT,
                        marks={
                            LOG10_MIN_SOURCE_COUNT: f"{MIN_SOURCE_COUNT}",
                            LOG10_DEFAULT_SOURCE_COUNT: f"{DEFAULT_SOURCE_COUNT}",
                            LOG10_MAX_SOURCE_COUNT: f"{MAX_SOURCE_COUNT:.0E}",
                        },
                    ),
                    dcc.Checklist(
                        id="options",
                        value=["parallel", "relative-error"],
                        options=[
                            {"label": "Parallel evaluation", "value": "parallel"},
                            {"label": "Show relative error", "value": "relative-error"},
                            {"label": "Calculate self-field", "value": "self-field"},
                            {"label": "Overlay source AABBs", "value": "aabbs"},
                            {"label": "Overlay AABB child links", "value": "aabb-links"},
                        ],
                    ),
                ],
                style=sidebar_style,
            ),
            html.Div(
                [
                    html.Div(
                        id="timing",
                        style={
                            "fontFamily": "monospace",
                            "padding": "16px 16px 0",
                            "whiteSpace": "pre-wrap",
                        },
                    ),
                    dcc.Graph(
                        id="field-figure",
                        config={"responsive": True},
                        style={"height": "620px", "marginBottom": "0"},
                    ),
                    dcc.Graph(
                        id="self-field-figure",
                        config={"responsive": True},
                        style={"height": "190px", "marginTop": "0"},
                    ),
                    dcc.Graph(
                        id="hierarchical-detail-figure",
                        config={"responsive": True},
                        style={"height": "620px", "marginTop": "8px"},
                    ),
                ],
                style=content_style,
            ),
        ],
        style={
            "fontFamily": "system-ui, sans-serif",
            "display": "flex",
            "minHeight": "100vh",
        },
    )

    @app.callback(
        Output("field-figure", "figure"),
        Output("self-field-figure", "figure"),
        Output("hierarchical-detail-figure", "figure"),
        Output("timing", "children"),
        Input("source-geometry", "value"),
        Input("geometry-layout", "value"),
        Input("twist-pitch", "value"),
        Input("helix-width", "value"),
        Input("bend-curvature", "value"),
        Input("loop-fraction", "value"),
        Input("source-radius", "value"),
        Input("construction", "value"),
        Input("field", "value"),
        Input("theta", "value"),
        Input("source-count", "value"),
        Input("options", "value"),
    )
    def update(
        source_geometry,
        geometry_layout,
        twist_pitch,
        helix_width,
        bend_curvature,
        loop_fraction,
        source_radius,
        construction,
        field,
        theta,
        log10_source_count,
        options,
    ):
        source_count = source_count_from_log10(float(log10_source_count))
        geometry = build_geometry(
            geometry_layout,
            int(source_count) + 1,
            GRID_N,
            float(twist_pitch),
            float(helix_width),
            float(bend_curvature),
            float(loop_fraction),
            float(source_radius),
        )
        opts = set(options or [])
        results = solve_fields(
            geometry,
            geometry_layout=geometry_layout,
            source_geometry=source_geometry,
            construction_method=str(construction),
            theta=float(theta),
            par="parallel" in opts,
            calc_self_field="self-field" in opts,
            overlay_aabbs="aabbs" in opts,
            overlay_tree_links="aabb-links" in opts,
        )
        fig = make_figure(
            geometry,
            results,
            field,
            "relative-error" in opts,
            "aabbs" in opts,
            "aabb-links" in opts,
        )
        self_fig = make_self_field_figure(results, field)
        detail_fig = make_hierarchical_detail_figure(
            geometry,
            results,
            field,
            "aabbs" in opts,
            "aabb-links" in opts,
        )
        self_field = results.get("self_field")
        self_text = ""
        if isinstance(self_field, dict):
            direct_self = (
                "skipped" if self_field["direct_skipped"] else f"{float(self_field['direct_time']):.3f}s"
            )
            self_speedup = (
                "n/a"
                if self_field["direct_skipped"]
                else speedup_text(
                    float(self_field["direct_time"]),
                    float(self_field["hierarchical_eval_time"]),
                )
            )
            self_text = (
                f"\nself-field source-source interactions={self_field['interactions']:.1E}\n"
                f"self direct:       evaluation={direct_self}\n"
                f"self hierarchical: construction={self_field['hierarchical_build_time']:.3f}s, "
                f"evaluation={self_field['hierarchical_eval_time']:.3f}s, "
                f"speedup={self_speedup}"
            )
        direct_build_time = result_float(results, "direct_build_time")
        direct_time = result_float(results, "direct_time")
        build_time = result_float(results, "build_time")
        eval_time = result_float(results, "eval_time")
        timing = (
            f"nsrc={results['source_count']}, nobs={geometry.obs[0].size}, "
            f"plane={geometry.obs_grid[0].shape[0]}x{geometry.obs_grid[0].shape[1]}\n"
            f"theta={float(theta):.2f}\n"
            f"geometry={geometry_layout}, kernel={source_geometry}, construction={construction}\n"
            f"twist_pitch={float(twist_pitch):.3f}, helix_width={float(helix_width):.3f}, "
            f"bend_curvature={float(bend_curvature):.3f}, loop_fraction={float(loop_fraction):.2f}, "
            f"source_radius={float(source_radius):.3f}\n"
            f"original source-target interactions={results['source_target_interactions']:.1E}\n"
            f"direct:       construction={direct_build_time:.3f}s, "
            f"evaluation={direct_time:.3f}s\n"
            f"hierarchical: construction={build_time:.3f}s, "
            f"evaluation={eval_time:.3f}s, "
            f"speedup={speedup_text(direct_time, eval_time)}"
            f"{self_text}"
        )
        return fig, self_fig, detail_fig, timing

    return app


if __name__ == "__main__":
    if not os.getenv("CFSEM_TESTING"):
        make_app().run(debug=True)
