# Pair correlation

In the previous example there was no mention on how the particles are distributed. This can be specified by choosing a pair correlation, and when no pair correlation is specified, the default is to assume that particles are uncorrelated, except they can not overlap.

Choosing the particle distribution only affects the effective wavenumbers and wavemodes, see references [1,2]

We note that the definition of the pair-correlation $g$, for a finite number of multi-species particles, is [1]

\[g(\mathbf r_1, \lambda_1; \mathbf r_2, \lambda_2) = \frac{p(\mathbf r_1, \lambda_1; \mathbf r_2, \lambda_2)}{p(\mathbf r_1, \lambda_1)p(\mathbf r_2, \lambda_2)}\frac{J -1}{J}\]

where $\mathbf r_j$ is the vector position of the centre of particle-$j$, $\lambda_j$ represents the size or other distinguishing properties of the type of particle, $p$ is the probability density function, and $J$ is the total number of particles. Note that when the particles become uncorrelated in the limit $|\mathbf r_1 - \mathbf r_2| \to \infty$ we have that

\[\lim_{|\mathbf r_1 - \mathbf r_2| \to \infty} g(\mathbf r_1, \lambda_1; \mathbf r_2, \lambda_2) = 1\]

## Bespoke pair-correlation

Here is an example of choosing your own pair-correlation for a material filled with only one type of particle

```
medium = Acoustic(3; ρ=1.2, c=1.0)
# Choose the species
r = 0.5
s = Specie(
Acoustic(3; ρ = 0.01, c = 0.01),
Sphere(r),
volume_fraction = 0.3,
separation_ratio = 1.0
);
```

Next we create a microstructure that has only this species, and has a specific pair-correlation

```
r = 1.0:0.1:10.0
my_pair_correlation = 1.0 .+ 0.2 .* sin.(r) ./ r.^3
Dim = 3
dpc = DiscretePairCorrelation(Dim, r |> collect, my_pair_correlation)
micro = Microstructure(medium, s, dpc);
```

Note that when specifying a pair-correlation, the minimal distance between particles will be taken to be `dpc.r[1]`

. This is stored in `dpc.minimal_distance`

. Previously when defining the `Specie`

we specified `separation_ratio = 1.01`

, which means the minimal distance between particles centres' is `2 * separation_ratio * r`

, this is used when no pair-correlation is specified, otherwise the value given in `dpc.minimal_distance`

will be used. In the future, we will phase out the use of `separation_ratio`

in the Specie.

## Percus-Yevick

Let us consider a material filled with only one type of particle and use the Percus-Yevick approximation to calculate the pair-correlation for 2D hard spheres. That is, sphere which do not attract of repel each other. For details see Notes on Percus-Yevick [2].

```
pair_type = PercusYevick(3; rtol = 1e-2, maxlength = 200)
micro = Microstructure(medium,s, pair_type);
# output
```

When calling the above, the pair correlation is calculated and stored in `micro.paircorrelations`

which is a matrix of the type `DiscretePairCorrelation`

. This type assume that the pair-correlation is only a function of the distance between particles, and the types of particles, like their size.

We can plot the result of the Percus-Yevick approximation with the package Plots:

```
using Plots
plot(micro.paircorrelations[1].r, micro.paircorrelations[1].g,
xlab = "distance", ylab = "P-Y"
)
```

which we can compare with Figure 8.3.1 from [1] below.

Note that for $x < 1$ the two particles of radius 0.5 would overlap, so the pair correlation should be zero. Also note that `dp`

is the variation from uncorrelated, which is why we add 1.0 to get the pair correlation.

## Calculate an effective wavenumber

The more points sampled within the pair correlation the longer it will take to calculate the effective wavenumber.

First we calculate the wavenumbers with the simplest pair correlation (hole correction), and then compare the results with Percus-Yevick.

```
micro = Microstructure(medium,s);
ω = 0.4
kps = wavenumbers(ω, micro;
basis_order = 1, num_wavenumbers = 4
)
pair_type = PercusYevick(3; meshsize = 0.1, maxlength = 50)
micro = Microstructure(medium, s, pair_type);
kps2 = wavenumbers(ω, micro;
basis_order = 1, num_wavenumbers = 4
)
```

we can then compare how the wavenumbers change with Percus-Yevick with a scatter plot

```
using Plots
scatter(kps, lab = "Hole correction")
scatter!(kps2, lab = "Percus-Yevick")
```

We can see that in this case, the effective wavenumbers with (`kps2`

) and without (`kps`

) Percus Yevick are quite different.

## Calculate reflection

To calculate the average reflection, or scattering, from a material is the same as before, except we just need to replace `Species`

with `Microstructure`

. For example to calculate reflection from a plate:

```
k_eff = kps2[1]
normal = [0.0,0.0,-1.0] # an outward normal to both surfaces of the
width = 150.0 # plate width
origin = [0.0,0.0,width/2] # the centre of the plate
plate = Plate(normal,width,origin)
# note below we use micro instead of species
material = Material(plate, micro)
source = PlaneSource(medium, [0.0,0.0,-1.0])
# Calculate the wavemode for the first wavenumber
# the WaveMode function calculates the types of waves and solves the needed boundary conditions
wavemodes = WaveMode(ω, k_eff, source, material; tol = 1e-6, basis_order = 1);
RTeff = reflection_transmission_coefficients(wavemodes, source, material)
abs.(RTeff)
```

We can compare the result of not using Percus-Yevick below.

```
k_eff = kps[1]
micro = Microstructure(medium,s);
material = Material(plate, micro)
wavemodes = WaveMode(ω, k_eff, source, material; tol = 1e-6, basis_order = 1);
RTeff = reflection_transmission_coefficients(wavemodes, source, material)
abs.(RTeff)
```

## References

[1] Kong, Jin Au, Leung Tsang, Kung-Hau Ding, and Chi On Ao. Scattering of electromagnetic waves: numerical simulations. John Wiley & Sons, 2004.

[2] Gerhard Kristensson. "The Percus-Yevick approximation". github.com/JuliaWaveScattering/EffectiveWaves.jl (2022).

[3] Gower, Artur L., and Gerhard Kristensson. "Effective waves for random three-dimensional particulate materials." New Journal of Physics 23.6 (2021): 063083.