qpms/tmatrices.py

   1 """
   2 Mostly legacy code, but still kept for reference.
   3 Might be removed in future versions.
   4 """
   5
   6 import numpy as np
   7 use_moble_quaternion = False
   8 try:
   9     import quaternion, spherical_functions as sf # because of the Wigner matrices. These imports are SLOW.
  10     use_moble_quaternion = True
  11 except ImportError:
  12     use_moble_quaternion = False
  13
  14 import re
  15 from .constants import hbar, eV, pi, c
  16 from .cycommon import get_mn_y, get_nelem
  17 from .cyquaternions import CQuat
  18 ň = np.newaxis
  19 from .types import NormalizationT, TMatrixSpec
  20
  21 # Transformations of spherical bases
  22 def WignerD_mm(l, quat):
  23     """
  24     Calculates Wigner D matrix (as an numpy (2*l+1,2*l+1)-shaped array)
  25     for order l, and a rotation given by quaternion quat.
  26
  27     This represents the rotation of spherical vector basis
  28     TODO doc
  29     """
  30
  31     if use_moble_quaternion:
  32         indices = np.array([ [l,i,j] for i in range(-l,l+1) for j in range(-l,l+1)])
  33         Delems = sf.Wigner_D_element(quat, indices).reshape(2*l+1,2*l+1)
  34         return Delems
  35     else:
  36         Delems = np.zeros((2*l+1, 2*l+1), dtype=complex)
  37         for i in range(-l,l+1):
  38             for j in range(-l,l+1):
  39                 Delems[i,j] = quat.wignerDelem(l, i, j)
  40         return Delems
  41
  42
  43 def WignerD_mm_fromvector(l, vect):
  44     """
  45     TODO doc
  46     """
  47     if use_moble_quaternion:
  48         return WignerD_mm(l, quaternion.from_rotation_vector(vect))
  49     else:
  50         return WignerD_mm(l, CQuat.from_rotvector(vect))
  51
  52
  53 def WignerD_yy(lmax, quat):
  54     """
  55     TODO doc
  56     """
  57     my, ny = get_mn_y(lmax)
  58     Delems = np.zeros((len(my),len(my)),dtype=complex)
  59     b_in = 0
  60     e_in =  None
  61     for l in range(1,lmax+1):
  62         e_in = b_in + 2*l+1
  63         Delems[b_in:e_in,b_in:e_in] = WignerD_mm(l, quat)
  64         b_in = e_in
  65     return Delems
  66
  67 def WignerD_yy_fromvector(lmax, vect):
  68     """
  69     TODO doc
  70     """
  71     if use_moble_quaternion:
  72         return WignerD_yy(lmax, quaternion.from_rotation_vector(vect))
  73     else:
  74         return WignerD_yy(lMax, CQuat.from_rotvector(vect))
  75
  76 def identity_yy(lmax):
  77     """
  78     TODO doc
  79     """
  80     return np.eye(lMax2nelem(lMax))
  81
  82 def identity_tyty(lmax):
  83     """
  84     TODO doc
  85     """
  86     nelem = lMax2nelem(lmax)
  87     return np.eye(2*nelem).reshape((2,nelem,2,nelem))
  88
  89 def xflip_yy(lmax):
  90     """
  91     TODO doc
  92     xflip = δ(m + m') δ(l - l')
  93     (i.e. ones on the (m' m) antidiagonal
  94     """
  95     my, ny = get_mn_y(lmax)
  96     elems = np.zeros((len(my),len(my)),dtype=int)
  97     b_in = 0
  98     e_in =  None
  99     for l in range(1,lmax+1):
 100         e_in = b_in + 2*l+1
 101         elems[b_in:e_in,b_in:e_in] = np.eye(2*l+1)[::-1,:]
 102         b_in = e_in
 103     return elems
 104
 105 def xflip_tyy(lmax):
 106     fl_yy = xflip_yy(lmax)
 107     return np.array([fl_yy,-fl_yy])
 108
 109 def xflip_tyty(lmax):
 110     fl_yy = xflip_yy(lmax)
 111     nelem = fl_yy.shape[0]
 112     fl_tyty = np.zeros((2,nelem,2,nelem),dtype=int)
 113     fl_tyty[0,:,0,:] = fl_yy
 114     fl_tyty[1,:,1,:] = -fl_yy
 115     return fl_tyty
 116
 117 def yflip_yy(lmax):
 118     """
 119     TODO doc
 120     yflip = rot(z,pi/2) * xflip * rot(z,-pi/2)
 121           = δ(m + m') δ(l - l') * (-1)**m
 122     """
 123     my, ny = get_mn_y(lmax)
 124     elems = xflip_yy(lmax)
 125     elems[(my % 2)==1] = elems[(my % 2)==1] * -1 # Obvious sign of tiredness (this is correct but ugly; FIXME)
 126     return elems
 127
 128 def yflip_tyy(lmax):
 129     fl_yy = yflip_yy(lmax)
 130     return np.array([fl_yy,-fl_yy])
 131
 132 def yflip_tyty(lmax):
 133     fl_yy = yflip_yy(lmax)
 134     nelem = fl_yy.shape[0]
 135     fl_tyty = np.zeros((2,nelem,2,nelem),dtype=int)
 136     fl_tyty[0,:,0,:] = fl_yy
 137     fl_tyty[1,:,1,:] = -fl_yy
 138     return fl_tyty
 139
 140 def zflip_yy(lmax):
 141     """
 142     TODO doc
 143     zflip = (-1)^(l+m)
 144     """
 145     my, ny = get_mn_y(lmax)
 146     elems = np.zeros((len(my), len(my)), dtype=int)
 147     b_in = 0
 148     e_in = None
 149     for l in range(1,lmax+1):
 150         e_in = b_in + 2*l+1
 151         elems[b_in:e_in,b_in:e_in] = np.diag([(-1)**i for i in range(e_in-b_in)])
 152         b_in = e_in
 153     return elems
 154
 155 def zflip_tyy(lmax):
 156     fl_yy = zflip_yy(lmax)
 157     return np.array([fl_yy,-fl_yy])
 158
 159 def zflip_tyty(lmax):
 160     fl_yy = zflip_yy(lmax)
 161     nelem = fl_yy.shape[0]
 162     fl_tyty = np.zeros((2,nelem,2,nelem),dtype=int)
 163     fl_tyty[0,:,0,:] = fl_yy
 164     fl_tyty[1,:,1,:] = -fl_yy
 165     return fl_tyty
 166
 167 def zrotN_yy(N, lMax):
 168     return WignerD_yy_fromvector(lMax, np.array([0,0,pi * (2/N)]))
 169
 170 def op_yy2tyty(yyop):
 171     '''
 172     Broadcasts an yy operator to tyty operator without considering mirroring effects.
 173     Good (maybe only) for rotations.
 174     '''
 175     return np.moveaxis(np.eye(2)[:,:,ň,ň] * yyop, 2,1)
 176
 177 def zrotN_tyty(N, lMax):
 178     return op_yy2tyty(zrotN_yy(N, lMax))
 179
 180 def parity_yy(lmax):
 181     """
 182     Parity operator (flip in x,y,z)
 183     parity = (-1)**l
 184     """
 185     my, ny = get_mn_y(lmax)
 186     return np.diag((-1)**ny)
 187
 188 # BTW parity (xyz-flip) is simply (-1)**ny
 189
 190
 191
 192 #----------------------------------------------------#
 193 # Loading T-matrices from scuff-tmatrix output files #
 194 #----------------------------------------------------#
 195
 196 # We don't really need this particular function anymore, but...
 197 def _scuffTMatrixConvert_EM_01(EM):
 198     #print(EM)
 199     if (EM == b'E'):
 200         return 1
 201     elif (EM == b'M'):
 202         return 0
 203     else:
 204         return None
 205
 206 def loadScuffTMatrices(fileName, normalisation = 1, version = 'old', freqscale = None, order = None):
 207     """
 208     TODO doc
 209
 210     version describes version of scuff-em. It is either 'old' or 'new'.
 211     default order is ('N','M') with 'old' version, ('M','N') with 'new'
 212     """
 213     oldversion = (version == 'old')
 214     μm = 1e-6
 215     table = np.genfromtxt(fileName,
 216                   converters={1: _scuffTMatrixConvert_EM_01, 4: _scuffTMatrixConvert_EM_01} if oldversion else None,
 217                   skip_header = 0 if oldversion else 5,
 218                   usecols = None if oldversion else (0, 2, 3, 4, 6, 7, 8, 9, 10),
 219                   dtype=[('freq', '<f8'),
 220                          ('outc_type', '<i8'), ('outc_l', '<i8'), ('outc_m', '<i8'),
 221                          ('inc_type', '<i8'), ('inc_l', '<i8'), ('inc_m', '<i8'),
 222                          ('Treal', '<f8'), ('Timag', '<f8')]
 223                    if oldversion else
 224                         [('freq', '<f8'),
 225                          ('outc_l', '<i8'), ('outc_m', '<i8'),  ('outc_type', '<i8'),
 226                          ('inc_l', '<i8'), ('inc_m', '<i8'), ('inc_type', '<i8'),
 227                          ('Treal', '<f8'), ('Timag', '<f8')]
 228                     )
 229     lMax=np.max(table['outc_l'])
 230     my,ny = get_mn_y(lMax)
 231     nelem = len(ny)
 232     TMatrix_sz = nelem**2 * 4 # number of rows for each frequency: nelem * nelem spherical incides, 2 * 2 E/M types
 233     freqs_weirdunits = table['freq'][::TMatrix_sz].copy()
 234     freqs = freqs_weirdunits * c / μm
 235     # The iteration in the TMatrix file goes in this order (the last one iterates fastest, i.e. in the innermost loop):
 236     # freq outc_l outc_m outc_type inc_l inc_m inc_type
 237     # The l,m mapping is the same as is given by my get_mn_y function, so no need to touch that
 238     TMatrices_tmp_real = table['Treal'].reshape(len(freqs), nelem, 2, nelem, 2)
 239     TMatrices_tmp_imag = table['Timag'].reshape(len(freqs), nelem, 2, nelem, 2)
 240     # There are two přoblems with the previous matrices. First, we want to have the
 241     # type indices first, so we want a shape (len(freqs), 2, nelem, 2, nelem) as in the older code.
 242     # Second, M-waves come first, so they have now 0-valued index, and E-waves have 1-valued index,
 243     # which we want to be inverted.
 244     TMatrices = np.zeros((len(freqs),2,nelem,2,nelem),dtype=complex)
 245     reorder = (0,1)
 246     if ((order == ('N', 'M')) and not oldversion): # reverse order for the new version
 247         reorder = (1,0)
 248     # TODO reverse order for the old version...
 249     for inc_type in (0,1):
 250         for outc_type in (0,1):
 251             TMatrices[:,reorder[outc_type],:,reorder[inc_type],:] = TMatrices_tmp_real[:,:,outc_type,:,inc_type]+1j*TMatrices_tmp_imag[:,:,outc_type,:,inc_type]
 252     # IMPORTANT: now we are going from Reid's/Kristensson's/Jackson's/whoseever convention to Taylor's convention
 253     # TODO make these consistent with what is defined in qpms_types.h (implement all possibilities)
 254     if normalisation == 1:
 255         TMatrices[:,:,:,:,:] = TMatrices[:,:,:,:,:] * np.sqrt(ny*(ny+1))[ň,ň,ň,ň,:] / np.sqrt(ny*(ny+1))[ň,ň,:,ň,ň]
 256     elif normalisation == 2: # Kristensson?
 257         pass
 258     if freqscale is not None:
 259         freqs *= freqscale
 260         freqs_weirdunits *= freqscale
 261     return (TMatrices, freqs, freqs_weirdunits, lMax)
 262
 263
 264 # misc tensor maniputalion
 265 def apply_matrix_left(matrix, tensor, axis):
 266     """
 267     TODO doc
 268     Apply square matrix to a given axis of a tensor, so that the result retains the shape
 269     of the original tensor. The summation goes over the second index of the matrix and the
 270     given tensor axis.
 271     """
 272     tmp = np.tensordot(matrix, tensor, axes=(-1,axis))
 273     return np.moveaxis(tmp, 0, axis)
 274
 275 def apply_ndmatrix_left(matrix,tensor,axes):
 276     """
 277     Generalized apply_matrix_left, the matrix can have more (2N) abstract dimensions,
 278     like M[i,j,k,...z,i,j,k,...,z]. N axes have to be specified in a tuple, corresponding
 279     to the axes 0,1,...N-1 of the matrix
 280     """
 281     N = len(axes)
 282     matrix = np.tensordot(matrix, tensor, axes=([-N+axn for axn in range(N)],axes))
 283     matrix = np.moveaxis(matrix, range(N), axes)
 284     return matrix
 285
 286 def apply_ndmatrix_right(tensor, matrix, axes):
 287     """
 288     Right-side analogue of apply_ndmatrix_lift.
 289     Multiplies a tensor with a 2N-dimensional matrix, conserving the axis order.
 290     """
 291     N = len(axes)
 292     matrix = np.tensordot(tensor, matrix, axes = (axes, range(N)))
 293     matrix = np.moveaxis(matrix, [-N+axn for axn in range(N)], axes)
 294     return matrix
 295
 296 def ndmatrix_Hconj(matrix):
 297     """
 298     For 2N-dimensional matrix, swap the first N and last N matrix, and complex conjugate.
 299     """
 300     twoN = len(matrix.shape)
 301     if not twoN % 2 == 0:
 302         raise ValueError("The matrix has to have even number of axes.")
 303     N = twoN//2
 304     matrix = np.moveaxis(matrix, range(N), range(N, 2*N))
 305     return matrix.conj()
 306
 307 def symz_indexarrays(lMax, npart = 1):
 308     """
 309     Returns indices that are used for separating the in-plane E ('TE' in the photonic crystal
 310     jargon) and perpendicular E ('TM' in the photonic crystal jargon) modes
 311     in the z-mirror symmetric systems.
 312
 313     Parameters
 314     ----------
 315     lMax : int
 316         The maximum degree cutoff for the T-matrix to which these indices will be applied.
 317
 318     npart : int
 319         Number of particles (TODO better description)
 320
 321     Returns
 322     -------
 323     TEč, TMč : (npart * 2 * nelem)-shaped bool ndarray
 324         Mask arrays corresponding to the 'TE' and 'TM' modes, respectively.
 325     """
 326     my, ny = get_mn_y(lMax)
 327     nelem = len(my)
 328     ž = np.arange(2*nelem) # single particle spherical wave indices
 329     tž = ž // nelem # tž == 0: electric waves, tž == 1: magnetic waves
 330     mž = my[ž%nelem]
 331     nž = ny[ž%nelem]
 332     TEž = ž[(mž+nž+tž) % 2 == 0]
 333     TMž = ž[(mž+nž+tž) % 2 == 1]
 334
 335     č = np.arange(npart*2*nelem) # spherical wave indices for multiple particles (e.g. in a unit cell)
 336     žč = č % (2* nelem)
 337     tč = tž[žč]
 338     mč = mž[žč]
 339     nč = nž[žč]
 340     TEč = č[(mč+nč+tč) % 2 == 0]
 341     TMč = č[(mč+nč+tč) % 2 == 1]
 342     return (TEč, TMč)
 343
 344
 345 """
 346 Processing T-matrix related operations from scripts
 347 ===================================================
 348
 349 see also scripts_common.py
 350
 351 The unit cell is defined by a dict particle_specs and a list TMatrix_specs.
 352 This particular module has to provide the T-matrices according to what is defined
 353 in TMatrix_specs
 354
 355 TMatrix_specs is a list of tuples (lMax_override, TMatrix_path, ops)
 356 where
 357   - TMatrix_path is path to the file generated by scuff-tmatrix
 358   - lMax_override is int or None; if it is int and less than the lMax found from the T-matrix file,
 359     lMax_override is used as the order cutoff for the output T-matrix.
 360   - ops is an iterable of tuples (optype, opargs) where currently optype can be 'sym' or 'tr'
 361     for symmetrization operation or some other transformation.
 362 """
 363
 364 #TODO FEATURE more basic group symmetry operations, cf. http://symmetry.otterbein.edu/tutorial/index.html
 365
 366 # This is for finite „fractional“ rotations along the z-axis (mCN means rotation of 2π*(m/N))
 367 reCN = re.compile('(\d*)C(\d+)') # TODO STYLE make this regexp also accept the 3*C_5-type input, eqiv. to 3C5
 368
 369 def get_TMatrix_fromspec(tmatrix_spec):
 370     ''' TODO doc
 371     returns (TMatrices, freqs, lMax)
 372     '''
 373     lMax_override, tmpath, ops = tmatrix_spec
 374     TMatrices, freqs, freqs_weirdunits, lMax = loadScuffTMatrices(tmpath)
 375     if lMax_override is not None and (lMax_override < lMax_orig):
 376         nelem = get_nelem(lMax_override)
 377         TMatrices = TMatrices[...,0:nelem,:,0:nelem]
 378         lMax = lMax_override
 379
 380     for (optype, opargs) in ops:
 381         if optype == 'sym':
 382             mCN = reCN.match(opargs)
 383             if opargs == 'C2' or opargs == 'C_2':
 384                 opmat = apply_matrix_left(yflip_yy(lMax), xflip_yy(lMax), -1)
 385                 TMatrices = (TMatrices + apply_matrix_left(opmat, apply_matrix_left(opmat, TMatrices, -3), -1))/2
 386             elif opargs == 'σ_x' or opargs == 'σx':
 387                 opmat = xflip_tyty(lMax)
 388                 TMatrices = (TMatrices + apply_ndmatrix_left(opmat, apply_ndmatrix_left(opmat, TMatrices, (-4,-3)),(-2,-1)))/2
 389             elif opargs == 'σ_y' or opargs == 'σy':
 390                 opmat = yflip_tyty(lMax)
 391                 TMatrices = (TMatrices + apply_ndmatrix_left(opmat, apply_ndmatrix_left(opmat, TMatrices, (-4,-3)),(-2,-1)))/2
 392             elif opargs == 'σ_z' or opargs == 'σz':
 393                 opmat = zflip_tyty(lMax)
 394                 TMatrices = (TMatrices + apply_ndmatrix_left(opmat, apply_ndmatrix_left(opmat, TMatrices, (-4,-3)),(-2,-1)))/2
 395             elif mCN:
 396                 rotN = int(mCN.group(2)) # the possible m is ignored
 397                 TMatrix_contribs = np.empty((rotN,)+TMatrices.shape, dtype=np.complex_)
 398                 for i in range(rotN):
 399                     rotangle = 2*np.pi*i / rotN
 400                     rot = WignerD_yy_fromvector(lMax, np.array([0,0,rotangle]))
 401                     rotinv = WignerD_yy_fromvector(lMax, np.array([0,0,-rotangle]))
 402                     TMatrix_contribs[i] = apply_matrix_left(rot, apply_matrix_left(rotinv, TMatrices, -3), -1)
 403                 TMatrices = np.sum(TMatrix_contribs, axis=0) / rotN
 404             elif opargs == 'C_inf' or opargs == 'Cinf' or opargs == 'C_∞' or opargs == 'C∞':
 405                 raise ValueError('not implemented: ', opargs)  # TODO FEATURE
 406             else:
 407                 raise ValueError('not implemented: ', opargs)
 408         elif optype == 'tr':
 409             mCN = reCN.match(opargs)
 410             if opargs == 'C2' or opargs == 'C_2':
 411                 opmat = apply_matrix_left(yflip_yy(lMax), xflip_yy(lMax), -1)
 412                 TMatrices = apply_matrix_left(opmat, apply_matrix_left(opmat, TMatrices, -3), -1)/2
 413             elif opargs == 'σ_x' or opargs == 'σx':
 414                 opmat = xflip_tyty(lMax)
 415                 TMatrices = apply_ndmatrix_left(opmat, apply_ndmatrix_left(opmat, TMatrices, (-4,-3)),(-2,-1))/2
 416             elif opargs == 'σ_y' or opargs == 'σy':
 417                 opmat = yflip_tyty(lMax)
 418                 TMatrices = apply_ndmatrix_left(opmat, apply_ndmatrix_left(opmat, TMatrices, (-4,-3)),(-2,-1))/2
 419             elif opargs == 'σ_z' or opargs == 'σz':
 420                 opmat = zflip_tyty(lMax)
 421                 TMatrices = apply_ndmatrix_left(opmat, apply_ndmatrix_left(opmat, TMatrices, (-4,-3)),(-2,-1))/2
 422             elif mCN:
 423                 rotN = int(mCN.group(2))
 424                 power = int(mCN.group(1)) if mCN.group(1) else 1
 425                 rotangle = 2*np.pi*power / rotN
 426                 rot = WignerD_yy_fromvector(lMax, np.array([0,0,rotangle]))
 427                 rotinv = WignerD_yy_fromvector(lMax, np.array([0,0,-rotangle]))
 428                 TMatrices = apply_matrix_left(rot, apply_matrix_left(rotinv, TMatrices, -3), -1)
 429             else:
 430                 raise ValueError('not implemented: ', opargs)
 431         else:
 432             raise ValueError('not implemented: ', optype)
 433     return (TMatrices, freqs, lMax)
 434
 435 class TMatrix(TMatrixSpec):
 436     '''
 437     TODO doc
 438
 439     TODO support for different/multiple interpolators
 440     '''
 441     def __init__(self, tmatrix_spec):
 442         from scipy import interpolate
 443         #self.specification = tmatrix_spec
 444         self.lMax_override = tmatrix_spec.lMax_override
 445         self.tmatrix_path = tmatrix_spec.tmatrix_path
 446         self.ops = tmatrix_spec.ops
 447         self.tmdata, self.freqs, self.lMax = get_TMatrix_fromspec(tmatrix_spec)
 448         self.nelem = get_nelem(self.lMax)
 449         #self._interpolators = dict()
 450         self.default_interpolator = interpolate.interp1d(self.freqs,
 451             self.tmdata, axis=0, kind='linear', fill_value='extrapolate')
 452         self.normalization = NormalizationT.TAYLOR # TODO others are not supported by the loading functions
 453
 454     def atfreq(self, freq):
 455         freqarray = np.array(freq, copy=False)
 456         if freqarray.shape: # not just a scalar
 457             tm_interp = np.empty(freqarray.shape + (2, self.nelem, 2, self.nelem), dtype=np.complex_)
 458             for i in np.ndindex(freqarray.shape):
 459                 tm_interp[i] = self.default_interpolator(freqarray[i])
 460             return tm_interp
 461         else: # scalar
 462             return self.default_interpolator(freq)
 463
 464     __getitem__ = atfreq # might be changed later, use atfreq to be sure
 465
 466 def perform_tmspecs(tmspecs):
 467     """Takes a sequence of TMatrixSpec or TMatrix instances and returns
 468     a list of corresponding TMatrix instances"""
 469     return [(tmspec if hasattr(tmspec, "tmdata") else TMatrix(tmspec))
 470                 for tmspec in tmspecs]
 471