From 22dd5ef87486506c6fe59a449b2656da188de027 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Thu, 13 Apr 2017 11:29:38 +0300 Subject: [PATCH] C implementation of the B translation coefficient. Former-commit-id: f7c8930de88c1d62ac99e29e9adda4ca78aa216e --- qpms/BUGS.rst | 6 ++++ qpms/test_translations.c | 16 +++++++-- qpms/translations.c | 72 ++++++++++++++++++++++++++++++++++++---- qpms/translations.h | 2 ++ 4 files changed, 86 insertions(+), 10 deletions(-) create mode 100644 qpms/BUGS.rst diff --git a/qpms/BUGS.rst b/qpms/BUGS.rst new file mode 100644 index 0000000..046dd66 --- /dev/null +++ b/qpms/BUGS.rst @@ -0,0 +1,6 @@ +gaunt.c +======= +abort při určitých vstupech, např. -6 10 8 10 +(ačkoliv fortran originál paří bez problémů) + + diff --git a/qpms/test_translations.c b/qpms/test_translations.c index a6210c3..c690523 100644 --- a/qpms/test_translations.c +++ b/qpms/test_translations.c @@ -16,9 +16,19 @@ testcase_single_trans_t testcases_Taylor[] = { int main() { for(testcase_single_trans_t *tc = testcases_Taylor; tc->J != QPMS_BESSEL_UNDEF; tc++) { - complex double A = qpms_trans_single_A_Taylor(tc->m, tc->n, tc->mu, tc->nu, tc->kdlj, true, tc->J); - printf("m=%d, n=%d, mu=%d, nu=%d, relerr=%.16f\n", tc->m,tc->n,tc->mu,tc->nu, - cabs(tc->result_A - A)/((cabs(A) < cabs(tc->result_A)) ? cabs(A) : cabs(tc->result_A))); + if (tc->n > 12 || tc->nu > 12 || !tc->n || !tc->nu ) continue; + + printf("m=%d, n=%d, mu=%d, nu=%d,\n", tc->m,tc->n,tc->mu,tc->nu); + complex double A = qpms_trans_single_A_Taylor(tc->m, tc->n, tc->mu, tc->nu, tc->kdlj, true, tc->J); + complex double B = qpms_trans_single_B_Taylor(tc->m, tc->n, tc->mu, tc->nu, tc->kdlj, true, tc->J); + printf("A = %.16f+%.16fj, relerr=%.16f, J=%d\n", + creal(A), cimag(A), + cabs(tc->result_A - A)/((cabs(A) < cabs(tc->result_A)) ? cabs(A) : cabs(tc->result_A)), + tc->J); + printf("B = %.16f+%.16fj, relerr=%.16f, J=%d\n", + creal(B), cimag(B), + cabs(tc->result_B - B)/((cabs(B) < cabs(tc->result_B)) ? cabs(B) : cabs(tc->result_B)), + tc->J); } } diff --git a/qpms/translations.c b/qpms/translations.c index 53a768e..8b35505 100644 --- a/qpms/translations.c +++ b/qpms/translations.c @@ -65,10 +65,9 @@ int qpms_sph_bessel_array(qpms_bessel_t typ, int lmax, double x, complex double 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 if(r_ge_d) J = QPMS_BESSEL_REGULAR; - 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)); + 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]; @@ -76,8 +75,7 @@ complex double qpms_trans_single_A_Taylor(int m, int n, int mu, int nu, sph_t kd gaunt_xu(-m,n,mu,nu,qmax,a1q,&err); double a1q0 = a1q[0]; if (err) abort(); - //double *leg = malloc(sizeof(double)*gsl_sf_legendre_array_n(n+nu)); - //if (!leg) 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]; @@ -94,8 +92,10 @@ complex double qpms_trans_single_A_Taylor(int m, int n, int mu, int nu, sph_t kd complex double summandq = (n*(n+1) + nu*(nu+1) - p*(p+1)) * min1pow(q) * a1q_n * zp * Pp; sum += summandq; } - //free(leg); - //free(bes); + + 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); @@ -104,3 +104,61 @@ 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_B_Taylor(int m, int n, int mu, int nu, sph_t kdlj, + bool r_ge_d, qpms_bessel_t J) { // TODO make J enum + 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 + 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))); + + return (presum / prenormratio) * sum; +} + + + diff --git a/qpms/translations.h b/qpms/translations.h index bc673ea..78cb325 100644 --- a/qpms/translations.h +++ b/qpms/translations.h @@ -20,5 +20,7 @@ int qpms_sph_bessel_array(qpms_bessel_t typ, int lmax, double x, complex double 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); +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); #endif // QPMS_TRANSLATIONS_H