wrote scattering from effective plate. Still needs testing.

Project.toml

@@ -19,7 +19,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-ApproxFun = "0.9, 0.10, ^0.11"
+ApproxFun = "0.9, 0.10, ^0.11,^0.12"
 GSL = "0.6, 0.7, 0.8, 0.9, 1.0"
 MDBM = "0.1"
 MultipleScattering = "^0.1.5"
@@ -28,7 +28,7 @@ Optim = "0.22, 1"
 Polynomials = "0.8"
 RecipesBase = "1.0"
 Reexport = "0.2"
-SpecialFunctions = "0.7, 0.8, 0.9, 10.0"
+SpecialFunctions = "0.3, 0.7, 0.8, 0.9, 0.10, 1.0"
 StaticArrays = "0.8,0.9,0.10,0.11,0.12"
 julia = "1"
 

docs/src/manual/reflection.md

@@ -119,30 +119,32 @@ s1 = Specie(
 );
 radius2 = 0.002
 s2 = Specie(
-    Acoustic(spatial_dim; ρ=0.2, c=4.1), radius2;
+    Acoustic(spatial_dim; ρ=0.00001, c=0.1), radius2;
+    #Acoustic(spatial_dim; ρ=0.2, c=4.1), radius2;
     volume_fraction=0.15
 
 );
 species = [s1,s2]
+species = [s2]
 
 # Choose the frequency
 ω = 1e-2
 k = ω / real(medium.c)
 
-# Check that we are indeed in a low frequency limit
-k * radius2 < 1e-4
-
+# For the limit of low frequencies we can define
 eff_medium = effective_medium(medium, species)
 
-normal = [0.0,0.0,-1.0] # an outward normal to the surface
-
-# Define the material region
-material = Material(Halfspace(normal),species)
+# Define a plate
+normal = [0.0,0.0,-1.0] # an outward normal to both surfaces of the plate
+width = 1.0 # plate width
+plate1 = Plate(normal,width)
 
 # define a plane wave source travelling at a 45 degree angle in relation to the material
 source = PlaneSource(medium, [cos(pi/4.0),0.0,sin(pi/4.0)])
 
-R = reflection_coefficient(ω,source, eff_medium, material.shape)
+amps = planewave_amplitudes(ω, source, eff_medium, plate1);
+Ramp = amps[1]
+Tamp = amps[2]
 
 # output
 
