Final adjustments to the two medium case and tests

src/effective_waves/plane_waves/boundary-condition.jl

@@ -76,6 +76,8 @@ function solve_boundary_condition(ω::T, k_eff::Complex{T}, eigvectors::Array{Co
     else
 
         # Setting parameters
+        ρ = psource.medium.ρ
+        ρ0 = material.microstructure.medium.ρ
         c = psource.medium.c
         c0 = material.microstructure.medium.c
         k = ω / c
@@ -102,19 +104,19 @@ function solve_boundary_condition(ω::T, k_eff::Complex{T}, eigvectors::Array{Co
         k0z = sqrt(k0^2 - (k^2 - kz^2))
 
         # Needs adjustments for complex wavespeeds
-        direction0 = real(k) * (psource.direction - dot(-conj(n), psource.direction) * psource.direction)
-        direction0 += real(k0z) * dot(-conj(n), psource.direction) * psource.direction
-        direction0 /= norm(direction0)
+        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 = k_eff * dot(-conj(n), direction)
 
-        γ0 = kz/k0z
-        ζR = (kz - k0z) / (kz + k0z)
-        ζT = 2k0z / (kz + k0z)
+        γ0 = (ρ0 * kz) /(ρ * k0z)
+        ζR = (ρ0 * kz - ρ * k0z) / (ρ0 * kz + ρ * k0z)
+        ζT = 2ρ * k0z / (ρ0 * kz + ρ * k0z)
 
-        rθφ = cartesian_to_radial_coordinates(direction0)
+        rθφ = cartesian_to_radial_coordinates(direction_k0)
         # spherical_harmonics() needs adjustments for complex wavespeeds
         Ys = spherical_harmonics(basis_order, rθφ[2], rθφ[3])
         lm_to_n = lm_to_spherical_harmonic_index
@@ -124,12 +126,9 @@ function solve_boundary_condition(ω::T, k_eff::Complex{T}, eigvectors::Array{Co
         for p = 1:M, dl = 0:basis_order for dm = -dl:dl, j in eachindex(species))
 
         particle_contribution *= (2pi * 1im) * (k_effz + k0z) / (k0 * k0z * (k0^2 - k_eff^2))
-
-        rθφ_n = cartesian_to_radial_coordinates(-conj(n))
-        Ys_n = spherical_harmonics(basis_order, rθφ_n[2], rθφ_n[3])
 
         wall_contribution = sum(
-            -Fs[lm_to_n(dl,dm),j,p] * nf[j] * 1im^dl * Ys_n[lm_to_n(dl, dm)] * exp(1im * (k_effz + k0z) * (Z0 + rs[j]))
+            -Fs[lm_to_n(dl,dm),j,p] * nf[j] * 1im^dl * Ys[lm_to_n(dl, dm)] * exp(1im * (k_effz + k0z) * (Z0 + rs[j]))
         for p = 1:M, dl = 0:basis_order for dm = -dl:dl, j in eachindex(species))
 
         wall_contribution *= ζR * 2 * pi / (1im * k0 * k0z * (k_effz + k0z))

src/effective_waves/plane_waves/reflection-transmission.jl

@@ -35,42 +35,89 @@ end
 function reflection_coefficient(wavemode::EffectivePlaneWaveMode{T,Dim}, psource::PlaneSource{T,3,1,Acoustic{T,3}}, material::Material{Halfspace{T,3}}) where {T<:AbstractFloat,Dim}
 
 
-    if psource.medium != material.microstructure.medium @error mismatched_medium end
+    if psource.medium == material.microstructure.medium
 
-    # Unpacking parameters
-    k = wavemode.ω / psource.medium.c
+        # Unpacking parameters
+        k = wavemode.ω / psource.medium.c
 
-    species = material.microstructure.species
-    S = length(species)
-    rs = outer_radius.(species)
+        species = material.microstructure.species
+        S = length(species)
+        rs = outer_radius.(species)
 
-    basis_order = wavemode.basis_order
-    ls, ms = spherical_harmonics_indices(basis_order)
+        basis_order = wavemode.basis_order
+        ls, ms = spherical_harmonics_indices(basis_order)
 
-    # make the normal outward pointing
-    plate = material.shape
-    n = plate.normal / norm(plate.normal)
-    if real(dot(n,psource.direction)) > 0
-        n = -n
-    end
+        # 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)
-    Yrefs = spherical_harmonics(basis_order, rθφ[2], rθφ[3]);
+        rθφ = cartesian_to_radial_coordinates(psource.direction)
+        Yrefs = spherical_harmonics(basis_order, rθφ[2], rθφ[3]);
 
-    Z0 = dot(-n,plate.origin)
-    kcos_in = k * dot(- conj(n), psource.direction)
+        Z0 = dot(-n,plate.origin)
+        kcos_in = k * dot(- conj(n), psource.direction)
 
-    kcos_eff = wavemode.wavenumber * dot(- conj(n), wavemode.direction)
+        kcos_eff = wavemode.wavenumber * dot(- conj(n), wavemode.direction)
 
-    Ramp = - sum(number_density(species[i[2]]) * wavemode.eigenvectors[i] * 2pi * (-1.0)^ms[i[1]] * (1.0im)^(ls[i[1]]-1) *
-        Yrefs[i[1]] * exp(im*(kcos_eff + kcos_in)*(Z0 + rs[i[2]])) / ((kcos_in + kcos_eff) * k * kcos_in)
-    for i in CartesianIndices(wavemode.eigenvectors))
-    Ramp = Ramp * exp(im * kcos_in * Z0)
+        Ramp = - sum(number_density(species[i[2]]) * wavemode.eigenvectors[i] * 2pi * (-1.0)^ms[i[1]] * (1.0im)^(ls[i[1]]-1) *
+            Yrefs[i[1]] * exp(im*(kcos_eff + kcos_in)*(Z0 + rs[i[2]])) / ((kcos_in + kcos_eff) * k * kcos_in)
+        for i in CartesianIndices(wavemode.eigenvectors))
+        Ramp = Ramp * exp(im * kcos_in * Z0)
 
-    # direction_ref = psource.direction - 2 * dot(n,psource.direction) * n
-    # reflected_wave = PlaneSource(psource.medium; direction = direction_ref, amplitude = Ramp)
+        return Ramp
+
+    else
+
+        # Unpacking parameters
+        ρ = psource.medium.ρ
+        ρ0 = material.microstructure.medium.ρ
+        c = psource.medium.c
+        c0 = material.microstructure.medium.c
+        ω = wavemode.ω
+        k = ω / c
+        k0 = ω / c0
+
+        species = material.microstructure.species
+        S = length(species)
+        rs = outer_radius.(species)
+
+        basis_order = wavemode.basis_order
+        ls, ms = spherical_harmonics_indices(basis_order)
+
+        # make the normal outward pointing
+        halfspace = material.shape
+        n = halfspace.normal / norm(halfspace.normal)
+        if real(dot(n,psource.direction)) > 0
+            n = -n
+        end
+
+        Z0 = dot(-n, halfspace.origin)
 
-    return Ramp
+        kz = k * dot(-conj(n), psource.direction)
+        k0z = sqrt(k0^2 - (k^2 - kz^2))
+        k_effz = wavemode.wavenumber * dot(-conj(n), wavemode.direction)
+
+        ζR = (ρ0 * kz - ρ * k0z) / (ρ0 * kz + ρ * k0z)
+        ζT = 2ρ * k0z / (ρ0 * kz + ρ * k0z)
+
+        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)
+
+        rθφ = cartesian_to_radial_coordinates(direction_k0)
+        Yrefs = spherical_harmonics(basis_order, rθφ[2], rθφ[3]);
+
+        B = - sum(number_density(species[i[2]]) * wavemode.eigenvectors[i] * 2pi * (1.0im)^(ls[i[1]]-1) *
+             Yrefs[i[1]] * exp(im * (k_effz + k0z) * (Z0 + rs[i[2]])) / (k0 * k0z * (k_effz + k0z))
+        for i in CartesianIndices(wavemode.eigenvectors))
+
+        Ramp = (ζR * field(psource, zeros(T, 3), ω) + ζT * B) * exp(im * kz * Z0)
+
+        return Ramp
+    end
 end
 
 function reflection_transmission_coefficients(wavemodes::Vector{E}, psource::PlaneSource{T,3,1,Acoustic{T,3}}, material::Material{Plate{T,3}}) where {T<:AbstractFloat,Dim, E<:EffectivePlaneWaveMode{T,Dim}}

src/effective_waves/wavemode.jl

@@ -57,7 +57,7 @@ end
 function WaveMode(ω::T, wavenumber::Complex{T}, psource::PlaneSource{T,Dim,1}, material::Material{Halfspace{T,Dim}};
     tol::T = 1e-6, kws...) where {T,Dim}
 
-    direction = transmission_direction(wavenumber, (ω / psource.medium.c) * psource.direction, material.shape.normal)
+    direction = transmission_direction(wavenumber, (ω / material.microstructure.medium.c) * psource.direction, material.shape.normal)
     eigvectors = eigenvectors(ω, wavenumber, psource, material; direction_eff = direction, kws...)
 
     α = solve_boundary_condition(ω, wavenumber, eigvectors, psource, material; kws...)

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

@@ -0,0 +1,96 @@
+using EffectiveWaves, Test, LinearAlgebra
+
+@testset "Compare one and two media halfspace" begin
+
+    # Low frequency 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; ρ=10.0, c=10.0), radius1;
+        volume_fraction=0.05
+    )
+    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)
+
+    # 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, halfspace, 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
+    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
+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
+    basis_order = 2
+
+    medium = Acoustic(spatial_dim; ρ = 2.0, c = 2.5)
+    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
+    );
+    species = [s1]
+
+    # For the limit of low frequencies we can define
+    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)
+
+    # 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)])
+    source = PlaneSource(medium, [0.0,0.0,1.0])
+
+    # reflection low freq approx
+    Rlows = [reflection_coefficient(ω, source, eff_medium, halfspace) for ω in ωs];
+
+    # 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,halfspace,species)
+    wavemodes = [WaveMode(ωs[i], k_effs[i], source, material; tol = 1e-6, basis_order = basis_order) for i in eachindex(ωs)];
+
+    Reffs = [reflection_coefficient(w, source, material) for w in wavemodes]
+
+    @test abs(Reffs[1] - Rlows[1]) / abs(Rlows[1]) < 0.01
+    @test norm(Reffs - Rlows) / norm(Rlows) < 0.1
+end
+

test/runtests.jl

@@ -29,6 +29,9 @@ include("complex.jl")
 
     include("pair-correlation.jl")
 
+# Eigensystems and fields for the case of two media
+    include("acoustics-3D/Two-media/planar-symmetry.jl")
+
 # test equivalence between methods for finding wavenumbers
     # test does not run on Julia version < 0.7 due to differences in Optim versions
     # include("path_mesh_wavenumbers.jl")