qpms/cytmatrices.pyx

   1 import numpy as np
   2 from .qpms_cdefs cimport *
   3 from .cybspec cimport BaseSpec
   4 from .cycommon import *
   5 from .cycommon cimport make_c_string
   6 from .cymaterials cimport EpsMuGenerator
   7 from .qpms_c cimport FinitePointGroup
   8 import warnings
   9 import os
  10 from libc.stdlib cimport free
  11
  12 cdef class TMatrixInterpolator:
  13     '''
  14     Wrapper over the qpms_tmatrix_interpolator_t structure.
  15     '''
  16     def __cinit__(self, filename, BaseSpec bspec,  *args, **kwargs):
  17         '''Creates a T-matrix interpolator object from a scuff-tmatrix output'''
  18         global qpms_load_scuff_tmatrix_crash_on_failure
  19         qpms_load_scuff_tmatrix_crash_on_failure = False
  20         self.spec = bspec
  21         cdef char * cpath = make_c_string(filename)
  22         retval = qpms_load_scuff_tmatrix(cpath, self.spec.rawpointer(),
  23                 &(self.nfreqs), &(self.freqs), &(self.freqs_su),
  24                 &(self.tmatrices_array), &(self.tmdata))
  25         if (retval != QPMS_SUCCESS):
  26             raise IOError("Could not read T-matrix from %s: %s" % (filename, os.strerror(retval)))
  27         if 'symmetrise' in kwargs:
  28             sym = kwargs['symmetrise']
  29             if isinstance(sym, FinitePointGroup):
  30                 if QPMS_SUCCESS != qpms_symmetrise_tmdata_finite_group(
  31                         self.tmdata, self.nfreqs, self.spec.rawpointer(),
  32                         (<FinitePointGroup?>sym).rawpointer()):
  33                     raise Exception("This should not happen.")
  34                 atol = kwargs['atol'] if 'atol' in kwargs else 1e-16
  35                 qpms_czero_roundoff_clean(self.tmdata, self.nfreqs * len(bspec)**2, atol)
  36             else:
  37                 warnings.warn('symmetrise argument type not supported; ignoring.')
  38         self.interp = qpms_tmatrix_interpolator_create(self.nfreqs,
  39                 self.freqs, self.tmatrices_array, gsl_interp_cspline)
  40         if not self.interp: raise Exception("Unexpected NULL at interpolator creation.")
  41     def __call__(self, double freq):
  42         '''Returns a TMatrix instance, corresponding to a given frequency.'''
  43         if freq < self.freqs[0] or freq > self.freqs[self.nfreqs-1]:# FIXME here I assume that the input is already sorted
  44             raise ValueError("input frequency %g is outside the interpolator domain (%g, %g)"
  45                     % (freq, self.freqs[0], self.freqs[self.nfreqs-1]))
  46         # This is a bit stupid, I should rethink the CTMatrix constuctors
  47         cdef qpms_tmatrix_t *t = qpms_tmatrix_interpolator_eval(self.interp, freq)
  48         cdef CTMatrix res = CTMatrix(self.spec, <cdouble[:len(self.spec),:len(self.spec)]>(t[0].m))
  49         qpms_tmatrix_free(t)
  50         return res
  51     def __dealloc__(self):
  52         qpms_tmatrix_interpolator_free(self.interp)
  53         free(self.tmatrices_array)
  54         free(self.tmdata)
  55         free(self.freqs_su)
  56         free(self.freqs)
  57     property freq_interval:
  58         def __get__(self):
  59             return [self.freqs[0], self.freqs[self.nfreqs-1]]
  60     property omega_table:
  61         def __get__(self):
  62             cdef size_t i
  63             omegas = np.empty((self.nfreqs,), dtype=float)
  64             cdef double[:] omegas_view = omegas
  65             for i in range(self.nfreqs):
  66                 omegas_view[i] = self.freqs[i]
  67             return omegas
  68
  69
  70 cdef class CTMatrix: # N.B. there is another type called TMatrix in tmatrices.py!
  71     '''
  72     Wrapper over the C qpms_tmatrix_t stucture.
  73     '''
  74
  75     def __cinit__(CTMatrix self, BaseSpec spec, matrix):
  76         self.spec = spec
  77         self.t.spec = self.spec.rawpointer();
  78         if (matrix is None) or not np.any(matrix):
  79             self.m = np.zeros((len(spec),len(spec)), dtype=complex, order='C')
  80         else:
  81             # The following will raise an exception if shape is wrong
  82             self.m = np.array(matrix, dtype=complex, copy=True, order='C').reshape((len(spec), len(spec)))
  83         #self.m.setflags(write=False) # checkme
  84         cdef cdouble[:,:] m_memview = self.m
  85         self.t.m = &(m_memview[0,0])
  86         self.t.owns_m = False # Memory in self.t.m is "owned" by self.m, not by self.t...
  87
  88     property rawpointer:
  89         def __get__(self):
  90             return <uintptr_t> &(self.t)
  91
  92     # Transparent access to the T-matrix elements.
  93     def __getitem__(self, key):
  94         return self.m[key]
  95     def __setitem__(self, key, value):
  96         self.m[key] = value
  97
  98     def as_ndarray(CTMatrix self):
  99         ''' Returns a copy of the T-matrix as a numpy array.'''
 100         # Maybe not totally needed after all, as np.array(T[...]) should be equivalent and not longer
 101         return np.array(self.m, copy=True)
 102
 103     def spherical_fill(CTMatrix self, double radius, cdouble k_int,
 104             cdouble k_ext, cdouble mu_int = 1, cdouble mu_ext = 1):
 105         ''' Replaces the contents of the T-matrix with those of a spherical particle.'''
 106         qpms_tmatrix_spherical_fill(&self.t, radius, k_int, k_ext, mu_int, mu_ext)
 107
 108     def spherical_perm_fill(CTMatrix self, double radius, double freq, cdouble epsilon_int,
 109             cdouble epsilon_ext):
 110         '''Replaces the contenst of the T-matrix with those of a spherical particle.'''
 111         qpms_tmatrix_spherical_mu0_fill(&self.t, radius, freq, epsilon_int, epsilon_ext)
 112
 113     @staticmethod
 114     def spherical(BaseSpec spec, double radius, cdouble k_int, cdouble k_ext,
 115             cdouble mu_int = 1, cdouble mu_ext = 1):
 116         ''' Creates a T-matrix of a spherical nanoparticle. '''
 117         tm = CTMatrix(spec, 0)
 118         tm.spherical_fill(radius, k_int, k_ext, mu_int, mu_ext)
 119         return tm
 120
 121     @staticmethod
 122     def spherical_perm(BaseSpec spec, double radius, double freq, cdouble epsilon_int, cdouble epsilon_ext):
 123         '''Creates a T-matrix of a spherical nanoparticle.'''
 124         tm = CTMatrix(spec, 0)
 125         tm.spherical_perm_fill(radius, freq, epsilon_int, epsilon_ext)
 126         return tm
 127
 128 cdef class __MieParams:
 129     # Not directly callable right now, serves just to be used by TMatrixGenerator.
 130     cdef qpms_tmatrix_generator_sphere_param_t cparam
 131     cdef EpsMuGenerator outside
 132     cdef EpsMuGenerator inside
 133     cdef inline  void *rawpointer(self):
 134         return <void *>&(self.cparam)
 135
 136     def __init__(self, outside, inside, r):
 137         self.inside = inside
 138         self.outside = outside
 139         self.cparam.inside = self.inside.raw();
 140         self.cparam.outside = self.outside.raw();
 141         self.cparam.radius = r
 142
 143     property r:
 144         def __get__(self):
 145             return self.cparam.radius
 146         def __set__(self, val):
 147             self.cparam.radius = val
 148
 149 cdef class __ArcCylinder:
 150     cdef qpms_arc_cylinder_params_t p
 151     cdef inline void *rawpointer(self):
 152         return <void *> &(self.p)
 153     def __init__(self, R, h):
 154         self.p.R = R
 155         self.p.h = h
 156
 157 cdef class __ArcSphere:
 158     cdef double r
 159     cdef inline void *rawpointer(self):
 160         return <void *> &(self.r)
 161     def __init__(self, r):
 162         self.r = r
 163
 164 cdef qpms_arc_function_retval_t userarc(double theta, const void *params):
 165     cdef object fun = <object> params
 166     cdef qpms_arc_function_retval_t retval
 167     retval.r, retval.beta = fun(theta)
 168     return retval
 169
 170
 171 cdef class ArcFunction:
 172     cdef qpms_arc_function_t g
 173     cdef object holder
 174     def __init__(self, what):
 175         if isinstance(what, __ArcCylinder):
 176             self.holder = what
 177             self.g.function = qpms_arc_cylinder
 178             self.g.params = (<__ArcCylinder?>self.holder).rawpointer()
 179         elif isinstance(what, __ArcSphere):
 180             self.holder = what
 181             self.g.function = qpms_arc_sphere
 182             self.g.params = (<__ArcSphere?>self.holder).rawpointer()
 183         elif callable(what):
 184             warnings.warn("Custom python (r, beta) arc functions are an experimental feature. Also expect it to be slow.")
 185             self.holder = what
 186             self.g.function = userarc
 187             self.g.params = <const void *> self.holder
 188         elif isinstance(what, ArcFunction): #Copy constructor
 189             self.holder = what.holder
 190             self.g = (<ArcFunction?>what).g
 191             self.holder.rawpointer()
 192
 193 cdef class __AxialSymParams:
 194     cdef qpms_tmatrix_generator_axialsym_param_t p
 195     cdef EpsMuGenerator outside
 196     cdef EpsMuGenerator inside
 197     cdef ArcFunction shape
 198     cdef void * rawpointer(self):
 199         return <void *> &(self.p)
 200     property lMax_extend:
 201         def __get__(self):
 202             return self.p.lMax_extend
 203         def __set__(self, val):
 204             self.p.lMax_extend = val
 205     def __init__(self, outside, inside, shape, *args, **kwargs):
 206         self.outside = outside
 207         self.p.outside = self.outside.g
 208         self.inside = inside
 209         self.p.inside = self.inside.g
 210         self.shape = shape
 211         self.p.shape = self.shape.g
 212         if len(args)>0:
 213             self.lMax_extend = args[0]
 214         if 'lMax_extend' in kwargs.keys():
 215             self.lMax_extend = kwargs['lMax_extend']
 216         if self.lMax_extend == 0:
 217             self.lMax_extend = 1
 218     def Q_transposed(self, cdouble omega, norm):
 219         cdef size_t n = 2*(self.p.lMax_extend*(self.p.lMax_extend +2))
 220         cdef np.ndarray[np.complex_t, ndim=2] arr = np.empty((n,n), dtype=complex, order='C')
 221         cdef cdouble[:,::1] arrview = arr
 222         qpms_tmatrix_generator_axialsym_RQ_transposed_fill(&arrview[0][0], omega, &self.p, norm, QPMS_HANKEL_PLUS)
 223         return arr
 224     def R_transposed(self, cdouble omega, norm):
 225         cdef size_t n = 2*(self.p.lMax_extend*(self.p.lMax_extend +2))
 226         cdef np.ndarray[np.complex_t, ndim=2] arr = np.empty((n,n), dtype=complex, order='C')
 227         cdef cdouble[:,::1] arrview = arr
 228         qpms_tmatrix_generator_axialsym_RQ_transposed_fill(&arrview[0][0], omega, &self.p, norm, QPMS_BESSEL_REGULAR)
 229         return arr
 230
 231 cdef class TMatrixFunction:
 232     '''
 233     Wrapper over qpms_tmatrix_function_t. The main functional difference between this
 234     and TMatrixGenerator is that this remembers a specific BaseSpec
 235     and its __call__ method takes only one mandatory argument (in addition to self).
 236     '''
 237     def __init__(self, TMatrixGenerator tmg, BaseSpec spec):
 238         self.generator = tmg
 239         self.spec = spec
 240         self.f.gen = self.generator.rawpointer()
 241         self.f.spec = self.spec.rawpointer()
 242
 243     def __call__(self, cdouble omega, fill = None):
 244         cdef CTMatrix tm
 245         if fill is None: # make a new CTMatrix
 246             tm = CTMatrix(self.spec, None)
 247         else: # TODO check whether fill has the same bspec as self?
 248             tm = fill
 249         if self.f.gen.function(tm.rawpointer(), omega, self.f.gen.params) != 0:
 250             raise ValueError("Something went wrong")
 251         else:
 252             return tm
 253
 254
 255 cdef class TMatrixGenerator:
 256     def __init__(self, what):
 257         if isinstance(what, __MieParams):
 258             self.holder = what
 259             self.g.function = qpms_tmatrix_generator_sphere
 260             self.g.params = (<__MieParams?>self.holder).rawpointer()
 261         elif isinstance(what,__AxialSymParams):
 262             self.holder = what
 263             self.g.function = qpms_tmatrix_generator_axialsym
 264             self.g.params = (<__AxialSymParams?>self.holder).rawpointer()
 265         elif isinstance(what, CTMatrix):
 266             self.holder = what
 267             self.g.function = qpms_tmatrix_generator_constant
 268             self.g.params = <void*>(<CTMatrix?>self.holder).rawpointer()
 269         elif isinstance(what, TMatrixInterpolator):
 270             self.holder = what
 271             self.g.function = qpms_tmatrix_generator_interpolator
 272             self.g.params = <void*>(<TMatrixInterpolator?>self.holder).rawpointer()
 273         else:
 274             raise TypeError("Can't construct TMatrixGenerator from that")
 275
 276     def __call__(self, arg, cdouble omega):
 277         cdef CTMatrix tm
 278         if isinstance(arg, CTMatrix): # fill the matrix
 279             tm = arg
 280             if self.g.function(tm.rawpointer(), omega, self.g.params) != 0:
 281                 raise ValueError("Something went wrong")
 282             return
 283         elif isinstance(arg, BaseSpec): # Make a new CTMatrix
 284             tm = CTMatrix(arg, None)
 285             if self.g.function(tm.rawpointer(), omega, self.g.params) != 0:
 286                 raise ValueError("Something went wrong")
 287             return tm
 288         else:
 289             raise ValueError("Must specify CTMatrix or BaseSpec")
 290
 291     def Q_transposed(self, cdouble omega, norm):
 292         if self.g.function != qpms_tmatrix_generator_axialsym:
 293             raise TypeError("Applicable only for axialsym generators")
 294         return self.holder.Q_transposed(omega, norm)
 295     def R_transposed(self, cdouble omega, norm):
 296         if self.g.function != qpms_tmatrix_generator_axialsym:
 297             raise TypeError("Applicable only for axialsym generators")
 298         return self.holder.R_transposed(omega, norm)
 299
 300     # Better "constructors":
 301     @staticmethod
 302     def sphere(outside, inside, r):
 303         return TMatrixGenerator(__MieParams(EpsMuGenerator(outside),
 304                     EpsMuGenerator(inside), r))
 305     @staticmethod
 306     def sphere_asarc(outside, inside, r, *args, **kwargs):
 307         return TMatrixGenerator(__AxialSymParams(
 308             EpsMuGenerator(outside), EpsMuGenerator(inside),
 309             ArcFunction(__ArcSphere(r)), *args, **kwargs))
 310     @staticmethod
 311     def cylinder(outside, inside, r, h, *args, **kwargs):
 312         return TMatrixGenerator(__AxialSymParams(
 313             EpsMuGenerator(outside), EpsMuGenerator(inside),
 314             ArcFunction(__ArcCylinder(r, h)), *args, **kwargs))
 315