diff --git a/qpms/cyewaldtest.pyx b/qpms/cyewaldtest.pyx index 4d545b1..c88215b 100644 --- a/qpms/cyewaldtest.pyx +++ b/qpms/cyewaldtest.pyx @@ -5,7 +5,7 @@ import numpy as np cdef extern from "ewald.h": void ewald3_2_sigma_long_Delta(cdouble *target, double *err, int maxn, cdouble x, int xbranch, cdouble z) void ewald3_2_sigma_long_Delta_series(cdouble *target, double *err, int maxn, cdouble x, int xbranch, cdouble z) - void ewald3_2_sigma_long_Delta_recurrent(cdouble *target, double *err, int maxn, cdouble x, int xbranch, cdouble z) + void ewald3_2_sigma_long_Delta_recurrent(cdouble *target, double *err, int maxn, cdouble x, int xbranch, cdouble z, bint bigimz) int complex_gamma_inc_e(double a, cdouble x, int xbranch, qpms_csf_result *result) def e32_Delta(int maxn, cdouble x, cdouble z, int xbranch = 0, get_err=True, method='auto'): @@ -17,7 +17,9 @@ def e32_Delta(int maxn, cdouble x, cdouble z, int xbranch = 0, get_err=True, met err_np = np.empty((maxn+1,), order='C') err_view = err_np if method == 'recurrent': - ewald3_2_sigma_long_Delta_recurrent(&target_view[0], &err_view[0] if get_err else NULL, maxn, x, xbranch, z) + ewald3_2_sigma_long_Delta_recurrent(&target_view[0], &err_view[0] if get_err else NULL, maxn, x, xbranch, z, False) + if method == 'recurrent_bigimz': + ewald3_2_sigma_long_Delta_recurrent(&target_view[0], &err_view[0] if get_err else NULL, maxn, x, xbranch, z, True) elif method == 'series': ewald3_2_sigma_long_Delta_series(&target_view[0], &err_view[0] if get_err else NULL, maxn, x, xbranch, z) else: diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index 0b5c6be..c0b17c8 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -1036,8 +1036,8 @@ cdef class _ScatteringSystemAtOmegaK: Returns ------- - array_like of complex, with the same shape as `evalpos` - Electric field at the positions given in `evalpos`. + ndarray of complex, with the same shape as `evalpos` + Electric field at the positions given in `evalpos`, in cartesian coordinates. """ if(btyp != QPMS_HANKEL_PLUS): raise NotImplementedError("Only first kind Bessel function-based fields are supported") @@ -1062,6 +1062,51 @@ cdef class _ScatteringSystemAtOmegaK: results[i,1] = res.y results[i,2] = res.z return results.reshape(evalpos.shape) + + def scattered_field_basis(self, evalpos, btyp=QPMS_HANKEL_PLUS): + # TODO examples + """Evaluate scattered field "basis" (periodic system) + + This function enables the evaluation of "scattered" fields + generated by the system for many different excitation + coefficients vectors, without the expensive re-evaluation of the respective + translation operators for each excitation coefficient vector. + + Parameters + ---------- + evalpos: array_like of floats and shape (..., 3) + Evaluation points in cartesian coordinates. + + Returns + ------- + ndarray of complex, with the shape `evalpos.shape[:-1] + (self.fecv_size, 3)` + "Basis" fields at the positions given in `evalpos`, in cartesian coordinates. + """ + if(btyp != QPMS_HANKEL_PLUS): + raise NotImplementedError("Only first kind Bessel function-based fields are supported") + cdef qpms_bessel_t btyp_c = BesselType(btyp) + cdef Py_ssize_t fecv_size = self.fecv_size + evalpos = np.array(evalpos, dtype=float, copy=False) + if evalpos.shape[-1] != 3: + raise ValueError("Last dimension of evalpos has to be 3") + cdef np.ndarray[double,ndim=2] evalpos_a = evalpos.reshape(-1,3) + cdef np.ndarray[complex, ndim=3] results = np.empty((evalpos_a.shape[0], fecv_size, 3), dtype=complex) + cdef ccart3_t *res + res = malloc(fecv_size*sizeof(ccart3_t)) + cdef cart3_t pos + cdef Py_ssize_t i, j + with nogil, wraparound(False), parallel(): + for i in prange(evalpos_a.shape[0]): + pos.x = evalpos_a[i,0] + pos.y = evalpos_a[i,1] + pos.z = evalpos_a[i,2] + qpms_scatsyswk_scattered_field_basis(res, &self.sswk, btyp_c, pos) + for j in range(fecv_size): + results[i,j,0] = res[j].x + results[i,j,1] = res[j].y + results[i,j,2] = res[j].z + free(res) + return results.reshape(evalpos.shape[:-1] + (self.fecv_size, 3)) property fecv_size: def __get__(self): return self.ssw_pyref.ss_pyref.fecv_size diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 8d43538..2a9f8c9 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -705,6 +705,8 @@ cdef extern from "scatsystem.h": const cdouble *f_excitation_vector_full, cart3_t where) nogil ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk, qpms_bessel_t btyp, const cdouble *f_excitation_vector_full, cart3_t where) nogil + qpms_errno_t qpms_scatsyswk_scattered_field_basis(ccart3_t *target, const qpms_scatsys_at_omega_k_t *sswk, + qpms_bessel_t btyp, cart3_t where) nogil double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, cdouble wavenumber, const double *wavevector) nogil diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 42e9b6a..f36c0de 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -2186,6 +2186,86 @@ ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk, return res; } +qpms_errno_t qpms_scatsyswk_scattered_field_basis( + ccart3_t *target, + const qpms_scatsys_at_omega_k_t *sswk, + const qpms_bessel_t btyp, + const cart3_t where + ) { + QPMS_UNTESTED; + if (btyp != QPMS_HANKEL_PLUS) + QPMS_NOT_IMPLEMENTED("Only scattered field with first kind Hankel functions currently implemented."); + const qpms_scatsys_t *ss = sswk->ssw->ss; + if (ss->lattice_dimension != 2) + QPMS_NOT_IMPLEMENTED("Only 2D-periodic lattices implemented"); + //ccart3_t res = {0,0,0}; + //ccart3_t res_kc = {0,0,0}; // kahan sum compensation + + static const int dipspecn = 3; // We have three basis vectors + // bspec containing only electric dipoles + const qpms_vswf_set_spec_t dipspec = { + .n = dipspecn, + .ilist = (qpms_uvswfi_t[]){ + qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, -1, 1), + qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, 0, 1), + qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, +1, 1), + }, + .lMax=1, .lMax_M=0, .lMax_N=1, .lMax_L=-1, + .capacity=0, + .norm = ss->c->normalisation, + }; + + ccart3_t regdipoles_0[dipspecn]; { + const sph_t origin_sph = {.r = 0, .theta = M_PI_2, .phi=0}; // Should work with any theta/phi (TESTWORTHY) + csphvec_t regdipoles_0_sph[dipspecn]; + QPMS_ENSURE_SUCCESS(qpms_uvswf_fill(regdipoles_0_sph, &dipspec, + sph2csph(origin_sph), QPMS_BESSEL_REGULAR)); + for(int i = 0; i < dipspecn; ++i) + regdipoles_0[i] = csphvec2ccart(regdipoles_0_sph[i], origin_sph); + } + + complex double *s; // Translation matrix + QPMS_CRASHING_MALLOC(s, ss->max_bspecn * sizeof(*s) * dipspec.n); + + memset(target, 0, ss->fecv_size * sizeof(*target)); + + for (qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) { + const qpms_vswf_set_spec_t *bspec = qpms_ss_bspec_pi(ss, pi); + const cart3_t particle_pos = ss->p[pi].pos; + //const complex double *particle_cv = cvf + ss->fecv_pstarts[pi]; + + const cart3_t origin_cart = {0, 0, 0}; + + QPMS_ASSERT(sswk->k[2] == 0); // At least not implemented now + QPMS_ASSERT(ss->per.lattice_basis[0].z == 0); + QPMS_ASSERT(ss->per.lattice_basis[1].z == 0); + + // Same choices as in qpms_ss_ppair_W32xy; TODO make it more dynamic + const double maxR = sqrt(ss->per.unitcell_volume) * 64; + const double maxK = 2048 * 2 * M_PI / maxR; + + QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_e32( + ss->c, s, NULL, + &dipspec, 1, bspec, dipspecn, + sswk->eta, sswk->ssw->wavenumber, + cart3xy2cart2(ss->per.lattice_basis[0]), cart3xy2cart2(ss->per.lattice_basis[1]), + cart2_from_double_array(sswk->k), cart3_substract(where, particle_pos) /*CHECKSIGN*/, + maxR, maxK)); + + for(size_t i = 0; i < bspec->n; ++i) + for(size_t j = 0; j < dipspecn; ++j){ + target[ss->fecv_pstarts[pi] + i] = ccart3_add(target[ss->fecv_pstarts[pi] + i], + ccart3_scale(s[dipspecn*i+j], regdipoles_0[j])); + //ccart3_t summand = ccart3_scale(particle_cv[i] * s[dipspecn*i+j], regdipoles_0[j]); + //ckahanadd(&(res.x), &(res_kc.x), summand.x); + //ckahanadd(&(res.y), &(res_kc.y), summand.y); + //ckahanadd(&(res.z), &(res_kc.z), summand.z); + } + } + free(s); + return QPMS_SUCCESS; +} + #if 0 ccart3_t qpms_scatsys_scattered_E_irrep(const qpms_scatsys_t *ss, qpms_iri_t iri, const complex double *cvr, cart3_t where) { diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 7325a6e..bcb3649 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -785,11 +785,24 @@ ccart3_t qpms_scatsysw_scattered_E__alt( */ ccart3_t qpms_scatsyswk_scattered_E( const qpms_scatsys_at_omega_k_t *sswk, - qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS). + qpms_bessel_t typ, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supported). const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$. cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated. ); +/// Evaluates "scattered" field basis functions in a periodic system. +/** + * + * \see qpms_uvswf_fill() evaluates a set of VSWF basis functions (for finite systems). + */ +qpms_errno_t qpms_scatsyswk_scattered_field_basis( + ccart3_t *target, ///< Target array of size sswk->ssw->ss->fecv_size + const qpms_scatsys_at_omega_k_t *sswk, + qpms_bessel_t typ, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supponted). + cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the basis 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]);