/
bari
/
HPLSE-poster-code
Обзор
Документация
Войти
/
bari
/
HPLSE-poster-code
Код
Запросы
0
Задачи
Вики
Пакеты
0
Релизы
0
CI/CD
Аналитика
Безопасность
Репозиторий помещен в архив 16 июня 2026. Вся функциональность доступна только для чтения.
main
tests/test_helpinn.py
464 строки
21 KB
bari
Add regression tests for spectral and time-domain propagation
13 апр 2026, 01:26
13 апр 2026, 01:26
74fd376
Код
Авторство
О чём код?
"""Tests for the public helpinn package.""" from __future__ import annotations import unittest import numpy as np from helpinn import GaussianPulse, GaussianTimeImpulse, MonochromaticPulse, SampledProfile, propagate_time_impulse from helpinn.plasma import ( DrudePlasmaStateEquation, ExponentialTailProfile, GaussianAbsorbingProfile, LayeredSlabProfile, LinearRampProfile, PlasmaProfile, PlasmaProfileBase, SinusoidalGratingProfile, StepInterfaceProfile, TanhTransitionProfile, ) from helpinn.pulse import ChirpedGaussianPulse, FlatTopPulse from helpinn.solver import RiccatiSolver, Result, SmmSolver, TimeResult, TmmSolver def fresnel_reflectance(eps_left: float, eps_right: float, angle_rad: float, polarization: str) -> float: """Return the Fresnel reflectance of a single interface.""" n_left = np.sqrt(eps_left + 0j) n_right = np.sqrt(eps_right + 0j) sin_theta_t = n_left * np.sin(angle_rad) / n_right cos_theta_i = np.cos(angle_rad) cos_theta_t = np.sqrt(1.0 - sin_theta_t**2 + 0j) if polarization == "TE": numerator = n_left * cos_theta_i - n_right * cos_theta_t denominator = n_left * cos_theta_i + n_right * cos_theta_t else: numerator = n_right * cos_theta_i - n_left * cos_theta_t denominator = n_right * cos_theta_i + n_left * cos_theta_t return float(np.abs(numerator / denominator) ** 2) def single_interface_profile(eps_left: float, eps_right: float) -> SampledProfile: """Build a one-interface sampled profile.""" return SampledProfile( z=np.array([0.0, 1.0]), eps=np.array([eps_left, eps_right]), eps_in=eps_left, eps_out=eps_right, discontinuity_mask=np.array([False, True]), ) def overdense_drude_slab_profile( wavelength: float = 0.25, omega_p: float = 180.0, collision_rate: float = 2.0, ) -> SampledProfile: """Build a strongly overdense plasma-like slab using Drude numbers.""" z = np.linspace(0.0, 1.0, 400) eps = np.ones_like(z, dtype=complex) omega = 2.0 * np.pi / wavelength state = DrudePlasmaStateEquation(collision_rate=collision_rate) eps_overdense = state.permittivity(np.array([omega]), np.array([omega_p]))[0] slab = (z >= 0.3) & (z <= 0.7) eps[slab] = eps_overdense jumps = np.zeros_like(z, dtype=bool) jumps[np.argmin(np.abs(z - 0.3))] = True jumps[np.argmin(np.abs(z - 0.7))] = True return SampledProfile(z=z, eps=eps, discontinuity_mask=jumps) def overdense_drude_gaussian_profile( wavelength: float = 0.25, omega_p_peak: float = 60.0, collision_rate: float = 2.0, width: float = 0.03, ) -> SampledProfile: """Build a smooth strongly overdense Gaussian plasma-like profile.""" z = np.linspace(0.0, 1.0, 500) omega = 2.0 * np.pi / wavelength state = DrudePlasmaStateEquation(collision_rate=collision_rate) omega_p = omega_p_peak * np.exp(-((z - 0.5) / width) ** 2) eps = state.permittivity(np.full_like(omega_p, omega), omega_p) return SampledProfile(z=z, eps=eps) class HelpinnCoreTests(unittest.TestCase): """Test the core data objects.""" def test_profile_requires_state_equation_for_w(self) -> None: """A w-profile should require a plasma state equation.""" with self.assertRaises(ValueError): SampledProfile(z=np.array([0.0, 1.0]), w=np.array([1.0, 2.0])) def test_profile_and_pulse_validation(self) -> None: """Construct a valid profile and pulse.""" profile = SampledProfile( z=np.array([0.0, 1.0, 2.0]), eps=np.array([1.0, 1.1, 1.2]), deps=np.array([0.1, 0.1, 0.1]), ) pulse = GaussianPulse( wavelength_grid=np.array([1.0, 2.0, 3.0]), center_wavelength=2.0, width=0.5, peak_amplitude=1.0 + 0.0j, theta_rad=np.array([0.0, 0.1, 0.2]), polarization="TE", ) self.assertEqual(profile.z.shape, (3,)) self.assertEqual(pulse.A.shape, (3,)) self.assertIsInstance(profile, PlasmaProfile) def test_gaussian_profile_is_readable(self) -> None: """Construct a concrete Gaussian absorbing profile.""" profile = GaussianAbsorbingProfile( z=np.array([0.0, 1.0, 2.0]), center=1.0, width=1.0, amplitude=0.2, ) self.assertIsInstance(profile, PlasmaProfileBase) self.assertEqual(profile.eps.shape, (3,)) def test_plot_methods_return_figures(self) -> None: """Plot methods should return Matplotlib figures and axes.""" profile = SampledProfile( z=np.array([0.0, 1.0, 2.0]), eps=np.array([1.0, 1.1, 1.2]), ) pulse = GaussianPulse( wavelength_grid=np.array([1.0, 2.0, 3.0]), center_wavelength=2.0, width=0.5, peak_amplitude=1.0 + 0.0j, theta_rad=np.array([0.0, 0.1, 0.2]), polarization="TE", ) result = TmmSolver(profile, pulse).solve() mode0 = result.mode_solution(0) pulse_solution = result.pulse_solution() metadata_pulse = result.export_metadata() metadata_mode0 = result.export_metadata(mode_index=0) self.assertIn("E", mode0) self.assertIn("E", pulse_solution) self.assertEqual(metadata_pulse["kind"], "pulse") self.assertEqual(metadata_mode0["kind"], "mode") self.assertEqual(metadata_mode0["mode_index"], 0) fig1, ax1 = profile.plot() fig_plane, ax_plane = result.plot_plane_wave_coefficients() fig_plane_mode, ax_plane_mode = result.plot_plane_wave_coefficients(mode_index=0) fig2, ax2 = result.plot_fields() fig2_mode, ax2_mode = result.plot_fields(mode_index=0) fig3, ax3 = result.plot_energy() fig3_mode, ax3_mode = result.plot_energy(mode_index=0) fig4, ax4 = result.plot_coefficients() fig4_mode, ax4_mode = result.plot_coefficients(mode_index=0) fig5, ax5 = result.plot_summary() fig5_mode, ax5_mode = result.plot_summary(mode_index=0) self.assertIsNotNone(fig1) self.assertIsNotNone(ax1) self.assertIsNotNone(fig_plane) self.assertIsNotNone(ax_plane) self.assertIsNotNone(fig_plane_mode) self.assertIsNotNone(ax_plane_mode) self.assertIsNotNone(fig2) self.assertIsNotNone(ax2) self.assertIsNotNone(fig2_mode) self.assertIsNotNone(ax2_mode) self.assertIsNotNone(fig3) self.assertIsNotNone(ax3) self.assertIsNotNone(fig3_mode) self.assertIsNotNone(ax3_mode) self.assertIsNotNone(fig4) self.assertIsNotNone(ax4) self.assertIsNotNone(fig4_mode) self.assertIsNotNone(ax4_mode) self.assertIsNotNone(fig5) self.assertIsNotNone(ax5) self.assertIsNotNone(fig5_mode) self.assertIsNotNone(ax5_mode) def test_density_to_plasma_frequency_conversion_is_explicit(self) -> None: """Density-based profiles should use the state-equation conversion.""" state_equation = DrudePlasmaStateEquation(density_scale=4.0) omega_p = state_equation.plasma_frequency_from_density(np.array([1.0, 9.0])) np.testing.assert_allclose(omega_p, np.array([2.0, 6.0])) def test_heatmap_plot_returns_figure(self) -> None: """Heatmap plotting should work for the pulse and a selected mode.""" profile = SampledProfile(z=np.linspace(0.0, 1.0, 100), eps=np.linspace(1.0, 1.5, 100)) pulse = GaussianPulse( wavelength_grid=np.linspace(0.08, 0.12, 5), center_wavelength=0.1, width=0.01, peak_amplitude=1.0 + 0.0j, theta_rad=0.1, polarization="TE", ) result = TmmSolver(profile, pulse).solve() fig1, ax1 = result.plot_heatmap(np.linspace(-0.2, 0.2, 50)) fig2, ax2 = result.plot_heatmap(np.linspace(-0.2, 0.2, 50), mode_index=0) self.assertIsNotNone(fig1) self.assertIsNotNone(ax1) self.assertIsNotNone(fig2) self.assertIsNotNone(ax2) def test_additional_profiles_construct_and_solve(self) -> None: """Additional plasma profiles should work with TMM and SMM.""" z = np.linspace(0.0, 1.0, 300) pulse = MonochromaticPulse(wavelength=0.08, angle_rad=0.1, polarization="TM") profiles = [ StepInterfaceProfile(z=z, interface_z=0.5, eps_left=1.0, eps_right=2.25), LinearRampProfile(z=z, eps_start=1.0, eps_end=2.0 + 0.1j), TanhTransitionProfile(z=z, center=0.5, width=0.08, eps_left=1.0, eps_right=2.2 + 0.1j), SinusoidalGratingProfile(z=z, mean_eps=1.5, amplitude=0.2 + 0.05j, period=0.2), LayeredSlabProfile(z=z, boundaries=(0.3, 0.6), eps_values=(1.0, 2.0 + 0.2j, 1.2)), ExponentialTailProfile(z=z, onset=0.2, scale=0.25, eps_base=1.0, eps_peak=1.8 + 0.1j), ] for profile in profiles: tmm = TmmSolver(profile, pulse).solve() smm = SmmSolver(profile, pulse).solve() self.assertFalse(np.isnan(tmm.cumulative_R)) self.assertFalse(np.isnan(smm.cumulative_R)) self.assertLess(abs(tmm.cumulative_R - smm.cumulative_R), 5e-3) def test_additional_pulse_shapes_construct(self) -> None: """Additional pulse-shape subclasses should validate.""" wavelength_grid = np.linspace(0.08, 0.12, 21) flat = FlatTopPulse(wavelength_grid, wavelength_min=0.09, wavelength_max=0.11) chirped = ChirpedGaussianPulse(wavelength_grid, center_wavelength=0.1, width=0.01, chirp=300.0) self.assertEqual(flat.A.shape, wavelength_grid.shape) self.assertEqual(chirped.A.shape, wavelength_grid.shape) def test_time_impulse_fft_pipeline(self) -> None: """A time-domain impulse should propagate through the FFT pipeline.""" t = np.linspace(-8.0, 8.0, 512, endpoint=False) impulse = GaussianTimeImpulse( t=t, carrier_wavelength=0.25, width=1.0, center=0.0, amplitude=1.0, angle_rad=0.1, polarization="TE", ) profile = SampledProfile(z=np.linspace(0.0, 1.0, 120), eps=np.ones(120)) result = propagate_time_impulse(TmmSolver, profile, impulse, threshold_ratio=1e-4, max_modes=80) self.assertIsInstance(result, TimeResult) self.assertEqual(result.input_signal.shape, t.shape) self.assertEqual(result.transmitted_signal.shape, t.shape) self.assertEqual(result.reflected_signal.shape, t.shape) self.assertEqual(result.E_time.shape, (profile.z.size, t.size)) self.assertGreater(result.spectral_result.mode_count, 0) self.assertIn("selected_mode_count", result.diagnostics) def test_time_impulse_runs_on_multiple_profiles(self) -> None: """The time-domain pipeline should work on several profile families.""" t = np.linspace(-6.0, 6.0, 256, endpoint=False) impulse = GaussianTimeImpulse( t=t, carrier_wavelength=0.18, width=0.8, center=0.0, angle_rad=0.15, polarization="TM", ) z = np.linspace(0.0, 1.0, 160) profiles = [ GaussianAbsorbingProfile(z=z, center=0.5, width=0.12, amplitude=0.5 + 0.1j), TanhTransitionProfile(z=z, center=0.5, width=0.08, eps_left=1.0, eps_right=2.0 + 0.15j), SinusoidalGratingProfile(z=z, mean_eps=1.4, amplitude=0.15 + 0.02j, period=0.18), ] for profile in profiles: result = propagate_time_impulse(TmmSolver, profile, impulse, threshold_ratio=1e-4, max_modes=64) self.assertTrue(np.all(np.isfinite(result.transmitted_signal))) self.assertTrue(np.all(np.isfinite(result.reflected_signal))) def test_overdense_plasma_like_monochromatic_solution(self) -> None: """A strongly overdense plasma-like slab should reflect almost everything.""" profile = overdense_drude_slab_profile() pulse = MonochromaticPulse(wavelength=0.25, angle_rad=0.05, polarization="TE") tmm = TmmSolver(profile, pulse).solve() smm = SmmSolver(profile, pulse).solve() riccati = RiccatiSolver(profile, pulse).solve() for result in (tmm, smm, riccati): self.assertTrue(np.isfinite(result.cumulative_R)) self.assertGreater(result.cumulative_R, 0.95) self.assertLess(smm.cumulative_T, 1e-6) self.assertLess(riccati.cumulative_T, 1e-6) def test_overdense_plasma_like_time_impulse(self) -> None: """A time-domain impulse should be strongly reflected by an overdense slab.""" t = np.linspace(-10.0, 10.0, 512, endpoint=False) impulse = GaussianTimeImpulse( t=t, carrier_wavelength=0.25, width=1.2, center=0.0, amplitude=1.0, angle_rad=0.05, polarization="TE", ) profile = overdense_drude_slab_profile() result = propagate_time_impulse(TmmSolver, profile, impulse, threshold_ratio=1e-4, max_modes=64) transmitted_max = float(np.max(np.abs(result.transmitted_signal))) reflected_max = float(np.max(np.abs(result.reflected_signal))) self.assertLess(transmitted_max, 0.2) self.assertGreater(reflected_max, 0.8) self.assertGreater(result.spectral_result.cumulative_R, 0.95) def test_overdense_gaussian_profile_monochromatic_both_polarizations(self) -> None: """A smooth overdense Gaussian profile should strongly reflect for TE and TM.""" profile = overdense_drude_gaussian_profile() for polarization in ("TE", "TM"): pulse = MonochromaticPulse(wavelength=0.25, angle_rad=0.05, polarization=polarization) tmm = TmmSolver(profile, pulse).solve() smm = SmmSolver(profile, pulse).solve() self.assertTrue(np.isfinite(tmm.cumulative_R)) self.assertTrue(np.isfinite(smm.cumulative_R)) self.assertGreater(tmm.cumulative_R, 0.9) self.assertGreater(smm.cumulative_R, 0.9) self.assertLess(tmm.cumulative_T, 0.05) self.assertLess(smm.cumulative_T, 0.05) self.assertLess(abs(tmm.cumulative_R - smm.cumulative_R), 5e-3) def test_overdense_gaussian_profile_time_impulse_both_polarizations(self) -> None: """Time-domain impulses should stay strongly reflected by an overdense Gaussian profile.""" t = np.linspace(-10.0, 10.0, 512, endpoint=False) profile = overdense_drude_gaussian_profile() for polarization in ("TE", "TM"): impulse = GaussianTimeImpulse( t=t, carrier_wavelength=0.25, width=1.2, center=0.0, amplitude=1.0, angle_rad=0.05, polarization=polarization, ) result = propagate_time_impulse(TmmSolver, profile, impulse, threshold_ratio=1e-4, max_modes=64) transmitted_max = float(np.max(np.abs(result.transmitted_signal))) reflected_max = float(np.max(np.abs(result.reflected_signal))) self.assertLess(transmitted_max, 0.25) self.assertGreater(reflected_max, 0.9) self.assertGreater(result.spectral_result.cumulative_R, 0.9) def test_magnus_diagnostic_is_reported(self) -> None: """Magnus metadata should report the aggregate check result.""" z = np.linspace(0.0, 1.0, 20) profile = SampledProfile(z=z, eps=1.0 + 0.6 * np.exp(-((z - 0.5) / 0.08) ** 2)) pulse = MonochromaticPulse(wavelength=0.01, angle_rad=0.3, polarization="TE") result = TmmSolver(profile, pulse, check_magnus=True).solve() self.assertIn("magnus_value", result.diagnostics) self.assertIn("magnus_ok", result.diagnostics) self.assertGreater(result.diagnostics["magnus_value"], 0.0) class HelpinnSolverTests(unittest.TestCase): """Test the solver backends.""" def setUp(self) -> None: """Create a reusable profile and pulse.""" self.profile = SampledProfile( z=np.array([0.0, 1.0, 2.0, 3.0]), eps=np.array([1.0, 1.05, 1.1, 1.15]), ) self.pulse = GaussianPulse( wavelength_grid=np.array([1.0, 1.5, 2.0, 2.5]), center_wavelength=1.5, width=0.5, peak_amplitude=1.0 + 0.0j, theta_rad=np.array([0.0, 0.05, 0.1, 0.15]), polarization="TM", ) def test_tmm_solver_returns_result(self) -> None: """The TMM solver should produce a result object.""" result = TmmSolver(self.profile, self.pulse).solve() self.assertEqual(result.z.shape, self.profile.z.shape) self.assertIn("solver", result.diagnostics) self.assertTrue(hasattr(TmmSolver, "solve")) self.assertIsInstance(result, Result) def test_smm_solver_returns_result(self) -> None: """The SMM solver should produce a result object.""" result = SmmSolver(self.profile, self.pulse).solve() self.assertEqual(result.z.shape, self.profile.z.shape) self.assertEqual(result.diagnostics.get("solver"), "smm") def test_riccati_solver_returns_result(self) -> None: """The Riccati solver should produce a result object.""" result = RiccatiSolver(self.profile, self.pulse, integration_method="rk4").solve() self.assertEqual(result.z.shape, self.profile.z.shape) self.assertEqual(result.diagnostics.get("integration_method"), "rk4") def test_monochromatic_pulse_is_single_component(self) -> None: """The monochromatic pulse should expose one sample.""" pulse = MonochromaticPulse(wavelength=0.8, polarization="TE") self.assertEqual(pulse.A.shape, (1,)) self.assertEqual(pulse.lambda_.shape, (1,)) def test_solvers_agree_on_lossy_case(self) -> None: """The three solvers should produce consistent lossy propagation.""" z = np.linspace(0.0, 5.0, 300) eps = 1.0 + 1.2 * np.exp(-((z - 2.5) / 0.45) ** 2) + 0.2j * np.exp(-((z - 2.5) / 0.6) ** 2) profile = SampledProfile(z=z, eps=eps) pulse = GaussianPulse( wavelength_grid=np.linspace(0.7, 0.95, 25), center_wavelength=0.8, width=0.04, peak_amplitude=1.0 + 0.0j, theta_rad=0.2, polarization="TM", ) tmm = TmmSolver(profile, pulse).solve() smm = SmmSolver(profile, pulse).solve() riccati = RiccatiSolver(profile, pulse, integration_method="rk4").solve() np.testing.assert_allclose(tmm.plane_wave_A[:, 0], 1.0, atol=1e-12) np.testing.assert_allclose(smm.plane_wave_A[:, 0], 1.0, atol=1e-12) np.testing.assert_allclose(riccati.plane_wave_A[:, 0], 1.0, atol=1e-12) self.assertLess(tmm.cumulative_T, 1.0) self.assertLess(smm.cumulative_T, 1.0) self.assertLess(riccati.cumulative_T, 1.0) self.assertAlmostEqual(tmm.cumulative_R, smm.cumulative_R, places=8) self.assertAlmostEqual(tmm.cumulative_T, smm.cumulative_T, places=8) self.assertLess(abs(smm.cumulative_R - riccati.cumulative_R), 1e-5) self.assertLess(abs(smm.cumulative_T - riccati.cumulative_T), 1e-2) def test_fresnel_normal_incidence_reference(self) -> None: """All solvers should reproduce the normal-incidence Fresnel limit.""" profile = single_interface_profile(1.0, 2.25) expected_R = fresnel_reflectance(1.0, 2.25, 0.0, "TE") pulse = MonochromaticPulse(wavelength=1.0, angle_rad=0.0, polarization="TE") for solver_cls in (TmmSolver, SmmSolver, RiccatiSolver): result = solver_cls(profile, pulse).solve() self.assertAlmostEqual(result.cumulative_R, expected_R, places=8) self.assertAlmostEqual(result.cumulative_T, 1.0 - expected_R, places=8) def test_brewster_angle_reference(self) -> None: """TM reflection should be minimal near Brewster angle.""" eps_left = 1.0 eps_right = 2.25 brewster = float(np.arctan(np.sqrt(eps_right / eps_left))) profile = single_interface_profile(eps_left, eps_right) pulse = MonochromaticPulse(wavelength=1.0, angle_rad=brewster, polarization="TM") for solver_cls in (TmmSolver, SmmSolver, RiccatiSolver): result = solver_cls(profile, pulse).solve() self.assertLess(result.cumulative_R, 1e-8) def test_critical_angle_reference(self) -> None: """TE reflection should approach unity above the critical angle.""" eps_left = 2.25 eps_right = 1.0 critical = float(np.arcsin(np.sqrt(eps_right / eps_left))) profile = single_interface_profile(eps_left, eps_right) pulse = MonochromaticPulse(wavelength=1.0, angle_rad=critical + 0.1, polarization="TE") for solver_cls in (TmmSolver, SmmSolver, RiccatiSolver): result = solver_cls(profile, pulse).solve() self.assertGreater(result.cumulative_R, 0.999999) if __name__ == "__main__": unittest.main()