Fixing reflection transmission from plate two media and adding test t…

…o compare with one medium case

src/effective_waves/plane_waves/boundary-condition.jl

@@ -252,7 +252,7 @@ function solve_boundary_condition(ω::T, k_eff::Complex{T}, eigvectors1::Array{C
         k_effz_m = -k_effz_p
 
         # Needed coefficients
-        Δ = 2 * k * k0 * ρ * ρ0 * cos(k0*Z) - 1im * ((k0 * ρ)^2 + (k * ρ0)^2 * sin(k0 * Z))
+        Δ = 2 * k * k0 * ρ * ρ0 * cos(k0*Z) - 1im * ((k0 * ρ)^2 + (k * ρ0)^2) * sin(k0 * Z)
         D_p = k0 * ρ * (k0 * ρ + k * ρ0) / Δ
         D_m = k0 * ρ * (k0 * ρ - k * ρ0) / Δ
         D_1 = ((k0 * ρ)^2 - (k * ρ0)^2) / 2Δ
@@ -264,7 +264,7 @@ function solve_boundary_condition(ω::T, k_eff::Complex{T}, eigvectors1::Array{C
         Ys = spherical_harmonics(basis_order, rθφ[2], rθφ[3])
         lm_to_n = lm_to_spherical_harmonic_index
 
-        Pr(x::Complex{T}, y::Complex{T}, r::Float64) = exp(2im * (x + y) / Z) * sin((x + y) * (Z / 2 - r)) / (x + y)
+        Pr(x::Complex{T}, y::Complex{T}, r::T) = exp(1im * (x + y) * Z / 2) * sin((x + y) * (Z / 2 - r)) / (x + y)
 
         Bp(F::Array{Complex{T}}, p_or_m::Int) = (4π / (k0^2)) * sum(
             1im^(-T(l)) * Ys[lm_to_n(l, m)] * Pr(p_or_m * k_eff, -k0, rs[j]) * F[lm_to_n(l, m), j, p] * nf[j]
@@ -274,22 +274,23 @@ function solve_boundary_condition(ω::T, k_eff::Complex{T}, eigvectors1::Array{C
             1im^(T(l)) * Ys[lm_to_n(l, m)] * Pr(p_or_m * k_eff, k0, rs[j]) * F[lm_to_n(l, m), j, p] * nf[j]
         for p = 1:size(F, 3), l = 0:basis_order for m = -l:l, j in eachindex(species))
 
-        Iu(F::Array{Complex{T}}, p_or_m::Int) = 4π * sum(
-            1im^(-T(dl)) * Ys[lm_to_n(dl, dm)] * exp(1im * (k0 + p_or_m * k_eff) * rs[j]) * (k_eff + p_or_m * k0) * F[lm_to_n(dl, dm), j, p] * nf[j]
+        Iu(F::Array{Complex{T}}, p_or_m::Int) = (2π / (k0^2)) * sum(
+            1im^(-T(dl)) * Ys[lm_to_n(dl, dm)] * exp(1im * (p_or_m * k_eff - k0) * rs[j]) * F[lm_to_n(dl, dm), j, p] * nf[j] / (k_eff - p_or_m * k0)
         for p = 1:size(F, 3), dl = 0:basis_order for dm = -dl:dl, j in eachindex(species))
 
-        Il(F::Array{Complex{T}}, p_or_m::Int) = 4π * sum(
-            1im^(T(dl)) * Ys[lm_to_n(dl, dm)] * exp(1im * (k0 + p_or_m * k_eff) * (Z - rs[j])) * (k_eff - p_or_m * k0) * F[lm_to_n(dl, dm), j, p] * nf[j]
+        Il(F::Array{Complex{T}}, p_or_m::Int) = (2π / (k0^2)) * sum(
+            1im^(T(dl)) * Ys[lm_to_n(dl, dm)] * exp(1im * (p_or_m * k_eff + k0) * (Z - rs[j])) * F[lm_to_n(dl, dm), j, p] * nf[j] / (k_eff + p_or_m * k0)
         for p = 1:size(F, 3), dl = 0:basis_order for dm = -dl:dl, j in eachindex(species))
 
-        MI11 = D_2 * Bp(F_p, 1) * exp(1im * k0 * Z) + D_1 * Bm(F_p, 1) * exp(-1im * k0 * Z) - Iu(F_p, 1)
-        MI12 = D_2 * Bp(F_m, -1) * exp(1im * k0 * Z) + D_1 * Bm(F_m, -1) * exp(-1im * k0 * Z) - Iu(F_m, -1)
-        MI21 = (D_1 * Bp(F_p, 1) + D_2 * Bm(F_p, 1)) * exp(1im * k0 * Z) + Il(F_p, 1)
-        MI22 = (D_1 * Bp(F_m, -1) + D_2 * Bm(F_m, -1)) * exp(1im * k0 * Z) + Il(F_m, -1)
+        MI11 = D_2 * Bp(F_p, 1) * exp(1im * k0 * Z) + D_1 * Bm(F_p, 1) * exp(-1im * k0 * Z) + 1im * Iu(F_p, 1)
+        MI12 = D_2 * Bp(F_m, -1) * exp(1im * k0 * Z) + D_1 * Bm(F_m, -1) * exp(-1im * k0 * Z) + 1im * Iu(F_m, -1)
+        MI21 = (D_1 * Bp(F_p, 1) + D_2 * Bm(F_p, 1)) * exp(1im * k0 * Z) - 1im * Il(F_p, 1)
+        MI22 = (D_1 * Bp(F_m, -1) + D_2 * Bm(F_m, -1)) * exp(1im * k0 * Z) - 1im * Il(F_m, -1)
 
-        MIs = [MI11 MI12; MI21 MI22]
+        MIs = [MI11 MI12;
+               MI21 MI22]
 
-        forcing = G * γ0 * [D_p * exp(-1im * k0 * Z), D_m * exp(1im * k0 * Z)]
+        forcing = -G * γ0 * [D_p * exp(-1im * k0 * Z), D_m * exp(1im * k0 * Z)]
 
         α = MIs \ forcing
 

src/effective_waves/plane_waves/reflection-transmission.jl

@@ -167,20 +167,20 @@ function reflection_transmission_coefficients(wavemodes::Vector{E}, psource::Pla
         k0z = sqrt(k0^2 - k^2 + kz^2)
 
         # Needed coefficients
-        Δ = 2 * k * k0 * ρ * ρ0 * cos(k0 * Z) - 1im * ((k0 * ρ)^2 + (k * ρ0)^2 * sin(k0 * Z))
+        Δ = 2 * k * k0 * ρ * ρ0 * cos(k0 * Z) - 1im * ((k0 * ρ)^2 + (k * ρ0)^2) * sin(k0 * Z)
         D_p = k0 * ρ * (k0 * ρ + k * ρ0) / Δ
         D_m = k0 * ρ * (k0 * ρ - k * ρ0) / Δ
         D_0 = 2 * k * k0 * ρ * ρ0 / Δ
         D_1 = ((k0 * ρ)^2 - (k * ρ0)^2) / 2Δ
 
-        Pr(x::Complex{T}, y::Complex{T}, r::T) = exp(2im * (x + y) / Z) * sin((x + y) * (Z / 2 - r)) / (x + y)
+        Pr(x::Complex{T}, y::Complex{T}, r::T) = exp(1im * (x + y) * Z / 2) * sin((x + y) * (Z / 2 - r)) / (x + y)
 
-        Bp(Fp::Array{Complex{T}}, Fm::Array{Complex{T}}) = (4π / (k0 * k0z)) * sum(
+        Bp(Fp::Array{Complex{T}}, Fm::Array{Complex{T}}) = (4π / (k0^2)) * sum(
             1im^(-T(l)) * Ys[lm_to_n(l, m)] * (Pr(k_eff, -k0, rs[j]) * Fp[lm_to_n(l, m), j, p] +
                                             Pr(-k_eff, -k0, rs[j]) * Fm[lm_to_n(l, m), j, p]) * nf[j]
         for p = 1:size(Fp, 3), l = 0:basis_order for m = -l:l, j in eachindex(species))
 
-        Bm(Fp::Array{Complex{T}}, Fm::Array{Complex{T}}) = (4π / (k0 * k0z)) * sum(
+        Bm(Fp::Array{Complex{T}}, Fm::Array{Complex{T}}) = (4π / (k0^2)) * sum(
             1im^(T(l)) * Ys[lm_to_n(l, m)] * (Pr(k_eff, k0, rs[j]) * Fp[lm_to_n(l, m), j, p] +
                                          Pr(-k_eff, k0, rs[j]) * Fm[lm_to_n(l, m), j, p]) * nf[j]
         for p = 1:size(Fp, 3), l = 0:basis_order for m = -l:l, j in eachindex(species))
@@ -193,10 +193,6 @@ function reflection_transmission_coefficients(wavemodes::Vector{E}, psource::Pla
                 D_p * Bp(wavemodes[1].eigenvectors, wavemodes[2].eigenvectors) +
                 D_0 * G) * exp(-1im * k * Z)
 
-        # direction_ref = psource.direction - 2 * dot(n,psource.direction) * n
-        # reflected_wave = PlaneSource(psource.medium; direction = direction_ref, amplitude = Ramp)
-        # transmitted_wave = PlaneSource(psource.medium; direction = psource.direction, amplitude = Tamp)
-
         return [Ramp, Tamp]
     else
         # Unpacking parameters

test/acoustics-3D/Two-media/planar-symmetry.jl

@@ -2,7 +2,7 @@ using EffectiveWaves, Test, LinearAlgebra
 
 @testset "Compare one and two media halfspace" begin
 
-    # Low frequency test
+    # Two media test
     spatial_dim = 3
     # Choose the frequency
     ωs = 0.001:0.5:3.001
@@ -15,7 +15,7 @@ using EffectiveWaves, Test, LinearAlgebra
     radius1 = 1.0
     s1 = Specie(
         Acoustic(spatial_dim; ρ=10.0, c=10.0), radius1;
-        volume_fraction=0.05
+        volume_fraction=0.2
     )
     species = [s1]
 
@@ -43,26 +43,74 @@ using EffectiveWaves, Test, LinearAlgebra
     Reffs1 = [reflection_coefficient(w, source1, material) for w in wavemodes1]
     Reffs2 = [reflection_coefficient(w, source2, material) for w in wavemodes2]
 
-    @test abs(Reffs2[1] - Reffs1[1]) / abs(Reffs1[1]) < 0.05
-    @test norm(Reffs2 - Reffs1) / norm(Reffs1) < 0.1
+    @test abs(Reffs2[1] - Reffs1[1]) / abs(Reffs1[1]) < 0.01
+    @test norm(Reffs2 - Reffs1) / norm(Reffs1) < 0.01
+end
+
+@testset "Compare one and two media plate" begin
+
+    # Two media test
+    spatial_dim = 3
+    # Choose the frequency
+    ωs = 0.001:0.5:3.001
+    basis_order = 2
+
+    medium = Acoustic(spatial_dim; ρ=1.001, c=1.001)
+    medium0 = Acoustic(spatial_dim; ρ=1.0, c=1.0)
+
+    # Choose the species
+    radius1 = 1.0
+    s1 = Specie(
+        Acoustic(spatial_dim; ρ=0.01, c=0.01), radius1;
+        volume_fraction=0.2
+    )
+    species = [s1]
+
+    # Define the halfspace
+    r = maximum(outer_radius.(species))
+    normal = [0.0, 0.0, -1.0] # an outward normal to both surfaces of the plate
+    width = 200.0;
+    plate = Plate(normal, width, [0.0, 0.0, width / 2])
+
+    # define a plane wave sources with slightly different material
+    source1 = PlaneSource(medium0, [0.0, 0.0, 1.0])
+    source2 = PlaneSource(medium, [0.0, 0.0, 1.0])
+
+    # Calculate the effective wavenumber and wavemode numerically from the general methods
+    kp_arr = [wavenumbers(ω, medium0, species; tol=1e-6, num_wavenumbers=2, basis_order=basis_order) for ω in ωs]
+    k_effs = [kps[1] for kps in kp_arr]
+
+    # Define material 
+    material = Material(medium0, plate, species)
+
+    # Compute both wavemodes, for one and two slightly different media
+    wavemodes1 = [WaveMode(ωs[i], k_effs[i], source1, material; tol=1e-6, basis_order=basis_order) for i in eachindex(ωs)]
+    wavemodes2 = [WaveMode(ωs[i], k_effs[i], source2, material; tol=1e-6, basis_order=basis_order) for i in eachindex(ωs)]
+
+    # Compute reflection coefficients for both cases
+    RTs1 = [reflection_transmission_coefficients(w, source1, material) for w in wavemodes1]
+    RTs2 = [reflection_transmission_coefficients(w, source2, material) for w in wavemodes2]
+
+    @test norm(RTs2[1] - RTs1[1]) / norm(RTs1[1]) < 0.001
+    @test sum(norm.(RTs2 - RTs1)) / sum(norm.(RTs1)) < 0.01
 end
 
 @testset "Compare low frequency halfspace two media" begin
 
     # Low frequency test
     spatial_dim = 3
     # Choose the frequency
-    ωs = 0.001:0.02:0.101
+    ωs = 0.0001:0.0002:0.0011
     basis_order = 2
 
-    medium = Acoustic(spatial_dim; ρ = 2.0, c = 2.5)
+    medium = Acoustic(spatial_dim; ρ = 3.0, c = 3.0)
     medium0 = Acoustic(spatial_dim; ρ = 1.0, c = 1.0)
 
     # Choose the species
     radius1 = 1.0
     s1 = Specie(
         Acoustic(spatial_dim; ρ=10.0, c=10.0), radius1;
-        volume_fraction=0.05
+        volume_fraction=0.2
     );
     species = [s1]
 
@@ -90,7 +138,8 @@ end
 
     Reffs = [reflection_coefficient(w, source, material) for w in wavemodes]
 
-    @test abs(Reffs[1] - Rlows[1]) / abs(Rlows[1]) < 0.0001
-    @test norm(Reffs - Rlows) / norm(Rlows) < 0.005
+    @test abs(Reffs[1] - Rlows[1]) / abs(Rlows[1]) < 0.00002
+    @test norm(Reffs - Rlows) / norm(Rlows) < 0.0002
 end
 
+