From 87887c3a8d3638fe89da1eb7a0cad1f63bae9033 Mon Sep 17 00:00:00 2001 From: Bernd Hofmann Date: Tue, 31 Jan 2023 18:04:17 +0100 Subject: [PATCH] implement general orientations of excitations --- Project.toml | 2 +- docs/src/details.md | 54 +++++++++- docs/src/dipoles.md | 25 +++-- docs/src/manual.md | 5 +- docs/src/planeWave.md | 21 ++-- docs/src/ringCurrents.md | 26 +++-- docs/src/uniformStatic.md | 8 +- src/SphericalScattering.jl | 2 +- src/UniformField/excitation.jl | 4 +- src/coordinateTransforms.jl | 182 ++++++++++++++++++++++++++++----- src/dipoles/excitation.jl | 40 +++++--- src/dipoles/incident.jl | 17 +-- src/dipoles/scattered.jl | 26 ++--- src/planeWave/excitation.jl | 15 ++- src/planeWave/scattered.jl | 7 +- src/ringCurrent/excitation.jl | 84 +++++++++------ src/ringCurrent/incident.jl | 7 +- src/ringCurrent/scattered.jl | 12 +-- test/coordinateTransforms.jl | 4 - test/dipoles.jl | 170 ++++++++++++++++++++++-------- test/planeWave.jl | 102 ++++++++++++------ test/ringCurrents.jl | 176 +++++++++++++++++++++---------- 22 files changed, 704 insertions(+), 285 deletions(-) diff --git a/Project.toml b/Project.toml index 61c9e6b..305198c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SphericalScattering" uuid = "1a9ea918-b599-4f1f-bd9a-d681e8bb5b3e" authors = ["Bernd Hofmann and contributors"] -version = "0.4.0" +version = "0.5.0" [deps] LegendrePolynomials = "3db4a2ba-fc88-11e8-3e01-49c72059a882" diff --git a/docs/src/details.md b/docs/src/details.md index 0af8f69..0c08fe7 100644 --- a/docs/src/details.md +++ b/docs/src/details.md @@ -1,6 +1,12 @@ # Further Details +--- +## Units + +SI units are employed everywhere. + + --- ## [Duality Relations](@id dualityRelations) @@ -17,6 +23,50 @@ are employed in order to compute, e.g., the field of the magnetic current counte --- -## Units +## [Rotations](@id rotationDetails) + +Due to the rotational symmetry the fields for different orientations of the excitations can be computed via rotations. + +#### Plane Wave + +A plane wave excitation with arbitrary direction ``\hat{\bm k}`` and polarization ``\hat{\bm p}`` (forming a valid combination) +can be related to the case ``\hat{\bm k} = \hat{\bm e}_z`` and polarization ``\hat{\bm p} = \hat{\bm e}_x`` by a rotation matrix +```math +\bm R = \begin{bmatrix} \hat{\bm p} & \hat{\bm k} \times \hat{\bm p} & \hat{\bm k} \end{bmatrix} \,. +``` +It constitutes a change from one right-handed orthogonal basis ``(\hat{\bm e}_x^\prime , \hat{\bm e}_y^\prime, \hat{\bm e}_z^\prime )`` to another right-handed orthogonal basis ``(\hat{\bm e}_x = \hat{\bm p}, \hat{\bm e}_y = \hat{\bm k} \times \hat{\bm p}, \hat{\bm e}_z = \hat{\bm k}``. +Hence, ``\bm R`` relates vectors (and points) ``\bm v'`` in the primed coordinate system ``(x',y',z')`` to vectors (and points) ``\bm v`` in the unprimed one ``(x,y,z)`` by +```math +\bm v = \bm R \bm v' \qquad \text{and} \qquad \bm v' = \bm R^\mathrm{T} \bm v +``` +leveraging that ``\bm R^{-1} = \bm R^\mathrm{T}``. +```@raw html +
+``` + +#### Ring-Currents and Dipoles -SI units are employed everywhere. \ No newline at end of file +Ring-currents and dipoles are characterized by only one vector ``\hat{\bm p}`` (defining the orientation). +An arbitrary ``\hat{\bm p}`` can be related to an initial one ``\hat{\bm p}_0`` (e.g., ``\hat{\bm p}_0 = \hat{\bm e}_z^\prime``) by a +right-handed rotation around a rotation axis ``\hat{\bm a} = (a_x, a_y, a_z)`` by an angle ``\alpha``. +The rotation axis is given by +```math +\hat{\bm a} = \cfrac{\hat{\bm p}_0 \times \hat{\bm p}}{|\hat{\bm p}_0 \times \hat{\bm p}|} +``` +and the angle ``\alpha`` is encoded in +```math +\cos(\alpha) = \hat{\bm p}_0 \cdot \hat{\bm p} \qquad \text{and} \qquad \sin(\alpha) = |\hat{\bm p}_0 \times \hat{\bm p}|\,. +``` +The rotation matrix is then found by the [Rodrigues' rotation formula](https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula) as +```math +\bm R = \mathbf{I} + \sin(\alpha)\bm{K} + (1 - \cos(\alpha)) \bm{K}^2 +``` +with +```math +\bm K = \begin{bmatrix} 0 & -a_z & a_y \\ a_z & 0 & -a_x \\ -a_y & a_x & 0 \end{bmatrix} \,. +``` +Vectors (and points) ``\bm v'`` in the initial coordinate system ``(x',y',z')`` (with, e.g., ``\hat{\bm e}_z^\prime = \hat{\bm p}_0``) are then related to vectors (and points) ``\bm v`` in the unprimed one ``(x,y,z)`` by +```math +\bm v = \bm R \bm v' \qquad \text{and} \qquad \bm v' = \bm R^\mathrm{T} \bm v +``` +leveraging that ``\bm R^{-1} = \bm R^\mathrm{T}``. diff --git a/docs/src/dipoles.md b/docs/src/dipoles.md index bdd6824..001ad3e 100644 --- a/docs/src/dipoles.md +++ b/docs/src/dipoles.md @@ -44,8 +44,8 @@ HertzianDipole FitzgeraldDipole ``` -!!! note - The `orientation` is automatically normalized. +!!! tip + The `orientation` vector is automatically normalized to a unit vector during the initialization. --- @@ -64,7 +64,7 @@ where \hat{\bm n} = \cfrac{\bm r - \bm r_0}{|\bm r - \bm r_0|}\,. ``` !!! note - The fields of the Fitzgerald dipole are computed via [duality relations](@ref dualityRelations). + The fields of the Fitzgerald dipole follow from the [duality relations](@ref dualityRelations). #### API @@ -81,15 +81,22 @@ FF = field(ex, FarField(point_cart)) --- ## Scattered Field -The scattered field computation is a generalization of the analysis in [[1, pp. 374ff]](@ref refs). For the magnetic ring current [duality relations](@ref dualityRelations) are employed. +The scattered field computation is a generalization of the analysis in [[1, pp. 374ff]](@ref refs). For the magnetic dipole [duality relations](@ref dualityRelations) are employed. -!!! note - Orientation and location of the dipole are restricted for the scattered field computation! +!!! tip + For the scattered field computation the orientation of the dipole is restricted: the dipole has to be perpendicular to the surface of the sphere. -The dipole has to be oriented such that it is normal to the sphere surface and located outside of the sphere. + This can be achieved, e.g., by defining the orientation first + ```julia + orientation = normalize(SVector(0.0,1.0,1.0)) + ``` + and then setting the position as a multiple of the orientaion: + ```julia + position = 2.0 * orientation + ``` -!!! warning - So far the dipole is assumed to be along the ``z``-axis! This is planned to be generalized. +!!! note + Internal details of the computations: Following [[1, pp. 347ff]](@ref refs) the dipoles are initially assumed to be aligned with the ``z``-axis. Arbitrary positions and orientations (forming a valid pair) are obtained via [rotations](@ref rotationDetails). #### API diff --git a/docs/src/manual.md b/docs/src/manual.md index 967f75a..f74a10a 100644 --- a/docs/src/manual.md +++ b/docs/src/manual.md @@ -11,7 +11,7 @@ The basic building blocks are introduced in the following simple example; more d ```julia using SphericalScattering, StaticArrays -# define excitation: plane wave travelling in negative z-direction with x-polarization +# define excitation: plane wave travelling in positive z-direction with x-polarization ex = planeWave(frequency=10e6) # Hz # define scatterer: PEC sphere @@ -24,7 +24,6 @@ point_cart = [SVector(2.0, 2.0, 3.2)] E = scatteredfield(sp, ex, ElectricField(point_cart)) H = scatteredfield(sp, ex, MagneticField(point_cart)) FF = scatteredfield(sp, ex, FarField(point_cart)) - ``` --- @@ -54,7 +53,7 @@ For all available excitations a simple constructor with keyword arguments and de - [Dipoles](@ref dipolesAPI) - [Ring currents](@ref rcAPI) - [Spherical modes](@ref modesAPI) -- [Uniform Static Field](@ref uniformAPI) +- [Uniform static field](@ref uniformAPI) --- diff --git a/docs/src/planeWave.md b/docs/src/planeWave.md index 30e3c1c..133b4ac 100644 --- a/docs/src/planeWave.md +++ b/docs/src/planeWave.md @@ -18,11 +18,15 @@ --- ## Definition -The plane wave with amplitude ``a``, wave vector ``\bm k``, and polarization ``\hat{\bm p}`` is assumed to have the field +A plane wave with amplitude ``a``, wave vector ``\bm k = k \hat{\bm k}``, and polarization ``\hat{\bm p}`` (vectors with a hat denote unit vectors) is defined by the field ```math -\bm e_\mathrm{PW}(\bm r) = a \hat{\bm p} \, \mathrm{e}^{-\mathrm{j} \bm k \cdot \bm r} \,. +\bm e_\mathrm{PW}(\bm r) = a \hat{\bm p} \, \mathrm{e}^{-\mathrm{j} \bm k \cdot \bm r} \,, ``` - +where the ploarization and wave vector are orthogonal, that is, +```math +\bm k \cdot \hat{\bm p} = 0 +``` +holds. --- ## [API](@id pwAPI) @@ -33,7 +37,10 @@ planeWave ``` !!! note - The `direction` and the `polarization` are each automatically normalized. + The provided `direction` and the `polarization` vectors have to be orthogonal. This is checked during initialization. + +!!! tip + The `direction` and the `polarization` vectors are each automatically normalized to unit vectors during the initialization. --- ## Incident Field @@ -42,7 +49,7 @@ The electric field of the plane wave is as given above. The magnetic field is gi ```math \bm h_\mathrm{PW}(\bm r) = \cfrac{a}{Z_\mathrm{F}} (\hat{\bm k} \times \hat{\bm p}) \mathrm{e}^{-\mathrm{j} \bm k \cdot \bm r} ``` -with ``Z_\mathrm{F} = \sqrt{\mu / \varepsilon}`` +with ``Z_\mathrm{F} = \sqrt{\mu / \varepsilon}``. #### API @@ -61,8 +68,8 @@ H = field(ex, MagneticField(point_cart)) The scattered field computation follows [[1, pp. 347ff]](@ref refs). -!!! warning - So far the plane wave is assumed to travel in positive ``z``-axis direction and to have a polarization along the ``x``-axis! This is planned to be generalized. +!!! note + Internal details of the computations: Following [[1, pp. 347ff]](@ref refs) the plane wave is initially assumed to travel in positive ``z``-axis direction and to have a polarization along the positive ``x``-axis. Arbitrary directions and orientations (forming a valid pair) are obtained via [rotations](@ref rotationDetails). #### API diff --git a/docs/src/ringCurrents.md b/docs/src/ringCurrents.md index c449c49..c79114e 100644 --- a/docs/src/ringCurrents.md +++ b/docs/src/ringCurrents.md @@ -23,7 +23,7 @@ The ring currents are defined by a circular loop of radius ``a`` carrying a unif #### Electric Ring Current -The electric ring current with amplitude ``I``, and orientation ``\hat{\bm e}_z`` at position ``\bm r_0`` is assumed to have the current density +The electric ring current with amplitude ``I``, and orientation ``\hat{\bm p} = \hat{\bm e}_z`` with its center at position ``\bm r_0 = (0,0,z_0)`` is assumed to have the current density ```math \bm{j}_\mathrm{rc} = I \hat{\bm e}_{\varphi'} \delta (r - R_0) \delta (\vartheta - \vartheta_0)\,, ``` @@ -31,11 +31,16 @@ where ``\delta`` denotes the Dirac delta distribution and ``R_0 = \sqrt{z_0^2 + #### Magnetic Ring Current -The magnetic ring current with amplitude ``M``, and orientation ``\hat{\bm e}_z`` at position ``\bm r_0`` is assumed to have the current density +The magnetic ring current with amplitude ``M``, and orientation ``\hat{\bm p} = \hat{\bm e}_z`` at position ``\bm r_0 = (0,0,z_0)`` is assumed to have the current density ```math \bm{m}_\mathrm{rc} = M \hat{\bm e}_{\varphi'} \delta (r - R_0) \delta (\vartheta - \vartheta_0)\,. ``` +!!! note + Arbitrary center positions ``\bm r_0`` and orientations ``\hat{\bm p}`` are computed by [rotations](@ref rotationDetails) and translations of the corresponding fields. + + + --- ## [API](@id rcAPI) @@ -73,13 +78,20 @@ FF = field(ex, FarField(point_cart)) The scattered field computation follows [[1, pp. 368ff]](@ref refs). For the magnetic ring current [duality relations](@ref dualityRelations) are employed. -!!! note - Orientation and location of the ring currents are restricted for the scattered field computation! +!!! tip + For the scattered field computation the orientation of the ring current is restricted: the orientation vector has to be perpendicular to the surface of the sphere. -The ring current has to be oriented such that it is normal to the sphere surface and located outside of the sphere. + This can be achieved, e.g., by defining the orientation first + ```julia + orientation = normalize(SVector(0.0,1.0,1.0)) + ``` + and then setting the center position as a multiple of the orientaion: + ```julia + center = 2.0 * orientation + ``` -!!! warning - So far the ring current is assumed to be along the ``z``-axis! This is planned to be generalized. +!!! note + Internal details of the computations: Following [[1, pp. 368ff]](@ref refs) the orientation vectors of the ring currents are initially assumed to be aligned with the ``z``-axis. Arbitrary centers and orientations (forming a valid pair) are obtained via [rotations](@ref rotationDetails). #### API diff --git a/docs/src/uniformStatic.md b/docs/src/uniformStatic.md index 4ae0233..9fd25b1 100644 --- a/docs/src/uniformStatic.md +++ b/docs/src/uniformStatic.md @@ -32,14 +32,14 @@ The API provides the following constructor with default values: UniformField ``` -!!! note - The direction is automatically normalized. +!!! tip + The `direction` vector is automatically normalized to a unit vector during the initialization. --- ## Incident Field -The electric field of the plane wave is as given above. +The electric field is as given above. #### API @@ -55,8 +55,6 @@ E = field(ex, ElectricField(point_cart)) The scattered field computation follows [[3]](@ref refs). -!!! warning - So far the static electric field is assumed to point in ``x``-direction. This is planned to be generalized. #### API diff --git a/src/SphericalScattering.jl b/src/SphericalScattering.jl index a0d509c..926d535 100644 --- a/src/SphericalScattering.jl +++ b/src/SphericalScattering.jl @@ -47,7 +47,6 @@ export field, scatteredfield # -------- included files -include("coordinateTransforms.jl") include("dataHandling.jl") include("sphere.jl") @@ -72,4 +71,5 @@ include("UniformField/incident.jl") include("UniformField/scattered.jl") include("totalFields.jl") +include("coordinateTransforms.jl") end diff --git a/src/UniformField/excitation.jl b/src/UniformField/excitation.jl index 86ae92c..f90ea0e 100644 --- a/src/UniformField/excitation.jl +++ b/src/UniformField/excitation.jl @@ -16,8 +16,8 @@ end ex = UniformField(; embedding=Medium(ε0, μ0), amplitude=1.0, - direction=SVector{3,Float64}(1.0, 0.0, 0.0) + direction=SVector{3,typeof(amplitude)}(1.0, 0.0, 0.0) ) """ -UniformField(; embedding=Medium(ε0, μ0), amplitude=1.0, direction=SVector{3,Float64}(1.0, 0.0, 0.0)) = +UniformField(; embedding=Medium(ε0, μ0), amplitude=1.0, direction=SVector{3,typeof(amplitude)}(1.0, 0.0, 0.0)) = UniformField(embedding, amplitude, direction) diff --git a/src/coordinateTransforms.jl b/src/coordinateTransforms.jl index 4b2528c..bf17c96 100644 --- a/src/coordinateTransforms.jl +++ b/src/coordinateTransforms.jl @@ -2,7 +2,7 @@ """ translate(points, translation::SVector{3,T}) -Translate the points in negative! direction of the translation vector. +Translate the points in the direction of the translation vector. All inputs are assumed do be in Cartesian coordinates. """ @@ -12,56 +12,163 @@ function translate(points, translation::SVector{3,T}) where {T} points_shifted = similar(points) for (ind, p) in enumerate(points) - points_shifted[ind] = p - translation + points_shifted[ind] = p + translation end return points_shifted end + """ - rotate!(points, rotation::SVector{2,T}) + rotate(excitation::Excitation, vectors_list; inverse=false) -Rotate points around x-axis and y-axis away from the z-axis by angles (in radian) specified in rotation[1] and rotation[2]. +Determine rotation matrix and perform rotation for general excitations. The points are assumed do be in Cartesian coordinates. + +The vectors_list is NOT modified. """ -function rotate!(points, rotation::SVector{2,T}) where {T} +function rotate(excitation::Excitation, vectors_list; inverse=false) - rotation == SVector{2,T}(0.0, 0.0) && return nothing # no rotation required + vectors_list_rot = deepcopy(vectors_list) + rotate!(excitation, vectors_list_rot; inverse=inverse) - # --- perform rotation - # wx = rotation[1] # rotation angle around x-axis (away from z-axis) - # wy = rotation[2] # rotation angle around y-axis (away from z-axis) + return vectors_list_rot +end + + +""" + rotate!(excitation::Excitation, vectors_list; inverse=false) + +Determine rotation matrix and perform rotation for a general excitation. - # cosy = cos(wy) - # siny = sin(wy) - # cosx = cos(wx) - # sinx = sin(wx) +The points are assumed do be in Cartesian coordinates. - # Ry = [ - # cosy 0 siny - # 0 1 0 - # -siny 0 cosy - # ] +The vectors_list IS modified (overwritten). +""" +function rotate!(excitation::Excitation, vectors_list; inverse=false) - # Rx = [ - # 1 0 0 - # 0 cosx sinx - # 0 -sinx cosx - # ] + # --- rotation matrix + R = rotationMatrix(excitation) + isnothing(R) && return nothing - # R = Rx * Ry + # --- inverse? + inverse && (R = R') # for inverse take transpose: inv(R) = R' - # for (ind, p) in enumerate(points) - # points[ind] = R * p - # end + # --- perform rotation + for ind in eachindex(vectors_list) + vectors_list[ind] = R * vectors_list[ind] + end return nothing end +""" + rotationMatrix(excitation::PlaneWave) + +Determine rotation matrix for a plane wave excitation. +""" +function rotationMatrix(excitation::PlaneWave) + + p = excitation.polarization # unit vector + k = excitation.direction # unit vector + T = eltype(p) + + # --- k == ez and p == ex => no ratation + k == SVector{3,T}(0, 0, 1) && p == SVector{3,T}(1, 0, 0) && return nothing + + # --- rotation matrix + R = SMatrix{3,3,T}([p k × p k]) + + return R +end + + +""" + rotationMatrix(excitation::Union{Dipole,RingCurrent}) + +Determine rotation matrix for a dipole or a ring-current excitation via the Rodrigues formula. +""" +function rotationMatrix(excitation::Union{Dipole,RingCurrent}) + + p = excitation.orientation # unit vector + T = eltype(p) + + # --- p == ez => no rotation + p == SVector{3,T}(0, 0, 1) && return nothing + + # --- rotation axis + aux = SVector{3,T}(0, 0, 1) × p + if norm(aux) == T(0) # rotation to -ez + rotAxis = SVector{3,T}(0, 1, 0) # rotate around y-axis + else # all other cases + rotAxis = normalize(aux) + end + + cosϑ = p[3] # inner product of unit vectors: ez ⋅ p = p[3] = cos(ϑ) + sinϑ = sqrt(p[1]^2 + p[2]^2) # cross product of unit vectors: |ez × p| = sin(ϑ) + + # --- auxiliary matrix + K = SMatrix{3,3,T}([ + 0 -rotAxis[3] rotAxis[2] + rotAxis[3] 0 -rotAxis[1] + -rotAxis[2] rotAxis[1] 0 + ]) + + # --- put it together + R = I + sinϑ * K + (1 - cosϑ) * K * K # Rodriguez rotation formula + + return R +end + + +# """ +# rotationMatrix(excitation::UniformField) + +# Determine rotation matrix for a uniform field excitation via the Rodrigues formula. +# """ +# function rotationMatrix(excitation::UniformField) + +# p = excitation.direction # unit vector +# T = eltype(p) + +# # --- p == ex => no rotation +# p == SVector{3,T}(1, 0, 0) && return nothing + +# # --- rotation axis +# aux = SVector{3,T}(1, 0, 0) × p +# if norm(aux) == T(0) # rotation to -ex +# rotAxis = SVector{3,T}(0, 1, 0) # rotate around y-axis +# else # all other cases +# rotAxis = normalize(aux) +# end + +# cosϑ = p[1] # inner product of unit vectors: ex ⋅ p = p[1] = cos(ϑ) +# sinϑ = sqrt(p[2]^2 + p[3]^2) # cross product of unit vectors: |ex × p| = sin(ϑ) + +# # --- auxiliary matrix +# K = SMatrix{3,3,T}([ +# 0 -rotAxis[3] rotAxis[2] +# rotAxis[3] 0 -rotAxis[1] +# -rotAxis[2] rotAxis[1] 0 +# ]) + +# # --- put it together +# R = I + sinϑ * K + (1 - cosϑ) * K * K # Rodriguez rotation formula + +# return R +# end + + +""" + convertSpherical2Cartesian(F_sph, point_sph) +Takes a 3D (3 entry) vector `F_sph` and converts it from its spherical basis to its Cartesian basis representation. + +The location of the vector has to be provided in spherical coordinates by `point_sph` ordered as ``(r, ϑ, φ)``. +""" function convertSpherical2Cartesian(F_sph, point_sph) T = eltype(F_sph) @@ -80,7 +187,14 @@ function convertSpherical2Cartesian(F_sph, point_sph) end -# ----------------------- Convert the computed field from Cartesian to spherical representation +""" +convertSpherical2Cartesian(F_sph, point_sph) + +Takes a 3D (3 entry) vector `F_cart` and converts it from its Cartesian basis to its spherical basis representation. + +The location of the vector has to be provided in spherical coordinates by `point_sph` ordered as ``(r, ϑ, φ)``. +""" + function convertCartesian2Spherical(F_cart, point_sph) T = eltype(F_cart) @@ -99,7 +213,13 @@ function convertCartesian2Spherical(F_cart, point_sph) end +""" + cart2sph(vec) + +Convert a 3D (3 entry) point from Cartesian to spherical coordinates. +""" function cart2sph(vec) + x = vec[1] y = vec[2] z = vec[3] @@ -122,7 +242,13 @@ function cart2sph(vec) end +""" + sph2cart(vec) + +Convert a 3D (e entry) point from spherical to Cartesian coordinates. +""" function sph2cart(vec) + r = vec[1] ϑ = vec[2] ϕ = vec[3] diff --git a/src/dipoles/excitation.jl b/src/dipoles/excitation.jl index 692d2a4..236b07f 100644 --- a/src/dipoles/excitation.jl +++ b/src/dipoles/excitation.jl @@ -5,15 +5,18 @@ struct HertzianDipole{T,R,C} <: Dipole embedding::Medium{C} frequency::R amplitude::T - center::SVector{3,R} + position::SVector{3,R} orientation::SVector{3,R} - # inner constructor: normalize orientation + # inner constructor: normalize orientation and check orientation function HertzianDipole( - embedding::Medium{C}, frequency::R, amplitude::T, center::SVector{3,R}, orientation::SVector{3,R} + embedding::Medium{C}, frequency::R, amplitude::T, position::SVector{3,R}, orientation::SVector{3,R} ) where {T,R,C} + or_normalized = normalize(orientation) - new{T,R,C}(embedding, frequency, amplitude, center, or_normalized) + orientation × position == SVector{3,R}(0, 0, 0) || @info "The dipole is not placed perpendicular to a sphere." + + new{T,R,C}(embedding, frequency, amplitude, position, or_normalized) end end @@ -21,23 +24,28 @@ struct FitzgeraldDipole{T,R,C} <: Dipole embedding::Medium{C} frequency::R amplitude::T - center::SVector{3,R} + position::SVector{3,R} orientation::SVector{3,R} + # inner constructor: normalize orientation and check orientation function FitzgeraldDipole( - embedding::Medium{C}, frequency::R, amplitude::T, center::SVector{3,R}, orientation::SVector{3,R} + embedding::Medium{C}, frequency::R, amplitude::T, position::SVector{3,R}, orientation::SVector{3,R} ) where {T,R,C} + or_normalized = normalize(orientation) - new{T,R,C}(embedding, frequency, amplitude, center, or_normalized) + orientation × position == SVector{3,R}(0, 0, 0) || @info "The dipole is not placed perpendicular to a sphere." + + new{T,R,C}(embedding, frequency, amplitude, position, or_normalized) end end + """ ex = HertzianDipole( embedding = Medium(ε0, μ0), frequency = error("missing argument `frequency`"), amplitude = 1.0, - center = SVector{3,typeof(frequency)}(0.0, 0.0, 0.0), + position = error("missing argument `position`"), orientation = SVector{3,typeof(frequency)}(0.0, 0.0, 1.0), ) """ @@ -45,9 +53,9 @@ HertzianDipole(; embedding = Medium(ε0, μ0), frequency = error("missing argument `frequency`"), amplitude = 1.0, - center = SVector{3,typeof(frequency)}(0.0, 0.0, 0.0), + position = error("missing argument `position`"), orientation = SVector{3,typeof(frequency)}(0.0, 0.0, 1.0), -) = HertzianDipole(embedding, frequency, amplitude, center, orientation) +) = HertzianDipole(embedding, frequency, amplitude, position, orientation) """ @@ -55,7 +63,7 @@ HertzianDipole(; embedding = Medium(ε0, μ0), frequency = error("missing argument `frequency`"), amplitude = 1.0, - center = SVector{3,typeof(frequency)}(0.0, 0.0, 0.0), + position = error("missing argument `position`"), orientation = SVector{3,typeof(frequency)}(0.0, 0.0, 1.0), ) """ @@ -63,9 +71,9 @@ FitzgeraldDipole(; embedding = Medium(ε0, μ0), frequency = error("missing argument `frequency`"), amplitude = 1.0, - center = SVector{3,typeof(frequency)}(0.0, 0.0, 0.0), + position = error("missing argument `position`"), orientation = SVector{3,typeof(frequency)}(0.0, 0.0, 1.0), -) = FitzgeraldDipole(embedding, frequency, amplitude, center, orientation) +) = FitzgeraldDipole(embedding, frequency, amplitude, position, orientation) @@ -92,15 +100,15 @@ function getFieldType(excitation::FitzgeraldDipole, quantity::Field) embedding = Medium(excitation.embedding.μ, excitation.embedding.ε) # exchange μ and ε (duality relations) if typeof(quantity) == ElectricField - exc = FitzgeraldDipole(embedding, excitation.frequency, -excitation.amplitude, excitation.center, excitation.orientation) # change sign (duality realations) + exc = FitzgeraldDipole(embedding, excitation.frequency, -excitation.amplitude, excitation.position, excitation.orientation) # change sign (duality realations) return MagneticField(quantity.locations), exc elseif typeof(quantity) == MagneticField - exc = FitzgeraldDipole(embedding, excitation.frequency, excitation.amplitude, excitation.center, excitation.orientation) + exc = FitzgeraldDipole(embedding, excitation.frequency, excitation.amplitude, excitation.position, excitation.orientation) return ElectricField(quantity.locations), exc else - exc = FitzgeraldDipole(embedding, excitation.frequency, -excitation.amplitude, excitation.center, excitation.orientation) # change sign (duality realations) + exc = FitzgeraldDipole(embedding, excitation.frequency, -excitation.amplitude, excitation.position, excitation.orientation) # change sign (duality realations) return quantity, exc end end diff --git a/src/dipoles/incident.jl b/src/dipoles/incident.jl index f6cb6c6..e155c38 100644 --- a/src/dipoles/incident.jl +++ b/src/dipoles/incident.jl @@ -13,18 +13,11 @@ function field(excitation::Dipole, quantity::Field; parameter::Parameter=Paramet # --- distinguish electric/magnetic current fieldType, exc = getFieldType(excitation, quantity) - # --- translate/rotate coordinates - #points = translate(quantity.locations, excitation.center) - # rotate!(points, -excitation.rotation) - # --- compute field in Cartesian representation for (ind, point) in enumerate(quantity.locations) F[ind] = field(exc, point, fieldType; parameter=parameter) end - # --- rotate resulting field - # rotate!(F, excitation.rotation) - return F end @@ -44,7 +37,7 @@ function field(excitation::Dipole, point, quantity::ElectricField; parameter::Pa ε = excitation.embedding.ε μ = excitation.embedding.μ - r0 = excitation.center + r0 = excitation.position p = excitation.orientation d = point - r0 @@ -70,7 +63,7 @@ function field(excitation::Dipole, point, quantity::MagneticField; parameter::Pa Il = excitation.amplitude k = wavenumber(excitation) - r0 = excitation.center + r0 = excitation.position p = excitation.orientation d = point - r0 @@ -96,7 +89,7 @@ function field(excitation::HertzianDipole, point, quantity::FarField; parameter: ε = excitation.embedding.ε μ = excitation.embedding.μ - r0 = excitation.center + r0 = excitation.position p = excitation.orientation n = point / norm(point) @@ -119,7 +112,7 @@ function field(excitation::FitzgeraldDipole, point, quantity::FarField; paramete Il = excitation.amplitude k = wavenumber(excitation) - r0 = excitation.center + r0 = excitation.position p = excitation.orientation n = point / norm(point) @@ -136,7 +129,7 @@ end # k = wavenumber(excitation) # I0 = excitation.amplitude -# R = norm(excitation.center) +# R = norm(excitation.position) # μ = excitation.embedding.μ # ε = excitation.embedding.ε diff --git a/src/dipoles/scattered.jl b/src/dipoles/scattered.jl index 17ad776..eb41098 100644 --- a/src/dipoles/scattered.jl +++ b/src/dipoles/scattered.jl @@ -6,18 +6,18 @@ Compute the field scattered by a PEC sphere excited by a dipole at some position """ function scatteredfield(sphere::PECSphere, excitation::Dipole, quantity::Field; parameter::Parameter=Parameter()) - sphere.embedding == excitation.embedding || error("Excitation and sphere are not in the same medium.") # verify excitation and sphere are in the same medium - T = typeof(excitation.frequency) + sphere.embedding == excitation.embedding || error("Excitation and sphere are not in the same medium.") # verify excitation and sphere are in the same medium + excitation.orientation × excitation.position == SVector{3,T}(0, 0, 0) || error("The dipole is not perpendicular to the sphere.") + F = zeros(SVector{3,Complex{T}}, size(quantity.locations)) # --- distinguish electric/magnetic current fieldType, exc = getFieldType(excitation, quantity) - # --- translate/rotate coordinates - points = quantity.locations # translate(quantity.locations, -excitation.center) - # rotate!(points, -excitation.rotation) + # --- rotate coordinates + points = rotate(excitation, quantity.locations; inverse=true) # --- compute field in Cartesian representation for (ind, point) in enumerate(points) @@ -25,7 +25,7 @@ function scatteredfield(sphere::PECSphere, excitation::Dipole, quantity::Field; end # --- rotate resulting field - # rotate!(F, excitation.rotation) + rotate!(excitation, F; inverse=false) return F end @@ -46,7 +46,7 @@ function scatteredfield(sphere::PECSphere, excitation::Dipole, point, quantity:: k = wavenumber(excitation) T = typeof(k) Il = excitation.amplitude - z0 = norm(excitation.center) # distance Dipole-origin + z0 = norm(excitation.position) # distance Dipole-origin μ = excitation.embedding.μ ε = excitation.embedding.ε @@ -96,8 +96,8 @@ function scatteredfield(sphere::PECSphere, excitation::Dipole, point, quantity:: ZF = sqrt(μ / ε) - Er *= -im * Il / 8 / k / z0 * sqrt(r / z0) / r / r * ZF - Eϑ *= im * Il / 8 / k / z0 / z0 * sinϑ / r * ZF + Er *= +im * Il / 8 / k / z0 * sqrt(r / z0) / r / r * ZF + Eϑ *= -im * Il / 8 / k / z0 / z0 * sinϑ / r * ZF return convertSpherical2Cartesian(SVector{3,Complex{T}}(Er, Eϑ, 0.0), point_sph) end @@ -118,7 +118,7 @@ function scatteredfield(sphere::PECSphere, excitation::Dipole, point, quantity:: k = wavenumber(excitation) T = typeof(k) Il = excitation.amplitude - z0 = norm(excitation.center) # distance Dipole-origin + z0 = norm(excitation.position) # distance Dipole-origin eps = parameter.relativeAccuracy @@ -177,7 +177,7 @@ function scatteredfield(sphere::PECSphere, excitation::HertzianDipole, point, qu k = wavenumber(excitation) T = typeof(k) Il = excitation.amplitude - z0 = norm(excitation.center) # distance dipole-origin + z0 = norm(excitation.position) # distance dipole-origin μ = excitation.embedding.μ ε = excitation.embedding.ε @@ -212,7 +212,7 @@ function scatteredfield(sphere::PECSphere, excitation::HertzianDipole, point, qu # print("did not converge: n=$n") # if Hankel function throws overflow error -> result still agrees with small argument approximation (verified) end - Eϑ *= Il / k * sinϑ / z0 / 4 * sqrt(k / z0 / 2 / π) * sqrt(μ / ε) + Eϑ *= -Il / k * sinϑ / z0 / 4 * sqrt(k / z0 / 2 / π) * sqrt(μ / ε) return SVector{3,Complex{T}}( Eϑ * cos(point_sph[2]) * cos(point_sph[3]), Eϑ * cos(point_sph[2]) * sin(point_sph[3]), -Eϑ * sin(point_sph[2]) @@ -235,7 +235,7 @@ function scatteredfield(sphere::PECSphere, excitation::FitzgeraldDipole, point, k = wavenumber(excitation) T = typeof(k) Ul = excitation.amplitude - z0 = norm(excitation.center) # distance dipole-origin + z0 = norm(excitation.position) # distance dipole-origin μ = excitation.embedding.μ ε = excitation.embedding.ε diff --git a/src/planeWave/excitation.jl b/src/planeWave/excitation.jl index 0e0c493..5af806f 100644 --- a/src/planeWave/excitation.jl +++ b/src/planeWave/excitation.jl @@ -6,12 +6,17 @@ struct PlaneWave{T,R,C} <: Excitation direction::SVector{3,R} polarization::SVector{3,R} - # inner constructor: normalize direction and polarization + # inner constructor: normalize direction and polarization and check their orthogonality function PlaneWave( embedding::Medium{C}, frequency::R, amplitude::T, direction::SVector{3,R}, polarization::SVector{3,R} ) where {T,R,C} + + dot(direction, polarization) == 0.0 || + error("No porper definition of a plane wave: direction and polarization vector are not perpendicular!") + dir_normalized = normalize(direction) pol_normalized = normalize(polarization) + new{T,R,C}(embedding, frequency, amplitude, dir_normalized, pol_normalized) end end @@ -22,7 +27,7 @@ end embedding = Medium(ε0, μ0), frequency = error("missing argument `frequency`"), amplitude = 1.0, - direction = SVector{3,typeof(frequency)}(0.0, 0.0, -1.0), + direction = SVector{3,typeof(frequency)}(0.0, 0.0, 1.0), polarization = SVector{3,typeof(frequency)}(1.0, 0.0, 0.0), ) """ @@ -30,7 +35,7 @@ planeWave(; embedding = Medium(ε0, μ0), frequency = error("missing argument `frequency`"), amplitude = 1.0, - direction = SVector{3,typeof(frequency)}(0.0, 0.0, -1.0), + direction = SVector{3,typeof(frequency)}(0.0, 0.0, 1.0), polarization = SVector{3,typeof(frequency)}(1.0, 0.0, 0.0), ) = PlaneWave(embedding, frequency, amplitude, direction, polarization) @@ -40,7 +45,7 @@ planeWave(; sp::Sphere; frequency = error("missing argument `frequency`"), amplitude = 1.0, - direction = SVector{3,typeof(frequency)}(0.0, 0.0, -1.0), + direction = SVector{3,typeof(frequency)}(0.0, 0.0, 1.0), polarization = SVector{3,typeof(frequency)}(1.0, 0.0, 0.0), ) """ @@ -48,6 +53,6 @@ planeWave( sp::Sphere; frequency = error("missing argument `frequency`"), amplitude = 1.0, - direction = SVector{3,typeof(frequency)}(0.0, 0.0, -1.0), + direction = SVector{3,typeof(frequency)}(0.0, 0.0, 1.0), polarization = SVector{3,typeof(frequency)}(1.0, 0.0, 0.0), ) = PlaneWave(sp.embedding, frequency, amplitude, direction, polarization) diff --git a/src/planeWave/scattered.jl b/src/planeWave/scattered.jl index eaddda7..3e6f636 100644 --- a/src/planeWave/scattered.jl +++ b/src/planeWave/scattered.jl @@ -9,19 +9,18 @@ function scatteredfield(sphere::Sphere, excitation::PlaneWave, quantity::Field; sphere.embedding == excitation.embedding || error("Excitation and sphere are not in the same medium.") # verify excitation and sphere are in the same medium T = typeof(excitation.frequency) - F = zeros(SVector{3,Complex{T}}, size(quantity.locations)) # --- rotate coordinates - # rotate!(points, -excitation.rotation) + points = rotate(excitation, quantity.locations; inverse=true) # --- compute field in Cartesian representation - for (ind, point) in enumerate(quantity.locations) + for (ind, point) in enumerate(points) F[ind] = scatteredfield(sphere, excitation, point, quantity; parameter=parameter) end # --- rotate resulting field - # rotate!(F, excitation.rotation) + rotate!(excitation, F; inverse=false) return F end diff --git a/src/ringCurrent/excitation.jl b/src/ringCurrent/excitation.jl index e7c01df..70ee237 100644 --- a/src/ringCurrent/excitation.jl +++ b/src/ringCurrent/excitation.jl @@ -7,7 +7,18 @@ struct ElectricRingCurrent{T,R,C} <: RingCurrent amplitude::T radius::R center::SVector{3,R} - rotation::SVector{2,R} + orientation::SVector{3,R} + + # inner constructor: normalize orientation and check orientation + function ElectricRingCurrent( + embedding::Medium{C}, frequency::R, amplitude::T, radius::R, center::SVector{3,R}, orientation::SVector{3,R} + ) where {T,R,C} + + or_normalized = normalize(orientation) + orientation × center == SVector{3,R}(0, 0, 0) || @info "The ring current is not placed perpendicular to a sphere." + + new{T,R,C}(embedding, frequency, amplitude, radius, center, or_normalized) + end end struct MagneticRingCurrent{T,R,C} <: RingCurrent @@ -16,48 +27,59 @@ struct MagneticRingCurrent{T,R,C} <: RingCurrent amplitude::T radius::R center::SVector{3,R} - rotation::SVector{2,R} + orientation::SVector{3,R} + + # inner constructor: normalize orientation and check orientation + function MagneticRingCurrent( + embedding::Medium{C}, frequency::R, amplitude::T, radius::R, center::SVector{3,R}, orientation::SVector{3,R} + ) where {T,R,C} + + or_normalized = normalize(orientation) + orientation × center == SVector{3,R}(0, 0, 0) || @info "The ring current is not placed perpendicular to a sphere." + + new{T,R,C}(embedding, frequency, amplitude, radius, center, or_normalized) + end end """ ex = electricRingCurrent(; - embedding = Medium(ε0, μ0), - frequency = error("missing argument `frequency`"), - amplitude = 1.0, - radius = error("missing argument `radius`"), - center = SVector{3,typeof(frequency)}(0.0, 0.0, 0.0), - rotation = SVector{2,typeof(frequency)}(0.0, 0.0), + embedding = Medium(ε0, μ0), + frequency = error("missing argument `frequency`"), + amplitude = 1.0, + radius = error("missing argument `radius`"), + center = error("missing argument `center`"), + orientation = SVector{3,typeof(frequency)}(0.0, 0.0, 1.0), ) """ electricRingCurrent(; - embedding = Medium(ε0, μ0), - frequency = error("missing argument `frequency`"), - amplitude = 1.0, - radius = error("missing argument `radius`"), - center = SVector{3,typeof(frequency)}(0.0, 0.0, 0.0), - rotation = SVector{2,typeof(frequency)}(0.0, 0.0), -) = ElectricRingCurrent(embedding, frequency, amplitude, radius, center, rotation) + embedding = Medium(ε0, μ0), + frequency = error("missing argument `frequency`"), + amplitude = 1.0, + radius = error("missing argument `radius`"), + center = error("missing argument `center`"), + orientation = SVector{3,typeof(frequency)}(0.0, 0.0, 1.0), +) = ElectricRingCurrent(embedding, frequency, amplitude, radius, center, orientation) """ ex = magneticRingCurrent(; - embedding = Medium(ε0, μ0), - frequency = error("missing argument `frequency`"), - amplitude = 1.0, - radius = error("missing argument `radius`"), - center = SVector{3,typeof(frequency)}(0.0, 0.0, 0.0), - rotation = SVector{2,typeof(frequency)}(0.0, 0.0), + embedding = Medium(ε0, μ0), + frequency = error("missing argument `frequency`"), + amplitude = 1.0, + radius = error("missing argument `radius`"), + center = error("missing argument `center`"), + orientation = SVector{3,typeof(frequency)}(0.0, 0.0, 1.0), ) """ magneticRingCurrent(; - embedding = Medium(ε0, μ0), - frequency = error("missing argument `frequency`"), - amplitude = 1.0, - radius = error("missing argument `radius`"), - center = SVector{3,typeof(frequency)}(0.0, 0.0, 0.0), - rotation = SVector{2,typeof(frequency)}(0.0, 0.0), -) = MagneticRingCurrent(embedding, frequency, amplitude, radius, center, rotation) + embedding = Medium(ε0, μ0), + frequency = error("missing argument `frequency`"), + amplitude = 1.0, + radius = error("missing argument `radius`"), + center = error("missing argument `center`"), + orientation = SVector{3,typeof(frequency)}(0.0, 0.0, 1.0), +) = MagneticRingCurrent(embedding, frequency, amplitude, radius, center, orientation) @@ -85,19 +107,19 @@ function getFieldType(excitation::MagneticRingCurrent, quantity::Field) if typeof(quantity) == ElectricField exc = MagneticRingCurrent( - embedding, excitation.frequency, -excitation.amplitude, excitation.radius, excitation.center, excitation.rotation + embedding, excitation.frequency, -excitation.amplitude, excitation.radius, excitation.center, excitation.orientation ) # change sign (duality realations) return MagneticField(quantity.locations), exc elseif typeof(quantity) == MagneticField exc = MagneticRingCurrent( - embedding, excitation.frequency, excitation.amplitude, excitation.radius, excitation.center, excitation.rotation + embedding, excitation.frequency, excitation.amplitude, excitation.radius, excitation.center, excitation.orientation ) return ElectricField(quantity.locations), exc else exc = MagneticRingCurrent( - embedding, excitation.frequency, -excitation.amplitude, excitation.radius, excitation.center, excitation.rotation + embedding, excitation.frequency, -excitation.amplitude, excitation.radius, excitation.center, excitation.orientation ) # change sign (duality realations) return quantity, exc end diff --git a/src/ringCurrent/incident.jl b/src/ringCurrent/incident.jl index f603361..1078ef3 100644 --- a/src/ringCurrent/incident.jl +++ b/src/ringCurrent/incident.jl @@ -1,5 +1,4 @@ - """ field(excitation::RingCurrent, quantity::ElectricField; parameter::Parameter=Parameter()) @@ -14,8 +13,8 @@ function field(excitation::RingCurrent, quantity::Field; parameter::Parameter=Pa fieldType, exc = getFieldType(excitation, quantity) # --- translate/rotate coordinates - points = translate(quantity.locations, excitation.center) - rotate!(points, -excitation.rotation) + points = translate(quantity.locations, -excitation.center) + rotate!(excitation, points; inverse=true) # --- compute field in Cartesian representation for (ind, point) in enumerate(points) @@ -23,7 +22,7 @@ function field(excitation::RingCurrent, quantity::Field; parameter::Parameter=Pa end # --- rotate resulting field - rotate!(F, excitation.rotation) + rotate!(excitation, F; inverse=false) return F end diff --git a/src/ringCurrent/scattered.jl b/src/ringCurrent/scattered.jl index 0375da7..29061e5 100644 --- a/src/ringCurrent/scattered.jl +++ b/src/ringCurrent/scattered.jl @@ -1,5 +1,4 @@ - """ field(excitation::ElectricRingCurrent, quantity::Field; parameter::Parameter=Parameter()) @@ -7,18 +6,19 @@ Compute the electric field radiated by an electric ring current at some position """ function scatteredfield(sphere::PECSphere, excitation::RingCurrent, quantity::Field; parameter::Parameter=Parameter()) - sphere.embedding == excitation.embedding || error("Excitation and sphere are not in the same medium.") # verify excitation and sphere are in the same medium - T = typeof(excitation.frequency) + sphere.embedding == excitation.embedding || error("Excitation and sphere are not in the same medium.") # verify excitation and sphere are in the same medium + excitation.orientation × excitation.center == SVector{3,T}(0, 0, 0) || + error("The ring current is not perpendicular to the sphere.") + F = zeros(SVector{3,Complex{T}}, size(quantity.locations)) # --- distinguish electric/magnetic current fieldType, exc = getFieldType(excitation, quantity) # --- translate/rotate coordinates - points = quantity.locations #translate(quantity.locations, -excitation.center) - #rotate!(points, -excitation.rotation) + points = rotate(excitation, quantity.locations; inverse=true) # --- compute field in Cartesian representation for (ind, point) in enumerate(points) @@ -26,7 +26,7 @@ function scatteredfield(sphere::PECSphere, excitation::RingCurrent, quantity::Fi end # --- rotate resulting field - #rotate!(F, excitation.rotation) + rotate!(excitation, F; inverse=false) return F end diff --git a/test/coordinateTransforms.jl b/test/coordinateTransforms.jl index 0e56c20..2ec12fd 100644 --- a/test/coordinateTransforms.jl +++ b/test/coordinateTransforms.jl @@ -22,8 +22,4 @@ @test F_sph[1] ≈ +1.0 @test F_sph[2] ≈ -1.0 @test F_sph[3] ≈ -1.0 - - # ----- rotate - point = SVector(2.0, 2.2, 3.1) - @test_nowarn SphericalScattering.rotate!(point, SVector(1.0, 2.3)) end diff --git a/test/dipoles.jl b/test/dipoles.jl index 5658a01..e65fe1c 100644 --- a/test/dipoles.jl +++ b/test/dipoles.jl @@ -9,7 +9,7 @@ T = assemble(𝑇, RT, RT) @testset "Hertzian dipole" begin - ex = HertzianDipole(; frequency=f, center=SVector(0.0, 0.0, 2.0)) + ex = HertzianDipole(; frequency=f, position=SVector(0.0, 0.0, 2.0)) @testset "Incident fields" begin @@ -35,35 +35,74 @@ T = assemble(𝑇, RT, RT) @testset "Scattered fields" begin - # ----- BEAST solution - 𝐸 = ex + @testset "Standard orientation" begin - 𝑒 = n × 𝐸 × n - e = assemble(𝑒, RT) + ex = HertzianDipole(; frequency=f, position=SVector(0.0, 0.0, 2.0)) - u = T \ e + # ----- BEAST solution + 𝐸 = ex - EF_MoM = +potential(MWSingleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) - HF_MoM = -potential(BEAST.MWDoubleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) - FF_MoM = -im * f / (2 * c) * potential(MWFarField3D(; gamma=𝑇.gamma), points_cartFF, u, RT) + 𝑒 = n × 𝐸 × n + e = -assemble(𝑒, RT) + u = T \ e - # ----- this package - sp = PECSphere(; radius=spRadius) + EF_MoM = potential(MWSingleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) + HF_MoM = 1 / (c * 𝜇) * potential(BEAST.MWDoubleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) + FF_MoM = -im * f / (2 * c) * potential(MWFarField3D(; gamma=𝑇.gamma), points_cartFF, u, RT) - EF = scatteredfield(sp, ex, ElectricField(points_cartNF)) - HF = scatteredfield(sp, ex, MagneticField(points_cartNF)) * c * 𝜇 - FF = scatteredfield(sp, ex, FarField(points_cartFF)) + # ----- this package + sp = PECSphere(; radius=spRadius) + EF = scatteredfield(sp, ex, ElectricField(points_cartNF)) + HF = scatteredfield(sp, ex, MagneticField(points_cartNF)) + FF = scatteredfield(sp, ex, FarField(points_cartFF)) - # ----- compare - diff_EF = norm.(EF - EF_MoM) ./ maximum(norm.(EF)) # worst case error - diff_HF = norm.(HF - HF_MoM) ./ maximum(norm.(HF)) # worst case error - diff_FF = norm.(FF - FF_MoM) ./ maximum(norm.(FF)) # worst case error + # ----- compare + diff_EF = norm.(EF - EF_MoM) ./ maximum(norm.(EF)) # worst case error + diff_HF = norm.(HF - HF_MoM) ./ maximum(norm.(HF)) # worst case error + diff_FF = norm.(FF - FF_MoM) ./ maximum(norm.(FF)) # worst case error - @test maximum(20 * log10.(abs.(diff_EF))) < -25 # dB - @test maximum(20 * log10.(abs.(diff_HF))) < -25 # dB - @test maximum(20 * log10.(abs.(diff_FF))) < -25 # dB + @test maximum(20 * log10.(abs.(diff_EF))) < -25 # dB + @test maximum(20 * log10.(abs.(diff_HF))) < -25 # dB + @test maximum(20 * log10.(abs.(diff_FF))) < -25 # dB + + end + + @testset "General orientation" begin + + ori = normalize(SVector(0.0, -1.0, -1.0)) + ex = HertzianDipole(; frequency=f, orientation=ori, position=ori * 2) + + # ----- BEAST solution + 𝐸 = ex + + 𝑒 = n × 𝐸 × n + e = -assemble(𝑒, RT) + + u = T \ e + + EF_MoM = potential(MWSingleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) + HF_MoM = 1 / (c * 𝜇) * potential(BEAST.MWDoubleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) + FF_MoM = -im * f / (2 * c) * potential(MWFarField3D(; gamma=𝑇.gamma), points_cartFF, u, RT) + + # ----- this package + sp = PECSphere(; radius=spRadius) + + EF = scatteredfield(sp, ex, ElectricField(points_cartNF)) + HF = scatteredfield(sp, ex, MagneticField(points_cartNF)) + FF = scatteredfield(sp, ex, FarField(points_cartFF)) + + # ----- compare + diff_EF = norm.(EF - EF_MoM) ./ maximum(norm.(EF)) # worst case error + diff_HF = norm.(HF - HF_MoM) ./ maximum(norm.(HF)) # worst case error + diff_FF = norm.(FF - FF_MoM) ./ maximum(norm.(FF)) # worst case error + + @test maximum(20 * log10.(abs.(diff_EF))) < -25 # dB + @test maximum(20 * log10.(abs.(diff_HF))) < -25 # dB + @test maximum(20 * log10.(abs.(diff_FF))) < -25 # dB + + end end end @@ -74,7 +113,7 @@ end #κ = 2π * f / c # Wavenumber - ex = FitzgeraldDipole(; frequency=f, center=SVector(0.0, 0.0, 2.0)) + ex = FitzgeraldDipole(; frequency=f, position=SVector(0.0, 0.0, 2.0)) @testset "Incident fields" begin @@ -100,36 +139,79 @@ end @testset "Scattered fields" begin - # ----- BEAST solution - 𝐸 = ex + @testset "Standard orientation" begin + + ex = FitzgeraldDipole(; frequency=f, position=SVector(0.0, 0.0, 2.0)) + + # ----- BEAST solution + 𝐸 = ex + + 𝑒 = n × 𝐸 × n + #𝑇 = Maxwell3D.singlelayer(; wavenumber=κ) + e = assemble(𝑒, RT) + #T = assemble(𝑇, RT, RT) + + u = T \ e + + EF_MoM = -potential(MWSingleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) + HF_MoM = -potential(BEAST.MWDoubleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) + FF_MoM = +im * f / (2 * c) * potential(MWFarField3D(; gamma=𝑇.gamma), points_cartFF, u, RT) + + + # ----- this package + sp = PECSphere(; radius=spRadius) + + EF = scatteredfield(sp, ex, ElectricField(points_cartNF)) + HF = scatteredfield(sp, ex, MagneticField(points_cartNF)) * c * 𝜇 + FF = scatteredfield(sp, ex, FarField(points_cartFF)) + + + # ----- compare + diff_EF = norm.(EF - EF_MoM) ./ maximum(norm.(EF)) # worst case error + diff_HF = norm.(HF - HF_MoM) ./ maximum(norm.(HF)) # worst case error + diff_FF = norm.(FF - FF_MoM) ./ maximum(norm.(FF)) # worst case error + + @test maximum(20 * log10.(abs.(diff_EF))) < -24 # dB + @test maximum(20 * log10.(abs.(diff_HF))) < -24 # dB + @test maximum(20 * log10.(abs.(diff_FF))) < -24 # dB + end + + @testset "General orientation" begin + + ori = normalize(SVector(0.0, -1.0, -1.0)) + ex = FitzgeraldDipole(; frequency=f, orientation=ori, position=ori * 2) + + # ----- BEAST solution + 𝐸 = ex - 𝑒 = n × 𝐸 × n - #𝑇 = Maxwell3D.singlelayer(; wavenumber=κ) - e = assemble(𝑒, RT) - #T = assemble(𝑇, RT, RT) + 𝑒 = n × 𝐸 × n + #𝑇 = Maxwell3D.singlelayer(; wavenumber=κ) + e = assemble(𝑒, RT) + #T = assemble(𝑇, RT, RT) - u = T \ e + u = T \ e - EF_MoM = -potential(MWSingleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) - HF_MoM = +potential(BEAST.MWDoubleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) - FF_MoM = +im * f / (2 * c) * potential(MWFarField3D(; gamma=𝑇.gamma), points_cartFF, u, RT) + EF_MoM = -potential(MWSingleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) + HF_MoM = -potential(BEAST.MWDoubleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) + FF_MoM = +im * f / (2 * c) * potential(MWFarField3D(; gamma=𝑇.gamma), points_cartFF, u, RT) - # ----- this package - sp = PECSphere(; radius=spRadius) + # ----- this package + sp = PECSphere(; radius=spRadius) - EF = scatteredfield(sp, ex, ElectricField(points_cartNF)) - HF = scatteredfield(sp, ex, MagneticField(points_cartNF)) * c * 𝜇 - FF = scatteredfield(sp, ex, FarField(points_cartFF)) + EF = scatteredfield(sp, ex, ElectricField(points_cartNF)) + HF = scatteredfield(sp, ex, MagneticField(points_cartNF)) * c * 𝜇 + FF = scatteredfield(sp, ex, FarField(points_cartFF)) - # ----- compare - diff_EF = norm.(EF - EF_MoM) ./ maximum(norm.(EF)) # worst case error - diff_HF = norm.(HF - HF_MoM) ./ maximum(norm.(HF)) # worst case error - diff_FF = norm.(FF - FF_MoM) ./ maximum(norm.(FF)) # worst case error + # ----- compare + diff_EF = norm.(EF - EF_MoM) ./ maximum(norm.(EF)) # worst case error + diff_HF = norm.(HF - HF_MoM) ./ maximum(norm.(HF)) # worst case error + diff_FF = norm.(FF - FF_MoM) ./ maximum(norm.(FF)) # worst case error - @test maximum(20 * log10.(abs.(diff_EF))) < -24 # dB - @test maximum(20 * log10.(abs.(diff_HF))) < -24 # dB - @test maximum(20 * log10.(abs.(diff_FF))) < -24 # dB + @test maximum(20 * log10.(abs.(diff_EF))) < -24 # dB + @test maximum(20 * log10.(abs.(diff_HF))) < -24 # dB + @test maximum(20 * log10.(abs.(diff_FF))) < -24 # dB + end end end diff --git a/test/planeWave.jl b/test/planeWave.jl index af4d336..4111289 100644 --- a/test/planeWave.jl +++ b/test/planeWave.jl @@ -2,6 +2,7 @@ @testset "PEC" begin f = 1e8 + κ = 2π * f / c # Wavenumber sp = PECSphere(; radius=spRadius, embedding=Medium(𝜀, 𝜇)) ex = planeWave(sp; frequency=f) @@ -24,47 +25,90 @@ @testset "Scattered fields" begin - # ----- BEAST solution - κ = 2π * f / c # Wavenumber + @testset "Standard orientation" begin - 𝐸 = Maxwell3D.planewave(; direction=ẑ, polarization=x̂, wavenumber=κ) + # ----- BEAST solution + 𝐸 = Maxwell3D.planewave(; direction=ẑ, polarization=x̂, wavenumber=κ) - 𝑒 = n × 𝐸 × n - 𝑇 = Maxwell3D.singlelayer(; wavenumber=κ, alpha=-im * 𝜇 * (2π * f), beta=1 / (-im * 𝜀 * (2π * f))) + 𝑒 = n × 𝐸 × n + 𝑇 = Maxwell3D.singlelayer(; wavenumber=κ, alpha=-im * 𝜇 * (2π * f), beta=1 / (-im * 𝜀 * (2π * f))) - e = -assemble(𝑒, RT) - T = assemble(𝑇, RT, RT) + e = -assemble(𝑒, RT) + T = assemble(𝑇, RT, RT) - u = T \ e + u = T \ e - EF_MoM₂ = potential(MWSingleLayerField3D(𝑇), points_cartNF, u, RT) - HF_MoM₂ = potential(BEAST.MWDoubleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) - FF_MoM = -im * f / (2 * c) * potential(MWFarField3D(𝑇), points_cartFF, u, RT) + EF_MoM₂ = potential(MWSingleLayerField3D(𝑇), points_cartNF, u, RT) + HF_MoM₂ = potential(BEAST.MWDoubleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) + FF_MoM = -im * f / (2 * c) * potential(MWFarField3D(𝑇), points_cartFF, u, RT) - # ----- this package - EF₂ = scatteredfield(sp, ex, ElectricField(points_cartNF)) - EF₁ = scatteredfield(sp, ex, ElectricField(points_cartNF_inside)) - HF₂ = scatteredfield(sp, ex, MagneticField(points_cartNF)) - HF₁ = scatteredfield(sp, ex, MagneticField(points_cartNF_inside)) - FF = scatteredfield(sp, ex, FarField(points_cartFF)) + # ----- this package + ex = planeWave(sp; frequency=f) + EF₂ = scatteredfield(sp, ex, ElectricField(points_cartNF)) + EF₁ = scatteredfield(sp, ex, ElectricField(points_cartNF_inside)) + HF₂ = scatteredfield(sp, ex, MagneticField(points_cartNF)) + HF₁ = scatteredfield(sp, ex, MagneticField(points_cartNF_inside)) + FF = scatteredfield(sp, ex, FarField(points_cartFF)) - # ----- compare - diff_EF₂ = norm.(EF₂ - EF_MoM₂) ./ maximum(norm.(EF₂)) # worst case error - diff_HF₂ = norm.(HF₂ - HF_MoM₂) ./ maximum(norm.(HF₂)) # worst case error - diff_FF = norm.(FF - FF_MoM) ./ maximum(norm.(FF)) # worst case error - #@show maximum(20 * log10.(abs.(diff_EF))) - #@show maximum(20 * log10.(abs.(diff_HF))) - #@show maximum(20 * log10.(abs.(diff_FF))) + # ----- compare + diff_EF₂ = norm.(EF₂ - EF_MoM₂) ./ maximum(norm.(EF₂)) # worst case error + diff_HF₂ = norm.(HF₂ - HF_MoM₂) ./ maximum(norm.(HF₂)) # worst case error + diff_FF = norm.(FF - FF_MoM) ./ maximum(norm.(FF)) # worst case error - @test maximum(20 * log10.(abs.(diff_EF₂))) < -25 # dB - @test norm(EF₁) == 0.0 - @test maximum(20 * log10.(abs.(diff_HF₂))) < -25 # dB - @test norm(HF₁) == 0.0 - @test maximum(20 * log10.(abs.(diff_FF))) < -25 # dB + @test maximum(20 * log10.(abs.(diff_EF₂))) < -25 # dB + @test norm(EF₁) == 0.0 + @test maximum(20 * log10.(abs.(diff_HF₂))) < -25 # dB + @test norm(HF₁) == 0.0 + @test maximum(20 * log10.(abs.(diff_FF))) < -25 # dB + end + + @testset "General orientation" begin + + # ----- BEAST solution + dir = normalize(SVector(0.0, 1.0, 1.0)) # normalization for BEAST + pol = normalize(SVector(-1.0, 0.0, 0.0)) + + + 𝐸 = Maxwell3D.planewave(; direction=dir, polarization=pol, wavenumber=κ) + + 𝑒 = n × 𝐸 × n + 𝑇 = Maxwell3D.singlelayer(; wavenumber=κ, alpha=-im * 𝜇 * (2π * f), beta=1 / (-im * 𝜀 * (2π * f))) + + e = -assemble(𝑒, RT) + T = assemble(𝑇, RT, RT) + + u = T \ e + + EF_MoM₂ = potential(MWSingleLayerField3D(𝑇), points_cartNF, u, RT) + HF_MoM₂ = potential(BEAST.MWDoubleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) + FF_MoM = -im * f / (2 * c) * potential(MWFarField3D(𝑇), points_cartFF, u, RT) + + # ----- this package + ex = planeWave(sp; frequency=f, direction=dir, polarization=pol) + + EF₂ = scatteredfield(sp, ex, ElectricField(points_cartNF)) + EF₁ = scatteredfield(sp, ex, ElectricField(points_cartNF_inside)) + HF₂ = scatteredfield(sp, ex, MagneticField(points_cartNF)) + HF₁ = scatteredfield(sp, ex, MagneticField(points_cartNF_inside)) + FF = scatteredfield(sp, ex, FarField(points_cartFF)) + + + # ----- compare + diff_EF₂ = norm.(EF₂ - EF_MoM₂) ./ maximum(norm.(EF₂)) # worst case error + diff_HF₂ = norm.(HF₂ - HF_MoM₂) ./ maximum(norm.(HF₂)) # worst case error + diff_FF = norm.(FF - FF_MoM) ./ maximum(norm.(FF)) # worst case error + + @test maximum(20 * log10.(abs.(diff_EF₂))) < -25 # dB + @test norm(EF₁) == 0.0 + @test maximum(20 * log10.(abs.(diff_HF₂))) < -25 # dB + @test norm(HF₁) == 0.0 + @test maximum(20 * log10.(abs.(diff_FF))) < -25 # dB + end end + @testset "Total fields" begin # define an observation point diff --git a/test/ringCurrents.jl b/test/ringCurrents.jl index f3ad200..282fadd 100644 --- a/test/ringCurrents.jl +++ b/test/ringCurrents.jl @@ -9,10 +9,10 @@ T = assemble(𝑇, RT, RT) @testset "Electric ring current" begin - ex = electricRingCurrent(; frequency=f, center=SVector(0.0, 0.0, 2.0), radius=0.5) - @testset "Incident fields" begin + ex = electricRingCurrent(; frequency=f, center=SVector(0.0, 0.0, 2.0), radius=0.5) + # define an observation point point_cart = [SVector(4.0, 2.0, 3.2), SVector(0.2, 0.1, 2.3)] @@ -35,52 +35,86 @@ T = assemble(𝑇, RT, RT) @testset "Scattered fields" begin - # ----- BEAST solution - 𝐸 = ex + @testset "Standard orientation" begin - 𝑒 = n × 𝐸 × n - #𝑇 = Maxwell3D.singlelayer(; wavenumber=κ) + ex = electricRingCurrent(; frequency=f, center=SVector(0.0, 0.0, 2.0), radius=0.5) - e = assemble(𝑒, RT) - #T = assemble(𝑇, RT, RT) + # ----- BEAST solution + 𝐸 = ex + 𝑒 = n × 𝐸 × n + #𝑇 = Maxwell3D.singlelayer(; wavenumber=κ) - u = T \ e + e = assemble(𝑒, RT) + #T = assemble(𝑇, RT, RT) - EF_MoM = +potential(MWSingleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) - HF_MoM = -potential(BEAST.MWDoubleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) - FF_MoM = -im * f / (2 * c) * potential(MWFarField3D(; gamma=𝑇.gamma), points_cartFF, u, RT) + u = T \ e + EF_MoM = +potential(MWSingleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) + HF_MoM = -potential(BEAST.MWDoubleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) + FF_MoM = -im * f / (2 * c) * potential(MWFarField3D(; gamma=𝑇.gamma), points_cartFF, u, RT) - # ----- this package - sp = PECSphere(; radius=spRadius) + # ----- this package + sp = PECSphere(; radius=spRadius) - EF = scatteredfield(sp, ex, ElectricField(points_cartNF)) - HF = scatteredfield(sp, ex, MagneticField(points_cartNF)) * c * 𝜇 - FF = scatteredfield(sp, ex, FarField(points_cartFF)) + EF = scatteredfield(sp, ex, ElectricField(points_cartNF)) + HF = scatteredfield(sp, ex, MagneticField(points_cartNF)) * c * 𝜇 + FF = scatteredfield(sp, ex, FarField(points_cartFF)) + # ----- compare + diff_EF = norm.(EF - EF_MoM) ./ maximum(norm.(EF)) # worst case error + diff_HF = norm.(HF - HF_MoM) ./ maximum(norm.(HF)) # worst case error + diff_FF = norm.(FF - FF_MoM) ./ maximum(norm.(FF)) # worst case error - # ----- compare - diff_EF = norm.(EF - EF_MoM) ./ maximum(norm.(EF)) # worst case error - diff_HF = norm.(HF - HF_MoM) ./ maximum(norm.(HF)) # worst case error - diff_FF = norm.(FF - FF_MoM) ./ maximum(norm.(FF)) # worst case error + @test maximum(20 * log10.(abs.(diff_EF))) < -24 # dB + @test maximum(20 * log10.(abs.(diff_HF))) < -24 # dB + @test maximum(20 * log10.(abs.(diff_FF))) < -24 # dB + end - @test maximum(20 * log10.(abs.(diff_EF))) < -24 # dB - @test maximum(20 * log10.(abs.(diff_HF))) < -24 # dB - @test maximum(20 * log10.(abs.(diff_FF))) < -24 # dB - end -end + @testset "General orientation" begin + ori = normalize(SVector(0.0, -1.0, -1.0)) + ex = electricRingCurrent(; frequency=f, radius=0.5, center=ori * 2, orientation=ori) -@testset "Magnetic ring current" begin + # ----- BEAST solution + 𝐸 = ex + 𝑒 = n × 𝐸 × n + #𝑇 = Maxwell3D.singlelayer(; wavenumber=κ) + + e = assemble(𝑒, RT) + #T = assemble(𝑇, RT, RT) - #f = 1e8 - #κ = 2π * f / c # Wavenumber + u = T \ e + + EF_MoM = +potential(MWSingleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) + HF_MoM = -potential(BEAST.MWDoubleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) + FF_MoM = -im * f / (2 * c) * potential(MWFarField3D(; gamma=𝑇.gamma), points_cartFF, u, RT) + + # ----- this package + sp = PECSphere(; radius=spRadius) + + EF = scatteredfield(sp, ex, ElectricField(points_cartNF)) + HF = scatteredfield(sp, ex, MagneticField(points_cartNF)) * c * 𝜇 + FF = scatteredfield(sp, ex, FarField(points_cartFF)) + + # ----- compare + diff_EF = norm.(EF - EF_MoM) ./ maximum(norm.(EF)) # worst case error + diff_HF = norm.(HF - HF_MoM) ./ maximum(norm.(HF)) # worst case error + diff_FF = norm.(FF - FF_MoM) ./ maximum(norm.(FF)) # worst case error + + @test maximum(20 * log10.(abs.(diff_EF))) < -24 # dB + @test maximum(20 * log10.(abs.(diff_HF))) < -24 # dB + @test maximum(20 * log10.(abs.(diff_FF))) < -24 # dB + end + end +end - ex = magneticRingCurrent(; frequency=f, center=SVector(0.0, 0.0, 2.0), radius=0.5) +@testset "Magnetic ring current" begin @testset "Incident fields" begin + ex = magneticRingCurrent(; frequency=f, center=SVector(0.0, 0.0, 2.0), radius=0.5) + # define an observation point point_cart = [SVector(4.0, 2.0, 3.2), SVector(0.2, 0.1, 2.3)] @@ -103,37 +137,75 @@ end @testset "Scattered fields" begin - # ----- BEAST solution - 𝐸 = ex + @testset "Standard orientation" begin + + ex = magneticRingCurrent(; frequency=f, center=SVector(0.0, 0.0, 2.0), radius=0.5) + + # ----- BEAST solution + 𝐸 = ex + 𝑒 = n × 𝐸 × n + #𝑇 = Maxwell3D.singlelayer(; wavenumber=κ) + + e = assemble(𝑒, RT) + #T = assemble(𝑇, RT, RT) + + u = T \ e + + EF_MoM = -potential(MWSingleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) + HF_MoM = +potential(BEAST.MWDoubleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) + FF_MoM = +im * f / (2 * c) * potential(MWFarField3D(; gamma=𝑇.gamma), points_cartFF, u, RT) + + # ----- this package + sp = PECSphere(; radius=spRadius) + + EF = scatteredfield(sp, ex, ElectricField(points_cartNF)) + HF = scatteredfield(sp, ex, MagneticField(points_cartNF)) * c * 𝜇 + FF = scatteredfield(sp, ex, FarField(points_cartFF)) + + # ----- compare + diff_EF = norm.(EF - EF_MoM) ./ maximum(norm.(EF)) # worst case error + diff_HF = norm.(HF - HF_MoM) ./ maximum(norm.(HF)) # worst case error + diff_FF = norm.(FF - FF_MoM) ./ maximum(norm.(FF)) # worst case error + + @test maximum(20 * log10.(abs.(diff_EF))) < -24 # dB + @test maximum(20 * log10.(abs.(diff_HF))) < -24 # dB + @test maximum(20 * log10.(abs.(diff_FF))) < -24 # dB + end - 𝑒 = n × 𝐸 × n - #𝑇 = Maxwell3D.singlelayer(; wavenumber=κ) + @testset "General orientation" begin - e = assemble(𝑒, RT) - #T = assemble(𝑇, RT, RT) + ori = normalize(SVector(0.0, -1.0, -1.0)) + ex = magneticRingCurrent(; frequency=f, radius=0.5, center=ori * 2, orientation=ori) - u = T \ e + # ----- BEAST solution + 𝐸 = ex + 𝑒 = n × 𝐸 × n + #𝑇 = Maxwell3D.singlelayer(; wavenumber=κ) - EF_MoM = -potential(MWSingleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) - HF_MoM = +potential(BEAST.MWDoubleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) - FF_MoM = +im * f / (2 * c) * potential(MWFarField3D(; gamma=𝑇.gamma), points_cartFF, u, RT) + e = assemble(𝑒, RT) + #T = assemble(𝑇, RT, RT) + u = T \ e - # ----- this package - sp = PECSphere(; radius=spRadius) + EF_MoM = -potential(MWSingleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) + HF_MoM = +potential(BEAST.MWDoubleLayerField3D(; wavenumber=κ), points_cartNF, u, RT) + FF_MoM = +im * f / (2 * c) * potential(MWFarField3D(; gamma=𝑇.gamma), points_cartFF, u, RT) - EF = scatteredfield(sp, ex, ElectricField(points_cartNF)) - HF = scatteredfield(sp, ex, MagneticField(points_cartNF)) * c * 𝜇 - FF = scatteredfield(sp, ex, FarField(points_cartFF)) + # ----- this package + sp = PECSphere(; radius=spRadius) + EF = scatteredfield(sp, ex, ElectricField(points_cartNF)) + HF = scatteredfield(sp, ex, MagneticField(points_cartNF)) * c * 𝜇 + FF = scatteredfield(sp, ex, FarField(points_cartFF)) - # ----- compare - diff_EF = norm.(EF - EF_MoM) ./ maximum(norm.(EF)) # worst case error - diff_HF = norm.(HF - HF_MoM) ./ maximum(norm.(HF)) # worst case error - diff_FF = norm.(FF - FF_MoM) ./ maximum(norm.(FF)) # worst case error + # ----- compare + diff_EF = norm.(EF - EF_MoM) ./ maximum(norm.(EF)) # worst case error + diff_HF = norm.(HF - HF_MoM) ./ maximum(norm.(HF)) # worst case error + diff_FF = norm.(FF - FF_MoM) ./ maximum(norm.(FF)) # worst case error - @test maximum(20 * log10.(abs.(diff_EF))) < -24 # dB - @test maximum(20 * log10.(abs.(diff_HF))) < -24 # dB - @test maximum(20 * log10.(abs.(diff_FF))) < -24 # dB + @test maximum(20 * log10.(abs.(diff_EF))) < -24 # dB + @test maximum(20 * log10.(abs.(diff_HF))) < -24 # dB + @test maximum(20 * log10.(abs.(diff_FF))) < -24 # dB + end end end