472 lines
17 KiB
Python
472 lines
17 KiB
Python
"""
|
||
Mostly legacy code, but still kept for reference.
|
||
Might be removed in future versions.
|
||
"""
|
||
|
||
import numpy as np
|
||
use_moble_quaternion = False
|
||
try:
|
||
import quaternion, spherical_functions as sf # because of the Wigner matrices. These imports are SLOW.
|
||
use_moble_quaternion = True
|
||
except ImportError:
|
||
use_moble_quaternion = False
|
||
|
||
import re
|
||
from .constants import hbar, eV, pi, c
|
||
from .cycommon import get_mn_y, get_nelem
|
||
from .cyquaternions import CQuat
|
||
ň = np.newaxis
|
||
from .types import NormalizationT, TMatrixSpec
|
||
|
||
# Transformations of spherical bases
|
||
def WignerD_mm(l, quat):
|
||
"""
|
||
Calculates Wigner D matrix (as an numpy (2*l+1,2*l+1)-shaped array)
|
||
for order l, and a rotation given by quaternion quat.
|
||
|
||
This represents the rotation of spherical vector basis
|
||
TODO doc
|
||
"""
|
||
|
||
if use_moble_quaternion:
|
||
indices = np.array([ [l,i,j] for i in range(-l,l+1) for j in range(-l,l+1)])
|
||
Delems = sf.Wigner_D_element(quat, indices).reshape(2*l+1,2*l+1)
|
||
return Delems
|
||
else:
|
||
Delems = np.zeros((2*l+1, 2*l+1), dtype=complex)
|
||
for i in range(-l,l+1):
|
||
for j in range(-l,l+1):
|
||
Delems[i,j] = quat.wignerDelem(l, i, j)
|
||
return Delems
|
||
|
||
|
||
def WignerD_mm_fromvector(l, vect):
|
||
"""
|
||
TODO doc
|
||
"""
|
||
if use_moble_quaternion:
|
||
return WignerD_mm(l, quaternion.from_rotation_vector(vect))
|
||
else:
|
||
return WignerD_mm(l, CQuat.from_rotvector(vect))
|
||
|
||
|
||
def WignerD_yy(lmax, quat):
|
||
"""
|
||
TODO doc
|
||
"""
|
||
my, ny = get_mn_y(lmax)
|
||
Delems = np.zeros((len(my),len(my)),dtype=complex)
|
||
b_in = 0
|
||
e_in = None
|
||
for l in range(1,lmax+1):
|
||
e_in = b_in + 2*l+1
|
||
Delems[b_in:e_in,b_in:e_in] = WignerD_mm(l, quat)
|
||
b_in = e_in
|
||
return Delems
|
||
|
||
def WignerD_yy_fromvector(lmax, vect):
|
||
"""
|
||
TODO doc
|
||
"""
|
||
if use_moble_quaternion:
|
||
return WignerD_yy(lmax, quaternion.from_rotation_vector(vect))
|
||
else:
|
||
return WignerD_yy(lMax, CQuat.from_rotvector(vect))
|
||
|
||
def identity_yy(lmax):
|
||
"""
|
||
TODO doc
|
||
"""
|
||
return np.eye(lMax2nelem(lMax))
|
||
|
||
def identity_tyty(lmax):
|
||
"""
|
||
TODO doc
|
||
"""
|
||
nelem = lMax2nelem(lmax)
|
||
return np.eye(2*nelem).reshape((2,nelem,2,nelem))
|
||
|
||
def xflip_yy(lmax):
|
||
"""
|
||
TODO doc
|
||
xflip = δ(m + m') δ(l - l')
|
||
(i.e. ones on the (m' m) antidiagonal
|
||
"""
|
||
my, ny = get_mn_y(lmax)
|
||
elems = np.zeros((len(my),len(my)),dtype=int)
|
||
b_in = 0
|
||
e_in = None
|
||
for l in range(1,lmax+1):
|
||
e_in = b_in + 2*l+1
|
||
elems[b_in:e_in,b_in:e_in] = np.eye(2*l+1)[::-1,:]
|
||
b_in = e_in
|
||
return elems
|
||
|
||
def xflip_tyy(lmax):
|
||
fl_yy = xflip_yy(lmax)
|
||
return np.array([fl_yy,-fl_yy])
|
||
|
||
def xflip_tyty(lmax):
|
||
fl_yy = xflip_yy(lmax)
|
||
nelem = fl_yy.shape[0]
|
||
fl_tyty = np.zeros((2,nelem,2,nelem),dtype=int)
|
||
fl_tyty[0,:,0,:] = fl_yy
|
||
fl_tyty[1,:,1,:] = -fl_yy
|
||
return fl_tyty
|
||
|
||
def yflip_yy(lmax):
|
||
"""
|
||
TODO doc
|
||
yflip = rot(z,pi/2) * xflip * rot(z,-pi/2)
|
||
= δ(m + m') δ(l - l') * (-1)**m
|
||
"""
|
||
my, ny = get_mn_y(lmax)
|
||
elems = xflip_yy(lmax)
|
||
elems[(my % 2)==1] = elems[(my % 2)==1] * -1 # Obvious sign of tiredness (this is correct but ugly; FIXME)
|
||
return elems
|
||
|
||
def yflip_tyy(lmax):
|
||
fl_yy = yflip_yy(lmax)
|
||
return np.array([fl_yy,-fl_yy])
|
||
|
||
def yflip_tyty(lmax):
|
||
fl_yy = yflip_yy(lmax)
|
||
nelem = fl_yy.shape[0]
|
||
fl_tyty = np.zeros((2,nelem,2,nelem),dtype=int)
|
||
fl_tyty[0,:,0,:] = fl_yy
|
||
fl_tyty[1,:,1,:] = -fl_yy
|
||
return fl_tyty
|
||
|
||
def zflip_yy(lmax):
|
||
"""
|
||
TODO doc
|
||
zflip = (-1)^(l+m)
|
||
"""
|
||
my, ny = get_mn_y(lmax)
|
||
elems = np.zeros((len(my), len(my)), dtype=int)
|
||
b_in = 0
|
||
e_in = None
|
||
for l in range(1,lmax+1):
|
||
e_in = b_in + 2*l+1
|
||
elems[b_in:e_in,b_in:e_in] = np.diag([(-1)**i for i in range(e_in-b_in)])
|
||
b_in = e_in
|
||
return elems
|
||
|
||
def zflip_tyy(lmax):
|
||
fl_yy = zflip_yy(lmax)
|
||
return np.array([fl_yy,-fl_yy])
|
||
|
||
def zflip_tyty(lmax):
|
||
fl_yy = zflip_yy(lmax)
|
||
nelem = fl_yy.shape[0]
|
||
fl_tyty = np.zeros((2,nelem,2,nelem),dtype=int)
|
||
fl_tyty[0,:,0,:] = fl_yy
|
||
fl_tyty[1,:,1,:] = -fl_yy
|
||
return fl_tyty
|
||
|
||
def zrotN_yy(N, lMax):
|
||
return WignerD_yy_fromvector(lMax, np.array([0,0,pi * (2/N)]))
|
||
|
||
def op_yy2tyty(yyop):
|
||
'''
|
||
Broadcasts an yy operator to tyty operator without considering mirroring effects.
|
||
Good (maybe only) for rotations.
|
||
'''
|
||
return np.moveaxis(np.eye(2)[:,:,ň,ň] * yyop, 2,1)
|
||
|
||
def zrotN_tyty(N, lMax):
|
||
return op_yy2tyty(zrotN_yy(N, lMax))
|
||
|
||
def parity_yy(lmax):
|
||
"""
|
||
Parity operator (flip in x,y,z)
|
||
parity = (-1)**l
|
||
"""
|
||
my, ny = get_mn_y(lmax)
|
||
return np.diag((-1)**ny)
|
||
|
||
# BTW parity (xyz-flip) is simply (-1)**ny
|
||
|
||
|
||
|
||
#----------------------------------------------------#
|
||
# Loading T-matrices from scuff-tmatrix output files #
|
||
#----------------------------------------------------#
|
||
|
||
# We don't really need this particular function anymore, but...
|
||
def _scuffTMatrixConvert_EM_01(EM):
|
||
#print(EM)
|
||
if (EM == b'E'):
|
||
return 1
|
||
elif (EM == b'M'):
|
||
return 0
|
||
else:
|
||
return None
|
||
|
||
def loadScuffTMatrices(fileName, normalisation = 1, version = 'old', freqscale = None, order = None):
|
||
"""
|
||
TODO doc
|
||
|
||
version describes version of scuff-em. It is either 'old' or 'new'.
|
||
default order is ('N','M') with 'old' version, ('M','N') with 'new'
|
||
"""
|
||
oldversion = (version == 'old')
|
||
μm = 1e-6
|
||
table = np.genfromtxt(fileName,
|
||
converters={1: _scuffTMatrixConvert_EM_01, 4: _scuffTMatrixConvert_EM_01} if oldversion else None,
|
||
skip_header = 0 if oldversion else 5,
|
||
usecols = None if oldversion else (0, 2, 3, 4, 6, 7, 8, 9, 10),
|
||
dtype=[('freq', '<f8'),
|
||
('outc_type', '<i8'), ('outc_l', '<i8'), ('outc_m', '<i8'),
|
||
('inc_type', '<i8'), ('inc_l', '<i8'), ('inc_m', '<i8'),
|
||
('Treal', '<f8'), ('Timag', '<f8')]
|
||
if oldversion else
|
||
[('freq', '<f8'),
|
||
('outc_l', '<i8'), ('outc_m', '<i8'), ('outc_type', '<i8'),
|
||
('inc_l', '<i8'), ('inc_m', '<i8'), ('inc_type', '<i8'),
|
||
('Treal', '<f8'), ('Timag', '<f8')]
|
||
)
|
||
lMax=np.max(table['outc_l'])
|
||
my,ny = get_mn_y(lMax)
|
||
nelem = len(ny)
|
||
TMatrix_sz = nelem**2 * 4 # number of rows for each frequency: nelem * nelem spherical incides, 2 * 2 E/M types
|
||
freqs_weirdunits = table['freq'][::TMatrix_sz].copy()
|
||
freqs = freqs_weirdunits * c / μm
|
||
# The iteration in the TMatrix file goes in this order (the last one iterates fastest, i.e. in the innermost loop):
|
||
# freq outc_l outc_m outc_type inc_l inc_m inc_type
|
||
# The l,m mapping is the same as is given by my get_mn_y function, so no need to touch that
|
||
TMatrices_tmp_real = table['Treal'].reshape(len(freqs), nelem, 2, nelem, 2)
|
||
TMatrices_tmp_imag = table['Timag'].reshape(len(freqs), nelem, 2, nelem, 2)
|
||
# There are two přoblems with the previous matrices. First, we want to have the
|
||
# type indices first, so we want a shape (len(freqs), 2, nelem, 2, nelem) as in the older code.
|
||
# Second, M-waves come first, so they have now 0-valued index, and E-waves have 1-valued index,
|
||
# which we want to be inverted.
|
||
TMatrices = np.zeros((len(freqs),2,nelem,2,nelem),dtype=complex)
|
||
reorder = (0,1)
|
||
if ((order == ('N', 'M')) and not oldversion): # reverse order for the new version
|
||
reorder = (1,0)
|
||
# TODO reverse order for the old version...
|
||
for inc_type in (0,1):
|
||
for outc_type in (0,1):
|
||
TMatrices[:,reorder[outc_type],:,reorder[inc_type],:] = TMatrices_tmp_real[:,:,outc_type,:,inc_type]+1j*TMatrices_tmp_imag[:,:,outc_type,:,inc_type]
|
||
# IMPORTANT: now we are going from Reid's/Kristensson's/Jackson's/whoseever convention to Taylor's convention
|
||
# TODO make these consistent with what is defined in qpms_types.h (implement all possibilities)
|
||
if normalisation == 1:
|
||
TMatrices[:,:,:,:,:] = TMatrices[:,:,:,:,:] * np.sqrt(ny*(ny+1))[ň,ň,ň,ň,:] / np.sqrt(ny*(ny+1))[ň,ň,:,ň,ň]
|
||
elif normalisation == 2: # Kristensson?
|
||
pass
|
||
if freqscale is not None:
|
||
freqs *= freqscale
|
||
freqs_weirdunits *= freqscale
|
||
return (TMatrices, freqs, freqs_weirdunits, lMax)
|
||
|
||
|
||
# misc tensor maniputalion
|
||
def apply_matrix_left(matrix, tensor, axis):
|
||
"""
|
||
TODO doc
|
||
Apply square matrix to a given axis of a tensor, so that the result retains the shape
|
||
of the original tensor. The summation goes over the second index of the matrix and the
|
||
given tensor axis.
|
||
"""
|
||
tmp = np.tensordot(matrix, tensor, axes=(-1,axis))
|
||
return np.moveaxis(tmp, 0, axis)
|
||
|
||
def apply_ndmatrix_left(matrix,tensor,axes):
|
||
"""
|
||
Generalized apply_matrix_left, the matrix can have more (2N) abstract dimensions,
|
||
like M[i,j,k,...z,i,j,k,...,z]. N axes have to be specified in a tuple, corresponding
|
||
to the axes 0,1,...N-1 of the matrix
|
||
"""
|
||
N = len(axes)
|
||
matrix = np.tensordot(matrix, tensor, axes=([-N+axn for axn in range(N)],axes))
|
||
matrix = np.moveaxis(matrix, range(N), axes)
|
||
return matrix
|
||
|
||
def apply_ndmatrix_right(tensor, matrix, axes):
|
||
"""
|
||
Right-side analogue of apply_ndmatrix_lift.
|
||
Multiplies a tensor with a 2N-dimensional matrix, conserving the axis order.
|
||
"""
|
||
N = len(axes)
|
||
matrix = np.tensordot(tensor, matrix, axes = (axes, range(N)))
|
||
matrix = np.moveaxis(matrix, [-N+axn for axn in range(N)], axes)
|
||
return matrix
|
||
|
||
def ndmatrix_Hconj(matrix):
|
||
"""
|
||
For 2N-dimensional matrix, swap the first N and last N matrix, and complex conjugate.
|
||
"""
|
||
twoN = len(matrix.shape)
|
||
if not twoN % 2 == 0:
|
||
raise ValueError("The matrix has to have even number of axes.")
|
||
N = twoN//2
|
||
matrix = np.moveaxis(matrix, range(N), range(N, 2*N))
|
||
return matrix.conj()
|
||
|
||
def symz_indexarrays(lMax, npart = 1):
|
||
"""
|
||
Returns indices that are used for separating the in-plane E ('TE' in the photonic crystal
|
||
jargon) and perpendicular E ('TM' in the photonic crystal jargon) modes
|
||
in the z-mirror symmetric systems.
|
||
|
||
Parameters
|
||
----------
|
||
lMax : int
|
||
The maximum degree cutoff for the T-matrix to which these indices will be applied.
|
||
|
||
npart : int
|
||
Number of particles (TODO better description)
|
||
|
||
Returns
|
||
-------
|
||
TEč, TMč : (npart * 2 * nelem)-shaped bool ndarray
|
||
Mask arrays corresponding to the 'TE' and 'TM' modes, respectively.
|
||
"""
|
||
my, ny = get_mn_y(lMax)
|
||
nelem = len(my)
|
||
ž = np.arange(2*nelem) # single particle spherical wave indices
|
||
tž = ž // nelem # tž == 0: electric waves, tž == 1: magnetic waves
|
||
mž = my[ž%nelem]
|
||
nž = ny[ž%nelem]
|
||
TEž = ž[(mž+nž+tž) % 2 == 0]
|
||
TMž = ž[(mž+nž+tž) % 2 == 1]
|
||
|
||
č = np.arange(npart*2*nelem) # spherical wave indices for multiple particles (e.g. in a unit cell)
|
||
žč = č % (2* nelem)
|
||
tč = tž[žč]
|
||
mč = mž[žč]
|
||
nč = nž[žč]
|
||
TEč = č[(mč+nč+tč) % 2 == 0]
|
||
TMč = č[(mč+nč+tč) % 2 == 1]
|
||
return (TEč, TMč)
|
||
|
||
|
||
"""
|
||
Processing T-matrix related operations from scripts
|
||
===================================================
|
||
|
||
see also scripts_common.py
|
||
|
||
The unit cell is defined by a dict particle_specs and a list TMatrix_specs.
|
||
This particular module has to provide the T-matrices according to what is defined
|
||
in TMatrix_specs
|
||
|
||
TMatrix_specs is a list of tuples (lMax_override, TMatrix_path, ops)
|
||
where
|
||
- TMatrix_path is path to the file generated by scuff-tmatrix
|
||
- lMax_override is int or None; if it is int and less than the lMax found from the T-matrix file,
|
||
lMax_override is used as the order cutoff for the output T-matrix.
|
||
- ops is an iterable of tuples (optype, opargs) where currently optype can be 'sym' or 'tr'
|
||
for symmetrization operation or some other transformation.
|
||
"""
|
||
|
||
#TODO FEATURE more basic group symmetry operations, cf. http://symmetry.otterbein.edu/tutorial/index.html
|
||
|
||
# This is for finite „fractional“ rotations along the z-axis (mCN means rotation of 2π*(m/N))
|
||
reCN = re.compile('(\d*)C(\d+)') # TODO STYLE make this regexp also accept the 3*C_5-type input, eqiv. to 3C5
|
||
|
||
def get_TMatrix_fromspec(tmatrix_spec):
|
||
''' TODO doc
|
||
returns (TMatrices, freqs, lMax)
|
||
'''
|
||
lMax_override, tmpath, ops = tmatrix_spec
|
||
TMatrices, freqs, freqs_weirdunits, lMax = loadScuffTMatrices(tmpath)
|
||
if lMax_override is not None and (lMax_override < lMax_orig):
|
||
nelem = get_nelem(lMax_override)
|
||
TMatrices = TMatrices[...,0:nelem,:,0:nelem]
|
||
lMax = lMax_override
|
||
|
||
for (optype, opargs) in ops:
|
||
if optype == 'sym':
|
||
mCN = reCN.match(opargs)
|
||
if opargs == 'C2' or opargs == 'C_2':
|
||
opmat = apply_matrix_left(yflip_yy(lMax), xflip_yy(lMax), -1)
|
||
TMatrices = (TMatrices + apply_matrix_left(opmat, apply_matrix_left(opmat, TMatrices, -3), -1))/2
|
||
elif opargs == 'σ_x' or opargs == 'σx':
|
||
opmat = xflip_tyty(lMax)
|
||
TMatrices = (TMatrices + apply_ndmatrix_left(opmat, apply_ndmatrix_left(opmat, TMatrices, (-4,-3)),(-2,-1)))/2
|
||
elif opargs == 'σ_y' or opargs == 'σy':
|
||
opmat = yflip_tyty(lMax)
|
||
TMatrices = (TMatrices + apply_ndmatrix_left(opmat, apply_ndmatrix_left(opmat, TMatrices, (-4,-3)),(-2,-1)))/2
|
||
elif opargs == 'σ_z' or opargs == 'σz':
|
||
opmat = zflip_tyty(lMax)
|
||
TMatrices = (TMatrices + apply_ndmatrix_left(opmat, apply_ndmatrix_left(opmat, TMatrices, (-4,-3)),(-2,-1)))/2
|
||
elif mCN:
|
||
rotN = int(mCN.group(2)) # the possible m is ignored
|
||
TMatrix_contribs = np.empty((rotN,)+TMatrices.shape, dtype=np.complex_)
|
||
for i in range(rotN):
|
||
rotangle = 2*np.pi*i / rotN
|
||
rot = WignerD_yy_fromvector(lMax, np.array([0,0,rotangle]))
|
||
rotinv = WignerD_yy_fromvector(lMax, np.array([0,0,-rotangle]))
|
||
TMatrix_contribs[i] = apply_matrix_left(rot, apply_matrix_left(rotinv, TMatrices, -3), -1)
|
||
TMatrices = np.sum(TMatrix_contribs, axis=0) / rotN
|
||
elif opargs == 'C_inf' or opargs == 'Cinf' or opargs == 'C_∞' or opargs == 'C∞':
|
||
raise ValueError('not implemented: ', opargs) # TODO FEATURE
|
||
else:
|
||
raise ValueError('not implemented: ', opargs)
|
||
elif optype == 'tr':
|
||
mCN = reCN.match(opargs)
|
||
if opargs == 'C2' or opargs == 'C_2':
|
||
opmat = apply_matrix_left(yflip_yy(lMax), xflip_yy(lMax), -1)
|
||
TMatrices = apply_matrix_left(opmat, apply_matrix_left(opmat, TMatrices, -3), -1)/2
|
||
elif opargs == 'σ_x' or opargs == 'σx':
|
||
opmat = xflip_tyty(lMax)
|
||
TMatrices = apply_ndmatrix_left(opmat, apply_ndmatrix_left(opmat, TMatrices, (-4,-3)),(-2,-1))/2
|
||
elif opargs == 'σ_y' or opargs == 'σy':
|
||
opmat = yflip_tyty(lMax)
|
||
TMatrices = apply_ndmatrix_left(opmat, apply_ndmatrix_left(opmat, TMatrices, (-4,-3)),(-2,-1))/2
|
||
elif opargs == 'σ_z' or opargs == 'σz':
|
||
opmat = zflip_tyty(lMax)
|
||
TMatrices = apply_ndmatrix_left(opmat, apply_ndmatrix_left(opmat, TMatrices, (-4,-3)),(-2,-1))/2
|
||
elif mCN:
|
||
rotN = int(mCN.group(2))
|
||
power = int(mCN.group(1)) if mCN.group(1) else 1
|
||
rotangle = 2*np.pi*power / rotN
|
||
rot = WignerD_yy_fromvector(lMax, np.array([0,0,rotangle]))
|
||
rotinv = WignerD_yy_fromvector(lMax, np.array([0,0,-rotangle]))
|
||
TMatrices = apply_matrix_left(rot, apply_matrix_left(rotinv, TMatrices, -3), -1)
|
||
else:
|
||
raise ValueError('not implemented: ', opargs)
|
||
else:
|
||
raise ValueError('not implemented: ', optype)
|
||
return (TMatrices, freqs, lMax)
|
||
|
||
class TMatrix(TMatrixSpec):
|
||
'''
|
||
TODO doc
|
||
|
||
TODO support for different/multiple interpolators
|
||
'''
|
||
def __init__(self, tmatrix_spec):
|
||
from scipy import interpolate
|
||
#self.specification = tmatrix_spec
|
||
self.lMax_override = tmatrix_spec.lMax_override
|
||
self.tmatrix_path = tmatrix_spec.tmatrix_path
|
||
self.ops = tmatrix_spec.ops
|
||
self.tmdata, self.freqs, self.lMax = get_TMatrix_fromspec(tmatrix_spec)
|
||
self.nelem = get_nelem(self.lMax)
|
||
#self._interpolators = dict()
|
||
self.default_interpolator = interpolate.interp1d(self.freqs,
|
||
self.tmdata, axis=0, kind='linear', fill_value='extrapolate')
|
||
self.normalization = NormalizationT.TAYLOR # TODO others are not supported by the loading functions
|
||
|
||
def atfreq(self, freq):
|
||
freqarray = np.array(freq, copy=False)
|
||
if freqarray.shape: # not just a scalar
|
||
tm_interp = np.empty(freqarray.shape + (2, self.nelem, 2, self.nelem), dtype=np.complex_)
|
||
for i in np.ndindex(freqarray.shape):
|
||
tm_interp[i] = self.default_interpolator(freqarray[i])
|
||
return tm_interp
|
||
else: # scalar
|
||
return self.default_interpolator(freq)
|
||
|
||
__getitem__ = atfreq # might be changed later, use atfreq to be sure
|
||
|
||
def perform_tmspecs(tmspecs):
|
||
"""Takes a sequence of TMatrixSpec or TMatrix instances and returns
|
||
a list of corresponding TMatrix instances"""
|
||
return [(tmspec if hasattr(tmspec, "tmdata") else TMatrix(tmspec))
|
||
for tmspec in tmspecs]
|
||
|