+001: r""" Cython-optimized mathematical utilities for the Poisson-Schrodinger problem.

  __pyx_t_4 = __Pyx_PyDict_NewPresized(0); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  if (PyDict_SetItem(__pyx_mstate_global->__pyx_d, __pyx_mstate_global->__pyx_n_u_test, __pyx_t_4) < (0)) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;

 002: 

 003: This module contains a couple useful high-performance math functions, including

 004: a tridiagonal matrix solver

 005: and an approximate Fermi-Dirac integral of order 1/2 and -1/2.

 006: These functions are compiled by Cython, and are tested for accuracy and speed by the

 007: :py:mod:`~pynitride.tests.core.test_cython_maths` module.

 008: """

 009: 

 010: cimport cython

 011: cimport numpy as cnp

+012: import numpy as np

  __pyx_t_1 = __Pyx_Import(__pyx_mstate_global->__pyx_n_u_numpy, 0, 0, NULL, 0); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 12, __pyx_L1_error)
  __pyx_t_4 = __pyx_t_1;
  __Pyx_GOTREF(__pyx_t_4);
  if (PyDict_SetItem(__pyx_mstate_global->__pyx_d, __pyx_mstate_global->__pyx_n_u_np, __pyx_t_4) < (0)) __PYX_ERR(0, 12, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;

 013: from libc.math cimport pow, exp, sqrt

+014: from scipy.sparse import lil_matrix

  {
    PyObject* const __pyx_imported_names[] = {__pyx_mstate_global->__pyx_n_u_lil_matrix};
    __pyx_t_1 = __Pyx_Import(__pyx_mstate_global->__pyx_n_u_scipy_sparse, __pyx_imported_names, 1, NULL, 0); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 14, __pyx_L1_error)
  }
  __pyx_t_4 = __pyx_t_1;
  __Pyx_GOTREF(__pyx_t_4);
  {
    PyObject* const __pyx_imported_names[] = {__pyx_mstate_global->__pyx_n_u_lil_matrix};
    __pyx_t_9 = 0; {
      __pyx_t_5 = __Pyx_ImportFrom(__pyx_t_4, __pyx_imported_names[__pyx_t_9]); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 14, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_5);
      if (PyDict_SetItem(__pyx_mstate_global->__pyx_d, __pyx_imported_names[__pyx_t_9], __pyx_t_5) < (0)) __PYX_ERR(0, 14, __pyx_L1_error)
      __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    }
  }
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;

+015: from functools import wraps

  {
    PyObject* const __pyx_imported_names[] = {__pyx_mstate_global->__pyx_n_u_wraps};
    __pyx_t_1 = __Pyx_Import(__pyx_mstate_global->__pyx_n_u_functools, __pyx_imported_names, 1, NULL, 0); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 15, __pyx_L1_error)
  }
  __pyx_t_4 = __pyx_t_1;
  __Pyx_GOTREF(__pyx_t_4);
  {
    PyObject* const __pyx_imported_names[] = {__pyx_mstate_global->__pyx_n_u_wraps};
    __pyx_t_9 = 0; {
      __pyx_t_5 = __Pyx_ImportFrom(__pyx_t_4, __pyx_imported_names[__pyx_t_9]); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 15, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_5);
      if (PyDict_SetItem(__pyx_mstate_global->__pyx_d, __pyx_imported_names[__pyx_t_9], __pyx_t_5) < (0)) __PYX_ERR(0, 15, __pyx_L1_error)
      __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    }
  }
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;

+016: cnp.import_array()

  __pyx_t_10 = __pyx_f_5numpy_import_array(); if (unlikely(__pyx_t_10 == ((int)-1))) __PYX_ERR(0, 16, __pyx_L1_error)

 017: from cpython cimport bool

 018: 

 019: ##########

 020: ### Tridiagonal matrix algorithm

 021: ##########

 022: 

 023: # c-function to implement tdma. See :py:func:`~pynitride.poissolve.maths.tmda` for more.

+024: @cython.boundscheck(False)

static PyArrayObject *__pyx_f_9pynitride_4core_12cython_maths_tdma_c(__Pyx_memviewslice __pyx_v_a, __Pyx_memviewslice __pyx_v_b, __Pyx_memviewslice __pyx_v_c, __Pyx_memviewslice __pyx_v_d) {
  int __pyx_v_N;
  int __pyx_v_k;
  double __pyx_v_m;
  PyArrayObject *__pyx_v_x_arr = 0;
  __Pyx_memviewslice __pyx_v_x = { 0, 0, { 0 }, { 0 }, { 0 } };
  PyArrayObject *__pyx_r = NULL;
  __Pyx_TraceDeclarationsFunc
  __Pyx_TraceFrameInit(((PyObject *)__pyx_mstate_global->__pyx_codeobj_tab[39]))
  __Pyx_TraceStartFunc("tdma_c", __pyx_f[0], 24, 0, 0, 0, __PYX_ERR(0, 24, __pyx_L1_error));
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_1);
  __Pyx_XDECREF(__pyx_t_2);
  __Pyx_XDECREF(__pyx_t_3);
  __Pyx_XDECREF(__pyx_t_4);
  __PYX_XCLEAR_MEMVIEW(&__pyx_t_6, 1);
  __Pyx_TraceException(__pyx_lineno, 0, 0);
  #if CYTHON_USE_SYS_MONITORING
  __Pyx_TraceExceptionUnwind(0, 0);
  #else
  __Pyx_TraceReturnValue(NULL, 0, 0, __PYX_ERR(0, 24, __pyx_L1_error));
  #endif
  __Pyx_AddTraceback("pynitride.core.cython_maths.tdma_c", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = 0;
  __pyx_L0:;
  __Pyx_XDECREF((PyObject *)__pyx_v_x_arr);
  __PYX_XCLEAR_MEMVIEW(&__pyx_v_x, 1);
  __Pyx_XGIVEREF((PyObject *)__pyx_r);
  __Pyx_PyMonitoring_ExitScope(0);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

 025: @cython.wraparound(False)

 026: cdef cnp.ndarray tdma_c(double[::1] a, double[::1] b, double[::1] c, double[::1] d):

 027:     cdef:

+028:         int N=b.shape[0]

  __pyx_v_N = (__pyx_v_b.shape[0]);

 029:         int k

 030:         double m

+031:         cnp.ndarray x_arr=np.empty(N)

  __pyx_t_2 = NULL;
  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 31, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_empty); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 31, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = __Pyx_PyLong_From_int(__pyx_v_N); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 31, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_5 = 1;
  #if CYTHON_UNPACK_METHODS
  if (unlikely(PyMethod_Check(__pyx_t_4))) {
    __pyx_t_2 = PyMethod_GET_SELF(__pyx_t_4);
    assert(__pyx_t_2);
    PyObject* __pyx__function = PyMethod_GET_FUNCTION(__pyx_t_4);
    __Pyx_INCREF(__pyx_t_2);
    __Pyx_INCREF(__pyx__function);
    __Pyx_DECREF_SET(__pyx_t_4, __pyx__function);
    __pyx_t_5 = 0;
  }
  #endif
  {
    PyObject *__pyx_callargs[2] = {__pyx_t_2, __pyx_t_3};
    __pyx_t_1 = __Pyx_PyObject_FastCall((PyObject*)__pyx_t_4, __pyx_callargs+__pyx_t_5, (2-__pyx_t_5) | (__pyx_t_5*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET));
    __Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 31, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
  }
  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_mstate_global->__pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 31, __pyx_L1_error)
  __pyx_v_x_arr = ((PyArrayObject *)__pyx_t_1);
  __pyx_t_1 = 0;

+032:         double[::1] x=x_arr

  __pyx_t_6 = __Pyx_PyObject_to_MemoryviewSlice_dc_double(((PyObject *)__pyx_v_x_arr), PyBUF_WRITABLE); if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 32, __pyx_L1_error)
  __pyx_v_x = __pyx_t_6;
  __pyx_t_6.memview = NULL;
  __pyx_t_6.data = NULL;

 033: 

 034:     # Implementation precisely follows

 035:     # https://www.cfd-online.com/Wiki/Tridiagonal_matrix_algorithm_-_TDMA_(Thomas_algorithm)

 036:     # Except that these arrays are zero-indexed!

 037:     # Makes the same assumption: a[0]=0, c[n-1]=0

+038:     for k in range(1,N):

  __pyx_t_7 = __pyx_v_N;
  __pyx_t_8 = __pyx_t_7;
  for (__pyx_t_9 = 1; __pyx_t_9 < __pyx_t_8; __pyx_t_9+=1) {
    __pyx_v_k = __pyx_t_9;

+039:         m=a[k]/b[k-1]

    __pyx_t_10 = __pyx_v_k;
    __pyx_t_11 = (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_a.data) + __pyx_t_10)) )));
    __pyx_t_10 = (__pyx_v_k - 1);
    __pyx_t_12 = (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_b.data) + __pyx_t_10)) )));
    if (unlikely(__pyx_t_12 == 0)) {
      PyErr_SetString(PyExc_ZeroDivisionError, "float division");
      __PYX_ERR(0, 39, __pyx_L1_error)
    }
    __pyx_v_m = (__pyx_t_11 / __pyx_t_12);

+040:         b[k]-=m*c[k-1]

    __pyx_t_10 = (__pyx_v_k - 1);
    __pyx_t_13 = __pyx_v_k;
    *((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_b.data) + __pyx_t_13)) )) -= (__pyx_v_m * (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_c.data) + __pyx_t_10)) ))));

+041:         d[k]-=m*d[k-1]

    __pyx_t_10 = (__pyx_v_k - 1);
    __pyx_t_13 = __pyx_v_k;
    *((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_d.data) + __pyx_t_13)) )) -= (__pyx_v_m * (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_d.data) + __pyx_t_10)) ))));
  }

+042:     for k in range(N):

  __pyx_t_7 = __pyx_v_N;
  __pyx_t_8 = __pyx_t_7;
  for (__pyx_t_9 = 0; __pyx_t_9 < __pyx_t_8; __pyx_t_9+=1) {
    __pyx_v_k = __pyx_t_9;

+043:         x[k]=d[k]/b[k]

    __pyx_t_10 = __pyx_v_k;
    __pyx_t_12 = (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_d.data) + __pyx_t_10)) )));
    __pyx_t_10 = __pyx_v_k;
    __pyx_t_11 = (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_b.data) + __pyx_t_10)) )));
    if (unlikely(__pyx_t_11 == 0)) {
      PyErr_SetString(PyExc_ZeroDivisionError, "float division");
      __PYX_ERR(0, 43, __pyx_L1_error)
    }
    __pyx_t_10 = __pyx_v_k;
    *((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_x.data) + __pyx_t_10)) )) = (__pyx_t_12 / __pyx_t_11);
  }

+044:     for k in range(N-2,-1,-1):

  for (__pyx_t_7 = (__pyx_v_N - 2); __pyx_t_7 > -1; __pyx_t_7-=1) {
    __pyx_v_k = __pyx_t_7;

+045:         x[k]=(d[k]-c[k]*x[k+1])/b[k]

    __pyx_t_10 = __pyx_v_k;
    __pyx_t_13 = __pyx_v_k;
    __pyx_t_14 = (__pyx_v_k + 1);
    __pyx_t_11 = ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_d.data) + __pyx_t_10)) ))) - ((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_c.data) + __pyx_t_13)) ))) * (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_x.data) + __pyx_t_14)) )))));
    __pyx_t_14 = __pyx_v_k;
    __pyx_t_12 = (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_b.data) + __pyx_t_14)) )));
    if (unlikely(__pyx_t_12 == 0)) {
      PyErr_SetString(PyExc_ZeroDivisionError, "float division");
      __PYX_ERR(0, 45, __pyx_L1_error)
    }
    __pyx_t_14 = __pyx_v_k;
    *((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_x.data) + __pyx_t_14)) )) = (__pyx_t_11 / __pyx_t_12);
  }

+046:     return x_arr

  __Pyx_XDECREF((PyObject *)__pyx_r);
  __Pyx_INCREF((PyObject *)__pyx_v_x_arr);
  __pyx_r = __pyx_v_x_arr;
  __Pyx_TraceReturnValue((PyObject *)__pyx_r, 93, 0, __PYX_ERR(0, 46, __pyx_L1_error));
  goto __pyx_L0;

 047: 

 048: 

+049: def tdma(a,b,c,d,copy=True):

