qpms/cytranslations.pyx

   1 import numpy as np
   2 import cmath
   3 from .qpms_cdefs cimport *
   4 from .cycommon import *
   5 from .cybspec cimport *
   6 cimport cython
   7 from cython.parallel cimport parallel, prange
   8 from libc.stdlib cimport malloc, free, calloc, abort
   9
  10 cdef int q_max(int m, int n, int mu, int nu):
  11     return min(n,nu,(n+nu-abs(m+mu)//2))
  12
  13 """
  14 Now we generate our own universal functions to be used with numpy.
  15
  16 Good way to see how this is done is to look at scipy/scipy/special/generate_ufuncs.py
  17 and scipy/scipy/special/generate_ufuncs.py
  18
  19 In simple words, it works like this:
  20 - Let's have a single element function. This can be function which returns or a "subroutine".
  21 - Then we need a loop function; this is a wrapper that gets bunch of pointers from numpy and
  22   has to properly call the single element function.
  23 - From those two, we build a python object using PyUFunc_FromFuncAndData.
  24   * If the ufunc is supposed to work on different kinds of input/output types,
  25     then a pair of single-element and loop functions is o be provided for
  26     each combination of types. However, the single-element function can be reused if
  27     the corresponding loop functions do the proper casting.
  28 """
  29
  30 ## as in scipy/special/_ufuncs_cxx.pyx
  31 ##-------------------------------------
  32 #cdef extern from "numpy/ufuncobject.h":
  33 #    int PyUFunc_getfperr() nogil
  34 #
  35 #cdef public int wrap_PyUFunc_getfperr() nogil:
  36 #    """
  37 #    Call PyUFunc_getfperr in a context where PyUFunc_API array is initialized;
  38 #
  39 #    """
  40 #    return PyUFunc_getfperr()
  41 #
  42 #cimport sf_error
  43 #-------------------------------------
  44
  45
  46 # This one is probably not used anymore an can perhaps be removed:
  47 cdef void loop_D_iiiidddii_As_D_lllldddbl(char **args, np.npy_intp *dims, np.npy_intp *steps, void *data) nogil:
  48     cdef np.npy_intp i, n = dims[0]
  49     cdef void *func = (<void**>data)#[0]
  50     #cdef char *func_name= <char*>(<void**>data)[1] # i am not using this, nor have I saved func_name to data
  51     cdef char *ip0 = args[0]
  52     cdef char *ip1 = args[1]
  53     cdef char *ip2 = args[2]
  54     cdef char *ip3 = args[3]
  55     cdef char *ip4 = args[4]
  56     cdef char *ip5 = args[5]
  57     cdef char *ip6 = args[6]
  58     cdef char *ip7 = args[7]
  59     cdef char *ip8 = args[8]
  60     cdef char *op0 = args[9]
  61     cdef cdouble ov0
  62     for i in range(n): # iterating over dimensions
  63         ov0 = (<double complex(*)(int, int, int, int, double, double, double, int, int) nogil>func)(
  64             <int>(<np.npy_long*>ip0)[0],
  65             <int>(<np.npy_long*>ip1)[0],
  66             <int>(<np.npy_long*>ip2)[0],
  67             <int>(<np.npy_long*>ip3)[0],
  68             <double>(<np.npy_double*>ip4)[0],
  69             <double>(<np.npy_double*>ip5)[0],
  70             <double>(<np.npy_double*>ip6)[0],
  71             <int>(<np.npy_bool*>ip7)[0],
  72             <int>(<np.npy_long*>ip8)[0],
  73         )
  74         (<cdouble *>op0)[0] = <cdouble>ov0
  75         ip0 += steps[0]
  76         ip1 += steps[1]
  77         ip2 += steps[2]
  78         ip3 += steps[3]
  79         ip4 += steps[4]
  80         ip5 += steps[5]
  81         ip6 += steps[6]
  82         ip7 += steps[7]
  83         ip8 += steps[8]
  84         op0 += steps[9]
  85 # FIXME ERROR HANDLING!!! requires correct import and different data passed (see scipy's generated ufuncs)
  86 #    sf_error.check_fpe(func_name)
  87
  88
  89
  90
  91 # Module initialisation
  92 # ---------------------
  93
  94 np.import_array()  # not sure whether this is really needed
  95 np.import_ufunc()
  96
  97 # Arrays passed to PyUFunc_FromFuncAndData()
  98 # ------------------------------------------
  99
 100 # BTW, aren't there anonymous arrays in cython?
 101
 102
 103 # types to be used for all of the single-type translation
 104 # coefficient retrieval ufuncs called like
 105 # coeff = func(m, n, mu, nu, r, theta, phi, r_ge_d, J)
 106 # currently supported signatures: (D_lllldddbl)
 107 cdef char ufunc__get_either_trans_coeff_types[10]
 108 ufunc__get_either_trans_coeff_types[0] = np.NPY_LONG
 109 ufunc__get_either_trans_coeff_types[1] = np.NPY_LONG
 110 ufunc__get_either_trans_coeff_types[2] = np.NPY_LONG
 111 ufunc__get_either_trans_coeff_types[3] = np.NPY_LONG
 112 ufunc__get_either_trans_coeff_types[4] = np.NPY_DOUBLE
 113 ufunc__get_either_trans_coeff_types[5] = np.NPY_DOUBLE
 114 ufunc__get_either_trans_coeff_types[6] = np.NPY_DOUBLE
 115 ufunc__get_either_trans_coeff_types[7] = np.NPY_BOOL
 116 ufunc__get_either_trans_coeff_types[8] = np.NPY_LONG
 117 ufunc__get_either_trans_coeff_types[9] = np.NPY_CDOUBLE
 118
 119 # types to be used for all of the both-type translation
 120 # coefficient retrieval ufuncs called like
 121 # errval = func(m, n, mu, nu, r, theta, phi, r_ge_d, J, &A, &B)
 122 # currently supported signatures: (lllldddbl_DD)
 123 cdef char ufunc__get_both_coeff_types[11]
 124 ufunc__get_both_coeff_types[0] = np.NPY_LONG
 125 ufunc__get_both_coeff_types[1] = np.NPY_LONG
 126 ufunc__get_both_coeff_types[2] = np.NPY_LONG
 127 ufunc__get_both_coeff_types[3] = np.NPY_LONG
 128 ufunc__get_both_coeff_types[4] = np.NPY_DOUBLE
 129 ufunc__get_both_coeff_types[5] = np.NPY_DOUBLE
 130 ufunc__get_both_coeff_types[6] = np.NPY_DOUBLE
 131 ufunc__get_both_coeff_types[7] = np.NPY_BOOL
 132 ufunc__get_both_coeff_types[8] = np.NPY_LONG
 133 ufunc__get_both_coeff_types[9] = np.NPY_CDOUBLE
 134 ufunc__get_both_coeff_types[10] = np.NPY_CDOUBLE
 135
 136
 137
 138 # ---------------------------------------------
 139 # Wrapper for the qpms_trans_calculator "class"
 140 # ---------------------------------------------
 141 ctypedef struct trans_calculator_get_X_data_t:
 142     qpms_trans_calculator* c
 143     void* cmethod
 144
 145 cdef void trans_calculator_loop_D_Ciiiidddii_As_D_lllldddbl(char **args, np.npy_intp *dims, np.npy_intp *steps, void *data) nogil:
 146     cdef np.npy_intp i, n = dims[0]
 147     cdef void *func = (<trans_calculator_get_X_data_t*>data)[0].cmethod
 148     #cdef cdouble (*func)(qpms_trans_calculator*, int, int, int, int, double, double, double, int, int) nogil = (<trans_calculator_get_X_data_t*>data)[0].cmethod
 149     cdef qpms_trans_calculator* c = (<trans_calculator_get_X_data_t*>data)[0].c
 150     #cdef char *func_name= <char*>(<void**>data)[1] # i am not using this, nor have I saved func_name to data
 151     cdef char *ip0 = args[0]
 152     cdef char *ip1 = args[1]
 153     cdef char *ip2 = args[2]
 154     cdef char *ip3 = args[3]
 155     cdef char *ip4 = args[4]
 156     cdef char *ip5 = args[5]
 157     cdef char *ip6 = args[6]
 158     cdef char *ip7 = args[7]
 159     cdef char *ip8 = args[8]
 160     cdef char *op0 = args[9]
 161     cdef cdouble ov0
 162     for i in range(n): # iterating over dimensions
 163         #ov0 = func(
 164         ov0 = (<double complex(*)(qpms_trans_calculator*, int, int, int, int, double, double, double, int, int) nogil>func)(
 165             c,
 166             <int>(<np.npy_long*>ip0)[0],
 167             <int>(<np.npy_long*>ip1)[0],
 168             <int>(<np.npy_long*>ip2)[0],
 169             <int>(<np.npy_long*>ip3)[0],
 170             <double>(<np.npy_double*>ip4)[0],
 171             <double>(<np.npy_double*>ip5)[0],
 172             <double>(<np.npy_double*>ip6)[0],
 173             <int>(<np.npy_bool*>ip7)[0],
 174             <int>(<np.npy_long*>ip8)[0],
 175         )
 176         (<cdouble *>op0)[0] = <cdouble>ov0
 177         ip0 += steps[0]
 178         ip1 += steps[1]
 179         ip2 += steps[2]
 180         ip3 += steps[3]
 181         ip4 += steps[4]
 182         ip5 += steps[5]
 183         ip6 += steps[6]
 184         ip7 += steps[7]
 185         ip8 += steps[8]
 186         op0 += steps[9]
 187 # FIXME ERROR HANDLING!!! requires correct import and different data passed (see scipy's generated ufuncs)
 188 #    sf_error.check_fpe(func_name)
 189
 190
 191 cdef void trans_calculator_loop_E_C_DD_iiiidddii_As_lllldddbl_DD(char **args, np.npy_intp *dims, np.npy_intp *steps, void *data) nogil:
 192     # E stands for error value (int), C for qpms_trans_calculator*
 193     cdef np.npy_intp i, n = dims[0]
 194     cdef void *func = (<trans_calculator_get_X_data_t*>data)[0].cmethod
 195     #cdef complex double (*func)(qpms_trans_calculator*, double complex *, double complex *, int, int, int, int, double, double, double, int, int) nogil = (<trans_calculator_get_X_data_t*>data)[0].cmethod
 196     cdef qpms_trans_calculator* c = (<trans_calculator_get_X_data_t*>data)[0].c
 197     #cdef char *func_name= <char*>(<void**>data)[1] # i am not using this, nor have I saved func_name to data
 198     cdef char *ip0 = args[0]
 199     cdef char *ip1 = args[1]
 200     cdef char *ip2 = args[2]
 201     cdef char *ip3 = args[3]
 202     cdef char *ip4 = args[4]
 203     cdef char *ip5 = args[5]
 204     cdef char *ip6 = args[6]
 205     cdef char *ip7 = args[7]
 206     cdef char *ip8 = args[8]
 207     cdef char *op0 = args[9]
 208     cdef char *op1 = args[10]
 209     cdef cdouble ov0
 210     cdef int errval
 211     for i in range(n): # iterating over dimensions
 212         #errval = func(
 213         errval = (<int(*)(qpms_trans_calculator*, double complex *, double complex *, int, int, int, int, double, double, double, int, int) nogil>func)(
 214             c,
 215             <cdouble *> op0,
 216             <cdouble *> op1,
 217             <int>(<np.npy_long*>ip0)[0],
 218             <int>(<np.npy_long*>ip1)[0],
 219             <int>(<np.npy_long*>ip2)[0],
 220             <int>(<np.npy_long*>ip3)[0],
 221             <double>(<np.npy_double*>ip4)[0],
 222             <double>(<np.npy_double*>ip5)[0],
 223             <double>(<np.npy_double*>ip6)[0],
 224             <int>(<np.npy_bool*>ip7)[0],
 225             <int>(<np.npy_long*>ip8)[0],
 226         )
 227         ip0 += steps[0]
 228         ip1 += steps[1]
 229         ip2 += steps[2]
 230         ip3 += steps[3]
 231         ip4 += steps[4]
 232         ip5 += steps[5]
 233         ip6 += steps[6]
 234         ip7 += steps[7]
 235         ip8 += steps[8]
 236         op0 += steps[9]
 237         op1 += steps[10]
 238         # TODO if (errval != 0): ...
 239 # FIXME ERROR HANDLING!!! requires correct import and different data passed (see scipy's generated ufuncs)
 240 #    sf_error.check_fpe(func_name)
 241
 242 @cython.boundscheck(False)
 243 @cython.wraparound(False)
 244 cdef void trans_calculator_parallel_loop_E_C_DD_iiiidddii_As_lllldddbl_DD(char **args, np.npy_intp *dims, np.npy_intp *steps, void *data) nogil:
 245     # E stands for error value (int), C for qpms_trans_calculator*
 246     cdef np.npy_intp i, n = dims[0]
 247     cdef void *func = (<trans_calculator_get_X_data_t*>data)[0].cmethod
 248     #cdef complex double (*func)(qpms_trans_calculator*, double complex *, double complex *, int, int, int, int, double, double, double, int, int) nogil = (<trans_calculator_get_X_data_t*>data)[0].cmethod
 249     cdef qpms_trans_calculator* c = (<trans_calculator_get_X_data_t*>data)[0].c
 250     #cdef char *func_name= <char*>(<void**>data)[1] # i am not using this, nor have I saved func_name to data
 251
 252     cdef char *ip0
 253     cdef char *ip1
 254     cdef char *ip2
 255     cdef char *ip3
 256     cdef char *ip4
 257     cdef char *ip5
 258     cdef char *ip6
 259     cdef char *ip7
 260     cdef char *ip8
 261     cdef char *op0
 262     cdef char *op1
 263     cdef int errval
 264     for i in prange(n): # iterating over dimensions
 265         ip0 = args[0] + i * steps[0]
 266         ip1 = args[1] + i * steps[1]
 267         ip2 = args[2] + i * steps[2]
 268         ip3 = args[3] + i * steps[3]
 269         ip4 = args[4] + i * steps[4]
 270         ip5 = args[5] + i * steps[5]
 271         ip6 = args[6] + i * steps[6]
 272         ip7 = args[7] + i * steps[7]
 273         ip8 = args[8] + i * steps[8]
 274         op0 = args[9] + i * steps[9]
 275         op1 = args[10] + i * steps[10]
 276         #errval = func(
 277         errval = (<int(*)(qpms_trans_calculator*, double complex *, double complex *, int, int, int, int, double, double, double, int, int) nogil>func)(
 278             c,
 279             <cdouble *> op0,
 280             <cdouble *> op1,
 281             <int>(<np.npy_long*>ip0)[0],
 282             <int>(<np.npy_long*>ip1)[0],
 283             <int>(<np.npy_long*>ip2)[0],
 284             <int>(<np.npy_long*>ip3)[0],
 285             <double>(<np.npy_double*>ip4)[0],
 286             <double>(<np.npy_double*>ip5)[0],
 287             <double>(<np.npy_double*>ip6)[0],
 288             <int>(<np.npy_bool*>ip7)[0],
 289             <int>(<np.npy_long*>ip8)[0],
 290         )
 291        # TODO if (errval != 0): ...
 292 # FIXME ERROR HANDLING!!! requires correct import and different data passed (see scipy's generated ufuncs)
 293 #    sf_error.check_fpe(func_name)
 294
 295
 296 cdef np.PyUFuncGenericFunction trans_calculator_get_X_loop_funcs[1]
 297 trans_calculator_get_X_loop_funcs[0] = trans_calculator_loop_D_Ciiiidddii_As_D_lllldddbl
 298
 299 cdef np.PyUFuncGenericFunction trans_calculator_get_AB_loop_funcs[1]
 300 #trans_calculator_get_AB_loop_funcs[0] = trans_calculator_parallel_loop_E_C_DD_iiiidddii_As_lllldddbl_DD
 301 trans_calculator_get_AB_loop_funcs[0] = trans_calculator_loop_E_C_DD_iiiidddii_As_lllldddbl_DD
 302 cdef void *trans_calculator_get_AB_elementwise_funcs[1]
 303 trans_calculator_get_AB_elementwise_funcs[0] = <void *>qpms_trans_calculator_get_AB_p_ext
 304
 305 '''
 306 cdef extern from "numpy/ndarrayobject.h":
 307     struct PyArrayInterface:
 308         int itemsize
 309         np.npy_uintp *shape
 310         np.npy_uintp *strides
 311         void *data
 312 '''
 313
 314
 315
 316
 317
 318 cdef class trans_calculator:
 319     cdef qpms_trans_calculator* c
 320     cdef trans_calculator_get_X_data_t get_A_data[1]
 321     cdef trans_calculator_get_X_data_t* get_A_data_p[1]
 322
 323     cdef trans_calculator_get_X_data_t get_B_data[1]
 324     cdef trans_calculator_get_X_data_t* get_B_data_p[1]
 325
 326     cdef trans_calculator_get_X_data_t get_AB_data[1]
 327     cdef trans_calculator_get_X_data_t* get_AB_data_p[1]
 328     cdef public: # TODO CHECK FOR CORRECT REFERENCE COUNTING AND LEAKS
 329         # have to be cdef public in order that __init__ can set these attributes
 330         object get_A, get_B, get_AB
 331
 332     def __cinit__(self, int lMax, int normalization = QPMS_NORMALISATION_DEFAULT):
 333         if (lMax <= 0):
 334             raise ValueError('lMax has to be greater than 0.')
 335         self.c = qpms_trans_calculator_init(lMax, normalization)
 336         if self.c is NULL:
 337             raise MemoryError
 338
 339     def __init__(self, int lMax, int normalization = QPMS_NORMALISATION_DEFAULT):
 340         if self.c is NULL:
 341             raise MemoryError()
 342         self.get_A_data[0].c = self.c
 343         self.get_A_data[0].cmethod = <void *>qpms_trans_calculator_get_A_ext
 344         self.get_A_data_p[0] = &(self.get_A_data[0])
 345         self.get_A = <object>np.PyUFunc_FromFuncAndData(# TODO CHECK FOR CORRECT REFERENCE COUNTING AND LEAKS
 346                 trans_calculator_get_X_loop_funcs, # func
 347                 <void **>self.get_A_data_p, #data
 348                 ufunc__get_either_trans_coeff_types, #types
 349                 1, # ntypes: number of supported input types
 350                 9, # nin: number of input args
 351                 1, # nout: number of output args
 352                 0, # identity element, unused
 353                 "get_A", #name
 354                 """
 355                 TODO doc
 356                 """, # doc
 357                 0 # unused
 358                 )
 359         self.get_B_data[0].c = self.c
 360         self.get_B_data[0].cmethod = <void *>qpms_trans_calculator_get_B_ext
 361         self.get_B_data_p[0] = &(self.get_B_data[0])
 362         self.get_B = <object>np.PyUFunc_FromFuncAndData(# TODO CHECK FOR CORRECT REFERENCE COUNTING AND LEAKS
 363                 trans_calculator_get_X_loop_funcs, # func
 364                 <void **>self.get_B_data_p, #data
 365                 ufunc__get_either_trans_coeff_types, #types
 366                 1, # ntypes: number of supported input types
 367                 9, # nin: number of input args
 368                 1, # nout: number of output args
 369                 0, # identity element, unused
 370                 "get_B", #name
 371                 """
 372                 TODO doc
 373                 """, # doc
 374                 0 # unused
 375                 )
 376         self.get_AB_data[0].c = self.c
 377         self.get_AB_data[0].cmethod = <void *>qpms_trans_calculator_get_AB_p_ext
 378         self.get_AB_data_p[0] = &(self.get_AB_data[0])
 379         self.get_AB = <object>np.PyUFunc_FromFuncAndData(# TODO CHECK FOR CORRECT REFERENCE COUNTING AND LEAKS
 380                 trans_calculator_get_AB_loop_funcs, # func
 381                 <void **>self.get_AB_data_p, #data
 382                 ufunc__get_both_coeff_types, #types
 383                 1, # ntypes: number of supported input types
 384                 9, # nin: number of input args
 385                 2, # nout: number of output args
 386                 0, # identity element, unused
 387                 "get_AB", #name
 388                 """
 389                 TODO doc
 390                 """, # doc
 391                 0 # unused
 392                 )
 393     def __dealloc__(self):
 394         if self.c is not NULL:
 395             qpms_trans_calculator_free(self.c)
 396         # TODO Reference counts to get_A, get_B, get_AB?
 397
 398     def lMax(self):
 399         return self.c[0].lMax
 400
 401     def nelem(self):
 402         return self.c[0].nelem
 403
 404     # THIS FUNCTION MIGHT BE OBSOLETE; NOT SURE WHETHER IT'S WORKING ANY MORE
 405     def get_AB_arrays(self, r, theta, phi, r_ge_d, int J,
 406             destaxis=None, srcaxis=None, expand=True):
 407         """
 408         Returns arrays of translation coefficients, inserting two new nelem-sized axes
 409         (corresponding to the destination and source axes of the translation matrix,
 410         respectively).
 411
 412         By default (expand==True), it inserts the new axes. or it can be provided with
 413         the resulting shape (with the corresponding axes dimensions equal to one).
 414         The provided axes positions are for the resulting array.
 415
 416         If none axis positions are provided, destaxis and srcaxis will be the second-to-last
 417         and last, respectively.
 418         """
 419         # TODO CHECK (and try to cast) INPUT ARRAY TYPES (now is done)
 420         # BIG FIXME: make skalars valid arguments, now r, theta, phi, r_ge_d have to be ndarrays
 421         cdef:
 422             int daxis, saxis, smallaxis, bigaxis, resnd, i, j, d, ax, errval
 423             np.npy_intp sstride, dstride, longi
 424             int *local_indices
 425             char *r_p
 426             char *theta_p
 427             char *phi_p
 428             char *r_ge_d_p
 429             char *a_p
 430             char *b_p
 431         # Process the array shapes
 432         baseshape = np.broadcast(r,theta,phi,r_ge_d).shape # nope, does not work as needed
 433         '''
 434         cdef int r_orignd = r.ndim if hasattr(r, "ndim") else 0
 435         cdef int theta_orignd = theta.ndim if hasattr(theta, "ndim") else 0
 436         cdef int phi_orignd = phi.ndim if hasattr(phi, "ndim") else 0
 437         cdef int r_ge_d_orignd = r_ge_d.ndim if hasattr(r_ge_d, "__len__") else 0
 438         cdef int basend = max(r_orignd, theta_orignd, phi_orignd, r_ge_d_orignd)
 439         baseshape = list()
 440         for d in range(basend):
 441             baseshape.append(max(
 442                     r.shape[d+r_orignd-basend] if d+r_orignd-basend >= 0 else 1,
 443                     theta.shape[d+theta_orignd-basend] if d+theta_orignd-basend >= 0 else 1,
 444                     phi.shape[d+phi_orignd-basend] if d+phi_orignd-basend >= 0 else 1,
 445                     r_ge_d.shape[d+r_ge_d_orignd-basend] if d+r_ge_d_orignd-basend >= 0 else 1,
 446                     ))
 447         baseshape = tuple(baseshape)
 448         '''
 449         if not expand:
 450             resnd = len(baseshape)
 451             if resnd < 2:
 452                 raise ValueError('Translation matrix arrays must have at least 2 dimensions!')
 453             daxis = (resnd-2) if destaxis is None else destaxis
 454             saxis = (resnd-1) if srcaxis is None else srcaxis
 455             if daxis < 0:
 456                 daxis = resnd + daxis
 457             if saxis < 0:
 458                 saxis = resnd + saxis
 459             if daxis < 0 or saxis < 0 or daxis >= resnd or saxis >= resnd or daxis == saxis:
 460                 raise ValueError('invalid axes provided (destaxis = %d, srcaxis = %d, # of axes: %d'
 461                         % (daxis, saxis, resnd))
 462             if baseshape[daxis] != 1 or baseshape[saxis] != 1:
 463                 raise ValueError('dimension mismatch (input argument dimensions have to be 1 both at'
 464                         'destaxis (==%d) and srcaxis (==%d) but are %d and %d' %
 465                         (daxis, saxis, baseshape[daxis], baseshape[saxis]))
 466             resultshape = list(baseshape)
 467         else:
 468             resnd = len(baseshape)+2
 469             daxis = (resnd-2) if destaxis is None else destaxis
 470             saxis = (resnd-1) if srcaxis is None else srcaxis
 471             if daxis < 0:
 472                 daxis = resnd + daxis
 473             if saxis < 0:
 474                 saxis = resnd + saxis
 475             if daxis < 0 or saxis < 0 or daxis >= resnd or saxis >= resnd or daxis == saxis:
 476                 raise ValueError('invalid axes provided') # TODO better error formulation
 477             resultshape = list(baseshape)
 478             if daxis > saxis:
 479                 smallaxis = saxis
 480                 bigaxis = daxis
 481             else:
 482                 smallaxis = daxis
 483                 bigaxis = saxis
 484             resultshape.insert(smallaxis,1)
 485             resultshape.insert(bigaxis,1)
 486             r = np.expand_dims(np.expand_dims(r.astype(np.float_, copy=False), smallaxis), bigaxis)
 487             theta = np.expand_dims(np.expand_dims(theta.astype(np.float_, copy=False), smallaxis), bigaxis)
 488             phi = np.expand_dims(np.expand_dims(phi.astype(np.float_, copy=False), smallaxis), bigaxis)
 489             r_ge_d = np.expand_dims(np.expand_dims(r_ge_d.astype(np.bool_, copy=False), smallaxis), bigaxis)
 490
 491         resultshape[daxis] = self.c[0].nelem
 492         resultshape[saxis] = self.c[0].nelem
 493         cdef np.ndarray r_c = np.broadcast_to(r,resultshape)
 494         cdef np.ndarray theta_c = np.broadcast_to(theta,resultshape)
 495         cdef np.ndarray phi_c = np.broadcast_to(phi,resultshape)
 496         cdef np.ndarray r_ge_d_c = np.broadcast_to(r_ge_d, resultshape)
 497         cdef np.ndarray a = np.empty(resultshape, dtype=complex)
 498         cdef np.ndarray b = np.empty(resultshape, dtype=complex)
 499         dstride = a.strides[daxis]
 500         sstride = a.strides[saxis]
 501         with nogil:
 502             errval = qpms_cython_trans_calculator_get_AB_arrays_loop(
 503                     self.c, J, resnd,
 504                     daxis, saxis,
 505                     a.data, a.shape, a.strides,
 506                     b.data, b.shape, b.strides,
 507                     r_c.data, r_c.shape, r_c.strides,
 508                     theta_c.data, theta_c.shape, theta_c.strides,
 509                     phi_c.data, phi_c.shape, phi_c.strides,
 510                     r_ge_d_c.data, r_ge_d_c.shape, r_ge_d_c.strides
 511                     )
 512         return a, b
 513
 514     def get_trans_array_bspec_sph(self, BaseSpec destspec, BaseSpec srcspec,
 515             kdlj, qpms_bessel_t J = QPMS_HANKEL_PLUS):
 516         kdlj = np.array(kdlj)
 517         if kdlj.shape != (3,):
 518             raise ValueError("Array of shape (3,) with spherical coordinates of the translation expected")
 519         cdef size_t destn = len(destspec)
 520         cdef size_t srcn = len(srcspec)
 521         cdef np.ndarray[np.complex_t, ndim=2] target = np.empty(
 522                 (destn, srcn), dtype=complex, order='C')
 523         cdef cdouble[:,::1] target_view = target
 524         cdef csph_t kdlj_sph
 525         kdlj_sph.r = kdlj[0]
 526         kdlj_sph.theta = kdlj[1].real
 527         kdlj_sph.phi = kdlj[2].real
 528         qpms_trans_calculator_get_trans_array(self.c, &target_view[0][0],
 529                 destspec.rawpointer(), srcn, srcspec.rawpointer(), 1,
 530                 kdlj_sph, False, J)
 531         return target
 532
 533     def get_trans_array_bspec_c3pos(self, BaseSpec destspec, BaseSpec srcspec,
 534             cdouble k, destpos, srcpos, qpms_bessel_t J = QPMS_HANKEL_PLUS):
 535         destpos = np.array(destpos)
 536         srcpos = np.array(srcpos)
 537         if destpos.shape != (3,) or srcpos.shape != (3,):
 538             raise ValueError("Array of shape (3,) with cartesian coordinates of the particle position expected")
 539         cdef size_t destn = len(destspec)
 540         cdef size_t srcn = len(srcspec)
 541         cdef np.ndarray[np.complex_t, ndim=2] target = np.empty(
 542                 (destn, srcn), dtype=complex, order='C')
 543         cdef cdouble[:,::1] target_view = target
 544         cdef cart3_t srcp, destp
 545         srcp.x = srcpos[0]
 546         srcp.y = srcpos[1]
 547         srcp.z = srcpos[2]
 548         destp.x = destpos[0]
 549         destp.y = destpos[1]
 550         destp.z = destpos[2]
 551         qpms_trans_calculator_get_trans_array_lc3p(self.c, &target_view[0][0],
 552                 destspec.rawpointer(), srcn, srcspec.rawpointer(), 1, k,
 553                 destp, srcp, J)
 554         return target
 555
 556
 557     # TODO make possible to access the attributes (to show normalization etc)
 558