[IN PRO6RESS] Xu's antinormalisation
Former-commit-id: 577e58118d7da93b927035b4e6e91c6c4cdab485
This commit is contained in:
Marek Nečada
parent
45edb30a62
commit
8f6181ced8
|
|
@ -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)];
|
||||
|
|
|
|||
|
|
@ -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 <math.h>
|
||||
#include <assert.h>
|
||||
// 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;
|
||||
|
|
|
|||
|
|
@ -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,
|
||||
|
|
|
|||
Loading…
Reference in New Issue