oldtests/modeproblem_implementations.py

   1 from qpms import Particle, CTMatrix, BaseSpec, FinitePointGroup, ScatteringSystem, TMatrixInterpolator, eV, hbar, c
   2 from qpms.symmetries import point_group_info
   3 import numpy as np
   4 import os
   5 import sys
   6 nm = 1e-9
   7
   8 sym = FinitePointGroup(point_group_info['D2h'])
   9 bspec = BaseSpec(lMax = 2)
  10 #tmfile = '/m/phys/project/qd/Marek/tmatrix-experiments/Cylinder/AaroBEC/cylinder_50nm_lMax4_cleaned.TMatrix'
  11 tmfile = '/home/mmn/repo/tmatrix-experiments/Cylinder/AaroBEC/cylinder_50nm_lMax4_cleaned.TMatrix'
  12 #outputdatadir = '/home/necadam1/wrkdir/AaroBECfinite_new'
  13 outputdatadir = '/tmp/u/46/necadam1/unix/project/AaroBECfinite_new'
  14 os.makedirs(outputdatadir, exist_ok = True)
  15 interp = TMatrixInterpolator(tmfile, bspec, symmetrise = sym, atol = 1e-8)
  16 # There is only one t-matrix in the system for each frequency. We initialize the matrix with the lowest frequency data.
  17 # Later, we can replace it using the tmatrix[...] = interp(freq) and s.update_tmatrices NOT YET; TODO
  18
  19 omega = 1.475 * eV/hbar
  20 sv_threshold = 0.5
  21
  22 # Now place the particles and set background index.
  23 px = 571*nm; py = 621*nm
  24 n = 1.51
  25 Nx = 5
  26 Ny = 7
  27
  28 orig_x = (np.arange(Nx/2) + (0 if (Nx % 2) else .5)) * px
  29 orig_y = (np.arange(Ny/2) + (0 if (Ny % 2) else .5)) * py
  30
  31 orig_xy = np.stack(np.meshgrid(orig_x, orig_y), axis = -1)
  32
  33 tmatrix = interp(omega)
  34 particles = [Particle(orig_xy[i], tmatrix) for i in np.ndindex(orig_xy.shape[:-1])]
  35
  36
  37 ss = ScatteringSystem(particles, sym)
  38 k = n * omega / c
  39
  40 for iri in range(ss.nirreps):
  41     mm_iri_orig = ss.modeproblem_matrix_packed(k, iri, version = None)
  42     mm_iri_alt = ss.modeproblem_matrix_packed(k, iri, version='R')
  43     mm_iri_paral = ss.modeproblem_matrix_packed(k, iri, version='pR')
  44     print(np.amax(abs(mm_iri_orig-mm_iri_alt)), np.amax(abs(mm_iri_orig-mm_iri_paral)))
  45