qpms/cymaterials.pyx

   1 # Cythonized parts of QPMS here
   2 # -----------------------------
   3
   4 import numpy as np
   5 import cmath
   6 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, cdouble
   7 from .cycommon cimport make_c_string
   8 cimport cython
   9 import enum
  10 import warnings
  11 import os
  12 from .constants import e as eV, hbar, c
  13 from libc.stdlib cimport malloc, free, calloc, abort
  14
  15 class EpsMuGeneratorType(enum.Enum):
  16     CONSTANT = 1
  17     PERMITTIVITY_INTERPOLATOR = 2
  18     LORENTZ_DRUDE = 3
  19     PYTHON_CALLABLE = 4
  20
  21 cdef class EpsMu:
  22     """Permittivity and permeability of an isotropic material.
  23
  24     This wraps the C structure qpms_epsmu_t.
  25
  26     See Also
  27     --------
  28     EpsMuGenerator : generates EpsMu objects as a function of frequency.
  29     """
  30     def __init__(self, *args ,**kwargs):
  31         """EpsMu object constructor
  32
  33         Parameters
  34         ----------
  35         eps : complex
  36             Relative electric permittivity.
  37         mu : complex
  38             Relative magnetic permeability.
  39         """
  40         self.em.eps = 1
  41         self.em.mu = 1
  42         if(len(args)>=1):
  43             self.em.eps = args[0]
  44         if(len(args)>=2):
  45             self.em.mu = args[1]
  46         if 'eps' in kwargs.keys():
  47             self.em.eps = kwargs['eps']
  48         if 'mu' in kwargs.keys():
  49             self.em.mu = kwargs['mu']
  50         return
  51
  52     def __repr__(self):
  53         return 'EpsMu(' + repr(self.em.eps) + ', ' + repr(self.em.mu) + ')'
  54
  55     def k(self, omega):
  56         """Wavenumber of the material at frequency omega
  57
  58         Parameters
  59         ----------
  60         omega : complex
  61             Angular frequency in 1/s.
  62
  63         Returns
  64         -------
  65         out : complex
  66             Wavenumber of the material at omega assuming permittivity and permeability
  67             from self.
  68         """
  69         return self.n * omega / c
  70
  71     property n:
  72         """Refractive index of a material specified by this permittivity and permeability."""
  73         def __get__(self):
  74             return (self.em.eps * self.em.mu)**.5
  75     property Z:
  76         """Wave impedance of a material specified by this permittivity and permeability."""
  77         def __get__(self):
  78             return (self.em.mu / self.em.eps)**.5
  79
  80 cdef class LorentzDrudeModel:
  81     """Lorentz-Drude model of material permittivity.
  82
  83     This wraps the C structure qpms_ldparams_t.
  84
  85     Some materials are available in the `lorentz_drude` dictionary.
  86     """
  87
  88     def __cinit__(self, eps_inf, omega_p, f_arr, omega_arr, gamma_arr):
  89         cdef size_t n = len(omega_arr)
  90         if (n != len(gamma_arr) or n != len(f_arr)):
  91             raise ValueError("omega_arr, gamma_arr and f_arr must have equal lengths!")
  92         cdef qpms_ldparams_t *p
  93         p = <qpms_ldparams_t *>malloc(sizeof(qpms_ldparams_t) + sizeof(qpms_ldparams_triple_t) * n)
  94         p[0].eps_inf = eps_inf
  95         p[0].omega_p = omega_p
  96         p[0].n = n
  97         cdef size_t i
  98         for i in range(0,n):
  99             p[0].data[i].f = f_arr[i]
 100             p[0].data[i].omega = omega_arr[i]
 101             p[0].data[i].gamma = gamma_arr[i]
 102         self.params = p
 103
 104     def __dealloc__(self):
 105         free(self.params)
 106         self.params = NULL
 107
 108     def __call__(self, omega):
 109         """Evaluates the permittivity at a given frequency
 110
 111         Parameters
 112         ----------
 113         omega : complex
 114             Angular frequency in 1/s.
 115
 116         Returns
 117         -------
 118         eps : complex
 119             Relative permittivity from the Lorentz-Drude model
 120         """
 121         return qpms_lorentzdrude_eps(omega, self.params)
 122
 123 cdef class _CLorentzDrudeModel:
 124     """Lorentz-Drude model of material permittivity.
 125
 126     This is an auxilliary class making the pre-compiled C Lorentz-Drude models accessible.
 127     For defining own Lorentz-Drude models from python, use LorentzDrudeModel class instead.
 128     """
 129     def __cinit__(self):
 130         "Do not use directly. Do not use from Python. Use the link() method instead."
 131         pass
 132
 133     @staticmethod
 134     cdef link(const qpms_ldparams_t *params):
 135         self = _CLorentzDrudeModel()
 136         self.params = params
 137         return self
 138
 139     def __call__(self, omega):
 140         """Evaluates the permittivity at a given frequency
 141
 142         Parameters
 143         ----------
 144         omega : complex
 145             Angular frequency in 1/s.
 146
 147         Returns
 148         -------
 149         eps : complex
 150             Relative permittivity from the Lorentz-Drude model
 151         """
 152         return qpms_lorentzdrude_eps(omega, self.params)
 153
 154 cdef double eh = eV/hbar
 155
 156 # Some basic Lorentz-Drude parameters
 157 lorentz_drude = {
 158         'Au_py' :  # This should give the same results as 'Au'; to be removed.
 159             LorentzDrudeModel(1, 9.03*eh,
 160                 (0.76, 0.024, 0.01, 0.071, 0.601, 4.384),
 161                 (0, 0.415*eh, 0.83*eh, 2.969*eh, 4.304*eh, 13.32*eh),
 162                 (0.053*eh, 0.241*eh, 0.345*eh, 0.87*eh, 2.494*eh, 2.214*eh)),
 163         'Ag_py' :  # This should give the same results as 'Ag'; to be removed.
 164             LorentzDrudeModel(1, 9.01*eh,
 165                 (0.845, 0.065,0.124, 0.011, 0.840, 5.646),
 166                 (0, 0.816*eh,4.481*eh, 8.185*eh, 9.083*eh, 20.29*eh),
 167                     (0.048*eh, 3.886*eh, 0.452*eh,0.065*eh, 0.916*eh, 2.419*eh)),
 168         'Au' : _CLorentzDrudeModel.link(QPMS_LDPARAMS_AU),
 169         'Ag' : _CLorentzDrudeModel.link(QPMS_LDPARAMS_AG),
 170         'Cu' : _CLorentzDrudeModel.link(QPMS_LDPARAMS_CU),
 171         'Al' : _CLorentzDrudeModel.link(QPMS_LDPARAMS_AL),
 172         'Cr' : _CLorentzDrudeModel.link(QPMS_LDPARAMS_CR),
 173         'Ti' : _CLorentzDrudeModel.link(QPMS_LDPARAMS_TI),
 174         'Be' : _CLorentzDrudeModel.link(QPMS_LDPARAMS_BE),
 175         'Ni' : _CLorentzDrudeModel.link(QPMS_LDPARAMS_NI),
 176         'Pd' : _CLorentzDrudeModel.link(QPMS_LDPARAMS_PD),
 177         'Pt' : _CLorentzDrudeModel.link(QPMS_LDPARAMS_PT),
 178         'W'  : _CLorentzDrudeModel.link(QPMS_LDPARAMS_W),
 179
 180 }
 181
 182 cdef qpms_epsmu_t python_epsmu_generator(cdouble omega, const void *params):
 183     cdef object fun = <object> params
 184     cdef qpms_epsmu_t em
 185     em.eps, em.mu = fun(omega)
 186     return em
 187
 188
 189
 190 cdef class EpsMuGenerator:
 191     def __init__(self, what):
 192         if isinstance(what, str):
 193             if what in lorentz_drude.keys():
 194                 what = lorentz_drude[what]
 195             else:
 196                 try: what = float(what)
 197                 except ValueError:
 198                     try: what = complex(what)
 199                     except ValueError as ve:
 200                         raise ValueError("Trying to create EpsMuGenerator from a string that can't be parsed as a number (constant refractive index) nor supported lorentz-drude key") from ve
 201         if isinstance(what, int):
 202             what = float(what)
 203         if isinstance(what, (float, complex)):
 204             what = EpsMu(what**2)
 205         if isinstance(what, EpsMu):
 206             self.holder = what
 207             self.g.function = qpms_epsmu_const_g
 208             self.g.params = (<EpsMu?>self.holder).rawpointer()
 209         elif isinstance(what, LorentzDrudeModel):
 210             self.holder = what
 211             self.g.function = qpms_lorentzdrude_epsmu_g
 212             self.g.params = (<LorentzDrudeModel?>self.holder).rawpointer()
 213         elif isinstance(what, _CLorentzDrudeModel):
 214             self.holder = what
 215             self.g.function = qpms_lorentzdrude_epsmu_g
 216             self.g.params = (<_CLorentzDrudeModel?>self.holder).rawpointer()
 217         elif isinstance(what, MaterialInterpolator):
 218             self.holder = what
 219             self.g.function = qpms_permittivity_interpolator_epsmu_g
 220             self.g.params = (<MaterialInterpolator?>self.holder).rawpointer()
 221         elif isinstance(what, EpsMuGenerator): # Copy constructor
 222             self.holder = (<EpsMuGenerator?>what).holder
 223             self.g = (<EpsMuGenerator?>what).g
 224         elif callable(what):
 225             warnings.warn("Custom python (eps,mu) generators are an experimental feature")
 226             self.holder = what
 227             self.g.function = python_epsmu_generator
 228             self.g.params = <const void *> what
 229         else:
 230             raise ValueError("Must be constructed from a number (refractive index), EpsMu, LorentzDrudeModel or MaterialInterpolator, or a python callable object that returns an (epsilon, mu) tuple.")
 231
 232     property typ:
 233         def __get__(self):
 234             if(self.g.function == qpms_epsmu_const_g):
 235                 return EpsMuGeneratorType.CONSTANT
 236             elif(self.g.function == qpms_lorentzdrude_epsmu_g):
 237                 return EpsMuGeneratorType.LORENTZ_DRUDE
 238             elif(self.g.function == qpms_permittivity_interpolator_epsmu_g):
 239                 return EpsMuGeneratorType.PERMITTIVITY_INTERPOLATOR
 240             elif(self.g.function == python_epsmu_generator):
 241                 return EpsMuGeneratorType.PYTHON_CALLABLE
 242             else:
 243                 raise ValueError("Damn, this is a bug.")
 244
 245     def __call__(self, omega):
 246         cdef qpms_epsmu_t em
 247         if self.g.function == qpms_permittivity_interpolator_epsmu_g:
 248             i = self.holder.freq_interval
 249             if(omega < i[0] or omega > i[1]):
 250                 raise ValueError("Input frequency %g is outside the interpolator domain (%g, %g)."
 251                     % (omega, i[0], i[1]))
 252         em = self.g.function(omega, self.g.params)
 253         return EpsMu(em.eps, em.mu)
 254
 255     def n(self, omega):
 256         return self(omega).n
 257     def Z(self, omega):
 258         return self(omega).Z
 259     def k(self, omega):
 260         return self(omega).k(omega)
 261
 262     cdef qpms_epsmu_generator_t raw(self):
 263         return self.g
 264
 265
 266 cdef class MaterialInterpolator:
 267     '''
 268     Wrapper over the qpms_permittivity_interpolator_t structure.
 269     '''
 270
 271     def __cinit__(self, filename, *args, **kwargs):
 272         '''Creates a permittivity interpolator.'''
 273         cdef char *cpath = make_c_string(filename)
 274         self.interp = qpms_permittivity_interpolator_from_yml(cpath, gsl_interp_cspline)
 275         if not self.interp:
 276             raise IOError("Could not load permittivity data from %s" % filename)
 277         self.omegamin = qpms_permittivity_interpolator_omega_min(self.interp)
 278         self.omegamax = qpms_permittivity_interpolator_omega_max(self.interp)
 279
 280     def __dealloc__(self):
 281         qpms_permittivity_interpolator_free(self.interp)
 282
 283     def __call__(self, double freq):
 284         '''Returns interpolated permittivity, corresponding to a given angular frequency.'''
 285         if freq < self.omegamin or freq > self.omegamax:
 286             raise ValueError("Input frequency %g is outside the interpolator domain (%g, %g)."
 287                     % (freq, self.minomega, self.freqs[self.maxomega]))
 288         return qpms_permittivity_interpolator_eps_at_omega(self.interp, freq)
 289
 290     property freq_interval:
 291         def __get__(self):
 292             return [self.omegamin, self.omegamax]
 293