Skip to content

Commit

Permalink
rewrote many functions for PhysicalMedium{dim,1} instead of Acoustics
Browse files Browse the repository at this point in the history
  • Loading branch information
@arturgower
arturgower committed Nov 30, 2023
1 parent a24e9f3 commit 4eb9c4c
Show file tree
Hide file tree
Showing 6 changed files with 165 additions and 148 deletions.
4 changes: 2 additions & 2 deletions docs/src/example/particles_in_circle/README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,8 +64,8 @@ big_result = run(big_particle_simulation, x, ωs)

#plot(result, lab = "scattering from particles")
#plot!(big_result,
# lab = "scattering from big particle",
# title="Compare scattered wave from one big particle, \n and a circle filled with small particles")
# lab = "scattering from big particle",
# title="Compare scattered wave from one big particle, \n and a circle filled with small particles")
```
![The response comparison](plot_response_compare.png)

2 changes: 1 addition & 1 deletion src/MultipleScattering.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ export boundary_functions, boundary_points, boundary_data, bounding_box, corners
export points_in_shape, bottomleft, topright

## Physical mediums
export PhysicalMedium, outgoing_basis_function, regular_basis_function, outgoing_radial_basis, regular_radial_basis, outgoing_translation_matrix, regular_translation_matrix, estimate_regular_basisorder, estimate_outgoing_basisorder, basisorder_to_basislength, basis_coefficients, basislength_to_basisorder, internal_field, check_material
export PhysicalMedium, ScalarMedium, outgoing_basis_function, regular_basis_function, outgoing_radial_basis, regular_radial_basis, outgoing_translation_matrix, regular_translation_matrix, estimate_regular_basisorder, estimate_outgoing_basisorder, basisorder_to_basislength, basis_coefficients, basislength_to_basisorder, internal_field, check_material
export spatial_dimension, field_dimension

## Electromagnetic
Expand Down
4 changes: 4 additions & 0 deletions src/particle.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -134,6 +134,10 @@ issubset(s::Shape, particle::AbstractParticle) = issubset(s, shape(particle))
issubset(particle::AbstractParticle, s::Shape) = issubset(shape(particle), s)
issubset(particle1::AbstractParticle, particle2::AbstractParticle) = issubset(shape(particle1), shape(particle2))


# NOTE that medium in both functions below is only used to get c and to identify typeof(medium)
regular_basis_function(p::Particle{Dim,P}, ω::T) where {Dim,T, P<:PhysicalMedium{Dim,1}} = regular_basis_function(ω / p.medium.c, p.medium)

"""
estimate_regular_basis_order(medium::P, ω::Number, radius::Number; tol = 1e-6)
"""
Expand Down
141 changes: 0 additions & 141 deletions src/physics/acoustics/acoustics.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,14 +11,6 @@ struct Acoustic{T,Dim} <: PhysicalMedium{Dim,1}
c::Complex{T} # Phase velocity
end

# basisorder_to_linearindices(::Type{Acoustic{T,3}}, order::Int) where T = (order^2 + 1):(order+1)^2
# basisorder_to_linearindices(::Type{Acoustic{T,2}}, order::Int) where T = 1:(2*order + 1)
basisorder_to_basislength(::Type{Acoustic{T,3}}, order::Int) where T = (order+1)^2
basisorder_to_basislength(::Type{Acoustic{T,2}}, order::Int) where T = 2*order + 1

basislength_to_basisorder(::Type{Acoustic{T,3}},len::Int) where T = Int(sqrt(len) - 1)
basislength_to_basisorder(::Type{Acoustic{T,2}},len::Int) where T = Int(T(len - 1) / T(2.0))

# Constructor which supplies the dimension without explicitly mentioning type
Acoustic(ρ::T,c::Union{T,Complex{T}},Dim::Integer) where {T<:Number} = Acoustic{T,Dim}(ρ,Complex{T}(c))
Acoustic(Dim::Integer; ρ::T = 0.0, c::Union{T,Complex{T}} = 0.0) where {T<:Number} = Acoustic{T,Dim}(ρ,Complex{T}(c))
Expand All @@ -43,139 +35,6 @@ Characteristic specific acoustic impedance (z₀) of medium
"""
impedance(medium::Acoustic) = medium.ρ * medium.c

function outgoing_radial_basis(medium::Acoustic{T,2}, ω::T, order::Integer, r::T) where {T<:Number}
k = ω/medium.c
return hankelh1.(-order:order,k*r)
end

function outgoing_basis_function(medium::Acoustic{T,2}, ω::T) where {T<:Number}
return function (order::Integer, x::AbstractVector{T})
r, θ = cartesian_to_radial_coordinates(x)
k = ω/medium.c
[hankelh1(m,k*r)*exp(im*θ*m) for m = -order:order]
end
end

function outgoing_radial_basis(medium::Acoustic{T,3}, ω::T, order::Integer, r::T) where {T<:Number}
k = ω/medium.c
hs = shankelh1.(0:order,k*r)
return [hs[l+1] for l = 0:order for m = -l:l]
end

function outgoing_basis_function(medium::Acoustic{T,3}, ω::T) where {T<:Number}
return function (order::Integer, x::AbstractVector{T})
r, θ, φ = cartesian_to_radial_coordinates(x)
k = ω/medium.c

Ys = spherical_harmonics(order, θ, φ)
hs = [shankelh1(l,k*r) for l = 0:order]

lm_to_n = lm_to_spherical_harmonic_index

return [hs[l+1] * Ys[lm_to_n(l,m)] for l = 0:order for m = -l:l]
end
end

function outgoing_translation_matrix(medium::Acoustic{T,2}, in_order::Integer, out_order::Integer, ω::T, x::AbstractVector{T}) where {T<:Number}
translation_vec = outgoing_basis_function(medium, ω)(in_order + out_order, x)
U = [
translation_vec[n-m + in_order + out_order + 1]
for n in -out_order:out_order, m in -in_order:in_order]

return U
end

function regular_translation_matrix(medium::Acoustic{T,2}, in_order::Integer, out_order::Integer, ω::T, x::AbstractVector{T}) where {T<:Number}
translation_vec = regular_basis_function(medium, ω)(in_order + out_order, x)
V = [
translation_vec[n-m + in_order + out_order + 1]
for n in -out_order:out_order, m in -in_order:in_order]

return V
end

function outgoing_translation_matrix(medium::Acoustic{T,3}, in_order::Integer, out_order::Integer, ω::T, x::AbstractVector{T}) where {T<:Number}

us = outgoing_basis_function(medium, ω)(in_order + out_order,x)
c = gaunt_coefficient

ind(order::Int) = basisorder_to_basislength(Acoustic{T,3},order)
U = [
begin
i1 = abs(l-dl) == 0 ? 1 : ind(abs(l-dl)-1) + 1
i2 = ind(l+dl)

cs = [c(T,l,m,dl,dm,l1,m1) for l1 = abs(l-dl):(l+dl) for m1 = -l1:l1]
sum(us[i1:i2] .* cs)
end
for dl = 0:in_order for dm = -dl:dl for l = 0:out_order for m = -l:l];
# U = [
# [(l,m),(dl,dm)]
# for dl = 0:order for dm = -dl:dl for l = 0:order for m = -l:l]

U = reshape(U, ((out_order+1)^2, (in_order+1)^2))

return U
end

# NOTE that medium in both functions below is only used to get c and to identify typeof(medium)
regular_basis_function(p::Particle{Dim,Acoustic{T,Dim}}, ω::T) where {Dim,T} = regular_basis_function(ω/p.medium.c, p.medium)
regular_basis_function(medium::Acoustic{T,Dim}, ω::Union{T,Complex{T}}) where {Dim,T} = regular_basis_function(ω/medium.c, medium)

function regular_radial_basis(medium::Acoustic{T,2}, ω::T, order::Integer, r::T) where {T<:Number}
k = ω/medium.c
return besselj.(-order:order,k*r)
end

function regular_basis_function(wavenumber::Union{T,Complex{T}}, ::Acoustic{T,2}) where T
return function (order::Integer, x::AbstractVector{T})
r, θ = cartesian_to_radial_coordinates(x)
k = wavenumber

return [besselj(m,k*r)*exp(im*θ*m) for m = -order:order]
end
end

function regular_radial_basis(medium::Acoustic{T,3}, ω::T, order::Integer, r::T) where {T<:Number}
k = ω / medium.c
js = sbesselj.(0:order,k*r)

return [js[l+1] for l = 0:order for m = -l:l]
end

function regular_basis_function(wavenumber::Union{T,Complex{T}}, ::Acoustic{T,3}) where T
return function (order::Integer, x::AbstractVector{T})
r, θ, φ = cartesian_to_radial_coordinates(x)

Ys = spherical_harmonics(order, θ, φ)
js = [sbesselj(l,wavenumber*r) for l = 0:order]

lm_to_n = lm_to_spherical_harmonic_index

return [js[l+1] * Ys[lm_to_n(l,m)] for l = 0:order for m = -l:l]
end
end

function regular_translation_matrix(medium::Acoustic{T,3}, in_order::Integer, out_order::Integer, ω::T, x::AbstractVector{T}) where {T<:Number}
vs = regular_basis_function(medium, ω)(in_order + out_order,x)
c = gaunt_coefficient

ind(order::Int) = basisorder_to_basislength(Acoustic{T,3},order)
V = [
begin
i1 = (abs(l-dl) == 0) ? 1 : ind(abs(l-dl)-1) + 1
i2 = ind(l+dl)

cs = [c(T,l,m,dl,dm,l1,m1) for l1 = abs(l-dl):(l+dl) for m1 = -l1:l1]
sum(vs[i1:i2] .* cs)
end
for dl = 0:in_order for dm = -dl:dl for l = 0:out_order for m = -l:l];

V = reshape(V, ((out_order+1)^2, (in_order+1)^2))

return V
end

# Check for material properties that don't make sense or haven't been implemented
"""
check_material(p::Particle, outer_medium::Acoustic)
Expand Down
159 changes: 157 additions & 2 deletions src/physics/physical_medium.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
PhysicalMedium{Dim,FieldDim}
An abstract type used to represent the physical medium, the dimension of
the field, and the number of spatial dimensions.
the field, and the number of spatial dimensions. We expect every sub-type of PhysicalMedium{Dim,1} to have a field c which is the complex wave speed.
"""
abstract type PhysicalMedium{Dim,FieldDim} end

Expand All @@ -18,11 +18,26 @@ spatial_dimension(p::PhysicalMedium{Dim,FieldDim}) where {Dim,FieldDim} = Dim
"Extract the dimension of the space that this type of physical property lives in"
spatial_dimension(p::Type{P}) where {Dim,FieldDim,P<:PhysicalMedium{Dim,FieldDim}} = Dim

"""
ScalarMedium{Dim}
An type used to represent a scalar wave which satisfies a Helmoholtz equations. Dim represents the number of spatial dimensions.
"""
struct ScalarMedium{T,Dim} <: PhysicalMedium{Dim,1}
c::Complex{T} # complex wavespeed
end


"""
A basis for regular functions, that is, smooth functions. A series expansion in this basis should converge to any regular function within a ball.
"""
regular_basis_function

"""
Basis of outgoing wave when the dependance on the angles has been removed.
"""
outgoing_radial_function

"""
Basis of outgoing wave. A series expansion in this basis should converge to any scattered field outside of a ball which contains the scatterer.
"""
Expand All @@ -38,8 +53,148 @@ A tuples of vectors of the field close to the boundary of the shape. The field i
"""
boundary_data

## Below are many of the essential functions for scalar wave potentials (i.e. PhysicalMedium{Dim,1}) which are used to build all other types of vector waves (as we strict ourselves to isotropy.)

# basisorder_to_linearindices(::Type{Acoustic{T,3}}, order::Int) where T = (order^2 + 1):(order+1)^2
# basisorder_to_linearindices(::Type{Acoustic{T,2}}, order::Int) where T = 1:(2*order + 1)
basisorder_to_basislength(::Type{P}, order::Int) where P <: PhysicalMedium{3,1} = (order+1)^2
basisorder_to_basislength(::Type{P}, order::Int) where P <: PhysicalMedium{2,1} = 2*order + 1

basislength_to_basisorder(::Type{P},len::Int) where P <: PhysicalMedium{3,1} = Int(sqrt(len) - 1)
basislength_to_basisorder(::Type{P},len::Int) where P <: PhysicalMedium{2,1} = Int((len - 1) / 2.0)

function outgoing_radial_basis(medium::PhysicalMedium{2,1}, ω::T, order::Integer, r::T) where {T<:Number}
k = ω/medium.c
return hankelh1.(-order:order,k*r)
end

function outgoing_basis_function(medium::PhysicalMedium{2,1}, ω::T) where {T<:Number}
return function (order::Integer, x::AbstractVector{T})
r, θ = cartesian_to_radial_coordinates(x)
k = ω/medium.c
[hankelh1(m,k*r)*exp(im*θ*m) for m = -order:order]
end
end

function outgoing_radial_basis(medium::PhysicalMedium{3,1}, ω::T, order::Integer, r::T) where {T<:Number}
k = ω/medium.c
hs = shankelh1.(0:order,k*r)
return [hs[l+1] for l = 0:order for m = -l:l]
end

function outgoing_basis_function(medium::PhysicalMedium{3,1}, ω::T) where {T<:Number}
return function (order::Integer, x::AbstractVector{T})
r, θ, φ = cartesian_to_radial_coordinates(x)
k = ω/medium.c

Ys = spherical_harmonics(order, θ, φ)
hs = [shankelh1(l,k*r) for l = 0:order]

lm_to_n = lm_to_spherical_harmonic_index

return [hs[l+1] * Ys[lm_to_n(l,m)] for l = 0:order for m = -l:l]
end
end

function outgoing_translation_matrix(medium::PhysicalMedium{2,1}, in_order::Integer, out_order::Integer, ω::T, x::AbstractVector{T}) where {T<:Number}
translation_vec = outgoing_basis_function(medium, ω)(in_order + out_order, x)
U = [
translation_vec[n-m + in_order + out_order + 1]
for n in -out_order:out_order, m in -in_order:in_order]

return U
end

function regular_translation_matrix(medium::PhysicalMedium{2,1}, in_order::Integer, out_order::Integer, ω::T, x::AbstractVector{T}) where {T<:Number}
translation_vec = regular_basis_function(medium, ω)(in_order + out_order, x)
V = [
translation_vec[n-m + in_order + out_order + 1]
for n in -out_order:out_order, m in -in_order:in_order]

return V
end

function outgoing_translation_matrix(medium::PhysicalMedium{3,1}, in_order::Integer, out_order::Integer, ω::T, x::AbstractVector{T}) where {T<:Number}

us = outgoing_basis_function(medium, ω)(in_order + out_order,x)
c = gaunt_coefficient

ind(order::Int) = basisorder_to_basislength(Acoustic{T,3},order)
U = [
begin
i1 = abs(l-dl) == 0 ? 1 : ind(abs(l-dl)-1) + 1
i2 = ind(l+dl)

cs = [c(T,l,m,dl,dm,l1,m1) for l1 = abs(l-dl):(l+dl) for m1 = -l1:l1]
sum(us[i1:i2] .* cs)
end
for dl = 0:in_order for dm = -dl:dl for l = 0:out_order for m = -l:l];
# U = [
# [(l,m),(dl,dm)]
# for dl = 0:order for dm = -dl:dl for l = 0:order for m = -l:l]

U = reshape(U, ((out_order+1)^2, (in_order+1)^2))

return U
end

# NOTE that medium in both functions below is only used to get c and to identify typeof(medium)
regular_basis_function(medium::PhysicalMedium{Dim,1}, ω::Union{T,Complex{T}}) where {Dim,T} = regular_basis_function(ω/medium.c, medium)

function regular_radial_basis(medium::PhysicalMedium{2,1}, ω::T, order::Integer, r::T) where {T<:Number}
k = ω / medium.c
return besselj.(-order:order,k*r)
end

function regular_basis_function(wavenumber::Union{T,Complex{T}}, ::PhysicalMedium{2,1}) where T
return function (order::Integer, x::AbstractVector{T})
r, θ = cartesian_to_radial_coordinates(x)
k = wavenumber

return [besselj(m,k*r)*exp(im*θ*m) for m = -order:order]
end
end

function regular_radial_basis(medium::PhysicalMedium{3,1}, ω::T, order::Integer, r::T) where {T<:Number}
k = ω / medium.c
js = sbesselj.(0:order,k*r)

return [js[l+1] for l = 0:order for m = -l:l]
end

function regular_basis_function(wavenumber::Union{T,Complex{T}}, ::PhysicalMedium{3,1}) where T
return function (order::Integer, x::AbstractVector{T})
r, θ, φ = cartesian_to_radial_coordinates(x)

Ys = spherical_harmonics(order, θ, φ)
js = [sbesselj(l,wavenumber*r) for l = 0:order]

lm_to_n = lm_to_spherical_harmonic_index

return [js[l+1] * Ys[lm_to_n(l,m)] for l = 0:order for m = -l:l]
end
end

function regular_translation_matrix(medium::PhysicalMedium{3,1}, in_order::Integer, out_order::Integer, ω::T, x::AbstractVector{T}) where {T<:Number}
vs = regular_basis_function(medium, ω)(in_order + out_order,x)
c = gaunt_coefficient

ind(order::Int) = basisorder_to_basislength(Acoustic{T,3},order)
V = [
begin
i1 = (abs(l-dl) == 0) ? 1 : ind(abs(l-dl)-1) + 1
i2 = ind(l+dl)

cs = [c(T,l,m,dl,dm,l1,m1) for l1 = abs(l-dl):(l+dl) for m1 = -l1:l1]
sum(vs[i1:i2] .* cs)
end
for dl = 0:in_order for dm = -dl:dl for l = 0:out_order for m = -l:l];

V = reshape(V, ((out_order+1)^2, (in_order+1)^2))

return V
end

# estimate_regular_basisorder(medium::P, ka) where P<:PhysicalMedium = estimate_regular_basisorder(P, ka)


"""
Expand Down
3 changes: 1 addition & 2 deletions test/end_to_end_tests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,8 +17,7 @@

result = run(sim, [1.0,2.0], 0.1)
result = run(particles, source, [1.0,2.0], 0.1)
result = 3.2*result + result*4.0im + 0.3+4.0im # changes result.field
3.2 + result
result = 3.2*result + result*4.0im + [0.3+4.0im] # changes result.field

@test field(result)[1] == result.field[1][1] # returns
@test field(result,1,1) == result.field[1][1] # returns
Expand Down

0 comments on commit 4eb9c4c

Please sign in to comment.