# Scattering from a sphere

Here we show how to calculate the average acoustic wave scattered from a sphere filled with particles. We will the analytic solution with an approximation using effective low frequency properties.

## Choose the microstructure

```
using EffectiveWaves, LinearAlgebra
spatial_dim = 3
medium = Acoustic(spatial_dim; ρ=1.0, c=1.0)
s1 = Specie(
Acoustic(spatial_dim; ρ=10.0, c=10.0), Sphere(spatial_dim, 1.0);
volume_fraction=0.15
);
s2 = Specie(
Acoustic(spatial_dim; ρ=0.22, c=0.22), Sphere(spatial_dim, 1.0);
volume_fraction=0.05
);
species = [s1,s2]
# calculate the low frequency effective properties
eff_medium = effective_medium(medium, species)
# choose incident angular frequency
ω = 0.22
# choose maximum order of the spherical harmonic for scattered waves from particles
basis_order = 2
# choose a numerical tolerance
tol = 1e-7
# calculate the effective wavenumbers which are valid for any frequency low frequency effective properties
opts = Dict(
:num_wavenumbers => 2, # calculate at least 3 wavenumbers
:basis_order => basis_order,
:tol => tol
)
# note this is the heaviest calculation in this example
k_effs = wavenumbers(ω, medium, species; opts...)
```

## Is one effective wavenumber enough?

To easily calculate scattered waves we need `k_effs[1]`

to have a smaller imaginary part than the other wavenumbers `k_effs[2:end]`

. This is easily checked with a scatter plot

```
using Plots
pyplot(size = (600, 350))
scatter(k_effs, lab = "kp", ylims = (0.0,Inf))
# savefig("../assets/sphere-wavenumbers.png")
```

When using only one wavenumber, and excluding the other wavenumber, will lead to an error that is roughly given by

```
λ = 2π * medium.c / ω
approximate_error = abs(exp(im * λ * k_effs[2]) / exp(im * λ * k_effs[1]))
# output
1.5188e-6 # approximately
```

The precise error depends on the application.

## Define the source and material

As we want to calculate effective scattered waves, and not just effective wavenumbers, we need to define a source and material domain. We will choose a plane-wave source, and a sphere for the domain.

```
# choose an incident wave
# NOTE: having the incident wave direction along the z-axis leads to azimuthal symmetry, which simplifies the next steps.
θ = 0.0
psource = PlaneSource(medium, [sin(θ),0.0,cos(θ)]);
# choose the size and position of the spherical domain of the material
R = 20.0
material = Material(medium,Sphere(spatial_dim, R),species);
# the larger the material region, the larger this basis_field_order needs to be
basis_field_order = 7
```

## Calculate the scattered waves

Below we calculate the outward going spherical wave coefficients. These coefficients, multiplied with a basis of outward going waves, will give the scattered field.

```
# create a sphere using the low frequency effective properties
effective_sphere = Particle(eff_medium, material.shape)
region_order = basis_order + basis_field_order
# calculate the T-matrix, also know as scattering matrix, and take the diagonal
t_mat = t_mat(effective_sphere, medium, ω, region_order);
# calculate the outward going spherical wave coefficients
source_coefficients = regular_spherical_coefficients(psource)(region_order,zeros(3),ω);
scat_low_coefficients = t_mat * source_coefficients
# for the analytic solution, we need the wavemode first
wave = WaveMode(ω, k_effs[1], psource, material;
tol=tol,
basis_order = basis_order,
basis_field_order = basis_field_order
)
# and then the outward spherical wave coefficients
scat_coefficients = material_scattering_coefficients(wave)
```

## Calculate the field

```
xlim = 3.0 * R
ylim = 4.4 * R
# a region to plot
region = Box([[-xlim,-ylim], [xlim,ylim]]);
# a quick fix to make it a 3D region
xs, xinds = points_in_shape(region; res = 180);
xs = [[x[1],0,x[2]] for x in xs];
# Calculate the field for the low frequency effective properties
particles = [Particle(eff_medium,material.shape)]
sim = FrequencySimulation(particles, psource)
fs = field(sim, ω, xs[xinds], scat_low_coefficients);
field_mat = zeros(typeof(fs[1]), length(xs), 1);
field_mat[xinds,:] = fs;
result_low = FrequencySimulationResult(field_mat, xs, [ω])
# Calculate the analytic scattered field
fs = field(sim, ω, xs[xinds], scat_coefficients);
field_mat = zeros(typeof(fs[1]), length(xs), 1);
field_mat[xinds,:] = fs;
result = FrequencySimulationResult(field_mat, xs, [ω])
```

## Plot the field

```
pyplot( frame = :axes,
size = (xlim * 5.0, ylim * 4.0),
flip = false,
xlims = :auto, ylims = :auto,
title="", xlab="ka",ylab="y"
)
plot(result_low, ω, seriestype = :contour, field_apply=abs
, color = :balance # for julia 1.1 change to color = :pu_or
, levels = 16, clims=(0.0,1.65)
, title = "abs low freq. scat"
, axis=false)
plot!(material.shape)
p_scat = plot!(leg = false, frame = :none
,axis=false, xlab="",ylab="")
# savefig("../assets/low-freq-scat-sphere.png")
```

