# flake8: noqa
from pkg_resources import resource_filename
from functools import lru_cache
import warnings
import numpy as np
from ...matlab_funcs import besselh, besselj, gammaln, lscov, quadl
from ...sci_funcs import legendrePlm
from ...core import stress2legendre
[docs]def boundary(costheta, a=1, epsilon=.1, nu=0):
"""Projected boundary of a prolate spheroid
Compute the boundary according to equation (4) in
:cite:`Boyde2009` with the addition of the
Poisson's ratio of the object.
.. math::
B(\\theta) = a (1+\\epsilon)
\\left[ (1+\\epsilon)^2 - \\epsilon (1+\\nu)
(2+\\epsilon (1-\\nu)) \\cos^2 \\theta \\right]^{-1/2}
This boundary function was derived for a prolate spheroid under
the assumption that the semi-major axis :math:`a` and the
semi-minor axes :math:`b=c` are defined as
.. math::
a = b \\cdot \\frac{1+ \\epsilon}{1- \\nu \\epsilon}
The boundary function :math:`B(\\theta)` can be derived with
the above relation using the equation for a prolate spheroid.
Parameters
----------
costheta: float or np.ndarray
Cosine of polar coordinates :math:`\\theta`
at which to compute the boundary.
a: float
Equatorial radii of prolate spheroid (semi-minor axis).
epsilon: float
Stretch ratio; defines size of semi-major axis:
:math:`a = (1+\\epsilon) b`. Note that this is not
the eccentricity of the prolate spheroid.
nu: float
Poisson's ratio :math:`\\nu` of the material.
Returns
-------
B: 1d ndarray
Radial object boundary in dependence of theta
:math:`B(\\theta)`.
Notes
-----
For :math:`\\nu=0`, the above equation becomes
equation (4) in :cite:`Boyde2009`.
"""
x = costheta
B = a * (1 + epsilon) \
/ ((1 + epsilon)**2
- epsilon * (1 + nu) * (2 + epsilon * (1 - nu)) * x**2)**.5
return B
[docs]@lru_cache(maxsize=32)
def get_hgc():
"""Load hypergeometric coefficients from *hypergeomdata2.dat*.
These coefficients were computed by Lars Boyde
using Wolfram Mathematica.
"""
hpath = resource_filename("ggf.stress.boyde2009", "hypergeomdata2.dat")
hgc = np.loadtxt(hpath)
return hgc
[docs]def stress(object_index=1.41, medium_index=1.3465, poisson_ratio=0.45,
semi_minor=2.8466e-6, stretch_ratio=0.1, wavelength=780e-9,
beam_waist=3, power_left=.6, power_right=.6, dist=100e-6,
n_points=100, theta_max=np.pi, field_approx="davis",
ret_legendre_decomp=False, verbose=False):
"""Compute the stress acting on a prolate spheroid
The prolate spheroid has semi-major axis :math:`a` and
semi-minor axis :math:`b=c`.
Parameters
----------
object_index: float
Refractive index of the spheroid
medium_index: float
Refractive index of the surrounding medium
poisson_ratio: float
Poisson's ratio of the spheroid material
semi_minor: float
Semi-minor axis (inner) radius of the stretched object
:math:`b=c`.
stretch_ratio: float
Measure of the deformation, defined as :math:`(a - b) / b`
wavelength: float
Wavelenth of the gaussian beam [m]
beam_waist: float
Beam waist radius of the gaussian beam [wavelengths]
power_left: float
Laser power of the left beam [W]
power_right: float
Laser power of the right beam [W]
dist: float
Distance between beam waist and object center [m]
n_points: int
Number of points to compute stresses for
theta_max: float
Maximum angle to compute stressed for
field_approx: str
TODO
ret_legendre_decomp: bool
If True, return coefficients of decomposition of stress
into Legendre polynomials
verbose: int
Increase verbosity
Returns
-------
theta: 1d ndarray
Angles for which stresses are computed
sigma_rr: 1d ndarray
Radial stress corresponding to angles
coeff: 1d ndarray
If `ret_legendre_decomp` is True, return compositions
of decomposition of stress into Legendre polynomials.
Notes
-----
- The angles `theta` are computed on a grid that does not
include zero and `theta_max`.
- This implementation was first presented in :cite:`Boyde2009`.
"""
if field_approx not in ["davis", "barton"]:
raise ValueError("`field_approx` must be 'davis' or 'barton'")
object_index = complex(object_index)
medium_index = complex(medium_index)
W0 = beam_waist * wavelength
epsilon = stretch_ratio
nu = poisson_ratio
# ZRL = 0.5*medium_index*2*np.pi/wavelength*W0**2 # Rayleigh range [m]
# WZ = W0*(1+(beam_pos+d)**2/ZRL**2)**0.5 # beam waist at specified
# position [m]
K0 = 2 * np.pi / wavelength # wave vector [m]
Alpha = semi_minor * K0 # size parameter
C = 3e8 # speed of light [m/s]
# maximum number of orders
lmax = np.int(np.round(2 + Alpha + 4 * (Alpha)**(1 / 3) + 10))
if lmax > 120:
msg = 'Required number of orders for accurate expansion exceeds allowed maximum!' \
+ 'Reduce size of trapped particle!'
raise ValueError(msg)
if epsilon == 0:
# spherical object, no point-matching needed (mmax = 0)
mmax = 3
else:
if (epsilon > 0.15):
warnings.warn(
'Warning, cell stretching ratio is high: {}'.format(epsilon))
# spheroidal object, point-matching required (mmax has to be divisible
# by 3)
mmax = 6 * lmax
# permittivity in surrounding medium [1]
EpsilonI = medium_index**2
EpsilonII = object_index**2 # permittivity in within cell [1]
MuI = 1.000 # permeability in surrounding medium [1]
MuII = 1.000 # permeability within cell [1]
# wave constant in Maxwell's equations (surrounding medium) [1/m]
K1I = 1j * K0 * EpsilonI
# wave constant in Maxwell's equations (within cell) [1/m]
K1II = 1j * K0 * EpsilonII
# wave constant in Maxwell's equations (surrounding medium) [1/m]
K2I = 1j * K0
# wave constant in Maxwell's equations (within cell) [1/m]
K2II = 1j * K0
KI = (-K1I * K2I)**0.5 # wave vector (surrounding medium) [1/m]
KII = (-K1II * K2II)**0.5 # wave vector (within cell) [1/m]
# dimensionless parameters
k0 = 1 # wave vector
a = semi_minor * K0 # internal radius of stretched cell
d = dist * K0 # distance from cell centre to optical stretcher
# ap = a*(1+stretch_ratio) # semi-major axis (after stretching)
# bp = a*(1-poisson_ratio*stretch_ratio) # semi-minor axis (after
# stretching)
w0 = W0 * K0 # Gaussian width
# wave constant in Maxwell's equations (surrounding medium)
k1I = K1I / K0
# wave constant in Maxwell's equations (within cell)
k1II = K1II / K0
# wave constant in Maxwell's equations (surrounding medium)
k2I = K2I / K0
# wave constant in Maxwell's equations (within cell)
k2II = K2II / K0
kI = KI / K0 # wave vector (surrounding medium)
kII = KII / K0 # wave vector (within cell)
beta = kI # wave vector of Gaussian beam
# other definitions
# amplitude of electric field of left laser [kg m/(s**2 C)]
EL = np.sqrt(power_left / (medium_index * C * W0**2))
# amplitude of electric field of right laser [kg m/(s**2 C)]
ER = np.sqrt(power_right / (medium_index * C * W0**2))
HL = beta / k0 * EL # left laser amplitude of magnetic field
HR = beta / k0 * ER # right laser amplitude of magnetic field
zR = beta * w0**2 / 2 # definition of Rayleigh range
S = (1 + 1j * d / zR)**(-1) # complex amplitude for Taylor expansion
s = 1 / (beta * w0) # expansion parameter for Gaussian (Barton)
# Functions
# object boundary function: r(th) = a*B1(x) x= cos(th)
def B1(x): return boundary(costheta=x, a=1, epsilon=epsilon, nu=nu)
# Riccati Bessel functions and their derivatives
# Riccati Bessel function (psi)
def psi(l, z): return (np.pi / 2 * z)**(1 / 2) * besselj(l + 1 / 2, z)
def psi1(l, z): return (np.pi / (2. * z))**(1 / 2) * \
(z * besselj(l - 1 / 2, z) - l *
besselj(l + 1 / 2, z)) # first derivative (psi')
def psi2(l, z): return (np.pi / 2)**(1 / 2) * (l + l**2 - z**2) * \
besselj(l + 1 / 2, z) * z**(-3 / 2) # second derivative (psi'')
# First order Taylor expansion of psi is too inaccurate for larger values of k*a*Eps.
# Hence, to match 1-st and higher order terms in Eps, subtract the 0-th order terms (no angular dependence)
# from the exact function psi (including angular dependence)
# Riccati Bessel function excluding angular dependence in 0-th order
# (psiex)
def psiex(l, z, x): return psi(l, z * B1(x)) - psi(l, z)
def psi1ex(l, z, x): return psi1(l, z * B1(x)) - \
psi1(l, z) # first derivative of psiex
def psi2ex(l, z, x): return psi2(l, z * B1(x)) - \
psi2(l, z) # second derivative of psi
# defined for abbreviation
def psixx(l, z, x): return psi(l, z * B1(x))
def psi1xx(l, z, x): return psi1(l, z * B1(x))
def psi2xx(l, z, x): return psi2(l, z * B1(x))
# Hankel function and its derivative
def xi(l, z): return (np.pi / 2 * z)**(1 / 2) * besselh(l + 1 / 2, z)
def xi1(l, z): return (np.pi / (2 * z))**(1 / 2) * \
((l + 1) * besselh(l + 1 / 2, z) - z * besselh(l + 3 / 2, z))
def xi2(l, z): return (np.pi / 2)**(1 / 2) / \
z**(3 / 2) * (l + l**2 - z**2) * besselh(l + 1 / 2, z)
# Comments: see above for psiex
def xiex(l, z, x): return xi(l, z * B1(x)) - xi(l, z)
def xi1ex(l, z, x): return xi1(l, z * B1(x)) - xi1(l, z)
def xi2ex(l, z, x): return xi2(l, z * B1(x)) - xi2(l, z)
def xixx(l, z, x): return xi(l, z * B1(x))
def xi1xx(l, z, x): return xi1(l, z * B1(x))
def xi2xx(l, z, x): return xi2(l, z * B1(x))
#% Associated Legendre functions P(m)_l(x) and their derivatives
#% select mth component of vector 'legendre' [P**(m)_l(x)]
# [zeros(m,1);1;zeros(l-m,1)].'*legendre(l,x)
#% legendre polynomial [P**(m)_l(x)]
def legendrePl(l, x): return legendrePlm(1, l, x)
#% legendre polynomial [P**(1)_l(x)]
def legendrePlm1(m, l, x): return (
(l - m + 1.) * legendrePlm(m, l + 1, x) - (l + 1.) * x * legendrePlm(m, l, x)) / (x**2 - 1)
#% derivative d/dx[P**(m)_l(x)]
def legendrePl1(l, x): return legendrePlm1(1, l, x)
# defined to avoid division by zero (which can occur for x=1 in
# legendrePl1...
def legendrePlmex1(m, l, x): return -((l - m + 1) *
legendrePlm(m, l + 1, x) - (l + 1) * x * legendrePlm(m, l, x))
def legendrePlex1(l, x): return legendrePlmex1(1, l, x)
# Hypergeometric and Gamma functions
hypergeomcoeff = get_hgc()
## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ##
# Gaussian beam (incident fields) - in Cartesian basis (either Davis first order or Barton fifth order fields)
# electric and magnetic fields according to Davis (first order)
if field_approx == "davis":
# left
def eExiL(r, th, phi): return EL * (1 + 1j * (r * np.cos(th) + d) / zR)**(-1) * np.exp(-r**2 *
np.sin(th)**2 / (w0**2 * (1 + 1j * (r * np.cos(th) + d) / zR))) * np.exp(1j * beta * (r * np.cos(th) + d))
def eEyiL(r, th, phi): return 0
def eEziL(r, th, phi): return -1j * (1 + 1j * (r * np.cos(th) + d) /
zR)**(-1) * r * np.sin(th) * np.cos(phi) / zR * eExiL(r, th, phi)
def eHxiL(r, th, phi): return 0
def eHyiL(r, th, phi): return HL * (1 + 1j * (r * np.cos(th) + d) / zR)**(-1) * np.exp(-r**2 *
np.sin(th)**2 / (w0**2 * (1 + 1j * (r * np.cos(th) + d) / zR))) * np.exp(1j * beta * (r * np.cos(th) + d))
def eHziL(r, th, phi): return -1j * (1 + 1j * (r * np.cos(th) + d) /
zR)**(-1) * r * np.sin(th) * np.sin(phi) / zR * eHyiL(r, th, phi)
# right
def eExiR(r, th, phi): return ER * (1 - 1j * (r * np.cos(th) - d) / zR)**(-1) * np.exp(-r**2 *
np.sin(th)**2 / (w0**2 * (1 - 1j * (r * np.cos(th) - d) / zR))) * np.exp(-1j * beta * (r * np.cos(th) - d))
def eEyiR(r, th, phi): return 0
def eEziR(r, th, phi): return +1j * (1 - 1j * (r * np.cos(th) - d) /
zR)**(-1) * r * np.sin(th) * np.cos(phi) / zR * eExiR(r, th, phi)
def eHxiR(r, th, phi): return 0
def eHyiR(r, th, phi): return -HR * (1 - 1j * (r * np.cos(th) - d) / zR)**(-1) * np.exp(-r**2 *
np.sin(th)**2 / (w0**2 * (1 - 1j * (r * np.cos(th) - d) / zR))) * np.exp(-1j * beta * (r * np.cos(th) - d))
def eHziR(r, th, phi): return +1j * (1 - 1j * (r * np.cos(th) - d) /
zR)**(-1) * r * np.sin(th) * np.sin(phi) / zR * eHyiR(r, th, phi)
else: # electric and magnetic fields according to Barton (fifth order)
# Note: in Barton propagation is: np.exp(i (omega*t-k*z)) and not np.exp(i (k*z-omega*t)).
# Hence, take complex conjugate of equations or make the changes z ->
# -z, Ez -> -Ez, Hz -> -Hz while x and y are not affected.
def Xi(r, th, phi): return r * np.sin(th) * np.cos(phi) / w0
def Eta(r, th, phi): return r * np.sin(th) * np.sin(phi) / w0
def ZetaL(r, th): return (r * np.cos(th) + d) / (beta * w0**2)
def Rho(r, th, phi): return np.sqrt(
Xi(r, th, phi)**2 + Eta(r, th, phi)**2)
def QL(r, th): return 1. / (1j + 2 * ZetaL(r, th))
def Psi0L(r, th, phi): return 1j * QL(r, th) * \
np.exp(-1j * (Rho(r, th, phi))**2 * QL(r, th))
def eExiL(r, th, phi): return np.conj(EL * Psi0L(r, th, phi) * np.exp(-1j * ZetaL(r, th) / s**2) *
(1 + s**2 * (-Rho(r, th, phi)**2 * QL(r, th)**2 + 1j * Rho(r, th, phi)**4 * QL(r, th)**3 - 2 * QL(r, th)**2 * Xi(r, th, phi)**2) +
s**4 * (2 * Rho(r, th, phi)**4 * QL(r, th)**4 - 3 * 1j * Rho(r, th, phi)**6 * QL(r, th)**5 - 0.5 * Rho(r, th, phi)**8 * QL(r, th)**6 +
(8 * Rho(r, th, phi)**2 * QL(r, th)**4 - 2 * 1j * Rho(r, th, phi)**4 * QL(r, th)**5) * Xi(r, th, phi)**2)))
def eEyiL(r, th, phi): return np.conj(EL * Psi0L(r, th, phi) * np.exp(-1j * ZetaL(r, th) / s**2) *
(s**2 * (-2 * QL(r, th)**2 * Xi(r, th, phi) * Eta(r, th, phi)) +
s**4 * ((8 * Rho(r, th, phi)**2 * QL(r, th)**4 - 2 * 1j * Rho(r, th, phi)**4 * QL(r, th)**5) * Xi(r, th, phi) * Eta(r, th, phi))))
def eEziL(r, th, phi): return np.conj(EL * Psi0L(r, th, phi) * np.exp(-1j * ZetaL(r, th) / s**2) *
(s * (-2 * QL(r, th) * Xi(r, th, phi)) + s**3 * (6 * Rho(r, th, phi)**2 * QL(r, th)**3 - 2 * 1j * Rho(r, th, phi)**4 * QL(r, th)**4) * Xi(r, th, phi) +
s**5 * (-20 * Rho(r, th, phi)**4 * QL(r, th)**5 + 10 * 1j * Rho(r, th, phi)**6 * QL(r, th)**6 + Rho(r, th, phi)**8 * QL(r, th)**7) * Xi(r, th, phi)))
def eHxiL(r, th, phi): return np.conj(+np.sqrt(EpsilonI) * EL * Psi0L(r, th, phi) * np.exp(-1j * ZetaL(r, th) / s**2) *
(s**2 * (-2 * QL(r, th)**2 * Xi(r, th, phi) * Eta(r, th, phi)) +
s**4 * ((8 * Rho(r, th, phi)**2 * QL(r, th)**4 - 2 * 1j * Rho(r, th, phi)**4 * QL(r, th)**5) * Xi(r, th, phi) * Eta(r, th, phi))))
def eHyiL(r, th, phi): return np.conj(+np.sqrt(EpsilonI) * EL * Psi0L(r, th, phi) * np.exp(-1j * ZetaL(r, th) / s**2) *
(1 + s**2 * (-Rho(r, th, phi)**2 * QL(r, th)**2 + 1j * Rho(r, th, phi)**4 * QL(r, th)**3 - 2 * QL(r, th)**2 * Eta(r, th, phi)**2) +
s**4 * (2 * Rho(r, th, phi)**4 * QL(r, th)**4 - 3 * 1j * Rho(r, th, phi)**6 * QL(r, th)**5 - 0.5 * Rho(r, th, phi)**8 * QL(r, th)**6 +
(8 * Rho(r, th, phi)**2 * QL(r, th)**4 - 2 * 1j * Rho(r, th, phi)**4 * QL(r, th)**5) * Eta(r, th, phi)**2)))
def eHziL(r, th, phi): return np.conj(np.sqrt(EpsilonI) * EL * Psi0L(r, th, phi) * np.exp(-1j * ZetaL(r, th) / s**2) *
(s * (-2 * QL(r, th) * Eta(r, th, phi)) + s**3 * (6 * Rho(r, th, phi)**2 * QL(r, th)**3 - 2 * 1j * Rho(r, th, phi)**4 * QL(r, th)**4) * Eta(r, th, phi) +
s**5 * (-20 * Rho(r, th, phi)**4 * QL(r, th)**5 + 10 * 1j * Rho(r, th, phi)**6 * QL(r, th)**6 + Rho(r, th, phi)**8 * QL(r, th)**7) * Eta(r, th, phi)))
# right
# Take left fiber fields and make coordinate changes (x,y,z,d) ->
# (x,-y,-z,d) and amplitude changes (Ex,Ey,Ez) -> (Ex,-Ey,-Ez).
def ZetaR(r, th): return -(r * np.cos(th) - d) / (beta * w0**2)
def QR(r, th): return 1. / (1j + 2 * ZetaR(r, th))
def Psi0R(r, th, phi): return 1j * QR(r, th) * \
np.exp(-1j * (Rho(r, th, phi))**2 * QR(r, th))
def eExiR(r, th, phi): return np.conj(ER * Psi0R(r, th, phi) * np.exp(-1j * ZetaR(r, th) / s**2) *
(1 + s**2 * (-Rho(r, th, phi)**2 * QR(r, th)**2 + 1j * Rho(r, th, phi)**4 * QR(r, th)**3 - 2 * QR(r, th)**2 * Xi(r, th, phi)**2) +
s**4 * (2 * Rho(r, th, phi)**4 * QR(r, th)**4 - 3 * 1j * Rho(r, th, phi)**6 * QR(r, th)**5 - 0.5 * Rho(r, th, phi)**8 * QR(r, th)**6 +
(8 * Rho(r, th, phi)**2 * QR(r, th)**4 - 2 * 1j * Rho(r, th, phi)**4 * QR(r, th)**5) * Xi(r, th, phi)**2)))
def eEyiR(r, th, phi): return np.conj(ER * Psi0R(r, th, phi) * np.exp(-1j * ZetaR(r, th) / s**2) *
(s**2 * (-2 * QR(r, th)**2 * Xi(r, th, phi) * Eta(r, th, phi)) +
s**4 * ((8 * Rho(r, th, phi)**2 * QR(r, th)**4 - 2 * 1j * Rho(r, th, phi)**4 * QR(r, th)**5) * Xi(r, th, phi) * Eta(r, th, phi))))
def eEziR(r, th, phi): return - np.conj(ER * Psi0R(r, th, phi) * np.exp(-1j * ZetaR(r, th) / s**2) *
(s * (-2 * QR(r, th) * Xi(r, th, phi)) + s**3 * (6 * Rho(r, th, phi)**2 * QR(r, th)**3 - 2 * 1j * Rho(r, th, phi)**4 * QR(r, th)**4) * Xi(r, th, phi) +
s**5 * (-20 * Rho(r, th, phi)**4 * QR(r, th)**5 + 10 * 1j * Rho(r, th, phi)**6 * QR(r, th)**6 + Rho(r, th, phi)**8 * QR(r, th)**7) * Xi(r, th, phi)))
def eHxiR(r, th, phi): return - np.conj(+np.sqrt(EpsilonI) * ER * Psi0R(r, th, phi) * np.exp(-1j * ZetaR(r, th) / s**2) *
(s**2 * (-2 * QR(r, th)**2 * Xi(r, th, phi) * Eta(r, th, phi)) +
s**4 * ((8 * Rho(r, th, phi)**2 * QR(r, th)**4 - 2 * 1j * Rho(r, th, phi)**4 * QR(r, th)**5) * Xi(r, th, phi) * Eta(r, th, phi))))
def eHyiR(r, th, phi): return - np.conj(+np.sqrt(EpsilonI) * ER * Psi0R(r, th, phi) * np.exp(-1j * ZetaR(r, th) / s**2) *
(1 + s**2 * (-Rho(r, th, phi)**2 * QR(r, th)**2 + 1j * Rho(r, th, phi)**4 * QR(r, th)**3 - 2 * QR(r, th)**2 * Eta(r, th, phi)**2) +
s**4 * (2 * Rho(r, th, phi)**4 * QR(r, th)**4 - 3 * 1j * Rho(r, th, phi)**6 * QR(r, th)**5 - 0.5 * Rho(r, th, phi)**8 * QR(r, th)**6 +
(8 * Rho(r, th, phi)**2 * QR(r, th)**4 - 2 * 1j * Rho(r, th, phi)**4 * QR(r, th)**5) * Eta(r, th, phi)**2)))
def eHziR(r, th, phi): return np.conj(np.sqrt(EpsilonI) * ER * Psi0R(r, th, phi) * np.exp(-1j * ZetaR(r, th) / s**2) *
(s * (-2 * QR(r, th) * Eta(r, th, phi)) + s**3 * (6 * Rho(r, th, phi)**2 * QR(r, th)**3 - 2 * 1j * Rho(r, th, phi)**4 * QR(r, th)**4) * Eta(r, th, phi) +
s**5 * (-20 * Rho(r, th, phi)**4 * QR(r, th)**5 + 10 * 1j * Rho(r, th, phi)**6 * QR(r, th)**6 + Rho(r, th, phi)**8 * QR(r, th)**7) * Eta(r, th, phi)))
# Gaussian beam (incident fields) - in spherical polar coordinates basis
# left
def eEriL(r, th, phi): return np.sin(th) * np.cos(phi) * eExiL(r, th, phi) + \
np.sin(th) * np.sin(phi) * eEyiL(r, th, phi) + \
np.cos(th) * eEziL(r, th, phi)
def eEthiL(r, th, phi): return np.cos(th) * np.cos(phi) * eExiL(r, th, phi) + \
np.cos(th) * np.sin(phi) * eEyiL(r, th, phi) - \
np.sin(th) * eEziL(r, th, phi)
def eEphiiL(r, th, phi): return - np.sin(phi) * \
eExiL(r, th, phi) + np.cos(phi) * eEyiL(r, th, phi)
def eHriL(r, th, phi): return np.sin(th) * np.cos(phi) * eHxiL(r, th, phi) + \
np.sin(th) * np.sin(phi) * eHyiL(r, th, phi) + \
np.cos(th) * eHziL(r, th, phi)
def eHthiL(r, th, phi): return np.cos(th) * np.cos(phi) * eHxiL(r, th, phi) + \
np.cos(th) * np.sin(phi) * eHyiL(r, th, phi) - \
np.sin(th) * eHziL(r, th, phi)
def eHphiiL(r, th, phi): return - np.sin(phi) * \
eHxiL(r, th, phi) + np.cos(phi) * eHyiL(r, th, phi)
# right
def eEriR(r, th, phi): return np.sin(th) * np.cos(phi) * eExiR(r, th, phi) + \
np.sin(th) * np.sin(phi) * eEyiR(r, th, phi) + \
np.cos(th) * eEziR(r, th, phi)
def eEthiR(r, th, phi): return np.cos(th) * np.cos(phi) * eExiR(r, th, phi) + \
np.cos(th) * np.sin(phi) * eEyiR(r, th, phi) - \
np.sin(th) * eEziR(r, th, phi)
def eEphiiR(r, th, phi): return - np.sin(phi) * \
eExiR(r, th, phi) + np.cos(phi) * eEyiR(r, th, phi)
def eHriR(r, th, phi): return np.sin(th) * np.cos(phi) * eHxiR(r, th, phi) + \
np.sin(th) * np.sin(phi) * eHyiR(r, th, phi) + \
np.cos(th) * eHziR(r, th, phi)
def eHthiR(r, th, phi): return np.cos(th) * np.cos(phi) * eHxiR(r, th, phi) + \
np.cos(th) * np.sin(phi) * eHyiR(r, th, phi) - \
np.sin(th) * eHziR(r, th, phi)
def eHphiiR(r, th, phi): return - np.sin(phi) * \
eHxiR(r, th, phi) + np.cos(phi) * eHyiR(r, th, phi)
eBL = np.zeros(lmax, dtype=complex)
mBL = np.zeros(lmax, dtype=complex)
eBR = np.zeros(lmax, dtype=complex)
# eBR_test = np.zeros(lmax, dtype=complex)
mBR = np.zeros(lmax, dtype=complex)
if field_approx == "davis": # Davis
tau = np.zeros(lmax, dtype=complex)
for ii in range(lmax):
l = ii + 1
tau[ii] = 0 # sum over m and n
niimax = int(np.floor((l - 1) / 3))
for n in range(niimax + 1):
miimax = int(np.floor((l - 1) / 2 - 3 / 2 * n))
for m in range(miimax + 1):
if l < (2 * m + 3 * n + 2):
Delta = 0
else:
Delta = 1
tau[ii] = tau[ii] + np.exp(gammaln(l / 2 - m - n) + gammaln(m + n + 2) - gammaln(l - 2 * m - 3 * n) - gammaln(l / 2 + 2) - gammaln(m + 1) - gammaln(n + 1)) * hypergeomcoeff[l - 1, m + n] \
* (-1)**(m + 1) / (S**m * (beta * zR)**n) * (S / (1j * beta * w0))**(2 * m + 2 * n) * (1 - Delta * 2 * S / (beta * zR) * (l - 1 - 2 * m - 3 * n))
# print(tau)
# calculate expansion coefficients of order l for electric and magnetic
# Debye potentials
def emB(l): return S * np.exp(1j * beta * d) * (1j * beta / (2 * kI))**(l - 1) * \
(l + 1 / 2)**2 / (l + 1) * \
np.exp(gammaln(2 * l) - gammaln(l + 1)) * tau[l - 1]
for ii in range(lmax):
l = ii + 1
# left
eBL[ii] = EL * emB(l)
mBL[ii] = HL * emB(l)
# right
# should include factor of (-1)**(l-1) for symmetry reasons, this
# is taken into account only after least square fitting (see below)
eBR[ii] = ER * emB(l)
# should include factor of (-1)**l for symmetry reasons, this is
# taken into account only after least square fitting (see
mBR[ii] = HR * emB(l)
else: # Barton
for ii in range(lmax):
l = ii + 1
eBL[ii] = np.sqrt(2) * quadl(lambda th: eEriL(a, th, np.pi / 4) * legendrePl(l, np.cos(
th)) * np.sin(th), 0, np.pi) * (2 * l + 1) / (2 * l * (l + 1)) * a**2 / psi(l, kI * a) * kI**2 / (l * (l + 1))
mBL[ii] = np.sqrt(2) * quadl(lambda th: eHriL(a, th, np.pi / 4) * legendrePl(l, np.cos(
th)) * np.sin(th), 0, np.pi) * (2 * l + 1) / (2 * l * (l + 1)) * a**2 / psi(l, kI * a) * kI**2 / (l * (l + 1))
eBR[ii] = np.sqrt(2) * quadl(lambda th: eEriR(a, th, np.pi / 4) * legendrePl(l, np.cos(
th)) * np.sin(th), 0, np.pi) * (2 * l + 1) / (2 * l * (l + 1)) * a**2 / psi(l, kI * a) * kI**2 / (l * (l + 1))
mBR[ii] = np.sqrt(2) * quadl(lambda th: eHriR(a, th, np.pi / 4) * legendrePl(l, np.cos(
th)) * np.sin(th), 0, np.pi) * (2 * l + 1) / (2 * l * (l + 1)) * a**2 / psi(l, kI * a) * kI**2 / (l * (l + 1))
# make symetrical with left expansion coefficients
eBR[ii] = eBR[ii] * (-1)**(l - 1)
# make symetrical with left expansion coefficients
mBR[ii] = mBR[ii] * (-1)**l
# coefficients for internal fields (eCl, mCl) and scattered fields (eDl,
# mDl)
eCL = np.zeros(lmax, dtype=complex)
mCL = np.zeros(lmax, dtype=complex)
eCR = np.zeros(lmax, dtype=complex)
mCR = np.zeros(lmax, dtype=complex)
eDL = np.zeros(lmax, dtype=complex)
mDL = np.zeros(lmax, dtype=complex)
eDR = np.zeros(lmax, dtype=complex)
mDR = np.zeros(lmax, dtype=complex)
for ii in range(lmax):
l = ii + 1
# internal (left and right)
eCL[ii] = k1I / kI * (kII)**2 * (xi(l, kI * a) * psi1(l, kI * a) - xi1(l, kI * a) * psi(l, kI * a)) / (
k1I * kII * xi(l, kI * a) * psi1(l, kII * a) - k1II * kI * xi1(l, kI * a) * psi(l, kII * a)) * eBL[ii]
mCL[ii] = k2I / kI * (kII)**2 * (xi(l, kI * a) * psi1(l, kI * a) - xi1(l, kI * a) * psi(l, kI * a)) / (
k2I * kII * xi(l, kI * a) * psi1(l, kII * a) - k2II * kI * xi1(l, kI * a) * psi(l, kII * a)) * mBL[ii]
eCR[ii] = k1I / kI * (kII)**2 * (xi(l, kI * a) * psi1(l, kI * a) - xi1(l, kI * a) * psi(l, kI * a)) / (
k1I * kII * xi(l, kI * a) * psi1(l, kII * a) - k1II * kI * xi1(l, kI * a) * psi(l, kII * a)) * eBR[ii]
mCR[ii] = k2I / kI * (kII)**2 * (xi(l, kI * a) * psi1(l, kI * a) - xi1(l, kI * a) * psi(l, kI * a)) / (
k2I * kII * xi(l, kI * a) * psi1(l, kII * a) - k2II * kI * xi1(l, kI * a) * psi(l, kII * a)) * mBR[ii]
# scattered (left and right)
eDL[ii] = (k1I * kII * psi(l, kI * a) * psi1(l, kII * a) - k1II * kI * psi1(l, kI * a) * psi(l, kII * a)) / \
(k1II * kI * xi1(l, kI * a) * psi(l, kII * a) - k1I *
kII * xi(l, kI * a) * psi1(l, kII * a)) * eBL[ii]
mDL[ii] = (k2I * kII * psi(l, kI * a) * psi1(l, kII * a) - k2II * kI * psi1(l, kI * a) * psi(l, kII * a)) / \
(k2II * kI * xi1(l, kI * a) * psi(l, kII * a) - k2I *
kII * xi(l, kI * a) * psi1(l, kII * a)) * mBL[ii]
eDR[ii] = (k1I * kII * psi(l, kI * a) * psi1(l, kII * a) - k1II * kI * psi1(l, kI * a) * psi(l, kII * a)) / \
(k1II * kI * xi1(l, kI * a) * psi(l, kII * a) - k1I *
kII * xi(l, kI * a) * psi1(l, kII * a)) * eBR[ii]
mDR[ii] = (k2I * kII * psi(l, kI * a) * psi1(l, kII * a) - k2II * kI * psi1(l, kI * a) * psi(l, kII * a)) / \
(k2II * kI * xi1(l, kI * a) * psi(l, kII * a) - k2I *
kII * xi(l, kI * a) * psi1(l, kII * a)) * mBR[ii]
# First Order Expansion Coefficients
# coefficients for internal fields (eCcl, mCcl) and scattered fields
# (eDdl, mDdl)
eLambda1L = {}
mLambda1L = {}
eLambda2L = {}
mLambda2L = {}
eLambda3L = {}
mLambda3L = {}
eLambda1R = {}
mLambda1R = {}
eLambda2R = {}
mLambda2R = {}
eLambda3R = {}
mLambda3R = {}
for jj in range(lmax):
l = jj + 1
# left
eLambda1L[l] = lambda x, l=l, jj=jj: (eBL[jj] / kI * psi1ex(l, kI * a, x) - eCL[jj] / kII * psi1ex(
l, kII * a, x) + eDL[jj] / kI * xi1ex(l, kI * a, x)) # electric parameter1 left
mLambda1L[l] = lambda x, l=l, jj=jj: (mBL[jj] / kI * psi1ex(l, kI * a, x) - mCL[jj] / kII * psi1ex(
l, kII * a, x) + mDL[jj] / kI * xi1ex(l, kI * a, x)) # magnetic parameter1 left
eLambda2L[l] = lambda x, l=l, jj=jj: (k1I / kI**2 * eBL[jj] * psiex(l, kI * a, x) - k1II / kII**2 * eCL[jj] * psiex(
l, kII * a, x) + k1I / kI**2 * eDL[jj] * xiex(l, kI * a, x)) # electric parameter2 left
mLambda2L[l] = lambda x, l=l, jj=jj: (k2I / kI**2 * mBL[jj] * psiex(l, kI * a, x) - k2II / kII**2 * mCL[jj] * psiex(
l, kII * a, x) + k2I / kI**2 * mDL[jj] * xiex(l, kI * a, x)) # magnetic parameter2 left
eLambda3L[l] = lambda x, l=l, jj=jj: (eBL[jj] * (psiex(l, kI * a, x) + psi2ex(l, kI * a, x)) * k1I - eCL[jj] * (psiex(l, kII * a, x) + psi2ex(l, kII * a, x)) * k1II
+ eDL[jj] * (xiex(l, kI * a, x) + xi2ex(l, kI * a, x)) * k1I) # electric parameter3 left
mLambda3L[l] = lambda x, l=l, jj=jj: (mBL[jj] * (psiex(l, kI * a, x) + psi2ex(l, kI * a, x)) * MuI - mCL[jj] * (psiex(l, kII * a, x) + psi2ex(l, kII * a, x)) * MuII
+ mDL[jj] * (xiex(l, kI * a, x) + xi2ex(l, kI * a, x)) * MuI) # magnetic parameter3 left
# right
eLambda1R[l] = lambda x, l=l, jj=jj: (eBR[jj] / kI * psi1ex(l, kI * a, x) - eCR[jj] / kII * psi1ex(
l, kII * a, x) + eDR[jj] / kI * xi1ex(l, kI * a, x)) # electric parameter1 right
mLambda1R[l] = lambda x, l=l, jj=jj: (mBR[jj] / kI * psi1ex(l, kI * a, x) - mCR[jj] / kII * psi1ex(
l, kII * a, x) + mDR[jj] / kI * xi1ex(l, kI * a, x)) # magnetic parameter1 right
eLambda2R[l] = lambda x, l=l, jj=jj: (k1I / kI**2 * eBR[jj] * psiex(l, kI * a, x) - k1II / kII**2 * eCR[jj] * psiex(
l, kII * a, x) + k1I / kI**2 * eDR[jj] * xiex(l, kI * a, x)) # electric parameter2 right
mLambda2R[l] = lambda x, l=l, jj=jj: (k2I / kI**2 * mBR[jj] * psiex(l, kI * a, x) - k2II / kII**2 * mCR[jj] * psiex(
l, kII * a, x) + k2I / kI**2 * mDR[jj] * xiex(l, kI * a, x)) # magnetic parameter2 right
eLambda3R[l] = lambda x, l=l, jj=jj: (eBR[jj] * (psiex(l, kI * a, x) + psi2ex(l, kI * a, x)) * k1I - eCR[jj] * (psiex(l, kII * a, x) + psi2ex(l, kII * a, x)) * k1II
+ eDR[jj] * (xiex(l, kI * a, x) + xi2ex(l, kI * a, x)) * k1I) # electric parameter3 right
mLambda3R[l] = lambda x, l=l, jj=jj: (mBR[jj] * (psiex(l, kI * a, x) + psi2ex(l, kI * a, x)) * MuI - mCR[jj] * (psiex(l, kII * a, x) + psi2ex(l, kII * a, x)) * MuII
+ mDR[jj] * (xiex(l, kI * a, x) + xi2ex(l, kI * a, x)) * MuI) # magnetic parameter3 right
# define points for least square fitting [similar to operating on
# equations with int_ {-1}**(+1} dx legendrePl(m,x)...]
x = {}
for m in range(1, int(mmax / 3) + 1):
x[int(m)] = (-1 + 0.1 * (m - 1) / (mmax / 3))
x[int(m + mmax / 3)] = (-0.9 + 1.8 * m / (mmax / 3))
x[int(m + 2 * mmax / 3)] = (+0.9 + 0.1 * m / (mmax / 3))
efun1aL = np.zeros((mmax, lmax), dtype=complex)
efun1bL = np.zeros((mmax, lmax), dtype=complex)
efun2aL = np.zeros((mmax, lmax), dtype=complex)
efun2bL = np.zeros((mmax, lmax), dtype=complex)
efun3L = np.zeros((mmax, lmax), dtype=complex)
mfun1aL = np.zeros((mmax, lmax), dtype=complex)
mfun1bL = np.zeros((mmax, lmax), dtype=complex)
mfun2aL = np.zeros((mmax, lmax), dtype=complex)
mfun2bL = np.zeros((mmax, lmax), dtype=complex)
mfun3L = np.zeros((mmax, lmax), dtype=complex)
efun1aR = np.zeros((mmax, lmax), dtype=complex)
efun1bR = np.zeros((mmax, lmax), dtype=complex)
efun2aR = np.zeros((mmax, lmax), dtype=complex)
efun2bR = np.zeros((mmax, lmax), dtype=complex)
efun3R = np.zeros((mmax, lmax), dtype=complex)
mfun1aR = np.zeros((mmax, lmax), dtype=complex)
mfun1bR = np.zeros((mmax, lmax), dtype=complex)
mfun2aR = np.zeros((mmax, lmax), dtype=complex)
mfun2bR = np.zeros((mmax, lmax), dtype=complex)
mfun3R = np.zeros((mmax, lmax), dtype=complex)
for ii in range(mmax):
m = ii + 1
# define points for least square fitting [similar to operating on equations with int_ {-1}**(+1} dx legendrePl(m,x)...]
# the points x[m] lie in the range [-1,1] or similarly th(m) is in the range [0,pi]
# x[m] = (-1 + 2*(m-1)/(mmax-1))
for jj in range(lmax):
l = jj + 1
# left
efun1aL[ii, jj] = eLambda1L[l](x[m]) * legendrePl(l, x[m])
efun1bL[ii, jj] = eLambda1L[l](x[m]) * legendrePlex1(l, x[m])
mfun1aL[ii, jj] = mLambda1L[l](x[m]) * legendrePl(l, x[m])
mfun1bL[ii, jj] = mLambda1L[l](x[m]) * legendrePlex1(l, x[m])
mfun2aL[ii, jj] = mLambda2L[l](x[m]) * legendrePl(l, x[m])
mfun2bL[ii, jj] = mLambda2L[l](x[m]) * legendrePlex1(l, x[m])
efun2aL[ii, jj] = eLambda2L[l](x[m]) * legendrePl(l, x[m])
efun2bL[ii, jj] = eLambda2L[l](x[m]) * legendrePlex1(l, x[m])
efun3L[ii, jj] = eLambda3L[l](x[m]) * legendrePl(l, x[m])
mfun3L[ii, jj] = mLambda3L[l](x[m]) * legendrePl(l, x[m])
# right
efun1aR[ii, jj] = eLambda1R[l](x[m]) * legendrePl(l, x[m])
efun1bR[ii, jj] = eLambda1R[l](x[m]) * legendrePlex1(l, x[m])
mfun1aR[ii, jj] = mLambda1R[l](x[m]) * legendrePl(l, x[m])
mfun1bR[ii, jj] = mLambda1R[l](x[m]) * legendrePlex1(l, x[m])
mfun2aR[ii, jj] = mLambda2R[l](x[m]) * legendrePl(l, x[m])
mfun2bR[ii, jj] = mLambda2R[l](x[m]) * legendrePlex1(l, x[m])
efun2aR[ii, jj] = eLambda2R[l](x[m]) * legendrePl(l, x[m])
efun2bR[ii, jj] = eLambda2R[l](x[m]) * legendrePlex1(l, x[m])
efun3R[ii, jj] = eLambda3R[l](x[m]) * legendrePl(l, x[m])
mfun3R[ii, jj] = mLambda3R[l](x[m]) * legendrePl(l, x[m])
# first order BC can be written in form: M11*eCc + M12*eDd + M13*mCc + M14*mDd = N1
# ...
# M61*eCc + M62*eDd + M63*mCc + M64*mDd = N6
# M11...M66 are matrices including the pre-factors for the individual terms in the sums of eCc[jj], etc
# eCc...mDd are vectors including all the expansion coefficients (eCc = [eCc(1),...,eCc(lmax)])
# N1...N6 include the 0-th order terms and angular-dependent Legendre and
# Bessel functions
N1L = np.zeros(mmax, dtype=complex)
N2L = np.zeros(mmax, dtype=complex)
N3L = np.zeros(mmax, dtype=complex)
N4L = np.zeros(mmax, dtype=complex)
N5L = np.zeros(mmax, dtype=complex)
N6L = np.zeros(mmax, dtype=complex)
N1R = np.zeros(mmax, dtype=complex)
N2R = np.zeros(mmax, dtype=complex)
N3R = np.zeros(mmax, dtype=complex)
N4R = np.zeros(mmax, dtype=complex)
N5R = np.zeros(mmax, dtype=complex)
N6R = np.zeros(mmax, dtype=complex)
for ii in range(mmax):
# left
N1L[ii] = np.sum(efun1aL[ii, :]) - np.sum(mfun2bL[ii, :])
N2L[ii] = np.sum(efun1bL[ii, :]) - np.sum(mfun2aL[ii, :])
N3L[ii] = np.sum(mfun1aL[ii, :]) - np.sum(efun2bL[ii, :])
N4L[ii] = np.sum(mfun1bL[ii, :]) - np.sum(efun2aL[ii, :])
N5L[ii] = np.sum(efun3L[ii, :])
N6L[ii] = np.sum(mfun3L[ii, :])
# right
N1R[ii] = np.sum(efun1aR[ii, :]) - np.sum(mfun2bR[ii, :])
N2R[ii] = np.sum(efun1bR[ii, :]) - np.sum(mfun2aR[ii, :])
N3R[ii] = np.sum(mfun1aR[ii, :]) - np.sum(efun2bR[ii, :])
N4R[ii] = np.sum(mfun1bR[ii, :]) - np.sum(efun2aR[ii, :])
N5R[ii] = np.sum(efun3R[ii, :])
N6R[ii] = np.sum(mfun3R[ii, :])
##
M11 = np.zeros((mmax, lmax), dtype=complex)
M12 = np.zeros((mmax, lmax), dtype=complex)
M13 = np.zeros((mmax, lmax), dtype=complex)
M14 = np.zeros((mmax, lmax), dtype=complex)
M21 = np.zeros((mmax, lmax), dtype=complex)
M22 = np.zeros((mmax, lmax), dtype=complex)
M23 = np.zeros((mmax, lmax), dtype=complex)
M24 = np.zeros((mmax, lmax), dtype=complex)
M31 = np.zeros((mmax, lmax), dtype=complex)
M32 = np.zeros((mmax, lmax), dtype=complex)
M33 = np.zeros((mmax, lmax), dtype=complex)
M34 = np.zeros((mmax, lmax), dtype=complex)
M41 = np.zeros((mmax, lmax), dtype=complex)
M42 = np.zeros((mmax, lmax), dtype=complex)
M43 = np.zeros((mmax, lmax), dtype=complex)
M44 = np.zeros((mmax, lmax), dtype=complex)
M51 = np.zeros((mmax, lmax), dtype=complex)
M52 = np.zeros((mmax, lmax), dtype=complex)
M53 = np.zeros((mmax, lmax), dtype=complex)
M54 = np.zeros((mmax, lmax), dtype=complex)
M61 = np.zeros((mmax, lmax), dtype=complex)
M62 = np.zeros((mmax, lmax), dtype=complex)
M63 = np.zeros((mmax, lmax), dtype=complex)
M64 = np.zeros((mmax, lmax), dtype=complex)
for ii in range(mmax):
m = ii + 1
for jj in range(lmax):
l = jj + 1
M11[ii, jj] = +1 / kII * \
psi1xx(l, kII * a, x[m]) * legendrePl(l, x[m])
M12[ii, jj] = -1 / kI * \
xi1xx(l, kI * a, x[m]) * legendrePl(l, x[m])
M13[ii, jj] = +1 / k1II * \
psixx(l, kII * a, x[m]) * legendrePlex1(l, x[m])
M14[ii, jj] = -1 / k1I * \
xixx(l, kI * a, x[m]) * legendrePlex1(l, x[m])
M21[ii, jj] = +1 / kII * \
psi1xx(l, kII * a, x[m]) * legendrePlex1(l, x[m])
M22[ii, jj] = -1 / kI * \
xi1xx(l, kI * a, x[m]) * legendrePlex1(l, x[m])
M23[ii, jj] = +1 / k1II * \
psixx(l, kII * a, x[m]) * legendrePl(l, x[m])
M24[ii, jj] = -1 / k1I * \
xixx(l, kI * a, x[m]) * legendrePl(l, x[m])
M31[ii, jj] = +1 / k2II * \
psixx(l, kII * a, x[m]) * legendrePlex1(l, x[m])
M32[ii, jj] = -1 / k2I * \
xixx(l, kI * a, x[m]) * legendrePlex1(l, x[m])
M33[ii, jj] = M11[ii, jj]
M34[ii, jj] = M12[ii, jj]
M41[ii, jj] = +1 / k2II * \
psixx(l, kII * a, x[m]) * legendrePl(l, x[m])
M42[ii, jj] = -1 / k2I * \
xixx(l, kI * a, x[m]) * legendrePl(l, x[m])
M43[ii, jj] = M21[ii, jj]
M44[ii, jj] = M22[ii, jj]
M51[ii, jj] = + k1II * (psixx(l, kII * a, x[m]) +
psi2xx(l, kII * a, x[m])) * legendrePl(l, x[m])
M52[ii, jj] = - k1I * (xixx(l, kI * a, x[m]) +
xi2xx(l, kI * a, x[m])) * legendrePl(l, x[m])
M53[ii, jj] = 0
M54[ii, jj] = 0
M61[ii, jj] = 0
M62[ii, jj] = 0
M63[ii, jj] = + MuII * (psixx(l, kII * a, x[m]) +
psi2xx(l, kII * a, x[m])) * legendrePl(l, x[m])
M64[ii, jj] = - MuI * (xixx(l, kI * a, x[m]) +
xi2xx(l, kI * a, x[m])) * legendrePl(l, x[m])
Matrix = np.zeros((6 * mmax, 6 * lmax), dtype=complex)
for ii in range(mmax):
m = ii + 1
for jj in range(lmax):
l = jj + 1
Matrix[ii, jj] = M11[ii, jj]
Matrix[ii, jj + lmax] = M12[ii, jj]
Matrix[ii, jj + 2 * lmax] = M13[ii, jj]
Matrix[ii, jj + 3 * lmax] = M14[ii, jj]
Matrix[ii + mmax, jj] = M21[ii, jj]
Matrix[ii + mmax, jj + lmax] = M22[ii, jj]
Matrix[ii + mmax, jj + 2 * lmax] = M23[ii, jj]
Matrix[ii + mmax, jj + 3 * lmax] = M24[ii, jj]
Matrix[ii + 2 * mmax, jj] = M31[ii, jj]
Matrix[ii + 2 * mmax, jj + lmax] = M32[ii, jj]
Matrix[ii + 2 * mmax, jj + 2 * lmax] = M33[ii, jj]
Matrix[ii + 2 * mmax, jj + 3 * lmax] = M34[ii, jj]
Matrix[ii + 3 * mmax, jj] = M41[ii, jj]
Matrix[ii + 3 * mmax, jj + lmax] = M42[ii, jj]
Matrix[ii + 3 * mmax, jj + 2 * lmax] = M43[ii, jj]
Matrix[ii + 3 * mmax, jj + 3 * lmax] = M44[ii, jj]
Matrix[ii + 4 * mmax, jj] = M51[ii, jj]
Matrix[ii + 4 * mmax, jj + lmax] = M52[ii, jj]
Matrix[ii + 4 * mmax, jj + 2 * lmax] = M53[ii, jj]
Matrix[ii + 4 * mmax, jj + 3 * lmax] = M54[ii, jj]
Matrix[ii + 5 * mmax, jj] = M61[ii, jj]
Matrix[ii + 5 * mmax, jj + lmax] = M62[ii, jj]
Matrix[ii + 5 * mmax, jj + 2 * lmax] = M63[ii, jj]
Matrix[ii + 5 * mmax, jj + 3 * lmax] = M64[ii, jj]
VectorL = np.zeros(6 * mmax, dtype=complex)
VectorR = np.zeros(6 * mmax, dtype=complex)
for ii in range(mmax):
# left and right
VectorL[ii] = N1L[ii]
VectorL[ii + mmax] = N2L[ii]
VectorL[ii + 2 * mmax] = N3L[ii]
VectorL[ii + 3 * mmax] = N4L[ii]
VectorL[ii + 4 * mmax] = N5L[ii]
VectorL[ii + 5 * mmax] = N6L[ii]
VectorR[ii] = N1R[ii]
VectorR[ii + mmax] = N2R[ii]
VectorR[ii + 2 * mmax] = N3R[ii]
VectorR[ii + 3 * mmax] = N4R[ii]
VectorR[ii + 4 * mmax] = N5R[ii]
VectorR[ii + 5 * mmax] = N6R[ii]
Weight = np.zeros(6 * mmax, dtype=complex)
# weights
for ii in range(mmax):
m = ii + 1
Weight[ii] = 1
Weight[ii + mmax] = 1
Weight[ii + 2 * mmax] = 1
Weight[ii + 3 * mmax] = 1
Weight[ii + 4 * mmax] = 1 / (100 * m**0.8)
Weight[ii + 5 * mmax] = 1 / (22 * m**0.20)
# Weights were commented out?:
# xL = lscov (Matrix,VectorL.',Weight.').'
xL = lscov(np.matrix(Matrix), np.matrix(VectorL).T, np.matrix(Weight).T)
xR = lscov(np.matrix(Matrix), np.matrix(VectorR).T, np.matrix(Weight).T)
if verbose:
print('If Eps=0, ignore previous error messages regarding rank deficiency!')
eCcL = np.zeros(lmax, dtype=complex)
eDdL = np.zeros(lmax, dtype=complex)
mCcL = np.zeros(lmax, dtype=complex)
mDdL = np.zeros(lmax, dtype=complex)
eCcR = np.zeros(lmax, dtype=complex)
eDdR = np.zeros(lmax, dtype=complex)
mCcR = np.zeros(lmax, dtype=complex)
mDdR = np.zeros(lmax, dtype=complex)
for jj in range(lmax):
l = jj + 1
# left and right first order expansion coefficients
eCcL[jj] = xL[jj]
eDdL[jj] = xL[jj + lmax]
mCcL[jj] = xL[jj + 2 * lmax]
mDdL[jj] = xL[jj + 3 * lmax]
eCcR[jj] = xR[jj] * (-1)**(l - 1)
eDdR[jj] = xR[jj + lmax] * (-1)**(l - 1)
mCcR[jj] = xR[jj + 2 * lmax] * (-1)**l
mDdR[jj] = xR[jj + 3 * lmax] * (-1)**l
# corrected expansion coefficients for right optical fiber*
eBR[jj] = eBR[jj] * (-1)**(l - 1)
eCR[jj] = eCR[jj] * (-1)**(l - 1)
eDR[jj] = eDR[jj] * (-1)**(l - 1)
mBR[jj] = mBR[jj] * (-1)**l
mCR[jj] = mCR[jj] * (-1)**l
mDR[jj] = mDR[jj] * (-1)**l
#*additional factors [(-1)**(l-1),(-1)**l] should have been included when eB, mB are calculated however due to least squares approach slight antisymmetry arises...
# hence include factors only after least square approach
# Field Expansions For Incident Fields
# sums of radial, zenithal and azimuthal field terms from l=1 to l=lmax
# initialisation of sums
# Paul: These are not used?!
#EriL = lambda r, th, phi: 0
#EriR = lambda r, th, phi: 0
#EthiL = lambda r, th, phi: 0
#EthiR = lambda r, th, phi: 0
#EphiiL = lambda r, th, phi: 0
#EphiiR = lambda r, th, phi: 0
#HriL = lambda r, th, phi: 0
#HriR = lambda r, th, phi: 0
#HthiL = lambda r, th, phi: 0
#HthiR = lambda r, th, phi: 0
#HphiiL = lambda r, th, phi: 0
#HphiiR = lambda r, th, phi: 0
# left
#EriL = lambda r, th, phi: EriL(r,th,phi) + eBL[l-1]*(psi(l,kI*r)+psi2(l,kI*r))*legendrePl(l,np.cos(th))*np.cos(phi)
# EthiL = lambda r, th, phi: EthiL(r,th,phi) -(eBL[l-1]*psi1 (l,kI*r)/(kI*r)*legendrePl1(l,np.cos(th))*np.sin(th) \
# + mBL[l-1]*psi (l,kI*r)/(k1I*r)*legendrePl (l,np.cos(th))/np.sin(th))*np.cos(phi)
# EphiiL = lambda r, th, phi: EphiiL(r,th,phi) -(eBL[l-1]*psi1 (l,kI*r)/(kI*r)*legendrePl (l,np.cos(th))/np.sin(th) \
# + mBL[l-1]*psi (l,kI*r)/(k1I*r)*legendrePl1(l,np.cos(th))*np.sin(th))*np.sin(phi)
#HriL = lambda r, th, phi: HriL(r,th,phi) + mBL[l-1]*(psi(l,kI*r)+psi2(l,kI*r))*legendrePl(l,np.cos(th))*np.sin(phi)
# HthiL = lambda r, th, phi: HthiL(r,th,phi) -(mBL[l-1]*psi1 (l,kI*r)/(kI*r)*legendrePl1(l,np.cos(th))*np.sin(th) \
# + eBL[l-1]*psi (l,kI*r)/(k2I*r)*legendrePl (l,np.cos(th))/np.sin(th))*np.sin(phi)
# HphiiL = lambda r, th, phi: HphiiL(r,th,phi) +(mBL[l-1]*psi1 (l,kI*r)/(kI*r)*legendrePl (l,np.cos(th))/np.sin(th) \
# + eBL[l-1]*psi (l,kI*r)/(k2I*r)*legendrePl1(l,np.cos(th))*np.sin(th))*np.cos(phi)
# right
#EriR = lambda r, th, phi: EriR(r,th,phi) + eBR[l-1]*(psi(l,kI*r)+psi2(l,kI*r))*legendrePl(l,np.cos(th))*np.cos(phi)
# EthiR = lambda r, th, phi: EthiR(r,th,phi) -(eBR[l-1]*psi1 (l,kI*r)/(kI*r)*legendrePl1(l,np.cos(th))*np.sin(th) \
# + mBR[l-1]*psi (l,kI*r)/(k1I*r)*legendrePl (l,np.cos(th))/np.sin(th))*np.cos(phi)
# EphiiR = lambda r, th, phi: EphiiR(r,th,phi) -(eBR[l-1]*psi1 (l,kI*r)/(kI*r)*legendrePl (l,np.cos(th))/np.sin(th) \
# + mBR[l-1]*psi (l,kI*r)/(k1I*r)*legendrePl1(l,np.cos(th))*np.sin(th))*np.sin(phi)
#HriR = lambda r, th, phi: HriR(r,th,phi) + mBR[l-1]*(psi(l,kI*r)+psi2(l,kI*r))*legendrePl(l,np.cos(th))*np.sin(phi)
# HthiR = lambda r, th, phi: HthiR(r,th,phi) -(mBR[l-1]*psi1 (l,kI*r)/(kI*r)*legendrePl1(l,np.cos(th))*np.sin(th) \
# + eBR[l-1]*psi (l,kI*r)/(k2I*r)*legendrePl (l,np.cos(th))/np.sin(th))*np.sin(phi)
# HphiiR = lambda r, th, phi: HphiiR(r,th,phi) +(mBR[l-1]*psi1 (l,kI*r)/(kI*r)*legendrePl (l,np.cos(th))/np.sin(th) \
# + eBR[l-1]*psi (l,kI*r)/(k2I*r)*legendrePl1(l,np.cos(th))*np.sin(th))*np.cos(phi)
# Field Expansions For Scattered Fields
# sums of radial, zenithal and azimuthal field terms from l=1 to l=lmax
# Paul: workaround to recursing into lambda functions
def wrapper_expansion(lambda_func):
def wrapped(r, th, ph):
result = 0
for jj in range(lmax):
l = jj + 1
result += lambda_func(r, th, ph, l)
return result
return wrapped
# left
def ErsL_it(r, th, phi, l): return (
eDL[l - 1] + eDdL[l - 1]) * (xi(l, kI * r) + xi2(l, kI * r)) * legendrePl(l, np.cos(th)) * np.cos(phi)
def EthsL_it(r, th, phi, l): return -((eDL[l - 1] + eDdL[l - 1]) * xi1(l, kI * r) / (kI * r) * legendrePl1(l, np.cos(th)) * np.sin(th)
+ (mDL[l - 1] + mDdL[l - 1]) * xi(l, kI * r) / (k1I * r) * legendrePl(l, np.cos(th)) / np.sin(th)) * np.cos(phi)
def EphisL_it(r, th, phi, l): return -((eDL[l - 1] + eDdL[l - 1]) * xi1(l, kI * r) / (kI * r) * legendrePl(l, np.cos(th)) / np.sin(th)
+ (mDL[l - 1] + mDdL[l - 1]) * xi(l, kI * r) / (k1I * r) * legendrePl1(l, np.cos(th)) * np.sin(th)) * np.sin(phi)
def HrsL_it(r, th, phi, l): return +(mDL[l - 1] + mDdL[l - 1]) * (
xi(l, kI * r) + xi2(l, kI * r)) * legendrePl(l, np.cos(th)) * np.sin(phi)
def HthsL_it(r, th, phi, l): return -((mDL[l - 1] + mDdL[l - 1]) * xi1(l, kI * r) / (kI * r) * legendrePl1(l, np.cos(th)) * np.sin(th)
+ (eDL[l - 1] + eDdL[l - 1]) * xi(l, kI * r) / (k2I * r) * legendrePl(l, np.cos(th)) / np.sin(th)) * np.sin(phi)
def HphisL_it(r, th, phi, l): return +((mDL[l - 1] + mDdL[l - 1]) * xi1(l, kI * r) / (kI * r) * legendrePl(l, np.cos(th)) / np.sin(th)
+ (eDL[l - 1] + eDdL[l - 1]) * xi(l, kI * r) / (k2I * r) * legendrePl1(l, np.cos(th)) * np.sin(th)) * np.cos(phi)
# right
def ErsR_it(r, th, phi, l): return +(eDR[l - 1] + eDdR[l - 1]) * (
xi(l, kI * r) + xi2(l, kI * r)) * legendrePl(l, np.cos(th)) * np.cos(phi)
def EthsR_it(r, th, phi, l): return -((eDR[l - 1] + eDdR[l - 1]) * xi1(l, kI * r) / (kI * r) * legendrePl1(l, np.cos(th)) * np.sin(th)
+ (mDR[l - 1] + mDdR[l - 1]) * xi(l, kI * r) / (k1I * r) * legendrePl(l, np.cos(th)) / np.sin(th)) * np.cos(phi)
def EphisR_it(r, th, phi, l): return -((eDR[l - 1] + eDdR[l - 1]) * xi1(l, kI * r) / (kI * r) * legendrePl(l, np.cos(th)) / np.sin(th)
+ (mDR[l - 1] + mDdR[l - 1]) * xi(l, kI * r) / (k1I * r) * legendrePl1(l, np.cos(th)) * np.sin(th)) * np.sin(phi)
def HrsR_it(r, th, phi, l): return +(mDR[l - 1] + mDdR[l - 1]) * (
xi(l, kI * r) + xi2(l, kI * r)) * legendrePl(l, np.cos(th)) * np.sin(phi)
def HthsR_it(r, th, phi, l): return -((mDR[l - 1] + mDdR[l - 1]) * xi1(l, kI * r) / (kI * r) * legendrePl1(l, np.cos(th)) * np.sin(th)
+ (eDR[l - 1] + eDdR[l - 1]) * xi(l, kI * r) / (k2I * r) * legendrePl(l, np.cos(th)) / np.sin(th)) * np.sin(phi)
def HphisR_it(r, th, phi, l): return +((mDR[l - 1] + mDdR[l - 1]) * xi1(l, kI * r) / (kI * r) * legendrePl(l, np.cos(th)) / np.sin(th)
+ (eDR[l - 1] + eDdR[l - 1]) * xi(l, kI * r) / (k2I * r) * legendrePl1(l, np.cos(th)) * np.sin(th)) * np.cos(phi)
ErsL = wrapper_expansion(ErsL_it)
EthsL = wrapper_expansion(EthsL_it)
EphisL = wrapper_expansion(EphisL_it)
HrsL = wrapper_expansion(HrsL_it)
HthsL = wrapper_expansion(HthsL_it)
HphisL = wrapper_expansion(HphisL_it)
ErsR = wrapper_expansion(ErsR_it)
EthsR = wrapper_expansion(EthsR_it)
EphisR = wrapper_expansion(EphisR_it)
HrsR = wrapper_expansion(HrsR_it)
HthsR = wrapper_expansion(HthsR_it)
HphisR = wrapper_expansion(HphisR_it)
# Field Expansions For Internal Fields (Within)
# sums of radial, zenithal and azimuthal field terms from l=1 to l=lmax
# initialisation of sums
def ErwL_it(r, th, phi, l): return +(eCL[l - 1] + eCcL[l - 1]) * (
psi(l, kII * r) + psi2(l, kII * r)) * legendrePl(l, np.cos(th)) * np.cos(phi)
def EthwL_it(r, th, phi, l): return -((eCL[l - 1] + eCcL[l - 1]) * psi1(l, kII * r) / (kII * r) * legendrePl1(l, np.cos(th)) * np.sin(th)
+ (mCL[l - 1] + mCcL[l - 1]) * psi(l, kII * r) / (k1II * r) * legendrePl(l, np.cos(th)) / np.sin(th)) * np.cos(phi)
def EphiwL_it(r, th, phi, l): return -((eCL[l - 1] + eCcL[l - 1]) * psi1(l, kII * r) / (kII * r) * legendrePl(l, np.cos(th)) / np.sin(th)
+ (mCL[l - 1] + mCcL[l - 1]) * psi(l, kII * r) / (k1II * r) * legendrePl1(l, np.cos(th)) * np.sin(th)) * np.sin(phi)
# Paul: These are not used?!
#HrwL_it = lambda r,th,phi,l: +(mCL[l-1] + mCcL[l-1])*(psi(l,kII*r)+psi2(l,kII*r))*legendrePl(l,np.cos(th))*np.sin(phi)
# HthwL_it = lambda r,th,phi,l: -((mCL[l-1] + mCcL[l-1])*psi1 (l,kII*r)/(kII*r)*legendrePl1(l,np.cos(th))*np.sin(th) \
# + (eCL[l-1] + eCcL[l-1])*psi (l,kII*r)/(k2II*r)*legendrePl (l,np.cos(th))/np.sin(th))*np.sin(phi)
# HphiwL_it = lambda r,th,phi,l: +((mCL[l-1] + mCcL[l-1])*psi1 (l,kII*r)/(kII*r)*legendrePl (l,np.cos(th))/np.sin(th) \
# + (eCL[l-1] + eCcL[l-1])*psi (l,kII*r)/(k2II*r)*legendrePl1(l,np.cos(th))*np.sin(th))*np.cos(phi)
# right
#ErwR_it = lambda r,th,phi,l: +(eCR[l-1] + eCcR[l-1])*(psi(l,kII*r)+psi2(l,kII*r))*legendrePl(l,np.cos(th))*np.cos(phi)
# EthwR_it = lambda r,th,phi,l: -((eCR[l-1] + eCcR[l-1])*psi1 (l,kII*r)/(kII*r)*legendrePl1(l,np.cos(th))*np.sin(th) \
# + (mCR[l-1] + mCcR[l-1])*psi (l,kII*r)/(k1II*r)*legendrePl (l,np.cos(th))/np.sin(th))*np.cos(phi)
# EphiwR_it = lambda r,th,phi,l: -((eCR[l-1] + eCcR[l-1])*psi1 (l,kII*r)/(kII*r)*legendrePl (l,np.cos(th))/np.sin(th) \
# + (mCR[l-1] + mCcR[l-1])*psi (l,kII*r)/(k1II*r)*legendrePl1(l,np.cos(th))*np.sin(th))*np.sin(phi)
#HrwR_it = lambda r,th,phi,l: +(mCR[l-1] + mCcR[l-1])*(psi(l,kII*r)+psi2(l,kII*r))*legendrePl(l,np.cos(th))*np.sin(phi)
# HthwR_it = lambda r,th,phi,l: -((mCR[l-1] + mCcR[l-1])*psi1 (l,kII*r)/(kII*r)*legendrePl1(l,np.cos(th))*np.sin(th) \
# + (eCR[l-1] + eCcR[l-1])*psi (l,kII*r)/(k2II*r)*legendrePl (l,np.cos(th))/np.sin(th))*np.sin(phi)
# HphiwR_it = lambda r,th,phi,l: +((mCR[l-1] + mCcR[l-1])*psi1 (l,kII*r)/(kII*r)*legendrePl (l,np.cos(th))/np.sin(th) \
# + (eCR[l-1] + eCcR[l-1])*psi (l,kII*r)/(k2II*r)*legendrePl1(l,np.cos(th))*np.sin(th))*np.cos(phi)
ErwL = wrapper_expansion(ErwL_it)
EthwL = wrapper_expansion(EthwL_it)
EphiwL = wrapper_expansion(EphiwL_it)
# Paul: These are not used?!
#HrwL = wrapper_expansion(HrwL_it)
#HthwL = wrapper_expansion(HthwL_it)
#HphiwL = wrapper_expansion(HphiwL_it)
#
#ErwR = wrapper_expansion(ErwR_it)
#EthwR = wrapper_expansion(EthwR_it)
#EphiwR = wrapper_expansion(EphiwR_it)
#HrwR = wrapper_expansion(HrwR_it)
#HthwR = wrapper_expansion(HthwR_it)
#HphiwR = wrapper_expansion(HphiwR_it)
# Overall Fields
# sum of the electric fields in region I (independent of time) [separated
# for left and right]
def ErL(r, th, phi): return ErsL(r, th, phi) + eEriL(r, th, phi)
def ErR(r, th, phi): return ErsR(r, th, phi) + eEriR(r, th, phi)
def EthL(r, th, phi): return EthsL(r, th, phi) + eEthiL(r, th, phi)
def EthR(r, th, phi): return EthsR(r, th, phi) + eEthiR(r, th, phi)
def EphiL(r, th, phi): return EphisL(r, th, phi) + eEphiiL(r, th, phi)
def EphiR(r, th, phi): return EphisR(r, th, phi) + eEphiiR(r, th, phi)
def HrL(r, th, phi): return HrsL(r, th, phi) + eHriL(r, th, phi)
def HrR(r, th, phi): return HrsR(r, th, phi) + eHriR(r, th, phi)
def HthL(r, th, phi): return HthsL(r, th, phi) + eHthiL(r, th, phi)
def HthR(r, th, phi): return HthsR(r, th, phi) + eHthiR(r, th, phi)
def HphiL(r, th, phi): return HphisL(r, th, phi) + eHphiiL(r, th, phi)
def HphiR(r, th, phi): return HphisR(r, th, phi) + eHphiiR(r, th, phi)
# Checking Boundary Conditions
# generate points on boundary
r = np.zeros(22, dtype=complex)
th = np.zeros(22, dtype=complex)
phi = np.pi / 4
for ii in range(22):
th[ii] = np.pi / 20.003 * (ii + 1 - 0.999)
r[ii] = a * B1(np.cos(th[ii]))
FEr1 = np.zeros(22, dtype=float)
FEr2 = np.zeros(22, dtype=float)
FEth1 = np.zeros(22, dtype=float)
FEth2 = np.zeros(22, dtype=float)
FEphi1 = np.zeros(22, dtype=float)
FEphi2 = np.zeros(22, dtype=float)
FEr = np.zeros(22, dtype=float)
FEth = np.zeros(22, dtype=float)
FEphi = np.zeros(22, dtype=float)
for ii in range(22):
# internal and external fields
FEr1[ii] = np.real(EpsilonI * ErL(r[ii], th[ii], phi))
FEr2[ii] = np.real(EpsilonII * ErwL(r[ii], th[ii], phi))
FEth1[ii] = np.real(EthL(r[ii], th[ii], phi))
FEth2[ii] = np.real(EthwL(r[ii], th[ii], phi))
FEphi1[ii] = np.real(EphiL(r[ii], th[ii], phi))
FEphi2[ii] = np.real(EphiwL(r[ii], th[ii], phi))
# ratios
FEr[ii] = FEr1[ii] / FEr2[ii]
FEth[ii] = FEth1[ii] / FEth2[ii]
FEphi[ii] = FEphi1[ii] / FEphi2[ii]
if verbose:
print('Following output indicates the accuracy of the matching of the boundary conditions.')
print('theta, Er, Er_ratio, Eth, Eth_ratio, Ephi, Ephi_ratio.')
print(np.array([th, FEr1, FEr, FEth1, FEth, FEphi1, FEphi]))
print('The ratio should be close to 1 for good matching of the BCs. However if the fields tend to zero the ratio is not that relevant.')
# Data for Plotting
# time-averaged, radial stress on infinitesimal area dA = r**2 np.sin(th)
# dr dth dphi
def sigmarr(r, th, phi): return -1 / (16 * np.pi) * np.real(
EpsilonI * (ErL(r, th, phi) * np.conj(ErL(r, th, phi)) - EthL(r, th, phi) *
np.conj(EthL(r, th, phi)) - EphiL(r, th, phi) * np.conj(EphiL(r, th, phi)))
+ EpsilonI * (ErR(r, th, phi) * np.conj(ErR(r, th, phi)) - EthR(r, th, phi) *
np.conj(EthR(r, th, phi)) - EphiR(r, th, phi) * np.conj(EphiR(r, th, phi)))
+ MuI * (HrL(r, th, phi) * np.conj(HrL(r, th, phi)) - HthL(r, th, phi) *
np.conj(HthL(r, th, phi)) - HphiL(r, th, phi) * np.conj(HphiL(r, th, phi)))
+ MuI * (HrR(r, th, phi) * np.conj(HrR(r, th, phi)) - HthR(r, th, phi) * np.conj(HthR(r, th, phi)) - HphiR(r, th, phi) * np.conj(HphiR(r, th, phi))))
th = np.zeros(n_points, dtype=float)
r = np.zeros(n_points, dtype=float)
sigma = np.zeros(n_points, dtype=float)
for ii in range(n_points):
th[ii] = theta_max / (n_points - 0.998) * (ii + 1 - 0.999)
r[ii] = a * B1(np.cos(th[ii]))
sigma[ii] = sigmarr(r[ii], th[ii], np.pi / 4)
res = [th, sigma]
if ret_legendre_decomp:
coeff = stress2legendre(stress=sigma, theta=th, n_poly=lmax)
res.append(coeff)
return res