Skip to main content

sim_lib_interference_core/
propagation.rs

1//! Pure `f64` Green-function contributions at one spatial point.
2
3use crate::{
4    Emitter, FieldAmplitude, InterferenceError, InterferenceProblem,
5    POINT_SOURCE_REFERENCE_DISTANCE_METRES, Point3M, Radians, UnitVector3,
6};
7
8fn require_finite(
9    source_id: &str,
10    name: &'static str,
11    value: f64,
12) -> Result<f64, InterferenceError> {
13    value
14        .is_finite()
15        .then_some(value)
16        .ok_or_else(|| InterferenceError::NonFinitePropagation {
17            source_id: source_id.to_owned(),
18            name,
19            value,
20        })
21}
22
23fn cartesian_components(
24    source_id: &str,
25    gain: f64,
26    angle: f64,
27) -> Result<(f64, f64), InterferenceError> {
28    let gain = require_finite(source_id, "gain", gain)?;
29    let angle = require_finite(source_id, "phase-radians", angle)?;
30    Ok((gain * angle.cos(), gain * angle.sin()))
31}
32
33/// Evaluates one emitter's complex field contribution at one point.
34///
35/// The tuple is `(real, imaginary)`. This is the dependency-free `f64`
36/// executable specification for a single source, not a grid solver or an
37/// accumulation policy.
38pub fn contribution_at(
39    problem: &InterferenceProblem,
40    source: &Emitter,
41    at: Point3M,
42) -> Result<(f64, f64), InterferenceError> {
43    match source {
44        Emitter::Point {
45            id,
46            position,
47            amplitude_at_reference,
48            phase,
49        } => point_contribution_at(problem, id, *position, *amplitude_at_reference, *phase, at),
50        Emitter::ForwardPlane {
51            id,
52            through,
53            direction,
54            amplitude,
55            phase,
56        } => {
57            forward_plane_contribution_at(problem, id, *through, *direction, *amplitude, *phase, at)
58        }
59    }
60}
61
62/// Evaluates the outgoing point-source Green function at one point.
63///
64/// For distance `r`, this returns the Cartesian components of
65/// `A * (R_ref / r) * exp(i * k_tilde * r + i * phase)`. A sample at or inside
66/// the problem's singularity radius is rejected rather than clamped.
67pub fn point_contribution_at(
68    problem: &InterferenceProblem,
69    source_id: &str,
70    position: Point3M,
71    amplitude_at_reference: FieldAmplitude,
72    phase: Radians,
73    at: Point3M,
74) -> Result<(f64, f64), InterferenceError> {
75    let distance = require_finite(source_id, "point-distance-metres", at.distance_to(position))?;
76    if distance <= problem.singularity_radius.get() {
77        return Err(InterferenceError::SingularPointSample {
78            source_id: source_id.to_owned(),
79            distance_metres: distance,
80            singularity_radius_metres: problem.singularity_radius.get(),
81        });
82    }
83
84    let wave_number = problem.wavenumber();
85    let phase_advance = require_finite(
86        source_id,
87        "propagation-phase-radians",
88        wave_number.real_radians_per_metre() * distance,
89    )?;
90    let attenuation_exponent = wave_number.imaginary_nepers_per_metre() * distance;
91    let attenuation =
92        if attenuation_exponent.is_infinite() && attenuation_exponent.is_sign_positive() {
93            0.0
94        } else {
95            (-require_finite(source_id, "attenuation-exponent", attenuation_exponent)?).exp()
96        };
97    let gain = amplitude_at_reference.get() * POINT_SOURCE_REFERENCE_DISTANCE_METRES * attenuation
98        / distance;
99
100    cartesian_components(source_id, gain, phase_advance + phase.get())
101}
102
103/// Evaluates a forward-plane Green function at one point.
104///
105/// For signed distance `s`, this returns the Cartesian components of
106/// `A * exp(i * k_tilde * s + i * phase)`. Only the forward half-space
107/// `s >= 0` belongs to this emitter model.
108pub fn forward_plane_contribution_at(
109    problem: &InterferenceProblem,
110    source_id: &str,
111    through: Point3M,
112    direction: UnitVector3,
113    amplitude: FieldAmplitude,
114    phase: Radians,
115    at: Point3M,
116) -> Result<(f64, f64), InterferenceError> {
117    let signed_distance = require_finite(
118        source_id,
119        "plane-signed-distance-metres",
120        direction.signed_distance_metres(through, at),
121    )?;
122    if signed_distance < 0.0 {
123        return Err(InterferenceError::BehindForwardPlane {
124            source_id: source_id.to_owned(),
125            signed_distance_metres: signed_distance,
126        });
127    }
128
129    let wave_number = problem.wavenumber();
130    let phase_advance = require_finite(
131        source_id,
132        "propagation-phase-radians",
133        wave_number.real_radians_per_metre() * signed_distance,
134    )?;
135    let attenuation_exponent = wave_number.imaginary_nepers_per_metre() * signed_distance;
136    let attenuation =
137        if attenuation_exponent.is_infinite() && attenuation_exponent.is_sign_positive() {
138            0.0
139        } else {
140            (-require_finite(source_id, "attenuation-exponent", attenuation_exponent)?).exp()
141        };
142    let gain = amplitude.get() * attenuation;
143
144    cartesian_components(source_id, gain, phase_advance + phase.get())
145}