2 #include "qpms_types.h"
3 #include "qpms_specfunc.h"
5 #include "translations.h"
6 #include "indexing.h" // TODO replace size_t and int with own index types here
8 #include <gsl/gsl_sf_legendre.h>
9 #include <gsl/gsl_sf_bessel.h>
10 #include "tiny_inlines.h"
11 #include "assert_cython_workaround.h"
13 #include <gsl/gsl_sf_coupling.h>
14 #include "qpms_error.h"
15 #include "normalisation.h"
16 #include "translations_inlines.h"
19 * Define macros with additional factors that "should not be there" according
20 * to the "original" formulae but are needed to work with my vswfs.
21 * (actually, I don't know whether the error is in using "wrong" implementation
22 * of vswfs, "wrong" implementation of Xu's translation coefficient formulae,
23 * error/inconsintency in Xu's paper or something else)
24 * Anyway, the zeroes give the correct _numerical_ values according to Xu's
25 * paper tables (without Xu's typos, of course), while
26 * the predefined macros give the correct translations of the VSWFs for the
27 * QPMS_NORMALIZATION_TAYLOR_CS norm.
29 #if !(defined AN0 || defined AN1 || defined AN2 || defined AN3)
30 #pragma message "using AN1 macro as default"
33 //#if !(defined AM0 || defined AM2)
36 #if !(defined BN0 || defined BN1 || defined BN2 || defined BN3)
37 #pragma message "using BN1 macro as default"
40 //#if !(defined BM0 || defined BM2)
43 //#if !(defined BF0 || defined BF1 || defined BF2 || defined BF3)
47 // if defined, the pointer B_multipliers[y] corresponds to the q = 1 element;
48 // otherwise, it corresponds to the q = 0 element, which should be identically zero
49 #ifdef QPMS_PACKED_B_MULTIPLIERS
56 // Translation operators for real sph. harm. based waves are not yet implemented...
57 static inline void TROPS_ONLY_EIMF_IMPLEMENTED(qpms_normalisation_t norm
) {
58 if (norm
& (QPMS_NORMALISATION_SPHARM_REAL
| QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE
))
59 QPMS_NOT_IMPLEMENTED("Translation operators for real or inverse complex spherical harmonics based waves are not implemented.");
62 // Use if only the symmetric form [[A, B], [B, A]] (without additional factors) of translation operator is allowed.
63 static inline void TROPS_ONLY_AB_SYMMETRIC_NORMS_IMPLEMENTED(qpms_normalisation_t norm
) {
64 switch (norm
& QPMS_NORMALISATION_NORM_BITS
) {
65 case QPMS_NORMALISATION_NORM_SPHARM
:
66 case QPMS_NORMALISATION_NORM_POWER
:
69 QPMS_NOT_IMPLEMENTED("Only spherical harmonic and power normalisation supported.");
72 ( !(norm
& QPMS_NORMALISATION_N_I
) != !(norm
& QPMS_NORMALISATION_M_I
) )
74 ( !(norm
& QPMS_NORMALISATION_N_MINUS
) != !(norm
& QPMS_NORMALISATION_M_MINUS
) )
76 QPMS_NOT_IMPLEMENTED("Only normalisations without a phase factors between M and N waves are supported.");
82 * [Xu_old] Yu-Lin Xu, Journal of Computational Physics 127, 285–298 (1996)
83 * [Xu] Yu-Lin Xu, Journal of Computational Physics 139, 137–165 (1998)
87 * GENERAL TODO: use normalised Legendre functions for Kristensson and Taylor conventions directly
88 * instead of normalising them here (the same applies for csphase).
91 static const double sqrtpi
= 1.7724538509055160272981674833411451827975494561223871;
92 //static const double ln2 = 0.693147180559945309417232121458176568075500134360255254120;
94 // Associated Legendre polynomial at zero argument (DLMF 14.5.1)
95 double qpms_legendre0(int m
, int n
) {
96 return pow(2,m
) * sqrtpi
/ tgamma(.5*n
- .5*m
+ .5) / tgamma(.5*n
-.5*m
);
99 // Derivative of associated Legendre polynomial at zero argument (DLMF 14.5.2)
100 double qpms_legendreD0(int m
, int n
) {
101 return -2 * qpms_legendre0(m
, n
);
105 static inline int imin(int x
, int y
) {
106 return x
> y
? y
: x
;
109 // The uppermost value of q index for the B coefficient terms from [Xu](60).
110 // N.B. this is different from [Xu_old](79) due to the n vs. n+1 difference.
111 // However, the trailing terms in [Xu_old] are analytically zero (although
112 // the numerical values will carry some non-zero rounding error).
113 static inline int gauntB_Q_max(int M
, int n
, int mu
, int nu
) {
114 return imin(n
, imin(nu
, (n
+nu
+1-abs(M
+mu
))/2));
117 static inline double qpms_trans_normlogfac(qpms_normalisation_t norm
,
118 int m
, int n
, int mu
, int nu
) {
119 return -0.5*(lgamma(n
+m
+1)-lgamma(n
-m
+1)+lgamma(nu
-mu
+1)-lgamma(nu
+mu
+1));
122 static inline double qpms_trans_normfac(qpms_normalisation_t norm
,
123 int m
, int n
, int mu
, int nu
) {
124 int csphase
= qpms_normalisation_t_csphase(norm
);
125 /* Account for csphase here. Alternatively, this could be done by
126 * using appropriate csphase in the legendre polynomials when calculating
127 * the translation operator.
129 double normfac
= (1 == csphase
) ? min1pow(m
-mu
) : 1.;
130 normfac
*= sqrt((n
*(n
+1.))/(nu
*(nu
+1.)));
131 normfac
*= sqrt((2.*n
+1)/(2.*nu
+1));
135 complex double qpms_trans_single_A(qpms_normalisation_t norm
,
136 int m
, int n
, int mu
, int nu
, csph_t kdlj
,
137 bool r_ge_d
, qpms_bessel_t J
) {
138 TROPS_ONLY_EIMF_IMPLEMENTED(norm
);
139 if(r_ge_d
) J
= QPMS_BESSEL_REGULAR
;
141 double costheta
= cos(kdlj
.theta
);
143 int qmax
= gaunt_q_max(-m
,n
,mu
,nu
); // nemá tu být +m?
147 gaunt_xu(-m
,n
,mu
,nu
,qmax
,a1q
,&err
);
148 QPMS_ENSURE_SUCCESS(err
);
149 double a1q0
= a1q
[0];
151 double leg
[gsl_sf_legendre_array_n(n
+nu
)];
152 QPMS_ENSURE_SUCCESS(gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE
,n
+nu
,costheta
,-1,leg
));
153 complex double bes
[n
+nu
+1];
154 QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(J
, n
+nu
, kdlj
.r
, bes
));
155 complex double sum
= 0;
156 for(int q
= 0; q
<= qmax
; ++q
) {
159 //if(p < abs(Pp_order)) continue; // FIXME raději nastav lépe meze
160 assert(p
>= abs(Pp_order
));
161 double a1q_n
= a1q
[q
] / a1q0
;
162 double Pp
= leg
[gsl_sf_legendre_array_index(p
, abs(Pp_order
))];
163 if (Pp_order
< 0) Pp
*= min1pow(mu
-m
) * exp(lgamma(1+p
+Pp_order
)-lgamma(1+p
-Pp_order
));
164 complex double zp
= bes
[p
];
165 complex double summandq
= (n
*(n
+1) + nu
*(nu
+1) - p
*(p
+1)) * min1pow(q
) * a1q_n
* zp
* Pp
;
166 sum
+= summandq
; // TODO KAHAN
169 double exponent
=(lgamma(2*n
+1)-lgamma(n
+2)+lgamma(2*nu
+3)-lgamma(nu
+2)
170 +lgamma(n
+nu
+m
-mu
+1)-lgamma(n
-m
+1)-lgamma(nu
+mu
+1)
171 +lgamma(n
+nu
+1) - lgamma(2*(n
+nu
)+1));
172 complex double presum
= exp(exponent
);
173 presum
*= cexp(I
*(mu
-m
)*kdlj
.phi
) * min1pow(m
) * ipow(nu
+n
) / (4*n
);
175 double normlogfac
= qpms_trans_normlogfac(norm
,m
,n
,mu
,nu
);
176 double normfac
= qpms_trans_normfac(norm
,m
,n
,mu
,nu
);
177 /// N<-N type coefficients w.r.t. Kristensson's convention. Csphase has been already taken into acct ^^^.
178 normfac
*= qpms_normalisation_factor_N_noCS(norm
, nu
, mu
)
179 / qpms_normalisation_factor_N_noCS(norm
, n
, m
);
180 // ipow(n-nu) is the difference from the Taylor formula!
181 presum
*= /*ipow(n-nu) * */
182 (normfac
* exp(normlogfac
))
204 complex double qpms_trans_single_B(qpms_normalisation_t norm
,
205 int m
, int n
, int mu
, int nu
, csph_t kdlj
,
206 bool r_ge_d
, qpms_bessel_t J
) {
207 TROPS_ONLY_EIMF_IMPLEMENTED(norm
);
208 if(r_ge_d
) J
= QPMS_BESSEL_REGULAR
;
209 double costheta
= cos(kdlj
.theta
);
211 int q2max
= gaunt_q_max(-m
-1,n
+1,mu
+1,nu
);
212 int Qmax
= gaunt_q_max(-m
,n
+1,mu
,nu
);
213 int realQmax
= gauntB_Q_max(-m
,n
,mu
,nu
);
214 double a2q
[q2max
+1], a3q
[Qmax
+1], a2q0
, a3q0
;
217 for (int q
= 0; q
<= q2max
; ++q
)
222 gaunt_xu(-m
-1,n
+1,mu
+1,nu
,q2max
,a2q
,&err
);
223 QPMS_ENSURE_SUCCESS(err
);
226 gaunt_xu(-m
,n
+1,mu
,nu
,Qmax
,a3q
,&err
);
227 QPMS_ENSURE_SUCCESS(err
);
230 double leg
[gsl_sf_legendre_array_n(n
+nu
+1)];
231 QPMS_ENSURE_SUCCESS(gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE
,n
+nu
+1,costheta
,-1,leg
));
232 complex double bes
[n
+nu
+2];
233 QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(J
, n
+nu
+1, kdlj
.r
, bes
));
235 complex double sum
= 0;
236 for (int q
= 0; q
<= realQmax
; ++q
) {
238 double a2q_n
= a2q
[q
]/a2q0
;
239 double a3q_n
= a3q
[q
]/a3q0
;
240 complex double zp_
= bes
[p
+1];
241 int Pp_order_
= mu
-m
;
242 //if(p+1 < abs(Pp_order_)) continue; // FIXME raději nastav lépe meze
243 assert(p
+1 >= abs(Pp_order_
));
244 double Pp_
= leg
[gsl_sf_legendre_array_index(p
+1, abs(Pp_order_
))];
245 if (Pp_order_
< 0) Pp_
*= min1pow(mu
-m
) * exp(lgamma(1+1+p
+Pp_order_
)-lgamma(1+1+p
-Pp_order_
));
246 complex double summandq
= ((2*(n
+1)*(nu
-mu
)*a2q_n
247 -(-nu
*(nu
+1) - n
*(n
+3) - 2*mu
*(n
+1)+p
*(p
+3))* a3q_n
)
248 *min1pow(q
) * zp_
* Pp_
);
249 sum
+= summandq
; //TODO KAHAN
252 double exponent
=(lgamma(2*n
+3)-lgamma(n
+2)+lgamma(2*nu
+3)-lgamma(nu
+2)
253 +lgamma(n
+nu
+m
-mu
+2)-lgamma(n
-m
+1)-lgamma(nu
+mu
+1)
254 +lgamma(n
+nu
+2) - lgamma(2*(n
+nu
)+3));
255 complex double presum
= exp(exponent
);
256 presum
*= cexp(I
*(mu
-m
)*kdlj
.phi
) * min1pow(m
) * ipow(nu
+n
+1) / (
257 (4*n
)*(n
+1)*(n
+m
+1));
259 double normlogfac
= qpms_trans_normlogfac(norm
,m
,n
,mu
,nu
);
260 double normfac
= qpms_trans_normfac(norm
,m
,n
,mu
,nu
);
261 /// N<-M type coefficients w.r.t. Kristensson's convention. Csphase has been already taken into acct ^^^.
262 normfac
*= qpms_normalisation_factor_M_noCS(norm
, nu
, mu
)
263 / qpms_normalisation_factor_N_noCS(norm
, n
, m
);
265 // ipow(n-nu) is the difference from the "old Taylor" formula
266 presum
*= /*ipow(n-nu) * */(exp(normlogfac
) * normfac
)
288 void qpms_trans_calculator_free(qpms_trans_calculator
*c
) {
289 free(c
->A_multipliers
[0]);
290 free(c
->A_multipliers
);
291 free(c
->B_multipliers
[0]);
292 free(c
->B_multipliers
);
294 qpms_ewald3_constants_free(c
->e3c
);
300 static inline size_t qpms_trans_calculator_index_mnmunu(const qpms_trans_calculator
*c
,
301 int m
, int n
, int mu
, int nu
){
302 return c
->nelem
* qpms_mn2y(m
,n
) + qpms_mn2y(mu
,nu
);
305 static inline size_t qpms_trans_calculator_index_yyu(const qpms_trans_calculator
*c
,
306 size_t y
, size_t yu
) {
307 return c
->nelem
* y
+ yu
;
311 #define SQ(x) ((x)*(x))
313 static inline double fsq(double x
) {return x
* x
; }
315 static void qpms_trans_calculator_multipliers_A(
316 qpms_normalisation_t norm
,
317 complex double *dest
, int m
, int n
, int mu
, int nu
, int qmax
) {
318 assert(qmax
== gaunt_q_max(-m
,n
,mu
,nu
));
321 gaunt_xu(-m
,n
,mu
,nu
,qmax
,a1q
,&err
);
322 QPMS_ENSURE_SUCCESS(err
);
323 double a1q0
= a1q
[0];
325 double normlogfac
= qpms_trans_normlogfac(norm
,m
,n
,mu
,nu
);
326 double normfac
= qpms_trans_normfac(norm
,m
,n
,mu
,nu
);
327 /// N<-N type coefficients w.r.t. Kristensson's convention. Csphase has been already taken into acct ^^^.
328 normfac
*= qpms_normalisation_factor_N_noCS(norm
, nu
, mu
)
329 / qpms_normalisation_factor_N_noCS(norm
, n
, m
);
331 normfac
*= min1pow(m
); //different from old Taylor
333 double exponent
=(lgamma(2*n
+1)-lgamma(n
+2)+lgamma(2*nu
+3)-lgamma(nu
+2)
334 +lgamma(n
+nu
+m
-mu
+1)-lgamma(n
-m
+1)-lgamma(nu
+mu
+1)
335 +lgamma(n
+nu
+1) - lgamma(2*(n
+nu
)+1))
337 complex double presum
= exp(exponent
);
338 presum
*= normfac
/ (4.*n
);
339 presum
*= ipow(n
+nu
); // different from old Taylor
341 for(int q
= 0; q
<= qmax
; q
++) {
343 int Pp_order
= mu
- m
;
344 assert(p
>= abs(Pp_order
));
345 double a1q_n
= a1q
[q
] / a1q0
;
346 // Assuming non_normalized legendre polynomials (normalisation done here by hand)!
347 double Ppfac
= (Pp_order
>= 0) ? 1 :
348 min1pow(mu
-m
) * exp(lgamma(1+p
+Pp_order
)-lgamma(1+p
-Pp_order
));
349 double summandfac
= (n
*(n
+1) + nu
*(nu
+1) - p
*(p
+1)) * min1pow(q
) * a1q_n
;
350 dest
[q
] = presum
* summandfac
* Ppfac
368 // FIXME I might not need complex here
374 static double cruzan_bfactor(int M
, int n
, int mu
, int nu
, int p
) {
375 double logprefac
= lgamma(n
+M
+1) - lgamma(n
-M
+1) + lgamma(nu
+mu
+1) - lgamma(nu
-mu
+1)
376 + lgamma(p
-M
-mu
+2) - lgamma(p
+M
+mu
+2);
378 return min1pow(mu
+M
) * (2*p
+3) * exp(logprefac
)
379 * gsl_sf_coupling_3j(2*n
, 2*nu
, 2*(p
+1), 2*M
, 2*mu
, 2*(-M
-mu
))
380 * gsl_sf_coupling_3j(2*n
, 2*nu
, 2*p
, 0, 0, 0);
384 void qpms_trans_calculator_multipliers_B(
385 qpms_normalisation_t norm
,
386 complex double *dest
, int m
, int n
, int mu
, int nu
, int Qmax
){
387 // This is according to the Cruzan-type formula [Xu](59)
388 assert(Qmax
== gauntB_Q_max(-m
,n
,mu
,nu
));
390 double normlogfac
= qpms_trans_normlogfac(norm
,m
,n
,mu
,nu
);
391 double normfac
= qpms_trans_normfac(norm
,m
,n
,mu
,nu
);
392 /// N<-M type coefficients w.r.t. Kristensson's convention. Csphase has been already taken into acct ^^^.
393 normfac
*= qpms_normalisation_factor_M_noCS(norm
, nu
, mu
)
394 / qpms_normalisation_factor_N_noCS(norm
, n
, m
);
396 double presum
= min1pow(1-m
) * (2*nu
+1)/(2.*(n
*(n
+1)))
397 * exp(lgamma(n
+m
+1) - lgamma(n
-m
+1) + lgamma(nu
-mu
+1) - lgamma(nu
+mu
+1)
401 for(int q
= BQ_OFFSET
; q
<= Qmax
; ++q
) {
403 int Pp_order
= mu
- m
;
404 // Assuming non-normalised Legendre polynomials, normalise here by hand.
405 // Ppfac_ differs from Ppfac in the A-case by the substitution p->p+1
406 double Ppfac_
= (Pp_order
>= 0)? 1 :
407 min1pow(mu
-m
) * exp(lgamma(1+1+p
+Pp_order
)-lgamma(1+1+p
-Pp_order
));
410 * (isq(n
+nu
+1)-isq(p
+1))
412 dest
[q
-BQ_OFFSET
] = presum
* t
* Ppfac_
413 * cruzan_bfactor(-m
,n
,mu
,nu
,p
) * ipow(p
+1)
441 qpms_trans_calculator
442 *qpms_trans_calculator_init (const int lMax
, const qpms_normalisation_t normalisation
) {
443 TROPS_ONLY_EIMF_IMPLEMENTED(normalisation
);
445 qpms_trans_calculator
*c
= malloc(sizeof(qpms_trans_calculator
));
447 c
->nelem
= lMax
* (lMax
+2);
448 c
->A_multipliers
= malloc((1+SQ(c
->nelem
)) * sizeof(complex double *));
449 c
->B_multipliers
= malloc((1+SQ(c
->nelem
)) * sizeof(complex double *));
450 c
->normalisation
= normalisation
;
451 size_t *qmaxes
= malloc(SQ(c
->nelem
) * sizeof(size_t));
453 for(size_t y
= 0; y
< c
->nelem
; y
++)
454 for(size_t yu
= 0; yu
< c
->nelem
; yu
++) {
456 qpms_y2mn_p(y
,&m
,&n
);
457 qpms_y2mn_p(yu
,&mu
,&nu
);
459 qmaxes
[qpms_trans_calculator_index_yyu(c
,y
,yu
)]
460 = gaunt_q_max(-m
,n
,mu
,nu
));
462 c
->A_multipliers
[0] = malloc(qmaxsum
* sizeof(complex double));
463 // calculate multiplier beginnings
464 for(size_t i
= 0; i
< SQ(c
->nelem
); ++i
)
465 c
->A_multipliers
[i
+1] = c
->A_multipliers
[i
] + qmaxes
[i
] + 1;
466 // calculate the multipliers
467 for(size_t y
= 0; y
< c
->nelem
; ++y
)
468 for(size_t yu
= 0; yu
< c
->nelem
; ++yu
) {
469 size_t i
= y
* c
->nelem
+ yu
;
471 qpms_y2mn_p(y
, &m
, &n
);
472 qpms_y2mn_p(yu
, &mu
, &nu
);
473 qpms_trans_calculator_multipliers_A(normalisation
,
474 c
->A_multipliers
[i
], m
, n
, mu
, nu
, qmaxes
[i
]);
478 for(size_t y
=0; y
< c
->nelem
; y
++)
479 for(size_t yu
= 0; yu
< c
->nelem
; yu
++) {
481 qpms_y2mn_p(y
,&m
,&n
);
482 qpms_y2mn_p(yu
,&mu
,&nu
);
483 qmaxsum
+= (1 - BQ_OFFSET
) + (
484 qmaxes
[qpms_trans_calculator_index_yyu(c
,y
,yu
)]
485 = gauntB_Q_max(-m
,n
,mu
,nu
));
487 c
->B_multipliers
[0] = malloc(qmaxsum
* sizeof(complex double));
488 // calculate multiplier beginnings
489 for(size_t i
= 0; i
< SQ(c
->nelem
); ++i
)
490 c
->B_multipliers
[i
+1] = c
->B_multipliers
[i
] + qmaxes
[i
] + (1 - BQ_OFFSET
);
491 // calculate the multipliers
492 for(size_t y
= 0; y
< c
->nelem
; ++y
)
493 for(size_t yu
= 0; yu
< c
->nelem
; ++yu
) {
494 size_t i
= y
* c
->nelem
+ yu
;
496 qpms_y2mn_p(y
, &m
, &n
);
497 qpms_y2mn_p(yu
, &mu
, &nu
);
498 qpms_trans_calculator_multipliers_B(normalisation
,
499 c
->B_multipliers
[i
], m
, n
, mu
, nu
, qmaxes
[i
]);
504 c
->e3c
= qpms_ewald3_constants_init(2 * lMax
+ 1, /*csphase*/ qpms_normalisation_t_csphase(normalisation
));
506 c
->legendre0
= malloc(gsl_sf_legendre_array_n(2*lMax
+1) * sizeof(double));
507 QPMS_ENSURE_SUCCESS(gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE
,2*lMax
+1,
508 0,-1,c
->legendre0
)); // TODO maybe use some "precise" analytical formula instead?
512 static inline complex double qpms_trans_calculator_get_A_precalcbuf(const qpms_trans_calculator
*c
,
513 int m
, int n
, int mu
, int nu
, double kdlj_phi
,
514 const complex double *bessel_buf
, const double *legendre_buf
) {
515 TROPS_ONLY_EIMF_IMPLEMENTED(c
->normalisation
);
516 size_t i
= qpms_trans_calculator_index_mnmunu(c
, m
, n
, mu
, nu
);
517 size_t qmax
= c
->A_multipliers
[i
+1] - c
->A_multipliers
[i
] - 1;
518 assert(qmax
== gaunt_q_max(-m
,n
,mu
,nu
));
519 complex double sum
, kahanc
;
520 ckahaninit(&sum
, &kahanc
);
521 for(size_t q
= 0; q
<= qmax
; ++q
) {
523 double Pp
= legendre_buf
[gsl_sf_legendre_array_index(p
, abs(mu
-m
))];
524 complex double zp
= bessel_buf
[p
];
525 complex double multiplier
= c
->A_multipliers
[i
][q
];
526 ckahanadd(&sum
, &kahanc
, Pp
* zp
* multiplier
);
528 complex double eimf
= cexp(I
*(mu
-m
)*kdlj_phi
);
532 complex double qpms_trans_calculator_get_A_buf(const qpms_trans_calculator
*c
,
533 int m
, int n
, int mu
, int nu
, csph_t kdlj
,
534 bool r_ge_d
, qpms_bessel_t J
,
535 complex double *bessel_buf
, double *legendre_buf
) {
536 // This functions gets preallocated memory for bessel and legendre functions, but computes them itself
537 if (r_ge_d
) J
= QPMS_BESSEL_REGULAR
;
538 if (0 == kdlj
.r
&& J
!= QPMS_BESSEL_REGULAR
)
541 int csphase
= qpms_normalisation_t_csphase(c
->normalisation
);
543 double costheta
= cos(kdlj
.theta
);
544 QPMS_ENSURE_SUCCESS(gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE
,n
+nu
,
545 costheta
,csphase
,legendre_buf
));
546 QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(J
, n
+nu
+1, kdlj
.r
, bessel_buf
));
547 return qpms_trans_calculator_get_A_precalcbuf(c
,m
,n
,mu
,nu
,
548 kdlj
.phi
,bessel_buf
,legendre_buf
);
551 static inline complex double qpms_trans_calculator_get_B_precalcbuf(const qpms_trans_calculator
*c
,
552 int m
, int n
, int mu
, int nu
, double kdlj_phi
,
553 const complex double *bessel_buf
, const double *legendre_buf
) {
554 TROPS_ONLY_EIMF_IMPLEMENTED(c
->normalisation
);
555 size_t i
= qpms_trans_calculator_index_mnmunu(c
, m
, n
, mu
, nu
);
556 size_t qmax
= c
->B_multipliers
[i
+1] - c
->B_multipliers
[i
] - (1 - BQ_OFFSET
);
557 assert(qmax
== gauntB_Q_max(-m
,n
,mu
,nu
));
558 complex double sum
, kahanc
;
559 ckahaninit(&sum
, &kahanc
);
560 for(int q
= BQ_OFFSET
; q
<= qmax
; ++q
) {
562 double Pp_
= legendre_buf
[gsl_sf_legendre_array_index(p
+1, abs(mu
-m
))];
563 complex double zp_
= bessel_buf
[p
+1];
564 complex double multiplier
= c
->B_multipliers
[i
][q
-BQ_OFFSET
];
565 ckahanadd(&sum
, &kahanc
, Pp_
* zp_
* multiplier
);
567 complex double eimf
= cexp(I
*(mu
-m
)*kdlj_phi
);
571 complex double qpms_trans_calculator_get_B_buf(const qpms_trans_calculator
*c
,
572 int m
, int n
, int mu
, int nu
, csph_t kdlj
,
573 bool r_ge_d
, qpms_bessel_t J
,
574 complex double *bessel_buf
, double *legendre_buf
) {
575 // This functions gets preallocated memory for bessel and legendre functions, but computes them itself
576 if (r_ge_d
) J
= QPMS_BESSEL_REGULAR
;
577 if (0 == kdlj
.r
&& J
!= QPMS_BESSEL_REGULAR
)
580 int csphase
= qpms_normalisation_t_csphase(c
->normalisation
);
581 double costheta
= cos(kdlj
.theta
);
582 QPMS_ENSURE_SUCCESS(gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE
,n
+nu
+1,
583 costheta
,csphase
,legendre_buf
));
584 QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(J
, n
+nu
+1, kdlj
.r
, bessel_buf
));
585 return qpms_trans_calculator_get_B_precalcbuf(c
,m
,n
,mu
,nu
,
586 kdlj
.phi
,bessel_buf
,legendre_buf
);
589 int qpms_trans_calculator_get_AB_buf_p(const qpms_trans_calculator
*c
,
590 complex double *Adest
, complex double *Bdest
,
591 int m
, int n
, int mu
, int nu
, csph_t kdlj
,
592 bool r_ge_d
, qpms_bessel_t J
,
593 complex double *bessel_buf
, double *legendre_buf
) {
594 if (r_ge_d
) J
= QPMS_BESSEL_REGULAR
;
595 if (0 == kdlj
.r
&& J
!= QPMS_BESSEL_REGULAR
) {
598 // TODO warn? different return value?
601 double costheta
= cos(kdlj
.theta
);
602 QPMS_ENSURE_SUCCESS(gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE
,n
+nu
+1,
603 costheta
,-1,legendre_buf
));
604 QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(J
, n
+nu
+1, kdlj
.r
, bessel_buf
));
605 *Adest
= qpms_trans_calculator_get_A_precalcbuf(c
,m
,n
,mu
,nu
,
606 kdlj
.phi
,bessel_buf
,legendre_buf
);
607 *Bdest
= qpms_trans_calculator_get_B_precalcbuf(c
,m
,n
,mu
,nu
,
608 kdlj
.phi
,bessel_buf
,legendre_buf
);
612 int qpms_trans_calculator_get_AB_arrays_precalcbuf(const qpms_trans_calculator
*c
,
613 complex double *Adest
, complex double *Bdest
,
614 size_t deststride
, size_t srcstride
, double kdlj_phi
,
615 const complex double *bessel_buf
, const double *legendre_buf
) {
616 size_t desti
= 0, srci
= 0;
617 for (int n
= 1; n
<= c
->lMax
; ++n
) for (int m
= -n
; m
<= n
; ++m
) {
618 for (int nu
= 1; nu
<= c
->lMax
; ++nu
) for (int mu
= -nu
; mu
<= nu
; ++mu
) {
620 size_t assertindex
= qpms_trans_calculator_index_mnmunu(c
,m
,n
,mu
,nu
);
622 assert(assertindex
== desti
*c
->nelem
+ srci
);
623 *(Adest
+ deststride
* desti
+ srcstride
* srci
) =
624 qpms_trans_calculator_get_A_precalcbuf(c
,m
,n
,mu
,nu
,
625 kdlj_phi
, bessel_buf
, legendre_buf
);
626 *(Bdest
+ deststride
* desti
+ srcstride
* srci
) =
627 qpms_trans_calculator_get_B_precalcbuf(c
,m
,n
,mu
,nu
,
628 kdlj_phi
,bessel_buf
,legendre_buf
);
637 int qpms_trans_calculator_get_AB_arrays_buf(const qpms_trans_calculator
*c
,
638 complex double *Adest
, complex double *Bdest
,
639 size_t deststride
, size_t srcstride
,
640 csph_t kdlj
, bool r_ge_d
, qpms_bessel_t J
,
641 complex double *bessel_buf
, double *legendre_buf
) {
642 if (r_ge_d
) J
= QPMS_BESSEL_REGULAR
;
643 if (0 == kdlj
.r
&& J
!= QPMS_BESSEL_REGULAR
) {
644 for (size_t i
= 0; i
< c
->nelem
; ++i
)
645 for (size_t j
= 0; j
< c
->nelem
; ++j
) {
646 *(Adest
+ i
*srcstride
+ j
*deststride
) = NAN
+I
*NAN
;
647 *(Bdest
+ i
*srcstride
+ j
*deststride
) = NAN
+I
*NAN
;
649 // TODO warn? different return value?
653 double costheta
= cos(kdlj
.theta
);
654 QPMS_ENSURE_SUCCESS(gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE
,2*c
->lMax
+1,
655 costheta
,-1,legendre_buf
));
656 QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(J
, 2*c
->lMax
+1, kdlj
.r
, bessel_buf
));
658 return qpms_trans_calculator_get_AB_arrays_precalcbuf(c
, Adest
, Bdest
,
659 deststride
, srcstride
, kdlj
.phi
, bessel_buf
, legendre_buf
);
662 complex double qpms_trans_calculator_get_A(const qpms_trans_calculator
*c
,
663 int m
, int n
, int mu
, int nu
, csph_t kdlj
,
664 bool r_ge_d
, qpms_bessel_t J
) {
665 double leg
[gsl_sf_legendre_array_n(n
+nu
)];
666 complex double bes
[n
+nu
+1]; // maximum order is 2n for A coeffs, plus the zeroth.
667 return qpms_trans_calculator_get_A_buf(c
,m
,n
,mu
,nu
,kdlj
,r_ge_d
,J
,
671 complex double qpms_trans_calculator_get_B(const qpms_trans_calculator
*c
,
672 int m
, int n
, int mu
, int nu
, csph_t kdlj
,
673 bool r_ge_d
, qpms_bessel_t J
) {
674 double leg
[gsl_sf_legendre_array_n(n
+nu
+1)];
675 complex double bes
[n
+nu
+2]; // maximum order is 2n+1 for B coeffs, plus the zeroth.
676 return qpms_trans_calculator_get_B_buf(c
,m
,n
,mu
,nu
,kdlj
,r_ge_d
,J
,
680 int qpms_trans_calculator_get_AB_p(const qpms_trans_calculator
*c
,
681 complex double *Adest
, complex double *Bdest
,
682 int m
, int n
, int mu
, int nu
, csph_t kdlj
,
683 bool r_ge_d
, qpms_bessel_t J
) {
684 double leg
[gsl_sf_legendre_array_n(2*c
->lMax
+1)];
685 complex double bes
[2*c
->lMax
+2]; // maximum order is 2n+1 for B coeffs, plus the zeroth.
686 return qpms_trans_calculator_get_AB_buf_p(c
,Adest
, Bdest
,m
,n
,mu
,nu
,kdlj
,r_ge_d
,J
,
690 int qpms_trans_calculator_get_AB_arrays(const qpms_trans_calculator
*c
,
691 complex double *Adest
, complex double *Bdest
,
692 size_t deststride
, size_t srcstride
,
693 csph_t kdlj
, bool r_ge_d
, qpms_bessel_t J
) {
694 double leg
[gsl_sf_legendre_array_n(c
->lMax
+c
->lMax
+1)];
695 complex double bes
[2*c
->lMax
+2]; // maximum order is 2n+1 for B coeffs, plus the zeroth.
696 return qpms_trans_calculator_get_AB_arrays_buf(c
,
697 Adest
, Bdest
, deststride
, srcstride
,
702 // Convenience functions using VSWF base specs
703 qpms_errno_t
qpms_trans_calculator_get_trans_array(const qpms_trans_calculator
*c
,
704 complex double *target
,
705 /// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm.
706 const qpms_vswf_set_spec_t
*destspec
, size_t deststride
,
707 /// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
708 const qpms_vswf_set_spec_t
*srcspec
, size_t srcstride
,
709 csph_t kdlj
, bool r_ge_d
, qpms_bessel_t J
)
711 TROPS_ONLY_AB_SYMMETRIC_NORMS_IMPLEMENTED(c
->normalisation
);
712 assert(c
->normalisation
== destspec
->norm
&& c
->normalisation
== srcspec
->norm
);
713 assert(c
->lMax
>= destspec
->lMax
&& c
->lMax
>= srcspec
->lMax
);
714 assert(destspec
->lMax_L
< 0 && srcspec
->lMax_L
< 0);
715 // TODO don't use c->lMax etc. if both destspec->lMax and srcspec->lMax are smaller
716 complex double A
[c
->nelem
][c
->nelem
];
717 complex double B
[c
->nelem
][c
->nelem
];
718 qpms_errno_t retval
= qpms_trans_calculator_get_AB_arrays(c
,
719 A
[0], B
[0], c
->nelem
, 1,
721 qpms_trans_array_from_AB(target
, destspec
, deststride
, srcspec
, srcstride
,
722 A
[0], B
[0], c
->lMax
);
726 qpms_errno_t
qpms_trans_calculator_get_trans_array_e32_e(const qpms_trans_calculator
*c
,
727 complex double *target
, double *err
,
728 /// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm.
729 const qpms_vswf_set_spec_t
*destspec
, size_t deststride
,
730 /// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
731 const qpms_vswf_set_spec_t
*srcspec
, size_t srcstride
,
732 const double eta
, const complex double k
,
733 cart2_t b1
, cart2_t b2
,
735 const cart2_t particle_shift
,
736 double maxR
, double maxK
,
737 const qpms_ewald_part parts
740 TROPS_ONLY_AB_SYMMETRIC_NORMS_IMPLEMENTED(c
->normalisation
);
741 QPMS_ENSURE(c
->normalisation
== destspec
->norm
&& c
->normalisation
== srcspec
->norm
,
742 "The normalisation conventions must be the same");
743 assert(c
->lMax
>= destspec
->lMax
&& c
->lMax
>= srcspec
->lMax
);
744 assert(destspec
->lMax_L
< 0 && srcspec
->lMax_L
< 0);
745 // TODO don't use c->lMax etc. if both destspec->lMax and srcspec->lMax are smaller
746 const ptrdiff_t ldAB
= c
->nelem
;
747 complex double *A
, *B
;
748 double *Aerr
= NULL
, *Berr
= NULL
;
749 QPMS_CRASHING_MALLOC(A
, c
->nelem
*c
->nelem
*sizeof(complex double));
750 QPMS_CRASHING_MALLOC(B
, c
->nelem
*c
->nelem
*sizeof(complex double));
752 QPMS_CRASHING_MALLOC(Aerr
, c
->nelem
*c
->nelem
*sizeof(double));
753 QPMS_CRASHING_MALLOC(Berr
, c
->nelem
*c
->nelem
*sizeof(double));
755 qpms_errno_t retval
= qpms_trans_calculator_get_AB_arrays_e32_e(c
,
756 A
, Aerr
, B
, Berr
, ldAB
, 1,
757 eta
, k
, b1
, b2
, beta
, particle_shift
, maxR
, maxK
, parts
);
758 for (size_t desti
= 0; desti
< destspec
->n
; ++desti
) {
759 // TODO replace with (modified) qpms_trans_array_from_AB()
760 qpms_y_t desty
; qpms_vswf_type_t destt
;
761 if(QPMS_SUCCESS
!= qpms_uvswfi2ty(destspec
->ilist
[desti
], &destt
, &desty
))
762 qpms_pr_error_at_flf(__FILE__
,__LINE__
,__func__
,
763 "Invalid u. vswf index %llx.", destspec
->ilist
[desti
]);
764 for (size_t srci
= 0; srci
< srcspec
->n
; ++srci
){
765 qpms_y_t srcy
; qpms_vswf_type_t srct
;
766 if(QPMS_SUCCESS
!= qpms_uvswfi2ty(srcspec
->ilist
[srci
], &srct
, &srcy
))
767 qpms_pr_error_at_flf(__FILE__
,__LINE__
,__func__
,
768 "Invalid u. vswf index %llx.", srcspec
->ilist
[srci
]);
769 target
[srci
* srcstride
+ desti
* deststride
]
770 = (srct
== destt
) ? A
[ldAB
*desty
+ srcy
] : B
[ldAB
*desty
+ srcy
];
771 if(err
) err
[srci
* srcstride
+ desti
* deststride
]
772 = (srct
== destt
) ? Aerr
[ldAB
*desty
+ srcy
] : Berr
[ldAB
*desty
+ srcy
];
776 if (err
) { free(Aerr
); free(Berr
); }
780 qpms_errno_t
qpms_trans_calculator_get_trans_array_e32(const qpms_trans_calculator
*c
,
781 complex double *target
, double *err
,
782 /// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm.
783 const qpms_vswf_set_spec_t
*destspec
, size_t deststride
,
784 /// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
785 const qpms_vswf_set_spec_t
*srcspec
, size_t srcstride
,
786 const double eta
, const complex double k
,
787 cart2_t b1
, cart2_t b2
,
789 const cart2_t particle_shift
,
790 double maxR
, double maxK
793 return qpms_trans_calculator_get_trans_array_e32_e(c
, target
, err
, destspec
, deststride
,
794 srcspec
, srcstride
, eta
, k
, b1
, b2
, beta
, particle_shift
, maxR
, maxK
, QPMS_EWALD_FULL
);
798 qpms_errno_t
qpms_trans_calculator_get_trans_array_lc3p(
799 const qpms_trans_calculator
*c
,
800 complex double *target
,
801 /// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm.
802 const qpms_vswf_set_spec_t
*destspec
, size_t deststride
,
803 /// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
804 const qpms_vswf_set_spec_t
*srcspec
, size_t srcstride
,
805 complex double k
, cart3_t destpos
, cart3_t srcpos
, qpms_bessel_t J
806 /// Workspace has to be at least 2 * c->neleme**2 long
809 csph_t kdlj
= cart2csph(cart3_substract(destpos
, srcpos
));
811 return qpms_trans_calculator_get_trans_array(c
, target
,
812 destspec
, deststride
, srcspec
, srcstride
, kdlj
,
817 int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_trans_calculator
*c
,
818 complex double * const Adest
, double * const Aerr
,
819 complex double * const Bdest
, double * const Berr
,
820 const ptrdiff_t deststride
, const ptrdiff_t srcstride
,
821 /* qpms_bessel_t J*/ // assume QPMS_HANKEL_PLUS
822 const double eta
, const double k
, const double unitcell_area
,
823 const size_t nRpoints
, const cart2_t
*Rpoints
, // n.b. can't contain 0; TODO automatic recognition and skip
824 const size_t nKpoints
, const cart2_t
*Kpoints
,
825 const double beta
,//DIFF21
826 const double particle_shift
//DIFF21
830 const qpms_y_t nelem2_sc
= qpms_lMax2nelem_sc(c
->e3c
->lMax
);
831 //const qpms_y_t nelem = qpms_lMax2nelem(c->lMax);
832 const bool doerr
= Aerr
|| Berr
;
833 const bool do_sigma0
= (particle_shift
== 0)//DIFF21((particle_shift.x == 0) && (particle_shift.y == 0)); // FIXME ignoring the case where particle_shift equals to lattice vector
835 complex double *sigmas_short
= malloc(sizeof(complex double)*nelem2_sc
);
836 complex double *sigmas_long
= malloc(sizeof(complex double)*nelem2_sc
);
837 complex double *sigmas_total
= malloc(sizeof(complex double)*nelem2_sc
);
838 double *serr_short
, *serr_long
, *serr_total
;
840 serr_short
= malloc(sizeof(double)*nelem2_sc
);
841 serr_long
= malloc(sizeof(double)*nelem2_sc
);
842 serr_total
= malloc(sizeof(double)*nelem2_sc
);
843 } else serr_short
= serr_long
= serr_total
= NULL
;
845 QPMS_ENSURE_SUCCESS(ewald31z_sigma_long_points_and_shift(sigmas_long
, serr_long
, //DIFF21
846 c
->e3c
, eta
, k
, unitcell_area
, nKpoints
, Kpoints
, beta
, particle_shift
));
848 QPMS_ENSURE_SUCCESS(ewald31z_sigma_short_points_and_shift(sigmas_short
, serr_short
, //DIFF21
849 c
->e3c
, eta
, k
, nRpoints
, Rpoints
, beta
, particle_shift
));
851 for(qpms_y_t y
= 0; y
< nelem2_sc
; ++y
)
852 sigmas_total
[y
] = sigmas_short
[y
] + sigmas_long
[y
];
853 if (doerr
) for(qpms_y_t y
= 0; y
< nelem2_sc
; ++y
)
854 serr_total
[y
] = serr_short
[y
] + serr_long
[y
];
856 complex double sigma0
= 0; double sigma0_err
= 0;
858 QPMS_ENSURE_SUCCESS(ewald31z_sigma0(&sigma0
, &sigma0_err
, c
->e3c
, eta
, k
));
859 const qpms_l_t y
= qpms_mn2y_sc(0,0);
860 sigmas_total
[y
] += sigma0
;
861 if(doerr
) serr_total
[y
] += sigma0_err
;
865 ptrdiff_t desti
= 0, srci
= 0;
866 for (qpms_l_t n
= 1; n
<= c
->lMax
; ++n
) for (qpms_m_t m
= -n
; m
<= n
; ++m
) {
867 for (qpms_l_t nu
= 1; nu
<= c
->lMax
; ++nu
) for (qpms_m_t mu
= -nu
; mu
<= nu
; ++mu
){
868 const size_t i
= qpms_trans_calculator_index_mnmunu(c
, m
, n
, mu
, nu
);
869 const size_t qmax
= c
->A_multipliers
[i
+1] - c
->A_multipliers
[i
] - 1;
870 complex double Asum
, Asumc
; ckahaninit(&Asum
, &Asumc
);
871 double Asumerr
, Asumerrc
; if(Aerr
) kahaninit(&Asumerr
, &Asumerrc
);
873 const qpms_m_t mu_m
= mu
- m
;
874 // TODO skip if ... (N.B. skip will be different for 31z and 32)
875 for(qpms_l_t q
= 0; q
<= qmax
; ++q
) {
876 const qpms_l_t p
= n
+ nu
- 2*q
;
877 const qpms_y_t y_sc
= qpms_mn2y_sc(mu_m
, p
);
878 const complex double multiplier
= c
->A_multipliers
[i
][q
];
879 complex double sigma
= sigmas_total
[y_sc
];
880 ckahanadd(&Asum
, &Asumc
, multiplier
* sigma
);
881 if (Aerr
) kahanadd(&Asumerr
, &Asumerrc
, multiplier
* serr_total
[y_sc
]);
884 *(Adest
+ deststride
* desti
+ srcstride
* srci
) = Asum
;
885 if (Aerr
) *(Aerr
+ deststride
* desti
+ srcstride
* srci
) = Asumerr
;
888 complex double Bsum
, Bsumc
; ckahaninit(&Bsum
, &Bsumc
);
889 double Bsumerr
, Bsumerrc
; if(Berr
) kahaninit(&Bsumerr
, &Bsumerrc
);
890 for(qpms_l_t q
= 0; q
<= qmax
; ++q
) {
891 const qpms_l_t p_
= n
+ nu
- 2*q
+ 1;
892 const qpms_y_t y_sc
= qpms_mn2y_sc(mu_m
, p_
);
893 const complex double multiplier
= c
->B_multipliers
[i
][q
-BQ_OFFSET
];
894 complex double sigma
= sigmas_total
[y_sc
];
895 ckahanadd(&Bsum
, &Bsumc
, multiplier
* sigma
);
896 if (Berr
) kahanadd(&Bsumerr
, &Bsumerrc
, multiplier
* serr_total
[y_sc
]);
899 *(Bdest
+ deststride
* desti
+ srcstride
* srci
) = Bsum
;
900 if (Berr
) *(Berr
+ deststride
* desti
+ srcstride
* srci
) = Bsumerr
;
919 #endif // LATTICESUMS_31
924 // N.B. alternative point generation strategy toggled by macro GEN_RSHIFTEDPOINTS
925 // and GEN_KSHIFTEDPOINTS.
926 // The results should be the same. The performance can slightly differ (especially
927 // if some optimizations in the point generators are implemented.)
928 int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator
*c
,
929 complex double * const Adest
, double * const Aerr
,
930 complex double * const Bdest
, double * const Berr
,
931 const ptrdiff_t deststride
, const ptrdiff_t srcstride
,
932 /* qpms_bessel_t J*/ // assume QPMS_HANKEL_PLUS
933 const double eta
, const complex double k
,
934 const cart2_t b1
, const cart2_t b2
,
936 const cart2_t particle_shift
,
937 double maxR
, double maxK
,
938 const qpms_ewald_part parts
942 const qpms_y_t nelem2_sc
= qpms_lMax2nelem_sc(c
->e3c
->lMax
);
943 //const qpms_y_t nelem = qpms_lMax2nelem(c->lMax);
944 const bool doerr
= Aerr
|| Berr
;
945 const bool do_sigma0
= ((particle_shift
.x
== 0) && (particle_shift
.y
== 0)); // FIXME ignoring the case where particle_shift equals to lattice vector
947 complex double *sigmas_short
= malloc(sizeof(complex double)*nelem2_sc
);
948 complex double *sigmas_long
= malloc(sizeof(complex double)*nelem2_sc
);
949 complex double *sigmas_total
= malloc(sizeof(complex double)*nelem2_sc
);
950 double *serr_short
, *serr_long
, *serr_total
;
952 serr_short
= malloc(sizeof(double)*nelem2_sc
);
953 serr_long
= malloc(sizeof(double)*nelem2_sc
);
954 serr_total
= malloc(sizeof(double)*nelem2_sc
);
955 } else serr_short
= serr_long
= serr_total
= NULL
;
957 const double unitcell_area
= l2d_unitcell_area(b1
, b2
);
958 cart2_t rb1
, rb2
; // reciprocal basis
959 QPMS_ENSURE_SUCCESS(l2d_reciprocalBasis2pi(b1
, b2
, &rb1
, &rb2
));
961 if (parts
& QPMS_EWALD_LONG_RANGE
) {
962 PGen Kgen
= PGen_xyWeb_new(rb1
, rb2
, BASIS_RTOL
,
963 #ifdef GEN_KSHIFTEDPOINTS
968 0, true, maxK
, false);
970 QPMS_ENSURE_SUCCESS(ewald3_sigma_long(sigmas_long
, serr_long
, c
->e3c
, eta
, k
,
971 unitcell_area
, LAT_2D_IN_3D_XYONLY
, &Kgen
,
972 #ifdef GEN_KSHIFTEDPOINTS
977 cart22cart3xy(beta
), cart22cart3xy(particle_shift
)));
978 if(Kgen
.stateData
) // PGen not consumed entirely (converged earlier)
982 if (parts
& QPMS_EWALD_SHORT_RANGE
) {
983 PGen Rgen
= PGen_xyWeb_new(b1
, b2
, BASIS_RTOL
,
984 #ifdef GEN_RSHIFTEDPOINTS
985 cart2_scale(-1 /*CHECKSIGN*/, particle_shift
),
989 0, !do_sigma0
, maxR
, false);
991 QPMS_ENSURE_SUCCESS(ewald3_sigma_short(sigmas_short
, serr_short
, c
->e3c
, eta
, k
,
992 LAT_2D_IN_3D_XYONLY
, &Rgen
,
993 #ifdef GEN_RSHIFTEDPOINTS
998 cart22cart3xy(beta
), cart22cart3xy(particle_shift
)));
1000 if(Rgen
.stateData
) // PGen not consumed entirely (converged earlier)
1001 PGen_destroy(&Rgen
);
1004 for(qpms_y_t y
= 0; y
< nelem2_sc
; ++y
)
1005 sigmas_total
[y
] = ((parts
& QPMS_EWALD_SHORT_RANGE
) ? sigmas_short
[y
] : 0)
1006 + ((parts
& QPMS_EWALD_LONG_RANGE
) ? sigmas_long
[y
] : 0);
1007 if (doerr
) for(qpms_y_t y
= 0; y
< nelem2_sc
; ++y
)
1008 serr_total
[y
] = ((parts
& QPMS_EWALD_SHORT_RANGE
) ? serr_short
[y
] : 0)
1009 + ((parts
& QPMS_EWALD_LONG_RANGE
) ? serr_long
[y
] : 0);
1011 complex double sigma0
= 0; double sigma0_err
= 0;
1012 if (do_sigma0
&& (parts
& QPMS_EWALD_0TERM
)) {
1013 QPMS_ENSURE_SUCCESS(ewald3_sigma0(&sigma0
, &sigma0_err
, c
->e3c
, eta
, k
));
1014 const qpms_l_t y
= qpms_mn2y_sc(0,0);
1015 sigmas_total
[y
] += sigma0
;
1016 if(doerr
) serr_total
[y
] += sigma0_err
;
1020 ptrdiff_t desti
= 0, srci
= 0;
1021 for (qpms_l_t n
= 1; n
<= c
->lMax
; ++n
) for (qpms_m_t m
= -n
; m
<= n
; ++m
) {
1022 for (qpms_l_t nu
= 1; nu
<= c
->lMax
; ++nu
) for (qpms_m_t mu
= -nu
; mu
<= nu
; ++mu
){
1023 const size_t i
= qpms_trans_calculator_index_mnmunu(c
, m
, n
, mu
, nu
);
1024 const size_t qmax
= c
->A_multipliers
[i
+1] - c
->A_multipliers
[i
] - 1;
1025 complex double Asum
, Asumc
; ckahaninit(&Asum
, &Asumc
);
1026 double Asumerr
, Asumerrc
; if(Aerr
) kahaninit(&Asumerr
, &Asumerrc
);
1028 const qpms_m_t mu_m
= mu
- m
;
1030 for(qpms_l_t q
= 0; q
<= qmax
; ++q
) {
1031 const qpms_l_t p
= n
+ nu
- 2*q
;
1032 const qpms_y_t y_sc
= qpms_mn2y_sc(mu_m
, p
);
1033 const complex double multiplier
= c
->A_multipliers
[i
][q
];
1034 complex double sigma
= sigmas_total
[y_sc
];
1035 ckahanadd(&Asum
, &Asumc
, multiplier
* sigma
);
1036 if (Aerr
) kahanadd(&Asumerr
, &Asumerrc
, multiplier
* serr_total
[y_sc
]);
1039 *(Adest
+ deststride
* desti
+ srcstride
* srci
) = Asum
;
1040 if (Aerr
) *(Aerr
+ deststride
* desti
+ srcstride
* srci
) = Asumerr
;
1043 complex double Bsum
, Bsumc
; ckahaninit(&Bsum
, &Bsumc
);
1044 double Bsumerr
, Bsumerrc
; if(Berr
) kahaninit(&Bsumerr
, &Bsumerrc
);
1045 for(qpms_l_t q
= 0; q
<= qmax
; ++q
) {
1046 const qpms_l_t p_
= n
+ nu
- 2*q
+ 1;
1047 const qpms_y_t y_sc
= qpms_mn2y_sc(mu_m
, p_
);
1048 const complex double multiplier
= c
->B_multipliers
[i
][q
-BQ_OFFSET
];
1049 complex double sigma
= sigmas_total
[y_sc
];
1050 ckahanadd(&Bsum
, &Bsumc
, multiplier
* sigma
);
1051 if (Berr
) kahanadd(&Bsumerr
, &Bsumerrc
, multiplier
* serr_total
[y_sc
]);
1054 *(Bdest
+ deststride
* desti
+ srcstride
* srci
) = Bsum
;
1055 if (Berr
) *(Berr
+ deststride
* desti
+ srcstride
* srci
) = Bsumerr
;
1075 int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator
*c
,
1076 complex double * const Adest
, double * const Aerr
,
1077 complex double * const Bdest
, double * const Berr
,
1078 const ptrdiff_t deststride
, const ptrdiff_t srcstride
,
1079 /* qpms_bessel_t J*/ // assume QPMS_HANKEL_PLUS
1080 const double eta
, const complex double k
,
1081 const cart2_t b1
, const cart2_t b2
,
1083 const cart2_t particle_shift
,
1084 double maxR
, double maxK
)
1086 return qpms_trans_calculator_get_AB_arrays_e32_e(
1087 c
, Adest
, Aerr
, Bdest
, Berr
, deststride
, srcstride
,
1088 eta
, k
, b1
, b2
, beta
, particle_shift
, maxR
, maxK
, QPMS_EWALD_FULL
);
1091 #endif // LATTICESUMS32
1094 complex double qpms_trans_calculator_get_A_ext(const qpms_trans_calculator
*c
,
1095 int m
, int n
, int mu
, int nu
,
1096 complex double kdlj_r
, double kdlj_theta
, double kdlj_phi
,
1097 int r_ge_d
, int J
) {
1098 csph_t kdlj
= {kdlj_r
, kdlj_theta
, kdlj_phi
};
1099 return qpms_trans_calculator_get_A(c
,m
,n
,mu
,nu
,kdlj
,r_ge_d
,J
);
1102 complex double qpms_trans_calculator_get_B_ext(const qpms_trans_calculator
*c
,
1103 int m
, int n
, int mu
, int nu
,
1104 complex double kdlj_r
, double kdlj_theta
, double kdlj_phi
,
1105 int r_ge_d
, int J
) {
1106 csph_t kdlj
= {kdlj_r
, kdlj_theta
, kdlj_phi
};
1107 return qpms_trans_calculator_get_B(c
,m
,n
,mu
,nu
,kdlj
,r_ge_d
,J
);
1110 int qpms_trans_calculator_get_AB_p_ext(const qpms_trans_calculator
*c
,
1111 complex double *Adest
, complex double *Bdest
,
1112 int m
, int n
, int mu
, int nu
,
1113 complex double kdlj_r
, double kdlj_theta
, double kdlj_phi
,
1114 int r_ge_d
, int J
) {
1115 csph_t kdlj
= {kdlj_r
, kdlj_theta
, kdlj_phi
};
1116 return qpms_trans_calculator_get_AB_p(c
,Adest
,Bdest
,m
,n
,mu
,nu
,kdlj
,r_ge_d
,J
);
1119 int qpms_trans_calculator_get_AB_arrays_ext(const qpms_trans_calculator
*c
,
1120 complex double *Adest
, complex double *Bdest
,
1121 size_t deststride
, size_t srcstride
,
1122 complex double kdlj_r
, double kdlj_theta
, double kdlj_phi
,
1123 int r_ge_d
, int J
) {
1124 csph_t kdlj
= {kdlj_r
, kdlj_theta
, kdlj_phi
};
1125 return qpms_trans_calculator_get_AB_arrays(c
,Adest
,Bdest
,deststride
,srcstride
,