diff --git a/qpms/bessel.c b/qpms/bessel.c index b5488b3..162a0f2 100644 --- a/qpms/bessel.c +++ b/qpms/bessel.c @@ -1,9 +1,9 @@ #include -#include "translations.h" +#include "qpms_specfunc.h" #include #include -int qpms_sph_bessel_fill(qpms_bessel_t typ, int lmax, double x, complex double *result_array) { +qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, double x, complex double *result_array) { int retval; double tmparr[lmax+1]; switch(typ) { diff --git a/qpms/legendre.c b/qpms/legendre.c new file mode 100644 index 0000000..71e2ae5 --- /dev/null +++ b/qpms/legendre.c @@ -0,0 +1,158 @@ +#include "qpms_specfunc.h" +#include "qpms_types.h" +#include +#include +#include +#include "indexing.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, + gsl_sf_legendre_t lnorm, double csphase) +{ + size_t n = gsl_sf_legendre_array_n(lMax); + double *legendre_tmp = malloc(n * sizeof(double)); + double *legendre_deriv_tmp = malloc(n * sizeof(double)); + int gsl_errno = gsl_sf_legendre_deriv_array_e( + lnorm, (size_t)lMax, x, csphase, legendre_tmp,legendre_deriv_tmp); + for (qpms_l_t l = 1; l <= lMax; ++l) + for (qpms_m_t m = 0; m <= l; ++m) { + qpms_y_t y = qpms_mn2y(m,l); + size_t i = gsl_sf_legendre_array_index(l,m); + target[y] = legendre_tmp[i]; + target_deriv[y] = legendre_deriv_tmp[i]; + } + switch(lnorm) { + case GSL_SF_LEGENDRE_NONE: + for (qpms_l_t l = 1; l <= lMax; ++l) + for (qpms_m_t m = 1; m <= l; ++m) { + 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(l+m+1))*((m%2)?-1:1); + target[y] = factor * legendre_tmp[i]; + 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 = 1; l <= lMax; ++l) + for (qpms_m_t m = 1; m <= l; ++m) { + 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 = ((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; + } + free(legendre_tmp); + free(legendre_deriv_tmp); + return QPMS_SUCCESS; +} + +qpms_errno_t qpms_legendre_deriv_y_get(double **target, double **dtarget, double x, qpms_l_t lMax, gsl_sf_legendre_t lnorm, + double csphase) +{ + + *target = malloc(sizeof(double)*qpms_lMax2nelem(lMax)); + *dtarget = malloc(sizeof(double)*qpms_lMax2nelem(lMax)); + return qpms_legendre_deriv_y_fill(*target, *dtarget, x, lMax, lnorm, csphase); +} + +qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, qpms_normalisation_t norm) +{ + const double csphase = qpms_normalisation_t_csphase(norm); + norm = qpms_normalisation_t_normonly(norm); + qpms_pitau_t res; + qpms_y_t nelem = qpms_lMax2nelem(lMax); + res.pi = malloc(nelem * sizeof(double)); + res.tau = malloc(nelem * sizeof(double)); + double ct = cos(theta), st = sin(theta); + if (1 == fabs(ct)) { // singular case, use DLMF 14.8.2 + memset(res.pi, 0, nelem*sizeof(double)); + memset(res.tau, 0, nelem*sizeof(double)); + res.leg = calloc(nelem, sizeof(double)); + switch(norm) { + 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)/2; + const double n = 0.5; + int lpar = (l%2)?-1:1; + res.pi [qpms_mn2y(+1, l)] = -((ct>0) ? -1 : lpar) * p * csphase; + res.pi [qpms_mn2y(-1, l)] = -((ct>0) ? -1 : lpar) * n * csphase; + res.tau[qpms_mn2y(+1, l)] = ((ct>0) ? +1 : lpar) * p * csphase; + res.tau[qpms_mn2y(-1, l)] = -((ct>0) ? +1 : lpar) * n * csphase; + } + break; + case QPMS_NORMALISATION_TAYLOR: + 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 * sqrt((2*l+1)*l*(l+1)*M_1_PI); + 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 * csphase; + res.tau[qpms_mn2y(-1, l)] = -((ct>0) ? +1 : lpar) * fl * csphase; + } + break; + case QPMS_NORMALISATION_POWER: + 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 * sqrt((2*l+1)*M_1_PI); + 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 * csphase; + res.tau[qpms_mn2y(-1, l)] = -((ct>0) ? +1 : lpar) * fl * csphase; + + } + break; + default: + abort(); + } + } + else { // cos(theta) in (-1,1), use normal calculation + double *legder = malloc(sizeof(double)*qpms_lMax2nelem(lMax)); + res.leg = malloc(sizeof(double)*qpms_lMax2nelem(lMax)); + if (qpms_legendre_deriv_y_fill(res.leg, legder, ct, lMax, + norm == QPMS_NORMALISATION_XU ? GSL_SF_LEGENDRE_NONE + : GSL_SF_LEGENDRE_SPHARM, csphase)) + abort(); + if (norm == QPMS_NORMALISATION_POWER) + /* for Xu (=non-normalized) and Taylor (=sph. harm. normalized) + * the correct normalisation is already obtained from gsl_sf_legendre_deriv_array_e(). + * However, Kristensson ("power") normalisation differs from Taylor + * by 1/sqrt(l*(l+1)) factor. + */ + for (qpms_l_t l = 1; l <= lMax; ++l) { + double prefac = 1./sqrt(l*(l+1)); + for (qpms_m_t m = -l; m <= l; ++m) { + res.leg[qpms_mn2y(m,l)] *= prefac; + legder[qpms_mn2y(m,l)] *= prefac; + } + } + for (qpms_l_t l = 1; l <= lMax; ++l) { + for (qpms_m_t m = -l; m <= l; ++m) { + res.pi [qpms_mn2y(m,l)] = m / st * res.leg[qpms_mn2y(m,l)]; + res.tau[qpms_mn2y(m,l)] = - st * legder[qpms_mn2y(m,l)]; + } + } + free(legder); + } + res.lMax = lMax; + return res; +} + +void qpms_pitau_free(qpms_pitau_t x) { + free(x.leg); + free(x.pi); + free(x.tau); +} + diff --git a/qpms/qpms_specfunc.h b/qpms/qpms_specfunc.h new file mode 100644 index 0000000..1f50262 --- /dev/null +++ b/qpms/qpms_specfunc.h @@ -0,0 +1,54 @@ +#ifndef QPMS_SPECFUNC_H +#define QPMS_SPECFUNC_H +#include "qpms_types.h" +#include + +/****************************************************************************** + * Spherical Bessel functions * + ******************************************************************************/ + +// TODO unify types +qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, double x, complex double *result_array); + + +/****************************************************************************** + * Legendre functions and their "angular derivatives" * + ******************************************************************************/ + +/* + * 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 + */ + +qpms_errno_t qpms_legendre_deriv_y_get(double **result, double **result_deriv, double x, qpms_l_t lMax, + gsl_sf_legendre_t lnorm, double csphase); // free() result and result_deriv yourself! +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); + + +double *qpms_legendre_y_get(double x, qpms_l_t lMax, qpms_normalisation_t norm);//NI +double *qpms_legendre0d_y_get(qpms_l_t lMax, qpms_normalisation_t norm); //NI +double *qpms_legendre_plus1d_y_get(qpms_l_t lMax, qpms_normalisation_t norm); //NI +double *qpms_legendre_minus1d_y_get(qpms_l_t lMax, qpms_normalisation_t norm); //NI + + + +// array of Legendre and pi, tau auxillary functions (see [1,(37)]) +// This should handle correct evaluation for theta -> 0 and theta -> pi +typedef struct { + //qpms_normalisation_t norm; + qpms_l_t lMax; + //qpms_y_t nelem; + double *leg, *pi, *tau; +} qpms_pitau_t; +qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, qpms_normalisation_t norm); +void qpms_pitau_free(qpms_pitau_t);//NI +void qpms_pitau_pfree(qpms_pitau_t*);//NI + +// Associated Legendre polynomial at zero argument (DLMF 14.5.1) DEPRECATED? +double qpms_legendre0(int m, int n); +// Associated Legendre polynomial derivative at zero argument (DLMF 14.5.2) +double qpms_legendred0(int m, int n); + +#endif // QPMS_SPECFUNC_H diff --git a/qpms/qpms_types.h b/qpms/qpms_types.h index e4b4712..921e80a 100644 --- a/qpms/qpms_types.h +++ b/qpms/qpms_types.h @@ -57,7 +57,9 @@ static inline int qpms_normalisation_t_normonly(qpms_normalisation_t norm) { } -// TODO move elsewhere +// TODO move the inlines elsewhere +#include +#include // relative to QPMS_NORMALISATION_KRISTENSSON_CS, i.e. // P_l^m[normtype] = P_l^m[Kristensson] static inline double qpms_normalisation_t_factor(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) { diff --git a/qpms/translations.h b/qpms/translations.h index 6e69e85..4904c61 100644 --- a/qpms/translations.h +++ b/qpms/translations.h @@ -5,65 +5,71 @@ #include #include #include -// Associated Legendre polynomial at zero argument (DLMF 14.5.1) -double qpms_legendre0(int m, int n); -// Associated Legendre polynomial derivative at zero argument (DLMF 14.5.2) -double qpms_legendred0(int m, int n); -// TODO unify types -int qpms_sph_bessel_fill(qpms_bessel_t typ, int lmax, double x, complex double *result_array); // TODO replace the xplicit "Taylor" functions with general, // taking qpms_bessel_t argument. -complex double qpms_trans_single_A_Taylor(int m, int n, int mu, int nu, sph_t kdlj, +complex double qpms_trans_single_A_Taylor(qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t 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, +complex double qpms_trans_single_B_Taylor(qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj, bool r_ge_d, qpms_bessel_t J); -complex double qpms_trans_single_A_Taylor_ext(int m, int n, int mu, int nu, double kdlj_r, +complex double qpms_trans_single_A_Taylor_ext(qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, double kdlj_r, double kdlj_th, double kdlj_phi, int r_ge_d, int J); -complex double qpms_trans_single_B_Taylor_ext(int m, int n, int mu, int nu, double kdlj_r, +complex double qpms_trans_single_B_Taylor_ext(qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, double kdlj_r, double kdlj_th, double kdlj_phi, int r_ge_d, int J); typedef struct qpms_trans_calculator { - int lMax; - size_t nelem; + qpms_normalisation_t normalisation; + qpms_l_t lMax; + qpms_y_t nelem; complex double **A_multipliers; complex double **B_multipliers; - qpms_normalisation_t normalisation; +#if 0 + // Normalised values of the Legendre functions and derivatives + // for θ == π/2, i.e. for the 2D case. + double *leg0; + double *pi0; + double *tau0; + // Spherical Bessel function coefficients: + // TODO +#endif } qpms_trans_calculator; -qpms_trans_calculator *qpms_trans_calculator_init(int lMax, qpms_normalisation_t nt); +qpms_trans_calculator *qpms_trans_calculator_init(qpms_l_t lMax, qpms_normalisation_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, + qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t 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, + qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj, bool r_ge_d, qpms_bessel_t J); +int qpms_trans_calculator_get_AB_p(const qpms_trans_calculator *c, + complex double *Adest, complex double *Bdest, + qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj, + bool r_ge_d, qpms_bessel_t J); +int qpms_trans_calculator_get_AB_arrays(const qpms_trans_calculator *c, + complex double *Adest, complex double *Bdest, + size_t deststride, size_t srcstride, + sph_t kdlj, bool r_ge_d, qpms_bessel_t J); + +// TODO update the types later complex double qpms_trans_calculator_get_A_ext(const qpms_trans_calculator *c, int m, int n, int mu, int nu, double kdlj_r, double kdlj_th, double kdlj_phi, int r_ge_d, int J); + complex double qpms_trans_calculator_get_B_ext(const qpms_trans_calculator *c, int m, int n, int mu, int nu, double kdlj_r, double kdlj_th, double kdlj_phi, int r_ge_d, int J); -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); int qpms_trans_calculator_get_AB_p_ext(const qpms_trans_calculator *c, complex double *Adest, complex double *Bdest, int m, int n, int mu, int nu, double kdlj_r, double kdlj_th, double kdlj_phi, int r_ge_d, int J); -int qpms_trans_calculator_get_AB_arrays(const qpms_trans_calculator *c, - complex double *Adest, complex double *Bdest, - size_t deststride, size_t srcstride, - sph_t kdlj, bool r_ge_d, qpms_bessel_t J); int qpms_trans_calculator_get_AB_arrays_ext(const qpms_trans_calculator *c, complex double *Adest, complex double *Bdest, size_t deststride, size_t srcstride, diff --git a/qpms/vswf.c b/qpms/vswf.c index c91cb8c..aa03ab9 100644 --- a/qpms/vswf.c +++ b/qpms/vswf.c @@ -4,162 +4,10 @@ #include "vswf.h" #include "indexing.h" #include "translations.h" // TODO move qpms_sph_bessel_fill elsewhere -#include +#include "qpms_specfunc.h" #include #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, - gsl_sf_legendre_t lnorm, double csphase) -{ - size_t n = gsl_sf_legendre_array_n(lMax); - double *legendre_tmp = malloc(n * sizeof(double)); - double *legendre_deriv_tmp = malloc(n * sizeof(double)); - int gsl_errno = gsl_sf_legendre_deriv_array_e( - lnorm, (size_t)lMax, x, csphase, legendre_tmp,legendre_deriv_tmp); - for (qpms_l_t l = 1; l <= lMax; ++l) - for (qpms_m_t m = 0; m <= l; ++m) { - qpms_y_t y = qpms_mn2y(m,l); - size_t i = gsl_sf_legendre_array_index(l,m); - target[y] = legendre_tmp[i]; - target_deriv[y] = legendre_deriv_tmp[i]; - } - switch(lnorm) { - case GSL_SF_LEGENDRE_NONE: - for (qpms_l_t l = 1; l <= lMax; ++l) - for (qpms_m_t m = 1; m <= l; ++m) { - 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(l+m+1))*((m%2)?-1:1); - target[y] = factor * legendre_tmp[i]; - 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 = 1; l <= lMax; ++l) - for (qpms_m_t m = 1; m <= l; ++m) { - 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 = ((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; - } - free(legendre_tmp); - free(legendre_deriv_tmp); - return QPMS_SUCCESS; -} - -qpms_errno_t qpms_legendre_deriv_y_get(double **target, double **dtarget, double x, qpms_l_t lMax, gsl_sf_legendre_t lnorm, - double csphase) -{ - - *target = malloc(sizeof(double)*qpms_lMax2nelem(lMax)); - *dtarget = malloc(sizeof(double)*qpms_lMax2nelem(lMax)); - return qpms_legendre_deriv_y_fill(*target, *dtarget, x, lMax, lnorm, csphase); -} - - -qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, qpms_normalisation_t norm) -{ - const double csphase = qpms_normalisation_t_csphase(norm); - norm = qpms_normalisation_t_normonly(norm); - qpms_pitau_t res; - qpms_y_t nelem = qpms_lMax2nelem(lMax); - res.pi = malloc(nelem * sizeof(double)); - res.tau = malloc(nelem * sizeof(double)); - double ct = cos(theta), st = sin(theta); - if (1 == fabs(ct)) { // singular case, use DLMF 14.8.2 - memset(res.pi, 0, nelem*sizeof(double)); - memset(res.tau, 0, nelem*sizeof(double)); - res.leg = calloc(nelem, sizeof(double)); - switch(norm) { - 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)/2; - const double n = 0.5; - int lpar = (l%2)?-1:1; - res.pi [qpms_mn2y(+1, l)] = -((ct>0) ? -1 : lpar) * p * csphase; - res.pi [qpms_mn2y(-1, l)] = -((ct>0) ? -1 : lpar) * n * csphase; - res.tau[qpms_mn2y(+1, l)] = ((ct>0) ? +1 : lpar) * p * csphase; - res.tau[qpms_mn2y(-1, l)] = -((ct>0) ? +1 : lpar) * n * csphase; - } - break; - case QPMS_NORMALISATION_TAYLOR: - 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 * sqrt((2*l+1)*l*(l+1)*M_1_PI); - 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 * csphase; - res.tau[qpms_mn2y(-1, l)] = -((ct>0) ? +1 : lpar) * fl * csphase; - } - break; - case QPMS_NORMALISATION_POWER: - 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 * sqrt((2*l+1)*M_1_PI); - 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 * csphase; - res.tau[qpms_mn2y(-1, l)] = -((ct>0) ? +1 : lpar) * fl * csphase; - - } - break; - default: - abort(); - } - } - else { // cos(theta) in (-1,1), use normal calculation - double *legder = malloc(sizeof(double)*qpms_lMax2nelem(lMax)); - res.leg = malloc(sizeof(double)*qpms_lMax2nelem(lMax)); - if (qpms_legendre_deriv_y_fill(res.leg, legder, ct, lMax, - norm == QPMS_NORMALISATION_XU ? GSL_SF_LEGENDRE_NONE - : GSL_SF_LEGENDRE_SPHARM, csphase)) - abort(); - if (norm == QPMS_NORMALISATION_POWER) - /* for Xu (=non-normalized) and Taylor (=sph. harm. normalized) - * the correct normalisation is already obtained from gsl_sf_legendre_deriv_array_e(). - * However, Kristensson ("power") normalisation differs from Taylor - * by 1/sqrt(l*(l+1)) factor. - */ - for (qpms_l_t l = 1; l <= lMax; ++l) { - double prefac = 1./sqrt(l*(l+1)); - for (qpms_m_t m = -l; m <= l; ++m) { - res.leg[qpms_mn2y(m,l)] *= prefac; - legder[qpms_mn2y(m,l)] *= prefac; - } - } - for (qpms_l_t l = 1; l <= lMax; ++l) { - for (qpms_m_t m = -l; m <= l; ++m) { - res.pi [qpms_mn2y(m,l)] = m / st * res.leg[qpms_mn2y(m,l)]; - res.tau[qpms_mn2y(m,l)] = - st * legder[qpms_mn2y(m,l)]; - } - } - free(legder); - } - res.lMax = lMax; - return res; -} - -void qpms_pitau_free(qpms_pitau_t x) { - free(x.leg); - free(x.pi); - free(x.tau); -} - - csphvec_t qpms_vswf_single_el(qpms_m_t m, qpms_l_t l, sph_t kdlj, qpms_bessel_t btyp, qpms_normalisation_t norm) { lmcheck(l,m); diff --git a/qpms/vswf.h b/qpms/vswf.h index eda78c2..778a3d0 100644 --- a/qpms/vswf.h +++ b/qpms/vswf.h @@ -60,23 +60,4 @@ qpms_vswfset_sph_t *qpms_vswfset_make(qpms_l_t lMax, sph_t kdlj, qpms_bessel_t btyp, qpms_normalisation_t norm);//NI void qpms_vswfset_sph_pfree(qpms_vswfset_sph_t *);//NI -double *qpms_legendre_y_get(double x, qpms_l_t lMax, qpms_normalisation_t norm);//NI -double *qpms_legendre0d_y_get(qpms_l_t lMax, qpms_normalisation_t norm); //NI -double *qpms_legendre_plus1d_y_get(qpms_l_t lMax, qpms_normalisation_t norm); //NI -double *qpms_legendre_minus1d_y_get(qpms_l_t lMax, qpms_normalisation_t norm); //NI - - - -// array of Legendre and pi, tau auxillary functions (see [1,(37)]) -// This should handle correct evaluation for theta -> 0 and theta -> pi -typedef struct { - //qpms_normalisation_t norm; - qpms_l_t lMax; - //qpms_y_t nelem; - double *leg, *pi, *tau; -} qpms_pitau_t; -qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, qpms_normalisation_t norm); -void qpms_pitau_free(qpms_pitau_t);//NI -void qpms_pitau_pfree(qpms_pitau_t*);//NI - #endif // QPMS_VSWF_H