/* Python wrapper */
static PyObject *__pyx_pw_9pynitride_4core_12cython_maths_1tdma(PyObject *__pyx_self, 
#if CYTHON_METH_FASTCALL
PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds
#else
PyObject *__pyx_args, PyObject *__pyx_kwds
#endif
); /*proto*/
PyDoc_STRVAR(__pyx_doc_9pynitride_4core_12cython_maths_tdma, "Implements the tridiagonal matrix algorithm.\n\n    Solves :math:`Mx=d` where :math:`M` is a tridiagonal matrix described by three diagonals :math:`a,b,c`.\n    A good discussion of the algorithm is available at\n    `CFD-Online <https://www.cfd-online.com/Wiki/Tridiagonal_matrix_algorithm_-_TDMA_(Thomas_algorithm)>`_.  All\n    arguments are given as 1-D numpy arrays of the same length.\n\n    :param a: the lower subdiagonal of :math:`M`. This should include a single leading zero so it is the same length as *b*.\n    :param b: the diagonal of :math:`M`.\n    :param c: the upper subdiagonal of :math:`M`. This should include a single trailing zero so it is the same length as *b*.\n    :param d: the desired result of :math:`Mx`.  Should be the same length as :math:`b`.\n    :param copy: If true, arrays will be copied (default); otherwise they will be overwritten.\n\n    :return: :math:`x`, the solution vector as a numpy array.\n    ");
static PyMethodDef __pyx_mdef_9pynitride_4core_12cython_maths_1tdma = {"tdma", (PyCFunction)(void(*)(void))(__Pyx_PyCFunction_FastCallWithKeywords)__pyx_pw_9pynitride_4core_12cython_maths_1tdma, __Pyx_METH_FASTCALL|METH_KEYWORDS, __pyx_doc_9pynitride_4core_12cython_maths_tdma};
static PyObject *__pyx_pw_9pynitride_4core_12cython_maths_1tdma(PyObject *__pyx_self, 
#if CYTHON_METH_FASTCALL
PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds
#else
PyObject *__pyx_args, PyObject *__pyx_kwds
#endif
) {
  PyObject *__pyx_v_a = 0;
  PyObject *__pyx_v_b = 0;
  PyObject *__pyx_v_c = 0;
  PyObject *__pyx_v_d = 0;
  PyObject *__pyx_v_copy = 0;
  #if !CYTHON_METH_FASTCALL
  CYTHON_UNUSED Py_ssize_t __pyx_nargs;
  #endif
  CYTHON_UNUSED PyObject *const *__pyx_kwvalues;
  PyObject *__pyx_r = 0;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("tdma (wrapper)", 0);
  #if !CYTHON_METH_FASTCALL
  #if CYTHON_ASSUME_SAFE_SIZE
  __pyx_nargs = PyTuple_GET_SIZE(__pyx_args);
  #else
  __pyx_nargs = PyTuple_Size(__pyx_args); if (unlikely(__pyx_nargs < 0)) return NULL;
  #endif
  #endif
  __pyx_kwvalues = __Pyx_KwValues_FASTCALL(__pyx_args, __pyx_nargs);
  {
    PyObject ** const __pyx_pyargnames[] = {&__pyx_mstate_global->__pyx_n_u_a,&__pyx_mstate_global->__pyx_n_u_b,&__pyx_mstate_global->__pyx_n_u_c,&__pyx_mstate_global->__pyx_n_u_d,&__pyx_mstate_global->__pyx_n_u_copy,0};
  PyObject* values[5] = {0,0,0,0,0};
    const Py_ssize_t __pyx_kwds_len = (__pyx_kwds) ? __Pyx_NumKwargs_FASTCALL(__pyx_kwds) : 0;
    if (unlikely(__pyx_kwds_len) < 0) __PYX_ERR(0, 49, __pyx_L3_error)
    if (__pyx_kwds_len > 0) {
      switch (__pyx_nargs) {
        case  5:
        values[4] = __Pyx_ArgRef_FASTCALL(__pyx_args, 4);
        if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[4])) __PYX_ERR(0, 49, __pyx_L3_error)
        CYTHON_FALLTHROUGH;
        case  4:
        values[3] = __Pyx_ArgRef_FASTCALL(__pyx_args, 3);
        if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[3])) __PYX_ERR(0, 49, __pyx_L3_error)
        CYTHON_FALLTHROUGH;
        case  3:
        values[2] = __Pyx_ArgRef_FASTCALL(__pyx_args, 2);
        if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[2])) __PYX_ERR(0, 49, __pyx_L3_error)
        CYTHON_FALLTHROUGH;
        case  2:
        values[1] = __Pyx_ArgRef_FASTCALL(__pyx_args, 1);
        if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[1])) __PYX_ERR(0, 49, __pyx_L3_error)
        CYTHON_FALLTHROUGH;
        case  1:
        values[0] = __Pyx_ArgRef_FASTCALL(__pyx_args, 0);
        if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[0])) __PYX_ERR(0, 49, __pyx_L3_error)
        CYTHON_FALLTHROUGH;
        case  0: break;
        default: goto __pyx_L5_argtuple_error;
      }
      const Py_ssize_t kwd_pos_args = __pyx_nargs;
      if (__Pyx_ParseKeywords(__pyx_kwds, __pyx_kwvalues, __pyx_pyargnames, 0, values, kwd_pos_args, __pyx_kwds_len, "tdma", 0) < (0)) __PYX_ERR(0, 49, __pyx_L3_error)
      if (!values[4]) values[4] = __Pyx_NewRef(((PyObject *)((PyObject*)Py_True)));
      for (Py_ssize_t i = __pyx_nargs; i < 4; i++) {
        if (unlikely(!values[i])) { __Pyx_RaiseArgtupleInvalid("tdma", 0, 4, 5, i); __PYX_ERR(0, 49, __pyx_L3_error) }
      }
    } else {
      switch (__pyx_nargs) {
        case  5:
        values[4] = __Pyx_ArgRef_FASTCALL(__pyx_args, 4);
        if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[4])) __PYX_ERR(0, 49, __pyx_L3_error)
        CYTHON_FALLTHROUGH;
        case  4:
        values[3] = __Pyx_ArgRef_FASTCALL(__pyx_args, 3);
        if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[3])) __PYX_ERR(0, 49, __pyx_L3_error)
        values[2] = __Pyx_ArgRef_FASTCALL(__pyx_args, 2);
        if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[2])) __PYX_ERR(0, 49, __pyx_L3_error)
        values[1] = __Pyx_ArgRef_FASTCALL(__pyx_args, 1);
        if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[1])) __PYX_ERR(0, 49, __pyx_L3_error)
        values[0] = __Pyx_ArgRef_FASTCALL(__pyx_args, 0);
        if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[0])) __PYX_ERR(0, 49, __pyx_L3_error)
        break;
        default: goto __pyx_L5_argtuple_error;
      }
      if (!values[4]) values[4] = __Pyx_NewRef(((PyObject *)((PyObject*)Py_True)));
    }
    __pyx_v_a = values[0];
    __pyx_v_b = values[1];
    __pyx_v_c = values[2];
    __pyx_v_d = values[3];
    __pyx_v_copy = values[4];
  }
  goto __pyx_L6_skip;
  __pyx_L5_argtuple_error:;
  __Pyx_RaiseArgtupleInvalid("tdma", 0, 4, 5, __pyx_nargs); __PYX_ERR(0, 49, __pyx_L3_error)
  __pyx_L6_skip:;
  goto __pyx_L4_argument_unpacking_done;
  __pyx_L3_error:;
  for (Py_ssize_t __pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) {
    Py_XDECREF(values[__pyx_temp]);
  }
  __Pyx_AddTraceback("pynitride.core.cython_maths.tdma", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __Pyx_RefNannyFinishContext();
  return NULL;
  __pyx_L4_argument_unpacking_done:;
  __pyx_r = __pyx_pf_9pynitride_4core_12cython_maths_tdma(__pyx_self, __pyx_v_a, __pyx_v_b, __pyx_v_c, __pyx_v_d, __pyx_v_copy);

  /* function exit code */
  for (Py_ssize_t __pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) {
    Py_XDECREF(values[__pyx_temp]);
  }
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

static PyObject *__pyx_pf_9pynitride_4core_12cython_maths_tdma(CYTHON_UNUSED PyObject *__pyx_self, PyObject *__pyx_v_a, PyObject *__pyx_v_b, PyObject *__pyx_v_c, PyObject *__pyx_v_d, PyObject *__pyx_v_copy) {
  Py_ssize_t __pyx_v_N;
  PyObject *__pyx_r = NULL;
  __Pyx_TraceDeclarationsFunc
  __Pyx_TraceFrameInit(((PyObject *)__pyx_mstate_global->__pyx_codeobj_tab[40]))
  __Pyx_TraceStartFunc("tdma", __pyx_f[0], 49, 0, 0, 0, __PYX_ERR(0, 49, __pyx_L1_error));
  __Pyx_INCREF(__pyx_v_a);
  __Pyx_INCREF(__pyx_v_b);
  __Pyx_INCREF(__pyx_v_c);
  __Pyx_INCREF(__pyx_v_d);
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_3);
  __Pyx_XDECREF(__pyx_t_4);
  __PYX_XCLEAR_MEMVIEW(&__pyx_t_6, 1);
  __PYX_XCLEAR_MEMVIEW(&__pyx_t_7, 1);
  __PYX_XCLEAR_MEMVIEW(&__pyx_t_8, 1);
  __PYX_XCLEAR_MEMVIEW(&__pyx_t_9, 1);
  __Pyx_TraceException(__pyx_lineno, 0, 0);
  #if CYTHON_USE_SYS_MONITORING
  __Pyx_TraceExceptionUnwind(0, 0);
  #else
  __Pyx_TraceReturnValue(NULL, 0, 0, __PYX_ERR(0, 49, __pyx_L1_error));
  #endif
  __Pyx_AddTraceback("pynitride.core.cython_maths.tdma", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = NULL;
  __pyx_L0:;
  __Pyx_XDECREF(__pyx_v_a);
  __Pyx_XDECREF(__pyx_v_b);
  __Pyx_XDECREF(__pyx_v_c);
  __Pyx_XDECREF(__pyx_v_d);
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_PyMonitoring_ExitScope(0);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
/* … */
  __pyx_t_4 = __Pyx_CyFunction_New(&__pyx_mdef_9pynitride_4core_12cython_maths_1tdma, 0, __pyx_mstate_global->__pyx_n_u_tdma, NULL, __pyx_mstate_global->__pyx_n_u_pynitride_core_cython_maths, __pyx_mstate_global->__pyx_d, ((PyObject *)__pyx_mstate_global->__pyx_codeobj_tab[40])); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 49, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  #if CYTHON_COMPILING_IN_CPYTHON && PY_VERSION_HEX >= 0x030E0000
  PyUnstable_Object_EnableDeferredRefcount(__pyx_t_4);
  #endif
  __Pyx_CyFunction_SetDefaultsTuple(__pyx_t_4, __pyx_mstate_global->__pyx_tuple[1]);
  if (PyDict_SetItem(__pyx_mstate_global->__pyx_d, __pyx_mstate_global->__pyx_n_u_tdma, __pyx_t_4) < (0)) __PYX_ERR(0, 49, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;

 050:     r"""Implements the tridiagonal matrix algorithm.

 051: 

 052:     Solves :math:`Mx=d` where :math:`M` is a tridiagonal matrix described by three diagonals :math:`a,b,c`.

 053:     A good discussion of the algorithm is available at

 054:     `CFD-Online <https://www.cfd-online.com/Wiki/Tridiagonal_matrix_algorithm_-_TDMA_(Thomas_algorithm)>`_.  All

 055:     arguments are given as 1-D numpy arrays of the same length.

 056: 

 057:     :param a: the lower subdiagonal of :math:`M`. This should include a single leading zero so it is the same length as *b*.

 058:     :param b: the diagonal of :math:`M`.

 059:     :param c: the upper subdiagonal of :math:`M`. This should include a single trailing zero so it is the same length as *b*.

 060:     :param d: the desired result of :math:`Mx`.  Should be the same length as :math:`b`.

 061:     :param copy: If true, arrays will be copied (default); otherwise they will be overwritten.

 062: 

 063:     :return: :math:`x`, the solution vector as a numpy array.

 064:     """

 065: 

+066:     N=len(b)

  __pyx_t_1 = PyObject_Length(__pyx_v_b); if (unlikely(__pyx_t_1 == ((Py_ssize_t)-1))) __PYX_ERR(0, 66, __pyx_L1_error)
  __pyx_v_N = __pyx_t_1;

+067:     assert len(a)==N, 'a should be the same length as b.'

  #ifndef CYTHON_WITHOUT_ASSERTIONS
  if (unlikely(__pyx_assertions_enabled())) {
    __pyx_t_1 = PyObject_Length(__pyx_v_a); if (unlikely(__pyx_t_1 == ((Py_ssize_t)-1))) __PYX_ERR(0, 67, __pyx_L1_error)
    __pyx_t_2 = (__pyx_t_1 == __pyx_v_N);
    if (unlikely(!__pyx_t_2)) {
      __Pyx_Raise(((PyObject *)(((PyTypeObject*)PyExc_AssertionError))), __pyx_mstate_global->__pyx_kp_u_a_should_be_the_same_length_as_b, 0, 0);
      __PYX_ERR(0, 67, __pyx_L1_error)
    }
  }
  #else
  if ((1)); else __PYX_ERR(0, 67, __pyx_L1_error)
  #endif

+068:     assert len(c)==N, 'c should be the same length as b.'

  #ifndef CYTHON_WITHOUT_ASSERTIONS
  if (unlikely(__pyx_assertions_enabled())) {
    __pyx_t_1 = PyObject_Length(__pyx_v_c); if (unlikely(__pyx_t_1 == ((Py_ssize_t)-1))) __PYX_ERR(0, 68, __pyx_L1_error)
    __pyx_t_2 = (__pyx_t_1 == __pyx_v_N);
    if (unlikely(!__pyx_t_2)) {
      __Pyx_Raise(((PyObject *)(((PyTypeObject*)PyExc_AssertionError))), __pyx_mstate_global->__pyx_kp_u_c_should_be_the_same_length_as_b, 0, 0);
      __PYX_ERR(0, 68, __pyx_L1_error)
    }
  }
  #else
  if ((1)); else __PYX_ERR(0, 68, __pyx_L1_error)
  #endif

+069:     assert len(d)==N, 'd should be the same length as b.'

  #ifndef CYTHON_WITHOUT_ASSERTIONS
  if (unlikely(__pyx_assertions_enabled())) {
    __pyx_t_1 = PyObject_Length(__pyx_v_d); if (unlikely(__pyx_t_1 == ((Py_ssize_t)-1))) __PYX_ERR(0, 69, __pyx_L1_error)
    __pyx_t_2 = (__pyx_t_1 == __pyx_v_N);
    if (unlikely(!__pyx_t_2)) {
      __Pyx_Raise(((PyObject *)(((PyTypeObject*)PyExc_AssertionError))), __pyx_mstate_global->__pyx_kp_u_d_should_be_the_same_length_as_b, 0, 0);
      __PYX_ERR(0, 69, __pyx_L1_error)
    }
  }
  #else
  if ((1)); else __PYX_ERR(0, 69, __pyx_L1_error)
  #endif

+070:     if copy:

  __pyx_t_2 = __Pyx_PyObject_IsTrue(__pyx_v_copy); if (unlikely((__pyx_t_2 < 0))) __PYX_ERR(0, 70, __pyx_L1_error)
  if (__pyx_t_2) {
/* … */
  }

+071:         a=a.copy()

    __pyx_t_4 = __pyx_v_a;
    __Pyx_INCREF(__pyx_t_4);
    __pyx_t_5 = 0;
    {
      PyObject *__pyx_callargs[2] = {__pyx_t_4, NULL};
      __pyx_t_3 = __Pyx_PyObject_FastCallMethod((PyObject*)__pyx_mstate_global->__pyx_n_u_copy, __pyx_callargs+__pyx_t_5, (1-__pyx_t_5) | (1*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET));
      __Pyx_XDECREF(__pyx_t_4); __pyx_t_4 = 0;
      if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 71, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_3);
    }
    __Pyx_DECREF_SET(__pyx_v_a, __pyx_t_3);
    __pyx_t_3 = 0;

+072:         b=b.copy()

    __pyx_t_4 = __pyx_v_b;
    __Pyx_INCREF(__pyx_t_4);
    __pyx_t_5 = 0;
    {
      PyObject *__pyx_callargs[2] = {__pyx_t_4, NULL};
      __pyx_t_3 = __Pyx_PyObject_FastCallMethod((PyObject*)__pyx_mstate_global->__pyx_n_u_copy, __pyx_callargs+__pyx_t_5, (1-__pyx_t_5) | (1*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET));
      __Pyx_XDECREF(__pyx_t_4); __pyx_t_4 = 0;
      if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 72, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_3);
    }
    __Pyx_DECREF_SET(__pyx_v_b, __pyx_t_3);
    __pyx_t_3 = 0;

+073:         c=c.copy()

    __pyx_t_4 = __pyx_v_c;
    __Pyx_INCREF(__pyx_t_4);
    __pyx_t_5 = 0;
    {
      PyObject *__pyx_callargs[2] = {__pyx_t_4, NULL};
      __pyx_t_3 = __Pyx_PyObject_FastCallMethod((PyObject*)__pyx_mstate_global->__pyx_n_u_copy, __pyx_callargs+__pyx_t_5, (1-__pyx_t_5) | (1*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET));
      __Pyx_XDECREF(__pyx_t_4); __pyx_t_4 = 0;
      if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 73, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_3);
    }
    __Pyx_DECREF_SET(__pyx_v_c, __pyx_t_3);
    __pyx_t_3 = 0;

+074:         d=d.copy()

    __pyx_t_4 = __pyx_v_d;
    __Pyx_INCREF(__pyx_t_4);
    __pyx_t_5 = 0;
    {
      PyObject *__pyx_callargs[2] = {__pyx_t_4, NULL};
      __pyx_t_3 = __Pyx_PyObject_FastCallMethod((PyObject*)__pyx_mstate_global->__pyx_n_u_copy, __pyx_callargs+__pyx_t_5, (1-__pyx_t_5) | (1*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET));
      __Pyx_XDECREF(__pyx_t_4); __pyx_t_4 = 0;
      if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 74, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_3);
    }
    __Pyx_DECREF_SET(__pyx_v_d, __pyx_t_3);
    __pyx_t_3 = 0;

+075:     return tdma_c(a,b,c,d)

  __Pyx_XDECREF(__pyx_r);
  __pyx_t_6 = __Pyx_PyObject_to_MemoryviewSlice_dc_double(__pyx_v_a, PyBUF_WRITABLE); if (unlikely(!__pyx_t_6.memview)) __PYX_ERR(0, 75, __pyx_L1_error)
  __pyx_t_7 = __Pyx_PyObject_to_MemoryviewSlice_dc_double(__pyx_v_b, PyBUF_WRITABLE); if (unlikely(!__pyx_t_7.memview)) __PYX_ERR(0, 75, __pyx_L1_error)
  __pyx_t_8 = __Pyx_PyObject_to_MemoryviewSlice_dc_double(__pyx_v_c, PyBUF_WRITABLE); if (unlikely(!__pyx_t_8.memview)) __PYX_ERR(0, 75, __pyx_L1_error)
  __pyx_t_9 = __Pyx_PyObject_to_MemoryviewSlice_dc_double(__pyx_v_d, PyBUF_WRITABLE); if (unlikely(!__pyx_t_9.memview)) __PYX_ERR(0, 75, __pyx_L1_error)
  __pyx_t_3 = ((PyObject *)__pyx_f_9pynitride_4core_12cython_maths_tdma_c(__pyx_t_6, __pyx_t_7, __pyx_t_8, __pyx_t_9)); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 75, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __PYX_XCLEAR_MEMVIEW(&__pyx_t_6, 1);
  __pyx_t_6.memview = NULL; __pyx_t_6.data = NULL;
  __PYX_XCLEAR_MEMVIEW(&__pyx_t_7, 1);
  __pyx_t_7.memview = NULL; __pyx_t_7.data = NULL;
  __PYX_XCLEAR_MEMVIEW(&__pyx_t_8, 1);
  __pyx_t_8.memview = NULL; __pyx_t_8.data = NULL;
  __PYX_XCLEAR_MEMVIEW(&__pyx_t_9, 1);
  __pyx_t_9.memview = NULL; __pyx_t_9.data = NULL;
  __pyx_r = __pyx_t_3;
  __pyx_t_3 = 0;
  __Pyx_TraceReturnValue(__pyx_r, 45, 0, __PYX_ERR(0, 75, __pyx_L1_error));
  goto __pyx_L0;

 076: 

 077: ##########

 078: ### Fermi-Dirac integrals

 079: ##########

 080: 

 081: 

 082: # These values come from Table 6 for use in the Fermi-Dirac 1/2 integrals

 083: # Wong et al Solid-State Electronics Vol. 37, No. I, pp. 61~64, 1994

 084: # http://dx.doi.org/10.1016/0038-1101(94)90105-8

 085: DEF ORDER=7

 086: cdef double[ORDER] a, b1, b2, c, ap, b1p, b2p, cp

+087: a[:]=[1,.353568,.192439,.122973,.077134,.036228,.008346]

  __pyx_t_11[0] = 1.0;
  __pyx_t_11[1] = .353568;
  __pyx_t_11[2] = .192439;
  __pyx_t_11[3] = .122973;
  __pyx_t_11[4] = .077134;
  __pyx_t_11[5] = .036228;
  __pyx_t_11[6] = .008346;
  memcpy(&(__pyx_v_9pynitride_4core_12cython_maths_a[0]), __pyx_t_11, sizeof(__pyx_v_9pynitride_4core_12cython_maths_a[0]) * (7));

+088: b1[:]=[.76514793,.60488667,.19003355,2.00193968e-2,-4.12643816e-3,-4.70958992e-4,1.50071469e-4]

  __pyx_t_12[0] = .76514793;
  __pyx_t_12[1] = .60488667;
  __pyx_t_12[2] = .19003355;
  __pyx_t_12[3] = 2.00193968e-2;
  __pyx_t_12[4] = -4.12643816e-3;
  __pyx_t_12[5] = -4.70958992e-4;
  __pyx_t_12[6] = 1.50071469e-4;
  memcpy(&(__pyx_v_9pynitride_4core_12cython_maths_b1[0]), __pyx_t_12, sizeof(__pyx_v_9pynitride_4core_12cython_maths_b1[0]) * (7));

+089: b2[:]=[.78095732,.57254453,.21419339,1.38432741e-2,-5.54949386e-3,6.48814900e-4,-2.84050520e-5]

  __pyx_t_13[0] = .78095732;
  __pyx_t_13[1] = .57254453;
  __pyx_t_13[2] = .21419339;
  __pyx_t_13[3] = 1.38432741e-2;
  __pyx_t_13[4] = -5.54949386e-3;
  __pyx_t_13[5] = 6.48814900e-4;
  __pyx_t_13[6] = -2.84050520e-5;
  memcpy(&(__pyx_v_9pynitride_4core_12cython_maths_b2[0]), __pyx_t_13, sizeof(__pyx_v_9pynitride_4core_12cython_maths_b2[0]) * (7));

+090: c[:]=[.752253,.928195,.680839,25.7829,-553.636,3531.43,-3254.65]

  __pyx_t_14[0] = .752253;
  __pyx_t_14[1] = .928195;
  __pyx_t_14[2] = .680839;
  __pyx_t_14[3] = 25.7829;
  __pyx_t_14[4] = -553.636;
  __pyx_t_14[5] = 3531.43;
  __pyx_t_14[6] = -3254.65;
  memcpy(&(__pyx_v_9pynitride_4core_12cython_maths_c[0]), __pyx_t_14, sizeof(__pyx_v_9pynitride_4core_12cython_maths_c[0]) * (7));

 091: 

 092: # These values will be used by the Fermi-Dirac 1/2 Derivative function

 093: # They come from modifying the above values by differentiating each expression in piecewise

 094: # from Equations 24-26 in the same reference.

