From 4475cffebae6de22e6cd3b97c6cd823a8a9bb148 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 18 Mar 2019 17:54:32 +0200 Subject: [PATCH] C Mie T-matrices but demand yet Bessel functions of complex arguments. Former-commit-id: 101201f5d72ab0f45eefc09c860b65b6841460bf --- qpms/bessel.c | 11 ++++- qpms/qpms_c.pyx | 20 ++++++++- qpms/qpms_cdefs.pxd | 12 ++++++ qpms/qpms_error.h | 3 ++ qpms/qpms_specfunc.h | 2 +- qpms/tmatrices.c | 91 ++++++++++++++++++++++++++++++++++++++-- qpms/tmatrices.h | 61 +++++++++++++++++++++++++++ tests/c_Mie_tmatrices.py | 61 +++++++++++++++++++++++++++ 8 files changed, 254 insertions(+), 7 deletions(-) create mode 100755 tests/c_Mie_tmatrices.py diff --git a/qpms/bessel.c b/qpms/bessel.c index e0600cf..825e30d 100644 --- a/qpms/bessel.c +++ b/qpms/bessel.c @@ -6,6 +6,7 @@ #include "kahansum.h" #include #include +#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; diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index be4a704..2a4cca1 100644 --- a/qpms/qpms_c.pyx +++ b/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 - self.m = np.array(matrix, dtype=complex, copy=True, order='C').reshape((len(spec), len(spec))) + if matrix is None or matrix == 0: + 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[] diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index be67b51..36ff55f 100644 --- a/qpms/qpms_cdefs.pxd +++ b/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: diff --git a/qpms/qpms_error.h b/qpms/qpms_error.h index 798b5d7..325a28a 100644 --- a/qpms/qpms_error.h +++ b/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 diff --git a/qpms/qpms_specfunc.h b/qpms/qpms_specfunc.h index 3a82d67..f95d819 100644 --- a/qpms/qpms_specfunc.h +++ b/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); diff --git a/qpms/tmatrices.c b/qpms/tmatrices.c index 4dfd6a1..6392ed5 100644 --- a/qpms/tmatrices.c +++ b/qpms/tmatrices.c @@ -16,6 +16,7 @@ #include #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; +} + diff --git a/qpms/tmatrices.h b/qpms/tmatrices.h index d8bdcec..3ca849d 100644 --- a/qpms/tmatrices.h +++ b/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 diff --git a/tests/c_Mie_tmatrices.py b/tests/c_Mie_tmatrices.py new file mode 100755 index 0000000..2a0389b --- /dev/null +++ b/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[ ]: + + + +