misc/finiterectlat-modes.py

   1 #!/usr/bin/env python3
   2
   3 import math
   4 from qpms.argproc import ArgParser, annotate_pdf_metadata
   5
   6
   7 ap = ArgParser(['rectlattice2d_finite', 'background_analytical', 'single_particle', 'single_lMax', ])
   8
   9 ap.add_argument("-t", "--rank-tolerance", type=float, default=1e11)
  10 ap.add_argument("-c", "--min-candidates", type=int, default=1, help='always try at least this many eigenvalue candidates, even if their SVs in the rank tests are lower than rank_tolerance')
  11 ap.add_argument("-T", "--residual-tolerance", type=float, default=0.)
  12
  13 ap.add_argument("-o", "--output", type=str, required=False, help='output path (if not provided, will be generated automatically)')
  14 #ap.add_argument("-O", "--plot-out", type=str, required=False, help="path to plot output (optional)")
  15 #ap.add_argument("-P", "--plot", action='store_true', help="if -p not given, plot to a default path")
  16 #ap.add_argument("-g", "--save-gradually", action='store_true', help="saves the partial result after computing each irrep")
  17
  18 ap.add_argument("-f", "--centre", type=complex, required=True, help='Contour centre in eV')
  19 ap.add_argument("--ai", type=float, default=0.05, help="Contour imaginary half-axis in eV")
  20 ap.add_argument("--ar", type=float, default=0.05, help="Contour real half-axis in eV")
  21 ap.add_argument("-N", type=int, default="150", help="Integration contour discretisation size")
  22 ap.add_argument("--D2", action='store_true', help="Use D2h symmetry even if the array has square symmetry")
  23
  24 ap.add_argument("--irrep", type=str, default="none", help="Irrep subspace (irrep index from 0 to 7 (9 for D4h), irrep label, or 'none' for no irrep decomposition")
  25
  26 a=ap.parse_args()
  27
  28 import logging
  29 logging.basicConfig(format='%(asctime)s %(message)s', level=logging.INFO)
  30
  31 Nx, Ny = a.size
  32 px, py = a.period
  33
  34 thegroup = 'D4h' if px == py and Nx == Ny and not a.D2 else 'D2h'
  35
  36 particlestr = ("sph" if a.height is None else "cyl") + ("_r%gnm" % (a.radius*1e9))
  37 if a.height is not None: particlestr += "_h%gnm" % (a.height * 1e9)
  38 defaultprefix = "%s_p%gnmx%gnm_%dx%d_m%s_B%s_L%d_c(%s±%g±%gj)eV_cn%d_%s" % (
  39     particlestr, px*1e9, py*1e9, Nx, Ny, str(a.material), str(a.background), a.lMax,
  40     str(a.centre), a.ar, a.ai, a.N,
  41     thegroup,
  42     )
  43 logging.info("Default file prefix: %s" % defaultprefix)
  44
  45 def inside_ellipse(point_xy, centre_xy, halfaxes_xy):
  46     x = point_xy[0] - centre_xy[0]
  47     y = point_xy[1] - centre_xy[1]
  48     ax = halfaxes_xy[0]
  49     ay = halfaxes_xy[1]
  50     return ((x/ax)**2 + (y/ay)**2) <= 1
  51
  52
  53 import numpy as np
  54 import qpms
  55 from qpms.cybspec import BaseSpec
  56 from qpms.cytmatrices import CTMatrix, TMatrixGenerator
  57 from qpms.qpms_c import Particle
  58 from qpms.cymaterials import EpsMu, EpsMuGenerator, LorentzDrudeModel, lorentz_drude
  59 from qpms.cycommon import DebugFlags, dbgmsg_enable
  60 from qpms import FinitePointGroup, ScatteringSystem, BesselType, eV, hbar
  61 from qpms.symmetries import point_group_info
  62 eh = eV/hbar
  63
  64 dbgmsg_enable(DebugFlags.INTEGRATION)
  65
  66 #Particle positions
  67 orig_x = (np.arange(Nx/2) + (0 if (Nx % 2) else .5)) * px
  68 orig_y = (np.arange(Ny/2) + (0 if (Ny % 2) else .5)) * py
  69
  70 orig_xy = np.stack(np.meshgrid(orig_x, orig_y), axis = -1)
  71
  72
  73 bspec = BaseSpec(lMax = a.lMax)
  74 medium = EpsMuGenerator(ap.background_epsmu)
  75 particles= [Particle(orig_xy[i], ap.tmgen, bspec) for i in np.ndindex(orig_xy.shape[:-1])]
  76
  77 sym = FinitePointGroup(point_group_info[thegroup])
  78 logging.info("Creating scattering system object")
  79 ss, ssw = ScatteringSystem.create(particles, medium, a.centre * eh, sym=sym)
  80
  81 if a.irrep == 'none':
  82     iri = None # no irrep decomposition
  83     irname = 'full'
  84     logging.info("Not performing irrep decomposition and working with the full space of dimension %d." % ss.fecv_size)
  85 else:
  86     try:
  87         iri = int(a.irrep)
  88     except ValueError:
  89         iri = ss.irrep_names.index(a.irrep)
  90     irname = ss.irrep_names[iri]
  91     logging.info("Using irrep subspace %s (iri = %d) of dimension %d." % (irname, iri, ss.saecv_sizes[iri]))
  92
  93 outfile_tmp = defaultprefix + ".tmp" if a.output is None else a.output + ".tmp"
  94
  95 logging.info("Starting Beyn's algorithm")
  96 results = ss.find_modes(iri=iri, omega_centre = a.centre*eh, omega_rr=a.ar*eh, omega_ri=a.ai*eh, contour_points=a.N,
  97         rank_tol=a.rank_tolerance, rank_min_sel=a.min_candidates, res_tol=a.residual_tolerance)
  98 results['inside_contour'] = inside_ellipse((results['eigval'].real, results['eigval'].imag),
  99         (a.centre.real*eh, a.centre.imag*eh), (a.ar*eh, a.ai*eh))
 100 results['refractive_index_internal'] = [medium(om).n for om in results['eigval']]
 101
 102 outfile = defaultprefix + (('_ir%s_%s.npz' % (str(iri), irname)) if iri is not None else '.npz') if a.output is None else a.output
 103 np.savez(outfile, meta={**vars(a), 'qpms_version' : qpms.__version__()}, **results)
 104 logging.info("Saved to %s" % outfile)
 105
 106 exit(0)
 107
 108 # TODO plots.
 109 if a.plot or (a.plot_out is not None):
 110     import matplotlib
 111     matplotlib.use('pdf')
 112     from matplotlib import pyplot as plt
 113     from matplotlib.backends.backend_pdf import PdfPages
 114     fig = plt.figure()
 115     ax = fig.add_subplot(111)
 116     ax.plot(sinalpha_list, σ_ext*1e12,label='$\sigma_\mathrm{ext}$')
 117     ax.plot(sinalpha_list, σ_scat*1e12, label='$\sigma_\mathrm{scat}$')
 118     ax.plot(sinalpha_list, σ_abs*1e12, label='$\sigma_\mathrm{abs}$')
 119     ax.legend()
 120     ax.set_xlabel('$\sin\\alpha$')
 121     ax.set_ylabel('$\sigma/\mathrm{\mu m^2}$')
 122
 123     plotfile = defaultprefix + ".pdf" if a.plot_out is None else a.plot_out
 124     with PdfPages(plotfile) as pdf:
 125         pdf.savefig(fig)
 126         annotate_pdf_metadata(pdf, scriptname='finiterectlat-modes.py')
 127
 128 exit(0)
 129