diff --git a/qpms/legendre.c b/qpms/legendre.c index 27f79cf..9338d87 100644 --- a/qpms/legendre.c +++ b/qpms/legendre.c @@ -79,7 +79,8 @@ qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, qpms_normalisation_t no memset(res.tau, 0, nelem*sizeof(double)); res.leg = calloc(nelem, sizeof(double)); switch(norm) { - case QPMS_NORMALISATION_XU: + /* FIXME get rid of multiplicating the five lines */ + case QPMS_NORMALISATION_NONE: 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; @@ -91,11 +92,22 @@ qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, qpms_normalisation_t no res.tau[qpms_mn2y(-1, l)] = -((ct>0) ? +1 : lpar) * n * csphase; } break; + case QPMS_NORMALISATION_XU: // Rather useless except for testing. + for (qpms_l_t l = 1; l <= lMax; ++l) { + res.leg[qpms_mn2y(0, l)] = ((l%2)?ct:1.)*sqrt(0.25*M_1_PI *(2*l+1)/(l*(l+1))); + double p = ..... HERE TODO ENDED CO + 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); + int lpar = (l%2)?-1:1; 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; @@ -105,13 +117,12 @@ qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, qpms_normalisation_t no 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); + int lpar = (l%2)?-1:1; 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: @@ -122,11 +133,11 @@ qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, qpms_normalisation_t no 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 + norm == QPMS_NORMALISATION_NONE ? GSL_SF_LEGENDRE_NONE : GSL_SF_LEGENDRE_SPHARM, csphase)) abort(); if (norm == QPMS_NORMALISATION_POWER) - /* for Xu (=non-normalized) and Taylor (=sph. harm. normalized) + /* for None (=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. @@ -138,6 +149,18 @@ qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, qpms_normalisation_t no legder[qpms_mn2y(m,l)] *= prefac; } } + else if (norm == QPMS_NORMALISATION_XU) + /* for Xu (anti-normalized), we start from spharm-normalized Legendre functions + * Do not use this normalisation except for testing + */ + for (qpms_l_t l = 1; l <= lMax; ++l) { + double prefac = (2*l + 1) / sqrt(4*M_PI / (l*(l+1))); + for (qpms_m_t m = -l; m <= l; ++m) { + double fac = prefac * exp(lgamma(l+m+1) - lgamma(l-m+1)); + res.leg[qpms_mn2y(m,l)] *= fac; + legder[qpms_mn2y(m,l)] *= fac; + } + } 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)]; diff --git a/qpms/qpms_types.h b/qpms/qpms_types.h index 707aee9..3878071 100644 --- a/qpms/qpms_types.h +++ b/qpms/qpms_types.h @@ -27,7 +27,7 @@ typedef enum { #define QPMS_NORMALISATION_T_CSBIT 128 typedef enum { // As in TODO - QPMS_NORMALISATION_XU = 4, // such that the numerical values in Xu's tables match + QPMS_NORMALISATION_XU = 4, // such that the numerical values in Xu's tables match, not recommended to use otherwise QPMS_NORMALISATION_NONE = 3, // genuine unnormalised waves (with unnormalised Legendre polynomials) // As in http://www.eit.lth.se/fileadmin/eit/courses/eit080f/Literature/book.pdf, power-normalised QPMS_NORMALISATION_KRISTENSSON = 2, @@ -58,6 +58,13 @@ static inline int qpms_normalisation_t_normonly(qpms_normalisation_t norm) { // TODO move the inlines elsewhere +/* Normalisation of the spherical waves is now scattered in at least three different files: + * here, we have the norm in terms of radiated power of outgoing wave. + * In file legendre.c, function qpms_pitau_get determines the norm used in the vswf.c + * spherical vector wave norms. The "dual" waves in vswf.c use the ..._abssquare function below. + * In file translations.c, the normalisations are again set by hand using the normfac and lognormfac + * functions. + */ #include #include // relative to QPMS_NORMALISATION_KRISTENSSON_CS, i.e. @@ -76,6 +83,9 @@ static inline double qpms_normalisation_t_factor(qpms_normalisation_t norm, qpms case QPMS_NORMALISATION_NONE: factor = sqrt(l*(l+1) * 4 * M_PI / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1))); break; + case QPMS_NORMALISATION_XU: + factor = sqrt(4 * M_PI) / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1)); + break; default: assert(0); } @@ -97,6 +107,10 @@ static inline double qpms_normalisation_t_factor_abssquare(qpms_normalisation_t case QPMS_NORMALISATION_NONE: return l*(l+1) * 4 * M_PI / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1)); break; + case QPMS_NORMALISATION_XU: + double fac = sqrt(4 * M_PI) / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1)); + return fac * fac; + break; default: assert(0); return NAN; diff --git a/qpms/translations.c b/qpms/translations.c index b7263f2..5bb5488 100644 --- a/qpms/translations.c +++ b/qpms/translations.c @@ -140,11 +140,9 @@ static inline double qpms_trans_normlogfac(qpms_normalisation_t norm, case QPMS_NORMALISATION_NONE: return -(lgamma(n+m+1)-lgamma(n-m+1)+lgamma(nu-mu+1)-lgamma(nu+mu+1)); break; - /* // TODO NONE and XU are going to be different after all... case QPMS_NORMALISATION_XU: return 0; break; - */ default: abort(); } @@ -166,10 +164,8 @@ static inline double qpms_trans_normfac(qpms_normalisation_t norm, case QPMS_NORMALISATION_NONE: normfac *= (2.*n+1)/(2.*nu+1); break; - /* // TODO NONE and XU are going to be different after all... case QPMS_NORMALISATION_XU: break; - */ default: abort(); } @@ -791,6 +787,7 @@ int qpms_trans_calculator_multipliers_A(qpms_normalisation_t norm, complex doubl break; #endif case QPMS_NORMALISATION_NONE: + case QPMS_NORMALISATION_XU: case QPMS_NORMALISATION_KRISTENSSON: qpms_trans_calculator_multipliers_A_general(norm, dest, m, n, mu, nu, qmax); return 0; @@ -810,6 +807,7 @@ int qpms_trans_calculator_multipliers_B(qpms_normalisation_t norm, complex doubl break; #endif case QPMS_NORMALISATION_NONE: + case QPMS_NORMALISATION_XU: case QPMS_NORMALISATION_KRISTENSSON: qpms_trans_calculator_multipliers_B_general(norm, dest, m, n, mu, nu, Qmax); return 0; @@ -919,6 +917,7 @@ complex double qpms_trans_calculator_get_A_buf(const qpms_trans_calculator *c, case QPMS_NORMALISATION_TAYLOR: case QPMS_NORMALISATION_KRISTENSSON: case QPMS_NORMALISATION_NONE: + case QPMS_NORMALISATION_XU: { double costheta = cos(kdlj.theta); if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu, @@ -968,6 +967,7 @@ complex double qpms_trans_calculator_get_B_buf(const qpms_trans_calculator *c, case QPMS_NORMALISATION_TAYLOR: case QPMS_NORMALISATION_KRISTENSSON: case QPMS_NORMALISATION_NONE: + case QPMS_NORMALISATION_XU: { double costheta = cos(kdlj.theta); if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1, @@ -999,6 +999,7 @@ int qpms_trans_calculator_get_AB_buf_p(const qpms_trans_calculator *c, case QPMS_NORMALISATION_TAYLOR: case QPMS_NORMALISATION_KRISTENSSON: case QPMS_NORMALISATION_NONE: + case QPMS_NORMALISATION_XU: { double costheta = cos(kdlj.theta); if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,