tests/test_vswf_transop_consistency.py

   1 #!/usr/bin/env python3
   2 """
   3 Tests whether direct evaluation of VSWFs gives the same result as their
   4 representation in terms of translation operators and regular electric dipole
   5 waves at origin
   6 """
   7
   8 from qpms import Particle, CTMatrix, lorentz_drude, EpsMuGenerator, TMatrixGenerator, BaseSpec, FinitePointGroup, ScatteringSystem, TMatrixInterpolator, EpsMu, dbgmsg_enable, dbgmsg_disable, dbgmsg_active, BesselType,eV, hbar
   9 from qpms.qpms_c import set_gsl_pythonic_error_handling
  10 from qpms.symmetries import point_group_info
  11 import numpy as np
  12 eh = eV/hbar
  13 np.random.seed(666)
  14 dbgmsg_enable(2)
  15
  16 part_radius = 80e-9
  17 p = 1580e-9
  18
  19 set_gsl_pythonic_error_handling()
  20
  21 sym = FinitePointGroup(point_group_info['D4h'])
  22 bspec1 = BaseSpec(lMax=3)
  23 medium=EpsMuGenerator(EpsMu(1.52**2))
  24 t1 = TMatrixGenerator.sphere(medium, EpsMuGenerator(lorentz_drude['Au']), r=part_radius)
  25 p1 = Particle((0,0,0),t1,bspec=bspec1)
  26 ss, ssw = ScatteringSystem.create([p1], EpsMuGenerator(EpsMu(1.52**2)), 1.4*eh, sym)
  27
  28
  29 points = np.random.random((100,3)) * p
  30 points = points[np.linalg.norm(points, axis=-1) > part_radius]
  31 t,l,m = bspec1.tlm()
  32
  33 fails=0
  34
  35 for i in range(ss.fecv_size):
  36     fvc = np.zeros((ss.fecv_size,), dtype=complex)
  37     fvc[i] = 1
  38
  39     E = ssw.scattered_E(fvc, points)
  40     E_alt = ssw.scattered_E(fvc, points,alt=True)
  41     diff = abs(E-E_alt)
  42     reldiffavg = np.nanmean(diff/(abs(E)+abs(E_alt)))
  43     fail = reldiffavg > 1e-3
  44     fails += fail
  45
  46     print('E' if t[i] == 2 else 'M', l[i], m[i], np.amax(diff), reldiffavg, 'FAIL!' if fail else 'OK')
  47
  48 exit(fails)