Merge pull request #91 from Timmmm/Timmmm-acoustic-wave-example

Add acoustic wave example

Author: 9prady9

Cargo.toml

@@ -46,3 +46,7 @@ path = "examples/snow.rs"
 [[example]]
 name = "histogram"
 path = "examples/histogram.rs"
+
+[[example]]
+name = "acoustic_wave"
+path = "examples/acoustic_wave.rs"

examples/acoustic_wave.rs

@@ -0,0 +1,81 @@
+extern crate arrayfire as af;
+
+use af::*;
+use std::f64::consts::*;
+
+fn main() {
+    set_device(0);
+    info();
+
+    acoustic_wave_simulation();
+}
+
+fn normalise(a: &Array) -> Array {
+    (a/(max_all(&abs(a)).0 as f32 * 2.0f32)) + 0.5f32
+}
+fn acoustic_wave_simulation() {
+    // Speed of sound
+    let c = 0.1;
+    // Distance step
+    let dx = 0.5;
+    // Time step
+    let dt = 1.0;
+
+    // Grid size.
+    let nx = 1500;
+    let ny = 1500;
+
+    // Grid dimensions.
+    let dims = Dim4::new(&[nx, ny, 1, 1]);
+
+    // Pressure field
+    let mut p = constant::<f32>(0.0, dims);
+    // d(pressure)/dt field
+    let mut p_dot = p.clone();
+
+    // Laplacian (Del^2) convolution kernel.
+    let laplacian_values = [0.0f32, 1.0, 0.0,
+                        1.0, -4.0, 1.0,
+                        0.0, 1.0, 0.0];
+    let laplacian_kernel = Array::new(&laplacian_values, Dim4::new(&[3, 3, 1, 1])) / (dx * dx);
+
+    // Create a window to show the waves.
+    let mut win = Window::new(1000, 1000, "Waves".to_string());
+
+    // Hann windowed pulse.
+    let pulse_time = 100.0f64;
+    let centre_freq = 0.05;
+
+    // Number of samples in pulse.
+    let pulse_n = (pulse_time/dt).floor() as u64;
+
+    let i = range::<f32>(Dim4::new(&[pulse_n, 1, 1, 1]), 0);
+    let t = i.clone() * dt;
+    let hamming_window = cos(&(i * (2.0 * PI / pulse_n as f64))) * -0.46 + 0.54;
+    let wave = sin(&(&t * centre_freq * 2.0 * PI));
+    let pulse = wave * hamming_window;
+
+    // Iteration count.
+    let mut it = 0;
+
+    while !win.is_closed() {
+        // Convole with laplacian to get spacial second derivative.
+        let lap_p = convolve2(&p, &laplacian_kernel, ConvMode::DEFAULT, ConvDomain::SPATIAL);
+        // Calculate the updated pressure and d(pressure)/dt fields.
+        p_dot += lap_p * (c * dt);
+        p += &p_dot * dt;
+
+        if it < pulse_n {
+            // Location of the source.
+            let seqs = &[Seq::new(700.0, 800.0, 1.0), Seq::new(800.0, 800.0, 1.0)];
+            // Set the pressure there.
+            p = assign_seq(&p, seqs, &index(&pulse, &[Seq::new(it as f64, it as f64, 1.0)]));
+        }
+
+        // Draw the image.
+        win.set_colormap(af::ColorMap::BLUE);
+        win.draw_image(&normalise(&p), None);
+
+        it += 1;
+    }
+}