juliawavescattering / multiplescattering.jl

A Julia library for simulating, processing, and plotting multiple scattering of waves.

MultipleScattering

A Julia library for simulating, processing, and plotting multiple scattering of waves.

The library focuses on multipole methods (addition translation theorems) to solve the inhomogeneous Helmholtz equation (time-harmonic waves). Multipole methods are particularly efficient at solving scattering from particles in an infinite domain. This library is configured to use T-matrices (also known as scattering matrices) to represent scattering from particles with any shape and properties (currently implemented for acoustics).

The package is setup to deal with different spatial dimensions and types of waves which satisfy Helmholtz equation's, e.g. acoustics, electromagnetism, elasticity. For details on some of the maths see Martin (1995) and Gower et al. (2017).

Installation

This package is available for Julia 1.0.5 and beyond. To get started, just add the package by typing

julia> ]
pkg> add MultipleScattering

then press the backspace key followed by

julia> using MultipleScattering

Documentation

• Latest Documentation — documentation of the in-development version.

Simple example

Define the properties of your host medium:

dimension = 2 # could also be 3, but then all 2D vectors need to be 3D vectors
host_medium = Acoustic(dimension; ρ = 1.0, c = 1.0); # 2D acoustic medium with density ρ = 1.0 and soundspeed c = 1.0

An acoustic medium in 2D with density 1 and wavespeed 1.

Next, define two dense, circular acoustic particles, the first centred at [-2,2] with radius 2 and the second at [-2,-2] with radius 0.5,

particle_medium =  Acoustic(dimension; ρ = 10.0, c = 2.0); # 2D acoustic particle with density ρ = 10.0 and soundspeed c = 2.0
p1 = Particle(particle_medium, Sphere([-2.0,2.0], 2.0));
p2 = Particle(particle_medium, Sphere([-2.0,-2.0], 0.5));
particles = [p1,p2];

Lastly we define the source, for example an incident plane wave () using a helper function.

source = plane_source(host_medium; direction = [1.0,0.0])

Once we have these three components, we can build our FrequencySimulation object

simulation = FrequencySimulation(particles, source)

To get numerical results, we run our simulation for specific positions and angular frequencies,

x = [[-10.0,0.0], [0.0,0.0]]
max_ω = 1.0
ω = 0.01:0.01:max_ω
result = run(simulation, x, ω)

Plot

The package also provides recipes to be used with the Plots package for plotting simulations after they have been run.

In our above simulation we ran the simulation for 100 different wavenumbers, and measured the response at the location (-10,0).

To plot the time-harmonic response across these wavenumbers type:

using Plots
plot(result)

For a better overview you can plot the whole field in space for a specific angular frequency by typing:

ω = 0.8
plot(simulation,ω)

This shows the field over a spatial domain for one particular angular frequency ω.

Note: most things in the package can be plotted by typing plot(thing) if you need an insight into a specific part of your simulation.

To calculate an incident plane wave pulse in time use:

time_result = frequency_to_time(result)
plot(time_result)

Or for a Gaussian impulse in time:

t_vec = LinRange(0.,700.,400)
time_result = frequency_to_time(result; t_vec = t_vec, impulse = GaussianImpulse(max_ω))
plot(time_result)

Examples

One way to learn the different features of this package is through the following examples.

Examples
Helmholtz resonator	Random particles
Choosing basis order	Particles in a circle
Times response lens	Statistical moments

Acknowledgements and contributing

This library was originally restructured from one written by Artur L Gower and Jonathan Deakin.

Contributions are welcome, and the aim is to expand this package to elasticity and electromagnetism.

multiplescattering.jl's Issues

Convergence example

We need an example showing that the response converges when we increase the number of Hankel functions and/or the box size. This may or may not lead to a good test, but it will be a good qualitative validation.

plot fields inside particles

The package already calculates the field inside the particles. So we should plot this field! Only graphical issue is that if the scatterers are impenetrable, then the fields inside will return 0 , which could make visualisation tricky.

I suggest we draw a grey outline around the particles, so that we can continue to see them.

Adding this feature will allow to simulate waveguides!

