Material generator cython wrappers (partly broken)

Former-commit-id: 881a182939bae3cccb545d51f1d138997bdf9c5c

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 *

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 <void *>&(self.interp)
 
+cdef class EpsMu:
+    cdef public qpms_epsmu_t em
+    cdef inline void *rawpointer(self):
+        return <void *>&(self.em)
+
+cdef class LorentzDrudeModel:
+    cdef const qpms_ldparams_t *params
+    #cdef bint owns_params
+    cdef inline void *rawpointer(self):
+        return <void *>&(self.params)
+
+cdef class EpsMuGenerator:
+    cdef qpms_epsmu_generator_t g
+    cdef object holder
+    cdef qpms_epsmu_generator_t raw(self)
 

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 = <qpms_ldparams_t *>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 = (<EpsMu?>self.holder).rawpointer()
+        elif isinstance(what, LorentzDrudeModel):
+            self.holder = what
+            self.g.function = qpms_lorentzdrude_epsmu_g
+            self.g.params = (<LorentzDrudeModel?>self.holder).rawpointer()
+        elif isinstance(what, MaterialInterpolator):
+            self.holder = what
+            self.g.function = qpms_permittivity_interpolator_epsmu_g
+            self.g.params = (<MaterialInterpolator?>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:
     '''

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