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, copy=True):
  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=copy, 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 cdef qpms_errno_t qpms_tmatrix_generator_pythoncallable(
 171         qpms_tmatrix_t *t, cdouble omega, const void *param):
 172     '''Use a python callable as a T-matrix generator
 173
 174     fun is expected to be callable as
 175     fun(CTMatrix tmatrix, cdouble omega)
 176     Therefore we must recreate the CTMatrix object
 177     and its BaseSpec object before passing it to fun
 178     '''
 179     cdef object fun = <object> param
 180     cdef size_t n = t[0].spec[0].n
 181     cdef BaseSpec bspec = BaseSpec.from_cpointer(t[0].spec)
 182     cdef CTMatrix tmatrix = CTMatrix(bspec,
 183             np.asarray(<np.complex[:n, :n]>(t[0].m)), copy=False)
 184     return fun(tmatrix, omega)
 185
 186
 187 cdef class ArcFunction:
 188     cdef qpms_arc_function_t g
 189     cdef object holder
 190     def __init__(self, what):
 191         if isinstance(what, __ArcCylinder):
 192             self.holder = what
 193             self.g.function = qpms_arc_cylinder
 194             self.g.params = (<__ArcCylinder?>self.holder).rawpointer()
 195         elif isinstance(what, __ArcSphere):
 196             self.holder = what
 197             self.g.function = qpms_arc_sphere
 198             self.g.params = (<__ArcSphere?>self.holder).rawpointer()
 199         elif callable(what):
 200             warnings.warn("Custom python (r, beta) arc functions are an experimental feature. Also expect it to be slow.")
 201             self.holder = what
 202             self.g.function = userarc
 203             self.g.params = <const void *> self.holder
 204         elif isinstance(what, ArcFunction): #Copy constructor
 205             self.holder = what.holder
 206             self.g = (<ArcFunction?>what).g
 207             self.holder.rawpointer()
 208
 209 cdef class __AxialSymParams:
 210     cdef qpms_tmatrix_generator_axialsym_param_t p
 211     cdef EpsMuGenerator outside
 212     cdef EpsMuGenerator inside
 213     cdef ArcFunction shape
 214     cdef void * rawpointer(self):
 215         return <void *> &(self.p)
 216     property lMax_extend:
 217         def __get__(self):
 218             return self.p.lMax_extend
 219         def __set__(self, val):
 220             self.p.lMax_extend = val
 221     def __init__(self, outside, inside, shape, *args, **kwargs):
 222         self.outside = outside
 223         self.p.outside = self.outside.g
 224         self.inside = inside
 225         self.p.inside = self.inside.g
 226         self.shape = shape
 227         self.p.shape = self.shape.g
 228         if len(args)>0:
 229             self.lMax_extend = args[0]
 230         if 'lMax_extend' in kwargs.keys() and kwargs['lMax_extend'] is not None:
 231             self.lMax_extend = kwargs['lMax_extend']
 232         if self.lMax_extend == 0:
 233             self.lMax_extend = 1
 234     def Q_transposed(self, cdouble omega, norm):
 235         cdef size_t n = 2*(self.p.lMax_extend*(self.p.lMax_extend +2))
 236         cdef np.ndarray[np.complex_t, ndim=2] arr = np.empty((n,n), dtype=complex, order='C')
 237         cdef cdouble[:,::1] arrview = arr
 238         qpms_tmatrix_generator_axialsym_RQ_transposed_fill(&arrview[0][0], omega, &self.p, norm, QPMS_HANKEL_PLUS)
 239         return arr
 240     def R_transposed(self, cdouble omega, norm):
 241         cdef size_t n = 2*(self.p.lMax_extend*(self.p.lMax_extend +2))
 242         cdef np.ndarray[np.complex_t, ndim=2] arr = np.empty((n,n), dtype=complex, order='C')
 243         cdef cdouble[:,::1] arrview = arr
 244         qpms_tmatrix_generator_axialsym_RQ_transposed_fill(&arrview[0][0], omega, &self.p, norm, QPMS_BESSEL_REGULAR)
 245         return arr
 246
 247 cdef class TMatrixFunction:
 248     '''
 249     Wrapper over qpms_tmatrix_function_t. The main functional difference between this
 250     and TMatrixGenerator is that this remembers a specific BaseSpec
 251     and its __call__ method takes only one mandatory argument (in addition to self).
 252     '''
 253     def __init__(self, TMatrixGenerator tmg, BaseSpec spec):
 254         self.generator = tmg
 255         self.spec = spec
 256         self.f.gen = self.generator.rawpointer()
 257         self.f.spec = self.spec.rawpointer()
 258
 259     def __call__(self, cdouble omega, fill = None):
 260         cdef CTMatrix tm
 261         if fill is None: # make a new CTMatrix
 262             tm = CTMatrix(self.spec, None)
 263         else: # TODO check whether fill has the same bspec as self?
 264             tm = fill
 265         if self.f.gen.function(tm.rawpointer(), omega, self.f.gen.params) != 0:
 266             raise ValueError("Something went wrong")
 267         else:
 268             return tm
 269
 270
 271 cdef class TMatrixGenerator:
 272     def __init__(self, what):
 273         if isinstance(what, __MieParams):
 274             self.holder = what
 275             self.g.function = qpms_tmatrix_generator_sphere
 276             self.g.params = (<__MieParams?>self.holder).rawpointer()
 277         elif isinstance(what,__AxialSymParams):
 278             self.holder = what
 279             self.g.function = qpms_tmatrix_generator_axialsym
 280             self.g.params = (<__AxialSymParams?>self.holder).rawpointer()
 281         elif isinstance(what, CTMatrix):
 282             self.holder = what
 283             self.g.function = qpms_tmatrix_generator_constant
 284             self.g.params = <void*>(<CTMatrix?>self.holder).rawpointer()
 285         elif isinstance(what, TMatrixInterpolator):
 286             self.holder = what
 287             self.g.function = qpms_tmatrix_generator_interpolator
 288             self.g.params = <void*>(<TMatrixInterpolator?>self.holder).rawpointer()
 289         elif callable(what):
 290             warnings.warn("Custom python T-matrix generators are an experimental feature. Also expect it to be slow.")
 291             self.holder = what
 292             self.g.function = qpms_tmatrix_generator_pythoncallable
 293             self.g.params = <void*>self.holder
 294         else:
 295             raise TypeError("Can't construct TMatrixGenerator from that")
 296
 297     def __call__(self, arg, cdouble omega):
 298         """Produces a T-matrix at a given frequency.
 299
 300         Parameters
 301         ----------
 302         arg : CTMatrix or BaseSpec
 303             If arg is a CTMatrix, its contents will be replaced.
 304             If arg is a BaseSpec, a new CTMatrix instance will be created.
 305         omega : complex
 306             Angular frequency at which the T-matrix shall be evaluated.
 307
 308         Returns
 309         -------
 310         t : CTMatrix
 311         """
 312         cdef CTMatrix tm
 313         if isinstance(arg, CTMatrix): # fill the matrix
 314             tm = arg
 315             if self.g.function(tm.rawpointer(), omega, self.g.params) != 0:
 316                 raise ValueError("Something went wrong")
 317             return
 318         elif isinstance(arg, BaseSpec): # Make a new CTMatrix
 319             tm = CTMatrix(arg, None)
 320             if self.g.function(tm.rawpointer(), omega, self.g.params) != 0:
 321                 raise ValueError("Something went wrong")
 322             return tm
 323         else:
 324             raise ValueError("Must specify CTMatrix or BaseSpec")
 325
 326     def Q_transposed(self, cdouble omega, norm):
 327         if self.g.function != qpms_tmatrix_generator_axialsym:
 328             raise TypeError("Applicable only for axialsym generators")
 329         return self.holder.Q_transposed(omega, norm)
 330     def R_transposed(self, cdouble omega, norm):
 331         if self.g.function != qpms_tmatrix_generator_axialsym:
 332             raise TypeError("Applicable only for axialsym generators")
 333         return self.holder.R_transposed(omega, norm)
 334
 335     # Better "constructors":
 336     @staticmethod
 337     def sphere(outside, inside, r):
 338         """Creates a T-matrix generator for a spherical scatterer.
 339
 340         This method uses analytical Mie-Lorentz formulae.
 341
 342         Parameters:
 343         -----------
 344         outside : EpsMuGenerator or EpsMu
 345             Optical properties of the surrounding medium.
 346         inside : EpsMuGenerator or EpsMu
 347             Optical properties of the material inside the sphere.
 348         r : double
 349             Sphere radius
 350         """
 351         return TMatrixGenerator(__MieParams(EpsMuGenerator(outside),
 352                     EpsMuGenerator(inside), r))
 353
 354     @staticmethod
 355     def sphere_asarc(outside, inside, r, *args, **kwargs):
 356         """Creates a T-matrix generator for a spherical scatterer.
 357
 358         This method uses numerical evaluation for generic axially-symmetric
 359         scatterer, and is intended for testing and benchmarking only.
 360         For regular use, see TMatrigGenerator.sphere() instead.
 361
 362         Parameters
 363         ----------
 364         outside : EpsMuGenerator or EpsMu
 365             Optical properties of the surrounding medium.
 366         inside : EpsMuGenerator or EpsMu
 367             Optical properties of the material inside the sphere.
 368         r : double
 369             Sphere radius
 370
 371         Returns
 372         -------
 373         tmgen : TMatrixGenerator
 374
 375         See Also
 376         --------
 377         TMatrigGenerator.sphere : Faster and more precise method.
 378         """
 379         return TMatrixGenerator(__AxialSymParams(
 380             EpsMuGenerator(outside), EpsMuGenerator(inside),
 381             ArcFunction(__ArcSphere(r)), *args, **kwargs))
 382
 383     @staticmethod
 384     def cylinder(outside, inside, r, h, *args, **kwargs):
 385         """Creates a T-matrix generator for a right circular cylinder.
 386
 387         Parameters:
 388         -----------
 389         outside : EpsMuGenerator or EpsMu
 390             Optical properties of the surrounding medium.
 391         inside : EpsMuGenerator or EpsMu
 392             Optical properties of the material inside the cylinder.
 393         r : double
 394             Cylinder base radius.
 395         h : double
 396             Cylinder height.
 397
 398         Returns
 399         -------
 400         tmgen : TMatrixGenerator
 401         """
 402         return TMatrixGenerator(__AxialSymParams(
 403             EpsMuGenerator(outside), EpsMuGenerator(inside),
 404             ArcFunction(__ArcCylinder(r, h)), *args, **kwargs))
 405