QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Rewrite vswf.c to use the new qpms_normalisation_t 

Former-commit-id: a623d0cf1c65b6134756e76e1739f189f3d6d53f
This commit is contained in:
Marek Nečada 2019-07-11 16:57:15 +03:00
parent 616095890f
commit bd40daccb7
5 changed files with 116 additions and 46 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

54
qpms/normalisation.h
View File

@ -7,6 +7,7 @@
#define NORMALISATION_H #define NORMALISATION_H
#include "qpms_types.h" #include "qpms_types.h"
#include "qpms_error.h"
#include <math.h> #include <math.h>
#include <complex.h> #include <complex.h>
@ -49,7 +50,7 @@ static inline complex double qpms_normalisation_factor_M_noCS(qpms_normalisation
* a `gsl_sf_legendre_*_e()` call. * a `gsl_sf_legendre_*_e()` call.
*/ */
static inline complex double qpms_normalisation_factor_M(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) { static inline complex double qpms_normalisation_factor_M(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
complex double fac = qn; complex double fac = qpms_normalisation_factor_M_noCS(norm, l, m);
return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac; return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac;
} }
@ -75,12 +76,11 @@ static inline complex double qpms_normalisation_factor_N_noCS(qpms_normalisation
* a `gsl_sf_legendre_*_e()` call. * a `gsl_sf_legendre_*_e()` call.
*/ */
static inline complex double qpms_normalisation_factor_N(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) { static inline complex double qpms_normalisation_factor_N(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
complex double fac = qn; complex double fac = qpms_normalisation_factor_N_noCS(norm, l, m);
return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac; return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac;
} }
#if 0
/// Returns the factors of a longitudinal VSWF of a given convention compared to the reference convention. /// Returns the factors of a longitudinal VSWF of a given convention compared to the reference convention.
/** /**
* This version ignores the Condon-Shortley phase bit (perhaps because the Condon-Shortley * This version ignores the Condon-Shortley phase bit (perhaps because the Condon-Shortley
@ -101,10 +101,9 @@ static inline complex double qpms_normalisation_factor_L_noCS(qpms_normalisation
* a `gsl_sf_legendre_*_e()` call. * a `gsl_sf_legendre_*_e()` call.
*/ */
static inline complex double qpms_normalisation_factor_L(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) { static inline complex double qpms_normalisation_factor_L(qpms_normalisation_t norm, qpms_l_t l, qpms_m_t m) {
complex double fac = qn; complex double fac = qpms_normalisation_factor_L_noCS(norm, l, m);
return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac; return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac;
} }
#endif
/// Returns normalisation flags corresponding to the dual spherical harmonics / waves. /// Returns normalisation flags corresponding to the dual spherical harmonics / waves.
/** /**
@ -122,20 +121,20 @@ static inline qpms_normalisation_t qpms_normalisation_dual(qpms_normalisation_t
/// Returns the asimuthal part of a spherical harmonic. /// Returns the asimuthal part of a spherical harmonic.
/** Returns \f[ e^{im\phi} \f] for standard complex spherical harmonics, /** Returns \f[ e^{im\phi} \f] for standard complex spherical harmonics,
* \f[ e^{-im\phi \f] for complex spherical harmonics * \f[ e^{-im\phi \f] for complex spherical harmonics
* and QPMS_NORMALISATION_REVERSE_ASIMUTHAL_PHASE set. * and QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE set.
* *
* For real spherical harmonics, this gives * For real spherical harmonics, this gives
* \f[ * \f[
* \sqrt{2}\cos{m \phi} \quad \mbox{if } m>0, \\ * \sqrt{2}\cos{m \phi} \quad \mbox{if } m>0, \\
* \sqrt{2}\sin{m \phi} \quad \mbox{if } m<0, \\ * \sqrt{2}\sin{m \phi} \quad \mbox{if } m<0, \\
* 1 \quad \mbox{if } m>0. \\ * 0 \quad \mbox{if } m>0. \\
* \f] * \f]
*/ */
static inline complex double qpms_spharm_azimuthal_part(qpms_normalisation_t norm, qpms_m_t m, double phi) { static inline complex double qpms_spharm_azimuthal_part(qpms_normalisation_t norm, qpms_m_t m, double phi) {
switch(norm & (QPMS_NORMALISATION_REVERSE_ASIMUTHAL_PHASE | QPMS_NORMALISATION_SPHARM_REAL)) { switch(norm & (QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE | QPMS_NORMALISATION_SPHARM_REAL)) {
case 0: case 0:
return cexp(I*m*phi); return cexp(I*m*phi);
case QPMS_NORMALISATION_REVERSE_ASIMUTHAL_PHASE: case QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE:
return cexp(-I*m*phi); return cexp(-I*m*phi);
case QPMS_NORMALISATION_SPHARM_REAL: case QPMS_NORMALISATION_SPHARM_REAL:
if (m > 0) return M_SQRT2 * cos(m*phi); if (m > 0) return M_SQRT2 * cos(m*phi);
@ -146,4 +145,41 @@ static inline complex double qpms_spharm_azimuthal_part(qpms_normalisation_t nor
} }
} }
/// Returns derivative of the asimuthal part of a spherical harmonic divided by \a m.
/**
*
* This is used to evaluate the VSWFs together with the \a pi member array of the
* qpms_pitau_t structure.
*
* Returns \f[ i e^{im\phi} \f] for standard complex spherical harmonics,
* \f[-i e^{-i\phi \f] for complex spherical harmonics
* and QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE set.
*
* For real spherical harmonics, this gives
* \f[
* -\sqrt{2}\sin{m \phi} \quad \mbox{if } m>0, \\
* \sqrt{2}\cos{m \phi} \quad \mbox{if } m<0, \\
* -1 \quad \mbox{if } \mbox{if }m=0. \\
* \f]
*
* (The value returned for \f$ m = 0 \f$ should not actually be used for
* anything except for multiplying by zero.)
*
*
*/
static inline complex double qpms_spharm_azimuthal_part_derivative_div_m(qpms_normalisation_t norm, qpms_m_t m, double phi) {
if(m==0) return 0;
switch(norm & (QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE | QPMS_NORMALISATION_SPHARM_REAL)) {
case 0:
return I*cexp(I*m*phi);
case QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE:
return -I*cexp(-I*m*phi);
case QPMS_NORMALISATION_SPHARM_REAL:
if (m > 0) return -M_SQRT2 * sin(m*phi);
else if (m < 0) return M_SQRT2 * cos(m*phi);
else return -1;
default:
QPMS_WTF;
}
}
#endif //NORMALISATION_H #endif //NORMALISATION_H

3
qpms/qpms_specfunc.h
View File

@ -1,3 +1,6 @@
/*! \file qpms_specfunc.h
* \brief Various special and auxillary functions.
*/
#ifndef QPMS_SPECFUNC_H #ifndef QPMS_SPECFUNC_H
#define QPMS_SPECFUNC_H #define QPMS_SPECFUNC_H
#include "qpms_types.h" #include "qpms_types.h"

2
qpms/qpms_types.h
View File

@ -119,10 +119,8 @@ typedef enum {
QPMS_NORMALISATION_M_MINUS = 32, ///< Include an additional \f$-1\f$ -factor into the magnetic waves. QPMS_NORMALISATION_M_MINUS = 32, ///< Include an additional \f$-1\f$ -factor into the magnetic waves.
QPMS_NORMALISATION_N_I = 64, ///< Include an additional \a i -factor into the electric waves. QPMS_NORMALISATION_N_I = 64, ///< Include an additional \a i -factor into the electric waves.
QPMS_NORMALISATION_N_MINUS = 128, ///< Include an additional \f$-1\f$ -factor into the magnetic waves. QPMS_NORMALISATION_N_MINUS = 128, ///< Include an additional \f$-1\f$ -factor into the magnetic waves.
#if 0
QPMS_NORMALISATION_L_I = 256, ///< Include an additional \a i -factor into the longitudinal waves. QPMS_NORMALISATION_L_I = 256, ///< Include an additional \a i -factor into the longitudinal waves.
QPMS_NORMALISATION_L_MINUS = 512, ///< Include an additional \f$-1\f$ -factor into the longitudinal waves. QPMS_NORMALISATION_L_MINUS = 512, ///< Include an additional \f$-1\f$ -factor into the longitudinal waves.
#endif
QPMS_NORMALISATION_NORM_BITSTART = 65536, QPMS_NORMALISATION_NORM_BITSTART = 65536,
/// The VSWFs shall be power-normalised. This is the "default". /// The VSWFs shall be power-normalised. This is the "default".
/** /**

92
qpms/vswf.c
View File

@ -8,6 +8,7 @@
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
#include "qpms_error.h" #include "qpms_error.h"
#include "normalisation.h"
qpms_vswf_set_spec_t *qpms_vswf_set_spec_init() { qpms_vswf_set_spec_t *qpms_vswf_set_spec_init() {
@ -98,35 +99,45 @@ 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) { qpms_bessel_t btyp, qpms_normalisation_t norm) {
lmcheck(l,m); lmcheck(l,m);
csphvec_t N; csphvec_t N;
complex double *bessel = malloc((l+1)*sizeof(complex double)); complex double *bessel;
if(qpms_sph_bessel_fill(btyp, l, kdlj.r, bessel)) abort(); QPMS_CRASHING_MALLOC(bessel,(l+1)*sizeof(complex double));
qpms_pitau_t pt = qpms_pitau_get(kdlj.theta, l, norm); QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(btyp, l, kdlj.r, bessel));
complex double eimf = cexp(m * kdlj.phi * I); qpms_pitau_t pt = qpms_pitau_get(kdlj.theta, l, qpms_normalisation_t_csphase(norm));
complex double eimf = qpms_spharm_azimuthal_part(norm, m, kdlj.phi);
complex double d_eimf_dmf = qpms_spharm_azimuthal_part_derivative_div_m(norm, m, kdlj.phi);
qpms_y_t y = qpms_mn2y(m,l); qpms_y_t y = qpms_mn2y(m,l);
N.rc = l*(l+1) * pt.leg[y] * bessel[l] / kdlj.r * eimf; N.rc = l*(l+1) * pt.leg[y] * bessel[l] / kdlj.r * eimf;
complex double besselfac = bessel[l-1] - l * bessel[l] / kdlj.r; complex double besselfac = bessel[l-1] - l * bessel[l] / kdlj.r;
N.thetac = pt.tau[y] * besselfac * eimf; N.thetac = pt.tau[y] * besselfac * eimf;
N.phic = pt.pi[y] * besselfac * I * eimf; N.phic = pt.pi[y] * besselfac * d_eimf_dmf;
N = csphvec_scale(qpms_normalisation_factor_N_noCS(norm, l, m), N);
qpms_pitau_free(pt); qpms_pitau_free(pt);
free(bessel); free(bessel);
return N; return N;
} }
csphvec_t qpms_vswf_single_mg(qpms_m_t m, qpms_l_t l, sph_t kdlj, csphvec_t qpms_vswf_single_mg(qpms_m_t m, qpms_l_t l, sph_t kdlj,
qpms_bessel_t btyp, qpms_normalisation_t norm) { qpms_bessel_t btyp, qpms_normalisation_t norm) {
lmcheck(l,m); lmcheck(l,m);
csphvec_t M; csphvec_t M;
complex double *bessel = malloc((l+1)*sizeof(complex double)); complex double *bessel;
if(qpms_sph_bessel_fill(btyp, l, kdlj.r, bessel)) abort(); QPMS_CRASHING_MALLOC(bessel,(l+1)*sizeof(complex double));
qpms_pitau_t pt = qpms_pitau_get(kdlj.theta, l, norm); QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(btyp, l, kdlj.r, bessel));
complex double eimf = cexp(m * kdlj.phi * I); qpms_pitau_t pt = qpms_pitau_get(kdlj.theta, l, qpms_normalisation_t_csphase(norm));
complex double eimf = qpms_spharm_azimuthal_part(norm, m, kdlj.phi);
complex double d_eimf_dmf = qpms_spharm_azimuthal_part_derivative_div_m(norm, m, kdlj.phi);
qpms_y_t y = qpms_mn2y(m,l); qpms_y_t y = qpms_mn2y(m,l);
M.rc = 0.; M.rc = 0.;
M.thetac = pt.pi[y] * bessel[l] * I * eimf; M.thetac = pt.pi[y] * bessel[l] * d_eimf_dmf;
M.phic = -pt.tau[y] * bessel[l] * eimf; M.phic = -pt.tau[y] * bessel[l] * eimf;
M = csphvec_scale(qpms_normalisation_factor_M_noCS(norm, l, m), M);
qpms_pitau_free(pt); qpms_pitau_free(pt);
free(bessel); free(bessel);
return M; return M;
@ -137,10 +148,9 @@ qpms_vswfset_sph_t *qpms_vswfset_make(qpms_l_t lMax, sph_t kdlj,
qpms_vswfset_sph_t *res = malloc(sizeof(qpms_vswfset_sph_t)); qpms_vswfset_sph_t *res = malloc(sizeof(qpms_vswfset_sph_t));
res->lMax = lMax; res->lMax = lMax;
qpms_y_t nelem = qpms_lMax2nelem(lMax); qpms_y_t nelem = qpms_lMax2nelem(lMax);
res->el = malloc(sizeof(csphvec_t)*nelem); QPMS_CRASHING_MALLOC(res->el, sizeof(csphvec_t)*nelem);
res->mg = malloc(sizeof(csphvec_t)*nelem); QPMS_CRASHING_MALLOC(res->mg, sizeof(csphvec_t)*nelem);
if(QPMS_SUCCESS != qpms_vswf_fill(NULL, res->mg, res->el, lMax, kdlj, btyp, norm)) QPMS_ENSURE_SUCCESS(qpms_vswf_fill(NULL, res->mg, res->el, lMax, kdlj, btyp, norm));
abort(); // or return NULL? or rather assert?
return res; return res;
} }
@ -163,11 +173,11 @@ csphvec_t qpms_vswf_L00(csph_t kr, qpms_bessel_t btyp,
qpms_errno_t qpms_vswf_fill_csph(csphvec_t *const longtarget, qpms_errno_t qpms_vswf_fill_csph(csphvec_t *const longtarget,
csphvec_t * const mgtarget, csphvec_t * const eltarget, qpms_l_t lMax, csphvec_t * const mgtarget, csphvec_t * const eltarget, qpms_l_t lMax,
csph_t kr, qpms_bessel_t btyp, qpms_normalisation_t norm) { csph_t kr, qpms_bessel_t btyp, const qpms_normalisation_t norm) {
assert(lMax >= 1); assert(lMax >= 1);
complex double *bessel = malloc((lMax+1)*sizeof(complex double)); complex double *bessel = malloc((lMax+1)*sizeof(complex double));
if(qpms_sph_bessel_fill(btyp, lMax, kr.r, bessel)) abort(); if(qpms_sph_bessel_fill(btyp, lMax, kr.r, bessel)) abort();
qpms_pitau_t pt = qpms_pitau_get(kr.theta, lMax, norm); qpms_pitau_t pt = qpms_pitau_get(kr.theta, lMax, qpms_normalisation_t_csphase(norm));
complex double const *pbes = bessel + 1; // starting from l = 1 complex double const *pbes = bessel + 1; // starting from l = 1
double const *pleg = pt.leg; double const *pleg = pt.leg;
double const *ppi = pt.pi; double const *ppi = pt.pi;
@ -183,28 +193,32 @@ qpms_errno_t qpms_vswf_fill_csph(csphvec_t *const longtarget,
} }
besderfac = *(pbes-1) - l * besfac; besderfac = *(pbes-1) - l * besfac;
for(qpms_m_t m = -l; m <= l; ++m) { for(qpms_m_t m = -l; m <= l; ++m) {
complex double eimf = cexp(m * kr.phi * I); complex double eimf = qpms_spharm_azimuthal_part(norm, m, kr.phi);
complex double d_eimf_dmf = qpms_spharm_azimuthal_part_derivative_div_m(norm, m, kr.phi);
if (longtarget) { QPMS_UNTESTED; if (longtarget) { QPMS_UNTESTED;
complex double longfac = sqrt(l*(l+1)) * eimf; double longfac = sqrt(l*(l+1));
plong->rc = // FATAL FIXME: I get wrong result here for plane wave re-expansion plong->rc = // FATAL FIXME: I get wrong result here for plane wave re-expansion
// whenever kr.r > 0 (for waves with longitudinal component, ofcoz) // whenever kr.r > 0 (for waves with longitudinal component, ofcoz)
/*(*(pbes-1) - (l+1)/kr.r* *pbes)*/ /*(*(pbes-1) - (l+1)/kr.r* *pbes)*/
(besderfac-besfac) (besderfac-besfac)
* (*pleg) * longfac; * (*pleg) * longfac * eimf;
plong->thetac = *ptau * besfac * longfac; plong->thetac = *ptau * besfac * longfac * eimf;
plong->phic = *ppi * I * besfac * longfac; plong->phic = *ppi * besfac * longfac * d_eimf_dmf;
*plong = csphvec_scale(qpms_normalisation_factor_L_noCS(norm, l, m), *plong);
++plong; ++plong;
} }
if (eltarget) { if (eltarget) {
pel->rc = l*(l+1) * (*pleg) * besfac * eimf; pel->rc = l*(l+1) * (*pleg) * besfac * eimf;
pel->thetac = *ptau * besderfac * eimf; pel->thetac = *ptau * besderfac * eimf;
pel->phic = *ppi * besderfac * I * eimf; pel->phic = *ppi * besderfac * d_eimf_dmf;
*pel = csphvec_scale(qpms_normalisation_factor_N_noCS(norm, l, m), *pel);
++pel; ++pel;
} }
if (mgtarget) { if (mgtarget) {
pmg->rc = 0.; pmg->rc = 0.;
pmg->thetac = *ppi * (*pbes) * I * eimf; pmg->thetac = *ppi * (*pbes) * d_eimf_dmf;
pmg->phic = - *ptau * (*pbes) * eimf; pmg->phic = - *ptau * (*pbes) * eimf;
*pmg = csphvec_scale(qpms_normalisation_factor_M_noCS(norm, l, m), *pmg);
++pmg; ++pmg;
} }
++pleg; ++ppi; ++ptau; ++pleg; ++ppi; ++ptau;
@ -234,7 +248,9 @@ qpms_errno_t qpms_vswf_fill_alternative(csphvec_t *const longtarget, csphvec_t *
complex double const *pbes = bessel + 1; // starting from l = 1 complex double const *pbes = bessel + 1; // starting from l = 1
qpms_y_t nelem = qpms_lMax2nelem(lMax); 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; csphvec_t *a;
QPMS_CRASHING_MALLOC(a, 3*nelem*sizeof(csphvec_t))
csphvec_t * const a1 = a, * const a2 = a1 + nelem, * const a3 = a2 + 2 * nelem;
if(qpms_vecspharm_fill(a1, a2, a3, lMax, kr, norm)) abort(); if(qpms_vecspharm_fill(a1, a2, a3, lMax, kr, norm)) abort();
const csphvec_t *p1 = a1; const csphvec_t *p1 = a1;
const csphvec_t *p2 = a2; const csphvec_t *p2 = a2;
@ -246,10 +262,12 @@ qpms_errno_t qpms_vswf_fill_alternative(csphvec_t *const longtarget, csphvec_t *
complex double besderfac = *(pbes-1) - l * besfac; complex double besderfac = *(pbes-1) - l * besfac;
double sqrtlfac = sqrt(l*(l+1)); double sqrtlfac = sqrt(l*(l+1));
for(qpms_m_t m = -l; m <= l; ++m) { for(qpms_m_t m = -l; m <= l; ++m) {
complex double eimf = cexp(m * kr.phi * I); // FIXME unused variable?!!!
if (longtarget) { if (longtarget) {
complex double L2Nfac = qpms_normalisation_factor_L_noCS(norm, l, m)
/ qpms_normalisation_factor_N_noCS(norm, l, m);
*plong = csphvec_add(csphvec_scale(besderfac-besfac, *p3), *plong = csphvec_add(csphvec_scale(besderfac-besfac, *p3),
csphvec_scale(sqrtlfac * besfac, *p2)); csphvec_scale(sqrtlfac * besfac, *p2));
*plong = csphvec_scale(L2Nfac, *plong);
++plong; ++plong;
} }
if (eltarget) { if (eltarget) {
@ -265,7 +283,7 @@ qpms_errno_t qpms_vswf_fill_alternative(csphvec_t *const longtarget, csphvec_t *
} }
++pbes; ++pbes;
} }
free(a1); free(a);
free(bessel); free(bessel);
return QPMS_SUCCESS; return QPMS_SUCCESS;
} }
@ -273,28 +291,31 @@ qpms_errno_t qpms_vswf_fill_alternative(csphvec_t *const longtarget, csphvec_t *
qpms_errno_t qpms_vecspharm_fill(csphvec_t *const a1target, csphvec_t *const a2target, csphvec_t *const a3target, 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_l_t lMax, sph_t dir, qpms_normalisation_t norm) {
assert(lMax >= 1); assert(lMax >= 1);
qpms_pitau_t pt = qpms_pitau_get(dir.theta, lMax, norm); qpms_pitau_t pt = qpms_pitau_get(dir.theta, lMax, qpms_normalisation_t_csphase(norm));
double const *pleg = pt.leg; double const *pleg = pt.leg;
double const *ppi = pt.pi; double const *ppi = pt.pi;
double const *ptau = pt.tau; double const *ptau = pt.tau;
csphvec_t *p1 = a1target, *p2 = a2target, *p3 = a3target; csphvec_t *p1 = a1target, *p2 = a2target, *p3 = a3target;
for (qpms_l_t l = 1; l <= lMax; ++l) { for (qpms_l_t l = 1; l <= lMax; ++l) {
for(qpms_m_t m = -l; m <= l; ++m) { for(qpms_m_t m = -l; m <= l; ++m) {
complex double eimf = cexp(m * dir.phi * I); const complex double Mfac = qpms_normalisation_factor_M_noCS(norm, l, m);
const complex double Nfac = qpms_normalisation_factor_N_noCS(norm, l, m);
const complex double eimf = qpms_spharm_azimuthal_part(norm, m, dir.phi);
const complex double deimf_dmf = qpms_spharm_azimuthal_part_derivative_div_m(norm, m, dir.phi);
if (a1target) { if (a1target) {
p1->rc = 0; p1->rc = 0;
p1->thetac = *ppi * I * eimf; p1->thetac = *ppi * deimf_dmf * Mfac;
p1->phic = -*ptau * eimf; p1->phic = -*ptau * eimf * Mfac;
++p1; ++p1;
} }
if (a2target) { if (a2target) {
p2->rc = 0; p2->rc = 0;
p2->thetac = *ptau * eimf; p2->thetac = *ptau * eimf * Nfac;
p2->phic = *ppi * I * eimf; p2->phic = *ppi * deimf_dmf * Nfac;
++p2; ++p2;
} }
if (a3target) { if (a3target) {
p3->rc = sqrt(l*(l+1)) * (*pleg) * eimf; p3->rc = sqrt(l*(l+1)) * (*pleg) * eimf * Nfac;
p3->thetac = 0; p3->thetac = 0;
p3->phic = 0; p3->phic = 0;
++p3; ++p3;
@ -308,6 +329,10 @@ qpms_errno_t qpms_vecspharm_fill(csphvec_t *const a1target, csphvec_t *const a2t
qpms_errno_t qpms_vecspharm_dual_fill(csphvec_t *const a1target, csphvec_t *const a2target, csphvec_t *const a3target, 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_l_t lMax, sph_t dir, qpms_normalisation_t norm) {
#if 1
return qpms_vecspharm_fill(a1target, a2target, a3target, lMax, dir,
qpms_normalisation_dual(norm));
#else
assert(lMax >= 1); assert(lMax >= 1);
qpms_pitau_t pt = qpms_pitau_get(dir.theta, lMax, norm); qpms_pitau_t pt = qpms_pitau_get(dir.theta, lMax, norm);
double const *pleg = pt.leg; double const *pleg = pt.leg;
@ -341,6 +366,7 @@ qpms_errno_t qpms_vecspharm_dual_fill(csphvec_t *const a1target, csphvec_t *cons
} }
qpms_pitau_free(pt); qpms_pitau_free(pt);
return QPMS_SUCCESS; return QPMS_SUCCESS;
#endif
} }

11
qpms/vswf.h
View File

@ -11,7 +11,8 @@
#include "qpms_types.h" #include "qpms_types.h"
#include <gsl/gsl_sf_legendre.h> #include <gsl/gsl_sf_legendre.h>
// Methods for qpms_vswf_spec_t // ---------------Methods for qpms_vswf_spec_t-----------------------
//
/// Creates a qpms_vswf_set_spec_t structure with an empty list of wave indices. /// Creates a qpms_vswf_set_spec_t structure with an empty list of wave indices.
qpms_vswf_set_spec_t *qpms_vswf_set_spec_init(void); qpms_vswf_set_spec_t *qpms_vswf_set_spec_init(void);
/// Appends a VSWF index to a \ref qpms_vswf_set_spec_t, also updating metadata. /// Appends a VSWF index to a \ref qpms_vswf_set_spec_t, also updating metadata.
@ -60,6 +61,8 @@ csphvec_t qpms_eval_uvswf(const qpms_vswf_set_spec_t *setspec,
csph_t kr, ///< Evaluation point. csph_t kr, ///< Evaluation point.
qpms_bessel_t btyp); qpms_bessel_t btyp);
// -----------------------------------------------------------------------
/// Electric wave N. /// Electric wave N.
csphvec_t qpms_vswf_single_el(int m, int n, sph_t kdlj, csphvec_t qpms_vswf_single_el(int m, int n, sph_t kdlj,
qpms_bessel_t btyp, qpms_normalisation_t norm); qpms_bessel_t btyp, qpms_normalisation_t norm);
@ -87,10 +90,14 @@ qpms_errno_t qpms_legendre_deriv_y_fill(double *where, double *where_deriv, doub
/// Evaluate the zeroth-degree longitudinal VSWF \f$ \mathbf{L}_0^0 \f$. /// Evaluate the zeroth-degree longitudinal VSWF \f$ \mathbf{L}_0^0 \f$.
/**
* Any `norm` is being ignored right now.
*/
csphvec_t qpms_vswf_L00( csphvec_t qpms_vswf_L00(
csph_t kdrj, //< VSWF evaluation point. csph_t kdrj, //< VSWF evaluation point.
qpms_bessel_t btyp, qpms_bessel_t btyp,
qpms_normalisation_t norm); qpms_normalisation_t norm //< Ignored!
);
/// Evaluate VSWFs at a given point from \a l = 1 up to a given degree \a lMax. /// Evaluate VSWFs at a given point from \a l = 1 up to a given degree \a lMax.
/** /**