implement general orientations of excitations

Project.toml

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

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
+<br/>
+```
+
+#### Ring-Currents and Dipoles
 
-SI units are employed everywhere.
+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}``.

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
 

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)
 
 
 ---

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
 

docs/src/ringCurrents.md

@@ -23,19 +23,24 @@ 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)\,,
 ```
 where ``\delta`` denotes the Dirac delta distribution and ``R_0 = \sqrt{z_0^2 + a^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
 

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
 

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

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)