From d2653e6a6d59a3687ce62428ea8c41de6899c75a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 10 May 2017 11:44:04 +0300 Subject: [PATCH] New symz_indexarray function, fix hexlatice_zsym_getSVD Former-commit-id: 2d42d93b8c7d2e7596127293c622197794fe7f20 --- misc/dispersion2D-SVD.py | 15 --------------- qpms/hexpoints.py | 22 ++++++++++++++-------- qpms/qpms_p.py | 38 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 52 insertions(+), 23 deletions(-) diff --git a/misc/dispersion2D-SVD.py b/misc/dispersion2D-SVD.py index b76c4cc..8bc93f9 100755 --- a/misc/dispersion2D-SVD.py +++ b/misc/dispersion2D-SVD.py @@ -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 diff --git a/qpms/hexpoints.py b/qpms/hexpoints.py index d53f7b7..03aba3a 100644 --- a/qpms/hexpoints.py +++ b/qpms/hexpoints.py @@ -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): @@ -473,7 +478,8 @@ 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) diff --git a/qpms/qpms_p.py b/qpms/qpms_p.py index 714c8de..d852a6c 100644 --- a/qpms/qpms_p.py +++ b/qpms/qpms_p.py @@ -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č) +