Rename Moments type

It's odd that Moment has a subtype which is called moments. It would make more sense to call it label+moments or moments with metadata. Another idea would be to borrow the word from Summary Statistics and call it a ModelSummary.

Fix time simulations

Some thought needs to go into the outer structures for the TimeSimulations again we can obviously reuse a lot of the code which does the mathsy heavy lifting, but we need some motivating examples for the structs. We also need to think about how the TimeSimulations will interact with the random machinery.

Remove parameter T from most all parametric types

The type T, used to indicate the type of numbers used, i.e. Float64, FLoat32, etc.. appear in many parametric types:

Shape{T<:AbstractFloat,Dim}
Particle{T<:AbstractFloat,Dim,P<:PhysicalMedium,S<:Shape}

and many more... Keeping track of T is pointless, as we do not dispatch based on this in any circumstance. It only currently helps enforce that a set of types all use the same T. Enforcing this has almost no purpose, and is at times annoying.

This T should be removed from almost all parametric types.

Needs pictures!

The repo readme needs some pictures of plots. Each example also needs pictures showing the results from when it is run. Should we put each example in its own folder with a readme? Or should each example be a Jupyter notebook?

RFC: Include example of vector waves (such as electromagnetism or elasticity)

Currently the machinery should work for acoustics, which is a scalar equation. But we have not yet tested for any physics with a FieldDim > 1. Without a test case it is hard to see what to do. The outer types and structure should be fine, but I expect the inner machinery will break. We also need to be careful that when we generalise the machinery fully that we don't make it slower.

Avoid strange parametric types

In FrequencySimulation, the field particles is defined as:

 particles::Vector{Particle{Dim,PP,S,T} where {PP<:PhysicalProperties,S<:Shape}}

but I think this leads to unintended behaviour:

    sound = Acoustic(0.1,0.1 + 0.0im,2)
    particles = [Particle(sound,Circle([x,0.0], .5)) for x =0.:2.:10.]

    #As expected
    Particles = Vector{Particle{D,P,S,T}} where {D,P<:PhysicalProperties,S<:Shape,T<:AbstractFloat}
    typeof(particles) <: Particles #true 

# However, when calling the constructor of `FrequencySimulation` it will perform the following type cast
    ParticleVec = Array{Particle{2,PP,S,Float64} where {PP<:PhysicalProperties,S<:Shape}}
    particles = ParticleVec(particles)

    # Which leads to:
    typeof(particles) <: Particles # False

This in turn makes it difficult to write functions for sim.particles. Instead we should use:

    particles = [Particle(sound,Circle([x,0.0], .5)) for x =0.:1.:2.]
    ParticleVec = Vector{Particle{2,PP,S,Float64}} where {PP<:PhysicalProperties,S<:Shape}
    particles = ParticleVec(particles)
    typeof(particles) <: Particles # true

Change the type of `FrequencySimulationResult.ω` from `RowVector` in to just `Vector`?

As of Julia 1.0, RowVector is no longer a thing. See this issue.

Instead there is a type Transpose, Adjoint, or just Array{2}. All in all, I think there is plenty of room for confusion now, and the intuitive idea of having a "lying down" vector is lost. So I suggest just using Vector for the type FrequencySimulationResult.ω.

Source.coef should return a vector

Similar to #27, Source.coef should be a function which returns a vector of the coefficients used for a regular basis function. This will lead to greater generalisation for different physics and more efficient code.

Allow the amplitude of a source to be a function

At the moment, the sources implement in physics/acoustics/source.jl take the follow arguments:
point_source(medium, source_position, amplitude), where amplitude is a constant.

Allowing amplitude to be a function of frequency ω has many uses: shift the source in time, allow for one source to use, say, a Gaussian filter, while another source can still use a sinc filter. There are almost no drawbacks to this generalisation, so I will implement.

Margins in field plotting

So at the moment plotting a field seems to only work for the least bounding box of the container and the source, so if I try to set the xlims or ylims beyond a spherical container's radius it just shows white:

It would be nice to see what the field looks like on the other side of the container as well.

