diff --git a/main.nf b/main.nf index f689482..d686546 100644 --- a/main.nf +++ b/main.nf @@ -26,6 +26,9 @@ params.radiation_model = "plane_wave" // plane_wave, flanged_piston, unflanged_ params.num_bands = 8 // Number of parallel jobs for the solver params.outdir = "./results" +// Directivity (opt-in, single mode only, requires BEM) +params.directivity = false + // Auto mode settings params.target_f_low = 500 params.target_f_high = 4000 @@ -177,6 +180,92 @@ process generate_dashboard { """ } +process render_horn_3d { + publishDir "${params.outdir}", mode: 'copy' + + input: + val throat_radius + val mouth_radius + val length + val profile + + output: + path "horn_3d.png" + + script: + """ + python3 -m horn_analysis.horn_render \ + horn_3d.png \ + --throat-radius ${throat_radius} \ + --mouth-radius ${mouth_radius} \ + --length ${length} \ + --profile ${profile} + """ +} + +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 { + publishDir "${params.outdir}", mode: 'copy' + + input: + path horn_step + + output: + path "directivity.csv" + + script: + """ + python3 -m horn_solver.solver \ + --step-file ${horn_step} \ + --output-file solver_directivity.csv \ + --min-freq ${params.min_freq} \ + --max-freq ${params.max_freq} \ + --num-intervals ${params.num_intervals} \ + --length ${params.length} \ + --mesh-size ${params.mesh_size} \ + --radiation-model bem \ + --compute-directivity \ + --directivity-file directivity.csv + """ +} + +process generate_directivity_plots { + publishDir "${params.outdir}/directivity", mode: 'copy' + + input: + path directivity_csv + + output: + path "polar_directivity.png" + path "directivity_contour.png" + path "beamwidth.png" + path "directivity_index.png" + + script: + """ + python3 -m horn_analysis.directivity_plot ${directivity_csv} --output-dir . + """ +} + // ======================================================================== // Auto mode processes // ======================================================================== @@ -414,6 +503,20 @@ workflow single { // 9. Combined dashboard generate_dashboard(ch_merged_results) + + // 10. 3D horn geometry render (runs in parallel with simulation) + render_horn_3d( + params.throat_radius, + params.mouth_radius, + params.length, + params.profile + ) + + // 11. Directivity (opt-in, requires BEM) + if (params.directivity) { + ch_directivity_csv = run_simulation_directivity(ch_step_file) + generate_directivity_plots(ch_directivity_csv) + } } workflow auto { @@ -468,6 +571,10 @@ workflow auto { ch_drivers_db, ch_prescreen, ) + + // 9. 3D horn geometry renders (one per profile, parallel) + ch_render_profiles = Channel.from("conical", "exponential", "hyperbolic") + render_auto_horn_3d(ch_render_profiles) } workflow { diff --git a/packages/horn-analysis/pyproject.toml b/packages/horn-analysis/pyproject.toml index 9d73c37..892d90b 100644 --- a/packages/horn-analysis/pyproject.toml +++ b/packages/horn-analysis/pyproject.toml @@ -23,6 +23,8 @@ horn-prescreen = "horn_analysis.prescreen:main" horn-rank = "horn_analysis.rank_pipeline:main" horn-auto-report = "horn_analysis.auto_report:main" horn-dashboard = "horn_analysis.dashboard:main" +horn-render = "horn_analysis.horn_render:main" +horn-directivity-plot = "horn_analysis.directivity_plot:main" [project.optional-dependencies] test = [ diff --git a/packages/horn-analysis/src/horn_analysis/directivity_plot.py b/packages/horn-analysis/src/horn_analysis/directivity_plot.py new file mode 100644 index 0000000..0238791 --- /dev/null +++ b/packages/horn-analysis/src/horn_analysis/directivity_plot.py @@ -0,0 +1,358 @@ +"""Akabak-style directivity / radiation pattern visualization. + +Reads a directivity CSV produced by the solver (columns: +``frequency``, ``theta_deg``, ``spl_db``) and generates four +diagnostic plots: + +1. **Polar directivity** — SPL vs angle at selected frequencies. +2. **Directivity contour** (sonogram) — SPL heat-map over frequency and angle. +3. **Beamwidth vs frequency** — coverage angle where SPL drops by a threshold. +4. **Directivity index vs frequency** — on-axis gain relative to omnidirectional. + +All plots use the shared :mod:`horn_analysis.plot_theme` styling. +""" + +import argparse +from pathlib import Path +from typing import List, Optional, Tuple + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + +from horn_analysis import plot_theme + + +# -- Data loading ------------------------------------------------------------- + +def load_directivity(csv_file: str) -> pd.DataFrame: + """Load and validate a directivity CSV file. + + Expected columns: ``frequency``, ``theta_deg``, ``spl_db``. + + Returns + ------- + pd.DataFrame + Sorted by (frequency, theta_deg). + """ + df = pd.read_csv(csv_file) + required = {"frequency", "theta_deg", "spl_db"} + missing = required - set(df.columns) + if missing: + raise ValueError(f"Directivity CSV missing columns: {missing}") + return df.sort_values(["frequency", "theta_deg"]).reset_index(drop=True) + + +# -- Polar directivity ------------------------------------------------------- + +def plot_polar_directivity( + df: pd.DataFrame, + frequencies: Optional[List[float]] = None, + output_file: str = "polar_directivity.png", + *, + n_auto: int = 6, +): + """Polar plot of SPL vs angle at selected frequencies. + + Parameters + ---------- + df : DataFrame + Directivity data with columns ``frequency``, ``theta_deg``, ``spl_db``. + frequencies : list of float, optional + Frequencies to plot. If *None*, auto-selects ``n_auto`` log-spaced + values from the data range. + output_file : str + Output image path. + n_auto : int + Number of frequencies to auto-select when *frequencies* is None. + """ + plot_theme.apply_theme() + + avail_freqs = np.sort(df["frequency"].unique()) + if frequencies is None: + indices = np.round(np.linspace(0, len(avail_freqs) - 1, n_auto)).astype(int) + frequencies = avail_freqs[indices] + else: + # Snap each requested frequency to the nearest available + frequencies = [avail_freqs[np.argmin(np.abs(avail_freqs - f))] for f in frequencies] + + fig, ax = plt.subplots(subplot_kw={"projection": "polar"}, figsize=(8, 8)) + ax.set_theta_zero_location("N") # 0 degrees at top (on-axis) + ax.set_theta_direction(-1) + ax.set_thetamin(0) + ax.set_thetamax(180) + + 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 + 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.set_title("Polar Directivity", pad=20) + ax.legend(loc="lower left", fontsize=7, bbox_to_anchor=(1.05, 0)) + + plot_theme.save_figure(fig, output_file) + + +# -- Directivity contour (sonogram) ------------------------------------------ + +def plot_directivity_contour( + df: pd.DataFrame, + output_file: str = "directivity_contour.png", + *, + db_range: float = 40, +): + """Frequency-angle heat-map of SPL (Akabak-style sonogram). + + Parameters + ---------- + df : DataFrame + Directivity data. + output_file : str + Output image path. + db_range : float + Dynamic range below peak to display (dB). + """ + plot_theme.apply_theme() + + freqs = np.sort(df["frequency"].unique()) + angles = np.sort(df["theta_deg"].unique()) + + spl_grid = np.full((len(angles), len(freqs)), np.nan) + for j, f in enumerate(freqs): + sub = df[df["frequency"] == f].sort_values("theta_deg") + for i, a in enumerate(angles): + row = sub[sub["theta_deg"] == a] + if not row.empty: + spl_grid[i, j] = row["spl_db"].values[0] + + vmax = np.nanmax(spl_grid) + vmin = vmax - db_range + + fig, ax = plot_theme.create_figure(figsize=(11, 6)) + mesh = ax.pcolormesh( + freqs, angles, spl_grid, + cmap="inferno", + vmin=vmin, + vmax=vmax, + shading="gouraud", + ) + plot_theme.setup_freq_axis(ax, freqs.min(), freqs.max()) + ax.set_ylabel("Angle (degrees)") + ax.set_title("Directivity Contour") + cbar = fig.colorbar(mesh, ax=ax, pad=0.02) + cbar.set_label("SPL (dB)") + plot_theme.setup_grid(ax) + + plot_theme.save_figure(fig, output_file) + + +# -- Beamwidth ---------------------------------------------------------------- + +def compute_beamwidth( + df: pd.DataFrame, + threshold_db: float = -6, +) -> Tuple[np.ndarray, np.ndarray]: + """Compute coverage angle vs frequency. + + The beamwidth is twice the angle at which the SPL drops by + ``threshold_db`` relative to the on-axis (0 degrees) value. + + Parameters + ---------- + df : DataFrame + Directivity data. + threshold_db : float + Drop from on-axis SPL (negative, e.g. -6). + + Returns + ------- + (frequencies, beamwidth_deg) : tuple of 1-D arrays + """ + freqs = np.sort(df["frequency"].unique()) + beamwidths = np.full(len(freqs), np.nan) + + for i, f in enumerate(freqs): + sub = df[df["frequency"] == f].sort_values("theta_deg") + theta = sub["theta_deg"].values + spl = sub["spl_db"].values + + # On-axis SPL (theta=0 or closest) + idx_on = np.argmin(np.abs(theta)) + spl_on = spl[idx_on] + target = spl_on + threshold_db # threshold_db is negative + + # Find first angle where SPL drops below target + below = np.where(spl < target)[0] + if len(below) > 0: + idx_cross = below[0] + if idx_cross > 0: + # Linear interpolation between adjacent points + t0, t1 = theta[idx_cross - 1], theta[idx_cross] + s0, s1 = spl[idx_cross - 1], spl[idx_cross] + if s1 != s0: + angle = t0 + (target - s0) * (t1 - t0) / (s1 - s0) + else: + angle = t0 + beamwidths[i] = 2 * angle # Full coverage angle + else: + beamwidths[i] = 0.0 + else: + beamwidths[i] = 2 * theta.max() # Never drops below threshold + + return freqs, beamwidths + + +def plot_beamwidth( + df: pd.DataFrame, + output_file: str = "beamwidth.png", + *, + threshold_db: float = -6, +): + """Plot beamwidth (coverage angle) vs frequency. + + Parameters + ---------- + df : DataFrame + Directivity data. + output_file : str + Output image path. + threshold_db : float + dB threshold for beamwidth calculation. + """ + plot_theme.apply_theme() + + freqs, bw = compute_beamwidth(df, threshold_db=threshold_db) + + fig, ax = plot_theme.create_figure(figsize=(10, 5)) + ax.plot(freqs, bw, color=plot_theme.COLORS["primary"], linewidth=1.4) + plot_theme.setup_freq_axis(ax, freqs.min(), freqs.max()) + ax.set_ylabel("Beamwidth (degrees)") + ax.set_title(f"Beamwidth ({threshold_db} dB) vs Frequency") + plot_theme.setup_grid(ax) + + plot_theme.save_figure(fig, output_file) + + +# -- Directivity index ------------------------------------------------------- + +def compute_directivity_index( + df: pd.DataFrame, +) -> Tuple[np.ndarray, np.ndarray]: + """Compute axisymmetric directivity index vs frequency. + + DI = 10 * log10( p_on_axis^2 / ) + + where is the power-weighted angular average: + = integral( p^2 * sin(theta) dtheta ) / integral( sin(theta) dtheta ) + + Returns + ------- + (frequencies, DI_dB) : tuple of 1-D arrays + """ + freqs = np.sort(df["frequency"].unique()) + di = np.full(len(freqs), np.nan) + + for i, f in enumerate(freqs): + sub = df[df["frequency"] == f].sort_values("theta_deg") + theta_deg = sub["theta_deg"].values + spl = sub["spl_db"].values + + # Convert SPL to linear pressure squared (relative) + p_sq = 10 ** (spl / 10) + + theta_rad = np.radians(theta_deg) + + # On-axis value (closest to 0 degrees) + idx_on = np.argmin(np.abs(theta_deg)) + p_sq_on = p_sq[idx_on] + + # Numerical integration using trapezoidal rule + integrand = p_sq * np.sin(theta_rad) + numerator = np.trapezoid(integrand, theta_rad) + denominator = np.trapezoid(np.sin(theta_rad), theta_rad) + + if denominator > 0 and numerator > 0: + p_sq_avg = numerator / denominator + di[i] = 10 * np.log10(p_sq_on / p_sq_avg) + + return freqs, di + + +def plot_directivity_index( + df: pd.DataFrame, + output_file: str = "directivity_index.png", +): + """Plot directivity index (DI) vs frequency. + + Parameters + ---------- + df : DataFrame + Directivity data. + output_file : str + Output image path. + """ + plot_theme.apply_theme() + + freqs, di = compute_directivity_index(df) + + fig, ax = plot_theme.create_figure(figsize=(10, 5)) + ax.plot(freqs, di, color=plot_theme.COLORS["primary"], linewidth=1.4) + plot_theme.setup_freq_axis(ax, freqs.min(), freqs.max()) + ax.set_ylabel("Directivity Index (dB)") + ax.set_title("Directivity Index vs Frequency") + plot_theme.setup_grid(ax) + + plot_theme.save_figure(fig, output_file) + + +# -- Combined report --------------------------------------------------------- + +def generate_directivity_report(csv_file: str, output_dir: str): + """Generate all four directivity plots into a directory. + + Creates: + - ``polar_directivity.png`` + - ``directivity_contour.png`` + - ``beamwidth.png`` + - ``directivity_index.png`` + + Parameters + ---------- + csv_file : str + Path to the directivity CSV. + output_dir : str + Directory for output images (created if needed). + """ + out = Path(output_dir) + out.mkdir(parents=True, exist_ok=True) + + df = load_directivity(csv_file) + + plot_polar_directivity(df, output_file=str(out / "polar_directivity.png")) + plot_directivity_contour(df, output_file=str(out / "directivity_contour.png")) + plot_beamwidth(df, output_file=str(out / "beamwidth.png")) + plot_directivity_index(df, output_file=str(out / "directivity_index.png")) + + +# -- CLI ---------------------------------------------------------------------- + +def main(): + """Command-line interface: ``horn-directivity-plot``.""" + parser = argparse.ArgumentParser( + description="Generate Akabak-style directivity plots from a CSV." + ) + parser.add_argument("csv_file", help="Directivity CSV (frequency, theta_deg, spl_db)") + parser.add_argument("--output-dir", default="directivity", + help="Output directory (default: directivity/)") + args = parser.parse_args() + + generate_directivity_report(args.csv_file, args.output_dir) + + +if __name__ == "__main__": + main() diff --git a/packages/horn-analysis/src/horn_analysis/horn_render.py b/packages/horn-analysis/src/horn_analysis/horn_render.py new file mode 100644 index 0000000..2cfbea8 --- /dev/null +++ b/packages/horn-analysis/src/horn_analysis/horn_render.py @@ -0,0 +1,283 @@ +"""3D horn geometry rendering with surface-of-revolution visualization. + +Produces a side-by-side figure: 3D surface-of-revolution plot (left) +and 2D wall profile cross-section (right). Uses the shared plot theme +for consistent styling across the analysis package. +""" + +import argparse +import base64 +import io +from pathlib import Path + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +from mpl_toolkits.mplot3d import Axes3D # noqa: F401 – registers 3D projection + +from horn_analysis import plot_theme + + +# -- Radius profile functions ------------------------------------------------ + +def _radius_profile( + z: np.ndarray, + throat_radius: float, + mouth_radius: float, + length: float, + profile: str, +) -> np.ndarray: + """Compute radius along the horn axis for a given flare profile. + + Parameters + ---------- + z : array-like + Axial positions (m), expected in [0, length]. + throat_radius, mouth_radius, length : float + Horn dimensions (m). + profile : str + One of ``"conical"``, ``"exponential"``, ``"hyperbolic"``. + + Returns + ------- + numpy array of radii at each z position. + """ + z = np.asarray(z, dtype=float) + if profile == "conical": + return throat_radius + (mouth_radius - throat_radius) * z / length + elif profile == "exponential": + m = np.log(mouth_radius / throat_radius) + return throat_radius * np.exp(m * z / length) + elif profile == "hyperbolic": + m = np.arccosh(mouth_radius / throat_radius) + return throat_radius * np.cosh(m * z / length) + else: + raise ValueError( + f"Unknown profile '{profile}'. " + f"Choose from: conical, exponential, hyperbolic" + ) + + +# -- 3D rendering ------------------------------------------------------------ + +def render_horn_3d( + throat_radius: float, + mouth_radius: float, + length: float, + profile: str, + output_file: str, + *, + n_z: int = 80, + n_theta: int = 60, + show_profile: bool = True, + figsize: tuple = (14, 6), +): + """Render a 3D surface-of-revolution horn with optional 2D profile panel. + + Parameters + ---------- + throat_radius, mouth_radius, length : float + Horn dimensions (m). + profile : str + Flare profile name. + output_file : str + Path to save the rendered image. + n_z : int + Number of axial samples for the surface mesh. + n_theta : int + Number of angular samples around the axis. + show_profile : bool + If True, add a 2D cross-section panel on the right. + figsize : tuple + Figure size (width, height) in inches. + """ + plot_theme.apply_theme() + + z_vals = np.linspace(0, length, n_z) + r_vals = _radius_profile(z_vals, throat_radius, mouth_radius, length, profile) + theta = np.linspace(0, 2 * np.pi, n_theta) + + Z, Theta = np.meshgrid(z_vals, theta) + R = np.meshgrid(r_vals, theta)[0] + X = R * np.cos(Theta) + Y = R * np.sin(Theta) + + ncols = 2 if show_profile else 1 + fig = plt.figure(figsize=figsize) + + # -- 3D surface panel -- + ax3d = fig.add_subplot(1, ncols, 1, projection="3d") + ax3d.plot_surface( + Z, X, Y, + cmap="coolwarm", + alpha=0.85, + edgecolor="none", + rstride=1, + cstride=1, + ) + 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.view_init(elev=20, azim=-60) + + # Equal aspect ratio for all three axes + max_range = max(length, 2 * mouth_radius) / 2 + mid_z = length / 2 + ax3d.set_xlim(mid_z - max_range, mid_z + max_range) + ax3d.set_ylim(-max_range, max_range) + ax3d.set_zlim(-max_range, max_range) + + # -- 2D profile panel -- + if show_profile: + ax2d = fig.add_subplot(1, 2, 2) + r_mm = r_vals * 1000 + z_mm = z_vals * 1000 + + ax2d.fill_between(z_mm, -r_mm, r_mm, alpha=0.15, color=plot_theme.COLORS["primary"]) + ax2d.plot(z_mm, r_mm, color=plot_theme.COLORS["primary"], linewidth=1.4, label="Wall") + ax2d.plot(z_mm, -r_mm, color=plot_theme.COLORS["primary"], linewidth=1.4) + + # Annotate throat and mouth radii + ax2d.annotate( + f"Throat: {throat_radius * 1000:.1f} mm", + xy=(0, throat_radius * 1000), + xytext=(length * 1000 * 0.15, mouth_radius * 1000 * 0.9), + fontsize=8, + arrowprops=dict(arrowstyle="->", color="#888888", lw=0.8), + color="#555555", + ) + ax2d.annotate( + f"Mouth: {mouth_radius * 1000:.1f} mm", + xy=(length * 1000, mouth_radius * 1000), + xytext=(length * 1000 * 0.65, mouth_radius * 1000 * 0.9), + fontsize=8, + arrowprops=dict(arrowstyle="->", color="#888888", lw=0.8), + color="#555555", + ) + + ax2d.set_xlabel("Axial position (mm)") + ax2d.set_ylabel("Radius (mm)") + ax2d.set_title(f"{profile.capitalize()} Horn — Wall Profile") + ax2d.set_aspect("equal", adjustable="datalim") + plot_theme.setup_grid(ax2d) + + plot_theme.save_figure(fig, output_file) + + +def fig_to_b64_3d( + throat_radius: float, + mouth_radius: float, + length: float, + profile: str, + **kwargs, +) -> str: + """Render a 3D horn to a base64 data-URI string for HTML embedding. + + Accepts the same keyword arguments as :func:`render_horn_3d` + (except ``output_file``). + + Returns + ------- + str + ``data:image/png;base64,...`` encoded image. + """ + plot_theme.apply_theme() + + buf = io.BytesIO() + # Render to a temporary buffer by writing to BytesIO + n_z = kwargs.pop("n_z", 80) + n_theta = kwargs.pop("n_theta", 60) + show_profile = kwargs.pop("show_profile", True) + figsize = kwargs.pop("figsize", (14, 6)) + + z_vals = np.linspace(0, length, n_z) + r_vals = _radius_profile(z_vals, throat_radius, mouth_radius, length, profile) + theta = np.linspace(0, 2 * np.pi, n_theta) + + Z, Theta = np.meshgrid(z_vals, theta) + R = np.meshgrid(r_vals, theta)[0] + X = R * np.cos(Theta) + Y = R * np.sin(Theta) + + ncols = 2 if show_profile else 1 + fig = plt.figure(figsize=figsize) + + ax3d = fig.add_subplot(1, ncols, 1, projection="3d") + ax3d.plot_surface( + Z, X, Y, + cmap="coolwarm", + alpha=0.85, + edgecolor="none", + rstride=1, + cstride=1, + ) + 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.view_init(elev=20, azim=-60) + + max_range = max(length, 2 * mouth_radius) / 2 + mid_z = length / 2 + ax3d.set_xlim(mid_z - max_range, mid_z + max_range) + ax3d.set_ylim(-max_range, max_range) + ax3d.set_zlim(-max_range, max_range) + + if show_profile: + ax2d = fig.add_subplot(1, 2, 2) + r_mm = r_vals * 1000 + z_mm = z_vals * 1000 + ax2d.fill_between(z_mm, -r_mm, r_mm, alpha=0.15, color=plot_theme.COLORS["primary"]) + ax2d.plot(z_mm, r_mm, color=plot_theme.COLORS["primary"], linewidth=1.4) + 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_aspect("equal", adjustable="datalim") + plot_theme.setup_grid(ax2d) + + fig.tight_layout() + fig.savefig(buf, format="png", dpi=150, bbox_inches="tight") + plt.close(fig) + buf.seek(0) + encoded = base64.b64encode(buf.read()).decode("ascii") + return f"data:image/png;base64,{encoded}" + + +# -- CLI ---------------------------------------------------------------------- + +def main(): + """Command-line interface: ``horn-render``.""" + parser = argparse.ArgumentParser( + description="Render a 3D horn geometry visualization." + ) + parser.add_argument("output", help="Output image file path (e.g. horn.png)") + parser.add_argument("--throat-radius", type=float, required=True, + help="Throat radius in metres") + parser.add_argument("--mouth-radius", type=float, required=True, + help="Mouth radius in metres") + parser.add_argument("--length", type=float, required=True, + help="Horn length in metres") + parser.add_argument("--profile", type=str, default="conical", + choices=["conical", "exponential", "hyperbolic"], + help="Horn flare profile (default: conical)") + parser.add_argument("--no-profile-panel", action="store_true", + help="Omit the 2D cross-section panel") + args = parser.parse_args() + + Path(args.output).parent.mkdir(parents=True, exist_ok=True) + + render_horn_3d( + throat_radius=args.throat_radius, + mouth_radius=args.mouth_radius, + length=args.length, + profile=args.profile, + output_file=args.output, + show_profile=not args.no_profile_panel, + ) + + +if __name__ == "__main__": + main() diff --git a/packages/horn-analysis/src/horn_analysis/html_report.py b/packages/horn-analysis/src/horn_analysis/html_report.py index 74338c8..9b42db4 100644 --- a/packages/horn-analysis/src/horn_analysis/html_report.py +++ b/packages/horn-analysis/src/horn_analysis/html_report.py @@ -17,6 +17,7 @@ from horn_core.parameters import DriverParameters from horn_analysis.scoring import TargetSpec from horn_analysis import plot_theme +from horn_analysis.horn_render import fig_to_b64_3d # -- Profile badge styles (HTML-only; matplotlib styles come from plot_theme) -- @@ -309,6 +310,8 @@ def _render_drivers_rows(drivers: Dict[str, DriverParameters]) -> str:
Lowest ripple
{best_ripple}
+{geometry_section} +

