From c4ed9852a6180ad87c5bdcb006e4fb61b958581e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Thu, 8 Aug 2019 17:16:07 +0300 Subject: [PATCH] Lorentz-Drude permittivity model (+LD params for Au, Ag) Former-commit-id: 76878ee0ec2dd9c7db06febe4266f5537e5d9376 --- qpms/CMakeLists.txt | 2 +- qpms/drudeparam_data.c | 39 ++++++++++++++++++++++++++++++ qpms/materials.c | 18 +++++++++++++- qpms/materials.h | 54 ++++++++++++++++++++++++++++++++++++++++++ qpms/tmatrices.c | 1 - 5 files changed, 111 insertions(+), 3 deletions(-) create mode 100644 qpms/drudeparam_data.c diff --git a/qpms/CMakeLists.txt b/qpms/CMakeLists.txt index 259eafd..b530588 100644 --- a/qpms/CMakeLists.txt +++ b/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}) diff --git a/qpms/drudeparam_data.c b/qpms/drudeparam_data.c new file mode 100644 index 0000000..1516dd7 --- /dev/null +++ b/qpms/drudeparam_data.c @@ -0,0 +1,39 @@ +#include "materials.h" +#include + +// 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; + diff --git a/qpms/materials.c b/qpms/materials.c index 0be5467..8905fc9 100644 --- a/qpms/materials.c +++ b/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; +} + diff --git a/qpms/materials.h b/qpms/materials.h index 55cc012..8fecb36 100644 --- a/qpms/materials.h +++ b/qpms/materials.h @@ -6,6 +6,60 @@ #include "qpms_types.h" #include +#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 { diff --git a/qpms/tmatrices.c b/qpms/tmatrices.c index 8925193..48fdb4d 100644 --- a/qpms/tmatrices.c +++ b/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))