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
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120

pub mod directivity;

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

use directivity::Directivity;

/// Calculate propagation of ultrasound wave
///
/// # Arguments
///
/// * `tr` - Source [Transducer]
/// * `attenuation` - Attenuation coefficient
/// * `wavenumber` - Wavenumber
/// * `target_pos` - Position of target
///
pub fn propagate<D: Directivity>(
    tr: &Transducer,
    attenuation: f64,
    wavenumber: f64,
    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(),
        -wavenumber * dist,
    )
}

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

    use rand::Rng;

    use crate::{
        defined::{mm, FREQ_40K},
        derive::Device,
        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, 1e-6 / mm);
            assert_approx_eq::assert_approx_eq!($a.im, $b.im, 1e-6 / mm);
        };
    }

    #[rstest::fixture]
    fn tr() -> Transducer {
        let mut rng = rand::thread_rng();
        Transducer::new(
            0,
            Vector3::new(
                rng.gen_range(-100.0..100.0),
                rng.gen_range(-100.0..100.0),
                rng.gen_range(-100.0..100.0),
            ),
            UnitQuaternion::from_axis_angle(
                &Vector3::x_axis(),
                rng.gen_range::<f64, _>(-180.0..180.0).to_radians(),
            ) * UnitQuaternion::from_axis_angle(
                &Vector3::y_axis(),
                rng.gen_range::<f64, _>(-180.0..180.0).to_radians(),
            ) * UnitQuaternion::from_axis_angle(
                &Vector3::z_axis(),
                rng.gen_range::<f64, _>(-180.0..180.0).to_radians(),
            ),
        )
    }

    #[rstest::fixture]
    fn target() -> Vector3 {
        let mut rng = rand::thread_rng();
        Vector3::new(
            rng.gen_range(-100.0..100.0),
            rng.gen_range(-100.0..100.0),
            rng.gen_range(-100.0..100.0),
        )
    }

    #[rstest::fixture]
    fn attenuation() -> f64 {
        let mut rng = rand::thread_rng();
        rng.gen_range(0.0..1e-6)
    }

    #[rstest::fixture]
    fn sound_speed() -> f64 {
        let mut rng = rand::thread_rng();
        rng.gen_range(300e3..400e3)
    }

    #[rstest::rstest]
    #[test]
    fn test_propagate(tr: Transducer, target: Vector3, attenuation: f64, sound_speed: f64) {
        let mut device = Device::new(0, vec![tr.clone()], FREQ_40K);
        device.sound_speed = sound_speed;
        let wavenumber = device.wavenumber();
        assert_complex_approx_eq!(
            {
                let diff = target - tr.position();
                let dist = diff.norm();
                let r = T4010A1_AMPLITUDE
                    * TestDirectivity::directivity_from_tr(&tr, &diff)
                    * (-dist * attenuation).exp()
                    / (4. * PI * dist);
                let phase = -wavenumber * dist;
                Complex::new(r * phase.cos(), r * phase.sin())
            },
            super::propagate::<TestDirectivity>(&tr, attenuation, wavenumber, &target)
        );
    }
}