WIP argproc multiparticle; update misc scripts
Former-commit-id: 50986e9e57873f8aca1ec7a8aafaa659434a839f
This commit is contained in:
Marek Nečada
parent
c5c148ca40
commit
3da6e5cfb2
|
|
@ -1,10 +1,11 @@
|
||||||
#!/usr/bin/env python3
|
#!/usr/bin/env python3
|
||||||
|
|
||||||
import math
|
import math
|
||||||
|
pi = math.pi
|
||||||
from qpms.argproc import ArgParser
|
from qpms.argproc import ArgParser
|
||||||
|
|
||||||
|
|
||||||
ap = ArgParser(['rectlattice2d_finite', 'single_particle', 'single_lMax', 'single_omega', 'planewave'])
|
ap = ArgParser(['rectlattice2d_finite', '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("-o", "--output", type=str, required=False, help='output path (if not provided, will be generated automatically)')
|
||||||
ap.add_argument("-O", "--plot-out", type=str, required=False, help="path to plot output (optional)")
|
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("-P", "--plot", action='store_true', help="if -p not given, plot to a default path")
|
||||||
|
|
@ -16,13 +17,8 @@ a=ap.parse_args()
|
||||||
import logging
|
import logging
|
||||||
logging.basicConfig(format='%(asctime)s %(message)s', level=logging.INFO)
|
logging.basicConfig(format='%(asctime)s %(message)s', level=logging.INFO)
|
||||||
|
|
||||||
Nx, Ny = a.size
|
|
||||||
px, py = a.period
|
|
||||||
|
|
||||||
|
|
||||||
import numpy as np
|
import numpy as np
|
||||||
import qpms
|
import qpms
|
||||||
import math
|
|
||||||
from qpms.qpms_p import cart2sph, sph2cart, sph_loccart2cart, sph_loccart_basis
|
from qpms.qpms_p import cart2sph, sph2cart, sph_loccart2cart, sph_loccart_basis
|
||||||
from qpms.cybspec import BaseSpec
|
from qpms.cybspec import BaseSpec
|
||||||
from qpms.cytmatrices import CTMatrix, TMatrixGenerator
|
from qpms.cytmatrices import CTMatrix, TMatrixGenerator
|
||||||
|
|
@ -32,15 +28,18 @@ from qpms.cycommon import DebugFlags, dbgmsg_enable
|
||||||
from qpms import FinitePointGroup, ScatteringSystem, BesselType, eV, hbar
|
from qpms import FinitePointGroup, ScatteringSystem, BesselType, eV, hbar
|
||||||
from qpms.symmetries import point_group_info
|
from qpms.symmetries import point_group_info
|
||||||
eh = eV/hbar
|
eh = eV/hbar
|
||||||
pi = math.pi
|
|
||||||
|
dbgmsg_enable(DebugFlags.INTEGRATION)
|
||||||
|
|
||||||
|
Nx, Ny = a.size
|
||||||
|
px, py = a.period
|
||||||
|
|
||||||
particlestr = ("sph" if a.height is None else "cyl") + ("_r%gnm" % (a.radius*1e9))
|
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)
|
if a.height is not None: particlestr += "_h%gnm" % (a.height * 1e9)
|
||||||
defaultprefix = "%s_p%gnmx%gnm_%dx%d_m%s_n%g_φ%gπ_θ(%g_%g)π_ψ%gπ_χ%gπ_f%geV_L%d" % (
|
defaultprefix = "%s_p%gnmx%gnm_%dx%d_m%s_bg%s_φ%gπ_θ(%g_%g)π_ψ%gπ_χ%gπ_f%s_L%d" % (
|
||||||
particlestr, px*1e9, py*1e9, Nx, Ny, 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, a.eV, a.lMax, )
|
particlestr, px*1e9, py*1e9, Nx, Ny, str(a.material), str(a.background), 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)
|
logging.info("Default file prefix: %s" % defaultprefix)
|
||||||
|
|
||||||
dbgmsg_enable(DebugFlags.INTEGRATION)
|
|
||||||
|
|
||||||
#Particle positions
|
#Particle positions
|
||||||
orig_x = (np.arange(Nx/2) + (0 if (Nx % 2) else .5)) * px
|
orig_x = (np.arange(Nx/2) + (0 if (Nx % 2) else .5)) * px
|
||||||
|
|
@ -48,91 +47,191 @@ orig_y = (np.arange(Ny/2) + (0 if (Ny % 2) else .5)) * py
|
||||||
|
|
||||||
orig_xy = np.stack(np.meshgrid(orig_x, orig_y), axis = -1)
|
orig_xy = np.stack(np.meshgrid(orig_x, orig_y), axis = -1)
|
||||||
|
|
||||||
|
|
||||||
omega = ap.omega
|
|
||||||
|
|
||||||
bspec = BaseSpec(lMax = a.lMax)
|
bspec = BaseSpec(lMax = a.lMax)
|
||||||
Tmatrix = ap.tmgen(bspec, ap.omega)
|
particles= [Particle(orig_xy[i], ap.tmgen, bspec=bspec) for i in np.ndindex(orig_xy.shape[:-1])]
|
||||||
particles= [Particle(orig_xy[i], Tmatrix) for i in np.ndindex(orig_xy.shape[:-1])]
|
|
||||||
|
|
||||||
sym = FinitePointGroup(point_group_info['D2h'])
|
sym = FinitePointGroup(point_group_info['D2h'])
|
||||||
ss, ssw = ScatteringSystem.create(particles, ap.background_emg, omega, sym=sym)
|
ss, ssw = ScatteringSystem.create(particles, ap.background_emg, ap.allomegas[0], sym=sym)
|
||||||
|
|
||||||
wavenumber = ap.background_epsmu.k(omega).real # Currently, ScatteringSystem does not "remember" frequency nor wavenumber
|
|
||||||
|
|
||||||
## Plane wave data
|
## Plane wave data
|
||||||
a.theta = np.array(a.theta)
|
a.theta = np.atleast_1d(np.array(a.theta))
|
||||||
k_sph_list = np.stack((np.broadcast_to(wavenumber, a.theta.shape), a.theta, np.broadcast_to(a.phi, a.theta.shape)), axis=-1)
|
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.psi), math.cos(a.psi)
|
||||||
sχ, cχ = math.sin(a.chi), math.cos(a.chi)
|
sχ, cχ = math.sin(a.chi), math.cos(a.chi)
|
||||||
E_sph = (0., cψ*cχ + 1j*sψ*sχ, sψ*cχ + 1j*cψ*sχ)
|
E_sph = (0., cψ*cχ + 1j*sψ*sχ, sψ*cχ + 1j*cψ*sχ)
|
||||||
|
|
||||||
k_cart_list = sph2cart(k_sph_list)
|
dir_cart_list = sph2cart(dir_sph_list)
|
||||||
E_cart_list = sph_loccart2cart(E_sph, k_sph_list)
|
E_cart_list = sph_loccart2cart(E_sph, dir_sph_list)
|
||||||
|
|
||||||
npoints = a.theta.shape[0]
|
nfreq = len(ap.allomegas)
|
||||||
|
ndir = a.theta.shape[0]
|
||||||
|
|
||||||
σ_ext_list_ir = np.empty((npoints, ss.nirreps), dtype=float)
|
k_cart_arr = np.empty((nfreq, ndir, 3), dtype=float)
|
||||||
σ_scat_list_ir = np.empty((npoints, ss.nirreps), dtype=float)
|
wavenumbers = np.empty((nfreq,), dtype=float)
|
||||||
|
|
||||||
|
σ_ext_arr_ir = np.empty((nfreq, ndir, ss.nirreps), dtype=float)
|
||||||
|
σ_scat_arr_ir = np.empty((nfreq, ndir, ss.nirreps), dtype=float)
|
||||||
|
|
||||||
outfile_tmp = defaultprefix + ".tmp" if a.output is None else a.output + ".tmp"
|
outfile_tmp = defaultprefix + ".tmp" if a.output is None else a.output + ".tmp"
|
||||||
|
|
||||||
for iri in range(ss.nirreps):
|
for i, omega in enumerate(ap.allomegas):
|
||||||
logging.info("processing irrep %d/%d" % (iri, ss.nirreps))
|
logging.info("Processing frequency %g eV" % (omega / eV,))
|
||||||
LU = None # to trigger garbage collection before the next call
|
if i != 0:
|
||||||
translation_matrix = None
|
ssw = ss(omega)
|
||||||
LU = ssw.scatter_solver(iri)
|
if ssw.wavenumber.imag != 0:
|
||||||
logging.info("LU solver created")
|
warnings.warn("The background medium wavenumber has non-zero imaginary part. Don't expect emaningful results for cross sections.")
|
||||||
translation_matrix = ss.translation_matrix_packed(wavenumber, iri, BesselType.REGULAR) + np.eye(ss.saecv_sizes[iri])
|
wavenumber = ssw.wavenumber.real
|
||||||
logging.info("auxillary translation matrix created")
|
wavenumbers[i] = wavenumber
|
||||||
|
|
||||||
for j in range(npoints):
|
k_sph_list = np.array(dir_sph_list, copy=True)
|
||||||
# the following two could be calculated only once, but probably not a big deal
|
k_sph_list[:,0] = wavenumber
|
||||||
ã = ss.planewave_full(k_cart=k_cart_list[j], E_cart=E_cart_list[j])
|
|
||||||
Tã = ssw.apply_Tmatrices_full(ã)
|
|
||||||
|
|
||||||
Tãi = ss.pack_vector(Tã, iri)
|
k_cart_arr[i] = sph2cart(k_sph_list)
|
||||||
ãi = ss.pack_vector(ã, iri)
|
|
||||||
fi = LU(Tãi)
|
|
||||||
σ_ext_list_ir[j, iri] = -np.vdot(ãi, fi).real/wavenumber**2
|
|
||||||
σ_scat_list_ir[j, iri] = np.vdot(fi,np.dot(translation_matrix, fi)).real/wavenumber**2
|
|
||||||
if a.save_gradually:
|
|
||||||
iriout = outfile_tmp + ".%d" % iri
|
|
||||||
np.savez(iriout, iri=iri, meta=vars(a), k_sph=k_sph_list, k_cart = k_cart_list, E_cart=E_cart_list, E_sph=np.array(E_sph),
|
|
||||||
omega=omega, wavenumber=wavenumber, σ_ext_list_ir=σ_ext_list_ir[:,iri], σ_scat_list_ir=σ_scat_list_ir[:,iri])
|
|
||||||
logging.info("partial results saved to %s"%iriout)
|
|
||||||
|
|
||||||
σ_abs_list_ir = σ_ext_list_ir - σ_scat_list_ir
|
for iri in range(ss.nirreps):
|
||||||
σ_abs= np.sum(σ_abs_list_ir, axis=-1)
|
logging.info("processing irrep %d/%d" % (iri, ss.nirreps))
|
||||||
σ_scat= np.sum(σ_scat_list_ir, axis=-1)
|
LU = None # to trigger garbage collection before the next call
|
||||||
σ_ext= np.sum(σ_ext_list_ir, axis=-1)
|
translation_matrix = None
|
||||||
|
LU = ssw.scatter_solver(iri)
|
||||||
|
logging.info("LU solver created")
|
||||||
|
translation_matrix = ssw.translation_matrix_packed(iri, BesselType.REGULAR) + np.eye(ss.saecv_sizes[iri])
|
||||||
|
logging.info("auxillary translation matrix created")
|
||||||
|
|
||||||
|
for j in range(ndir):
|
||||||
|
k_cart = k_cart_arr[i,j]
|
||||||
|
# the following two could be calculated only once, but probably not a big deal
|
||||||
|
ã = ss.planewave_full(k_cart=k_cart_arr[i,j], E_cart=E_cart_list[j])
|
||||||
|
Tã = ssw.apply_Tmatrices_full(ã)
|
||||||
|
|
||||||
|
Tãi = ss.pack_vector(Tã, iri)
|
||||||
|
ãi = ss.pack_vector(ã, iri)
|
||||||
|
fi = LU(Tãi)
|
||||||
|
σ_ext_arr_ir[i, j, iri] = -np.vdot(ãi, fi).real/wavenumber**2
|
||||||
|
σ_scat_arr_ir[i, j, iri] = np.vdot(fi,np.dot(translation_matrix, fi)).real/wavenumber**2
|
||||||
|
if a.save_gradually:
|
||||||
|
iriout = outfile_tmp + ".%d.%d" % (i, iri)
|
||||||
|
np.savez(iriout, omegai=i, iri=iri, meta=vars(a), omega=omega, k_sph=k_sph_list, k_cart = k_cart_arr, E_cart=E_cart_list, E_sph=np.array(E_sph),
|
||||||
|
wavenumber=wavenumber, σ_ext_list_ir=σ_ext_arr_ir[i,:,iri], σ_scat_list_ir=σ_scat_list_ir[i,:,iri])
|
||||||
|
logging.info("partial results saved to %s"%iriout)
|
||||||
|
|
||||||
|
σ_abs_arr_ir = σ_ext_arr_ir - σ_scat_arr_ir
|
||||||
|
σ_abs_arr = np.sum(σ_abs_arr_ir, axis=-1)
|
||||||
|
σ_scat_arr = np.sum(σ_scat_arr_ir, axis=-1)
|
||||||
|
σ_ext_arr = np.sum(σ_ext_arr_ir, axis=-1)
|
||||||
|
|
||||||
|
|
||||||
outfile = defaultprefix + ".npz" if a.output is None else a.output
|
outfile = defaultprefix + ".npz" if a.output is None else a.output
|
||||||
np.savez(outfile, meta=vars(a), k_sph=k_sph_list, k_cart = k_cart_list, E_cart=E_cart_list, E_sph=np.array(E_sph),
|
np.savez(outfile, meta=vars(a), k_sph=k_sph_list, k_cart = k_cart_arr, E_cart=E_cart_list, E_sph=np.array(E_sph),
|
||||||
σ_ext=σ_ext,σ_abs=σ_abs,σ_scat=σ_scat,
|
σ_ext=σ_ext_arr,σ_abs=σ_abs_arr,σ_scat=σ_scat_arr,
|
||||||
σ_ext_ir=σ_ext_list_ir,σ_abs_ir=σ_abs_list_ir,σ_scat_ir=σ_scat_list_ir, omega=omega, wavenumber=wavenumber
|
σ_ext_ir=σ_ext_arr_ir,σ_abs_ir=σ_abs_arr_ir,σ_scat_ir=σ_scat_arr_ir, omega=ap.allomegas, wavenumbers=wavenumbers
|
||||||
)
|
)
|
||||||
logging.info("Saved to %s" % outfile)
|
logging.info("Saved to %s" % outfile)
|
||||||
|
|
||||||
|
|
||||||
if a.plot or (a.plot_out is not None):
|
if a.plot or (a.plot_out is not None):
|
||||||
import matplotlib
|
import matplotlib
|
||||||
|
from matplotlib.backends.backend_pdf import PdfPages
|
||||||
matplotlib.use('pdf')
|
matplotlib.use('pdf')
|
||||||
from matplotlib import pyplot as plt
|
from matplotlib import pyplot as plt
|
||||||
|
from scipy.interpolate import griddata
|
||||||
|
|
||||||
fig = plt.figure()
|
|
||||||
ax = fig.add_subplot(111)
|
|
||||||
sintheta = np.sin(a.theta)
|
|
||||||
ax.plot(sintheta, σ_ext*1e12,label='$\sigma_\mathrm{ext}$')
|
|
||||||
ax.plot(sintheta, σ_scat*1e12, label='$\sigma_\mathrm{scat}$')
|
|
||||||
ax.plot(sintheta, σ_abs*1e12, label='$\sigma_\mathrm{abs}$')
|
|
||||||
ax.legend()
|
|
||||||
ax.set_xlabel('$\sin\\theta$')
|
|
||||||
ax.set_ylabel('$\sigma/\mathrm{\mu m^2}$')
|
|
||||||
|
|
||||||
plotfile = defaultprefix + ".pdf" if a.plot_out is None else a.plot_out
|
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)
|
exit(0)
|
||||||
|
|
||||||
|
|
|
||||||
|
|
@ -34,8 +34,8 @@ px, py = a.period
|
||||||
|
|
||||||
particlestr = ("sph" if a.height is None else "cyl") + ("_r%gnm" % (a.radius*1e9))
|
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)
|
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" % (
|
defaultprefix = "%s_p%gnmx%gnm_m%s_bg%s_φ%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)
|
particlestr, px*1e9, py*1e9, str(a.material), str(a.background), 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)
|
logging.info("Default file prefix: %s" % defaultprefix)
|
||||||
|
|
||||||
|
|
||||||
|
|
|
||||||
|
|
@ -3,7 +3,7 @@
|
||||||
import math
|
import math
|
||||||
from qpms.argproc import ArgParser
|
from qpms.argproc import ArgParser
|
||||||
|
|
||||||
ap = ArgParser(['rectlattice2d', 'single_particle', 'single_lMax'])
|
ap = ArgParser(['rectlattice2d', 'const_real_background', 'single_particle', 'single_lMax']) # const_real_background needed for calculation of the diffracted orders
|
||||||
ap.add_argument("-k", nargs=2, type=float, required=True, help='k vector', metavar=('K_X', 'K_Y'))
|
ap.add_argument("-k", nargs=2, type=float, required=True, help='k vector', metavar=('K_X', 'K_Y'))
|
||||||
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.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.add_argument("--rank-tol", type=float, required=False)
|
ap.add_argument("--rank-tol", type=float, required=False)
|
||||||
|
|
|
||||||
291
qpms/argproc.py
291
qpms/argproc.py
|
|
@ -4,6 +4,14 @@ Common snippets for argument processing in command line scripts; legacy scripts
|
||||||
|
|
||||||
import argparse
|
import argparse
|
||||||
import sys
|
import sys
|
||||||
|
import warnings
|
||||||
|
|
||||||
|
def flatten(S):
|
||||||
|
if S == []:
|
||||||
|
return S
|
||||||
|
if isinstance(S[0], list):
|
||||||
|
return flatten(S[0]) + flatten(S[1:])
|
||||||
|
return S[:1] + flatten(S[1:])
|
||||||
|
|
||||||
def make_action_sharedlist(opname, listname):
|
def make_action_sharedlist(opname, listname):
|
||||||
class opAction(argparse.Action):
|
class opAction(argparse.Action):
|
||||||
|
|
@ -13,6 +21,40 @@ def make_action_sharedlist(opname, listname):
|
||||||
getattr(args, listname).append((opname, values))
|
getattr(args, listname).append((opname, values))
|
||||||
return opAction
|
return opAction
|
||||||
|
|
||||||
|
def make_dict_action(argtype=None, postaction='store', first_is_key=True):
|
||||||
|
class DictAction(argparse.Action):
|
||||||
|
#def __init__(self, option_strings, dest, nargs=None, **kwargs):
|
||||||
|
# if nargs is not None:
|
||||||
|
# raise ValueError("nargs not allowed")
|
||||||
|
# super(DictAction, self).__init__(option_strings, dest, **kwargs)
|
||||||
|
def __call__(self, parser, namespace, values, option_string=None):
|
||||||
|
if first_is_key: # For the labeled versions
|
||||||
|
key = values[0]
|
||||||
|
vals = values[1:]
|
||||||
|
else: # For the default values
|
||||||
|
key = None
|
||||||
|
vals = values
|
||||||
|
if argtype is not None:
|
||||||
|
if (first_is_key and self.nargs == 2) or (not first_is_key and self.nargs == 1):
|
||||||
|
vals = argtype(vals[0]) # avoid having lists in this case
|
||||||
|
else:
|
||||||
|
vals = [argtype(val) for val in vals]
|
||||||
|
ledict = getattr(namespace, self.dest, {})
|
||||||
|
if ledict is None:
|
||||||
|
ledict = {}
|
||||||
|
if postaction=='store':
|
||||||
|
ledict[key] = vals
|
||||||
|
elif postaction=='append':
|
||||||
|
lelist = ledict.get(key, [])
|
||||||
|
lelist.append(vals)
|
||||||
|
ledict[key] = lelist
|
||||||
|
setattr(namespace, self.dest, ledict)
|
||||||
|
return DictAction
|
||||||
|
|
||||||
|
|
||||||
|
class ArgumentProcessingError(Exception):
|
||||||
|
pass
|
||||||
|
|
||||||
class AppendTupleAction(argparse.Action):
|
class AppendTupleAction(argparse.Action):
|
||||||
''' A variation on the 'append' builtin action from argparse, but uses tuples for the internal groupings instead of lists '''
|
''' A variation on the 'append' builtin action from argparse, but uses tuples for the internal groupings instead of lists '''
|
||||||
def __call__(self, parser, args, values, option_string=None):
|
def __call__(self, parser, args, values, option_string=None):
|
||||||
|
|
@ -39,17 +81,14 @@ def float_range(string):
|
||||||
steps = None
|
steps = None
|
||||||
match = re.match(r's?([^:]+):([^|]+)\|(.+)', string)
|
match = re.match(r's?([^:]+):([^|]+)\|(.+)', string)
|
||||||
if match:
|
if match:
|
||||||
#print("first:last|steps", match.group(1,2,3))
|
|
||||||
steps = int(match.group(3))
|
steps = int(match.group(3))
|
||||||
else:
|
else:
|
||||||
match = re.match(r's?([^:]+):([^:]+):(.+)', string)
|
match = re.match(r's?([^:]+):([^:]+):(.+)', string)
|
||||||
if match:
|
if match:
|
||||||
#print("first:last:increment", match.group(1,2,3))
|
|
||||||
increment = float(match.group(3))
|
increment = float(match.group(3))
|
||||||
else:
|
else:
|
||||||
match = re.match(r's?([^:]+):(.+)', string)
|
match = re.match(r's?([^:]+):(.+)', string)
|
||||||
if match:
|
if match:
|
||||||
#print("first:last", match.group(1,2))
|
|
||||||
steps = 50
|
steps = 50
|
||||||
else:
|
else:
|
||||||
argparse.ArgumentTypeError('Invalid float/sequence format: "%s"' % string)
|
argparse.ArgumentTypeError('Invalid float/sequence format: "%s"' % string)
|
||||||
|
|
@ -61,6 +100,24 @@ def float_range(string):
|
||||||
else:
|
else:
|
||||||
return np.arange(first, last, increment)
|
return np.arange(first, last, increment)
|
||||||
|
|
||||||
|
def material_spec(string):
|
||||||
|
"""Tries to parse a string as a material specification, i.e. a
|
||||||
|
real or complex number or one of the string in built-in Lorentz-Drude models.
|
||||||
|
|
||||||
|
Tries to interpret the string as 1) float, 2) complex, 3) Lorentz-Drude key.
|
||||||
|
Raises argparse.ArgumentTypeError on failure.
|
||||||
|
"""
|
||||||
|
from .cymaterials import lorentz_drude
|
||||||
|
if string in lorentz_drude.keys():
|
||||||
|
return string
|
||||||
|
else:
|
||||||
|
try: lemat = float(string)
|
||||||
|
except ValueError:
|
||||||
|
try: lemat = complex(string)
|
||||||
|
except ValueError as ve:
|
||||||
|
raise argpares.ArgumentTypeError("Material specification must be a supported material name %s, or a number" % (str(lorentz_drude.keys()),)) from ve
|
||||||
|
return lemat
|
||||||
|
|
||||||
class ArgParser:
|
class ArgParser:
|
||||||
''' Common argument parsing engine for QPMS python CLI scripts. '''
|
''' Common argument parsing engine for QPMS python CLI scripts. '''
|
||||||
|
|
||||||
|
|
@ -87,56 +144,143 @@ class ArgParser:
|
||||||
pwgrp.add_argument("-χ", "--chi", type=float, default=0,
|
pwgrp.add_argument("-χ", "--chi", type=float, default=0,
|
||||||
help='Polarisation parameter χ in multiples of π/2. 0 for linear, 0.5 for circular pol.')
|
help='Polarisation parameter χ in multiples of π/2. 0 for linear, 0.5 for circular pol.')
|
||||||
|
|
||||||
|
def __add_manyparticle_argparse_group(ap):
|
||||||
|
mpgrp = ap.add_argument_group('Many particle specification', "TODO DOC")
|
||||||
|
mpgrp.add_argument("-p", "--position", nargs='+', action=make_dict_action(argtype=float, postaction='append',
|
||||||
|
first_is_key=False), help="Particle positions, cartesion coordinates (default particle properties)")
|
||||||
|
mpgrp.add_argument("+p", "++position", nargs='+', action=make_dict_action(argtype=float, postaction='append',
|
||||||
|
first_is_key=True), help="Particle positions, cartesian coordinates (labeled)")
|
||||||
|
mpgrp.add_argument("-L", "--lMax", nargs=1, default={},
|
||||||
|
action=make_dict_action(argtype=int, postaction='store', first_is_key=False,),
|
||||||
|
help="Cutoff multipole degree (default)")
|
||||||
|
mpgrp.add_argument("+L", "++lMax", nargs=2,
|
||||||
|
action=make_dict_action(argtype=int, postaction='store', first_is_key=True,),
|
||||||
|
help="Cutoff multipole degree (labeled)")
|
||||||
|
mpgrp.add_argument("-m", "--material", nargs=1, default={},
|
||||||
|
action=make_dict_action(argtype=material_spec, postaction='store', first_is_key=False,),
|
||||||
|
help='particle material (Au, Ag, ... for Lorentz-Drude or number for constant refractive index) (default)')
|
||||||
|
mpgrp.add_argument("+m", "++material", nargs=2,
|
||||||
|
action=make_dict_action(argtype=material_spec, postaction='store', first_is_key=True,),
|
||||||
|
help='particle material (Au, Ag, ... for Lorentz-Drude or number for constant refractive index) (labeled)')
|
||||||
|
mpgrp.add_argument("-r", "--radius", nargs=1, default={},
|
||||||
|
action=make_dict_action(argtype=float, postaction='store', first_is_key=False,),
|
||||||
|
help='particle radius (sphere or cylinder; default)')
|
||||||
|
mpgrp.add_argument("+r", "++radius", nargs=2,
|
||||||
|
action=make_dict_action(argtype=float, postaction='store', first_is_key=True,),
|
||||||
|
help='particle radius (sphere or cylinder; labeled)')
|
||||||
|
mpgrp.add_argument("-H", "--height", nargs=1, default={},
|
||||||
|
action=make_dict_action(argtype=float, postaction='store', first_is_key=False,),
|
||||||
|
help='particle radius (cylinder; default)')
|
||||||
|
mpgrp.add_argument("+H", "++height", nargs=2,
|
||||||
|
action=make_dict_action(argtype=float, postaction='store', first_is_key=True,),
|
||||||
|
help='particle radius (cylinder; labeled)')
|
||||||
|
|
||||||
atomic_arguments = {
|
atomic_arguments = {
|
||||||
'rectlattice2d_periods': lambda ap: ap.add_argument("-p", "--period", type=float, nargs='+', required=True, help='square/rectangular lattice periods', metavar=('px','[py]')),
|
'rectlattice2d_periods': lambda ap: ap.add_argument("-p", "--period", type=float, nargs='+', required=True, help='square/rectangular lattice periods', metavar=('px','[py]')),
|
||||||
'rectlattice2d_counts': lambda ap: ap.add_argument("--size", type=int, nargs=2, required=True, help='rectangular array size (particle column, row count)', metavar=('NCOLS', 'NROWS')),
|
'rectlattice2d_counts': lambda ap: ap.add_argument("--size", type=int, nargs=2, required=True, help='rectangular array size (particle column, row count)', metavar=('NCOLS', 'NROWS')),
|
||||||
'single_frequency_eV': lambda ap: ap.add_argument("-f", "--eV", type=float, required=True, help='radiation angular frequency in eV'),
|
'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)'),
|
'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')),
|
'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'),
|
'real_frequencies_eV_ng': lambda ap: ap.add_argument("-f", "--eV", type=float_range, nargs=1, action='append', required=True, help='Angular frequency (or angular frequency range) in eV'), # nargs='+', action='extend' would be better, but action='extend' requires python>=3.8
|
||||||
'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_material': lambda ap: ap.add_argument("-m", "--material", help='particle material (Au, Ag, ... for Lorentz-Drude or number for constant refractive index)', type=material_spec, required=True),
|
||||||
'single_radius': lambda ap: ap.add_argument("-r", "--radius", type=float, required=True, help='particle radius (sphere or cylinder)'),
|
'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'),
|
'single_height': lambda ap: ap.add_argument("-H", "--height", type=float, help='cylindrical particle height; if not provided, particle is assumed to be spherical'),
|
||||||
'single_kvec2': lambda ap: ap.add_argument("-k", '--kx-lim', nargs=2, type=float, required=True, help='k vector', metavar=('KX_MIN', 'KX_MAX')),
|
'single_kvec2': lambda ap: ap.add_argument("-k", '--kx-lim', nargs=2, type=float, required=True, help='k vector', metavar=('KX_MIN', 'KX_MAX')),
|
||||||
'kpi': lambda ap: 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)"),
|
'kpi': lambda ap: 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)"),
|
||||||
'bg_refractive_index': lambda ap: ap.add_argument("-n", "--refractive-index", type=float, default=1.52, help='background medium refractive index'),
|
'bg_real_refractive_index': lambda ap: ap.add_argument("-n", "--refractive-index", type=float, default=1., help='background medium strictly real refractive index'),
|
||||||
|
'bg_analytical': lambda ap: ap.add_argument("-B", "--background", type=material_spec, default=1., help="Background medium specification (constant real or complex refractive index, or supported material label)"),
|
||||||
'single_lMax': lambda ap: ap.add_argument("-L", "--lMax", type=int, required=True, default=3, help='multipole degree cutoff'),
|
'single_lMax': lambda ap: ap.add_argument("-L", "--lMax", type=int, required=True, default=3, help='multipole degree cutoff'),
|
||||||
'single_lMax_extend': lambda ap: ap.add_argument("--lMax-extend", type=int, required=False, default=6, help='multipole degree cutoff for T-matrix calculation (cylindrical particles only'),
|
'single_lMax_extend': lambda ap: ap.add_argument("--lMax-extend", type=int, required=False, default=6, help='multipole degree cutoff for T-matrix calculation (cylindrical particles only'),
|
||||||
'outfile': lambda ap: ap.add_argument("-o", "--output", type=str, required=False, help='output path (if not provided, will be generated automatically)'), # TODO consider type=argparse.FileType('w')
|
'outfile': lambda ap: ap.add_argument("-o", "--output", type=str, required=False, help='output path (if not provided, will be generated automatically)'), # TODO consider type=argparse.FileType('w')
|
||||||
'plot_out': lambda ap: ap.add_argument("-O", "--plot-out", type=str, required=False, help="path to plot output (optional)"),
|
'plot_out': lambda ap: ap.add_argument("-O", "--plot-out", type=str, required=False, help="path to plot output (optional)"),
|
||||||
'plot_do': lambda ap: ap.add_argument("-P", "--plot", action='store_true', help="if -p not given, plot to a default path"),
|
'plot_do': lambda ap: ap.add_argument("-P", "--plot", action='store_true', help="if -p not given, plot to a default path"),
|
||||||
'lattice2d_basis': lambda ap: ap.add_argument("-b", "--basis-vector", action=AppendTupleAction, help="basis vector in xy-cartesian coordinates (two required)", dest='basis_vectors', metavar=('X', 'Y')),
|
'lattice2d_basis': lambda ap: ap.add_argument("-b", "--basis-vector", nargs='+', action=AppendTupleAction, help="basis vector in xy-cartesian coordinates (two required)", required=True, dest='basis_vectors', metavar=('X', 'Y')),
|
||||||
'planewave_pol_angles': __add_planewave_argparse_group,
|
'planewave_pol_angles': __add_planewave_argparse_group,
|
||||||
|
'multi_particle': __add_manyparticle_argparse_group,
|
||||||
}
|
}
|
||||||
|
|
||||||
feature_sets_available = { # name : (description, dependencies, atoms not in other dependencies, methods called after parsing)
|
feature_sets_available = { # name : (description, dependencies, atoms not in other dependencies, methods called after parsing, "virtual" features provided)
|
||||||
'background': ("Background medium definition (currently only constant epsilon supported)", (), ('bg_refractive_index',), ('_eval_background_epsmu',)),
|
'const_real_background': ("Background medium with constant real refractive index", (), ('bg_real_refractive_index',), ('_eval_const_background_epsmu',), ('background', 'background_analytical')),
|
||||||
'single_particle': ("Single particle definition (shape [currently spherical or cylindrical]) and materials, incl. background)", ('background',), ('single_material', 'single_radius', 'single_height', 'single_lMax_extend'), ('_eval_single_tmgen',)),
|
'background' : ("Most general background medium specification currently supported", ('background_analytical',), (), (), ()),
|
||||||
'single_lMax': ("Single particle lMax definition", (), ('single_lMax',), ()),
|
'background_analytical' : ("Background medium model holomorphic for 'reasonably large' complex frequency areas", (), ('bg_analytical',), ('_eval_analytical_background_epsmugen',), ('background',)),
|
||||||
'single_omega': ("Single angular frequency", (), ('single_frequency_eV',), ('_eval_single_omega',)),
|
'single_particle': ("Single particle definition (shape [currently spherical or cylindrical]) and materials, incl. background)", ('background',), ('single_material', 'single_radius', 'single_height', 'single_lMax_extend'), ('_eval_single_tmgen',), ()),
|
||||||
'omega_seq': ("Equidistant real frequency range with possibility of adding individual frequencies", (), ('seq_frequency_eV', 'multiple_frequency_eV_optional',), ('_eval_omega_seq',)),
|
'multi_particle': ("One or more particle definition (shape [curently spherical or cylindrical]), materials, and positions)", ('background',), ('multi_particle',), ('_process_multi_particle',), ()),
|
||||||
'omega_seq_real_ng': ("Equidistant real frequency ranges or individual frequencies (new syntax)", (), ('real_frequencies_eV_ng',), ('_eval_omega_seq_real_ng',)),
|
'single_lMax': ("Single particle lMax definition", (), ('single_lMax',), (), ()),
|
||||||
'lattice2d': ("Specification of a generic 2d lattice (spanned by the x,y axes)", (), ('lattice2d_basis',), ('_eval_lattice2d',)),
|
'single_omega': ("Single angular frequency", (), ('single_frequency_eV',), ('_eval_single_omega',), ()),
|
||||||
'rectlattice2d': ("Specification of a rectangular 2d lattice; conflicts with lattice2d", (), ('rectlattice2d_periods',), ('_eval_rectlattice2d',)),
|
'omega_seq': ("Equidistant real frequency range with possibility of adding individual frequencies", (), ('seq_frequency_eV', 'multiple_frequency_eV_optional',), ('_eval_omega_seq',), ()),
|
||||||
'rectlattice2d_finite': ("Specification of a rectangular 2d lattice; conflicts with lattice2d", ('rectlattice2d',), ('rectlattice2d_counts',), ()),
|
'omega_seq_real_ng': ("Equidistant real frequency ranges or individual frequencies (new syntax)", (), ('real_frequencies_eV_ng',), ('_eval_omega_seq_real_ng',), ()),
|
||||||
'planewave': ("Specification of a normalised plane wave (typically used for scattering) with a full polarisation state", (), ('planewave_pol_angles',), ("_process_planewave_angles",)),
|
'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',), (), ()),
|
||||||
|
'planewave': ("Specification of a normalised plane wave (typically used for scattering) with a full polarisation state", (), ('planewave_pol_angles',), ("_process_planewave_angles",), ()),
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
def __init__(self, features=[]):
|
def __init__(self, features=[]):
|
||||||
self.ap = argparse.ArgumentParser()
|
prefix_chars = '+-' if 'multi_particle' in features else '-'
|
||||||
|
self.ap = argparse.ArgumentParser(prefix_chars=prefix_chars)
|
||||||
self.features_enabled = set()
|
self.features_enabled = set()
|
||||||
self.call_at_parse_list = []
|
self.call_at_parse_list = []
|
||||||
self.parsed = False
|
self.parsed = False
|
||||||
for feat in features:
|
for feat in features:
|
||||||
self.add_feature(feat)
|
self.add_feature(feat)
|
||||||
|
self._emg_register = {} # EpsMuGenerator dictionary to avoid recreating equivalent instances; filled by _add_emg()
|
||||||
|
self._tmg_register = {} # TMatrixGenerator dictionary to avoid recreating equivalent instances; filled by _add_tmg()
|
||||||
|
self._bspec_register = {} # Dictionary of used BaseSpecs to keep the equivalent instances unique; filled by _add_bspec()
|
||||||
|
|
||||||
|
def _add_emg(self, emgspec):
|
||||||
|
"""Looks up whether if an EpsMuGenerator from given material_spec has been already registered, and if not, creates a new one"""
|
||||||
|
from .cymaterials import EpsMu, EpsMuGenerator, lorentz_drude
|
||||||
|
if emgspec in self._emg_register.keys():
|
||||||
|
return self._emg_register[emgspec]
|
||||||
|
else:
|
||||||
|
if isinstance(emgspec, (float, complex)):
|
||||||
|
emg = EpsMuGenerator(EpsMu(emgspec**2))
|
||||||
|
else:
|
||||||
|
emg = EpsMuGenerator(lorentz_drude[emgspec])
|
||||||
|
self._emg_register[emgspec] = emg
|
||||||
|
return emg
|
||||||
|
|
||||||
|
def _add_tmg(self, tmgspec):
|
||||||
|
"""Looks up whether if a T-matrix from given T-matrix specification tuple has been already registered, and if not, creates a new one
|
||||||
|
|
||||||
|
T-matrix specification shall be of the form
|
||||||
|
(bg_material_spec, fg_material_spec, shape_spec) where shape_spec is
|
||||||
|
(radius, height, lMax_extend)
|
||||||
|
"""
|
||||||
|
if tmgspec in self._tmg_register.keys():
|
||||||
|
return self._tmg_register[tmgspec]
|
||||||
|
else:
|
||||||
|
from .cytmatrices import TMatrixGenerator
|
||||||
|
bgspec, fgspec, (radius, height, lMax_extend) = tmgspec
|
||||||
|
bg = self._add_emg(bgspec)
|
||||||
|
fg = self._add_emg(fgspec)
|
||||||
|
if height is None:
|
||||||
|
tmgen = TMatrixGenerator.sphere(bg, fg, radius)
|
||||||
|
else:
|
||||||
|
tmgen = TMatrixGenerator.cylinder(bg, fg, radius, height, lMax_extend=lMax_extend)
|
||||||
|
self._tmg_register[tmgspec] = tmgen
|
||||||
|
return tmgen
|
||||||
|
|
||||||
|
def _add_bspec(self, key):
|
||||||
|
if key in self._bspec_register.keys():
|
||||||
|
return self._bspec_register[key]
|
||||||
|
else:
|
||||||
|
from .cybspec import BaseSpec
|
||||||
|
if isinstance(key, BaseSpec):
|
||||||
|
bspec = key
|
||||||
|
elif isinstance(key, int):
|
||||||
|
bspec = self._add_bspec(BaseSpec(lMax=key))
|
||||||
|
else: raise TypeError("Can't register this as a BaseSpec")
|
||||||
|
self._bspec_register[key] = bspec
|
||||||
|
return bspec
|
||||||
|
|
||||||
def add_feature(self, feat):
|
def add_feature(self, feat):
|
||||||
if feat not in self.features_enabled:
|
if feat not in self.features_enabled:
|
||||||
if feat not in ArgParser.feature_sets_available:
|
if feat not in ArgParser.feature_sets_available:
|
||||||
raise ValueError("Unknown ArgParser feature: %s" % feat)
|
raise ValueError("Unknown ArgParser feature: %s" % feat)
|
||||||
#resolve dependencies
|
#resolve dependencies
|
||||||
_, deps, atoms, atparse = ArgParser.feature_sets_available[feat]
|
_, deps, atoms, atparse, provides_virtual = ArgParser.feature_sets_available[feat]
|
||||||
for dep in deps:
|
for dep in deps:
|
||||||
self.add_feature(dep)
|
self.add_feature(dep)
|
||||||
for atom in atoms: # maybe check whether that atom has already been added sometimes in the future?
|
for atom in atoms: # maybe check whether that atom has already been added sometimes in the future?
|
||||||
|
|
@ -144,6 +288,8 @@ class ArgParser:
|
||||||
for methodname in atparse:
|
for methodname in atparse:
|
||||||
self.call_at_parse_list.append(methodname)
|
self.call_at_parse_list.append(methodname)
|
||||||
self.features_enabled.add(feat)
|
self.features_enabled.add(feat)
|
||||||
|
for feat_virt in provides_virtual:
|
||||||
|
self.features_enabled.add(feat_virt)
|
||||||
|
|
||||||
def add_argument(self, *args, **kwargs):
|
def add_argument(self, *args, **kwargs):
|
||||||
'''Add a custom argument directly to the standard library ArgParser object'''
|
'''Add a custom argument directly to the standard library ArgParser object'''
|
||||||
|
|
@ -153,7 +299,11 @@ class ArgParser:
|
||||||
self.args = self.ap.parse_args(*args, **kwargs)
|
self.args = self.ap.parse_args(*args, **kwargs)
|
||||||
if process_data:
|
if process_data:
|
||||||
for method in self.call_at_parse_list:
|
for method in self.call_at_parse_list:
|
||||||
getattr(self, method)()
|
try:
|
||||||
|
getattr(self, method)()
|
||||||
|
except ArgumentProcessingError:
|
||||||
|
err = sys.exc_info()[1]
|
||||||
|
self.ap.error(str(err))
|
||||||
return self.args
|
return self.args
|
||||||
|
|
||||||
def __getattr__(self, name):
|
def __getattr__(self, name):
|
||||||
|
|
@ -162,30 +312,24 @@ class ArgParser:
|
||||||
|
|
||||||
# Methods to initialise the related data structures:
|
# Methods to initialise the related data structures:
|
||||||
|
|
||||||
def _eval_background_epsmu(self): # feature: background
|
def _eval_const_background_epsmu(self): # feature: const_real_background
|
||||||
from .cymaterials import EpsMu, EpsMuGenerator
|
self.args.background = self.args.refractive_index
|
||||||
self.background_epsmu = EpsMu(self.args.refractive_index**2)
|
self._eval_analytical_background_epsmugen()
|
||||||
self.background_emg = EpsMuGenerator(self.background_epsmu)
|
|
||||||
|
def _eval_analytical_background_epsmugen(self): # feature: background_analytical
|
||||||
|
a = self.args
|
||||||
|
from .cymaterials import EpsMu
|
||||||
|
if isinstance(a.background, (float, complex)):
|
||||||
|
self.background_epsmu = EpsMu(a.background**2)
|
||||||
|
self.background_emg = self._add_emg(a.background)
|
||||||
|
|
||||||
def _eval_single_tmgen(self): # feature: single_particle
|
def _eval_single_tmgen(self): # feature: single_particle
|
||||||
a = self.args
|
a = self.args
|
||||||
from .cymaterials import EpsMuGenerator, lorentz_drude
|
from .cymaterials import EpsMuGenerator, lorentz_drude
|
||||||
from .cytmatrices import TMatrixGenerator
|
from .cytmatrices import TMatrixGenerator
|
||||||
if a.material in lorentz_drude.keys():
|
self.foreground_emg = self._add_emg(a.material)
|
||||||
self.foreground_emg = EpsMuGenerator(lorentz_drude[a.material])
|
self.tmgen = self._add_tmg((a.background, a.material, (a.radius, a.height, a.lMax_extend)))
|
||||||
else:
|
self.bspec = self._add_bspec(a.lMax)
|
||||||
try: lemat = float(a.material)
|
|
||||||
except ValueError:
|
|
||||||
try: lemat = complex(a.material)
|
|
||||||
except ValueError as ve:
|
|
||||||
raise ValueError("--material must be either a label such as 'Ag', 'Au', or a number") from ve
|
|
||||||
a.material = lemat
|
|
||||||
self.foreground_emg = EpsMuGenerator(EpsMu(a.material**2))
|
|
||||||
|
|
||||||
if a.height is None:
|
|
||||||
self.tmgen = TMatrixGenerator.sphere(self.background_emg, self.foreground_emg, a.radius)
|
|
||||||
else:
|
|
||||||
self.tmgen = TMatrixGenerator.cylinder(self.background_emg, self.foreground_emg, a.radius, a.height, lMax_extend = a.lMax_extend)
|
|
||||||
|
|
||||||
def _eval_single_omega(self): # feature: single_omega
|
def _eval_single_omega(self): # feature: single_omega
|
||||||
from .constants import eV, hbar
|
from .constants import eV, hbar
|
||||||
|
|
@ -205,7 +349,7 @@ class ArgParser:
|
||||||
import numpy as np
|
import numpy as np
|
||||||
from .constants import eV, hbar
|
from .constants import eV, hbar
|
||||||
eh = eV / hbar
|
eh = eV / hbar
|
||||||
self.omegas = [omega_eV * eh for omega_eV in self.args.eV]
|
self.omegas = [omega_eV * eh for omega_eV in flatten(self.args.eV)]
|
||||||
self.omega_max = max(om if isinstance(om, float) else max(om) for om in self.omegas)
|
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_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_singles = [om for om in self.omegas if isinstance(om, float)]
|
||||||
|
|
@ -225,7 +369,7 @@ class ArgParser:
|
||||||
l = len(self.args.basis_vectors)
|
l = len(self.args.basis_vectors)
|
||||||
if l != 2: raise ValueError('Two basis vectors must be specified (have %d)' % l)
|
if l != 2: raise ValueError('Two basis vectors must be specified (have %d)' % l)
|
||||||
from .qpms_c import lll_reduce
|
from .qpms_c import lll_reduce
|
||||||
self.direct_basis = lll_reduce(self.args.basis_vector, delta=1.)
|
self.direct_basis = lll_reduce(self.args.basis_vectors, delta=1.)
|
||||||
import numpy as np
|
import numpy as np
|
||||||
self.reciprocal_basis1 = np.linalg.inv(self.direct_basis)
|
self.reciprocal_basis1 = np.linalg.inv(self.direct_basis)
|
||||||
self.reciprocal_basis2pi = 2 * np.pi * self.reciprocal_basis1
|
self.reciprocal_basis2pi = 2 * np.pi * self.reciprocal_basis1
|
||||||
|
|
@ -255,3 +399,66 @@ class ArgParser:
|
||||||
a.theta = a.theta * pi2
|
a.theta = a.theta * pi2
|
||||||
a.phi = a.phi * pi2
|
a.phi = a.phi * pi2
|
||||||
|
|
||||||
|
def _process_multi_particle(self): # feature: multi_particle
|
||||||
|
a = self.args
|
||||||
|
self.tmspecs = {}
|
||||||
|
self.tmgens = {}
|
||||||
|
self.bspecs = {}
|
||||||
|
self.positions = {}
|
||||||
|
pos13, pos23, pos33 = False, False, False # used to
|
||||||
|
if len(a.position.keys()) == 0:
|
||||||
|
warnings.warn("No particle position (-p or +p) specified, assuming single particle in the origin / single particle per unit cell!")
|
||||||
|
a.position[None] = [(0.,0.,0.)]
|
||||||
|
for poslabel in a.position.keys():
|
||||||
|
try:
|
||||||
|
lMax = a.lMax.get(poslabel, False) or a.lMax[None]
|
||||||
|
radius = a.radius.get(poslabel, False) or a.radius[None]
|
||||||
|
# Height is "inherited" only together with radius
|
||||||
|
height = a.height.get(poslabel, None) if poslabel in a.radius.keys() else a.height.get(None, None)
|
||||||
|
if hasattr(a, 'lMax_extend'):
|
||||||
|
lMax_extend = a.lMax_extend.get(poslabel, False) or a.lMax_extend.get(None, False) or None
|
||||||
|
else:
|
||||||
|
lMax_extend = None
|
||||||
|
material = a.material.get(poslabel, False) or a.material[None]
|
||||||
|
except (TypeError, KeyError) as exc:
|
||||||
|
if poslabel is None:
|
||||||
|
raise ArgumentProcessingError("Unlabeled particles' positions (-p) specified, but some default particle properties are missing (--lMax, --radius, and --material have to be specified)") from exc
|
||||||
|
else:
|
||||||
|
raise ArgumentProcessingError(("Incomplete specification of '%s'-labeled particles: you must"
|
||||||
|
"provide at least ++lMax, ++radius, ++material arguments with the label, or the fallback arguments"
|
||||||
|
"--lMax, --radius, --material.")%(str(poslabel),)) from exc
|
||||||
|
tmspec = (a.background, material, (radius, height, lMax_extend))
|
||||||
|
self.tmspecs[poslabel] = tmspec
|
||||||
|
self.tmgens[poslabel] = self._add_tmg(tmspec)
|
||||||
|
self.bspecs[poslabel] = self._add_bspec(lMax)
|
||||||
|
poslist_cured = []
|
||||||
|
for pos in a.position[poslabel]:
|
||||||
|
if len(pos) == 1:
|
||||||
|
pos_cured = (0., 0., pos[0])
|
||||||
|
pos13 = True
|
||||||
|
elif len(pos) == 2:
|
||||||
|
pos_cured = (pos[0], pos[1], 0.)
|
||||||
|
pos23 = True
|
||||||
|
elif len(pos) == 3:
|
||||||
|
pos_cured = pos
|
||||||
|
pos33 = True
|
||||||
|
else:
|
||||||
|
raise argparse.ArgumentTypeError("Each -p / +p argument requires 1 to 3 cartesian coordinates")
|
||||||
|
poslist_cured.append(pos_cured)
|
||||||
|
self.positions[poslabel] = poslist_cured
|
||||||
|
if pos13 and pos23:
|
||||||
|
warnings.warn("Both 1D and 2D position specifications used. The former are interpreted as z coordinates while the latter as x, y coordinates")
|
||||||
|
|
||||||
|
def get_particles(self):
|
||||||
|
"""Creates a list of Particle instances that can be directly used in ScatteringSystem.create().
|
||||||
|
|
||||||
|
Assumes that self._process_multi_particle() has been already called.
|
||||||
|
"""
|
||||||
|
from .qpms_c import Particle
|
||||||
|
plist = []
|
||||||
|
for poslabel, poss in self.positions.items():
|
||||||
|
t = self.tmgens[poslabel]
|
||||||
|
bspec = self.bspecs[poslabel]
|
||||||
|
plist.extend([Particle(pos, t, bspec=bspec) for pos in poss])
|
||||||
|
return plist
|
||||||
|
|
||||||
|
|
|
||||||
|
|
@ -726,14 +726,14 @@ cdef class ScatteringSystem:
|
||||||
qpms_scatsys_periodic_build_translation_matrix_full(&target_view[0][0], self.s, wavenumber, &blochvector_c)
|
qpms_scatsys_periodic_build_translation_matrix_full(&target_view[0][0], self.s, wavenumber, &blochvector_c)
|
||||||
return target
|
return target
|
||||||
|
|
||||||
def translation_matrix_packed(self, double k, qpms_iri_t iri, J = QPMS_HANKEL_PLUS):
|
def translation_matrix_packed(self, cdouble wavenumber, qpms_iri_t iri, J = QPMS_HANKEL_PLUS):
|
||||||
self.check_s()
|
self.check_s()
|
||||||
cdef size_t rlen = self.saecv_sizes[iri]
|
cdef size_t rlen = self.saecv_sizes[iri]
|
||||||
cdef np.ndarray[np.complex_t, ndim=2] target = np.empty(
|
cdef np.ndarray[np.complex_t, ndim=2] target = np.empty(
|
||||||
(rlen,rlen),dtype=complex, order='C')
|
(rlen,rlen),dtype=complex, order='C')
|
||||||
cdef cdouble[:,::1] target_view = target
|
cdef cdouble[:,::1] target_view = target
|
||||||
qpms_scatsys_build_translation_matrix_e_irrep_packed(&target_view[0][0],
|
qpms_scatsys_build_translation_matrix_e_irrep_packed(&target_view[0][0],
|
||||||
self.s, iri, k, J)
|
self.s, iri, wavenumber, J)
|
||||||
return target
|
return target
|
||||||
|
|
||||||
property fullvec_psizes:
|
property fullvec_psizes:
|
||||||
|
|
@ -1020,6 +1020,9 @@ cdef class _ScatteringSystemAtOmega:
|
||||||
def translation_matrix_full(self, blochvector = None):
|
def translation_matrix_full(self, blochvector = None):
|
||||||
return self.ss_pyref.translation_matrix_full(wavenumber=self.wavenumber, blochvector=blochvector)
|
return self.ss_pyref.translation_matrix_full(wavenumber=self.wavenumber, blochvector=blochvector)
|
||||||
|
|
||||||
|
def translation_matrix_packed(self, iri, J = QPMS_HANKEL_PLUS):
|
||||||
|
return self.ss_pyref.translation_matrix_packed(wavenumber=self.wavenumber, iri=iri, J=J)
|
||||||
|
|
||||||
def scattered_E(self, scatcoeffvector_full, evalpos, bint alt=False, btyp=QPMS_HANKEL_PLUS):
|
def scattered_E(self, scatcoeffvector_full, evalpos, bint alt=False, btyp=QPMS_HANKEL_PLUS):
|
||||||
cdef qpms_bessel_t btyp_c = BesselType(btyp)
|
cdef qpms_bessel_t btyp_c = BesselType(btyp)
|
||||||
evalpos = np.array(evalpos, dtype=float, copy=False)
|
evalpos = np.array(evalpos, dtype=float, copy=False)
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue