oldtests/scatsys.py

   1 from qpms import Particle, CTMatrix, BaseSpec, FinitePointGroup, ScatteringSystem
   2 from qpms.symmetries import point_group_info
   3 import numpy as np
   4
   5 np.random.seed(444)
   6
   7 sym = FinitePointGroup(point_group_info['D3h'])
   8 bspec2 = BaseSpec(lMax=2)
   9 bspec1 = BaseSpec(lMax=1)
  10 t1 = CTMatrix(bspec1, np.diag(np.random.random(len(bspec1))))
  11 t2 = CTMatrix(bspec2, np.diag(np.random.random(len(bspec2))))
  12 p1 = Particle((1,2,),t1)
  13 p2 = Particle((1,2,3),t1)
  14 p3 = Particle((0.1,2),t2)
  15 ss = ScatteringSystem([p1, p2, p3], sym)
  16
  17 #print(ss.fecv_size, ss.saecv_sizes, ss.nirreps, ss.nirrep_names)
  18
  19 fullvector = np.random.rand(ss.fecv_size) + 1j*np.random.rand(ss.fecv_size)
  20 packedvectors = [(iri, ss.pack_vector(fullvector, iri)) for iri in range(ss.nirreps)]
  21 unpackedvectors = np.array([ss.unpack_vector(v[1], v[0]) for v in packedvectors])
  22 rec_fullvector = np.sum(unpackedvectors, axis=0)
  23 thediff = np.amax(abs(rec_fullvector-fullvector))
  24 assert(thediff < 1e-8)
  25
  26 packedmatrices = list()
  27 for iri in range(ss.nirreps):
  28     d = ss.saecv_sizes[iri]
  29     m = np.random.rand(d,d)+1j*np.random.rand(d,d)
  30     packedmatrices.append((iri,m))
  31
  32 fullmatrix = np.zeros((ss.fecv_size, ss.fecv_size), dtype=complex)
  33 for iri, m in packedmatrices:
  34     fullmatrix += ss.unpack_matrix(m, iri)
  35
  36 for iri, m in packedmatrices:
  37     print (m.shape)
  38     repackedmatrix = ss.pack_matrix(fullmatrix,iri)
  39     print(np.amax(abs(repackedmatrix-m)))
  40
  41 k = 1.7
  42 modematrix_full = ss.modeproblem_matrix_full(k)
  43 modematrix_packed_list = [(iri, ss.pack_matrix(modematrix_full,iri)) for iri in range(ss.nirreps)]
  44 modematrix_full_rec = np.empty((ss.fecv_size, ss.fecv_size), dtype=complex)
  45 for iri, m in modematrix_packed_list:
  46     modematrix_full_rec += ss.unpack_matrix(m,iri)
  47 print(np.amax(abs(modematrix_full-modematrix_full_rec)))