fix #5

shinolab · Jan 9, 2024 · 4cd462f · 4cd462f
1 parent 0197ae3
commit 4cd462f
Show file tree

Hide file tree

Showing 3 changed files with 16 additions and 8 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -1,3 +1,7 @@
+# 20.0.2
+
+- Fix [#5](https://github.com/shinolab/autd3-rs/issues/5): missing 1/4π factor in propagate
+
 # 20.0.1
 
 - Fix for `sync` feature

diff --git a/autd3-driver/src/acoustics/mod.rs b/autd3-driver/src/acoustics/mod.rs
@@ -4,7 +4,7 @@
  * Created Date: 04/10/2023
  * Author: Shun Suzuki
  * -----
- * Last Modified: 01/12/2023
+ * Last Modified: 09/01/2024
  * Modified By: Shun Suzuki (suzuki@hapis.k.u-tokyo.ac.jp)
  * -----
  * Copyright (c) 2023 Shun Suzuki. All rights reserved.
@@ -14,7 +14,7 @@
 pub mod directivity;
 
 use crate::{
-    defined::{float, Complex},
+    defined::{float, Complex, PI, T4010A1_AMPLITUDE},
     geometry::{Transducer, Vector3},
 };
 
@@ -37,7 +37,9 @@ pub fn propagate<D: Directivity>(
 ) -> Complex {
     let diff = target_pos - tr.position();
     let dist = diff.norm();
-    let r = D::directivity_from_tr(tr, &diff) * (-dist * attenuation).exp() / dist;
+    let r = T4010A1_AMPLITUDE / (4. * PI) / dist
+        * D::directivity_from_tr(tr, &diff)
+        * (-dist * attenuation).exp();
     let phase = -tr.wavenumber(sound_speed) * dist;
     Complex::new(r * phase.cos(), r * phase.sin())
 }
@@ -64,7 +66,7 @@ mod tests {
 
         let tr = crate::geometry::Transducer::new(0, Vector3::zeros(), UnitQuaternion::identity());
 
-        let atten = rng.gen_range(0.0..1.0);
+        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),
@@ -74,8 +76,10 @@ mod tests {
 
         let expect = {
             let dist = target.norm();
-            let r =
-                TestDirectivity::directivity_from_tr(&tr, &target) * (-dist * atten).exp() / dist;
+            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())
         };

diff --git a/autd3-driver/src/defined.rs b/autd3-driver/src/defined.rs
@@ -4,7 +4,7 @@
  * Created Date: 05/12/2022
  * Author: Shun Suzuki
  * -----
- * Last Modified: 19/12/2023
+ * Last Modified: 09/01/2024
  * Modified By: Shun Suzuki (suzuki@hapis.k.u-tokyo.ac.jp)
  * -----
  * Copyright (c) 2022-2023 Shun Suzuki. All rights reserved.
@@ -42,4 +42,4 @@ pub type Complex = nalgebra::Complex<float>;
 
 pub const ABSOLUTE_THRESHOLD_OF_HEARING: float = 20e-6;
 
-pub const T4010A1_AMPLITUDE: float = 21.9295017696 * 200.0 * MILLIMETER; // [Pa*mm]
+pub const T4010A1_AMPLITUDE: float = 275.574246625 * 200.0 * MILLIMETER; // [Pa*mm]