From 7a0386086a0b44234972c2bac54bba954c87e856 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 13 Aug 2019 13:30:48 +0300 Subject: [PATCH] Functions to calculate the R,Q matrices for axially symmetric particles separately. Former-commit-id: 0861a9014f293153b786b879bb3803d5ee002b56 --- qpms/cytmatrices.pyx | 27 ++++++++++++++- qpms/qpms_cdefs.pxd | 2 ++ qpms/tmatrices.c | 82 +++++++++++++++++++++++++++++++++++++++++--- qpms/tmatrices.h | 16 +++++++++ 4 files changed, 122 insertions(+), 5 deletions(-) diff --git a/qpms/cytmatrices.pyx b/qpms/cytmatrices.pyx index aa4bd4a..d8fa27b 100644 --- a/qpms/cytmatrices.pyx +++ b/qpms/cytmatrices.pyx @@ -204,6 +204,22 @@ cdef class __AxialSymParams: self.lMax_extend = args[0] if 'lMax_extend' in kwargs.keys(): self.lMax_extend = kwargs['lMax_extend'] + if self.lMax_extend == 0: + self.lMax_extend = 1 + def Q_transposed(self, cdouble omega, norm): + cdef size_t n = 2*(self.p.lMax_extend*(self.p.lMax_extend +2)) + cdef np.ndarray[np.complex_t, ndim=2] arr = np.empty((n,n), dtype=complex, order='C') + cdef cdouble[:,::1] arrview = arr + qpms_tmatrix_generator_axialsym_RQ_transposed_fill(&arrview[0][0], omega, &self.p, norm, QPMS_HANKEL_PLUS) + return arr + def R_transposed(self, cdouble omega, norm): + cdef size_t n = 2*(self.p.lMax_extend*(self.p.lMax_extend +2)) + cdef np.ndarray[np.complex_t, ndim=2] arr = np.empty((n,n), dtype=complex, order='C') + cdef cdouble[:,::1] arrview = arr + qpms_tmatrix_generator_axialsym_RQ_transposed_fill(&arrview[0][0], omega, &self.p, norm, QPMS_BESSEL_REGULAR) + return arr + + cdef class TMatrixGenerator: cdef qpms_tmatrix_generator_t g @@ -221,7 +237,7 @@ cdef class TMatrixGenerator: self.g.params = (<__AxialSymParams?>self.holder).rawpointer() # TODO INTERPOLATOR else: - raise ValueError("Can't construct TMatrixGenerator from that") + raise TypeError("Can't construct TMatrixGenerator from that") def __call__(self, arg, cdouble omega): cdef CTMatrix tm @@ -238,6 +254,15 @@ cdef class TMatrixGenerator: else: raise ValueError("Must specify CTMatrix or BaseSpec") + def Q_transposed(self, cdouble omega, norm): + if self.g.function != qpms_tmatrix_generator_axialsym: + raise TypeError("Applicable only for axialsym generators") + return self.holder.Q_transposed(omega, norm) + def R_transposed(self, cdouble omega, norm): + if self.g.function != qpms_tmatrix_generator_axialsym: + raise TypeError("Applicable only for axialsym generators") + return self.holder.R_transposed(omega, norm) + # Better "constructors": @staticmethod def sphere(outside, inside, r): diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index f43f45f..0aee2a3 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -424,6 +424,8 @@ cdef extern from "tmatrices.h": cdouble epsilon_fg, cdouble epsilon_bg) qpms_tmatrix_t *qpms_tmatrix_spherical_mu0(const qpms_vswf_set_spec_t *bspec, double a, double omega, cdouble epsilon_fg, cdouble epsilon_bg) + qpms_errno_t qpms_tmatrix_generator_axialsym_RQ_transposed_fill(cdouble *target, cdouble omega, + const qpms_tmatrix_generator_axialsym_param_t *param, qpms_normalisation_t norm, qpms_bessel_t J) cdef extern from "pointgroups.h": bint qpms_pg_is_finite_axial(qpms_pointgroup_class cls) diff --git a/qpms/tmatrices.c b/qpms/tmatrices.c index 6d9c309..95d6298 100644 --- a/qpms/tmatrices.c +++ b/qpms/tmatrices.c @@ -724,7 +724,6 @@ qpms_errno_t qpms_tmatrix_axialsym_fill( for(size_t i2 = 0; i2 < bspecQR->n; ++i2) { // NOTE that the Q', R' matrices are !TRANSPOSED! here because of zgetrs const size_t iQR = i1 + i2 * bspecQR->n; - double abserr; // gsl_integration_qag() needs this; thrown away for now qpms_m_t m1, m2; qpms_l_t l1, l2; qpms_vswf_type_t t1, t2; @@ -739,7 +738,8 @@ qpms_errno_t qpms_tmatrix_axialsym_fill( p.l = l1; p.l_in = l2; p.t = t1; p.t_in = t2; - double result; // We throw the quadrature error estimate away. + double result; + double abserr; // gsl_integration_qag() needs this; thrown away for now // Re(R') p.btype = QPMS_BESSEL_REGULAR; p.realpart = true; @@ -794,7 +794,82 @@ qpms_errno_t qpms_tmatrix_axialsym_fill( return QPMS_SUCCESS; } -qpms_errno_t qpms_tmatrix_generator_axialsym(qpms_tmatrix_t *t, complex double omega, const void *param) { +qpms_errno_t qpms_tmatrix_axialsym_RQ_transposed_fill(complex double *target, + complex double omega, qpms_epsmu_t outside, qpms_epsmu_t inside, + qpms_arc_function_t shape, qpms_l_t lMaxQR, qpms_normalisation_t norm,qpms_bessel_t J) +{ + QPMS_UNTESTED; // TODO also check whether this is valid for all norms. + struct tmatrix_axialsym_integral_param_t p; + + qpms_vswf_set_spec_t *bspecQR = qpms_vswf_set_spec_from_lMax(lMaxQR, norm); + + p.k = qpms_wavenumber(omega, outside); + p.k_in = qpms_wavenumber(omega, inside); + p.z = qpms_waveimpedance(outside); + p.z_in = qpms_waveimpedance(inside); + p.f = shape; + p.bspec = bspecQR; + p.btype = J; + const gsl_function f = {tmatrix_axialsym_integrand, (void *) &p}; + + + // Not sure what the size should be, but this should be more than enough. + const size_t intlimit = 1024; + const double epsrel = TMATRIX_AXIALSYM_INTEGRAL_EPSREL; + const double epsabs = TMATRIX_AXIALSYM_INTEGRAL_EPSABS; + gsl_integration_workspace *w = gsl_integration_workspace_alloc(intlimit); + for(size_t i1 = 0; i1 < bspecQR->n; ++i1) + for(size_t i2 = 0; i2 < bspecQR->n; ++i2) { + // NOTE that the Q', R' matrices are !TRANSPOSED! here because of zgetrs + const size_t iQR = i1 + i2 * bspecQR->n; + qpms_m_t m1, m2; + qpms_l_t l1, l2; + qpms_vswf_type_t t1, t2; + const qpms_uvswfi_t u1 = bspecQR->ilist[i1], u2 = bspecQR->ilist[i2]; + QPMS_ENSURE_SUCCESS(qpms_uvswfi2tmn(u1, &t1, &m1, &l1)); + QPMS_ENSURE_SUCCESS(qpms_uvswfi2tmn(u2, &t2, &m2, &l2)); + if (m1 + m2) { + target[iQR] = 0; + } else { + p.m = m1; + p.l = l1; p.l_in = l2; + p.t = t1; p.t_in = t2; + + double result; + double abserr; // gsl_integration_qag() needs this; thrown away for now + // Re(R' or Q') + p.realpart = true; + QPMS_ENSURE_SUCCESS(gsl_integration_qag(&f, 0, M_PI, epsabs, epsrel, + intlimit, 2, w, &result, &abserr)); + target[iQR] = result; + + // Im(R' or Q') + p.realpart = false; + QPMS_ENSURE_SUCCESS(gsl_integration_qag(&f, 0, M_PI, epsabs, epsrel, + intlimit, 2, w, &result, &abserr)); + target[iQR] += I*result; + } + } + gsl_integration_workspace_free(w); + + qpms_vswf_set_spec_free(bspecQR); + return QPMS_SUCCESS; +} + +qpms_errno_t qpms_tmatrix_generator_axialsym_RQ_transposed_fill(complex double *target, + complex double omega, const qpms_tmatrix_generator_axialsym_param_t *param, + qpms_normalisation_t norm, qpms_bessel_t J) +{ + const qpms_tmatrix_generator_axialsym_param_t *p = param; + return qpms_tmatrix_axialsym_RQ_transposed_fill(target, omega, + p->outside.function(omega, p->outside.params), + p->inside.function(omega, p->inside.params), + p->shape, + p->lMax_extend, norm, J); +} + +qpms_errno_t qpms_tmatrix_generator_axialsym(qpms_tmatrix_t *t, complex double omega, const void *param) +{ const qpms_tmatrix_generator_axialsym_param_t *p = param; return qpms_tmatrix_axialsym_fill(t, omega, p->outside.function(omega, p->outside.params), @@ -814,4 +889,3 @@ qpms_errno_t qpms_tmatrix_generator_constant(qpms_tmatrix_t *t, return QPMS_SUCCESS; } - diff --git a/qpms/tmatrices.h b/qpms/tmatrices.h index b39b799..3178c3c 100644 --- a/qpms/tmatrices.h +++ b/qpms/tmatrices.h @@ -488,6 +488,22 @@ qpms_errno_t qpms_tmatrix_generator_axialsym(qpms_tmatrix_t *t, ///< T-matrix to ); +/// Computes the (reduced) transposed R or Q matrix for axially symmetric particle (useful for debugging). +qpms_errno_t qpms_tmatrix_generator_axialsym_RQ_transposed_fill(complex double *target, + complex double omega, + const qpms_tmatrix_generator_axialsym_param_t *param, + qpms_normalisation_t norm, + qpms_bessel_t J + ); + +/// Computes the (reduced) transposed R or Q matrix for axially symmetric particle (useful mostly for debugging). +qpms_errno_t qpms_tmatrix_axialsym_RQ_transposed_fill(complex double *target, + complex double omega, qpms_epsmu_t outside, qpms_epsmu_t inside, + qpms_arc_function_t shape, qpms_l_t lMaxQR, qpms_normalisation_t norm, + qpms_bessel_t J ///< Use QPMS_BESSEL_REGULAR to calculate \f$ R^T\f$ or QPMS_HANKEL_PLUS to calculate \f$ Q^T\f$. + ); + + #if 0 // Abstract types that describe T-matrix/particle/scatsystem symmetries // To be implemented later. See also the thoughts in the beginning of groups.h.