Cython evaluation of "normal" and Ewald-summed SSWF
This commit is contained in:
Marek Nečada
parent
210ddaa527
commit
0078f4db56
|
|
@ -16,7 +16,7 @@ add_library (qpms SHARED translations.c tmatrices.c vecprint.c vswf.c wigner.c e
|
||||||
lattices2d.c gaunt.c error.c legendre.c symmetries.c vecprint.c
|
lattices2d.c gaunt.c error.c legendre.c symmetries.c vecprint.c
|
||||||
bessel.c own_zgemm.c parsing.c scatsystem.c materials.c drudeparam_data.c
|
bessel.c own_zgemm.c parsing.c scatsystem.c materials.c drudeparam_data.c
|
||||||
lll.c beyn.c trivialgroup.c
|
lll.c beyn.c trivialgroup.c
|
||||||
translations_dbg.c
|
translations_dbg.c scatsystem_dbg.c
|
||||||
)
|
)
|
||||||
use_c99()
|
use_c99()
|
||||||
|
|
||||||
|
|
|
||||||
|
|
@ -144,6 +144,46 @@ def get_y_mn_unsigned(int nmax):
|
||||||
i = i + 1
|
i = i + 1
|
||||||
return(ymn_plus, ymn_minus)
|
return(ymn_plus, ymn_minus)
|
||||||
|
|
||||||
|
def get_nelem_scalar(lMax):
|
||||||
|
# TODO DOC, exceptions etc.
|
||||||
|
#return qpms_lMax2nelem(lMax)
|
||||||
|
return lMax * (lMax + 2) + 1 # = (lMax + 1)**2
|
||||||
|
|
||||||
|
## Auxillary function for retrieving the "meshgrid-like" indices; inc. nmax
|
||||||
|
@cython.boundscheck(False)
|
||||||
|
def get_mn_y_scalar(int nmax):
|
||||||
|
"""
|
||||||
|
Auxillary function for retreiving the 'meshgrid-like' indices from the flat indexing;
|
||||||
|
inc. nmax.
|
||||||
|
('y to mn' conversion)
|
||||||
|
|
||||||
|
Parameters
|
||||||
|
----------
|
||||||
|
|
||||||
|
nmax : int
|
||||||
|
The maximum order to which the SSWFs / Legendre functions etc. will be evaluated.
|
||||||
|
|
||||||
|
Returns
|
||||||
|
-------
|
||||||
|
|
||||||
|
output : (m, n)
|
||||||
|
Tuple of two arrays of type np.array(shape=(nmax*nmax + 2*nmax + 1), dtype=np.int),
|
||||||
|
where [(m[y],n[y]) for y in range(nmax*nmax + 2*nmax + 1)] covers all possible
|
||||||
|
integer pairs n >= 0, -n <= m <= n.
|
||||||
|
"""
|
||||||
|
cdef Py_ssize_t nelems = (nmax + 1)**2
|
||||||
|
cdef np.ndarray[np.int_t,ndim=1] m_arr = np.empty([nelems], dtype=np.int)
|
||||||
|
cdef np.ndarray[np.int_t,ndim=1] n_arr = np.empty([nelems], dtype=np.int)
|
||||||
|
cdef Py_ssize_t i = 0
|
||||||
|
cdef np.int_t n, m
|
||||||
|
for n in range(0,nmax+1):
|
||||||
|
for m in range(-n,n+1):
|
||||||
|
m_arr[i] = m
|
||||||
|
n_arr[i] = n
|
||||||
|
i = i + 1
|
||||||
|
return (m_arr, n_arr)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
def tlm2uvswfi(t, l, m):
|
def tlm2uvswfi(t, l, m):
|
||||||
''' TODO doc
|
''' TODO doc
|
||||||
|
|
|
||||||
|
|
@ -1,6 +1,10 @@
|
||||||
from .qpms_cdefs cimport *
|
from .qpms_cdefs cimport *
|
||||||
|
from .qpms_c cimport ScatteringSystem, _ScatteringSystemAtOmega, _ScatteringSystemAtOmegaK
|
||||||
|
from .cycommon import EwaldPart, BesselType
|
||||||
from libc.stdlib cimport malloc, free, calloc
|
from libc.stdlib cimport malloc, free, calloc
|
||||||
import numpy as np
|
import numpy as np
|
||||||
|
from cython.parallel cimport prange, parallel
|
||||||
|
from cython cimport boundscheck, wraparound
|
||||||
|
|
||||||
cdef extern from "ewald.h":
|
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(cdouble *target, double *err, int maxn, cdouble x, int xbranch, cdouble z)
|
||||||
|
|
@ -10,6 +14,24 @@ cdef extern from "ewald.h":
|
||||||
extern int ewald_factor_ipow_l
|
extern int ewald_factor_ipow_l
|
||||||
extern int ewald_factor_ipow_m
|
extern int ewald_factor_ipow_m
|
||||||
|
|
||||||
|
cdef extern from "translations_dbg.h":
|
||||||
|
int qpms_trans_calculator_test_sswf(const qpms_trans_calculator *c,
|
||||||
|
cdouble *dest, const csph_t kdlj, const qpms_bessel_t J) nogil
|
||||||
|
int qpms_trans_calculator_test_sswf_lc3p(const qpms_trans_calculator *c,
|
||||||
|
cdouble *dest, const cdouble k, const cart3_t r, const qpms_bessel_t J) nogil
|
||||||
|
int qpms_trans_calculator_test_sswf_e32(const qpms_trans_calculator *c,
|
||||||
|
cdouble * const sigmas_total, double * serr_total,
|
||||||
|
const double eta, const cdouble k,
|
||||||
|
const cart2_t b1, const cart2_t b2,
|
||||||
|
const cart2_t beta,
|
||||||
|
const cart3_t particle_shift,
|
||||||
|
double maxR, double maxK,
|
||||||
|
const qpms_ewald_part parts) nogil
|
||||||
|
|
||||||
|
cdef extern from "scatsystem_dbg.h":
|
||||||
|
int qpms_scatsyswk_test_sswf_basis_e(cdouble *target, const qpms_scatsys_at_omega_k_t *sswk,
|
||||||
|
qpms_bessel_t typ, cart3_t evalpoint, qpms_ewald_part parts) nogil
|
||||||
|
|
||||||
def e32_Delta(int maxn, cdouble x, cdouble z, int xbranch = 0, get_err=True, method='auto'):
|
def e32_Delta(int maxn, cdouble x, cdouble z, int xbranch = 0, get_err=True, method='auto'):
|
||||||
cdef np.ndarray[double, ndim=1] err_np
|
cdef np.ndarray[double, ndim=1] err_np
|
||||||
cdef double[::1] err_view
|
cdef double[::1] err_view
|
||||||
|
|
@ -31,6 +53,104 @@ def e32_Delta(int maxn, cdouble x, cdouble z, int xbranch = 0, get_err=True, met
|
||||||
else:
|
else:
|
||||||
return target_np
|
return target_np
|
||||||
|
|
||||||
|
def scatsyswk_sswf_basis(_ScatteringSystemAtOmegaK sswk, evalpos, btyp=QPMS_HANKEL_PLUS, ewaldparts=EwaldPart.FULL):
|
||||||
|
"""Evaluate Ewald-summed scalar SWFs that are used to compose the tranlation operators.
|
||||||
|
|
||||||
|
Parameters
|
||||||
|
----------
|
||||||
|
evalpos: array_like of floats and shape (..., 3)
|
||||||
|
Evaluation points in cartesian coordinates.
|
||||||
|
|
||||||
|
Returns
|
||||||
|
-------
|
||||||
|
ndarray of complex, with the shape `evalpos.shape[:-1] + (particle_count, pnelem)`
|
||||||
|
"""
|
||||||
|
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 qpms_ewald_part ewaldparts_c = EwaldPart(ewaldparts)
|
||||||
|
evalpos = np.array(evalpos, dtype=float, copy=False)
|
||||||
|
if evalpos.shape[-1] != 3:
|
||||||
|
raise ValueError("Last dimension of evalpos has to be 3")
|
||||||
|
|
||||||
|
cdef ScatteringSystem ss = sswk.ssw_pyref.ss_pyref
|
||||||
|
cdef qpms_scatsys_t *ss_c = ss.s
|
||||||
|
cdef qpms_scatsys_at_omega_k_t *sswk_c = sswk.rawpointer()
|
||||||
|
cdef const qpms_trans_calculator *c = ss_c[0].c
|
||||||
|
cdef qpms_l_t sswf_lMax = 2*c[0].lMax+1
|
||||||
|
cdef qpms_y_t pnelem = qpms_lMax2nelem_sc(sswf_lMax)
|
||||||
|
cdef qpms_ss_pi_t particle_count = sswk.ssw_pyref.ss_pyref.particle_count
|
||||||
|
|
||||||
|
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], particle_count, pnelem), dtype=complex)
|
||||||
|
cdef cdouble *res
|
||||||
|
cdef cart3_t pos
|
||||||
|
cdef Py_ssize_t i, pi, j
|
||||||
|
with nogil, wraparound(False), parallel():
|
||||||
|
res = <cdouble *> malloc(particle_count * pnelem * sizeof(cdouble))
|
||||||
|
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_test_sswf_basis_e(res, sswk_c, btyp_c, pos, ewaldparts_c)
|
||||||
|
for pi in range(particle_count):
|
||||||
|
for j in range(pnelem):
|
||||||
|
results[i,pi,j] = res[pi * pnelem + j]
|
||||||
|
free(res)
|
||||||
|
return results.reshape(evalpos.shape[:-1] + (particle_count, pnelem))
|
||||||
|
|
||||||
|
def scatsysw_eval_sswf(_ScatteringSystemAtOmega ssw, evalpos, btyp=QPMS_HANKEL_PLUS):
|
||||||
|
"""Evaluate sswfs by internal qpms_trans_calculator at given locations.
|
||||||
|
|
||||||
|
Parameters
|
||||||
|
----------
|
||||||
|
ssw: _ScatteringSystemAtOmega
|
||||||
|
We take the wavenumber and qpms_trans_calculator from this, nothing else.
|
||||||
|
(In a way, this is an overkill, but there is no up-to-date standalone cython
|
||||||
|
wrapper for qpms_trans_calculator.)
|
||||||
|
evalpos: array_like of floats and shape (..., 3)
|
||||||
|
Evaluation points in cartesian coordinates.
|
||||||
|
btyp: BesselType
|
||||||
|
Kind of Bessel functions to be used
|
||||||
|
|
||||||
|
Returns
|
||||||
|
-------
|
||||||
|
ndarray of complex, with the shape `evalpos.shape[:-1] + (pnelem,)`
|
||||||
|
where `pnelem` is TODO
|
||||||
|
"""
|
||||||
|
cdef qpms_bessel_t btyp_c = BesselType(btyp)
|
||||||
|
evalpos = np.array(evalpos, dtype=float, copy=False)
|
||||||
|
if evalpos.shape[-1] != 3:
|
||||||
|
raise ValueError("Last dimension of evalpos has to be 3")
|
||||||
|
|
||||||
|
cdef ScatteringSystem ss = ssw.ss_pyref
|
||||||
|
cdef qpms_scatsys_t *ss_c = ss.s
|
||||||
|
cdef const qpms_trans_calculator *c = ss_c[0].c
|
||||||
|
cdef qpms_l_t sswf_lMax = 2*c[0].lMax+1
|
||||||
|
cdef qpms_y_t pnelem = qpms_lMax2nelem_sc(sswf_lMax)
|
||||||
|
|
||||||
|
cdef cdouble k = ssw.wavenumber
|
||||||
|
|
||||||
|
cdef np.ndarray[double,ndim=2] evalpos_a = evalpos.reshape(-1,3)
|
||||||
|
cdef np.ndarray[complex,ndim=2] results = np.empty((evalpos_a.shape[0], pnelem), dtype=complex)
|
||||||
|
cdef cdouble *res
|
||||||
|
cdef cart3_t pos
|
||||||
|
cdef csph_t kdlj
|
||||||
|
cdef Py_ssize_t i, j
|
||||||
|
with nogil, wraparound(False), parallel():
|
||||||
|
res = <cdouble *> malloc(pnelem * sizeof(cdouble))
|
||||||
|
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]
|
||||||
|
kdlj = cart2csph(pos)
|
||||||
|
kdlj.r *= k
|
||||||
|
qpms_trans_calculator_test_sswf(c, res, kdlj, btyp_c)
|
||||||
|
for j in range(pnelem):
|
||||||
|
results[i,j] = res[j]
|
||||||
|
free(res)
|
||||||
|
return results.reshape(evalpos.shape[:-1] + (pnelem,))
|
||||||
|
|
||||||
def gamma_inc(double a, cdouble x, int xbranch=0):
|
def gamma_inc(double a, cdouble x, int xbranch=0):
|
||||||
cdef qpms_csf_result res
|
cdef qpms_csf_result res
|
||||||
complex_gamma_inc_e(a, x, xbranch, &res)
|
complex_gamma_inc_e(a, x, xbranch, &res)
|
||||||
|
|
|
||||||
|
|
@ -190,10 +190,11 @@ cdef extern from "vswf.h":
|
||||||
csphvec_t qpms_vswf_single_mg_csph(qpms_m_t m, qpms_l_t n, csph_t kdlj, qpms_bessel_t btyp, qpms_normalisation_t norm)
|
csphvec_t qpms_vswf_single_mg_csph(qpms_m_t m, qpms_l_t n, csph_t kdlj, qpms_bessel_t btyp, qpms_normalisation_t norm)
|
||||||
|
|
||||||
cdef extern from "indexing.h":
|
cdef extern from "indexing.h":
|
||||||
qpms_y_t qpms_lMax2nelem(qpms_l_t lMax)
|
qpms_y_t qpms_lMax2nelem(qpms_l_t lMax) nogil
|
||||||
qpms_uvswfi_t qpms_tmn2uvswfi(qpms_vswf_type_t t, qpms_m_t m, qpms_l_t n)
|
qpms_uvswfi_t qpms_tmn2uvswfi(qpms_vswf_type_t t, qpms_m_t m, qpms_l_t n) nogil
|
||||||
qpms_errno_t qpms_uvswfi2tmn(qpms_uvswfi_t u, qpms_vswf_type_t* t, qpms_m_t* m, qpms_l_t* n)
|
qpms_errno_t qpms_uvswfi2tmn(qpms_uvswfi_t u, qpms_vswf_type_t* t, qpms_m_t* m, qpms_l_t* n) nogil
|
||||||
qpms_m_t qpms_uvswfi2m(qpms_uvswfi_t u)
|
qpms_m_t qpms_uvswfi2m(qpms_uvswfi_t u) nogil
|
||||||
|
qpms_y_t qpms_lMax2nelem_sc(qpms_l_t lmax) nogil
|
||||||
# maybe more if needed
|
# maybe more if needed
|
||||||
|
|
||||||
# Point generators from lattices.h
|
# Point generators from lattices.h
|
||||||
|
|
@ -256,95 +257,95 @@ cdef extern from "lattices.h":
|
||||||
|
|
||||||
|
|
||||||
cdef extern from "vectors.h":
|
cdef extern from "vectors.h":
|
||||||
cart2_t cart2_add(const cart2_t a, const cart2_t b)
|
cart2_t cart2_add(const cart2_t a, const cart2_t b) nogil
|
||||||
cart2_t cart2_substract(const cart2_t a, const cart2_t b)
|
cart2_t cart2_substract(const cart2_t a, const cart2_t b) nogil
|
||||||
cart2_t cart2_scale(const double c, const cart2_t v)
|
cart2_t cart2_scale(const double c, const cart2_t v) nogil
|
||||||
double cart2_dot(const cart2_t a, const cart2_t b)
|
double cart2_dot(const cart2_t a, const cart2_t b) nogil
|
||||||
double cart2_normsq(const cart2_t a)
|
double cart2_normsq(const cart2_t a) nogil
|
||||||
double cart2norm(const cart2_t v)
|
double cart2norm(const cart2_t v) nogil
|
||||||
pol_t cart2pol(const cart2_t cart)
|
pol_t cart2pol(const cart2_t cart) nogil
|
||||||
sph_t pol2sph_equator(const pol_t pol)
|
sph_t pol2sph_equator(const pol_t pol) nogil
|
||||||
sph_t cart22sph(const cart2_t cart)
|
sph_t cart22sph(const cart2_t cart) nogil
|
||||||
sph_t cart12sph_zaxis(double z)
|
sph_t cart12sph_zaxis(double z) nogil
|
||||||
cart3_t cart12cart3z(double z)
|
cart3_t cart12cart3z(double z) nogil
|
||||||
cart3_t cart22cart3xy(const cart2_t a)
|
cart3_t cart22cart3xy(const cart2_t a) nogil
|
||||||
cart2_t cart3xy2cart2(const cart3_t a)
|
cart2_t cart3xy2cart2(const cart3_t a) nogil
|
||||||
double cart3_dot(const cart3_t a, const cart3_t b)
|
double cart3_dot(const cart3_t a, const cart3_t b) nogil
|
||||||
cart3_t cart3_vectorprod(const cart3_t a, const cart3_t b)
|
cart3_t cart3_vectorprod(const cart3_t a, const cart3_t b) nogil
|
||||||
double cart3_tripleprod(const cart3_t a, const cart3_t b, const cart3_t c)
|
double cart3_tripleprod(const cart3_t a, const cart3_t b, const cart3_t c) nogil
|
||||||
double cart3_normsq(const cart3_t a)
|
double cart3_normsq(const cart3_t a) nogil
|
||||||
double cart3norm(const cart3_t v)
|
double cart3norm(const cart3_t v) nogil
|
||||||
sph_t cart2sph(const cart3_t cart)
|
sph_t cart2sph(const cart3_t cart) nogil
|
||||||
cart3_t sph2cart(const sph_t sph)
|
cart3_t sph2cart(const sph_t sph) nogil
|
||||||
cart2_t pol2cart(const pol_t pol)
|
cart2_t pol2cart(const pol_t pol) nogil
|
||||||
cart3_t pol2cart3_equator(const pol_t pol)
|
cart3_t pol2cart3_equator(const pol_t pol) nogil
|
||||||
cart3_t cart3_add(const cart3_t a, const cart3_t b)
|
cart3_t cart3_add(const cart3_t a, const cart3_t b) nogil
|
||||||
cart3_t cart3_substract(const cart3_t a, const cart3_t b)
|
cart3_t cart3_substract(const cart3_t a, const cart3_t b) nogil
|
||||||
cart3_t cart3_scale(const double c, const cart3_t v)
|
cart3_t cart3_scale(const double c, const cart3_t v) nogil
|
||||||
cart3_t cart3_divscale( const cart3_t v, const double c)
|
cart3_t cart3_divscale( const cart3_t v, const double c) nogil
|
||||||
double cart3_dist(const cart3_t a, const cart3_t b)
|
double cart3_dist(const cart3_t a, const cart3_t b) nogil
|
||||||
bint cart3_isclose(const cart3_t a, const cart3_t b, double rtol, double atol)
|
bint cart3_isclose(const cart3_t a, const cart3_t b, double rtol, double atol) nogil
|
||||||
ccart3_t ccart3_scale(const cdouble c, const ccart3_t v)
|
ccart3_t ccart3_scale(const cdouble c, const ccart3_t v) nogil
|
||||||
ccart3_t ccart3_add(const ccart3_t a, const ccart3_t b)
|
ccart3_t ccart3_add(const ccart3_t a, const ccart3_t b) nogil
|
||||||
ccart3_t ccart3_substract(const ccart3_t a, const ccart3_t b)
|
ccart3_t ccart3_substract(const ccart3_t a, const ccart3_t b) nogil
|
||||||
cdouble ccart3_dotnc(const ccart3_t a, const ccart3_t b)
|
cdouble ccart3_dotnc(const ccart3_t a, const ccart3_t b) nogil
|
||||||
ccart3_t cart32ccart3(cart3_t c)
|
ccart3_t cart32ccart3(cart3_t c) nogil
|
||||||
csphvec_t csphvec_add(const csphvec_t a, const csphvec_t b)
|
csphvec_t csphvec_add(const csphvec_t a, const csphvec_t b) nogil
|
||||||
csphvec_t csphvec_substract(const csphvec_t a, const csphvec_t b)
|
csphvec_t csphvec_substract(const csphvec_t a, const csphvec_t b) nogil
|
||||||
csphvec_t csphvec_scale(cdouble c, const csphvec_t v)
|
csphvec_t csphvec_scale(cdouble c, const csphvec_t v) nogil
|
||||||
cdouble csphvec_dotnc(const csphvec_t a, const csphvec_t b)
|
cdouble csphvec_dotnc(const csphvec_t a, const csphvec_t b) nogil
|
||||||
sph_t sph_scale(double c, const sph_t s)
|
sph_t sph_scale(double c, const sph_t s) nogil
|
||||||
csph_t sph_cscale(cdouble c, const sph_t s)
|
csph_t sph_cscale(cdouble c, const sph_t s) nogil
|
||||||
ccart3_t csphvec2ccart(const csphvec_t sphvec, const sph_t at)
|
ccart3_t csphvec2ccart(const csphvec_t sphvec, const sph_t at) nogil
|
||||||
ccart3_t csphvec2ccart_csph(const csphvec_t sphvec, const csph_t at)
|
ccart3_t csphvec2ccart_csph(const csphvec_t sphvec, const csph_t at) nogil
|
||||||
csphvec_t ccart2csphvec(const ccart3_t cartvec, const sph_t at)
|
csphvec_t ccart2csphvec(const ccart3_t cartvec, const sph_t at) nogil
|
||||||
csph_t sph2csph(sph_t s)
|
csph_t sph2csph(sph_t s) nogil
|
||||||
sph_t csph2sph(csph_t s)
|
sph_t csph2sph(csph_t s) nogil
|
||||||
csph_t ccart2csph(const ccart3_t cart)
|
csph_t ccart2csph(const ccart3_t cart) nogil
|
||||||
csph_t cart2csph(const cart3_t cart)
|
csph_t cart2csph(const cart3_t cart) nogil
|
||||||
ccart3_t csph2ccart(const csph_t sph)
|
ccart3_t csph2ccart(const csph_t sph) nogil
|
||||||
void ccart3_kahaninit(ccart3_t *sum, ccart3_t *compensation)
|
void ccart3_kahaninit(ccart3_t *sum, ccart3_t *compensation) nogil
|
||||||
void csphvec_kahaninit(csphvec_t *sum, csphvec_t *compensation)
|
void csphvec_kahaninit(csphvec_t *sum, csphvec_t *compensation) nogil
|
||||||
void ccart3_kahanadd(ccart3_t *sum, ccart3_t *compensation, const ccart3_t input)
|
void ccart3_kahanadd(ccart3_t *sum, ccart3_t *compensation, const ccart3_t input) nogil
|
||||||
void csphvec_kahanadd(csphvec_t *sum, csphvec_t *compensation, const csphvec_t input)
|
void csphvec_kahanadd(csphvec_t *sum, csphvec_t *compensation, const csphvec_t input) nogil
|
||||||
double csphvec_norm(const csphvec_t a)
|
double csphvec_norm(const csphvec_t a) nogil
|
||||||
double csphvec_reldiff_abstol(const csphvec_t a, const csphvec_t b, double tolerance)
|
double csphvec_reldiff_abstol(const csphvec_t a, const csphvec_t b, double tolerance) nogil
|
||||||
double csphvec_reldiff(const csphvec_t a, const csphvec_t b)
|
double csphvec_reldiff(const csphvec_t a, const csphvec_t b) nogil
|
||||||
sph_t anycoord2sph(anycoord_point_t p, qpms_coord_system_t t)
|
sph_t anycoord2sph(anycoord_point_t p, qpms_coord_system_t t) nogil
|
||||||
cart3_t anycoord2cart3(anycoord_point_t p, qpms_coord_system_t t)
|
cart3_t anycoord2cart3(anycoord_point_t p, qpms_coord_system_t t) nogil
|
||||||
double anycoord_norm(anycoord_point_t p, qpms_coord_system_t t)
|
double anycoord_norm(anycoord_point_t p, qpms_coord_system_t t) nogil
|
||||||
pol_t anycoord2pol(anycoord_point_t p, qpms_coord_system_t t)
|
pol_t anycoord2pol(anycoord_point_t p, qpms_coord_system_t t) nogil
|
||||||
cart2_t anycoord2cart2(anycoord_point_t p, qpms_coord_system_t t)
|
cart2_t anycoord2cart2(anycoord_point_t p, qpms_coord_system_t t) nogil
|
||||||
double anycoord2cart1(anycoord_point_t p, qpms_coord_system_t t)
|
double anycoord2cart1(anycoord_point_t p, qpms_coord_system_t t) nogil
|
||||||
void qpms_array_coord_transform(void *dest, qpms_coord_system_t tdest,
|
void qpms_array_coord_transform(void *dest, qpms_coord_system_t tdest,
|
||||||
const void *src, qpms_coord_system_t tsrc, size_t nmemb)
|
const void *src, qpms_coord_system_t tsrc, size_t nmemb) nogil
|
||||||
void anycoord_arr2something(void *dest, qpms_coord_system_t tdest,
|
void anycoord_arr2something(void *dest, qpms_coord_system_t tdest,
|
||||||
const void *src, qpms_coord_system_t tsrc, size_t nmemb)
|
const void *src, qpms_coord_system_t tsrc, size_t nmemb) nogil
|
||||||
void cart3_to_double_array(double *a, cart3_t b)
|
void cart3_to_double_array(double *a, cart3_t b) nogil
|
||||||
cart3_t cart3_from_double_array(const double *a)
|
cart3_t cart3_from_double_array(const double *a) nogil
|
||||||
void cart2_to_double_array(double *a, cart2_t b)
|
void cart2_to_double_array(double *a, cart2_t b) nogil
|
||||||
cart2_t cart2_from_double_array(const double *a)
|
cart2_t cart2_from_double_array(const double *a) nogil
|
||||||
const cart2_t CART2_ZERO
|
const cart2_t CART2_ZERO
|
||||||
|
|
||||||
|
|
||||||
cdef extern from "quaternions.h":
|
cdef extern from "quaternions.h":
|
||||||
qpms_quat_t qpms_quat_2c_from_4d(qpms_quat4d_t q)
|
qpms_quat_t qpms_quat_2c_from_4d(qpms_quat4d_t q) nogil
|
||||||
qpms_quat4d_t qpms_quat_4d_from_2c(qpms_quat_t q)
|
qpms_quat4d_t qpms_quat_4d_from_2c(qpms_quat_t q) nogil
|
||||||
qpms_quat_t qpms_quat_mult(qpms_quat_t p, qpms_quat_t q)
|
qpms_quat_t qpms_quat_mult(qpms_quat_t p, qpms_quat_t q) nogil
|
||||||
qpms_quat_t qpms_quat_add(qpms_quat_t p, qpms_quat_t q)
|
qpms_quat_t qpms_quat_add(qpms_quat_t p, qpms_quat_t q) nogil
|
||||||
qpms_quat_t qpms_quat_rscale(double s, qpms_quat_t q)
|
qpms_quat_t qpms_quat_rscale(double s, qpms_quat_t q) nogil
|
||||||
qpms_quat_t qpms_quat_conj(qpms_quat_t q)
|
qpms_quat_t qpms_quat_conj(qpms_quat_t q) nogil
|
||||||
double qpms_quat_norm(qpms_quat_t q)
|
double qpms_quat_norm(qpms_quat_t q) nogil
|
||||||
double qpms_quat_imnorm(qpms_quat_t q)
|
double qpms_quat_imnorm(qpms_quat_t q) nogil
|
||||||
qpms_quat_t qpms_quat_normalise(qpms_quat_t q)
|
qpms_quat_t qpms_quat_normalise(qpms_quat_t q) nogil
|
||||||
qpms_quat_t qpms_quat_log(qpms_quat_t q)
|
qpms_quat_t qpms_quat_log(qpms_quat_t q) nogil
|
||||||
qpms_quat_t qpms_quat_exp(qpms_quat_t q)
|
qpms_quat_t qpms_quat_exp(qpms_quat_t q) nogil
|
||||||
qpms_quat_t qpms_quat_pow(qpms_quat_t q, double exponent)
|
qpms_quat_t qpms_quat_pow(qpms_quat_t q, double exponent) nogil
|
||||||
cdouble qpms_wignerD_elem(qpms_quat_t q, qpms_l_t l,
|
cdouble qpms_wignerD_elem(qpms_quat_t q, qpms_l_t l,
|
||||||
qpms_m_t mp, qpms_m_t m)
|
qpms_m_t mp, qpms_m_t m) nogil
|
||||||
qpms_irot3_t qpms_irot3_mult(qpms_irot3_t p, qpms_irot3_t q)
|
qpms_irot3_t qpms_irot3_mult(qpms_irot3_t p, qpms_irot3_t q) nogil
|
||||||
qpms_irot3_t qpms_irot3_pow(qpms_irot3_t p, int n)
|
qpms_irot3_t qpms_irot3_pow(qpms_irot3_t p, int n) nogil
|
||||||
qpms_quat_t qpms_quat_from_rotvector(cart3_t v)
|
qpms_quat_t qpms_quat_from_rotvector(cart3_t v) nogil
|
||||||
|
|
||||||
cdef extern from "groups.h":
|
cdef extern from "groups.h":
|
||||||
struct qpms_finite_group_irrep_t:
|
struct qpms_finite_group_irrep_t:
|
||||||
|
|
@ -625,6 +626,7 @@ cdef extern from "scatsystem.h":
|
||||||
size_t *saecv_sizes
|
size_t *saecv_sizes
|
||||||
const qpms_finite_group_t *sym
|
const qpms_finite_group_t *sym
|
||||||
qpms_scatsys_periodic_info_t per
|
qpms_scatsys_periodic_info_t per
|
||||||
|
qpms_trans_calculator *c
|
||||||
|
|
||||||
# per[] and other stuff not currently needed in cython
|
# per[] and other stuff not currently needed in cython
|
||||||
void qpms_scatsys_free(qpms_scatsys_t *s)
|
void qpms_scatsys_free(qpms_scatsys_t *s)
|
||||||
|
|
|
||||||
|
|
@ -0,0 +1,48 @@
|
||||||
|
#include "scatsystem_dbg.h"
|
||||||
|
#include "translations_dbg.h"
|
||||||
|
#include "indexing.h"
|
||||||
|
|
||||||
|
|
||||||
|
// analogue of qpms_scatsyswk_scattered_field_basis_e()
|
||||||
|
qpms_errno_t qpms_scatsyswk_test_sswf_basis_e(
|
||||||
|
complex double *target, ///< Target array of size (2 * sswk->ssw->ss->c->lMax + 1) * sswk->ssw->ss->p_count
|
||||||
|
const qpms_scatsys_at_omega_k_t *sswk,
|
||||||
|
qpms_bessel_t btyp, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supponted).
|
||||||
|
cart3_t where, ///< A point \f$ \vect r \f$, at which the basis is evaluated.
|
||||||
|
qpms_ewald_part parts
|
||||||
|
) {
|
||||||
|
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
|
||||||
|
|
||||||
|
const size_t particle_nelem = qpms_lMax2nelem_sc(2*ss->c->lMax + 1);
|
||||||
|
|
||||||
|
|
||||||
|
for (qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) {
|
||||||
|
const cart3_t particle_pos = ss->p[pi].pos;
|
||||||
|
|
||||||
|
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_test_sswf_e32(ss->c,
|
||||||
|
target + pi * particle_nelem, NULL,
|
||||||
|
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, parts));
|
||||||
|
}
|
||||||
|
|
||||||
|
return QPMS_SUCCESS;
|
||||||
|
}
|
||||||
|
|
@ -0,0 +1,10 @@
|
||||||
|
#include "scatsystem.h"
|
||||||
|
|
||||||
|
qpms_errno_t qpms_scatsyswk_test_sswf_basis_e(
|
||||||
|
complex double *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.
|
||||||
|
qpms_ewald_part parts
|
||||||
|
);
|
||||||
|
|
||||||
Loading…
Reference in New Issue