diff --git a/notes/cylinderT.lyx b/notes/cylinderT.lyx index 5e4b7f9..a4dc377 100644 --- a/notes/cylinderT.lyx +++ b/notes/cylinderT.lyx @@ -131,9 +131,9 @@ where refers to the particle inside; then \begin_inset Formula -\[ -T_{nn'}=-\sum_{n''}R_{nn''}Q_{n''n}^{-1}. -\] +\begin{equation} +T_{nn'}=-\sum_{n''}R_{nn''}Q_{n''n}^{-1}.\label{eq:T matrix from R and Q} +\end{equation} \end_inset @@ -196,6 +196,52 @@ For fully axially symmetric particles the integrals vanish for \end_inset asimuthal factor in the integrand. + One then has +\begin_inset Formula +\begin{equation} +T_{nn'}=-\sum_{n''}R'_{nn''}Q'_{n''n}^{-1}\label{eq:T-matrix from reduced R and Q} +\end{equation} + +\end_inset + +where +\begin_inset Formula +\begin{align*} +R'_{nn'} & =\int_{0}^{\pi}\left(\frac{\eta}{\eta_{1}}\wfkcreg_{n}\left(k\vect r\right)\times\wfkcreg_{\overline{n'}}\left(k_{1}\vect r\right)+\wfkcreg_{\overline{n}}\left(k\vect r\right)\times\wfkcreg_{n'}\left(k_{1}\vect r\right)\right)\cdot\left(\uvec r\cos\beta\left(\theta\right)+\uvec{\theta}\sin\beta\left(\theta\right)\right)\frac{\left(r\left(\theta\right)\right)^{2}\sin\theta}{\cos\beta\left(\theta\right)}\ud\theta,\\ +Q'_{nn'} & =\int_{0}^{\pi}\left(\frac{\eta}{\eta_{1}}\wfkcreg_{n}\left(k\vect r\right)\times\wfkcreg_{\overline{n'}}\left(k_{1}\vect r\right)+\wfkcreg_{\overline{n}}\left(k\vect r\right)\times\wfkcreg_{n'}\left(k_{1}\vect r\right)\right)\cdot\left(\uvec r\cos\beta\left(\theta\right)+\uvec{\theta}\sin\beta\left(\theta\right)\right)\frac{\left(r\left(\theta\right)\right)^{2}\sin\theta}{\cos\beta\left(\theta\right)}\ud\theta +\end{align*} + +\end_inset + +where +\begin_inset Formula $\vect r=\vect r\left(\theta\right)=\left(r\left(\theta\right),\theta,0\right)$ +\end_inset + +. + Matrices +\begin_inset Formula $Q',R'$ +\end_inset + + differ from the original +\begin_inset Formula $R,Q$ +\end_inset + + matrices in +\begin_inset CommandInset ref +LatexCommand eqref +reference "eq:T matrix from R and Q" +plural "false" +caps "false" +noprefix "false" + +\end_inset + + by a factor of +\begin_inset Formula $2\pi ik^{2}$ +\end_inset + +, but this cancels out in the matrix product. + \end_layout \begin_layout Standard diff --git a/qpms/tmatrices.c b/qpms/tmatrices.c index 8b0be37..8dac691 100644 --- a/qpms/tmatrices.c +++ b/qpms/tmatrices.c @@ -1,4 +1,9 @@ #define _POSIX_C_SOURCE 200809L // for getline() +#define lapack_int int +#define lapack_complex_double complex double +#define lapack_complex_double_real(z) (creal(z)) +#define lapack_complex_double_imag(z) (cimag(z)) +#include #include #include #include @@ -18,12 +23,19 @@ #include "qpms_specfunc.h" #include "normalisation.h" #include +#include + +// These are quite arbitrarily chosen constants for the quadrature in qpms_tmatrix_axialsym_fill() +#define TMATRIX_AXIALSYM_INTEGRAL_EPSREL (1e-6) +#define TMATRIX_AXIALSYM_INTEGRAL_EPSABS (1e-10) #define HBAR (1.05457162825e-34) #define ELECTRONVOLT (1.602176487e-19) #define SCUFF_OMEGAUNIT (3e14) #define SQ(x) ((x)*(x)) +#define MAX(x,y) ((x) < (y) ? (y) : (x)) + qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec) { qpms_tmatrix_t *t = malloc(sizeof(qpms_tmatrix_t)); if (!t) abort(); @@ -612,27 +624,148 @@ qpms_arc_function_retval_t qpms_arc_cylinder(double theta, const void *param) { struct tmatrix_axialsym_integral_param_t { const qpms_vswf_set_spec_t *bspec; - qpms_l_t l1, l2; - qpms_m_t m1; // m2 = -m1 - qpms_vswf_type_t t1; // t2 = 2 - t1 + qpms_l_t l, l_in; + qpms_m_t m; // m_in = -m + qpms_vswf_type_t t, t_in; qpms_arc_function_t f; complex double k_in, k, z_in, z; bool realpart; // Otherwise imaginary part - bool Q; // Otherwise R + qpms_bessel_t btype; // For Q QPMS_HANKEL_PLUS, for R QPMS_BESSEL_REGULAR }; -#if 0 static double tmatrix_axialsym_integrand(double theta, void *param) { + // Pretty inefficient; either real or imaginary part is thrown away etc. struct tmatrix_axialsym_integral_param_t *p = param; + qpms_arc_function_retval_t rb = p->f.function(theta, p->f.params); + csph_t kr = {rb.r * p->k, theta, 0}, + kr_in = {rb.r * p->k_in, theta, 0}; + csphvec_t y_el = qpms_vswf_single_el_csph(p->m, p->l, kr, p->btype, p->bspec->norm); + csphvec_t y_mg = qpms_vswf_single_mg_csph(p->m, p->l, kr, p->btype, p->bspec->norm); + csphvec_t v_in_el = qpms_vswf_single_el_csph(-p->m, p->l_in, kr_in, QPMS_BESSEL_REGULAR, p->bspec->norm); + csphvec_t v_in_mg = qpms_vswf_single_mg_csph(-p->m, p->l_in, kr_in, QPMS_BESSEL_REGULAR, p->bspec->norm); + csphvec_t y1, y2, v_in1, v_in2; + switch(p->t) { + case QPMS_VSWF_ELECTRIC: y1 = y_el; y2 = y_mg; break; + case QPMS_VSWF_MAGNETIC: y1 = y_mg; y2 = y_el; break; + default: QPMS_WTF; + } + switch(p->t_in) { + case QPMS_VSWF_ELECTRIC: v_in1 = v_in_mg; v_in2 = v_in_el; break; + case QPMS_VSWF_MAGNETIC: v_in1 = v_in_el; v_in2 = v_in_mg; break; + default: QPMS_WTF; + } + // Normal vector components + double nrc = cos(rb.beta), nthetac = sin(rb.beta); + // First triple product + complex double tp1 = nrc * (y1.thetac * v_in1.phic - y1.phic * v_in1.thetac) + + nthetac * (y1.phic * v_in1.rc - y1.rc * v_in1.phic); + // Second thiple product + complex double tp2 = nrc * (y2.thetac * v_in2.phic - y2.phic * v_in2.thetac) + + nthetac * (y2.phic * v_in2.rc - y2.rc * v_in2.phic); + double jac = SQ(rb.r) * sin(theta) / nrc; // Jacobian + complex double res = p->z/p->z_in * tp1 + tp2; + return p->realpart ? creal(res) : cimag(res); } qpms_errno_t qpms_tmatrix_axialsym_fill( - qpms_tmatrix_t *t, complex double omega, qpms_epsmu_generator_t outside, - qpms_epsmu_generator_t inside,qpms_arc_function_t shape) + qpms_tmatrix_t *t, complex double omega, qpms_epsmu_t outside, + qpms_epsmu_t inside,qpms_arc_function_t shape, qpms_l_t lMaxQR) { + QPMS_UNTESTED; + const qpms_vswf_set_spec_t *bspec = t->spec; + struct tmatrix_axialsym_integral_param_t p; + 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; + const gsl_function f = {tmatrix_axialsym_integrand, (void *) &p}; + if (lMaxQR < bspec->lMax) lMaxQR = bspec->lMax; + qpms_vswf_set_spec_t *bspecQR = qpms_vswf_set_spec_from_lMax(lMaxQR, bspec->norm); + size_t *reindex = qpms_vswf_set_reindex(bspec, bspecQR); + // Q' and R' matrices + complex double *Q, *R; + QPMS_CRASHING_MALLOC(Q, sizeof(complex double) * SQ(bspecQR->n)); + QPMS_CRASHING_MALLOC(R, sizeof(complex double) * SQ(bspecQR->n)); + + // 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) { + Q[iQR] = 0; + R[iQR] = 0; + } else { + p.m = m1; + p.l = l1; p.l_in = l2; + p.t = t1; p.t_in = t2; + + double result; // We throw the quadrature error estimate away. + // Re(R') + p.btype = QPMS_BESSEL_REGULAR; + p.realpart = true; + QPMS_ENSURE_SUCCESS(gsl_integration_qag(&f, 0, M_PI, epsabs, epsrel, + intlimit, 2, w, &result, NULL)); + R[iQR] = result; + + // Im(R') + p.realpart = false; + QPMS_ENSURE_SUCCESS(gsl_integration_qag(&f, 0, M_PI, epsabs, epsrel, + intlimit, 2, w, &result, NULL)); + R[iQR] += I*result; + + // Re(Q') + p.btype = QPMS_HANKEL_PLUS; + p.realpart = true; + QPMS_ENSURE_SUCCESS(gsl_integration_qag(&f, 0, M_PI, epsabs, epsrel, + intlimit, 2, w, &result, NULL)); + Q[iQR] = result; + + // Im(Q') + p.realpart = false; + QPMS_ENSURE_SUCCESS(gsl_integration_qag(&f, 0, M_PI, epsabs, epsrel, + intlimit, 2, w, &result, NULL)); + Q[iQR] += I*result; + } + } + gsl_integration_workspace_free(w); + + // Compute T-matrix; maybe it would be better solved with some sparse matrix mechanism, + // but fukkit. + const size_t n = bspecQR->n; + lapack_int *pivot; + QPMS_CRASHING_MALLOC(pivot, sizeof(lapack_int) * n); + QPMS_ENSURE_SUCCESS(LAPACKE_zgetrf(LAPACK_ROW_MAJOR, n, n, Q, n, pivot)); + // Solve Q'^T X = R'^T where X will be -T^T + // Note that Q'^T, R'^T are already (transposed) in Q, R. + const complex double minus1 = -1; + QPMS_ENSURE_SUCCESS(LAPACKE_zgetrs(LAPACK_ROW_MAJOR, 'N', n, n /*nrhs*/, + Q, n, pivot, R, n)); + // R now contains -T^T. + for(size_t i1 = 0; i1 < bspec->n; ++i1) + for(size_t i2 = 0; i2 < bspec->n; ++i2) { + if (reindex[i1] == ~(size_t) 0 && reindex[i2] == ~(size_t) 0) QPMS_WTF; + const size_t it = i1 * bspec->n + i2; + const size_t iQR = reindex[i1] + reindex[i2] * bspecQR->n; + t->m[it] = -R[iQR]; + } + + free(pivot); free(R); free(Q); free(reindex); + qpms_vswf_set_spec_free(bspecQR); + return QPMS_SUCCESS; } -#endif diff --git a/qpms/tmatrices.h b/qpms/tmatrices.h index 3bba2e0..1094a22 100644 --- a/qpms/tmatrices.h +++ b/qpms/tmatrices.h @@ -392,21 +392,32 @@ qpms_arc_function_retval_t qpms_arc_sphere(double theta, qpms_errno_t qpms_tmatrix_axialsym_fill( qpms_tmatrix_t *t, ///< T-matrix whose contents are to be replaced. Not NULL. complex double omega, ///< Angular frequency. - qpms_epsmu_generator_t outside, ///< Optical properties of the outside medium. - qpms_epsmu_generator_t inside, ///< Optical properties of the particle's material. - qpms_arc_function_t shape ///< Particle surface parametrisation. + qpms_epsmu_t outside, ///< Optical properties of the outside medium. + qpms_epsmu_t inside, ///< Optical properties of the particle's material. + qpms_arc_function_t shape, ///< Particle surface parametrisation. + /** If `lMax_extend > t->bspec->lMax`, then the internal \a Q, \a R matrices will be + * trimmed at this larger lMax; the final T-matrix will then be trimmed + * according to bspec. + */ + qpms_l_t lMax_extend ); /// Creates a new T-matrix of a particle with \f$ C_\infty \f$ symmetry. static inline qpms_tmatrix_t *qpms_tmatrix_axialsym( const qpms_vswf_set_spec_t *bspec, complex double omega, ///< Angular frequency. - qpms_epsmu_generator_t outside, ///< Optical properties of the outside medium. - qpms_epsmu_generator_t inside, ///< Optical properties of the particle's material. - qpms_arc_function_t shape ///< Particle surface parametrisation. + qpms_epsmu_t outside, ///< Optical properties of the outside medium. + qpms_epsmu_t inside, ///< Optical properties of the particle's material. + qpms_arc_function_t shape, ///< Particle surface parametrisation. + /// Internal lMax to improve precision of the result. + /** If `lMax_extend > bspec->lMax`, then the internal \a Q, \a R matrices will be + * trimmed at this larger lMax; the final T-matrix will then be trimmed + * according to bspec. + */ + qpms_l_t lMax_extend ) { qpms_tmatrix_t *t = qpms_tmatrix_init(bspec); - qpms_tmatrix_axialsym_fill(t, omega, outside, inside, shape); + qpms_tmatrix_axialsym_fill(t, omega, outside, inside, shape, lMax_extend); return t; } diff --git a/qpms/vswf.c b/qpms/vswf.c index 2775962..c319091 100644 --- a/qpms/vswf.c +++ b/qpms/vswf.c @@ -95,7 +95,52 @@ void qpms_vswf_set_spec_free(qpms_vswf_set_spec_t *s) { free(s); } -csphvec_t qpms_vswf_single_el(qpms_m_t m, qpms_l_t l, sph_t kdlj, +struct bspec_reindex_pair { + qpms_uvswfi_t ui; + size_t i_orig; +}; + +static int cmp_bspec_reindex_pair(const void *aa, const void *bb) { + const struct bspec_reindex_pair *a = aa, *b = bb; + if (a->ui < b->ui) return -1; + else if (a->ui == b->ui) return 0; + else return 1; +} + +size_t *qpms_vswf_set_reindex(const qpms_vswf_set_spec_t *small, const qpms_vswf_set_spec_t *big) { + QPMS_UNTESTED; + struct bspec_reindex_pair *small_pairs, *big_pairs; + size_t *r; + QPMS_CRASHING_MALLOC(small_pairs, sizeof(struct bspec_reindex_pair) * small->n); + QPMS_CRASHING_MALLOC(big_pairs, sizeof(struct bspec_reindex_pair) * big->n); + QPMS_CRASHING_MALLOC(r, sizeof(size_t) * small->n); + for(size_t i = 0; i < small->n; ++i) { + small_pairs[i].ui = small->ilist[i]; + small_pairs[i].i_orig = i; + } + for(size_t i = 0 ; i < big->n; ++i) { + big_pairs[i].ui = big->ilist[i]; + big_pairs[i].i_orig = i; + } + qsort(small_pairs, small->n, sizeof(struct bspec_reindex_pair), cmp_bspec_reindex_pair); + qsort(big_pairs, big->n, sizeof(struct bspec_reindex_pair), cmp_bspec_reindex_pair); + + size_t bi = 0; + for(size_t si = 0; si < small->n; ++si) { + while(big_pairs[bi].ui < small_pairs[si].ui) + ++bi; + if(big_pairs[bi].ui == small_pairs[si].ui) + r[small_pairs[si].i_orig] = big_pairs[si].i_orig; + else + r[small_pairs[si].i_orig] = ~(size_t)0; + } + + free(small_pairs); + free(big_pairs); + return r; +} + +csphvec_t qpms_vswf_single_el_csph(qpms_m_t m, qpms_l_t l, csph_t kdlj, qpms_bessel_t btyp, qpms_normalisation_t norm) { lmcheck(l,m); csphvec_t N; @@ -120,7 +165,7 @@ csphvec_t qpms_vswf_single_el(qpms_m_t m, qpms_l_t l, sph_t kdlj, 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_csph(qpms_m_t m, qpms_l_t l, csph_t kdlj, qpms_bessel_t btyp, qpms_normalisation_t norm) { lmcheck(l,m); csphvec_t M; @@ -143,6 +188,16 @@ csphvec_t qpms_vswf_single_mg(qpms_m_t m, qpms_l_t l, sph_t kdlj, return M; } +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) { + return qpms_vswf_single_el_csph(m, l, sph2csph(kdlj), btyp, norm); +} + +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) { + return qpms_vswf_single_mg_csph(m, l, sph2csph(kdlj), btyp, norm); +} + qpms_vswfset_sph_t *qpms_vswfset_make(qpms_l_t lMax, sph_t kdlj, qpms_bessel_t btyp, qpms_normalisation_t norm) { qpms_vswfset_sph_t *res = malloc(sizeof(qpms_vswfset_sph_t)); diff --git a/qpms/vswf.h b/qpms/vswf.h index b4a6f21..334b001 100644 --- a/qpms/vswf.h +++ b/qpms/vswf.h @@ -23,7 +23,7 @@ typedef qpms_errno_t (*qpms_incfield_t)( bool add ///< If true, add to target; rewrite target if false. ); -// ---------------Methods for qpms_vswf_spec_t----------------------- +// ---------------Methods for qpms_vswf_set_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); @@ -55,6 +55,19 @@ static inline ssize_t qpms_vswf_set_spec_find_uvswfi(const qpms_vswf_set_spec_t return -1; } +/// Creates an index mapping between two bspecs. +/** + * Creates an array r such that small->ilist[i] == big->ilist[r[i]]. + * It's not lossless if the two bspecs contain different combinations of waves. + * + * Preferably, big->ilist contains everything small->ilist does. + * If small->ilist[i] is not found in big->ilist, r[i] will be set to ~(size_t)0. + * + * Discard with free() after use. + */ +size_t *qpms_vswf_set_reindex(const qpms_vswf_set_spec_t *small, const qpms_vswf_set_spec_t *big); + + /// Evaluates a set of VSWF basis functions at a given point. /** The list of basis wave indices is specified in \a setspec; * \a setspec->norm must be set as well. @@ -120,6 +133,12 @@ csphvec_t qpms_vswf_single_el(int m, int n, sph_t kdlj, /// Magnetic wave M. csphvec_t qpms_vswf_single_mg(int m, int n, sph_t kdlj, qpms_bessel_t btyp, qpms_normalisation_t norm); +/// Electric wave N, complex wave number version. +csphvec_t qpms_vswf_single_el_csph(int m, int n, csph_t kdlj, + qpms_bessel_t btyp, qpms_normalisation_t norm); +/// Magnetic wave M, complex wave number version.. +csphvec_t qpms_vswf_single_mg_csph(int m, int n, csph_t kdlj, + qpms_bessel_t btyp, qpms_normalisation_t norm); /// Set of electric and magnetic VSWF values in spherical coordinate basis. /** This is supposed to contain all the waves up to $l = lMax$.