diff --git a/README.md b/README.md
index fd9f22ec..4755f041 100644
--- a/README.md
+++ b/README.md
@@ -173,6 +173,8 @@ Contributors (0.9.0):
 - Oscar Benjamin (OB)
 - Daniel Simmons-Marengo (DSM)
 - ForeverHaibara (FH)
+- Jeroen Hanselman (JH)
+- Sam Schiavone (SS)
 
 Changes (0.9.0):
 
@@ -198,6 +200,7 @@ Changes (0.9.0):
 - [gh-359](https://github.com/flintlib/python-flint/pull/359),
   Sort factorisations of all mpoly types. (OB)
 - [gh-374](https://github.com/flintlib/python-flint/pull/374), Fixed a bug in `nmod.__hash__`. (FH)
+- [gh-392](https://github.com/flintlib/python-flint/pull/392), Add interface to `acb_theta_jet` (JH, SS)
 
 0.8.0
 -----
diff --git a/doc/source/acb_theta.rst b/doc/source/acb_theta.rst
index f9ccfdd2..b89cd759 100644
--- a/doc/source/acb_theta.rst
+++ b/doc/source/acb_theta.rst
@@ -3,3 +3,5 @@
 
 .. autofunction :: flint.types.acb_theta.acb_theta
 
+.. autofunction :: flint.types.acb_theta.acb_theta_jets
+
diff --git a/meson.build b/meson.build
index fbabce65..61b9d036 100644
--- a/meson.build
+++ b/meson.build
@@ -73,9 +73,6 @@ endif
 # flint.pc was missing -lflint until Flint 3.1.0
 if flint_dep.version().version_compare('<3.1')
   flint_dep = cc.find_library('flint')
-  have_acb_theta = false
-else
-  have_acb_theta = true
 endif
 
 pyflint_deps = [gmp_dep, mpfr_dep, flint_dep]
diff --git a/setup.py b/setup.py
index 669bcfe4..c9e5e22d 100644
--- a/setup.py
+++ b/setup.py
@@ -119,6 +119,7 @@
     ("flint.types.acb_poly", ["src/flint/types/acb_poly.pyx"]),
     ("flint.types.acb_mat", ["src/flint/types/acb_mat.pyx"]),
     ("flint.types.acb_series", ["src/flint/types/acb_series.pyx"]),
+    ("flint.types.acb_theta", ["src/flint/types/acb_theta.pyx"]),
 
     ("flint.types.dirichlet", ["src/flint/types/dirichlet.pyx"]),
 
diff --git a/src/flint/__init__.py b/src/flint/__init__.py
index e7811272..87155e8f 100644
--- a/src/flint/__init__.py
+++ b/src/flint/__init__.py
@@ -51,22 +51,6 @@
 __version__ = "0.8.0"
 
 
-def _flint_version_at_least(major: int, minor: int) -> bool:
-    version_parts = __FLINT_VERSION__.split(".")
-    if len(version_parts) < 2:
-        return False
-    try:
-        current_major = int(version_parts[0])
-        current_minor = int(version_parts[1])
-    except ValueError:
-        return False
-    return (current_major, current_minor) >= (major, minor)
-
-
-def _has_acb_theta() -> bool:
-    return _flint_version_at_least(3, 1)
-
-
 __all__ = [
     "ctx",
     "fmpz",
diff --git a/src/flint/test/test_acb_mat.py b/src/flint/test/test_acb_mat.py
index d13ab095..9751006e 100644
--- a/src/flint/test/test_acb_mat.py
+++ b/src/flint/test/test_acb_mat.py
@@ -1,10 +1,13 @@
 from __future__ import annotations
 
-import sys
-from unittest.mock import patch
-
-from flint import _has_acb_theta, acb, acb_mat, arb_mat, ctx, fmpq_mat, fmpz_mat
-from flint.test.helpers import is_close_acb, is_close_acb_mat as is_close, is_close_arb_mat, raises
+from flint import acb, acb_mat, arb_mat, ctx, fmpq_mat, fmpz_mat
+from flint.test.helpers import (
+    is_close_acb,
+    is_close_acb_mat as is_close,
+    is_close_arb_mat,
+    raises,
+)
+from flint.types.acb_theta import acb_theta
 
 
 class _DummyMatrix:
@@ -304,7 +307,7 @@ def test_acb_mat_eig_theta_and_helper() -> None:
     assert acb_mat(0, 0).eig() == []
 
     tau = acb_mat([[1j]])
-    if _has_acb_theta():
+    if acb_theta is not None:
         theta_vals = acb_mat.theta(tau, acb_mat([[0]]))
         assert isinstance(theta_vals, acb_mat)
         assert theta_vals.nrows() == 1
@@ -312,9 +315,6 @@ def test_acb_mat_eig_theta_and_helper() -> None:
     else:
         assert raises(lambda: acb_mat.theta(tau, acb_mat([[0]])), NotImplementedError)
 
-    with patch.dict(sys.modules, {"flint.types.acb_theta": None}):
-        assert raises(lambda: acb_mat.theta(tau, acb_mat([[0]])), NotImplementedError)
-
     assert is_close(a, [[1, 0], [0, 2]]) is True
     assert is_close(a, acb_mat([[1, 0], [0, 2]])) is True
     assert is_close(a, [[1, 0, 0], [0, 2, 0]]) is False
diff --git a/src/flint/test/test_acb_theta.py b/src/flint/test/test_acb_theta.py
index 4ea6237a..4446002c 100644
--- a/src/flint/test/test_acb_theta.py
+++ b/src/flint/test/test_acb_theta.py
@@ -1,37 +1,35 @@
 from __future__ import annotations
 
-from flint import _has_acb_theta, acb, acb_mat
+from flint import acb, acb_mat
 from flint.test.helpers import is_close_acb_mat as is_close, raises
+from flint.types.acb_theta import acb_theta, acb_theta_jets
 
 
 def test_acb_theta_basic() -> None:
-    if not _has_acb_theta():
+    if acb_theta is None:
         return
-
-    from flint.types.acb_theta import acb_theta
+    theta = acb_theta
 
     z = acb(1 + 1j)
     tau = acb(1.25 + 3j)
     zmat = acb_mat([[z]])
     taumat = acb_mat([[tau]])
 
-    direct = acb_theta(zmat, taumat)
+    direct = theta(zmat, taumat)
     via_method = taumat.theta(zmat)
     assert is_close(direct, via_method, tol=1e-12, rel_tol=1e-12, max_width=1e-12)
     assert direct.nrows() == 1
     assert direct.ncols() == 4
 
-    squared = acb_theta(zmat, taumat, square=True)
+    squared = theta(zmat, taumat, square=True)
     assert squared.nrows() == 1
     assert squared.ncols() == 4
 
 
 def test_acb_theta_shape_assertions() -> None:
-    if not _has_acb_theta():
+    if acb_theta is None:
         return
 
-    from flint.types.acb_theta import acb_theta
-
     z = acb_mat([[0]])
     tau_bad = acb_mat([[1j, 0]])
     assert raises(lambda: acb_theta(z, tau_bad), AssertionError)
@@ -45,3 +43,20 @@ def test_acb_theta_shape_assertions() -> None:
 
     assert raises(lambda: acb_theta(object(), tau), TypeError)  # type: ignore[arg-type]
     assert raises(lambda: acb_theta(z, object()), TypeError)  # type: ignore[arg-type]
+
+
+def test_acb_theta_jets_basic() -> None:
+    if acb_theta_jets is None:
+        return
+
+    z = acb(1 + 1j)
+    tau = acb(1.25 + 3j)
+    zmat = acb_mat([[z]])
+    taumat = acb_mat([[tau]])
+    ord = 2
+
+    direct = acb_theta_jets(zmat, taumat, ord)
+    via_method = taumat.theta_jets(zmat, ord)
+    assert is_close(direct, via_method, tol=1e-12, rel_tol=1e-12, max_width=1e-12)
+    assert direct.nrows() == 4
+    assert direct.ncols() == 3
diff --git a/src/flint/test/test_docstrings.py b/src/flint/test/test_docstrings.py
index 66de7210..8c2b8c56 100644
--- a/src/flint/test/test_docstrings.py
+++ b/src/flint/test/test_docstrings.py
@@ -10,7 +10,7 @@
 test_flint_at_least = {
     "flint.types._gr.gr_ctx.gens": 30100,
     "flint.types._gr.gr_ctx.neg": 30100,
-    "flint.types.acb_theta.acb_theta": 30300,
+    "flint.types.acb_theta.acb_theta": 30100,
 }
 
 
diff --git a/src/flint/types/acb_mat.pyi b/src/flint/types/acb_mat.pyi
index 87a8736b..75a4e523 100644
--- a/src/flint/types/acb_mat.pyi
+++ b/src/flint/types/acb_mat.pyi
@@ -96,6 +96,7 @@ class acb_mat(flint_mat[acb]):
     ) -> Any: ...
 
     def theta(self, z: acb_mat, square: bool = False) -> acb_mat: ...
+    def theta_jets(self, z: acb_mat, ord: int) -> acb_mat: ...
 
     def str(self, n: int = 0, radius: bool = True, more: bool = False, condense: int = 0) -> _str: ...
     def repr(self) -> _str: ...
diff --git a/src/flint/types/acb_mat.pyx b/src/flint/types/acb_mat.pyx
index ae65fe71..dec64eb7 100644
--- a/src/flint/types/acb_mat.pyx
+++ b/src/flint/types/acb_mat.pyx
@@ -5,9 +5,8 @@ from flint.types.arb_mat cimport arb_mat
 from flint.types.fmpz_mat cimport fmpz_mat
 from flint.types.fmpq_mat cimport fmpq_mat
 from flint.types.arb cimport arb
-from flint.types.acb cimport acb
+from flint.types.acb cimport acb, any_as_acb
 from flint.types.acb_poly cimport acb_poly
-from flint.types.acb cimport any_as_acb
 from flint.types.fmpz cimport fmpz
 from flint.types.fmpq cimport fmpq
 
@@ -38,6 +37,8 @@ from flint.flintlib.types.acb cimport (
 
 cimport cython
 
+from flint.types.acb_theta import acb_theta, acb_theta_jets
+
 cdef acb_mat_coerce_operands(x, y):
     if isinstance(y, (fmpz_mat, fmpq_mat, arb_mat)):
         return x, acb_mat(y)
@@ -834,8 +835,25 @@ cdef class acb_mat(flint_mat):
         for the ordering of the theta characteristics.
 
         """
-        try:
-            from .acb_theta import acb_theta
-        except ImportError:
+        if acb_theta is None:
             raise NotImplementedError("acb_mat.theta needs Flint >= 3.1.0")
         return acb_theta(z, tau, square=square)
+
+    def theta_jets(tau, z, ord):
+        r"""
+        Computes Taylor approximation for the vector-valued Riemann theta function
+        `(\theta_{a,b}(z, \tau) : a, b \in \{0,1\}^{g})` or its squares,
+        where `\tau` is given by ``self``.
+
+        This is a wrapper for :func:`.acb_theta_jet.acb_theta_jet`; see the
+        documentation for that method for details and examples.
+        This follows the same conventions of the C-function
+        `acb_theta_jet <https://flintlib.org/doc/acb_theta.html#c.acb_theta_jet>`_
+        for the ordering of the theta characteristics.
+        The result is an ``acb_mat`` with one row for each theta characteristic and
+        one column for each jet coefficient.
+
+        """
+        if acb_theta_jets is None:
+            raise NotImplementedError("acb_mat.theta_jets needs Flint >= 3.3.0")
+        return acb_theta_jets(z, tau, ord)
diff --git a/src/flint/types/acb_theta.pyi b/src/flint/types/acb_theta.pyi
index 37f1e1e3..e5bc0239 100644
--- a/src/flint/types/acb_theta.pyi
+++ b/src/flint/types/acb_theta.pyi
@@ -1,6 +1,18 @@
 from __future__ import annotations
 
+from typing import Protocol
+
 from flint.types.acb_mat import acb_mat
 
 
-def acb_theta(z: acb_mat, tau: acb_mat, square: bool | int = False) -> acb_mat: ...
+class _AcbThetaFunc(Protocol):
+    def __call__(self, z: acb_mat, tau: acb_mat, square: bool | int = False) -> acb_mat: ...
+
+
+class _AcbThetaJetsFunc(Protocol):
+    def __call__(self, z: acb_mat, tau: acb_mat, ord: int) -> acb_mat: ...
+
+
+acb_theta: _AcbThetaFunc | None
+
+acb_theta_jets: _AcbThetaJetsFunc | None
diff --git a/src/flint/types/acb_theta.pyx b/src/flint/types/acb_theta.pyx
index 80ed5df1..240faaf6 100644
--- a/src/flint/types/acb_theta.pyx
+++ b/src/flint/types/acb_theta.pyx
@@ -1,12 +1,68 @@
 from flint.flint_base.flint_context cimport getprec
+from flint.flint_base.flint_base import FLINT_RELEASE
 from flint.types.acb cimport acb
 from flint.types.acb_mat cimport acb_mat
 from flint.flintlib.functions.acb cimport *
 from flint.flintlib.types.acb cimport (
+    acb_mat_t,
     acb_mat_entry,
+    acb_ptr,
+    acb_srcptr,
 )
 from flint.flintlib.functions.acb_mat cimport *
-from flint.flintlib.functions.acb_theta cimport *
+from flint.flintlib.types.flint cimport slong, ulong
+
+
+cdef extern from *:
+    """
+    #include "flint/flint.h"
+    #include "flint/acb.h"
+    #include "flint/acb_mat.h"
+
+    #if __FLINT_RELEASE >= 30100 /* Flint 3.1.0 or later */
+    #include "flint/acb_theta.h"
+    static inline void
+    compat_acb_theta_all(acb_ptr th, acb_srcptr z, const acb_mat_t tau, int sqr, slong prec)
+    {
+        acb_theta_all(th, z, tau, sqr, prec);
+    }
+    #else
+    static inline void
+    compat_acb_theta_all(acb_ptr th, acb_srcptr z, const acb_mat_t tau, int sqr, slong prec)
+    {
+    }
+    #endif
+
+    #if __FLINT_RELEASE >= 30300 /* Flint 3.3.0 or later */
+    static inline slong
+    compat_acb_theta_jet_nb(slong ord, slong g)
+    {
+        return acb_theta_jet_nb(g, ord);
+    }
+
+    static inline void
+    compat_acb_theta_jet(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau,
+            slong ord, ulong ab, int all, int sqr, slong prec)
+    {
+        acb_theta_jet(th, zs, nb, tau, ord, ab, all, sqr, prec);
+    }
+    #else
+    static inline slong
+    compat_acb_theta_jet_nb(slong ord, slong g)
+    {
+        return 0;
+    }
+
+    static inline void
+    compat_acb_theta_jet(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau,
+            slong ord, ulong ab, int all, int sqr, slong prec)
+    {
+    }
+    #endif
+    """
+    void compat_acb_theta_all(acb_ptr th, acb_srcptr z, const acb_mat_t tau, int sqr, slong prec)
+    slong compat_acb_theta_jet_nb(slong ord, slong g)
+    void compat_acb_theta_jet(acb_ptr th, acb_srcptr zs, slong nb, const acb_mat_t tau, slong ord, ulong ab, int all, int sqr, slong prec)
 
 
 def acb_theta(acb_mat z, acb_mat tau, ulong square=False):
@@ -84,7 +140,7 @@ def acb_theta(acb_mat z, acb_mat tau, ulong square=False):
     cdef slong nb = 1 << (2 * g)
     cdef acb_ptr theta = _acb_vec_init(nb)
 
-    acb_theta_all(theta, zvec, tau.val, square, getprec())
+    compat_acb_theta_all(theta, zvec, tau.val, square, getprec())
     _acb_vec_clear(zvec, g)
     # copy the output
     res = []
@@ -95,3 +151,85 @@ def acb_theta(acb_mat z, acb_mat tau, ulong square=False):
         res.append(r)
     _acb_vec_clear(theta, nb)
     return acb_mat([res])
+
+
+def acb_theta_jets(acb_mat z, acb_mat tau, slong ord):
+    r"""
+    Computes the coefficients of the Taylor expansion of the vector valued Riemann
+    theta function `(\theta_{a,b}(z, \tau) : a, b \in \{0,1\}^{g})` or its squares.
+
+    This is a wrapper for the C-function
+    `acb_theta_jet <https://flintlib.org/doc/acb_theta.html#c.acb_theta_jet>`_
+    and it follows the same conventions for the ordering of the theta characteristics.
+    The result is an ``acb_mat`` with one row for each theta characteristic and
+    one column for each jet coefficient.
+
+    This should be used via the method :meth:`.acb_mat.theta_jets`, explicitly ``tau.theta_jets(z, ord)``.
+
+        >>> from flint import acb, acb_mat, showgood, ctx
+        >>> z = acb(1+1j); tau = acb(1.25+3j)
+        >>> acb_mat([[tau]]).theta_jets(acb_mat([[z]]),2)  # doctest: +SKIP
+        [[0.969443038779670 +/- 5.67e-16] + [-0.0305569612081680 +/- 5.13e-17]j,
+        [-0.191993710594950 +/- 4.89e-16] + [0.191993710747776 +/- 7.42e-16]j,
+        [0.60317023860834 +/- 2.93e-15] + [0.60317023764810 +/- 4.86e-15]j]
+        [[1.03055696119601 +/- 3.89e-15] + [0.0305569612081680 +/- 5.13e-17]j,
+        [0.191993710594950 +/- 4.89e-16] + [-0.191993710442123 +/- 4.62e-16]j,
+        [-0.60317023668787 +/- 5.90e-15] + [-0.60317023764810 +/- 4.86e-15]j]
+        [[-1.22079026757697 +/- 4.36e-15] + [-1.82705551679115 +/- 5.17e-15]j,
+        [-5.71849316258739 +/- 7.02e-15] + [3.82088827346268 +/- 5.75e-15]j,
+        [6.0241074288587 +/- 3.88e-14] + [9.0163253443780 +/- 2.05e-14]j]
+        [[-1.82023591012499 +/- 2.67e-15] + [1.21625195015448 +/- 4.14e-15]j,
+        [3.8353056542516 +/- 4.99e-14] + [5.73981078971270 +/- 6.74e-15]j,
+        [8.9823364151977 +/- 2.44e-14] + [-6.0022138700195 +/- 3.72e-14]j]
+    """
+    g = tau.nrows()
+    if g == 0:
+        return acb_mat(0, 0)
+
+    # Calculate the length of the jet for one characteristic
+    # This is the number of multi-indices (alpha) such that |alpha| < ord
+    cdef slong nj = compat_acb_theta_jet_nb(ord, g)
+
+    # Total number of characteristics
+    cdef slong nb = 1 << (2 * g)
+
+    # Total number of acb elements to allocate
+    # FLINT stores nj coefficients for each of the nb characteristics
+    cdef slong total_size = nb * nj
+
+    cdef acb_ptr zvec = _acb_vec_init(g)
+    cdef slong i, j
+    for i in range(g):
+        acb_set(zvec + i, acb_mat_entry(z.val, i, 0))
+
+    # Parameters for characteristics
+    cdef slong nb_in = 1  # Number of input z vectors
+    cdef ulong ab = 0     # Base characteristic
+    cdef ulong all = True  # Compute all 2^2g characteristics
+    cdef ulong square = False  # Don't compute the squares of the thetas.
+
+    # Initialize the output buffer
+    cdef acb_ptr theta = _acb_vec_init(total_size)
+
+    # Call the FLINT C function
+    # Note: Computes all partial derivatives up to total order 'ord'
+    compat_acb_theta_jet(theta, zvec, nb_in, tau.val, ord, ab, all, square, getprec())
+
+    # Copy the output into a structured format
+    res_mat = acb_mat(nb, nj)
+    for i in range(nb):
+        for j in range(nj):
+            acb_set(acb_mat_entry(res_mat.val, i, j), theta + (i * nj + j))
+
+    # Cleanup
+    _acb_vec_clear(zvec, g)
+    _acb_vec_clear(theta, total_size)
+
+    return res_mat
+
+
+if FLINT_RELEASE < 30100:
+    acb_theta = None
+
+if FLINT_RELEASE < 30300:
+    acb_theta_jets = None
diff --git a/src/flint/types/meson.build b/src/flint/types/meson.build
index b894ffd9..1e527091 100644
--- a/src/flint/types/meson.build
+++ b/src/flint/types/meson.build
@@ -83,6 +83,7 @@ exts = [
   'arb',
   'arb_poly',
   'arb_mat',
+  'acb_theta',
   'arb_series',
 
   'acb',
@@ -95,10 +96,6 @@ exts = [
   '_gr',
 ]
 
-if have_acb_theta
-  exts += ['acb_theta']
-endif
-
 py.install_sources(
   pyfiles,
   pure: false,