diff --git a/qpms/lattices.py b/qpms/lattices.py index 0e45436..960fc24 100644 --- a/qpms/lattices.py +++ b/qpms/lattices.py @@ -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