From 1bf9c0b44e600d7bff1a7866d882205cfeb5d9c1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 13 Aug 2019 17:18:39 +0300 Subject: [PATCH] Cython bindings for "pi, tau" functions Former-commit-id: 74edc5a7a431a190150c4e1cbb3563809c6b9c01 --- qpms/__init__.py | 2 +- qpms/legendre.c | 49 +++++++++++++++++++++++++++++--------------- qpms/qpms_c.pyx | 13 +++++++++++- qpms/qpms_cdefs.pxd | 13 ++++++++++++ qpms/qpms_specfunc.h | 12 ++++++++++- 5 files changed, 70 insertions(+), 19 deletions(-) diff --git a/qpms/__init__.py b/qpms/__init__.py index cc43bed..5f73da3 100644 --- a/qpms/__init__.py +++ b/qpms/__init__.py @@ -7,7 +7,7 @@ from sys import platform as __platform import warnings as __warnings try: - from .qpms_c import PointGroup, FinitePointGroup, FinitePointGroupElement, Particle, scatsystem_set_nthreads, ScatteringSystem, ScatteringMatrix + from .qpms_c import PointGroup, FinitePointGroup, FinitePointGroupElement, Particle, scatsystem_set_nthreads, ScatteringSystem, ScatteringMatrix, pitau except ImportError as ex: if __platform == "linux" or __platform == "linux2": if 'LD_LIBRARY_PATH' not in __os.environ.keys(): diff --git a/qpms/legendre.c b/qpms/legendre.c index 9d3f19a..a440897 100644 --- a/qpms/legendre.c +++ b/qpms/legendre.c @@ -53,50 +53,67 @@ qpms_errno_t qpms_legendre_deriv_y_get(double **target, double **dtarget, double qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, const double csphase) { - QPMS_ENSURE(fabs(csphase) == 1, "The csphase argument must be either 1 or -1, not %g.", csphase); qpms_pitau_t res; const qpms_y_t nelem = qpms_lMax2nelem(lMax); QPMS_CRASHING_MALLOC(res.leg, nelem * sizeof(double)); QPMS_CRASHING_MALLOC(res.pi, nelem * sizeof(double)); QPMS_CRASHING_MALLOC(res.tau, nelem * sizeof(double)); + qpms_pitau_fill(res.leg, res.pi, res.tau, theta, lMax, csphase); + return res; +} + +qpms_errno_t qpms_pitau_fill(double *target_leg, double *pi, double *tau, double theta, qpms_l_t lMax, double csphase) +{ + QPMS_ENSURE(fabs(csphase) == 1, "The csphase argument must be either 1 or -1, not %g.", csphase); + const qpms_y_t nelem = qpms_lMax2nelem(lMax); + double ct = cos(theta), st = sin(theta); if (1 == fabs(ct)) { // singular case, use DLMF 14.8.2 - memset(res.pi, 0, nelem * sizeof(double)); - memset(res.tau, 0, nelem * sizeof(double)); - memset(res.leg, 0, nelem * sizeof(double)); + if(pi) memset(pi, 0, nelem * sizeof(double)); + if(tau) memset(tau, 0, nelem * sizeof(double)); + if(target_leg) memset(target_leg, 0, nelem * sizeof(double)); 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))); + if(target_leg) target_leg[qpms_mn2y(0, l)] = ((l%2)?ct:1.)*sqrt((2*l+1)/(4*M_PI *l*(l+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; + if(pi) { + pi [qpms_mn2y(+1, l)] = -((ct>0) ? -1 : lpar) * fl * csphase; + pi [qpms_mn2y(-1, l)] = -((ct>0) ? -1 : lpar) * fl * csphase; + } + if(tau) { + tau[qpms_mn2y(+1, l)] = ((ct>0) ? +1 : lpar) * fl * csphase; + tau[qpms_mn2y(-1, l)] = -((ct>0) ? +1 : lpar) * fl * csphase; + } } } else { // cos(theta) in (-1,1), use normal calculation - double *legder; + double *leg, *legder; + if (target_leg) + leg = target_leg; + else + QPMS_CRASHING_MALLOC(leg, nelem*sizeof(double)); QPMS_CRASHING_MALLOC(legder, nelem * sizeof(double)); - QPMS_ENSURE_SUCCESS(qpms_legendre_deriv_y_fill(res.leg, legder, ct, lMax, + QPMS_ENSURE_SUCCESS(qpms_legendre_deriv_y_fill(leg, legder, ct, lMax, GSL_SF_LEGENDRE_SPHARM, csphase)); // Multiply by the "power normalisation" factor for (qpms_l_t l = 1; l <= lMax; ++l) { double prefac = 1./sqrt(l*(l+1)); for (qpms_m_t m = -l; m <= l; ++m) { - res.leg[qpms_mn2y(m,l)] *= prefac; + leg[qpms_mn2y(m,l)] *= prefac; legder[qpms_mn2y(m,l)] *= prefac; } } 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)]; - res.tau[qpms_mn2y(m,l)] = - st * legder[qpms_mn2y(m,l)]; + if(pi) pi [qpms_mn2y(m,l)] = m / st * leg[qpms_mn2y(m,l)]; + if(tau) tau[qpms_mn2y(m,l)] = - st * legder[qpms_mn2y(m,l)]; } } free(legder); + if(!target_leg) + free(leg); } - res.lMax = lMax; - return res; + return QPMS_SUCCESS; } void qpms_pitau_free(qpms_pitau_t x) { diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index cd23411..58941cf 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -488,4 +488,15 @@ cdef class ScatteringMatrix: qpms_scatsys_scatter_solve(&f_view[0], &a_view[0], self.lu) return f - +def pitau(double theta, qpms_l_t lMax, double csphase = 1): + if(abs(csphase) != 1): + raise ValueError("csphase must be 1 or -1, is %g" % csphase) + cdef size_t nelem = qpms_lMax2nelem(lMax) + cdef np.ndarray[np.float_t, ndim=1] lega = np.empty((nelem,), dtype=float) + cdef np.ndarray[np.float_t, ndim=1] pia = np.empty((nelem,), dtype=float) + cdef np.ndarray[np.float_t, ndim=1] taua = np.empty((nelem,), dtype=float) + cdef double[::1] leg = lega + cdef double[::1] pi = pia + cdef double[::1] tau = taua + qpms_pitau_fill(&leg[0], &pi[0], &tau[0], theta, lMax, csphase) + return (lega, pia, taua) diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 0aee2a3..fc7d7b6 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -65,6 +65,7 @@ cdef extern from "qpms_types.h": ctypedef int qpms_lm_t ctypedef int qpms_l_t ctypedef int qpms_m_t + ctypedef size_t qpms_y_t struct qpms_quat_t: cdouble a cdouble b @@ -161,6 +162,7 @@ cdef extern from "vswf.h": const cart3_t evalpoint, const void *args, bint add) cdef extern from "indexing.h": + qpms_y_t qpms_lMax2nelem(qpms_l_t lMax) qpms_uvswfi_t qpms_tmn2uvswfi(qpms_vswf_type_t t, qpms_m_t m, qpms_l_t n) qpms_errno_t qpms_uvswfi2tmn(qpms_uvswfi_t u, qpms_vswf_type_t* t, qpms_m_t* m, qpms_l_t* n) qpms_m_t qpms_uvswfi2m(qpms_uvswfi_t u) @@ -313,6 +315,17 @@ cdef extern from "translations.h": qpms_bessel_t J ); +cdef extern from "qpms_specfunc.h": + struct qpms_pitau_t: + qpms_l_t lMax + double *leg + double *pi + double *tau + qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, double csphase) + qpms_errno_t qpms_pitau_fill(double *leg, double *pi, double *tau, double theta, qpms_l_t lMax, double csphase) + void qpms_pitau_free(qpms_pitau_t pitau) + + cdef extern from "gsl/gsl_interp.h": struct gsl_interp_type: pass diff --git a/qpms/qpms_specfunc.h b/qpms/qpms_specfunc.h index 9db24dd..3c07a91 100644 --- a/qpms/qpms_specfunc.h +++ b/qpms/qpms_specfunc.h @@ -68,7 +68,7 @@ typedef struct { double *leg, *pi, *tau; } qpms_pitau_t; -/// Returns an array of normalised Legendre and and auxillary \f$\pi_{lm}, \tau_{lm} \f$ functions. +/// Returns an array of normalised Legendre and auxillary \f$\pi_{lm}, \tau_{lm} \f$ functions. /** * The normalised Legendre function here is defined as * \f[ @@ -88,8 +88,18 @@ typedef struct { * */ qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, double csphase); + +/// Directly fills (pre-allocated) arrays of normalised Legendre and auxillary \f$\pi_{lm}, \tau_{lm} \f$ functions. +/** + * Arrays must be preallocated for `lMax * (lMax + 2)` elements. `NULL` targets are skipped. + * For details, see qpms_pitau_get(). + */ +qpms_errno_t qpms_pitau_fill(double *target_leg, double *target_pi, double *target_tau, + double theta, qpms_l_t lMax, double csphase); + /// Frees the dynamically allocated arrays from qpms_pitau_t. void qpms_pitau_free(qpms_pitau_t); + //void qpms_pitau_pfree(qpms_pitau_t*);//NI // Associated Legendre polynomial at zero argument (DLMF 14.5.1) DEPRECATED?