implement total fields computation

Project.toml

@@ -1,7 +1,7 @@
 name = "SphericalScattering"
 uuid = "1a9ea918-b599-4f1f-bd9a-d681e8bb5b3e"
 authors = ["Bernd Hofmann <[email protected]> and contributors"]
-version = "0.3.0"
+version = "0.4.0"
 
 [deps]
 LegendrePolynomials = "3db4a2ba-fc88-11e8-3e01-49c72059a882"

docs/src/dipoles.md

@@ -105,9 +105,6 @@ FF = scatteredfield(sp, ex, FarField(point_cart))
 ---
 ## Total Field
 
-!!! warning
-    Not fully implemented yet.
-
 #### API
 
 The general API is employed:

docs/src/manual.md

@@ -12,7 +12,7 @@ The basic building blocks are introduced in the following simple example; more d
 using SphericalScattering, StaticArrays
 
 # define excitation: plane wave travelling in negative z-direction with x-polarization
-ex = planeWave(wavenumber=30.0) # ≈ 10 MHz
+ex = planeWave(frequency=10e6) # Hz
 
 # define scatterer: PEC sphere
 sp = PECSphere(radius = 1.0)

docs/src/planeWave.md

@@ -51,10 +51,11 @@ The general API is employed:
 E  = field(ex, ElectricField(point_cart))
 
 H  = field(ex, MagneticField(point_cart))
-
-FF = field(ex, FarField(point_cart))
 ```
 
+!!! note
+    The far-field of a plane wave is not defined.
+
 ---
 ## Scattered Field
 
@@ -77,16 +78,14 @@ FF = scatteredfield(sp, ex, FarField(point_cart))
 ---
 ## Total Field
 
-!!! warning
-    Not fully implemented yet.
-
 #### API
 
 The general API is employed:
 ```julia
 E  = field(sp, ex, ElectricField(point_cart))
 
 H  = field(sp, ex, MagneticField(point_cart))
+```
 
-FF = field(sp, ex, FarField(point_cart))
-```
+!!! note
+    The total far-field is not defined (since the incident far-field is not defined).

docs/src/ringCurrents.md

@@ -95,9 +95,6 @@ FF = scatteredfield(sp, ex, FarField(point_cart))
 ---
 ## Total Field
 
-!!! warning
-    Not fully implemented yet.
-
 #### API
 
 The general API is employed:

docs/src/sphModes.md

@@ -86,7 +86,6 @@ and
 
 
 
-
 #### API
 
 The general API is employed:
@@ -98,6 +97,10 @@ H  = field(ex, MagneticField(point_cart))
 FF = field(ex, FarField(point_cart))
 ```
 
+!!! note
+    The far-field for a spherical mode makes only sense when travelling outwards.
+
+
 ---
 ## Scattered Field
 
@@ -134,16 +137,14 @@ FF = scatteredfield(sp, ex, FarField(point_cart))
 ---
 ## Total Field
 
-!!! warning
-    Not fully implemented yet.
-
 #### API
 
 The general API is employed:
 ```julia
 E  = field(sp, ex, ElectricField(point_cart))
 
 H  = field(sp, ex, MagneticField(point_cart))
+```
 
-FF = field(sp, ex, FarField(point_cart))
-```
+!!! note
+    The total far-field for a spherical mode excitation is not defined.

docs/src/uniformStatic.md

@@ -70,14 +70,11 @@ E = scatteredfield(sp, ex, ElectricField(point_cart))
 ---
 ## Total Field
 
