Jdu do práce
Former-commit-id: 3be0f5adda5334dae7397d21b6d48835a190632d
This commit is contained in:
Marek Nečada
parent
2a7015b80b
commit
00f9f5234a
|
|
@ -39,6 +39,9 @@ typedef enum {
|
|||
QPMS_NORMALISATION_UNDEF = 0
|
||||
} qpms_normalisation_t;
|
||||
|
||||
|
||||
|
||||
|
||||
static inline int qpms_normalisation_t_csphase(qpms_normalisation_t norm) {
|
||||
return (norm & QPMS_NORMALISATION_T_CSBIT)? -1 : 1;
|
||||
}
|
||||
|
|
@ -48,6 +51,46 @@ static inline int qpms_normalisation_t_normonly(qpms_normalisation_t norm) {
|
|||
}
|
||||
|
||||
|
||||
// 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) {
|
||||
int csphase = qpms_normalisation_t_csphase(norm);
|
||||
norm = qpms_normalisation_t_normonly(norm);
|
||||
double factor;
|
||||
switch (norm) {
|
||||
case QPMS_NORMALISATION_KRISTENSSON:
|
||||
factor = 1.;
|
||||
break;
|
||||
case QPMS_NORMALISATION_TAYLOR:
|
||||
factor = sqrt(l*(l+1));
|
||||
break;
|
||||
case QPMS_NORMALISATION_XU:
|
||||
factor = sqrt(l*(l+1) * 4 * M_PI / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1)));
|
||||
break;
|
||||
default:
|
||||
assert(0);
|
||||
}
|
||||
factor *= (m%2)?(-csphase):1;
|
||||
return factor;
|
||||
}
|
||||
|
||||
static inline double qpms_normalisation_t_factor_abssquare(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
|
||||
norm = qpms_normalisation_t_normonly(norm);
|
||||
switch (norm) {
|
||||
case QPMS_NORMALISATION_KRISTENSSON:
|
||||
return 1.;
|
||||
break;
|
||||
case QPMS_NORMALISATION_TAYLOR:
|
||||
return l*(l+1);
|
||||
break;
|
||||
case QPMS_NORMALISATION_XU:
|
||||
return l*(l+1) * 4 * M_PI / (2*l+1) * exp(lgamma(l+m+1)-lgamma(l-m+1));
|
||||
break;
|
||||
default:
|
||||
assert(0);
|
||||
return NAN;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
typedef enum {
|
||||
|
|
@ -75,6 +118,11 @@ typedef struct {
|
|||
double r, theta, phi;
|
||||
} sph_t;
|
||||
|
||||
typedef struct { // Do I really need this???
|
||||
complex double r;
|
||||
double theta, phi;
|
||||
} csph_t;
|
||||
|
||||
// complex vector components in local spherical basis
|
||||
typedef struct {
|
||||
complex double rc, thetac, phic;
|
||||
|
|
|
|||
|
|
@ -27,6 +27,21 @@ static inline cart3_t sph2cart(const sph_t sph) {
|
|||
return cart;
|
||||
}
|
||||
|
||||
static inline cart3_t cart3_add(const cart3_t a, const cart3_t b) {
|
||||
cart3_t res = {a.x+b.x, a.y+b.y, a.z+b.z};
|
||||
return res;
|
||||
}
|
||||
|
||||
static inline csphvec_t csphvec_add(const csphvec_t a, const csphvec_t b) {
|
||||
csphvec_t res = {a.rc + b.rc, a.thetac + b.thetac, a.phic + b.phic};
|
||||
return res;
|
||||
}
|
||||
|
||||
static inline csphvec_t scphvec_scale(complex double c, const scphvec_t v) {
|
||||
csphvec_t res = {c * v.rc, c * v.thetac, c * v.phic};
|
||||
return res;
|
||||
}
|
||||
|
||||
// equivalent to sph_loccart2cart in qpms_p.py
|
||||
static inline ccart3_t csphvec2cart(const csphvec_t sphvec, const sph_t at) {
|
||||
ccart3_t res = {0, 0, 0};
|
||||
|
|
|
|||
176
qpms/vswf.c
176
qpms/vswf.c
|
|
@ -217,7 +217,7 @@ void qpms_vswfset_sph_pfree(qpms_vswfset_sph_t *w) {
|
|||
free(w);
|
||||
}
|
||||
|
||||
qpms_errno_t qpms_vswf_fill(csphvec_t *mgtarget, csphvec_t *eltarget,
|
||||
qpms_errno_t qpms_vswf_fill(csphvec_t *const longtarget, csphvec_t * const mgtarget, csphvec_t * const eltarget,
|
||||
qpms_l_t lMax, sph_t kr,
|
||||
qpms_bessel_t btyp, qpms_normalisation_t norm) {
|
||||
assert(lMax >= 1);
|
||||
|
|
@ -228,19 +228,32 @@ qpms_errno_t qpms_vswf_fill(csphvec_t *mgtarget, csphvec_t *eltarget,
|
|||
double const *pleg = pt.leg;
|
||||
double const *ppi = pt.pi;
|
||||
double const *ptau = pt.tau;
|
||||
csphvec_t *pmg = mgtarget, *pel = eltarget;
|
||||
csphvec_t *plong = longtarget, *pmg = mgtarget, *pel = eltarget;
|
||||
for(qpms_l_t l = 1; l <= lMax; ++l) {
|
||||
complex double besfac = *pbes / kr.r;
|
||||
complex double besderfac = *(pbes-1) - l * besfac;
|
||||
for(qpms_m_t m = -l; m <= l; ++m) {
|
||||
complex double eimf = cexp(m * kr.phi * I);
|
||||
pel->rc = l*(l+1) * (*pleg) * besfac * eimf;
|
||||
pel->thetac = *ptau * besderfac * eimf;
|
||||
pel->phic = *ppi * besderfac * I * eimf;
|
||||
pmg->rc = 0.;
|
||||
pmg->thetac = *ppi * (*pbes) * I * eimf;
|
||||
pmg->phic = - *ptau * (*pbes) * eimf;
|
||||
++pleg; ++ppi; ++ptau; ++pel; ++pmg;
|
||||
if (longtarget) {
|
||||
complex double longfac = sqrt(l*(l+1)) * eimf;
|
||||
plong->rc = (besderfac-besfac) * (*pleg) * longfac;
|
||||
plong->thetac = *ptau * besfac * longfac;
|
||||
plong->phic = *ppi * I * besfac * longfac;
|
||||
++plong;
|
||||
}
|
||||
if (eltarget) {
|
||||
pel->rc = l*(l+1) * (*pleg) * besfac * eimf;
|
||||
pel->thetac = *ptau * besderfac * eimf;
|
||||
pel->phic = *ppi * besderfac * I * eimf;
|
||||
++pel;
|
||||
}
|
||||
if (mgtarget) {
|
||||
pmg->rc = 0.;
|
||||
pmg->thetac = *ppi * (*pbes) * I * eimf;
|
||||
pmg->phic = - *ptau * (*pbes) * eimf;
|
||||
++pmg;
|
||||
}
|
||||
++pleg; ++ppi; ++ptau;
|
||||
}
|
||||
++pbes;
|
||||
}
|
||||
|
|
@ -248,3 +261,148 @@ qpms_errno_t qpms_vswf_fill(csphvec_t *mgtarget, csphvec_t *eltarget,
|
|||
qpms_pitau_free(pt);
|
||||
return QPMS_SUCCESS;
|
||||
}
|
||||
|
||||
// consistency check: this should give the same results as the above function (up to rounding errors)
|
||||
qpms_errno_t qpms_vswf_fill_alternative(csphvec_t *const longtarget, csphvec_t * const mgtarget, csphvec_t * const eltarget,
|
||||
qpms_l_t lMax, sph_t kr,
|
||||
qpms_bessel_t btyp, qpms_normalisation_t norm) {
|
||||
assert(lMax >= 1);
|
||||
complex double *bessel = malloc((lMax+1)*sizeof(complex double));
|
||||
if(qpms_sph_bessel_fill(btyp, lMax, kr.r, bessel)) abort();
|
||||
complex double const *pbes = bessel + 1; // starting from l = 1
|
||||
|
||||
qpms_y_t nelem = qpms_lMax2nelem(lMax);
|
||||
csphvec_t * const a1 = malloc(3*nelem*sizeof(csphvec_t)), * const a2 = a1 + nelem, * const a3 = a2 + nelem;
|
||||
if(qpms_vecspharm_fill(a1, a2, a3, lMax, kr, norm)) abort();
|
||||
const sphvec_t *p1 = a1, const *p2 = a2, const *p3 = a3;
|
||||
|
||||
csphvec_t *plong = longtarget, *pmg = mgtarget, *pel = eltarget;
|
||||
for(qpms_l_t l = 1; l <= lMax; ++l) {
|
||||
complex double besfac = *pbes / kr.r;
|
||||
complex double besderfac = *(pbes-1) - l * besfac;
|
||||
double sqrtlfac = sqrt(l*(l+1));
|
||||
for(qpms_m_t m = -l; m <= l; ++m) {
|
||||
complex double eimf = cexp(m * kr.phi * I);
|
||||
if (longtarget) {
|
||||
*plong = csphvec_add(csphvec_scale(besderfac-besfac, *p3),
|
||||
csphvec_scale(sqrtlfac * besfac, *p2));
|
||||
++plong;
|
||||
}
|
||||
if (eltarget) {
|
||||
*pel = csphvec_add(csphvec_scale(besderfac, *p2),
|
||||
csphvec_scale(sqrtlfac * besfac, *p3));
|
||||
++pel;
|
||||
}
|
||||
if (mgtarget) {
|
||||
*pmg = csphvec_scale(*pbes, *p1);
|
||||
++pmg;
|
||||
}
|
||||
++p1; ++p2; ++p3;
|
||||
}
|
||||
++pbes;
|
||||
}
|
||||
free(a1);
|
||||
free(bessel);
|
||||
return QPMS_SUCCESS;
|
||||
}
|
||||
|
||||
qpms_errno_t qpms_vecspharm_fill(csphvec_t *const a1target, csphvec_t *const a2target, csphvec_t *const a3target,
|
||||
qpms_l_t lMax, sph_t dir, qpms_normalisation_t norm) {
|
||||
assert(lMax >= 1);
|
||||
qpms_pitau_t pt = qpms_pitau_get(dir.theta, lMax, norm);
|
||||
double const *pleg = pt.leg;
|
||||
double const *ppi = pt.pi;
|
||||
double const *ptau = pt.tau;
|
||||
csphvec_t *p1 = a1target, *p2 = a2target, *p3 = a3target;
|
||||
for(qpms_l_t = 1; l <= lMax; ++l) {
|
||||
for(qpms_m_t m = -l; m <= l; ++m) {
|
||||
complex double eimf = cexp(m * dir.phi * I);
|
||||
if (a1target) {
|
||||
p1->rc = 0;
|
||||
p1->thetac = *ppi * I * eimf;
|
||||
p1->phic = -*ptau * eimf;
|
||||
++p1;
|
||||
}
|
||||
if (a2target) {
|
||||
p2->rc = 0;
|
||||
p2->thetac = *ptau * eimf;
|
||||
p2->phic = *ppi * I * eimf;
|
||||
++p2;
|
||||
}
|
||||
if (a3target) {
|
||||
p3->rc = sqrt(l*(l+1)) * (*pleg) * eimf;
|
||||
p3->thetac = 0;
|
||||
p3->phic = 0;
|
||||
++p3;
|
||||
}
|
||||
}
|
||||
++pleg; ++ppi; ++ptau;
|
||||
}
|
||||
qpms_pitau_free(pt);
|
||||
return QPMS_SUCCESS;
|
||||
}
|
||||
|
||||
qpms_errno_t qpms_vecspharm_dual_fill(csphvec_t *const a1target, csphvec_t *const a2target, csphvec_t *const a3target,
|
||||
qpms_l_t lMax, sph_t dir, qpms_normalisation_t norm) {
|
||||
assert(lMax >= 1);
|
||||
qpms_pitau_t pt = qpms_pitau_get(dir.theta, lMax, norm);
|
||||
double const *pleg = pt.leg;
|
||||
double const *ppi = pt.pi;
|
||||
double const *ptau = pt.tau;
|
||||
csphvec_t *p1 = a1target, *p2 = a2target, *p3 = a3target;
|
||||
for(qpms_l_t = 1; l <= lMax; ++l) {
|
||||
for(qpms_m_t m = -l; m <= l; ++m) {
|
||||
double normfac = 1./qpms_normalisation_t_factor_abssquare(norm, l, m); // factor w.r.t. Kristensson
|
||||
complex double eimf = cexp(m * dir.phi * I);
|
||||
if (a1target) {
|
||||
p1->rc = 0;
|
||||
p1->thetac = conj(*ppi * normfac * I * eimf);
|
||||
p1->phic = conj(-*ptau * normfac * eimf);
|
||||
++p1;
|
||||
}
|
||||
if (a2target) {
|
||||
p2->rc = 0;
|
||||
p2->thetac = conj(*ptau * normfac * eimf);
|
||||
p2->phic = conj(*ppi * normfac * I * eimf);
|
||||
++p2;
|
||||
}
|
||||
if (a3target) {
|
||||
p3->rc = conj(sqrt(l*(l+1)) * (*pleg) * normfac * eimf);
|
||||
p3->thetac = 0;
|
||||
p3->phic = 0;
|
||||
++p3;
|
||||
}
|
||||
}
|
||||
++pleg; ++ppi; ++ptau;
|
||||
}
|
||||
qpms_pitau_free(pt);
|
||||
return QPMS_SUCCESS;
|
||||
}
|
||||
|
||||
|
||||
qpms_errno_t qpms_planewave2vswf_fill_sph(sph_t wavedir, csphvec_t amplitude,
|
||||
complex double *targt_longcoeff, complex double *target_mgcoeff,
|
||||
complex double *target_elcoeff, qpms_l_t lMax, qpms_normalisation_t norm) {
|
||||
abort(); //NI
|
||||
qpms_y_t nelem = qpms_lMax2nelem(lMax);
|
||||
csphvec_t * const dual_A1 = malloc(3*nelem*sizeof(csphvec_t)), *const dual_A2 = dual_A1 + nelem,
|
||||
* const dual_A3 = dual_A2 + nelem;
|
||||
if (QPMS_SUCCESS != qpms_vecspharm_dual_fill(dual_A1, dual_A2, dual_A3, lMax, wavedir, norm))
|
||||
abort();
|
||||
for (qpms_l_t l = 1; l <= lMax; ++l) {
|
||||
...
|
||||
}
|
||||
|
||||
free(dual_A1);
|
||||
return QPMS_SUCCESS;
|
||||
}
|
||||
|
||||
qpms_errno_t qpms_planewave2vswf_fill_cart(cart3_t wavedir /*allow complex k?*/, ccart3_t amplitude,
|
||||
complex double * const longcoeff, complex double * const mgcoeff,
|
||||
complex double * const elcoeff, qpms_normalisation_t norm)
|
||||
{
|
||||
abort(); //NI
|
||||
|
||||
return QPMS_SUCCESS;
|
||||
}
|
||||
|
||||
|
|
|
|||
21
qpms/vswf.h
21
qpms/vswf.h
|
|
@ -31,8 +31,27 @@ qpms_errno_t qpms_legendre_deriv_y_get(double **result, double **result_deriv, d
|
|||
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);
|
||||
|
||||
qpms_errno_t qpms_vswf_fill(csphvec_t *resultM, csphvec_t *resultN, qpms_l_t lMax, sph_t kdrj,
|
||||
/* some of the result targets may be NULL */
|
||||
qpms_errno_t qpms_vswf_fill(csphvec_t *resultL, csphvec_t *resultM, csphvec_t *resultN, qpms_l_t lMax, sph_t kdrj,
|
||||
qpms_bessel_t btyp, qpms_normalisation_t norm);
|
||||
// Should give the same results: for consistency checks
|
||||
qpms_errno_t qpms_vswf_fill_alternative(csphvec_t *resultL, csphvec_t *resultM, csphvec_t *resultN, qpms_l_t lMax, sph_t kdrj,
|
||||
qpms_bessel_t btyp, qpms_normalisation_t norm);
|
||||
|
||||
qpms_errno_t qpms_vecspharm_fill(csphvec_t *const a1target, csphvec_t *const a2target, csphvec_t *const a3target,
|
||||
qpms_l_t lMax, sph_t dir, qpms_normalisation_t norm);
|
||||
qpms_errno_t qpms_vecspharm_dual_fill(csphvec_t *const a1target, csphvec_t *const a2target, csphvec_t *const a3target,
|
||||
qpms_l_t lMax, sph_t dir, qpms_normalisation_t norm);
|
||||
|
||||
qpms_errno_t qpms_planewave2vswf_fill_cart(ccart3_t amplitude, ccart3_t wavedir /* real or complex? */,
|
||||
complex double *targt_longcoeff, complex double *target_mgcoeff, complex double *target_elcoeff,
|
||||
qpms_normalisation_t norm);
|
||||
qpms_errno_t qpms_planewave2vswf_fill_cart(sph_t wavedir, csphvec_t amplitude,
|
||||
complex double *targt_longcoeff, complex double *target_mgcoeff, complex double *target_elcoeff,
|
||||
qpms_normalisation_t norm);
|
||||
|
||||
|
||||
|
||||
|
||||
qpms_vswfset_sph_t *qpms_vswfset_make(qpms_l_t lMax, sph_t kdlj,
|
||||
qpms_bessel_t btyp, qpms_normalisation_t norm);//NI
|
||||
|
|
|
|||
Loading…
Reference in New Issue