diff --git a/qpms/test_translations2.c b/qpms/test_translations2.c new file mode 100644 index 0000000..ff7f958 --- /dev/null +++ b/qpms/test_translations2.c @@ -0,0 +1,54 @@ +#include "translations.h" +#include +//#include +#include + +typedef struct { + int m, n, mu, nu; + sph_t kdlj; + qpms_bessel_t J; + complex double result_A, result_B; +} testcase_single_trans_t; + +testcase_single_trans_t testcases_Taylor[] = { +#include "testcases_taylor" +}; + +int lMax=10; + +int main() { + qpms_trans_calculator *c = qpms_trans_calculator_init(lMax, QPMS_NORMALIZATION_TAYLOR); + + for(testcase_single_trans_t *tc = testcases_Taylor; tc->J != QPMS_BESSEL_UNDEF; tc++) { + if (!tc->n || !tc->nu || tc->n > lMax || tc->nu > lMax ) 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); + complex double A2 = qpms_trans_calculator_get_A(c, tc->m, tc->n, tc->mu, tc->nu, tc->kdlj, true, tc->J); + complex double B2 = qpms_trans_calculator_get_B(c, 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), (0 == cabs(tc->result_A - A)) ? 0 : + cabs(tc->result_A - A)/((cabs(A) < cabs(tc->result_A)) ? cabs(A) : cabs(tc->result_A)), + tc->J); + printf("A' = %.16f+%.16fj, relerr=%.16f, relerr2=%.3e\n", + creal(A2), cimag(A2), (0 == cabs(tc->result_A - A2)) ? 0 : + cabs(tc->result_A - A2)/((cabs(A2) < cabs(tc->result_A)) ? cabs(A2) : cabs(tc->result_A)), + (0 == cabs(A - A2)) ? 0 : + cabs(A - A2)/((cabs(A2) < cabs(A)) ? cabs(A2) : cabs(A)) + ); + printf("B = %.16f+%.16fj, relerr=%.16f, J=%d\n", + creal(B), cimag(B), (0 == cabs(tc->result_B - B)) ? 0 : + cabs(tc->result_B - B)/((cabs(B) < cabs(tc->result_B)) ? cabs(B) : cabs(tc->result_B)), + tc->J); + printf("B' = %.16f+%.16fj, relerr=%.16f, relerr2=%.3e\n", + creal(B2), cimag(B2), (0 == cabs(tc->result_B - B2)) ? 0 : + cabs(tc->result_B - B2)/((cabs(B2) < cabs(tc->result_B)) ? cabs(B2) : cabs(tc->result_B)), + (0 == cabs(B - B2)) ? 0 : + cabs(B - B2)/((cabs(B2) < cabs(B)) ? cabs(B2) : cabs(B)) + ); + } + qpms_trans_calculator_free(c); +} + + diff --git a/qpms/translations.c b/qpms/translations.c index 9821706..9ded01e 100644 --- a/qpms/translations.c +++ b/qpms/translations.c @@ -84,7 +84,8 @@ complex double qpms_trans_single_A_Taylor(int m, int n, int mu, int nu, sph_t kd 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 + //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)); @@ -105,7 +106,7 @@ complex double qpms_trans_single_A_Taylor(int m, int n, int mu, int nu, sph_t kd } 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 + bool r_ge_d, qpms_bessel_t J) { if(r_ge_d) J = QPMS_BESSEL_REGULAR; double costheta = cos(kdlj.theta); @@ -186,7 +187,8 @@ static inline size_t qpms_mn2y(int m, int n) { } 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(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) { @@ -207,8 +209,148 @@ static inline size_t qpms_trans_calculator_index_yyu(const qpms_trans_calculator return c->nelem * y + yu; } + #define SQ(x) ((x)*(x)) +//#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(Qmax == gaunt_q_max(-m,n+1,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[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 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 = 0; 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] = Ppfac * summandq * presum; + } +} + +int qpms_trans_calculator_multipliers_A(qpms_normalization_t norm, complex double *dest, int m, int n, int mu, int nu, int qmax) { + switch (norm) { + case QPMS_NORMALIZATION_TAYLOR: + qpms_trans_calculator_multipliers_A_Taylor(dest,m,n,mu,nu,qmax); + return 0; + break; + default: + abort(); + } + assert(0); +} + +int qpms_trans_calculator_multipliers_B(qpms_normalization_t norm, complex double *dest, int m, int n, int mu, int nu, int qmax) { + switch (norm) { + case QPMS_NORMALIZATION_TAYLOR: + qpms_trans_calculator_multipliers_B_Taylor(dest,m,n,mu,nu,qmax); + return 0; + break; + default: + abort(); + } + assert(0); +} + qpms_trans_calculator *qpms_trans_calculator_init (int lMax, qpms_normalization_t normalization) { qpms_trans_calculator *c = malloc(sizeof(qpms_trans_calculator)); @@ -224,14 +366,24 @@ qpms_trans_calculator int m,n, mu, nu; qpms_y2mn_p(y,&m,&n); qpms_y2mn_p(yu,&mu,&nu); - qmaxsum += (qmaxes[qpms_trans_calculator_index_yyu(c,y,yu)] = gaunt_q_max(-m,n,mu,nu)); + qmaxsum += 1 + ( + qmaxes[qpms_trans_calculator_index_yyu(c,y,yu)] + = gaunt_q_max(-m,n,mu,nu)); } c->A_multipliers[0] = malloc(qmaxsum * sizeof(complex double)); - for(size_t i = 0; i < SQ(c->nelem); ++i) - c->A_multipliers[i+1] = c->A_multipliers[i] + qmaxes[i]; - // TODO here comes the filling of A_multipliers - - + // calculate multiplier beginnings + for(size_t i = 0; i < SQ(c->nelem); ++i) + c->A_multipliers[i+1] = c->A_multipliers[i] + qmaxes[i] + 1; + // calculate the multipliers + for(size_t y = 0; y < c->nelem; ++y) + for(size_t yu = 0; yu < c->nelem; ++yu) { + size_t i = y * c->nelem + yu; + int m, n, mu, nu; + qpms_y2mn_p(y, &m, &n); + qpms_y2mn_p(yu, &mu, &nu); + qpms_trans_calculator_multipliers_A(normalization, + c->A_multipliers[i], m, n, mu, nu, qmaxes[i]); + } qmaxsum = 0; for(size_t y=0; y < c->nelem; y++) @@ -239,16 +391,122 @@ qpms_trans_calculator int m, n, mu, nu; qpms_y2mn_p(y,&m,&n); qpms_y2mn_p(yu,&mu,&nu); - qmaxsum += (qmaxes[qpms_trans_calculator_index_yyu(c,y,yu)] = gaunt_q_max(-m,n+1,mu,nu)); + qmaxsum += 1 + ( + qmaxes[qpms_trans_calculator_index_yyu(c,y,yu)] + = gaunt_q_max(-m,n+1,mu,nu)); } c->B_multipliers[0] = malloc(qmaxsum * sizeof(complex double)); - for(size_t i = 0; i < SQ(c->nelem); ++i) - c->B_multipliers[i+1] = c->B_multipliers[i] + qmaxes[i]; - - // TODO here comes the filling of B_multipliers + // calculate multiplier beginnings + for(size_t i = 0; i < SQ(c->nelem); ++i) + c->B_multipliers[i+1] = c->B_multipliers[i] + qmaxes[i] + 1; + // calculate the multipliers + for(size_t y = 0; y < c->nelem; ++y) + for(size_t yu = 0; yu < c->nelem; ++yu) { + size_t i = y * c->nelem + yu; + int m, n, mu, nu; + qpms_y2mn_p(y, &m, &n); + qpms_y2mn_p(yu, &mu, &nu); + qpms_trans_calculator_multipliers_B(normalization, + c->B_multipliers[i], m, n, mu, nu, qmaxes[i]); + } free(qmaxes); return c; } +complex double qpms_trans_calculator_get_A_buf(const qpms_trans_calculator *c, + int m, int n, int mu, int nu, sph_t kdlj, + bool r_ge_d, qpms_bessel_t J, + complex double *bessel_buf, double *legendre_buf) { + if (r_ge_d) J = QPMS_BESSEL_REGULAR; + switch(c->normalization) { + case QPMS_NORMALIZATION_TAYLOR: + { + double costheta = cos(kdlj.theta); + if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu,costheta,-1,legendre_buf)) abort(); + if (qpms_sph_bessel_array(J, n+nu, kdlj.r, bessel_buf)) abort(); + 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 = 0; + for(size_t q = 0; q <= qmax; ++q) { + int p = n+nu-2*q; + double Pp = legendre_buf[gsl_sf_legendre_array_index(p, abs(mu-m))]; + complex double zp = bessel_buf[p]; + complex double multiplier = c->A_multipliers[i][q]; + sum += Pp * zp * multiplier; + } + complex double eimf = cexp(I*(mu-m)*kdlj.phi); + return sum * eimf; + } + break; + default: + abort(); + } + assert(0); +} +complex double qpms_trans_calculator_get_B_buf(const qpms_trans_calculator *c, + int m, int n, int mu, int nu, sph_t kdlj, + bool r_ge_d, qpms_bessel_t J, + complex double *bessel_buf, double *legendre_buf) { + if (r_ge_d) J = QPMS_BESSEL_REGULAR; + switch(c->normalization) { + case QPMS_NORMALIZATION_TAYLOR: + { + 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_array(J, n+nu+2, kdlj.r, bessel_buf)) abort(); + 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; + assert(qmax == gaunt_q_max(-m,n+1,mu,nu)); + complex double sum = 0; + for(int q = 0; q <= qmax; ++q) { + int p = n+nu-2*q; + double Pp_ = legendre_buf[gsl_sf_legendre_array_index(p+1, abs(mu-m))]; + complex double eimf = cexp(I * kdlj.phi); + complex double zp_ = bessel_buf[p+1]; + complex double multiplier = c->B_multipliers[i][q]; + sum += Pp_ * zp_ * multiplier; + } + complex double eimf = cexp(I*(mu-m)*kdlj.phi); + return sum * eimf; + } + break; + default: + abort(); + } + assert(0); +} + +complex double qpms_trans_calculator_get_A(const qpms_trans_calculator *c, + int m, int n, int mu, int nu, sph_t kdlj, + bool r_ge_d, qpms_bessel_t J) { + double leg[gsl_sf_legendre_array_n(n+nu)]; + complex double bes[n+nu+1]; + return qpms_trans_calculator_get_A_buf(c,m,n,mu,nu,kdlj,r_ge_d,J, + bes,leg); +} + +complex double qpms_trans_calculator_get_B(const qpms_trans_calculator *c, + int m, int n, int mu, int nu, sph_t kdlj, + bool r_ge_d, qpms_bessel_t J) { + double leg[gsl_sf_legendre_array_n(n+nu+1)]; + complex double bes[n+nu+2]; + return qpms_trans_calculator_get_B_buf(c,m,n,mu,nu,kdlj,r_ge_d,J, + bes,leg); +} + +#if 0 +int qpms_trans_calculator_get_AB_p(const qpms_trans_calculator *c, + complex double *Adest, complex double *Bdest, + int m, int n, int mu, int nu, sph_t kdlj, + bool r_ge_d, qpms_bessel_t J) { + double leg[gsl_sf_legendre_array_n(n+nu+1)]; + complex double bes[n+nu+2]; + return qpms_trans_calculator_get_AB_buf(c,Adest,Bdest, + m,n,mu,nu,kdlj,r_ge_d,J, + bes,leg); +} +#endif diff --git a/qpms/translations.h b/qpms/translations.h index 3ac1eb9..3ca425b 100644 --- a/qpms/translations.h +++ b/qpms/translations.h @@ -46,4 +46,18 @@ typedef struct qpms_trans_calculator { qpms_trans_calculator *qpms_trans_calculator_init(int lMax, qpms_normalization_t nt); void qpms_trans_calculator_free(qpms_trans_calculator *); +complex double qpms_trans_calculator_get_A(const qpms_trans_calculator *c, + int m, int n, int mu, int nu, sph_t kdlj, + bool r_ge_d, qpms_bessel_t J); +complex double qpms_trans_calculator_get_B(const qpms_trans_calculator *c, + int m, int n, int mu, int nu, sph_t kdlj, + bool r_ge_d, qpms_bessel_t J); +#if 0 +int qpms_trans_calculator_get_AB_p(const qpms_trans_calculator *c, + complex double *Adest, complex double *Bdest, + int m, int n, int mu, int nu, sph_t kdlj, + bool r_ge_d, qpms_bessel_t J); +#endif + + #endif // QPMS_TRANSLATIONS_H