extending time response manual, correcting examples

jondea · Sep 16, 2019 · 124573c · 124573c
1 parent 51c34bd
commit 124573c
Show file tree

Hide file tree

Showing 24 changed files with 200 additions and 154 deletions.
diff --git a/docs/make.jl b/docs/make.jl
@@ -24,17 +24,17 @@ makedocs(
             "library/acoustics.md",
             "library/random.md"
         ],
+        "Theory" => "maths/README.md",
         "Examples" =>
         [
             "example/README.md",
             "example/box_size/README.md",
             "example/hankel_convergence/README.md",
-            "example/lens/README.md",
             "example/moments/README.md",
             "example/near_surface_backscattering/README.md",
             "example/random_particles/README.md",
             "example/time_response_single_particle/README.md",
-            "example/two_particles/README.md",
+            "example/two_particles/README.md"
         ]
     ]
 )

diff --git a/docs/src/assets/backscatter_harmonic.gif b/docs/src/assets/backscatter_harmonic.gif
diff --git a/docs/src/assets/lens-particles.png b/docs/src/assets/lens-particles.png
diff --git a/docs/src/assets/lens-response.png b/docs/src/assets/lens-response.png
diff --git a/docs/src/example/README.md b/docs/src/example/README.md
@@ -1,22 +1,11 @@
-# Examples
+# Overview
 
 ## [Simple random particles](random_particles/README.md)
 - Simple example of random particles in rectangle
 
-## [Two particles](two_particles/README.md)
-- Specify particle positions manually
-- Give different particles different radii and different material properties
-
-## [Lens](lens/README.md)
-- Use the shape functions to create particles in a new geometry
-- See time response from a semi-circular wall
-
 ## [Random particles in a circle](particles_in_circle/README.md)
 - Generate random particles in a circle
 
-## [Make a gif of harmonic scattering](plot/README.md)
-- Generate random particles in a circle
-
 ## [StatisticalMoments](moments/README.md)
 - How to extract statistical information from a batch of simulations, in this
   case: mean, standard deviation, skew and kurtosis (also known as moments).

diff --git a/docs/src/example/lens/README.md b/docs/src/example/lens/README.md
diff --git a/docs/src/example/lens/lens.jl b/docs/src/example/lens/lens.jl
diff --git a/docs/src/example/lens/plot_lens.png b/docs/src/example/lens/plot_lens.png
diff --git a/docs/src/example/moments/README.md b/docs/src/example/moments/README.md
@@ -1,64 +1,99 @@
 # StatisticalMoments
 
