QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

New symz_indexarray function, fix hexlatice_zsym_getSVD 

Former-commit-id: 2d42d93b8c7d2e7596127293c622197794fe7f20
This commit is contained in:
Marek Nečada 2017-05-10 11:44:04 +03:00
parent 3a1f7c95aa
commit d2653e6a6d
3 changed files with 52 additions and 23 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
misc/dispersion2D-SVD.py
View File

@ -139,21 +139,6 @@ nelem = len(my)
if pargs.lMax: #force commandline specified lMax
TMatrices_orig = TMatrices_orig[...,0:nelem,:,0:nelem]
ž = np.arange(2*nelem)
tž = ž // nelem
mž = my[ž%nelem]
nž = ny[ž%nelem]
TEž = ž[(mž+nž+tž) % 2 == 0]
TMž = ž[(mž+nž+tž) % 2 == 1]
č = np.arange(2*2*nelem)
žč = č % (2* nelem)
tč = tž[žč]
mč = mž[žč]
nč = nž[žč]
TEč = č[(mč+nč+tč) % 2 == 0]
TMč = č[(mč+nč+tč) % 2 == 1]
TMatrices = np.array(np.broadcast_to(TMatrices_orig[:,nx,:,:,:,:],(len(freqs_orig),2,2,nelem,2,nelem)) )
#TMatrices[:,:,:,:,:,ny==3] *= factor13inc

20
qpms/hexpoints.py
View File

@ -436,8 +436,11 @@ def hexlattice_get_AB(lMax, k_hexside, maxlayer, circular=True, return_points =
return d
from scipy.constants import c
from .timetrack import _time_b, _time_e
from .qpms_p import symz_indexarrays
def hexlattice_zsym_getSVD(lMax, TMatrices_om, epsilon_b, hexside, maxlayer, omega, klist, gaussianSigma=False, onlyNmin=0):
def hexlattice_zsym_getSVD(lMax, TMatrices_om, epsilon_b, hexside, maxlayer, omega, klist, gaussianSigma=False, onlyNmin=0, verbose=False):
btime = _time_b(verbose)
nelem = lMax * (lMax + 2)
n2id = np.identity(2*nelem)
n2id.shape = (2,nelem,2,nelem)
@ -455,9 +458,11 @@ def hexlattice_zsym_getSVD(lMax, TMatrices_om, epsilon_b, hexside, maxlayer, ome
d2u_translations = tdic['d2u_tr']*hexside*_s3
if gaussianSigma:
sbtime = _time_b(verbose, step='Calculating gaussian envelope')
unitcell_envelope = np.exp(-np.sum(tdic['self_tr']**2,axis=-1)/(2*gaussianSigma**2))
u2d_envelope = np.exp(-np.sum(tdic['u2d_tr']**2,axis=-1)/(2*gaussianSigma**2))
d2u_envelope = np.exp(-np.sum(tdic['d2u_tr']**2,axis=-1)/(2*gaussianSigma**2))
_time_e(sbtime, verbose, step='Calculating gaussian envelope')
#TMatrices_om = TMatrices_interp(omega)
if(not onlyNmin):
@ -474,6 +479,7 @@ def hexlattice_zsym_getSVD(lMax, TMatrices_om, epsilon_b, hexside, maxlayer, ome
leftmatrixlist = np.full((klist.shape[0],2,2,nelem,2,2,nelem),np.nan,dtype=complex)
isNaNlist = np.zeros((klist.shape[0]), dtype=bool)
sbtime = _time_b(verbose, step='Initializing matrices for SVD for a given list of k\'s.')
# sem nějaká rozumná smyčka
for ki in range(klist.shape[0]):
k = klist[ki]
@ -516,21 +522,21 @@ def hexlattice_zsym_getSVD(lMax, TMatrices_om, epsilon_b, hexside, maxlayer, ome
leftmatrixlist[ki] = leftmatrix
nnlist = np.logical_not(isNaNlist)
leftmatrixlist_s = np.reshape(leftmatrixlist,(klist.shape[0], 2*2*nelem,2*2*nelem))[nnlist]
TEč, TMč = symz_indexarrays(lMax, 2)
leftmatrixlist_TE = leftmatrixlist_s[np.ix_(np.arange(leftmatrixlist_s.shape[0]),TEč,TEč)]
leftmatrixlist_TM = leftmatrixlist_s[np.ix_(np.arange(leftmatrixlist_s.shape[0]),TMč,TMč)]
_time_e(sbtime, verbose, step='Initializing matrices for SVD for a given list of k\'s.')
sbtime = _time_b(verbose, step='Calculating SVDs for a given list of k\'s.')
if(not onlyNmin):
svUfullTElist[nnlist], svSfullTElist[nnlist], svVfullTElist[nnlist] = np.linalg.svd(leftmatrixlist_TE, compute_uv=True)
svUfullTMlist[nnlist], svSfullTMlist[nnlist], svVfullTMlist[nnlist] = np.linalg.svd(leftmatrixlist_TM, compute_uv=True)
_time_e(sbtime, verbose, step='Calculating SVDs for a given list of k\'s.')
return ((svUfullTElist, svSfullTElist, svVfullTElist), (svUfullTMlist, svSfullTMlist, svVfullTMlist))
else:
minsvTElist[nnlist] = np.linalg.svd(leftmatrixlist_TE, compute_uv=False)[...,-onlyNmin:]
minsvTMlist[nnlist] = np.linalg.svd(leftmatrixlist_TM, compute_uv=False)[...,-onlyNmin:]
_time_e(sbtime, verbose, step='Calculating SVDs for a given list of k\'s.')
return (minsvTElist, minsvTMlist)

38
qpms/qpms_p.py
View File

@ -930,3 +930,41 @@ def apply_ndmatrix_left(matrix,tensor,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 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.
"""
nelem = lMax * (lMax + 2)
ž = 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č)