Make infinite rect lat. / real freq. SVD work with ScatteringSystem.

Former-commit-id: 56beaeca0b44ca9087208c5e27c007d492572302

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 
-#from qpms.symmetries import point_group_info # TODO
 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)
+pp = Particle(orig_xy[0][0], ap.tmgen, bspec=bspec)
-par = [pp]
 
-u = unitcell(a1, a2, par, refractive_index=a.refractive_index)
+ss, ssw = ScatteringSystem.create([pp], ap.background_emg, omegas[0], latticebasis=ap.direct_basis)
-eta = (np.pi / u.Area)**.5
+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,))
+SVs = [None] * ss1.nirreps
-beta = np.array([0.,0.])
+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)
+        ImTW = ssw.modeproblem_matrix_full(k)
-    SVs[i] = np.linalg.svd(ImTW, compute_uv = False)
+    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):
+    cc = plt.rcParams['axes.prop_cycle']()
-        ax.plot(omegas/eh, SVs[:,-1-i])
+    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)

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

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)
+      retval = gsl_sf_bessel_jl_steed_array(lmax, x, tmparr);
-                                     : gsl_sf_bessel_jl_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)
+      retval = gsl_sf_bessel_jl_steed_array(lmax, x, tmparr);
-                                     : gsl_sf_bessel_jl_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
+    else {
-      QPMS_PR_ERROR("Amos's zbes%c failed with ierr == %d.",
+      if (retry_counter < 5) {
-         kindchar, (int) ierr);
+        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));
+    if ( NULL == (c->akn = realloc(c->akn, sizeof(double) * akn_index(lMax + 2, 0))))
-    // QPMS_CRASHING_REALLOC(c->bkn, sizeof(complex double) * bkn_index(lMax + 1, 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);

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 = self._sswk(k)
-            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)
 
+    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'):

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,