Rankings

@@ -365,6 +368,8 @@ def generate_html_report( target: TargetSpec, csv_pairs: List[Tuple[str, str]], top_n: int = 5, + mouth_radius: Optional[float] = None, + length: Optional[float] = None, ) -> str: """Generate a self-contained HTML report string. @@ -376,6 +381,8 @@ def generate_html_report( target: Target frequency specification. csv_pairs: List of (csv_path, label) for coupled SPL plots. top_n: Number of top candidates to include. + mouth_radius: Horn mouth radius in metres (enables 3D geometry renders). + length: Horn length in metres (enables 3D geometry renders). Returns: Complete HTML document as a string. @@ -394,6 +401,24 @@ 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) + geometry_html = "" + if mouth_radius is not None and length is not None: + geom_imgs = [] + for profile in sorted(solver_csvs.keys()): + b64 = fig_to_b64_3d( + throat_radius=throat_radius, + mouth_radius=mouth_radius, + length=length, + profile=profile, + figsize=(12, 5), + ) + geom_imgs.append( + f'
' + ) + geometry_html = "\n".join(geom_imgs) + # Generate plots plot_coupled_spl = _plot_coupled_spl_comparison(csv_pairs, target) plot_raw_spl = _plot_raw_profile_spl(solver_csvs, target) @@ -406,6 +431,17 @@ def generate_html_report( timestamp = datetime.now(timezone.utc).strftime("%Y-%m-%d %H:%M UTC") + # Build geometry section HTML + if geometry_html: + geometry_section = ( + '

Horn Geometry

\n' + '
\n' + f'{geometry_html}\n' + '
' + ) + else: + geometry_section = "" + return _HTML_TEMPLATE.format_map({ "target_low": target.f_low_hz, "target_high": target.f_high_hz, @@ -418,6 +454,7 @@ def generate_html_report( "best_bw": best_bw, "best_sens": best_sens, "best_ripple": best_ripple, + "geometry_section": geometry_section, "rankings_rows": rankings_rows, "plot_coupled_spl": plot_coupled_spl, "plot_raw_spl": plot_raw_spl, diff --git a/packages/horn-analysis/tests/test_analysis.py b/packages/horn-analysis/tests/test_analysis.py index 2733951..0a4906f 100644 --- a/packages/horn-analysis/tests/test_analysis.py +++ b/packages/horn-analysis/tests/test_analysis.py @@ -240,3 +240,163 @@ def test_dashboard_creates_image(self, solver_csv_with_impedance_phase, tmp_path generate_dashboard(str(solver_csv_with_impedance_phase), str(output)) assert output.exists() assert output.stat().st_size > 0 + + +# -- Horn 3D Render Tests --------------------------------------------------- + +class TestHornRender: + def test_render_conical_creates_image(self, tmp_path): + from horn_analysis.horn_render import render_horn_3d + + output = tmp_path / "conical.png" + render_horn_3d(0.05, 0.2, 0.5, "conical", str(output)) + assert output.exists() + assert output.stat().st_size > 0 + + def test_render_exponential_creates_image(self, tmp_path): + from horn_analysis.horn_render import render_horn_3d + + output = tmp_path / "exponential.png" + render_horn_3d(0.05, 0.2, 0.5, "exponential", str(output)) + assert output.exists() + assert output.stat().st_size > 0 + + def test_render_hyperbolic_creates_image(self, tmp_path): + from horn_analysis.horn_render import render_horn_3d + + output = tmp_path / "hyperbolic.png" + render_horn_3d(0.05, 0.2, 0.5, "hyperbolic", str(output)) + assert output.exists() + assert output.stat().st_size > 0 + + def test_render_without_profile_panel(self, tmp_path): + from horn_analysis.horn_render import render_horn_3d + + output = tmp_path / "no_profile.png" + render_horn_3d(0.05, 0.2, 0.5, "conical", str(output), show_profile=False) + assert output.exists() + assert output.stat().st_size > 0 + + def test_radius_profile_values(self): + from horn_analysis.horn_render import _radius_profile + + r_t, r_m, L = 0.05, 0.2, 0.5 + for profile in ("conical", "exponential", "hyperbolic"): + r = _radius_profile(np.array([0.0, L]), r_t, r_m, L, profile) + assert r[0] == pytest.approx(r_t, rel=1e-10) + assert r[-1] == pytest.approx(r_m, rel=1e-10) + + def test_radius_profile_unknown_raises(self): + from horn_analysis.horn_render import _radius_profile + + with pytest.raises(ValueError, match="Unknown profile"): + _radius_profile(np.array([0.0]), 0.05, 0.2, 0.5, "parabolic") + + def test_fig_to_b64_3d_returns_data_uri(self): + from horn_analysis.horn_render import fig_to_b64_3d + + result = fig_to_b64_3d(0.05, 0.2, 0.5, "conical") + assert result.startswith("data:image/png;base64,") + assert len(result) > 100 + + +# -- Directivity Plot Tests -------------------------------------------------- + +@pytest.fixture +def directivity_csv(tmp_path): + """Create a synthetic directivity CSV with cos^2 pattern narrowing with frequency.""" + csv_path = tmp_path / "directivity.csv" + freqs = [500, 1000, 2000, 4000, 8000] + angles = np.arange(0, 181, 5) + rows = [] + for f in freqs: + # Narrower beam at higher frequency + n = f / 500 # exponent increases with frequency + for theta in angles: + # cos^n pattern gives narrowing beam + spl = 90 + 10 * np.log10(max(np.cos(np.radians(theta)) ** n, 1e-10)) + rows.append({"frequency": f, "theta_deg": theta, "spl_db": spl}) + df = pd.DataFrame(rows) + df.to_csv(csv_path, index=False) + return csv_path + + +class TestDirectivityPlot: + def test_load_directivity(self, directivity_csv): + from horn_analysis.directivity_plot import load_directivity + + df = load_directivity(str(directivity_csv)) + assert "frequency" in df.columns + assert "theta_deg" in df.columns + assert "spl_db" in df.columns + + def test_polar_directivity(self, directivity_csv, tmp_path): + from horn_analysis.directivity_plot import load_directivity, plot_polar_directivity + + df = load_directivity(str(directivity_csv)) + output = tmp_path / "polar.png" + plot_polar_directivity(df, output_file=str(output)) + assert output.exists() + assert output.stat().st_size > 0 + + def test_directivity_contour(self, directivity_csv, tmp_path): + from horn_analysis.directivity_plot import load_directivity, plot_directivity_contour + + df = load_directivity(str(directivity_csv)) + output = tmp_path / "contour.png" + plot_directivity_contour(df, output_file=str(output)) + assert output.exists() + assert output.stat().st_size > 0 + + def test_beamwidth(self, directivity_csv, tmp_path): + from horn_analysis.directivity_plot import load_directivity, plot_beamwidth + + df = load_directivity(str(directivity_csv)) + output = tmp_path / "beamwidth.png" + plot_beamwidth(df, output_file=str(output)) + assert output.exists() + assert output.stat().st_size > 0 + + def test_directivity_index(self, directivity_csv, tmp_path): + from horn_analysis.directivity_plot import load_directivity, plot_directivity_index + + df = load_directivity(str(directivity_csv)) + output = tmp_path / "di.png" + plot_directivity_index(df, output_file=str(output)) + assert output.exists() + assert output.stat().st_size > 0 + + def test_compute_beamwidth_values(self, directivity_csv): + from horn_analysis.directivity_plot import load_directivity, compute_beamwidth + + df = load_directivity(str(directivity_csv)) + freqs, bw = compute_beamwidth(df) + # Higher frequency should have narrower beamwidth + assert len(freqs) == 5 + # Compare lowest and highest frequency beamwidths + assert bw[-1] < bw[0], "Higher frequency should have narrower beam" + + def test_compute_directivity_index_values(self, directivity_csv): + from horn_analysis.directivity_plot import load_directivity, compute_directivity_index + + df = load_directivity(str(directivity_csv)) + freqs, di = compute_directivity_index(df) + # Horn concentrates sound -> DI should be positive + assert np.all(di > 0), "Directivity index should be positive for a horn" + + def test_generate_directivity_report(self, directivity_csv, tmp_path): + from horn_analysis.directivity_plot import generate_directivity_report + + output_dir = tmp_path / "directivity_report" + generate_directivity_report(str(directivity_csv), str(output_dir)) + + expected_files = [ + "polar_directivity.png", + "directivity_contour.png", + "beamwidth.png", + "directivity_index.png", + ] + for fname in expected_files: + fpath = output_dir / fname + assert fpath.exists(), f"Missing: {fname}" + assert fpath.stat().st_size > 0, f"Empty: {fname}" diff --git a/packages/horn-solver/src/horn_solver/bem_coupling.py b/packages/horn-solver/src/horn_solver/bem_coupling.py index ee500d9..5cab566 100644 --- a/packages/horn-solver/src/horn_solver/bem_coupling.py +++ b/packages/horn-solver/src/horn_solver/bem_coupling.py @@ -138,6 +138,7 @@ def coupled_solve( outlet_tag: int, k: float, bcs=None, + return_trace_data: bool = False, ): """Solve the coupled FEM-BEM system for exterior radiation at the outlet. @@ -173,11 +174,17 @@ def coupled_solve( Wavenumber. bcs : list Dirichlet boundary conditions. + return_trace_data : bool + If True, return ``(p_h, trace_data)`` where *trace_data* is a dict + with keys ``trace_space``, ``p_trace``, ``dpdn_trace`` needed for + far-field (directivity) computation. Returns ------- p_h : dolfinx.fem.Function Solution pressure field. + trace_data : dict (only when *return_trace_data* is True) + ``{"trace_space": ..., "p_trace": ..., "dpdn_trace": ...}`` """ check_bempp_available() @@ -270,6 +277,13 @@ def coupled_solve( solver.destroy() x_vec.destroy() + if return_trace_data: + trace_data = { + "trace_space": trace_space, + "p_trace": p_trace, + "dpdn_trace": neumann_outlet, + } + return p_h, trace_data return p_h diff --git a/packages/horn-solver/src/horn_solver/solver.py b/packages/horn-solver/src/horn_solver/solver.py index 373a9c8..231e666 100644 --- a/packages/horn-solver/src/horn_solver/solver.py +++ b/packages/horn-solver/src/horn_solver/solver.py @@ -30,6 +30,7 @@ extract_outlet_trace, build_bem_operators, coupled_solve, + compute_far_field, ) except ImportError: BEMPP_AVAILABLE = False @@ -176,6 +177,9 @@ def run_simulation( throat_area: Optional[float] = None, z_horn_initial: Optional[Dict[str, np.ndarray]] = None, radiation_model: str = "plane_wave", + compute_directivity: bool = False, + directivity_file: Optional[str] = None, + directivity_angles: Optional[np.ndarray] = None, ) -> Path: """Run the FEM simulation for the Helmholtz equation. @@ -197,6 +201,12 @@ def run_simulation( ``"flanged_piston"`` (analytical piston in infinite baffle), ``"unflanged_piston"`` (Levine-Schwinger approximation), ``"bem"`` (nonlocal BEM coupling via bempp-cl). + compute_directivity: If True, compute far-field directivity at each + frequency (requires ``radiation_model="bem"``). + directivity_file: Output CSV path for directivity data. Defaults to + ``directivity.csv`` next to the main output file. + directivity_angles: Array of polar angles in degrees for directivity + computation. Defaults to ``np.arange(0, 181, 5)``. Returns: Path to the output CSV file. @@ -213,6 +223,25 @@ def run_simulation( "Run with: mpirun -n 1 python ..." ) + if compute_directivity: + if radiation_model != "bem": + raise ValueError( + "Directivity computation requires --radiation-model bem" + ) + if directivity_angles is None: + directivity_angles = np.arange(0, 181, 5) + if directivity_file is None: + directivity_file = str(Path(output_file).with_name("directivity.csv")) + # Build unit direction vectors for far-field evaluation + # Axisymmetric: directions in the xz-plane, theta from z-axis + theta_rad = np.radians(directivity_angles) + ff_directions = np.column_stack([ + np.sin(theta_rad), + np.zeros(len(theta_rad)), + np.cos(theta_rad), + ]) + directivity_results = [] + if bc_mode == "neumann": if driver is None or throat_area is None or z_horn_initial is None: raise ValueError( @@ -308,7 +337,7 @@ def run_simulation( # --- Solve --- if radiation_model == "bem": - p_h = coupled_solve( + solve_result = coupled_solve( A_fem=a, b_fem=L, V=V, @@ -316,7 +345,28 @@ def run_simulation( outlet_tag=OUTLET_TAG, k=k, bcs=bcs or None, + return_trace_data=compute_directivity, ) + if compute_directivity: + p_h, trace_data = solve_result + # Compute far-field directivity at this frequency + p_far = compute_far_field( + trace_data["trace_space"], + trace_data["p_trace"], + trace_data["dpdn_trace"], + k, + ff_directions, + ) + p_ref = 20e-6 + spl_far = 20 * np.log10(np.abs(p_far) / p_ref + 1e-12) + for angle_deg, spl_val in zip(directivity_angles, spl_far): + directivity_results.append({ + "frequency": frequency, + "theta_deg": float(angle_deg), + "spl_db": float(spl_val), + }) + else: + p_h = solve_result else: problem = LinearProblem( a, L, bcs=bcs, @@ -383,6 +433,13 @@ def run_simulation( results_df = pd.DataFrame(results) results_df.to_csv(output_path, index=False) + # Write directivity results if computed + if compute_directivity and directivity_results: + dir_df = pd.DataFrame(directivity_results) + dir_path = Path(directivity_file) + dir_df.to_csv(dir_path, index=False) + print(f"Directivity data written to: {dir_path}") + print(f"\nSuccessfully ran simulation and generated results: {output_path}") return output_path @@ -414,10 +471,20 @@ def main(): parser.add_argument("--radiation-model", type=str, default="plane_wave", choices=["plane_wave", "flanged_piston", "unflanged_piston", "bem"], help="Radiation impedance model at the outlet (default: plane_wave).") + parser.add_argument("--compute-directivity", action="store_true", + help="Compute far-field directivity (requires --radiation-model bem).") + parser.add_argument("--directivity-file", type=str, default=None, + help="Output CSV path for directivity data.") args = parser.parse_args() # Build extra kwargs for neumann mode - extra_kwargs = {"bc_mode": args.bc_mode, "radiation_model": args.radiation_model} + extra_kwargs = { + "bc_mode": args.bc_mode, + "radiation_model": args.radiation_model, + "compute_directivity": args.compute_directivity, + } + if args.directivity_file: + extra_kwargs["directivity_file"] = args.directivity_file if args.bc_mode == "neumann": if not all([args.driver_json, args.driver_id, args.throat_area, args.phase_a_csv]):