From bb289a46b9fa1f5e1b68637b53fb0265dbf21d45 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Sun, 17 Dec 2017 17:42:18 +0200 Subject: [PATCH] C indexing functions to separate file, new normalizations in progress Former-commit-id: 30b4ab302fe69fedc8c12bdb7b3b35f532cfc35d --- qpms/indexing.h | 27 +++++++ qpms/qpms_types.h | 2 + qpms/translations.c | 175 ++++++++++++++++++++++++++++++++++++++++---- 3 files changed, 188 insertions(+), 16 deletions(-) create mode 100644 qpms/indexing.h diff --git a/qpms/indexing.h b/qpms/indexing.h new file mode 100644 index 0000000..8cabaf1 --- /dev/null +++ b/qpms/indexing.h @@ -0,0 +1,27 @@ +#ifndef QPMS_INDEXING_H +#define QPMS_INDEXING_H + +#include "qpms_types.h" + +static inline qpms_y_t qpms_mn2y(qpms_m_t m, qpms_l_t n) { + return (qpms_y_t) n * (n + 1) + m - 1; +} + +static inline qpms_lm_t qpms_y2n(qpms_y_t y) { + //return (sqrt(5+y)-2)/2; // the cast will truncate the fractional part, which is what we want + return sqrt(y+1); +} + +static inline qpms_m_t qpms_yn2m(qpms_y_t y, qpms_l_t n) { + return y-qpms_mn2y(0,n); +} + +static inline void qpms_y2mn_p(qpms_y_t y, qpms_m_t *m, qpms_l_t *n){ + *m=qpms_yn2m(y,*n=qpms_y2n(y)); +} + +static inline qpms_y_t qpms_lMax2nelem(qpms_l_t lmax){ + return lmax * ((qpms_y_t)lmax + 2); +} + +#endif //QPMS_INDEXING_H diff --git a/qpms/qpms_types.h b/qpms/qpms_types.h index 4ce5e14..28e5352 100644 --- a/qpms/qpms_types.h +++ b/qpms/qpms_types.h @@ -7,6 +7,8 @@ // integer index types typedef int qpms_lm_t; +typedef unsigned int qpms_l_t; +typedef qpms_lm_t qpms_m_t; typedef size_t qpms_y_t; // Normalisations diff --git a/qpms/translations.c b/qpms/translations.c index 07c493e..8bc2422 100644 --- a/qpms/translations.c +++ b/qpms/translations.c @@ -1,11 +1,17 @@ #include #include "gaunt.h" #include "translations.h" +#include "indexing.h" // TODO replace size_t and int with own index types here #include #include #include #include "assert_cython_workaround.h" +#include //abort() +/* + * References: + * [1] Yu-Lin Xu, Journal of Computational Physics 127, 285–298 (1996) + */ static const double sqrtpi = 1.7724538509055160272981674833411451827975494561223871; //static const double ln2 = 0.693147180559945309417232121458176568075500134360255254120; @@ -61,9 +67,10 @@ int qpms_sph_bessel_array(qpms_bessel_t typ, int lmax, double x, complex double assert(0); } - -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) { // TODO make J enum +// [1], eq. (82) +complex double qpms_trans_single_A_Xu(int m, int n, int mu, int nu, sph_t kdlj, + bool r_ge_d, qpms_bessel_t J) { + abort(); // FIXME, THIS IS STILL TAYLOR if(r_ge_d) J = QPMS_BESSEL_REGULAR; double costheta = cos(kdlj.theta); @@ -105,6 +112,155 @@ complex double qpms_trans_single_A_Taylor(int m, int n, int mu, int nu, sph_t kd return (presum / prenormratio) * sum; } + + +complex double qpms_trans_single_A_Kristensson(int m, int n, int mu, int nu, sph_t kdlj, + bool r_ge_d, qpms_bessel_t J) { + abort();// FIXME, THIS IS STILL TAYLOR + 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_array(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; + } + + 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; +} + + + +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_array(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; + } + + 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; +} + +// [1], 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) { + abort(); // FIXME, this is still Taylor + 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); + 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_array(J, n+nu+1, kdlj.r, bes)) abort(); + + complex double sum = 0; + for (int q = 0; q <= Qmax; ++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; + } + + 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) * sqrt(((2.*nu+1)/(2.*n+1))* exp( + lgamma(n+m+1)-lgamma(n-m+1)+lgamma(nu-mu+1)-lgamma(nu+mu+1))); + + + 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) { if(r_ge_d) J = QPMS_BESSEL_REGULAR; @@ -186,19 +342,6 @@ static inline size_t qpms_mn2y(int m, int n) { return (size_t) n * (n + 1) + m - 1; } -static inline int qpms_y2n(size_t y) { - //return (sqrt(5+y)-2)/2; // the cast will truncate the fractional part, which is what we want - return sqrt(y+1); -} - -static inline int qpms_yn2m(size_t y, int n) { - return y-qpms_mn2y(0,n); -} - -static inline void qpms_y2mn_p(size_t y, int *m, int *n){ - *m=qpms_yn2m(y,*n=qpms_y2n(y)); -} - static inline size_t qpms_trans_calculator_index_mnmunu(const qpms_trans_calculator *c, int m, int n, int mu, int nu){ return c->nelem * qpms_mn2y(m,n) + qpms_mn2y(mu,nu);