I'm guessing an option would be to make separate shapes for domain=Rectangle(), then container=Shape(), make particles=random_particles(shape=container), and then do a simulation with FrequencySimulation(particles, shape=domain), but some kind of margin around a container might be a nice addition, and easier to work with? Or some kind of way to attach a container shape to a domain shape, so the user wouldn't need to keep track of so many things?

Frequency or wavenumber

We should really use one of frequency or wavenumber almost everywhere. Although wavenumber is important, I suspect frequency (angular or not?) is more natural to use. We already call the main type a FrequencyModel, even though we parameterise by wavenumber. Given that wavenumber (in the medium) is still important to the validity, we could always print the wavenumber range when calculating the model. What do you think?

Implement and test a physics type with dim != 2

The capability exists but we have not yet implemented acoustics (or any other physics) in a space of dimension == 2. 1D would be easy, although perhaps useless; 3D looks much more complicated and may raise some errors however.

Specialise function frequency_to_time for real time signal

We should assume all time signals are real and specialise the inverse fourier for that scenario. In which case we can use only positive frequencies, but need to take the real part and multiple by two:

Note that any frequency ωi not given is the same as taking F(ωi) =0. For very large ωi, we expect F(ωi) to decay and therefore become less significant, but this is not true for ωi =0. If F(0) is not given, one simple solution is to interpolate F(0) ~= (ω2 F(ω1) - ω1 F(ω2))/(ω2 - ω1).

Thoughts?

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

An easy to use type impulse

Currently the struct impulse is very elegant: it has two fields, which are functions representing the impulse in time and frequency (one being the analytic Fourier Transform of the other). However, I think the struct should be easy to use and understand for the most common cases of use. These are:

1. simple analytic functions like a discrete delta and Gaussian
2. providing a discrete impulse in either time (most common) or frequency. This is always the case when comparing with experiments.

In both cases, it should be easy to plot both the impulse in time and frequency, even if the user only provided the impulse in, say, time. Currently I think it is not easy to provide 2, and even the inability to plot the discrete delta function in 1 will confuse most. For example, after sampling the discrete delta in frequency, it can be useful to look at the resulting impulse in time, and vice-versa (a common exercise in signal processing). It helps decide how to sample the frequency.

I think it will make our lives easier to use instead

struct Impulse{T<:AbstractFloat}
    in_time::Vector{T}
    in_freq:::Vector{Complex{T}}
    t::Vector{T}
    ω::Vector{T}
end

One draw back is that it seems that both the vectors ω and t need to be supplied. But this is not really an issue in practice, and, for example, the impulse in time can be easily re-sampled from the impulse in frequency for a different t vector.

Although the struct impulse is elegant, and is more flexible than above, there is no clear example where this flexibility is needed. And I think it is not worth the added complications of trying to convert impulse vectors to functions.

Did you method consider the absorption boundaries?

Hi,

Nice work!

Curious that if your method has considered the absorption boundaries, such as the "perfect matched layer".

For my application, I would be more interested in the ones with open boundaries.

Would like to hear your kind reply!

Thanks
Guangyuan

FrequencySimulationResult, FrequencySimulation, and basis order

This is just to debate out current types.

Am not sure how useful is the type FrequencySimulation now. Essentially it is not just an added step were the user always has to call:

sim = FrequencySimulation(particles, source)
result = run(sim, x, ω)

Why not just have:

result = run(particles, source, x, ω)

how is FrequencySimulation useful?

Also for debate, I think FrequencySimulationResult should have a field basis_order::Vector{Int}. I think by having the field basis_order it makes it easier for new users to discover this option, rather than just being keyword of run. Second, even if the used doesn't specify basis_order it can still be very important to know and keep a record of what the basis order was. On the other hand, if FrequencySimulationResult didn't have the field basis_order that would allow us more freedom to, for example, have a different basis order for each particle, depending on it's radius....

Allow for any smooth particle

I think we should: use the T-matrix to allow for any (smooth) particle.
For some details on the maths, see T-matrix.pdf, which also shows how to use a generic incident wave #10.

For the particle shall we just:

abstract Particle

