Lorentz-Drude permittivity model (+LD params for Au, Ag)

Former-commit-id: 76878ee0ec2dd9c7db06febe4266f5537e5d9376

qpms/CMakeLists.txt

@@ -12,7 +12,7 @@ include_directories(${DIRS})
 
 add_library (qpms translations.c tmatrices.c vecprint.c vswf.c wigner.c
 	lattices2d.c gaunt.c error.c legendre.c symmetries.c vecprint.c
-	bessel.c own_zgemm.c parsing.c scatsystem.c materials.c)
+	bessel.c own_zgemm.c parsing.c scatsystem.c materials.c drudeparam_data.c)
 use_c99()
 
 set(LIBS ${LIBS} ${GSL_LIBRARIES} ${GSLCBLAS_LIBRARIES})

qpms/drudeparam_data.c

@@ -0,0 +1,39 @@
+#include "materials.h"
+#include <gsl/gsl_const_mksa.h>
+
+// Data from drudelorentz.com
+
+#define EH (GSL_CONST_MKSA_ELECTRON_VOLT / GSL_CONST_MKSA_PLANCKS_CONSTANT_HBAR)
+
+static const qpms_ldparams_t LDPARAMS_AU = {
+  1, // eps_inf
+  9.03*EH, // omega_p
+  6, // n
+  {
+    {0.76 *EH,  0.   *EH, 0.053*EH},
+    {0.024*EH,  0.415*EH, 0.241*EH},
+    {0.01 *EH,  0.83 *EH, 0.345*EH},
+    {0.071*EH,  2.969*EH, 0.87 *EH},
+    {0.601*EH,  4.304*EH, 2.494*EH},
+    {4.384*EH, 13.32 *EH, 2.214*EH}
+  }
+};
+
+const qpms_ldparams_t *const QPMS_LDPARAMS_AU = &LDPARAMS_AU;
+
+static const qpms_ldparams_t LDPARAMS_AG = {
+  1, // eps_inf
+  9.01*EH, // omega_p
+  6, // n
+  {
+    {0.84 *EH,  0.   *EH, 0.053*EH},
+    {0.065*EH,  0.816*EH, 3.886*EH},
+    {0.124*EH,  4.481*EH, 0.452*EH},
+    {0.111*EH,  8.185*EH, 0.065*EH},
+    {0.84 *EH,  9.083*EH, 0.916*EH},
+    {5.646*EH, 20.29 *EH, 2.419*EH}
+  }
+};
+
+const qpms_ldparams_t *const QPMS_LDPARAMS_AG = &LDPARAMS_AG;
+

qpms/materials.c

@@ -6,7 +6,7 @@
 #include "materials.h"
 #include "qpms_error.h"
 
-#define SPEED_OF_LIGHT (2.99792458e8)
+#define SQ(x) ((x)*(x))
 
 qpms_permittivity_interpolator_t *qpms_permittivity_interpolator_create(
     const size_t incount, const double *wavelen_m, const double *n, const double *k,
@@ -140,3 +140,19 @@ double qpms_permittivity_interpolator_omega_min(
   return 2*M_PI*SPEED_OF_LIGHT / ip->wavelength_m[ip->size-1];
 }
 
+complex double qpms_lorentzdrude_eps(const qpms_ldparams_t *p, complex double omega) {
+  complex double eps = 0;
+  for(size_t j = 0; j < p->n; ++j) {
+    const qpms_ldparams_triple_t d = p->data[j];
+    eps += d.f * SQ(p->omega_p) / (SQ(d.omega) - SQ(omega) + I*omega*d.gamma );
+  }
+  return eps;
+}
+
+qpms_epsmu_t qpms_lorentzdrude_epsmu(const qpms_ldparams_t *p, complex double omega) {
+  qpms_epsmu_t em;
+  em.eps = qpms_lorentzdrude_eps(p, omega);
+  em.mu = 1;
+  return em;
+}
+

qpms/materials.h

@@ -6,6 +6,60 @@
 #include "qpms_types.h"
 #include <gsl/gsl_spline.h>
 
+#ifndef SPEED_OF_LIGHT
+/// Speed of light in m/s.
+#define SPEED_OF_LIGHT (2.99792458e8)
+#endif
+
+
+/// Gets refractive index of a material from its permeability and permittivity.
+/** \f[ n = \sqrt{\mu_r \varepsilon_r} \f] */
+static inline complex double qpms_refindex(qpms_epsmu_t em) {
+	return csqrt(em.eps * em.mu);
+}
+
+/// Gets wave number \a k from angular frequency and material permeability and permittivity.
+/** \f[ k = \frac{n\omega}{c_0}  = \frac{\omega\sqrt{\mu_r \varepsilon_r}}{c_0} \f] */
+static inline complex double qpms_wavenumber(complex double omega, qpms_epsmu_t em) {
+	return qpms_refindex(em)*omega/SPEED_OF_LIGHT;
+}
+
+/// Gets (relative) wave impedance \f$ \eta_r \f$ from material permeability and permittivity.
+/** \eta_r = \sqrt{\mu_r  / \varepsilon_r} \f] */
+static inline complex double qpms_waveimpedance(qpms_epsmu_t em) {
+	return csqrt(em.mu / em.eps);
+}
+
+/// A \f$ (f_j, \omega_j, \gamma_j) \f for qpms_ldparams_t.
+typedef struct qpms_ldparams_triple_t {
+	double f;
+	double omega;
+	double gamma;
+} qpms_ldparams_triple_t;
+
+/// Structure holding Lorentz-Drude model parameters of a material.
+/** \f[
+ * \varepsilon = \varepsilon_\infty + \sum_j=0^{n-1}
+ *    \frac{f_j \omega_p^2}{\omega_j^2-\omega^2+i\omega\gamma_j} 
+ *  \f]
+ */
+typedef struct qpms_ldparams_t {
+	complex double eps_inf; ///< Permittivity at infinity.
+	double omega_p; ///< Plasma frequency.
+	size_t n; ///< Number of "oscillators".
+	qpms_ldparams_triple_t data[]; ///< "Oscillator" parameters.
+} qpms_ldparams_t;
+
+extern const qpms_ldparams_t *const QPMS_LDPARAMS_AG; ///< Lorentz-Drude parameters for silver.
+extern const qpms_ldparams_t *const QPMS_LDPARAMS_AU; ///< Lorentz-Drude parameters for gold.
+
+/// Lorentz-Drude permittivity.
+complex double qpms_lorentzdrude_eps(const qpms_ldparams_t *, complex double omega);
+
+/// Lorentz-Drude optical properties, with relative permeability set always to one.
+qpms_epsmu_t qpms_lorentzdrude_epsmu(const qpms_ldparams_t *, complex double omega);
+
+
 /// Interpolator of tabulated optical properties.
 // TODO use gsl_interp instead of gsl_spline.
 typedef struct qpms_permittivity_interpolator_t {

qpms/tmatrices.c

@@ -21,7 +21,6 @@
 
 #define HBAR (1.05457162825e-34)
 #define ELECTRONVOLT (1.602176487e-19)
-#define SPEED_OF_LIGHT (2.99792458e8)
 #define SCUFF_OMEGAUNIT (3e14)
 
 #define SQ(x) ((x)*(x))