diff --git a/.github/workflows/test_with_conda.yml b/.github/workflows/test_with_conda.yml
index 19c0600d..30a20df0 100644
--- a/.github/workflows/test_with_conda.yml
+++ b/.github/workflows/test_with_conda.yml
@@ -39,6 +39,7 @@ jobs:
           matplotlib
           jupyter
           utm
+          numba
           pytest
           pytest-cov
           sphinx
@@ -87,6 +88,7 @@ jobs:
           matplotlib
           jupyter
           utm
+          numba
           pytest
           pytest-cov
           sphinx
diff --git a/.gitignore b/.gitignore
index e54e277a..0ed7dd09 100644
--- a/.gitignore
+++ b/.gitignore
@@ -98,7 +98,11 @@ docs/api/generated/*
 
 # Jupyter
 *.ipynb
+
+#Cython generated files
 geoana/kernels/_extensions/rTE.cpp
+geoana/kernels/_extensions/potential_field_prism.c
+geoana/kernels/_extensions/potential_field_prism_api.h
 
 # setuptools_scm
 geoana/version.py
diff --git a/geoana/kernels/_extensions/__init__.py b/geoana/kernels/_extensions/__init__.py
index e69de29b..0b1da38d 100644
--- a/geoana/kernels/_extensions/__init__.py
+++ b/geoana/kernels/_extensions/__init__.py
@@ -0,0 +1,54 @@
+try:
+    # register numba jitable versions of the prism functions
+    # if numba is available (and this module is installed).
+    from numba.extending import (
+        overload,
+        get_cython_function_address
+    )
+    from numba import types
+    import ctypes
+
+    from .potential_field_prism import (
+        prism_f,
+        prism_fz,
+        prism_fzz,
+        prism_fzx,
+        prism_fzy,
+        prism_fzzz,
+        prism_fxxy,
+        prism_fxxz,
+        prism_fxyz,
+    )
+    funcs = [
+        prism_f,
+        prism_fz,
+        prism_fzz,
+        prism_fzx,
+        prism_fzy,
+        prism_fzzz,
+        prism_fxxy,
+        prism_fxxz,
+        prism_fxyz,
+    ]
+
+    def _numba_register_prism_func(prism_func):
+        module = 'geoana.kernels._extensions.potential_field_prism'
+        name = prism_func.__name__
+
+        func_address = get_cython_function_address(module, name)
+        func_type = ctypes.CFUNCTYPE(ctypes.c_double, ctypes.c_double, ctypes.c_double, ctypes.c_double)
+        c_func = func_type(func_address)
+
+        @overload(prism_func)
+        def numba_func(x, y, z):
+            if isinstance(x, types.Float):
+                if isinstance(y, types.Float):
+                    if isinstance(z, types.Float):
+                        def f(x, y, z):
+                            return c_func(x, y, z)
+                        return f
+    for func in funcs:
+        _numba_register_prism_func(func)
+
+except ImportError as err:
+    pass
\ No newline at end of file
diff --git a/geoana/kernels/_extensions/potential_field_prism.c b/geoana/kernels/_extensions/potential_field_prism.c
deleted file mode 100644
index 081d423c..00000000
--- a/geoana/kernels/_extensions/potential_field_prism.c
+++ /dev/null
@@ -1,579 +0,0 @@
-#define PY_SSIZE_T_CLEAN
-#include <Python.h>
-#include "numpy/ndarraytypes.h"
-#include "numpy/ufuncobject.h"
-#include "numpy/npy_3kcompat.h"
-#include <math.h>
-
-static PyMethodDef PrismMethods[] = {
-    {NULL, NULL, 0, NULL}
-};
-
-/* The loop definition must precede the PyMODINIT_FUNC. */
-
-static void double_prism_f(char **args, const npy_intp *dimensions,
-                         const npy_intp *steps, void *data)
-{
-    npy_intp i;
-    npy_intp n = dimensions[0];
-    char *inx = args[0], *iny = args[1], *inz = args[2];
-    char *out = args[3];
-    npy_intp inx_step = steps[0], iny_step = steps[1], inz_step = steps[2];
-    npy_intp out_step = steps[3];
-
-    double x, y, z, v, r, temp;
-
-    for (i = 0; i < n; i++) {
-        /* BEGIN main ufunc computation */
-        x = *(double *)inx;
-        y = *(double *)iny;
-        z = *(double *)inz;
-        v = 0.0;
-        r = sqrt(x * x + y * y + z * z);
-        if (x != 0.0){
-            if (y != 0.0){
-                temp = z + r;
-                // check if x & y were small relative to -z
-                if (temp > 0){
-                    v -= x * y * log(temp);
-                }
-            }
-            v += 0.5 * x * x * atan( y * z / (x * r));
-        }
-        if (y != 0.0){
-            if (z != 0.0){
-                temp = x + r;
-                // check if y & z were small relative to -x
-                if (temp > 0){
-                    v -= y *z * log(temp);
-                }
-            }
-            v += 0.5 * y * y * atan(z * x / (y * r));
-        }
-        if (z != 0.0){
-            if (x != 0.0){
-                temp = y + r;
-                // check if x & z were small relative to -y
-                if (temp > 0){
-                    v -= z * x * log(temp);
-                }
-            }
-            v += 0.5 * z * z * atan(x * y / (z * r));
-        }
-        *((double *)out) = v;
-        /* END main ufunc computation */
-
-        inx += inx_step;
-        iny += iny_step;
-        inz += inz_step;
-        out += out_step;
-    }
-}
-
-/* This a pointer to the above function */
-PyUFuncGenericFunction funcs_prism_f[1] = {&double_prism_f};
-
-static void double_prism_fz(char **args, const npy_intp *dimensions,
-                         const npy_intp *steps, void *data)
-{
-    npy_intp i;
-    npy_intp n = dimensions[0];
-    char *inx = args[0], *iny = args[1], *inz = args[2];
-    char *out = args[3];
-    npy_intp inx_step = steps[0], iny_step = steps[1], inz_step = steps[2];
-    npy_intp out_step = steps[3];
-
-    double x, y, z, v, r, temp;
-
-    for (i = 0; i < n; i++) {
-        /* BEGIN main ufunc computation */
-        x = *(double *)inx;
-        y = *(double *)iny;
-        z = *(double *)inz;
-        v = 0.0;
-        r = sqrt(x * x + y * y + z * z);
-        if (x != 0.0){
-            temp = y + r;
-            // check if x & z were small relative to -y
-            if (temp > 0.0){
-                v += x * log(temp);
-            }
-        }
-        if (y != 0.0){
-            temp = x + r;
-            // check if y & z were small relative to -x
-            if (temp > 0.0){
-                v += y * log(temp);
-            }
-        }
-        if (z != 0.0){
-            v -= z * atan(x * y / (z * r));
-        }
-        *((double *)out) = v;
-        /* END main ufunc computation */
-
-        inx += inx_step;
-        iny += iny_step;
-        inz += inz_step;
-        out += out_step;
-    }
-}
-
-/* This a pointer to the above function */
-PyUFuncGenericFunction funcs_prism_fz[1] = {&double_prism_fz};
-
-static void double_prism_fzz(char **args, const npy_intp *dimensions,
-                         const npy_intp *steps, void *data)
-{
-    npy_intp i;
-    npy_intp n = dimensions[0];
-    char *inx = args[0], *iny = args[1], *inz = args[2];
-    char *out = args[3];
-    npy_intp inx_step = steps[0], iny_step = steps[1], inz_step = steps[2];
-    npy_intp out_step = steps[3];
-
-    double x, y, z, v, r;
-
-    for (i = 0; i < n; i++) {
-        /* BEGIN main ufunc computation */
-        x = *(double *)inx;
-        y = *(double *)iny;
-        z = *(double *)inz;
-        if (z != 0.0){
-            r = sqrt(x * x + y * y + z * z);
-            v = atan(x * y / (z * r));
-        }else{
-            v = 0.0;
-        }
-        *((double *)out) = v;
-        /* END main ufunc computation */
-
-        inx += inx_step;
-        iny += iny_step;
-        inz += inz_step;
-        out += out_step;
-    }
-}
-
-/* This a pointer to the above function */
-PyUFuncGenericFunction funcs_prism_fzz[1] = {&double_prism_fzz};
-
-static void double_prism_fzx(char **args, const npy_intp *dimensions,
-                         const npy_intp *steps, void *data)
-{
-    npy_intp i;
-    npy_intp n = dimensions[0];
-    char *inx = args[0], *iny = args[1], *inz = args[2];
-    char *out = args[3];
-    npy_intp inx_step = steps[0], iny_step = steps[1], inz_step = steps[2];
-    npy_intp out_step = steps[3];
-
-    double x, y, z, v, r;
-    for (i = 0; i < n; i++) {
-        /* BEGIN main ufunc computation */
-        x = *(double *)inx;
-        y = *(double *)iny;
-        z = *(double *)inz;
-        r = sqrt(x * x + y * y + z * z);
-        v = y + r;
-        if (v == 0.0){
-            if (y < 0){
-                v = log(-2 * y);
-            }
-        }else{
-            v = -log(v);
-        }
-        *((double *)out) = v;
-        /* END main ufunc computation */
-
-        inx += inx_step;
-        iny += iny_step;
-        inz += inz_step;
-        out += out_step;
-    }
-}
-
-/* This a pointer to the above function */
-PyUFuncGenericFunction funcs_prism_fzx[1] = {&double_prism_fzx};
-
-static void double_prism_fzy(char **args, const npy_intp *dimensions,
-                         const npy_intp *steps, void *data)
-{
-    npy_intp i;
-    npy_intp n = dimensions[0];
-    char *inx = args[0], *iny = args[1], *inz = args[2];
-    char *out = args[3];
-    npy_intp inx_step = steps[0], iny_step = steps[1], inz_step = steps[2];
-    npy_intp out_step = steps[3];
-
-    double x, y, z, v, r;
-    for (i = 0; i < n; i++) {
-        /* BEGIN main ufunc computation */
-        x = *(double *)inx;
-        y = *(double *)iny;
-        z = *(double *)inz;
-        r = sqrt(x * x + y * y + z * z);
-        v = x + r;
-        if (v == 0.0){
-            if (x < 0){
-                v = log(-2 * x);
-            }
-        }else{
-            v = -log(v);
-        }
-        *((double *)out) = v;
-        /* END main ufunc computation */
-
-        inx += inx_step;
-        iny += iny_step;
-        inz += inz_step;
-        out += out_step;
-    }
-}
-
-/* This a pointer to the above function */
-PyUFuncGenericFunction funcs_prism_fzy[1] = {&double_prism_fzy};
-
-static void double_prism_fzzz(char **args, const npy_intp *dimensions,
-                         const npy_intp *steps, void *data)
-{
-    npy_intp i;
-    npy_intp n = dimensions[0];
-    char *inx = args[0], *iny = args[1], *inz = args[2];
-    char *out = args[3];
-    npy_intp inx_step = steps[0], iny_step = steps[1], inz_step = steps[2];
-    npy_intp out_step = steps[3];
-
-    double x, y, z, v, r, v1, v2;
-    for (i = 0; i < n; i++) {
-        /* BEGIN main ufunc computation */
-        x = *(double *)inx;
-        y = *(double *)iny;
-        z = *(double *)inz;
-        r = sqrt(x * x + y * y + z * z);
-        v1 = x * x + z * z;
-        v2 = y * y + z * z;
-        v = 0.0;
-        if (v1 != 0.0){
-            v += 1.0/v1;
-        }
-        if (v2 != 0.0){
-            v += 1.0/v2;
-        }
-        if (r != 0.0){
-            v *= x * y / r;
-        }
-
-        *((double *)out) = v;
-        /* END main ufunc computation */
-
-        inx += inx_step;
-        iny += iny_step;
-        inz += inz_step;
-        out += out_step;
-    }
-}
-
-/* This a pointer to the above function */
-PyUFuncGenericFunction funcs_prism_fzzz[1] = {&double_prism_fzzz};
-
-static void double_prism_fxxy(char **args, const npy_intp *dimensions,
-                         const npy_intp *steps, void *data)
-{
-    npy_intp i;
-    npy_intp n = dimensions[0];
-    char *inx = args[0], *iny = args[1], *inz = args[2];
-    char *out = args[3];
-    npy_intp inx_step = steps[0], iny_step = steps[1], inz_step = steps[2];
-    npy_intp out_step = steps[3];
-
-    double x, y, z, v, r;
-    for (i = 0; i < n; i++) {
-        /* BEGIN main ufunc computation */
-        x = *(double *)inx;
-        y = *(double *)iny;
-        z = *(double *)inz;
-        if (x != 0.0){
-            v = x * x + y * y;
-            r = sqrt(x * x + y * y + z * z);
-            v = - x * z / (v * r);
-        }else{
-            v = 0.0;
-        }
-        *((double *)out) = v;
-        /* END main ufunc computation */
-
-        inx += inx_step;
-        iny += iny_step;
-        inz += inz_step;
-        out += out_step;
-    }
-}
-
-/* This a pointer to the above function */
-PyUFuncGenericFunction funcs_prism_fxxy[1] = {&double_prism_fxxy};
-
-static void double_prism_fxxz(char **args, const npy_intp *dimensions,
-                         const npy_intp *steps, void *data)
-{
-    npy_intp i;
-    npy_intp n = dimensions[0];
-    char *inx = args[0], *iny = args[1], *inz = args[2];
-    char *out = args[3];
-    npy_intp inx_step = steps[0], iny_step = steps[1], inz_step = steps[2];
-    npy_intp out_step = steps[3];
-
-    double x, y, z, v, r;
-    for (i = 0; i < n; i++) {
-        /* BEGIN main ufunc computation */
-        x = *(double *)inx;
-        y = *(double *)iny;
-        z = *(double *)inz;
-        if (x != 0.0){
-            v = x * x + z * z;
-            r = sqrt(x * x + y * y + z * z);
-            v = - x * y / (v * r);
-        }else{
-            v = 0.0;
-        }
-        *((double *)out) = v;
-        /* END main ufunc computation */
-
-        inx += inx_step;
-        iny += iny_step;
-        inz += inz_step;
-        out += out_step;
-    }
-}
-
-/* This a pointer to the above function */
-PyUFuncGenericFunction funcs_prism_fxxz[1] = {&double_prism_fxxz};
-
-static void double_prism_fxyz(char **args, const npy_intp *dimensions,
-                         const npy_intp *steps, void *data)
-{
-    npy_intp i;
-    npy_intp n = dimensions[0];
-    char *inx = args[0], *iny = args[1], *inz = args[2];
-    char *out = args[3];
-    npy_intp inx_step = steps[0], iny_step = steps[1], inz_step = steps[2];
-    npy_intp out_step = steps[3];
-
-    double x, y, z, v, r;
-    for (i = 0; i < n; i++) {
-        /* BEGIN main ufunc computation */
-        x = *(double *)inx;
-        y = *(double *)iny;
-        z = *(double *)inz;
-        r = sqrt(x * x + y * y + z * z);
-        if (r != 0.0){
-            v = 1.0/r;
-        }else{
-            v = 0.0;
-        }
-        *((double *)out) = v;
-        /* END main ufunc computation */
-
-        inx += inx_step;
-        iny += iny_step;
-        inz += inz_step;
-        out += out_step;
-    }
-}
-
-/* This a pointer to the above function */
-PyUFuncGenericFunction funcs_prism_fxyz[1] = {&double_prism_fxyz};
-
-/* These are the input and return dtypes of the prism functions.*/
-static char types[4] = {NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE};
-
-static struct PyModuleDef moduledef = {
-    PyModuleDef_HEAD_INIT,
-    "potential_field_prism",
-    NULL,
-    -1,
-    PrismMethods,
-    NULL,
-    NULL,
-    NULL,
-    NULL
-};
-
-static void *data[1] = {NULL,};
-
-PyMODINIT_FUNC PyInit_potential_field_prism(void)
-{
-    PyObject *m, *d;
-    PyObject *prism_f, *prism_fz, *prism_fzz, *prism_fzx, *prism_fzy;
-    PyObject *prism_fzzz, *prism_fxxy, *prism_fxxz, *prism_fxyz;
-    m = PyModule_Create(&moduledef);
-    if (!m) {
-        return NULL;
-    }
-
-    import_array();
-    import_umath();
-
-    prism_f = PyUFunc_FromFuncAndData(funcs_prism_f, data, types, 1, 3, 1,
-                                    PyUFunc_None, "prism_f",
-"Evaluates the indefinite volume integral for the 1/r kernel.\n\n"
-"This is used to evaluate the gravitational potential of dense prisms.\n\n"
-"Parameters\n"
-"----------\n"
-"x, y, z : (...) numpy.ndarray\n"
-"    The nodal locations to evaluate the function at\n\n"
-"Returns\n"
-"-------\n"
-"(...) numpy.ndarray", 0);
-
-    prism_fz = PyUFunc_FromFuncAndData(funcs_prism_fz, data, types, 1, 3, 1,
-                                    PyUFunc_None, "prism_fz",
-"Evaluates the indefinite volume integral for the d/dz * 1/r kernel.\n\n"
-"This is used to evaluate the gravitational field of dense prisms.\n\n"
-"Parameters\n"
-"----------\n"
-"x, y, z : (...) numpy.ndarray\n"
-"    The nodal locations to evaluate the function at\n\n"
-"Returns\n"
-"-------\n"
-"(...) numpy.ndarray\n\n"
-"Notes\n"
-"-----\n"
-"Can be used to compute other components by cycling the inputs", 0);
-
-    prism_fzz = PyUFunc_FromFuncAndData(funcs_prism_fzz, data, types, 1, 3, 1,
-                                    PyUFunc_None, "prism_fzz",
-"Evaluates the indefinite volume integral for the d**2/dz**2 * 1/r kernel.\n\n"
-"This is used to evaluate the gravitational potential of dense prisms.\n\n"
-"Parameters\n"
-"----------\n"
-"x, y, z : (...) numpy.ndarray\n"
-"    The nodal locations to evaluate the function at\n\n"
-"Returns\n"
-"-------\n"
-"  (...) numpy.ndarray\n\n"
-"Notes\n"
-"-----\n"
-"Can be used to compute other components by cycling the inputs", 0);
-
-    prism_fzx = PyUFunc_FromFuncAndData(funcs_prism_fzx, data, types, 1, 3, 1,
-                                    PyUFunc_None, "prism_fzx",
-"Evaluates the indefinite volume integral for the d**2/(dz*dx) * 1/r kernel.\n\n"
-"This is used to evaluate the gravitational potential of dense prisms.\n\n"
-"Parameters\n"
-"----------\n"
-"x, y, z : (...) numpy.ndarray\n"
-"    The nodal locations to evaluate the function at\n\n"
-"Returns\n"
-"-------\n"
-"(...) numpy.ndarray\n\n"
-"Notes\n"
-"-----\n"
-"Can be used to compute other components by cycling the inputs", 0);
-
-    prism_fzy = PyUFunc_FromFuncAndData(funcs_prism_fzy, data, types, 1, 3, 1,
-                                    PyUFunc_None, "prism_fzy",
-"Evaluates the indefinite volume integral for the d**2/(dz*dx) * 1/r kernel.\n\n"
-"This is used to evaluate the gravitational potential of dense prisms.\n\n"
-"Parameters\n"
-"----------\n"
-"x, y, z : (...) numpy.ndarray\n"
-"    The nodal locations to evaluate the function at\n\n"
-"Returns\n"
-"-------\n"
-"(...) numpy.ndarray\n\n"
-"Notes\n"
-"-----\n"
-"Can be used to compute other components by cycling the inputs", 0);
-
-    prism_fzzz = PyUFunc_FromFuncAndData(funcs_prism_fzzz, data, types, 1, 3, 1,
-                                    PyUFunc_None, "prism_fzzz",
-"Evaluates the indefinite volume integral for the d**3/(dz**3) * 1/r kernel.\n\n"
-"This is used to evaluate the magnetic gradient of susceptible prisms.\n\n"
-"Parameters\n"
-"----------\n"
-"x, y, z : (...) numpy.ndarray\n"
-"    The nodal locations to evaluate the function at\n\n"
-"Returns\n"
-"-------\n"
-"(...) numpy.ndarray\n\n"
-"Notes\n"
-"-----\n"
-"Can be used to compute other components by cycling the inputs", 0);
-
-    prism_fxxy = PyUFunc_FromFuncAndData(funcs_prism_fxxy, data, types, 1, 3, 1,
-                                    PyUFunc_None, "prism_fxxy",
-"Evaluates the indefinite volume integral for the d**3/(dx**2 * dy) * 1/r kernel.\n\n"
-"This is used to evaluate the magnetic gradient of susceptible prisms.\n\n"
-"Parameters\n"
-"----------\n"
-"x, y, z : (...) numpy.ndarray\n"
-"    The nodal locations to evaluate the function at\n\n"
-"Returns\n"
-"-------\n"
-"(...) numpy.ndarray\n\n"
-"Notes\n"
-"-----\n"
-"Can be used to compute other components by cycling the inputs", 0);
-
-    prism_fxxz = PyUFunc_FromFuncAndData(funcs_prism_fxxz, data, types, 1, 3, 1,
-                                    PyUFunc_None, "prism_fxxz",
-"Evaluates the indefinite volume integral for the d**3/(dx**2 * dz) * 1/r kernel.\n\n"
-"This is used to evaluate the magnetic gradient of susceptible prisms.\n\n"
-"Parameters\n"
-"----------\n"
-"x, y, z : (...) numpy.ndarray\n"
-"    The nodal locations to evaluate the function at\n\n"
-"Returns\n"
-"-------\n"
-"(...) numpy.ndarray\n\n"
-"Notes\n"
-"-----\n"
-"Can be used to compute other components by cycling the inputs", 0);
-
-    prism_fxyz = PyUFunc_FromFuncAndData(funcs_prism_fxyz, data, types, 1, 3, 1,
-                                    PyUFunc_None, "prism_fxyz",
-"Evaluates the indefinite volume integral for the d**3/(dx * dy * dz) * 1/r kernel.\n\n"
-"This is used to evaluate the magnetic gradient of susceptible prisms.\n\n"
-"Parameters\n"
-"----------\n"
-"x, y, z : (...) numpy.ndarray\n"
-"    The nodal locations to evaluate the function at\n\n"
-"Returns\n"
-"-------\n"
-"(...) numpy.ndarray\n\n"
-"Notes\n"
-"-----\n"
-"Can be used to compute other components by cycling the inputs", 0);
-
-    d = PyModule_GetDict(m);
-
-    PyDict_SetItemString(d, "prism_f", prism_f);
-    Py_DECREF(prism_f);
-
-    PyDict_SetItemString(d, "prism_fz", prism_fz);
-    Py_DECREF(prism_fz);
-
-    PyDict_SetItemString(d, "prism_fzz", prism_fzz);
-    Py_DECREF(prism_fzz);
-
-    PyDict_SetItemString(d, "prism_fzx", prism_fzx);
-    Py_DECREF(prism_fzx);
-
-    PyDict_SetItemString(d, "prism_fzy", prism_fzy);
-    Py_DECREF(prism_fzy);
-
-    PyDict_SetItemString(d, "prism_fzzz", prism_fzzz);
-    Py_DECREF(prism_fzzz);
-
-    PyDict_SetItemString(d, "prism_fxxy", prism_fxxy);
-    Py_DECREF(prism_fxxy);
-
-    PyDict_SetItemString(d, "prism_fxxz", prism_fxxz);
-    Py_DECREF(prism_fxxz);
-
-    PyDict_SetItemString(d, "prism_fxyz", prism_fxyz);
-    Py_DECREF(prism_fxyz);
-
-    return m;
-}
diff --git a/geoana/kernels/_extensions/potential_field_prism.pyx b/geoana/kernels/_extensions/potential_field_prism.pyx
new file mode 100644
index 00000000..07134a39
--- /dev/null
+++ b/geoana/kernels/_extensions/potential_field_prism.pyx
@@ -0,0 +1,298 @@
+cimport cython
+
+from libc.math cimport sqrt, log, atan
+
+@cython.cdivision
+@cython.ufunc
+cdef api double prism_f(double x, double y, double z) nogil:
+    """Evaluates the indefinite volume integral for the 1/r kernel.
+
+    This is used to evaluate the gravitational potential of dense prisms.
+
+    Parameters
+    ----------
+    x, y, z : (...) numpy.ndarray
+        The nodal locations to evaluate the function at
+
+    Returns
+    -------
+    (...) numpy.ndarray
+    """
+    cdef:
+        double v = 0.0
+        double r, temp
+
+    r = sqrt(x * x + y * y + z * z)
+    if x != 0.0 and y != 0.0:
+        temp = z + r
+        if temp > 0.0:
+            v -= x * y * log(temp)
+        v += 0.5 * x * x * atan( y * z / (x * r))
+    if y != 0.0 and z != 0.0:
+        temp = x + r
+        if temp > 0.0:
+            v -= y * z * log(temp)
+        v += 0.5 * y * y * atan(z * x / (y * r))
+    if z != 0.0 and x != 0.0:
+        temp = y + r
+        if temp > 0.0:
+            v -= z * x * log(temp)
+        v += 0.5 * z * z * atan(x * y / (z * r))
+    return v
+
+@cython.cdivision
+@cython.ufunc
+cdef api double prism_fz(double x, double y, double z) nogil:
+    """Evaluates the indefinite volume integral for the d/dz * 1/r kernel.
+
+    This is used to evaluate the gravitational field of dense prisms.
+
+    Parameters
+    ----------
+    x, y, z : (...) numpy.ndarray
+        The nodal locations to evaluate the function at
+
+    Returns
+    -------
+    (...) numpy.ndarray
+
+    Notes
+    -----
+    Can be used to compute other components by cycling the inputs
+    """
+    cdef:
+        double v = 0.0
+        double r, temp
+
+    r = sqrt(x * x + y * y + z * z)
+    if x != 0.0:
+        temp = y + r
+        if temp > 0.0:
+            v += x * log(temp)
+    if y != 0.0:
+        temp = x + r
+        if temp > 0.0:
+            v += y * log(temp)
+    if z != 0.0:
+        v -= z * atan(x * y / (z * r))
+    return v
+
+
+@cython.ufunc
+@cython.cdivision
+cdef api double prism_fzz(double x, double y, double z) nogil:
+    """Evaluates the indefinite volume integral for the d**2/dz**2 * 1/r kernel.
+
+    This is used to evaluate the gravitational gradient of dense prisms.
+
+    Parameters
+    ----------
+    x, y, z : (...) numpy.ndarray
+        The nodal locations to evaluate the function at
+
+    Returns
+    -------
+    (...) numpy.ndarray
+
+    Notes
+    -----
+    Can be used to compute other components by cycling the inputs
+    """
+    cdef:
+        double v = 0.0
+        double r
+
+    if z != 0.0:
+        r = sqrt(x * x + y * y + z * z)
+        v = atan(x * y / (z * r))
+    return v
+
+
+@cython.ufunc
+@cython.cdivision
+cdef api double prism_fzx(double x, double y, double z) nogil:
+    """Evaluates the indefinite volume integral for the d**2/(dz*dx) * 1/r kernel.
+
+    This is used to evaluate the gravitational gradient of dense prisms.
+
+    Parameters
+    ----------
+    x, y, z : (...) numpy.ndarray
+        The nodal locations to evaluate the function at
+
+    Returns
+    -------
+    (...) numpy.ndarray
+
+    Notes
+    -----
+    Can be used to compute other components by cycling the inputs
+    """
+    cdef:
+        double v = 0.0
+        double r
+    r = sqrt(x * x + y * y + z * z)
+    v = y + r
+    if v == 0.0:
+        if y < 0:
+            v = log(-2 * y)
+    else:
+        v = -log(v)
+    return v
+
+
+@cython.ufunc
+@cython.cdivision
+cdef api double prism_fzy(double x, double y, double z) nogil:
+    """Evaluates the indefinite volume integral for the d**2/(dz*dy) * 1/r kernel.
+
+    This is used to evaluate the gravitational gradient of dense prisms.
+
+    Parameters
+    ----------
+    x, y, z : (...) numpy.ndarray
+        The nodal locations to evaluate the function at
+
+    Returns
+    -------
+    (...) numpy.ndarray
+
+    Notes
+    -----
+    Can be used to compute other components by cycling the inputs
+    """
+    cdef:
+        double v = 0.0
+        double r, temp
+    r = sqrt(x * x + y * y + z * z)
+    v = x + r
+    if v == 0.0:
+        if x < 0:
+            v = log(-2 * x)
+    else:
+        v = -log(v)
+    return v
+
+
+@cython.ufunc
+@cython.cdivision
+cdef api double prism_fzzz(double x, double y, double z) nogil:
+    """Evaluates the indefinite volume integral for the d**3/(dz**3) * 1/r kernel.
+
+    This is used to evaluate the magnetic gradient of susceptible prisms.
+
+    Parameters
+    ----------
+    x, y, z : (...) numpy.ndarray
+        The nodal locations to evaluate the function at
+
+    Returns
+    -------
+    (...) numpy.ndarray
+
+    Notes
+    -----
+    Can be used to compute other components by cycling the inputs
+    """
+    cdef:
+        double v = 0.0
+        double r, v1, v2
+    r = sqrt(x * x + y * y + z * z)
+    v1 = x * x + z * z
+    v2 = y * y + z * z
+    if v1 != 0.0:
+        v += 1.0/v1
+    if v2 != 0.0:
+        v += 1.0/v2
+    if r != 0.0:
+        v *= x * y / r
+    return v
+
+
+@cython.ufunc
+@cython.cdivision
+cdef api double prism_fxxy(double x, double y, double z) nogil:
+    """Evaluates the indefinite volume integral for the d**3/(dx**2 * dy) * 1/r kernel.
+
+    This is used to evaluate the magnetic gradient of susceptible prisms.
+
+    Parameters
+    ----------
+    x, y, z : (...) numpy.ndarray
+        The nodal locations to evaluate the function at
+
+    Returns
+    -------
+    (...) numpy.ndarray
+
+    Notes
+    -----
+    Can be used to compute other components by cycling the inputs
+    """
+    cdef:
+        double v = 0.0
+        double r
+    if x != 0.0:
+        v = x * x + y * y
+        r = sqrt(x * x + y * y + z * z)
+        v = - x * z / (v * r)
+    return v
+
+
+@cython.ufunc
+@cython.cdivision
+cdef api double prism_fxxz(double x, double y, double z) nogil:
+    """Evaluates the indefinite volume integral for the d**3/(dx**2 * dz) * 1/r kernel.
+
+    This is used to evaluate the magnetic gradient of susceptible prisms.
+
+    Parameters
+    ----------
+    x, y, z : (...) numpy.ndarray
+        The nodal locations to evaluate the function at
+
+    Returns
+    -------
+    (...) numpy.ndarray
+
+    Notes
+    -----
+    Can be used to compute other components by cycling the inputs
+    """
+    cdef:
+        double v = 0.0
+        double r
+    if x != 0.0:
+        v = x * x + z * z
+        r = sqrt(x * x + y * y + z * z)
+        v = - x * y / (v * r)
+    return v
+
+
+@cython.ufunc
+@cython.cdivision
+cdef api double prism_fxyz(double x, double y, double z) nogil:
+    """Evaluates the indefinite volume integral for the d**3/(dx * dy * dz) * 1/r kernel.
+
+    This is used to evaluate the magnetic gradient of susceptible prisms.
+
+    Parameters
+    ----------
+    x, y, z : (...) numpy.ndarray
+        The nodal locations to evaluate the function at
+
+    Returns
+    -------
+    (...) numpy.ndarray
+
+    Notes
+    -----
+    Can be used to compute other components by cycling the inputs
+    """
+    cdef:
+        double v = 0.0
+        double r
+    r = sqrt(x * x + y * y + z * z)
+    if r != 0.0:
+        v = 1.0/r
+    return v
\ No newline at end of file
diff --git a/geoana/kernels/_extensions/setup.py b/geoana/kernels/_extensions/setup.py
index 674461f0..9209270f 100644
--- a/geoana/kernels/_extensions/setup.py
+++ b/geoana/kernels/_extensions/setup.py
@@ -13,6 +13,8 @@ def configuration(parent_package="", top_path=None):
         from Cython.Build import cythonize
 
         cythonize(os.path.join(base_path, "rTE.pyx"))
