diff --git a/main.nf b/main.nf index f75c625..75878b2 100644 --- a/main.nf +++ b/main.nf @@ -11,7 +11,7 @@ params.mode = "single" params.throat_radius = 0.05 // Radius of the horn's throat in meters params.mouth_radius = 0.2 // Radius of the horn's mouth in meters params.length = 0.5 // Length of the horn in meters -params.profile = "conical" // Horn flare profile: conical, exponential, hyperbolic +params.profile = "conical" // Horn flare profile: conical, exponential, hyperbolic, tractrix, os, lecleach, cd params.num_sections = 20 // Number of cross-sections for lofting // Simulation Settings @@ -862,7 +862,7 @@ workflow auto { ) // 2. Generate geometry for each profile using throat radius from prescreen - ch_profiles = Channel.from("conical", "exponential", "hyperbolic") + ch_profiles = Channel.from("conical", "exponential", "hyperbolic", "tractrix", "os", "lecleach", "cd") // Extract throat radius from prescreen result ch_throat_radius = ch_prescreen.map { json_file -> @@ -904,7 +904,7 @@ workflow auto { ) // 9. 3D horn geometry renders (one per profile, parallel) - ch_render_profiles = Channel.from("conical", "exponential", "hyperbolic") + ch_render_profiles = Channel.from("conical", "exponential", "hyperbolic", "tractrix", "os", "lecleach", "cd") render_auto_horn_3d(ch_render_profiles) } diff --git a/packages/horn-analysis/src/horn_analysis/horn_render.py b/packages/horn-analysis/src/horn_analysis/horn_render.py index 2cfbea8..0da4d91 100644 --- a/packages/horn-analysis/src/horn_analysis/horn_render.py +++ b/packages/horn-analysis/src/horn_analysis/horn_render.py @@ -52,10 +52,39 @@ def _radius_profile( elif profile == "hyperbolic": m = np.arccosh(mouth_radius / throat_radius) return throat_radius * np.cosh(m * z / length) + elif profile == "tractrix": + t = np.linspace(np.pi - 1e-6, np.pi / 2, 500) + y, x = np.sin(t), np.log(np.tan(t / 2)) + np.cos(t) + x -= x[0] + x_n, y_n = x / x[-1], y / y[-1] + return throat_radius + (mouth_radius - throat_radius) * np.interp(z, x_n * length, y_n) + elif profile == "os": + theta = np.arctan2(np.sqrt(mouth_radius**2 - throat_radius**2), length) + return np.sqrt(throat_radius**2 + (z * np.tan(theta)) ** 2) + elif profile == "lecleach": + t = np.linspace(np.pi - 1e-6, np.pi / 2, 500) + y, x = np.sin(t), np.log(np.tan(t / 2)) + np.cos(t) + x -= x[0] + idx = np.searchsorted(y, throat_radius / mouth_radius) + x_c, y_c = x[idx:] - x[idx], y[idx:] + return np.interp(z, x_c / x_c[-1] * length, y_c / y_c[-1] * mouth_radius) + elif profile == "cd": + frac = 0.3 + z_t = frac * length + r_trans = throat_radius * (mouth_radius / throat_radius) ** frac + result = np.empty_like(z) + mask_exp = z <= z_t + result[mask_exp] = throat_radius * np.exp( + np.log(r_trans / throat_radius) * z[mask_exp] / z_t + ) + result[~mask_exp] = r_trans + (mouth_radius - r_trans) * ( + z[~mask_exp] - z_t + ) / (length - z_t) + return result else: raise ValueError( f"Unknown profile '{profile}'. " - f"Choose from: conical, exponential, hyperbolic" + f"Choose from: conical, exponential, hyperbolic, tractrix, os, lecleach, cd" ) @@ -261,7 +290,8 @@ def main(): 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"], + choices=["conical", "exponential", "hyperbolic", + "tractrix", "os", "lecleach", "cd"], help="Horn flare profile (default: conical)") parser.add_argument("--no-profile-panel", action="store_true", help="Omit the 2D cross-section panel") diff --git a/packages/horn-analysis/src/horn_analysis/html_report.py b/packages/horn-analysis/src/horn_analysis/html_report.py index 9ad3aff..e72d8e2 100644 --- a/packages/horn-analysis/src/horn_analysis/html_report.py +++ b/packages/horn-analysis/src/horn_analysis/html_report.py @@ -26,6 +26,10 @@ "conical": {"badge_bg": "#dbeafe", "badge_fg": "#1e40af"}, "exponential": {"badge_bg": "#dcfce7", "badge_fg": "#166534"}, "hyperbolic": {"badge_bg": "#fef3c7", "badge_fg": "#92400e"}, + "tractrix": {"badge_bg": "#ede9fe", "badge_fg": "#5b21b6"}, + "os": {"badge_bg": "#ccfbf1", "badge_fg": "#115e59"}, + "lecleach": {"badge_bg": "#fce7f3", "badge_fg": "#9d174d"}, + "cd": {"badge_bg": "#e2e8f0", "badge_fg": "#334155"}, } _DEFAULT_BADGE = {"badge_bg": "#f3f4f6", "badge_fg": "#374151"} diff --git a/packages/horn-analysis/src/horn_analysis/plot_theme.py b/packages/horn-analysis/src/horn_analysis/plot_theme.py index 8be7887..c8c77cf 100644 --- a/packages/horn-analysis/src/horn_analysis/plot_theme.py +++ b/packages/horn-analysis/src/horn_analysis/plot_theme.py @@ -42,6 +42,10 @@ "conical": {"color": "#1f4e79", "linestyle": "-"}, "exponential": {"color": "#27864e", "linestyle": "--"}, "hyperbolic": {"color": "#d4831a", "linestyle": "-."}, + "tractrix": {"color": "#7b4ea3", "linestyle": ":"}, + "os": {"color": "#0d9488", "linestyle": "-"}, + "lecleach": {"color": "#be185d", "linestyle": "--"}, + "cd": {"color": "#475569", "linestyle": "-."}, } _DEFAULT_PROFILE_STYLE = {"color": "#6b7280", "linestyle": "-"} diff --git a/packages/horn-analysis/src/horn_analysis/single_report.py b/packages/horn-analysis/src/horn_analysis/single_report.py index 7f4ea08..215411e 100644 --- a/packages/horn-analysis/src/horn_analysis/single_report.py +++ b/packages/horn-analysis/src/horn_analysis/single_report.py @@ -20,6 +20,10 @@ "conical": {"badge_bg": "#dbeafe", "badge_fg": "#1e40af"}, "exponential": {"badge_bg": "#dcfce7", "badge_fg": "#166534"}, "hyperbolic": {"badge_bg": "#fef3c7", "badge_fg": "#92400e"}, + "tractrix": {"badge_bg": "#ede9fe", "badge_fg": "#5b21b6"}, + "os": {"badge_bg": "#ccfbf1", "badge_fg": "#115e59"}, + "lecleach": {"badge_bg": "#fce7f3", "badge_fg": "#9d174d"}, + "cd": {"badge_bg": "#e2e8f0", "badge_fg": "#334155"}, } _DEFAULT_BADGE = {"badge_bg": "#f3f4f6", "badge_fg": "#374151"} diff --git a/packages/horn-analysis/tests/test_analysis.py b/packages/horn-analysis/tests/test_analysis.py index 0a4906f..bd19ec3 100644 --- a/packages/horn-analysis/tests/test_analysis.py +++ b/packages/horn-analysis/tests/test_analysis.py @@ -281,10 +281,15 @@ 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"): + for profile in ("conical", "exponential", "hyperbolic", "os", "cd"): 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) + # Tractrix and Le Cléac'h use parametric interpolation; endpoints are approximate + for profile in ("tractrix", "lecleach"): + r = _radius_profile(np.array([0.0, L]), r_t, r_m, L, profile) + assert r[0] == pytest.approx(r_t, rel=0.02) + assert r[-1] == pytest.approx(r_m, rel=0.02) def test_radius_profile_unknown_raises(self): from horn_analysis.horn_render import _radius_profile diff --git a/packages/horn-core/src/horn_core/candidates.py b/packages/horn-core/src/horn_core/candidates.py index 32fcc61..241c3f9 100644 --- a/packages/horn-core/src/horn_core/candidates.py +++ b/packages/horn-core/src/horn_core/candidates.py @@ -17,7 +17,7 @@ # Standard 1" compression driver throat radii DEFAULT_THROAT_RADII = [0.0125, 0.0175, 0.025] # metres -DEFAULT_PROFILES = ["conical", "exponential", "hyperbolic"] +DEFAULT_PROFILES = ["conical", "exponential", "hyperbolic", "tractrix", "os", "lecleach", "cd"] @dataclass diff --git a/packages/horn-core/src/horn_core/geometry_designer.py b/packages/horn-core/src/horn_core/geometry_designer.py index 10e6d6a..c736f6a 100644 --- a/packages/horn-core/src/horn_core/geometry_designer.py +++ b/packages/horn-core/src/horn_core/geometry_designer.py @@ -26,7 +26,7 @@ # Speed of sound in air at ~20C C0 = 343.0 -DEFAULT_PROFILES = ["conical", "exponential", "hyperbolic"] +DEFAULT_PROFILES = ["conical", "exponential", "hyperbolic", "tractrix", "os", "lecleach", "cd"] @dataclass diff --git a/packages/horn-core/src/horn_core/parameters.py b/packages/horn-core/src/horn_core/parameters.py index dc1afd4..c4acf17 100644 --- a/packages/horn-core/src/horn_core/parameters.py +++ b/packages/horn-core/src/horn_core/parameters.py @@ -10,6 +10,10 @@ class FlareProfile(str, Enum): EXPONENTIAL = "exponential" HYPERBOLIC = "hyperbolic" CONICAL = "conical" + TRACTRIX = "tractrix" + OS = "os" + LECLEACH = "lecleach" + CD = "cd" @dataclass diff --git a/packages/horn-core/tests/test_candidates.py b/packages/horn-core/tests/test_candidates.py index 4110165..3c59315 100644 --- a/packages/horn-core/tests/test_candidates.py +++ b/packages/horn-core/tests/test_candidates.py @@ -47,7 +47,7 @@ def test_all_profiles_represented(self): max_length=0.5, max_mouth_radius=0.2, ) profiles = {c.profile for c in candidates} - assert profiles == {"conical", "exponential", "hyperbolic"} + assert profiles == {"conical", "exponential", "hyperbolic", "tractrix", "os", "lecleach", "cd"} def test_custom_profiles(self): candidates = generate_candidates( @@ -67,15 +67,15 @@ def test_unique_candidate_ids(self): assert len(ids) == len(set(ids)) def test_expected_grid_size(self): - """With default throat radii (3) and 4x4 grid across 3 profiles, - maximum would be 3*3*4*4 = 144 but invalid combos are skipped.""" + """With default throat radii (3) and 4x4 grid across 7 profiles, + maximum would be 7*3*4*4 = 336 but invalid combos are skipped.""" candidates = generate_candidates( target_f_low=500, target_f_high=4000, max_length=0.5, max_mouth_radius=0.2, ) - # Some combos filtered (mouth <= throat), so count should be < 144 - assert len(candidates) <= 144 - assert len(candidates) > 50 # but we should still have plenty + # Some combos filtered (mouth <= throat), so count should be < 336 + assert len(candidates) <= 336 + assert len(candidates) > 100 # but we should still have plenty class TestWriteCandidatesCsv: diff --git a/packages/horn-core/tests/test_geometry_designer.py b/packages/horn-core/tests/test_geometry_designer.py index ac7df79..c3013ca 100644 --- a/packages/horn-core/tests/test_geometry_designer.py +++ b/packages/horn-core/tests/test_geometry_designer.py @@ -99,15 +99,15 @@ def test_sim_range_brackets_target(self): class TestGenerateFullautoCandidates: def test_default_grid_size(self): - """3 profiles x 1 throat x 3 mouth x 3 lengths = 27 max.""" + """7 profiles x 1 throat x 3 mouth x 3 lengths = 63 max.""" candidates, derived = generate_fullauto_candidates( target_f_low=500, target_f_high=4000, throat_radii=[0.025], ) - # All 27 should pass since derived mouth radii >> 0.025 - assert len(candidates) == 27 - assert derived.candidate_count == 27 + # All 63 should pass since derived mouth radii >> 0.025 + assert len(candidates) == 63 + assert derived.candidate_count == 63 def test_all_profiles_represented(self): candidates, _ = generate_fullauto_candidates( @@ -116,7 +116,7 @@ def test_all_profiles_represented(self): throat_radii=[0.025], ) profiles = {c.profile for c in candidates} - assert profiles == {"conical", "exponential", "hyperbolic"} + assert profiles == {"conical", "exponential", "hyperbolic", "tractrix", "os", "lecleach", "cd"} def test_mouth_exceeds_throat(self): """Every candidate must have mouth_radius > throat_radius.""" @@ -163,7 +163,7 @@ def test_length_within_derived_range(self): assert lo - 1e-9 <= c.length <= hi + 1e-9 def test_custom_grid_size(self): - """num_mouth_radii=2, num_lengths=2 -> 3*1*2*2=12 max.""" + """num_mouth_radii=2, num_lengths=2 -> 7*1*2*2=28 max.""" candidates, _ = generate_fullauto_candidates( target_f_low=500, target_f_high=4000, @@ -171,7 +171,7 @@ def test_custom_grid_size(self): num_mouth_radii=2, num_lengths=2, ) - assert len(candidates) == 12 + assert len(candidates) == 28 def test_multiple_throat_radii(self): """Two throat radii should roughly double the candidates.""" @@ -199,8 +199,8 @@ def test_large_throat_filters_small_mouths(self): ) # Should still produce some valid candidates assert len(candidates) > 0 - # But fewer than the full grid - assert len(candidates) < 27 + # But fewer than the full grid (7 profiles x 1 throat x 3 mouth x 3 lengths = 63) + assert len(candidates) < 63 def test_derived_geometry_populated(self): _, derived = generate_fullauto_candidates( diff --git a/packages/horn-geometry/src/horn_geometry/generator.py b/packages/horn-geometry/src/horn_geometry/generator.py index 36218db..5a56a69 100644 --- a/packages/horn-geometry/src/horn_geometry/generator.py +++ b/packages/horn-geometry/src/horn_geometry/generator.py @@ -1,3 +1,5 @@ +import math + import gmsh import numpy as np from pathlib import Path @@ -111,10 +113,105 @@ def radius_func(z: float) -> float: return _loft_horn_profile(radius_func, length, num_sections, output_file, "hyperbolic_horn") +def create_tractrix_horn( + throat_radius: float, + mouth_radius: float, + length: float, + output_file: Path, + num_sections: int = 20, +) -> Path: + """Generate a tractrix horn STEP file. + + Scaled tractrix parametric curve: rapid initial expansion, decelerating toward mouth. + """ + t = np.linspace(np.pi - 1e-6, np.pi / 2, 500) + y, x = np.sin(t), np.log(np.tan(t / 2)) + np.cos(t) + x -= x[0] + x_n, y_n = x / x[-1], y / y[-1] + + def radius_func(z: float) -> float: + return throat_radius + (mouth_radius - throat_radius) * float( + np.interp(z, x_n * length, y_n) + ) + + return _loft_horn_profile(radius_func, length, num_sections, output_file, "tractrix_horn") + + +def create_os_horn( + throat_radius: float, + mouth_radius: float, + length: float, + output_file: Path, + num_sections: int = 20, +) -> Path: + """Generate an oblate spheroidal (OS / Geddes) horn STEP file. + + Simple closed-form: r(z) = sqrt(r_t^2 + (z * tan(theta))^2). + """ + theta = math.atan2(math.sqrt(mouth_radius**2 - throat_radius**2), length) + + def radius_func(z: float) -> float: + return math.sqrt(throat_radius**2 + (z * math.tan(theta)) ** 2) + + return _loft_horn_profile(radius_func, length, num_sections, output_file, "os_horn") + + +def create_lecleach_horn( + throat_radius: float, + mouth_radius: float, + length: float, + output_file: Path, + num_sections: int = 20, +) -> Path: + """Generate a Le Cléac'h horn STEP file. + + Tractrix curve with constant a = mouth_radius, clipped where radius >= throat_radius. + Gentle expansion at throat, aggressive flare at mouth. + """ + t = np.linspace(np.pi - 1e-6, np.pi / 2, 500) + y, x = np.sin(t), np.log(np.tan(t / 2)) + np.cos(t) + x -= x[0] + idx = np.searchsorted(y, throat_radius / mouth_radius) + x_c, y_c = x[idx:] - x[idx], y[idx:] + + def radius_func(z: float) -> float: + return float(np.interp(z, x_c / x_c[-1] * length, y_c / y_c[-1] * mouth_radius)) + + return _loft_horn_profile(radius_func, length, num_sections, output_file, "lecleach_horn") + + +def create_cd_horn( + throat_radius: float, + mouth_radius: float, + length: float, + output_file: Path, + num_sections: int = 30, +) -> Path: + """Generate a constant directivity (CD) horn STEP file. + + Compound: exponential throat (30% of length) transitioning to conical body. + """ + frac = 0.3 + z_t = frac * length + r_trans = throat_radius * (mouth_radius / throat_radius) ** frac + + def radius_func(z: float) -> float: + if z <= z_t: + return throat_radius * math.exp(math.log(r_trans / throat_radius) * z / z_t) + else: + return r_trans + (mouth_radius - r_trans) * (z - z_t) / (length - z_t) + + return _loft_horn_profile(radius_func, length, num_sections, output_file, "cd_horn") + + _PROFILE_DISPATCH = { "conical": create_conical_horn, "exponential": create_exponential_horn, "hyperbolic": create_hyperbolic_horn, + "tractrix": create_tractrix_horn, + "os": create_os_horn, + "lecleach": create_lecleach_horn, + "cd": create_cd_horn, } @@ -160,7 +257,7 @@ def main(): parser.add_argument( "--profile", type=str, - choices=["conical", "exponential", "hyperbolic"], + choices=list(_PROFILE_DISPATCH.keys()), default="conical", help="Horn flare profile (default: conical).", ) diff --git a/packages/horn-geometry/tests/test_profiles.py b/packages/horn-geometry/tests/test_profiles.py index 158744a..4beb0fb 100644 --- a/packages/horn-geometry/tests/test_profiles.py +++ b/packages/horn-geometry/tests/test_profiles.py @@ -13,6 +13,10 @@ create_conical_horn, create_exponential_horn, create_hyperbolic_horn, + create_tractrix_horn, + create_os_horn, + create_lecleach_horn, + create_cd_horn, create_horn, ) @@ -21,7 +25,7 @@ class TestProfileSmoke: """Smoke tests: each profile creates a non-empty STEP file.""" - @pytest.fixture(params=["conical", "exponential", "hyperbolic"]) + @pytest.fixture(params=["conical", "exponential", "hyperbolic", "tractrix", "os", "lecleach", "cd"]) def profile_name(self, request): return request.param @@ -117,9 +121,89 @@ def test_hyperbolic_volume(self, tmp_path): f"Hyperbolic volume mismatch: expected={expected:.6e}, actual={actual:.6e}" ) + def _numerical_volume(self, radius_func, length, n=2000): + """Compute volume by numerical integration: V = integral pi*r(z)^2 dz.""" + z = np.linspace(0, length, n) + r = np.array([radius_func(zi) for zi in z]) + return np.trapz(np.pi * r**2, z) + + def test_tractrix_volume(self, tmp_path): + step = create_tractrix_horn( + self.THROAT_R, self.MOUTH_R, self.LENGTH, + tmp_path / "tractrix.step", num_sections=40, + ) + from horn_geometry.generator import create_tractrix_horn as _ + # Build the same radius_func used internally + t = np.linspace(np.pi - 1e-6, np.pi / 2, 500) + y, x = np.sin(t), np.log(np.tan(t / 2)) + np.cos(t) + x -= x[0] + x_n, y_n = x / x[-1], y / y[-1] + def rf(z): + return self.THROAT_R + (self.MOUTH_R - self.THROAT_R) * float( + np.interp(z, x_n * self.LENGTH, y_n) + ) + expected = self._numerical_volume(rf, self.LENGTH) + actual = self._get_step_volume(step) + assert np.isclose(actual, expected, rtol=0.05), ( + f"Tractrix volume mismatch: expected={expected:.6e}, actual={actual:.6e}" + ) + + def test_os_volume(self, tmp_path): + import math + step = create_os_horn( + self.THROAT_R, self.MOUTH_R, self.LENGTH, + tmp_path / "os.step", num_sections=40, + ) + theta = math.atan2(math.sqrt(self.MOUTH_R**2 - self.THROAT_R**2), self.LENGTH) + def rf(z): + return math.sqrt(self.THROAT_R**2 + (z * math.tan(theta))**2) + expected = self._numerical_volume(rf, self.LENGTH) + actual = self._get_step_volume(step) + assert np.isclose(actual, expected, rtol=0.02), ( + f"OS volume mismatch: expected={expected:.6e}, actual={actual:.6e}" + ) + + def test_lecleach_volume(self, tmp_path): + step = create_lecleach_horn( + self.THROAT_R, self.MOUTH_R, self.LENGTH, + tmp_path / "lecleach.step", num_sections=40, + ) + t = np.linspace(np.pi - 1e-6, np.pi / 2, 500) + y, x = np.sin(t), np.log(np.tan(t / 2)) + np.cos(t) + x -= x[0] + idx = np.searchsorted(y, self.THROAT_R / self.MOUTH_R) + x_c, y_c = x[idx:] - x[idx], y[idx:] + def rf(z): + return float(np.interp(z, x_c / x_c[-1] * self.LENGTH, y_c / y_c[-1] * self.MOUTH_R)) + expected = self._numerical_volume(rf, self.LENGTH) + actual = self._get_step_volume(step) + assert np.isclose(actual, expected, rtol=0.02), ( + f"Le Cléac'h volume mismatch: expected={expected:.6e}, actual={actual:.6e}" + ) + + def test_cd_volume(self, tmp_path): + import math + step = create_cd_horn( + self.THROAT_R, self.MOUTH_R, self.LENGTH, + tmp_path / "cd.step", num_sections=40, + ) + frac = 0.3 + z_t = frac * self.LENGTH + r_trans = self.THROAT_R * (self.MOUTH_R / self.THROAT_R) ** frac + def rf(z): + if z <= z_t: + return self.THROAT_R * math.exp(math.log(r_trans / self.THROAT_R) * z / z_t) + else: + return r_trans + (self.MOUTH_R - r_trans) * (z - z_t) / (self.LENGTH - z_t) + expected = self._numerical_volume(rf, self.LENGTH) + actual = self._get_step_volume(step) + assert np.isclose(actual, expected, rtol=0.02), ( + f"CD volume mismatch: expected={expected:.6e}, actual={actual:.6e}" + ) + def test_throat_mouth_radii_match(self, tmp_path): """Verify that throat (z=0) and mouth (z=L) surfaces have correct areas.""" - for profile in ["exponential", "hyperbolic"]: + for profile in ["exponential", "hyperbolic", "tractrix", "os", "lecleach", "cd"]: step = create_horn( profile=profile, throat_radius=self.THROAT_R, @@ -150,10 +234,12 @@ def test_throat_mouth_radii_match(self, tmp_path): expected_mouth = np.pi * self.MOUTH_R**2 assert throat_area is not None, f"{profile}: no throat face found at z=0" assert mouth_area is not None, f"{profile}: no mouth face found at z=L" - assert np.isclose(throat_area, expected_throat, rtol=0.01), ( + # Le Cléac'h clips the tractrix curve, so throat radius is approximate + tol = 0.02 if profile in ("lecleach", "tractrix") else 0.01 + assert np.isclose(throat_area, expected_throat, rtol=tol), ( f"{profile} throat area: expected={expected_throat:.6e}, actual={throat_area:.6e}" ) - assert np.isclose(mouth_area, expected_mouth, rtol=0.01), ( + assert np.isclose(mouth_area, expected_mouth, rtol=tol), ( f"{profile} mouth area: expected={expected_mouth:.6e}, actual={mouth_area:.6e}" )