type CircularParticle{T <: AbstractFloat} <: Particle
    x::Vector{T} #position of the centre
    r::T # radius
    c::Complex{T} # phase wave speed inside the particle
    ρ::T # the particles density
end

function tmatrix{T}(sim::FrequencySimulation{T}, p::CircularRectangle{T}, k::T, hankel_order::Int, exciting_hankel_order::Int=hankel_order)
    if hankel_order == exciting_hankel_order
        Zn{T}(sim, p, k, hankel_order)
   else 0
   end
end

Do you envisage any efficiency problems with this?

Change Plots submodule to use RecipesBase

Because the Plots package is a dependency of MultipleScattering.Plots, we cannot precompile the whole MultipleScattering package even though it would make loadup much faster. Instead of using the Plots package, we are really supposed to use the RecipesBase which adds recipes (or methods) to the plots function. It would also allow us to just call plot(model) and the recipe would be used.

Recipe docs

Basis function should return a vector

Currently the outgoing_basis_function and regular_basis_function return a function, say basis, where basis(m, x) is the basis function evaluated at the order m and position x. This is awkward because the different ways to add basis(m, x) over m depends on the physics and spatial dimension. Things need to be rewritten so that, for example, outgoing_basis_function should return either a vector evaluated over all m smaller than a certain order, or return a function basis, where basis(max_order, x) is a vector.

For current code examples see: general and acoustic regular waves.

Reintegrate old tests

On the version 0.2 branch, there are a lot of old tests lying around unused. We need to adapt them to the new structures and get it tested!

Allow different incident waves

Clearly having the option to use a point source (or collection of point sources) or an incident plane wave, will greatly increase the usefulness of this library. So we need a type for all Incident waves.

The essential data for this type is a scattering function (x,y) -> scattering coefficients, i.e. given any point returns the scattering coefficients of a Bessel J expansion. One major issue is that we are not allowed to save functions. So we can not have this scattering function as a field of, say, FrequencySimulation. This is why I removed the impulse function and replaced with an impulse vector for the Fourier transforms. One possible solution is to allow for only plane waves or sum of point sources, then by recording which one was used, we can then recreate scattering function.

Simple shape functions not counted in coverage when precompiled

When calling name on a shape, like this:

name(shape::Circle) = "Circle"

coverage says that it is not covered but it is called in tests. I also have the same problem when using coverage locally. However, we did not have this problem before we enabled precompile.

I think this is a Julia bug, and so should be reported there, but I'm a bit nervous that I'm wrong...

Time of flight shape is wrong

Time of flight has been written as constant time of flight from a point source, however it was intended (from the old package) to be the time of flight from a planar wave (noticed by Art). This means that all the functions related to it are wrong. Suggesting moving current shape to type called TimeOfFlightFromPoint and making TimeOfFlight the planar version (to make it work with old data). Could perhaps change this to TimeOfFlightFromPlanar in a later release to avoid ambiguity.

Simple interface for plot field for both time and frequency

I think the following should work:

• plot(SimulationResult) should return several lines plots for each x
• plot(SimulationResult, linetype=:contour) should attempt a contour plot. This should work for both TimeSimulationResult and FrequencySimulationResult. Note the function run(Simulation,shape) can return both TimeSimulationResult and FrequencySimulationResult, which can then be plotted
• plot(FrequencySimulation, ω) should generate a contour plot.

There might be a dispatch issue for the first two, as I'm struggling to get plot recipes to do this...

suggestion: meshes.jl interop

that would really increase the number of supported shapes and would be in line with what seems to be the most supported shape library in julia

Boundary conditions test fails randomly

The boundary conditions test uses randomly generated particles, and this sometimes cause the test to fail (example). We need to make the test more permissive (but still strong), or alternatively use a seed which is known to pass.

Error while running example: ERROR: UndefVarError: bounding_rectangle not defined

Issue

The "Simple Example" at https://github.com/JuliaWaveScattering/MultipleScattering.jl doesn't work.

Procedure

Run the "Simple Example":

using MultipleScattering
using Plots

