adding more properties to SphericalHelmholtz

docs/src/example/random_particles/README.md

@@ -1,5 +1,7 @@
 # Simple random particles example
 
+```
+
 ## Define particle properties
 Define the volume fraction of particles, the region to place the particles, and their radius
 ```julia
@@ -33,25 +35,25 @@ result = run(simulation, x, ωs)
 We use the `Plots` package to plot both the response at the listener position x
 
 ```julia
-using Plots; #pyplot(linewidth = 2.0)
-plot(result, field_apply=real) # plot result
-plot!(result, field_apply=imag)
-#savefig("plot_result.png")
+    using Plots; #pyplot(linewidth = 2.0)
+    plot(result, field_apply=real) # plot result
+    plot!(result, field_apply=imag)
+    #savefig("plot_result.png")
 ```
 ![Plot of response against wavenumber](plot_result.png)
 
 And plot the whole field inside the region_shape `bounds` for a specific wavenumber (`ω=0.8`)
 ```julia
-bottomleft = [-15.,-max_width]
-topright = [max_width,max_width]
-bounds = Box([bottomleft,topright])
+    bottomleft = [-15.,-max_width]
+    topright = [max_width,max_width]
+    bounds = Box([bottomleft,topright])
 
-#plot(simulation,0.8; res=80, bounds=bounds)
-#plot!(region_shape, linecolor=:red)
-#plot!(simulation)
-#scatter!([x[1]],[x[2]], lab="receiver")
+    #plot(simulation,0.8; res=80, bounds=bounds)
+    #plot!(region_shape, linecolor=:red)
+    #plot!(simulation)
+    #scatter!([x[1]],[x[2]], lab="receiver")
 
-#savefig("plot_field.png")
+    #savefig("plot_field.png")
 ```
 ![Plot real part of acoustic field](plot_field.png)
 ## Things to try

docs/src/maths/capsule-boundary-conditions.wl

@@ -255,13 +255,47 @@
 
 
 (* ::Input:: *)
