Compare commits
No commits in common. 'master' and 'fdtd' have entirely different histories.
@ -1,68 +0,0 @@
|
|||||||
import numpy, scipy, gridlock, fdfd_tools
|
|
||||||
from fdfd_tools import bloch
|
|
||||||
from numpy.linalg import norm
|
|
||||||
import logging
|
|
||||||
|
|
||||||
logging.basicConfig(level=logging.DEBUG)
|
|
||||||
logger = logging.getLogger(__name__)
|
|
||||||
|
|
||||||
|
|
||||||
dx = 40
|
|
||||||
x_period = 400
|
|
||||||
y_period = z_period = 2000
|
|
||||||
g = gridlock.Grid([numpy.arange(-x_period/2, x_period/2, dx),
|
|
||||||
numpy.arange(-1000, 1000, dx),
|
|
||||||
numpy.arange(-1000, 1000, dx)],
|
|
||||||
shifts=numpy.array([[0,0,0]]),
|
|
||||||
initial=1.445**2,
|
|
||||||
periodic=True)
|
|
||||||
|
|
||||||
g.draw_cuboid([0,0,0], [200e8, 220, 220], eps=3.47**2)
|
|
||||||
|
|
||||||
#x_period = y_period = z_period = 13000
|
|
||||||
#g = gridlock.Grid([numpy.arange(3), ]*3,
|
|
||||||
# shifts=numpy.array([[0, 0, 0]]),
|
|
||||||
# initial=2.0**2,
|
|
||||||
# periodic=True)
|
|
||||||
|
|
||||||
g2 = g.copy()
|
|
||||||
g2.shifts = numpy.zeros((6,3))
|
|
||||||
g2.grids = [numpy.zeros(g.shape) for _ in range(6)]
|
|
||||||
|
|
||||||
epsilon = [g.grids[0],] * 3
|
|
||||||
reciprocal_lattice = numpy.diag(1e6/numpy.array([x_period, y_period, z_period])) #cols are vectors
|
|
||||||
|
|
||||||
#print('Finding k at 1550nm')
|
|
||||||
#k, f = bloch.find_k(frequency=1/1550,
|
|
||||||
# tolerance=(1/1550 - 1/1551),
|
|
||||||
# direction=[1, 0, 0],
|
|
||||||
# G_matrix=reciprocal_lattice,
|
|
||||||
# epsilon=epsilon,
|
|
||||||
# band=0)
|
|
||||||
#
|
|
||||||
#print("k={}, f={}, 1/f={}, k/f={}".format(k, f, 1/f, norm(reciprocal_lattice @ k) / f ))
|
|
||||||
|
|
||||||
print('Finding f at [0.25, 0, 0]')
|
|
||||||
for k0x in [.25]:
|
|
||||||
k0 = numpy.array([k0x, 0, 0])
|
|
||||||
|
|
||||||
kmag = norm(reciprocal_lattice @ k0)
|
|
||||||
tolerance = (1e6/1550) * 1e-4/1.5 # df = f * dn_eff / n
|
|
||||||
logger.info('tolerance {}'.format(tolerance))
|
|
||||||
|
|
||||||
n, v = bloch.eigsolve(4, k0, G_matrix=reciprocal_lattice, epsilon=epsilon, tolerance=tolerance)
|
|
||||||
v2e = bloch.hmn_2_exyz(k0, G_matrix=reciprocal_lattice, epsilon=epsilon)
|
|
||||||
v2h = bloch.hmn_2_hxyz(k0, G_matrix=reciprocal_lattice, epsilon=epsilon)
|
|
||||||
ki = bloch.generate_kmn(k0, reciprocal_lattice, g.shape)
|
|
||||||
|
|
||||||
z = 0
|
|
||||||
e = v2e(v[0])
|
|
||||||
for i in range(3):
|
|
||||||
g2.grids[i] += numpy.real(e[i])
|
|
||||||
g2.grids[i+3] += numpy.imag(e[i])
|
|
||||||
|
|
||||||
f = numpy.sqrt(numpy.real(numpy.abs(n))) # TODO
|
|
||||||
print('k0x = {:3g}\n eigval = {}\n f = {}\n'.format(k0x, n, f))
|
|
||||||
n_eff = norm(reciprocal_lattice @ k0) / f
|
|
||||||
print('kmag/f = n_eff = {} \n wl = {}\n'.format(n_eff, 1/f ))
|
|
||||||
|
|
@ -1,513 +0,0 @@
|
|||||||
'''
|
|
||||||
Bloch eigenmode solver/operators
|
|
||||||
|
|
||||||
This module contains functions for generating and solving the
|
|
||||||
3D Bloch eigenproblem. The approach is to transform the problem
|
|
||||||
into the (spatial) fourier domain, transforming the equation
|
|
||||||
1/mu * curl(1/eps * curl(H)) = (w/c)^2 H
|
|
||||||
into
|
|
||||||
conv(1/mu_k, ik x conv(1/eps_k, ik x H_k)) = (w/c)^2 H_k
|
|
||||||
where:
|
|
||||||
- the _k subscript denotes a 3D fourier transformed field
|
|
||||||
- each component of H_k corresponds to a plane wave with wavevector k
|
|
||||||
- x is the cross product
|
|
||||||
- conv denotes convolution
|
|
||||||
|
|
||||||
Since k and H are orthogonal for each plane wave, we can use each
|
|
||||||
k to create an orthogonal basis (k, m, n), with k x m = n, and
|
|
||||||
|m| = |n| = 1. The cross products are then simplified with
|
|
||||||
|
|
||||||
k @ h = kx hx + ky hy + kz hz = 0 = hk
|
|
||||||
h = hk + hm + hn = hm + hn
|
|
||||||
k = kk + km + kn = kk = |k|
|
|
||||||
|
|
||||||
k x h = (ky hz - kz hy,
|
|
||||||
kz hx - kx hz,
|
|
||||||
kx hy - ky hx)
|
|
||||||
= ((k x h) @ k, (k x h) @ m, (k x h) @ n)_kmn
|
|
||||||
= (0, (m x k) @ h, (n x k) @ h)_kmn # triple product ordering
|
|
||||||
= (0, kk (-n @ h), kk (m @ h))_kmn # (m x k) = -|k| n, etc.
|
|
||||||
= |k| (0, -h @ n, h @ m)_kmn
|
|
||||||
|
|
||||||
k x h = (km hn - kn hm,
|
|
||||||
kn hk - kk hn,
|
|
||||||
kk hm - km hk)_kmn
|
|
||||||
= (0, -kk hn, kk hm)_kmn
|
|
||||||
= (-kk hn)(mx, my, mz) + (kk hm)(nx, ny, nz)
|
|
||||||
= |k| (hm * (nx, ny, nz) - hn * (mx, my, mz))
|
|
||||||
|
|
||||||
where h is shorthand for H_k, (...)_kmn deontes the (k, m, n) basis,
|
|
||||||
and e.g. hm is the component of h in the m direction.
|
|
||||||
|
|
||||||
We can also simplify conv(X_k, Y_k) as fftn(X * ifftn(Y_k)).
|
|
||||||
|
|
||||||
Using these results and storing H_k as h = (hm, hn), we have
|
|
||||||
e_xyz = fftn(1/eps * ifftn(|k| (hm * n - hn * m)))
|
|
||||||
b_mn = |k| (-e_xyz @ n, e_xyz @ m)
|
|
||||||
h_mn = fftn(1/mu * ifftn(b_m * m + b_n * n))
|
|
||||||
which forms the operator from the left side of the equation.
|
|
||||||
|
|
||||||
We can then use a preconditioned block Rayleigh iteration algorithm, as in
|
|
||||||
SG Johnson and JD Joannopoulos, Block-iterative frequency-domain methods
|
|
||||||
for Maxwell's equations in a planewave basis, Optics Express 8, 3, 173-190 (2001)
|
|
||||||
(similar to that used in MPB) to find the eigenvectors for this operator.
|
|
||||||
|
|
||||||
===
|
|
||||||
|
|
||||||
Typically you will want to do something like
|
|
||||||
|
|
||||||
recip_lattice = numpy.diag(1/numpy.array(epsilon[0].shape * dx))
|
|
||||||
n, v = bloch.eigsolve(5, k0, recip_lattice, epsilon)
|
|
||||||
f = numpy.sqrt(-numpy.real(n[0]))
|
|
||||||
n_eff = norm(recip_lattice @ k0) / f
|
|
||||||
|
|
||||||
v2e = bloch.hmn_2_exyz(k0, recip_lattice, epsilon)
|
|
||||||
e_field = v2e(v[0])
|
|
||||||
|
|
||||||
k, f = find_k(frequency=1/1550,
|
|
||||||
tolerance=(1/1550 - 1/1551),
|
|
||||||
direction=[1, 0, 0],
|
|
||||||
G_matrix=recip_lattice,
|
|
||||||
epsilon=epsilon,
|
|
||||||
band=0)
|
|
||||||
|
|
||||||
'''
|
|
||||||
|
|
||||||
from typing import List, Tuple, Callable, Dict
|
|
||||||
import logging
|
|
||||||
import numpy
|
|
||||||
from numpy.fft import fftn, ifftn, fftfreq
|
|
||||||
import scipy
|
|
||||||
import scipy.optimize
|
|
||||||
from scipy.linalg import norm
|
|
||||||
import scipy.sparse.linalg as spalg
|
|
||||||
|
|
||||||
from .eigensolvers import rayleigh_quotient_iteration
|
|
||||||
from . import field_t
|
|
||||||
|
|
||||||
logger = logging.getLogger(__name__)
|
|
||||||
|
|
||||||
|
|
||||||
def generate_kmn(k0: numpy.ndarray,
|
|
||||||
G_matrix: numpy.ndarray,
|
|
||||||
shape: numpy.ndarray
|
|
||||||
) -> Tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray]:
|
|
||||||
"""
|
|
||||||
Generate a (k, m, n) orthogonal basis for each k-vector in the simulation grid.
|
|
||||||
|
|
||||||
:param k0: [k0x, k0y, k0z], Bloch wavevector, in G basis.
|
|
||||||
:param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
|
|
||||||
:param shape: [nx, ny, nz] shape of the simulation grid.
|
|
||||||
:return: (|k|, m, n) where |k| has shape tuple(shape) + (1,)
|
|
||||||
and m, n have shape tuple(shape) + (3,).
|
|
||||||
All are given in the xyz basis (e.g. |k|[0,0,0] = norm(G_matrix @ k0)).
|
|
||||||
"""
|
|
||||||
k0 = numpy.array(k0)
|
|
||||||
|
|
||||||
Gi_grids = numpy.meshgrid(*(fftfreq(n, 1/n) for n in shape[:3]), indexing='ij')
|
|
||||||
Gi = numpy.stack(Gi_grids, axis=3)
|
|
||||||
|
|
||||||
k_G = k0[None, None, None, :] - Gi
|
|
||||||
k_xyz = numpy.rollaxis(G_matrix @ numpy.rollaxis(k_G, 3, 2), 3, 2)
|
|
||||||
|
|
||||||
m = numpy.broadcast_to([0, 1, 0], tuple(shape[:3]) + (3,)).astype(float)
|
|
||||||
n = numpy.broadcast_to([0, 0, 1], tuple(shape[:3]) + (3,)).astype(float)
|
|
||||||
|
|
||||||
xy_non0 = numpy.any(k_xyz[:, :, :, 0:1] != 0, axis=3)
|
|
||||||
if numpy.any(xy_non0):
|
|
||||||
u = numpy.cross(k_xyz[xy_non0], [0, 0, 1])
|
|
||||||
m[xy_non0, :] = u / norm(u, axis=1)[:, None]
|
|
||||||
|
|
||||||
z_non0 = numpy.any(k_xyz != 0, axis=3)
|
|
||||||
if numpy.any(z_non0):
|
|
||||||
v = numpy.cross(k_xyz[z_non0], m[z_non0])
|
|
||||||
n[z_non0, :] = v / norm(v, axis=1)[:, None]
|
|
||||||
|
|
||||||
k_mag = norm(k_xyz, axis=3)[:, :, :, None]
|
|
||||||
return k_mag, m, n
|
|
||||||
|
|
||||||
|
|
||||||
def maxwell_operator(k0: numpy.ndarray,
|
|
||||||
G_matrix: numpy.ndarray,
|
|
||||||
epsilon: field_t,
|
|
||||||
mu: field_t = None
|
|
||||||
) -> Callable[[numpy.ndarray], numpy.ndarray]:
|
|
||||||
"""
|
|
||||||
Generate the Maxwell operator
|
|
||||||
conv(1/mu_k, ik x conv(1/eps_k, ik x ___))
|
|
||||||
which is the spatial-frequency-space representation of
|
|
||||||
1/mu * curl(1/eps * curl(___))
|
|
||||||
|
|
||||||
The operator is a function that acts on a vector h_mn of size (2 * epsilon[0].size)
|
|
||||||
|
|
||||||
See the module-level docstring for more information.
|
|
||||||
|
|
||||||
:param k0: Bloch wavevector, [k0x, k0y, k0z].
|
|
||||||
:param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
|
|
||||||
:param epsilon: Dielectric constant distribution for the simulation.
|
|
||||||
All fields are sampled at cell centers (i.e., NOT Yee-gridded)
|
|
||||||
:param mu: Magnetic permability distribution for the simulation.
|
|
||||||
Default None (1 everywhere).
|
|
||||||
:return: Function which applies the maxwell operator to h_mn.
|
|
||||||
"""
|
|
||||||
|
|
||||||
shape = epsilon[0].shape + (1,)
|
|
||||||
k_mag, m, n = generate_kmn(k0, G_matrix, shape)
|
|
||||||
|
|
||||||
epsilon = numpy.stack(epsilon, 3)
|
|
||||||
if mu is not None:
|
|
||||||
mu = numpy.stack(mu, 3)
|
|
||||||
|
|
||||||
def operator(h: numpy.ndarray):
|
|
||||||
"""
|
|
||||||
Maxwell operator for Bloch eigenmode simulation.
|
|
||||||
|
|
||||||
h is complex 2-field in (m, n) basis, vectorized
|
|
||||||
|
|
||||||
:param h: Raveled h_mn; size (2 * epsilon[0].size).
|
|
||||||
:return: Raveled conv(1/mu_k, ik x conv(1/eps_k, ik x h_mn)).
|
|
||||||
"""
|
|
||||||
hin_m, hin_n = [hi.reshape(shape) for hi in numpy.split(h, 2)]
|
|
||||||
|
|
||||||
#{d,e,h}_xyz fields are complex 3-fields in (1/x, 1/y, 1/z) basis
|
|
||||||
|
|
||||||
# cross product and transform into xyz basis
|
|
||||||
d_xyz = (n * hin_m -
|
|
||||||
m * hin_n) * k_mag
|
|
||||||
|
|
||||||
# divide by epsilon
|
|
||||||
e_xyz = fftn(ifftn(d_xyz, axes=range(3)) / epsilon, axes=range(3))
|
|
||||||
|
|
||||||
# cross product and transform into mn basis
|
|
||||||
b_m = numpy.sum(e_xyz * n, axis=3)[:, :, :, None] * -k_mag
|
|
||||||
b_n = numpy.sum(e_xyz * m, axis=3)[:, :, :, None] * +k_mag
|
|
||||||
|
|
||||||
if mu is None:
|
|
||||||
h_m, h_n = b_m, b_n
|
|
||||||
else:
|
|
||||||
# transform from mn to xyz
|
|
||||||
b_xyz = (m * b_m[:, :, :, None] +
|
|
||||||
n * b_n[:, :, :, None])
|
|
||||||
|
|
||||||
# divide by mu
|
|
||||||
h_xyz = fftn(ifftn(b_xyz, axes=range(3)) / mu, axes=range(3))
|
|
||||||
|
|
||||||
# transform back to mn
|
|
||||||
h_m = numpy.sum(h_xyz * m, axis=3)
|
|
||||||
h_n = numpy.sum(h_xyz * n, axis=3)
|
|
||||||
return numpy.hstack((h_m.ravel(), h_n.ravel()))
|
|
||||||
|
|
||||||
return operator
|
|
||||||
|
|
||||||
|
|
||||||
def hmn_2_exyz(k0: numpy.ndarray,
|
|
||||||
G_matrix: numpy.ndarray,
|
|
||||||
epsilon: field_t,
|
|
||||||
) -> Callable[[numpy.ndarray], field_t]:
|
|
||||||
"""
|
|
||||||
Generate an operator which converts a vectorized spatial-frequency-space
|
|
||||||
h_mn into an E-field distribution, i.e.
|
|
||||||
ifft(conv(1/eps_k, ik x h_mn))
|
|
||||||
|
|
||||||
The operator is a function that acts on a vector h_mn of size (2 * epsilon[0].size)
|
|
||||||
|
|
||||||
See the module-level docstring for more information.
|
|
||||||
|
|
||||||
:param k0: Bloch wavevector, [k0x, k0y, k0z].
|
|
||||||
:param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
|
|
||||||
:param epsilon: Dielectric constant distribution for the simulation.
|
|
||||||
All fields are sampled at cell centers (i.e., NOT Yee-gridded)
|
|
||||||
:return: Function for converting h_mn into E_xyz
|
|
||||||
"""
|
|
||||||
shape = epsilon[0].shape + (1,)
|
|
||||||
epsilon = numpy.stack(epsilon, 3)
|
|
||||||
|
|
||||||
k_mag, m, n = generate_kmn(k0, G_matrix, shape)
|
|
||||||
|
|
||||||
def operator(h: numpy.ndarray) -> field_t:
|
|
||||||
hin_m, hin_n = [hi.reshape(shape) for hi in numpy.split(h, 2)]
|
|
||||||
d_xyz = (n * hin_m -
|
|
||||||
m * hin_n) * k_mag
|
|
||||||
|
|
||||||
# divide by epsilon
|
|
||||||
return [ei for ei in numpy.rollaxis(ifftn(d_xyz, axes=range(3)) / epsilon, 3)]
|
|
||||||
|
|
||||||
return operator
|
|
||||||
|
|
||||||
|
|
||||||
def hmn_2_hxyz(k0: numpy.ndarray,
|
|
||||||
G_matrix: numpy.ndarray,
|
|
||||||
epsilon: field_t
|
|
||||||
) -> Callable[[numpy.ndarray], field_t]:
|
|
||||||
"""
|
|
||||||
Generate an operator which converts a vectorized spatial-frequency-space
|
|
||||||
h_mn into an H-field distribution, i.e.
|
|
||||||
ifft(h_mn)
|
|
||||||
|
|
||||||
The operator is a function that acts on a vector h_mn of size (2 * epsilon[0].size)
|
|
||||||
|
|
||||||
See the module-level docstring for more information.
|
|
||||||
|
|
||||||
:param k0: Bloch wavevector, [k0x, k0y, k0z].
|
|
||||||
:param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
|
|
||||||
:param epsilon: Dielectric constant distribution for the simulation.
|
|
||||||
Only epsilon[0].shape is used.
|
|
||||||
:return: Function for converting h_mn into H_xyz
|
|
||||||
"""
|
|
||||||
shape = epsilon[0].shape + (1,)
|
|
||||||
k_mag, m, n = generate_kmn(k0, G_matrix, shape)
|
|
||||||
|
|
||||||
def operator(h: numpy.ndarray):
|
|
||||||
hin_m, hin_n = [hi.reshape(shape) for hi in numpy.split(h, 2)]
|
|
||||||
h_xyz = (m * hin_m +
|
|
||||||
n * hin_n)
|
|
||||||
return [ifftn(hi) for hi in numpy.rollaxis(h_xyz, 3)]
|
|
||||||
|
|
||||||
return operator
|
|
||||||
|
|
||||||
|
|
||||||
def inverse_maxwell_operator_approx(k0: numpy.ndarray,
|
|
||||||
G_matrix: numpy.ndarray,
|
|
||||||
epsilon: field_t,
|
|
||||||
mu: field_t = None
|
|
||||||
) -> Callable[[numpy.ndarray], numpy.ndarray]:
|
|
||||||
"""
|
|
||||||
Generate an approximate inverse of the Maxwell operator,
|
|
||||||
ik x conv(eps_k, ik x conv(mu_k, ___))
|
|
||||||
which can be used to improve the speed of ARPACK in shift-invert mode.
|
|
||||||
|
|
||||||
See the module-level docstring for more information.
|
|
||||||
|
|
||||||
:param k0: Bloch wavevector, [k0x, k0y, k0z].
|
|
||||||
:param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
|
|
||||||
:param epsilon: Dielectric constant distribution for the simulation.
|
|
||||||
All fields are sampled at cell centers (i.e., NOT Yee-gridded)
|
|
||||||
:param mu: Magnetic permability distribution for the simulation.
|
|
||||||
Default None (1 everywhere).
|
|
||||||
:return: Function which applies the approximate inverse of the maxwell operator to h_mn.
|
|
||||||
"""
|
|
||||||
shape = epsilon[0].shape + (1,)
|
|
||||||
epsilon = numpy.stack(epsilon, 3)
|
|
||||||
|
|
||||||
k_mag, m, n = generate_kmn(k0, G_matrix, shape)
|
|
||||||
|
|
||||||
if mu is not None:
|
|
||||||
mu = numpy.stack(mu, 3)
|
|
||||||
|
|
||||||
def operator(h: numpy.ndarray):
|
|
||||||
"""
|
|
||||||
Approximate inverse Maxwell operator for Bloch eigenmode simulation.
|
|
||||||
|
|
||||||
h is complex 2-field in (m, n) basis, vectorized
|
|
||||||
|
|
||||||
:param h: Raveled h_mn; size (2 * epsilon[0].size).
|
|
||||||
:return: Raveled ik x conv(eps_k, ik x conv(mu_k, h_mn))
|
|
||||||
"""
|
|
||||||
hin_m, hin_n = [hi.reshape(shape) for hi in numpy.split(h, 2)]
|
|
||||||
|
|
||||||
#{d,e,h}_xyz fields are complex 3-fields in (1/x, 1/y, 1/z) basis
|
|
||||||
|
|
||||||
if mu is None:
|
|
||||||
b_m, b_n = hin_m, hin_n
|
|
||||||
else:
|
|
||||||
# transform from mn to xyz
|
|
||||||
h_xyz = (m * hin_m[:, :, :, None] +
|
|
||||||
n * hin_n[:, :, :, None])
|
|
||||||
|
|
||||||
# multiply by mu
|
|
||||||
b_xyz = fftn(ifftn(h_xyz, axes=range(3)) * mu, axes=range(3))
|
|
||||||
|
|
||||||
# transform back to mn
|
|
||||||
b_m = numpy.sum(b_xyz * m, axis=3)
|
|
||||||
b_n = numpy.sum(b_xyz * n, axis=3)
|
|
||||||
|
|
||||||
# cross product and transform into xyz basis
|
|
||||||
e_xyz = (n * b_m -
|
|
||||||
m * b_n) / k_mag
|
|
||||||
|
|
||||||
# multiply by epsilon
|
|
||||||
d_xyz = fftn(ifftn(e_xyz, axes=range(3)) * epsilon, axes=range(3))
|
|
||||||
|
|
||||||
# cross product and transform into mn basis crossinv_t2c
|
|
||||||
h_m = numpy.sum(e_xyz * n, axis=3)[:, :, :, None] / +k_mag
|
|
||||||
h_n = numpy.sum(e_xyz * m, axis=3)[:, :, :, None] / -k_mag
|
|
||||||
|
|
||||||
return numpy.hstack((h_m.ravel(), h_n.ravel()))
|
|
||||||
|
|
||||||
return operator
|
|
||||||
|
|
||||||
|
|
||||||
def eigsolve(num_modes: int,
|
|
||||||
k0: numpy.ndarray,
|
|
||||||
G_matrix: numpy.ndarray,
|
|
||||||
epsilon: field_t,
|
|
||||||
mu: field_t = None,
|
|
||||||
tolerance = 1e-8,
|
|
||||||
) -> Tuple[numpy.ndarray, numpy.ndarray]:
|
|
||||||
"""
|
|
||||||
Find the first (lowest-frequency) num_modes eigenmodes with Bloch wavevector
|
|
||||||
k0 of the specified structure.
|
|
||||||
|
|
||||||
:param k0: Bloch wavevector, [k0x, k0y, k0z].
|
|
||||||
:param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
|
|
||||||
:param epsilon: Dielectric constant distribution for the simulation.
|
|
||||||
All fields are sampled at cell centers (i.e., NOT Yee-gridded)
|
|
||||||
:param mu: Magnetic permability distribution for the simulation.
|
|
||||||
Default None (1 everywhere).
|
|
||||||
:return: (eigenvalues, eigenvectors) where eigenvalues[i] corresponds to the
|
|
||||||
vector eigenvectors[i, :]
|
|
||||||
"""
|
|
||||||
h_size = 2 * epsilon[0].size
|
|
||||||
|
|
||||||
kmag = norm(G_matrix @ k0)
|
|
||||||
|
|
||||||
'''
|
|
||||||
Generate the operators
|
|
||||||
'''
|
|
||||||
mop = maxwell_operator(k0=k0, G_matrix=G_matrix, epsilon=epsilon, mu=mu)
|
|
||||||
imop = inverse_maxwell_operator_approx(k0=k0, G_matrix=G_matrix, epsilon=epsilon, mu=mu)
|
|
||||||
|
|
||||||
scipy_op = spalg.LinearOperator(dtype=complex, shape=(h_size, h_size), matvec=mop)
|
|
||||||
scipy_iop = spalg.LinearOperator(dtype=complex, shape=(h_size, h_size), matvec=imop)
|
|
||||||
|
|
||||||
y_shape = (h_size, num_modes)
|
|
||||||
|
|
||||||
def rayleigh_quotient(Z: numpy.ndarray, approx_grad: bool = True):
|
|
||||||
"""
|
|
||||||
Absolute value of the block Rayleigh quotient, and the associated gradient.
|
|
||||||
|
|
||||||
See Johnson and Joannopoulos, Opt. Expr. 8, 3 (2001) for details (full
|
|
||||||
citation in module docstring).
|
|
||||||
|
|
||||||
===
|
|
||||||
|
|
||||||
Notes on my understanding of the procedure:
|
|
||||||
|
|
||||||
Minimize f(Y) = |trace((Y.H @ A @ Y)|, making use of Y = Z @ inv(Z.H @ Z)^(1/2)
|
|
||||||
(a polar orthogonalization of Y). This gives f(Z) = |trace(Z.H @ A @ Z @ U)|,
|
|
||||||
where U = inv(Z.H @ Z). We minimize the absolute value to find the eigenvalues
|
|
||||||
with smallest magnitude.
|
|
||||||
|
|
||||||
The gradient is P @ (A @ Z @ U), where P = (1 - Z @ U @ Z.H) is a projection
|
|
||||||
onto the space orthonormal to Z. If approx_grad is True, the approximate
|
|
||||||
inverse of the maxwell operator is used to precondition the gradient.
|
|
||||||
"""
|
|
||||||
z = Z.view(dtype=complex).reshape(y_shape)
|
|
||||||
U = numpy.linalg.inv(z.conj().T @ z)
|
|
||||||
zU = z @ U
|
|
||||||
AzU = scipy_op @ zU
|
|
||||||
zTAzU = z.conj().T @ AzU
|
|
||||||
f = numpy.real(numpy.trace(zTAzU))
|
|
||||||
if approx_grad:
|
|
||||||
df_dy = scipy_iop @ (AzU - zU @ zTAzU)
|
|
||||||
else:
|
|
||||||
df_dy = (AzU - zU @ zTAzU)
|
|
||||||
|
|
||||||
df_dy_flat = df_dy.view(dtype=float).ravel()
|
|
||||||
return numpy.abs(f), numpy.sign(f) * df_dy_flat
|
|
||||||
|
|
||||||
'''
|
|
||||||
Use the conjugate gradient method and the approximate gradient calculation to
|
|
||||||
quickly find approximate eigenvectors.
|
|
||||||
'''
|
|
||||||
result = scipy.optimize.minimize(rayleigh_quotient,
|
|
||||||
numpy.random.rand(*y_shape, 2),
|
|
||||||
jac=True,
|
|
||||||
method='L-BFGS-B',
|
|
||||||
tol=1e-20,
|
|
||||||
options={'maxiter': 2000, 'gtol':0, 'ftol':1e-20 , 'disp':True})#, 'maxls':80, 'm':30})
|
|
||||||
|
|
||||||
|
|
||||||
result = scipy.optimize.minimize(lambda y: rayleigh_quotient(y, True),
|
|
||||||
result.x,
|
|
||||||
jac=True,
|
|
||||||
method='L-BFGS-B',
|
|
||||||
tol=1e-20,
|
|
||||||
options={'maxiter': 2000, 'gtol':0, 'disp':True})
|
|
||||||
|
|
||||||
result = scipy.optimize.minimize(lambda y: rayleigh_quotient(y, False),
|
|
||||||
result.x,
|
|
||||||
jac=True,
|
|
||||||
method='L-BFGS-B',
|
|
||||||
tol=1e-20,
|
|
||||||
options={'maxiter': 2000, 'gtol':0, 'disp':True})
|
|
||||||
|
|
||||||
for i in range(20):
|
|
||||||
result = scipy.optimize.minimize(lambda y: rayleigh_quotient(y, False),
|
|
||||||
result.x,
|
|
||||||
jac=True,
|
|
||||||
method='L-BFGS-B',
|
|
||||||
tol=1e-20,
|
|
||||||
options={'maxiter': 70, 'gtol':0, 'disp':True})
|
|
||||||
if result.nit == 0:
|
|
||||||
# We took 0 steps, so re-running won't help
|
|
||||||
break
|
|
||||||
|
|
||||||
|
|
||||||
z = result.x.view(dtype=complex).reshape(y_shape)
|
|
||||||
|
|
||||||
'''
|
|
||||||
Recover eigenvectors from Z
|
|
||||||
'''
|
|
||||||
U = numpy.linalg.inv(z.conj().T @ z)
|
|
||||||
y = z @ scipy.linalg.sqrtm(U)
|
|
||||||
w = y.conj().T @ (scipy_op @ y)
|
|
||||||
|
|
||||||
eigvals, w_eigvecs = numpy.linalg.eig(w)
|
|
||||||
eigvecs = y @ w_eigvecs
|
|
||||||
|
|
||||||
for i in range(len(eigvals)):
|
|
||||||
v = eigvecs[:, i]
|
|
||||||
n = eigvals[i]
|
|
||||||
v /= norm(v)
|
|
||||||
eigness = norm(scipy_op @ v - (v.conj() @ (scipy_op @ v)) * v )
|
|
||||||
f = numpy.sqrt(-numpy.real(n))
|
|
||||||
df = numpy.sqrt(-numpy.real(n + eigness))
|
|
||||||
neff_err = kmag * (1/df - 1/f)
|
|
||||||
logger.info('eigness {}: {}\n neff_err: {}'.format(i, eigness, neff_err))
|
|
||||||
|
|
||||||
order = numpy.argsort(numpy.abs(eigvals))
|
|
||||||
return eigvals[order], eigvecs.T[order]
|
|
||||||
|
|
||||||
|
|
||||||
def find_k(frequency: float,
|
|
||||||
tolerance: float,
|
|
||||||
direction: numpy.ndarray,
|
|
||||||
G_matrix: numpy.ndarray,
|
|
||||||
epsilon: field_t,
|
|
||||||
mu: field_t = None,
|
|
||||||
band: int = 0,
|
|
||||||
k_min: float = 0,
|
|
||||||
k_max: float = 0.5,
|
|
||||||
) -> Tuple[numpy.ndarray, float]:
|
|
||||||
"""
|
|
||||||
Search for a bloch vector that has a given frequency.
|
|
||||||
|
|
||||||
:param frequency: Target frequency.
|
|
||||||
:param tolerance: Target frequency tolerance.
|
|
||||||
:param direction: k-vector direction to search along.
|
|
||||||
:param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns.
|
|
||||||
:param epsilon: Dielectric constant distribution for the simulation.
|
|
||||||
All fields are sampled at cell centers (i.e., NOT Yee-gridded)
|
|
||||||
:param mu: Magnetic permability distribution for the simulation.
|
|
||||||
Default None (1 everywhere).
|
|
||||||
:param band: Which band to search in. Default 0 (lowest frequency).
|
|
||||||
return: (k, actual_frequency) The found k-vector and its frequency
|
|
||||||
"""
|
|
||||||
|
|
||||||
direction = numpy.array(direction) / norm(direction)
|
|
||||||
|
|
||||||
def get_f(k0_mag: float, band: int = 0):
|
|
||||||
k0 = direction * k0_mag
|
|
||||||
n, _v = eigsolve(band + 1, k0, G_matrix=G_matrix, epsilon=epsilon)
|
|
||||||
f = numpy.sqrt(numpy.abs(numpy.real(n[band])))
|
|
||||||
return f
|
|
||||||
|
|
||||||
res = scipy.optimize.minimize_scalar(lambda x: abs(get_f(x, band) - frequency),
|
|
||||||
(k_min + k_max) / 2,
|
|
||||||
method='Bounded',
|
|
||||||
bounds=(k_min, k_max),
|
|
||||||
options={'xatol': abs(tolerance)})
|
|
||||||
return res.x * direction, res.fun + frequency
|
|
||||||
|
|
||||||
|
|
@ -1,120 +0,0 @@
|
|||||||
"""
|
|
||||||
Solvers for eigenvalue / eigenvector problems
|
|
||||||
"""
|
|
||||||
from typing import Tuple, List
|
|
||||||
import numpy
|
|
||||||
from numpy.linalg import norm
|
|
||||||
from scipy import sparse
|
|
||||||
import scipy.sparse.linalg as spalg
|
|
||||||
|
|
||||||
|
|
||||||
def power_iteration(operator: sparse.spmatrix,
|
|
||||||
guess_vector: numpy.ndarray = None,
|
|
||||||
iterations: int = 20,
|
|
||||||
) -> Tuple[complex, numpy.ndarray]:
|
|
||||||
"""
|
|
||||||
Use power iteration to estimate the dominant eigenvector of a matrix.
|
|
||||||
|
|
||||||
:param operator: Matrix to analyze.
|
|
||||||
:param guess_vector: Starting point for the eigenvector. Default is a randomly chosen vector.
|
|
||||||
:param iterations: Number of iterations to perform. Default 20.
|
|
||||||
:return: (Largest-magnitude eigenvalue, Corresponding eigenvector estimate)
|
|
||||||
"""
|
|
||||||
if numpy.any(numpy.equal(guess_vector, None)):
|
|
||||||
v = numpy.random.rand(operator.shape[0])
|
|
||||||
else:
|
|
||||||
v = guess_vector
|
|
||||||
|
|
||||||
for _ in range(iterations):
|
|
||||||
v = operator @ v
|
|
||||||
v /= norm(v)
|
|
||||||
|
|
||||||
lm_eigval = v.conj() @ (operator @ v)
|
|
||||||
return lm_eigval, v
|
|
||||||
|
|
||||||
|
|
||||||
def rayleigh_quotient_iteration(operator: sparse.spmatrix or spalg.LinearOperator,
|
|
||||||
guess_vector: numpy.ndarray,
|
|
||||||
iterations: int = 40,
|
|
||||||
tolerance: float = 1e-13,
|
|
||||||
solver=None,
|
|
||||||
) -> Tuple[complex, numpy.ndarray]:
|
|
||||||
"""
|
|
||||||
Use Rayleigh quotient iteration to refine an eigenvector guess.
|
|
||||||
|
|
||||||
:param operator: Matrix to analyze.
|
|
||||||
:param guess_vector: Eigenvector to refine.
|
|
||||||
:param iterations: Maximum number of iterations to perform. Default 40.
|
|
||||||
:param tolerance: Stop iteration if (A - I*eigenvalue) @ v < tolerance.
|
|
||||||
Default 1e-13.
|
|
||||||
:param solver: Solver function of the form x = solver(A, b).
|
|
||||||
By default, use scipy.sparse.spsolve for sparse matrices and
|
|
||||||
scipy.sparse.bicgstab for general LinearOperator instances.
|
|
||||||
:return: (eigenvalue, eigenvector)
|
|
||||||
"""
|
|
||||||
try:
|
|
||||||
_test = operator - sparse.eye(operator.shape[0])
|
|
||||||
shift = lambda eigval: eigval * sparse.eye(operator.shape[0])
|
|
||||||
if solver is None:
|
|
||||||
solver = spalg.spsolve
|
|
||||||
except TypeError:
|
|
||||||
shift = lambda eigval: spalg.LinearOperator(shape=operator.shape,
|
|
||||||
dtype=operator.dtype,
|
|
||||||
matvec=lambda v: eigval * v)
|
|
||||||
if solver is None:
|
|
||||||
solver = lambda A, b: spalg.bicgstab(A, b)[0]
|
|
||||||
|
|
||||||
v = guess_vector
|
|
||||||
v /= norm(v)
|
|
||||||
for _ in range(iterations):
|
|
||||||
eigval = v.conj() @ (operator @ v)
|
|
||||||
if norm(operator @ v - eigval * v) < tolerance:
|
|
||||||
break
|
|
||||||
|
|
||||||
shifted_operator = operator - shift(eigval)
|
|
||||||
v = solver(shifted_operator, v)
|
|
||||||
v /= norm(v)
|
|
||||||
return eigval, v
|
|
||||||
|
|
||||||
|
|
||||||
def signed_eigensolve(operator: sparse.spmatrix or spalg.LinearOperator,
|
|
||||||
how_many: int,
|
|
||||||
negative: bool = False,
|
|
||||||
) -> Tuple[numpy.ndarray, numpy.ndarray]:
|
|
||||||
"""
|
|
||||||
Find the largest-magnitude positive-only (or negative-only) eigenvalues and
|
|
||||||
eigenvectors of the provided matrix.
|
|
||||||
|
|
||||||
:param operator: Matrix to analyze.
|
|
||||||
:param how_many: How many eigenvalues to find.
|
|
||||||
:param negative: Whether to find negative-only eigenvalues.
|
|
||||||
Default False (positive only).
|
|
||||||
:return: (sorted list of eigenvalues, 2D ndarray of corresponding eigenvectors)
|
|
||||||
eigenvectors[:, k] corresponds to the k-th eigenvalue
|
|
||||||
"""
|
|
||||||
# Use power iteration to estimate the dominant eigenvector
|
|
||||||
lm_eigval, _ = power_iteration(operator)
|
|
||||||
|
|
||||||
'''
|
|
||||||
Shift by the absolute value of the largest eigenvalue, then find a few of the
|
|
||||||
largest-magnitude (shifted) eigenvalues. A positive shift ensures that we find the
|
|
||||||
largest _positive_ eigenvalues, since any negative eigenvalues will be shifted to the
|
|
||||||
range 0 >= neg_eigval + abs(lm_eigval) > abs(lm_eigval)
|
|
||||||
'''
|
|
||||||
shift = numpy.abs(lm_eigval)
|
|
||||||
if negative:
|
|
||||||
shift *= -1
|
|
||||||
|
|
||||||
# Try to combine, use general LinearOperator if we fail
|
|
||||||
try:
|
|
||||||
shifted_operator = operator + shift * sparse.eye(operator.shape[0])
|
|
||||||
except TypeError:
|
|
||||||
shifted_operator = operator + spalg.LinearOperator(shape=operator.shape,
|
|
||||||
matvec=lambda v: shift * v)
|
|
||||||
|
|
||||||
shifted_eigenvalues, eigenvectors = spalg.eigs(shifted_operator, which='LM', k=how_many, ncv=50)
|
|
||||||
eigenvalues = shifted_eigenvalues - shift
|
|
||||||
|
|
||||||
k = eigenvalues.argsort()
|
|
||||||
return eigenvalues[k], eigenvectors[:, k]
|
|
||||||
|
|
@ -1,220 +0,0 @@
|
|||||||
"""
|
|
||||||
Functions for performing near-to-farfield transformation (and the reverse).
|
|
||||||
"""
|
|
||||||
from typing import Dict, List
|
|
||||||
import numpy
|
|
||||||
from numpy.fft import fft2, fftshift, fftfreq, ifft2, ifftshift
|
|
||||||
from numpy import pi
|
|
||||||
|
|
||||||
|
|
||||||
def near_to_farfield(E_near: List[numpy.ndarray],
|
|
||||||
H_near: List[numpy.ndarray],
|
|
||||||
dx: float,
|
|
||||||
dy: float,
|
|
||||||
padded_size: List[int] = None
|
|
||||||
) -> Dict[str]:
|
|
||||||
"""
|
|
||||||
Compute the farfield, i.e. the distribution of the fields after propagation
|
|
||||||
through several wavelengths of uniform medium.
|
|
||||||
|
|
||||||
The input fields should be complex phasors.
|
|
||||||
|
|
||||||
:param E_near: List of 2 ndarrays containing the 2D phasor field slices for the transverse
|
|
||||||
E fields (e.g. [Ex, Ey] for calculating the farfield toward the z-direction).
|
|
||||||
:param H_near: List of 2 ndarrays containing the 2D phasor field slices for the transverse
|
|
||||||
H fields (e.g. [Hx, hy] for calculating the farfield towrad the z-direction).
|
|
||||||
:param dx: Cell size along x-dimension, in units of wavelength.
|
|
||||||
:param dy: Cell size along y-dimension, in units of wavelength.
|
|
||||||
:param padded_size: Shape of the output. A single integer `n` will be expanded to `(n, n)`.
|
|
||||||
Powers of 2 are most efficient for FFT computation.
|
|
||||||
Default is the smallest power of 2 larger than the input, for each axis.
|
|
||||||
:returns: Dict with keys
|
|
||||||
'E_far': Normalized E-field farfield; multiply by
|
|
||||||
(i k exp(-i k r) / (4 pi r)) to get the actual field value.
|
|
||||||
'H_far': Normalized H-field farfield; multiply by
|
|
||||||
(i k exp(-i k r) / (4 pi r)) to get the actual field value.
|
|
||||||
'kx', 'ky': Wavevector values corresponding to the x- and y- axes in E_far and H_far,
|
|
||||||
normalized to wavelength (dimensionless).
|
|
||||||
'dkx', 'dky': step size for kx and ky, normalized to wavelength.
|
|
||||||
'theta': arctan2(ky, kx) corresponding to each (kx, ky).
|
|
||||||
This is the angle in the x-y plane, counterclockwise from above, starting from +x.
|
|
||||||
'phi': arccos(kz / k) corresponding to each (kx, ky).
|
|
||||||
This is the angle away from +z.
|
|
||||||
"""
|
|
||||||
|
|
||||||
if not len(E_near) == 2:
|
|
||||||
raise Exception('E_near must be a length-2 list of ndarrays')
|
|
||||||
if not len(H_near) == 2:
|
|
||||||
raise Exception('H_near must be a length-2 list of ndarrays')
|
|
||||||
|
|
||||||
s = E_near[0].shape
|
|
||||||
if not all(s == f.shape for f in E_near + H_near):
|
|
||||||
raise Exception('All fields must be the same shape!')
|
|
||||||
|
|
||||||
if padded_size is None:
|
|
||||||
padded_size = (2**numpy.ceil(numpy.log2(s))).astype(int)
|
|
||||||
if not hasattr(padded_size, '__len__'):
|
|
||||||
padded_size = (padded_size, padded_size)
|
|
||||||
|
|
||||||
En_fft = [fftshift(fft2(fftshift(Eni), s=padded_size)) for Eni in E_near]
|
|
||||||
Hn_fft = [fftshift(fft2(fftshift(Hni), s=padded_size)) for Hni in H_near]
|
|
||||||
|
|
||||||
# Propagation vectors kx, ky
|
|
||||||
k = 2 * pi
|
|
||||||
kxx = 2 * pi * fftshift(fftfreq(padded_size[0], dx))
|
|
||||||
kyy = 2 * pi * fftshift(fftfreq(padded_size[1], dy))
|
|
||||||
|
|
||||||
kx, ky = numpy.meshgrid(kxx, kyy, indexing='ij')
|
|
||||||
kxy2 = kx * kx + ky * ky
|
|
||||||
kxy = numpy.sqrt(kxy2)
|
|
||||||
kz = numpy.sqrt(k * k - kxy2)
|
|
||||||
|
|
||||||
sin_th = ky / kxy
|
|
||||||
cos_th = kx / kxy
|
|
||||||
cos_phi = kz / k
|
|
||||||
|
|
||||||
sin_th[numpy.logical_and(kx == 0, ky == 0)] = 0
|
|
||||||
cos_th[numpy.logical_and(kx == 0, ky == 0)] = 1
|
|
||||||
|
|
||||||
# Normalized vector potentials N, L
|
|
||||||
N = [-Hn_fft[1] * cos_phi * cos_th + Hn_fft[0] * cos_phi * sin_th,
|
|
||||||
Hn_fft[1] * sin_th + Hn_fft[0] * cos_th]
|
|
||||||
L = [ En_fft[1] * cos_phi * cos_th - En_fft[0] * cos_phi * sin_th,
|
|
||||||
-En_fft[1] * sin_th - En_fft[0] * cos_th]
|
|
||||||
|
|
||||||
E_far = [-L[1] - N[0],
|
|
||||||
L[0] - N[1]]
|
|
||||||
H_far = [-E_far[1],
|
|
||||||
E_far[0]]
|
|
||||||
|
|
||||||
theta = numpy.arctan2(ky, kx)
|
|
||||||
phi = numpy.arccos(cos_phi)
|
|
||||||
theta[numpy.logical_and(kx == 0, ky == 0)] = 0
|
|
||||||
phi[numpy.logical_and(kx == 0, ky == 0)] = 0
|
|
||||||
|
|
||||||
# Zero fields beyond valid (phi, theta)
|
|
||||||
invalid_ind = kxy2 >= k * k
|
|
||||||
theta[invalid_ind] = 0
|
|
||||||
phi[invalid_ind] = 0
|
|
||||||
for i in range(2):
|
|
||||||
E_far[i][invalid_ind] = 0
|
|
||||||
H_far[i][invalid_ind] = 0
|
|
||||||
|
|
||||||
outputs = {
|
|
||||||
'E': E_far,
|
|
||||||
'H': H_far,
|
|
||||||
'dkx': kx[1]-kx[0],
|
|
||||||
'dky': ky[1]-ky[0],
|
|
||||||
'kx': kx,
|
|
||||||
'ky': ky,
|
|
||||||
'theta': theta,
|
|
||||||
'phi': phi,
|
|
||||||
}
|
|
||||||
|
|
||||||
return outputs
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
def far_to_nearfield(E_far: List[numpy.ndarray],
|
|
||||||
H_far: List[numpy.ndarray],
|
|
||||||
dkx: float,
|
|
||||||
dky: float,
|
|
||||||
padded_size: List[int] = None
|
|
||||||
) -> Dict[str]:
|
|
||||||
"""
|
|
||||||
Compute the farfield, i.e. the distribution of the fields after propagation
|
|
||||||
through several wavelengths of uniform medium.
|
|
||||||
|
|
||||||
The input fields should be complex phasors.
|
|
||||||
|
|
||||||
:param E_far: List of 2 ndarrays containing the 2D phasor field slices for the transverse
|
|
||||||
E fields (e.g. [Ex, Ey] for calculating the nearfield toward the z-direction).
|
|
||||||
Fields should be normalized so that
|
|
||||||
E_far = E_far_actual / (i k exp(-i k r) / (4 pi r))
|
|
||||||
:param H_far: List of 2 ndarrays containing the 2D phasor field slices for the transverse
|
|
||||||
H fields (e.g. [Hx, hy] for calculating the nearfield toward the z-direction).
|
|
||||||
Fields should be normalized so that
|
|
||||||
H_far = H_far_actual / (i k exp(-i k r) / (4 pi r))
|
|
||||||
:param dkx: kx discretization, in units of wavelength.
|
|
||||||
:param dky: ky discretization, in units of wavelength.
|
|
||||||
:param padded_size: Shape of the output. A single integer `n` will be expanded to `(n, n)`.
|
|
||||||
Powers of 2 are most efficient for FFT computation.
|
|
||||||
Default is the smallest power of 2 larger than the input, for each axis.
|
|
||||||
:returns: Dict with keys
|
|
||||||
'E': E-field nearfield
|
|
||||||
'H': H-field nearfield
|
|
||||||
'dx', 'dy': spatial discretization, normalized to wavelength (dimensionless)
|
|
||||||
"""
|
|
||||||
|
|
||||||
if not len(E_far) == 2:
|
|
||||||
raise Exception('E_far must be a length-2 list of ndarrays')
|
|
||||||
if not len(H_far) == 2:
|
|
||||||
raise Exception('H_far must be a length-2 list of ndarrays')
|
|
||||||
|
|
||||||
s = E_far[0].shape
|
|
||||||
if not all(s == f.shape for f in E_far + H_far):
|
|
||||||
raise Exception('All fields must be the same shape!')
|
|
||||||
|
|
||||||
if padded_size is None:
|
|
||||||
padded_size = (2**numpy.ceil(numpy.log2(s))).astype(int)
|
|
||||||
if not hasattr(padded_size, '__len__'):
|
|
||||||
padded_size = (padded_size, padded_size)
|
|
||||||
|
|
||||||
|
|
||||||
k = 2 * pi
|
|
||||||
kxs = fftshift(fftfreq(s[0], 1/(s[0] * dkx)))
|
|
||||||
kys = fftshift(fftfreq(s[0], 1/(s[1] * dky)))
|
|
||||||
|
|
||||||
kx, ky = numpy.meshgrid(kxs, kys, indexing='ij')
|
|
||||||
kxy2 = kx * kx + ky * ky
|
|
||||||
kxy = numpy.sqrt(kxy2)
|
|
||||||
|
|
||||||
kz = numpy.sqrt(k * k - kxy2)
|
|
||||||
|
|
||||||
sin_th = ky / kxy
|
|
||||||
cos_th = kx / kxy
|
|
||||||
cos_phi = kz / k
|
|
||||||
|
|
||||||
sin_th[numpy.logical_and(kx == 0, ky == 0)] = 0
|
|
||||||
cos_th[numpy.logical_and(kx == 0, ky == 0)] = 1
|
|
||||||
|
|
||||||
# Zero fields beyond valid (phi, theta)
|
|
||||||
invalid_ind = kxy2 >= k * k
|
|
||||||
theta[invalid_ind] = 0
|
|
||||||
phi[invalid_ind] = 0
|
|
||||||
for i in range(2):
|
|
||||||
E_far[i][invalid_ind] = 0
|
|
||||||
H_far[i][invalid_ind] = 0
|
|
||||||
|
|
||||||
|
|
||||||
# Normalized vector potentials N, L
|
|
||||||
L = [0.5 * E_far[1],
|
|
||||||
-0.5 * E_far[0]]
|
|
||||||
N = [L[1],
|
|
||||||
-L[0]]
|
|
||||||
|
|
||||||
En_fft = [-( L[0] * sin_th + L[1] * cos_phi * cos_th)/cos_phi,
|
|
||||||
-(-L[0] * cos_th + L[1] * cos_phi * sin_th)/cos_phi]
|
|
||||||
|
|
||||||
Hn_fft = [( N[0] * sin_th + N[1] * cos_phi * cos_th)/cos_phi,
|
|
||||||
(-N[0] * cos_th + N[1] * cos_phi * sin_th)/cos_phi]
|
|
||||||
|
|
||||||
for i in range(2):
|
|
||||||
En_fft[i][cos_phi == 0] = 0
|
|
||||||
Hn_fft[i][cos_phi == 0] = 0
|
|
||||||
|
|
||||||
E_near = [ifftshift(ifft2(ifftshift(Ei), s=padded_size)) for Ei in En_fft]
|
|
||||||
H_near = [ifftshift(ifft2(ifftshift(Hi), s=padded_size)) for Hi in Hn_fft]
|
|
||||||
|
|
||||||
dx = 2 * pi / (s[0] * dkx)
|
|
||||||
dy = 2 * pi / (s[0] * dky)
|
|
||||||
|
|
||||||
outputs = {
|
|
||||||
'E': E_near,
|
|
||||||
'H': H_near,
|
|
||||||
'dx': dx,
|
|
||||||
'dy': dy,
|
|
||||||
}
|
|
||||||
|
|
||||||
return outputs
|
|
||||||
|
|
Loading…
Reference in New Issue