Merge branch 'piva'

Merging branch with master after changes

src/physics/acoustics/helmholtz_circle.jl

@@ -1,7 +1,7 @@
 """
     t_matrix(Particle{2,Acoustic{T,2},Sphere{T,2}}, Acoustic{T,2}, ω, order)
 
-The T-matrix for a 2D circlular Helmholtz resonator in a 2D acoustic medium.
+The T-matrix for a 2D isotropic Helmholtz resonator in a 2D acoustic medium.
 """
 function t_matrix(p::Particle{2,Acoustic{T,2},IsotropicHelmholtz{T,2}}, outer_medium::Acoustic{T,2}, ω::T, basis_order::Integer)::Diagonal{Complex{T}} where T <: AbstractFloat
 
@@ -37,3 +37,58 @@ function t_matrix(p::Particle{2,Acoustic{T,2},IsotropicHelmholtz{T,2}}, outer_me
 
     return - Diagonal(vcat(reverse(Zns), Zns[2:end]))
 end
+
+"""
+    t_matrix(Particle{2,Acoustic{T,2},Sphere{T,2}}, Acoustic{T,2}, ω, order)
+
+The T-matrix for a 2D non-isotropic Helmholtz resonator in a 2D acoustic medium.
+"""
+function t_matrix(p::Particle{2,Acoustic{T,2},Helmholtz{T,2}}, outer_medium::Acoustic{T,2}, ω::T, basis_order::Integer)::Matrix{Complex{T}} where T <: AbstractFloat
+
+    M = basis_order
+    l = p.shape.aperture
+    γe = Base.MathConstants.eulergamma
+    k = ω/outer_medium.c
+    b = p.shape.radius
+    a = p.shape.inner_radius
+    θ0 = p.shape.orientation
+
+    if a != b
+        @warn "Thick-walled Helmholtz resonator not implemented, using thin-walled instead."
+    end
+
+    kb = k*b
+    Qsum = sum([(besselj(m, kb)*diffhankelh1(m, kb) + diffbesselj(m, kb)*hankelh1(m, kb))^2/(diffbesselj(m, kb)*diffhankelh1(m, kb)) for m in -2M:2M])
+    hε = (4im / pi) * (γe - 1im*pi/2 + log(k*l/4)) - Qsum / 2
+
+    # Check for material properties that don't mkbe sense or haven't been implemented
+    check_material(p, outer_medium)
+
+    if p.medium.ρ < Inf || real(p.medium.c) < Inf
+        @warn "Theory not done for penetrable Helmholtz resonators, using Newmann boundaries instead."
+    end
+
+    function δ(x,y)
+        if x == y
+            return 1
+        else
+            return 0
+        end
+    end
+
+    "Returns a ratio used in multiple scattering which reflects the material properties of the particles"
+    function Zn(m::Integer,n::Integer)::Complex{T}
+
+        # set the scattering strength and type
+        numer = diffbesselj(m, kb)
+        denom = diffhankelh1(m, kb)
+        resonance_term = 2 * exp(-1im * (m - n) * θ0) / (hε * (pi * kb)^2 * diffhankelh1(m, kb) * diffhankelh1(n, kb))
+
+        return -(numer / denom) * δ(n,m) + resonance_term
+    end
+
+    # Get Zns
+    Zns = [Zn(i, j) for i in -M:M, j in -M:M]
+
+    return Zns
+end

src/random/random_particles.jl

@@ -93,7 +93,7 @@ function random_particles(particle_medium::P, particle_shape::S, region_shape::S
         max_attempts_to_place_particle::Int = 3000, # Maximum number of attempts to place a particle
         current_particles::AbstractParticles{Dim} = AbstractParticle{Dim}[] # Particles already present.
 ) where {Dim,P<:PhysicalMedium{Dim},S<:Shape{Dim}}
-
+    
     # Check volume fraction is not impossible
     volfrac = N * volume(particle_shape) / volume(region_shape)
     max_packing = if length(current_particles) > 0
@@ -135,7 +135,12 @@ function random_particles(particle_medium::P, particle_shape::S, region_shape::S
             outside_box = true
             while outside_box
                 x = (box.dimensions ./ 2) .* (1 .- 2 .* rand(randgen,typeof(origin(box)))) + origin(box)
-                particles[n] = Particle(particle_medium, congruent(particle_shape, x))
+                if S == Helmholtz{Float64, 2}
+                    θ = 2π*rand(randgen)
+                    particles[n] = Particle(particle_medium, congruent(particle_shape, x, θ))
+                else
+                    particles[n] = Particle(particle_medium, congruent(particle_shape, x))
+                end
                 outside_box = !(particles[n] ⊆ region_shape)
             end
 

src/shapes/sphere.jl

@@ -173,6 +173,10 @@ function congruent(s::IsotropicHelmholtz, x)
     IsotropicHelmholtz(x, s.radius, s.aperture)
 end
 
+function congruent(s::Helmholtz, x, θ)
+    Helmholtz(x, s.radius, s.inner_radius, s.aperture, θ)
+end
+
 function Circle(sphere::AbstractSphere; y = sphere.origin[2])
     if abs(y - sphere.origin[2]) > sphere.radius
         return EmptyShape(sphere)