diff --git a/qpms/vswf.c b/qpms/vswf.c index 06bb210..5195dee 100644 --- a/qpms/vswf.c +++ b/qpms/vswf.c @@ -1,4 +1,5 @@ #include "vswf.h" +#include // Legendre functions also for negative m, see DLMF 14.9.3 qpms_errno_t qpms_legendre_deriv_y_fill(double *target, double *target_deriv, double x, qpms_l_t lMax, @@ -25,12 +26,22 @@ qpms_errno_t qpms_legendre_deriv_y_fill(double *target, double *target_deriv, do // 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); target[y] = factor * legendre_tmp[i]; - target_deriv[y] = factor * legendre_deriv_tmp; + target_deriv[y] = factor * legendre_deriv_tmp[i]; } break; case GSL_SF_LEGENDRE_SCHMIDT: case GSL_SF_LEGENDRE_SPHARM: case GSL_SF_LEGENDRE_FULL: + for (qpms_l_t l = 0; l <= lMax; ++l) + for (qpms_m_t m = 1; m <= l; ++m) { + qpms_y_t y = gpms_mn2y(-m,l); + size_t i = gsl_sf_legenre_array_index(l,m); + // viz DLMF 14.9.3, čert ví, jak je to s cs fasí. + double factor = ((m%2)?-1:1); // this is the difference from the unnormalised case + target[y] = factor * legendre_tmp[i]; + target_deriv[y] = factor * legendre_deriv_tmp[i]; + } + break; default: abort(); //NI break; diff --git a/qpms/vswf.h b/qpms/vswf.h index 3615669..8f989bc 100644 --- a/qpms/vswf.h +++ b/qpms/vswf.h @@ -20,9 +20,15 @@ typedef struct { } qpms_vswfset_sph_t; -double *qpms_legenre_deriv_y_get(double x, qpms_l_t lMax, - gsle_sf_legendre_t lnorm, double csphase); -qpms_errno_t qpms_legenre_deriv_y_fill(double *where, double x, +/* + * N.B. for the norm definitions, see + * https://www.gnu.org/software/gsl/manual/html_node/Associated-Legendre-Polynomials-and-Spherical-Harmonics.html + * ( gsl/specfunc/legendre_source.c and 7.24.2 of gsl docs + * + +double *qpms_legendre_deriv_y_get(double x, qpms_l_t lMax, + gsle_sf_legendre_t lnorm, double csphase); //FIXME what about derivative? +qpms_errno_t qpms_legendre_deriv_y_fill(double *where, double *where_deriv, double x, qpms_l_t lMax, gsl_sf_legendre_t lnorm, double csphase); //NI qpms_vswfset_sph_t *qpms_vswfset_make(qpms_l_t lMax, sph_t kdlj,