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

Fixes 

Former-commit-id: 332fdd586e784b734f245b0d149cc63d6a2a4d68
This commit is contained in:
Marek Nečada 2016-12-22 09:14:58 +02:00
parent 9e390e0d30
commit 1f995be463
1 changed files with 11 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
qpms/lattices.py
View File

@ -200,12 +200,14 @@ class Scattering_2D_zsym(Scattering):
self.TE_NMz = (self.my + self.ny) % 2
self.TM_NMz = 1 - self.TE_NMz
# TODO možnost zadávat T-matice rovnou ve zhuštěné podobě
TMatrices_TE = TMatrices[...,TE_NMz[:,nx],TE_yz[:,nx],TE_NMZ[nx,:],TE_yz[nx,:]]
TMatrices_TM = TMatrices[...,TM_NMz[:,nx],TM_yz[:,nx],TM_NMZ[nx,:],TM_yz[nx,:]]
TMatrices_TE = TMatrices[...,self.TE_NMz[:,nx],self.TE_yz[:,nx],self.TE_NMz[nx,:],self.TE_yz[nx,:]]
TMatrices_TM = TMatrices[...,self.TM_NMz[:,nx],self.TM_yz[:,nx],self.TM_NMz[nx,:],self.TM_yz[nx,:]]
self.TMatrices_TE = np.broadcast_to(TMatrices_TE, (self.N, self.nelem, self.nelem))
self.TMatrices_TM = np.broadcast_to(TMatrices_TM, (self.N, self.nelem, self.nelem))
self.prepared_TE = False
self.prepared_TM = False
self.interaction_matrix_TE = None
self.interaction_matrix_TM= None
def prepare_partial(self, TE_or_TM, keep_interaction_matrix = False, verbose=False):
'''
@ -219,7 +221,7 @@ class Scattering_2D_zsym(Scattering):
self.lupiv_TE = scipy.linalg.lu_factor(self.interaction_matrix_TE, overwrite_a = not keep_interaction_matrix)
self.prepared_TE = True
if (TE_or_TM == 1): #TM
in not self.prepared_TM:
if not self.prepared_TM:
if not self.interaction_matrix_TM:
self.build_interaction_matrix(1, verbose)
self.lupiv_TM = scipy.linalg.lu_factor(self.interaction_matrix_TM, overwrite_a = not keep_interaction_matrix)
@ -261,11 +263,12 @@ class Scattering_2D_zsym(Scattering):
tr = B̃(my[yj],ny[yj],my[yi],ny[yi],kdlj=d_i2j[0]*self.k_0,θlj=d_i2j[1],φlj=d_i2j[2],r_ge_d=False,J=self.J_scat)
leftmatrix[j,yj,i,yi] = tr
leftmatrix[i,yi,j,yj] = tr if (0 == (my[yj]+my[yi]) % 2) else -tr
_ time_e(sbtime, verbose, step = 'Calculating interparticle translation coefficients, T%s part' % ('M' if EoM else 'E'))
_time_e(sbtime, verbose, step = 'Calculating interparticle translation coefficients, T%s part' % ('M' if EoM else 'E'))
for j in range(N):
leftmatrix = - np.tensordot(self.TMatrices_TM[j] if EoM else self.TMatrices_TE[j],leftmatrix[j],
leftmatrix[j] = - np.tensordot(self.TMatrices_TM[j] if EoM else self.TMatrices_TE[j],leftmatrix[j],
axes = ([-1],[0]))
leftmatrix[j,:,j,:] += idm
leftmatrix.shape = (self.N*self.nelem, self.N*self.nelem)
if EoM == 0:
self.interaction_matrix_TE = leftmatrix
if EoM == 1:
@ -285,9 +288,11 @@ class Scattering_2D_zsym(Scattering):
MP_0[j] = np.tensordot(self.TMatrices_TM[j], pq_0[j], axes=([-1],[-1]))
else: #TE
MP_0[j] = np.tensordot(self.TMatrices_TE[j], pq_0[j], axes=([-1],[-1]))
MP_0.shape = (self.N*self.nelem,)
solvebtime = _time_b(verbose,step='Solving the linear equation')
ab = scipy.linalg.lu_solve(self.lupiv_TM if TE_or_EM else self.lupiv.TE, MP_0)
ab = scipy.linalg.lu_solve(self.lupiv_TM if TE_or_TM else self.lupiv_TE, MP_0)
_time_e(solvebtime, verbose, step='Solving the linear equation')
ab.shape = (self.N, self.nelem)
_time_e(btime,verbose)
return ab