Merge pull request #23 from sbadrian/feature/improved_layered_sphere_…

…support_2

Extend permeability and permittivity function to LayeredSphere

src/SphericalScattering.jl

@@ -45,6 +45,8 @@ export DielectricSphereThinImpedanceLayer
 export field, scatteredfield
 export rcs
 export sphericalGridPoints, phiCutPoints, thetaCutPoints
+export numlayers, layer
+export permittivity, permeability, medium
 
 
 # -------- extensions

src/sphere.jl

@@ -108,8 +108,18 @@ end
 
 Constructor for the layered dielectric sphere.
 """
-LayeredSphere(; radii=error("Missing argument `radii`"), filling=error("`missing argument `filling`")) = LayeredSphere(radii, filling)
+function LayeredSphere(; radii=error("Missing argument `radii`"), filling=error("`missing argument `filling`"))
 
+    if sort(radii) != radii
+        error("Radii are not ordered ascendingly.")
+    end
+
+    if length(radii) != length(filling)
+        error("Number of fillings does not match number of radii.")
+    end
+
+    LayeredSphere(radii, filling)
+end
 
 
 
@@ -126,9 +136,32 @@ end
 
 Constructor for the layered dielectric sphere.
 """
-LayeredSpherePEC(; radii=error("Missing argument `radii`"), filling=error("Missing argument `filling`")) =
+function LayeredSpherePEC(; radii=error("Missing argument `radii`"), filling=error("Missing argument `filling`"))
+
+    if sort(radii) != radii
+        error("Radii are not ordered ascendingly.")
+    end
+
+    if length(radii) != length(filling) + 1
+        error("Number of fillings does not match number of radii.")
+    end
+
     LayeredSpherePEC(radii, filling)
+end
+
 
+"""
+    numlayers(sp::Sphere)
+
+Returns the number of layers.
+"""
+function numlayers(sp::Sphere)
+    return 2
+end
+
+function numlayers(sp::Union{LayeredSphere,LayeredSpherePEC})
+    return length(sp.radii) + 1
+end
 
 
 
@@ -149,6 +182,38 @@ function layer(sp::Sphere, r)
     end
 end
 
+
+
+"""
+    layer(sp::LayeredSphere, r)
+
+Returns the index of the layer, `r` is located, where `1` denotes the inner most layer. 
+"""
+function layer(sp::Union{LayeredSphere,LayeredSpherePEC}, r)
+    # Using Jin's numbering from the multi-layered cartesian
+    # in anticipation. For PEC and dielectric sphere;
+    # 1 = interior, 2 = exterior
+
+    r < 0.0 && error("The radius must be a positive number.")
+
+    N = numlayers(sp)
+
+    for i in 1:(N - 1)
+        if r < sp.radii[i]
+            return i # Convention: Boundary belongs to outer layer
+        end
+    end
+
+    return N
+end
+
+
+"""
+wavenumber(sp::Sphere, ex::Excitation, r)
+
+Returns the wavenumber at radius `r` in the sphere `sp`.
+If this part is PEC, k=0 is rweturned.
+"""
 function wavenumber(sp::PECSphere, ex::Excitation, r)
     ε = ex.embedding.ε
     μ = ex.embedding.μ
@@ -178,6 +243,12 @@ function wavenumber(sp::DielectricSphere, ex::Excitation, r)
     return k
 end
 
+"""
+    impedance(sp::Sphere, ex::Excitation, r)
+
+Returns the wavenumber at radius `r` in the sphere `sp`.
+If this part is PEC, a zero medium is returned.
+"""
 function impedance(sp::DielectricSphere, ex::Excitation, r)
     if layer(sp, r) == 2
         ε = ex.embedding.ε
@@ -190,6 +261,12 @@ function impedance(sp::DielectricSphere, ex::Excitation, r)
     return sqrt(μ / ε)
 end
 
+"""
+    impedance(sp::Sphere, r)
+
+Returns the impedance of the sphere `sp` at radius `r`.
+If this part is PEC, a zero medium is returned.
+"""
 function impedance(sp::PECSphere, ex::Excitation, r)
     if layer(sp, r) == 2
         ε = ex.embedding.ε
@@ -200,3 +277,106 @@ function impedance(sp::PECSphere, ex::Excitation, r)
 
     return sqrt(μ / ε)
 end