-(*subN = {a[1]->2.0,a[0]->1.5,\[Rho]o->0.00001,\[Rho][1]->2.0,ko-> 0.0001,k[1] -> 2.0,Subscript[J, n][ko a[1]]->1.0,fo[n]->1, Subscript[H, n_][x_] ->HankelH1[n,x], Subscript[J, n_][x_] ->BesselJ[n,x], Derivative[1][Subscript[H, n_]][x_] ->D[HankelH1[n,x],x],Derivative[1][Subscript[J, n_]][x_] ->D[BesselJ[n,x],x]};*)
-(*subsol//.subN//Flatten;*)
+(*subN = {a[1]->2.0,a[0]->1.5,\[Rho]o->0.00001,\[Rho][1]->2.0,ko-> 0.0001,k[1] -> 2.0,Subscript[J, n][ko a[1]]->1.0, Subscript[H, n_][x_] ->HankelH1[n,x], Subscript[J, n_][x_] ->BesselJ[n,x], Derivative[1][Subscript[H, n_]][x_] ->D[HankelH1[n,x],x],Derivative[1][Subscript[J, n_]][x_] ->D[BesselJ[n,x],x]};*)
+(**)
+(**)
+(*\[Epsilon] = 10^-12;*)
+(*ns = Range[0,30,3];*)
+(**)
+(*subNsol1 = subsol//.subN/.{fo[n] -> 1 + RandomReal[\[Epsilon]] +RandomReal[\[Epsilon]] I} //Flatten;*)
+(*subNsol2 = subsol//.subN/.{fo[n] -> 1 + RandomReal[\[Epsilon]] +RandomReal[\[Epsilon]] I} //Flatten;*)
+(**)
+(**)
+(*Nsol1= Transpose[#/.Rule-> List][[2]]&/@(subNsol1/.n->#&/@ns);*)
+(*Nsol2= Transpose[#/.Rule-> List][[2]]&/@(subNsol2/.n->#&/@ns);*)
+(**)
+(*(*The difference appears to be substantial, but it is only the coefficient multiplying Jn which gets larger errors. And as Jn becomes very small the result may not be a large error in the field*)*)
+(*(Nsol1 - Nsol2)*)
 
 
 (* ::Input:: *)
-(*%/.n->#&/@Range[1,20,3]*)
+(*(*The scattered wave appears to be stable*)*)
+(**)
+(*#[[1]]&/@(Nsol1 - Nsol2);*)
+(*Abs[% *( HankelH1[#,ko a[1]//.subN]&/@ns)]*)
+(**)
+(**)
+(*(*The internal besselJ wave*)*)
+(*#[[2]]&/@(Nsol1 - Nsol2);*)
+(*Abs[% *( BesselJ[#,k[1] a[1]//.subN]&/@ns)]*)
 (**)
+(*(*The internal HankelH1 wave*)*)
+(*#[[3]]&/@(Nsol1 - Nsol2);*)
+(*Abs[% *( HankelH1[#,k[1] a[1]//.subN]&/@ns)]*)
+(**)
+
+
+(* ::Input:: *)
+(*Norm/@(Nsol1 - Nsol2)*)
+(**)
+
+
+(* ::Input:: *)
+(*subsol*)
 
 
 (* ::Input:: *)
@@ -281,16 +315,15 @@
 
 (* ::Input:: *)
 (*subN = {a[1]->2.0,a[0]->1.5,k[1] -> 2.0,Subscript[J, n][ko a[1]]->1.0,fo[n]->1, Subscript[H, n_][x_] ->HankelH1[n,x], Subscript[J, n_][x_] ->BesselJ[n,x]};*)
-(*ns = Range[1,30,3];*)
+(*ns = Range[1,40,3];*)
 (*subsol//.subN//Flatten;*)
 (*subNsol=Flatten[%/.n->#&/@ns];*)
 (**)
 (*eqs//.subN/.n->#&/@ns;*)
-(*%/.subNsol*)
+(*Norm/@(%/.subNsol)*)
 
 
 (* ::Input:: *)
-(**)
 (*{Abs@f[1,n]Abs@BesselJ[n,a[1] k[1]],Abs@A[1,n]Abs@HankelH1[n,a[1] k[1]]}//.subN//Flatten;*)
 (*%/.n->#&/@ns;*)
 (*%/.Flatten@subNsol*)
@@ -301,12 +334,10 @@
 (*(*Is the matrix system ill-posed?*)*)
 (*NM = M //.subN;*)
 (*Nb = b //.subN;*)
-(*ns = Range[1,70,6];*)
 (*{Abs@Det[NM/.n->#],SingularValueList[NM/.n->#]}&/@ns*)
 (**)
 (*(*If we numerically solve the system:*)*)
-(*sols =Flatten[ Thread[(vars/.subN/.n->#)->LinearSolve[NM/.n->#,Nb/.n->#]]&/@ns];*)
-(**)
+(*sols =Flatten[ Thread[(vars/.subN/.n->#)->LinearSolve[NM/.n->#,Nb/.n->#]]&/@ns]//Quiet;*)
 (**)
 
 
@@ -317,22 +348,37 @@
 
 
 (* ::Input:: *)
-(*subsol*)
-(**)
+(*invM = Inverse[M]*)
+(*invM . b*)
+(*NinvM =invM //.subN;*)
+(*Nb = b //.subN;*)
 
 
-(* ::InheritFromParent:: *)
+(* ::Input:: *)
+(*invsols =Flatten[ Thread[(vars/.subN/.n->#)->(-NinvM . (Nb + RandomReal[{10^-18,10^-15}])/.n->#)]&/@ns];*)
+(*Norm/@(eqs//.subN/.n->#&/@ns/.invsols)*)
 (**)
 
 
-(* ::InheritFromParent:: *)
+(* ::Input:: *)
 (**)
 
 
-Test a capsule filled with air
+(* ::Input:: *)
+(*eqs//.subN/.n->#&/@ns/.invsols*)
+(**)
+
 
+(* ::Input:: *)
+(*Nx = -NinvM . Nb/.n-> 1*)
+(*M . Nx//.subN/.n-> 1*)
+(*b//.subN/.n-> 1*)
+(**)
+(**)
 
 
+(* ::Input:: *)
+(*#->#&/@ns*)
 
 
 (* ::Input:: *)
@@ -494,6 +540,10 @@ the background material.
 
 
 
+
+
+
+
 (* ::Section:: *)
 (*Addition Theorems*)
 

src/shapes/sphere.jl

@@ -24,16 +24,28 @@ A [`Shape`](@ref) which represents a 2D thin-walled isotropic Helmholtz resonato
 struct SphericalHelmholtz{T,Dim} <: AbstractSphere{Dim}
     origin::SVector{Dim,T}
     radius::T
+    inner_radius::T
     aperture::T
+    orientation::T
 end
 
 # Alternate constructors, where type is inferred naturally
 Sphere(origin::NTuple{Dim}, radius::T) where {T,Dim} = Sphere{T,Dim}(origin, radius)
 Sphere(origin::AbstractVector, radius::T) where {T} = Sphere{T,length(origin)}(origin, radius)
-SphericalHelmholtz(origin::NTuple{2}, radius::T, aperture::T) where {T} = SphericalHelmholtz{T,2}(origin, radius, aperture)
+SphericalHelmholtz(origin::NTuple{Dim}, radius::T, inner_radius::T, aperture::T, orientation::T) where {T, Dim} = SphericalHelmholtz{T,Dim}(origin, radius, inner_radius, aperture, orientation)
+
+SphericalHelmholtz(origin::NTuple{Dim}, radius::T; inner_radius::T = zero(T), aperture::T = zero(T), orientation::T = zero(T)) where {T, Dim} = SphericalHelmholtz{T,Dim}(origin, radius, inner_radius, aperture, orientation)
 
 Sphere(Dim, radius::T) where {T} = Sphere{T,Dim}(zeros(T,Dim), radius)
-SphericalHelmholtz(radius::T, aperture::T) where {T} = SphericalHelmholtz{T,2}(zeros(T,2), radius, aperture)
+
+function SphericalHelmholtz(Dim::Int, radius::T; 
+        origin::T = zeros(T,Dim), 
+        aperture::T = zero(T), 
+        inner_radius::T = zero(T),
+        orientation::T = zero(T)) where {T} 
+
+    return SphericalHelmholtz{T,Dim}(origin, radius, inner_radius, aperture, orientation)
+end
 
 Circle(origin::Union{AbstractVector{T},NTuple{2,T}}, radius::T) where T <: AbstractFloat = Sphere{T,2}(origin, radius::T)
 Circle(radius::T) where T <: AbstractFloat = Sphere(2, radius::T)
@@ -87,7 +99,7 @@ function ==(c1::Sphere, c2::Sphere)
 end
 
 function ==(c1::SphericalHelmholtz, c2::SphericalHelmholtz)
-    c1.origin == c2.origin && c1.radius == c2.radius && c1.aperture == c2.aperture
+    c1.origin == c2.origin && c1.radius == c2.radius && c1.inner_radius == c2.inner_radius && c1.aperture == c2.aperture && c1.orientation == c2.orientation
 end
 
 import Base.isequal
@@ -112,7 +124,7 @@ function iscongruent(s1::Sphere, s2::Sphere)
 end
 
 function iscongruent(s1::SphericalHelmholtz, s2::SphericalHelmholtz)
-    s1.radius == s2.radius && s1.aperture == s2.aperture
+    s1.radius == s2.radius && s1.aperture == s2.aperture && s1.inner_radius == s2.inner_radius && s1.orientation == s2.orientation
 end
 
 function congruent(s::Sphere, x)
@@ -136,7 +148,7 @@ end
 
 
 
-function boundary_functions(sphere::AbstractSphere{3})
+function boundary_functions(sphere::Sphere{T,3}) where T
 
     function x(t,s=0)
         check_boundary_coord_range(t)
@@ -162,7 +174,22 @@ function boundary_functions(sphere::AbstractSphere{3})
     return x, y, z
 end
 
-function boundary_functions(circle::AbstractSphere{2})
+function boundary_functions(circle::Sphere{T,2}) where T
+
+    function x(t)
+        check_boundary_coord_range(t)
+        circle.radius * cos(2π * t) + origin(circle)[1]
+    end
+
+    function y(t)
+        check_boundary_coord_range(t)
+        circle.radius * sin(2π * t) + origin(circle)[2]
+    end
+
+    return x, y
+end
+
+function boundary_functions(circle::SphericalHelmholtz{T,2}) where T
 
     function x(t)
         check_boundary_coord_range(t)

test/particle_tests.jl

@@ -53,10 +53,12 @@
         circle_identical = Sphere((1.0,3.0),2.0)
         circle_congruent = Sphere((4.0,7.0),2.0)
         rect = Box((1.0,2.0),(2.0,4.0))
-        resonator = SphericalHelmholtz((1.0,3.0),2.0, 0.2)
-        resonator_dif_aperture = SphericalHelmholtz((1.0,3.0),2.0, 0.1)
-        resonator_identical = SphericalHelmholtz((1.0,3.0), 2.0, 0.2)
-        resonator_congruent = SphericalHelmholtz((4.0,7.0), 2.0, 0.2)
+
+        resonator = SphericalHelmholtz((1.0,3.0),2.0, 0.2, 0.01, -1.3)
+        resonator_kws = SphericalHelmholtz((1.0,3.0),2.0; inner_radius = 0.2, aperture = 0.1, orientation = -1.3)
+        resonator_dif_aperture = SphericalHelmholtz((1.0,3.0),2.0; aperture = 0.1)
+        resonator_identical = SphericalHelmholtz((1.0,3.0), 2.0; orientation = -1.3, inner_radius = 0.2, aperture = 0.01)
+        resonator_congruent = SphericalHelmholtz((4.0,7.0), 2.0; orientation = -1.3, inner_radius = 0.2, aperture = 0.01)
 
         # Construct four particles, with two the same
         p = Particle(a2,circle)