Fixing two media plate and adding more tests

src/effective_waves/plane_waves/boundary-condition.jl

@@ -237,26 +237,14 @@ function solve_boundary_condition(ω::T, k_eff::Complex{T}, eigvectors1::Array{C
         Z1 = Z0 - Z / 2
         Z2 = Z0 + Z / 2
 
-        # next the components of the wavenumbers in the direction of the inward normal
-        kz = k * dot(-conj(n), psource.direction)
-        k0z = sqrt(k0^2 - (k^2 - kz^2))
-
-        # Needs adjustments for complex wavespeeds
-        direction_k0 = real(k) * (psource.direction - dot(-conj(n), psource.direction) * psource.direction)
-        direction_k0 = direction_k0 + real(k0z) * dot(-conj(n), psource.direction) * psource.direction
-        direction_k0 = direction_k0 / norm(direction_k0)
-
-        direction = transmission_direction(k_eff, k * psource.direction, n)
-
-        k_effz_p = k_eff * dot(-conj(n), direction)
-        k_effz_m = -k_effz_p
+        direction = [0.0, 0.0, 1.0]
 
         # 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Δ
-        D_2 = (k0 * ρ - k * ρ0)^2 / 2Δ
+        D_1 = ((k0 * ρ)^2 - (k * ρ0)^2) / (2 * Δ)
+        D_2 = (k0 * ρ - k * ρ0)^2 / (2 * Δ)
         γ0 = k * ρ0 / (k0 * ρ)
 
         rθφ = cartesian_to_radial_coordinates(direction)

src/effective_waves/plane_waves/reflection-transmission.jl

@@ -171,7 +171,7 @@ function reflection_transmission_coefficients(wavemodes::Vector{E}, psource::Pla
         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Δ
+        D_1 = ((k0 * ρ)^2 - (k * ρ0)^2) / (2 * Δ)
 
         Pr(x::Complex{T}, y::Complex{T}, r::T) = exp(1im * (x + y) * Z / 2) * sin((x + y) * (Z / 2 - r)) / (x + y)
 
@@ -193,7 +193,7 @@ 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)
 
-        return [Ramp, Tamp]
+        return [Ramp, Tamp * exp(im * kz * Z2)]
     else
         # Unpacking parameters
         ω = wavemodes[1].ω

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

@@ -20,7 +20,6 @@ using EffectiveWaves, Test, LinearAlgebra
     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
     halfspace = Halfspace(normal)
 
@@ -55,24 +54,23 @@ end
     ωs = 0.001:0.5:3.001
     basis_order = 2
 
-    medium = Acoustic(spatial_dim; ρ=1.001, c=1.001)
+    medium = Acoustic(spatial_dim; ρ=1.00001, c=1.00001)
     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
+        Acoustic(spatial_dim; ρ=10.0, c=10.0), radius1;
+        volume_fraction=0.1
     )
     species = [s1]
 
-    # Define the halfspace
-    r = maximum(outer_radius.(species))
+    # Define the plate
     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])
+    Z = 200.0;
+    plate = Plate(normal, Z, [0.0, 0.0, Z / 2])
 
-    # define a plane wave sources with slightly different material
+    # define 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])
 
@@ -91,8 +89,10 @@ end
     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
+    @test abs(RTs2[1][1] - RTs1[1][1]) / abs(RTs2[1][1]) < 0.02
+    @test abs(RTs1[1][2] - RTs2[1][2]) / abs(RTs2[1][2]) < 0.003
+    errorsT = [norm([abs(RTs1[i][2] - RTs2[i][2])]) for i in eachindex(ωs)]
+    @test norm(errorsT) / norm(norm.(RTs2)) < 0.005
 end
 
 @testset "Compare low frequency halfspace two media" begin
@@ -118,7 +118,6 @@ end
     eff_medium = effective_medium(medium0, species)
 
     # Define the halfspace
-    r = maximum(outer_radius.(species))
     normal = [0.0,0.0,-1.0] # an outward normal to both surfaces of the plate
     halfspace = Halfspace(normal)
 
@@ -142,4 +141,61 @@ end
     @test norm(Reffs - Rlows) / norm(Rlows) < 0.0002
 end
 
+@testset "Compare low frequency plate two media" begin
+
+    # Low frequency test
+    spatial_dim = 3
+    # Choose the frequency
+    ωs = 0.0001:0.0002:0.0011
+    basis_order = 2
+
+    # Define all properties
+    c = 3.0
+    ρ = 3.0
+    c0 = 1.0
+    ρ0 = 1.0
+    cs = 10.0
+    ρs = 10.0
+    ϕ = 0.1
+
+    medium = Acoustic(spatial_dim; ρ = ρ, c = c)
+    medium0 = Acoustic(spatial_dim; ρ = ρ0, c = c0)
+
+    # Choose the species
+    r = 1.0
+    s1 = Specie(
+        Acoustic(spatial_dim; ρ = ρs, c = cs), r;
+        volume_fraction = ϕ
+    );
+    species = [s1]
+
+    # Define the plate
+    normal = [0.0, 0.0, -1.0] # an outward normal to both surfaces of the plate
+    Z = 200.0
+    plate = Plate(normal, Z, [0.0, 0.0, Z / 2])
+
+    # define a plane wave source travelling to the material
+    source = PlaneSource(medium, [0.0,0.0,1.0])
+
+    # Reflection and transmittion low frequency approximation
+    eff_medium = effective_medium(medium0, species)
+    amps = [planewave_amplitudes(ω, source, eff_medium, plate) for ω in ωs]
+    Ramp = [amps[i][1] for i in eachindex(amps)]
+    Tamp = [amps[i][2] for i in eachindex(amps)]
+    RTlows = [[Ramp[i], Tamp[i]] for i in eachindex(amps)]
+
+    # 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]
+
+    material = Material(medium0,plate,species)
+    wavemodes = [WaveMode(ωs[i], k_effs[i], source, material; tol = 1e-6, basis_order = basis_order) for i in eachindex(ωs)];
+
+    RTeffs = [reflection_transmission_coefficients(w, source, material) for w in wavemodes]
 
+    @test abs(RTeffs[1][1] - RTlows[1][1]) / abs(RTlows[1][1]) < 0.09
+    @test abs(RTeffs[1][2] - RTlows[1][2]) / abs(RTlows[1][2]) < 0.009
+    errors = [norm([abs(RTeffs[i][1] - RTlows[i][1]), abs(RTeffs[i][2] - RTlows[i][2])]) for i in eachindex(ωs)]
+    @test norm(errors) / norm(norm.(RTlows)) < 0.06
+
+end