+
+function medium(sp::LayeredSphere, ex::Excitation, r)
+    N = numlayers(sp) # Number of interior layers
+
+    if layer(sp, r) == N # Outer layer has largest index
+        return ex.embedding
+    else
+        return sp.filling[layer(sp, r)]
+    end
+end
+
+
+
+function medium(sp::LayeredSpherePEC, ex::Excitation, r)
+    N = numlayers(sp) # Number of interior layers
+
+    if layer(sp, r) == N # Outer layer has largest index
+        return ex.embedding
+    elseif layer(sp, r) == 1
+        Z = promote_type(typeof(sp.filling[end].ε), typeof(sp.filling[end].μ), typeof(ex.embedding.ε), typeof(ex.embedding.μ))(0.0)
+        return Medium(Z, Z)
+    else
+        return sp.filling[layer(sp, r) - 1]
+    end
+end
+
+
+
+function medium(sp::Union{DielectricSphere,DielectricSphereThinImpedanceLayer}, ex::Excitation, r)
+    if layer(sp, r) == 2
+        return ex.embedding
+    else
+        return sp.filling
+    end
+end
+
+
+
+"""
+    medium(sp::Sphere, ex::Excitation, r)
+
+Returns the medium of the sphere `sp` at radius `r`.
+If this part is PEC, a zero medium is returned. The
+argument `r` may be a vector of positions.
+"""
+function medium(sp::PECSphere, ex::Excitation, r)
+    if layer(sp, r) == 2
+        return ex.embedding
+    else
+        Z = promote_type(typeof(ex.embedding.ε), typeof(ex.embedding.μ))(0.0)
+        return Medium(Z, Z)
+    end
+end
+
+"""
+    permittivity(sp::Sphere, ex::Excitation, r)
+
+Returns the permittivity of the sphere `sp` at radius `r`.
+If this part is PEC, a zero permittivity is returned. The
+argument `r` may be a vector of positions.
+"""
+function permittivity(sp::Sphere, ex::Excitation, r)
+    md = medium(sp, ex, r)
+    return md.ε
+end
+
+"""
+    permeability(sp::Sphere, ex::Excitation, r)
+
+Returns the permeability of the sphere `sp` at radius `r`.
+If this part is PEC, a zero permeability is returned. The
+argument `r` may be a vector of positions.
+"""
+function permeability(sp::Sphere, ex::Excitation, r)
+    md = medium(sp, ex, r)
+    return md.μ
+end
+
+function permittivity(sp::Sphere, ex::Excitation, pts::AbstractVecOrMat)
+
+    ε = permittivity(sp, ex, norm(first(pts)))
+    F = zeros(typeof(ε), size(pts))
+
+    # --- compute field in Cartesian representation
+    for (ind, point) in enumerate(pts)
+        F[ind] = permittivity(sp, ex, norm(point))
+    end
+
+    return F
+end
+
+function permeability(sp::Sphere, ex::Excitation, pts::AbstractVecOrMat)
+
+    ε = permeability(sp, ex, norm(first(pts)))
+    F = zeros(typeof(ε), size(pts))
+
+    # --- compute field in Cartesian representation
+    for (ind, point) in enumerate(pts)
+        F[ind] = permeability(sp, ex, norm(point))
+    end
+
+    return F
+end

test/sphere.jl

@@ -11,12 +11,73 @@
 
     @test Medium{Float64}(mdf) isa Medium{Float64}
 
+    pecsp = PECSphere(; radius=Float32(1.0))
+
+    @test numlayers(pecsp) == 2
+
+    ex_md = planeWave(; frequency=1e8, embedding=md)
+
+    @test medium(pecsp, ex_md, 3.0) == ex_md.embedding
+
+    @test medium(pecsp, ex_md, 0.5) == Medium(0.0, 0.0)
+
+    @test SphericalScattering.wavenumber(pecsp, ex_md, 0.0) == 0.0
+
     sp = DielectricSphere(; radius=Float32(1.0), filling=md)
 
+    @test medium(sp, ex_md, 0.5) == md
+    @test medium(sp, ex_md, 1.5) == md
+    @test medium(sp, planeWave(; frequency=1e8, embedding=mdf), 1.5) == mdf
+
     @test sp isa DielectricSphere{Float64,Float32}
 
     @test SphericalScattering.impedance(sp, planeWave(; frequency=1e8, embedding=mdf), 2.0) ≈ 0.707106 rtol = 1e-5
 
     @test SphericalScattering.wavenumber(sp, planeWave(; frequency=1e8, embedding=mdf), 2.0) ≈ 8.8857658e8 rtol = 1e-5
 
+    md1 = Medium(15.0, -2.1) # innermost medium
+    md2 = Medium(-10, -2.0) # innermost medium
+
+    @test_throws ErrorException("Radii are not ordered ascendingly.") LayeredSphere(;
+        radii=SVector(0.25, 0.5, 0.3), filling=SVector(md1, md2, md)
+    )
+    @test_throws ErrorException("Number of fillings does not match number of radii.") LayeredSphere(;
+        radii=SVector(0.25, 0.3, 0.5), filling=SVector(md1, md2)
+    )
+
+    @test_throws ErrorException("Radii are not ordered ascendingly.") LayeredSpherePEC(;
+        radii=SVector(0.25, 0.5, 0.3), filling=SVector(md1, md2)
+    )
+    @test_throws ErrorException("Number of fillings does not match number of radii.") LayeredSpherePEC(;
+        radii=SVector(0.25, 0.3, 0.5), filling=SVector(md1)
+    )
+
+    spl = LayeredSphere(; radii=SVector(0.25, 0.3, 0.5), filling=SVector(md1, md2, md))
+
+    @test numlayers(spl) == 4
+
+    @test layer(spl, 0.1) == 1
+    @test layer(spl, 0.24999) == 1
+    @test layer(spl, 0.25) == 2
+    @test layer(spl, 0.4) == 3
+    @test layer(spl, 0.5) == 4
+    @test layer(spl, 1.6) == 4
+
+    ex = planeWave(; frequency=1e8, embedding=Medium(-3.0, 2.0))
+
+    @test permittivity(spl, ex, 0.1) == 15.0
+    @test permittivity(spl, ex, 0.24999) == 15.0
+    @test permeability(spl, ex, 0.25) == -2.0
+    @test permeability(spl, ex, 1.25) == 2.0
+
+    spl = LayeredSpherePEC(; radii=SVector(0.25, 0.3, 0.5), filling=SVector(md1, md2))
+    @test medium(spl, ex_md, 1.5) == md
+
+    @test permittivity(spl, ex, 0.1) == 0.0
+    @test permittivity(spl, ex, 0.25) == 15.0
+
+    pts = [SVector(0.1, 0.0, 0.0), SVector(0.25, 0.0, 0.0)]
+
+    @test permittivity(spl, ex, pts) == [0.0, 15.0]
+    @test permeability(spl, ex, pts) == [0.0, -2.1]
 end