T-matrix for axially symmetric particles; untested

Former-commit-id: dfc7a0771a52df097a2afba172bb286369f085a2

notes/cylinderT.lyx

@@ -131,9 +131,9 @@ where
 
  refers to the particle inside; then
 \begin_inset Formula 
-\[
+\begin{equation}
-T_{nn'}=-\sum_{n''}R_{nn''}Q_{n''n}^{-1}.
+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

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 <lapacke.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <cblas.h>
@@ -18,12 +23,19 @@
 #include "qpms_specfunc.h"
 #include "normalisation.h"
 #include <errno.h>
+#include <gsl/gsl_integration.h>
+
+// 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_l_t l, l_in; 
-  qpms_m_t m1; // m2 = -m1
+  qpms_m_t m; // m_in = -m
-  qpms_vswf_type_t t1; // t2 = 2 - t1
+  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_tmatrix_t *t, complex double omega, qpms_epsmu_t outside,
-		qpms_epsmu_generator_t inside,qpms_arc_function_t shape)
+		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

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_t outside, ///< Optical properties of the outside medium.
-		qpms_epsmu_generator_t inside, ///< Optical properties of the particle's material.
+		qpms_epsmu_t inside, ///< Optical properties of the particle's material.
-		qpms_arc_function_t shape ///< Particle surface parametrisation.
+		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_t outside, ///< Optical properties of the outside medium.
-		qpms_epsmu_generator_t inside, ///< Optical properties of the particle's material.
+		qpms_epsmu_t inside, ///< Optical properties of the particle's material.
-		qpms_arc_function_t shape ///< Particle surface parametrisation.
+		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;
 }
 

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));

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$.