-!!! warning
-    Not fully implemented yet.
-
 #### API
 
 The following API is employed:
 ```julia
 Φ = field(sp, ex, ScalarPotential(point_cart))
 
 E = field(sp, ex, ElectricField(point_cart))
-```
+```

src/SphericalScattering.jl

@@ -70,4 +70,6 @@ include("sphericalModes/scattered.jl")
 include("UniformField/excitation.jl")
 include("UniformField/incident.jl")
 include("UniformField/scattered.jl")
+
+include("totalFields.jl")
 end

src/UniformField/incident.jl

@@ -28,12 +28,12 @@ The point and the returned field are in Cartesian coordinates.
 function field(excitation::UniformField, point, quantity::ElectricField; parameter::Parameter=Parameter())
 
     a = excitation.amplitude
-
     p = excitation.direction
 
     return a * p
 end
 
+
 """
     field(excitation::UniformField, point, quantity::ScalarPotential; parameter::Parameter=Parameter())
 

src/dipoles/scattered.jl

@@ -2,7 +2,7 @@
 """
     scatteredfield(sphere::PECSphere, excitation::Dipole, quantity::Field; parameter::Parameter=Parameter())
 
-Compute the electric field scattered by a dipole at some position and orientation.
+Compute the field scattered by a PEC sphere excited by a dipole at some position and orientation.
 """
 function scatteredfield(sphere::PECSphere, excitation::Dipole, quantity::Field; parameter::Parameter=Parameter())
 
@@ -278,9 +278,9 @@ end
 
 
 """
-    scatterCoeff(sphere::PECSphere, excitation::ElectricRingCurrent, n::Int, ka)
+    scatterCoeff(sphere::PECSphere, excitation::FitzgeraldDipole, n::Int, ka)
 
-Compute scattering coefficient for electric ring current.
+Compute scattering coefficient for Fitzgerald dipole.
 """
 function scatterCoeff(sphere::PECSphere, excitation::FitzgeraldDipole, n::Int, ka)
 
@@ -294,9 +294,9 @@ end
 
 
 """
-    scatterCoeff(sphere::PECSphere, excitation::MagneticRingCurrent, n::Int, ka)
+    scatterCoeff(sphere::PECSphere, excitation::HertzianDipole, n::Int, ka)
 
-Compute scattering coefficient for magnetic ring current.
+Compute scattering coefficient for Hertzian dipole.
 """
 function scatterCoeff(sphere::PECSphere, excitation::HertzianDipole, n::Int, ka)
 

src/planeWave/incident.jl

@@ -60,3 +60,15 @@ function field(excitation::PlaneWave, point, quantity::MagneticField; parameter:
 
     return a * sqrt(ε / μ) * exp(-im * k * dot(d, point)) * (d × p)
 end
+
+
+
+"""
+    field(excitation::PlaneWave, point, quantity::FarField; parameter::Parameter=Parameter())
+
+Throw error since the far-field of a plane wave is not defined.
+"""
+function field(excitation::PlaneWave, point, quantity::FarField; parameter::Parameter=Parameter())
+
+    error("The far-field of a plane wave is not defined.")
+end

src/totalFields.jl

@@ -0,0 +1,37 @@
+
+"""
+    field(sphere::Sphere, excitation::Excitation, quantity::Field; parameter::Parameter=Parameter())
+
+Compute the total field in the presence of a sphere for a given excitaion.
+"""
+function field(sphere::Sphere, excitation::Excitation, quantity::Field; parameter::Parameter=Parameter())
+
+    # incident and scattered field
+    Finc = field(excitation, quantity; parameter=parameter)
+    Fsca = scatteredfield(sphere, excitation, quantity; parameter=parameter)
+
+    return Finc .+ Fsca
+end
+
+
+
+"""
+    field(sphere::Sphere, excitation::PlaneWave, quantity::Field; parameter::Parameter=Parameter())
+
+Descriptive error for the total far-field in the presence of a sphere for an incident plane wave.
+"""
+function field(sphere::Sphere, excitation::PlaneWave, quantity::FarField; parameter::Parameter=Parameter())
+
+    return error("The total far-field for a plane-wave excitation is not defined")
+end
+
+
+"""
+    field(sphere::Sphere, excitation::SphericalMode, quantity::Field; parameter::Parameter=Parameter())
+
+Descriptive error for the total far-field in the presence of a sphere for an incident spherical mode.
+"""
+function field(sphere::Sphere, excitation::SphericalMode, quantity::FarField; parameter::Parameter=Parameter())
+
+    return error("The total far-field for a spherical mode excitation is not defined")
+end

test/dipoles.jl

@@ -1,11 +1,15 @@
 
-@testset "Hertzian dipole" begin
+f = 1e8
+κ = 2π * f / c   # Wavenumber
 
-    f = 1e8
-    κ = 2π * f / c   # Wavenumber
+# BEAST impedance matrix
+𝑇 = Maxwell3D.singlelayer(; wavenumber=κ)
+T = assemble(𝑇, RT, RT)
 
-    ex = HertzianDipole(; frequency=f, center=SVector(0.0, 0.0, 2.0))
 
+@testset "Hertzian dipole" begin
+
+    ex = HertzianDipole(; frequency=f, center=SVector(0.0, 0.0, 2.0))
 
     @testset "Incident fields" begin
 
@@ -35,10 +39,7 @@
         𝐸 = ex
 
         𝑒 = n × 𝐸 × n
-        𝑇 = Maxwell3D.singlelayer(; wavenumber=κ)
-
         e = assemble(𝑒, RT)
-        T = assemble(𝑇, RT, RT)
 
         u = T \ e
 
@@ -69,8 +70,8 @@ end
 
 @testset "Fitzgerald dipole" begin
 
-    f = 1e8
-    κ = 2π * f / c   # Wavenumber
+    #f = 1e8
+    #κ = 2π * f / c   # Wavenumber
 
 
     ex = FitzgeraldDipole(; frequency=f, center=SVector(0.0, 0.0, 2.0))
@@ -103,10 +104,9 @@ end
         𝐸 = ex
 
         𝑒 = n × 𝐸 × n
-        𝑇 = Maxwell3D.singlelayer(; wavenumber=κ)
-
+        #𝑇 = Maxwell3D.singlelayer(; wavenumber=κ)
         e = assemble(𝑒, RT)
-        T = assemble(𝑇, RT, RT)
+        #T = assemble(𝑇, RT, RT)
 
         u = T \ e
 

test/planeWave.jl

@@ -7,6 +7,7 @@
     ex = planeWave(sp; frequency=f)
 
     @testset "Planewave excitation" begin
+        @test planeWave(; frequency=f) isa PlaneWave{Float64,Float64,Float64}
         @test planeWave(sp; frequency=f) isa PlaneWave{Float64,Float64,Float64}
     end
 
@@ -17,6 +18,8 @@
         @test_nowarn E = field(ex, ElectricField(point_cart))
         @test_nowarn H = field(ex, MagneticField(point_cart))
 
+        @test_throws ErrorException("The far-field of a plane wave is not defined.") field(ex, FarField(point_cart))
+
     end
 
     @testset "Scattered fields" begin
@@ -61,4 +64,30 @@
         @test norm(HF₁) == 0.0
         @test maximum(20 * log10.(abs.(diff_FF))) < -25 # dB
     end
+
+    @testset "Total fields" begin
+
+        # define an observation point
+        point_cart = [SVector(2.0, 2.0, 3.2), SVector(3.1, 4, 2)]
+
+        # compute scattered fields
+        Es = scatteredfield(sp, ex, ElectricField(point_cart))
+        Hs = scatteredfield(sp, ex, MagneticField(point_cart))
+        #FFs = scatteredfield(sp, ex, FarField(point_cart))
+
+        Ei = field(ex, ElectricField(point_cart))
+        Hi = field(ex, MagneticField(point_cart))
+        #FFi = field(ex, FarField(point_cart))
+
+        # total field
+        E = field(sp, ex, ElectricField(point_cart))
+        H = field(sp, ex, MagneticField(point_cart))
+        @test_throws ErrorException("The total far-field for a plane-wave excitation is not defined") field(
+            sp, ex, FarField(point_cart)
+        )
+
+        # is it the sum?
+        @test E[1] == Es[1] .+ Ei[1]
+        @test H[1] == Hs[1] .+ Hi[1]
+    end
 end

test/ringCurrents.jl

@@ -1,9 +1,13 @@
 
-@testset "Electric ring current" begin
+f = 1e8
+κ = 2π * f / c   # Wavenumber
+
+# BEAST impedance matrix
+𝑇 = Maxwell3D.singlelayer(; wavenumber=κ)
+T = assemble(𝑇, RT, RT)
 
-    f = 1e8
-    κ = 2π * f / c   # Wavenumber
 
+@testset "Electric ring current" begin
 
     ex = electricRingCurrent(; frequency=f, center=SVector(0.0, 0.0, 2.0), radius=0.5)
 
@@ -35,10 +39,10 @@
         𝐸 = ex
 
         𝑒 = n × 𝐸 × n
-        𝑇 = Maxwell3D.singlelayer(; wavenumber=κ)
+        #𝑇 = Maxwell3D.singlelayer(; wavenumber=κ)
 
         e = assemble(𝑒, RT)
-        T = assemble(𝑇, RT, RT)
+        #T = assemble(𝑇, RT, RT)
 
         u = T \ e
 
@@ -69,8 +73,8 @@ end
 
 @testset "Magnetic ring current" begin
 
-    f = 1e8
-    κ = 2π * f / c   # Wavenumber
+    #f = 1e8
+    #κ = 2π * f / c   # Wavenumber
 
 
     ex = magneticRingCurrent(; frequency=f, center=SVector(0.0, 0.0, 2.0), radius=0.5)
@@ -103,10 +107,10 @@ end
         𝐸 = ex
 
         𝑒 = n × 𝐸 × n
-        𝑇 = Maxwell3D.singlelayer(; wavenumber=κ)
+        #𝑇 = Maxwell3D.singlelayer(; wavenumber=κ)
 
         e = assemble(𝑒, RT)
-        T = assemble(𝑇, RT, RT)
+        #T = assemble(𝑇, RT, RT)
 
         u = T \ e