From bd40daccb7e7e4093612a0d8f39e26b73d34a9b6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Thu, 11 Jul 2019 16:57:15 +0300 Subject: [PATCH] Rewrite vswf.c to use the new qpms_normalisation_t Former-commit-id: a623d0cf1c65b6134756e76e1739f189f3d6d53f --- qpms/normalisation.h | 54 +++++++++++++++++++++----- qpms/qpms_specfunc.h | 3 ++ qpms/qpms_types.h | 2 - qpms/vswf.c | 92 ++++++++++++++++++++++++++++---------------- qpms/vswf.h | 11 +++++- 5 files changed, 116 insertions(+), 46 deletions(-) diff --git a/qpms/normalisation.h b/qpms/normalisation.h index 02970df..8dad9da 100644 --- a/qpms/normalisation.h +++ b/qpms/normalisation.h @@ -7,6 +7,7 @@ #define NORMALISATION_H #include "qpms_types.h" +#include "qpms_error.h" #include #include @@ -49,7 +50,7 @@ static inline complex double qpms_normalisation_factor_M_noCS(qpms_normalisation * 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) { - complex double fac = qn; + complex double fac = qpms_normalisation_factor_M_noCS(norm, l, m); 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. */ 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; } -#if 0 /// 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 @@ -101,10 +101,9 @@ static inline complex double qpms_normalisation_factor_L_noCS(qpms_normalisation * 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) { - complex double fac = qn; + complex double fac = qpms_normalisation_factor_L_noCS(norm, l, m); return ((norm & QPMS_NORMALISATION_CSPHASE) && (m % 2)) ? -fac : fac; } -#endif /// 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 \f[ e^{im\phi} \f] for standard 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 * \f[ * \sqrt{2}\cos{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] */ 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: return cexp(I*m*phi); - case QPMS_NORMALISATION_REVERSE_ASIMUTHAL_PHASE: + case QPMS_NORMALISATION_REVERSE_AZIMUTHAL_PHASE: return cexp(-I*m*phi); case QPMS_NORMALISATION_SPHARM_REAL: 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 diff --git a/qpms/qpms_specfunc.h b/qpms/qpms_specfunc.h index 2ca9ad6..3db4b1b 100644 --- a/qpms/qpms_specfunc.h +++ b/qpms/qpms_specfunc.h @@ -1,3 +1,6 @@ +/*! \file qpms_specfunc.h + * \brief Various special and auxillary functions. + */ #ifndef QPMS_SPECFUNC_H #define QPMS_SPECFUNC_H #include "qpms_types.h" diff --git a/qpms/qpms_types.h b/qpms/qpms_types.h index d297a32..69bc50b 100644 --- a/qpms/qpms_types.h +++ b/qpms/qpms_types.h @@ -119,10 +119,8 @@ typedef enum { 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_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_MINUS = 512, ///< Include an additional \f$-1\f$ -factor into the longitudinal waves. -#endif QPMS_NORMALISATION_NORM_BITSTART = 65536, /// The VSWFs shall be power-normalised. This is the "default". /** diff --git a/qpms/vswf.c b/qpms/vswf.c index 86afc0a..a77925f 100644 --- a/qpms/vswf.c +++ b/qpms/vswf.c @@ -8,6 +8,7 @@ #include #include #include "qpms_error.h" +#include "normalisation.h" 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) { lmcheck(l,m); csphvec_t N; - complex double *bessel = malloc((l+1)*sizeof(complex double)); - if(qpms_sph_bessel_fill(btyp, l, kdlj.r, bessel)) abort(); - qpms_pitau_t pt = qpms_pitau_get(kdlj.theta, l, norm); - complex double eimf = cexp(m * kdlj.phi * I); + complex double *bessel; + QPMS_CRASHING_MALLOC(bessel,(l+1)*sizeof(complex double)); + QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(btyp, l, kdlj.r, bessel)); + 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); N.rc = l*(l+1) * pt.leg[y] * bessel[l] / kdlj.r * eimf; complex double besselfac = bessel[l-1] - l * bessel[l] / kdlj.r; 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); free(bessel); return N; } + 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) { lmcheck(l,m); csphvec_t M; - complex double *bessel = malloc((l+1)*sizeof(complex double)); - if(qpms_sph_bessel_fill(btyp, l, kdlj.r, bessel)) abort(); - qpms_pitau_t pt = qpms_pitau_get(kdlj.theta, l, norm); - complex double eimf = cexp(m * kdlj.phi * I); + complex double *bessel; + QPMS_CRASHING_MALLOC(bessel,(l+1)*sizeof(complex double)); + QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(btyp, l, kdlj.r, bessel)); + 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); 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 = csphvec_scale(qpms_normalisation_factor_M_noCS(norm, l, m), M); + qpms_pitau_free(pt); free(bessel); 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)); res->lMax = lMax; qpms_y_t nelem = qpms_lMax2nelem(lMax); - res->el = malloc(sizeof(csphvec_t)*nelem); - res->mg = malloc(sizeof(csphvec_t)*nelem); - if(QPMS_SUCCESS != qpms_vswf_fill(NULL, res->mg, res->el, lMax, kdlj, btyp, norm)) - abort(); // or return NULL? or rather assert? + QPMS_CRASHING_MALLOC(res->el, sizeof(csphvec_t)*nelem); + QPMS_CRASHING_MALLOC(res->mg, sizeof(csphvec_t)*nelem); + QPMS_ENSURE_SUCCESS(qpms_vswf_fill(NULL, res->mg, res->el, lMax, kdlj, btyp, norm)); 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, 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); complex double *bessel = malloc((lMax+1)*sizeof(complex double)); 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 double const *pleg = pt.leg; 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; 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; - 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 // whenever kr.r > 0 (for waves with longitudinal component, ofcoz) /*(*(pbes-1) - (l+1)/kr.r* *pbes)*/ (besderfac-besfac) - * (*pleg) * longfac; - plong->thetac = *ptau * besfac * longfac; - plong->phic = *ppi * I * besfac * longfac; + * (*pleg) * longfac * eimf; + plong->thetac = *ptau * besfac * longfac * eimf; + plong->phic = *ppi * besfac * longfac * d_eimf_dmf; + *plong = csphvec_scale(qpms_normalisation_factor_L_noCS(norm, l, m), *plong); ++plong; } if (eltarget) { pel->rc = l*(l+1) * (*pleg) * besfac * 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; } if (mgtarget) { pmg->rc = 0.; - pmg->thetac = *ppi * (*pbes) * I * eimf; + pmg->thetac = *ppi * (*pbes) * d_eimf_dmf; pmg->phic = - *ptau * (*pbes) * eimf; + *pmg = csphvec_scale(qpms_normalisation_factor_M_noCS(norm, l, m), *pmg); ++pmg; } ++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 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(); const csphvec_t *p1 = a1; 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; double sqrtlfac = sqrt(l*(l+1)); for(qpms_m_t m = -l; m <= l; ++m) { - complex double eimf = cexp(m * kr.phi * I); // FIXME unused variable?!!! 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), csphvec_scale(sqrtlfac * besfac, *p2)); + *plong = csphvec_scale(L2Nfac, *plong); ++plong; } if (eltarget) { @@ -265,7 +283,7 @@ qpms_errno_t qpms_vswf_fill_alternative(csphvec_t *const longtarget, csphvec_t * } ++pbes; } - free(a1); + free(a); free(bessel); 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_l_t lMax, sph_t dir, qpms_normalisation_t norm) { 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 *ppi = pt.pi; double const *ptau = pt.tau; csphvec_t *p1 = a1target, *p2 = a2target, *p3 = a3target; for (qpms_l_t l = 1; l <= lMax; ++l) { 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) { p1->rc = 0; - p1->thetac = *ppi * I * eimf; - p1->phic = -*ptau * eimf; + p1->thetac = *ppi * deimf_dmf * Mfac; + p1->phic = -*ptau * eimf * Mfac; ++p1; } if (a2target) { p2->rc = 0; - p2->thetac = *ptau * eimf; - p2->phic = *ppi * I * eimf; + p2->thetac = *ptau * eimf * Nfac; + p2->phic = *ppi * deimf_dmf * Nfac; ++p2; } if (a3target) { - p3->rc = sqrt(l*(l+1)) * (*pleg) * eimf; + p3->rc = sqrt(l*(l+1)) * (*pleg) * eimf * Nfac; p3->thetac = 0; p3->phic = 0; ++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_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); qpms_pitau_t pt = qpms_pitau_get(dir.theta, lMax, norm); 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); return QPMS_SUCCESS; +#endif } diff --git a/qpms/vswf.h b/qpms/vswf.h index ca5da8c..23d1224 100644 --- a/qpms/vswf.h +++ b/qpms/vswf.h @@ -11,7 +11,8 @@ #include "qpms_types.h" #include -// 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. 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. @@ -60,6 +61,8 @@ csphvec_t qpms_eval_uvswf(const qpms_vswf_set_spec_t *setspec, csph_t kr, ///< Evaluation point. qpms_bessel_t btyp); +// ----------------------------------------------------------------------- + /// Electric wave N. csphvec_t qpms_vswf_single_el(int m, int n, sph_t kdlj, 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$. +/** + * Any `norm` is being ignored right now. + */ csphvec_t qpms_vswf_L00( csph_t kdrj, //< VSWF evaluation point. 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. /**