translations legacy code removal
Former-commit-id: 98e883f9c7586208623acca05585f747dc57f911
This commit is contained in:
Marek Nečada
parent
0b0ebf2bce
commit
8b9c5a794c
|
|
@ -43,7 +43,7 @@ In simple words, it works like this:
|
||||||
#-------------------------------------
|
#-------------------------------------
|
||||||
|
|
||||||
|
|
||||||
|
# This one is probably not used anymore an can perhaps be removed:
|
||||||
cdef void loop_D_iiiidddii_As_D_lllldddbl(char **args, np.npy_intp *dims, np.npy_intp *steps, void *data) nogil:
|
cdef void loop_D_iiiidddii_As_D_lllldddbl(char **args, np.npy_intp *dims, np.npy_intp *steps, void *data) nogil:
|
||||||
cdef np.npy_intp i, n = dims[0]
|
cdef np.npy_intp i, n = dims[0]
|
||||||
cdef void *func = (<void**>data)#[0]
|
cdef void *func = (<void**>data)#[0]
|
||||||
|
|
@ -99,11 +99,6 @@ np.import_ufunc()
|
||||||
|
|
||||||
# BTW, aren't there anonymous arrays in cython?
|
# BTW, aren't there anonymous arrays in cython?
|
||||||
|
|
||||||
cdef np.PyUFuncGenericFunction trans_X_taylor_loop_func[1]
|
|
||||||
cdef void *trans_A_taylor_elementwise_funcs[1]
|
|
||||||
cdef void *trans_B_taylor_elementwise_funcs[1]
|
|
||||||
|
|
||||||
trans_X_taylor_loop_func[0] = loop_D_iiiidddii_As_D_lllldddbl
|
|
||||||
|
|
||||||
# types to be used for all of the single-type translation
|
# types to be used for all of the single-type translation
|
||||||
# coefficient retrieval ufuncs called like
|
# coefficient retrieval ufuncs called like
|
||||||
|
|
@ -139,38 +134,6 @@ ufunc__get_both_coeff_types[9] = np.NPY_CDOUBLE
|
||||||
ufunc__get_both_coeff_types[10] = np.NPY_CDOUBLE
|
ufunc__get_both_coeff_types[10] = np.NPY_CDOUBLE
|
||||||
|
|
||||||
|
|
||||||
trans_A_taylor_elementwise_funcs[0] = <void*> qpms_trans_single_A_Taylor_ext
|
|
||||||
trans_B_taylor_elementwise_funcs[0] = <void*> qpms_trans_single_B_Taylor_ext
|
|
||||||
|
|
||||||
trans_A_Taylor = np.PyUFunc_FromFuncAndData(
|
|
||||||
trans_X_taylor_loop_func, # func
|
|
||||||
trans_A_taylor_elementwise_funcs, #data
|
|
||||||
ufunc__get_either_trans_coeff_types, # types
|
|
||||||
1, # ntypes: number of supported input types
|
|
||||||
9, # nin: number of input args
|
|
||||||
1, # nout: number of output args
|
|
||||||
0, # identity element, unused
|
|
||||||
"trans_A_Taylor", # name
|
|
||||||
"""
|
|
||||||
TODO computes the E-E or M-M translation coefficient in Taylor's normalisation
|
|
||||||
""", # doc
|
|
||||||
0 # unused, for backward compatibility of numpy c api
|
|
||||||
)
|
|
||||||
|
|
||||||
trans_B_Taylor = np.PyUFunc_FromFuncAndData(
|
|
||||||
trans_X_taylor_loop_func,
|
|
||||||
trans_B_taylor_elementwise_funcs,
|
|
||||||
ufunc__get_either_trans_coeff_types,
|
|
||||||
1, # number of supported input types
|
|
||||||
9, # number of input args
|
|
||||||
1, # number of output args
|
|
||||||
0, # identity element, unused
|
|
||||||
"trans_B_Taylor",
|
|
||||||
"""
|
|
||||||
TODO computes the E-E or M-M translation coefficient in Taylor's normalisation
|
|
||||||
""",
|
|
||||||
0 # unused
|
|
||||||
)
|
|
||||||
|
|
||||||
# ---------------------------------------------
|
# ---------------------------------------------
|
||||||
# Wrapper for the qpms_trans_calculator "class"
|
# Wrapper for the qpms_trans_calculator "class"
|
||||||
|
|
|
||||||
|
|
@ -265,10 +265,6 @@ cdef extern from "symmetries.h":
|
||||||
|
|
||||||
|
|
||||||
cdef extern from "translations.h":
|
cdef extern from "translations.h":
|
||||||
cdouble qpms_trans_single_A_Taylor_ext(int m, int n, int mu, int nu,
|
|
||||||
double r, double th, double ph, int r_ge_d, int J) nogil
|
|
||||||
cdouble qpms_trans_single_B_Taylor_ext(int m, int n, int mu, int nu,
|
|
||||||
double r, double th, double ph, int r_ge_d, int J) nogil
|
|
||||||
struct qpms_trans_calculator:
|
struct qpms_trans_calculator:
|
||||||
int lMax
|
int lMax
|
||||||
size_t nelem
|
size_t nelem
|
||||||
|
|
|
||||||
|
|
@ -114,41 +114,6 @@ static inline int gauntB_Q_max(int M, int n, int mu, int nu) {
|
||||||
return imin(n, imin(nu, (n+nu+1-abs(M+mu))/2));
|
return imin(n, imin(nu, (n+nu+1-abs(M+mu))/2));
|
||||||
}
|
}
|
||||||
|
|
||||||
#if 0 //Apparently, this is duplicit to that in bessel.c (which also support complex x)
|
|
||||||
int qpms_sph_bessel_fill(qpms_bessel_t typ, int lmax, double x, complex double *result_array) {
|
|
||||||
int retval;
|
|
||||||
double tmparr[lmax+1];
|
|
||||||
switch(typ) {
|
|
||||||
case QPMS_BESSEL_REGULAR:
|
|
||||||
retval = gsl_sf_bessel_jl_steed_array(lmax, x, tmparr);
|
|
||||||
for (int l = 0; l <= lmax; ++l) result_array[l] = tmparr[l];
|
|
||||||
return retval;
|
|
||||||
break;
|
|
||||||
case QPMS_BESSEL_SINGULAR: //FIXME: is this precise enough? Would it be better to do it one-by-one?
|
|
||||||
retval = gsl_sf_bessel_yl_array(lmax,x,tmparr);
|
|
||||||
for (int l = 0; l <= lmax; ++l) result_array[l] = tmparr[l];
|
|
||||||
return retval;
|
|
||||||
break;
|
|
||||||
case QPMS_HANKEL_PLUS:
|
|
||||||
case QPMS_HANKEL_MINUS:
|
|
||||||
retval = gsl_sf_bessel_jl_steed_array(lmax, x, tmparr);
|
|
||||||
for (int l = 0; l <= lmax; ++l) result_array[l] = tmparr[l];
|
|
||||||
if(retval) return retval;
|
|
||||||
retval = gsl_sf_bessel_yl_array(lmax, x, tmparr);
|
|
||||||
if (typ==QPMS_HANKEL_PLUS)
|
|
||||||
for (int l = 0; l <= lmax; ++l) result_array[l] += I * tmparr[l];
|
|
||||||
else
|
|
||||||
for (int l = 0; l <= lmax; ++l) result_array[l] +=-I * tmparr[l];
|
|
||||||
return retval;
|
|
||||||
break;
|
|
||||||
default:
|
|
||||||
abort();
|
|
||||||
//return GSL_EDOM;
|
|
||||||
}
|
|
||||||
assert(0);
|
|
||||||
}
|
|
||||||
#endif
|
|
||||||
|
|
||||||
static inline double qpms_trans_normlogfac(qpms_normalisation_t norm,
|
static inline double qpms_trans_normlogfac(qpms_normalisation_t norm,
|
||||||
int m, int n, int mu, int nu) {
|
int m, int n, int mu, int nu) {
|
||||||
return -0.5*(lgamma(n+m+1)-lgamma(n-m+1)+lgamma(nu-mu+1)-lgamma(nu+mu+1));
|
return -0.5*(lgamma(n+m+1)-lgamma(n-m+1)+lgamma(nu-mu+1)-lgamma(nu+mu+1));
|
||||||
|
|
@ -237,109 +202,6 @@ complex double qpms_trans_single_A(qpms_normalisation_t norm,
|
||||||
return presum * sum;
|
return presum * sum;
|
||||||
}
|
}
|
||||||
|
|
||||||
complex double qpms_trans_single_A_Taylor(int m, int n, int mu, int nu, sph_t kdlj,
|
|
||||||
bool r_ge_d, qpms_bessel_t J) {
|
|
||||||
if(r_ge_d) J = QPMS_BESSEL_REGULAR;
|
|
||||||
|
|
||||||
double costheta = cos(kdlj.theta);
|
|
||||||
|
|
||||||
int qmax = gaunt_q_max(-m,n,mu,nu); // nemá tu být +m?
|
|
||||||
// N.B. -m !!!!!!
|
|
||||||
double a1q[qmax+1];
|
|
||||||
int err;
|
|
||||||
gaunt_xu(-m,n,mu,nu,qmax,a1q,&err);
|
|
||||||
double a1q0 = a1q[0];
|
|
||||||
if (err) abort();
|
|
||||||
|
|
||||||
double leg[gsl_sf_legendre_array_n(n+nu)];
|
|
||||||
if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu,costheta,-1,leg)) abort();
|
|
||||||
complex double bes[n+nu+1];
|
|
||||||
if (qpms_sph_bessel_fill(J, n+nu, kdlj.r, bes)) abort();
|
|
||||||
complex double sum = 0;
|
|
||||||
for(int q = 0; q <= qmax; ++q) {
|
|
||||||
int p = n+nu-2*q;
|
|
||||||
int Pp_order = mu-m;
|
|
||||||
//if(p < abs(Pp_order)) continue; // FIXME raději nastav lépe meze
|
|
||||||
assert(p >= abs(Pp_order));
|
|
||||||
double a1q_n = a1q[q] / a1q0;
|
|
||||||
double Pp = leg[gsl_sf_legendre_array_index(p, abs(Pp_order))];
|
|
||||||
if (Pp_order < 0) Pp *= min1pow(mu-m) * exp(lgamma(1+p+Pp_order)-lgamma(1+p-Pp_order));
|
|
||||||
complex double zp = bes[p];
|
|
||||||
complex double summandq = (n*(n+1) + nu*(nu+1) - p*(p+1)) * min1pow(q) * a1q_n * zp * Pp;
|
|
||||||
sum += summandq; // TODO KAHAN
|
|
||||||
}
|
|
||||||
|
|
||||||
double exponent=(lgamma(2*n+1)-lgamma(n+2)+lgamma(2*nu+3)-lgamma(nu+2)
|
|
||||||
+lgamma(n+nu+m-mu+1)-lgamma(n-m+1)-lgamma(nu+mu+1)
|
|
||||||
+lgamma(n+nu+1) - lgamma(2*(n+nu)+1));
|
|
||||||
complex double presum = exp(exponent);
|
|
||||||
presum *= cexp(I*(mu-m)*kdlj.phi) * min1pow(m) * ipow(nu+n) / (4*n);
|
|
||||||
|
|
||||||
// N.B. ipow(nu-n) is different from the general formula!
|
|
||||||
complex double prenormratio = ipow(nu-n) * sqrt(((2.*nu+1)/(2.*n+1))* exp(
|
|
||||||
lgamma(n+m+1)-lgamma(n-m+1)+lgamma(nu-mu+1)-lgamma(nu+mu+1)));
|
|
||||||
return (presum / prenormratio) * sum;
|
|
||||||
}
|
|
||||||
|
|
||||||
// [Xu_old], eq. (83)
|
|
||||||
complex double qpms_trans_single_B_Xu(int m, int n, int mu, int nu, sph_t kdlj,
|
|
||||||
bool r_ge_d, qpms_bessel_t J) {
|
|
||||||
assert(0); // FIXME probably gives wrong values, do not use.
|
|
||||||
if(r_ge_d) J = QPMS_BESSEL_REGULAR;
|
|
||||||
double costheta = cos(kdlj.theta);
|
|
||||||
|
|
||||||
// TODO Qmax cleanup: can I replace Qmax with realQmax???
|
|
||||||
int q2max = gaunt_q_max(-m-1,n+1,mu+1,nu);
|
|
||||||
int Qmax = gaunt_q_max(-m,n+1,mu,nu);
|
|
||||||
int realQmax = gauntB_Q_max(-m, n, mu, nu);
|
|
||||||
double a2q[q2max+1], a3q[Qmax+1], a2q0, a3q0;
|
|
||||||
int err;
|
|
||||||
if (mu == nu) {
|
|
||||||
for (int q = 0; q <= q2max; ++q)
|
|
||||||
a2q[q] = 0;
|
|
||||||
a2q0 = 1;
|
|
||||||
}
|
|
||||||
else {
|
|
||||||
gaunt_xu(-m-1,n+1,mu+1,nu,q2max,a2q,&err); if (err) abort();
|
|
||||||
a2q0 = a2q[0];
|
|
||||||
}
|
|
||||||
gaunt_xu(-m,n+1,mu,nu,Qmax,a3q,&err); if (err) abort();
|
|
||||||
a3q0 = a3q[0];
|
|
||||||
|
|
||||||
double leg[gsl_sf_legendre_array_n(n+nu+1)];
|
|
||||||
if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,costheta,-1,leg)) abort();
|
|
||||||
complex double bes[n+nu+2];
|
|
||||||
if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bes)) abort();
|
|
||||||
|
|
||||||
complex double sum = 0;
|
|
||||||
for (int q = 0; q <= realQmax; ++q) {
|
|
||||||
int p = n+nu-2*q;
|
|
||||||
double a2q_n = a2q[q]/a2q0;
|
|
||||||
double a3q_n = a3q[q]/a3q0;
|
|
||||||
complex double zp_ = bes[p+1];
|
|
||||||
int Pp_order_ = mu-m;
|
|
||||||
//if(p+1 < abs(Pp_order_)) continue; // FIXME raději nastav lépe meze
|
|
||||||
assert(p+1 >= abs(Pp_order_));
|
|
||||||
double Pp_ = leg[gsl_sf_legendre_array_index(p+1, abs(Pp_order_))];
|
|
||||||
if (Pp_order_ < 0) Pp_ *= min1pow(mu-m) * exp(lgamma(1+1+p+Pp_order_)-lgamma(1+1+p-Pp_order_));
|
|
||||||
complex double summandq = ((2*(n+1)*(nu-mu)*a2q_n
|
|
||||||
-(-nu*(nu+1) - n*(n+3) - 2*mu*(n+1)+p*(p+3))* a3q_n)
|
|
||||||
*min1pow(q) * zp_ * Pp_);
|
|
||||||
sum += summandq; // TODO KAHAN
|
|
||||||
}
|
|
||||||
|
|
||||||
double exponent=(lgamma(2*n+3)-lgamma(n+2)+lgamma(2*nu+3)-lgamma(nu+2)
|
|
||||||
+lgamma(n+nu+m-mu+2)-lgamma(n-m+1)-lgamma(nu+mu+1)
|
|
||||||
+lgamma(n+nu+2) - lgamma(2*(n+nu)+3));
|
|
||||||
complex double presum = exp(exponent);
|
|
||||||
presum *= cexp(I*(mu-m)*kdlj.phi) * min1pow(m) * ipow(nu+n+1) / (
|
|
||||||
(4*n)*(n+1)*(n+m+1));
|
|
||||||
|
|
||||||
// Taylor normalisation v2, proven to be equivalent
|
|
||||||
complex double prenormratio = ipow(nu-n);
|
|
||||||
|
|
||||||
return (presum / prenormratio) * sum;
|
|
||||||
}
|
|
||||||
|
|
||||||
complex double qpms_trans_single_B(qpms_normalisation_t norm,
|
complex double qpms_trans_single_B(qpms_normalisation_t norm,
|
||||||
int m, int n, int mu, int nu, sph_t kdlj,
|
int m, int n, int mu, int nu, sph_t kdlj,
|
||||||
|
|
@ -428,93 +290,6 @@ complex double qpms_trans_single_B(qpms_normalisation_t norm,
|
||||||
return presum * sum;
|
return presum * sum;
|
||||||
}
|
}
|
||||||
|
|
||||||
complex double qpms_trans_single_B_Taylor(int m, int n, int mu, int nu, sph_t kdlj,
|
|
||||||
bool r_ge_d, qpms_bessel_t J) {
|
|
||||||
assert(0); // FIXME probably gives wrong values, do not use.
|
|
||||||
if(r_ge_d) J = QPMS_BESSEL_REGULAR;
|
|
||||||
double costheta = cos(kdlj.theta);
|
|
||||||
|
|
||||||
int q2max = gaunt_q_max(-m-1,n+1,mu+1,nu);
|
|
||||||
int Qmax = gaunt_q_max(-m,n+1,mu,nu);
|
|
||||||
int realQmax = gauntB_Q_max(-m,n,mu,nu);
|
|
||||||
double a2q[q2max+1], a3q[Qmax+1], a2q0, a3q0;
|
|
||||||
int err;
|
|
||||||
if (mu == nu) {
|
|
||||||
for (int q = 0; q <= q2max; ++q)
|
|
||||||
a2q[q] = 0;
|
|
||||||
a2q0 = 1;
|
|
||||||
}
|
|
||||||
else {
|
|
||||||
gaunt_xu(-m-1,n+1,mu+1,nu,q2max,a2q,&err); if (err) abort();
|
|
||||||
a2q0 = a2q[0];
|
|
||||||
}
|
|
||||||
gaunt_xu(-m,n+1,mu,nu,Qmax,a3q,&err); if (err) abort();
|
|
||||||
a3q0 = a3q[0];
|
|
||||||
|
|
||||||
double leg[gsl_sf_legendre_array_n(n+nu+1)];
|
|
||||||
if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,costheta,-1,leg)) abort();
|
|
||||||
complex double bes[n+nu+2];
|
|
||||||
if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bes)) abort();
|
|
||||||
|
|
||||||
complex double sum = 0;
|
|
||||||
for (int q = 0; q <= realQmax; ++q) {
|
|
||||||
int p = n+nu-2*q;
|
|
||||||
double a2q_n = a2q[q]/a2q0;
|
|
||||||
double a3q_n = a3q[q]/a3q0;
|
|
||||||
complex double zp_ = bes[p+1];
|
|
||||||
int Pp_order_ = mu-m;
|
|
||||||
//if(p+1 < abs(Pp_order_)) continue; // FIXME raději nastav lépe meze
|
|
||||||
assert(p+1 >= abs(Pp_order_));
|
|
||||||
double Pp_ = leg[gsl_sf_legendre_array_index(p+1, abs(Pp_order_))];
|
|
||||||
if (Pp_order_ < 0) Pp_ *= min1pow(mu-m) * exp(lgamma(1+1+p+Pp_order_)-lgamma(1+1+p-Pp_order_));
|
|
||||||
complex double summandq = ((2*(n+1)*(nu-mu)*a2q_n
|
|
||||||
-(-nu*(nu+1) - n*(n+3) - 2*mu*(n+1)+p*(p+3))* a3q_n)
|
|
||||||
*min1pow(q) * zp_ * Pp_);
|
|
||||||
sum += summandq; //TODO KAHAN
|
|
||||||
}
|
|
||||||
|
|
||||||
double exponent=(lgamma(2*n+3)-lgamma(n+2)+lgamma(2*nu+3)-lgamma(nu+2)
|
|
||||||
+lgamma(n+nu+m-mu+2)-lgamma(n-m+1)-lgamma(nu+mu+1)
|
|
||||||
+lgamma(n+nu+2) - lgamma(2*(n+nu)+3));
|
|
||||||
complex double presum = exp(exponent);
|
|
||||||
presum *= cexp(I*(mu-m)*kdlj.phi) * min1pow(m) * ipow(nu+n+1) / (
|
|
||||||
(4*n)*(n+1)*(n+m+1));
|
|
||||||
|
|
||||||
// Taylor normalisation v2, proven to be equivalent
|
|
||||||
// ipow(nu-n) is different from the new general formula!!!
|
|
||||||
complex double prenormratio = ipow(nu-n) * sqrt(((2.*nu+1)/(2.*n+1))* exp(
|
|
||||||
lgamma(n+m+1)-lgamma(n-m+1)+lgamma(nu-mu+1)-lgamma(nu+mu+1)));
|
|
||||||
|
|
||||||
return (presum / prenormratio) * sum;
|
|
||||||
}
|
|
||||||
|
|
||||||
complex double qpms_trans_single_A_Taylor_ext(int m, int n, int mu, int nu,
|
|
||||||
double kdlj_r, double kdlj_theta, double kdlj_phi, int r_ge_d, int J) {
|
|
||||||
sph_t kdlj = {kdlj_r, kdlj_theta, kdlj_phi};
|
|
||||||
return qpms_trans_single_A_Taylor(m,n,mu,nu,kdlj,r_ge_d,J);
|
|
||||||
}
|
|
||||||
|
|
||||||
complex double qpms_trans_single_B_Taylor_ext(int m, int n, int mu, int nu,
|
|
||||||
double kdlj_r, double kdlj_theta, double kdlj_phi, int r_ge_d, int J) {
|
|
||||||
sph_t kdlj = {kdlj_r, kdlj_theta, kdlj_phi};
|
|
||||||
return qpms_trans_single_B_Taylor(m,n,mu,nu,kdlj,r_ge_d,J);
|
|
||||||
}
|
|
||||||
|
|
||||||
void qpms_trans_calculator_free(qpms_trans_calculator *c) {
|
|
||||||
free(c->A_multipliers[0]);
|
|
||||||
free(c->A_multipliers);
|
|
||||||
free(c->B_multipliers[0]);
|
|
||||||
free(c->B_multipliers);
|
|
||||||
#ifdef LATTICESUMS
|
|
||||||
qpms_ewald3_constants_free(e3c);
|
|
||||||
#endif
|
|
||||||
#ifdef LATTICESUMS_OLD
|
|
||||||
free(c->hct);
|
|
||||||
#endif
|
|
||||||
free(c->legendre0);
|
|
||||||
free(c);
|
|
||||||
}
|
|
||||||
|
|
||||||
static inline size_t qpms_trans_calculator_index_mnmunu(const qpms_trans_calculator *c,
|
static inline size_t qpms_trans_calculator_index_mnmunu(const qpms_trans_calculator *c,
|
||||||
int m, int n, int mu, int nu){
|
int m, int n, int mu, int nu){
|
||||||
return c->nelem * qpms_mn2y(m,n) + qpms_mn2y(mu,nu);
|
return c->nelem * qpms_mn2y(m,n) + qpms_mn2y(mu,nu);
|
||||||
|
|
@ -530,7 +305,7 @@ static inline size_t qpms_trans_calculator_index_yyu(const qpms_trans_calculator
|
||||||
|
|
||||||
static inline double fsq(double x) {return x * x; }
|
static inline double fsq(double x) {return x * x; }
|
||||||
|
|
||||||
static void qpms_trans_calculator_multipliers_A_general(
|
static void qpms_trans_calculator_multipliers_A(
|
||||||
qpms_normalisation_t norm,
|
qpms_normalisation_t norm,
|
||||||
complex double *dest, int m, int n, int mu, int nu, int qmax) {
|
complex double *dest, int m, int n, int mu, int nu, int qmax) {
|
||||||
assert(qmax == gaunt_q_max(-m,n,mu,nu));
|
assert(qmax == gaunt_q_max(-m,n,mu,nu));
|
||||||
|
|
@ -589,7 +364,7 @@ static void qpms_trans_calculator_multipliers_A_general(
|
||||||
|
|
||||||
|
|
||||||
// as in [Xu](61)
|
// as in [Xu](61)
|
||||||
double cruzan_bfactor(int M, int n, int mu, int nu, int p) {
|
static double cruzan_bfactor(int M, int n, int mu, int nu, int p) {
|
||||||
double logprefac = lgamma(n+M+1) - lgamma(n-M+1) + lgamma(nu+mu+1) - lgamma(nu-mu+1)
|
double logprefac = lgamma(n+M+1) - lgamma(n-M+1) + lgamma(nu+mu+1) - lgamma(nu-mu+1)
|
||||||
+ lgamma(p-M-mu+2) - lgamma(p+M+mu+2);
|
+ lgamma(p-M-mu+2) - lgamma(p+M+mu+2);
|
||||||
logprefac *= 0.5;
|
logprefac *= 0.5;
|
||||||
|
|
@ -599,7 +374,7 @@ double cruzan_bfactor(int M, int n, int mu, int nu, int p) {
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
void qpms_trans_calculator_multipliers_B_general(
|
void qpms_trans_calculator_multipliers_B(
|
||||||
qpms_normalisation_t norm,
|
qpms_normalisation_t norm,
|
||||||
complex double *dest, int m, int n, int mu, int nu, int Qmax){
|
complex double *dest, int m, int n, int mu, int nu, int Qmax){
|
||||||
// This is according to the Cruzan-type formula [Xu](59)
|
// This is according to the Cruzan-type formula [Xu](59)
|
||||||
|
|
@ -658,192 +433,6 @@ void qpms_trans_calculator_multipliers_B_general(
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
/*static*/ void qpms_trans_calculator_multipliers_B_general_oldXu(
|
|
||||||
qpms_normalisation_t norm,
|
|
||||||
complex double *dest, int m, int n, int mu, int nu, int Qmax) {
|
|
||||||
assert(0); // FIXME probably gives wrong values, do not use.
|
|
||||||
assert(Qmax == gauntB_Q_max(-m,n,mu,nu));
|
|
||||||
int q2max = gaunt_q_max(-m-1,n+1,mu+1,nu);
|
|
||||||
// assert(Qmax == q2max);
|
|
||||||
// FIXME is it safe to replace q2max with Qmax in gaunt_xu??
|
|
||||||
double a2q[q2max+1], a3q[Qmax+1], a2q0, a3q0;
|
|
||||||
int err;
|
|
||||||
if (mu == nu) {
|
|
||||||
for (int q = 0; q <= q2max; ++q)
|
|
||||||
a2q[q] = 0;
|
|
||||||
a2q0 = 1;
|
|
||||||
}
|
|
||||||
else {
|
|
||||||
gaunt_xu(-m-1,n+1,mu+1,nu,q2max,a2q,&err); if (err) abort();
|
|
||||||
a2q0 = a2q[0];
|
|
||||||
}
|
|
||||||
gaunt_xu(-m,n+1,mu,nu,q2max,a3q,&err); if (err) abort(); // FIXME this should probably go away
|
|
||||||
a3q0 = a3q[0];
|
|
||||||
|
|
||||||
|
|
||||||
int csphase = qpms_normalisation_t_csphase(norm); //TODO FIXME use this
|
|
||||||
double normlogfac= qpms_trans_normlogfac(norm,m,n,mu,nu);
|
|
||||||
double normfac = qpms_trans_normfac(norm,m,n,mu,nu);
|
|
||||||
/// N<-M type coefficients w.r.t. Kristensson's convention. Csphase has been already taken into acct ^^^.
|
|
||||||
normfac *= qpms_normalisation_factor_M_noCS(norm, nu, mu)
|
|
||||||
/ qpms_normalisation_factor_N_noCS(norm, n, m);
|
|
||||||
// TODO use csphase to modify normfac here!!!!
|
|
||||||
// normfac = xxx ? -normfac : normfac;
|
|
||||||
normfac *= min1pow(m);//different from old taylor
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
double exponent=(lgamma(2*n+3)-lgamma(n+2)+lgamma(2*nu+3)-lgamma(nu+2)
|
|
||||||
+lgamma(n+nu+m-mu+2)-lgamma(n-m+1)-lgamma(nu+mu+1)
|
|
||||||
+lgamma(n+nu+2) - lgamma(2*(n+nu)+3))
|
|
||||||
+normlogfac;
|
|
||||||
complex double presum = exp(exponent);
|
|
||||||
presum *= I * ipow(nu+n) /*different from old Taylor */ * normfac / (
|
|
||||||
(4*n)*(n+1)*(n+m+1));
|
|
||||||
|
|
||||||
for (int q = BQ_OFFSET; q <= Qmax; ++q) {
|
|
||||||
int p = n+nu-2*q;
|
|
||||||
double a2q_n = a2q[q]/a2q0;
|
|
||||||
double a3q_n = a3q[q]/a3q0;
|
|
||||||
int Pp_order_ = mu-m;
|
|
||||||
//if(p+1 < abs(Pp_order_)) continue; // FIXME raději nastav lépe meze
|
|
||||||
assert(p+1 >= abs(Pp_order_));
|
|
||||||
double Ppfac = (Pp_order_ >= 0) ? 1 :
|
|
||||||
|
|
||||||
min1pow(mu-m) * exp(lgamma(1+1+p+Pp_order_)-lgamma(1+1+p-Pp_order_));
|
|
||||||
double summandq = ((2*(n+1)*(nu-mu)*a2q_n
|
|
||||||
-(-nu*(nu+1) - n*(n+3) - 2*mu*(n+1)+p*(p+3))* a3q_n)
|
|
||||||
*min1pow(q));
|
|
||||||
dest[q-BQ_OFFSET] = Ppfac * summandq * presum;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
//#if 0
|
|
||||||
static void qpms_trans_calculator_multipliers_A_Taylor(
|
|
||||||
complex double *dest, int m, int n, int mu, int nu, int qmax) {
|
|
||||||
assert(qmax == gaunt_q_max(-m,n,mu,nu));
|
|
||||||
double a1q[qmax+1];
|
|
||||||
int err;
|
|
||||||
gaunt_xu(-m,n,mu,nu,qmax,a1q,&err);
|
|
||||||
if (err) abort();
|
|
||||||
double a1q0 = a1q[0];
|
|
||||||
|
|
||||||
double exponent=(lgamma(2*n+1)-lgamma(n+2)+lgamma(2*nu+3)-lgamma(nu+2)
|
|
||||||
+lgamma(n+nu+m-mu+1)-lgamma(n-m+1)-lgamma(nu+mu+1)
|
|
||||||
+lgamma(n+nu+1) - lgamma(2*(n+nu)+1)) - 0.5*( // ex-prenormratio
|
|
||||||
lgamma(n+m+1)-lgamma(n-m+1)+lgamma(nu-mu+1)-lgamma(nu+mu+1));
|
|
||||||
double presum = exp(exponent);
|
|
||||||
presum *= min1pow(m+n) * sqrt((2.*n+1)/(2.*nu+1)) / (4*n);
|
|
||||||
|
|
||||||
for(int q = 0; q <= qmax; q++) {
|
|
||||||
int p = n+nu-2*q;
|
|
||||||
int Pp_order = mu - m;
|
|
||||||
assert(p >= abs(Pp_order));
|
|
||||||
double a1q_n = a1q[q] / a1q0;
|
|
||||||
// Assuming non_normalized legendre polynomials!
|
|
||||||
double Ppfac = (Pp_order >= 0) ? 1 :
|
|
||||||
min1pow(mu-m) * exp(lgamma(1+p+Pp_order)-lgamma(1+p-Pp_order));
|
|
||||||
double summandfac = (n*(n+1) + nu*(nu+1) - p*(p+1)) * min1pow(q) * a1q_n;
|
|
||||||
dest[q] = presum * summandfac * Ppfac;
|
|
||||||
// FIXME I might not need complex here
|
|
||||||
}
|
|
||||||
}
|
|
||||||
//#endif
|
|
||||||
#if 0
|
|
||||||
static void qpms_trans_calculator_multipliers_A_Taylor(
|
|
||||||
complex double *dest, int m, int n, int mu, int nu, int qmax) {
|
|
||||||
assert(qmax == gaunt_q_max(-m,n,mu,nu));
|
|
||||||
double a1q[qmax+1];
|
|
||||||
int err;
|
|
||||||
gaunt_xu(-m,n,mu,nu,qmax,a1q,&err);
|
|
||||||
if (err) abort();
|
|
||||||
double a1q0 = a1q[0];
|
|
||||||
for(int q = 0; q <= qmax; ++q) {
|
|
||||||
int p = n+nu-2*q;
|
|
||||||
int Pp_order = mu-m;
|
|
||||||
//if(p < abs(Pp_order)) continue; // FIXME raději nastav lépe meze
|
|
||||||
assert(p >= abs(Pp_order));
|
|
||||||
double a1q_n = a1q[q] / a1q0;
|
|
||||||
//double Pp = leg[gsl_sf_legendre_array_index(p, abs(Pp_order))];
|
|
||||||
//complex double zp = bes[p];
|
|
||||||
dest[q] = (n*(n+1) + nu*(nu+1) - p*(p+1)) * min1pow(q) * a1q_n /* * zp * Pp*/;
|
|
||||||
if (Pp_order < 0) dest[q] *= min1pow(mu-m) * exp(lgamma(1+p+Pp_order)-lgamma(1+p-Pp_order));
|
|
||||||
//sum += summandq;
|
|
||||||
}
|
|
||||||
|
|
||||||
double exponent=(lgamma(2*n+1)-lgamma(n+2)+lgamma(2*nu+3)-lgamma(nu+2)
|
|
||||||
+lgamma(n+nu+m-mu+1)-lgamma(n-m+1)-lgamma(nu+mu+1)
|
|
||||||
+lgamma(n+nu+1) - lgamma(2*(n+nu)+1));
|
|
||||||
complex double presum = exp(exponent);
|
|
||||||
presum *=/* cexp(I*(mu-m)*kdlj.phi) * */ min1pow(m) * ipow(nu+n) / (4*n);
|
|
||||||
|
|
||||||
complex double prenormratio = ipow(nu-n) * sqrt(((2.*nu+1)/(2.*n+1))* exp(
|
|
||||||
lgamma(n+m+1)-lgamma(n-m+1)+lgamma(nu-mu+1)-lgamma(nu+mu+1)));
|
|
||||||
//return (presum / prenormratio) * sum;
|
|
||||||
for(int q=0;q<=qmax;++q) dest[q] *= presum / prenormratio;
|
|
||||||
}
|
|
||||||
#endif
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
static void qpms_trans_calculator_multipliers_B_Taylor(
|
|
||||||
complex double *dest, int m, int n, int mu, int nu, int Qmax) {
|
|
||||||
assert(0); // FIXME probably gives wrong values, do not use.
|
|
||||||
assert(Qmax == gauntB_Q_max(-m,n,mu,nu));
|
|
||||||
int q2max = gaunt_q_max(-m-1,n+1,mu+1,nu);
|
|
||||||
//assert(Qmax == q2max);
|
|
||||||
// FIXME remove the q2max variable altogether, as it is probably equal
|
|
||||||
// to Qmax
|
|
||||||
double a2q[q2max+1], a3q[q2max+1], a2q0, a3q0;
|
|
||||||
int err;
|
|
||||||
if (mu == nu) {
|
|
||||||
for (int q = 0; q <= q2max; ++q)
|
|
||||||
a2q[q] = 0;
|
|
||||||
a2q0 = 1;
|
|
||||||
}
|
|
||||||
else {
|
|
||||||
gaunt_xu(-m-1,n+1,mu+1,nu,q2max,a2q,&err); if (err) abort();
|
|
||||||
a2q0 = a2q[0];
|
|
||||||
}
|
|
||||||
gaunt_xu(-m,n+1,mu,nu,q2max,a3q,&err); if (err) abort();
|
|
||||||
a3q0 = a3q[0];
|
|
||||||
|
|
||||||
double exponent=(lgamma(2*n+3)-lgamma(n+2)+lgamma(2*nu+3)-lgamma(nu+2)
|
|
||||||
+lgamma(n+nu+m-mu+2)-lgamma(n-m+1)-lgamma(nu+mu+1)
|
|
||||||
+lgamma(n+nu+2) - lgamma(2*(n+nu)+3)) - 0.5 * (
|
|
||||||
lgamma(n+m+1)-lgamma(n-m+1)+lgamma(nu-mu+1)
|
|
||||||
-lgamma(nu+mu+1));
|
|
||||||
complex double presum = exp(exponent);
|
|
||||||
presum *= I * (min1pow(m+n) *sqrt((2.*n+1)/(2.*nu+1)) / (
|
|
||||||
(4*n)*(n+1)*(n+m+1)));
|
|
||||||
|
|
||||||
for (int q = BQ_OFFSET; q <= Qmax; ++q) {
|
|
||||||
int p = n+nu-2*q;
|
|
||||||
double a2q_n = a2q[q]/a2q0;
|
|
||||||
double a3q_n = a3q[q]/a3q0;
|
|
||||||
int Pp_order_ = mu-m;
|
|
||||||
//if(p+1 < abs(Pp_order_)) continue; // FIXME raději nastav lépe meze
|
|
||||||
assert(p+1 >= abs(Pp_order_));
|
|
||||||
double Ppfac = (Pp_order_ >= 0) ? 1 :
|
|
||||||
|
|
||||||
min1pow(mu-m) * exp(lgamma(1+1+p+Pp_order_)-lgamma(1+1+p-Pp_order_));
|
|
||||||
double summandq = ((2*(n+1)*(nu-mu)*a2q_n
|
|
||||||
-(-nu*(nu+1) - n*(n+3) - 2*mu*(n+1)+p*(p+3))* a3q_n)
|
|
||||||
*min1pow(q));
|
|
||||||
dest[q-BQ_OFFSET] = Ppfac * summandq * presum;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
int qpms_trans_calculator_multipliers_A(qpms_normalisation_t norm, complex double *dest, int m, int n, int mu, int nu, int qmax) {
|
|
||||||
qpms_trans_calculator_multipliers_A_general(norm, dest, m, n, mu, nu, qmax);
|
|
||||||
return 0;
|
|
||||||
}
|
|
||||||
|
|
||||||
int qpms_trans_calculator_multipliers_B(qpms_normalisation_t norm, complex double *dest, int m, int n, int mu, int nu, int Qmax) {
|
|
||||||
qpms_trans_calculator_multipliers_B_general(norm, dest, m, n, mu, nu, Qmax);
|
|
||||||
return 0;
|
|
||||||
}
|
|
||||||
|
|
||||||
qpms_trans_calculator
|
qpms_trans_calculator
|
||||||
*qpms_trans_calculator_init (const int lMax, const qpms_normalisation_t normalisation) {
|
*qpms_trans_calculator_init (const int lMax, const qpms_normalisation_t normalisation) {
|
||||||
TROPS_ONLY_EIMF_IMPLEMENTED(normalisation);
|
TROPS_ONLY_EIMF_IMPLEMENTED(normalisation);
|
||||||
|
|
@ -1556,249 +1145,6 @@ int qpms_trans_calculator_e32_short_points_and_shift(const qpms_trans_calculator
|
||||||
#endif // 0
|
#endif // 0
|
||||||
#endif // LATTICESUMS32
|
#endif // LATTICESUMS32
|
||||||
|
|
||||||
#ifdef LATTICESUMS_OLD
|
|
||||||
|
|
||||||
int qpms_trans_calculator_get_shortrange_AB_arrays_buf(const qpms_trans_calculator *c,
|
|
||||||
complex double *Adest, complex double *Bdest,
|
|
||||||
size_t deststride, size_t srcstride,
|
|
||||||
sph_t kdlj, qpms_bessel_t J,
|
|
||||||
qpms_l_t lrcutoff, unsigned kappa, double cc, // regularisation params
|
|
||||||
complex double *bessel_buf, double *legendre_buf
|
|
||||||
) {
|
|
||||||
assert(J == QPMS_HANKEL_PLUS); // support only J == 3 for now
|
|
||||||
if (0 == kdlj.r && J != QPMS_BESSEL_REGULAR) {
|
|
||||||
for (size_t i = 0; i < c->nelem; ++i)
|
|
||||||
for (size_t j = 0; j < c->nelem; ++j) {
|
|
||||||
*(Adest + i*srcstride + j*deststride) = NAN+I*NAN;
|
|
||||||
*(Bdest + i*srcstride + j*deststride) = NAN+I*NAN;
|
|
||||||
}
|
|
||||||
// TODO warn? different return value?
|
|
||||||
return 0;
|
|
||||||
}
|
|
||||||
{
|
|
||||||
double costheta = cos(kdlj.theta);
|
|
||||||
if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,2*c->lMax+1,
|
|
||||||
costheta,-1,legendre_buf)) abort();
|
|
||||||
// if (qpms_sph_bessel_fill(J, 2*c->lMax+1, kdlj.r, bessel_buf)) abort(); // original
|
|
||||||
hankelparts_fill(NULL, bessel_buf, 2*c->lMax+1, lrcutoff, c->hct, kappa, cc, kdlj.r);
|
|
||||||
size_t desti = 0, srci = 0;
|
|
||||||
for (int n = 1; n <= c->lMax; ++n) for (int m = -n; m <= n; ++m) {
|
|
||||||
for (int nu = 1; nu <= c->lMax; ++nu) for (int mu = -nu; mu <= nu; ++mu) {
|
|
||||||
size_t assertindex = qpms_trans_calculator_index_mnmunu(c,m,n,mu,nu);
|
|
||||||
assert(assertindex == desti*c->nelem + srci);
|
|
||||||
*(Adest + deststride * desti + srcstride * srci) =
|
|
||||||
qpms_trans_calculator_get_A_precalcbuf(c,m,n,mu,nu,
|
|
||||||
kdlj,false,J,bessel_buf,legendre_buf);
|
|
||||||
*(Bdest + deststride * desti + srcstride * srci) =
|
|
||||||
qpms_trans_calculator_get_B_precalcbuf(c,m,n,mu,nu,
|
|
||||||
kdlj,false,J,bessel_buf,legendre_buf);
|
|
||||||
++srci;
|
|
||||||
}
|
|
||||||
++desti;
|
|
||||||
srci = 0;
|
|
||||||
}
|
|
||||||
return 0;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
int qpms_trans_calculator_get_shortrange_AB_buf_p(const qpms_trans_calculator *c,
|
|
||||||
complex double *Adest, complex double *Bdest,
|
|
||||||
int m, int n, int mu, int nu, sph_t kdlj,
|
|
||||||
qpms_bessel_t J,
|
|
||||||
qpms_l_t lrcutoff, unsigned kappa, double cc, // regularisation params
|
|
||||||
complex double *bessel_buf, double *legendre_buf) {
|
|
||||||
assert(J == QPMS_HANKEL_PLUS); // support only J == 3 for now
|
|
||||||
if (0 == kdlj.r && J != QPMS_BESSEL_REGULAR) {
|
|
||||||
*Adest = NAN+I*NAN;
|
|
||||||
*Bdest = NAN+I*NAN;
|
|
||||||
// TODO warn? different return value?
|
|
||||||
return 0;
|
|
||||||
}
|
|
||||||
{
|
|
||||||
double costheta = cos(kdlj.theta);
|
|
||||||
if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,
|
|
||||||
costheta,-1,legendre_buf)) abort();
|
|
||||||
//if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bessel_buf)) abort(); // original
|
|
||||||
hankelparts_fill(NULL, bessel_buf, 2*c->lMax+1, lrcutoff, c->hct, kappa, cc, kdlj.r);
|
|
||||||
|
|
||||||
*Adest = qpms_trans_calculator_get_A_precalcbuf(c,m,n,mu,nu,
|
|
||||||
kdlj,false,J,bessel_buf,legendre_buf);
|
|
||||||
*Bdest = qpms_trans_calculator_get_B_precalcbuf(c,m,n,mu,nu,
|
|
||||||
kdlj,false,J,bessel_buf,legendre_buf);
|
|
||||||
return 0;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// Short-range parts of the translation coefficients
|
|
||||||
int qpms_trans_calculator_get_shortrange_AB_p(const qpms_trans_calculator *c,
|
|
||||||
complex double *Adest, complex double *Bdest,
|
|
||||||
qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj,
|
|
||||||
qpms_bessel_t J /* Only J=3 valid for now */,
|
|
||||||
qpms_l_t lrcutoff, unsigned kappa, double cc) {
|
|
||||||
double leg[gsl_sf_legendre_array_n(2*c->lMax+1)];
|
|
||||||
complex double bes[2*c->lMax+2]; // maximum order is 2n+1 for B coeffs, plus the zeroth.
|
|
||||||
return qpms_trans_calculator_get_shortrange_AB_buf_p(c,Adest, Bdest,m,n,mu,nu,kdlj,J,
|
|
||||||
lrcutoff, kappa, cc,
|
|
||||||
bes, leg);
|
|
||||||
}
|
|
||||||
|
|
||||||
int qpms_trans_calculator_get_shortrange_AB_arrays(const qpms_trans_calculator *c,
|
|
||||||
complex double *Adest, complex double *Bdest,
|
|
||||||
size_t deststride, size_t srcstride,
|
|
||||||
sph_t kdlj, qpms_bessel_t J /* Only J=3 valid for now */,
|
|
||||||
qpms_l_t lrcutoff, unsigned kappa, double cc) {
|
|
||||||
double leg[gsl_sf_legendre_array_n(c->lMax+c->lMax+1)];
|
|
||||||
complex double bes[2*c->lMax+2]; // maximum order is 2n+1 for B coeffs, plus the zeroth.
|
|
||||||
return qpms_trans_calculator_get_shortrange_AB_arrays_buf(c,
|
|
||||||
Adest, Bdest, deststride, srcstride,
|
|
||||||
kdlj, J,
|
|
||||||
lrcutoff, kappa, cc,
|
|
||||||
bes, leg);
|
|
||||||
}
|
|
||||||
|
|
||||||
// Long-range parts
|
|
||||||
static inline complex double qpms_trans_calculator_get_2DFT_longrange_A_precalcbuf(const qpms_trans_calculator *c,
|
|
||||||
int m, int n, int mu, int nu, sph_t k_sph /* theta must be M_PI_2 */,
|
|
||||||
qpms_bessel_t J /* must be 3 for now */,
|
|
||||||
const complex double *lrhankel_recparts_buf) {
|
|
||||||
assert(J == QPMS_HANKEL_PLUS);
|
|
||||||
//assert(k_sph.theta == M_PI_2); CHECK IN ADVANCE INSTEAD
|
|
||||||
//assert(k_sph.r > 0);
|
|
||||||
size_t i = qpms_trans_calculator_index_mnmunu(c, m, n, mu, nu);
|
|
||||||
size_t qmax = c->A_multipliers[i+1] - c->A_multipliers[i] - 1;
|
|
||||||
assert(qmax == gaunt_q_max(-m,n,mu,nu));
|
|
||||||
complex double sum, kahanc;
|
|
||||||
ckahaninit(&sum, &kahanc);
|
|
||||||
for(size_t q = 0; q <= qmax; ++q) {
|
|
||||||
int p = n+nu-2*q;
|
|
||||||
double Pp = c->legendre0[gsl_sf_legendre_array_index(p, abs(mu-m))];
|
|
||||||
complex double zp = trindex_cd(lrhankel_recparts_buf, p)[abs(mu-m)]; // orig: bessel_buf[p];
|
|
||||||
if (mu - m < 0) zp *= min1pow(mu-m); // DLMF 10.4.1
|
|
||||||
complex double multiplier = c->A_multipliers[i][q];
|
|
||||||
ckahanadd(&sum, &kahanc, Pp * zp * multiplier);
|
|
||||||
}
|
|
||||||
complex double eimf = cexp(I*(mu-m)*k_sph.phi);
|
|
||||||
return sum * eimf * ipow(mu-m);
|
|
||||||
}
|
|
||||||
|
|
||||||
static inline complex double qpms_trans_calculator_get_2DFT_longrange_B_precalcbuf(const qpms_trans_calculator *c,
|
|
||||||
int m, int n, int mu, int nu, sph_t k_sph /* theta must be M_PI_2 */,
|
|
||||||
qpms_bessel_t J /* must be 3 for now */,
|
|
||||||
const complex double *lrhankel_recparts_buf) {
|
|
||||||
assert(J == QPMS_HANKEL_PLUS);
|
|
||||||
size_t i = qpms_trans_calculator_index_mnmunu(c, m, n, mu, nu);
|
|
||||||
size_t qmax = c->B_multipliers[i+1] - c->B_multipliers[i] - (1 - BQ_OFFSET);
|
|
||||||
assert(qmax == gauntB_Q_max(-m,n,mu,nu));
|
|
||||||
complex double sum, kahanc;
|
|
||||||
ckahaninit(&sum, &kahanc);
|
|
||||||
for(int q = BQ_OFFSET; q <= qmax; ++q) {
|
|
||||||
int p = n+nu-2*q;
|
|
||||||
double Pp_ = c->legendre0[gsl_sf_legendre_array_index(p+1, abs(mu-m))];
|
|
||||||
complex double zp_ = trindex_cd(lrhankel_recparts_buf, p+1)[abs(mu-m)]; // orig: bessel_buf[p+1];
|
|
||||||
if (mu - m < 0) zp_ *= min1pow(mu-m); // DLMF 10.4.1
|
|
||||||
complex double multiplier = c->B_multipliers[i][q-BQ_OFFSET];
|
|
||||||
ckahanadd(&sum, &kahanc, Pp_ * zp_ * multiplier);
|
|
||||||
}
|
|
||||||
complex double eimf = cexp(I*(mu-m)*k_sph.phi);
|
|
||||||
return sum * eimf * ipow(mu-m);
|
|
||||||
}
|
|
||||||
|
|
||||||
int qpms_trans_calculator_get_2DFT_longrange_AB_buf_p(const qpms_trans_calculator *c,
|
|
||||||
complex double *Adest, complex double *Bdest,
|
|
||||||
int m, int n, int mu, int nu, sph_t k_sph,
|
|
||||||
qpms_bessel_t J,
|
|
||||||
qpms_l_t lrk_cutoff, unsigned kappa, double cv, double k0,
|
|
||||||
complex double *lrhankel_recparts_buf) {
|
|
||||||
|
|
||||||
assert (J == QPMS_HANKEL_PLUS);
|
|
||||||
assert(k_sph.theta == M_PI_2);
|
|
||||||
|
|
||||||
{
|
|
||||||
//double costheta = cos(kdlj.theta);
|
|
||||||
//if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,
|
|
||||||
// costheta,-1,legendre_buf)) abort();
|
|
||||||
//if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bessel_buf)) abort();
|
|
||||||
lrhankel_recpart_fill(lrhankel_recparts_buf, 2*c->lMax+1 /* TODO n+nu+1 might be enough */,
|
|
||||||
lrk_cutoff, c->hct, kappa, cv, k0, k_sph.r);
|
|
||||||
*Adest = qpms_trans_calculator_get_2DFT_longrange_A_precalcbuf(c,m,n,mu,nu,
|
|
||||||
k_sph,J,lrhankel_recparts_buf);
|
|
||||||
*Bdest = qpms_trans_calculator_get_2DFT_longrange_B_precalcbuf(c,m,n,mu,nu,
|
|
||||||
k_sph,J,lrhankel_recparts_buf);
|
|
||||||
return 0;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// Fourier transforms of the long-range parts of the translation coefficients
|
|
||||||
int qpms_trans_calculator_get_2DFT_longrange_AB_p(const qpms_trans_calculator *c,
|
|
||||||
complex double *Adest, complex double *Bdest,
|
|
||||||
qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t k_sph,
|
|
||||||
qpms_bessel_t J /* Only J=3 valid for now */,
|
|
||||||
qpms_l_t lrcutoff, unsigned kappa, double cv, double k0) {
|
|
||||||
int maxp = 2*c->lMax+1; // TODO this may not be needed here, n+nu+1 could be enough instead
|
|
||||||
complex double lrhankel_recpart[maxp * (maxp+1) / 2];
|
|
||||||
return qpms_trans_calculator_get_2DFT_longrange_AB_buf_p(c, Adest, Bdest,m,n,mu,nu,k_sph,
|
|
||||||
J, lrcutoff, kappa, cv, k0, lrhankel_recpart);
|
|
||||||
}
|
|
||||||
|
|
||||||
int qpms_trans_calculator_get_2DFT_longrange_AB_arrays_buf(const qpms_trans_calculator *c,
|
|
||||||
complex double *Adest, complex double *Bdest,
|
|
||||||
size_t deststride, size_t srcstride,
|
|
||||||
sph_t k_sph, qpms_bessel_t J /* must be 3 for now */,
|
|
||||||
qpms_l_t lrk_cutoff, unsigned kappa, double cv, double k0,
|
|
||||||
complex double *lrhankel_recparts_buf) {
|
|
||||||
assert(J == QPMS_HANKEL_PLUS);
|
|
||||||
assert(k_sph.theta == M_PI_2);
|
|
||||||
#if 0
|
|
||||||
if (0 == kdlj.r && J != QPMS_BESSEL_REGULAR) {
|
|
||||||
for (size_t i = 0; i < c->nelem; ++i)
|
|
||||||
for (size_t j = 0; j < c->nelem; ++j) {
|
|
||||||
*(Adest + i*srcstride + j*deststride) = NAN+I*NAN;
|
|
||||||
*(Bdest + i*srcstride + j*deststride) = NAN+I*NAN;
|
|
||||||
}
|
|
||||||
// TODO warn? different return value?
|
|
||||||
return 0;
|
|
||||||
}
|
|
||||||
#endif
|
|
||||||
{
|
|
||||||
lrhankel_recpart_fill(lrhankel_recparts_buf, 2*c->lMax+1,
|
|
||||||
lrk_cutoff, c->hct, kappa, cv, k0, k_sph.r);
|
|
||||||
// if (qpms_sph_bessel_fill(J, 2*c->lMax+1, kdlj.r, bessel_buf)) abort();
|
|
||||||
size_t desti = 0, srci = 0;
|
|
||||||
for (int n = 1; n <= c->lMax; ++n) for (int m = -n; m <= n; ++m) {
|
|
||||||
for (int nu = 1; nu <= c->lMax; ++nu) for (int mu = -nu; mu <= nu; ++mu) {
|
|
||||||
size_t assertindex = qpms_trans_calculator_index_mnmunu(c,m,n,mu,nu);
|
|
||||||
assert(assertindex == desti*c->nelem + srci);
|
|
||||||
*(Adest + deststride * desti + srcstride * srci) =
|
|
||||||
qpms_trans_calculator_get_2DFT_longrange_A_precalcbuf(c,m,n,mu,nu,
|
|
||||||
k_sph,J,lrhankel_recparts_buf);
|
|
||||||
*(Bdest + deststride * desti + srcstride * srci) =
|
|
||||||
qpms_trans_calculator_get_2DFT_longrange_B_precalcbuf(c,m,n,mu,nu,
|
|
||||||
k_sph,J,lrhankel_recparts_buf);
|
|
||||||
++srci;
|
|
||||||
}
|
|
||||||
++desti;
|
|
||||||
srci = 0;
|
|
||||||
}
|
|
||||||
return 0;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// FIXME i get stack smashing error inside the following function if compiled with optimizations
|
|
||||||
int qpms_trans_calculator_get_2DFT_longrange_AB_arrays(const qpms_trans_calculator *c,
|
|
||||||
complex double *Adest, complex double *Bdest,
|
|
||||||
size_t deststride, size_t srcstride,
|
|
||||||
sph_t k_sph, qpms_bessel_t J /* Only J=3 valid for now */,
|
|
||||||
qpms_l_t lrcutoff, unsigned kappa, double cv, double k0) {
|
|
||||||
int maxp = 2*c->lMax+1;
|
|
||||||
complex double lrhankel_recpart[maxp * (maxp+1) / 2];
|
|
||||||
return qpms_trans_calculator_get_2DFT_longrange_AB_arrays_buf(c,
|
|
||||||
Adest, Bdest, deststride, srcstride, k_sph, J,
|
|
||||||
lrcutoff, kappa, cv, k0,
|
|
||||||
lrhankel_recpart);
|
|
||||||
}
|
|
||||||
|
|
||||||
#endif // LATTICESUMS_OLD
|
|
||||||
|
|
||||||
|
|
||||||
complex double qpms_trans_calculator_get_A_ext(const qpms_trans_calculator *c,
|
complex double qpms_trans_calculator_get_A_ext(const qpms_trans_calculator *c,
|
||||||
int m, int n, int mu, int nu,
|
int m, int n, int mu, int nu,
|
||||||
|
|
|
||||||
|
|
@ -40,19 +40,6 @@
|
||||||
#include "ewald.h"
|
#include "ewald.h"
|
||||||
#endif
|
#endif
|
||||||
|
|
||||||
// TODO replace the xplicit "Taylor" functions with general,
|
|
||||||
// taking qpms_normalisation_t argument.
|
|
||||||
complex double qpms_trans_single_A_Taylor(qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj,
|
|
||||||
bool r_ge_d, qpms_bessel_t J);
|
|
||||||
|
|
||||||
complex double qpms_trans_single_B_Taylor(qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj,
|
|
||||||
bool r_ge_d, qpms_bessel_t J);
|
|
||||||
|
|
||||||
complex double qpms_trans_single_A_Taylor_ext(qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, double kdlj_r,
|
|
||||||
double kdlj_th, double kdlj_phi, int r_ge_d, int J);
|
|
||||||
|
|
||||||
complex double qpms_trans_single_B_Taylor_ext(qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, double kdlj_r,
|
|
||||||
double kdlj_th, double kdlj_phi, int r_ge_d, int J);
|
|
||||||
|
|
||||||
complex double qpms_trans_single_A(qpms_normalisation_t norm, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj,
|
complex double qpms_trans_single_A(qpms_normalisation_t norm, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj,
|
||||||
bool r_ge_d, qpms_bessel_t J);
|
bool r_ge_d, qpms_bessel_t J);
|
||||||
|
|
@ -168,34 +155,6 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p(
|
||||||
/// Workspace has to be at least 2 * c->neleme**2 long
|
/// Workspace has to be at least 2 * c->neleme**2 long
|
||||||
);
|
);
|
||||||
|
|
||||||
#ifdef LATTICESUMS_OLD
|
|
||||||
// Short-range parts of the translation coefficients
|
|
||||||
int qpms_trans_calculator_get_shortrange_AB_p(const qpms_trans_calculator *c,
|
|
||||||
complex double *Adest, complex double *Bdest,
|
|
||||||
qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj,
|
|
||||||
qpms_bessel_t J /* Only J=3 valid for now */,
|
|
||||||
qpms_l_t longrange_order_cutoff, unsigned kappa, double cc);
|
|
||||||
int qpms_trans_calculator_get_shortrange_AB_arrays(const qpms_trans_calculator *c,
|
|
||||||
complex double *Adest, complex double *Bdest,
|
|
||||||
size_t deststride, size_t srcstride,
|
|
||||||
sph_t kdlj, qpms_bessel_t J /* Only J=3 valid for now */,
|
|
||||||
qpms_l_t longrange_order_cutoff, unsigned kappa, double cc);
|
|
||||||
|
|
||||||
// Fourier transforms of the long-range parts of the translation coefficients
|
|
||||||
int qpms_trans_calculator_get_2DFT_longrange_AB_p(const qpms_trans_calculator *c,
|
|
||||||
complex double *Adest, complex double *Bdest,
|
|
||||||
qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t k_sph,
|
|
||||||
qpms_bessel_t J /* Only J=3 valid for now */,
|
|
||||||
qpms_l_t longrange_order_cutoff, unsigned kappa, double cv, double k0);
|
|
||||||
|
|
||||||
int qpms_trans_calculator_get_2DFT_longrange_AB_arrays(const qpms_trans_calculator *c,
|
|
||||||
complex double *Adest, complex double *Bdest,
|
|
||||||
size_t deststride, size_t srcstride,
|
|
||||||
sph_t k_sph, qpms_bessel_t J /* Only J=3 valid for now */,
|
|
||||||
qpms_l_t longrange_order_cutoff, unsigned kappa, double cv, double k0);
|
|
||||||
#endif // LATTICESUMS_OLD
|
|
||||||
|
|
||||||
|
|
||||||
#ifdef LATTICESUMS32
|
#ifdef LATTICESUMS32
|
||||||
// for the time being there are only those relatively low-level quick&dirty functions
|
// for the time being there are only those relatively low-level quick&dirty functions
|
||||||
// according to what has been implemented from ewald.h;
|
// according to what has been implemented from ewald.h;
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue