QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Legendre, pi, tau calculation fixes 

Former-commit-id: 98e265cc0c6d73bd3eff3a601096ce84d63cc7a2
This commit is contained in:
Marek Nečada 2018-01-29 10:13:13 +00:00
parent 92e9622b64
commit bd5aacd7d9
2 changed files with 6 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
qpms/test_vswf.c
View File

@ -8,7 +8,7 @@ const qpms_l_t lMax = 3;
int main() { int main() {
qpms_y_t nelem = qpms_lMax2nelem(lMax); qpms_y_t nelem = qpms_lMax2nelem(lMax);
for (double theta = 0.; theta <= M_PI; theta += dtheta) { 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 for(qpms_normalisation_t norm = 1; norm <= 3; ++norm) {//fujka :D
qpms_pitau_t pt = qpms_pitau_get(theta, lMax, norm); qpms_pitau_t pt = qpms_pitau_get(theta, lMax, norm);
for (qpms_y_t y = 0; y < nelem; ++y) for (qpms_y_t y = 0; y < nelem; ++y)

10
qpms/vswf.c
View File

@ -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); qpms_y_t y = qpms_mn2y(-m,l);
size_t i = gsl_sf_legendre_array_index(l,m); size_t i = gsl_sf_legendre_array_index(l,m);
// viz DLMF 14.9.3, čert ví, jak je to s cs fasí. // 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[y] = factor * legendre_tmp[i];
target_deriv[y] = factor * legendre_deriv_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: case QPMS_NORMALISATION_XU:
for (qpms_l_t l = 1; l <= lMax; ++l) { for (qpms_l_t l = 1; l <= lMax; ++l) {
res.leg[qpms_mn2y(0, l)] = (l%2)?ct:1.; res.leg[qpms_mn2y(0, l)] = (l%2)?ct:1.;
double p = l*(l+1)*(l+2)/2; double p = l*(l+1)/2;
double n = (l+2)/2.; const double n = 0.5;
int lpar = (l%2)?-1:1; 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) * p;
res.pi [qpms_mn2y(-1, l)] = ((ct>0) ? -1 : lpar) * n * CSPHASE; 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) { 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); 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; 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;
res.pi [qpms_mn2y(-1, l)] = ((ct>0) ? -1 : lpar) * fl * CSPHASE; res.pi [qpms_mn2y(-1, l)] = ((ct>0) ? -1 : lpar) * fl * CSPHASE;
res.tau[qpms_mn2y(+1, l)] = ((ct>0) ? +1 : lpar) * fl; 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) { 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))); 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; 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;
res.pi [qpms_mn2y(-1, l)] = ((ct>0) ? -1 : lpar) * fl * CSPHASE; res.pi [qpms_mn2y(-1, l)] = ((ct>0) ? -1 : lpar) * fl * CSPHASE;
res.tau[qpms_mn2y(+1, l)] = ((ct>0) ? +1 : lpar) * fl; res.tau[qpms_mn2y(+1, l)] = ((ct>0) ? +1 : lpar) * fl;