Adjust eta for each (omega, k) pair to prevent hi-freq breakdown.

This effectively reverts 0b9129, because it does not make to define
Ewald parameter for just each frequency after all...


Former-commit-id: 6a3df5ecc1eecd6c120a74c70df5b747d593aae3

qpms/qpms_c.pyx

@@ -963,6 +963,13 @@ cdef class _ScatteringSystemAtOmegaK:
 
     cdef qpms_scatsys_at_omega_k_t *rawpointer(self):
         return &self.sswk
+    
+    property eta:
+        """Ewald parameter η"""
+        def __get__(self):
+            return self.sswk.eta
+        def __set__(self, double eta):
+            self.sswk.eta = eta
 
 cdef class _ScatteringSystemAtOmega:
     '''
@@ -1021,6 +1028,7 @@ cdef class _ScatteringSystemAtOmega:
         sswk.sswk.k[0] = k[0]
         sswk.sswk.k[1] = k[1]
         sswk.sswk.k[2] = k[2]
+        sswk.eta = qpms_ss_adjusted_eta(self.ssw[0].ss, self.ssw[0].wavenumber, sswk.sswk.k)
         return sswk
 
     property fecv_size: 
@@ -1033,20 +1041,6 @@ 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,7 +630,6 @@ 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,6 +687,7 @@ cdef extern from "scatsystem.h":
     struct qpms_scatsys_at_omega_k_t:
         const qpms_scatsys_at_omega_t *ssw
         double k[3]
+        double eta
     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, 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)
@@ -703,6 +703,8 @@ cdef extern from "scatsystem.h":
             const cdouble *f_excitation_vector_full, cart3_t where)
     ccart3_t qpms_scatsysw_scattered_E__alt(const qpms_scatsys_at_omega_t *ssw, qpms_bessel_t btyp,
             const cdouble *f_excitation_vector_full, cart3_t where)
+    double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, cdouble wavenumber, const double *wavevector);
+
 
 cdef extern from "ewald.h":
     struct qpms_csf_result:

qpms/scatsystem.c

@@ -23,6 +23,7 @@
 #include "kahansum.h"
 #include "tolerances.h"
 #include "beyn.h"
+#include "tiny_inlines.h"
 
 #ifdef QPMS_SCATSYSTEM_USE_OWN_BLAS
 #include "qpmsblas.h"
@@ -31,11 +32,15 @@
 #define SERIAL_ZGEMM cblas_zgemm
 #endif
 
-#define SQ(x) ((x)*(x))
 #define QPMS_SCATSYS_LEN_RTOL 1e-13
 #define QPMS_SCATSYS_TMATRIX_ATOL 1e-12
 #define QPMS_SCATSYS_TMATRIX_RTOL 1e-12
 
+// This is used in adjustment of Ewald parameter to avoid high frequency breakdown.
+// Very roughly, the value of 16 should lead to allowing terms containing incomplete Gamma
+// functions with magnitudes around exp(16) == 8.9e6
+static const double QPMS_SCATSYS_EWALD_MAX_EXPONENT = 16.;
+
 long qpms_scatsystem_nthreads_default = 4;
 long qpms_scatsystem_nthreads_override = 0;
 
@@ -59,18 +64,23 @@ 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; 
+// Adjust Ewald parameter to avoid high-frequency breakdown
-// TODO make this actually do something (other than just copying ss's eta.)
+double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, complex double wavenumber, const double k[3]) {
-static inline double ss_adjusted_eta(const qpms_scatsys_t *ss, complex double omega) {
 	qpms_ss_ensure_periodic(ss);
-	return ss->per.eta;
+	const double eta_default = ss->per.eta;
+	// FIXME here we silently assume that k lies in the first Brillioun zone, we should ensure that.
+	const double k2 = k[0]*k[0] + k[1]*k[1] + k[2] * k[2]; 
+	const double kappa2 = SQ(cabs(wavenumber)); // maybe creal would be enough
+  if(kappa2 < k2) // This should happen only for pretty low frequencies
+    return eta_default;
+	const qpms_l_t maxj = ss->c->lMax; // Based on ewald.c:301, note that VSWF (c's) lMax is already half of corresponding translation matrix Ewald factors' (c->e3c's) lMax
+	const double mina = 0.5 * (ss->lattice_dimension - 1) - maxj; // minimum incomplete Gamma first argument, based on ewald.c:301; CHECKME whether this is fine also for 3D lattice
+  const double eta_min = sqrt(fabs((kappa2 - k2) * (mina - 1.) / QPMS_SCATSYS_EWALD_MAX_EXPONENT));
+  return MAX(eta_default, eta_min);
 }
 
 // ------------ Stupid implementation of qpms_scatsys_apply_symmetry() -------------
 
-#define MIN(x,y) (((x)<(y))?(x):(y))
-#define MAX(x,y) (((x)>(y))?(x):(y))
-
 // The following functions are just to make qpms_scatsys_apply_symmetry more readable.
 // They are not to be used elsewhere, as they do not perform any array capacity checks etc.
 
@@ -549,8 +559,6 @@ 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;
 }
 
@@ -604,7 +612,6 @@ 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));
@@ -1167,7 +1174,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, ssw->eta);
+  return qpms_scatsys_periodic_build_translation_matrix_full(target, ss, wavenumber, &k_cart3, sswk->eta);
 }
 
 complex double *qpms_scatsys_build_translation_matrix_e_full(
@@ -1247,8 +1254,12 @@ complex double *qpms_scatsys_periodic_build_translation_matrix_full(
     complex double wavenumber, const cart3_t *wavevector, double eta) {
   QPMS_UNTESTED;
   qpms_ss_ensure_periodic(ss);
-  if (eta == 0 || isnan(eta))
+  if (eta == 0 || isnan(eta)) {
-    eta = ss->per.eta;
+    double tmp[3];
+    cart3_to_double_array(tmp, *wavevector);
+    eta = qpms_ss_adjusted_eta(ss, wavenumber, tmp);
+  }
+
   const size_t full_len = ss->fecv_size;
   if(!target)
     QPMS_CRASHING_MALLOC(target, SQ(full_len) * sizeof(complex double));
@@ -1273,7 +1284,8 @@ static inline complex double *qpms_scatsysw_scatsyswk_build_modeproblem_matrix_f
     /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated.
     complex double *target,
     const qpms_scatsys_at_omega_t *ssw,
-    const double k[] // NULL if non-periodic
+    const double k[], // NULL if non-periodic
+    const double eta // ignored if non-periodic
     )
 {
   const complex double wavenumber = ssw->wavenumber;
@@ -1307,7 +1319,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, ssw->eta));
+                QPMS_EWALD_FULL, eta));
         }
 
         cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
@@ -1334,14 +1346,14 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_full(
     complex double *target, const qpms_scatsys_at_omega_t *ssw) {
   qpms_ss_ensure_nonperiodic_a(ssw->ss, "qpms_scatsyswk_build_modeproblem_matrix_full()");
   return qpms_scatsysw_scatsyswk_build_modeproblem_matrix_full(
-      target, ssw, NULL);
+      target, ssw, NULL, NAN);
 }
 
 complex double *qpms_scatsyswk_build_modeproblem_matrix_full(
     complex double *target, const qpms_scatsys_at_omega_k_t *sswk) 
 {
   qpms_ss_ensure_periodic_a(sswk->ssw->ss, "qpms_scatsysw_build_modeproblem_matrix_full()");
-  return qpms_scatsysw_scatsyswk_build_modeproblem_matrix_full(target, sswk->ssw, sswk->k);
+  return qpms_scatsysw_scatsyswk_build_modeproblem_matrix_full(target, sswk->ssw, sswk->k, sswk->eta);
 }
 
 
@@ -2246,7 +2258,8 @@ static int qpms_scatsys_periodic_eval_Beyn_ImTW(complex double *target,
   QPMS_ENSURE(ssw != NULL, "qpms_scatsys_at_omega() returned NULL");
   qpms_scatsys_at_omega_k_t sswk = {
     .ssw = ssw,
-    .k = {p->k[0], p->k[1], p->k[2]}
+    .k = {p->k[0], p->k[1], p->k[2]},
+    .eta = qpms_ss_adjusted_eta(p->ss, ssw->wavenumber, p->k)
   };
   QPMS_ASSERT(m == p->ss->fecv_size);
   QPMS_ENSURE(NULL != 

qpms/scatsystem.h

@@ -246,7 +246,6 @@ 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;
 
 
@@ -489,6 +488,7 @@ complex double *qpms_scatsys_scatter_solve(
 typedef struct qpms_scatsys_at_omega_k_t {
 	const qpms_scatsys_at_omega_t *ssw;
 	double k[3]; ///< The k-vector's cartesian coordinates.
+	double eta; ///< Ewald parameter η.
 } qpms_scatsys_at_omega_k_t;
 
 /// Creates the full \f$ (I - WS) \f$ matrix of the periodic scattering system. 
@@ -769,6 +769,10 @@ ccart3_t qpms_scatsysw_scattered_E__alt(
 		cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
 		);
 
+/// Adjusted Ewadl parameter to avoid high-frequency breakdown.
+// TODO DOC
+double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, complex double wavenumber, const double wavevector[3]);
+
 #if 0
 /** Evaluates partial scattered fields (corresponding to a given irrep-reduced excitation vector)
  * at a given point.