# 2D acoustic medium with density ρ = 1.0 and soundspeed c = 1.0
host_medium = Acoustic(2; ρ=1.0, c=1.0);  

# 2D acoustic particle with density ρ = 10.0 and soundspeed c = 2.0
particle_medium =  Acoustic(2; ρ=10.0, c=2.0); 

p1 = Particle(particle_medium, Sphere([-2.0,2.0], 2.0));
p2 = Particle(particle_medium, Sphere([-2.0,-2.0], 0.5));
particles = [p1,p2];

source = plane_source(host_medium; direction = [1.0,0.0])

simulation = FrequencySimulation(particles, source)

x = [[-10.0,0.0], [0.0,0.0]]
max_ω = 1.0
ω = 0.01:0.01:max_ω
result = run(simulation, x, ω)

time_result = frequency_to_time(result)

# plot(time_result)
plot(simulation,0.8)

gui()

Expected

No error

What I get:

ERROR: UndefVarError: bounding_rectangle not defined
Stacktrace:
 [1] macro expansion
   @ ~/.julia/packages/MultipleScattering/mcGmd/src/plot/plot.jl:62 [inlined]
 [2] apply_recipe(plotattributes::AbstractDict{Symbol, Any}, sim::FrequencySimulation{Float64, 2, Acoustic{Float64, 2}}, ω::Float64)
   @ MultipleScattering ~/.julia/packages/RecipesBase/92zOw/src/RecipesBase.jl:282
 [3] _process_userrecipes!(plt::Any, plotattributes::Any, args::Any)
   @ RecipesPipeline ~/.julia/packages/RecipesPipeline/BY2Dd/src/user_recipe.jl:36
 [4] recipe_pipeline!(plt::Any, plotattributes::Any, args::Any)
   @ RecipesPipeline ~/.julia/packages/RecipesPipeline/BY2Dd/src/RecipesPipeline.jl:70
 [5] _plot!(plt::Plots.Plot, plotattributes::Any, args::Any)
   @ Plots ~/.julia/packages/Plots/SVksJ/src/plot.jl:172
 [6] plot(::Any, ::Vararg{Any, N} where N; kw::Any)
   @ Plots ~/.julia/packages/Plots/SVksJ/src/plot.jl:58
 [7] plot(::Any, ::Any)
   @ Plots ~/.julia/packages/Plots/SVksJ/src/plot.jl:52
 [8] top-level scope
   @ REPL[51]:1

Environment

Julia version: Version 1.6.2

Project.toml

[deps]
MultipleScattering = "80a8ab25-5750-5d93-a6d7-4adc97cdd5fb"

Manifest.toml

# This file is machine-generated - editing it directly is not advised

[[Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.1"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[GSL]]
deps = ["GSL_jll", "Libdl", "Markdown"]
git-tree-sha1 = "3ebd07d519f5ec318d5bc1b4971e2472e14bd1f0"
uuid = "92c85e6c-cbff-5e0c-80f7-495c94daaecd"
version = "1.0.1"

[[GSL_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "9a6d2b6c33aa7061654d2f1b31d9f7e1c2bd2dbc"
uuid = "1b77fbbe-d8ee-58f0-85f9-836ddc23a7a4"
version = "2.7.1+0"

[[HalfIntegers]]
git-tree-sha1 = "8aeed192be8b368c9f5c149806db90f4451c8bb1"
uuid = "f0d1745a-41c9-11e9-1dd9-e5d34d218721"
version = "1.3.3"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[MultipleScattering]]
deps = ["GSL", "LinearAlgebra", "OffsetArrays", "Printf", "ProgressMeter", "Random", "RecipesBase", "SpecialFunctions", "StaticArrays", "Statistics", "WignerSymbols"]
git-tree-sha1 = "505fbc04c05ebbb01baf06deb5024a871936470b"
uuid = "80a8ab25-5750-5d93-a6d7-4adc97cdd5fb"
version = "0.1.8"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OffsetArrays]]
deps = ["Adapt"]
git-tree-sha1 = "c0f4a4836e5f3e0763243b8324200af6d0e0f90c"
uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
version = "1.10.5"

[[OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[Primes]]
deps = ["Test"]
git-tree-sha1 = "ff1a2323cb468ec5f201838fcbe3c232266b1f95"
uuid = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
version = "0.4.0"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[ProgressMeter]]
deps = ["Distributed", "Printf"]
git-tree-sha1 = "afadeba63d90ff223a6a48d2009434ecee2ec9e8"
uuid = "92933f4c-e287-5a05-a399-4b506db050ca"
version = "1.7.1"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[RationalRoots]]
git-tree-sha1 = "4165f49c451e6fc4bc0a2eeb1250abdf8c6e550e"
uuid = "308eb6b3-cc68-5ff3-9e97-c3c4da4fa681"
version = "0.1.0"