+095: for j in range(ORDER):

  for (__pyx_t_9 = 0; __pyx_t_9 < 7; __pyx_t_9+=1) {
    __pyx_t_4 = PyLong_FromSsize_t(__pyx_t_9); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 95, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_4);
    if (PyDict_SetItem(__pyx_mstate_global->__pyx_d, __pyx_mstate_global->__pyx_n_u_j, __pyx_t_4) < (0)) __PYX_ERR(0, 95, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;

+096:     ap[j] =a[j]*(j+1)

    __Pyx_GetModuleGlobalName(__pyx_t_4, __pyx_mstate_global->__pyx_n_u_j); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 96, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_4);
    __pyx_t_15 = __Pyx_PyIndex_AsSsize_t(__pyx_t_4); if (unlikely((__pyx_t_15 == (Py_ssize_t)-1) && PyErr_Occurred())) __PYX_ERR(0, 96, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    __pyx_t_4 = PyFloat_FromDouble((__pyx_v_9pynitride_4core_12cython_maths_a[__pyx_t_15])); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 96, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_4);
    __Pyx_GetModuleGlobalName(__pyx_t_5, __pyx_mstate_global->__pyx_n_u_j); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 96, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_5);
    __pyx_t_16 = __Pyx_PyLong_AddObjC(__pyx_t_5, __pyx_mstate_global->__pyx_int_1, 1, 0, 0); if (unlikely(!__pyx_t_16)) __PYX_ERR(0, 96, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_16);
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    __pyx_t_5 = PyNumber_Multiply(__pyx_t_4, __pyx_t_16); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 96, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_5);
    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    __Pyx_DECREF(__pyx_t_16); __pyx_t_16 = 0;
    __pyx_t_17 = __Pyx_PyFloat_AsDouble(__pyx_t_5); if (unlikely((__pyx_t_17 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 96, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    __Pyx_GetModuleGlobalName(__pyx_t_5, __pyx_mstate_global->__pyx_n_u_j); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 96, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_5);
    __pyx_t_15 = __Pyx_PyIndex_AsSsize_t(__pyx_t_5); if (unlikely((__pyx_t_15 == (Py_ssize_t)-1) && PyErr_Occurred())) __PYX_ERR(0, 96, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    (__pyx_v_9pynitride_4core_12cython_maths_ap[__pyx_t_15]) = __pyx_t_17;

+097:     b1p[ORDER-j-1]=(ORDER-j-1)*b1[ORDER-j-1]

    __Pyx_GetModuleGlobalName(__pyx_t_5, __pyx_mstate_global->__pyx_n_u_j); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 97, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_5);
    __pyx_t_16 = __Pyx_PyLong_SubtractCObj(__pyx_mstate_global->__pyx_int_7, __pyx_t_5, 7, 0, 0); if (unlikely(!__pyx_t_16)) __PYX_ERR(0, 97, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_16);
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    __pyx_t_5 = __Pyx_PyLong_SubtractObjC(__pyx_t_16, __pyx_mstate_global->__pyx_int_1, 1, 0, 0); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 97, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_5);
    __Pyx_DECREF(__pyx_t_16); __pyx_t_16 = 0;
    __Pyx_GetModuleGlobalName(__pyx_t_16, __pyx_mstate_global->__pyx_n_u_j); if (unlikely(!__pyx_t_16)) __PYX_ERR(0, 97, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_16);
    __pyx_t_4 = __Pyx_PyLong_SubtractCObj(__pyx_mstate_global->__pyx_int_7, __pyx_t_16, 7, 0, 0); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 97, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_4);
    __Pyx_DECREF(__pyx_t_16); __pyx_t_16 = 0;
    __pyx_t_16 = __Pyx_PyLong_SubtractObjC(__pyx_t_4, __pyx_mstate_global->__pyx_int_1, 1, 0, 0); if (unlikely(!__pyx_t_16)) __PYX_ERR(0, 97, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_16);
    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    __pyx_t_15 = __Pyx_PyIndex_AsSsize_t(__pyx_t_16); if (unlikely((__pyx_t_15 == (Py_ssize_t)-1) && PyErr_Occurred())) __PYX_ERR(0, 97, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_16); __pyx_t_16 = 0;
    __pyx_t_16 = PyFloat_FromDouble((__pyx_v_9pynitride_4core_12cython_maths_b1[__pyx_t_15])); if (unlikely(!__pyx_t_16)) __PYX_ERR(0, 97, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_16);
    __pyx_t_4 = PyNumber_Multiply(__pyx_t_5, __pyx_t_16); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 97, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_4);
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    __Pyx_DECREF(__pyx_t_16); __pyx_t_16 = 0;
    __pyx_t_17 = __Pyx_PyFloat_AsDouble(__pyx_t_4); if (unlikely((__pyx_t_17 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 97, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    __Pyx_GetModuleGlobalName(__pyx_t_4, __pyx_mstate_global->__pyx_n_u_j); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 97, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_4);
    __pyx_t_16 = __Pyx_PyLong_SubtractCObj(__pyx_mstate_global->__pyx_int_7, __pyx_t_4, 7, 0, 0); if (unlikely(!__pyx_t_16)) __PYX_ERR(0, 97, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_16);
    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    __pyx_t_4 = __Pyx_PyLong_SubtractObjC(__pyx_t_16, __pyx_mstate_global->__pyx_int_1, 1, 0, 0); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 97, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_4);
    __Pyx_DECREF(__pyx_t_16); __pyx_t_16 = 0;
    __pyx_t_15 = __Pyx_PyIndex_AsSsize_t(__pyx_t_4); if (unlikely((__pyx_t_15 == (Py_ssize_t)-1) && PyErr_Occurred())) __PYX_ERR(0, 97, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    (__pyx_v_9pynitride_4core_12cython_maths_b1p[__pyx_t_15]) = __pyx_t_17;

+098:     b2p[ORDER-j-1]=(ORDER-j-1)*b2[ORDER-j-1]

    __Pyx_GetModuleGlobalName(__pyx_t_4, __pyx_mstate_global->__pyx_n_u_j); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 98, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_4);
    __pyx_t_16 = __Pyx_PyLong_SubtractCObj(__pyx_mstate_global->__pyx_int_7, __pyx_t_4, 7, 0, 0); if (unlikely(!__pyx_t_16)) __PYX_ERR(0, 98, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_16);
    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    __pyx_t_4 = __Pyx_PyLong_SubtractObjC(__pyx_t_16, __pyx_mstate_global->__pyx_int_1, 1, 0, 0); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 98, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_4);
    __Pyx_DECREF(__pyx_t_16); __pyx_t_16 = 0;
    __Pyx_GetModuleGlobalName(__pyx_t_16, __pyx_mstate_global->__pyx_n_u_j); if (unlikely(!__pyx_t_16)) __PYX_ERR(0, 98, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_16);
    __pyx_t_5 = __Pyx_PyLong_SubtractCObj(__pyx_mstate_global->__pyx_int_7, __pyx_t_16, 7, 0, 0); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 98, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_5);
    __Pyx_DECREF(__pyx_t_16); __pyx_t_16 = 0;
    __pyx_t_16 = __Pyx_PyLong_SubtractObjC(__pyx_t_5, __pyx_mstate_global->__pyx_int_1, 1, 0, 0); if (unlikely(!__pyx_t_16)) __PYX_ERR(0, 98, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_16);
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    __pyx_t_15 = __Pyx_PyIndex_AsSsize_t(__pyx_t_16); if (unlikely((__pyx_t_15 == (Py_ssize_t)-1) && PyErr_Occurred())) __PYX_ERR(0, 98, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_16); __pyx_t_16 = 0;
    __pyx_t_16 = PyFloat_FromDouble((__pyx_v_9pynitride_4core_12cython_maths_b2[__pyx_t_15])); if (unlikely(!__pyx_t_16)) __PYX_ERR(0, 98, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_16);
    __pyx_t_5 = PyNumber_Multiply(__pyx_t_4, __pyx_t_16); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 98, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_5);
    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    __Pyx_DECREF(__pyx_t_16); __pyx_t_16 = 0;
    __pyx_t_17 = __Pyx_PyFloat_AsDouble(__pyx_t_5); if (unlikely((__pyx_t_17 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 98, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    __Pyx_GetModuleGlobalName(__pyx_t_5, __pyx_mstate_global->__pyx_n_u_j); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 98, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_5);
    __pyx_t_16 = __Pyx_PyLong_SubtractCObj(__pyx_mstate_global->__pyx_int_7, __pyx_t_5, 7, 0, 0); if (unlikely(!__pyx_t_16)) __PYX_ERR(0, 98, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_16);
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    __pyx_t_5 = __Pyx_PyLong_SubtractObjC(__pyx_t_16, __pyx_mstate_global->__pyx_int_1, 1, 0, 0); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 98, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_5);
    __Pyx_DECREF(__pyx_t_16); __pyx_t_16 = 0;
    __pyx_t_15 = __Pyx_PyIndex_AsSsize_t(__pyx_t_5); if (unlikely((__pyx_t_15 == (Py_ssize_t)-1) && PyErr_Occurred())) __PYX_ERR(0, 98, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    (__pyx_v_9pynitride_4core_12cython_maths_b2p[__pyx_t_15]) = __pyx_t_17;

+099:     cp[j] =(1.5-2*j)*c[j]

    __Pyx_GetModuleGlobalName(__pyx_t_5, __pyx_mstate_global->__pyx_n_u_j); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 99, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_5);
    __pyx_t_16 = __Pyx_PyLong_MultiplyCObj(__pyx_mstate_global->__pyx_int_2, __pyx_t_5, 2, 0, 0); if (unlikely(!__pyx_t_16)) __PYX_ERR(0, 99, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_16);
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    __pyx_t_5 = __Pyx_PyFloat_SubtractCObj(__pyx_mstate_global->__pyx_float_1_5, __pyx_t_16, 1.5, 0, 0); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 99, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_5);
    __Pyx_DECREF(__pyx_t_16); __pyx_t_16 = 0;
    __Pyx_GetModuleGlobalName(__pyx_t_16, __pyx_mstate_global->__pyx_n_u_j); if (unlikely(!__pyx_t_16)) __PYX_ERR(0, 99, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_16);
    __pyx_t_15 = __Pyx_PyIndex_AsSsize_t(__pyx_t_16); if (unlikely((__pyx_t_15 == (Py_ssize_t)-1) && PyErr_Occurred())) __PYX_ERR(0, 99, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_16); __pyx_t_16 = 0;
    __pyx_t_16 = PyFloat_FromDouble((__pyx_v_9pynitride_4core_12cython_maths_c[__pyx_t_15])); if (unlikely(!__pyx_t_16)) __PYX_ERR(0, 99, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_16);
    __pyx_t_4 = PyNumber_Multiply(__pyx_t_5, __pyx_t_16); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 99, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_4);
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    __Pyx_DECREF(__pyx_t_16); __pyx_t_16 = 0;
    __pyx_t_17 = __Pyx_PyFloat_AsDouble(__pyx_t_4); if (unlikely((__pyx_t_17 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 99, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    __Pyx_GetModuleGlobalName(__pyx_t_4, __pyx_mstate_global->__pyx_n_u_j); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 99, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_4);
    __pyx_t_15 = __Pyx_PyIndex_AsSsize_t(__pyx_t_4); if (unlikely((__pyx_t_15 == (Py_ssize_t)-1) && PyErr_Occurred())) __PYX_ERR(0, 99, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    (__pyx_v_9pynitride_4core_12cython_maths_cp[__pyx_t_15]) = __pyx_t_17;
  }

 100: 

 101: 

 102: # c-function to implement Fermi-Dirac 1/2 integral.  See :py:func:`~pynitride.poissolve.maths.fd12` for more.

+103: cdef double fd12_scalar(double x):

static double __pyx_f_9pynitride_4core_12cython_maths_fd12_scalar(double __pyx_v_x) {
  double __pyx_v_partial_sum;
  long __pyx_v_j;
  long __pyx_v_s;
  double __pyx_r;
  __Pyx_TraceDeclarationsFunc
  __Pyx_TraceFrameInit(((PyObject *)__pyx_mstate_global->__pyx_codeobj_tab[41]))
  __Pyx_TraceStartFunc("fd12_scalar", __pyx_f[0], 103, 0, 0, 0, __PYX_ERR(0, 103, __pyx_L1_error));
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_TraceException(__pyx_lineno, 0, 0);
  #if CYTHON_USE_SYS_MONITORING
  __Pyx_TraceExceptionUnwind(0, 0);
  #else
  __Pyx_TraceReturnValue(NULL, 0, 0, __PYX_ERR(0, 103, __pyx_L1_error));
  #endif
  __Pyx_AddTraceback("pynitride.core.cython_maths.fd12_scalar", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = -1;
  __pyx_L0:;
  __Pyx_PyMonitoring_ExitScope(0);
  return __pyx_r;
}

 104: 

 105:     cdef:

 106:         double partial_sum

 107:         long j

 108:         long s

 109: 

 110:     # Follows the sums from Wong et al (see documentation for fd12)

+111:     partial_sum=0

  __pyx_v_partial_sum = 0.0;

+112:     if x<=-10:

  __pyx_t_1 = (__pyx_v_x <= -10.0);
  if (__pyx_t_1) {
/* … */
    goto __pyx_L3;
  }

+113:         partial_sum=exp(x)

    __pyx_v_partial_sum = exp(__pyx_v_x);

+114:     elif x<=0:

  __pyx_t_1 = (__pyx_v_x <= 0.0);
  if (__pyx_t_1) {
/* … */
    goto __pyx_L3;
  }

+115:         s=1

    __pyx_v_s = 1;

+116:         for j in range(ORDER):

    for (__pyx_t_2 = 0; __pyx_t_2 < 7; __pyx_t_2+=1) {
      __pyx_v_j = __pyx_t_2;

+117:             partial_sum+=s*a[j]*exp((j+1)*x)

      __pyx_v_partial_sum = (__pyx_v_partial_sum + ((__pyx_v_s * (__pyx_v_9pynitride_4core_12cython_maths_a[__pyx_v_j])) * exp(((__pyx_v_j + 1) * __pyx_v_x))));

+118:             s*=-1

      __pyx_v_s = (__pyx_v_s * -1L);
    }

+119:     elif x<=2:

  __pyx_t_1 = (__pyx_v_x <= 2.0);
  if (__pyx_t_1) {
/* … */
    goto __pyx_L3;
  }

+120:         partial_sum=b1[ORDER-1]

    __pyx_v_partial_sum = (__pyx_v_9pynitride_4core_12cython_maths_b1[0x6]);

+121:         for j in range(1,ORDER):

    for (__pyx_t_2 = 1; __pyx_t_2 < 7; __pyx_t_2+=1) {
      __pyx_v_j = __pyx_t_2;

+122:             partial_sum=partial_sum*x + b1[ORDER-j-1]

      __pyx_v_partial_sum = ((__pyx_v_partial_sum * __pyx_v_x) + (__pyx_v_9pynitride_4core_12cython_maths_b1[((7 - __pyx_v_j) - 1)]));
    }

+123:     elif x<=5:

  __pyx_t_1 = (__pyx_v_x <= 5.0);
  if (__pyx_t_1) {
/* … */
    goto __pyx_L3;
  }

+124:         partial_sum=b2[ORDER-1]

    __pyx_v_partial_sum = (__pyx_v_9pynitride_4core_12cython_maths_b2[0x6]);

+125:         for j in range(1,ORDER):

    for (__pyx_t_2 = 1; __pyx_t_2 < 7; __pyx_t_2+=1) {
      __pyx_v_j = __pyx_t_2;

+126:             partial_sum=partial_sum*x + b2[ORDER-j-1]

      __pyx_v_partial_sum = ((__pyx_v_partial_sum * __pyx_v_x) + (__pyx_v_9pynitride_4core_12cython_maths_b2[((7 - __pyx_v_j) - 1)]));
    }

 127:     else:

+128:         for j in range(ORDER):

  /*else*/ {
    for (__pyx_t_2 = 0; __pyx_t_2 < 7; __pyx_t_2+=1) {
      __pyx_v_j = __pyx_t_2;

+129:             partial_sum+=c[j]*pow(x,1.5-2*j)

      __pyx_v_partial_sum = (__pyx_v_partial_sum + ((__pyx_v_9pynitride_4core_12cython_maths_c[__pyx_v_j]) * pow(__pyx_v_x, (1.5 - (2 * __pyx_v_j)))));
    }
  }
  __pyx_L3:;

+130:     return partial_sum

  __pyx_r = __pyx_v_partial_sum;
  __Pyx_TraceReturnCValue(__pyx_r, PyFloat_FromDouble, 106, 0, __PYX_ERR(0, 130, __pyx_L1_error));
  goto __pyx_L0;

 131: 

 132: 

 133: # c-function to implement Fermi-Dirac 1/2 Integral Derivative.  See :py:func:`~pynitride.poissolve.maths.fd12p` for more.

+134: cdef double fd12p_scalar(double x):

static double __pyx_f_9pynitride_4core_12cython_maths_fd12p_scalar(double __pyx_v_x) {
  double __pyx_v_partial_sum;
  long __pyx_v_j;
  long __pyx_v_s;
  double __pyx_r;
  __Pyx_TraceDeclarationsFunc
  __Pyx_TraceFrameInit(((PyObject *)__pyx_mstate_global->__pyx_codeobj_tab[42]))
  __Pyx_TraceStartFunc("fd12p_scalar", __pyx_f[0], 134, 0, 0, 0, __PYX_ERR(0, 134, __pyx_L1_error));
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_TraceException(__pyx_lineno, 0, 0);
  #if CYTHON_USE_SYS_MONITORING
  __Pyx_TraceExceptionUnwind(0, 0);
  #else
  __Pyx_TraceReturnValue(NULL, 0, 0, __PYX_ERR(0, 134, __pyx_L1_error));
  #endif
  __Pyx_AddTraceback("pynitride.core.cython_maths.fd12p_scalar", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = -1;
  __pyx_L0:;
  __Pyx_PyMonitoring_ExitScope(0);
  return __pyx_r;
}

 135:     cdef:

 136:         double partial_sum

 137:         long j

 138:         long s

+139:     partial_sum=0

  __pyx_v_partial_sum = 0.0;

 140: 

 141:     # Code is the same as in fd12_scalar except using the primed coefficients,

 142:     # and changing the power in the final sum

+143:     if x<=-10:

  __pyx_t_1 = (__pyx_v_x <= -10.0);
  if (__pyx_t_1) {
/* … */
    goto __pyx_L3;
  }

+144:         partial_sum=exp(x)

    __pyx_v_partial_sum = exp(__pyx_v_x);

+145:     elif x<=0:

  __pyx_t_1 = (__pyx_v_x <= 0.0);
  if (__pyx_t_1) {
/* … */
    goto __pyx_L3;
  }

+146:         s=1

    __pyx_v_s = 1;

+147:         for j in range(ORDER):

    for (__pyx_t_2 = 0; __pyx_t_2 < 7; __pyx_t_2+=1) {
      __pyx_v_j = __pyx_t_2;

+148:             partial_sum+=s*ap[j]*exp((j+1)*x)

      __pyx_v_partial_sum = (__pyx_v_partial_sum + ((__pyx_v_s * (__pyx_v_9pynitride_4core_12cython_maths_ap[__pyx_v_j])) * exp(((__pyx_v_j + 1) * __pyx_v_x))));

+149:             s*=-1

      __pyx_v_s = (__pyx_v_s * -1L);
    }

+150:     elif x<=2:

  __pyx_t_1 = (__pyx_v_x <= 2.0);
  if (__pyx_t_1) {
/* … */
    goto __pyx_L3;
  }

+151:         partial_sum=b1p[ORDER-1]

    __pyx_v_partial_sum = (__pyx_v_9pynitride_4core_12cython_maths_b1p[0x6]);

+152:         for j in range(1,ORDER-1):

    for (__pyx_t_2 = 1; __pyx_t_2 < 0x6; __pyx_t_2+=1) {
      __pyx_v_j = __pyx_t_2;

+153:             partial_sum=partial_sum*x + b1p[ORDER-j-1]

      __pyx_v_partial_sum = ((__pyx_v_partial_sum * __pyx_v_x) + (__pyx_v_9pynitride_4core_12cython_maths_b1p[((7 - __pyx_v_j) - 1)]));
    }

+154:     elif x<=5:

  __pyx_t_1 = (__pyx_v_x <= 5.0);
  if (__pyx_t_1) {
/* … */
    goto __pyx_L3;
  }

+155:         partial_sum=b2p[ORDER-1]

    __pyx_v_partial_sum = (__pyx_v_9pynitride_4core_12cython_maths_b2p[0x6]);

+156:         for j in range(1,ORDER-1):

    for (__pyx_t_2 = 1; __pyx_t_2 < 0x6; __pyx_t_2+=1) {
      __pyx_v_j = __pyx_t_2;

+157:             partial_sum=partial_sum*x + b2p[ORDER-j-1]

      __pyx_v_partial_sum = ((__pyx_v_partial_sum * __pyx_v_x) + (__pyx_v_9pynitride_4core_12cython_maths_b2p[((7 - __pyx_v_j) - 1)]));
    }

 158:     else:

+159:         for j in range(ORDER):

  /*else*/ {
    for (__pyx_t_2 = 0; __pyx_t_2 < 7; __pyx_t_2+=1) {
      __pyx_v_j = __pyx_t_2;

+160:             partial_sum+=cp[j]*pow(x,.5-2*j)

      __pyx_v_partial_sum = (__pyx_v_partial_sum + ((__pyx_v_9pynitride_4core_12cython_maths_cp[__pyx_v_j]) * pow(__pyx_v_x, (.5 - (2 * __pyx_v_j)))));
    }
  }
  __pyx_L3:;

+161:     return partial_sum

  __pyx_r = __pyx_v_partial_sum;
  __Pyx_TraceReturnCValue(__pyx_r, PyFloat_FromDouble, 106, 0, __PYX_ERR(0, 161, __pyx_L1_error));
  goto __pyx_L0;

 162: 

+163: @cython.boundscheck(False)

static PyArrayObject *__pyx_f_9pynitride_4core_12cython_maths_map1(PyArrayObject *__pyx_v_inarr, double (*__pyx_v_func)(double), struct __pyx_opt_args_9pynitride_4core_12cython_maths_map1 *__pyx_optional_args) {
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_1);
  __Pyx_XDECREF(__pyx_t_2);
  __Pyx_XDECREF(__pyx_t_3);
  __Pyx_XDECREF(__pyx_t_5);
  __Pyx_XDECREF(__pyx_t_7);
  __Pyx_XDECREF(__pyx_t_8);
  __Pyx_XDECREF(__pyx_t_9);
  __Pyx_XDECREF(__pyx_t_10);
  { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
    __Pyx_PyThreadState_declare
    __Pyx_PyThreadState_assign
    __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_subinarr.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_suboutarr.rcbuffer->pybuffer);
  __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
  __Pyx_TraceException(__pyx_lineno, 0, 0);
  #if CYTHON_USE_SYS_MONITORING
  __Pyx_TraceExceptionUnwind(0, 0);
  #else
  __Pyx_TraceReturnValue(NULL, 0, 0, __PYX_ERR(0, 163, __pyx_L1_error));
  #endif
  __Pyx_AddTraceback("pynitride.core.cython_maths.map1", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = 0;
  goto __pyx_L2;
  __pyx_L0:;
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_subinarr.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_suboutarr.rcbuffer->pybuffer);
  __pyx_L2:;
  __Pyx_XDECREF((PyObject *)__pyx_v_subinarr);
  __Pyx_XDECREF((PyObject *)__pyx_v_suboutarr);
  __Pyx_XDECREF(__pyx_v_it);
  __Pyx_XDECREF((PyObject *)__pyx_v_outarr);
  __Pyx_XGIVEREF((PyObject *)__pyx_r);
  __Pyx_PyMonitoring_ExitScope(0);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

+164: cdef cnp.ndarray map1(  cnp.ndarray inarr,

struct __pyx_opt_args_9pynitride_4core_12cython_maths_map1 {
  int __pyx_n;
  PyArrayObject *outarr;
};

 165:                         double (*func)(double),

+166:                         cnp.ndarray outarr=None):

  PyArrayObject *__pyx_v_outarr = ((PyArrayObject *)Py_None);
  int __pyx_v_j;
  PyArrayObject *__pyx_v_subinarr = 0;
  PyArrayObject *__pyx_v_suboutarr = 0;
  PyObject *__pyx_v_it = NULL;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_subinarr;
  __Pyx_Buffer __pyx_pybuffer_subinarr;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_suboutarr;
  __Pyx_Buffer __pyx_pybuffer_suboutarr;
  PyArrayObject *__pyx_r = NULL;
  __Pyx_TraceDeclarationsFunc
  __Pyx_TraceFrameInit(((PyObject *)__pyx_mstate_global->__pyx_codeobj_tab[43]))
  __Pyx_TraceStartFunc("map1", __pyx_f[0], 163, 0, 0, 0, __PYX_ERR(0, 163, __pyx_L1_error));
  if (__pyx_optional_args) {
    if (__pyx_optional_args->__pyx_n > 0) {
      __pyx_v_outarr = __pyx_optional_args->outarr;
    }
  }
  __Pyx_INCREF((PyObject *)__pyx_v_outarr);
  __pyx_pybuffer_subinarr.pybuffer.buf = NULL;
  __pyx_pybuffer_subinarr.refcount = 0;
  __pyx_pybuffernd_subinarr.data = NULL;
  __pyx_pybuffernd_subinarr.rcbuffer = &__pyx_pybuffer_subinarr;
  __pyx_pybuffer_suboutarr.pybuffer.buf = NULL;
  __pyx_pybuffer_suboutarr.refcount = 0;
  __pyx_pybuffernd_suboutarr.data = NULL;
  __pyx_pybuffernd_suboutarr.rcbuffer = &__pyx_pybuffer_suboutarr;

 167:     r""" Map the function `func` over `inarr` to fill `outarr`.

 168: 

 169:     Args:

 170:         inarr: double array of inputs

 171:         func: a function to apply to each element of the input

 172:             takes a single `double` as an argument

 173:         outarr: optional double array of outputs to fill

 174:     Returns:

 175:         the double array of outputs

 176:     """

 177:     cdef int j

 178:     cdef cnp.ndarray[double] subinarr, suboutarr

 179: 

+180:     assert inarr.dtype==np.float64

  #ifndef CYTHON_WITHOUT_ASSERTIONS
  if (unlikely(__pyx_assertions_enabled())) {
    __pyx_t_1 = __Pyx_PyObject_GetAttrStr(((PyObject *)__pyx_v_inarr), __pyx_mstate_global->__pyx_n_u_dtype); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 180, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 180, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_2);
    __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_mstate_global->__pyx_n_u_float64); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 180, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_3);
    __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
    __pyx_t_2 = PyObject_RichCompare(__pyx_t_1, __pyx_t_3, Py_EQ); __Pyx_XGOTREF(__pyx_t_2); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 180, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    __pyx_t_4 = __Pyx_PyObject_IsTrue(__pyx_t_2); if (unlikely((__pyx_t_4 < 0))) __PYX_ERR(0, 180, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
    if (unlikely(!__pyx_t_4)) {
      __Pyx_Raise(((PyObject *)(((PyTypeObject*)PyExc_AssertionError))), 0, 0, 0);
      __PYX_ERR(0, 180, __pyx_L1_error)
    }
  }
  #else
  if ((1)); else __PYX_ERR(0, 180, __pyx_L1_error)
  #endif

+181:     if outarr is None:

  __pyx_t_4 = (((PyObject *)__pyx_v_outarr) == Py_None);
  if (__pyx_t_4) {
/* … */
  }

+182:         outarr=np.empty_like(inarr)

    __pyx_t_3 = NULL;
    __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 182, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    __pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_mstate_global->__pyx_n_u_empty_like); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 182, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_5);
    __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
    __pyx_t_6 = 1;
    #if CYTHON_UNPACK_METHODS
    if (unlikely(PyMethod_Check(__pyx_t_5))) {
      __pyx_t_3 = PyMethod_GET_SELF(__pyx_t_5);
      assert(__pyx_t_3);
      PyObject* __pyx__function = PyMethod_GET_FUNCTION(__pyx_t_5);
      __Pyx_INCREF(__pyx_t_3);
      __Pyx_INCREF(__pyx__function);
      __Pyx_DECREF_SET(__pyx_t_5, __pyx__function);
      __pyx_t_6 = 0;
    }
    #endif
    {
      PyObject *__pyx_callargs[2] = {__pyx_t_3, ((PyObject *)__pyx_v_inarr)};
      __pyx_t_2 = __Pyx_PyObject_FastCall((PyObject*)__pyx_t_5, __pyx_callargs+__pyx_t_6, (2-__pyx_t_6) | (__pyx_t_6*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET));
      __Pyx_XDECREF(__pyx_t_3); __pyx_t_3 = 0;
      __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
      if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 182, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_2);
    }
    if (!(likely(((__pyx_t_2) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_2, __pyx_mstate_global->__pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 182, __pyx_L1_error)
    __Pyx_DECREF_SET(__pyx_v_outarr, ((PyArrayObject *)__pyx_t_2));
    __pyx_t_2 = 0;

 183: 

+184:     it = np.nditer([inarr,outarr], flags=['external_loop','buffered'],

  __pyx_t_5 = NULL;
  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 184, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_nditer); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 184, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = PyList_New(2); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 184, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_INCREF((PyObject *)__pyx_v_inarr);
  __Pyx_GIVEREF((PyObject *)__pyx_v_inarr);
  if (__Pyx_PyList_SET_ITEM(__pyx_t_3, 0, ((PyObject *)__pyx_v_inarr)) != (0)) __PYX_ERR(0, 184, __pyx_L1_error);
  __Pyx_INCREF((PyObject *)__pyx_v_outarr);
  __Pyx_GIVEREF((PyObject *)__pyx_v_outarr);
  if (__Pyx_PyList_SET_ITEM(__pyx_t_3, 1, ((PyObject *)__pyx_v_outarr)) != (0)) __PYX_ERR(0, 184, __pyx_L1_error);
  __pyx_t_7 = PyList_New(2); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 184, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_7);
  __Pyx_INCREF(__pyx_mstate_global->__pyx_n_u_external_loop);
  __Pyx_GIVEREF(__pyx_mstate_global->__pyx_n_u_external_loop);
  if (__Pyx_PyList_SET_ITEM(__pyx_t_7, 0, __pyx_mstate_global->__pyx_n_u_external_loop) != (0)) __PYX_ERR(0, 184, __pyx_L1_error);
  __Pyx_INCREF(__pyx_mstate_global->__pyx_n_u_buffered);
  __Pyx_GIVEREF(__pyx_mstate_global->__pyx_n_u_buffered);
  if (__Pyx_PyList_SET_ITEM(__pyx_t_7, 1, __pyx_mstate_global->__pyx_n_u_buffered) != (0)) __PYX_ERR(0, 184, __pyx_L1_error);

+185:                    op_flags=[['readonly'], ['writeonly']])

  __pyx_t_8 = PyList_New(1); if (unlikely(!__pyx_t_8)) __PYX_ERR(0, 185, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_8);
  __Pyx_INCREF(__pyx_mstate_global->__pyx_n_u_readonly);
  __Pyx_GIVEREF(__pyx_mstate_global->__pyx_n_u_readonly);
  if (__Pyx_PyList_SET_ITEM(__pyx_t_8, 0, __pyx_mstate_global->__pyx_n_u_readonly) != (0)) __PYX_ERR(0, 185, __pyx_L1_error);
  __pyx_t_9 = PyList_New(1); if (unlikely(!__pyx_t_9)) __PYX_ERR(0, 185, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_9);
  __Pyx_INCREF(__pyx_mstate_global->__pyx_n_u_writeonly);
  __Pyx_GIVEREF(__pyx_mstate_global->__pyx_n_u_writeonly);
  if (__Pyx_PyList_SET_ITEM(__pyx_t_9, 0, __pyx_mstate_global->__pyx_n_u_writeonly) != (0)) __PYX_ERR(0, 185, __pyx_L1_error);
  __pyx_t_10 = PyList_New(2); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 185, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_10);
  __Pyx_GIVEREF(__pyx_t_8);
  if (__Pyx_PyList_SET_ITEM(__pyx_t_10, 0, __pyx_t_8) != (0)) __PYX_ERR(0, 185, __pyx_L1_error);
  __Pyx_GIVEREF(__pyx_t_9);
  if (__Pyx_PyList_SET_ITEM(__pyx_t_10, 1, __pyx_t_9) != (0)) __PYX_ERR(0, 185, __pyx_L1_error);
  __pyx_t_8 = 0;
  __pyx_t_9 = 0;
  __pyx_t_6 = 1;
  #if CYTHON_UNPACK_METHODS
  if (unlikely(PyMethod_Check(__pyx_t_1))) {
    __pyx_t_5 = PyMethod_GET_SELF(__pyx_t_1);
    assert(__pyx_t_5);
    PyObject* __pyx__function = PyMethod_GET_FUNCTION(__pyx_t_1);
    __Pyx_INCREF(__pyx_t_5);
    __Pyx_INCREF(__pyx__function);
    __Pyx_DECREF_SET(__pyx_t_1, __pyx__function);
    __pyx_t_6 = 0;
  }
  #endif
  {
    PyObject *__pyx_callargs[2 + ((CYTHON_VECTORCALL) ? 2 : 0)] = {__pyx_t_5, __pyx_t_3};
    __pyx_t_9 = __Pyx_MakeVectorcallBuilderKwds(2); if (unlikely(!__pyx_t_9)) __PYX_ERR(0, 184, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_9);
    if (__Pyx_VectorcallBuilder_AddArg(__pyx_mstate_global->__pyx_n_u_flags, __pyx_t_7, __pyx_t_9, __pyx_callargs+2, 0) < (0)) __PYX_ERR(0, 184, __pyx_L1_error)
    if (__Pyx_VectorcallBuilder_AddArg(__pyx_mstate_global->__pyx_n_u_op_flags, __pyx_t_10, __pyx_t_9, __pyx_callargs+2, 1) < (0)) __PYX_ERR(0, 184, __pyx_L1_error)
    __pyx_t_2 = __Pyx_Object_Vectorcall_CallFromBuilder((PyObject*)__pyx_t_1, __pyx_callargs+__pyx_t_6, (2-__pyx_t_6) | (__pyx_t_6*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET), __pyx_t_9);
    __Pyx_XDECREF(__pyx_t_5); __pyx_t_5 = 0;
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
    __Pyx_DECREF(__pyx_t_10); __pyx_t_10 = 0;
    __Pyx_DECREF(__pyx_t_9); __pyx_t_9 = 0;
    __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
    if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 184, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_2);
  }
  __pyx_v_it = __pyx_t_2;
  __pyx_t_2 = 0;

