Skip to content

Commit

Permalink
Final adjustments in halfspace case with two different acoustic media
Browse files Browse the repository at this point in the history
  • Loading branch information
@pivaps
pivaps committed Mar 6, 2024
1 parent c3194d7 commit 7d713ab
Show file tree
Hide file tree
Showing 2 changed files with 44 additions and 61 deletions.
74 changes: 44 additions & 30 deletions src/effective_waves/plane_waves/boundary-condition.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,39 +43,38 @@ function solve_boundary_condition(ω::T, k_eff::Complex{T}, eigvectors::Array{Co
@warn "The effective wavenumber: $k_eff has more than one eigenvector. For plane-waves this case has not been fully implemented"
end

rθφ = cartesian_to_radial_coordinates(psource.direction);

Ys = spherical_harmonics(basis_order, rθφ[2], rθφ[3]);
lm_to_n = lm_to_spherical_harmonic_index

# if psource.medium == material.microstructure.medium
if psource.medium == material.microstructure.medium

# First we calculate the outward pointing normal
# n = material.shape.normal;
# n = - n .* sign(real(dot(n,psource.direction)));
n = material.shape.normal;
n = - n .* sign(real(dot(n,psource.direction)));

# calculate the distances from the origin to the face of the halfspace
# Z0 = dot(-n, material.shape.origin)
Z0 = dot(-n, material.shape.origin)

# next the components of the wavenumbers in the direction of the inward normal
# direction = transmission_direction(k_eff, (ω / psource.medium.c) * psource.direction, n)
direction = transmission_direction(k_eff, (ω / psource.medium.c) * psource.direction, n)

k = (ω / psource.medium.c)
kz = k * dot(-conj(n), psource.direction)
k_effz = k_eff * dot(-conj(n), direction)
rθφ = cartesian_to_radial_coordinates(psource.direction);

# k = (ω / psource.medium.c)
# kz = k * dot(-conj(n), psource.direction)
# k_effz = k_eff * dot(-conj(n), 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,p] * Ys[lm_to_n(dl,dm)] * im^T(dl+1) * T(2π) * T(-1)^dl *
# exp(im * (k_effz - kz)*(Z0 + outer_radius(material.microstructure.species[j]))) * number_density(material.microstructure.species[j]) / (k*kz * (kz - k_effz))
# for p = 1:size(F,3), dl = 0:basis_order for dm = -dl:dl, j in eachindex(material.microstructure.species))
I1(F::Array{Complex{T}}, k_effz::Complex{T}) = sum(
- F[lm_to_n(dl,dm),j,p] * Ys[lm_to_n(dl,dm)] * im^T(dl+1) * T(2π) * T(-1)^dl *
exp(im * (k_effz - kz)*(Z0 + outer_radius(material.microstructure.species[j]))) * number_density(material.microstructure.species[j]) / (k*kz * (kz - k_effz))
for p = 1:size(F,3), dl = 0:basis_order for dm = -dl:dl, j in eachindex(material.microstructure.species))

# forcing = - field(psource,zeros(T,3),ω)
forcing = - field(psource,zeros(T,3),ω)

# α = I1(eigvectors,k_effz) \ forcing
α = I1(eigvectors,k_effz) \ forcing

# return α
return α

# else
else

# Setting parameters
ρ = psource.medium.ρ
Expand All @@ -90,6 +89,10 @@ function solve_boundary_condition(ω::T, k_eff::Complex{T}, eigvectors::Array{Co
Fs = eigvectors
M = size(Fs,3)

if (c0 - real(c0) != 0)
@warn "The case of complex wavespeed has not been fully implemented for two medium case, using real part of wavespeeds instead."
end

# First we calculate the outward pointing normal
n = material.shape.normal;
n = - n .* sign(real(dot(n,psource.direction)));
Expand All @@ -98,27 +101,38 @@ function solve_boundary_condition(ω::T, k_eff::Complex{T}, eigvectors::Array{Co
Z0 = dot(-n, material.shape.origin)

# next the components of the wavenumbers in the direction of the inward normal
direction = transmission_direction(k_eff, k * psource.direction, ρ, ρ0, n)

k = (ω / c)
kz = k * dot(-conj(n), psource.direction)
k0z = kz * ρ0 / ρ
in_direction = k * (psource.direction - dot(-conj(n), psource.direction) * psource.direction)
in_direction += k0z * dot(-conj(n), psource.direction) * psource.direction

# 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 = 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)

# Needs adjustments for complex wavespeeds
rθφ = cartesian_to_radial_coordinates(direction0)
Ys = spherical_harmonics(basis_order, rθφ[2], rθφ[3])
lm_to_n = lm_to_spherical_harmonic_index

particle_contribution = sum(
-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]))
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))

particle_contribution *= (2pi * 1im) * (k_effz + k0z) / (k0 * k0z * (k_eff^2 - k0^2))
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[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_n[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))
Expand All @@ -130,7 +144,7 @@ function solve_boundary_condition(ω::T, k_eff::Complex{T}, eigvectors::Array{Co
α = M_dot \ forcing

return α
#end
end

end

Expand Down
31 changes: 0 additions & 31 deletions src/specialfunctions.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,37 +47,6 @@ function transmission_direction(k_eff::Complex{T}, psource::PlaneSource{T}, mat
transmission_direction(k_eff, psource.direction, material.shape.normal; tol = tol)
end

function transmission_direction(k_eff::Complex{T}, incident_wavevector::AbstractArray{CT} where {CT<:Union{T,Complex{T}}}, ρ::T,ρ0::T, surface_normal::AbstractArray{T};
tol::T=sqrt(eps(T))) where {T<:AbstractFloat,Dim}

ratioρ = ρ0 / ρ

# Let us first make the surface normal a unit vector which is outward pointing
surface_normal = surface_normal / norm(surface_normal)
if real(dot(surface_normal, incident_wavevector)) > 0
surface_normal = -surface_normal
end

wnp = incident_wavevector - dot(surface_normal, incident_wavevector) .* surface_normal
α = sqrt(k_eff^2 - sum(x^2 for x in wnp)) * ratioρ

# The conditions below guarantee (in order of priority) that either: 1) the wave attenuates when travelling into the material (imag(α) < 0), or 2) the wave does not attenuate, but travels into the material.
if tol * real(α) > abs(imag(α))
α = -α # leads to real(α) < 0
elseif imag(α) > 0
α = -α # leads to imag(α) < 0
elseif isnan(α)
α = -one(T) # covers the cases k_eff = Inf and k_eff = 0.0
return wnp + α .* surface_normal
end

return (wnp + α .* surface_normal) ./ k_eff
end

function transmission_direction(k_eff::Complex{T}, psource::PlaneSource{T}, ρ::T,ρ0::T, material::Material; tol::T=sqrt(eps(T))) where {T}
transmission_direction(k_eff, psource.direction, ρ::T, ρ0::T, material.shape.normal; tol=tol)
end

transmission_angle(pwave::Union{PlaneSource,EffectivePlaneWaveMode}, material::Material) = transmission_angle(pwave, material.shape)

transmission_angle(pwave::Union{PlaneSource,EffectivePlaneWaveMode}, shape::Union{Plate,Halfspace}) = transmission_angle(pwave.direction, shape.normal)
Expand Down

0 comments on commit 7d713ab

Please sign in to comment.