@@ -155,41 +157,18 @@ Using only one plane wave mode we can calculate both reflection and transmission
 
 ```julia 2
 k_effs = wavenumbers(ω, medium, species; tol = 1e-6, num_wavenumbers = 1, basis_order = 1)
-
-# Define a plate
-normal = [0.0,0.0,-1.0] # an outward normal to both surfaces of the plate
-width = 1.0 # plate width
-
-# Define the material region
-material = Material(Plate(normal,width),species)
-
-basis_order = 1;
-kws = Dict(:basis_order => basis_order)
-
-MM = eigensystem(ω, source, material; kws...)
-
-# calculate eigenvectors
-using LinearAlgebra
-
 k_eff = k_effs[1]
-MM_svd = svd(MM(k_eff))
-MM_svd.S[end]
-v1 = MM_svd.V[:,end]
-
-k_eff = -k_effs[1]
-MM_svd = svd(MM(k_eff))
-MM_svd.S[end]
-v2 = MM_svd.V[:,end]
-
 
+abs(k_eff - ω / eff_medium.c) < 1e-12
 
-
-vecs = eigenvectors(ω, k_effs[1], source, material; kws...)
-
+plate = Plate(normal,width)
+material = Material(plate,species)
 
 # Calculate the wavemode for the first wavenumber
-wave1 = WaveMode(ω, k_effs[1], source, material; tol = 1e-6, basis_order = 1)
+# 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);
 
+material_scattering_coefficients(wavemodes, source, material)
 
 ```
 Currently implementing... formulas from [Gower & Kristensson 2020](https://arxiv.org/pdf/2010.00934.pdf).

src/EffectiveWaves.jl

@@ -20,7 +20,7 @@ export dispersion_equation, dispersion_complex, eigensystem # supplies a matrix
 
 export scattering_amplitudes_average
 
-export  reflection_coefficient, reflection_coefficients
+export  reflection_coefficient, reflection_coefficients, planewave_amplitudes
 export  effective_medium
 
 # List of shorthand for some materials

src/acoustics/effective_wave/planewaves/boundary_condition.jl

@@ -43,38 +43,38 @@ function solve_boundary_condition(ω::T, k_eff::Complex{T}, eigvectors1::Array{C
 
     scale_number_density = one(T) - one(T) / material.numberofparticles
 
-    # First we calculate the innerward pointing normal
+    # First we calculate the outward pointing normal
     n = material.shape.normal;
-    n = n .* sign(real(dot(n,psource.direction)));
+    n = - n .* sign(real(dot(n,psource.direction)));
 
     # calculate the distances from the origin to the two faces of the plate
-    Z0 = dot(material.shape.normal, material.shape.origin)
-    Z1 = Z0 - material.shape.width / T(2)
-    Z2 = Z0 + material.shape.width / T(2)
+    Z0 = dot(-n, material.shape.origin)
+    Z1 = Z0 - material.shape.width / 2
+    Z2 = Z0 + material.shape.width / 2
 
     # next the components of the wavenumbers in the direction of the inward normal
-    direction1 = transmission_direction(k_eff, (ω / psource.medium.c) * psource.direction, -n)
-    direction2 = transmission_direction(- k_eff, (ω / psource.medium.c) * psource.direction, -n)
+    direction1 = transmission_direction(k_eff, (ω / psource.medium.c) * psource.direction, n)
+    direction2 = transmission_direction(- k_eff, (ω / psource.medium.c) * psource.direction, n)
 
     k = (ω / psource.medium.c)
-    kz = k * dot(conj(n), psource.direction)
-    k1_effz = k_eff * dot(conj(n), direction1)
-    k2_effz = k_eff * dot(conj(n), direction2)
+    kz = k * dot(-conj(n), psource.direction)
+    k1_effz = k_eff * dot(-conj(n), direction1)
+    k2_effz = k_eff * dot(-conj(n), direction2)
 
     rθφ = cartesian_to_radial_coordinates(psource.direction);
 
     Ys = spherical_harmonics(basis_order, rθφ[2], rθφ[3]);
     lm_to_n = lm_to_spherical_harmonic_index
 
     I1(F::Array{Complex{T}}, k_effz::Complex{T}) = sum(
-        - F[lm_to_n(dl,dm),j,1] * Ys[lm_to_n(dl,dm)] * im^T(dl+1) * T(2π) * T(1)^dl *
+        - F[lm_to_n(dl,dm),j,1] * Ys[lm_to_n(dl,dm)] * im^T(dl+1) * T(2π) * T(-1)^dl *
         exp(im * (k_effz - kz)*(Z1 + outer_radius(material.species[j]))) *
         scale_number_density * number_density(material.species[j]) / (k*kz * (kz - k_effz))
     for dl = 0:basis_order for dm = -dl:dl, j in eachindex(material.species))
 
     I2(F::Array{Complex{T}}, k_effz::Complex{T}) = sum(
-        - F[lm_to_n(dl,dm),j,1] * Ys[lm_to_n(dl,dm)] * im^T(dl+1) * T(2π) * T(1)^dm *
-        exp(im * (k_effz - kz)*(Z2 - outer_radius(material.species[j]))) *
+        - F[lm_to_n(dl,dm),j,1] * Ys[lm_to_n(dl,dm)] * im^T(dl+1) * T(2π) * T(-1)^dm *
+        exp(im * (k_effz + kz)*(Z2 - outer_radius(material.species[j]))) *
         scale_number_density * number_density(material.species[j]) / (k*kz * (kz + k_effz))
     for dl = 0:basis_order for dm = -dl:dl, j in eachindex(material.species))
 

src/acoustics/effective_wave/planewaves/reflection_coefficient.jl → src/acoustics/effective_wave/planewaves/material_scattering_coefficients.jl

@@ -30,6 +30,61 @@ function reflection_coefficient(ω::T, wave_eff::EffectivePlaneWaveMode{T}, psou
     return R
 end
 
+function material_scattering_coefficients(wavemodes::Vector{E}, psource::PlaneSource{T,3,1,Acoustic{T,3}}, material::Material{3,Plate{T,3}}) where {T<:AbstractFloat,Dim, E<:EffectivePlaneWaveMode{T,Dim}}
+
+    # Unpacking parameters
+    k = real(wavemodes[1].ω / psource.medium.c)
+
+    species = material.species
+    S = length(species)
+    rs = outer_radius.(species)
+
+    ho = maximum(w.basis_order for w in wavemodes)
+    ls, ms = spherical_harmonics_indices(ho)
+
+    # make the normal outward pointing
+    plate = material.shape
+    n = plate.normal / norm(plate.normal)
+    if real(dot(n,psource.direction)) > 0
+        n = -n
+    end
+
+    rθφ = cartesian_to_radial_coordinates(psource.direction)
+    Ys = spherical_harmonics(ho, rθφ[2], rθφ[3]);
+
+    direction_ref = psource.direction - 2 * dot(n,psource.direction) * n
+    rθφ = cartesian_to_radial_coordinates(direction_ref)
+    Yrefs = spherical_harmonics(ho, rθφ[2], rθφ[3]);
+
+    Z0 = dot(-n,plate.origin - psource.position)
+    Z1 = Z0 - plate.width / 2
+    Z2 = Z0 + plate.width / 2
+
+    kcos_in = k * dot(- conj(n), psource.direction)
+
+    RTs = map(wavemodes) do w
+        kcos_eff = w.wavenumber * dot(- conj(n), w.direction)
+
+        Rp = sum(
+            number_density(species[i[2]]) * w.eigenvectors[i] * 2pi * (1.0im)^(ls[i[1]]-1) *
+            Yrefs[i[1]] * (exp(im*(kcos_eff + kcos_in)*(Z2 - rs[i[2]])) - exp(im*(kcos_eff + kcos_in)*(Z1 + rs[i[2]]))) /
+            ((kcos_in + kcos_eff) * k * kcos_in)
+        for i in CartesianIndices(w.eigenvectors))
+
+        Tp = sum(
+            number_density(species[i[2]]) * w.eigenvectors[i] * 2pi * (-1.0)^ls[i[1]] * (1.0im)^(ls[i[1]]+1) *
+            Ys[i[1]] * (exp(im*(kcos_eff - kcos_in)*(Z2 - rs[i[2]])) - exp(im*(kcos_eff - kcos_in)*(Z1 + rs[i[2]]))) /
+            ((kcos_in - kcos_eff) * k * kcos_in)
+        for i in CartesianIndices(w.eigenvectors))
+
+        [Rp, Tp]
+    end
+
+    Ramp, Tamp = sum(RTs) + [0.0,1.0]
+
+    return [Ramp, Tamp]
+end
+
 "The average reflection coefficient"
 function wienerhopf_reflection_coefficient(ω::T, psource::PlaneSource{T,2,1,Acoustic{T,2}}, material::Material{2,Halfspace{T,2}};
         tol::T = T(1e-7),

src/acoustics/export.jl

@@ -11,7 +11,7 @@ include("low_frequency.jl")
 include("effective_wave/asymptotics.jl")
 
 include("effective_wave/planewaves/eigensystem.jl")
-include("effective_wave/planewaves/reflection_coefficient.jl")
+include("effective_wave/planewaves/material_scattering_coefficients.jl")
 include("effective_wave/planewaves/boundary_condition.jl")
 include("effective_wave/planewaves/wavemode-wienerhopf.jl")
 

src/acoustics/low_frequency.jl

@@ -54,11 +54,11 @@ function reflection_coefficient(ω::T, psource::PlaneSource{T,Dim,1,Acoustic{T,D
 end
 
 """
-    planewave_coefficients(psource::PlaneSource, reflect_medium::Acoustic, plate::Plate)
+    planewave_amplitudes(psource::PlaneSource, reflect_medium::Acoustic, plate::Plate)
 
-Calculates the coefficients, or amplityudes, of the plane waves scattering from and transmitted inside the plate, for the plane wave source `psource` and the `plate`. The function returns `[R, T, P1, P2]` where `R` is the reflection coefficient, `T` is the coefficient of the transmitted wave, and `P1` (`P2`) are the amplitudes of the wave travling forward (backward) inside the plate.
+Calculates the coefficients, or amplitudes, of the plane waves scattering from and transmitted inside the plate, for the plane wave source `psource` and the `plate`. The function returns `[R, T, P1, P2]` where `R` is the reflection coefficient, `T` is the coefficient of the transmitted wave, and `P1` (`P2`) are the amplitudes of the wave travling forward (backward) inside the plate.
 """
-function planewave_coefficients(ω::T, psource::PlaneSource{T,Dim,1,Acoustic{T,Dim}}, plate_medium::Acoustic{T,Dim}, plate::Plate{T,Dim}) where {T<:AbstractFloat,Dim}
+function planewave_amplitudes(ω::T, psource::PlaneSource{T,Dim,1,Acoustic{T,Dim}}, plate_medium::Acoustic{T,Dim}, plate::Plate{T,Dim}) where {T<:AbstractFloat,Dim}
 
     k0 = ω / psource.medium.c
     k1 = if abs(plate_medium.c) == zero(T)
@@ -85,17 +85,12 @@ function planewave_coefficients(ω::T, psource::PlaneSource{T,Dim,1,Acoustic{T,D
     Z1 = Z0 - plate.width / T(2)
     Z2 = Z0 + plate.width / T(2)
 
-    UR = (C0 - C1)*(C0 + C1) * exp(2im*k0*Z1*cos(θ0)) * (exp(2im*k1*Z1*cos(θ1)) - exp(2im*k1*Z2*cos(θ1))) /
-        ((C0 + C1)^2 * exp(2im*k1*Z1*cos(θ1)) - (C0 - C1)^2 * exp(2im*k1*Z2*cos(θ1)))
+    denom = ((C0 + C1)^2 * exp(2im*k1*Z1*cos(θ1)) - (C0 - C1)^2 * exp(2im*k1*Z2*cos(θ1)));
 
-    UP1 = 2C0*(C0 + C1) * exp(im*Z1*(k0*cos(θ0) + k1*cos(θ1))) /
-        ((C0 + C1)^2*exp(2im*k1*Z1*cos(θ1)) - (C0 - C1)^2*exp(2im*k1*Z2*cos(θ1)))
-
-    UP2 = 2*C0*(C0 - C1)*exp(im*(k0*Z1*cos(θ0) + k1*(Z1 + 2*Z2)*cos(θ1))) /
-        (-((C0 + C1)^2*exp((2*im)*k1*Z1*cos(θ1))) + (C0 - C1)^2*exp(2im*k1*Z2*cos(θ1)))
-
-    UT = 4*C0*C1*exp(im*(k0*(Z1 - Z2)*cos(θ0) + k1*(Z1 + Z2)*cos(θ1))) /
-    ((C0 + C1)^2*exp(2im*k1*Z1*cos(θ1)) - (C0 - C1)^2*exp(2im*k1*Z2*cos(θ1)))
+    UR = (C0 - C1)*(C0 + C1) * exp(2im*k0*Z1*cos(θ0)) * (exp(2im*k1*Z1*cos(θ1)) - exp(2im*k1*Z2*cos(θ1))) / denom
+    UP1 = 2C0*(C0 + C1) * exp(im*Z1*(k0*cos(θ0) + k1*cos(θ1))) / denom
+    UP2 = - 2*C0*(C0 - C1)*exp(im*(k0*Z1*cos(θ0) + k1*(Z1 + 2*Z2)*cos(θ1))) / denom
+    UT = 4*C0*C1*exp(im*(k0*(Z1 - Z2)*cos(θ0) + k1*(Z1 + Z2)*cos(θ1))) / denom 
 
     return [UR, UT, UP1, UP2]
 end

src/effective_wave/wavemode.jl

@@ -52,7 +52,7 @@ function WaveMode(ω::T, wavenumber::Complex{T}, psource::PlaneSource{T,Dim,1},
     direction2 = transmission_direction(- wavenumber, (ω / psource.medium.c) * psource.direction, material.shape.normal)
     eigvectors2 = eigenvectors(ω, - wavenumber, psource, material; direction_eff = direction2, kws...)
 
-    α = solve_boundary_condition(ω, k_eff, eigvectors1, eigvectors2, psource, material; kws...)
+    α = solve_boundary_condition(ω, wavenumber, eigvectors1, eigvectors2, psource, material; kws...)
 
     # apply normalisation
     eigvectors1[:,:,1] = eigvectors1[:,:,1] .* α[1]