From bd5aacd7d9c82eddc5a47bba695aae1c3c8ab6d0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 29 Jan 2018 10:13:13 +0000 Subject: [PATCH] Legendre, pi, tau calculation fixes Former-commit-id: 98e265cc0c6d73bd3eff3a601096ce84d63cc7a2 --- qpms/test_vswf.c | 2 +- qpms/vswf.c | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/qpms/test_vswf.c b/qpms/test_vswf.c index 5a307e0..12d145a 100644 --- a/qpms/test_vswf.c +++ b/qpms/test_vswf.c @@ -8,7 +8,7 @@ const qpms_l_t lMax = 3; int main() { qpms_y_t nelem = qpms_lMax2nelem(lMax); for (double theta = 0.; theta <= M_PI; theta += dtheta) { - printf("%.5e ", theta); + printf("%.5e %.5e ", theta, cos(theta)); for(qpms_normalisation_t norm = 1; norm <= 3; ++norm) {//fujka :D qpms_pitau_t pt = qpms_pitau_get(theta, lMax, norm); for (qpms_y_t y = 0; y < nelem; ++y) diff --git a/qpms/vswf.c b/qpms/vswf.c index 2ac55bd..c58ed19 100644 --- a/qpms/vswf.c +++ b/qpms/vswf.c @@ -33,7 +33,7 @@ qpms_errno_t qpms_legendre_deriv_y_fill(double *target, double *target_deriv, do qpms_y_t y = qpms_mn2y(-m,l); size_t i = gsl_sf_legendre_array_index(l,m); // viz DLMF 14.9.3, čert ví, jak je to s cs fasí. - double factor = exp(lgamma(l-m+1)-lgamma(n+m+1))*((m%2)?-1:1); + double factor = exp(lgamma(l-m+1)-lgamma(l+m+1))*((m%2)?-1:1); target[y] = factor * legendre_tmp[i]; target_deriv[y] = factor * legendre_deriv_tmp[i]; } @@ -85,8 +85,8 @@ qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, qpms_normalisation_t no case QPMS_NORMALISATION_XU: for (qpms_l_t l = 1; l <= lMax; ++l) { res.leg[qpms_mn2y(0, l)] = (l%2)?ct:1.; - double p = l*(l+1)*(l+2)/2; - double n = (l+2)/2.; + double p = l*(l+1)/2; + const double n = 0.5; int lpar = (l%2)?-1:1; res.pi [qpms_mn2y(+1, l)] = ((ct>0) ? -1 : lpar) * p; res.pi [qpms_mn2y(-1, l)] = ((ct>0) ? -1 : lpar) * n * CSPHASE; @@ -98,7 +98,7 @@ qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, qpms_normalisation_t no for (qpms_l_t l = 1; l <= lMax; ++l) { res.leg[qpms_mn2y(0, l)] = ((l%2)?ct:1.)*sqrt((2*l+1)*0.25*M_1_PI); int lpar = (l%2)?-1:1; - double fl = 0.25 * (l+2) * sqrt((2*l+1)*l*(l+1)*M_1_PI); + double fl = 0.25 * sqrt((2*l+1)*l*(l+1)*M_1_PI); res.pi [qpms_mn2y(+1, l)] = ((ct>0) ? -1 : lpar) * fl; res.pi [qpms_mn2y(-1, l)] = ((ct>0) ? -1 : lpar) * fl * CSPHASE; res.tau[qpms_mn2y(+1, l)] = ((ct>0) ? +1 : lpar) * fl; @@ -109,7 +109,7 @@ qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, qpms_normalisation_t no for (qpms_l_t l = 1; l <= lMax; ++l) { res.leg[qpms_mn2y(0, l)] = ((l%2)?ct:1.)*sqrt((2*l+1)/(4*M_PI *l*(l+1))); int lpar = (l%2)?-1:1; - double fl = 0.25 * (l+2) * sqrt((2*l+1)*M_1_PI); + double fl = 0.25 * sqrt((2*l+1)*M_1_PI); res.pi [qpms_mn2y(+1, l)] = ((ct>0) ? -1 : lpar) * fl; res.pi [qpms_mn2y(-1, l)] = ((ct>0) ? -1 : lpar) * fl * CSPHASE; res.tau[qpms_mn2y(+1, l)] = ((ct>0) ? +1 : lpar) * fl;