Separate Ewald parameter for different frequencies

Former-commit-id: 0b9129fda0224411c2a1d8372fe715db4e071ecd

qpms/qpms_c.pyx

@@ -568,7 +568,10 @@ cdef class ScatteringSystem:
                 return None
 
     property eta:
-        """Ewald parameter η"""
+        """Ewald parameter η (only relevant for periodic systems)
+        
+        This is the default value that will be copied into
+        _ScatteringSystemAtOmega by __call__"""
         def __get__(self):
             self.check_s()
             if self.lattice_dimension:
@@ -723,7 +726,7 @@ cdef class ScatteringSystem:
                 self.s, iri, 0)
         return target_np
 
-    def translation_matrix_full(self, cdouble wavenumber, blochvector = None, J = QPMS_HANKEL_PLUS):
+    def translation_matrix_full(self, cdouble wavenumber, blochvector = None, J = QPMS_HANKEL_PLUS, double eta = 0):
         """Constructs full translation matrix of a scattering system.
 
         This method enables to use any wave number for the background medium ignoring the
@@ -741,6 +744,11 @@ cdef class ScatteringSystem:
         J : BesselType
             Optionally, one can replace Hankel functions of the first kind with different
             Bessel functions.
+
+        eta : float
+            Ewald parameter η, only relevant for periodic lattices; 
+            if set to 0 (default), the actual value is determined 
+            automatically.
         
         See Also
         --------
@@ -763,7 +771,7 @@ cdef class ScatteringSystem:
                     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]}
                 with pgsl_ignore_error(15):
-                    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, eta)
         return target
 
     def translation_matrix_packed(self, cdouble wavenumber, qpms_iri_t iri, J = QPMS_HANKEL_PLUS):
@@ -1025,6 +1033,20 @@ cdef class _ScatteringSystemAtOmega:
         def __get__(self): return self.ss_pyref.nirreps
     property wavenumber:
         def __get__(self): return self.ssw[0].wavenumber
+    property eta:
+        """Ewald parameter η (only relevant for periodic systems)"""
+        def __get__(self):
+            self.check_s()
+            if self.lattice_dimension:
+                return self.ssw[0].eta
+            else:
+                return None
+        def __set__(self, double eta):
+            self.check_s()
+            if self.lattice_dimension:
+                self.ssw[0].eta = eta
+            else:
+                raise AttributeError("Cannot set Ewald parameter for finite system") # different exception?
         
 
     def modeproblem_matrix_full(self, k=None):

qpms/qpms_cdefs.pxd

@@ -630,6 +630,7 @@ cdef extern from "scatsystem.h":
         cdouble omega
         qpms_epsmu_t medium
         cdouble wavenumber
+        double eta
     qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qpms_finite_group_t *sym,
             cdouble omega, const qpms_tolerance_spec_t *tol)
     qpms_scatsys_at_omega_t *qpms_scatsys_at_omega(const qpms_scatsys_t *ss, cdouble omega)
@@ -688,7 +689,7 @@ cdef extern from "scatsystem.h":
         const qpms_scatsys_at_omega_t *ssw
         double k[3]
     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)
+    cdouble *qpms_scatsys_periodic_build_translation_matrix_full(cdouble *target, const qpms_scatsys_t *ss, cdouble wavenumber, const cart3_t *wavevector, double eta)
     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, 

qpms/scatsystem.c

@@ -59,6 +59,13 @@ static inline void qpms_ss_ensure_nonperiodic_a(const qpms_scatsys_t *ss, const
   QPMS_ENSURE(ss->lattice_dimension == 0, "This method is applicable only to nonperiodic systems. Use %s instead.", s);
 }
 
+// Adjust Ewald parameter to avoid high-frequency breakdown; 
+// TODO make this actually do something (other than just copying ss's eta.)
+static inline double ss_adjusted_eta(const qpms_scatsys_t *ss, complex double omega) {
+	qpms_ss_ensure_periodic(ss);
+	return ss->per.eta;
+}
+
 // ------------ Stupid implementation of qpms_scatsys_apply_symmetry() -------------
 
 #define MIN(x,y) (((x)<(y))?(x):(y))
@@ -541,6 +548,9 @@ qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig,
   }
   
   ss->c = qpms_trans_calculator_init(lMax, normalisation);
+
+  ssw->eta = ss->lattice_dimension ? ss_adjusted_eta(ss, omega) : NAN;
+
   return ssw;
 }
 
@@ -594,6 +604,7 @@ qpms_scatsys_at_omega_t *qpms_scatsys_at_omega(const qpms_scatsys_t *ss,
   ssw->ss = ss;
   ssw->medium = qpms_epsmu_generator_eval(ss->medium, omega);
   ssw->wavenumber = qpms_wavenumber(omega, ssw->medium);
+  ssw->eta = ss->lattice_dimension ? ss_adjusted_eta(ss, omega) : NAN;
   QPMS_CRASHING_CALLOC(ssw->tm, ss->tm_count, sizeof(*ssw->tm));
   qpms_tmatrix_t **tmatrices_preop;
   QPMS_CRASHING_CALLOC(tmatrices_preop, ss->tmg_count, sizeof(*tmatrices_preop));
@@ -1156,7 +1167,7 @@ complex double *qpms_scatsyswk_build_translation_matrix_full(
   const qpms_scatsys_t *ss = ssw->ss;
   qpms_ss_ensure_periodic(ss);
   const cart3_t k_cart3 = cart3_from_double_array(sswk->k);
-  return qpms_scatsys_periodic_build_translation_matrix_full(target, ss, wavenumber, &k_cart3);
+  return qpms_scatsys_periodic_build_translation_matrix_full(target, ss, wavenumber, &k_cart3, ssw->eta);
 }
 
 complex double *qpms_scatsys_build_translation_matrix_e_full(
@@ -1201,17 +1212,17 @@ complex double *qpms_scatsys_build_translation_matrix_e_full(
 static inline int qpms_ss_ppair_W32xy(const qpms_scatsys_t *ss,
     qpms_ss_pi_t pdest, qpms_ss_pi_t psrc, complex double wavenumber, const cart2_t kvector,
     complex double *target, const ptrdiff_t deststride, const ptrdiff_t srcstride,
-    qpms_ewald_part parts) {
+    qpms_ewald_part parts, double eta) {
   const qpms_vswf_set_spec_t *srcspec = qpms_ss_bspec_pi(ss, psrc);
   const qpms_vswf_set_spec_t *destspec = qpms_ss_bspec_pi(ss, pdest);
 
   // This might be a bit arbitrary, roughly "copied" from Unitcell constructor. TODO review.
-  const double maxR = sqrt(ss->per.unitcell_volume) * 32;
+  const double maxR = sqrt(ss->per.unitcell_volume) * 64;
-  const double maxK = 1024 * 2 * M_PI / maxR;
+  const double maxK = 2048 * 2 * M_PI / maxR;
 
   return qpms_trans_calculator_get_trans_array_e32_e(ss->c,
       target, NULL /*err*/, destspec, deststride, srcspec, srcstride,
-      ss->per.eta, wavenumber, 
+      eta, wavenumber, 
       cart3xy2cart2(ss->per.lattice_basis[0]), cart3xy2cart2(ss->per.lattice_basis[1]),
       kvector, 
       cart2_substract(cart3xy2cart2(ss->p[pdest].pos), cart3xy2cart2(ss->p[psrc].pos)),
@@ -1221,21 +1232,23 @@ static inline int qpms_ss_ppair_W32xy(const qpms_scatsys_t *ss,
 static inline int qpms_ss_ppair_W(const qpms_scatsys_t *ss,
     qpms_ss_pi_t pdest, qpms_ss_pi_t psrc, complex double wavenumber, const double wavevector[],
     complex double *target, const ptrdiff_t deststride, const ptrdiff_t srcstride,
-    qpms_ewald_part parts) {
+    qpms_ewald_part parts, double eta) {
   if(ss->lattice_dimension == 2 && // Currently, we can only the xy-plane
       !ss->per.lattice_basis[0].z && !ss->per.lattice_basis[1].z &&
       !wavevector[2]) 
     return qpms_ss_ppair_W32xy(ss, pdest, psrc, wavenumber, cart2_from_double_array(wavevector), 
-              target, deststride, srcstride, parts);
+              target, deststride, srcstride, parts, eta);
   else 
     QPMS_NOT_IMPLEMENTED("Only 2D xy-lattices currently supported");
 }
 
 complex double *qpms_scatsys_periodic_build_translation_matrix_full(
     complex double *target, const qpms_scatsys_t *ss,
-    complex double wavenumber, const cart3_t *wavevector) {
+    complex double wavenumber, const cart3_t *wavevector, double eta) {
   QPMS_UNTESTED;
   qpms_ss_ensure_periodic(ss);
+  if (eta == 0 || isnan(eta))
+    eta = ss->per.eta;
   const size_t full_len = ss->fecv_size;
   if(!target)
     QPMS_CRASHING_MALLOC(target, SQ(full_len) * sizeof(complex double));
@@ -1248,7 +1261,7 @@ complex double *qpms_scatsys_periodic_build_translation_matrix_full(
       for (qpms_ss_pi_t ps = 0; ps < ss->p_count; ++ps) {
         QPMS_ENSURE_SUCCESS(qpms_ss_ppair_W32xy(ss, pd, ps, wavenumber, cart3xy2cart2(*wavevector), 
               target + deststride * ss->fecv_pstarts[pd] + srcstride * ss->fecv_pstarts[ps],
-              deststride, srcstride, QPMS_EWALD_FULL));
+              deststride, srcstride, QPMS_EWALD_FULL, eta));
       }
   } else 
     QPMS_NOT_IMPLEMENTED("Only 2D xy-lattices currently supported");
@@ -1294,7 +1307,7 @@ static inline complex double *qpms_scatsysw_scatsyswk_build_modeproblem_matrix_f
         } else { // periodic case 
           QPMS_ENSURE_SUCCESS(qpms_ss_ppair_W(ss, piR, piC, wavenumber, k,
                 tmp /*target*/, bspecC->n /*deststride*/, 1 /*srcstride*/,
-                QPMS_EWALD_FULL));
+                QPMS_EWALD_FULL, ssw->eta));
         }
 
         cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,

qpms/scatsystem.h

@@ -147,7 +147,12 @@ typedef struct qpms_scatsys_periodic_info_t {
 	 */
 	double unitcell_volume;
 
-	/// Ewald parameter \f$ \eta \f$.
+	/// Default Ewald parameter \f$ \eta \f$.
+	/** Normally, this just gets copied into qpms_scatsys_at_omega_t,
+	 * which is then used in the Ewald sums.
+	 * However, for higher frequencies it must be adjusted to avoid
+	 * numerical instability.
+	 */
 	double eta;
 } qpms_scatsys_periodic_info_t;
 
@@ -241,6 +246,7 @@ typedef struct qpms_scatsys_at_omega_t {
 	complex double omega; ///< Angular frequency
 	qpms_epsmu_t medium; ///< Background medium optical properties at the given frequency
 	complex double wavenumber; ///< Background medium wave number
+	double eta; ///< Ewald parameter \f$ \eta \f$.
 } qpms_scatsys_at_omega_t;
 
 
@@ -501,7 +507,8 @@ complex double *qpms_scatsys_periodic_build_translation_matrix_full(
 		complex double *target,
 		const qpms_scatsys_t *ss,
 		complex double wavenumber, ///< Wave number to use in the translation matrix.
-		const cart3_t *wavevector ///< Wavevector / pseudomomentum in cartesian coordinates.
+		const cart3_t *wavevector, ///< Wavevector / pseudomomentum in cartesian coordinates.
+		double eta ///< Ewald parameter eta. Pass 0 or NaN to use the default value in \a ss.
 		);
 
 /// Global translation matrix.