From 00f9f5234a32ba90c7edb7923c30e7c575b9edc9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Thu, 1 Feb 2018 04:40:45 +0000 Subject: [PATCH] =?UTF-8?q?Jdu=20do=20pr=C3=A1ce?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Former-commit-id: 3be0f5adda5334dae7397d21b6d48835a190632d --- qpms/qpms_types.h | 48 +++++++++++++ qpms/vectors.h | 15 ++++ qpms/vswf.c | 176 +++++++++++++++++++++++++++++++++++++++++++--- qpms/vswf.h | 21 +++++- 4 files changed, 250 insertions(+), 10 deletions(-) diff --git a/qpms/qpms_types.h b/qpms/qpms_types.h index e041b73..c37cf90 100644 --- a/qpms/qpms_types.h +++ b/qpms/qpms_types.h @@ -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; diff --git a/qpms/vectors.h b/qpms/vectors.h index 1724f23..5f965fa 100644 --- a/qpms/vectors.h +++ b/qpms/vectors.h @@ -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}; diff --git a/qpms/vswf.c b/qpms/vswf.c index f49e700..5b69477 100644 --- a/qpms/vswf.c +++ b/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; +} + diff --git a/qpms/vswf.h b/qpms/vswf.h index d76d8cb..b159972 100644 --- a/qpms/vswf.h +++ b/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