diff --git a/qpms/__init__.py b/qpms/__init__.py index fb7026d..4595639 100644 --- a/qpms/__init__.py +++ b/qpms/__init__.py @@ -7,7 +7,7 @@ from .cyquaternions import CQuat, IRot3 from .cybspec import VSWFNorm, BaseSpec from .cytmatrices import CTMatrix, TMatrixInterpolator from .cytranslations import trans_calculator -from .cymaterials import MaterialInterpolator +from .cymaterials import MaterialInterpolator, EpsMu, LorentzDrudeModel, lorentz_drude, EpsMuGenerator from .cycommon import dbgmsg_enable, dbgmsg_disable, dbgmsg_active, BesselType from .lattices2d import * from .hexpoints import * diff --git a/qpms/cymaterials.pxd b/qpms/cymaterials.pxd index 3a295b6..873adf6 100644 --- a/qpms/cymaterials.pxd +++ b/qpms/cymaterials.pxd @@ -1,8 +1,25 @@ -from .qpms_cdefs cimport qpms_permittivity_interpolator_t +from .qpms_cdefs cimport qpms_permittivity_interpolator_t, qpms_epsmu_generator_t, qpms_epsmu_t, qpms_ldparams_t cdef class MaterialInterpolator: cdef qpms_permittivity_interpolator_t *interp cdef readonly double omegamin cdef readonly double omegamax + cdef inline void *rawpointer(self): + return &(self.interp) +cdef class EpsMu: + cdef public qpms_epsmu_t em + cdef inline void *rawpointer(self): + return &(self.em) + +cdef class LorentzDrudeModel: + cdef const qpms_ldparams_t *params + #cdef bint owns_params + cdef inline void *rawpointer(self): + return &(self.params) + +cdef class EpsMuGenerator: + cdef qpms_epsmu_generator_t g + cdef object holder + cdef qpms_epsmu_generator_t raw(self) diff --git a/qpms/cymaterials.pyx b/qpms/cymaterials.pyx index 26c0af4..2769f5e 100644 --- a/qpms/cymaterials.pyx +++ b/qpms/cymaterials.pyx @@ -3,14 +3,118 @@ import numpy as np import cmath -from .qpms_cdefs cimport qpms_permittivity_interpolator_from_yml, qpms_permittivity_interpolator_free, qpms_permittivity_interpolator_omega_min, qpms_permittivity_interpolator_omega_max, gsl_interp_type, qpms_permittivity_interpolator_t, gsl_interp_cspline, qpms_permittivity_interpolator_eps_at_omega +from .qpms_cdefs cimport qpms_permittivity_interpolator_from_yml, qpms_permittivity_interpolator_free, qpms_permittivity_interpolator_omega_min, qpms_permittivity_interpolator_omega_max, gsl_interp_type, qpms_permittivity_interpolator_t, gsl_interp_cspline, qpms_permittivity_interpolator_eps_at_omega, qpms_epsmu_const_g, qpms_permittivity_interpolator_epsmu_g, qpms_epsmu_const_g, qpms_lorentzdrude_epsmu_g, qpms_ldparams_triple_t, qpms_lorentzdrude_eps from .cycommon cimport make_c_string cimport cython import enum import warnings import os +from scipy.constants import e as eV, hbar from libc.stdlib cimport malloc, free, calloc, abort +class EpsMuGeneratorType(enum.IntEnum): + CONSTANT = 1 + PERMITTIVITY_INTERPOLATOR = 2 + LORENTZ_DRUDE = 3 + +cdef class EpsMu: + def __init__(self, *args ,**kwargs): + self.em.eps = 1 + self.em.mu = 1 + if(len(args)>=1): + self.em.eps = args[0] + if(len(args)>=2): + self.em.mu = args[1] + if 'eps' in kwargs.keys(): + self.em.eps = kwargs['eps'] + if 'mu' in kwargs.keys(): + self.em.mu = kwargs['mu'] + return + + def __repr__(self): + return 'EpsMu(' + repr(self.em.eps) + ', ' + repr(self.em.mu) + ')' + +cdef class LorentzDrudeModel: + def __cinit__(self, eps_inf, omega_p, f_arr, omega_arr, gamma_arr): + cdef size_t n = len(omega_arr) + if (n != len(gamma_arr) or n != len(f_arr)): + raise ValueError("omega_arr, gamma_arr and f_arr must have equal lengths!") + cdef qpms_ldparams_t *p + p = malloc(sizeof(qpms_ldparams_t) + sizeof(qpms_ldparams_triple_t) * n) + p[0].eps_inf = eps_inf + p[0].omega_p = omega_p + p[0].n = n + cdef size_t i + for i in range(0,n): + p[0].data[i].f = f_arr[i] + p[0].data[i].omega = omega_arr[i] + p[0].data[i].gamma = gamma_arr[i] + self.params = p + + def __dealloc__(self): + free(self.params) + self.params = NULL + + def __call__(self, omega): + return qpms_lorentzdrude_eps(omega, self.params) + +cdef double eh = eV/hbar + +# Some basic Lorentz-Drude parameters +lorentz_drude = { + 'Au' : + LorentzDrudeModel(1, 9.03*eh, + (0.76, 0.024, 0.01, 0.071, 0.601, 4.384), + (0, 0.415*eh, 0.83*eh, 2.969*eh, 4.304*eh, 13.32*eh), + (0.053*eh, 0.241*eh, 0.345*eh, 0.87*eh, 2.494*eh, 2.214*eh)), + 'Ag' : + LorentzDrudeModel(1, 9.01*eh, + (0.84, 0.065,0.124, 0.111, 0.840, 5.646), + (0, 0.816*eh,4.481*eh, 8.185*eh, 9.083*eh, 20.29*eh), + (0.053*eh, 3.886*eh, 0.452*eh,0.065*eh, 0.916*eh, 2.419*eh)), +} + +cdef class EpsMuGenerator: + def __init__(self, what): + if isinstance(what, EpsMu): + self.holder = what + self.g.function = qpms_epsmu_const_g + self.g.params = (self.holder).rawpointer() + elif isinstance(what, LorentzDrudeModel): + self.holder = what + self.g.function = qpms_lorentzdrude_epsmu_g + self.g.params = (self.holder).rawpointer() + elif isinstance(what, MaterialInterpolator): + self.holder = what + self.g.function = qpms_permittivity_interpolator_epsmu_g + self.g.params = (self.holder).rawpointer() + else: + raise ValueError("Must be constructed from EpsMu, LorentzDrudeModel or MaterialInterpolator") + + property typ: + def __get__(self): + if(self.g.function == qpms_epsmu_const_g): + return EpsMuGeneratorType.CONSTANT + elif(self.g.function == qpms_lorentzdrude_epsmu_g): + return EpsMuGeneratorType.LORENTZ_DRUDE + elif(self.g.function == qpms_permittivity_interpolator_epsmu_g): + return EpsMuGeneratorType.PERMITTIVITY_INTERPOLATOR + else: + raise ValueError("Damn, this is a bug.") + + def __call__(self, omega): + cdef qpms_epsmu_t em + if self.g.function == qpms_permittivity_interpolator_epsmu_g: + i = self.holder.freq_interval + if(omega < i[0] or omega > i[1]): + raise ValueError("Input frequency %g is outside the interpolator domain (%g, %g)." + % (omega, i[0], i[1])) + em = self.g.function(omega, self.g.params) + return EpsMu(em.eps, em.mu) + + cdef qpms_epsmu_generator_t raw(self): + return self.g + cdef class MaterialInterpolator: ''' diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 9c9d010..f43f45f 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -338,10 +338,18 @@ cdef extern from "materials.h": double qpms_permittivity_interpolator_omega_max(const qpms_permittivity_interpolator_t *interp) double qpms_permittivity_interpolator_omega_min(const qpms_permittivity_interpolator_t *interp) void qpms_permittivity_interpolator_free(qpms_permittivity_interpolator_t *interp) + struct qpms_ldparams_triple_t: + double f + double omega + double gamma struct qpms_ldparams_t: - pass + cdouble eps_inf + double omega_p + size_t n + qpms_ldparams_triple_t data[0] const qpms_ldparams_t *const QPMS_LDPARAMS_AG const qpms_ldparams_t *const QPMS_LDPARAMS_AU + cdouble qpms_lorentzdrude_eps(cdouble, const qpms_ldparams_t *) cdef extern from "tmatrices.h": bint qpms_load_scuff_tmatrix_crash_on_failure