1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79

pub mod directivity;

use crate::{
    defined::{float, Complex, PI, T4010A1_AMPLITUDE},
    geometry::{Transducer, Vector3},
};

use directivity::Directivity;

/// Calculate propagation of ultrasound wave
///
/// # Arguments
///
/// * `tr` - Source [Transducer]
/// * `attenuation` - Attenuation coefficient
/// * `sound_speed` - Speed of sound
/// * `target_pos` - Position of target
///
pub fn propagate<D: Directivity>(
    tr: &Transducer,
    attenuation: float,
    sound_speed: float,
    target_pos: &Vector3,
) -> Complex {
    let diff = target_pos - tr.position();
    let dist = diff.norm();
    Complex::from_polar(
        T4010A1_AMPLITUDE / (4. * PI) / dist
            * D::directivity_from_tr(tr, &diff)
            * (-dist * attenuation).exp(),
        -tr.wavenumber(sound_speed) * dist,
    )
}

#[cfg(test)]
mod tests {
    use super::*;

    use rand::Rng;

    use crate::geometry::UnitQuaternion;
    use directivity::tests::TestDirectivity;

    macro_rules! assert_complex_approx_eq {
        ($a:expr, $b:expr) => {
            assert_approx_eq::assert_approx_eq!($a.re, $b.re);
            assert_approx_eq::assert_approx_eq!($a.im, $b.im);
        };
    }

    #[test]
    fn propagate() {
        let mut rng = rand::thread_rng();

        let tr = crate::geometry::Transducer::new(0, Vector3::zeros(), UnitQuaternion::identity());

        let atten = rng.gen_range(0.0..1e-6);
        let c = rng.gen_range(300e3..400e3);
        let target = Vector3::new(
            rng.gen_range(-100.0..100.0),
            rng.gen_range(-100.0..100.0),
            rng.gen_range(-100.0..100.0),
        );

        let expect = {
            let dist = target.norm();
            let r = T4010A1_AMPLITUDE
                * TestDirectivity::directivity_from_tr(&tr, &target)
                * (-dist * atten).exp()
                / (4. * PI * dist);
            let phase = -tr.wavenumber(c) * dist;
            Complex::new(r * phase.cos(), r * phase.sin())
        };
        assert_complex_approx_eq!(
            expect,
            super::propagate::<TestDirectivity>(&tr, atten, c, &target)
        );
    }
}