From 3458acca16120b52810c16b92611283bd5c1a9e8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Sun, 5 Apr 2020 21:41:15 +0300 Subject: [PATCH] Infinite rect lattices script multiple freqs etc. Former-commit-id: 8db37dfaf051abf1d6a542eb6ca9b40317848469 --- misc/infiniterectlat-scatter.py | 206 ++++++++++++++++++++++++-------- qpms/argproc.py | 22 ++++ 2 files changed, 175 insertions(+), 53 deletions(-) diff --git a/misc/infiniterectlat-scatter.py b/misc/infiniterectlat-scatter.py index 72c6b57..b7111b6 100755 --- a/misc/infiniterectlat-scatter.py +++ b/misc/infiniterectlat-scatter.py @@ -1,13 +1,11 @@ #!/usr/bin/env python3 import math +pi = math.pi from qpms.argproc import ArgParser -ap = ArgParser(['rectlattice2d', 'single_particle', 'single_lMax', 'single_omega']) -ap.add_argument("-k", '--kx-lim', nargs=2, type=float, required=True, help='k vector', metavar=('KX_MIN', 'KX_MAX')) -# ap.add_argument("--kpi", action='store_true', help="Indicates that the k vector is given in natural units instead of SI, i.e. the arguments given by -k shall be automatically multiplied by pi / period (given by -p argument)") +ap = ArgParser(['rectlattice2d', 'single_particle', 'single_lMax', 'omega_seq_real_ng', 'planewave']) ap.add_argument("-o", "--output", type=str, required=False, help='output path (if not provided, will be generated automatically)') -ap.add_argument("-N", type=int, default="151", help="Number of angles") ap.add_argument("-O", "--plot-out", type=str, required=False, help="path to plot output (optional)") ap.add_argument("-P", "--plot", action='store_true', help="if -p not given, plot to a default path") #ap.add_argument("-g", "--save-gradually", action='store_true', help="saves the partial result after computing each irrep") @@ -18,17 +16,9 @@ a=ap.parse_args() import logging logging.basicConfig(format='%(asctime)s %(message)s', level=logging.INFO) -px, py = a.period - -particlestr = ("sph" if a.height is None else "cyl") + ("_r%gnm" % (a.radius*1e9)) -if a.height is not None: particlestr += "_h%gnm" % (a.height * 1e9) -defaultprefix = "%s_p%gnmx%gnm_m%s_n%g_angles(%g_%g)_Ey_f%geV_L%d_cn%d" % ( - particlestr, px*1e9, py*1e9, str(a.material), a.refractive_index, a.kx_lim[0], a.kx_lim[1], a.eV, a.lMax, a.N) -logging.info("Default file prefix: %s" % defaultprefix) - - import numpy as np import qpms +from qpms.qpms_p import cart2sph, sph2cart, sph_loccart2cart, sph_loccart_basis import warnings from qpms.cybspec import BaseSpec from qpms.cytmatrices import CTMatrix, TMatrixGenerator @@ -40,6 +30,15 @@ eh = eV/hbar dbgmsg_enable(DebugFlags.INTEGRATION) +px, py = a.period + +particlestr = ("sph" if a.height is None else "cyl") + ("_r%gnm" % (a.radius*1e9)) +if a.height is not None: particlestr += "_h%gnm" % (a.height * 1e9) +defaultprefix = "%s_p%gnmx%gnm_m%s_n%g_φ%g_θ(%g_%g)π_ψ%gπ_χ%gπ_f%s_L%d" % ( + particlestr, px*1e9, py*1e9, str(a.material), a.refractive_index, a.phi/pi, np.amin(a.theta)/pi, np.amax(a.theta)/pi, a.psi/pi, a.chi/pi, ap.omega_descr, a.lMax) +logging.info("Default file prefix: %s" % defaultprefix) + + a1 = ap.direct_basis[0] a2 = ap.direct_basis[1] @@ -48,72 +47,173 @@ orig_x = [0] orig_y = [0] orig_xy = np.stack(np.meshgrid(orig_x,orig_y),axis=-1) -omega = ap.omega - bspec = BaseSpec(lMax = a.lMax) # The parameters here should probably be changed (needs a better qpms_c.Particle implementation) pp = Particle(orig_xy[0][0], ap.tmgen, bspec=bspec) par = [pp] -ss, ssw = ScatteringSystem.create(par, ap.background_emg, omega, latticebasis = ap.direct_basis) +ss, ssw = ScatteringSystem.create(par, ap.background_emg, ap.allomegas[0], latticebasis = ap.direct_basis) -if ssw.wavenumber.imag != 0: - warnings.warn("The background medium wavenumber has non-zero imaginary part. Don't expect meaningful results for cross sections.") -wavenumber = ssw.wavenumber.real -sinalpha_list = np.linspace(a.kx_lim[0],a.kx_lim[1],a.N) +## Plane wave data +a.theta = np.array(a.theta) +dir_sph_list = np.stack((np.broadcast_to(1, a.theta.shape), a.theta, np.broadcast_to(a.phi, a.theta.shape)), axis=-1) +sψ, cψ = math.sin(a.psi), math.cos(a.psi) +sχ, cχ = math.sin(a.chi), math.cos(a.chi) +E_sph = (0., cψ*cχ + 1j*sψ*sχ, sψ*cχ + 1j*cψ*sχ) -# Plane wave data -E_cart_list = np.empty((a.N,3), dtype=complex) -E_cart_list[:,:] = np.array((0,1,0))[None,:] -k_cart_list = np.empty((a.N,3), dtype=float) -k_cart_list[:,0] = sinalpha_list -k_cart_list[:,1] = 0 -k_cart_list[:,2] = np.sqrt(1-sinalpha_list**2) -k_cart_list *= wavenumber +dir_cart_list = sph2cart(dir_sph_list) +E_cart_list = sph_loccart2cart(E_sph, dir_sph_list) -σ_ext_list = np.empty((a.N,), dtype=float) -σ_scat_list = np.empty((a.N,), dtype=float) +nfreq = len(ap.allomegas) +ndir = len(dir_sph_list) + +k_cart_arr = np.empty((nfreq, ndir, 3), dtype=float) +wavenumbers = np.empty((nfreq,), dtype=float) + +σ_ext_arr = np.empty((nfreq,ndir), dtype=float) +σ_scat_arr = np.empty((nfreq,ndir), dtype=float) with pgsl_ignore_error(15): # avoid gsl crashing on underflow - for j in range(a.N): - k_cart = k_cart_list[j] - blochvector = (k_cart[0], k_cart[1], 0) - # the following two could be calculated only once, but probably not a big deal - LU = ssw.scatter_solver(k=blochvector) - ã = ss.planewave_full(k_cart=k_cart, E_cart=E_cart_list[j]) - Tã = ssw.apply_Tmatrices_full(ã) - f = LU(Tã) + for i, omega in enumerate(ap.allomegas): + if i != 0: + ssw = ss(omega) + if ssw.wavenumber.imag != 0: + warnings.warn("The background medium wavenumber has non-zero imaginary part. Don't expect meaningful results for cross sections.") + wavenumber = ssw.wavenumber.real + wavenumbers[i] = wavenumber - σ_ext_list[j] = -np.vdot(ã, f).real/wavenumber**2 - translation_matrix = ssw.translation_matrix_full(blochvector=blochvector) + np.eye(ss.fecv_size) - σ_scat_list[j] = np.vdot(f,np.dot(translation_matrix, f)).real/wavenumber**2 + k_sph_list = np.array(dir_sph_list, copy=True) + k_sph_list[:,0] = wavenumber -σ_abs_list = σ_ext_list - σ_scat_list + k_cart_arr[i] = sph2cart(k_sph_list) + + + for j in range(ndir): + k_cart = k_cart_arr[i,j] + blochvector = (k_cart[0], k_cart[1], 0) + # the following two could be calculated only once, but probably not a big deal + LU = ssw.scatter_solver(k=blochvector) + ã = ss.planewave_full(k_cart=k_cart, E_cart=E_cart_list[j]) + Tã = ssw.apply_Tmatrices_full(ã) + f = LU(Tã) + + σ_ext_arr[i,j] = -np.vdot(ã, f).real/wavenumber**2 + translation_matrix = ssw.translation_matrix_full(blochvector=blochvector) + np.eye(ss.fecv_size) + σ_scat_arr[i,j] = np.vdot(f,np.dot(translation_matrix, f)).real/wavenumber**2 + +σ_abs_arr = σ_ext_arr - σ_scat_arr outfile = defaultprefix + ".npz" if a.output is None else a.output -np.savez(outfile, meta=vars(a), sinalpha=sinalpha_list, k_cart = k_cart_list, E_cart=E_cart_list, σ_ext=σ_ext_list,σ_abs=σ_abs_list,σ_scat=σ_scat_list, omega=omega, wavenumber=wavenumber, unitcell_area=ss.unitcell_volume +np.savez(outfile, meta=vars(a), dir_sph=dir_sph_list, k_cart = k_cart_arr, omega = ap.allomegas, E_cart = E_cart_list, wavenumbers= wavenumbers, σ_ext=σ_ext_arr,σ_abs=σ_abs_arr,σ_scat=σ_scat_arr, unitcell_area=ss.unitcell_volume ) logging.info("Saved to %s" % outfile) if a.plot or (a.plot_out is not None): import matplotlib + from matplotlib.backends.backend_pdf import PdfPages matplotlib.use('pdf') from matplotlib import pyplot as plt + from scipy.interpolate import griddata - fig = plt.figure() - ax = fig.add_subplot(111) - ax.plot(sinalpha_list, σ_ext_list*1e12,label='$\sigma_\mathrm{ext}$') - ax.plot(sinalpha_list, σ_scat_list*1e12, label='$\sigma_\mathrm{scat}$') - ax.plot(sinalpha_list, σ_abs_list*1e12, label='$\sigma_\mathrm{abs}$') - ax.legend() - ax.set_xlabel('$\sin\\alpha$') - ax.set_ylabel('$\sigma/\mathrm{\mu m^2}$') - plotfile = defaultprefix + ".pdf" if a.plot_out is None else a.plot_out - fig.savefig(plotfile) + with PdfPages(plotfile) as pdf: + ipm = 'nearest' + sintheta = np.sin(a.theta) + if False: #len(ap.omega_ranges) != 0: + # angle plot --------------------------------- + fig = plt.figure(figsize=(210/25.4, 297/25.4)) + vmax = max(np.amax(σ_ext_arr), np.amax(σ_scat_arr), np.amax(σ_abs_arr)) + vmin = min(np.amin(σ_ext_arr), np.amin(σ_scat_arr), np.amin(σ_abs_arr)) + + ax = fig.add_subplot(311) + ax.pcolormesh(a.theta, ap.allomegas/eh, σ_ext_arr, vmin=vmin, vmax=vmax) + ax.set_xlabel('$\\theta$') + ax.set_ylabel('$\\hbar\\omega / \\mathrm{eV}$') + ax.set_title('$\\sigma_\\mathrm{ext}$') + + ax = fig.add_subplot(312) + ax.pcolormesh(a.theta, ap.allomegas/eh, σ_scat_arr, vmin=vmin, vmax=vmax) + ax.set_xlabel('$\\theta$') + ax.set_ylabel('$\\hbar\\omega / \\mathrm{eV}$') + ax.set_title('$\\sigma_\\mathrm{scat}$') + + ax = fig.add_subplot(313) + im = ax.pcolormesh(a.theta, ap.allomegas/eh, σ_abs_arr, vmin=vmin, vmax=vmax) + ax.set_xlabel('$\\theta$') + ax.set_ylabel('$\\hbar\\omega / \\mathrm{eV}$') + ax.set_title('$\\sigma_\\mathrm{abs}$') + + + fig.subplots_adjust(right=0.8) + fig.colorbar(im, cax = fig.add_axes([0.85, 0.15, 0.05, 0.7])) + + pdf.savefig(fig) + plt.close(fig) + + if len(ap.omega_ranges) != 0: + # "k-space" plot ----------------------------- + domega = np.amin(np.diff(ap.allomegas)) + dsintheta = np.amin(abs(np.diff(sintheta))) + dk = dsintheta * wavenumbers[0] + + # target image grid + grid_y, grid_x = np.mgrid[ap.allomegas[0] : ap.allomegas[-1] : domega, np.amin(sintheta) * wavenumbers[-1] : np.amax(sintheta) * wavenumbers[-1] : dk] + imextent = (np.amin(sintheta) * wavenumbers[-1] / 1e6, np.amax(sintheta) * wavenumbers[-1] / 1e6, ap.allomegas[0] / eh, ap.allomegas[-1] / eh) + + # source coordinates for griddata + ktheta = sintheta[None, :] * wavenumbers[:, None] + omegapoints = np.broadcast_to(ap.allomegas[:, None], ktheta.shape) + points = np.stack( (ktheta.flatten(), omegapoints.flatten()), axis = -1) + + fig = plt.figure(figsize=(210/25.4, 297/25.4)) + vmax = np.amax(σ_ext_arr) + + ax = fig.add_subplot(311) + grid_z = griddata(points, σ_ext_arr.flatten(), (grid_x, grid_y), method = ipm) + ax.imshow(grid_z, extent = imextent, origin = 'lower', vmin = 0, vmax = vmax, aspect = 'auto', interpolation='none') + ax.set_xlabel('$k_\\theta / \\mathrm{\\mu m^{-1}}$') + ax.set_ylabel('$\\hbar\\omega / \\mathrm{eV}$') + ax.set_title('$\\sigma_\\mathrm{ext}$') + + ax = fig.add_subplot(312) + grid_z = griddata(points, σ_scat_arr.flatten(), (grid_x, grid_y), method = ipm) + ax.imshow(grid_z, extent = imextent, origin = 'lower', vmin = 0, vmax = vmax, aspect = 'auto', interpolation='none') + ax.set_xlabel('$k_\\theta / \\mathrm{\\mu m^{-1}}$') + ax.set_ylabel('$\\hbar\\omega / \\mathrm{eV}$') + ax.set_title('$\\sigma_\\mathrm{scat}$') + + ax = fig.add_subplot(313) + grid_z = griddata(points, σ_abs_arr.flatten(), (grid_x, grid_y), method = ipm) + im = ax.imshow(grid_z, extent = imextent, origin = 'lower', vmin = 0, vmax = vmax, aspect = 'auto', interpolation='none') + ax.set_xlabel('$k_\\theta / \\mathrm{\\mu m^{-1}}$') + ax.set_ylabel('$\\hbar\\omega / \\mathrm{eV}$') + ax.set_title('$\\sigma_\\mathrm{abs}$') + + fig.subplots_adjust(right=0.8) + fig.colorbar(im, cax = fig.add_axes([0.85, 0.15, 0.05, 0.7])) + + pdf.savefig(fig) + plt.close(fig) + + for omega in ap.omega_singles: + i = np.searchsorted(ap.allomegas, omega) + + fig = plt.figure() + fig.suptitle("%g eV" % (omega / eh)) + ax = fig.add_subplot(111) + sintheta = np.sin(a.theta) + ax.plot(sintheta, σ_ext_arr[i]*1e12,label='$\sigma_\mathrm{ext}$') + ax.plot(sintheta, σ_scat_arr[i]*1e12, label='$\sigma_\mathrm{scat}$') + ax.plot(sintheta, σ_abs_arr[i]*1e12, label='$\sigma_\mathrm{abs}$') + ax.legend() + ax.set_xlabel('$\sin\\theta$') + ax.set_ylabel('$\sigma/\mathrm{\mu m^2}$') + + pdf.savefig(fig) + plt.close(fig) exit(0) diff --git a/qpms/argproc.py b/qpms/argproc.py index 1502cda..6942dbd 100644 --- a/qpms/argproc.py +++ b/qpms/argproc.py @@ -93,6 +93,7 @@ class ArgParser: 'single_frequency_eV': lambda ap: ap.add_argument("-f", "--eV", type=float, required=True, help='radiation angular frequency in eV'), 'multiple_frequency_eV_optional': lambda ap: ap.add_argument("-f", "--eV", type=float, nargs='*', help='radiation angular frequency in eV (additional)'), 'seq_frequency_eV': lambda ap: ap.add_argument("-F", "--eV-seq", type=float, nargs=3, required=True, help='uniform radiation angular frequency sequence in eV', metavar=('FIRST', 'INCREMENT', 'LAST')), + 'real_frequencies_eV_ng': lambda ap: ap.add_argument("-f", "--eV", type=float_range, nargs='+', required=True, help='Angular frequency (or angular frequency range) in eV'), 'single_material': lambda ap: ap.add_argument("-m", "--material", help='particle material (Au, Ag, ... for Lorentz-Drude or number for constant refractive index)', default='Au', required=True), 'single_radius': lambda ap: ap.add_argument("-r", "--radius", type=float, required=True, help='particle radius (sphere or cylinder)'), 'single_height': lambda ap: ap.add_argument("-H", "--height", type=float, help='cylindrical particle height; if not provided, particle is assumed to be spherical'), @@ -114,6 +115,7 @@ class ArgParser: 'single_lMax': ("Single particle lMax definition", (), ('single_lMax',), ()), 'single_omega': ("Single angular frequency", (), ('single_frequency_eV',), ('_eval_single_omega',)), 'omega_seq': ("Equidistant real frequency range with possibility of adding individual frequencies", (), ('seq_frequency_eV', 'multiple_frequency_eV_optional',), ('_eval_omega_seq',)), + 'omega_seq_real_ng': ("Equidistant real frequency ranges or individual frequencies (new syntax)", (), ('real_frequencies_eV_ng',), ('_eval_omega_seq_real_ng',)), 'lattice2d': ("Specification of a generic 2d lattice (spanned by the x,y axes)", (), ('lattice2d_basis',), ('_eval_lattice2d',)), 'rectlattice2d': ("Specification of a rectangular 2d lattice; conflicts with lattice2d", (), ('rectlattice2d_periods',), ('_eval_rectlattice2d',)), 'rectlattice2d_finite': ("Specification of a rectangular 2d lattice; conflicts with lattice2d", ('rectlattice2d',), ('rectlattice2d_counts',), ()), @@ -199,6 +201,26 @@ class ArgParser: self.omegas.sort() self.omegas *= eV/hbar + def _eval_omega_seq_real_ng(self): # feature: omega_seq_real_ng + import numpy as np + from .constants import eV, hbar + eh = eV / hbar + self.omegas = [omega_eV * eh for omega_eV in self.args.eV] + self.omega_max = max(om if isinstance(om, float) else max(om) for om in self.omegas) + self.omega_min = min(om if isinstance(om, float) else min(om) for om in self.omegas) + self.omega_singles = [om for om in self.omegas if isinstance(om, float)] + self.omega_ranges = [om for om in self.omegas if not isinstance(om, float)] + self.omega_descr = ("%geV" % (self.omega_max / eh)) if (self.omega_max == self.omega_min) else ( + "%g–%geV" % (self.omega_min / eh, self.omega_max / eh)) + self.allomegas = [] + for om in self.omegas: + if isinstance(om, float): + self.allomegas.append(om) + else: + self.allomegas.extend(om) + self.allomegas = np.unique(self.allomegas) + + def _eval_lattice2d(self): # feature: lattice2d l = len(self.args.basis_vectors) if l != 2: raise ValueError('Two basis vectors must be specified (have %d)' % l)