C Mie T-matrices but demand yet Bessel functions of complex arguments.

Former-commit-id: 101201f5d72ab0f45eefc09c860b65b6841460bf

qpms/bessel.c

@@ -6,6 +6,7 @@
 #include "kahansum.h"
 #include <gsl/gsl_sf_bessel.h>
 #include <complex.h>
+#include "qpms_error.h"
 
 #ifndef M_LN2
 #define M_LN2   0.69314718055994530942  /* log_e 2 */
@@ -16,7 +17,7 @@ static inline complex double ipow(int x) {
 }
 
 // There is a big issue with gsl's precision of spherical bessel function; these have to be implemented differently
-qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, double x, complex double *result_array) {
+qpms_errno_t qpms_sph_bessel_realx_fill(qpms_bessel_t typ, qpms_l_t lmax, double x, complex double *result_array) {
   int retval;
   double tmparr[lmax+1];
   switch(typ) {
@@ -49,6 +50,14 @@ qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, double x, co
   assert(0);
 }
 
+qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, complex double x, complex double *result_array) {
+  int retval = qpms_sph_bessel_realx_fill(typ, lmax, creal(x), result_array);
+  if(!cimag(x))
+    return retval;
+  else
+    QPMS_NOT_IMPLEMENTED("complex argument of Bessel function.");
+}
+
 static inline ptrdiff_t akn_index(qpms_l_t n, qpms_l_t k) {
   assert(k <= n);
   return ((ptrdiff_t) n + 1) * n / 2 + k;

qpms/qpms_c.pyx

@@ -1131,8 +1131,11 @@ cdef class CTMatrix: # N.B. there is another type called TMatrix in tmatrices.py
     def __cinit__(CTMatrix self, BaseSpec spec, matrix):
         self.spec = spec
         self.t.spec = self.spec.rawpointer();
-        # The following will raise an exception if shape is wrong
+        if matrix is None or matrix == 0:
-        self.m = np.array(matrix, dtype=complex, copy=True, order='C').reshape((len(spec), len(spec)))
+            self.m = np.zeros((len(spec),len(spec)), dtype=complex, order='C')
+        else:
+            # The following will raise an exception if shape is wrong
+            self.m = np.array(matrix, dtype=complex, copy=True, order='C').reshape((len(spec), len(spec)))
         #self.m.setflags(write=False) # checkme
         cdef cdouble[:,:] m_memview = self.m
         self.t.m = &(m_memview[0,0])
@@ -1158,6 +1161,19 @@ cdef class CTMatrix: # N.B. there is another type called TMatrix in tmatrices.py
         # Maybe not totally needed after all, as np.array(T[...]) should be equivalent and not longer
         return np.array(self.m, copy=True)
 
+    def spherical_fill(CTMatrix self, double radius, cdouble k_int,
+            cdouble k_ext, cdouble mu_int = 1, cdouble mu_ext = 1):
+        ''' Replaces the contents of the T-matrix with those of a spherical particle. '''
+        qpms_tmatrix_spherical_fill(&self.t, radius, k_int, k_ext, mu_int, mu_ext)
+    
+    @staticmethod
+    def spherical(BaseSpec spec, double radius, cdouble k_int, cdouble k_ext, 
+            cdouble mu_int = 1, cdouble mu_ext = 1):
+        ''' Creates a T-matrix of a spherical nanoparticle. '''
+        tm = CTMatrix(spec, 0)
+        tm.spherical_fill(radius, k_int, k_ext, mu_int, mu_ext)
+        return tm
+
 cdef char *make_c_string(pythonstring):
     '''
     Copies contents of a python string into a char[]

qpms/qpms_cdefs.pxd

@@ -38,6 +38,12 @@ cdef extern from "qpms_types.h":
         QPMS_NORMALISATION_SPHARM
         QPMS_NORMALISATION_SPHARM_CS
         QPMS_NORMALISATION_UNDEF
+    ctypedef enum qpms_bessel_t:
+        QPMS_BESSEL_REGULAR
+        QPMS_BESSEL_SINGULAR
+        QPMS_HANKEL_PLUS
+        QPMS_HANKEL_MINUS
+        QPMS_HANKEL_UNDEF
     ctypedef int qpms_lm_t
     ctypedef int qpms_l_t
     ctypedef int qpms_m_t
@@ -258,6 +264,12 @@ cdef extern from "tmatrices.h":
     qpms_errno_t qpms_load_scuff_tmatrix(const char *path, const qpms_vswf_set_spec_t *bspec,
             size_t *n, double **freqs, double **freqs_su, qpms_tmatrix_t **tmatrices_array,
             cdouble **tmdata)
+    cdouble *qpms_mie_coefficients_reflection(cdouble *target, const qpms_vswf_set_spec_t *bspec,
+            double a, cdouble k_i, cdouble k_e, cdouble mu_i, cdouble mu_e, qpms_bessel_t J_ext, qpms_bessel_t J_scat)
+    qpms_errno_t qpms_tmatrix_spherical_fill(qpms_tmatrix_t *t, double a, 
+            cdouble k_i, cdouble k_e, cdouble mu_i, cdouble mu_e)
+    qpms_tmatrix_t *qpms_tmatrix_spherical(const qpms_vswf_set_spec_t *bspec,
+            double a, cdouble k_i, cdouble k_e, cdouble mu_i, cdouble mu_e)
 
 cdef extern from "scatsystem.h":
     struct qpms_particle_t:

qpms/qpms_error.h

@@ -29,4 +29,7 @@ void qpms_pr_debug_at_flf(const char *filename, unsigned int linenum,
 
 #define QPMS_ENSURE(x, msg, ...) {if(!(x)) qpms_pr_error_at_flf(__FILE__,__LINE__,__func__,msg, ##__VA_ARGS__); }
 
+#define QPMS_NOT_IMPLEMENTED(msg, ...) qpms_pr_error_at_flf(__FILE__,__LINE__,__func__, \
+		"Not implemented:" msg, ##__VA_ARGS__)
+
 #endif 

qpms/qpms_specfunc.h

@@ -8,7 +8,7 @@
  ******************************************************************************/
 
 // TODO unify types
-qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, double x, complex double *result_array);
+qpms_errno_t qpms_sph_bessel_fill(qpms_bessel_t typ, qpms_l_t lmax, complex double x, complex double *result_array);
 
 
 

qpms/tmatrices.c

@@ -16,6 +16,7 @@
 #include <string.h>
 #include "qpms_error.h"
 #include "tmatrices.h"
+#include "qpms_specfunc.h"
 
 #define HBAR (1.05457162825e-34)
 #define ELECTRONVOLT (1.602176487e-19)
@@ -314,16 +315,21 @@ void qpms_tmatrix_interpolator_free(qpms_tmatrix_interpolator_t *ip) {
 
 qpms_tmatrix_t *qpms_tmatrix_interpolator_eval(const qpms_tmatrix_interpolator_t *ip, double freq) {
   qpms_tmatrix_t *t = qpms_tmatrix_init(ip->bspec);
+  QPMS_ENSURE_SUCCESS(qpms_tmatrix_interpolator_eval_fill(t, ip, freq));
+  return t;
+}
+
+qpms_errno_t qpms_tmatrix_interpolator_eval_fill(qpms_tmatrix_t *t,
+    const qpms_tmatrix_interpolator_t *ip, double freq) {
+  QPMS_ENSURE(qpms_vswf_set_spec_isidentical(t->spec, ip->bspec), "Tried to fill a T-matrix with an incompatible interpolator!");
   const size_t n = ip->bspec->n;
   for (size_t i = 0; i < n*n; ++i){
     if (ip->splines_real[i]) t->m[i] = gsl_spline_eval(ip->splines_real[i], freq, NULL /*does this work?*/);
     if (ip->splines_imag[i]) t->m[i] += I* gsl_spline_eval(ip->splines_imag[i], freq, NULL /*does this work?*/);
   }
-  return t;
+  return QPMS_SUCCESS;
 }
 
-
-
 double qpms_SU2eV(double e_SU) {
   return e_SU * SCUFF_OMEGAUNIT / (ELECTRONVOLT / HBAR);
 }
@@ -454,3 +460,82 @@ qpms_errno_t qpms_load_scuff_tmatrix(
   return retval;
 }
 
+complex double *qpms_mie_coefficients_reflection(
+		complex double *target, ///< Target array of length bspec->n. If NULL, a new one will be allocated.
+		const qpms_vswf_set_spec_t *bspec, ///< Defines which of the coefficients are calculated.
+		double a, ///< Radius of the sphere.
+		complex double k_i, ///< Wave number of the internal material of the sphere.
+		complex double k_e, ///< Wave number of the surrounding medium.
+		complex double mu_i, ///< Relative permeability of the sphere material.
+		complex double mu_e, ///< Relative permeability of the surrounding medium.
+    qpms_bessel_t J_ext,
+    qpms_bessel_t J_scat
+		// TODO J_ext, J_scat?
+		) {
+  /*
+   *  This implementation pretty much copies mie_coefficients() from qpms_p.py, so 
+   *  any bugs there should affect this function as well and perhaps vice versa.
+   */
+  QPMS_ENSURE(J_ext != J_scat, "J_ext and J_scat must not be equal. Perhaps you want J_ext = 1 and J_scat = 3 anyways.");
+  if (!target) QPMS_CRASHING_MALLOC(target, bspec->n * sizeof(complex double));
+  qpms_l_t lMax = bspec->lMax;
+  complex double x_i = k_i * a;
+  complex double x_e = k_e * a;
+  complex double m = k_i / k_e;
+  complex double eta_inv_i = k_i / mu_i;
+  complex double eta_inv_e = k_e / mu_e;
+
+  complex double zi[lMax + 2];
+  complex double ze[lMax + 2];
+  complex double zs[lMax + 2];
+  complex double RH[lMax + 1] /* MAGNETIC */, RV[lMax+1] /* ELECTRIC */;
+  QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(QPMS_BESSEL_REGULAR, lMax+1, x_i, zi));
+  QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(J_ext, lMax+1, x_e, ze));
+  QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(J_scat, lMax+1, x_e, zs));
+  for (qpms_l_t l = 0; l <= lMax; ++l) {
+    // Bessel function derivatives as in DLMF 10.51.2
+    complex double dzi = -zi[l+1] + l/x_i*zi[l];
+    complex double dze = -ze[l+1] + l/x_e*ze[l];
+    complex double dzs = -zs[l+1] + l/x_e*zs[l];
+    complex double fi = zi[l]/x_i+dzi;
+    complex double fs = zs[l]/x_e+dzs;
+    complex double fe = ze[l]/x_e+dze;
+    RV[l] = -((-eta_inv_i * fe * zi[l] + eta_inv_e * ze[l] * fi)/(-eta_inv_e * fi * zs[l] + eta_inv_i * zi[l] * fs));
+    RH[l] = -((-eta_inv_e * fe * zi[l] + eta_inv_i * ze[l] * fi)/(-eta_inv_i * fi * zs[l] + eta_inv_e * zi[l] * fs));
+  }
+
+  for (size_t i = 0; i < bspec->n; ++i) {
+    qpms_l_t l; qpms_m_t m; qpms_vswf_type_t t;
+    QPMS_ENSURE_SUCCESS(qpms_uvswfi2tmn(bspec->ilist[i], &t, &m, &l));
+    assert(l <= lMax);
+    switch(t) {
+      case QPMS_VSWF_ELECTRIC:
+        target[i] = RV[l];
+        break;
+      case QPMS_VSWF_MAGNETIC:
+        target[i] = RH[l];
+        break;
+      default: QPMS_WTF;
+    }
+  }
+  return target;
+}
+
+/// Replaces the contents of an existing T-matrix with that of a spherical nanoparticle calculated using the Lorentz-mie theory.
+qpms_errno_t qpms_tmatrix_spherical_fill(qpms_tmatrix_t *t, ///< T-matrix whose contents are to be replaced. Not NULL.
+		double a, ///< Radius of the sphere.
+		complex double k_i, ///< Wave number of the internal material of the sphere.
+		complex double k_e, ///< Wave number of the surrounding medium.
+		complex double mu_i, ///< Relative permeability of the sphere material.
+		complex double mu_e ///< Relative permeability of the surrounding medium.
+		) {
+  qpms_l_t lMax = t->spec->lMax;
+  complex double *miecoeffs = qpms_mie_coefficients_reflection(NULL, t->spec, a, k_i, k_e, mu_i, mu_e,
+    QPMS_BESSEL_REGULAR, QPMS_HANKEL_PLUS);
+  memset(t->m, 0, SQ(t->spec->n));
+  for(size_t i = 0; i < t->spec->n; ++i) 
+    t->m[t->spec->n*i + i] = miecoeffs[i];
+  free(miecoeffs);
+  return QPMS_SUCCESS;
+}
+

qpms/tmatrices.h

@@ -233,6 +233,10 @@ typedef struct qpms_tmatrix_interpolator_t {
 /// Free a T-matrix interpolator.
 void qpms_tmatrix_interpolator_free(qpms_tmatrix_interpolator_t *interp);
 
+/// Fills an existing T-matrix with new interpolated values.
+qpms_errno_t qpms_tmatrix_interpolator_eval_fill(qpms_tmatrix_t *target, ///< T-matrix to be updated, not NULL.
+		const qpms_tmatrix_interpolator_t *interp, double freq);
+
 /// Evaluate a T-matrix interpolated value.
 /** The result is to be freed using qpms_tmatrix_free().*/
 qpms_tmatrix_t *qpms_tmatrix_interpolator_eval(const qpms_tmatrix_interpolator_t *interp, double freq);
@@ -244,6 +248,63 @@ qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, ///< Num
 	       //, bool copy_bspec ///< if true, copies its own copy of basis spec from the first T-matrix.
        	       /*, ...? */);
 
+/// Calculates the reflection Mie-Lorentz coefficients for a spherical particle.
+/**
+ * This function is based on the previous python implementation mie_coefficients() from qpms_p.py,
+ * so any bugs therein should affect this function as well and perhaps vice versa.
+ *
+ * Most importantly, the case of magnetic material, \a mu_i != 0 or \a mu_e != 0 has never been tested
+ * and might give wrong results.
+ *
+ * \return Array with the Mie-Lorentz reflection coefficients in the order determined by bspec.
+ * If \a target was not NULL, this is target, otherwise a newly allocated array.
+ *
+ * TODO better doc.
+ */
+complex double *qpms_mie_coefficients_reflection(
+		complex double *target, ///< Target array of length bspec->n. If NULL, a new one will be allocated.
+		const qpms_vswf_set_spec_t *bspec, ///< Defines which of the coefficients are calculated.
+		double a, ///< Radius of the sphere.
+		complex double k_i, ///< Wave number of the internal material of the sphere.
+		complex double k_e, ///< Wave number of the surrounding medium.
+		complex double mu_i, ///< Relative permeability of the sphere material.
+		complex double mu_e, ///< Relative permeability of the surrounding medium.
+		qpms_bessel_t J_ext, ///< Kind of the "incoming" waves. Most likely QPMS_BESSEL_REGULAR.
+		qpms_bessel_t J_scat ///< Kind of the "scattered" waves. Most likely QPMS_HANKEL_PLUS.
+		);
+
+/// Replaces the contents of an existing T-matrix with that of a spherical nanoparticle calculated using the Lorentz-mie theory.
+qpms_errno_t qpms_tmatrix_spherical_fill(qpms_tmatrix_t *t, ///< T-matrix whose contents are to be replaced. Not NULL.
+		double a, ///< Radius of the sphere.
+		complex double k_i, ///< Wave number of the internal material of the sphere.
+		complex double k_e, ///< Wave number of the surrounding medium.
+		complex double mu_i, ///< Relative permeability of the sphere material.
+		complex double mu_e ///< Relative permeability of the surrounding medium.
+		);
+
+/// Creates a new T-matrix of a spherical particle using the Lorentz-Mie theory.
+static inline qpms_tmatrix_t *qpms_tmatrix_spherical(
+		const qpms_vswf_set_spec_t *bspec,
+		double a, ///< Radius of the sphere.
+		complex double k_i, ///< Wave number of the internal material of the sphere.
+		complex double k_e, ///< Wave number of the surrounding medium.
+		complex double mu_i, ///< Relative permeability of the sphere material.
+		complex double mu_e ///< Relative permeability of the surrounding medium.
+		) {
+	qpms_tmatrix_t *t = qpms_tmatrix_init(bspec);
+	qpms_tmatrix_spherical_fill(t, a, k_i, k_e, mu_i, mu_e);
+	return t;
+}
+
+/// Relative permittivity from the Drude model.
+static inline complex double qpms_drude_epsilon(
+		complex double eps_inf, ///< Permittivity at 
+		complex double omega_p, ///< Plasma frequency \f$ \omega_p \f$ of the material.
+		complex double gamma_p, ///< Decay constant \f$ \gamma_p \f$ of the material.
+		complex double omega ///< Frequency \f$ \omega \f$ at which the permittivity is evaluated.
+		) {
+	return eps_inf - omega_p*omega_p/(omega*(omega+I*gamma_p));
+}
 
 #if 0
 // Abstract types that describe T-matrix/particle/scatsystem symmetries

tests/c_Mie_tmatrices.py

@@ -0,0 +1,61 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+# In[1]:
+
+
+from qpms import *
+
+
+# In[6]:
+
+
+R = 40e-9
+ω_p = 9*eV/ℏ #9*eV/ℏ
+ε_inf = 4.6
+γ_p = 0.1*eV/ℏ
+ε_b = 2.13
+lMax = 3
+
+ω = 1.5*eV/ℏ
+
+
+# In[7]:
+
+
+ε_m = ε_drude(ε_inf, ω_p, γ_p, ω)
+
+
+# In[9]:
+
+
+k_i = cmath.sqrt(ε_m)*ω/c
+k_e = cmath.sqrt(ε_b)*ω/c
+RH, RV, TH, TV = mie_coefficients(a=R, nmax=lMax, k_i=k_i, k_e=k_e, J_ext=1, J_scat=3)
+
+
+# In[11]:
+
+
+spec = BaseSpec(lMax=lMax)
+cT = CTMatrix.spherical(spec, R, k_i, k_e, 1, 1)
+
+
+# In[16]:
+
+
+print(np.diag(cT[...]))
+
+
+# In[18]:
+
+
+print(RV)
+print(RH)
+
+
+# In[ ]:
+
+
+
+