Infinite periodic system scattered field "bases"
This commit is contained in:
Marek Nečada
parent
bf297c11c3
commit
868e603f1c
|
|
@ -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:
|
||||
|
|
|
|||
|
|
@ -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")
|
||||
|
|
@ -1063,6 +1063,51 @@ cdef class _ScatteringSystemAtOmegaK:
|
|||
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 = <ccart3_t *> 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
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -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) {
|
||||
|
|
|
|||
|
|
@ -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]);
|
||||
|
|
|
|||
Loading…
Reference in New Issue