```
plot(result,ω, seriestype = :contour, field_apply=abs
, levels = 16, clims=(0.0,1.65)
, title = "abs scattered field"
, color = :balance
, colorbar = false
)
plot!(material.shape)
plot!(leg = false
, axis=false, xlab="",ylab="")
# savefig("../assets/scat-sphere.png")
```

# Plot a gif

Finally, for fun, we plot a GIF of the real part of the full field.

```
max_c = maximum(abs.(field(result)))
min_c = - max_c
ts = LinRange(0.,2pi/ω,30)
anim = @animate for t in ts
plot(result, ω, seriestype=:contour
, levels = 20
, phase_time=t
, leg = false
, field_apply = real
, clim = (min_c, max_c)
, colorbar = false
, color = :balance
, frame = :none
, title = "Re. full field"
)
plot!(material.shape)
plot!(leg = false, frame = :none
, axis=false, xlab="",ylab="")
end
# gif(anim,"../assets/real-sphere-scatter.gif", fps = 7)
anim = @animate for t in ts
plot(result_low, ω, seriestype=:contour
, levels = 20
, phase_time=t
, leg = false
, field_apply = real
, clim = (min_c, max_c)
, colorbar = false
, color = :balance
, frame = :none
, title = "Re. low freq full field"
)
plot!(material.shape)
plot!(leg = false, frame = :none
, axis=false, xlab="",ylab="")
end
# gif(anim,"../assets/real-sphere-low-scat-scatter.gif", fps = 7)
```

# A sphere with radial symmetry

The simplest possible scenario with a finite number of particles is radial symmetry. That is, to excite a sphere filled with uniformly distributed particles with a radially symmetric incident field.

At first it may appear that the scenario of radial symmetry is a very specific case and perhaps not that useful. However, the paper ? shows that the same effective wavenumbers appear for both the radial symmetry case and all other cases including plane wave incidence on a plate.

Below we show two different methods: using just one effective wavenumber and another which solves a 1D integral equation which includes all effective wavenumbers

## Effective wave with radial symmetry

Let us first choose the material parameters again

```
using EffectiveWaves
ω = 0.1
medium = Acoustic(3; ρ=1.0, c=1.0);
# Choose the particle properties
particle_medium = Acoustic(3; ρ=8.0, c=8.0);
r = 1.0;
s1 = Specie(
particle_medium, Sphere(r);
volume_fraction = 0.1
)
k = ω / medium.c;
ka = k * r;
# here we use the Percus-Yevick approximation for the inter particle pair-correlation
py_pair = PercusYevick(3; meshsize = 0.12, maxlength = 30, rtol = 1e-4)
microstructure = Microstructure(medium, [s1], py_pair)
# Choose the material properties
R = 5.0
material = Material(Sphere(R),microstructure);
Symmetry(material,source) == RadialSymmetry{3}()
# Choose a radially symmetric source
source = regular_spherical_source(medium, [1.0+0.0im];
position = [0.0,0.0,0.0],
symmetry = RadialSymmetry{3}()
);
# Calculate the effective wavenumbers and wavemode
basis_order = Int(round(2.0 * ka)) .+ 1
keffs = wavenumbers(ω, microstructure;
basis_order = basis_order,
symmetry = PlanarAzimuthalSymmetry()
)
wave = WaveMode(ω, keffs[1], source, material;
basis_order = basis_order
)
# Calculate the scattering coefficients from the whole material
scattering_coefficients = material_scattering_coefficients(wave);
# Due to radial symmetry scattering_coefficients has only one component
# Calculate the far-field patterns
Ys = spherical_harmonics(0, 0.0, 0.0);
u∞ = sum((1/k) .* Ys[1] .* scattering_coefficients[1] .* exp.(-(pi*im/2)))
```

## Integral method with radial symmetry

Below shows how to use the same parameters defined above to solve a 1D integral equation.

```
# The integral methods requires a number to determine how fine the mesh should be
basis_field_order = Int(round(0.5 * k .* R)) .+ 1
basis_field_order = max(2,basis_field_order)
# This depends on the pair-correlation chosen. It depends how many Fourier series to use to represent the pair-correlation
polynomial_order = 10;
# This solves the integral equation and may take a few minutes
discretewave = discrete_system_radial(ω, source, material, Symmetry(source,material);
basis_order = basis_order,
basis_field_order = basis_field_order,
polynomial_order = polynomial_order
);
scattering_coefficients2 = material_scattering_coefficients(discretewave)
u∞2 = (1/k) * Ys[1] * scattering_coefficients2[1] * exp(-(pi*im/2))
abs(u∞ - u∞2) / abs(u∞) < 1e-3
```

## Full Integral method with radial symmetry

Below is a method which takes the integral equation and does not use radial symmetry to simplify it. This removes the Gibb's phenomena which comes with using a Fourier series to approximate the pair-correlation, but requires more mesh points to reach the same level of precision.

```
rtol = 1e-2; maxevals = Int(1e5);
dwave = discrete_system(ω, source, material;
basis_order = basis_order,
basis_field_order = basis_field_order,
rtol = rtol, maxevals = maxevals,
);
scattering_coefficients_discrete = material_scattering_coefficients(dwave)
u∞_discrete = (1/k) * Ys[1] * scattering_coefficients_discrete[1] * exp(-(pi*im/2))
abs(u∞ - u∞_discrete) / abs(u∞) < 1e-3
abs(u∞2 - u∞_discrete) / abs(u∞) < 3e-4
```