+186:     for subinarr, suboutarr in it:

  if (likely(PyList_CheckExact(__pyx_v_it)) || PyTuple_CheckExact(__pyx_v_it)) {
    __pyx_t_2 = __pyx_v_it; __Pyx_INCREF(__pyx_t_2);
    __pyx_t_11 = 0;
    __pyx_t_12 = NULL;
  } else {
    __pyx_t_11 = -1; __pyx_t_2 = PyObject_GetIter(__pyx_v_it); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 186, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_2);
    __pyx_t_12 = (CYTHON_COMPILING_IN_LIMITED_API) ? PyIter_Next : __Pyx_PyObject_GetIterNextFunc(__pyx_t_2); if (unlikely(!__pyx_t_12)) __PYX_ERR(0, 186, __pyx_L1_error)
  }
  for (;;) {
    if (likely(!__pyx_t_12)) {
      if (likely(PyList_CheckExact(__pyx_t_2))) {
        {
          Py_ssize_t __pyx_temp = __Pyx_PyList_GET_SIZE(__pyx_t_2);
          #if !CYTHON_ASSUME_SAFE_SIZE
          if (unlikely((__pyx_temp < 0))) __PYX_ERR(0, 186, __pyx_L1_error)
          #endif
          if (__pyx_t_11 >= __pyx_temp) break;
        }
        __pyx_t_1 = __Pyx_PyList_GetItemRefFast(__pyx_t_2, __pyx_t_11, __Pyx_ReferenceSharing_OwnStrongReference);
        ++__pyx_t_11;
      } else {
        {
          Py_ssize_t __pyx_temp = __Pyx_PyTuple_GET_SIZE(__pyx_t_2);
          #if !CYTHON_ASSUME_SAFE_SIZE
          if (unlikely((__pyx_temp < 0))) __PYX_ERR(0, 186, __pyx_L1_error)
          #endif
          if (__pyx_t_11 >= __pyx_temp) break;
        }
        #if CYTHON_ASSUME_SAFE_MACROS && !CYTHON_AVOID_BORROWED_REFS
        __pyx_t_1 = __Pyx_NewRef(PyTuple_GET_ITEM(__pyx_t_2, __pyx_t_11));
        #else
        __pyx_t_1 = __Pyx_PySequence_ITEM(__pyx_t_2, __pyx_t_11);
        #endif
        ++__pyx_t_11;
      }
      if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 186, __pyx_L1_error)
    } else {
      __pyx_t_1 = __pyx_t_12(__pyx_t_2);
      if (unlikely(!__pyx_t_1)) {
        PyObject* exc_type = PyErr_Occurred();
        if (exc_type) {
          if (unlikely(!__Pyx_PyErr_GivenExceptionMatches(exc_type, PyExc_StopIteration))) __PYX_ERR(0, 186, __pyx_L1_error)
          PyErr_Clear();
        }
        break;
      }
    }
    __Pyx_GOTREF(__pyx_t_1);
    if ((likely(PyTuple_CheckExact(__pyx_t_1))) || (PyList_CheckExact(__pyx_t_1))) {
      PyObject* sequence = __pyx_t_1;
      Py_ssize_t size = __Pyx_PySequence_SIZE(sequence);
      if (unlikely(size != 2)) {
        if (size > 2) __Pyx_RaiseTooManyValuesError(2);
        else if (size >= 0) __Pyx_RaiseNeedMoreValuesError(size);
        __PYX_ERR(0, 186, __pyx_L1_error)
      }
      #if CYTHON_ASSUME_SAFE_MACROS && !CYTHON_AVOID_BORROWED_REFS
      if (likely(PyTuple_CheckExact(sequence))) {
        __pyx_t_9 = PyTuple_GET_ITEM(sequence, 0);
        __Pyx_INCREF(__pyx_t_9);
        __pyx_t_10 = PyTuple_GET_ITEM(sequence, 1);
        __Pyx_INCREF(__pyx_t_10);
      } else {
        __pyx_t_9 = __Pyx_PyList_GetItemRefFast(sequence, 0, __Pyx_ReferenceSharing_SharedReference);
        if (unlikely(!__pyx_t_9)) __PYX_ERR(0, 186, __pyx_L1_error)
        __Pyx_XGOTREF(__pyx_t_9);
        __pyx_t_10 = __Pyx_PyList_GetItemRefFast(sequence, 1, __Pyx_ReferenceSharing_SharedReference);
        if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 186, __pyx_L1_error)
        __Pyx_XGOTREF(__pyx_t_10);
      }
      #else
      __pyx_t_9 = __Pyx_PySequence_ITEM(sequence, 0); if (unlikely(!__pyx_t_9)) __PYX_ERR(0, 186, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_9);
      __pyx_t_10 = __Pyx_PySequence_ITEM(sequence, 1); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 186, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_10);
      #endif
      __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
    } else {
      Py_ssize_t index = -1;
      __pyx_t_7 = PyObject_GetIter(__pyx_t_1); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 186, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_7);
      __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
      __pyx_t_13 = (CYTHON_COMPILING_IN_LIMITED_API) ? PyIter_Next : __Pyx_PyObject_GetIterNextFunc(__pyx_t_7);
      index = 0; __pyx_t_9 = __pyx_t_13(__pyx_t_7); if (unlikely(!__pyx_t_9)) goto __pyx_L6_unpacking_failed;
      __Pyx_GOTREF(__pyx_t_9);
      index = 1; __pyx_t_10 = __pyx_t_13(__pyx_t_7); if (unlikely(!__pyx_t_10)) goto __pyx_L6_unpacking_failed;
      __Pyx_GOTREF(__pyx_t_10);
      if (__Pyx_IternextUnpackEndCheck(__pyx_t_13(__pyx_t_7), 2) < (0)) __PYX_ERR(0, 186, __pyx_L1_error)
      __pyx_t_13 = NULL;
      __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
      goto __pyx_L7_unpacking_done;
      __pyx_L6_unpacking_failed:;
      __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
      __pyx_t_13 = NULL;
      if (__Pyx_IterFinish() == 0) __Pyx_RaiseNeedMoreValuesError(index);
      __PYX_ERR(0, 186, __pyx_L1_error)
      __pyx_L7_unpacking_done:;
    }
    if (!(likely(((__pyx_t_9) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_9, __pyx_mstate_global->__pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 186, __pyx_L1_error)
    if (!(likely(((__pyx_t_10) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_10, __pyx_mstate_global->__pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 186, __pyx_L1_error)
    {
      __Pyx_BufFmt_StackElem __pyx_stack[1];
      __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_subinarr.rcbuffer->pybuffer);
      __pyx_t_14 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_subinarr.rcbuffer->pybuffer, (PyObject*)((PyArrayObject *)__pyx_t_9), &__Pyx_TypeInfo_double, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack);
      if (unlikely(__pyx_t_14 < 0)) {
        PyErr_Fetch(&__pyx_t_15, &__pyx_t_16, &__pyx_t_17);
        if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_subinarr.rcbuffer->pybuffer, (PyObject*)__pyx_v_subinarr, &__Pyx_TypeInfo_double, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {
          Py_XDECREF(__pyx_t_15); Py_XDECREF(__pyx_t_16); Py_XDECREF(__pyx_t_17);
          __Pyx_RaiseBufferFallbackError();
        } else {
          PyErr_Restore(__pyx_t_15, __pyx_t_16, __pyx_t_17);
        }
        __pyx_t_15 = __pyx_t_16 = __pyx_t_17 = 0;
      }
      __pyx_pybuffernd_subinarr.diminfo[0].strides = __pyx_pybuffernd_subinarr.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_subinarr.diminfo[0].shape = __pyx_pybuffernd_subinarr.rcbuffer->pybuffer.shape[0];
      if (unlikely((__pyx_t_14 < 0))) __PYX_ERR(0, 186, __pyx_L1_error)
    }
    __Pyx_XDECREF_SET(__pyx_v_subinarr, ((PyArrayObject *)__pyx_t_9));
    __pyx_t_9 = 0;
    {
      __Pyx_BufFmt_StackElem __pyx_stack[1];
      __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_suboutarr.rcbuffer->pybuffer);
      __pyx_t_14 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_suboutarr.rcbuffer->pybuffer, (PyObject*)((PyArrayObject *)__pyx_t_10), &__Pyx_TypeInfo_double, PyBUF_FORMAT| PyBUF_STRIDES| PyBUF_WRITABLE, 1, 0, __pyx_stack);
      if (unlikely(__pyx_t_14 < 0)) {
        PyErr_Fetch(&__pyx_t_17, &__pyx_t_16, &__pyx_t_15);
        if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_suboutarr.rcbuffer->pybuffer, (PyObject*)__pyx_v_suboutarr, &__Pyx_TypeInfo_double, PyBUF_FORMAT| PyBUF_STRIDES| PyBUF_WRITABLE, 1, 0, __pyx_stack) == -1)) {
          Py_XDECREF(__pyx_t_17); Py_XDECREF(__pyx_t_16); Py_XDECREF(__pyx_t_15);
          __Pyx_RaiseBufferFallbackError();
        } else {
          PyErr_Restore(__pyx_t_17, __pyx_t_16, __pyx_t_15);
        }
        __pyx_t_17 = __pyx_t_16 = __pyx_t_15 = 0;
      }
      __pyx_pybuffernd_suboutarr.diminfo[0].strides = __pyx_pybuffernd_suboutarr.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_suboutarr.diminfo[0].shape = __pyx_pybuffernd_suboutarr.rcbuffer->pybuffer.shape[0];
      if (unlikely((__pyx_t_14 < 0))) __PYX_ERR(0, 186, __pyx_L1_error)
    }
    __Pyx_XDECREF_SET(__pyx_v_suboutarr, ((PyArrayObject *)__pyx_t_10));
    __pyx_t_10 = 0;
/* … */
  }
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;

+187:         for j in range(subinarr.shape[0]):

    __pyx_t_18 = (__pyx_f_5numpy_7ndarray_5shape_shape(((PyArrayObject *)__pyx_v_subinarr))[0]);
    __pyx_t_19 = __pyx_t_18;
    for (__pyx_t_14 = 0; __pyx_t_14 < __pyx_t_19; __pyx_t_14+=1) {
      __pyx_v_j = __pyx_t_14;

+188:             suboutarr[j]=func(subinarr[j])

      __pyx_t_20 = __pyx_v_j;
      if (__pyx_t_20 < 0) __pyx_t_20 += __pyx_pybuffernd_subinarr.diminfo[0].shape;
      __pyx_t_21 = __pyx_v_func((*__Pyx_BufPtrStrided1d(double *, __pyx_pybuffernd_subinarr.rcbuffer->pybuffer.buf, __pyx_t_20, __pyx_pybuffernd_subinarr.diminfo[0].strides))); if (unlikely(PyErr_Occurred())) __PYX_ERR(0, 188, __pyx_L1_error)
      __pyx_t_20 = __pyx_v_j;
      if (__pyx_t_20 < 0) __pyx_t_20 += __pyx_pybuffernd_suboutarr.diminfo[0].shape;
      *__Pyx_BufPtrStrided1d(double *, __pyx_pybuffernd_suboutarr.rcbuffer->pybuffer.buf, __pyx_t_20, __pyx_pybuffernd_suboutarr.diminfo[0].strides) = __pyx_t_21;
    }