+
+        cythonize(os.path.join(base_path, "potential_field_prism.pyx"))
     except ImportError:
         pass
 
@@ -22,6 +24,8 @@ def configuration(parent_package="", top_path=None):
     )
 
     config.add_extension(
-        'potential_field_prism', sources=['potential_field_prism.c']
+        'potential_field_prism',
+        sources=['potential_field_prism.c'],
+        include_dirs=[get_numpy_include_dirs()]
     )
     return config
diff --git a/requirements_dev.txt b/requirements_dev.txt
index 56c838dc..501ab7fb 100644
--- a/requirements_dev.txt
+++ b/requirements_dev.txt
@@ -1,5 +1,4 @@
 sphinx
-sphinx_rtd_theme
 sphinx-gallery
 sphinx-toolbox
 sphinxcontrib-apidoc
diff --git a/tests/test_potential_prism.py b/tests/test_potential_prism.py
index 949bfe7c..ff602bd3 100644
--- a/tests/test_potential_prism.py
+++ b/tests/test_potential_prism.py
@@ -6,6 +6,7 @@
 import geoana.kernels.potential_field_prism as pf
 import geoana.gravity as grav
 from geoana.em.static import MagneticPrism, MagneticDipoleWholeSpace
+from numba import njit
 
 
 class TestCompiledVsNumpy():
@@ -282,3 +283,35 @@ def test_grav_init_and_errors():
 
     with pytest.raises(TypeError):
         prism.rho = 'abc'
+
+
+@pytest.mark.parametrize('function', [
+    pf.prism_f,
+    pf.prism_fz,
+    pf.prism_fzx,
+    pf.prism_fzy,
+    pf.prism_fzz,
+    pf.prism_fzzz,
+    pf.prism_fxxy,
+    pf.prism_fxxz,
+    pf.prism_fxyz
+])
+def test_numba_jitting_nopython(function):
+
+    # create a vectorized jit function of it:
+
+    x = np.random.rand(10)
+    y = np.random.rand(10)
+    z = np.random.rand(10)
+    @njit
+    def jitted_func(x, y, z):
+        n = len(x)
+        out = np.empty_like(x)
+        for i in range(n):
+            out[i] = function(x[i], y[i], z[i])
+        return out
+
+    v1 = jitted_func(x, y, z)
+    v2 = function(x, y, z)
+
+    np.allclose(v1, v2)