+```@meta
+DocTestSetup = quote
+    using MultipleScattering
+end
+```
+
 Here we are going to simulate the scattered wave for many different configurations of particles. We can then take the average and standard deviation (the moments) of the scattered wave. In statistical mechanics this process is called [ensemble average](https://en.wikipedia.org/wiki/Ensemble_average_(statistical_mechanics)).
 
-## Choose the type of particles
-```julia
-using MultipleScattering
+## Region and particles properties
 
-volfrac = 0.01
-radius = 1.0
-num_particles = 10
+First we choose the region to place particles and the receiver position:
+```jldoctest moments; output = false
+bottomleft = [0.0;-25.0]
+topright = [50.0;25.0]
+shape = Rectangle(bottomleft, topright)
+x = [-10.0,0.0]
+
+# output
 
-# region to place the particles
-shape = Rectangle(volfrac, radius, num_particles)
+2-element Array{Float64,1}:
+ -10.0
+   0.0
 ```
-To see the region where the particle will be placed, and the receiver position:
 ```julia
 using Plots
-pyplot()
-listener = [-10.0, 0.0]
 plot(shape);
-scatter!([listener[1]],[listener[2]]);
-plot_shape = annotate!([(listener_position[1], listener_position[2] -2., "Receiver")])
+scatter!([x[1]],[x[2]], label="");
+plot_shape = annotate!([(x[1], x[2] -2., "Receiver")])
 ```
 ![Plot of shape and receiver](shape_receiver.png)
 
-## Calculate the moments of the scattered wave
-The code below chooses a random (uniform distribution) configuration of particles inside `shape` and calculates the received signal at `listener` for wavenumbers `k_arr`,
-```julia
-k_arr = collect(LinRange(0.01,1.0,100))
-simulation = FrequencySimulation(volfrac,radius,k_arr; shape=shape, listener_positions = listener, seed = 1)
-plot(simulation)
-```
-![Plot of response against wavenumber](plot_simulation.png)
+Next we fill this `shape` with a random (uniform distribution) configuration of particles:
+```jldoctest moments
+volfrac = 0.05
+radius = 1.0
+
+particles = random_particles(Acoustic(2; ρ=0.0, c=0.0), Circle(radius);
+        region_shape = shape,
+        volume_fraction = volfrac,
+        seed=2
+);
 
+typeof(particles)
+
+# output
+Array{AbstractParticle{Float64,2},1}
+```
 To see the position of the chosen particles:
 ```julia
 plot(plot_shape)
-plot!.(simulation.particles);
+plot!(particles);
 plot!()
 ```
-![Plot of the position of the particles](plot_particles.png)
+![Plot particles](plot_particles.png)
 
+Scattering a plane-wave from these particles  
+```julia
+ωs = LinRange(0.01,1.0,100)
+plane_wave = plane_source(Acoustic(1.0, 1.0, 2);
+    direction = [1.0,0.0], position = x);
+sim = FrequencySimulation(particles, plane_wave);
+```
+```julia
+plot(run(sim,x,ωs))
+```
+![](plot_result.png)
+
+
+## The moments of the scattered wave
 Now we will do simulations for particles placed in many different configurations and take the moments:
 ```julia
-simulations = [
-    FrequencySimulation(volfrac,radius,k_arr; shape=shape, listener_positions = listener, seed = i)
-for i = 1:20]
-real_moments = StatisticalMoments(simulations; response_apply=real) # moments of the real part
-plot(real_moments);
-plot!(xlabel="wavenumbers", title="Moments of the real part")
+results = map(1:20) do i
+    particles = random_particles(Acoustic(2; ρ=0.0, c=0.0), Circle(radius);
+            region_shape = shape,
+            volume_fraction = volfrac,
+            seed=i
+    )
+    run(FrequencySimulation(particles, plane_wave), x, ωs)
+end
+
+# package Plots changed it's argument, the below no longer works..
+# num_moments = 3
+# plot(results; field_apply = real, num_moments = num_moments)
+# plot!(xlabel="wavenumbers", title="Moments of the real part")
 ```
 ![Moments of the real part the scattered waves](plot_moments.png)
 
 ## Calculate the moments of the scattered wave in time
 ```julia
-time_simulations = TimeSimulation.(simulations)
-time_simulations[1].time_arr # the time_arr chosen will be based on the discrete Fourier transform of simulations[1].k_arr
-real_time_moments = StatisticalMoments(time_simulations; response_apply=real) # moments of the real part
-plot(real_time_moments,xlims=(0,300));
-plot!(xlabel="time", title="Moments of the real part of the time wave")
+time_simulations = frequency_to_time.(results)
+time_simulations[1].t # the time_arr chosen will be based on the discrete Fourier transform of simulations[1].k_arr
+# real_time_moments = StatisticalMoments(time_simulations; response_apply=real) # moments of the real part
+# plot(real_time_moments,xlims=(0,300));
+# plot!(xlabel="time", title="Moments of the real part of the time wave")
 ```
 ![Moments of the real part the scattered waves in time](plot_time_moments.png)
 

diff --git a/docs/src/example/moments/plot_particles.png b/docs/src/example/moments/plot_particles.png
diff --git a/docs/src/example/moments/plot_result.png b/docs/src/example/moments/plot_result.png
diff --git a/docs/src/example/moments/shape_receiver.png b/docs/src/example/moments/shape_receiver.png
diff --git a/docs/src/example/plot/README.md b/docs/src/example/plot/README.md
diff --git a/docs/src/example/plot/backscatter_harmonic.gif b/docs/src/example/plot/backscatter_harmonic.gif
diff --git a/docs/src/example/plot/harmonic_field.png b/docs/src/example/plot/harmonic_field.png
diff --git a/docs/src/example/two_particles/README.md b/docs/src/example/two_particles/README.md
@@ -1,3 +1,5 @@
+# OUT-DATED needs to be updated
+
 # Two particles
 
 Define two particles with the first centred at [1.,-2.], with radius 1.0, sound speed 2.0 and density 10.0

diff --git a/docs/src/example/two_particles/two_particles.jl b/docs/src/example/two_particles/two_particles.jl
@@ -24,7 +24,7 @@ simulation = FrequencySimulation(particles, w_arr;
 
 # to plot the frequency response over a region that includes the particles
 w = 3.2
-plot(simulation,w; res=80, resp_fnc=abs)
+plot(simulation,w; res=80, field_apply=abs)
 # the green circle in the plot is the reciever position
 
 TimeSimulation = TimeSimulation(simulation)

diff --git a/docs/src/manual/plot.md b/docs/src/manual/plot.md
@@ -1,3 +1,47 @@
 # Plotting
 
 Here will be examples of videos and more...
+
+## GIF - Harmonic from random particles
+
+```julia
+using MultipleScattering
+using Plots; pyplot()
+
+num_particles = 70
+radius = 1.0
+ω = 1.0
+
+host_medium = Acoustic(1.0, 1.0, 2)
+particle_medium = Acoustic(0.2, 0.3, 2)
+particle_shape = Circle(radius)
+
+max_width = 50*radius
+bottomleft = [0.,-max_width]
+topright = [max_width,max_width]
+shape = Rectangle(bottomleft,topright)
+
+particles = random_particles(particle_medium, particle_shape; region_shape = shape, num_particles = num_particles)
+
+source =  plane_source(host_medium; direction = [1.0,0.5])
+
+simulation = FrequencySimulation(particles, source)
+
+bottomleft = [-25.,-max_width]
+bounds = Rectangle(bottomleft,topright)
+result = run(simulation, bounds, [ω]; res=100)
+
+ts = LinRange(0.,2pi/ω,30)
+
+maxc = round(10*maximum(real.(field(result))))/10
+minc = round(10*minimum(real.(field(result))))/10
+
+anim = @animate for t in ts
+    plot(result,ω; seriestype = :contour, phase_time=t, clim=(minc,maxc), c=:balance)
+    plot!(simulation)
+    plot!(colorbar=false, title="",axis=false, xlab="",ylab="")
+end
+#
+gif(anim,"backscatter_harmonic.gif", fps = 7)
+```
+![backscattering from harmonic wave](../assets/backscatter_harmonic.gif)
diff --git a/docs/src/manual/time_response.md b/docs/src/manual/time_response.md
@@ -122,6 +122,52 @@ time_response = frequency_to_time(freq_response; t_vec = t_vec, discrete_impulse
 ```
 ![](../assets/triangle_time_response.png)
 
+## Scattering example
+As an example, we will make a reflective lens out of particles. To achieve this we will place the particles into a region with the shape [`TimeOfFlight`](@ref).
+
+First we choose the properties of the lens:
+```julia
+p_radius = 0.1
+volfrac = 0.3
+
+x = [-10.0;0.0]
+outertime = 34.8
+innertime = 34.0
+
+# Generate particles which are at most outertime away from our listener
+outershape = TimeOfFlight(x, outertime)
+outerparticles = random_particles(Acoustic(2; ρ=0.0, c=0.0), Circle(p_radius);
+        region_shape = outershape,
+        volume_fraction = volfrac,
+        seed=2
+);
+
+# Filter out particles which are less than innertime away
+innershape = TimeOfFlight(x, innertime + 4*p_radius); # the + 4*p_radius is to account for double the particle diameter
+particles = filter(p -> p⊈innershape, outerparticles);
+
+plot(particles)
+```
+![](../assets/lens-particles.png)
+
+Next we simulate an impulse plane-wave starting at $x = -10$:
+```julia
+ωs = LinRange(0.01,2.0,100)
+
+plane_wave = plane_source(Acoustic(1.0, 1.0, 2); direction = [1.0,0.0], position = x);
+sim = FrequencySimulation(particles, plane_wave);
+
+freq_response = run(sim, x, ωs);
+
+t_vec = -10.:0.2:81.
+time_response = frequency_to_time(freq_response; t_vec=t_vec, impulse = GaussianImpulse(1.5, 1.0))
+
+xticks = [0.,20.,34.,40.0,60.,80.]
+plot(time_response, title="Time response from lens", label="", xticks=xticks)
+```
+![](../assets/lens-response.png)
+
+The first peak is the incident wave, and the next peak is the wave scattered from the lens which should arrive close to `t = 34`.
 
 ## [Technical details](@id impulse_details)
 
@@ -145,4 +191,4 @@ where $\omega_M = \Omega$ and $\Delta \omega_m$ depends on the scheme used, with
 To learn more see the notes [Discrete Fourier Transform](../maths/DiscreteFourier.pdf) or the tests in the folder test of the source code.
 
 !!! tip
-    The most standard way to sample the frequencies is to take $\omega_m = m \Delta \omega$ and $\Delta \omega_m = \Delta \omega$, for some fixed $\Delta \omega$. If we substitute this sampling into $\phi(\mathbf x, t)$ we find that $\phi$ becomes periodic in time with period $T = \frac{2\pi n}{\Delta \omega}$, that is $\phi(\mathbf x, t + T) = \phi(\mathbf x, t)$ for every $t$. Suppose you wanted to calculate a scattered wave that arrives at time $t = T + \tau$, what would happen? The answer is you would see this scattered wave arrive at time $t = \tau$ if $\tau < T$. Note, this occurs often in this package because it is easy to get trapped waves when there is strong multiple scattering.
+    The standard way to sample the frequencies is to take $\omega_m = m \Delta \omega$ and $\Delta \omega_m = \Delta \omega$ for some constant $\Delta \omega$. If we substitute this sampling into the approximation for $\phi(\mathbf x, t)$, shown above, we find that $\phi(\mathbf x,t)$ becomes periodic in time $t$ with period $T = 2\pi / \Delta \omega$. That is $\phi(\mathbf x, t + T) = \phi(\mathbf x, t)$ for every $t$. Suppose you were calculating a scattered wave that arrives at time $t = T + \tau$, what would happen? The answer is you would see this scattered wave arrive at time $t = \tau$, assuming $\tau < T$. This wrong arrival time occurs often when waves get trapped due to strong multiple scattering.