+189:     return outarr

  __Pyx_XDECREF((PyObject *)__pyx_r);
  __Pyx_INCREF((PyObject *)__pyx_v_outarr);
  __pyx_r = __pyx_v_outarr;
  __Pyx_TraceReturnValue((PyObject *)__pyx_r, 54, 0, __PYX_ERR(0, 189, __pyx_L1_error));
  goto __pyx_L0;

 190: 

+191: @cython.boundscheck(False)

static PyArrayObject *__pyx_f_9pynitride_4core_12cython_maths_map2(PyArrayObject *__pyx_v_inarr1, PyArrayObject *__pyx_v_inarr2, double (*__pyx_v_func)(double, double), struct __pyx_opt_args_9pynitride_4core_12cython_maths_map2 *__pyx_optional_args) {
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_1);
  __Pyx_XDECREF(__pyx_t_2);
  __Pyx_XDECREF(__pyx_t_3);
  __Pyx_XDECREF(__pyx_t_5);
  __Pyx_XDECREF(__pyx_t_7);
  __Pyx_XDECREF(__pyx_t_8);
  __Pyx_XDECREF(__pyx_t_9);
  __Pyx_XDECREF(__pyx_t_10);
  __Pyx_XDECREF(__pyx_t_11);
  { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
    __Pyx_PyThreadState_declare
    __Pyx_PyThreadState_assign
    __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_subinarr1.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_subinarr2.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_suboutarr.rcbuffer->pybuffer);
  __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
  __Pyx_TraceException(__pyx_lineno, 0, 0);
  #if CYTHON_USE_SYS_MONITORING
  __Pyx_TraceExceptionUnwind(0, 0);
  #else
  __Pyx_TraceReturnValue(NULL, 0, 0, __PYX_ERR(0, 191, __pyx_L1_error));
  #endif
  __Pyx_AddTraceback("pynitride.core.cython_maths.map2", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = 0;
  goto __pyx_L2;
  __pyx_L0:;
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_subinarr1.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_subinarr2.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_suboutarr.rcbuffer->pybuffer);
  __pyx_L2:;
  __Pyx_XDECREF((PyObject *)__pyx_v_subinarr1);
  __Pyx_XDECREF((PyObject *)__pyx_v_subinarr2);
  __Pyx_XDECREF((PyObject *)__pyx_v_suboutarr);
  __Pyx_XDECREF(__pyx_v_it);
  __Pyx_XDECREF((PyObject *)__pyx_v_outarr);
  __Pyx_XGIVEREF((PyObject *)__pyx_r);
  __Pyx_PyMonitoring_ExitScope(0);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

+192: cdef cnp.ndarray map2(  cnp.ndarray inarr1,

struct __pyx_opt_args_9pynitride_4core_12cython_maths_map2 {
  int __pyx_n;
  PyArrayObject *outarr;
};

 193:                         cnp.ndarray inarr2,

 194:                         double (*func)(double,double),

+195:                         cnp.ndarray outarr=None):

  PyArrayObject *__pyx_v_outarr = ((PyArrayObject *)Py_None);
  int __pyx_v_j;
  PyArrayObject *__pyx_v_subinarr1 = 0;
  PyArrayObject *__pyx_v_subinarr2 = 0;
  PyArrayObject *__pyx_v_suboutarr = 0;
  PyObject *__pyx_v_it = NULL;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_subinarr1;
  __Pyx_Buffer __pyx_pybuffer_subinarr1;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_subinarr2;
  __Pyx_Buffer __pyx_pybuffer_subinarr2;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_suboutarr;
  __Pyx_Buffer __pyx_pybuffer_suboutarr;
  PyArrayObject *__pyx_r = NULL;
  __Pyx_TraceDeclarationsFunc
  __Pyx_TraceFrameInit(((PyObject *)__pyx_mstate_global->__pyx_codeobj_tab[44]))
  __Pyx_TraceStartFunc("map2", __pyx_f[0], 191, 0, 0, 0, __PYX_ERR(0, 191, __pyx_L1_error));
  if (__pyx_optional_args) {
    if (__pyx_optional_args->__pyx_n > 0) {
      __pyx_v_outarr = __pyx_optional_args->outarr;
    }
  }
  __Pyx_INCREF((PyObject *)__pyx_v_outarr);
  __pyx_pybuffer_subinarr1.pybuffer.buf = NULL;
  __pyx_pybuffer_subinarr1.refcount = 0;
  __pyx_pybuffernd_subinarr1.data = NULL;
  __pyx_pybuffernd_subinarr1.rcbuffer = &__pyx_pybuffer_subinarr1;
  __pyx_pybuffer_subinarr2.pybuffer.buf = NULL;
  __pyx_pybuffer_subinarr2.refcount = 0;
  __pyx_pybuffernd_subinarr2.data = NULL;
  __pyx_pybuffernd_subinarr2.rcbuffer = &__pyx_pybuffer_subinarr2;
  __pyx_pybuffer_suboutarr.pybuffer.buf = NULL;
  __pyx_pybuffer_suboutarr.refcount = 0;
  __pyx_pybuffernd_suboutarr.data = NULL;
  __pyx_pybuffernd_suboutarr.rcbuffer = &__pyx_pybuffer_suboutarr;

 196:     r""" Map the function `func` over `inarr1` and `inarr2` to fill `outarr`.

 197: 

 198:     Args:

 199:         inarr1: double array of inputs

 200:         inarr2: double array of inputs

 201:         func: a function to apply to each element of the inputs

 202:             takes a two `double`s as arguments

 203:         outarr: optional double array of outputs to fill

 204:     Returns:

 205:         the double array of outputs

 206:     """

 207:     cdef int j

 208:     cdef cnp.ndarray[double] subinarr1, subinarr2, suboutarr

 209: 

+210:     assert inarr1.dtype==np.float64

  #ifndef CYTHON_WITHOUT_ASSERTIONS
  if (unlikely(__pyx_assertions_enabled())) {
    __pyx_t_1 = __Pyx_PyObject_GetAttrStr(((PyObject *)__pyx_v_inarr1), __pyx_mstate_global->__pyx_n_u_dtype); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 210, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 210, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_2);
    __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_mstate_global->__pyx_n_u_float64); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 210, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_3);
    __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
    __pyx_t_2 = PyObject_RichCompare(__pyx_t_1, __pyx_t_3, Py_EQ); __Pyx_XGOTREF(__pyx_t_2); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 210, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    __pyx_t_4 = __Pyx_PyObject_IsTrue(__pyx_t_2); if (unlikely((__pyx_t_4 < 0))) __PYX_ERR(0, 210, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
    if (unlikely(!__pyx_t_4)) {
      __Pyx_Raise(((PyObject *)(((PyTypeObject*)PyExc_AssertionError))), 0, 0, 0);
      __PYX_ERR(0, 210, __pyx_L1_error)
    }
  }
  #else
  if ((1)); else __PYX_ERR(0, 210, __pyx_L1_error)
  #endif

+211:     assert inarr2.dtype==np.float64

  #ifndef CYTHON_WITHOUT_ASSERTIONS
  if (unlikely(__pyx_assertions_enabled())) {
    __pyx_t_2 = __Pyx_PyObject_GetAttrStr(((PyObject *)__pyx_v_inarr2), __pyx_mstate_global->__pyx_n_u_dtype); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 211, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_2);
    __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 211, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_3);
    __pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_float64); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 211, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    __pyx_t_3 = PyObject_RichCompare(__pyx_t_2, __pyx_t_1, Py_EQ); __Pyx_XGOTREF(__pyx_t_3); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 211, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
    __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
    __pyx_t_4 = __Pyx_PyObject_IsTrue(__pyx_t_3); if (unlikely((__pyx_t_4 < 0))) __PYX_ERR(0, 211, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    if (unlikely(!__pyx_t_4)) {
      __Pyx_Raise(((PyObject *)(((PyTypeObject*)PyExc_AssertionError))), 0, 0, 0);
      __PYX_ERR(0, 211, __pyx_L1_error)
    }
  }
  #else
  if ((1)); else __PYX_ERR(0, 211, __pyx_L1_error)
  #endif

+212:     if outarr is None:

  __pyx_t_4 = (((PyObject *)__pyx_v_outarr) == Py_None);
  if (__pyx_t_4) {
/* … */
  }

+213:         outarr=np.empty_like(inarr1)

    __pyx_t_1 = NULL;
    __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 213, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_2);
    __pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_mstate_global->__pyx_n_u_empty_like); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 213, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_5);
    __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
    __pyx_t_6 = 1;
    #if CYTHON_UNPACK_METHODS
    if (unlikely(PyMethod_Check(__pyx_t_5))) {
      __pyx_t_1 = PyMethod_GET_SELF(__pyx_t_5);
      assert(__pyx_t_1);
      PyObject* __pyx__function = PyMethod_GET_FUNCTION(__pyx_t_5);
      __Pyx_INCREF(__pyx_t_1);
      __Pyx_INCREF(__pyx__function);
      __Pyx_DECREF_SET(__pyx_t_5, __pyx__function);
      __pyx_t_6 = 0;
    }
    #endif
    {
      PyObject *__pyx_callargs[2] = {__pyx_t_1, ((PyObject *)__pyx_v_inarr1)};
      __pyx_t_3 = __Pyx_PyObject_FastCall((PyObject*)__pyx_t_5, __pyx_callargs+__pyx_t_6, (2-__pyx_t_6) | (__pyx_t_6*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET));
      __Pyx_XDECREF(__pyx_t_1); __pyx_t_1 = 0;
      __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
      if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 213, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_3);
    }
    if (!(likely(((__pyx_t_3) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_3, __pyx_mstate_global->__pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 213, __pyx_L1_error)
    __Pyx_DECREF_SET(__pyx_v_outarr, ((PyArrayObject *)__pyx_t_3));
    __pyx_t_3 = 0;

 214: 

+215:     it = np.nditer([inarr1,inarr2,outarr], flags=['external_loop','buffered'],

  __pyx_t_5 = NULL;
  __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 215, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_mstate_global->__pyx_n_u_nditer); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 215, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_1 = PyList_New(3); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 215, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_INCREF((PyObject *)__pyx_v_inarr1);
  __Pyx_GIVEREF((PyObject *)__pyx_v_inarr1);
  if (__Pyx_PyList_SET_ITEM(__pyx_t_1, 0, ((PyObject *)__pyx_v_inarr1)) != (0)) __PYX_ERR(0, 215, __pyx_L1_error);
  __Pyx_INCREF((PyObject *)__pyx_v_inarr2);
  __Pyx_GIVEREF((PyObject *)__pyx_v_inarr2);
  if (__Pyx_PyList_SET_ITEM(__pyx_t_1, 1, ((PyObject *)__pyx_v_inarr2)) != (0)) __PYX_ERR(0, 215, __pyx_L1_error);
  __Pyx_INCREF((PyObject *)__pyx_v_outarr);
  __Pyx_GIVEREF((PyObject *)__pyx_v_outarr);
  if (__Pyx_PyList_SET_ITEM(__pyx_t_1, 2, ((PyObject *)__pyx_v_outarr)) != (0)) __PYX_ERR(0, 215, __pyx_L1_error);
  __pyx_t_7 = PyList_New(2); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 215, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_7);
  __Pyx_INCREF(__pyx_mstate_global->__pyx_n_u_external_loop);
  __Pyx_GIVEREF(__pyx_mstate_global->__pyx_n_u_external_loop);
  if (__Pyx_PyList_SET_ITEM(__pyx_t_7, 0, __pyx_mstate_global->__pyx_n_u_external_loop) != (0)) __PYX_ERR(0, 215, __pyx_L1_error);
  __Pyx_INCREF(__pyx_mstate_global->__pyx_n_u_buffered);
  __Pyx_GIVEREF(__pyx_mstate_global->__pyx_n_u_buffered);
  if (__Pyx_PyList_SET_ITEM(__pyx_t_7, 1, __pyx_mstate_global->__pyx_n_u_buffered) != (0)) __PYX_ERR(0, 215, __pyx_L1_error);

+216:                    op_flags=[['readonly'], ['readonly'], ['writeonly']])

  __pyx_t_8 = PyList_New(1); if (unlikely(!__pyx_t_8)) __PYX_ERR(0, 216, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_8);
  __Pyx_INCREF(__pyx_mstate_global->__pyx_n_u_readonly);
  __Pyx_GIVEREF(__pyx_mstate_global->__pyx_n_u_readonly);
  if (__Pyx_PyList_SET_ITEM(__pyx_t_8, 0, __pyx_mstate_global->__pyx_n_u_readonly) != (0)) __PYX_ERR(0, 216, __pyx_L1_error);
  __pyx_t_9 = PyList_New(1); if (unlikely(!__pyx_t_9)) __PYX_ERR(0, 216, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_9);
  __Pyx_INCREF(__pyx_mstate_global->__pyx_n_u_readonly);
  __Pyx_GIVEREF(__pyx_mstate_global->__pyx_n_u_readonly);
  if (__Pyx_PyList_SET_ITEM(__pyx_t_9, 0, __pyx_mstate_global->__pyx_n_u_readonly) != (0)) __PYX_ERR(0, 216, __pyx_L1_error);
  __pyx_t_10 = PyList_New(1); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 216, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_10);
  __Pyx_INCREF(__pyx_mstate_global->__pyx_n_u_writeonly);
  __Pyx_GIVEREF(__pyx_mstate_global->__pyx_n_u_writeonly);
  if (__Pyx_PyList_SET_ITEM(__pyx_t_10, 0, __pyx_mstate_global->__pyx_n_u_writeonly) != (0)) __PYX_ERR(0, 216, __pyx_L1_error);
  __pyx_t_11 = PyList_New(3); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 216, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_11);
  __Pyx_GIVEREF(__pyx_t_8);
  if (__Pyx_PyList_SET_ITEM(__pyx_t_11, 0, __pyx_t_8) != (0)) __PYX_ERR(0, 216, __pyx_L1_error);
  __Pyx_GIVEREF(__pyx_t_9);
  if (__Pyx_PyList_SET_ITEM(__pyx_t_11, 1, __pyx_t_9) != (0)) __PYX_ERR(0, 216, __pyx_L1_error);
  __Pyx_GIVEREF(__pyx_t_10);
  if (__Pyx_PyList_SET_ITEM(__pyx_t_11, 2, __pyx_t_10) != (0)) __PYX_ERR(0, 216, __pyx_L1_error);
  __pyx_t_8 = 0;
  __pyx_t_9 = 0;
  __pyx_t_10 = 0;
  __pyx_t_6 = 1;
  #if CYTHON_UNPACK_METHODS
  if (unlikely(PyMethod_Check(__pyx_t_2))) {
    __pyx_t_5 = PyMethod_GET_SELF(__pyx_t_2);
    assert(__pyx_t_5);
    PyObject* __pyx__function = PyMethod_GET_FUNCTION(__pyx_t_2);
    __Pyx_INCREF(__pyx_t_5);
    __Pyx_INCREF(__pyx__function);
    __Pyx_DECREF_SET(__pyx_t_2, __pyx__function);
    __pyx_t_6 = 0;
  }
  #endif
  {
    PyObject *__pyx_callargs[2 + ((CYTHON_VECTORCALL) ? 2 : 0)] = {__pyx_t_5, __pyx_t_1};
    __pyx_t_10 = __Pyx_MakeVectorcallBuilderKwds(2); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 215, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_10);
    if (__Pyx_VectorcallBuilder_AddArg(__pyx_mstate_global->__pyx_n_u_flags, __pyx_t_7, __pyx_t_10, __pyx_callargs+2, 0) < (0)) __PYX_ERR(0, 215, __pyx_L1_error)
    if (__Pyx_VectorcallBuilder_AddArg(__pyx_mstate_global->__pyx_n_u_op_flags, __pyx_t_11, __pyx_t_10, __pyx_callargs+2, 1) < (0)) __PYX_ERR(0, 215, __pyx_L1_error)
    __pyx_t_3 = __Pyx_Object_Vectorcall_CallFromBuilder((PyObject*)__pyx_t_2, __pyx_callargs+__pyx_t_6, (2-__pyx_t_6) | (__pyx_t_6*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET), __pyx_t_10);
    __Pyx_XDECREF(__pyx_t_5); __pyx_t_5 = 0;
    __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
    __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
    __Pyx_DECREF(__pyx_t_11); __pyx_t_11 = 0;
    __Pyx_DECREF(__pyx_t_10); __pyx_t_10 = 0;
    __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
    if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 215, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_3);
  }
  __pyx_v_it = __pyx_t_3;
  __pyx_t_3 = 0;

+217:     for subinarr1, subinarr2, suboutarr in it:

  if (likely(PyList_CheckExact(__pyx_v_it)) || PyTuple_CheckExact(__pyx_v_it)) {
    __pyx_t_3 = __pyx_v_it; __Pyx_INCREF(__pyx_t_3);
    __pyx_t_12 = 0;
    __pyx_t_13 = NULL;
  } else {
    __pyx_t_12 = -1; __pyx_t_3 = PyObject_GetIter(__pyx_v_it); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 217, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_3);
    __pyx_t_13 = (CYTHON_COMPILING_IN_LIMITED_API) ? PyIter_Next : __Pyx_PyObject_GetIterNextFunc(__pyx_t_3); if (unlikely(!__pyx_t_13)) __PYX_ERR(0, 217, __pyx_L1_error)
  }
  for (;;) {
    if (likely(!__pyx_t_13)) {
      if (likely(PyList_CheckExact(__pyx_t_3))) {
        {
          Py_ssize_t __pyx_temp = __Pyx_PyList_GET_SIZE(__pyx_t_3);
          #if !CYTHON_ASSUME_SAFE_SIZE
          if (unlikely((__pyx_temp < 0))) __PYX_ERR(0, 217, __pyx_L1_error)
          #endif
          if (__pyx_t_12 >= __pyx_temp) break;
        }
        __pyx_t_2 = __Pyx_PyList_GetItemRefFast(__pyx_t_3, __pyx_t_12, __Pyx_ReferenceSharing_OwnStrongReference);
        ++__pyx_t_12;
      } else {
        {
          Py_ssize_t __pyx_temp = __Pyx_PyTuple_GET_SIZE(__pyx_t_3);
          #if !CYTHON_ASSUME_SAFE_SIZE
          if (unlikely((__pyx_temp < 0))) __PYX_ERR(0, 217, __pyx_L1_error)
          #endif
          if (__pyx_t_12 >= __pyx_temp) break;
        }
        #if CYTHON_ASSUME_SAFE_MACROS && !CYTHON_AVOID_BORROWED_REFS
        __pyx_t_2 = __Pyx_NewRef(PyTuple_GET_ITEM(__pyx_t_3, __pyx_t_12));
        #else
        __pyx_t_2 = __Pyx_PySequence_ITEM(__pyx_t_3, __pyx_t_12);
        #endif
        ++__pyx_t_12;
      }
      if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 217, __pyx_L1_error)
    } else {
      __pyx_t_2 = __pyx_t_13(__pyx_t_3);
      if (unlikely(!__pyx_t_2)) {
        PyObject* exc_type = PyErr_Occurred();
        if (exc_type) {
          if (unlikely(!__Pyx_PyErr_GivenExceptionMatches(exc_type, PyExc_StopIteration))) __PYX_ERR(0, 217, __pyx_L1_error)
          PyErr_Clear();
        }
        break;
      }
    }
    __Pyx_GOTREF(__pyx_t_2);
    if ((likely(PyTuple_CheckExact(__pyx_t_2))) || (PyList_CheckExact(__pyx_t_2))) {
      PyObject* sequence = __pyx_t_2;
      Py_ssize_t size = __Pyx_PySequence_SIZE(sequence);
      if (unlikely(size != 3)) {
        if (size > 3) __Pyx_RaiseTooManyValuesError(3);
        else if (size >= 0) __Pyx_RaiseNeedMoreValuesError(size);
        __PYX_ERR(0, 217, __pyx_L1_error)
      }
      #if CYTHON_ASSUME_SAFE_MACROS && !CYTHON_AVOID_BORROWED_REFS
      if (likely(PyTuple_CheckExact(sequence))) {
        __pyx_t_10 = PyTuple_GET_ITEM(sequence, 0);
        __Pyx_INCREF(__pyx_t_10);
        __pyx_t_11 = PyTuple_GET_ITEM(sequence, 1);
        __Pyx_INCREF(__pyx_t_11);
        __pyx_t_7 = PyTuple_GET_ITEM(sequence, 2);
        __Pyx_INCREF(__pyx_t_7);
      } else {
        __pyx_t_10 = __Pyx_PyList_GetItemRefFast(sequence, 0, __Pyx_ReferenceSharing_SharedReference);
        if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 217, __pyx_L1_error)
        __Pyx_XGOTREF(__pyx_t_10);
        __pyx_t_11 = __Pyx_PyList_GetItemRefFast(sequence, 1, __Pyx_ReferenceSharing_SharedReference);
        if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 217, __pyx_L1_error)
        __Pyx_XGOTREF(__pyx_t_11);
        __pyx_t_7 = __Pyx_PyList_GetItemRefFast(sequence, 2, __Pyx_ReferenceSharing_SharedReference);
        if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 217, __pyx_L1_error)
        __Pyx_XGOTREF(__pyx_t_7);
      }
      #else
      __pyx_t_10 = __Pyx_PySequence_ITEM(sequence, 0); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 217, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_10);
      __pyx_t_11 = __Pyx_PySequence_ITEM(sequence, 1); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 217, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_11);
      __pyx_t_7 = __Pyx_PySequence_ITEM(sequence, 2); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 217, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_7);
      #endif
      __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
    } else {
      Py_ssize_t index = -1;
      __pyx_t_1 = PyObject_GetIter(__pyx_t_2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 217, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_1);
      __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
      __pyx_t_14 = (CYTHON_COMPILING_IN_LIMITED_API) ? PyIter_Next : __Pyx_PyObject_GetIterNextFunc(__pyx_t_1);
      index = 0; __pyx_t_10 = __pyx_t_14(__pyx_t_1); if (unlikely(!__pyx_t_10)) goto __pyx_L6_unpacking_failed;
      __Pyx_GOTREF(__pyx_t_10);
      index = 1; __pyx_t_11 = __pyx_t_14(__pyx_t_1); if (unlikely(!__pyx_t_11)) goto __pyx_L6_unpacking_failed;
      __Pyx_GOTREF(__pyx_t_11);
      index = 2; __pyx_t_7 = __pyx_t_14(__pyx_t_1); if (unlikely(!__pyx_t_7)) goto __pyx_L6_unpacking_failed;
      __Pyx_GOTREF(__pyx_t_7);
      if (__Pyx_IternextUnpackEndCheck(__pyx_t_14(__pyx_t_1), 3) < (0)) __PYX_ERR(0, 217, __pyx_L1_error)
      __pyx_t_14 = NULL;
      __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
      goto __pyx_L7_unpacking_done;
      __pyx_L6_unpacking_failed:;
      __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
      __pyx_t_14 = NULL;
      if (__Pyx_IterFinish() == 0) __Pyx_RaiseNeedMoreValuesError(index);
      __PYX_ERR(0, 217, __pyx_L1_error)
      __pyx_L7_unpacking_done:;
    }
    if (!(likely(((__pyx_t_10) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_10, __pyx_mstate_global->__pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 217, __pyx_L1_error)
    if (!(likely(((__pyx_t_11) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_11, __pyx_mstate_global->__pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 217, __pyx_L1_error)
    if (!(likely(((__pyx_t_7) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_7, __pyx_mstate_global->__pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 217, __pyx_L1_error)
    {
      __Pyx_BufFmt_StackElem __pyx_stack[1];
      __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_subinarr1.rcbuffer->pybuffer);
      __pyx_t_15 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_subinarr1.rcbuffer->pybuffer, (PyObject*)((PyArrayObject *)__pyx_t_10), &__Pyx_TypeInfo_double, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack);
      if (unlikely(__pyx_t_15 < 0)) {
        PyErr_Fetch(&__pyx_t_16, &__pyx_t_17, &__pyx_t_18);
        if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_subinarr1.rcbuffer->pybuffer, (PyObject*)__pyx_v_subinarr1, &__Pyx_TypeInfo_double, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {
          Py_XDECREF(__pyx_t_16); Py_XDECREF(__pyx_t_17); Py_XDECREF(__pyx_t_18);
          __Pyx_RaiseBufferFallbackError();
        } else {
          PyErr_Restore(__pyx_t_16, __pyx_t_17, __pyx_t_18);
        }
        __pyx_t_16 = __pyx_t_17 = __pyx_t_18 = 0;
      }
      __pyx_pybuffernd_subinarr1.diminfo[0].strides = __pyx_pybuffernd_subinarr1.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_subinarr1.diminfo[0].shape = __pyx_pybuffernd_subinarr1.rcbuffer->pybuffer.shape[0];
      if (unlikely((__pyx_t_15 < 0))) __PYX_ERR(0, 217, __pyx_L1_error)
    }
    __Pyx_XDECREF_SET(__pyx_v_subinarr1, ((PyArrayObject *)__pyx_t_10));
    __pyx_t_10 = 0;
    {
      __Pyx_BufFmt_StackElem __pyx_stack[1];
      __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_subinarr2.rcbuffer->pybuffer);
      __pyx_t_15 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_subinarr2.rcbuffer->pybuffer, (PyObject*)((PyArrayObject *)__pyx_t_11), &__Pyx_TypeInfo_double, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack);
      if (unlikely(__pyx_t_15 < 0)) {
        PyErr_Fetch(&__pyx_t_18, &__pyx_t_17, &__pyx_t_16);
        if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_subinarr2.rcbuffer->pybuffer, (PyObject*)__pyx_v_subinarr2, &__Pyx_TypeInfo_double, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {
          Py_XDECREF(__pyx_t_18); Py_XDECREF(__pyx_t_17); Py_XDECREF(__pyx_t_16);
          __Pyx_RaiseBufferFallbackError();
        } else {
          PyErr_Restore(__pyx_t_18, __pyx_t_17, __pyx_t_16);
        }
        __pyx_t_18 = __pyx_t_17 = __pyx_t_16 = 0;
      }
      __pyx_pybuffernd_subinarr2.diminfo[0].strides = __pyx_pybuffernd_subinarr2.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_subinarr2.diminfo[0].shape = __pyx_pybuffernd_subinarr2.rcbuffer->pybuffer.shape[0];
      if (unlikely((__pyx_t_15 < 0))) __PYX_ERR(0, 217, __pyx_L1_error)
    }
    __Pyx_XDECREF_SET(__pyx_v_subinarr2, ((PyArrayObject *)__pyx_t_11));
    __pyx_t_11 = 0;
    {
      __Pyx_BufFmt_StackElem __pyx_stack[1];
      __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_suboutarr.rcbuffer->pybuffer);
      __pyx_t_15 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_suboutarr.rcbuffer->pybuffer, (PyObject*)((PyArrayObject *)__pyx_t_7), &__Pyx_TypeInfo_double, PyBUF_FORMAT| PyBUF_STRIDES| PyBUF_WRITABLE, 1, 0, __pyx_stack);
      if (unlikely(__pyx_t_15 < 0)) {
        PyErr_Fetch(&__pyx_t_16, &__pyx_t_17, &__pyx_t_18);
        if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_suboutarr.rcbuffer->pybuffer, (PyObject*)__pyx_v_suboutarr, &__Pyx_TypeInfo_double, PyBUF_FORMAT| PyBUF_STRIDES| PyBUF_WRITABLE, 1, 0, __pyx_stack) == -1)) {
          Py_XDECREF(__pyx_t_16); Py_XDECREF(__pyx_t_17); Py_XDECREF(__pyx_t_18);
          __Pyx_RaiseBufferFallbackError();
        } else {
          PyErr_Restore(__pyx_t_16, __pyx_t_17, __pyx_t_18);
        }
        __pyx_t_16 = __pyx_t_17 = __pyx_t_18 = 0;
      }
      __pyx_pybuffernd_suboutarr.diminfo[0].strides = __pyx_pybuffernd_suboutarr.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_suboutarr.diminfo[0].shape = __pyx_pybuffernd_suboutarr.rcbuffer->pybuffer.shape[0];
      if (unlikely((__pyx_t_15 < 0))) __PYX_ERR(0, 217, __pyx_L1_error)
    }
    __Pyx_XDECREF_SET(__pyx_v_suboutarr, ((PyArrayObject *)__pyx_t_7));
    __pyx_t_7 = 0;
/* … */
  }
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;

+218:         for j in range(subinarr1.shape[0]):

    __pyx_t_19 = (__pyx_f_5numpy_7ndarray_5shape_shape(((PyArrayObject *)__pyx_v_subinarr1))[0]);
    __pyx_t_20 = __pyx_t_19;
    for (__pyx_t_15 = 0; __pyx_t_15 < __pyx_t_20; __pyx_t_15+=1) {
      __pyx_v_j = __pyx_t_15;

+219:             suboutarr[j]=func(subinarr1[j],subinarr2[j])

      __pyx_t_21 = __pyx_v_j;
      if (__pyx_t_21 < 0) __pyx_t_21 += __pyx_pybuffernd_subinarr1.diminfo[0].shape;
      __pyx_t_22 = __pyx_v_j;
      if (__pyx_t_22 < 0) __pyx_t_22 += __pyx_pybuffernd_subinarr2.diminfo[0].shape;
      __pyx_t_23 = __pyx_v_func((*__Pyx_BufPtrStrided1d(double *, __pyx_pybuffernd_subinarr1.rcbuffer->pybuffer.buf, __pyx_t_21, __pyx_pybuffernd_subinarr1.diminfo[0].strides)), (*__Pyx_BufPtrStrided1d(double *, __pyx_pybuffernd_subinarr2.rcbuffer->pybuffer.buf, __pyx_t_22, __pyx_pybuffernd_subinarr2.diminfo[0].strides))); if (unlikely(PyErr_Occurred())) __PYX_ERR(0, 219, __pyx_L1_error)
      __pyx_t_22 = __pyx_v_j;
      if (__pyx_t_22 < 0) __pyx_t_22 += __pyx_pybuffernd_suboutarr.diminfo[0].shape;
      *__Pyx_BufPtrStrided1d(double *, __pyx_pybuffernd_suboutarr.rcbuffer->pybuffer.buf, __pyx_t_22, __pyx_pybuffernd_suboutarr.diminfo[0].strides) = __pyx_t_23;
    }

+220:     return outarr

  __Pyx_XDECREF((PyObject *)__pyx_r);
  __Pyx_INCREF((PyObject *)__pyx_v_outarr);
  __pyx_r = __pyx_v_outarr;
  __Pyx_TraceReturnValue((PyObject *)__pyx_r, 67, 0, __PYX_ERR(0, 220, __pyx_L1_error));
  goto __pyx_L0;

 221: 

 222: # Table 6 and Equations 24-26

+223: def fd12(x):

/* Python wrapper */
static PyObject *__pyx_pw_9pynitride_4core_12cython_maths_3fd12(PyObject *__pyx_self, 
#if CYTHON_METH_FASTCALL
PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds
#else
PyObject *__pyx_args, PyObject *__pyx_kwds
#endif
); /*proto*/
PyDoc_STRVAR(__pyx_doc_9pynitride_4core_12cython_maths_2fd12, "Implements the Fermi-Dirac integral of order 1/2 (Van Halen and Pulfrey approximation).\n\n    This is essentially as given in Table 6 and Equations 24-26 of\n    `(Wong et al Solid-State Electronics 1994) <http://dx.doi.org/10.1016/0038-1101(94)90105-8>`_, except that\n    for very negative arguments (:math:`x < -10`), this function will replace the sum over integrals by simply\n    the first term (:math:`e^x`).  Since, in a typical simulation, the Fermi-Dirac integral is often computed for\n    arguments with very negative :math:`x` (ie Fermi-level in midgap), this often results in a several-times\n    speedup.\n\n    Args:\n        x: the argument to the Fermi-Dirac 1/2 integral, as a numpy array.\n\n    Returns:\n        the evaluation, as a numpy array.\n    ");
static PyMethodDef __pyx_mdef_9pynitride_4core_12cython_maths_3fd12 = {"fd12", (PyCFunction)(void(*)(void))(__Pyx_PyCFunction_FastCallWithKeywords)__pyx_pw_9pynitride_4core_12cython_maths_3fd12, __Pyx_METH_FASTCALL|METH_KEYWORDS, __pyx_doc_9pynitride_4core_12cython_maths_2fd12};
static PyObject *__pyx_pw_9pynitride_4core_12cython_maths_3fd12(PyObject *__pyx_self, 
#if CYTHON_METH_FASTCALL
PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds
#else
PyObject *__pyx_args, PyObject *__pyx_kwds
#endif
) {
  PyObject *__pyx_v_x = 0;
  #if !CYTHON_METH_FASTCALL
  CYTHON_UNUSED Py_ssize_t __pyx_nargs;
  #endif
  CYTHON_UNUSED PyObject *const *__pyx_kwvalues;
  PyObject *__pyx_r = 0;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("fd12 (wrapper)", 0);
  #if !CYTHON_METH_FASTCALL
  #if CYTHON_ASSUME_SAFE_SIZE
  __pyx_nargs = PyTuple_GET_SIZE(__pyx_args);
  #else
  __pyx_nargs = PyTuple_Size(__pyx_args); if (unlikely(__pyx_nargs < 0)) return NULL;
  #endif
  #endif
  __pyx_kwvalues = __Pyx_KwValues_FASTCALL(__pyx_args, __pyx_nargs);
  {
    PyObject ** const __pyx_pyargnames[] = {&__pyx_mstate_global->__pyx_n_u_x,0};
  PyObject* values[1] = {0};
    const Py_ssize_t __pyx_kwds_len = (__pyx_kwds) ? __Pyx_NumKwargs_FASTCALL(__pyx_kwds) : 0;
    if (unlikely(__pyx_kwds_len) < 0) __PYX_ERR(0, 223, __pyx_L3_error)
    if (__pyx_kwds_len > 0) {
      switch (__pyx_nargs) {
        case  1:
        values[0] = __Pyx_ArgRef_FASTCALL(__pyx_args, 0);
        if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[0])) __PYX_ERR(0, 223, __pyx_L3_error)
        CYTHON_FALLTHROUGH;
        case  0: break;
        default: goto __pyx_L5_argtuple_error;
      }
      const Py_ssize_t kwd_pos_args = __pyx_nargs;
      if (__Pyx_ParseKeywords(__pyx_kwds, __pyx_kwvalues, __pyx_pyargnames, 0, values, kwd_pos_args, __pyx_kwds_len, "fd12", 0) < (0)) __PYX_ERR(0, 223, __pyx_L3_error)
      for (Py_ssize_t i = __pyx_nargs; i < 1; i++) {
        if (unlikely(!values[i])) { __Pyx_RaiseArgtupleInvalid("fd12", 1, 1, 1, i); __PYX_ERR(0, 223, __pyx_L3_error) }
      }
    } else if (unlikely(__pyx_nargs != 1)) {
      goto __pyx_L5_argtuple_error;
    } else {
      values[0] = __Pyx_ArgRef_FASTCALL(__pyx_args, 0);
      if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[0])) __PYX_ERR(0, 223, __pyx_L3_error)
    }
    __pyx_v_x = values[0];
  }
  goto __pyx_L6_skip;
  __pyx_L5_argtuple_error:;
  __Pyx_RaiseArgtupleInvalid("fd12", 1, 1, 1, __pyx_nargs); __PYX_ERR(0, 223, __pyx_L3_error)
  __pyx_L6_skip:;
  goto __pyx_L4_argument_unpacking_done;
  __pyx_L3_error:;
  for (Py_ssize_t __pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) {
    Py_XDECREF(values[__pyx_temp]);
  }
  __Pyx_AddTraceback("pynitride.core.cython_maths.fd12", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __Pyx_RefNannyFinishContext();
  return NULL;
  __pyx_L4_argument_unpacking_done:;
  __pyx_r = __pyx_pf_9pynitride_4core_12cython_maths_2fd12(__pyx_self, __pyx_v_x);
  int __pyx_lineno = 0;
  const char *__pyx_filename = NULL;
  int __pyx_clineno = 0;

  /* function exit code */
  for (Py_ssize_t __pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) {
    Py_XDECREF(values[__pyx_temp]);
  }
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

static PyObject *__pyx_pf_9pynitride_4core_12cython_maths_2fd12(CYTHON_UNUSED PyObject *__pyx_self, PyObject *__pyx_v_x) {
  PyObject *__pyx_r = NULL;
  __Pyx_TraceDeclarationsFunc
  __Pyx_TraceFrameInit(((PyObject *)__pyx_mstate_global->__pyx_codeobj_tab[45]))
  __Pyx_TraceStartFunc("fd12", __pyx_f[0], 223, 0, 0, 0, __PYX_ERR(0, 223, __pyx_L1_error));
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_1);
  __Pyx_XDECREF(__pyx_t_2);
  __Pyx_XDECREF(__pyx_t_3);
  __Pyx_XDECREF(__pyx_t_4);
  __Pyx_XDECREF(__pyx_t_5);
  __Pyx_TraceException(__pyx_lineno, 0, 0);
  #if CYTHON_USE_SYS_MONITORING
  __Pyx_TraceExceptionUnwind(0, 0);
  #else
  __Pyx_TraceReturnValue(NULL, 0, 0, __PYX_ERR(0, 223, __pyx_L1_error));
  #endif
  __Pyx_AddTraceback("pynitride.core.cython_maths.fd12", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = NULL;
  __pyx_L0:;
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_PyMonitoring_ExitScope(0);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
/* … */
  __pyx_t_4 = __Pyx_CyFunction_New(&__pyx_mdef_9pynitride_4core_12cython_maths_3fd12, 0, __pyx_mstate_global->__pyx_n_u_fd12, NULL, __pyx_mstate_global->__pyx_n_u_pynitride_core_cython_maths, __pyx_mstate_global->__pyx_d, ((PyObject *)__pyx_mstate_global->__pyx_codeobj_tab[45])); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 223, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  #if CYTHON_COMPILING_IN_CPYTHON && PY_VERSION_HEX >= 0x030E0000
  PyUnstable_Object_EnableDeferredRefcount(__pyx_t_4);
  #endif
  if (PyDict_SetItem(__pyx_mstate_global->__pyx_d, __pyx_mstate_global->__pyx_n_u_fd12, __pyx_t_4) < (0)) __PYX_ERR(0, 223, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;

 224:     r"""Implements the Fermi-Dirac integral of order 1/2 (Van Halen and Pulfrey approximation).

 225: 

 226:     This is essentially as given in Table 6 and Equations 24-26 of

 227:     `(Wong et al Solid-State Electronics 1994) <http://dx.doi.org/10.1016/0038-1101(94)90105-8>`_, except that

 228:     for very negative arguments (:math:`x < -10`), this function will replace the sum over integrals by simply

 229:     the first term (:math:`e^x`).  Since, in a typical simulation, the Fermi-Dirac integral is often computed for

 230:     arguments with very negative :math:`x` (ie Fermi-level in midgap), this often results in a several-times

 231:     speedup.

 232: 

 233:     Args:

 234:         x: the argument to the Fermi-Dirac 1/2 integral, as a numpy array.

 235: 

 236:     Returns:

 237:         the evaluation, as a numpy array.

 238:     """

+239:     return map1(np.asarray(x,dtype=np.double),fd12_scalar)

  __Pyx_XDECREF(__pyx_r);
  __pyx_t_2 = NULL;
  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 239, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_asarray); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 239, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 239, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_double); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 239, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_6 = 1;
  #if CYTHON_UNPACK_METHODS
  if (unlikely(PyMethod_Check(__pyx_t_4))) {
    __pyx_t_2 = PyMethod_GET_SELF(__pyx_t_4);
    assert(__pyx_t_2);
    PyObject* __pyx__function = PyMethod_GET_FUNCTION(__pyx_t_4);
    __Pyx_INCREF(__pyx_t_2);
    __Pyx_INCREF(__pyx__function);
    __Pyx_DECREF_SET(__pyx_t_4, __pyx__function);
    __pyx_t_6 = 0;
  }
  #endif
  {
    PyObject *__pyx_callargs[2 + ((CYTHON_VECTORCALL) ? 1 : 0)] = {__pyx_t_2, __pyx_v_x};
    __pyx_t_3 = __Pyx_MakeVectorcallBuilderKwds(1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 239, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_3);
    if (__Pyx_VectorcallBuilder_AddArg(__pyx_mstate_global->__pyx_n_u_dtype, __pyx_t_5, __pyx_t_3, __pyx_callargs+2, 0) < (0)) __PYX_ERR(0, 239, __pyx_L1_error)
    __pyx_t_1 = __Pyx_Object_Vectorcall_CallFromBuilder((PyObject*)__pyx_t_4, __pyx_callargs+__pyx_t_6, (2-__pyx_t_6) | (__pyx_t_6*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET), __pyx_t_3);
    __Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 239, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
  }
  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_mstate_global->__pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 239, __pyx_L1_error)
  __pyx_t_4 = ((PyObject *)__pyx_f_9pynitride_4core_12cython_maths_map1(((PyArrayObject *)__pyx_t_1), __pyx_f_9pynitride_4core_12cython_maths_fd12_scalar, NULL)); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 239, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_r = __pyx_t_4;
  __pyx_t_4 = 0;
  __Pyx_TraceReturnValue(__pyx_r, 1, 0, __PYX_ERR(0, 239, __pyx_L1_error));
  goto __pyx_L0;

 240: 

+241: def fd12p(x):

/* Python wrapper */
static PyObject *__pyx_pw_9pynitride_4core_12cython_maths_5fd12p(PyObject *__pyx_self, 
#if CYTHON_METH_FASTCALL
PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds
#else
PyObject *__pyx_args, PyObject *__pyx_kwds
#endif
); /*proto*/
PyDoc_STRVAR(__pyx_doc_9pynitride_4core_12cython_maths_4fd12p, "Implements the derivative of the :py:func:`~pynitride.poissolve.maths.fd12`.\n\n    Computes :math:`\\mathcal{F}_{1/2}'(x)`, which is equal to :math:`\\mathcal{F}_{-1/2}(x)`.  The appropriate sums\n    were obtained by simply differentiating the expressions used to form :py:func:`~pynitride.poissolve.maths.fd12`.\n\n    Args:\n        the argument to the Fermi-Dirac -1/2 integral, as a numpy array.\n    Returns:\n        the evaluation, as a numpy array.\n    ");
static PyMethodDef __pyx_mdef_9pynitride_4core_12cython_maths_5fd12p = {"fd12p", (PyCFunction)(void(*)(void))(__Pyx_PyCFunction_FastCallWithKeywords)__pyx_pw_9pynitride_4core_12cython_maths_5fd12p, __Pyx_METH_FASTCALL|METH_KEYWORDS, __pyx_doc_9pynitride_4core_12cython_maths_4fd12p};
static PyObject *__pyx_pw_9pynitride_4core_12cython_maths_5fd12p(PyObject *__pyx_self, 
#if CYTHON_METH_FASTCALL
PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds
#else
PyObject *__pyx_args, PyObject *__pyx_kwds
#endif
) {
  PyObject *__pyx_v_x = 0;
  #if !CYTHON_METH_FASTCALL
  CYTHON_UNUSED Py_ssize_t __pyx_nargs;
  #endif
  CYTHON_UNUSED PyObject *const *__pyx_kwvalues;
  PyObject *__pyx_r = 0;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("fd12p (wrapper)", 0);
  #if !CYTHON_METH_FASTCALL
  #if CYTHON_ASSUME_SAFE_SIZE
  __pyx_nargs = PyTuple_GET_SIZE(__pyx_args);
  #else
  __pyx_nargs = PyTuple_Size(__pyx_args); if (unlikely(__pyx_nargs < 0)) return NULL;
  #endif
  #endif
  __pyx_kwvalues = __Pyx_KwValues_FASTCALL(__pyx_args, __pyx_nargs);
  {
    PyObject ** const __pyx_pyargnames[] = {&__pyx_mstate_global->__pyx_n_u_x,0};
  PyObject* values[1] = {0};
    const Py_ssize_t __pyx_kwds_len = (__pyx_kwds) ? __Pyx_NumKwargs_FASTCALL(__pyx_kwds) : 0;
    if (unlikely(__pyx_kwds_len) < 0) __PYX_ERR(0, 241, __pyx_L3_error)
    if (__pyx_kwds_len > 0) {
      switch (__pyx_nargs) {
        case  1:
        values[0] = __Pyx_ArgRef_FASTCALL(__pyx_args, 0);
        if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[0])) __PYX_ERR(0, 241, __pyx_L3_error)
        CYTHON_FALLTHROUGH;
        case  0: break;
        default: goto __pyx_L5_argtuple_error;
      }
      const Py_ssize_t kwd_pos_args = __pyx_nargs;
      if (__Pyx_ParseKeywords(__pyx_kwds, __pyx_kwvalues, __pyx_pyargnames, 0, values, kwd_pos_args, __pyx_kwds_len, "fd12p", 0) < (0)) __PYX_ERR(0, 241, __pyx_L3_error)
      for (Py_ssize_t i = __pyx_nargs; i < 1; i++) {
        if (unlikely(!values[i])) { __Pyx_RaiseArgtupleInvalid("fd12p", 1, 1, 1, i); __PYX_ERR(0, 241, __pyx_L3_error) }
      }
    } else if (unlikely(__pyx_nargs != 1)) {
      goto __pyx_L5_argtuple_error;
    } else {
      values[0] = __Pyx_ArgRef_FASTCALL(__pyx_args, 0);
      if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[0])) __PYX_ERR(0, 241, __pyx_L3_error)
    }
    __pyx_v_x = values[0];
  }
  goto __pyx_L6_skip;
  __pyx_L5_argtuple_error:;
  __Pyx_RaiseArgtupleInvalid("fd12p", 1, 1, 1, __pyx_nargs); __PYX_ERR(0, 241, __pyx_L3_error)
  __pyx_L6_skip:;
  goto __pyx_L4_argument_unpacking_done;
  __pyx_L3_error:;
  for (Py_ssize_t __pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) {
    Py_XDECREF(values[__pyx_temp]);
  }
  __Pyx_AddTraceback("pynitride.core.cython_maths.fd12p", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __Pyx_RefNannyFinishContext();
  return NULL;
  __pyx_L4_argument_unpacking_done:;
  __pyx_r = __pyx_pf_9pynitride_4core_12cython_maths_4fd12p(__pyx_self, __pyx_v_x);
  int __pyx_lineno = 0;
  const char *__pyx_filename = NULL;
  int __pyx_clineno = 0;

  /* function exit code */
  for (Py_ssize_t __pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) {
    Py_XDECREF(values[__pyx_temp]);
  }
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

static PyObject *__pyx_pf_9pynitride_4core_12cython_maths_4fd12p(CYTHON_UNUSED PyObject *__pyx_self, PyObject *__pyx_v_x) {
  PyObject *__pyx_r = NULL;
  __Pyx_TraceDeclarationsFunc
  __Pyx_TraceFrameInit(((PyObject *)__pyx_mstate_global->__pyx_codeobj_tab[46]))
  __Pyx_TraceStartFunc("fd12p", __pyx_f[0], 241, 0, 0, 0, __PYX_ERR(0, 241, __pyx_L1_error));
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_1);
  __Pyx_XDECREF(__pyx_t_2);
  __Pyx_XDECREF(__pyx_t_3);
  __Pyx_XDECREF(__pyx_t_4);
  __Pyx_XDECREF(__pyx_t_5);
  __Pyx_TraceException(__pyx_lineno, 0, 0);
  #if CYTHON_USE_SYS_MONITORING
  __Pyx_TraceExceptionUnwind(0, 0);
  #else
  __Pyx_TraceReturnValue(NULL, 0, 0, __PYX_ERR(0, 241, __pyx_L1_error));
  #endif
  __Pyx_AddTraceback("pynitride.core.cython_maths.fd12p", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = NULL;
  __pyx_L0:;
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_PyMonitoring_ExitScope(0);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
/* … */
  __pyx_t_4 = __Pyx_CyFunction_New(&__pyx_mdef_9pynitride_4core_12cython_maths_5fd12p, 0, __pyx_mstate_global->__pyx_n_u_fd12p, NULL, __pyx_mstate_global->__pyx_n_u_pynitride_core_cython_maths, __pyx_mstate_global->__pyx_d, ((PyObject *)__pyx_mstate_global->__pyx_codeobj_tab[46])); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 241, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  #if CYTHON_COMPILING_IN_CPYTHON && PY_VERSION_HEX >= 0x030E0000
  PyUnstable_Object_EnableDeferredRefcount(__pyx_t_4);
  #endif
  if (PyDict_SetItem(__pyx_mstate_global->__pyx_d, __pyx_mstate_global->__pyx_n_u_fd12p, __pyx_t_4) < (0)) __PYX_ERR(0, 241, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;

 242:     r"""Implements the derivative of the :py:func:`~pynitride.poissolve.maths.fd12`.

 243: 

 244:     Computes :math:`\mathcal{F}_{1/2}'(x)`, which is equal to :math:`\mathcal{F}_{-1/2}(x)`.  The appropriate sums

 245:     were obtained by simply differentiating the expressions used to form :py:func:`~pynitride.poissolve.maths.fd12`.

 246: 

 247:     Args:

 248:         the argument to the Fermi-Dirac -1/2 integral, as a numpy array.

 249:     Returns:

 250:         the evaluation, as a numpy array.

 251:     """

+252:     return map1(np.asarray(x,dtype=np.double),fd12p_scalar)

  __Pyx_XDECREF(__pyx_r);
  __pyx_t_2 = NULL;
  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 252, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_asarray); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 252, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 252, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_double); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 252, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_6 = 1;
  #if CYTHON_UNPACK_METHODS
  if (unlikely(PyMethod_Check(__pyx_t_4))) {
    __pyx_t_2 = PyMethod_GET_SELF(__pyx_t_4);
    assert(__pyx_t_2);
    PyObject* __pyx__function = PyMethod_GET_FUNCTION(__pyx_t_4);
    __Pyx_INCREF(__pyx_t_2);
    __Pyx_INCREF(__pyx__function);
    __Pyx_DECREF_SET(__pyx_t_4, __pyx__function);
    __pyx_t_6 = 0;
  }
  #endif
  {
    PyObject *__pyx_callargs[2 + ((CYTHON_VECTORCALL) ? 1 : 0)] = {__pyx_t_2, __pyx_v_x};
    __pyx_t_3 = __Pyx_MakeVectorcallBuilderKwds(1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 252, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_3);
    if (__Pyx_VectorcallBuilder_AddArg(__pyx_mstate_global->__pyx_n_u_dtype, __pyx_t_5, __pyx_t_3, __pyx_callargs+2, 0) < (0)) __PYX_ERR(0, 252, __pyx_L1_error)
    __pyx_t_1 = __Pyx_Object_Vectorcall_CallFromBuilder((PyObject*)__pyx_t_4, __pyx_callargs+__pyx_t_6, (2-__pyx_t_6) | (__pyx_t_6*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET), __pyx_t_3);
    __Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 252, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
  }
  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_mstate_global->__pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 252, __pyx_L1_error)
  __pyx_t_4 = ((PyObject *)__pyx_f_9pynitride_4core_12cython_maths_map1(((PyArrayObject *)__pyx_t_1), __pyx_f_9pynitride_4core_12cython_maths_fd12p_scalar, NULL)); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 252, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_r = __pyx_t_4;
  __pyx_t_4 = 0;
  __Pyx_TraceReturnValue(__pyx_r, 1, 0, __PYX_ERR(0, 252, __pyx_L1_error));
  goto __pyx_L0;

 253: 

 254: 

 255: 

 256: # See comments for idd

+257: cdef double idd_scalar(double eta, double g):

static double __pyx_f_9pynitride_4core_12cython_maths_idd_scalar(double __pyx_v_eta, double __pyx_v_g) {
  double __pyx_r;
  __Pyx_TraceDeclarationsFunc
  __Pyx_TraceFrameInit(((PyObject *)__pyx_mstate_global->__pyx_codeobj_tab[47]))
  __Pyx_TraceStartFunc("idd_scalar", __pyx_f[0], 257, 0, 0, 0, __PYX_ERR(0, 257, __pyx_L1_error));
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_TraceException(__pyx_lineno, 0, 0);
  #if CYTHON_USE_SYS_MONITORING
  __Pyx_TraceExceptionUnwind(0, 0);
  #else
  __Pyx_TraceReturnValue(NULL, 0, 0, __PYX_ERR(0, 257, __pyx_L1_error));
  #endif
  __Pyx_AddTraceback("pynitride.core.cython_maths.idd_scalar", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = -1;
  __pyx_L0:;
  __Pyx_PyMonitoring_ExitScope(0);
  return __pyx_r;
}

+258:     if eta>500:

  __pyx_t_1 = (__pyx_v_eta > 500.0);
  if (__pyx_t_1) {
/* … */
  }

+259:         return 0

    __pyx_r = 0.0;
    __Pyx_TraceReturnCValue(__pyx_r, PyFloat_FromDouble, 6, 0, __PYX_ERR(0, 259, __pyx_L1_error));
    goto __pyx_L0;

 260:     else:

+261:         return 1/(1+g*exp(eta))

  /*else*/ {
    __pyx_t_2 = (1.0 + (__pyx_v_g * exp(__pyx_v_eta)));
    if (unlikely(__pyx_t_2 == 0)) {
      PyErr_SetString(PyExc_ZeroDivisionError, "float division");
      __PYX_ERR(0, 261, __pyx_L1_error)
    }
    __pyx_r = (1.0 / __pyx_t_2);
    __Pyx_TraceReturnCValue(__pyx_r, PyFloat_FromDouble, 8, 0, __PYX_ERR(0, 261, __pyx_L1_error));
    goto __pyx_L0;
  }

 262: 

+263: def idd(eta,g):

/* Python wrapper */
static PyObject *__pyx_pw_9pynitride_4core_12cython_maths_7idd(PyObject *__pyx_self, 
#if CYTHON_METH_FASTCALL
PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds
#else
PyObject *__pyx_args, PyObject *__pyx_kwds
#endif
); /*proto*/
PyDoc_STRVAR(__pyx_doc_9pynitride_4core_12cython_maths_6idd, "\n    Computes the ionized dopant density given the normalized Fermi position and degeneracy factor\n\n\n    See Tiwari Compound Semiconductor Devices pg 31-32.  Prevents underflow by zeroing at large eta.\n\n    Args:\n        eta: For donors,   this should be `\\eta=(E_F-E_c + E)/kT`\n                where `E` is the dopant energy down from the conduction edge\n            For acceptors, this should be `\\eta=(E_v-E_F + E)/kT`\n                where `E` is the dopant energy down from the valence edge\n        g: degeneracy factor (as given by Tiwari)\n\n    Returns:\n        the evaluation as a numpy array\n    ");
static PyMethodDef __pyx_mdef_9pynitride_4core_12cython_maths_7idd = {"idd", (PyCFunction)(void(*)(void))(__Pyx_PyCFunction_FastCallWithKeywords)__pyx_pw_9pynitride_4core_12cython_maths_7idd, __Pyx_METH_FASTCALL|METH_KEYWORDS, __pyx_doc_9pynitride_4core_12cython_maths_6idd};
static PyObject *__pyx_pw_9pynitride_4core_12cython_maths_7idd(PyObject *__pyx_self, 
#if CYTHON_METH_FASTCALL
PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds
#else
PyObject *__pyx_args, PyObject *__pyx_kwds
#endif
) {
  PyObject *__pyx_v_eta = 0;
  PyObject *__pyx_v_g = 0;
  #if !CYTHON_METH_FASTCALL
  CYTHON_UNUSED Py_ssize_t __pyx_nargs;
  #endif
  CYTHON_UNUSED PyObject *const *__pyx_kwvalues;
  PyObject *__pyx_r = 0;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("idd (wrapper)", 0);
  #if !CYTHON_METH_FASTCALL
  #if CYTHON_ASSUME_SAFE_SIZE
  __pyx_nargs = PyTuple_GET_SIZE(__pyx_args);
  #else
  __pyx_nargs = PyTuple_Size(__pyx_args); if (unlikely(__pyx_nargs < 0)) return NULL;
  #endif
  #endif
  __pyx_kwvalues = __Pyx_KwValues_FASTCALL(__pyx_args, __pyx_nargs);
  {
    PyObject ** const __pyx_pyargnames[] = {&__pyx_mstate_global->__pyx_n_u_eta,&__pyx_mstate_global->__pyx_n_u_g,0};
  PyObject* values[2] = {0,0};
    const Py_ssize_t __pyx_kwds_len = (__pyx_kwds) ? __Pyx_NumKwargs_FASTCALL(__pyx_kwds) : 0;
    if (unlikely(__pyx_kwds_len) < 0) __PYX_ERR(0, 263, __pyx_L3_error)
    if (__pyx_kwds_len > 0) {
      switch (__pyx_nargs) {
        case  2:
        values[1] = __Pyx_ArgRef_FASTCALL(__pyx_args, 1);
        if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[1])) __PYX_ERR(0, 263, __pyx_L3_error)
        CYTHON_FALLTHROUGH;
        case  1:
        values[0] = __Pyx_ArgRef_FASTCALL(__pyx_args, 0);
        if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[0])) __PYX_ERR(0, 263, __pyx_L3_error)
        CYTHON_FALLTHROUGH;
        case  0: break;
        default: goto __pyx_L5_argtuple_error;
      }
      const Py_ssize_t kwd_pos_args = __pyx_nargs;
      if (__Pyx_ParseKeywords(__pyx_kwds, __pyx_kwvalues, __pyx_pyargnames, 0, values, kwd_pos_args, __pyx_kwds_len, "idd", 0) < (0)) __PYX_ERR(0, 263, __pyx_L3_error)
      for (Py_ssize_t i = __pyx_nargs; i < 2; i++) {
        if (unlikely(!values[i])) { __Pyx_RaiseArgtupleInvalid("idd", 1, 2, 2, i); __PYX_ERR(0, 263, __pyx_L3_error) }
      }
    } else if (unlikely(__pyx_nargs != 2)) {
      goto __pyx_L5_argtuple_error;
    } else {
      values[0] = __Pyx_ArgRef_FASTCALL(__pyx_args, 0);
      if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[0])) __PYX_ERR(0, 263, __pyx_L3_error)
      values[1] = __Pyx_ArgRef_FASTCALL(__pyx_args, 1);
      if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[1])) __PYX_ERR(0, 263, __pyx_L3_error)
    }
    __pyx_v_eta = values[0];
    __pyx_v_g = values[1];
  }
  goto __pyx_L6_skip;
  __pyx_L5_argtuple_error:;
  __Pyx_RaiseArgtupleInvalid("idd", 1, 2, 2, __pyx_nargs); __PYX_ERR(0, 263, __pyx_L3_error)
  __pyx_L6_skip:;
  goto __pyx_L4_argument_unpacking_done;
  __pyx_L3_error:;
  for (Py_ssize_t __pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) {
    Py_XDECREF(values[__pyx_temp]);
  }
  __Pyx_AddTraceback("pynitride.core.cython_maths.idd", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __Pyx_RefNannyFinishContext();
  return NULL;
  __pyx_L4_argument_unpacking_done:;
  __pyx_r = __pyx_pf_9pynitride_4core_12cython_maths_6idd(__pyx_self, __pyx_v_eta, __pyx_v_g);
  int __pyx_lineno = 0;
  const char *__pyx_filename = NULL;
  int __pyx_clineno = 0;

  /* function exit code */
  for (Py_ssize_t __pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) {
    Py_XDECREF(values[__pyx_temp]);
  }
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

static PyObject *__pyx_pf_9pynitride_4core_12cython_maths_6idd(CYTHON_UNUSED PyObject *__pyx_self, PyObject *__pyx_v_eta, PyObject *__pyx_v_g) {
  PyObject *__pyx_r = NULL;
  __Pyx_TraceDeclarationsFunc
  __Pyx_TraceFrameInit(((PyObject *)__pyx_mstate_global->__pyx_codeobj_tab[48]))
  __Pyx_TraceStartFunc("idd", __pyx_f[0], 263, 0, 0, 0, __PYX_ERR(0, 263, __pyx_L1_error));
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_1);
  __Pyx_XDECREF(__pyx_t_2);
  __Pyx_XDECREF(__pyx_t_3);
  __Pyx_XDECREF(__pyx_t_4);
  __Pyx_XDECREF(__pyx_t_5);
  __Pyx_XDECREF(__pyx_t_7);
  __Pyx_TraceException(__pyx_lineno, 0, 0);
  #if CYTHON_USE_SYS_MONITORING
  __Pyx_TraceExceptionUnwind(0, 0);
  #else
  __Pyx_TraceReturnValue(NULL, 0, 0, __PYX_ERR(0, 263, __pyx_L1_error));
  #endif
  __Pyx_AddTraceback("pynitride.core.cython_maths.idd", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = NULL;
  __pyx_L0:;
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_PyMonitoring_ExitScope(0);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
/* … */
  __pyx_t_4 = __Pyx_CyFunction_New(&__pyx_mdef_9pynitride_4core_12cython_maths_7idd, 0, __pyx_mstate_global->__pyx_n_u_idd, NULL, __pyx_mstate_global->__pyx_n_u_pynitride_core_cython_maths, __pyx_mstate_global->__pyx_d, ((PyObject *)__pyx_mstate_global->__pyx_codeobj_tab[48])); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 263, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  #if CYTHON_COMPILING_IN_CPYTHON && PY_VERSION_HEX >= 0x030E0000
  PyUnstable_Object_EnableDeferredRefcount(__pyx_t_4);
  #endif
  if (PyDict_SetItem(__pyx_mstate_global->__pyx_d, __pyx_mstate_global->__pyx_n_u_idd, __pyx_t_4) < (0)) __PYX_ERR(0, 263, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;

 264:     """

 265:     Computes the ionized dopant density given the normalized Fermi position and degeneracy factor

 266: 

 267: 

 268:     See Tiwari Compound Semiconductor Devices pg 31-32.  Prevents underflow by zeroing at large eta.

 269: 

 270:     Args:

 271:         eta: For donors,   this should be `\eta=(E_F-E_c + E)/kT`

 272:                 where `E` is the dopant energy down from the conduction edge

 273:             For acceptors, this should be `\eta=(E_v-E_F + E)/kT`

 274:                 where `E` is the dopant energy down from the valence edge

 275:         g: degeneracy factor (as given by Tiwari)

 276: 

 277:     Returns:

 278:         the evaluation as a numpy array

 279:     """

+280:     return map2(np.asarray(eta,dtype=np.double),np.asarray(g,dtype=np.double),idd_scalar)

  __Pyx_XDECREF(__pyx_r);
  __pyx_t_2 = NULL;
  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 280, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_asarray); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 280, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 280, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_double); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 280, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_6 = 1;
  #if CYTHON_UNPACK_METHODS
  if (unlikely(PyMethod_Check(__pyx_t_4))) {
    __pyx_t_2 = PyMethod_GET_SELF(__pyx_t_4);
    assert(__pyx_t_2);
    PyObject* __pyx__function = PyMethod_GET_FUNCTION(__pyx_t_4);
    __Pyx_INCREF(__pyx_t_2);
    __Pyx_INCREF(__pyx__function);
    __Pyx_DECREF_SET(__pyx_t_4, __pyx__function);
    __pyx_t_6 = 0;
  }
  #endif
  {
    PyObject *__pyx_callargs[2 + ((CYTHON_VECTORCALL) ? 1 : 0)] = {__pyx_t_2, __pyx_v_eta};
    __pyx_t_3 = __Pyx_MakeVectorcallBuilderKwds(1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 280, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_3);
    if (__Pyx_VectorcallBuilder_AddArg(__pyx_mstate_global->__pyx_n_u_dtype, __pyx_t_5, __pyx_t_3, __pyx_callargs+2, 0) < (0)) __PYX_ERR(0, 280, __pyx_L1_error)
    __pyx_t_1 = __Pyx_Object_Vectorcall_CallFromBuilder((PyObject*)__pyx_t_4, __pyx_callargs+__pyx_t_6, (2-__pyx_t_6) | (__pyx_t_6*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET), __pyx_t_3);
    __Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 280, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
  }
  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_mstate_global->__pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 280, __pyx_L1_error)
  __pyx_t_3 = NULL;
  __Pyx_GetModuleGlobalName(__pyx_t_5, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 280, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_5, __pyx_mstate_global->__pyx_n_u_asarray); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 280, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __Pyx_GetModuleGlobalName(__pyx_t_5, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 280, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __pyx_t_7 = __Pyx_PyObject_GetAttrStr(__pyx_t_5, __pyx_mstate_global->__pyx_n_u_double); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 280, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_7);
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __pyx_t_6 = 1;
  #if CYTHON_UNPACK_METHODS
  if (unlikely(PyMethod_Check(__pyx_t_2))) {
    __pyx_t_3 = PyMethod_GET_SELF(__pyx_t_2);
    assert(__pyx_t_3);
    PyObject* __pyx__function = PyMethod_GET_FUNCTION(__pyx_t_2);
    __Pyx_INCREF(__pyx_t_3);
    __Pyx_INCREF(__pyx__function);
    __Pyx_DECREF_SET(__pyx_t_2, __pyx__function);
    __pyx_t_6 = 0;
  }
  #endif
  {
    PyObject *__pyx_callargs[2 + ((CYTHON_VECTORCALL) ? 1 : 0)] = {__pyx_t_3, __pyx_v_g};
    __pyx_t_5 = __Pyx_MakeVectorcallBuilderKwds(1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 280, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_5);
    if (__Pyx_VectorcallBuilder_AddArg(__pyx_mstate_global->__pyx_n_u_dtype, __pyx_t_7, __pyx_t_5, __pyx_callargs+2, 0) < (0)) __PYX_ERR(0, 280, __pyx_L1_error)
    __pyx_t_4 = __Pyx_Object_Vectorcall_CallFromBuilder((PyObject*)__pyx_t_2, __pyx_callargs+__pyx_t_6, (2-__pyx_t_6) | (__pyx_t_6*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET), __pyx_t_5);
    __Pyx_XDECREF(__pyx_t_3); __pyx_t_3 = 0;
    __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
    if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 280, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_4);
  }
  if (!(likely(((__pyx_t_4) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_4, __pyx_mstate_global->__pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 280, __pyx_L1_error)
  __pyx_t_2 = ((PyObject *)__pyx_f_9pynitride_4core_12cython_maths_map2(((PyArrayObject *)__pyx_t_1), ((PyArrayObject *)__pyx_t_4), __pyx_f_9pynitride_4core_12cython_maths_idd_scalar, NULL)); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 280, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_r = __pyx_t_2;
  __pyx_t_2 = 0;
  __Pyx_TraceReturnValue(__pyx_r, 1, 0, __PYX_ERR(0, 280, __pyx_L1_error));
  goto __pyx_L0;

 281: 

 282: 

 283: # See comments for iddd

+284: cdef double iddd_scalar(double eta, double g):

static double __pyx_f_9pynitride_4core_12cython_maths_iddd_scalar(double __pyx_v_eta, double __pyx_v_g) {
  double __pyx_r;
  __Pyx_TraceDeclarationsFunc
  __Pyx_TraceFrameInit(((PyObject *)__pyx_mstate_global->__pyx_codeobj_tab[49]))
  __Pyx_TraceStartFunc("iddd_scalar", __pyx_f[0], 284, 0, 0, 0, __PYX_ERR(0, 284, __pyx_L1_error));
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_TraceException(__pyx_lineno, 0, 0);
  #if CYTHON_USE_SYS_MONITORING
  __Pyx_TraceExceptionUnwind(0, 0);
  #else
  __Pyx_TraceReturnValue(NULL, 0, 0, __PYX_ERR(0, 284, __pyx_L1_error));
  #endif
  __Pyx_AddTraceback("pynitride.core.cython_maths.iddd_scalar", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = -1;
  __pyx_L0:;
  __Pyx_PyMonitoring_ExitScope(0);
  return __pyx_r;
}

+285:     if eta>500:

  __pyx_t_1 = (__pyx_v_eta > 500.0);
  if (__pyx_t_1) {
/* … */
  }

+286:         return 0

    __pyx_r = 0.0;
    __Pyx_TraceReturnCValue(__pyx_r, PyFloat_FromDouble, 6, 0, __PYX_ERR(0, 286, __pyx_L1_error));
    goto __pyx_L0;

 287:     else:

+288:         return g*exp(eta)/(1+g*exp(eta))**2

  /*else*/ {
    __pyx_t_2 = (__pyx_v_g * exp(__pyx_v_eta));
    __pyx_t_3 = pow((1.0 + (__pyx_v_g * exp(__pyx_v_eta))), 2.0);
    if (unlikely(__pyx_t_3 == 0)) {
      PyErr_SetString(PyExc_ZeroDivisionError, "float division");
      __PYX_ERR(0, 288, __pyx_L1_error)
    }
    __pyx_r = (__pyx_t_2 / __pyx_t_3);
    __Pyx_TraceReturnCValue(__pyx_r, PyFloat_FromDouble, 8, 0, __PYX_ERR(0, 288, __pyx_L1_error));
    goto __pyx_L0;
  }

 289: 

+290: def iddd(eta,g):

/* Python wrapper */
static PyObject *__pyx_pw_9pynitride_4core_12cython_maths_9iddd(PyObject *__pyx_self, 
#if CYTHON_METH_FASTCALL
PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds
#else
PyObject *__pyx_args, PyObject *__pyx_kwds
#endif
); /*proto*/
PyDoc_STRVAR(__pyx_doc_9pynitride_4core_12cython_maths_8iddd, "\n    Computes the derivative of the ionized dopant density (with respect to eta),\n    given the normalized Fermi position and degeneracy factor\n\n\n    See Tiwari Compound Semiconductor Devices pg 31-32.  Prevents underflow by zeroing at large eta.\n\n    Args:\n        eta: For donors,   this should be `\\eta=(E_F-E_c + E)/kT`\n                where `E` is the dopant energy down from the conduction edge\n            For acceptors, this should be `\\eta=(E_v-E_F + E)/kT`\n                where `E` is the dopant energy down from the valence edge\n        g: degeneracy factor (as given by Tiwari)\n\n    Returns:\n        the evaluation as a numpy array\n    ");
static PyMethodDef __pyx_mdef_9pynitride_4core_12cython_maths_9iddd = {"iddd", (PyCFunction)(void(*)(void))(__Pyx_PyCFunction_FastCallWithKeywords)__pyx_pw_9pynitride_4core_12cython_maths_9iddd, __Pyx_METH_FASTCALL|METH_KEYWORDS, __pyx_doc_9pynitride_4core_12cython_maths_8iddd};
static PyObject *__pyx_pw_9pynitride_4core_12cython_maths_9iddd(PyObject *__pyx_self, 
#if CYTHON_METH_FASTCALL
PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds
#else
PyObject *__pyx_args, PyObject *__pyx_kwds
#endif
) {
  PyObject *__pyx_v_eta = 0;
  PyObject *__pyx_v_g = 0;
  #if !CYTHON_METH_FASTCALL
  CYTHON_UNUSED Py_ssize_t __pyx_nargs;
  #endif
  CYTHON_UNUSED PyObject *const *__pyx_kwvalues;
  PyObject *__pyx_r = 0;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("iddd (wrapper)", 0);
  #if !CYTHON_METH_FASTCALL
  #if CYTHON_ASSUME_SAFE_SIZE
  __pyx_nargs = PyTuple_GET_SIZE(__pyx_args);
  #else
  __pyx_nargs = PyTuple_Size(__pyx_args); if (unlikely(__pyx_nargs < 0)) return NULL;
  #endif
  #endif
  __pyx_kwvalues = __Pyx_KwValues_FASTCALL(__pyx_args, __pyx_nargs);
  {
    PyObject ** const __pyx_pyargnames[] = {&__pyx_mstate_global->__pyx_n_u_eta,&__pyx_mstate_global->__pyx_n_u_g,0};
  PyObject* values[2] = {0,0};
    const Py_ssize_t __pyx_kwds_len = (__pyx_kwds) ? __Pyx_NumKwargs_FASTCALL(__pyx_kwds) : 0;
    if (unlikely(__pyx_kwds_len) < 0) __PYX_ERR(0, 290, __pyx_L3_error)
    if (__pyx_kwds_len > 0) {
      switch (__pyx_nargs) {
        case  2:
        values[1] = __Pyx_ArgRef_FASTCALL(__pyx_args, 1);
        if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[1])) __PYX_ERR(0, 290, __pyx_L3_error)
        CYTHON_FALLTHROUGH;
        case  1:
        values[0] = __Pyx_ArgRef_FASTCALL(__pyx_args, 0);
        if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[0])) __PYX_ERR(0, 290, __pyx_L3_error)
        CYTHON_FALLTHROUGH;
        case  0: break;
        default: goto __pyx_L5_argtuple_error;
      }
      const Py_ssize_t kwd_pos_args = __pyx_nargs;
      if (__Pyx_ParseKeywords(__pyx_kwds, __pyx_kwvalues, __pyx_pyargnames, 0, values, kwd_pos_args, __pyx_kwds_len, "iddd", 0) < (0)) __PYX_ERR(0, 290, __pyx_L3_error)
      for (Py_ssize_t i = __pyx_nargs; i < 2; i++) {
        if (unlikely(!values[i])) { __Pyx_RaiseArgtupleInvalid("iddd", 1, 2, 2, i); __PYX_ERR(0, 290, __pyx_L3_error) }
      }
    } else if (unlikely(__pyx_nargs != 2)) {
      goto __pyx_L5_argtuple_error;
    } else {
      values[0] = __Pyx_ArgRef_FASTCALL(__pyx_args, 0);
      if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[0])) __PYX_ERR(0, 290, __pyx_L3_error)
      values[1] = __Pyx_ArgRef_FASTCALL(__pyx_args, 1);
      if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[1])) __PYX_ERR(0, 290, __pyx_L3_error)
    }
    __pyx_v_eta = values[0];
    __pyx_v_g = values[1];
  }
  goto __pyx_L6_skip;
  __pyx_L5_argtuple_error:;
  __Pyx_RaiseArgtupleInvalid("iddd", 1, 2, 2, __pyx_nargs); __PYX_ERR(0, 290, __pyx_L3_error)
  __pyx_L6_skip:;
  goto __pyx_L4_argument_unpacking_done;
  __pyx_L3_error:;
  for (Py_ssize_t __pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) {
    Py_XDECREF(values[__pyx_temp]);
  }
  __Pyx_AddTraceback("pynitride.core.cython_maths.iddd", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __Pyx_RefNannyFinishContext();
  return NULL;
  __pyx_L4_argument_unpacking_done:;
  __pyx_r = __pyx_pf_9pynitride_4core_12cython_maths_8iddd(__pyx_self, __pyx_v_eta, __pyx_v_g);
  int __pyx_lineno = 0;
  const char *__pyx_filename = NULL;
  int __pyx_clineno = 0;

  /* function exit code */
  for (Py_ssize_t __pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) {
    Py_XDECREF(values[__pyx_temp]);
  }
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

static PyObject *__pyx_pf_9pynitride_4core_12cython_maths_8iddd(CYTHON_UNUSED PyObject *__pyx_self, PyObject *__pyx_v_eta, PyObject *__pyx_v_g) {
  PyObject *__pyx_r = NULL;
  __Pyx_TraceDeclarationsFunc
  __Pyx_TraceFrameInit(((PyObject *)__pyx_mstate_global->__pyx_codeobj_tab[50]))
  __Pyx_TraceStartFunc("iddd", __pyx_f[0], 290, 0, 0, 0, __PYX_ERR(0, 290, __pyx_L1_error));
/* … */
  __pyx_t_4 = __Pyx_CyFunction_New(&__pyx_mdef_9pynitride_4core_12cython_maths_9iddd, 0, __pyx_mstate_global->__pyx_n_u_iddd, NULL, __pyx_mstate_global->__pyx_n_u_pynitride_core_cython_maths, __pyx_mstate_global->__pyx_d, ((PyObject *)__pyx_mstate_global->__pyx_codeobj_tab[50])); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 290, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  #if CYTHON_COMPILING_IN_CPYTHON && PY_VERSION_HEX >= 0x030E0000
  PyUnstable_Object_EnableDeferredRefcount(__pyx_t_4);
  #endif
  if (PyDict_SetItem(__pyx_mstate_global->__pyx_d, __pyx_mstate_global->__pyx_n_u_iddd, __pyx_t_4) < (0)) __PYX_ERR(0, 290, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;

 291:     """

 292:     Computes the derivative of the ionized dopant density (with respect to eta),

 293:     given the normalized Fermi position and degeneracy factor

 294: 

 295: 

 296:     See Tiwari Compound Semiconductor Devices pg 31-32.  Prevents underflow by zeroing at large eta.

 297: 

 298:     Args:

 299:         eta: For donors,   this should be `\eta=(E_F-E_c + E)/kT`

 300:                 where `E` is the dopant energy down from the conduction edge

 301:             For acceptors, this should be `\eta=(E_v-E_F + E)/kT`

 302:                 where `E` is the dopant energy down from the valence edge

 303:         g: degeneracy factor (as given by Tiwari)

 304: 

 305:     Returns:

 306:         the evaluation as a numpy array

 307:     """

+308:     return map2(np.asarray(eta,dtype=np.double),np.asarray(g,dtype=np.double),iddd_scalar)

  __Pyx_XDECREF(__pyx_r);
  __pyx_t_2 = NULL;
  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 308, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_asarray); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 308, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 308, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_double); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 308, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_6 = 1;
  #if CYTHON_UNPACK_METHODS
  if (unlikely(PyMethod_Check(__pyx_t_4))) {
    __pyx_t_2 = PyMethod_GET_SELF(__pyx_t_4);
    assert(__pyx_t_2);
    PyObject* __pyx__function = PyMethod_GET_FUNCTION(__pyx_t_4);
    __Pyx_INCREF(__pyx_t_2);
    __Pyx_INCREF(__pyx__function);
    __Pyx_DECREF_SET(__pyx_t_4, __pyx__function);
    __pyx_t_6 = 0;
  }
  #endif
  {
    PyObject *__pyx_callargs[2 + ((CYTHON_VECTORCALL) ? 1 : 0)] = {__pyx_t_2, __pyx_v_eta};
    __pyx_t_3 = __Pyx_MakeVectorcallBuilderKwds(1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 308, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_3);
    if (__Pyx_VectorcallBuilder_AddArg(__pyx_mstate_global->__pyx_n_u_dtype, __pyx_t_5, __pyx_t_3, __pyx_callargs+2, 0) < (0)) __PYX_ERR(0, 308, __pyx_L1_error)
    __pyx_t_1 = __Pyx_Object_Vectorcall_CallFromBuilder((PyObject*)__pyx_t_4, __pyx_callargs+__pyx_t_6, (2-__pyx_t_6) | (__pyx_t_6*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET), __pyx_t_3);
    __Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 308, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
  }
  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_mstate_global->__pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 308, __pyx_L1_error)
  __pyx_t_3 = NULL;
  __Pyx_GetModuleGlobalName(__pyx_t_5, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 308, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_5, __pyx_mstate_global->__pyx_n_u_asarray); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 308, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __Pyx_GetModuleGlobalName(__pyx_t_5, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 308, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __pyx_t_7 = __Pyx_PyObject_GetAttrStr(__pyx_t_5, __pyx_mstate_global->__pyx_n_u_double); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 308, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_7);
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __pyx_t_6 = 1;
  #if CYTHON_UNPACK_METHODS
  if (unlikely(PyMethod_Check(__pyx_t_2))) {
    __pyx_t_3 = PyMethod_GET_SELF(__pyx_t_2);
    assert(__pyx_t_3);
    PyObject* __pyx__function = PyMethod_GET_FUNCTION(__pyx_t_2);
    __Pyx_INCREF(__pyx_t_3);
    __Pyx_INCREF(__pyx__function);
    __Pyx_DECREF_SET(__pyx_t_2, __pyx__function);
    __pyx_t_6 = 0;
  }
  #endif
  {
    PyObject *__pyx_callargs[2 + ((CYTHON_VECTORCALL) ? 1 : 0)] = {__pyx_t_3, __pyx_v_g};
    __pyx_t_5 = __Pyx_MakeVectorcallBuilderKwds(1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 308, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_5);
    if (__Pyx_VectorcallBuilder_AddArg(__pyx_mstate_global->__pyx_n_u_dtype, __pyx_t_7, __pyx_t_5, __pyx_callargs+2, 0) < (0)) __PYX_ERR(0, 308, __pyx_L1_error)
    __pyx_t_4 = __Pyx_Object_Vectorcall_CallFromBuilder((PyObject*)__pyx_t_2, __pyx_callargs+__pyx_t_6, (2-__pyx_t_6) | (__pyx_t_6*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET), __pyx_t_5);
    __Pyx_XDECREF(__pyx_t_3); __pyx_t_3 = 0;
    __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
    if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 308, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_4);
  }
  if (!(likely(((__pyx_t_4) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_4, __pyx_mstate_global->__pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 308, __pyx_L1_error)
  __pyx_t_2 = ((PyObject *)__pyx_f_9pynitride_4core_12cython_maths_map2(((PyArrayObject *)__pyx_t_1), ((PyArrayObject *)__pyx_t_4), __pyx_f_9pynitride_4core_12cython_maths_iddd_scalar, NULL)); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 308, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_r = __pyx_t_2;
  __pyx_t_2 = 0;
  __Pyx_TraceReturnValue(__pyx_r, 1, 0, __PYX_ERR(0, 308, __pyx_L1_error));
  goto __pyx_L0;

 309: 

 310: