From 8683c69c0e104e1df0aeff289e098603d1cdadd2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Sun, 24 Mar 2019 14:47:42 +0000
Subject: [PATCH] Cython wrapper over qpms_permittivity_interpolator_t.

Former-commit-id: 38b362b52a8e8f5b6defca688b0589b25a2dfe22
---
 qpms/qpms_c.pyx     | 31 +++++++++++++++++++++++++++++++
 qpms/qpms_cdefs.pxd | 18 ++++++++++++++++++
 2 files changed, 49 insertions(+)

diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx
index 2a4cca1..ff1f3f4 100644
--- a/qpms/qpms_c.pyx
+++ b/qpms/qpms_c.pyx
@@ -1065,6 +1065,37 @@ cdef class IRot3:
         qpms_irot3_uvswfi_dense(&view[0,0], bspec.rawpointer(), self.qd)
         return m
 
+cdef class MaterialInterpolator:
+    '''
+    Wrapper over the qpms_permittivity_interpolator_t structure.
+    '''
+    cdef qpms_permittivity_interpolator_t *interp
+    cdef readonly double omegamin
+    cdef readonly double omegamax
+
+    def __cinit__(self, filename, *args, **kwargs):
+        '''Creates a permittivity interpolator.'''
+        cdef char *cpath = make_c_string(filename)
+        self.interp = qpms_permittivity_interpolator_from_yml(cpath, gsl_interp_cspline)
+        if not self.interp: 
+            raise IOError("Could not load permittivity data from %s" % filename)
+        self.omegamin = qpms_permittivity_interpolator_omega_min(self.interp)
+        self.omegamax = qpms_permittivity_interpolator_omega_max(self.interp)
+
+    def __dealloc__(self):
+        qpms_permittivity_interpolator_free(self.interp)
+
+    def __call__(self, double freq):
+        '''Returns interpolated permittivity, corresponding to a given angular frequency.'''
+        if freq < self.omegamin or freq > self.omegamax:
+            raise ValueError("Input frequency %g is outside the interpolator domain (%g, %g)."
+                    % (freq, self.minomega, self.freqs[self.maxomega]))
+        return qpms_permittivity_interpolator_eps_at_omega(self.interp, freq)
+
+    property freq_interval:
+        def __get__(self):
+            return [self.omegamin, self.omegamax]
+
 cdef class TMatrixInterpolator:
     '''
     Wrapper over the qpms_tmatrix_interpolator_t structure.
diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd
index 36ff55f..5116dd7 100644
--- a/qpms/qpms_cdefs.pxd
+++ b/qpms/qpms_cdefs.pxd
@@ -234,6 +234,8 @@ cdef extern from "gsl/gsl_interp.h":
 cdef extern from "tmatrices.h":
     struct qpms_tmatrix_interpolator_t:
         const qpms_vswf_set_spec_t *bspec
+    struct qpms_permittivity_interpolator_t:
+        pass
     void qpms_tmatrix_interpolator_free(qpms_tmatrix_interpolator_t *interp)
     qpms_tmatrix_t *qpms_tmatrix_interpolator_eval(const qpms_tmatrix_interpolator_t *interp, double freq)
     qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, double *freqs, 
@@ -266,10 +268,26 @@ cdef extern from "tmatrices.h":
             cdouble **tmdata)
     cdouble *qpms_mie_coefficients_reflection(cdouble *target, const qpms_vswf_set_spec_t *bspec,
             double a, cdouble k_i, cdouble k_e, cdouble mu_i, cdouble mu_e, qpms_bessel_t J_ext, qpms_bessel_t J_scat)
+    qpms_tmatrix_t *qpms_tmatrix_spherical(const qpms_vswf_set_spec_t *bspec, double a, 
+            cdouble k_i, cdouble k_e, cdouble mu_i, cdouble mu_e)
     qpms_errno_t qpms_tmatrix_spherical_fill(qpms_tmatrix_t *t, double a, 
             cdouble k_i, cdouble k_e, cdouble mu_i, cdouble mu_e)
     qpms_tmatrix_t *qpms_tmatrix_spherical(const qpms_vswf_set_spec_t *bspec,
             double a, cdouble k_i, cdouble k_e, cdouble mu_i, cdouble mu_e)
+    cdouble qpms_drude_epsilon(cdouble eps_inf, cdouble omega_p, cdouble gamma_p, cdouble omega)
+    qpms_errno_t qpms_tmatrix_spherical_mu0_fill(qpms_tmatrix_t *t, double a, double omega,
+            cdouble epsilon_fg, cdouble epsilon_bg)
+    qpms_tmatrix_t *qpms_tmatrix_spherical_mu0(const qpms_vswf_set_spec_t *bspec, double a,
+            double omega, cdouble epsilon_fg, cdouble epsilon_bg)
+    qpms_permittivity_interpolator_t *qpms_permittivity_interpolator_create(const size_t incount,
+            cdouble *wavelength_m, cdouble *n, cdouble *k, const gsl_interp_type *iptype)
+    qpms_permittivity_interpolator_t *qpms_permittivity_interpolator_from_yml(const char *path,
+            const gsl_interp_type *iptype)
+    cdouble qpms_permittivity_interpolator_eps_at_omega(const qpms_permittivity_interpolator_t *interp, double omega_SI)
+    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)
+
 
 cdef extern from "scatsystem.h":
     struct qpms_particle_t: