diff --git a/qpms/cytmatrices.pxd b/qpms/cytmatrices.pxd index c955257..51641d9 100644 --- a/qpms/cytmatrices.pxd +++ b/qpms/cytmatrices.pxd @@ -11,6 +11,8 @@ cdef class TMatrixInterpolator: cdef double *freqs_su cdef size_t nfreqs cdef qpms_tmatrix_interpolator_t *interp + cdef inline qpms_tmatrix_interpolator_t *rawpointer(self): + return self.interp cdef class CTMatrix: # N.B. there is another type called TMatrix in tmatrices.py! cdef readonly np.ndarray m # Numpy array holding the matrix data @@ -23,3 +25,7 @@ cdef class CTMatrix: # N.B. there is another type called TMatrix in tmatrices.py Don't forget to reference the BaseSpec object itself when storing the pointer anywhere!!! ''' return &(self.t) + + + + diff --git a/qpms/cytmatrices.pyx b/qpms/cytmatrices.pyx index fb19961..a15dc33 100644 --- a/qpms/cytmatrices.pyx +++ b/qpms/cytmatrices.pyx @@ -3,6 +3,7 @@ from .qpms_cdefs cimport * from .cybspec cimport BaseSpec from .cycommon import * from .cycommon cimport make_c_string +from .cymaterials cimport EpsMuGenerator from .qpms_c cimport FinitePointGroup import warnings import os @@ -114,3 +115,107 @@ cdef class CTMatrix: # N.B. there is another type called TMatrix in tmatrices.py tm = CTMatrix(spec, 0) tm.spherical_perm_fill(radius, freq, epsilon_int, epsilon_ext) return tm + +cdef class __MieParams: + # Not directly callable right now, serves just to be used by TMatrixGenerator. + cdef qpms_tmatrix_generator_sphere_param_t cparam + cdef EpsMuGenerator outside + cdef EpsMuGenerator inside + cdef inline void *rawpointer(self): + return &(self.cparam) + + def __init__(self, outside, inside, r): + self.inside = inside + self.outside = outside + self.cparam.inside = self.inside.raw(); + self.cparam.outside = self.outside.raw(); + self.cparam.radius = r + + property r: + def __get__(self): + return self.cparam.radius + def __set__(self, val): + self.cparam.radius = val + +cdef class __ArcCylinder: + cdef qpms_arc_cylinder_params_t p + cdef inline void *rawpointer(self): + return &(self.p) + +cdef class __ArcSphere: + cdef double r + cdef inline void *rawpointer(self): + return &(self.r) + +cdef qpms_arc_function_retval_t userarc(double theta, const void *params): + cdef object fun = params + cdef qpms_arc_function_retval_t retval + retval.r, retval.beta = fun(theta) + return retval + + +cdef class ArcFunction: + cdef qpms_arc_function_t g + cdef object holder + def __init__(self, what): + if isinstance(what, __ArcCylinder): + self.holder = what + self.g.function = qpms_arc_cylinder + self.g.params = (<__ArcCylinder?>self.holder).rawpointer() + elif isinstance(what, __ArcSphere): + self.holder = what + self.g.function = qpms_arc_sphere + self.g.params = (<__ArcSphere?>self.holder).rawpointer() + elif callable(what): + warnings.warn("Custom python (r, beta) arc functions are an experimental feature. Also expect it to be slow.") + self.holder = what + self.g.function = userarc + self.g.params = self.holder + elif isinstance(what, ArcFunction): #Copy constructor + self.holder = what.holder + self.g = (what).g + self.holder.rawpointer() + +cdef class __AxialSymParams: + cdef qpms_tmatrix_generator_axialsym_param_t p + cdef EpsMuGenerator outside + cdef EpsMuGenerator inside + cdef ArcFunction shape + cdef void * rawpointer(self): + return &(self.p) + property lMax_extend: + def __get__(self): + return self.p.lMax_extend + def __set__(self, val): + self.p.lMax_extend = val + def __init__(self, outside, inside, shape, *args, **kwargs): + self.outside = outside + self.p.outside = self.outside.g + self.inside = inside + self.p.inside = self.inside.g + self.shape = shape + self.p.shape = self.shape.g + if len(args)>0: + self.lMax_extend = args[0] + if 'lMax_extend' in kwargs.keys(): + self.lMax_extend = kwargs['lMax_extend'] + +cdef class TMatrixGenerator: + cdef qpms_tmatrix_generator_t g + cdef object holder + cdef qpms_tmatrix_generator_t raw(self): + return self.g + def __init__(self, what): + if isinstance(what, __MieParams): + self.holder = what + self.g.function = qpms_tmatrix_generator_sphere + self.g.params = (<__MieParams?>self.holder).rawpointer() + elif isinstance(what,__AxialSymParams): + self.holder = what + self.g.function = qpms_tmatrix_generator_axialsym + self.g.params = (<__AxialSymParams?>self.holder).rawpointer() + else: + raise ValueError("Can't construct TMatrixGenerator from that") + + + diff --git a/qpms/tmatrices.c b/qpms/tmatrices.c index 6b27ad5..d1ce6ee 100644 --- a/qpms/tmatrices.c +++ b/qpms/tmatrices.c @@ -35,6 +35,7 @@ #define SQ(x) ((x)*(x)) #define MAX(x,y) ((x) < (y) ? (y) : (x)) +#define MIN(x,y) ((x) > (y) ? (y) : (x)) qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec) { qpms_tmatrix_t *t = malloc(sizeof(qpms_tmatrix_t)); @@ -692,7 +693,7 @@ qpms_errno_t qpms_tmatrix_axialsym_fill( qpms_tmatrix_t *t, complex double omega, qpms_epsmu_t outside, qpms_epsmu_t inside,qpms_arc_function_t shape, qpms_l_t lMaxQR) { - QPMS_UNTESTED; + QPMS_UNTESTED; // TODO also check whether this is valid for all norms. const qpms_vswf_set_spec_t *bspec = t->spec; struct tmatrix_axialsym_integral_param_t p; p.k = qpms_wavenumber(omega, outside); @@ -797,3 +798,16 @@ qpms_errno_t qpms_tmatrix_generator_axialsym(qpms_tmatrix_t *t, complex double o p->shape, p->lMax_extend); } + +qpms_errno_t qpms_tmatrix_generator_constant(qpms_tmatrix_t *t, + complex double omega, const void *tmorig) { + QPMS_UNTESTED; + // TODO more thorough check on bspec! + const qpms_tmatrix_t *orig = tmorig; + const size_t tocopy = SQ(MIN(t->spec->n, orig->spec->n)); + memcpy(t->m, orig->m, tocopy*sizeof(complex double)); + memset(t->m+tocopy, 0, (SQ(t->spec->n)-tocopy)*sizeof(complex double)); + return QPMS_SUCCESS; +} + + diff --git a/qpms/tmatrices.h b/qpms/tmatrices.h index ab21a56..a642011 100644 --- a/qpms/tmatrices.h +++ b/qpms/tmatrices.h @@ -241,6 +241,7 @@ complex double *qpms_apply_tmatrix( * qpms_tmatrix_generator_axialsym() * qpms_tmatrix_generator_interpolator() * qpms_tmatrix_generator_sphere() + * qpms_tmatrix_generator_constant() */ typedef struct qpms_tmatrix_generator_t { qpms_errno_t (*function) (qpms_tmatrix_t *t, ///< T-matrix to fill. @@ -250,6 +251,16 @@ typedef struct qpms_tmatrix_generator_t { const void *params; ///< Parameter pointer passed to the function. } qpms_tmatrix_generator_t; +/// Implementation of qpms_matrix_generator_t that just copies a constant matrix. +/** N.B. this does almost no checks at all, so use it only for t-matrices with + * the same base spec. + */ +qpms_errno_t qpms_tmatrix_generator_constant(qpms_tmatrix_t *t, + complex double omega, + /// Source T-matrix, real type is (const qpms_tmatrix_t*). + const void *tmatrix_orig + ); + /* Fuck this, include the whole typedef struct gsl_spline gsl_spline; // Forward declaration for the interpolator struct typedef struct gsl_interp_type gsl_interp_type;