[[RecipesBase]]
git-tree-sha1 = "b3fb709f3c97bfc6e948be68beeecb55a0b340ae"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.1.1"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[SpecialFunctions]]
deps = ["OpenSpecFun_jll"]
git-tree-sha1 = "d8d8b8a9f4119829410ecd706da4cc8594a1e020"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "0.10.3"

[[StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "da4cf579416c81994afd6322365d00916c79b8ae"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "0.12.5"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[WignerSymbols]]
deps = ["HalfIntegers", "Primes", "RationalRoots"]
git-tree-sha1 = "dfb49a51c4c8e205d1d7ddc7e6a1a29da871ef52"
uuid = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b"
version = "1.1.0"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"

Numerical instability related to hankel order

This is accentuated when particles are small and there is at least two particles far apart. This primarily affects the internal waves, which makes verifying the boundary conditions difficult.

This does not happen in the master version! The way that internal waves are calculated are identical in master and version-0.2, so I think it is related to solving the scattering_matrix.

import StaticArrays: SVector
using MultipleScattering

ω = 0.1

medium = Acoustic(1.,1.,2)
sound = Acoustic(medium.ρ, 4. + 0.0im,2)
p1 = Particle(sound,Circle([-5.0,0.0], .2))
p2 = Particle(sound,Circle([5.0,0.0], .2))

source = TwoDimAcousticPlanarSource(medium, SVector(0.0,0.0), SVector(1.0,0.0), 1.)

particles = [p1, p2]
sim = FrequencySimulation(medium, particles, source)

widths = -5.4:0.01:-4.4
x_vec = [ SVector(x,.0) for x in widths]
result = run(sim, ω, x_vec; basis_order = 7)

using Plots; pyplot()
plot(widths, abs.(field(result)[:]), ylim = (0.9,2.))

But by lowering the Hankel order by one

result = run(sim, ω, x_vec; basis_order = 6)
plot(widths, abs.(field(result)[:]))

Note that as the density of the particle is the same as the background, the field should be smooth everywhere. This instability goes away if we increase the frequency or bring the particles closer together. However, it does not go away by just adding more particle between them. For example:

particles = [Particle(sound,Circle([x,0.0], .2)) for x = -5.:0.6:5 ]
sim = FrequencySimulation(medium, particles, source)

widths = -5.4:0.01:0.4
x_vec = [ SVector(x,.0) for x in widths]
result = run(sim, ω, x_vec; basis_order = 7)
plot(widths, abs.(field(result)[:]))

Add docstring with signatures

Quite a lot of functions are missing docstrings, and most with docstrings are missing the function signatures. See Julia documentation for recommendations on how to write docstring.

Source should have a field medium

At present Source only has the fields coef and field. This leads to some problems, first it means that the constructor is used in the form FrequencySimulation(medium, particles, source), because Source has forgotten its medium. This then also allows us to do nonsense like this:

source_medium = Acoustic(2.0, 2.0, 2)
medium = Acoustic(1.0, 1.0, 2)
source = plane_source(medium)
simulation = FrequencySimulation(medium, source)

If we had particles, then this would lead to non-sense results. The constructor FrequencySimulation should never have the freedom to choose a medium which is different from the Source.

Fix random simulations

The ability to do random simulations has not yet been finished in version 0.2. Many of the inner functions can be unchanged, but the outer structures (SimulationDistribution and StatisticalMoments) need some more careful thought and motivating examples.