diff --git a/misc/infiniterectlat-k0realfreqsvd.py b/misc/infiniterectlat-k0realfreqsvd.py index 251a720..4e8f8d6 100755 --- a/misc/infiniterectlat-k0realfreqsvd.py +++ b/misc/infiniterectlat-k0realfreqsvd.py @@ -8,6 +8,7 @@ ap.add_argument("-o", "--output", type=str, required=False, help='output path (i 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("-s", "--singular_values", type=int, default=10, help="Number of singular values to plot") +ap.add_argument("--D2", action='store_true', help="Use D2h symmetry even if the x and y periods are equal") a=ap.parse_args() @@ -17,7 +18,7 @@ logging.basicConfig(format='%(asctime)s %(message)s', level=logging.INFO) px, py = a.period #Important! The particles are supposed to be of D2h/D4h symmetry -thegroup = 'D4h' if px == py else 'D2h' +thegroup = 'D4h' if px == py and not a.D2 else 'D2h' 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) @@ -28,14 +29,14 @@ 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 -from qpms.cyunitcell import unitcell -#from qpms.symmetries import point_group_info # TODO +from qpms.symmetries import point_group_info eh = eV/hbar # not used; TODO: @@ -46,7 +47,9 @@ irrep_labels = {"B2''":"$B_2''$", "A2''":"$A_2''$", "B1''":"$B_1''$", "A2'":"$A_2'$", - "B1'":"$B_1'$"} + "B1'":"$B_1'$", + "E'":"$E'$", + "E''":"$E''$",} dbgmsg_enable(DebugFlags.INTEGRATION) @@ -65,23 +68,30 @@ logging.info("%d frequencies from %g to %g eV" % (len(omegas), omegas[0]/eh, ome bspec = BaseSpec(lMax = a.lMax) nelem = len(bspec) # The parameters here should probably be changed (needs a better qpms_c.Particle implementation) -pp = Particle(orig_xy[0][0], tmgen=ap.tmgen, bspec=bspec) -par = [pp] +pp = Particle(orig_xy[0][0], ap.tmgen, bspec=bspec) -u = unitcell(a1, a2, par, refractive_index=a.refractive_index) -eta = (np.pi / u.Area)**.5 +ss, ssw = ScatteringSystem.create([pp], ap.background_emg, omegas[0], latticebasis=ap.direct_basis) +k = np.array([0.,0.,0]) +# Auxillary finite scattering system for irrep decomposition, quite a hack +ss1, ssw1 = ScatteringSystem.create([pp], ap.background_emg, omegas[0],sym=FinitePointGroup(point_group_info[thegroup])) wavenumbers = np.empty(omegas.shape) -SVs = np.empty(omegas.shape+(nelem,)) -beta = np.array([0.,0.]) +SVs = [None] * ss1.nirreps +for iri in range(ss1.nirreps): + SVs[iri] = np.empty(omegas.shape+(ss1.saecv_sizes[iri],)) for i, omega in enumerate(omegas): - wavenumbers[i] = ap.background_epsmu.k(omega).real # Currently, ScatteringSystem does not "remember" frequency nor wavenumber + ssw = ss(omega) + wavenumbers[i] = ssw.wavenumber.real + if ssw.wavenumber.imag: + warnings.warn("Non-zero imaginary wavenumber encountered") with pgsl_ignore_error(15): # avoid gsl crashing on underflow; maybe not needed - ImTW = u.evaluate_ImTW(eta, omega, beta) - SVs[i] = np.linalg.svd(ImTW, compute_uv = False) + ImTW = ssw.modeproblem_matrix_full(k) + for iri in range(ss1.nirreps): + ImTW_packed = ss1.pack_matrix(ImTW, iri) + SVs[iri][i] = np.linalg.svd(ImTW_packed, compute_uv = False) outfile = defaultprefix + ".npz" if a.output is None else a.output -np.savez(outfile, meta=vars(a), omegas=omegas, wavenumbers=wavenumbers, SVs=SVs, unitcell_area=u.Area +np.savez(outfile, meta=vars(a), omegas=omegas, wavenumbers=wavenumbers, SVs=np.concatenate(SVs, axis=-1), irrep_names=ss1.irrep_names, irrep_sizes=ss1.saecv_sizes, unitcell_area=ss.unitcell_volume ) logging.info("Saved to %s" % outfile) @@ -93,10 +103,18 @@ if a.plot or (a.plot_out is not None): fig = plt.figure() ax = fig.add_subplot(111) - for i in range(a.singular_values): - ax.plot(omegas/eh, SVs[:,-1-i]) + cc = plt.rcParams['axes.prop_cycle']() + for iri in range(ss1.nirreps): + cargs = next(cc) + nlines = min(a.singular_values, ss1.saecv_sizes[iri]) + for i in range(nlines): + ax.plot(omegas/eh, SVs[iri][:,-1-i], + label= None if i else irrep_labels[ss1.irrep_names[iri]], + **cargs) + ax.set_ylim([0,1.1]) ax.set_xlabel('$\hbar \omega / \mathrm{eV}$') ax.set_ylabel('Singular values') + ax.legend() plotfile = defaultprefix + ".pdf" if a.plot_out is None else a.plot_out fig.savefig(plotfile) diff --git a/qpms/argproc.py b/qpms/argproc.py index 5f4ccf5..ce5b12d 100644 --- a/qpms/argproc.py +++ b/qpms/argproc.py @@ -26,7 +26,8 @@ class ArgParser: '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')), 'single_frequency_eV': lambda ap: ap.add_argument("-f", "--eV", type=float, required=True, help='radiation angular frequency in eV'), - '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')), + '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')), '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'), @@ -46,7 +47,7 @@ class ArgParser: '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',)), '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", (), ('seq_frequency_eV',), ('_eval_omega_seq',)), + 'omega_seq': ("Equidistant real frequency range with possibility of adding individual frequencies", (), ('seq_frequency_eV', 'multiple_frequency_eV_optional',), ('_eval_omega_seq',)), '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',), ()), @@ -126,6 +127,9 @@ class ArgParser: from .constants import eV, hbar start, step, stop = self.args.eV_seq self.omegas = np.arange(start, stop, step) + if self.args.eV: + self.omegas = np.concatenate((self.omegas, np.array(self.args.eV))) + self.omegas.sort() self.omegas *= eV/hbar def _eval_lattice2d(self): # feature: lattice2d diff --git a/qpms/bessel.c b/qpms/bessel.c index 6c8dc80..b22d4be 100644 --- a/qpms/bessel.c +++ b/qpms/bessel.c @@ -48,42 +48,27 @@ static inline complex double spherical_yn(qpms_l_t l, complex double z) { } #endif -// Don't use Stead algorithm from GSL above this threshold (it fails for x's around 1000000 and higher) -static const double stead_threshold = 1000; - // There is a big issue with gsl's precision of spherical bessel function; these have to be implemented differently qpms_errno_t qpms_sph_bessel_realx_fill(qpms_bessel_t typ, qpms_l_t lmax, double x, complex double *result_array) { int retval; double tmparr[lmax+1]; switch(typ) { case QPMS_BESSEL_REGULAR: - retval = (x < stead_threshold) ? gsl_sf_bessel_jl_steed_array(lmax, x, tmparr) - : gsl_sf_bessel_jl_array(lmax, x, tmparr); + retval = gsl_sf_bessel_jl_steed_array(lmax, x, tmparr); for (int l = 0; l <= lmax; ++l) result_array[l] = tmparr[l]; return retval; break; case QPMS_BESSEL_SINGULAR: //FIXME: is this precise enough? Would it be better to do it one-by-one? - if(QPMS_UNLIKELY(x == 0)) { - for (int l = 0; l <= lmax; ++l) - result_array[l] = NAN; - return QPMS_ESING; // GSL would have returned GSL_EDOM without setting NANs. - } retval = gsl_sf_bessel_yl_array(lmax,x,tmparr); for (int l = 0; l <= lmax; ++l) result_array[l] = tmparr[l]; return retval; break; case QPMS_HANKEL_PLUS: case QPMS_HANKEL_MINUS: - retval = (x < stead_threshold) ? gsl_sf_bessel_jl_steed_array(lmax, x, tmparr) - : gsl_sf_bessel_jl_array(lmax, x, tmparr); + retval = gsl_sf_bessel_jl_steed_array(lmax, x, tmparr); for (int l = 0; l <= lmax; ++l) result_array[l] = tmparr[l]; if(retval) return retval; retval = gsl_sf_bessel_yl_array(lmax, x, tmparr); - if(QPMS_UNLIKELY(x == 0)) { - for (int l = 0; l <= lmax; ++l) - result_array[l] += I * NAN; - return QPMS_ESING; // GSL would have returned GSL_EDOM without setting NANs. - } if (typ==QPMS_HANKEL_PLUS) for (int l = 0; l <= lmax; ++l) result_array[l] += I * tmparr[l]; else @@ -91,9 +76,10 @@ qpms_errno_t qpms_sph_bessel_realx_fill(qpms_bessel_t typ, qpms_l_t lmax, double return retval; break; default: - QPMS_INVALID_ENUM(typ); + abort(); + //return GSL_EDOM; } - QPMS_WTF; + assert(0); } // TODO DOC @@ -106,6 +92,8 @@ qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, complex doub else if (fpclassify(creal(x)) == FP_INFINITE) for(qpms_l_t l = 0; l <= lmax; ++l) res[l] = INFINITY + I * INFINITY; else { +try_again: ; + int retry_counter = 0; const DOUBLE_PRECISION_t zr = creal(x), zi = cimag(x), fnu = 0.5; const INTEGER_t n = lmax + 1, kode = 1 /* No exponential scaling */; DOUBLE_PRECISION_t cyr[n], cyi[n]; @@ -133,7 +121,7 @@ qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, complex doub } break; default: - QPMS_INVALID_ENUM(typ); + QPMS_WTF; } // TODO check for underflows? (nz != 0) if (ierr == 0 || ierr == 3) { @@ -145,9 +133,15 @@ qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, complex doub kindchar); return QPMS_SUCCESS; //TODO maybe something else if ierr == 3 } - else - QPMS_PR_ERROR("Amos's zbes%c failed with ierr == %d.", - kindchar, (int) ierr); + else { + if (retry_counter < 5) { + QPMS_WARN("Amos's zbes%c failed with ierr == %d (lMax = %d, x = %+.16g%+.16gi). Retrying.\n", + kindchar, (int) ierr, lmax, creal(x), cimag(x)); + ++retry_counter; + goto try_again; + } else QPMS_PR_ERROR("Amos's zbes%c failed with ierr == %d (lMax = %d, x = %+.16g%+.16gi).", + kindchar, (int) ierr, lmax, creal(x), cimag(x)); + } } return QPMS_SUCCESS; } @@ -165,8 +159,10 @@ static inline qpms_errno_t qpms_sbessel_calculator_ensure_lMax(qpms_sbessel_calc if (lMax <= c->lMax) return QPMS_SUCCESS; else { - QPMS_CRASHING_REALLOC(c->akn, sizeof(double) * akn_index(lMax + 2, 0)); - // QPMS_CRASHING_REALLOC(c->bkn, sizeof(complex double) * bkn_index(lMax + 1, 0)); + if ( NULL == (c->akn = realloc(c->akn, sizeof(double) * akn_index(lMax + 2, 0)))) + abort(); + //if ( NULL == (c->bkn = realloc(c->bkn, sizeof(complex double) * bkn_index(lMax + 1, 0)))) + // abort(); for(qpms_l_t n = c->lMax+1; n <= lMax + 1; ++n) for(qpms_l_t k = 0; k <= n; ++k) c->akn[akn_index(n,k)] = exp(lgamma(n + k + 1) - k*M_LN2 - lgamma(k + 1) - lgamma(n - k + 1)); @@ -177,7 +173,8 @@ static inline qpms_errno_t qpms_sbessel_calculator_ensure_lMax(qpms_sbessel_calc } complex double qpms_sbessel_calc_h1(qpms_sbessel_calculator_t *c, qpms_l_t n, complex double x) { - QPMS_ENSURE_SUCCESS(qpms_sbessel_calculator_ensure_lMax(c, n)); + if(QPMS_SUCCESS != qpms_sbessel_calculator_ensure_lMax(c, n)) + abort(); complex double z = I/x; complex double result = 0; for (qpms_l_t k = n; k >= 0; --k) @@ -190,7 +187,8 @@ complex double qpms_sbessel_calc_h1(qpms_sbessel_calculator_t *c, qpms_l_t n, co qpms_errno_t qpms_sbessel_calc_h1_fill(qpms_sbessel_calculator_t * const c, const qpms_l_t lMax, const complex double x, complex double * const target) { - QPMS_ENSURE_SUCCESS(qpms_sbessel_calculator_ensure_lMax(c, lMax)); + if(QPMS_SUCCESS != qpms_sbessel_calculator_ensure_lMax(c, lMax)) + abort(); memset(target, 0, sizeof(complex double) * lMax); complex double kahancomp[lMax]; memset(kahancomp, 0, sizeof(complex double) * lMax); diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index a9c97d2..77c98c2 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -821,13 +821,17 @@ cdef class _ScatteringSystemAtOmega: 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] + sswk = self._sswk(k) return ScatteringMatrix(ssw=self, sswk=sswk, iri=None) + def _sswk(self, k): + cdef _ScatteringSystemAtOmegaK 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 sswk + property fecv_size: def __get__(self): return self.ss_pyref.fecv_size property saecv_sizes: @@ -846,7 +850,12 @@ cdef class _ScatteringSystemAtOmega: cdef np.ndarray[np.complex_t, ndim=2] target = np.empty( (flen,flen),dtype=complex, order='C') cdef cdouble[:,::1] target_view = target - qpms_scatsysw_build_modeproblem_matrix_full(&target_view[0][0], self.ssw) + cdef _ScatteringSystemAtOmegaK sswk + if k is not None: + sswk = self._sswk(k) + qpms_scatsyswk_build_modeproblem_matrix_full(&target_view[0][0], &sswk.sswk) + else: + qpms_scatsysw_build_modeproblem_matrix_full(&target_view[0][0], self.ssw) return target def modeproblem_matrix_packed(self, qpms_iri_t iri, version='pR'): diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index a08f41d..cb5ff75 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -618,7 +618,7 @@ cdef extern from "scatsystem.h": struct qpms_scatsys_at_omega_k_t: const qpms_scatsys_at_omega_t *ssw double k[3] - cdouble *qpms_scatsyswk_build_modeproblem_motrix_full(cdouble *target, const qpms_scatsys_at_omega_k_t *sswk) + cdouble *qpms_scatsyswk_build_modeproblem_matrix_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) 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,