diff --git a/misc/infiniterectlat-scatter.py b/misc/infiniterectlat-scatter.py new file mode 100755 index 0000000..72c6b57 --- /dev/null +++ b/misc/infiniterectlat-scatter.py @@ -0,0 +1,119 @@ +#!/usr/bin/env python3 + +import math +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.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") + + +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 +import warnings +from qpms.cybspec import BaseSpec +from qpms.cytmatrices import CTMatrix, TMatrixGenerator +from qpms.qpms_c import Particle, pgsl_ignore_error +from qpms.cymaterials import EpsMu, EpsMuGenerator, LorentzDrudeModel, lorentz_drude +from qpms.cycommon import DebugFlags, dbgmsg_enable +from qpms import FinitePointGroup, ScatteringSystem, BesselType, eV, hbar +eh = eV/hbar + +dbgmsg_enable(DebugFlags.INTEGRATION) + +a1 = ap.direct_basis[0] +a2 = ap.direct_basis[1] + +#Particle positions +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) + +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 +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 + +σ_ext_list = np.empty((a.N,), dtype=float) +σ_scat_list = np.empty((a.N,), 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ã) + + σ_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 + +σ_abs_list = σ_ext_list - σ_scat_list + +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 + ) +logging.info("Saved to %s" % outfile) + + +if a.plot or (a.plot_out is not None): + import matplotlib + matplotlib.use('pdf') + from matplotlib import pyplot as plt + + 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) + +exit(0) + diff --git a/qpms/CMakeLists.txt b/qpms/CMakeLists.txt index 525ca52..4a7f0cc 100644 --- a/qpms/CMakeLists.txt +++ b/qpms/CMakeLists.txt @@ -15,7 +15,7 @@ add_library (qpms SHARED translations.c tmatrices.c vecprint.c vswf.c wigner.c e ewaldsf.c pointgroups.c latticegens.c lattices2d.c gaunt.c error.c legendre.c symmetries.c vecprint.c bessel.c own_zgemm.c parsing.c scatsystem.c materials.c drudeparam_data.c - lll.c beyn.c + lll.c beyn.c trivialgroup.c ) use_c99() diff --git a/qpms/argproc.py b/qpms/argproc.py index ba0363f..adf8b9c 100644 --- a/qpms/argproc.py +++ b/qpms/argproc.py @@ -91,8 +91,9 @@ class ArgParser: # Methods to initialise the related data structures: def _eval_background_epsmu(self): # feature: background - from .cymaterials import EpsMu + from .cymaterials import EpsMu, EpsMuGenerator self.background_epsmu = EpsMu(self.args.refractive_index**2) + self.background_emg = EpsMuGenerator(self.background_epsmu) def _eval_single_tmgen(self): # feature: single_particle a = self.args @@ -110,9 +111,9 @@ class ArgParser: self.foreground_emg = EpsMuGenerator(EpsMu(a.material**2)) if a.height is None: - self.tmgen = TMatrixGenerator.sphere(self.background_epsmu, self.foreground_emg, a.radius) + self.tmgen = TMatrixGenerator.sphere(self.background_emg, self.foreground_emg, a.radius) else: - self.tmgen = TMatrixGenerator.cylinder(self.background_epsmu, self.foreground_emg, a.radius, a.height, lMax_extend = a.lMax_extend) + 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 from .constants import eV, hbar diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index 1143a15..1b51bd5 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -438,7 +438,7 @@ cdef class ScatteringSystem: assert(len(latticebasis) <= 3 and len(latticebasis) > 0) orig.lattice_dimension = len(latticebasis) for d in range(len(latticebasis)): - orig.per.lattice_basis[d] = {'x' : latticebasis[d][0], 'y' : latticebasis[d][1], 'z' : latticebasis[d][2]} + orig.per.lattice_basis[d] = {'x' : latticebasis[d][0], 'y' : latticebasis[d][1], 'z' : latticebasis[d][2] if len(latticebasis[d]) >= 3 else 0} else: orig.lattice_dimension = 0 ssw = qpms_scatsys_apply_symmetry(&orig, sym.rawpointer(), omega, &QPMS_TOLERANCE_DEFAULT) ss = ssw[0].ss @@ -495,6 +495,17 @@ cdef class ScatteringSystem: def __get__(self): self.check_s() return self.s[0].sym[0].nirreps + property lattice_dimension: + def __get__(self): + return self.s[0].lattice_dimension + + property unitcell_volume: + def __get__(self): + self.check_s() + if self.lattice_dimension: + return self.s[0].per.unitcell_volume + else: + return None def pack_vector(self, vect, iri): self.check_s() @@ -552,13 +563,24 @@ cdef class ScatteringSystem: self.s, iri, 0) return target_np - def translation_matrix_full(self, double k, J = QPMS_HANKEL_PLUS): + def translation_matrix_full(self, cdouble wavenumber, blochvector = None, J = QPMS_HANKEL_PLUS): self.check_s() cdef size_t flen = self.s[0].fecv_size cdef np.ndarray[np.complex_t, ndim=2] target = np.empty( (flen,flen),dtype=complex, order='C') cdef cdouble[:,::1] target_view = target - qpms_scatsys_build_translation_matrix_e_full(&target_view[0][0], self.s, k, J) + cdef cart3_t blochvector_c + if self.lattice_dimension == 0: + if blochvector is None: + qpms_scatsys_build_translation_matrix_e_full(&target_view[0][0], self.s, wavenumber, J) + else: raise ValueError("Can't use blochvector with non-periodic system") + else: + if blochvector is None: raise ValueError("Valid blochvector must be specified for periodic system") + else: + if J != QPMS_HANKEL_PLUS: + raise NotImplementedError("Translation operators based on other than Hankel+ functions not supperted in periodic systems") + blochvector_c = {'x': blochvector[0], 'y': blochvector[1], 'z': blochvector[2]} + qpms_scatsys_periodic_build_translation_matrix_full(&target_view[0][0], self.s, wavenumber, &blochvector_c) return target def translation_matrix_packed(self, double k, qpms_iri_t iri, J = QPMS_HANKEL_PLUS): @@ -637,7 +659,7 @@ cdef class ScatteringSystem: rank_tol, rank_min_sel, res_tol) if res == NULL: raise RuntimeError - cdef size_t neig = res[0].neig + cdef size_t neig = res[0].neig, i, j cdef size_t vlen = res[0].vlen # should be equal to self.s.fecv_size cdef np.ndarray[complex, ndim=1] eigval = np.empty((neig,), dtype=complex) @@ -676,6 +698,16 @@ cdef class ScatteringSystem: return retdict +cdef class _ScatteringSystemAtOmegaK: + ''' + Wrapper over the C qpms_scatsys_at_omega_k_t structure + ''' + cdef qpms_scatsys_at_omega_k_t sswk + cdef _ScatteringSystemAtOmega ssw_pyref + + cdef qpms_scatsys_at_omega_k_t *rawpointer(self): + return &self.sswk + cdef class _ScatteringSystemAtOmega: ''' Wrapper over the C qpms_scatsys_at_omega_t structure @@ -691,6 +723,11 @@ cdef class _ScatteringSystemAtOmega: self.ss_pyref.check_s() #TODO is there a way to disable the constructor outside this module? + + def ensure_finite(self): + if self.ssw[0].ss[0].lattice_dimension != 0: + raise NotImplementedError("Operation not supported for periodic systems") + def __dealloc__(self): if (self.ssw): qpms_scatsys_at_omega_free(self.ssw) @@ -711,9 +748,20 @@ cdef class _ScatteringSystemAtOmega: cdef qpms_scatsys_at_omega_t *rawpointer(self): return self.ssw - def scatter_solver(self, iri=None): + def scatter_solver(self, iri=None, k=None): self.check() - return ScatteringMatrix(self, iri) + cdef _ScatteringSystemAtOmegaK sswk # used only for periodic systems + if(self.ssw[0].ss[0].lattice_dimension == 0): + return ScatteringMatrix(self, iri=iri) + else: + if iri is not None: + raise NotImplementedError("Irrep decomposition not (yet) supported for periodic systems") + sswk = _ScatteringSystemAtOmegaK() + sswk.sswk.ssw = self.ssw + sswk.sswk.k[0] = k[0] + sswk.sswk.k[1] = k[1] + sswk.sswk.k[2] = k[2] + return ScatteringMatrix(ssw=self, sswk=sswk, iri=None) property fecv_size: def __get__(self): return self.ss_pyref.fecv_size @@ -723,8 +771,11 @@ cdef class _ScatteringSystemAtOmega: def __get__(self): return self.ss_pyref.irrep_names property nirreps: def __get__(self): return self.ss_pyref.nirreps + property wavenumber: + def __get__(self): return self.ssw[0].wavenumber + - def modeproblem_matrix_full(self): + def modeproblem_matrix_full(self, k=None): self.check() cdef size_t flen = self.ss_pyref.s[0].fecv_size cdef np.ndarray[np.complex_t, ndim=2] target = np.empty( @@ -735,6 +786,7 @@ cdef class _ScatteringSystemAtOmega: def modeproblem_matrix_packed(self, qpms_iri_t iri, version='pR'): self.check() + self.ensure_finite() cdef size_t rlen = self.saecv_sizes[iri] cdef np.ndarray[np.complex_t, ndim=2] target = np.empty( (rlen,rlen),dtype=complex, order='C') @@ -748,24 +800,35 @@ cdef class _ScatteringSystemAtOmega: qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial(&target_view[0][0], self.ssw, iri) return target + def translation_matrix_full(self, blochvector = None): + return self.ss_pyref.translation_matrix_full(wavenumber=self.wavenumber, blochvector=blochvector) + + cdef class ScatteringMatrix: ''' Wrapper over the C qpms_ss_LU structure that keeps the factorised mode problem matrix. ''' cdef _ScatteringSystemAtOmega ssw # Here we keep the reference to the parent scattering system + cdef _ScatteringSystemAtOmegaK sswk cdef qpms_ss_LU lu - def __cinit__(self, _ScatteringSystemAtOmega ssw, iri=None): + def __cinit__(self, _ScatteringSystemAtOmega ssw, sswk=None, iri=None): ssw.check() self.ssw = ssw - # TODO? pre-allocate the matrix with numpy to make it transparent? - if iri is None: - self.lu = qpms_scatsysw_build_modeproblem_matrix_full_LU( - NULL, NULL, ssw.rawpointer()) + if sswk is None: + ssw.ensure_finite() + # TODO? pre-allocate the matrix with numpy to make it transparent? + if iri is None: + self.lu = qpms_scatsysw_build_modeproblem_matrix_full_LU( + NULL, NULL, ssw.rawpointer()) + else: + self.lu = qpms_scatsysw_build_modeproblem_matrix_irrep_packed_LU( + NULL, NULL, ssw.rawpointer(), iri) else: - self.lu = qpms_scatsysw_build_modeproblem_matrix_irrep_packed_LU( - NULL, NULL, ssw.rawpointer(), iri) + # TODO check sswk validity + self.sswk = sswk + self.lu = qpms_scatsyswk_build_modeproblem_matrix_full_LU(NULL, NULL, self.sswk.rawpointer()) def __dealloc__(self): qpms_ss_LU_free(self.lu) diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 1f39ebe..1687d73 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -531,6 +531,7 @@ cdef extern from "scatsystem.h": qpms_tmatrix_operation_t op struct qpms_scatsys_periodic_info_t: cart3_t lattice_basis[3] + double unitcell_volume #etc. struct qpms_scatsys_t: int lattice_dimension @@ -616,8 +617,7 @@ cdef extern from "scatsystem.h": double k[3] cdouble *qpms_scatsyswk_build_modeproblem_motrix_full(cdouble *target, const qpms_scatsys_at_omega_k_t *sswk) cdouble *qpms_scatsys_periodic_build_translation_matrix_full(cdouble *target, const qpms_scatsys_t *ss, cdouble wavenumber, const cart3_t *wavevector) - cdouble *qpms_scatsyswk_build_translation_matrix_full(cdouble *target, const qpms_scatsys_at_omega_k_t *sswk) - qpms_ss_LU qpms_scatsyswk_build_modeproblem_matrix_full_LU(cdouble *target, int *target_piv, const qpms_scatsys_at_omega_t *sswk) + qpms_ss_LU qpms_scatsyswk_build_modeproblem_matrix_full_LU(cdouble *target, int *target_piv, const qpms_scatsys_at_omega_k_t *sswk) beyn_result_t *qpms_scatsys_periodic_find_eigenmodes(const qpms_scatsys_t *ss, const double *k, cdouble omega_centre, double omega_rr, double omega_ri, size_t contour_npoints, double rank_tol, size_t rank_sel_min, double res_tol) diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 73a5281..e3dc4d6 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -1266,7 +1266,7 @@ static inline complex double *qpms_scatsysw_scatsyswk_build_modeproblem_matrix_f { const complex double wavenumber = ssw->wavenumber; const qpms_scatsys_t *ss = ssw->ss; - qpms_ss_ensure_nonperiodic(ss); + qpms_ss_ensure_periodic(ss); const size_t full_len = ss->fecv_size; if(!target) QPMS_CRASHING_MALLOC(target, SQ(full_len) * sizeof(complex double));