Repository: rodluger/starry (language: python)

Issue:
Fixes gradients
Fixes a bug with gradients and adds a test for the oblate gradients. Should work with pym3 now!


Discussion / hints:

Fix patch:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index ec99e69c1..51d1094ed 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -14,11 +14,10 @@ repos:
   rev: 19.3b0
   hooks:
   - id: black
-    language_version: python3.7
     exclude: lib
 
 - repo: https://github.com/dfm/black_nbconvert
-  rev: stable
+  rev: v0.3.0
   hooks:
   - id: black_nbconvert
     exclude: lib
diff --git a/docs/MapFactory.rst b/docs/MapFactory.rst
index 26cb847eb..f536fe1da 100644
--- a/docs/MapFactory.rst
+++ b/docs/MapFactory.rst
@@ -15,3 +15,4 @@ arguments passed to it.
     LimbDarkenedMap
     RadialVelocityMap
     ReflectedLightMap
+    OblateMap
diff --git a/docs/OblateMap.rst b/docs/OblateMap.rst
new file mode 100644
index 000000000..da862ee26
--- /dev/null
+++ b/docs/OblateMap.rst
@@ -0,0 +1,7 @@
+Oblate Maps
+===========
+
+.. autoclass:: starry._OblateMap()
+    :noindex:
+    :members:
+    :inherited-members:
diff --git a/docs/hacks.py b/docs/hacks.py
index 15a8edcdf..b2401d4ad 100644
--- a/docs/hacks.py
+++ b/docs/hacks.py
@@ -64,17 +64,22 @@ class _RadialVelocityMap(CustomBase, starry.maps.RVBase, starry.maps.MapBase):
     __doc__ = starry.maps.RVBase.__doc__
 
 
+class _OblateMap(CustomBase, starry.maps.OblateBase, starry.maps.MapBase):
+    __doc__ = starry.maps.OblateBase.__doc__
+
+
 starry._SphericalHarmonicMap = _SphericalHarmonicMap
 starry._LimbDarkenedMap = _LimbDarkenedMap
 starry._ReflectedLightMap = _ReflectedLightMap
 starry._RadialVelocityMap = _RadialVelocityMap
-
+starry._OblateMap = _OblateMap
 
 replace = {
     "_SphericalHarmonicMap": "Map",
     "_LimbDarkenedMap": "Map",
     "_ReflectedLightMap": "Map",
     "_RadialVelocityMap": "Map",
+    "_OblateMap": "Map",
 }
 
 
diff --git a/notebooks/todo/EasySampling.ipynb b/notebooks/todo/EasySampling.ipynb
new file mode 100644
index 000000000..d384e32e2
--- /dev/null
+++ b/notebooks/todo/EasySampling.ipynb
@@ -0,0 +1,203 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%run ../notebook_setup.py"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import starry\n",
+    "import pymc3 as pm\n",
+    "import pymc3_ext as pmx"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Generate"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "map = starry.Map(oblate=True, lazy=False)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "truths = {\n",
+    "    \"f\": 0.3,\n",
+    "    \"omega\": 0.5,\n",
+    "    \"inc\": 60,\n",
+    "    \"obl\": 30,\n",
+    "    \"b0\": 0.5,\n",
+    "    \"ro\": 0.1\n",
+    "}\n",
+    "\n",
+    "ferr = 3e-4"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "map.f = truths[\"f\"]\n",
+    "map.omega = truths[\"omega\"]\n",
+    "map.inc = truths[\"inc\"]\n",
+    "map.obl = truths[\"obl\"]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "xo = np.linspace(-1, 1, 300)\n",
+    "yo = truths[\"b0\"] * np.ones_like(xo)\n",
+    "ro = truths[\"ro\"]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(1)\n",
+    "ax.plot(xo, yo, \"k-\")\n",
+    "ax.plot(xo, yo - ro, \"k--\")\n",
+    "ax.plot(xo, yo + ro, \"k--\")\n",
+    "map.show(ax=ax)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "flux0 = map.flux(xo=xo, yo=yo, ro=ro)\n",
+    "flux = flux0 + ferr * np.random.randn(len(flux0))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.plot(flux0);\n",
+    "plt.plot(flux, \"k.\", ms=3, alpha=0.75);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Infer"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "with pm.Model() as model:\n",
+    "    f = pm.Uniform(\"f\", lower=0.0, upper=0.9, testval=truths[\"f\"])\n",
+    "    omega = pm.Uniform(\"omega\", lower=0.0, upper=0.9, testval=truths[\"omega\"])\n",
+    "    inc = pm.Uniform(\"inc\", lower=0.0, upper=90.0, testval=truths[\"inc\"])\n",
+    "    obl = pm.Uniform(\"obl\", lower=-90, upper=90.0, testval=truths[\"obl\"])\n",
+    "    b0 = pm.Uniform(\"b0\", lower=-1.0, upper=1.0, testval=truths[\"b0\"])\n",
+    "    ro = pm.Uniform(\"ro\", lower=1e-5, upper=0.5, testval=truths[\"ro\"])\n",
+    "    map = starry.Map(oblate=True)\n",
+    "    map.f = f\n",
+    "    map.omega = omega\n",
+    "    map.inc = inc\n",
+    "    map.obl = obl\n",
+    "    yo = b0 * np.ones_like(xo)\n",
+    "    flux_model = map.flux(xo=xo, yo=yo, ro=ro)\n",
+    "    pm.Normal(\"obs\", mu=flux_model, sd=ferr, observed=flux)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "np.random.seed(42)\n",
+    "with model:\n",
+    "    trace = pmx.sample(\n",
+    "        tune=1000,\n",
+    "        draws=1000,\n",
+    "        chains=2,\n",
+    "        cores=1,\n",
+    "        target_accept=0.9,\n",
+    "        return_inferencedata=True,\n",
+    "    )"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pm.summary(trace)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pm.plot_trace(trace);"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/notebooks/todo/SmoothGradients.ipynb b/notebooks/todo/SmoothGradients.ipynb
new file mode 100644
index 000000000..eef076f9f
--- /dev/null
+++ b/notebooks/todo/SmoothGradients.ipynb
@@ -0,0 +1,184 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import starry\n",
+    "import theano\n",
+    "import theano.tensor as tt"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "starry.config.quiet = True"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "map = starry.Map(oblate=True, lazy=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "map.f = 0.3\n",
+    "map.omega = 0.5\n",
+    "map.inc = 60\n",
+    "map.obl = 30\n",
+    "xo = np.linspace(-1, 1, 300)\n",
+    "yo = 0.5 * np.ones_like(xo)\n",
+    "ro = 0.1"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(1)\n",
+    "ax.plot(xo, yo, \"k-\")\n",
+    "ax.plot(xo, yo - ro, \"k--\")\n",
+    "ax.plot(xo, yo + ro, \"k--\")\n",
+    "map.show(ax=ax)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def plot_grad(wrt, res=1000):\n",
+    "    var = tt.dscalar()\n",
+    "    n = tt.iscalar()\n",
+    "    map.f = 0.5\n",
+    "    map.omega = 0.5\n",
+    "    map.inc = 60\n",
+    "    map.obl = 30\n",
+    "    xo = np.linspace(-1, 1, res)\n",
+    "    yo = 0.5 * np.ones_like(xo)\n",
+    "    ro = 0.1\n",
+    "    if wrt == \"f\":\n",
+    "        map.f = var\n",
+    "    elif wrt == \"omega\":\n",
+    "        map.omega = var\n",
+    "    elif wrt == \"inc\":\n",
+    "        map.inc = var\n",
+    "    elif wrt == \"obl\":\n",
+    "        map.obl = var\n",
+    "    elif wrt == \"b0\":\n",
+    "        yo = var * np.ones_like(xo)\n",
+    "    elif wrt == \"ro\":\n",
+    "        ro = var\n",
+    "    else:\n",
+    "        raise ValueError(wrt)\n",
+    "    grad = theano.function([var, n], tt.grad(map.flux(xo=xo, yo=yo, ro=ro)[n], var))\n",
+    "    arr = np.zeros(res)\n",
+    "    for n in range(res):\n",
+    "        arr[n] = grad(0.2, n)\n",
+    "    plt.plot(arr)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plot_grad(\"f\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plot_grad(\"omega\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plot_grad(\"inc\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plot_grad(\"obl\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plot_grad(\"b0\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plot_grad(\"ro\")"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/setup.cfg b/setup.cfg
index 73621c3ee..337eaca2e 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -5,3 +5,6 @@ setup_requires =
 
 [metadata]
 license_files = LICENSE
+
+[tool:pytest]
+norecursedirs = starry/_core/ops/lib/include/oblate/*
diff --git a/setup.py b/setup.py
index a1161ec57..173b24d9c 100644
--- a/setup.py
+++ b/setup.py
@@ -12,22 +12,10 @@
 on_rtd = os.getenv("READTHEDOCS") == "True"
 
 # Custom compiler flags
-macros = dict(
-    STARRY_NDIGITS=16,
-    STARRY_ELLIP_MAX_ITER=200,
-    STARRY_MAX_LMAX=50,
-    STARRY_BCUT=1.0e-3,
-    STARRY_MN_MAX_ITER=100,
-    STARRY_IJ_MAX_ITER=200,
-    STARRY_REFINE_J_AT=25,
-    STARRY_USE_INCOMPLETE_INTEGRALS=0,
-    STARRY_QUAD_POINTS=100,
-    STARRY_EL2_MAX_ITER=100,
-)
-
-# Override with user values
-for key, value in macros.items():
-    macros[key] = os.getenv(key, value)
+macros = {}
+for key in os.environ.keys():
+    if key.startswith("STARRY_"):
+        macros[key] = os.getenv(key)
 
 # Compiler optimization flag -O
 optimize = int(os.getenv("STARRY_O", 2))
@@ -47,6 +35,10 @@
     macros["STARRY_O"] = 0
     macros["STARRY_DEBUG"] = 1
 
+# Unit tests for C++ functions?
+if bool(int(os.getenv("STARRY_UNIT_TESTS", 0))):
+    macros["STARRY_UNIT_TESTS"] = 1
+
 # Numerical override at high l?
 if bool(int(os.getenv("STARRY_KL_NUMERICAL", 0))):
     macros["STARRY_KL_NUMERICAL"] = 1
@@ -170,7 +162,7 @@ def build_extensions(self):
         build_ext.build_extensions(self)
 
 
-if int(macros["STARRY_NDIGITS"]) > 16:
+if int(os.getenv("STARRY_NDIGITS", 16)) > 16:
     include_dirs = ["starry/_core/ops/lib/vendor/boost_1_66_0"]
 else:
     include_dirs = []
@@ -236,6 +228,7 @@ def build_extensions(self):
             "scikit-image",
             "pillow",
             "tqdm",
+            "mpmath",
         ],
         "docs": [
             "sphinx>=1.7.5",
diff --git a/starry/_core/core.py b/starry/_core/core.py
index 45cee6fd2..3960cb31a 100644
--- a/starry/_core/core.py
+++ b/starry/_core/core.py
@@ -6,6 +6,7 @@
     sTOp,
     rTReflectedOp,
     sTReflectedOp,
+    sTOblateOp,
     dotROp,
     tensordotRzOp,
     FOp,
@@ -22,6 +23,8 @@
 from .utils import logger, autocompile, is_tensor, clear_cache
 from .math import lazy_math as math
 from scipy.special import legendre as LegendreP
+from scipy.sparse import csc_matrix
+from scipy.sparse.linalg import inv as sparse_inv
 import numpy as np
 from astropy import units
 import os
@@ -35,13 +38,22 @@
     from .. import _c_ops
 
 
-__all__ = ["OpsYlm", "OpsLD", "OpsReflected", "OpsRV", "OpsSystem"]
+__all__ = [
+    "OpsYlm",
+    "OpsLD",
+    "OpsReflected",
+    "OpsRV",
+    "OpsOblate",
+    "OpsSystem",
+]
 
 
 class OpsYlm(object):
     """Class housing Theano operations for spherical harmonics maps."""
 
-    def __init__(self, ydeg, udeg, fdeg, nw, reflected=False, **kwargs):
+    def __init__(
+        self, ydeg, udeg, fdeg, nw, reflected=False, oblate=False, **kwargs
+    ):
         # Ingest kwargs
         self.ydeg = ydeg
         self.udeg = udeg
@@ -49,6 +61,7 @@ def __init__(self, ydeg, udeg, fdeg, nw, reflected=False, **kwargs):
         self.deg = ydeg + udeg + fdeg
         self.filter = (fdeg > 0) or (udeg > 0)
         self.nw = nw
+        self._oblate = oblate
         self._reflected = reflected
         self.Ny = (self.ydeg + 1) ** 2
 
@@ -126,11 +139,17 @@ def sT(self, b, r):
 
     @autocompile
     def tensordotRz(self, matrix, theta):
-        return self._tensordotRz(matrix, theta)
+        if self.ydeg == 0:
+            return matrix
+        else:
+            return self._tensordotRz(matrix, theta)
 
     @autocompile
     def dotR(self, matrix, ux, uy, uz, theta):
-        return self._dotR(matrix, ux, uy, uz, theta)
+        if self.ydeg == 0:
+            return matrix
+        else:
+            return self._dotR(matrix, ux, uy, uz, theta)
 
     @autocompile
     def F(self, u, f):
@@ -427,10 +446,6 @@ def right_project(self, M, inc, obl, theta):
         that transforms a spherical harmonic coefficient vector in the
         input frame to a vector in the observer's frame.
         """
-        # Trivial case
-        if self.ydeg == 0:
-            return M
-
         # Rotate to the sky frame
         # TODO: Do this in a single compound rotation
         M = self.dotR(
@@ -487,10 +502,6 @@ def left_project(self, M, inc, obl, theta):
         input frame to a vector in the observer's frame.
 
         """
-        # Trivial case
-        if self.ydeg == 0:
-            return M
-
         # Note that here we are using the fact that R . M = (M^T . R^T)^T
         MT = tt.transpose(M)
 
@@ -667,10 +678,9 @@ def spot(self, contrast, radius, lat, lon):
 class OpsLD(object):
     """Class housing Theano operations for limb-darkened maps."""
 
-    def __init__(self, ydeg, udeg, fdeg, nw, reflected=False, **kwargs):
+    def __init__(self, ydeg, udeg, fdeg, nw, **kwargs):
         # Sanity checks
         assert ydeg == fdeg == 0
-        assert reflected is False
 
         # Ingest kwargs
         self.udeg = udeg
@@ -1357,6 +1367,340 @@ def compute_illumination(self, xyz, xs, ys, zs, Rs, sigr, on94_exact):
             return tt.sum(I, axis=2) / self.source_npts
 
 
+class OpsOblate(OpsYlm):
+    def __init__(self, *args, **kwargs):
+        super(OpsOblate, self).__init__(*args, oblate=True, **kwargs)
+        self._sT = sTOblateOp(self._c_ops.sTOblate, self._c_ops.N)
+
+        # Spherical harmonic transform operator
+        smoothing = 0.0
+        eps4 = 1e-9
+        npts = 4 * (self.fdeg + 1) ** 2
+        z = np.linspace(-1, 1, npts)
+        B = np.hstack(
+            [
+                np.sqrt(2 * l + 1) * LegendreP(l)(z).reshape(-1, 1)
+                for l in range(self.fdeg + 1)
+            ]
+        )
+        A = np.linalg.solve(B.T @ B + eps4 * np.eye(self.fdeg + 1), B.T)
+        l = np.arange(self.fdeg + 1)
+        idx = l * (l + 1)
+        S = np.exp(-0.5 * idx * smoothing ** 2)
+        self.SHT = tt.as_tensor_variable(S[:, None] * A)
+        self.z = tt.as_tensor_variable(z)
+        self.idx = idx
+
+        # We need to instantiate custom versions of the C++ ops
+        # so we can apply and manipulate the gravity darkening filter
+        # Several steps in `render` and `X` require operators with
+        # different dimensions. Some day it would make sense to
+        # handle this logic on the C++ side so we don't have to
+        # instantiate so many Ops objects.
+        if self.fdeg == 0 and self.udeg == 0:
+            ops_y_0_f = self._c_ops
+            ops_yf_0_0 = self._c_ops
+            ops_yu_0_0 = self._c_ops
+            ops_f_0_0 = self._c_ops
+            ops_yf_u_0 = self._c_ops
+            ops_yuf_0_0 = self._c_ops
+        elif self.udeg == 0:
+            ops_y_0_f = self._c_ops
+            ops_yf_0_0 = _c_ops.Ops(self.ydeg + self.fdeg, 0, 0)
+            ops_yu_0_0 = _c_ops.Ops(self.ydeg, 0, 0)
+            ops_f_0_0 = _c_ops.Ops(self.fdeg, 0, 0)
+            ops_yf_u_0 = ops_yf_0_0
+            ops_yuf_0_0 = ops_yf_0_0
+        else:
+            ops_y_0_f = _c_ops.Ops(self.ydeg, 0, self.fdeg)
+            ops_yf_0_0 = _c_ops.Ops(self.ydeg + self.fdeg, 0, 0)
+            ops_yu_0_0 = _c_ops.Ops(self.ydeg + self.udeg, 0, 0)
+            ops_f_0_0 = _c_ops.Ops(self.fdeg, 0, 0)
+            ops_yf_u_0 = _c_ops.Ops(self.ydeg + self.fdeg, self.udeg, 0)
+            ops_yuf_0_0 = _c_ops.Ops(self.ydeg + self.udeg + self.fdeg, 0, 0)
+
+        # Filter operator (gravity darkening)
+        self.Fg = FOp(
+            ops_y_0_f.F, (self.ydeg + self.fdeg + 1) ** 2, (self.ydeg + 1) ** 2
+        )
+
+        # Change of basis matrices
+        if self.ydeg == 0 and self.fdeg == 0:
+            self.A1Inv_Nyf_x_Nyf = csc_matrix([[np.pi]])
+        else:
+            self.A1Inv_Nyf_x_Nyf = ts.as_sparse_variable(
+                sparse_inv(ops_y_0_f.A1Big)
+            )
+        self.A1_Nyf_x_Nyf = ts.as_sparse_variable(ops_y_0_f.A1Big)
+
+        # Rotation operator (gravity-darkened map)
+        self._dotR = dotROp(ops_yf_0_0.dotR)
+
+        # Rotation operator (limb-darkened map)
+        self._dotR_Nyu_x_Nyu = dotROp(ops_yu_0_0.dotR)
+
+        # Rotation operator (gravity darkening filter)
+        self._dotR_Nf_x_Nf = dotROp(ops_f_0_0.dotR)
+
+        # Filter operator (limb darkening)
+        self.Fu = FOp(
+            ops_yf_u_0.F,
+            (self.ydeg + self.udeg + self.fdeg + 1) ** 2,
+            (self.ydeg + self.fdeg + 1) ** 2,
+        )
+
+        # Greens-to-poly change of basis
+        self.A2_Nyuf_x_Nyuf = ts.as_sparse_variable(ops_yuf_0_0.A2)
+
+    @autocompile
+    def tensordotRz(self, matrix, theta):
+        if self.ydeg + self.fdeg == 0:
+            return matrix
+        else:
+            return self._tensordotRz(matrix, theta)
+
+    @autocompile
+    def dotR(self, matrix, ux, uy, uz, theta):
+        if self.ydeg + self.fdeg == 0:
+            return matrix
+        else:
+            return self._dotR(matrix, ux, uy, uz, theta)
+
+    @autocompile
+    def render(self, res, projection, theta, inc, obl, fproj, y, u, f):
+        """Render the map on a Cartesian grid."""
+        # Compute the Cartesian grid
+        xyz = ifelse(
+            tt.eq(projection, STARRY_RECTANGULAR_PROJECTION),
+            self.compute_rect_grid(res)[-1],
+            ifelse(
+                tt.eq(projection, STARRY_MOLLWEIDE_PROJECTION),
+                self.compute_moll_grid(res)[-1],
+                self.compute_ortho_grid_fproj(res, obl, fproj)[-1],
+            ),
+        )
+
+        # Compute the polynomial basis
+        pT = self.pT(xyz[0], xyz[1], xyz[2])
+
+        # Apply the grav dark filter to the map
+        if self.fdeg > 0:
+            if self.ydeg == 0:
+                # Trivial!
+                y = f
+            else:
+                if self.nw is None:
+                    F = self.Fg(tt.as_tensor_variable([-1.0]), f)
+                    A1InvFA1 = ts.dot(ts.dot(self.A1Inv_Nyf_x_Nyf, F), self.A1)
+                    y = tt.dot(A1InvFA1, y)
+                else:
+                    # We need to compute the filter for each wavelength bin.
+                    # We could speed this up by vectorization on the C++ side
+                    A1InvFA1 = tt.zeros(
+                        (
+                            self.nw,
+                            (self.ydeg + self.fdeg + 1) ** 2,
+                            (self.ydeg + 1) ** 2,
+                        )
+                    )
+                    for i in range(self.nw):
+                        F = self.Fg(tt.as_tensor_variable([-1.0]), f[:, i])
+                        A1InvFA1 = tt.set_subtensor(
+                            A1InvFA1[i],
+                            ts.dot(ts.dot(self.A1Inv_Nyf_x_Nyf, F), self.A1),
+                        )
+                    y = tt.transpose(tt.batched_dot(A1InvFA1, tt.transpose(y)))
+
+        # If orthographic, rotate the map to the correct frame
+        if self.nw is None:
+            Ry = ifelse(
+                tt.eq(projection, STARRY_ORTHOGRAPHIC_PROJECTION),
+                self.left_project(
+                    tt.transpose(tt.tile(y, [theta.shape[0], 1])),
+                    inc,
+                    obl * 0,  # We already account for the obliquity in `pT`
+                    theta,
+                ),
+                tt.transpose(tt.tile(y, [theta.shape[0], 1])),
+            )
+        else:
+            Ry = ifelse(
+                tt.eq(projection, STARRY_ORTHOGRAPHIC_PROJECTION),
+                tt.reshape(
+                    self.left_project(
+                        y,
+                        inc,
+                        obl
+                        * 0,  # We already account for the obliquity in `pT`
+                        tt.tile(theta[0], self.nw),
+                    ),
+                    (-1, self.nw),
+                ),
+                y,
+            )
+
+        # Change basis to polynomials
+        A1Ry = ts.dot(self.A1_Nyf_x_Nyf, Ry)
+
+        # Apply the limb darkening *only if orthographic*
+        if self.udeg > 0:
+            A1Ry = ifelse(
+                tt.eq(projection, STARRY_ORTHOGRAPHIC_PROJECTION),
+                tt.dot(self.Fu(u, tt.as_tensor_variable([np.pi])), A1Ry),
+                A1Ry,
+            )
+
+        # Dot the polynomial into the basis
+        res = tt.reshape(tt.dot(pT, A1Ry), [res, res, -1])
+
+        # We need the shape to be (nframes, npix, npix)
+        return res.dimshuffle(2, 0, 1)
+
+    @autocompile
+    def compute_ortho_grid_fproj(self, res, obl, fproj):
+        dx = 2.0 / (res - 0.01)
+        y, x = tt.mgrid[-1:1:dx, -1:1:dx]
+        xp = x * tt.cos(obl) + y * tt.sin(obl)
+        yp = -x * tt.sin(obl) + y * tt.cos(obl)
+        x = xp
+        y = yp / (1 - fproj)
+        z = tt.sqrt(1 - x ** 2 - y ** 2)
+        y = tt.set_subtensor(y[tt.isnan(z)], np.nan)
+        x = tt.reshape(x, [1, -1])
+        y = tt.reshape(y, [1, -1])
+        z = tt.reshape(z, [1, -1])
+        lat = tt.reshape(0.5 * np.pi - tt.arccos(y), [1, -1])
+        lon = tt.reshape(tt.arctan2(x, z), [1, -1])
+        return tt.concatenate((lat, lon)), tt.concatenate((x, y, z))
+
+    @autocompile
+    def sT(self, f, theta, bo, ro):
+        return self._sT(f, theta, bo, ro)
+
+    @autocompile
+    def X(self, theta, xo, yo, zo, ro, inc, obl, fproj, u, f):
+        """Compute the light curve design matrix."""
+        # Determine shapes
+        rows = theta.shape[0]
+        cols = self.Ny
+        X = tt.zeros((rows, cols))
+
+        # Compute the occultation mask
+        bo = tt.sqrt(xo ** 2 + yo ** 2)
+        thetao = tt.arctan2(xo, yo)
+
+        # Occultation + phase curve operator
+        sT = self.sT(fproj, thetao, bo, ro)
+
+        # Limb darkening
+        if self.udeg > 0:
+            sTA2 = ts.dot(sT, self.A2_Nyuf_x_Nyuf)
+            F = self.Fu(u, tt.as_tensor_variable([np.pi]))
+            FA1 = ts.dot(F, self.A1_Nyf_x_Nyf)
+            sTA = tt.dot(sTA2, FA1)
+        else:
+            sTA = ts.dot(sT, self.A)
+
+        # Projection onto the sky
+        sTAR = self.right_project(sTA, inc, obl, theta)
+
+        # Gravity darkening
+        if self.nw is None and self.fdeg > 0:
+            F = self.Fg(tt.as_tensor_variable([-1.0]), f)
+            A1InvFA1 = ts.dot(ts.dot(self.A1Inv_Nyf_x_Nyf, F), self.A1)
+            sTAR = tt.dot(sTAR, A1InvFA1)
+
+        return sTAR
+
+    @autocompile
+    def flux(self, theta, xo, yo, zo, ro, inc, obl, fproj, y, u, f):
+        """Compute the light curve."""
+        if self.nw is None or self.fdeg == 0:
+            return tt.dot(
+                self.X(theta, xo, yo, zo, ro, inc, obl, fproj, u, f), y
+            )
+        else:
+
+            if self.ydeg == 0:
+                # Trivial!
+                return tt.dot(
+                    self.X(theta, xo, yo, zo, ro, inc, obl, fproj, u, f), f
+                )
+            else:
+
+                # We need to pre-weight the spherical harmonic vector by
+                # the filter in each wavelength bin
+                y = self.weight_ylms_by_grav_dark_filter(y, f)
+                return tt.dot(
+                    self.X(theta, xo, yo, zo, ro, inc, obl, fproj, u, f), y
+                )
+
+    @autocompile
+    def weight_ylms_by_grav_dark_filter(self, y, f):
+        # We need to compute the filter for each wavelength bin.
+        # We could speed this up by vectorization on the C++ side
+        A1InvFA1 = tt.zeros(
+            (self.nw, (self.ydeg + self.fdeg + 1) ** 2, (self.ydeg + 1) ** 2)
+        )
+        for i in range(self.nw):
+            F = self.Fg(tt.as_tensor_variable([-1.0]), f[:, i])
+            A1InvFA1 = tt.set_subtensor(
+                A1InvFA1[i], ts.dot(ts.dot(self.A1Inv_Nyf_x_Nyf, F), self.A1)
+            )
+        y = tt.transpose(tt.batched_dot(A1InvFA1, tt.transpose(y)))
+        return y
+
+    @autocompile
+    def grav_dark(self, z, wavnorm, omega, fobl, beta, tpole):
+        b = 1.0 - fobl
+        z2 = z * z
+        term = (-(omega ** 2) * (z2 * b ** 2 - z2 + 1) ** 1.5 + 1.0) ** 2.0
+        temp = (
+            tpole
+            * b ** (2 * beta)
+            * ((-z2 * b ** 2 + (z2 - 1) * term) / (-z2 * b ** 2 + z2 - 1) ** 3)
+            ** (0.5 * beta)
+        )
+        return 1.0 / (tt.exp(1.0 / (wavnorm * temp)) - 1.0)
+
+    @autocompile
+    def compute_grav_dark_filter(self, wavnorm, omega, fobl, beta, tpole):
+
+        # Compute the map expansion in the polar frame
+        nw = 1 if self.nw is None else self.nw
+        y = tt.zeros(((self.fdeg + 1) ** 2, nw))
+        y = tt.set_subtensor(
+            y[self.idx],
+            tt.dot(
+                self.SHT,
+                self.grav_dark(
+                    tt.reshape(self.z, (-1, 1)),
+                    tt.reshape(wavnorm, (1, -1)),
+                    omega,
+                    fobl,
+                    beta,
+                    tpole,
+                ),
+            ),
+        )
+
+        # Rotate down to the standard frame
+        y = tt.reshape(
+            tt.transpose(
+                self._dotR_Nf_x_Nf(
+                    tt.transpose(y),
+                    math.to_tensor(1.0),
+                    math.to_tensor(0.0),
+                    math.to_tensor(0.0),
+                    0.5 * np.pi,
+                )
+            ),
+            (-1,) if self.nw is None else (-1, self.nw),
+        )
+
+        return y
+
+
 class OpsSystem(object):
     """Class housing ops for modeling Keplerian systems."""
 
@@ -1366,6 +1710,7 @@ def __init__(
         secondaries,
         reflected=False,
         rv=False,
+        oblate=False,
         light_delay=False,
         texp=None,
         oversample=7,
@@ -1375,6 +1720,7 @@ def __init__(
         self.primary = primary
         self.secondaries = secondaries
         self._reflected = reflected
+        self._oblate = oblate
         self._rv = rv
         self.nw = self.primary._map.nw
         self.light_delay = light_delay
@@ -1448,6 +1794,7 @@ def X(
         pri_amp,
         pri_inc,
         pri_obl,
+        pri_fproj,
         pri_u,
         pri_f,
         sec_r,
@@ -1533,17 +1880,31 @@ def X(
         ) + tt.shape_padright(sec_theta0)
 
         # Compute all the phase curves
-        phase_pri = pri_amp * self.primary.map.ops.X(
-            theta_pri,
-            tt.zeros_like(t),
-            tt.zeros_like(t),
-            tt.zeros_like(t),
-            math.to_tensor(0.0),
-            pri_inc,
-            pri_obl,
-            pri_u,
-            pri_f,
-        )
+        if self._oblate:
+            phase_pri = pri_amp * self.primary.map.ops.X(
+                theta_pri,
+                tt.zeros_like(t),
+                tt.zeros_like(t),
+                tt.zeros_like(t),
+                math.to_tensor(0.0),
+                pri_inc,
+                pri_obl,
+                pri_fproj,
+                pri_u,
+                pri_f,
+            )
+        else:
+            phase_pri = pri_amp * self.primary.map.ops.X(
+                theta_pri,
+                tt.zeros_like(t),
+                tt.zeros_like(t),
+                tt.zeros_like(t),
+                math.to_tensor(0.0),
+                pri_inc,
+                pri_obl,
+                pri_u,
+                pri_f,
+            )
         if self._reflected:
             phase_sec = [
                 pri_amp
@@ -1606,23 +1967,43 @@ def X(
                 tt.ge(b, 1.0 + ro) | tt.le(zo, 0.0) | tt.eq(ro, 0.0)
             )
             idx = tt.arange(b.shape[0])[b_occ]
-            occ_pri = tt.set_subtensor(
-                occ_pri[idx],
-                occ_pri[idx]
-                + pri_amp
-                * self.primary.map.ops.X(
-                    theta_pri[idx],
-                    xo[idx],
-                    yo[idx],
-                    zo[idx],
-                    ro,
-                    pri_inc,
-                    pri_obl,
-                    pri_u,
-                    pri_f,
+            if self._oblate:
+                occ_pri = tt.set_subtensor(
+                    occ_pri[idx],
+                    occ_pri[idx]
+                    + pri_amp
+                    * self.primary.map.ops.X(
+                        theta_pri[idx],
+                        xo[idx],
+                        yo[idx],
+                        zo[idx],
+                        ro,
+                        pri_inc,
+                        pri_obl,
+                        pri_fproj,
+                        pri_u,
+                        pri_f,
+                    )
+                    - phase_pri[idx],
+                )
+            else:
+                occ_pri = tt.set_subtensor(
+                    occ_pri[idx],
+                    occ_pri[idx]
+                    + pri_amp
+                    * self.primary.map.ops.X(
+                        theta_pri[idx],
+                        xo[idx],
+                        yo[idx],
+                        zo[idx],
+                        ro,
+                        pri_inc,
+                        pri_obl,
+                        pri_u,
+                        pri_f,
+                    )
+                    - phase_pri[idx],
                 )
-                - phase_pri[idx],
-            )
 
         # Compute occultations by the primary
         for i, sec in enumerate(self.secondaries):
@@ -1636,7 +2017,29 @@ def X(
                 tt.ge(b, 1.0 + ro) | tt.le(zo, 0.0) | tt.eq(ro, 0.0)
             )
             idx = tt.arange(b.shape[0])[b_occ]
-            if self._reflected:
+            if self._oblate:
+                # TODO: Occultations *by* an oblate occultor are not
+                # currently supported. The following code ignores any
+                # oblateness and instead treats the body as a spherical
+                # occultor with radius equal to its equatorial radius.
+                occ_sec[i] = tt.set_subtensor(
+                    occ_sec[i][idx],
+                    occ_sec[i][idx]
+                    + sec_amp[i]
+                    * sec.map.ops.X(
+                        theta_sec[i, idx],
+                        xo[idx],
+                        yo[idx],
+                        zo[idx],
+                        ro,
+                        sec_inc[i],
+                        sec_obl[i],
+                        sec_u[i],
+                        sec_f[i],
+                    )
+                    - phase_sec[i][idx],
+                )
+            elif self._reflected:
                 occ_sec[i] = tt.set_subtensor(
                     occ_sec[i][idx],
                     occ_sec[i][idx]
@@ -1769,6 +2172,7 @@ def rv(
         pri_amp,
         pri_inc,
         pri_obl,
+        pri_fproj,
         pri_y,
         pri_u,
         pri_alpha,
@@ -1828,6 +2232,7 @@ def rv(
             pri_amp,
             pri_inc,
             pri_obl,
+            pri_fproj,
             pri_u,
             pri_f,
             sec_r,
@@ -1858,6 +2263,7 @@ def rv(
             pri_amp,
             pri_inc,
             pri_obl,
+            pri_fproj,
             pri_u,
             pri_f0,
             sec_r,
@@ -1945,6 +2351,7 @@ def render(
         pri_theta0,
         pri_inc,
         pri_obl,
+        pri_fproj,
         pri_y,
         pri_u,
         pri_f,
@@ -2002,16 +2409,29 @@ def render(
         ) + tt.shape_padright(sec_theta0)
 
         # Compute all the maps
-        img_pri = self.primary.map.ops.render(
-            res,
-            STARRY_ORTHOGRAPHIC_PROJECTION,
-            theta_pri,
-            pri_inc,
-            pri_obl,
-            pri_y,
-            pri_u,
-            pri_f,
-        )
+        if self._oblate:
+            img_pri = self.primary.map.ops.render(
+                res,
+                STARRY_ORTHOGRAPHIC_PROJECTION,
+                theta_pri,
+                pri_inc,
+                pri_obl,
+                pri_fproj,
+                pri_y,
+                pri_u,
+                pri_f,
+            )
+        else:
+            img_pri = self.primary.map.ops.render(
+                res,
+                STARRY_ORTHOGRAPHIC_PROJECTION,
+                theta_pri,
+                pri_inc,
+                pri_obl,
+                pri_y,
+                pri_u,
+                pri_f,
+            )
         if self._reflected:
             img_sec = tt.as_tensor_variable(
                 [
diff --git a/starry/_core/math.py b/starry/_core/math.py
index 85a2c89e4..393afd874 100644
--- a/starry/_core/math.py
+++ b/starry/_core/math.py
@@ -11,8 +11,14 @@
 
 
 # Cholesky solve
-_solve_lower = slinalg.Solve(A_structure="lower_triangular", lower=True)
-_solve_upper = slinalg.Solve(A_structure="upper_triangular", lower=False)
+try:
+    # theano
+    _solve_lower = slinalg.Solve(A_structure="lower_triangular", lower=True)
+    _solve_upper = slinalg.Solve(A_structure="upper_triangular", lower=False)
+except TypeError:
+    # aesara (TODO: Check me)
+    _solve_lower = slinalg.Solve(assume_a="sym", lower=True)
+    _solve_upper = slinalg.Solve(assume_a="sym", lower=False)
 
 
 def _cho_solve(cho_A, b):
diff --git a/starry/_core/ops/integration.py b/starry/_core/ops/integration.py
index 763e28bda..d0459420d 100644
--- a/starry/_core/ops/integration.py
+++ b/starry/_core/ops/integration.py
@@ -1,9 +1,9 @@
 # -*- coding: utf-8 -*-
-from ...compat import Apply, Op, tt
+from ...compat import Apply, Op, tt, floatX
 import numpy as np
 
 
-__all__ = ["sTOp", "rTReflectedOp", "sTReflectedOp"]
+__all__ = ["sTOp", "rTReflectedOp", "sTReflectedOp", "sTOblateOp"]
 
 
 class sTOp(Op):
@@ -204,3 +204,56 @@ def perform(self, node, inputs, outputs):
         outputs[2][0] = np.reshape(bbo, np.shape(bo))
         outputs[3][0] = np.array(np.reshape(bro, np.shape(ro)))
         outputs[4][0] = np.array(np.reshape(bsigr, np.shape(sigr)))
+
+
+class sTOblateOp(Op):
+    def __init__(self, func, N):
+        self.func = func
+        self.N = N
+        self._grad_op = sTOblateGradientOp(func)
+
+    def make_node(self, *inputs):
+        inputs = [tt.as_tensor_variable(i) for i in inputs]
+        outputs = [tt.TensorType(inputs[-1].dtype, (False, False))()]
+        return Apply(self, inputs, outputs)
+
+    def infer_shape(self, *args):
+        shapes = args[-1]
+        return [shapes[1] + (tt.as_tensor(self.N),)]
+
+    def R_op(self, inputs, eval_points):
+        if eval_points[0] is None:
+            return eval_points
+        return self.grad(inputs, eval_points)
+
+    def perform(self, node, inputs, outputs):
+        f, theta, bo, ro = inputs
+        sT = self.func(f, theta, bo, ro)
+        outputs[0][0] = sT
+
+    def grad(self, inputs, gradients):
+        return self._grad_op(*(inputs + gradients))
+
+
+class sTOblateGradientOp(Op):
+    def __init__(self, func):
+        self.func = func
+
+    def make_node(self, *inputs):
+        inputs = [tt.as_tensor_variable(i) for i in inputs]
+        outputs = [
+            i.astype(floatX).type() for i in inputs[:4]
+        ]  # TODO: Investigate why this is necessary sometimes
+        return Apply(self, inputs, outputs)
+
+    def infer_shape(self, *args):
+        shapes = args[-1]
+        return shapes[:4]
+
+    def perform(self, node, inputs, outputs):
+        f, theta, bo, ro, bsT = inputs
+        bf, btheta, bbo, bro = self.func(f, theta, bo, ro, bsT)
+        outputs[0][0] = np.reshape(bf, np.shape(f))
+        outputs[1][0] = np.reshape(btheta, np.shape(theta))
+        outputs[2][0] = np.reshape(bbo, np.shape(bo))
+        outputs[3][0] = np.array(np.reshape(bro, np.shape(ro)))
diff --git a/starry/_core/ops/lib/include/interface.cpp b/starry/_core/ops/lib/include/interface.cpp
index c460020d8..e8c10d818 100644
--- a/starry/_core/ops/lib/include/interface.cpp
+++ b/starry/_core/ops/lib/include/interface.cpp
@@ -113,6 +113,15 @@ PYBIND11_MODULE(_c_ops, m) {
 #endif
   });
 
+  // Change of basis matrix: poly to green's
+  Ops.def_property_readonly("A2", [](starry::Ops<Scalar> &ops) {
+#ifdef STARRY_MULTI
+    return (ops.B.A2.template cast<double>()).eval();
+#else
+    return ops.B.A2;
+#endif
+  });
+
   // Augmented change of basis matrix: Ylm to poly
   Ops.def_property_readonly("A1Big", [](starry::Ops<Scalar> &ops) {
 #ifdef STARRY_MULTI
@@ -252,6 +261,87 @@ PYBIND11_MODULE(_c_ops, m) {
     return py::make_tuple(result, ddb, ddtheta, ddbo, ddro, ddsigr);
   });
 
+  // Occultation of an oblate spheroid
+  Ops.def("sTOblate", [](starry::Ops<Scalar> &ops, const double &f_,
+                         const Vector<double> &theta_,
+                         const Vector<double> &bo_, const double &ro_) {
+    // Total number of terms in `s^T`
+    int K = theta_.size();
+
+    // Cast to ADScalar
+    ADScalar<Scalar, 0> f, theta, bo, ro;
+    f.value() = static_cast<Scalar>(f_);
+    ro.value() = static_cast<Scalar>(ro_);
+
+    // The output
+    Matrix<double> result(K, ops.N);
+
+    // Loop through the timeseries
+    for (int k = 0; k < K; ++k) {
+
+      // Compute sT for this timestep
+      theta.value() = static_cast<Scalar>(theta_(k));
+      bo.value() = static_cast<Scalar>(bo_(k));
+      ops.OBL.compute(bo, ro, f, theta);
+
+      // Process the ADScalar
+      for (int n = 0; n < ops.N; ++n) {
+        result(k, n) = static_cast<double>(ops.OBL.sT(n).value());
+      }
+    }
+    return result;
+  });
+
+  // Occultation in reflected light (backprop)
+  Ops.def("sTOblate", [](starry::Ops<Scalar> &ops, const double &f_,
+                         const Vector<double> &theta_,
+                         const Vector<double> &bo_, const double &ro_,
+                         const Matrix<double> &bsT) {
+    // Total number of terms in `s^T`
+    int K = theta_.size();
+
+    // Seed the derivatives
+    ADScalar<Scalar, 4> f, theta, bo, ro;
+    f.derivatives() = Vector<Scalar>::Unit(4, 0);
+    theta.derivatives() = Vector<Scalar>::Unit(4, 1);
+    bo.derivatives() = Vector<Scalar>::Unit(4, 2);
+    ro.derivatives() = Vector<Scalar>::Unit(4, 3);
+    f.value() = f_ < STARRY_MIN_F ? STARRY_MIN_F : static_cast<Scalar>(f_);
+    ro.value() = static_cast<Scalar>(ro_);
+
+    // The output
+    Matrix<double> ddf(K, ops.N);
+    Matrix<double> ddtheta(K, ops.N);
+    Matrix<double> ddbo(K, ops.N);
+    Matrix<double> ddro(K, ops.N);
+
+    // Loop through the timeseries
+    for (int k = 0; k < K; ++k) {
+
+      // Compute sT for this timestep
+      theta.value() = static_cast<Scalar>(theta_(k));
+      bo.value() = static_cast<Scalar>(bo_(k));
+      ops.OBLAD.compute(bo, ro, f, theta);
+
+      // Process the ADScalar
+      for (int n = 0; n < ops.N; ++n) {
+        ddf(k, n) = static_cast<double>(ops.OBLAD.sT(n).derivatives()(0));
+        ddtheta(k, n) = static_cast<double>(ops.OBLAD.sT(n).derivatives()(1));
+        ddbo(k, n) = static_cast<double>(ops.OBLAD.sT(n).derivatives()(2));
+        ddro(k, n) = static_cast<double>(ops.OBLAD.sT(n).derivatives()(3));
+      }
+    }
+
+    // Chain rule
+    double bf = bsT.cwiseProduct(ddf).sum();
+    Vector<double> btheta = bsT.cwiseProduct(ddtheta).rowwise().sum();
+    Vector<double> bbo = bsT.cwiseProduct(ddbo).rowwise().sum();
+    double bro = bsT.cwiseProduct(ddro).sum();
+
+    // Return the backprop grads
+    return py::make_tuple(bf, btheta, bbo, bro);
+  });
+
   // Rotation solution in emitted light dotted into Ylm space
   Ops.def_property_readonly("rTA1", [](starry::Ops<Scalar> &ops) {
     return ops.B.rTA1.template cast<double>();
@@ -415,4 +505,16 @@ PYBIND11_MODULE(_c_ops, m) {
                                                static_cast<Scalar>(a),
                                                static_cast<Scalar>(b));
         });
+
+#ifdef STARRY_UNIT_TESTS
+
+  m.attr("STARRY_UNIT_TESTS") = py::bool_(1);
+
+  // TODO: Unit tests for C++ functions here.
+
+#else
+
+  m.attr("STARRY_UNIT_TESTS") = py::bool_(0);
+
+#endif
 }
\ No newline at end of file
diff --git a/starry/_core/ops/lib/include/macros.h b/starry/_core/ops/lib/include/macros.h
new file mode 100644
index 000000000..f9d5df6b1
--- /dev/null
+++ b/starry/_core/ops/lib/include/macros.h
@@ -0,0 +1,233 @@
+/**
+\file macros.h
+\brief Compile-time-tunable macros
+
+*/
+
+#ifndef _STARRY_CONSTANTS_H_
+#define _STARRY_CONSTANTS_H_
+
+//! Number of digits of precision (16 = double)
+#ifndef STARRY_NDIGITS
+#define STARRY_NDIGITS 16
+#endif
+
+//! Number of Gaussian-Legendre quadrature points for numerical integration
+#ifndef STARRY_QUAD_POINTS
+#define STARRY_QUAD_POINTS 100
+#endif
+
+//! Max iterations in elliptic integrals
+#ifndef STARRY_ELLIP_MAX_ITER
+#define STARRY_ELLIP_MAX_ITER 200
+#endif
+
+//! Max iterations in computing the M & N integrals
+#ifndef STARRY_MN_MAX_ITER
+#define STARRY_MN_MAX_ITER 100
+#endif
+
+//! Max iterations in computing the I & J integrals
+#ifndef STARRY_IJ_MAX_ITER
+#define STARRY_IJ_MAX_ITER 200
+#endif
+
+//! Refine the downward recursion in the J integral at this index
+#ifndef STARRY_REFINE_J_AT
+#define STARRY_REFINE_J_AT 25
+#endif
+
+//! Cutoff value for `b` below which we reparametrize LD evaluation
+#ifndef STARRY_BCUT
+#define STARRY_BCUT 1.0e-3
+#endif
+
+//! Use the incomplete elliptic integrals to compute P?
+#ifndef STARRY_USE_INCOMPLETE_INTEGRALS
+#define STARRY_USE_INCOMPLETE_INTEGRALS 0
+#endif
+
+//! Maximum number of iterations when computing `el2`
+#ifndef STARRY_EL2_MAX_ITER
+#define STARRY_EL2_MAX_ITER 100
+#endif
+
+//! Square root of the desired precision in `el2`
+#ifndef STARRY_EL2_CA
+#define STARRY_EL2_CA 1e-8
+#endif
+
+//! Replace `inf` with this value in argument to `el2`
+#ifndef STARRY_HUGE_TAN
+#define STARRY_HUGE_TAN 1e15
+#endif
+
+//! Nudge k^2 away from 1 when it gets this close
+#ifndef STARRY_K2_ONE_TOL
+#define STARRY_K2_ONE_TOL 1e-12
+#endif
+
+//! Maximum number of iterations when computing 2F1
+#ifndef STARRY_2F1_MAXITER
+#define STARRY_2F1_MAXITER 200
+#endif
+
+//! Tolerance when computing 2F1
+#ifndef STARRY_2F1_MAXTOL
+#define STARRY_2F1_MAXTOL 1e-15
+#endif
+
+//! Tolerance when computing 2F1
+#ifndef STARRY_2F1_MINTOL
+#define STARRY_2F1_MINTOL 1e-12
+#endif
+
+//! Things currently go numerically unstable in our bases for high `l`
+#ifndef STARRY_MAX_LMAX
+#define STARRY_MAX_LMAX 50
+#endif
+
+//! If |sin(theta)| or |cos(theta)| is less than this, set  0
+#ifndef STARRY_T_TOL
+#define STARRY_T_TOL 1e-12
+#endif
+
+//! Tolerance for quartic root polishing
+#ifndef STARRY_ROOT_TOL_LOW
+#define STARRY_ROOT_TOL_LOW 1e-2
+#endif
+
+//! Tolerance for theta = +/- pi/2 in oblate case
+#ifndef STARRY_ROOT_TOL_THETA_PI_TWO
+#define STARRY_ROOT_TOL_THETA_PI_TWO 1e-5
+#endif
+
+//! Tolerance for quartic root polishing
+#ifndef STARRY_ROOT_TOL_MED
+#define STARRY_ROOT_TOL_MED 1e-10
+#endif
+
+//! Tolerance for quartic root polishing
+#ifndef STARRY_ROOT_TOL_HIGH
+#define STARRY_ROOT_TOL_HIGH 1e-15
+#endif
+
+//! Tolerance for duplicate roots
+#ifndef STARRY_ROOT_TOL_DUP
+#define STARRY_ROOT_TOL_DUP 1e-7
+#endif
+
+//! Tolerance for quartic root polishing
+#ifndef STARRY_ROOT_TOL_FINAL
+#define STARRY_ROOT_TOL_FINAL 1e-13
+#endif
+
+//! Maximum number of root polishing iterations
+#ifndef STARRY_ROOT_MAX_ITER
+#define STARRY_ROOT_MAX_ITER 30
+#endif
+
+//! If |b| is less than this value, set equal to 0
+#define STARRY_B_ZERO_TOL 1e-8
+
+//! If |b| is within this value of 1, set equal to 1
+#ifndef STARRY_B_ONE_TOL
+#define STARRY_B_ONE_TOL 1e-6
+#endif
+
+//! Minimum value of the oblateness `f`
+#ifndef STARRY_MIN_F
+#define STARRY_MIN_F 1e-6
+#endif
+
+//! Tolerance at singular point `bo ~ 1 - ro`
+#ifndef STARRY_BO_EQUALS_ONE_MINUS_RO_TOL
+#define STARRY_BO_EQUALS_ONE_MINUS_RO_TOL 1e-6
+#endif
+
+//! Tolerance for various functions that calculate phi, xi, and lam
+#ifndef STARRY_ANGLE_TOL
+#define STARRY_ANGLE_TOL 1e-13
+#endif
+
+//! Nudge bo away from ro when it gets this close
+#ifndef STARRY_BO_EQUALS_RO_TOL
+#define STARRY_BO_EQUALS_RO_TOL 1e-8
+#endif
+
+//! Nudge bo away from ro when bo ~ ro ~ 0.5
+#ifndef STARRY_BO_EQUALS_RO_EQUALS_HALF_TOL
+#define STARRY_BO_EQUALS_RO_EQUALS_HALF_TOL 1e-5
+#endif
+
+//! Nudge bo away from zero when it gets this close (oblate case)
+#ifndef STARRY_BO_EQUALS_ZERO_TOL
+#define STARRY_BO_EQUALS_ZERO_TOL 1e-12
+#endif
+
+//! Nudge ro away from zero when it gets this close (oblate case)
+#ifndef STARRY_RO_EQUALS_ZERO_TOL
+#define STARRY_RO_EQUALS_ZERO_TOL 1e-12
+#endif
+
+//! Nudge theta away from pi / 2 when ro = 1
+#ifndef STARRY_THETA_UNIT_RADIUS_TOL
+#define STARRY_THETA_UNIT_RADIUS_TOL 1e-5
+#endif
+
+//! Smallest value of sin(kappa / 2) in computing W
+#ifndef STARRY_MIN_SIN_ALPHA
+#define STARRY_MIN_SIN_ALPHA 1e-8
+#endif
+
+//! Tolerance for complete occultations
+#ifndef STARRY_COMPLETE_OCC_TOL
+#define STARRY_COMPLETE_OCC_TOL 1e-8
+#endif
+
+//! Tolerance for no occultation
+#ifndef STARRY_NO_OCC_TOL
+#define STARRY_NO_OCC_TOL 1e-8
+#endif
+
+//! Tolerance for grazing occultations
+#ifndef STARRY_GRAZING_TOL
+#define STARRY_GRAZING_TOL 1e-7
+#endif
+
+//! Tolerance for the Pal (2012) solver
+#ifndef STARRY_CRJ_MAX_ITER
+#define STARRY_CRJ_MAX_ITER 100
+#endif
+
+//! Tolerance for the Pal (2012) solver
+#ifndef STARRY_CRJ_LO_LIM
+#define STARRY_CRJ_LO_LIM 2e-26
+#endif
+
+//! Tolerance for the Pal (2012) solver
+#ifndef STARRY_CRJ_HI_LIM
+#define STARRY_CRJ_HI_LIM 3e24
+#endif
+
+//! Tolerance for the Pal (2012) solver
+#ifndef STARRY_CRJ_TOL
+#define STARRY_CRJ_TOL 2e-2
+#endif
+
+//! Tolerance for the Pal (2012) solver
+#ifndef STARRY_PAL_BO_EQUALS_RO_TOL
+#define STARRY_PAL_BO_EQUALS_RO_TOL 1e-3
+#endif
+
+//! Tolerance for the Pal (2012) solver
+#ifndef STARRY_PAL_BO_EQUALS_RO_MINUS_ONE_TOL
+#define STARRY_PAL_BO_EQUALS_RO_MINUS_ONE_TOL 1e-3
+#endif
+
+//! Tolerance for the Pal (2012) solver
+#ifndef STARRY_PAL_BO_EQUALS_ONE_MINUS_RO_TOL
+#define STARRY_PAL_BO_EQUALS_ONE_MINUS_RO_TOL 1e-3
+#endif
+
+#endif
\ No newline at end of file
diff --git a/starry/_core/ops/lib/include/oblate/geometry.h b/starry/_core/ops/lib/include/oblate/geometry.h
new file mode 100644
index 000000000..7a52d3ad8
--- /dev/null
+++ b/starry/_core/ops/lib/include/oblate/geometry.h
@@ -0,0 +1,417 @@
+/**
+\file geometry.h
+\brief Circle-ellipse intersection stuff.
+
+*/
+
+#ifndef _STARRY_OBLATE_GEOMETRY_H_
+#define _STARRY_OBLATE_GEOMETRY_H_
+
+#include "../utils.h"
+#include <Eigen/Eigenvalues>
+
+namespace starry {
+namespace oblate {
+namespace geometry {
+
+using namespace utils;
+
+/**
+    Polynomial root finder using an eigensolver.
+    `coeffs` is a vector of coefficients, highest power first.
+
+    Adapted from http://www.sgh1.net/posts/cpp-root-finder.md
+*/
+template <class Scalar>
+inline std::vector<std::complex<Scalar>>
+eigen_roots(const std::vector<Scalar> &coeffs, bool &success) {
+
+  int matsz = coeffs.size() - 1;
+  std::vector<std::complex<Scalar>> vret;
+  Matrix<Scalar> companion_mat(matsz, matsz);
+  companion_mat.setZero();
+
+  for (int n = 0; n < matsz; ++n) {
+    for (int m = 0; m < matsz; ++m) {
+
+      if (n == m + 1)
+        companion_mat(n, m) = 1.0;
+
+      if (m == matsz - 1)
+        companion_mat(n, m) = -coeffs[matsz - n] / coeffs[0];
+    }
+  }
+
+  Eigen::EigenSolver<Matrix<Scalar>> solver(companion_mat);
+  if (solver.info() == Eigen::Success) {
+    success = true;
+  } else {
+    success = false;
+  }
+
+  Matrix<std::complex<Scalar>> eig = solver.eigenvalues();
+  for (int i = 0; i < matsz; ++i)
+    vret.push_back(eig(i));
+
+  return vret;
+}
+
+/**
+    Compute the points of intersection between a circle and an ellipse
+    in the frame where the ellipse is centered at the origin,
+    the semi-major axis of the ellipse is aligned with the x axis,
+    and the circle is centered at `(xo, yo)`.
+
+*/
+template <class Scalar, int N>
+inline void
+get_roots(const ADScalar<Scalar, N> &b_, const ADScalar<Scalar, N> &theta_,
+          const ADScalar<Scalar, N> &costheta_,
+          const ADScalar<Scalar, N> &sintheta_, const ADScalar<Scalar, N> &bo_,
+          const ADScalar<Scalar, N> &ro_, Vector<ADScalar<Scalar, N>> &x,
+          Vector<ADScalar<Scalar, N>> &y) {
+
+  // Initialize
+  x.resize(4);
+  y.resize(4);
+
+  // Get the *values*
+  using Complex = std::complex<Scalar>;
+  Scalar b = b_.value();
+  Scalar costheta = costheta_.value();
+  Scalar sintheta = sintheta_.value();
+  Scalar bo = bo_.value();
+  Scalar ro = ro_.value();
+
+  // We'll solve for circle-ellipse intersections
+  // in the frame where the ellipse is centered at the origin,
+  // the semi-major axis of the ellipse is aligned with the x axis,
+  // and the circle is centered at `(xo, yo)`.
+  Scalar xo = bo * sintheta;
+  Scalar yo = bo * costheta;
+
+  // Useful quantities
+  Scalar b2 = b * b;
+  Scalar b4 = b2 * b2;
+  Scalar ro2 = ro * ro;
+  Scalar xo2 = xo * xo;
+  Scalar yo2 = yo * yo;
+
+  // Get the roots (eigenvalue problem)
+  std::vector<Scalar> coeffs;
+  coeffs.push_back((1 - b2) * (1 - b2));
+  coeffs.push_back(-4 * xo * (1 - b2));
+  coeffs.push_back(-2 *
+                   (b4 + ro2 - 3 * xo2 - yo2 - b2 * (1 + ro2 - xo2 + yo2)));
+  coeffs.push_back(-4 * xo * (b2 - ro2 + xo2 + yo2));
+  coeffs.push_back(b4 - 2 * b2 * (ro2 - xo2 + yo2) +
+                   (ro2 - xo2 - yo2) * (ro2 - xo2 - yo2));
+  bool success = false;
+  std::vector<Complex> roots = eigen_roots(coeffs, success);
+  if (!success) {
+    std::stringstream args;
+    args << "b_ = " << b_ << ", "
+         << "theta_ = " << theta_ << ", "
+         << "costheta_ = " << costheta_ << ", "
+         << "sintheta_ = " << sintheta_ << ", "
+         << "bo_ = " << bo_ << ", "
+         << "ro_ = " << ro_;
+    throw StarryException("Root eigensolver did not converge.",
+                          "oblate/geometry.h", "get_roots", args.str());
+  }
+
+  // Apply Newton's method to polish the roots
+  int nroots = 0;
+  for (int n = 0; n < 4; ++n) {
+
+    Complex root = roots[n];
+    Complex minxc = root;
+    Complex root2, root3, root4;
+    Complex f, df;
+    Scalar absf, minerr = INFINITY;
+
+    // Polish
+    for (int k = 0; k < STARRY_ROOT_MAX_ITER; ++k) {
+
+      // Compute the error
+      root2 = root * root;
+      root3 = root2 * root;
+      root4 = root3 * root;
+      f = coeffs[0] * root4 + coeffs[1] * root3 + coeffs[2] * root2 +
+          coeffs[3] * root + coeffs[4];
+      absf = abs(f);
+      if (absf <= minerr) {
+        minerr = absf;
+        minxc = root;
+        if (minerr <= STARRY_ROOT_TOL_HIGH)
+          break;
+      }
+
+      // Take a step
+      df = 4.0 * coeffs[0] * root3 + 3.0 * coeffs[1] * root2 +
+           2.0 * coeffs[2] * root + coeffs[3];
+      if (df == 0.0)
+        break;
+      root -= f / df;
+    }
+    root = minxc;
+
+    // Keep the root if it's real
+    if (abs(root.imag()) < STARRY_ROOT_TOL_HIGH) {
+
+      // Nudge the root away from the endpoints
+      if (root.real() > 1)
+        root = 1.0 - STARRY_ROOT_TOL_HIGH;
+      else if (root.real() < -1)
+        root = -1.0 + STARRY_ROOT_TOL_HIGH;
+      else if (root.real() < xo - ro)
+        root = xo - ro + STARRY_ROOT_TOL_HIGH;
+      else if (root.real() > xo + ro)
+        root = xo + ro - STARRY_ROOT_TOL_HIGH;
+
+      // Determine the y value of the point on the ellipse
+      // corresponding to each root and the signs of the
+      // functions describing the intersecting circle &
+      // ellipse segments.
+      Scalar s0, s1;
+      Complex fA, fB;
+      fA = b * sqrt(Scalar(1.0) - root * root);
+      fB = sqrt(ro2 - (root - xo) * (root - xo));
+      Vector<Scalar> diff(4);
+      diff <<                   //
+          abs(fA - (yo + fB)),  //
+          abs(fA - (yo - fB)),  //
+          abs(-fA - (yo + fB)), //
+          abs(-fA - (yo - fB)); //
+      typename Vector<Scalar>::Index idx(4);
+      diff.minCoeff(&idx);
+      if (idx < 2) {
+        s0 = 1.0;
+      } else {
+        s0 = -1.0;
+      }
+      if (is_even(idx)) {
+        s1 = 1.0;
+      } else {
+        s1 = -1.0;
+      }
+
+      // Save the root
+      x(nroots).value() = root.real();
+      y(nroots).value() = s0 * fA.real();
+
+      // Compute its derivatives
+      if (N > 0) {
+        Scalar dxdb, dxdtheta, dxdbo, dxdro;
+        Scalar p, q, v, w, t;
+        q = sqrt(ro2 - (root.real() - xo) * (root.real() - xo));
+        p = sqrt(1 - root.real() * root.real());
+        v = (root.real() - xo) / q;
+        w = b / p;
+        t = 1.0 / (w * root.real() - (s1 * s0) * v);
+        dxdb = t * p;
+        dxdtheta =
+            -(s1 * costheta * v - sintheta) * (bo * t * s0); // TODO: Check me
+        dxdbo = -(costheta + s1 * sintheta * v) * (t * s0);
+        dxdro = -ro * t / q * s1 * s0;
+        x(nroots).derivatives() =
+            dxdb * b_.derivatives() + dxdtheta * theta_.derivatives() +
+            dxdbo * bo_.derivatives() + dxdro * ro_.derivatives();
+        y(nroots) = s0 * b_ * sqrt(1.0 - x(nroots) * x(nroots));
+      }
+
+      ++nroots;
+    }
+  }
+
+  // We're done!
+  x.conservativeResize(nroots);
+  y.conservativeResize(nroots);
+}
+
+/**
+    Compute the angles at which the circle intersects the ellipse.
+
+*/
+template <class Scalar, int N>
+inline void
+get_angles(const ADScalar<Scalar, N> &bo_, const ADScalar<Scalar, N> &ro_,
+           const ADScalar<Scalar, N> &f_, const ADScalar<Scalar, N> &theta_,
+           ADScalar<Scalar, N> &phi1, ADScalar<Scalar, N> &phi2,
+           ADScalar<Scalar, N> &xi1, ADScalar<Scalar, N> &xi2) {
+
+  using A = ADScalar<Scalar, N>;
+
+  // We may need to adjust these, so make a copy
+  A bo = bo_;
+  A ro = ro_;
+  A f = f_;
+  A b = 1 - f_;
+  A theta = theta_;
+
+  // Enforce bo >= 0
+  if (bo < 0) {
+    bo = -bo;
+    theta -= pi<Scalar>();
+  }
+
+  A costheta = cos(theta);
+  A sintheta = sin(theta);
+
+  // Trivial cases
+  if (bo <= ro - 1 + STARRY_COMPLETE_OCC_TOL) {
+
+    // Complete occultation
+    phi1 = phi2 = xi1 = 0.0;
+    xi2 = 2 * pi<Scalar>();
+    return;
+
+  } else if (bo + ro + f <= 1 + STARRY_GRAZING_TOL) {
+
+    // Regular occultation, but occultor doesn't touch the limb
+    phi1 = xi1 = xi2 = 0.0;
+    phi2 = 2 * pi<Scalar>();
+    return;
+
+  } else if (bo >= 1 + ro - STARRY_NO_OCC_TOL) {
+
+    // No occultation
+    phi1 = phi2 = xi1 = xi2 = 0.0;
+    return;
+  }
+
+  // This grazing configuration leads to instabilities
+  // in the root solver. Let's avoid it.
+  if ((1 - ro - STARRY_GRAZING_TOL <= bo) &&
+      (bo <= 1 - ro + STARRY_GRAZING_TOL))
+    bo = 1 - ro + STARRY_GRAZING_TOL;
+
+  // The eigensolver doesn't converge when ro = 1 and theta = pi / 2.
+  if ((abs(1 - ro) < STARRY_THETA_UNIT_RADIUS_TOL) &&
+      (abs(costheta) < STARRY_THETA_UNIT_RADIUS_TOL)) {
+    costheta += (costheta > 0 ? STARRY_THETA_UNIT_RADIUS_TOL
+                              : -STARRY_THETA_UNIT_RADIUS_TOL);
+  }
+
+  // Get the points of intersection between the circle & ellipse
+  // These are the roots to a quartic equation.
+  A xo = bo * sintheta;
+  A yo = bo * costheta;
+  Vector<A> x, y;
+  get_roots(b, theta, costheta, sintheta, bo, ro, x, y);
+
+  int nroots = x.size();
+
+  if (nroots == 0) {
+
+    // There are no intersections between the circle
+    // and the ellipse. There are 3 possibilies.
+
+    // Is the center of the circle outside the ellipse?
+    if ((abs(xo) > 1) || (abs(yo) > b * sqrt(1 - xo * xo))) {
+
+      // Is the center of the ellipse outside the circle?
+      if (bo > ro) {
+
+        // No occultation
+        phi1 = phi2 = xi1 = xi2 = 0.0;
+
+      } else {
+
+        // Complete occultation
+        phi1 = phi2 = xi1 = 0.0;
+        xi2 = 2 * pi<Scalar>();
+      }
+
+    } else {
+
+      // Regular occultation, but occultor doesn't touch the limb
+      phi1 = xi1 = xi2 = 0.0;
+      phi2 = 2 * pi<Scalar>();
+    }
+
+  } else if (nroots == 1) {
+
+    // Grazing configuration?
+    // TODO: Perturb and repeat.
+    std::stringstream args;
+    args << "bo_ = " << bo_ << ", "
+         << "ro_ = " << ro_ << ", "
+         << "f_ = " << f_ << ", "
+         << "theta_ = " << theta_;
+    throw StarryException("Unexpected branch.", "oblate/geometry.h",
+                          "get_angles", args.str());
+
+  } else if (nroots == 2) {
+
+    // Regular occultation
+    A xm, ym, mid, term;
+
+    // Get the angles
+    phi1 = theta + atan2(y(0) - yo, x(0) - xo);
+    phi2 = theta + atan2(y(1) - yo, x(1) - xo);
+    term = 1 - x(0) * x(0);
+    xi1 = atan2(sqrt(term < 0 ? 0 : term), x(0));
+    if (y(0) < 0)
+      xi1 *= -1;
+    term = 1 - x(1) * x(1);
+    xi2 = atan2(sqrt(term < 0 ? 0 : term), x(1));
+    if (y(1) < 0)
+      xi2 *= -1;
+
+    // Wrap and sort the angles
+    phi1 = angle(phi1);
+    phi2 = angle(phi2);
+    xi1 = angle(xi1);
+    xi2 = angle(xi2);
+
+    // xi is always counter-clockwise
+    if (xi1 > xi2) {
+      std::swap(xi1, xi2);
+      std::swap(phi1, phi2);
+    }
+
+    // Ensure the T integral takes us through the inside of the occultor
+    mid = 0.5 * (xi1 + xi2);
+    xm = cos(mid);
+    ym = b * sin(mid);
+    if ((xm - xo) * (xm - xo) + (ym - yo) * (ym - yo) >= ro * ro) {
+      std::swap(xi1, xi2);
+      xi2 += 2 * pi<Scalar>();
+    }
+
+    // Ensure the P integral takes us through the inside of the star
+    mid = 0.5 * (phi1 + phi2);
+    xm = xo + ro * cos(theta - mid);
+    ym = yo - ro * sin(theta - mid);
+    if (ym * ym > b * b * (1 - xm * xm)) {
+      if (phi1 < phi2) {
+        phi1 += 2 * pi<Scalar>();
+      } else {
+        phi2 += 2 * pi<Scalar>();
+      }
+    }
+
+    // phi is always counter-clockwise
+    if (phi1 > phi2) {
+      std::swap(phi1, phi2);
+    }
+
+  } else {
+
+    // Pathological case?
+    std::stringstream args;
+    args << "bo_ = " << bo_ << ", "
+         << "ro_ = " << ro_ << ", "
+         << "f_ = " << f_ << ", "
+         << "theta_ = " << theta_;
+    throw StarryException("Unexpected branch.", "oblate/geometry.h",
+                          "get_angles", args.str());
+  }
+}
+
+} // namespace geometry
+} // namespace oblate
+} // namespace starry
+
+#endif
diff --git a/starry/_core/ops/lib/include/oblate/numerical.h b/starry/_core/ops/lib/include/oblate/numerical.h
new file mode 100644
index 000000000..b0b6c7ea7
--- /dev/null
+++ b/starry/_core/ops/lib/include/oblate/numerical.h
@@ -0,0 +1,209 @@
+/**
+\file special.h
+\brief Numerical integration.
+
+*/
+
+#ifndef _STARRY_OBLATE_NUMERICAL_H_
+#define _STARRY_OBLATE_NUMERICAL_H_
+
+#include "../utils.h"
+
+namespace starry {
+namespace oblate {
+namespace numerical {
+
+using std::abs;
+using namespace utils;
+
+template <typename T>
+using IntegralVector = Eigen::Matrix<T, STARRY_QUAD_POINTS, 1>;
+
+template <typename T>
+using IntegralRowVector = Eigen::Matrix<T, 1, STARRY_QUAD_POINTS>;
+
+template <typename T>
+using IntegralArray = Eigen::Array<T, STARRY_QUAD_POINTS, 1>;
+
+template <typename T>
+using IntegralMatrix = Eigen::Matrix<T, STARRY_QUAD_POINTS, Eigen::Dynamic>;
+
+/*! Implementation of Gauss-Legendre quadrature based on
+ *  http://en.wikipedia.org/wiki/Gaussian_quadrature
+ *  http://rosettacode.org/wiki/Numerical_integration/Gauss-Legendre_Quadrature
+ */
+template <typename Scalar, int N> class Quad {
+
+  using A = ADScalar<Scalar, N>;
+
+private:
+  A p;
+  IntegralVector<A> q;
+  IntegralRowVector<A> vT;
+
+public:
+  template <typename Function>
+  inline Vector<A> integrate(const A &a, const A &b, Function f,
+                             const int &nint) {
+    const LegendrePolynomial &L = s_LegendrePolynomial;
+
+    // Initialize
+    IntegralMatrix<A> integrand;
+    integrand.setZero(STARRY_QUAD_POINTS, nint);
+
+    // Compute the function on the grid
+    p = 0.5 * (b - a);
+    vT = p * L.wT;
+    q.setConstant(0.5 * (b + a));
+    f(p * L.r + q, integrand);
+    return vT * integrand;
+  }
+
+private:
+  class LegendrePolynomial {
+  public:
+    IntegralVector<Scalar> r;
+    IntegralRowVector<Scalar> wT;
+
+    LegendrePolynomial() {
+      for (int i = 0; i < STARRY_QUAD_POINTS; ++i) {
+        Scalar dr = 1;
+        Evaluation eval(cos(M_PI * (i + 0.75) / (STARRY_QUAD_POINTS + 0.5)));
+        do {
+          dr = eval.v() / eval.d();
+          eval.evaluate(eval.x() - dr);
+        } while (abs(dr) > 2e-16);
+        this->r(i) = eval.x();
+        this->wT(i) = 2 / ((1 - eval.x() * eval.x()) * eval.d() * eval.d());
+      }
+    }
+
+    Scalar root(int i) const { return this->r(i); }
+    Scalar weight(int i) const { return this->wT(i); }
+
+  private:
+    class Evaluation {
+    public:
+      explicit Evaluation(Scalar x) : _x(x), _v(1), _d(0) { this->evaluate(x); }
+
+      void evaluate(Scalar x) {
+        this->_x = x;
+        Scalar vsub1 = x;
+        Scalar vsub2 = 1;
+        Scalar f = 1 / (x * x - 1);
+        for (int i = 2; i <= STARRY_QUAD_POINTS; ++i) {
+          this->_v = ((2 * i - 1) * x * vsub1 - (i - 1) * vsub2) / i;
+          this->_d = i * f * (x * this->_v - vsub1);
+          vsub2 = vsub1;
+          vsub1 = this->_v;
+        }
+      }
+
+      Scalar v() const { return this->_v; }
+      Scalar d() const { return this->_d; }
+      Scalar x() const { return this->_x; }
+
+    private:
+      Scalar _x;
+      Scalar _v;
+      Scalar _d;
+    };
+  };
+
+  static LegendrePolynomial s_LegendrePolynomial;
+};
+
+template <typename Scalar, int N>
+typename Quad<Scalar, N>::LegendrePolynomial
+    Quad<Scalar, N>::s_LegendrePolynomial;
+
+template <typename Scalar, int N>
+inline void pTodd_integrand(const int &deg, const ADScalar<Scalar, N> &bo,
+                            const ADScalar<Scalar, N> &ro,
+                            const ADScalar<Scalar, N> &f,
+                            const ADScalar<Scalar, N> &theta,
+                            const ADScalar<Scalar, N> &costheta,
+                            const ADScalar<Scalar, N> &sintheta,
+                            const IntegralVector<ADScalar<Scalar, N>> &phi,
+                            IntegralMatrix<ADScalar<Scalar, N>> &integrand) {
+
+  using Array = IntegralArray<ADScalar<Scalar, N>>;
+
+  // Pre-compute common terms
+  int ncoeff = (deg + 1) * (deg + 1);
+  Array cosphi = cos(phi.array());
+  Array sinphi = sin(phi.array());
+  Array cosvphi = costheta * cosphi + sintheta * sinphi;
+  Array sinvphi = -sintheta * cosphi + costheta * sinphi;
+  Array x = ro * cosvphi + bo * sintheta;
+  Array y = (ro * sinvphi + bo * costheta) / (1.0 - f);
+  Array z2 = 1.0 - x * x - y * y;
+  Array z = (z2 <= 0).select(Array::Zero(), sqrt(z2));
+  Array z3 = z2 * z;
+  Array z3x = -ro * (1.0 - f) * sinvphi * z3;
+  Array z3y = ro * cosvphi * z3;
+
+  // Case 2
+  integrand.col(2) = ro * (ro + bo * sinphi) * (1.0 - z3) / (3.0 * (1.0 - z2));
+
+  // Cases 3-5
+  Array xi, yj;
+  int n;
+  xi.setOnes();
+  for (int i = 0; i < deg - 1; ++i) {
+    if (is_even(i)) {
+      // Case 3
+      n = i * i + 6 * i + 7;
+      if (n < ncoeff) {
+        integrand.col(n) = xi * z3x;
+      }
+      // Case 4
+      n += 2 * i + 7;
+      if (n < ncoeff) {
+        integrand.col(n) = xi * y * z3x;
+      }
+    }
+    // Case 5
+    yj.setOnes();
+    for (int j = 0; j < deg - 1 - i; ++j) {
+      n = (i + j) * (i + j) + 4 * i + 6 * j + 5;
+      integrand.col(n) = xi * yj * z3y;
+      yj = yj * y;
+    }
+    xi = xi * x;
+  }
+}
+
+/**
+  Numerical evaluation of the `pT` integral for all odd `mu`.
+
+*/
+template <typename Scalar, int N>
+inline Vector<ADScalar<Scalar, N>>
+pTodd(const int &deg, const ADScalar<Scalar, N> &bo,
+      const ADScalar<Scalar, N> &ro, const ADScalar<Scalar, N> &f,
+      const ADScalar<Scalar, N> &theta, const ADScalar<Scalar, N> &costheta,
+      const ADScalar<Scalar, N> &sintheta, const ADScalar<Scalar, N> &phi1,
+      const ADScalar<Scalar, N> &phi2, Quad<Scalar, N> &QUAD) {
+  std::function<void( //
+      const IntegralVector<ADScalar<Scalar, N>> &,
+      IntegralMatrix<ADScalar<Scalar, N>> & //
+      )>
+      func =                                              //
+      [deg, bo, ro, f, theta, costheta, sintheta]         //
+      (                                                   //
+          const IntegralVector<ADScalar<Scalar, N>> &phi, //
+          IntegralMatrix<ADScalar<Scalar, N>> &integrand  //
+      ) {
+        pTodd_integrand(deg, bo, ro, f, theta, costheta, sintheta, phi,
+                        integrand);
+      };
+  int ncoeff = (deg + 1) * (deg + 1);
+  return QUAD.integrate(phi1, phi2, func, ncoeff);
+}
+
+} // namespace numerical
+} // namespace oblate
+} // namespace starry
+
+#endif
diff --git a/starry/_core/ops/lib/include/oblate/oblate.cpp b/starry/_core/ops/lib/include/oblate/oblate.cpp
new file mode 100644
index 000000000..66bc1c6e9
--- /dev/null
+++ b/starry/_core/ops/lib/include/oblate/oblate.cpp
@@ -0,0 +1,85 @@
+// Enable debug mode?
+#ifdef STARRY_DEBUG
+#undef NDEBUG
+#endif
+
+// Includes
+#include "../utils.h"
+#include "geometry.h"
+#include "occultation.h"
+#include <iostream>
+#include <pybind11/eigen.h>
+#include <pybind11/embed.h>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+namespace py = pybind11;
+
+// Register the Python module
+PYBIND11_MODULE(_c_ops, m) {
+  using namespace starry::utils;
+  using namespace starry::oblate::occultation;
+  using namespace starry::oblate::geometry;
+
+  m.def("sT", [](const int &deg, const double &bo_, const double &ro_,
+                 const double &f_, const double &theta_, const int &nruns) {
+    using A = ADScalar<double, 0>;
+    A bo, ro, f, theta;
+    bo.value() = bo_;
+    ro.value() = ro_;
+    f.value() = f_;
+    theta.value() = theta_;
+    auto occ = Occultation<double, 0>(deg);
+
+    for (int n = 0; n < nruns; ++n)
+      occ.compute(bo, ro, f, theta);
+
+    Vector<double> sT_value(occ.sT.size());
+    for (int n = 0; n < occ.sT.size(); ++n)
+      sT_value(n) = occ.sT(n).value();
+    return sT_value;
+  });
+
+  m.def("dsT", [](const int &deg, const double &bo_, const double &ro_,
+                  const double &f_, const double &theta_, const int &nruns) {
+    using A = ADScalar<double, 4>;
+    A bo, ro, f, theta;
+    bo.value() = bo_;
+    bo.derivatives() = Vector<double>::Unit(4, 0);
+    ro.value() = ro_;
+    ro.derivatives() = Vector<double>::Unit(4, 1);
+    f.value() = f_;
+    f.derivatives() = Vector<double>::Unit(4, 2);
+    theta.value() = theta_;
+    theta.derivatives() = Vector<double>::Unit(4, 3);
+    auto occ = Occultation<double, 4>(deg);
+
+    for (int n = 0; n < nruns; ++n)
+      occ.compute(bo, ro, f, theta);
+
+    Vector<double> sT_value(occ.sT.size());
+    Vector<double> Dbo(occ.sT.size());
+    Vector<double> Dro(occ.sT.size());
+    Vector<double> Df(occ.sT.size());
+    Vector<double> Dtheta(occ.sT.size());
+    for (int n = 0; n < occ.sT.size(); ++n) {
+      Dbo(n) = occ.sT(n).derivatives()(0);
+      Dro(n) = occ.sT(n).derivatives()(1);
+      Df(n) = occ.sT(n).derivatives()(2);
+      Dtheta(n) = occ.sT(n).derivatives()(3);
+    }
+    return py::make_tuple(Dbo, Dro, Df, Dtheta);
+  });
+
+  m.def("angles", [](const double &bo_, const double &ro_, const double &f_,
+                     const double &theta_) {
+    using A = ADScalar<double, 0>;
+    A bo, ro, f, theta, phi1, phi2, xi1, xi2;
+    bo.value() = bo_;
+    ro.value() = ro_;
+    f.value() = f_;
+    theta.value() = theta_;
+    get_angles(bo, ro, f, theta, phi1, phi2, xi1, xi2);
+    return py::make_tuple(phi1.value(), phi2.value(), xi1.value(), xi2.value());
+  });
+}
\ No newline at end of file
diff --git a/starry/_core/ops/lib/include/oblate/oblate/__init__.py b/starry/_core/ops/lib/include/oblate/oblate/__init__.py
new file mode 100644
index 000000000..9631f2446
--- /dev/null
+++ b/starry/_core/ops/lib/include/oblate/oblate/__init__.py
@@ -0,0 +1,1 @@
+from .oblate import *
diff --git a/starry/_core/ops/lib/include/oblate/oblate/oblate.py b/starry/_core/ops/lib/include/oblate/oblate/oblate.py
new file mode 100644
index 000000000..871e9c486
--- /dev/null
+++ b/starry/_core/ops/lib/include/oblate/oblate/oblate.py
@@ -0,0 +1,1111 @@
+from ._c_ops import sT, dsT, angles
+import numpy as np
+from scipy.integrate import quad
+import mpmath
+from mpmath import ellipe, ellipf, findroot, mpf, mpc
+from scipy.special import comb, hyp2f1
+import matplotlib.pyplot as plt
+from matplotlib.patches import Circle, Ellipse
+
+
+__all__ = [
+    "PythonSolver",
+    "NumericalSolver",
+    "CppSolver",
+    "BruteSolver",
+    "draw",
+]
+
+
+def quartic_coeffs(f, xo, yo, ro):
+    b = 1 - f
+    A = (1 - b ** 2) ** 2
+    B = -4 * xo * (1 - b ** 2)
+    C = -2 * (
+        b ** 4
+        + ro ** 2
+        - 3 * xo ** 2
+        - yo ** 2
+        - b ** 2 * (1 + ro ** 2 - xo ** 2 + yo ** 2)
+    )
+    D = -4 * xo * (b ** 2 - ro ** 2 + xo ** 2 + yo ** 2)
+    E = (
+        b ** 4
+        - 2 * b ** 2 * (ro ** 2 - xo ** 2 + yo ** 2)
+        + (ro ** 2 - xo ** 2 - yo ** 2) ** 2
+    )
+    return [A, B, C, D, E]
+
+
+def find_intersections(bo, ro, f, theta):
+    coeff = np.array(
+        quartic_coeffs(f, bo * np.sin(theta), bo * np.cos(theta), ro)
+    )
+    return np.roots(coeff), np.polyval(np.abs(coeff), np.abs(np.roots(coeff)))
+
+
+def circle_err(x, y, bo, ro, theta):
+    return np.abs(
+        (y - bo * np.cos(theta)) ** 2
+        - (ro ** 2 - (x - bo * np.sin(theta)) ** 2)
+    )
+
+
+def get_angles(bo, ro, f, theta):
+
+    # HACK to avoid instabilities
+    if f < 1e-10:
+        f = 0
+
+    # Domain adjustment
+    if bo < 0:
+        bo = -bo
+        theta -= np.pi
+
+    # HACK to avoid singular point (python impl. only!)
+    if np.abs(theta - 0.5 * np.pi) < 1e-3:
+        theta = 0.5 * np.pi + 1e-3
+    elif np.abs(theta - 1.5 * np.pi) < 1e-3:
+        theta = 1.5 * np.pi + 1e-3
+
+    # Get roots of quartic
+    roots, err = find_intersections(bo, ro, f, theta)
+
+    # Polish the roots with 30 decimal places of precision
+    try:
+        polished_roots = []
+        with mpmath.workdps(30):
+            xo = bo * mpmath.sin(theta)
+            yo = bo * mpmath.cos(theta)
+            A, B, C, D, E = quartic_coeffs(mpf(f), xo, yo, mpf(ro))
+            poly = lambda x: A * x ** 4 + B * x ** 3 + C * x ** 2 + D * x + E
+            for i, root in enumerate(roots):
+                polished_root = findroot(poly, root)
+                if np.abs(float(polished_root.imag)) < 1e-15:
+                    polished_roots += [float(polished_root.real)]
+        roots = np.array(polished_roots)
+    except:
+        pass
+
+    # Center of occultor in this frame
+    xo = bo * np.sin(theta)
+    yo = bo * np.cos(theta)
+
+    # Compute the corresponding angles
+    phi = []
+    xi = []
+    if len(roots) == 2:
+
+        # Circle intersects the edge of the ellipse twice
+        for x_root in roots:
+            x0, y0 = bo * np.sin(theta), bo * np.cos(theta)
+            y_int = (1 - f) * np.sqrt(1 - x_root ** 2)
+            phi_pos = theta + np.arctan2(y_int - y0, x_root - x0)
+            xi_pos = np.arctan2(np.sqrt(1 - x_root ** 2), x_root)
+            phi_neg = theta - np.arctan2(y_int + y0, x_root - x0)
+            xi_neg = np.arctan2(-np.sqrt(1 - x_root ** 2), x_root)
+            if circle_err(x_root, y_int, bo, ro, theta) > circle_err(
+                x_root, -y_int, bo, ro, theta
+            ):
+                phi.append(phi_neg)
+                xi.append(xi_neg)
+            else:
+                phi.append(phi_pos)
+                xi.append(xi_pos)
+
+        # Wrap the angles
+        phi = np.array(phi) % (2 * np.pi)
+        xi = np.array(xi) % (2 * np.pi)
+
+        # xi is always counter-clockwise
+        if xi[0] > xi[1]:
+            xi = xi[::-1]
+            phi = phi[::-1]
+
+        # Ensure the T integral takes us
+        # through the inside of the occultor
+        xm = np.cos(xi.mean())
+        ym = (1 - f) * np.sin(xi.mean())
+        if (xm - xo) ** 2 + (ym - yo) ** 2 >= ro ** 2:
+            xi = xi[::-1]
+            xi[1] += 2 * np.pi
+
+        # Ensure the P integral takes us through
+        # the inside of the star
+        xm = xo + ro * np.cos(theta - phi.mean())
+        ym = yo - ro * np.sin(theta - phi.mean())
+        if ym ** 2 > (1 - f) ** 2 * (1 - xm ** 2):
+            if phi[0] < phi[1]:
+                phi[0] += 2 * np.pi
+            else:
+                phi[1] += 2 * np.pi
+
+        # phi is always counter-clockwise
+        if phi[1] < phi[0]:
+            phi = phi[::-1]
+
+    elif len(roots) < 2:
+
+        # Is the center of the circle outside the ellipse?
+        if (np.abs(xo) > 1) or abs(yo) > (1 - f) * np.sqrt(1 - xo * xo):
+
+            # Is the center of the ellipse outside the circle?
+            if bo > ro:
+
+                # No occultation
+                phi = np.array([0.0, 0.0])
+                xi = np.array([0.0, 0.0])
+
+            else:
+
+                # Complete occultation
+                phi = np.array([0.0, 0.0])
+                xi = np.array([0.0, 2 * np.pi])
+
+        elif bo <= ro - 1:
+
+            # Complete occultation
+            phi = np.array([0.0, 0.0])
+            xi = np.array([0.0, 2 * np.pi])
+
+        else:
+
+            # Regular occultation, but occultor doesn't touch the limb
+            phi = np.array([0.0, 2 * np.pi])
+            xi = np.array([0.0, 0.0])
+
+    else:
+
+        raise ValueError("Unexpected branch.")
+
+    return phi, xi
+
+
+def G(n, f):
+
+    # Get the mu, nu indices
+    l = int(np.floor(np.sqrt(n)))
+    m = n - l * l - l
+    mu = l - m
+    nu = l + m
+
+    # NOTE: The abs prevents NaNs when the argument of the sqrt is
+    # zero but floating point error causes it to be ~ -eps.
+    z = lambda x, y: np.maximum(
+        1e-12, np.sqrt(np.abs(1 - x ** 2 - (y / (1 - f)) ** 2))
+    )
+
+    if nu % 2 == 0:
+
+        G = [
+            lambda x, y: 0,
+            lambda x, y: x ** (0.5 * (mu + 2)) * (y / (1 - f)) ** (0.5 * nu),
+        ]
+
+    elif (l == 1) and (m == 0):
+
+        def G0(x, y):
+            z_ = z(x, y)
+            if z_ > 1 - 1e-8:
+                return -0.5 * y
+            else:
+                return (1 - z_ ** 3) / (3 * (1 - z_ ** 2)) * (-y)
+
+        def G1(x, y):
+            z_ = z(x, y)
+            if z_ > 1 - 1e-8:
+                return 0.5 * x
+            else:
+                return (1 - z_ ** 3) / (3 * (1 - z_ ** 2)) * x
+
+        G = [G0, G1]
+
+    elif (mu == 1) and (l % 2 == 0):
+
+        G = [
+            lambda x, y: (1 - f) * x ** (l - 2) * z(x, y) ** 3,
+            lambda x, y: 0,
+        ]
+
+    elif (mu == 1) and (l % 2 != 0):
+
+        G = [lambda x, y: x ** (l - 3) * y * z(x, y) ** 3, lambda x, y: 0]
+
+    else:
+
+        G = [
+            lambda x, y: 0,
+            lambda x, y: x ** (0.5 * (mu - 3))
+            * (y / (1 - f)) ** (0.5 * (nu - 1))
+            * z(x, y) ** 3,
+        ]
+
+    return G
+
+
+def primitive(x, y, dx, dy, phi1, phi2, f, theta, n=0):
+    def func(phi):
+        Gx, Gy = G(n, f)
+        xr = lambda phi: x(phi) * np.cos(theta) + y(phi) * np.sin(theta)
+        yr = lambda phi: -x(phi) * np.sin(theta) + y(phi) * np.cos(theta)
+        dxr = lambda phi: dx(phi) * np.cos(theta) + dy(phi) * np.sin(theta)
+        dyr = lambda phi: -dx(phi) * np.sin(theta) + dy(phi) * np.cos(theta)
+        return Gx(xr(phi), yr(phi)) * dxr(phi) + Gy(xr(phi), yr(phi)) * dyr(
+            phi
+        )
+
+    res, _ = quad(func, phi1, phi2, epsabs=1e-12, epsrel=1e-12)
+    return res
+
+
+def pTbar_numerical(lmax, phi, b, r, f, theta, linear=False):
+    pT = np.zeros((lmax + 1) ** 2)
+    for n in range((lmax + 1) ** 2):
+        for phi1, phi2 in np.reshape(phi, (-1, 2)):
+            l = int(np.floor(np.sqrt(n)))
+            m = n - l * l - l
+            mu = l - m
+            nu = l + m
+            if (not linear) or (nu % 2 == 0) or (l == 1 and m == 0):
+                # Exact integral
+                x = lambda phi: r * np.cos(phi)
+                y = lambda phi: (b + r * np.sin(phi))
+                dx = lambda phi: -r * np.sin(phi)
+                dy = lambda phi: r * np.cos(phi)
+                pT[n] += primitive(x, y, dx, dy, phi1, phi2, f, theta, n)
+            else:
+                # Linear approximation in `f`
+                x = lambda phi: r * np.cos(theta - phi) + b * np.sin(theta)
+                y = lambda phi: b * np.cos(theta) - r * np.sin(theta - phi)
+                z = lambda phi: np.sqrt(np.abs(1 - x(phi) ** 2 - y(phi) ** 2))
+                dz3df = lambda phi: -3 * y(phi) ** 2 * z(phi)
+                if (mu == 1) and (l % 2 == 0):
+                    func = lambda phi: (
+                        r
+                        * (1 - f)
+                        * x(phi) ** (l - 2)
+                        * (z(phi) ** 3 + dz3df(phi) * f)
+                        * np.sin(theta - phi)
+                    )
+                elif (mu == 1) and (l % 2 != 0):
+                    func = lambda phi: (
+                        r
+                        * x(phi) ** (l - 3)
+                        * y(phi)
+                        * (z(phi) ** 3 + dz3df(phi) * f)
+                        * np.sin(theta - phi)
+                    )
+                else:
+                    func = lambda phi: (
+                        r
+                        * (1 - f) ** ((1 - nu) // 2)
+                        * x(phi) ** ((mu - 3) // 2)
+                        * y(phi) ** ((nu - 1) // 2)
+                        * (z(phi) ** 3 + dz3df(phi) * f)
+                        * np.cos(theta - phi)
+                    )
+                res, _ = quad(func, phi1, phi2, epsabs=1e-12, epsrel=1e-12)
+                pT[n] += res
+
+    return pT
+
+
+def tTbar_numerical(lmax, xi, b, r, f, theta):
+    tT = np.zeros((lmax + 1) ** 2)
+    for n in range((lmax + 1) ** 2):
+        for xi1, xi2 in np.reshape(xi, (-1, 2)):
+            x = lambda xi: np.cos(theta) * np.cos(xi) - (1 - f) * np.sin(
+                theta
+            ) * np.sin(xi)
+            y = lambda xi: np.sin(theta) * np.cos(xi) + (1 - f) * np.cos(
+                theta
+            ) * np.sin(xi)
+            dx = lambda xi: -np.cos(theta) * np.sin(xi) - (1 - f) * np.sin(
+                theta
+            ) * np.cos(xi)
+            dy = lambda xi: -np.sin(theta) * np.sin(xi) + (1 - f) * np.cos(
+                theta
+            ) * np.cos(xi)
+            tT[n] += primitive(x, y, dx, dy, xi1, xi2, f, theta, n)
+    return tT
+
+
+def sTbar_numerical(lmax, b, r, f, theta, linear=False):
+    phi, xi = get_angles(b, r, f, theta)
+    return pTbar_numerical(
+        lmax, phi, b, r, f, theta, linear=linear
+    ) + tTbar_numerical(lmax, xi, b, r, f, theta)
+
+
+def integrate(f, a, b, epsabs=1e-12, epsrel=1e-12, **kwargs):
+    return quad(f, a, b, epsabs=epsabs, epsrel=epsrel, **kwargs)[0]
+
+
+def get_sT0(lmax):
+    # Compute the solution vector for no occultation
+    L = np.zeros((lmax + 3, lmax + 3))
+    L[0, 0] = 2 * np.pi
+    for u in range(0, lmax + 3, 2):
+        for v in range(2, lmax + 3, 2):
+            fac = (v - 1.0) / (u + v)
+            L[u, v] = fac * L[u, v - 2]
+            L[v, u] = fac * L[v - 2, u]
+    sT0 = np.zeros((lmax + 1) ** 2)
+    n = 0
+    for l in range(lmax + 1):
+        for m in range(-l, l + 1):
+            mu = l - m
+            nu = l + m
+            if nu % 2 == 0:
+                # Case 1
+                sT0[n] = L[mu // 2 + 2, nu // 2]
+            elif (l == 1) and (m == 0):
+                # Case 2
+                sT0[n] = 2 * np.pi / 3
+            n += 1
+    return sT0
+
+
+class PythonSolver:
+    def __init__(self, lmax=5):
+
+        # Pre-compute case conditions
+        self.lmax = lmax
+        self.cases = np.zeros((self.lmax + 1) ** 2, dtype=int)
+        self.l = np.zeros((self.lmax + 1) ** 2, dtype=int)
+        self.m = np.zeros((self.lmax + 1) ** 2, dtype=int)
+        for n in range((self.lmax + 1) ** 2):
+            self.cases[n] = self.get_case(n)
+            self.l[n] = int(np.floor(np.sqrt(n)))
+            self.m[n] = n - self.l[n] ** 2 - self.l[n]
+        self.sT0 = get_sT0(self.lmax)
+
+    def get_case(self, n):
+        l = int(np.floor(np.sqrt(n)))
+        m = n - l ** 2 - l
+        mu = l - m
+        nu = l + m
+        if (nu % 2) == 0:
+            return 1
+        elif (l == 1) and (m == 0):
+            return 2
+        elif ((nu % 2) == 1) and (mu == 1) and ((l % 2) == 0):
+            return 3
+        elif ((nu % 2) == 1) and (mu == 1) and ((l % 2) == 1):
+            return 4
+        else:
+            return 5
+
+    def solve(self, f0, fN, A, B, C, N):
+        """
+        Solve a recursion relation of the form
+        
+            f(j) = A(j) * f(j - 1) + B(j) * f (j - 2) + C(j)
+        
+        given a lower boundary condition f(0) and 
+        an upper boundary condition f(N).
+        
+        """
+        # Set up the tridiagonal problem
+        a = np.empty(N - 1)
+        b = np.empty(N - 1)
+        c = np.empty(N - 1)
+        for i, j in enumerate(range(2, N + 1)):
+            a[i] = -A(j)
+            b[i] = -B(j)
+            c[i] = C(j)
+
+        # Add the boundary conditions
+        c[0] -= b[0] * f0
+        c[-1] -= fN
+
+        # Construct the tridiagonal matrix
+        A = np.diag(a, 0) + np.diag(b[1:], -1) + np.diag(np.ones(N - 2), 1)
+
+        # Solve using a dense solver for stability
+        soln = np.linalg.solve(A, c)
+        return np.concatenate(([f0], soln, [fN]))
+
+    def get_W(self, k2, phi1, phi2):
+        """
+        Return the vector of `W` integrals, computed
+        recursively given a lower boundary condition
+        (analytic in terms of elliptic integrals) and an upper
+        boundary condition (computed numerically).
+
+        The term `W_j` is the solution to the integral of
+        
+            sin(u)^(2 * j) * sqrt(1 - sin(u)^2 / (1 - k^2))
+
+        from u = u1 to u = u2, where 
+        
+            u = (pi - 2 * phi) / 4
+
+        """
+        # Useful quantities
+        N = 2 * self.lmax + 4
+        kc2 = 1 - k2
+        u1 = 0.25 * (np.pi - 2 * phi1)
+        u2 = 0.25 * (np.pi - 2 * phi2)
+        u = np.array([u1, u2])
+        sinu = np.sin(u)
+        cosu = np.cos(u)
+        diff = lambda x: np.diff(x)[0]
+        D2 = 1 - sinu ** 2 / kc2
+        if k2 >= 1:
+            sgn = 1
+        else:
+            sgn = -1
+        D3 = sgn * abs(D2) ** 1.5
+
+        # The two boundary conditions
+        if sgn == 1:
+            f0 = float((ellipe(u[1], 1 / kc2) - ellipe(u[0], 1 / kc2)).real)
+        else:
+            f0 = float((ellipe(u[1], 1 / kc2) - ellipe(u[0], 1 / kc2)).imag)
+            """
+            For reference, this expression can be written in terms of the
+            real parts of elliptic integrals as follows:
+
+                n = np.floor((u + 0.5 * np.pi) / np.pi)
+                term = np.sqrt((sinu ** 2 / kc2 - 1) * (1 / k2 - 1))
+                A = n * np.pi + (-1) ** n * np.arcsin(term / sinu)
+                F = float((ellipf(A[1], k2) - ellipf(A[0], k2)).real)
+                E = float((ellipe(A[1], k2) - ellipe(A[0], k2)).real)
+                C = diff(
+                    (np.sin(A) * np.cos(A)) / np.sqrt(1 - k2 * np.sin(A) ** 2)
+                )
+                f0 = -np.sqrt(kc2) * (F - E / kc2 + k2 / kc2 * C)
+            """
+
+        fN = quad(
+            lambda u: np.sin(u) ** (2 * N)
+            * abs(1 - np.sin(u) ** 2 / kc2) ** 0.5,
+            u1,
+            u2,
+            epsabs=1e-12,
+            epsrel=1e-12,
+        )[0]
+
+        # The recursion coefficients
+        A = lambda j: 2 * (j + (j - 1) * kc2) / (2 * j + 1)
+        B = lambda j: -(2 * j - 3) / (2 * j + 1) * kc2
+        C = lambda j: diff(kc2 * sinu ** (2 * j - 3) * cosu * D3) / (2 * j + 1)
+
+        # Solve the linear system
+        return self.solve(f0, fN, A, B, C, N)
+
+    def get_V(self, k2, phi1, phi2):
+        """
+        Return the vector of `V` integrals, computed
+        recursively given a lower boundary condition
+        (trivial) and an upper boundary condition 
+        (computed from `2F1`). 
+        
+        The term `V_i` is the solution to the integral of
+        
+            cos(phi) sin(phi)^i sqrt(1 - w^2 sin(phi))
+        
+        from phi=phi1 to phi=phi2.
+        
+        """
+        N = 2 * self.lmax + 4
+        invw2 = 2 * k2 - 1
+        if k2 < 0.5:
+            sgn = -1
+        else:
+            sgn = 1
+
+        V = np.zeros((N + 1, 2))
+        for i, phi in enumerate([phi1, phi2]):
+
+            # Useful quantities
+            sinphi = np.sin(phi)
+            x = sinphi / invw2
+
+            # The two boundary conditions
+            if sgn == 1:
+                f0 = (2.0 / 3.0) * (1 - (1 - x) ** 1.5) * invw2
+                fN = (
+                    sinphi ** (N + 1) / (N + 1) * hyp2f1(-0.5, N + 1, N + 2, x)
+                )
+            else:
+                # NOTE: These are the original expressions:
+                #
+                # f0 = ((2.0 / 3.0) * (1 - (1 - x + 0j) ** 1.5) * invw2).imag
+                # fN = -(
+                #    sinphi ** (N + 1)
+                #    / (N + 1)
+                #    * hyp2f1(-0.5, N + 1, N + 2, x + 0j)
+                # ).imag
+                #
+                term = (2.0 / 3.0) * (x - 1) ** 1.5
+                f0 = term * invw2
+                # This is equivalent:
+                # fN = f0 * invw2 ** N * hyp2f1(1.5, -N, 2.5, 1 - x)
+                fN = (
+                    term * sinphi ** (N + 1) * hyp2f1(1.0, N + 2.5, 2.5, 1 - x)
+                )
+
+            # The recursion coefficients
+            A = lambda j: (2 * j + (2 * j + 1) * x) * invw2 / (2 * j + 3)
+            B = lambda j: -2 * (j - 1) * sinphi * invw2 / (2 * j + 3)
+            C = lambda j: 0.0
+
+            # Solve the linear system
+            V[:, i] = self.solve(f0, fN, A, B, C, N)
+
+        return np.diff(V, axis=1)
+
+    def get_J(self, k2, phi1, phi2):
+        """
+        Return the `J` matrix of integrals. The term `J_{p,q}`
+        is equal to the integral of
+        
+            cos(phi)^p sin(phi)^q sqrt(1 - w^2 sin(phi))
+            
+        from phi=phi1 to phi=phi2.
+        
+        """
+        W = self.get_W(k2, phi1, phi2)
+        V = self.get_V(k2, phi1, phi2)
+        w2 = 1 / (2 * k2 - 1)
+        c1 = -2.0 * abs(1 - w2) ** 0.5
+        J = np.zeros((self.lmax + 3, self.lmax + 3))
+        for s in range((self.lmax + 3) // 2):
+            for q in range(self.lmax + 3):
+                fac = 1
+                for i in range(s + 1):
+                    term = 0
+                    amp = 1
+                    for j in range(2 * i + q + 1):
+                        term += amp * W[j]
+                        amp *= -2 * (2 * i + q - j) / (j + 1)
+                    J[2 * s, q] += c1 * fac * term
+                    J[2 * s + 1, q] += fac * V[2 * i + q]
+                    fac *= (i - s) / (i + 1)
+        return J
+
+    def get_L(self, phip1, phip2):
+        """
+        Return the `L` matrix of integrals. The term `L_{i, j}`
+        is equal to the integral of
+        
+            cos(phi)^i sin(phi)^j
+            
+        from phi=phip1 to phi=phip2, where
+        
+            phip = phi - theta
+            
+        (pT integral) or
+        
+            phip = xi
+            
+        (tT integral).
+        
+        """
+        # Lower boundary
+        L = np.zeros((self.lmax + 3, self.lmax + 3))
+        L[0, 0] = phip2 - phip1
+        L[1, 0] = np.sin(phip2) - np.sin(phip1)
+        L[0, 1] = np.cos(phip1) - np.cos(phip2)
+        L[1, 1] = 0.5 * (np.cos(phip1) ** 2 - np.cos(phip2) ** 2)
+
+        # Recursion coeffs
+        cp1 = np.cos(phip1)
+        cp2 = np.cos(phip2)
+        sp1 = np.sin(phip1)
+        sp2 = np.sin(phip2)
+        A0 = cp1 * sp1
+        B0 = cp2 * sp2
+        C0 = cp2 * sp2
+        D0 = cp1 * sp1
+
+        # Recurse
+        for u in range(self.lmax + 3):
+            A = A0
+            B = B0
+            C = C0
+            D = D0
+            for v in range(2, self.lmax + 3):
+                fac = 1.0 / (u + v)
+                L[u, v] = fac * (A - B + (v - 1) * L[u, v - 2])
+                L[v, u] = fac * (C - D + (v - 1) * L[v - 2, u])
+                A *= sp1
+                B *= sp2
+                C *= cp2
+                D *= cp1
+            A0 *= cp1
+            B0 *= cp2
+            C0 *= sp2
+            D0 *= sp1
+
+        return L
+
+    def get_pT2(self, r, theta, f, b, phi1, phi2):
+        """
+        Return the `pT[2]` integral, computed by numerical integration.
+        
+        """
+        bc = b * np.cos(theta)
+        bs = b * np.sin(theta)
+        p = (1 - f) ** (-2)
+
+        def func(phi):
+            z = np.sqrt(
+                1
+                - (r * np.cos(phi - theta) + bs) ** 2
+                - p * (r * np.sin(phi - theta) + bc) ** 2
+            )
+            z = max(1e-12, min(z, 1 - 1e-8))
+            return (
+                r * (r + b * np.sin(phi)) * (1 - z ** 3) / (3 * (1 - z ** 2))
+            )
+
+        res, _ = quad(func, phi1, phi2, epsabs=1e-12, epsrel=1e-12)
+        return res
+
+    def _get_sTbar(
+        self, bo=0.58, ro=0.4, f=0.2, theta=0.5,
+    ):
+        """
+        Return the complement of the `sT` solution vector.
+
+        """
+        # Stability hacks
+        costheta = np.cos(theta)
+        sintheta = np.sin(theta)
+        if np.abs(bo) < 1e-15:
+            bo = 1e-15
+        if np.abs(sintheta) < 1e-15:
+            sintheta = 1e-15
+        elif np.abs(costheta) < 1e-15:
+            costheta = 1e-15
+        if np.abs(bo - ro) < 1e-12:
+            bo = ro + 1e-12
+        if np.abs(bo - ro) < 1e-5 and np.abs(ro - 0.5) < 1e-5:
+            bo = ro + 1e-5
+        if np.abs(1 - bo - ro) < 1e-6:
+            bo = 1 - ro + 1e-6
+
+        # Get the angles
+        (phi1, phi2), (xi1, xi2) = get_angles(bo, ro, f, theta)
+
+        # Special cases
+        if (phi1 == phi2 == xi1 == 0) and (xi2 == 2 * np.pi):
+            # No occultation?
+            sT = np.zeros((self.lmax + 1) ** 2)
+            Lxi = (1 - f) * self.get_L(xi1, xi2)
+            for n in range((self.lmax + 1) ** 2):
+                l = self.l[n]
+                m = self.m[n]
+                mu = l - m
+                nu = l + m
+                if self.cases[n] == 1:
+                    tT = Lxi[mu // 2 + 2, nu // 2]
+                    sT[n] = tT
+                elif self.cases[n] == 2:
+                    tT = (1 - f) * (xi2 - xi1) / 3
+                    sT[n] = tT
+            return sT
+        elif phi1 == phi2 == xi1 == xi2 == 0:
+            # Complete occultation
+            return np.zeros((self.lmax + 1) ** 2)
+
+        # Useful variables
+        b = bo
+        r = ro
+        tantheta = sintheta / costheta
+        invtantheta = 1.0 / tantheta
+        gamma = 1 - b ** 2 - r ** 2
+        sqrtgamma = np.sqrt(np.abs(gamma))
+        w2 = 2 * b * r / gamma
+        k2 = (gamma + 2 * b * r) / (4 * b * r)
+
+        # Compute the S, C matrices
+        S = np.zeros((self.lmax + 3, self.lmax + 3))
+        C = np.zeros((self.lmax + 3, self.lmax + 3))
+        fac0 = 1
+        for i in range(self.lmax + 3):
+            facs = r * fac0
+            facc = fac0
+            for j in range(i + 1):
+                S[i, self.lmax + 2 + j - i] = facs
+                C[self.lmax + 2 + j - i, i] = facc
+                fac = b * (i - j) / (r * (j + 1))
+                facs *= fac * sintheta
+                facc *= fac * costheta
+            fac0 *= r
+
+        # Compute the B matrices
+        A = np.empty((self.lmax + 3, self.lmax + 3))
+        B = np.empty((self.lmax + 3, self.lmax + 3))
+        fac0l = 1
+        fac0r = 1
+        for i in range(self.lmax + 3):
+            facl = fac0l
+            facr = fac0r
+            for j in range(self.lmax + 3):
+                A[i, j] = facl
+                B[i, j] = facr
+                fac = (i - j) / (j + 1)
+                facl *= fac * tantheta
+                facr *= -fac * invtantheta
+            fac0l *= costheta
+            fac0r *= -sintheta
+
+        # Compute the first f derivative matrix
+        D = -3 * np.array(
+            [
+                [
+                    0.5
+                    * (
+                        b ** 2 + r ** 2 + (b + r) * (b - r) * np.cos(2 * theta)
+                    ),
+                    2 * b * r * costheta ** 2,
+                    r ** 2 * np.cos(2 * theta),
+                ],
+                [
+                    -2 * b * r * costheta * sintheta,
+                    -2 * r ** 2 * costheta * sintheta,
+                    0,
+                ],
+            ]
+        )
+
+        # Compute the L and M integrals
+        L = self.get_L(phi1 - theta, phi2 - theta)
+        Lxi = (1 - f) * self.get_L(xi1, xi2)
+        M = S[:-1, 1:] @ L[::-1, ::-1][:-1, :] @ C[:, :-1]
+
+        # Compute the J integrals
+        J = self.get_J(k2, phi1, phi2)
+        J32 = gamma * sqrtgamma * (J[:-1, :-1] - w2 * J[:-1, 1:])
+        J12 = sqrtgamma * J[:-1, :-1]
+
+        # Compute the I integral
+        I = np.zeros((self.lmax + 2, self.lmax + 2))
+        T = np.zeros((self.lmax + 1, self.lmax))
+        for i in range(self.lmax + 1):
+            for j in range(min(self.lmax, i + 1)):
+
+                # Taylor series in `f`
+                J0 = J32[i - j, j]
+                J1 = np.sum(D * (J12[i - j :, j:])[:2, :3])
+                T[i, j] = J0 + f * J1
+
+                for k in range(i + 1):
+                    imk = i - k
+                    for l in range(max(0, j - imk), min(k, j) + 1):
+                        I[self.lmax + 1 - imk, self.lmax + 1 - k] += (
+                            A[imk, j - l] * T[i, j] * B[k, l]
+                        )
+
+        # Compute the K, H integrals
+        K = S[:-1, 2:] @ I[:-1, :] @ C[1:, :-1]
+        H = S[:-1, 1:] @ I[:, :-1] @ C[2:, :-1]
+
+        # Compute sT on a case-by-case basis
+        sT = np.zeros((self.lmax + 1) ** 2)
+        for n in range((self.lmax + 1) ** 2):
+            l = self.l[n]
+            m = self.m[n]
+            mu = l - m
+            nu = l + m
+            if self.cases[n] == 1:
+                pT = (1 - f) ** (-nu // 2) * M[(mu + 2) // 2, nu // 2]
+                tT = Lxi[mu // 2 + 2, nu // 2]
+                sT[n] = pT + tT
+            elif self.cases[n] == 2:
+                pT = self.get_pT2(r, theta, f, b, phi1, phi2)
+                tT = (1 - f) * (xi2 - xi1) / 3
+                sT[n] = pT + tT
+            elif self.cases[n] == 3:
+                sT[n] = -(1 - f) * H[l - 2, 0]
+            elif self.cases[n] == 4:
+                sT[n] = -H[l - 3, 1]
+            else:
+                sT[n] = (1 - f) ** ((1 - nu) // 2) * K[
+                    (mu - 3) // 2, (nu - 1) // 2
+                ]
+        return sT
+
+    def get_sT(self, nruns=1, f=0.2, **kwargs):
+        for n in range(nruns):
+            res = (1 - f) * self.sT0 - self._get_sTbar(f=f, **kwargs)
+        return res
+
+
+class NumericalSolver:
+    def __init__(self, lmax=5, linear=False):
+        self.lmax = lmax
+        self.linear = linear
+        self.sT0 = get_sT0(self.lmax)
+
+    def get_sT(self, bo=0.58, ro=0.4, f=0.2, theta=0.5, nruns=1):
+        for n in range(nruns):
+            res = (1 - f) * self.sT0 - sTbar_numerical(
+                self.lmax, bo, ro, f, theta, linear=self.linear
+            )
+        return res
+
+
+class CppSolver:
+    def __init__(self, lmax=5):
+        self.lmax = lmax
+
+    def get_sT(
+        self, bo=0.58, ro=0.4, f=0.2, theta=0.5, nruns=1,
+    ):
+        return sT(self.lmax, bo, ro, f, theta, nruns)
+
+
+class BruteSolver:
+    def __init__(self, lmax=5, res=999, grid_star=False):
+        self.lmax = lmax
+        self.res = res
+        self.grid_star = grid_star
+        self.sT0 = get_sT0(self.lmax)
+
+    def g(self, n, x, y, z=None):
+        """
+        Return the nth term of the Green's basis.
+        
+        """
+        if z is None:
+            z2 = 1 - x ** 2 - y ** 2
+            z = np.sqrt(np.abs(z2))
+            on_star = z2 > 0
+            z[~on_star] = np.nan
+        l = int(np.floor(np.sqrt(n)))
+        m = n - l * l - l
+        mu = l - m
+        nu = l + m
+        if nu % 2 == 0:
+            I = [mu // 2]
+            J = [nu // 2]
+            K = [0]
+            C = [(mu + 2) // 2]
+        elif (l == 1) and (m == 0):
+            I = [0]
+            J = [0]
+            K = [1]
+            C = [1]
+        elif (mu == 1) and (l % 2 == 0):
+            I = [l - 2]
+            J = [1]
+            K = [1]
+            C = [3]
+        elif mu == 1:
+            I = [l - 3, l - 1, l - 3]
+            J = [0, 0, 2]
+            K = [1, 1, 1]
+            C = [-1, 1, 4]
+        else:
+            I = [(mu - 5) // 2, (mu - 5) // 2, (mu - 1) // 2]
+            J = [(nu - 1) // 2, (nu + 3) // 2, (nu - 1) // 2]
+            K = [1, 1, 1]
+            C = [(mu - 3) // 2, -(mu - 3) // 2, -(mu + 3) // 2]
+        res = z * 0
+        for i, j, k, c in zip(I, J, K, C):
+            if c != 0:
+                res += c * x ** i * y ** j * z ** k
+
+        return res
+
+    def get_sT(
+        self, bo=0.58, ro=0.4, f=0.2, theta=0.5, nruns=1,
+    ):
+        if self.grid_star:
+
+            # Grid up the star (accurate for large `ro`)
+
+            for n in range(nruns):
+
+                # Semi-minor axis of the ellipse
+                b = 1 - f
+
+                # There are two equivalent ways of doing this:
+                if False:
+
+                    # Grid up the ellipse
+                    x = np.linspace(-1 - ro, 1 + ro, self.res)
+                    y = np.linspace(-b - ro, b + ro, self.res)
+                    x, y = np.meshgrid(x, y)
+                    on_star = 1 - x ** 2 - (y / b) ** 2 > 0
+
+                    # Compute the Green's basis in (x, y') where
+                    #
+                    #     y' = y / b
+                    #
+                    # is the transformed y coordinate
+                    g = np.array(
+                        [
+                            self.g(n, x, y / b)
+                            for n in range((self.lmax + 1) ** 2)
+                        ]
+                    )
+                    g /= np.count_nonzero(on_star)
+
+                    # Compute the visible pixels
+                    xo = bo * np.sin(theta)
+                    yo = bo * np.cos(theta)
+                    under_occultor = (x - xo) ** 2 + (y - yo) ** 2 <= ro ** 2
+                    inds = on_star & ~under_occultor
+
+                else:
+
+                    # Grid up the unit disk
+                    x = np.linspace(-1 - ro, 1 + ro, self.res)
+                    y = np.linspace(-1 - ro, 1 + ro, self.res)
+                    x, y = np.meshgrid(x, y)
+                    on_star = 1 - x ** 2 - y ** 2 > 0
+
+                    # Compute the Green's basis as usual
+                    g = np.array(
+                        [self.g(n, x, y) for n in range((self.lmax + 1) ** 2)]
+                    )
+                    g /= np.count_nonzero(on_star)
+
+                    # Compute the visible pixels in a stretched frame where
+                    #
+                    #     y' = y * b
+                    #
+                    # is the transformed y coordinate
+                    xo = bo * np.sin(theta)
+                    yo = bo * np.cos(theta)
+                    under_occultor = (x - xo) ** 2 + (
+                        y * b - yo
+                    ) ** 2 <= ro ** 2
+                    inds = on_star & ~under_occultor
+
+            # Integrate
+            sT = b * np.pi * np.sum(g[:, inds], axis=1)
+
+            return sT
+
+        else:
+
+            # Grid up the planet and subtract the integral from
+            # the solution over the entire disk (accurate for small `ro`)
+
+            for n in range(nruns):
+
+                # Semi-minor axis of the ellipse
+                b = 1 - f
+
+                # Location of the occultor
+                xo = bo * np.sin(theta)
+                yo = bo * np.cos(theta)
+
+                # Grid up the occultor
+                x = np.linspace(xo - ro, xo + ro, self.res)
+                y = np.linspace(yo - ro, yo + ro, self.res)
+                x, y = np.meshgrid(x, y)
+                on_star = 1 - x ** 2 - (y / b) ** 2 > 0
+
+                # Compute the Green's basis in (x, y') where
+                #
+                #     y' = y / b
+                #
+                # is the transformed y coordinate
+                g = np.array(
+                    [self.g(n, x, y / b) for n in range((self.lmax + 1) ** 2)]
+                )
+
+                # Compute the masked pixels
+                under_occultor = (x - xo) ** 2 + (y - yo) ** 2 <= ro ** 2
+                inds = on_star & under_occultor
+                g /= np.count_nonzero(under_occultor)
+
+            # The solution vector for the complement of the visible region
+            sTbar = np.pi * ro ** 2 * np.sum(g[:, inds], axis=1)
+
+            return b * self.sT0 - sTbar
+
+
+def draw(bo=0.58, ro=0.4, f=0.2, theta=0.5, file=None):
+
+    # Domain adjustment
+    if bo < 0:
+        bo = -bo
+        theta -= np.pi
+
+    # Set up the figure
+    fig, ax = plt.subplots(1, figsize=(8, 8))
+    ax.set_xlim(-1.01 - 2 * ro, 1.01 + 2 * ro)
+    ax.set_ylim(-1.01 - 2 * ro, 1.01 + 2 * ro)
+    ax.set_aspect(1)
+    ax.axis("off")
+
+    # Draw the two bodies
+    occulted = Ellipse((0, 0), 1.0 * 2, (1 - f) * 2, 0, fill=False, color="k")
+    ax.add_artist(occulted)
+    occulted_max = Circle(
+        (0, 0), 1.0, fill=False, color="k", alpha=0.1, ls="--", lw=1
+    )
+    ax.add_artist(occulted_max)
+    occulted_min = Circle(
+        (0, 0), 1.0 - f, fill=False, color="k", alpha=0.1, ls="--", lw=1
+    )
+    ax.add_artist(occulted_min)
+    occultor = Circle(
+        (bo * np.sin(theta), bo * np.cos(theta)), ro, fill=False, color="r"
+    )
+    ax.add_artist(occultor)
+    ax.plot(0, 0, "ko", ms=3)
+    ax.plot(bo * np.sin(theta), bo * np.cos(theta), "ro", ms=3)
+
+    # Get angles
+    phi, xi = get_angles(bo, ro, f, theta)
+    ax.set_title(
+        r"$\phi: [{:.0f}^\circ, {:.0f}^\circ] \,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \xi: [{:.0f}^\circ, {:.0f}^\circ]$".format(
+            *(phi * 180 / np.pi), *(xi * 180 / np.pi)
+        ),
+        fontsize=12,
+    )
+
+    # Integration path along occultor (white --> red)
+    x = np.zeros(10000)
+    y = np.zeros(10000)
+    for k, phi_k in enumerate(np.linspace(phi[0], phi[1], 10000)):
+        x0 = ro * np.cos(phi_k)
+        y0 = bo + ro * np.sin(phi_k)
+        x[k] = x0 * np.cos(theta) + y0 * np.sin(theta)
+        y[k] = -x0 * np.sin(theta) + y0 * np.cos(theta)
+    plt.plot(x, y, "r", zorder=99)
+    plt.plot(x[0], y[0], "ro", zorder=100, ms=5)
+    plt.plot(x[-1], y[-1], "ro", zorder=100, ms=5)
+    for i in [len(x) // 4, len(x) // 2, (3 * len(x)) // 4, len(x) - 2]:
+        plt.annotate(
+            "",
+            xytext=(x[i], y[i]),
+            xy=(x[i + 1], y[i + 1]),
+            arrowprops=dict(arrowstyle="->", color="r"),
+            size=15,
+        )
+
+    # Integration path along star (white --> black)
+    x = np.zeros(10000)
+    y = np.zeros(10000)
+    for k, xi_k in enumerate(np.linspace(xi[0], xi[1], 10000)):
+        x[k] = np.cos(xi_k)
+        y[k] = (1 - f) * np.sin(xi_k)
+    plt.plot(x, y, "k", zorder=99)
+    plt.plot(x[0], y[0], "ko", zorder=101, ms=3)
+    plt.plot(x[-1], y[-1], "ko", zorder=101, ms=3)
+    for i in [len(x) // 4, len(x) // 2, (3 * len(x)) // 4, len(x) - 2]:
+        plt.annotate(
+            "",
+            xytext=(x[i], y[i]),
+            xy=(x[i + 1], y[i + 1]),
+            arrowprops=dict(arrowstyle="->", color="k"),
+            size=15,
+        )
+
+    # Save or display
+    if file is not None:
+        fig.savefig(file, bbox_inches="tight")
+        plt.close()
+    else:
+        plt.show()
diff --git a/starry/_core/ops/lib/include/oblate/occultation.h b/starry/_core/ops/lib/include/oblate/occultation.h
new file mode 100644
index 000000000..8a3c843ae
--- /dev/null
+++ b/starry/_core/ops/lib/include/oblate/occultation.h
@@ -0,0 +1,343 @@
+
+/**
+\file occultation.h
+\brief Solver for occultations of oblate bodies.
+
+*/
+
+#ifndef _STARRY_OBLATE_OCCULTATION_H_
+#define _STARRY_OBLATE_OCCULTATION_H_
+
+#include "../quad.h"
+#include "../utils.h"
+#include "geometry.h"
+#include "numerical.h"
+
+namespace starry {
+namespace oblate {
+namespace occultation {
+
+using namespace utils;
+using namespace numerical;
+using namespace geometry;
+using std::min, std::max;
+
+template <class Scalar, int N> class Occultation {
+
+  using A = ADScalar<Scalar, N>;
+
+protected:
+  // Misc
+  int deg;
+  int ncoeff;
+
+  // Inputs
+  A bo;
+  A ro;
+  A f;
+  A theta;
+  A costheta;
+  A sintheta;
+
+  // Integrals
+  Matrix<A> M;
+  Matrix<A> Lp;
+  Matrix<A> Lt;
+  Vector<A> K;
+
+  // Helper matrices
+  Matrix<A> S;
+  Matrix<A> C;
+
+  // Numerical integration
+  Quad<Scalar, N> QUAD;
+
+  /**
+   *
+   * Compute the matrix of `L` integrals.
+   *
+   * The term `L_{i,j}` is the solution to the integral of
+   *
+   *    cos(phi)^i * sin(phi)^j
+   *
+   * from phi = phip1 to phi = phip2, where
+   *
+   *    phip = phi - theta
+   *
+   * (pT integral) or
+   *
+   *    phip = xi
+   *
+   * (tT integral).
+   *
+   */
+  inline void compute_L(const A &phip1, const A &phip2, Matrix<A> &L) {
+
+    //
+    int nmax = deg + 3;
+    A cp1 = cos(phip1);
+    A cp2 = cos(phip2);
+    A sp1 = sin(phip1);
+    A sp2 = sin(phip2);
+    L.resize(nmax, nmax);
+
+    // Lower boundary
+    L(0, 0) = phip2 - phip1;
+    L(1, 0) = sp2 - sp1;
+    L(0, 1) = cp1 - cp2;
+    L(1, 1) = 0.5 * (cp1 * cp1 - cp2 * cp2);
+
+    // Recursion coeffs
+    A fac, A1, B1, C1, D1;
+    A A0 = cp1 * sp1;
+    A B0 = cp2 * sp2;
+    A C0 = cp2 * sp2;
+    A D0 = cp1 * sp1;
+
+    // Recurse
+    for (int u = 0; u < nmax; ++u) {
+      A1 = A0;
+      B1 = B0;
+      C1 = C0;
+      D1 = D0;
+      for (int v = 2; v < nmax; ++v) {
+        fac = 1.0 / (u + v);
+        L(u, v) = fac * (A1 - B1 + (v - 1) * L(u, v - 2));
+        L(v, u) = fac * (C1 - D1 + (v - 1) * L(v - 2, u));
+        A1 *= sp1;
+        B1 *= sp2;
+        C1 *= cp2;
+        D1 *= cp1;
+      }
+      A0 *= cp1;
+      B0 *= cp2;
+      C0 *= sp2;
+      D0 *= sp1;
+    }
+  }
+
+  /**
+   *
+   * Compute the matrix of `L` integrals for the special case
+   * phi1 = 0, phi2 = 2 * pi (no occultation).
+   *
+   */
+  inline void compute_L0(Matrix<A> &L) {
+
+    //
+    int nmax = deg + 3;
+    L.setZero(nmax, nmax);
+
+    // Lower boundary
+    L(0, 0) = 2 * pi<Scalar>();
+
+    // Recurse
+    A fac;
+    for (int u = 0; u < nmax; u += 2) {
+      for (int v = 2; v < nmax; v += 2) {
+        fac = (v - 1.0) / (u + v);
+        L(u, v) = fac * L(u, v - 2);
+        L(v, u) = fac * L(v - 2, u);
+      }
+    }
+  }
+
+  /**
+   * Nudge the inputs away from singular points
+   */
+  inline void nudge_inputs() {
+    if (abs(bo - ro) < STARRY_BO_EQUALS_RO_TOL)
+      bo = ro + (bo > ro ? STARRY_BO_EQUALS_RO_TOL : -STARRY_BO_EQUALS_RO_TOL);
+    if ((abs(bo - ro) < STARRY_BO_EQUALS_RO_EQUALS_HALF_TOL) &&
+        (abs(ro - 0.5) < STARRY_BO_EQUALS_RO_EQUALS_HALF_TOL))
+      bo = ro + STARRY_BO_EQUALS_RO_EQUALS_HALF_TOL;
+    if (abs(bo) < STARRY_BO_EQUALS_ZERO_TOL)
+      bo = STARRY_BO_EQUALS_ZERO_TOL;
+    if ((ro > 0) && (ro < STARRY_RO_EQUALS_ZERO_TOL))
+      ro = STARRY_RO_EQUALS_ZERO_TOL;
+    if (abs(1 - bo - ro) < STARRY_BO_EQUALS_ONE_MINUS_RO_TOL)
+      bo = 1 - ro + STARRY_BO_EQUALS_ONE_MINUS_RO_TOL;
+    if (abs(theta - 0.5 * pi<Scalar>()) < STARRY_ROOT_TOL_THETA_PI_TWO) {
+      theta += (theta > 0.5 * pi<Scalar>() ? 1.0 : -1.0) *
+               STARRY_ROOT_TOL_THETA_PI_TWO;
+      costheta = cos(theta);
+      sintheta = sin(theta);
+    } else if (abs(theta + 0.5 * pi<Scalar>()) < STARRY_ROOT_TOL_THETA_PI_TWO) {
+      theta += (theta > -0.5 * pi<Scalar>() ? 1.0 : -1.0) *
+               STARRY_ROOT_TOL_THETA_PI_TWO;
+      costheta = cos(theta);
+      sintheta = sin(theta);
+    } else if (abs(sintheta) < STARRY_T_TOL) {
+      sintheta = sintheta > 0 ? STARRY_T_TOL : -STARRY_T_TOL;
+      theta = costheta > 0 ? 0 : pi<Scalar>();
+    }
+    if (f < STARRY_MIN_F)
+      f = STARRY_MIN_F;
+  }
+
+public:
+  RowVector<A> sT;
+  RowVector<A> sTbar;
+  RowVector<A> sT0;
+  A phi1, phi2, xi1, xi2;
+
+  explicit Occultation(int deg) : deg(deg), ncoeff((deg + 1) * (deg + 1)) {
+    compute_phase();
+  }
+
+  /**
+      Compute the solution vector in the case
+      of no occultation. The phase curve solution
+      vector for any value of `f` is just
+
+        sT = (1 - f) * sT0
+
+  */
+  inline void compute_phase() {
+    compute_L0(Lt);
+    sT0.setZero(ncoeff);
+    int mu, nu, n = 0;
+    for (int l = 0; l < deg + 1; ++l) {
+      for (int m = -l; m < l + 1; ++m) {
+        mu = l - m;
+        nu = l + m;
+        if (is_even(nu)) {
+          // Case 1
+          sT0(n) = Lt(mu / 2 + 2, nu / 2);
+        } else if ((l == 1) && (m == 0)) {
+          // Case 2
+          sT0(n) = 2 * pi<Scalar>() / 3;
+        }
+        ++n;
+      }
+    }
+  }
+
+  /**
+      Compute the full solution vector s^T.
+
+
+  */
+  inline void compute(const A &bo, const A &ro, const A &f, const A &theta) {
+    compute_complement(bo, ro, f, theta);
+    sT = (1 - f) * sT0 - sTbar;
+  }
+
+  /**
+      Compute the complement of the solution vector.
+      This is the integral *under* the occultor.
+
+  */
+  inline void compute_complement(const A &bo_, const A &ro_, const A &f_,
+                                 const A &theta_) {
+    // Make local copies of the inputs
+    // Nudge them away from singular points
+    bo = bo_;
+    ro = ro_;
+    f = f_;
+    theta = theta_;
+    costheta = cos(theta);
+    sintheta = sin(theta);
+    nudge_inputs();
+
+    // Compute the angles of intersection
+    if (ro == 0) {
+
+      // No occultation
+      sTbar.setZero(ncoeff);
+      return;
+
+    } else {
+
+      get_angles(bo, ro, f, theta, phi1, phi2, xi1, xi2);
+    }
+
+    // Special cases
+    if ((phi1 == 0.0) && (phi2 == 0.0) && (xi1 == 0.0) &&
+        (xi2 == 2 * pi<Scalar>())) {
+
+      // Complete occultation
+      sTbar = (1 - f) * sT0;
+      return;
+
+    } else if ((phi1 == 0.0) && (phi2 == 0.0) && (xi1 == 0.0) && (xi2 == 0.0)) {
+
+      // No occultation
+      sTbar.setZero(ncoeff);
+      return;
+    }
+
+    // Compute the M integral analytically (even `mu`)
+    S.setZero(deg + 3, deg + 3);
+    C.setZero(deg + 3, deg + 3);
+    A fac, fac0, facc, facs;
+    fac0 = 1.0;
+    for (int i = 0; i < deg + 3; ++i) {
+      facs = ro * fac0;
+      facc = fac0;
+      for (int j = 0; j < deg + 3; ++j) {
+        if (j < i + 1) {
+          S(i, deg + 2 + j - i) = facs;
+          C(deg + 2 + j - i, i) = facc;
+          fac = bo * (i - j) / (ro * (j + 1.0));
+          facs *= fac * sintheta;
+          facc *= fac * costheta;
+        }
+      }
+      fac0 *= ro;
+    }
+    compute_L(phi1 - theta, phi2 - theta, Lp);
+    M = S.topRightCorner(deg + 2, deg + 2) * Lp.reverse().topRows(deg + 2) *
+        C.leftCols(deg + 2);
+
+    // Compute the K integral numerically (odd `mu`)
+    K = pTodd(deg, bo, ro, f, theta, costheta, sintheta, phi1, phi2, QUAD);
+
+    // Compute L (t integral)
+    compute_L(xi1, xi2, Lt);
+
+    // Go through the cases
+    A pT, tT;
+    sTbar.resize(ncoeff);
+    int mu, nu, n = 0;
+    for (int l = 0; l < deg + 1; ++l) {
+      for (int m = -l; m < l + 1; ++m) {
+        mu = l - m;
+        nu = l + m;
+
+        // Compute the pT and tT integrals
+        if (is_even(nu)) {
+          // Case 1
+          pT = pow(1 - f, -(nu / 2)) * M((mu + 2) / 2, nu / 2);
+          tT = (1 - f) * Lt(mu / 2 + 2, nu / 2);
+        } else if ((l == 1) && (m == 0)) {
+          // Case 2
+          pT = K(n);
+          tT = (1 - f) * (xi2 - xi1) / 3;
+        } else if (is_even(l) && (mu == 1)) {
+          // Case 3
+          pT = K(n);
+          tT = 0.0;
+        } else if (!is_even(l) && (mu == 1)) {
+          // Case 4
+          pT = K(n);
+          tT = 0.0;
+        } else {
+          // Case 5
+          pT = K(n);
+          tT = 0.0;
+        }
+
+        // The surface integral is just their sum
+        sTbar(n) = pT + tT;
+        ++n;
+      }
+    }
+  }
+};
+
+} // namespace occultation
+} // namespace oblate
+} // namespace starry
+
+#endif
\ No newline at end of file
diff --git a/starry/_core/ops/lib/include/oblate/setup.py b/starry/_core/ops/lib/include/oblate/setup.py
new file mode 100644
index 000000000..81befbab5
--- /dev/null
+++ b/starry/_core/ops/lib/include/oblate/setup.py
@@ -0,0 +1,141 @@
+"""Install script for `starry`."""
+from setuptools import setup, Extension, find_packages
+from setuptools.command.build_ext import build_ext
+import warnings
+import setuptools
+import subprocess
+import sys
+import os
+import glob
+
+# Debug mode?
+debug = bool(int(os.getenv("STARRY_DEBUG", 0)))
+if debug:
+    optimize = 0
+    macros = {}
+    macros["STARRY_O"] = 0
+    macros["STARRY_DEBUG"] = 1
+else:
+    optimize = 2
+    macros = {}
+
+
+class get_pybind_include(object):
+    """
+    Helper class to determine the pybind11 include path.
+
+    The purpose of this class is to postpone importing pybind11
+    until it is actually installed, so that the ``get_include()``
+    method can be invoked.
+    """
+
+    def __init__(self, user=False):
+        self.user = user
+
+    def __str__(self):
+        import pybind11
+
+        return pybind11.get_include(self.user)
+
+
+# As of Python 3.6, CCompiler has a `has_flag` method.
+# cf http://bugs.python.org/issue26689
+def has_flag(compiler, flagname):
+    """Return a boolean indicating whether a flag name is supported on
+    the specified compiler.
+    """
+    import tempfile
+
+    with tempfile.NamedTemporaryFile("w", suffix=".cpp") as f:
+        f.write("int main (int argc, char **argv) { return 0; }")
+        try:
+            compiler.compile([f.name], extra_postargs=["-w", flagname])
+        except setuptools.distutils.errors.CompileError:
+            return False
+    return True
+
+
+def cpp_flag(compiler):
+    """Return the -std=c++[11/14/17] compiler flag.
+
+    The newer version is prefered over c++11 (when it is available).
+    """
+    flags = ["-std=c++17", "-std=c++14", "-std=c++11"]
+
+    for flag in flags:
+        if has_flag(compiler, flag):
+            return flag
+
+    raise RuntimeError(
+        "Unsupported compiler -- at least C++11 support is needed!"
+    )
+
+
+class BuildExt(build_ext):
+    """A custom build extension for adding compiler-specific options."""
+
+    c_opts = {"msvc": ["/EHsc"], "unix": []}
+    l_opts = {"msvc": [], "unix": []}
+
+    if sys.platform == "darwin":
+        darwin_opts = ["-stdlib=libc++", "-mmacosx-version-min=10.14"]
+        c_opts["unix"] += darwin_opts
+        l_opts["unix"] += darwin_opts
+
+    def build_extensions(self):
+        ct = self.compiler.compiler_type
+        opts = self.c_opts.get(ct, [])
+        link_opts = self.l_opts.get(ct, [])
+        if ct == "unix":
+            opts.append(
+                '-DVERSION_INFO="%s"' % self.distribution.get_version()
+            )
+            opts.append(cpp_flag(self.compiler))
+            if has_flag(self.compiler, "-fcolor-diagnostics"):
+                opts.append("-fcolor-diagnostics")
+        elif ct == "msvc":
+            opts.append(
+                '/DVERSION_INFO=\\"%s\\"' % self.distribution.get_version()
+            )
+            opts.append("/Zm10")  # debug for C1060
+        extra_args = ["-O%d" % optimize]
+        if debug:
+            extra_args += ["-g", "-Wall", "-fno-lto"]
+        else:
+            if has_flag(self.compiler, "-fvisibility=hidden"):
+                opts.append("-fvisibility=hidden")
+            extra_args += ["-g0"]
+        for ext in self.extensions:
+            ext.extra_compile_args = opts + extra_args
+            ext.extra_link_args = link_opts + extra_args
+
+        build_ext.build_extensions(self)
+
+
+ext_modules = [
+    Extension(
+        "oblate._c_ops",
+        ["oblate.cpp"],
+        include_dirs=[
+            # Path to pybind11 headers
+            get_pybind_include(),
+            get_pybind_include(user=True),
+            ".",
+            "../../vendor/eigen_3.3.5",
+            "../../vendor/boost_1_75_0",
+        ],
+        language="c++",
+        define_macros=[(key, value) for key, value in macros.items()],
+    )
+]
+cmdclass = {"build_ext": BuildExt}
+
+
+setup(
+    name="oblate",
+    packages=find_packages(),
+    ext_modules=ext_modules,
+    cmdclass=cmdclass,
+    include_package_data=True,
+    zip_safe=False,
+)
diff --git a/starry/_core/ops/lib/include/ops.h b/starry/_core/ops/lib/include/ops.h
index 0df34a4e8..191695e31 100644
--- a/starry/_core/ops/lib/include/ops.h
+++ b/starry/_core/ops/lib/include/ops.h
@@ -7,6 +7,7 @@
 #include "basis.h"
 #include "filter.h"
 #include "misc.h"
+#include "oblate/occultation.h"
 #include "reflected/occultation.h"
 #include "reflected/phasecurve.h"
 #include "solver.h"
@@ -29,12 +30,19 @@ template <class Scalar> class Ops {
   const int deg;
   const int N;
 
+  // Standard starry
   basis::Basis<Scalar> B;
   wigner::Wigner<Scalar> W;
-  solver::Greens<Scalar> G; /**< The occultation integral solver class */
+  solver::Greens<Scalar> G;
+  filter::Filter<Scalar> F;
+
+  // Reflected light starry
   reflected::phasecurve::PhaseCurve<ADScalar<Scalar, 2>> RP;
   reflected::occultation::Occultation<ADScalar<Scalar, 5>> RO;
-  filter::Filter<Scalar> F;
+
+  // Oblate starry
+  oblate::occultation::Occultation<Scalar, 0> OBL;
+  oblate::occultation::Occultation<Scalar, 4> OBLAD;
 
   // Spot gradients
   RowVector<Scalar> bamp;
@@ -47,7 +55,7 @@ template <class Scalar> class Ops {
       : ydeg(ydeg), Ny((ydeg + 1) * (ydeg + 1)), udeg(udeg), Nu(udeg + 1),
         fdeg(fdeg), Nf((fdeg + 1) * (fdeg + 1)), deg(ydeg + udeg + fdeg),
         N((deg + 1) * (deg + 1)), B(ydeg, udeg, fdeg), W(ydeg, udeg, fdeg),
-        G(deg), RP(deg, B), RO(deg, B), F(B) {
+        G(deg), F(B), RP(deg, B), RO(deg, B), OBL(deg), OBLAD(deg) {
     // Bounds checks
     if ((ydeg < 0) || (ydeg > STARRY_MAX_LMAX))
       throw std::out_of_range("Spherical harmonic degree out of range.");
diff --git a/starry/_core/ops/lib/include/quad.h b/starry/_core/ops/lib/include/quad.h
index 9bf1c8095..fa90c0183 100644
--- a/starry/_core/ops/lib/include/quad.h
+++ b/starry/_core/ops/lib/include/quad.h
@@ -16,45 +16,60 @@ namespace quad {
 using namespace utils;
 
 /*! Implementation of Gauss-Legendre quadrature
-*  http://en.wikipedia.org/wiki/Gaussian_quadrature
-*  http://rosettacode.org/wiki/Numerical_integration/Gauss-Legendre_Quadrature
-*
-*/
+ *  http://en.wikipedia.org/wiki/Gaussian_quadrature
+ *  http://rosettacode.org/wiki/Numerical_integration/Gauss-Legendre_Quadrature
+ *
+ */
 template <typename T> class Quad {
 public:
   enum { eDEGREE = STARRY_QUAD_POINTS };
 
   /*! Compute the integral of a functor
-  *
-  *   @param a    lower limit of integration
-  *   @param b    upper limit of integration
-  *   @param f    the function to integrate
-  *   @param err  callback in case of problems
-  */
+   *
+   *   @param a    lower limit of integration
+   *   @param b    upper limit of integration
+   *   @param f    the function to integrate
+   */
   template <typename Function> inline T integrate(T a, T b, Function f) {
     T p = (b - a) / 2;
     T q = (b + a) / 2;
     const LegendrePolynomial &legpoly = s_LegendrePolynomial;
 
     T sum = 0;
-    for (int i = 1; i <= eDEGREE; ++i) {
+    for (int i = 0; i < eDEGREE; ++i) {
       sum += legpoly.weight(i) * f(p * legpoly.root(i) + q);
     }
 
     return p * sum;
   }
 
+  /*! Compute the vectorized integral of a functor
+   *
+   *   @param a    lower limit of integration
+   *   @param b    upper limit of integration
+   *   @param f    the function to integrate
+   *   @param npts the dimension of the function output
+   */
+  template <typename Function>
+  inline Vector<T> integrate(T a, T b, Function f, int npts) {
+    T p = 0.5 * (b - a);
+    Vector<T> q((int)eDEGREE);
+    q.setConstant(0.5 * (b + a));
+    const LegendrePolynomial &legpoly = s_LegendrePolynomial;
+    return p * legpoly.wT * f(p * legpoly.r + q);
+  }
+
   /*! Print out roots and weights for information
-  */
+   */
   void print_roots_and_weights(std::ostream &out) const {
     const LegendrePolynomial &legpoly = s_LegendrePolynomial;
     out << "Roots:  ";
-    for (int i = 0; i <= eDEGREE; ++i) {
+    for (int i = 0; i < eDEGREE; ++i) {
       out << ' ' << legpoly.root(i);
     }
     out << '\n';
     out << "Weights:";
-    for (int i = 0; i <= eDEGREE; ++i) {
+    for (int i = 0; i < eDEGREE; ++i) {
       out << ' ' << legpoly.weight(i);
     }
     out << '\n';
@@ -62,37 +77,38 @@ template <typename T> class Quad {
 
 private:
   /*! Implementation of the Legendre polynomials that form
-  *   the basis of this quadrature
-  */
+   *   the basis of this quadrature
+   */
   class LegendrePolynomial {
   public:
+    Eigen::Matrix<T, eDEGREE, 1> r;
+    Eigen::Matrix<T, 1, eDEGREE> wT;
+
     LegendrePolynomial() {
       // Solve roots and weights
-      for (int i = 0; i <= eDEGREE; ++i) {
+
+      for (int i = 0; i < eDEGREE; ++i) {
         T dr = 1;
 
         // Find zero
-        Evaluation eval(cos(M_PI * (i - 0.25) / (eDEGREE + 0.5)));
+        Evaluation eval(cos(M_PI * (i + 0.75) / (eDEGREE + 0.5)));
         do {
           dr = eval.v() / eval.d();
           eval.evaluate(eval.x() - dr);
         } while (abs(dr) > 2e-16);
 
-        this->_r[i] = eval.x();
-        this->_w[i] = 2 / ((1 - eval.x() * eval.x()) * eval.d() * eval.d());
+        this->r(i) = eval.x();
+        this->wT(i) = 2 / ((1 - eval.x() * eval.x()) * eval.d() * eval.d());
       }
     }
 
-    T root(int i) const { return this->_r[i]; }
-    T weight(int i) const { return this->_w[i]; }
+    T root(int i) const { return this->r(i); }
+    T weight(int i) const { return this->wT(i); }
 
   private:
-    T _r[eDEGREE + 1];
-    T _w[eDEGREE + 1];
-
     /*! Evaluate the value *and* derivative of the
-    *   Legendre polynomial
-    */
+     *   Legendre polynomial
+     */
     class Evaluation {
     public:
       explicit Evaluation(T x) : _x(x), _v(1), _d(0) { this->evaluate(x); }
@@ -125,7 +141,7 @@ template <typename T> class Quad {
   };
 
   /*! Pre-compute the weights and abscissae of the Legendre polynomials
-  */
+   */
   static LegendrePolynomial s_LegendrePolynomial;
 };
 
diff --git a/starry/_core/ops/lib/include/reflected/constants.h b/starry/_core/ops/lib/include/reflected/constants.h
deleted file mode 100644
index e9e73f328..000000000
--- a/starry/_core/ops/lib/include/reflected/constants.h
+++ /dev/null
@@ -1,109 +0,0 @@
-/**
-\file utils.h
-\brief Miscellaneous utilities and definitions used throughout the code.
-
-*/
-
-#ifndef _STARRY_REFLECTED_CONSTANTS_H_
-#define _STARRY_REFLECTED_CONSTANTS_H_
-
-// --------------------------
-// ------ User-tunable ------
-// --------------------------
-
-//! Use the incomplete elliptic integrals to compute P?
-#ifndef STARRY_USE_INCOMPLETE_INTEGRALS
-#define STARRY_USE_INCOMPLETE_INTEGRALS 0
-#endif
-
-//! Maximum number of iterations when computing `el2`
-#ifndef STARRY_EL2_MAX_ITER
-#define STARRY_EL2_MAX_ITER 100
-#endif
-
-// --------------------------
-// ---------- Fixed ---------
-// --------------------------
-
-// Integration codes
-#define FLUX_ZERO 0
-#define FLUX_DAY_OCC 1
-#define FLUX_DAY_VIS 2
-#define FLUX_NIGHT_OCC 3
-#define FLUX_NIGHT_VIS 4
-#define FLUX_SIMPLE_OCC 5
-#define FLUX_SIMPLE_REFL 6
-#define FLUX_SIMPLE_OCC_REFL 7
-#define FLUX_TRIP_DAY_OCC 8
-#define FLUX_TRIP_NIGHT_OCC 9
-#define FLUX_QUAD_DAY_VIS 10
-#define FLUX_QUAD_NIGHT_VIS 11
-#define FLUX_NOON 12
-
-// Maximum number of iterations & tolerance when computing 2F1
-#define STARRY_2F1_MAXITER 200
-#define STARRY_2F1_MAXTOL 1e-15
-#define STARRY_2F1_MINTOL 1e-12
-
-// Square root of the desired precision in `el2`
-#define STARRY_EL2_CA 1e-8
-
-// Replace `inf` with this value in argument to `el2`
-#define STARRY_HUGE_TAN 1e15
-
-// If |sin(theta)| or |cos(theta)| is less than this, set  0
-#define STARRY_T_TOL 1e-12
-
-// Low, medium, and high tolerance for root polishing
-#define STARRY_ROOT_TOL_LOW 1e-2
-#define STARRY_ROOT_TOL_MED 1e-10
-#define STARRY_ROOT_TOL_HIGH 1e-15
-
-// Tolerance for duplicate roots
-#define STARRY_ROOT_TOL_DUP 1e-7
-
-// Maximum number of root polishing iterations
-#define STARRY_ROOT_MAX_ITER 50
-
-// If |b| is less than this value, set equal to 0
-#define STARRY_B_ZERO_TOL 1e-8
-
-// If |b| is within this value of 1, set equal to 1
-#define STARRY_B_ONE_TOL 1e-6
-
-// Tolerance for various functions that calculate phi, xi, and lam
-#define STARRY_ANGLE_TOL 1e-13
-
-// Nudge k^2 away from 1 when it gets this close
-#define STARRY_K2_ONE_TOL 1e-12
-
-// Nudge bo away from ro when it gets this close
-#define STARRY_BO_EQUALS_RO_TOL 1e-8
-
-// Nudge theta away from pi / 2 when ro = 1
-#define STARRY_THETA_UNIT_RADIUS_TOL 1e-5
-
-// Smallest value of sin(kappa / 2) in computing W
-#define STARRY_MIN_SIN_ALPHA 1e-8
-
-/*
-Determining the integration paths close to the singular
-points of the occultation is quite hard, and the solution can
-often oscillate between two regimes. These tolerances prevent us
-from entering those regimes, at the cost of precision loss near
-these singular points.
-*/
-#define STARRY_COMPLETE_OCC_TOL 1e-8
-#define STARRY_NO_OCC_TOL 1e-8
-#define STARRY_GRAZING_TOL 1e-7
-
-// Tolerances for the Pal (2012) solver, which is very unstable
-#define STARRY_CRJ_MAX_ITER 100
-#define STARRY_CRJ_LO_LIM 2e-26
-#define STARRY_CRJ_HI_LIM 3e24
-#define STARRY_CRJ_TOL 2e-2
-#define STARRY_PAL_BO_EQUALS_RO_TOL 1e-3
-#define STARRY_PAL_BO_EQUALS_RO_MINUS_ONE_TOL 1e-3
-#define STARRY_PAL_BO_EQUALS_ONE_MINUS_RO_TOL 1e-3
-
-#endif
\ No newline at end of file
diff --git a/starry/_core/ops/lib/include/reflected/ellip.h b/starry/_core/ops/lib/include/reflected/ellip.h
index feaef2770..299890b37 100644
--- a/starry/_core/ops/lib/include/reflected/ellip.h
+++ b/starry/_core/ops/lib/include/reflected/ellip.h
@@ -15,7 +15,6 @@ Elliptic integrals computed following
 
 #include "../ellip.h"
 #include "../utils.h"
-#include "constants.h"
 
 namespace starry {
 
diff --git a/starry/_core/ops/lib/include/reflected/geometry.h b/starry/_core/ops/lib/include/reflected/geometry.h
index ec3dad690..e3f14a534 100644
--- a/starry/_core/ops/lib/include/reflected/geometry.h
+++ b/starry/_core/ops/lib/include/reflected/geometry.h
@@ -8,7 +8,6 @@
 #define _STARRY_GEOMETRY_H_
 
 #include "../utils.h"
-#include "constants.h"
 #include <Eigen/Eigenvalues>
 
 namespace starry {
@@ -213,12 +212,10 @@ inline Vector<T> get_roots(const T &b_, const T &theta_, const T &costheta_,
       dxdb(n) = s *
                 sqrt((1 - x(n) * x(n)) * (ro2 - (x(n) - xo) * (x(n) - xo))) /
                 (x(n) - xo);
-      dxdtheta(n) = bo * (costheta -
-                          s * sqrt(ro2 - (x(n) - xo) * (x(n) - xo)) /
-                              (x(n) - xo) * sintheta);
-      dxdbo(n) =
-          sintheta +
-          s * sqrt(ro2 - (x(n) - xo) * (x(n) - xo)) / (x(n) - xo) * costheta;
+      dxdtheta(n) = bo * (costheta - s * sqrt(ro2 - (x(n) - xo) * (x(n) - xo)) /
+                                         (x(n) - xo) * sintheta);
+      dxdbo(n) = sintheta + s * sqrt(ro2 - (x(n) - xo) * (x(n) - xo)) /
+                                (x(n) - xo) * costheta;
       dxdro(n) = ro / (x(n) - xo);
     }
 
diff --git a/starry/_core/ops/lib/include/reflected/occultation.h b/starry/_core/ops/lib/include/reflected/occultation.h
index 6c879b1cd..1ff56d64c 100644
--- a/starry/_core/ops/lib/include/reflected/occultation.h
+++ b/starry/_core/ops/lib/include/reflected/occultation.h
@@ -13,7 +13,6 @@ light).
 #include "../quad.h"
 #include "../solver.h"
 #include "../utils.h"
-#include "constants.h"
 #include "geometry.h"
 #include "phasecurve.h"
 #include "primitive.h"
@@ -99,7 +98,7 @@ template <class T> class Occultation {
 
   /**
    *
-  */
+   */
   inline RowVector<T> sTr(const T &b, const T &theta, const T &sigr) {
 
     // Compute the reflection solution in the terminator frame
diff --git a/starry/_core/ops/lib/include/reflected/primitive.h b/starry/_core/ops/lib/include/reflected/primitive.h
index 4c3ca60d0..a716c3865 100644
--- a/starry/_core/ops/lib/include/reflected/primitive.h
+++ b/starry/_core/ops/lib/include/reflected/primitive.h
@@ -8,7 +8,6 @@
 #define _STARRY_PRIMITIVE_H_
 
 #include "../utils.h"
-#include "constants.h"
 #include "ellip.h"
 #include "special.h"
 
@@ -246,11 +245,11 @@ inline Vector<T> W_indef(const int nmax, const T &s2_, const T &q2,
     T x = q2 * q3 * s2nmax;
 
     // Upper boundary condition
-    result(nmax) =
-        (s2 * s2nmax * (3.0 / (nmax + 1.0) *
-                            hyp2f1(-0.5, nmax + 1.0, nmax + 2.0, T(1.0 - q2)) +
-                        2 * q3) /
-         (2.0 * nmax + 5.0));
+    result(nmax) = (s2 * s2nmax *
+                    (3.0 / (nmax + 1.0) *
+                         hyp2f1(-0.5, nmax + 1.0, nmax + 2.0, T(1.0 - q2)) +
+                     2 * q3) /
+                    (2.0 * nmax + 5.0));
 
     // Recurse down
     T f, A, B;
diff --git a/starry/_core/ops/lib/include/reflected/scatter.h b/starry/_core/ops/lib/include/reflected/scatter.h
index 61595da94..517c8720f 100644
--- a/starry/_core/ops/lib/include/reflected/scatter.h
+++ b/starry/_core/ops/lib/include/reflected/scatter.h
@@ -9,7 +9,6 @@
 
 #include "../basis.h"
 #include "../utils.h"
-#include "constants.h"
 #include "oren_nayar.h"
 
 namespace starry {
@@ -130,9 +129,8 @@ inline Vector<T> OrenNayarPolynomial(const T &b, const T &theta, const T &sigr,
     // coefficients in the polynomial basis
     for (int n = 0; n < STARRY_OREN_NAYAR_N; ++n) {
       for (int m = 0; m < STARRY_OREN_NAYAR_NB * STARRY_OREN_NAYAR_NB; ++m) {
-        f(n) += STARRY_OREN_NAYAR_COEFFS[m +
-                                         STARRY_OREN_NAYAR_NB *
-                                             STARRY_OREN_NAYAR_NB * n] *
+        f(n) += STARRY_OREN_NAYAR_COEFFS[m + STARRY_OREN_NAYAR_NB *
+                                                 STARRY_OREN_NAYAR_NB * n] *
                 bbc(m);
       }
     }
diff --git a/starry/_core/ops/lib/include/reflected/special.h b/starry/_core/ops/lib/include/reflected/special.h
index a95291f98..b0f73185e 100644
--- a/starry/_core/ops/lib/include/reflected/special.h
+++ b/starry/_core/ops/lib/include/reflected/special.h
@@ -9,7 +9,6 @@
 
 #include "../quad.h"
 #include "../utils.h"
-#include "constants.h"
 
 namespace starry {
 namespace reflected {
@@ -133,13 +132,13 @@ inline T P2_numerical(const T &bo, const T &ro, const Vector<T> &kappa,
     res.derivatives() -=
         f(kappa(i).value() - pi<Scalar>()) * kappa(i).derivatives();
     res.derivatives() +=
-        bo.derivatives() *
-        QUAD.integrate(kappa(i).value() - pi<Scalar>(),
-                       kappa(i + 1).value() - pi<Scalar>(), dfdbo);
+        bo.derivatives() * QUAD.integrate(kappa(i).value() - pi<Scalar>(),
+                                          kappa(i + 1).value() - pi<Scalar>(),
+                                          dfdbo);
     res.derivatives() +=
-        ro.derivatives() *
-        QUAD.integrate(kappa(i).value() - pi<Scalar>(),
-                       kappa(i + 1).value() - pi<Scalar>(), dfdro);
+        ro.derivatives() * QUAD.integrate(kappa(i).value() - pi<Scalar>(),
+                                          kappa(i + 1).value() - pi<Scalar>(),
+                                          dfdro);
   }
 
   return res;
@@ -299,9 +298,9 @@ inline T J_numerical(const int N, const T &k2, const Vector<T> &kappa,
         0.5 * f(0.5 * kappa(i + 1).value()) * kappa(i + 1).derivatives();
     res.derivatives() -=
         0.5 * f(0.5 * kappa(i).value()) * kappa(i).derivatives();
-    res.derivatives() += k2.derivatives() *
-                         QUAD.integrate(0.5 * kappa(i).value(),
-                                        0.5 * kappa(i + 1).value(), dfdk2);
+    res.derivatives() +=
+        k2.derivatives() * QUAD.integrate(0.5 * kappa(i).value(),
+                                          0.5 * kappa(i + 1).value(), dfdk2);
   }
 
   return res;
diff --git a/starry/_core/ops/lib/include/utils.h b/starry/_core/ops/lib/include/utils.h
index fc4d77d7e..52ae3da35 100644
--- a/starry/_core/ops/lib/include/utils.h
+++ b/starry/_core/ops/lib/include/utils.h
@@ -8,6 +8,7 @@
 #define _STARRY_UTILS_H_
 
 // Includes
+#include "macros.h"
 #include <Eigen/Core>
 #include <Eigen/Dense>
 #include <Eigen/SparseLU>
@@ -22,14 +23,14 @@
 #include <unsupported/Eigen/AutoDiff>
 #include <vector>
 
-//! Number of digits (16 = double)
-#ifndef STARRY_NDIGITS
-#define STARRY_NDIGITS 16
-#else
+//! Number of digits of precision (16 = double)
+#ifdef STARRY_NDIGITS
 #if STARRY_NDIGITS > 16
 #define STARRY_ENABLE_BOOST
 #endif
 #endif
+
+//! Use the boost library?
 #ifdef STARRY_ENABLE_BOOST
 #include <boost/math/special_functions/gamma.hpp>
 #include <boost/multiprecision/cpp_dec_float.hpp>
@@ -49,41 +50,6 @@ typedef void Multi;
 #define unlikely(x) __builtin_expect(!!(x), 0)
 #endif
 
-//! Number of Gaussian-Legendre quadrature points for numerical integration
-#ifndef STARRY_QUAD_POINTS
-#define STARRY_QUAD_POINTS 100
-#endif
-
-//! Max iterations in elliptic integrals
-#ifndef STARRY_ELLIP_MAX_ITER
-#define STARRY_ELLIP_MAX_ITER 200
-#endif
-
-//! Max iterations in computing the M & N integrals
-#ifndef STARRY_MN_MAX_ITER
-#define STARRY_MN_MAX_ITER 100
-#endif
-
-//! Max iterations in computing the I & J integrals
-#ifndef STARRY_IJ_MAX_ITER
-#define STARRY_IJ_MAX_ITER 200
-#endif
-
-//! Refine the downward recursion in the J integral at this index
-#ifndef STARRY_REFINE_J_AT
-#define STARRY_REFINE_J_AT 25
-#endif
-
-//! Cutoff value for `b` below which we reparametrize LD evaluation
-#ifndef STARRY_BCUT
-#define STARRY_BCUT 1.0e-3
-#endif
-
-//! Things currently go numerically unstable in our bases for high `l`
-#ifndef STARRY_MAX_LMAX
-#define STARRY_MAX_LMAX 50
-#endif
-
 //! The value of `pi` in double precision
 #ifndef M_PI
 #define M_PI 3.14159265358979323846264338328
@@ -94,6 +60,21 @@ typedef void Multi;
 #define M_SQRTPI 1.77245385090551602729816748334
 #endif
 
+// Integration codes (reflected light / oblate)
+#define FLUX_ZERO 0
+#define FLUX_DAY_OCC 1
+#define FLUX_DAY_VIS 2
+#define FLUX_NIGHT_OCC 3
+#define FLUX_NIGHT_VIS 4
+#define FLUX_SIMPLE_OCC 5
+#define FLUX_SIMPLE_REFL 6
+#define FLUX_SIMPLE_OCC_REFL 7
+#define FLUX_TRIP_DAY_OCC 8
+#define FLUX_TRIP_NIGHT_OCC 9
+#define FLUX_QUAD_DAY_VIS 10
+#define FLUX_QUAD_NIGHT_VIS 11
+#define FLUX_NOON 12
+
 // Bounds checks
 #ifdef STARRY_DEBUG
 //! Check matrix shape (debug mode only)
@@ -127,7 +108,6 @@ typedef void Multi;
 #define CHECK_BOUNDS(INDEX, IMIN, IMAX)                                        \
   do {                                                                         \
   } while (0)
-
 #endif
 
 namespace starry {
@@ -156,6 +136,7 @@ template <typename T> using Vector = Eigen::Matrix<T, Eigen::Dynamic, 1>;
 template <typename T> using UnitVector = Eigen::Matrix<T, 3, 1>;
 template <typename T> using RowVector = Eigen::Matrix<T, 1, Eigen::Dynamic>;
 template <typename T> using OneByOne = Eigen::Matrix<T, 1, 1>;
+template <typename T> using Pair = Eigen::Matrix<T, 2, 1>;
 template <typename T, int StorageOrder = ColMajor>
 using Matrix = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, StorageOrder>;
 template <typename T, int N>
diff --git a/starry/_plotting.py b/starry/_plotting.py
index d242065ae..6309d0638 100644
--- a/starry/_plotting.py
+++ b/starry/_plotting.py
@@ -69,7 +69,9 @@ def get_moll_longitude_lines(dlon=np.pi / 6, npts=1000, niter=100):
     return res
 
 
-def get_ortho_latitude_lines(inc=np.pi / 2, obl=0, dlat=np.pi / 6, npts=1000):
+def get_ortho_latitude_lines(
+    inc=np.pi / 2, obl=0, fproj=0, dlat=np.pi / 6, npts=1000
+):
     """
 
     """
@@ -107,15 +109,15 @@ def get_ortho_latitude_lines(inc=np.pi / 2, obl=0, dlat=np.pi / 6, npts=1000):
 
         # Rotate them
         for y in (y1, y2):
-            xr = -x * co - y * so
-            yr = -x * so + y * co
+            xr = -x * co - (y * (1 - fproj)) * so
+            yr = -x * so + (y * (1 - fproj)) * co
             res.append((xr, yr))
 
     return res
 
 
 def get_ortho_longitude_lines(
-    inc=np.pi / 2, obl=0, theta=0, dlon=np.pi / 6, npts=1000
+    inc=np.pi / 2, obl=0, fproj=0, theta=0, dlon=np.pi / 6, npts=1000
 ):
     """
 
@@ -185,8 +187,8 @@ def get_ortho_longitude_lines(
                         y[imax:] = np.nan
 
             # Rotate by the obliquity
-            xr = -x * co - y * so
-            yr = -x * so + y * co
+            xr = -x * co - (y * (1 - fproj)) * so
+            yr = -x * so + (y * (1 - fproj)) * co
             res.append((xr, yr))
 
     return res
diff --git a/starry/compat.py b/starry/compat.py
index 5eeb6cf40..339a5a19b 100644
--- a/starry/compat.py
+++ b/starry/compat.py
@@ -7,6 +7,7 @@
 from aesara_theano_fallback import change_flags, ifelse, USE_AESARA
 from aesara_theano_fallback.tensor import slinalg
 from aesara_theano_fallback.graph import basic, op, params_type, fg
+from inspect import getmro
 
 __all__ = [
     "theano",
@@ -54,8 +55,8 @@ def evaluator(**kwargs):
 
     """
     # Store the kwargs
-    kwargs_point = kwargs.pop("point", None)
-    kwargs_model = kwargs.pop("model", None)
+    kwargs_point = kwargs.get("point", None)
+    kwargs_model = kwargs.get("model", None)
 
     if kwargs_point is not None:
 
@@ -68,7 +69,12 @@ def evaluator(**kwargs):
         model = kwargs_model
         if model is None:
             model = pm.Model.get_context()
-        get_val = lambda x: pmx.eval_in_model(x, model=model, point=point)
+
+        def get_val(x):
+            if Node not in getmro(type(x)):
+                return x
+            else:
+                return pmx.eval_in_model(x, model=model, point=point)
 
     else:
 
@@ -76,6 +82,12 @@ def evaluator(**kwargs):
 
         def get_val(x):
 
+            if Node not in getmro(type(x)):
+
+                # If it's not a tensor type, we don't
+                # have to do anything!
+                return x
+
             try:
 
                 # Try to directly evaluate it
diff --git a/starry/kepler.py b/starry/kepler.py
index 67427b668..a5b73db8a 100644
--- a/starry/kepler.py
+++ b/starry/kepler.py
@@ -9,6 +9,7 @@
 from inspect import getmro
 import matplotlib.pyplot as plt
 from matplotlib.animation import FuncAnimation
+from matplotlib.patches import Ellipse
 from IPython.display import HTML
 import os
 import logging
@@ -449,6 +450,7 @@ def __init__(
         ), "Reflected light map not allowed for the primary body."
         self._primary = primary
         self._rv = primary._map.__props__["rv"]
+        self._oblate = primary._map.__props__["oblate"]
         self._lazy = primary._lazy
         self._math = primary._math
         if self._lazy:
@@ -473,6 +475,9 @@ def __init__(
             assert (
                 sec._lazy == self._lazy
             ), "All bodies must have the same evaluation mode (lazy/greedy)."
+            assert not sec._map.__props__[
+                "oblate"
+            ], "Oblate secondary bodies are not currently supported."
 
         reflected = [sec._map.__props__["reflected"] for sec in secondaries]
         if np.all(reflected):
@@ -505,6 +510,7 @@ def __init__(
             self._secondaries,
             reflected=self._reflected,
             rv=self._rv,
+            oblate=self._oblate,
             light_delay=self._light_delay,
             texp=self._texp,
             oversample=self._oversample,
@@ -642,8 +648,17 @@ def show(
             self._primary._prot,
             self._primary._t0,
             self._primary._theta0,
-            self._primary._map._inc,
-            self._primary._map._obl,
+            getattr(
+                self._primary._map,
+                "_inc",
+                self._math.to_array_or_tensor(0.5 * np.pi),
+            ),
+            getattr(
+                self._primary._map, "_obl", self._math.to_array_or_tensor(0.0)
+            ),
+            getattr(
+                self._primary._map, "fproj", self._math.to_array_or_tensor(0.0)
+            ),
             self._primary._map._y,
             self._primary._map._u,
             self._primary._map._f,
@@ -676,10 +691,22 @@ def show(
                 [sec._inc for sec in self._secondaries]
             ),
             self._math.to_array_or_tensor(
-                [sec._map._inc for sec in self._secondaries]
+                [
+                    getattr(
+                        sec._map,
+                        "_inc",
+                        self._math.to_array_or_tensor(0.5 * np.pi),
+                    )
+                    for sec in self._secondaries
+                ]
             ),
             self._math.to_array_or_tensor(
-                [sec._map._obl for sec in self._secondaries]
+                [
+                    getattr(
+                        sec._map, "_obl", self._math.to_array_or_tensor(0.0)
+                    )
+                    for sec in self._secondaries
+                ]
             ),
             self._math.to_array_or_tensor(
                 [sec._map._y for sec in self._secondaries]
@@ -691,11 +718,16 @@ def show(
                 [sec._map._f for sec in self._secondaries]
             ),
             self._math.to_array_or_tensor(
-                [sec._map._sigr for sec in self._secondaries]
+                [
+                    getattr(
+                        sec._map, "_sigr", self._math.to_array_or_tensor(0.0)
+                    )
+                    for sec in self._secondaries
+                ]
             ),
         )
 
-        # Convert to units of the primary radiu
+        # Convert to units of the primary radius
         x, y, z = (
             x / self._primary._r,
             y / self._primary._r,
@@ -707,13 +739,12 @@ def show(
         r = r / self._primary._r
 
         # Evaluate if needed
-        if self._lazy:
-            img_pri = get_val(img_pri)
-            img_sec = get_val(img_sec)
-            x = get_val(x)
-            y = get_val(y)
-            z = get_val(z)
-            r = get_val(r)
+        img_pri = get_val(img_pri)
+        img_sec = get_val(img_sec)
+        x = get_val(x)
+        y = get_val(y)
+        z = get_val(z)
+        r = get_val(r)
 
         # We need this to be of shape (nplanet, nframe)
         x = x.T
@@ -750,9 +781,21 @@ def show(
             animated=animated,
             zorder=0.0,
         )
-        circ[0] = plt.Circle(
-            (0, 0), 1, color="k", fill=False, zorder=1e-3, lw=2
-        )
+        if hasattr(self._primary._map, "fproj"):
+            circ[0] = Ellipse(
+                (0, 0),
+                2,
+                2 * (1 - get_val(self._primary._map.fproj)),
+                angle=180 / np.pi * get_val(self._primary._map._obl),
+                color="k",
+                fill=False,
+                zorder=1e-3,
+                lw=2,
+            )
+        else:
+            circ[0] = plt.Circle(
+                (0, 0), 1, color="k", fill=False, zorder=1e-3, lw=2
+            )
         ax.add_artist(circ[0])
         for i, _ in enumerate(self._secondaries):
             extent = np.array([x[i, 0], x[i, 0], y[i, 0], y[i, 0]]) + (
@@ -876,9 +919,20 @@ def design_matrix(self, t):
             self._primary._prot,
             self._primary._t0,
             self._primary._theta0,
-            self._math.to_array_or_tensor(1.0),
-            self._primary._map._inc,
-            self._primary._map._obl,
+            self._math.to_array_or_tensor(
+                1.0
+            ),  # we treat `amp` seprately in `flux()`
+            getattr(
+                self._primary._map,
+                "_inc",
+                self._math.to_array_or_tensor(0.5 * np.pi),
+            ),
+            getattr(
+                self._primary._map, "_obl", self._math.to_array_or_tensor(0.0)
+            ),
+            getattr(
+                self._primary._map, "fproj", self._math.to_array_or_tensor(0.0)
+            ),
             self._primary._map._u,
             self._primary._map._f,
             self._math.to_array_or_tensor(
@@ -909,17 +963,28 @@ def design_matrix(self, t):
             self._math.to_array_or_tensor(
                 [sec._inc for sec in self._secondaries]
             ),
+            self._math.to_array_or_tensor(
+                np.ones(
+                    len(self._secondaries)
+                )  # we treat `amp` seprately in `flux()`
+            ),
             self._math.to_array_or_tensor(
                 [
-                    self._math.to_array_or_tensor(1.0)
+                    getattr(
+                        sec._map,
+                        "_inc",
+                        self._math.to_array_or_tensor(0.5 * np.pi),
+                    )
                     for sec in self._secondaries
                 ]
             ),
             self._math.to_array_or_tensor(
-                [sec._map._inc for sec in self._secondaries]
-            ),
-            self._math.to_array_or_tensor(
-                [sec._map._obl for sec in self._secondaries]
+                [
+                    getattr(
+                        sec._map, "_obl", self._math.to_array_or_tensor(0.0)
+                    )
+                    for sec in self._secondaries
+                ]
             ),
             self._math.to_array_or_tensor(
                 [sec._map._u for sec in self._secondaries]
@@ -928,11 +993,16 @@ def design_matrix(self, t):
                 [sec._map._f for sec in self._secondaries]
             ),
             self._math.to_array_or_tensor(
-                [sec._map._sigr for sec in self._secondaries]
+                [
+                    getattr(
+                        sec._map, "_sigr", self._math.to_array_or_tensor(0.0)
+                    )
+                    for sec in self._secondaries
+                ]
             ),
         )
 
-    def flux(self, t, total=True):
+    def flux(self, t, total=True, integrated=False):
         """Compute the system flux at times ``t``.
 
         Args:
@@ -953,16 +1023,57 @@ def flux(self, t, total=True):
                 self._primary.map.amp * body.map.amp * body._map._y
                 for body in self._secondaries
             ]
+        elif self._oblate and self.primary.map.nw is not None:
+            # The problem is not strictly linear, since the design matrix
+            # is now 3D (time, ylm index, wavelength). So we pre-weight
+            # the spherical harmonic vector by the filter at each wavelength
+            # to get a matrix (ylm index, wavelength), to which we can apply
+            # the design matrix (time, ylm index) to get the spectral
+            # light curve out.
+            if self._primary.map.ydeg == 0:
+                y = self._primary.map._f
+            else:
+                y = self._primary.map.ops.weight_ylms_by_grav_dark_filter(
+                    self._primary.map.y, self._primary.map._f
+                )
+            ay = [self._primary.map.amp * y] + [
+                body.map.amp * body._map._y for body in self._secondaries
+            ]
+
         else:
             ay = [body.map.amp * body._map._y for body in self._bodies]
 
         if total:
-            return self._math.dot(X, self._math.concatenate(ay))
+            ay = self._math.concatenate(ay)
+            if integrated and self.primary.map.nw is not None:
+                ay = self._math.sum(ay, axis=1)
+            return self._math.dot(X, ay)
         else:
-            return [
-                self._math.dot(X[:, idx], ay[i])
-                for i, idx in enumerate(self._inds)
-            ]
+            if self._oblate:
+                # Because of our weighting hack above, the indices of
+                # the primary's coefficients in X changed, so let's
+                # re-compute all indices
+                Npri = (
+                    self._primary._map.ydeg + self._primary._map.fdeg + 1
+                ) ** 2
+                Ny = [Npri] + [sec._map.Ny for sec in self._secondaries]
+                inds = []
+                cur = 0
+                for N in Ny:
+                    inds.append(cur + np.arange(N))
+                    cur += N
+            else:
+                inds = self._inds
+            if integrated and self.primary.map.nw is not None:
+                return [
+                    self._math.dot(X[:, idx], self._math.sum(ay[i], axis=1))
+                    for i, idx in enumerate(inds)
+                ]
+            else:
+                return [
+                    self._math.dot(X[:, idx], ay[i])
+                    for i, idx in enumerate(inds)
+                ]
 
     def rv(self, t, keplerian=True, total=True):
         """Compute the observed radial velocity of the system at times ``t``.
@@ -1051,8 +1162,17 @@ def rv(self, t, keplerian=True, total=True):
             self._primary._t0,
             self._primary._theta0,
             self._primary._map._amp,
-            self._primary._map._inc,
-            self._primary._map._obl,
+            getattr(
+                self._primary._map,
+                "_inc",
+                self._math.to_array_or_tensor(0.5 * np.pi),
+            ),
+            getattr(
+                self._primary._map, "_obl", self._math.to_array_or_tensor(0.0)
+            ),
+            getattr(
+                self._primary._map, "fproj", self._math.to_array_or_tensor(0.0)
+            ),
             self._primary._map._y,
             self._primary._map._u,
             self._primary._map._alpha,
@@ -1089,10 +1209,22 @@ def rv(self, t, keplerian=True, total=True):
                 [sec._map._amp for sec in self._secondaries]
             ),
             self._math.to_array_or_tensor(
-                [sec._map._inc for sec in self._secondaries]
+                [
+                    getattr(
+                        sec._map,
+                        "_inc",
+                        self._math.to_array_or_tensor(0.5 * np.pi),
+                    )
+                    for sec in self._secondaries
+                ]
             ),
             self._math.to_array_or_tensor(
-                [sec._map._obl for sec in self._secondaries]
+                [
+                    getattr(
+                        sec._map, "_obl", self._math.to_array_or_tensor(0.0)
+                    )
+                    for sec in self._secondaries
+                ]
             ),
             self._math.to_array_or_tensor(
                 [sec._map._y for sec in self._secondaries]
@@ -1104,7 +1236,12 @@ def rv(self, t, keplerian=True, total=True):
                 [sec._map._alpha for sec in self._secondaries]
             ),
             self._math.to_array_or_tensor(
-                [sec._map._sigr for sec in self._secondaries]
+                [
+                    getattr(
+                        sec._map, "_sigr", self._math.to_array_or_tensor(0.0)
+                    )
+                    for sec in self._secondaries
+                ]
             ),
             self._math.to_array_or_tensor(
                 [sec._map._veq for sec in self._secondaries]
diff --git a/starry/maps.py b/starry/maps.py
index 944d61530..b89103cf4 100644
--- a/starry/maps.py
+++ b/starry/maps.py
@@ -1,7 +1,7 @@
 # -*- coding: utf-8 -*-
 from . import config, legacy
 from ._constants import *
-from ._core import OpsYlm, OpsLD, OpsReflected, OpsRV, math
+from ._core import OpsYlm, OpsLD, OpsReflected, OpsRV, OpsOblate, math
 from ._core.utils import is_tensor
 from ._indices import integers, get_ylm_inds, get_ul_inds, get_ylmw_inds
 from ._plotting import (
@@ -35,6 +35,7 @@
     "LimbDarkenedBase",
     "RVBase",
     "ReflectedBase",
+    "OblateBase",
 ]
 
 
@@ -73,20 +74,16 @@ def __init__(self, ydeg, udeg, fdeg, nw, **kwargs):
         self._N = (ydeg + udeg + fdeg + 1) ** 2
         self._nw = nw
 
-        # Basic properties
-        self._inc = self._math.cast(0.5 * np.pi)
-        self._obl = self._math.cast(0.0)
-        self._sigr = self._math.cast(0.0)
-
         # Units
-        self.angle_unit = kwargs.pop("angle_unit", units.degree)
+        self.angle_unit = kwargs.get("angle_unit", units.degree)
+        self.wav_unit = kwargs.get("wav_unit", units.nm)
 
         # Initialize
         self.reset(**kwargs)
 
     @property
     def angle_unit(self):
-        """An ``astropy.units`` unit defining the angle metric for this map."""
+        """An ``astropy.units`` quantity defining the angle unit for this map."""
         return self._angle_unit
 
     @angle_unit.setter
@@ -95,6 +92,17 @@ def angle_unit(self, value):
         self._angle_unit = value
         self._angle_factor = value.in_units(units.radian)
 
+    @property
+    def wav_unit(self):
+        """An ``astropy.units`` quantity defining the wavelength unit for this map."""
+        return self._wav_unit
+
+    @angle_unit.setter
+    def wav_unit(self, value):
+        assert value.physical_type == "length"
+        self._wav_unit = value
+        self._wav_factor = value.in_units(units.m)
+
     @property
     def ydeg(self):
         """Spherical harmonic degree of the map. *Read-only*"""
@@ -177,6 +185,24 @@ def u(self):
         """
         return self._u
 
+    @property
+    def wav(self):
+        """
+        The wavelength(s) at which the flux is measured in units of ``wav_unit``.
+
+        """
+        return self._wav / self._wav_factor
+
+    @wav.setter
+    def wav(self, value):
+        value = self._math.cast(value) * self._wav_factor
+        if self._nw is None:
+            assert value.ndim == 0, "Wavelength must be a scalar."
+            self._wav = value
+        else:
+            assert value.ndim == 1, "Wavelength must be a vector."
+            self._wav = value * self._math.ones(self._nw)
+
     def __getitem__(self, idx):
         if isinstance(idx, integers) or isinstance(idx, slice):
             # User is accessing a limb darkening index
@@ -253,23 +279,57 @@ def __setitem__(self, idx, val):
         else:
             raise ValueError("Invalid map index.")
 
-    def _check_kwargs(self, method, kwargs):
-        if not config.quiet:
-            for key in kwargs.keys():
-                message = "Invalid keyword `{0}` in call to `{1}()`. Ignoring."
-                message = message.format(key, method)
-                logger.warning(message)
-
-    def _get_flux_kwargs(self, kwargs):
-        xo = kwargs.pop("xo", 0.0)
-        yo = kwargs.pop("yo", 0.0)
-        zo = kwargs.pop("zo", 1.0)
-        ro = kwargs.pop("ro", 0.0)
-        theta = kwargs.pop("theta", 0.0)
-        theta, xo, yo, zo = self._math.vectorize(theta, xo, yo, zo)
-        theta, xo, yo, zo, ro = self._math.cast(theta, xo, yo, zo, ro)
-        theta *= self._angle_factor
-        return theta, xo, yo, zo, ro
+    def _get_norm(self, image_arr, **kwargs):
+        norm = kwargs.get("norm", None)
+        if norm is None:
+            vmin = np.nanmin(image_arr)
+            vmax = np.nanmax(image_arr)
+            if np.abs(vmin - vmax) < 1e-9:
+                vmin -= 1e-9
+                vmax += 1e-9
+            norm = colors.Normalize(vmin=vmin, vmax=vmax)
+        return norm
+
+    def _get_projection(self, **kwargs):
+        return get_projection(kwargs.get("projection", "ortho"))
+
+    def _get_ortho_borders(self, **kwargs):
+        xp = np.linspace(-1, 1, 10000)
+        yp = np.sqrt(1 - xp ** 2)
+        xm = xp
+        ym = -yp
+        return xp, yp, xm, ym
+
+    def _updatefig(
+        self,
+        i,
+        img,
+        image,
+        img_overlay,
+        overlay,
+        projection,
+        grid,
+        lonlines,
+        latlines,
+        borders,
+        kwargs,
+    ):
+        img.set_array(image[i])
+        images = [img]
+        if overlay is not None:
+            img_overlay.set_array(overlay[i])
+            images += [img_overlay]
+        if (
+            projection == STARRY_ORTHOGRAPHIC_PROJECTION
+            and grid
+            and len(image) > 1
+            and self.nw is None
+        ):
+            lons = self._get_ortho_longitude_lines(i=i, **kwargs)
+            for n, l in enumerate(lons):
+                lonlines[n].set_xdata(l[0])
+                lonlines[n].set_ydata(l[1])
+        return tuple(images + lonlines + latlines + borders)
 
     def reset(self, **kwargs):
         """Reset all map coefficients and attributes.
@@ -281,10 +341,16 @@ def reset(self, **kwargs):
         if self.nw is None:
             y = np.zeros(self.Ny)
             y[0] = 1.0
+            wav = kwargs.get("wav", 650.0)
         else:
             y = np.zeros((self.Ny, self.nw))
             y[0, :] = 1.0
+            if self.nw == 1:
+                wav = kwargs.get("wav", np.array([650.0]))
+            else:
+                wav = kwargs.get("wav", np.linspace(400.0, 900.0, self.nw))
         self._y = self._math.cast(y)
+        self.wav = wav
 
         u = np.zeros(self.Nu)
         u[0] = -1.0
@@ -294,7 +360,7 @@ def reset(self, **kwargs):
         f[0] = np.pi
         self._f = self._math.cast(f)
 
-        self._amp = self._math.cast(kwargs.pop("amp", np.ones(self.nw)))
+        self._amp = self._math.cast(kwargs.get("amp", np.ones(self.nw)))
 
         # Reset data and priors
         self._flux = None
@@ -303,12 +369,6 @@ def reset(self, **kwargs):
         self._L = None
         self._solution = None
 
-        # Check for bad kwargs, with the following exceptions
-        kwargs.pop("source_npts", None)
-        kwargs.pop("dr_oversample", None)
-        kwargs.pop("dr_lam", None)
-        self._check_kwargs("reset", kwargs)
-
     def show(self, **kwargs):
         """
         Display an image of the map, with optional animation. See the
@@ -369,108 +429,37 @@ def show(self, **kwargs):
 
         """
         # Get kwargs
-        get_val = evaluator(**kwargs)
-        cmap = kwargs.pop("cmap", "plasma")
-        grid = kwargs.pop("grid", True)
-        interval = kwargs.pop("interval", 75)
-        file = kwargs.pop("file", None)
-        html5_video = kwargs.pop("html5_video", True)
-        norm = kwargs.pop("norm", None)
-        dpi = kwargs.pop("dpi", None)
-        figsize = kwargs.pop("figsize", None)
-        bitrate = kwargs.pop("bitrate", None)
-        colorbar = kwargs.pop("colorbar", False)
-        ax = kwargs.pop("ax", None)
+        cmap = kwargs.get("cmap", "plasma")
+        grid = kwargs.get("grid", True)
+        interval = kwargs.get("interval", 75)
+        file = kwargs.get("file", None)
+        html5_video = kwargs.get("html5_video", True)
+        dpi = kwargs.get("dpi", None)
+        figsize = kwargs.get("figsize", None)
+        bitrate = kwargs.get("bitrate", None)
+        colorbar = kwargs.get("colorbar", False)
+        ax = kwargs.get("ax", None)
         if ax is None:
             custom_ax = False
         else:
             custom_ax = True
 
-        # Ylm-base maps only
-        if not self.__props__["limbdarkened"]:
-
-            projection = get_projection(kwargs.get("projection", "ortho"))
-
-            # Get the map orientation
-            if self.lazy:
-                inc = get_val(self._inc)
-                obl = get_val(self._obl)
-            else:
-                inc = self._inc
-                obl = self._obl
-
-            # Get the rotational phase
-            if self.lazy:
-                theta = get_val(
-                    self._math.vectorize(
-                        self._math.cast(kwargs.pop("theta", 0.0))
-                        * self._angle_factor
-                    )
-                )
-            else:
-                theta = np.atleast_1d(
-                    np.array(kwargs.pop("theta", 0.0)) * self._angle_factor
-                )
-
-        else:
-
-            inc = np.array(0.5 * np.pi)
-            obl = np.array(0)
-            theta = np.array([0])
-
         # Render the map if needed
-        image = kwargs.pop("image", None)  # undocumented, used internally
-        illum = kwargs.pop("illum", None)  # undocumented, used internally
+        image = kwargs.get("image", None)  # undocumented, used internally
+        overlay = kwargs.get("overlay", None)  # undocumented, used internally
         if image is None:
-
-            # We need to evaluate the variables so we can plot the map!
-            if self.lazy:
-
-                # Get kwargs
-                res = kwargs.pop("res", 300)
-
-                # Evaluate the variables
-                u = get_val(self._u)
-
-                if not self.__props__["limbdarkened"]:
-
-                    inc = get_val(self._inc)
-                    obl = get_val(self._obl)
-                    y = get_val(self._y)
-                    f = get_val(self._f)
-
-                    # Explicitly call the compiled version of `render`
-                    image = get_val(self._amp).reshape(
-                        -1, 1, 1
-                    ) * self.ops.render(
-                        res, projection, theta, inc, obl, y, u, f
-                    )
-
-                else:
-
-                    # Explicitly call the compiled version of `render`
-                    image = get_val(self._amp).reshape(
-                        -1, 1, 1
-                    ) * self.ops.render_ld(res, u)
-
-            else:
-
-                # Easy!
-                if not self.__props__["limbdarkened"]:
-                    image = self.render(
-                        theta=theta / self._angle_factor, **kwargs
-                    )
-                else:
-                    image = self.render(**kwargs)
-                kwargs.pop("res", None)
-
+            image = self._render_greedy(**kwargs)
         if len(image.shape) == 3:
             nframes = image.shape[0]
         else:
             nframes = 1
             image = np.reshape(image, (1,) + image.shape)
-            if illum is not None:
-                illum = np.reshape(illum, (1,) + illum.shape)
+            if overlay is not None:
+                overlay = np.reshape(overlay, (1,) + overlay.shape)
+
+        # Additional kwargs
+        projection = self._get_projection(**kwargs)
+        norm = self._get_norm(image, **kwargs)
 
         # Animation
         animated = nframes > 1
@@ -478,10 +467,9 @@ def show(self, **kwargs):
         latlines = []
         lonlines = []
 
-        if (
-            not self.__props__["limbdarkened"]
-            and projection != STARRY_ORTHOGRAPHIC_PROJECTION
-        ):
+        # Set up the figure
+        if projection != STARRY_ORTHOGRAPHIC_PROJECTION:
+
             # Set up the plot
             if figsize is None:
                 figsize = (7, 3.75)
@@ -586,53 +574,31 @@ def show(self, **kwargs):
             extent = (-1 - dx, 1, -1 - dx, 1)
 
             # Anti-aliasing at the edges
-            x = np.linspace(-1, 1, 10000)
-            y = np.sqrt(1 - x ** 2)
-            borders += [ax.fill_between(x, 1.1 * y, y, color="w", zorder=-1)]
-            borders += [ax.fill_betweenx(x, 1.1 * y, y, color="w", zorder=-1)]
-            borders += [ax.fill_between(x, -1.1 * y, -y, color="w", zorder=-1)]
+            xp, yp, xm, ym = self._get_ortho_borders()
             borders += [
-                ax.fill_betweenx(x, -1.1 * y, -y, color="w", zorder=-1)
+                ax.fill_between(xp, 1.1 * yp, yp, color="w", zorder=-1)
+            ]
+            borders += [
+                ax.fill_between(xm, 1.1 * ym, ym, color="w", zorder=-1)
             ]
 
             # Grid lines
             if grid:
-                x = np.linspace(-1, 1, 10000)
-                y = np.sqrt(1 - x ** 2)
-                borders += ax.plot(x, y, "k-", alpha=1, lw=1.5, zorder=0)
-                borders += ax.plot(x, -y, "k-", alpha=1, lw=1.5, zorder=0)
-                lats = get_ortho_latitude_lines(inc=inc, obl=obl)
+                borders += ax.plot(xp, yp, "k-", alpha=1, lw=1.5, zorder=0)
+                borders += ax.plot(xm, ym, "k-", alpha=1, lw=1.5, zorder=0)
+                lats = self._get_ortho_latitude_lines(**kwargs)
                 latlines = [None for n in lats]
                 for n, l in enumerate(lats):
                     (latlines[n],) = ax.plot(
                         l[0], l[1], "k-", lw=0.5, alpha=0.5, zorder=0
                     )
-                lons = get_ortho_longitude_lines(
-                    inc=inc, obl=obl, theta=theta[0]
-                )
+                lons = self._get_ortho_longitude_lines(**kwargs)
                 lonlines = [None for n in lons]
                 for n, l in enumerate(lons):
                     (lonlines[n],) = ax.plot(
                         l[0], l[1], "k-", lw=0.5, alpha=0.5, zorder=0
                     )
 
-        # Plot the first frame of the image
-        if norm is None or norm == "rv":
-            vmin = np.nanmin(image)
-            vmax = np.nanmax(image)
-            # Set a minimum contrast
-            if np.abs(vmin - vmax) < 1e-12:
-                vmin -= 1e-12
-                vmax += 1e-12
-            if norm is None:
-                norm = colors.Normalize(vmin=vmin, vmax=vmax)
-            elif norm == "rv":
-                try:
-                    norm = colors.TwoSlopeNorm(vmin=vmin, vcenter=0, vmax=vmax)
-                except AttributeError:  # pragma: no cover
-                    # TwoSlopeNorm was introduced in matplotlib 3.2
-                    norm = colors.Normalize(vmin=vmin, vmax=vmax)
-
         img = ax.imshow(
             image[0],
             origin="lower",
@@ -644,19 +610,19 @@ def show(self, **kwargs):
             zorder=-3,
         )
 
-        if illum is not None:
-
-            # Apply the transparency filter for the illumination
+        if overlay is None:
+            img_overlay = None
+        else:
             cmapI = colors.LinearSegmentedColormap.from_list(
-                "illum", ["k", "k"], 256
+                "overlay", ["k", "k"], 256
             )
             cmapI._init()
             alphas = np.linspace(1.0, 0.0, cmapI.N + 3)
             cmapI._lut[:, -1] = alphas
             cmapI.set_under((0, 0, 0, 1))
             cmapI.set_over((0, 0, 0, 0))
-            img_illum = ax.imshow(
-                illum[0],
+            img_overlay = ax.imshow(
+                overlay[0],
                 origin="lower",
                 extent=extent,
                 cmap=cmapI,
@@ -678,30 +644,21 @@ def show(self, **kwargs):
         # Display or save the image / animation
         if animated:
 
-            def updatefig(i):
-                img.set_array(image[i])
-                images = [img]
-                if illum is not None:
-                    img_illum.set_array(illum[i])
-                    images += [img_illum]
-                if (
-                    not self.__props__["limbdarkened"]
-                    and projection == STARRY_ORTHOGRAPHIC_PROJECTION
-                    and grid
-                    and len(theta) > 1
-                    and self.nw is None
-                ):
-                    lons = get_ortho_longitude_lines(
-                        inc=inc, obl=obl, theta=theta[i]
-                    )
-                    for n, l in enumerate(lons):
-                        lonlines[n].set_xdata(l[0])
-                        lonlines[n].set_ydata(l[1])
-                return tuple(images + lonlines + latlines + borders)
-
             ani = FuncAnimation(
                 fig,
-                updatefig,
+                self._updatefig,
+                fargs=(
+                    img,
+                    image,
+                    img_overlay,
+                    overlay,
+                    projection,
+                    grid,
+                    lonlines,
+                    latlines,
+                    borders,
+                    kwargs,
+                ),
                 interval=interval,
                 blit=True,
                 frames=image.shape[0],
@@ -766,19 +723,6 @@ def updatefig(i):
             elif not custom_ax:
                 plt.show()
 
-        # Check for invalid kwargs
-        kwargs.pop("point", None)
-        kwargs.pop("model", None)
-        if self.__props__["rv"]:
-            kwargs.pop("rv", None)
-        if not self.__props__["limbdarkened"]:
-            kwargs.pop("projection", None)
-        if self.__props__["reflected"]:
-            kwargs.pop("xs", None)
-            kwargs.pop("ys", None)
-            kwargs.pop("zs", None)
-        self._check_kwargs("show", kwargs)
-
     def limbdark_is_physical(self):
         """Check whether the limb darkening profile (if any) is physical.
 
@@ -1030,14 +974,56 @@ class YlmBase(legacy.YlmBase):
 
     _ops_class_ = OpsYlm
 
+    def _render_greedy(self, **kwargs):
+        get_val = evaluator(**kwargs)
+        amp = get_val(self.amp).reshape(-1, 1, 1)
+        return amp * self.ops.render(
+            kwargs.get("res", 300),
+            get_projection(kwargs.get("projection", "ortho")),
+            np.atleast_1d(get_val(kwargs.get("theta", 0.0)))
+            * self._angle_factor,
+            get_val(self._inc),
+            get_val(self._obl),
+            get_val(self._y),
+            get_val(self._u),
+            get_val(self._f),
+        )
+
+    def _get_ortho_latitude_lines(self, **kwargs):
+        get_val = evaluator(**kwargs)
+        inc = get_val(self._inc)
+        obl = get_val(self._obl)
+        return get_ortho_latitude_lines(inc=inc, obl=obl)
+
+    def _get_ortho_longitude_lines(self, i=0, **kwargs):
+        get_val = evaluator(**kwargs)
+        inc = get_val(self._inc)
+        obl = get_val(self._obl)
+        theta = (
+            np.atleast_1d(get_val(kwargs.get("theta", 0.0)))
+            * self._angle_factor
+        )
+        return get_ortho_longitude_lines(inc=inc, obl=obl, theta=theta[i])
+
+    def _get_flux_kwargs(self, kwargs):
+        xo = kwargs.get("xo", 0.0)
+        yo = kwargs.get("yo", 0.0)
+        zo = kwargs.get("zo", 1.0)
+        ro = kwargs.get("ro", 0.0)
+        theta = kwargs.get("theta", 0.0)
+        theta, xo, yo, zo = self._math.vectorize(theta, xo, yo, zo)
+        theta, xo, yo, zo, ro = self._math.cast(theta, xo, yo, zo, ro)
+        theta *= self._angle_factor
+        return theta, xo, yo, zo, ro
+
     def reset(self, **kwargs):
         if kwargs.get("inc", None) is not None:
-            self.inc = kwargs.pop("inc")
+            self.inc = kwargs.get("inc")
         else:
             self._inc = self._math.cast(0.5 * np.pi)
 
         if kwargs.get("obl", None) is not None:
-            self.obl = kwargs.pop("obl")
+            self.obl = kwargs.get("obl")
         else:
             self._obl = self._math.cast(0.0)
 
@@ -1083,9 +1069,6 @@ def design_matrix(self, **kwargs):
         # Orbital kwargs
         theta, xo, yo, zo, ro = self._get_flux_kwargs(kwargs)
 
-        # Check for invalid kwargs
-        self._check_kwargs("design_matrix", kwargs)
-
         # Compute & return
         return self.ops.X(
             theta, xo, yo, zo, ro, self._inc, self._obl, self._u, self._f
@@ -1133,15 +1116,15 @@ def flux(self, **kwargs):
                 this body's radius.
             theta (scalar or vector, optional): Angular phase of the body
                 in units of :py:attr:`angle_unit`.
+            integrated (bool, optional): If True, dots the flux with the
+                amplitude. Default False, in which case this returns a
+                2d array (wavelength-dependent maps only).
         """
         # Orbital kwargs
         theta, xo, yo, zo, ro = self._get_flux_kwargs(kwargs)
 
-        # Check for invalid kwargs
-        self._check_kwargs("flux", kwargs)
-
         # Compute & return
-        return self.amp * self.ops.flux(
+        flux = self.ops.flux(
             theta,
             xo,
             yo,
@@ -1154,6 +1137,11 @@ def flux(self, **kwargs):
             self._f,
         )
 
+        if kwargs.get("integrated", False):
+            return self._math.dot(flux, self.amp)
+        else:
+            return self.amp * flux
+
     def intensity(self, lat=0, lon=0, **kwargs):
         """
         Compute and return the intensity of the map.
@@ -1180,14 +1168,11 @@ def intensity(self, lat=0, lon=0, **kwargs):
         theta *= self._angle_factor
 
         # If limb-darkened, allow a `limbdarken` keyword
-        if kwargs.pop("limbdarken", True) and self.udeg > 0:
+        if kwargs.get("limbdarken", True) and self.udeg > 0:
             ld = np.array(True)
         else:
             ld = np.array(False)
 
-        # Check for invalid kwargs
-        self._check_kwargs("intensity", kwargs)
-
         # Compute & return
         return self.amp * self.ops.intensity(
             lat, lon, self._y, self._u, self._f, theta, ld
@@ -1415,7 +1400,7 @@ def load(
 
             # Scale the coeffs?
             if I < 0:
-                fac = self._amp / (self._amp - np.pi * I)
+                fac = self.amp / (self.amp - np.pi * I)
                 if self.lazy:
                     self._y *= fac
                     self._y = self.ops.set_map_vector(self._y, 0, 1.0)
@@ -1579,14 +1564,14 @@ def minimize(self, oversample=1, ntries=1, bounds=None, return_info=False):
             return (
                 lat / self._angle_factor,
                 lon / self._angle_factor,
-                self._amp * I,
+                self.amp * I,
                 self.ops._minimize.result,
             )
         else:
             return (
                 lat / self._angle_factor,
                 lon / self._angle_factor,
-                self._amp * I,
+                self.amp * I,
             )
 
     def get_pixel_transforms(self, oversample=2, lam=1e-6, eps=1e-6):
@@ -1784,6 +1769,38 @@ class LimbDarkenedBase(object):
 
     _ops_class_ = OpsLD
 
+    def _render_greedy(self, **kwargs):
+        get_val = evaluator(**kwargs)
+        amp = get_val(self.amp).reshape(-1, 1, 1)
+        return amp.reshape(-1, 1, 1) * self.ops.render_ld(
+            kwargs.get("res", 300), get_val(self._u)
+        )
+
+    def _get_projection(self, **kwargs):
+        projection = get_projection(kwargs.get("projection", "ortho"))
+        assert (
+            projection == STARRY_ORTHOGRAPHIC_PROJECTION
+        ), "Only orthographic projection is supported for limb-darkened maps."
+        return projection
+
+    def _get_ortho_latitude_lines(self, **kwargs):
+        return get_ortho_latitude_lines()
+
+    def _get_ortho_longitude_lines(self, i=0, **kwargs):
+        return get_ortho_longitude_lines()
+
+    def _updatefig(i):
+        return ()
+
+    def _get_flux_kwargs(self, kwargs):
+        xo = kwargs.get("xo", 0.0)
+        yo = kwargs.get("yo", 0.0)
+        zo = kwargs.get("zo", 1.0)
+        ro = kwargs.get("ro", 0.0)
+        xo, yo, zo = self._math.vectorize(xo, yo, zo)
+        xo, yo, zo, ro = self._math.cast(xo, yo, zo, ro)
+        return xo, yo, zo, ro
+
     def flux(self, **kwargs):
         """
         Compute and return the light curve.
@@ -1799,14 +1816,7 @@ def flux(self, **kwargs):
                 this body's radius.
         """
         # Orbital kwargs
-        theta = kwargs.pop("theta", None)
-        _, xo, yo, zo, ro = self._get_flux_kwargs(kwargs)
-
-        # Check for invalid kwargs
-        if theta is not None:
-            # If the user passed in `theta`, make sure a warning is raised
-            kwargs["theta"] = theta
-        self._check_kwargs("flux", kwargs)
+        xo, yo, zo, ro = self._get_flux_kwargs(kwargs)
 
         # Compute & return
         return self.amp * self.ops.flux(xo, yo, zo, ro, self._u)
@@ -1886,10 +1896,26 @@ class RVBase(object):
 
     _ops_class_ = OpsRV
 
+    def _get_norm(self, image, **kwargs):
+        norm = kwargs.get("norm", None)
+        if norm is None:
+            vmin = np.nanmin(image)
+            vmax = np.nanmax(image)
+            if np.abs(vmin - vmax) < 1e-9:
+                vmin -= 1e-9
+                vmax += 1e-9
+            try:
+                norm = colors.TwoSlopeNorm(vmin=vmin, vcenter=0, vmax=vmax)
+            except (AttributeError, ValueError):  # pragma: no cover
+                # TwoSlopeNorm was introduced in matplotlib 3.2
+                norm = colors.Normalize(vmin=vmin, vmax=vmax)
+        else:
+            return norm
+
     def reset(self, **kwargs):
-        self.velocity_unit = kwargs.pop("velocity_unit", units.m / units.s)
-        self.veq = kwargs.pop("veq", 0.0)
-        self.alpha = kwargs.pop("alpha", 0.0)
+        self.velocity_unit = kwargs.get("velocity_unit", units.m / units.s)
+        self.veq = kwargs.get("veq", 0.0)
+        self.alpha = kwargs.get("alpha", 0.0)
         super(RVBase, self).reset(**kwargs)
 
     @property
@@ -1983,9 +2009,6 @@ def rv(self, **kwargs):
         # Orbital kwargs
         theta, xo, yo, zo, ro = self._get_flux_kwargs(kwargs)
 
-        # Check for invalid kwargs
-        self._check_kwargs("rv", kwargs)
-
         # Compute
         return self.ops.rv(
             theta,
@@ -2020,7 +2043,7 @@ def intensity(self, **kwargs):
 
         """
         # Compute the velocity-weighted intensity if `rv==True`
-        rv = kwargs.pop("rv", True)
+        rv = kwargs.get("rv", True)
         if rv:
             self._set_RV_filter()
         res = super(RVBase, self).intensity(**kwargs)
@@ -2056,7 +2079,7 @@ def render(self, **kwargs):
         # Render the velocity map if `rv==True`
         # Override the `projection` kwarg if we're
         # plotting the radial velocity.
-        rv = kwargs.pop("rv", True)
+        rv = kwargs.get("rv", True)
         if rv:
             kwargs.pop("projection", None)
             self._set_RV_filter()
@@ -2072,8 +2095,7 @@ def show(self, rv=True, **kwargs):
         if rv:
             kwargs.pop("projection", None)
             self._set_RV_filter()
-            kwargs["cmap"] = kwargs.pop("cmap", "RdBu_r")
-            kwargs["norm"] = kwargs.pop("norm", "rv")
+            kwargs["cmap"] = kwargs.get("cmap", "RdBu_r")
         res = super(RVBase, self).show(rv=rv, **kwargs)
         if rv:
             self._unset_RV_filter()
@@ -2120,7 +2142,7 @@ class ReflectedBase(object):
     _ops_class_ = OpsReflected
 
     def reset(self, **kwargs):
-        self.roughness = kwargs.pop("roughness", self._math.cast(0.0))
+        self.roughness = kwargs.get("roughness", self._math.cast(0.0))
         super(ReflectedBase, self).reset(**kwargs)
 
     @property
@@ -2146,15 +2168,15 @@ def roughness(self, value):
         self._sigr = self._math.cast(value) * self._angle_factor
 
     def _get_flux_kwargs(self, kwargs):
-        xo = kwargs.pop("xo", 0.0)
-        yo = kwargs.pop("yo", 0.0)
-        zo = kwargs.pop("zo", 1.0)
-        ro = kwargs.pop("ro", 0.0)
-        xs = kwargs.pop("xs", 0.0)
-        ys = kwargs.pop("ys", 0.0)
-        zs = kwargs.pop("zs", 1.0)
-        Rs = kwargs.pop("rs", 0.0)
-        theta = kwargs.pop("theta", 0.0)
+        xo = kwargs.get("xo", 0.0)
+        yo = kwargs.get("yo", 0.0)
+        zo = kwargs.get("zo", 1.0)
+        ro = kwargs.get("ro", 0.0)
+        xs = kwargs.get("xs", 0.0)
+        ys = kwargs.get("ys", 0.0)
+        zs = kwargs.get("zs", 1.0)
+        Rs = kwargs.get("rs", 0.0)
+        theta = kwargs.get("theta", 0.0)
         theta, xs, ys, zs, xo, yo, zo = self._math.vectorize(
             theta, xs, ys, zs, xo, yo, zo
         )
@@ -2196,9 +2218,6 @@ def design_matrix(self, **kwargs):
         # Orbital kwargs
         theta, xs, ys, zs, Rs, xo, yo, zo, ro = self._get_flux_kwargs(kwargs)
 
-        # Check for invalid kwargs
-        self._check_kwargs("X", kwargs)
-
         # Compute & return
         return self.ops.X(
             theta,
@@ -2240,16 +2259,16 @@ def flux(self, **kwargs):
                 this body's radius.
             theta (scalar or vector, optional): Angular phase of the body
                 in units of :py:attr:`angle_unit`.
+            integrated (bool, optional): If True, dots the flux with the
+                amplitude. Default False, in which case this returns a
+                2d array (wavelength-dependent maps only).
 
         """
         # Orbital kwargs
         theta, xs, ys, zs, Rs, xo, yo, zo, ro = self._get_flux_kwargs(kwargs)
 
-        # Check for invalid kwargs
-        self._check_kwargs("flux", kwargs)
-
         # Compute & return
-        return self.amp * self.ops.flux(
+        flux = self.ops.flux(
             theta,
             xs,
             ys,
@@ -2267,6 +2286,11 @@ def flux(self, **kwargs):
             self._sigr,
         )
 
+        if kwargs.get("integrated", False):
+            return self._math.dot(flux, self.amp)
+        else:
+            return self.amp * flux
+
     def intensity(
         self,
         lat=0,
@@ -2458,17 +2482,17 @@ def show(self, **kwargs):
 
         # Get kwargs
         get_val = evaluator(**kwargs)
-        res = kwargs.pop("res", 300)
+        res = kwargs.get("res", 300)
         projection = get_projection(kwargs.get("projection", "ortho"))
-        theta = self._math.cast(kwargs.pop("theta", 0.0)) * self._angle_factor
-        xs = self._math.cast(kwargs.pop("xs", 0))
-        ys = self._math.cast(kwargs.pop("ys", 0))
-        zs = self._math.cast(kwargs.pop("zs", 1))
-        Rs = self._math.cast(kwargs.pop("rs", 0))
+        theta = self._math.cast(kwargs.get("theta", 0.0)) * self._angle_factor
+        xs = self._math.cast(kwargs.get("xs", 0))
+        ys = self._math.cast(kwargs.get("ys", 0))
+        zs = self._math.cast(kwargs.get("zs", 1))
+        Rs = self._math.cast(kwargs.get("rs", 0))
         theta, xs, ys, zs = self._math.vectorize(theta, xs, ys, zs)
-        illuminate = int(kwargs.pop("illuminate", True))
-        on94_exact = int(kwargs.pop("on94_exact", False))
-        screen = bool(kwargs.pop("screen", True))
+        illuminate = int(kwargs.get("illuminate", True))
+        on94_exact = int(kwargs.get("on94_exact", False))
+        screen = bool(kwargs.get("screen", True))
 
         if self.nw is None:
             amp = self.amp
@@ -2477,27 +2501,19 @@ def show(self, **kwargs):
             # so we must reshape `amp` to take the product correctly
             amp = self.amp[:, np.newaxis, np.newaxis]
 
-        if self.lazy:
-            # Evaluate the variables
-            theta = get_val(theta)
-            xs = get_val(xs)
-            ys = get_val(ys)
-            zs = get_val(zs)
-            Rs = get_val(Rs)
-            inc = get_val(self._inc)
-            obl = get_val(self._obl)
-            y = get_val(self._y)
-            u = get_val(self._u)
-            f = get_val(self._f)
-            sigr = get_val(self._sigr)
-            amp = get_val(amp)
-        else:
-            inc = self._inc
-            obl = self._obl
-            y = self._y
-            u = self._u
-            f = self._f
-            sigr = self._sigr
+        # Evaluate the variables
+        theta = get_val(theta)
+        xs = get_val(xs)
+        ys = get_val(ys)
+        zs = get_val(zs)
+        Rs = get_val(Rs)
+        inc = get_val(self._inc)
+        obl = get_val(self._obl)
+        y = get_val(self._y)
+        u = get_val(self._u)
+        f = get_val(self._f)
+        sigr = get_val(self._sigr)
+        amp = get_val(amp)
 
         if screen and illuminate:
 
@@ -2542,7 +2558,7 @@ def show(self, **kwargs):
             )
             if np.nanmax(illum) > 0:
                 illum /= np.nanmax(illum)
-            kwargs["illum"] = illum
+            kwargs["overlay"] = illum
 
         else:
 
@@ -2570,13 +2586,345 @@ def show(self, **kwargs):
         return super(ReflectedBase, self).show(**kwargs)
 
 
+class OblateAmplitude(Amplitude):
+    def __get__(self, instance, owner):
+        return instance._amp * instance._norm
+
+
+class OblateBase(object):
+    """The oblate ``starry`` map class.
+
+    This class handles maps of oblate spheroids, such as fast-rotating
+    stars and planets or bodies with (moderate) tidal distortion.
+    It has all the same attributes and methods as :py:class:`starry.maps.YlmBase`,
+    with the additions and modifications listed below.
+
+    .. note::
+        Instantiate this class by calling :py:func:`starry.Map` with
+        ``oblate`` set to True.
+    """
+
+    _ops_class_ = OpsOblate
+
+    # Special handling for the map amplitude
+    amp = OblateAmplitude()
+
+    # Constants for Planck law (mks)
+    _twohcsq = 1.19104295e-16
+    _hcdivkb = 1.43877735e-2
+
+    @property
+    def _norm(self):
+        # Planck's law normalization
+        pnorm = self._twohcsq / self._wav ** 5
+
+        if self.__props__.get("normalized", True):
+
+            # Normalize by the flux from a uniform map
+            if self.nw is None:
+                uniform_y = np.zeros(self.Ny)
+                uniform_y[0] = 1.0
+                uniform_y = self._math.cast(uniform_y)
+            else:
+                uniform_y = np.zeros((self.Ny, self.nw))
+                uniform_y[0, :] = 1.0
+                uniform_y = self._math.cast(uniform_y)
+            flux0 = self.ops.flux(
+                self._math.cast([0.0]),  # theta
+                self._math.cast([0.0]),  # xo
+                self._math.cast([0.0]),  # yo
+                self._math.cast([0.0]),  # zo
+                self._math.cast(0.0),  # ro
+                self._inc,
+                self._obl,
+                self.fproj,
+                uniform_y,
+                self._u,
+                self._f,
+            )[0]
+            return pnorm / self._math.dot(pnorm * flux0, self._amp)
+
+        else:
+
+            # Trivial
+            return pnorm
+
+    def render(self, res=300, projection="ortho", theta=0.0):
+        # Multiple frames?
+        if self.nw is not None:
+            animated = True
+        else:
+            if is_tensor(theta):
+                animated = hasattr(theta, "ndim") and theta.ndim > 0
+            else:
+                animated = hasattr(theta, "__len__")
+
+        # Convert
+        projection = get_projection(projection)
+        theta = self._math.vectorize(
+            self._math.cast(theta) * self._angle_factor
+        )
+
+        # Compute
+        if self.nw is None or self.lazy:
+            amp = self.amp
+        else:
+            # The intensity has shape `(nw, res, res)`
+            # so we must reshape `amp` to take the product correctly
+            amp = self.amp[:, np.newaxis, np.newaxis]
+        image = amp * self.ops.render(
+            res,
+            projection,
+            theta,
+            self._inc,
+            self._obl,
+            self.fproj,
+            self._y,
+            self._u,
+            self._f,
+        )
+
+        # Squeeze?
+        if animated:
+            return image
+        else:
+            return self._math.reshape(image, [res, res])
+
+    def _render_greedy(self, **kwargs):
+        get_val = evaluator(**kwargs)
+        amp = get_val(self.amp).reshape(-1, 1, 1)
+        image = self.ops.render(
+            kwargs.get("res", 300),
+            get_projection(kwargs.get("projection", "ortho")),
+            np.atleast_1d(get_val(kwargs.get("theta", 0.0)))
+            * self._angle_factor,
+            get_val(self._inc),
+            get_val(self._obl),
+            get_val(self.fproj),
+            get_val(self._y),
+            get_val(self._u),
+            get_val(self._f),
+        )
+        return amp * image
+
+    def _get_ortho_borders(self, **kwargs):
+        get_val = evaluator(**kwargs)
+        obl = get_val(self._obl)
+        fproj = get_val(self.fproj)
+        xp = np.linspace(-1, 1, 10000)
+        yp = np.sqrt(1 - xp ** 2)
+        xm = xp
+        ym = -yp
+        xpr = xp * np.cos(obl) - yp * (1 - fproj) * np.sin(obl)
+        ypr = xp * np.sin(obl) + yp * (1 - fproj) * np.cos(obl)
+        xmr = xm * np.cos(obl) - ym * (1 - fproj) * np.sin(obl)
+        ymr = xm * np.sin(obl) + ym * (1 - fproj) * np.cos(obl)
+        return xpr, ypr, xmr, ymr
+
+    def _get_ortho_latitude_lines(self, **kwargs):
+        get_val = evaluator(**kwargs)
+        inc = get_val(self._inc)
+        obl = get_val(self._obl)
+        fproj = get_val(self.fproj)
+        return get_ortho_latitude_lines(inc=inc, obl=obl, fproj=fproj)
+
+    def _get_ortho_longitude_lines(self, i=0, **kwargs):
+        get_val = evaluator(**kwargs)
+        inc = get_val(self._inc)
+        obl = get_val(self._obl)
+        fproj = get_val(self.fproj)
+        theta = (
+            np.atleast_1d(get_val(kwargs.get("theta", 0.0)))
+            * self._angle_factor
+        )
+        return get_ortho_longitude_lines(
+            inc=inc, obl=obl, fproj=fproj, theta=theta[i]
+        )
+
+    def reset(self, **kwargs):
+        self._fobl = self._math.cast(kwargs.get("f", 0.0))
+        self._omega = self._math.cast(kwargs.get("omega", 0.0))
+        self._tpole = self._math.cast(kwargs.get("tpole", 5800.0))
+        self._beta = self._math.cast(kwargs.get("beta", 0.23))
+        super(OblateBase, self).reset(**kwargs)
+        self._set_grav_dark_filter()
+
+    @property
+    def f(self):
+        """The oblateness of the spheroid.
+
+        This is the ratio of the difference between the equatorial and
+        polar radii to the equatorial radius, and must be in the range
+        ``[0, 1)``.
+        """
+        return self._fobl
+
+    @f.setter
+    def f(self, value):
+        self._fobl = self._math.cast(value)
+        self._set_grav_dark_filter()
+
+    @property
+    def omega(self):
+        """
+        The angular rotational velocity of the body as a fraction of the breakup velocity.
+
+        """
+        return self._omega
+
+    @omega.setter
+    def omega(self, value):
+        self._omega = self._math.cast(value)
+        self._set_grav_dark_filter()
+
+    @property
+    def fproj(self):
+        """The projected oblateness at the current inclination."""
+        return 1 - self._math.sqrt(
+            1 - self._fobl * (2 - self._fobl) * self._math.sin(self._inc) ** 2
+        )
+
+    @property
+    def tpole(self):
+        """The polar temperature of the star in K."""
+        return self._tpole
+
+    @tpole.setter
+    def tpole(self, value):
+        self._tpole = self._math.cast(value)
+        self._set_grav_dark_filter()
+
+    @property
+    def beta(self):
+        """The von Zeipel law coefficient."""
+        return self._beta
+
+    @beta.setter
+    def beta(self, value):
+        self._beta = self._math.cast(value)
+        self._set_grav_dark_filter()
+
+    @property
+    def wav(self):
+        """
+        The wavelength(s) at which the flux is measured in units of ``wav_unit``.
+
+        """
+        return self._wav / self._wav_factor
+
+    @wav.setter
+    def wav(self, value):
+        value = self._math.cast(value) * self._wav_factor
+        if self._nw is None:
+            assert value.ndim == 0, "Wavelength must be a scalar."
+            self._wav = value
+        else:
+            assert value.ndim == 1, "Wavelength must be a vector."
+            self._wav = value * self._math.ones(self._nw)
+        self._set_grav_dark_filter()
+
+    def _set_grav_dark_filter(self):
+        if self._fdeg > 0:
+            self._f = self.ops.compute_grav_dark_filter(
+                self._wav / self._hcdivkb,
+                self._omega,
+                self._fobl,
+                self._beta,
+                self._tpole,
+            )
+
+    def design_matrix(self, **kwargs):
+        r"""Compute and return the light curve design matrix :math:`A`.
+
+        This matrix is used to compute the flux :math:`f` from a vector of spherical
+        harmonic coefficients :math:`y` and the map amplitude :math:`a`:
+        :math:`f = a A y`.
+
+        Args:
+            xo (scalar or vector, optional): x coordinate of the occultor
+                relative to this body in units of this body's radius.
+            yo (scalar or vector, optional): y coordinate of the occultor
+                relative to this body in units of this body's radius.
+            zo (scalar or vector, optional): z coordinate of the occultor
+                relative to this body in units of this body's radius.
+            ro (scalar, optional): Radius of the occultor in units of
+                this body's radius.
+            theta (scalar or vector, optional): Angular phase of the body
+                in units of :py:attr:`angle_unit`.
+        """
+        # The problem isn't quite linear for spectral maps
+        self._no_spectral()
+
+        # Orbital kwargs
+        theta, xo, yo, zo, ro = self._get_flux_kwargs(kwargs)
+
+        # Compute & return
+        return self.ops.X(
+            theta,
+            xo,
+            yo,
+            zo,
+            ro,
+            self._inc,
+            self._obl,
+            self.fproj,
+            self._u,
+            self._f,
+        )
+
+    def flux(self, **kwargs):
+        """
+        Compute and return the flux from the map.
+
+        Args:
+            xo (scalar or vector, optional): x coordinate of the occultor
+                relative to this body in units of this body's radius.
+            yo (scalar or vector, optional): y coordinate of the occultor
+                relative to this body in units of this body's radius.
+            zo (scalar or vector, optional): z coordinate of the occultor
+                relative to this body in units of this body's radius.
+            ro (scalar, optional): Radius of the occultor in units of
+                this body's radius.
+            theta (scalar or vector, optional): Angular phase of the body
+                in units of :py:attr:`angle_unit`.
+            integrated (bool, optional): If True, dots the flux with the
+                amplitude. Default False, in which case this returns a
+                2d array (wavelength-dependent maps only).
+        """
+        # Orbital kwargs
+        theta, xo, yo, zo, ro = self._get_flux_kwargs(kwargs)
+
+        # Compute & return
+        flux = self.ops.flux(
+            theta,
+            xo,
+            yo,
+            zo,
+            ro,
+            self._inc,
+            self._obl,
+            self.fproj,
+            self._y,
+            self._u,
+            self._f,
+        )
+
+        if kwargs.get("integrated", False):
+            return self._math.dot(flux, self.amp)
+        else:
+            return self.amp * flux
+
+
 def Map(
     ydeg=0,
     udeg=0,
     nw=None,
     rv=False,
     reflected=False,
+    oblate=False,
     source_npts=1,
+    gdeg=4,
+    normalized=True,
     lazy=None,
     **kwargs
 ):
@@ -2585,26 +2933,38 @@ def Map(
     This function is a class factory that returns either
     a :doc:`spherical harmonic map <SphericalHarmonicMap>`,
     a :doc:`limb darkened map <LimbDarkenedMap>`,
-    a :doc:`radial velocity map <RadialVelocityMap>`, or
-    a :doc:`reflected light map <ReflectedLightMap>`,
+    a :doc:`radial velocity map <RadialVelocityMap>`,
+    a :doc:`reflected light map <ReflectedLightMap>`, or
+    an :doc:`oblate map <OblateMap>`
     depending on the arguments provided by the user. The default is
     a :doc:`spherical harmonic map <SphericalHarmonicMap>`. If ``rv`` is True,
-    instantiates a :doc:`radial velocity map <RadialVelocityMap>` map, and
+    instantiates a :doc:`radial velocity map <RadialVelocityMap>` map,
     if ``reflected`` is True, instantiates a :doc:`reflected light map
-    <ReflectedLightMap>`. Otherwise, if ``ydeg`` is zero, instantiates a
-    :doc:`limb darkened map <LimbDarkenedMap>`.
+    <ReflectedLightMap>`, or if ``oblate`` is True, instantiates an
+    :doc:`oblate map <OblateMap>`. Otherwise, if ``ydeg`` is zero,
+    instantiates a :doc:`limb darkened map <LimbDarkenedMap>`.
 
     Args:
         ydeg (int, optional): Degree of the spherical harmonic map.
             Defaults to 0.
         udeg (int, optional): Degree of the limb darkening filter.
             Defaults to 0.
+        gdeg (int, optional): Degree of the gravity darkening filter.
+            Defaults to 4. Enabled only if `oblate` is True.
         nw (int, optional): Number of wavelength bins. Defaults to None
             (for monochromatic light curves).
         rv (bool, optional): If True, enable computation of radial velocities
             for modeling the Rossiter-McLaughlin effect. Defaults to False.
         reflected (bool, optional): If True, models light curves in reflected
             light. Defaults to False.
+        oblate (bool, optional): If True, models the object as an oblate
+            spheroid. Defaults to False.
+        normalized (bool, optional): Whether the flux from a uniform map should
+            be normalized to unity regardless of the projected oblateness of
+            the object (oblate maps only). Default is True. If False, the flux
+            from a uniform map is proportional to its projected area, which
+            scales as `1 / (1 - fproj)` where `fproj` is the projected
+            oblateness.
         source_npts (int, optional): Number of points used to approximate the
             finite illumination source size. Default is 1. Valid only if
             `reflected` is True.
@@ -2617,6 +2977,8 @@ def Map(
     if nw is not None:
         nw = int(nw)
         assert nw > 0, "Number of wavelength bins must be positive."
+    gdeg = int(gdeg)
+    assert gdeg >= 0, "Keyword `gdeg` must be positive."
     source_npts = int(source_npts)
     if source_npts < 1:
         source_npts = 1
@@ -2631,7 +2993,13 @@ def Map(
         )
 
     # Limb-darkened?
-    if (ydeg == 0) and (rv is False) and (reflected is False):
+    if (
+        (ydeg == 0)
+        and (udeg > 0)
+        and (rv is False)
+        and (reflected is False)
+        and (oblate is False)
+    ):
 
         # TODO: Add support for wavelength-dependent limb darkening
         if nw is not None:
@@ -2650,13 +3018,16 @@ def Map(
     elif reflected:
         Bases = (ReflectedBase,) + Bases
         fdeg = 0
+    elif oblate:
+        Bases = (OblateBase,) + Bases
+        fdeg = gdeg
     else:
         fdeg = 0
 
-    # Ensure we're not doing both
-    if rv and reflected:
+    # Ensure we're not doing a combo
+    if np.count_nonzero([rv, reflected, oblate]) > 1:
         raise NotImplementedError(
-            "Radial velocity maps not implemented in reflected light."
+            "Combinations of `rv`, `reflected`, and `oblate` not yet implemented."
         )
 
     # Construct the class
@@ -2667,8 +3038,10 @@ class Map(*Bases):
             limbdarkened=LimbDarkenedBase in Bases,
             reflected=ReflectedBase in Bases,
             rv=RVBase in Bases,
+            oblate=OblateBase in Bases,
             spectral=nw is not None,
             source_npts=source_npts,
+            normalized=normalized,
         )
 
         def __init__(self, *args, **kwargs):


Test patch:
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index e99f76130..44c04ab75 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -28,9 +28,10 @@ jobs:
           - python-version: "3.9"
             pymc-version: "pymc3==3.11"
             arviz-version: ""
-          - python-version: "3.8"
-            pymc-version: "https://github.com/pymc-devs/pymc3/archive/main.zip"
-            arviz-version: ""
+          # TODO:
+          # - python-version: "3.8"
+          #  pymc-version: "https://github.com/pymc-devs/pymc3/archive/main.zip"
+          #  arviz-version: ""
 
     steps:
       - uses: actions/checkout@v2
@@ -61,27 +62,12 @@ jobs:
         run: |
           python -m pip install -U pip
           python -m pip install ${{ matrix.pymc-version }} ${{ matrix.arviz-version }}
-          python -m pip install -e ".[tests]"
-
-      - name: Run tests (greedy)
-        if: steps.install.outcome == 'success'
-        shell: bash -l {0}
-        run: python -m pytest -v tests/greedy --junitxml=junit/test-results-greedy.xml --cov=starry --cov-append
-
-      - name: Run tests (lazy)
-        if: steps.install.outcome == 'success'
-        shell: bash -l {0}
-        run: python -m pytest -v tests/lazy --junitxml=junit/test-results-lazy.xml --cov=starry --cov-append
-
-      - name: Run extension tests (greedy)
-        if: steps.install.outcome == 'success'
-        shell: bash -l {0}
-        run: python -m pytest -v starry/extensions/tests/greedy --junitxml=junit/test-results-extensions-greedy.xml --cov=starry --cov-append
+          STARRY_UNIT_TESTS=1 python -m pip install -e ".[tests]"
 
-      - name: Run extension tests (lazy)
+      - name: Run tests
         if: steps.install.outcome == 'success'
         shell: bash -l {0}
-        run: python -m pytest -v starry/extensions/tests/lazy --junitxml=junit/test-results-extensions-lazy.xml --cov=starry --cov-append
+        run: python -m pytest -v --junitxml=junit/test-results.xml --cov=starry
 
       - name: Get unique id
         id: unique-id
@@ -97,7 +83,7 @@ jobs:
         if: always()
         with:
           files: junit/test-*.xml
-          comment_on_pr: false
+          comment_mode: off
 
       - name: Coveralls
         uses: AndreMiras/coveralls-python-action@develop
diff --git a/starry/_core/ops/lib/include/oblate/tests/conftest.py b/starry/_core/ops/lib/include/oblate/tests/conftest.py
new file mode 100644
index 000000000..91d88ac5a
--- /dev/null
+++ b/starry/_core/ops/lib/include/oblate/tests/conftest.py
@@ -0,0 +1,18 @@
+import oblate
+import pytest
+
+
+@pytest.fixture(scope="session", autouse=True)
+def solvers(lmax=5):
+    cppsolver = oblate.CppSolver(lmax)
+    numsolver = oblate.NumericalSolver(lmax, linear=False)
+    linsolver = oblate.NumericalSolver(lmax, linear=True)
+    pysolver = oblate.PythonSolver(lmax)
+    brutesolver = oblate.BruteSolver(lmax)
+    return {
+        "cpp": cppsolver,
+        "num": numsolver,
+        "lin": linsolver,
+        "py": pysolver,
+        "brute": brutesolver,
+    }
diff --git a/starry/_core/ops/lib/include/oblate/tests/derivs.py b/starry/_core/ops/lib/include/oblate/tests/derivs.py
new file mode 100644
index 000000000..cf4ab2259
--- /dev/null
+++ b/starry/_core/ops/lib/include/oblate/tests/derivs.py
@@ -0,0 +1,53 @@
+import numpy as np
+from oblate.oblate import sT, dsT
+
+deg = 2
+nruns = 1
+bo = 0.95
+ro = 0.25
+f = 0.2
+theta = 0.5
+eps = 1e-8
+
+
+Dbo, Dro, Df, Dtheta = dsT(deg, bo, ro, f, theta, nruns)
+
+Dbo_n = (
+    sT(deg, bo + eps, ro, f, theta, nruns)
+    - sT(deg, bo - eps, ro, f, theta, nruns)
+) / (2 * eps)
+
+Dro_n = (
+    sT(deg, bo, ro + eps, f, theta, nruns)
+    - sT(deg, bo, ro - eps, f, theta, nruns)
+) / (2 * eps)
+
+Df_n = (
+    sT(deg, bo, ro, f + eps, theta, nruns)
+    - sT(deg, bo, ro, f - eps, theta, nruns)
+) / (2 * eps)
+
+Dtheta_n = (
+    sT(deg, bo, ro, f, theta + eps, nruns)
+    - sT(deg, bo, ro, f, theta - eps, nruns)
+) / (2 * eps)
+
+# Case 1
+for n in [0, 1, 2, 3, 4, 5]:
+    print("n = {}".format(n))
+    for name, d, d_n in zip(
+        ["bo", "ro", "f", "theta"],
+        [Dbo, Dro, Df, Dtheta],
+        [Dbo_n, Dro_n, Df_n, Dtheta_n],
+    ):
+        diff = np.abs(d[n] - d_n[n])
+        if diff > 10 * eps:
+            print(
+                "{:7s} {:11.8f}   {:11.8f}   {:4.1e}".format(
+                    name, d[n], d_n[n], diff
+                )
+            )
+        else:
+            print("{:7s} {:11.8f}   {:11.8f}".format(name, d[n], d_n[n]))
+    print("")
+
diff --git a/starry/_core/ops/lib/include/oblate/tests/test_derivs.py b/starry/_core/ops/lib/include/oblate/tests/test_derivs.py
new file mode 100644
index 000000000..0a06508cf
--- /dev/null
+++ b/starry/_core/ops/lib/include/oblate/tests/test_derivs.py
@@ -0,0 +1,5 @@
+import oblate
+import numpy as np
+import pytest
+
+# TODO!
diff --git a/starry/_core/ops/lib/include/oblate/tests/test_oblate.py b/starry/_core/ops/lib/include/oblate/tests/test_oblate.py
new file mode 100644
index 000000000..760a913b4
--- /dev/null
+++ b/starry/_core/ops/lib/include/oblate/tests/test_oblate.py
@@ -0,0 +1,322 @@
+import oblate
+import numpy as np
+import pytest
+import matplotlib.pyplot as plt
+
+
+class Compare:
+    @pytest.mark.parametrize(
+        "compare",
+        [
+            ("cpp", "py"),
+            ("cpp", "lin"),
+            ("cpp", "num"),
+            ("cpp", "brute"),
+            ("brute", "num"),
+        ],
+        ids=["cpp-py", "cpp-lin", "cpp-num", "cpp-brute", "brute-num"],
+    )
+    def test_compare(self, solvers, kwargs, compare, request):
+        kwargs = dict(kwargs)
+        max_error_cpp_py = kwargs.pop("max_error_cpp_py", 1e-12)
+        max_error_cpp_lin = kwargs.pop("max_error_cpp_lin", 1e-12)
+        max_error_cpp_num = kwargs.pop("max_error_cpp_num", 2e-3)
+        max_error_cpp_brute = kwargs.pop("max_error_cpp_brute", 2e-3)
+        max_error_brute_num = kwargs.pop("max_error_brute_num", 2e-3)
+        atol = {
+            ("cpp", "py"): max_error_cpp_py,
+            ("cpp", "lin"): max_error_cpp_lin,
+            ("cpp", "num"): max_error_cpp_num,
+            ("cpp", "brute"): max_error_cpp_brute,
+            ("brute", "num"): max_error_brute_num,
+        }
+        try:
+            sT0 = solvers[compare[0]].get_sT(**kwargs)
+            sT1 = solvers[compare[1]].get_sT(**kwargs)
+            diff = np.abs(sT0 - sT1)
+            print(
+                "Med: {:.3e} | Max: {:.3e} | Tol: {:.3e}".format(
+                    np.median(diff), np.max(diff), atol[compare]
+                )
+            )
+            assert np.allclose(sT0, sT1, atol=atol[compare])
+        except Exception as e:
+            file = "{}.pdf".format(request.node.location[-1])
+            oblate.draw(file=file, **kwargs)
+            raise e
+
+
+@pytest.mark.parametrize("kwargs", [dict(bo=0.5, ro=0.1, f=0.1, theta=0.5)])
+class TestOccultorDoesntTouchLimb(Compare):
+    pass
+
+
+@pytest.mark.parametrize("kwargs", [dict(bo=0.85, ro=0.1, f=0.1, theta=0.5)])
+class TestKsqGreaterThanOne(Compare):
+    pass
+
+
+@pytest.mark.parametrize(
+    "kwargs",
+    [
+        dict(bo=0.91, ro=0.1, f=0.1, theta=0.8),
+        dict(bo=0.85, ro=0.1, f=0.1, theta=0.1),
+    ],
+)
+class TestKsqLessThanOne(Compare):
+    pass
+
+
+@pytest.mark.parametrize(
+    "kwargs",
+    [
+        dict(bo=0.95, ro=0.1, f=0.0, theta=0.8, max_error_cpp_num=1e-12),
+        dict(bo=0.95, ro=0.1, f=1e-15, theta=0.8, max_error_cpp_num=1e-12),
+        dict(bo=0.95, ro=0.1, f=1e-12, theta=0.8, max_error_cpp_num=1e-12),
+    ],
+)
+class TestNearlySpherical(Compare):
+    pass
+
+
+@pytest.mark.parametrize(
+    "kwargs",
+    [
+        dict(bo=0.0, ro=0.1, f=0.2, theta=0.8, max_error_cpp_num=5e-2),
+        dict(bo=1e-15, ro=0.1, f=0.1, theta=0.8),
+        dict(bo=1e-12, ro=0.1, f=0.1, theta=0.8),
+    ],
+)
+class TestSmallImpactParameter(Compare):
+    pass
+
+
+@pytest.mark.parametrize(
+    "kwargs",
+    [
+        dict(bo=0.5, ro=1e-6, f=0.1, theta=0.5, max_error_cpp_num=1e-12),
+        dict(bo=0.5, ro=1e-4, f=0.1, theta=0.5, max_error_cpp_num=1e-9),
+    ],
+)
+class TestSmallRadius(Compare):
+    pass
+
+
+@pytest.mark.parametrize(
+    "kwargs",
+    [
+        dict(bo=0.85, ro=0.1, f=0.1, theta=0.0),
+        dict(
+            bo=0.85,
+            ro=0.1,
+            f=0.1,
+            theta=0.5 * np.pi,
+            max_error_cpp_py=1e-6,
+            max_error_cpp_lin=1e-6,
+        ),
+    ],
+)
+class TestSinOrCosThetaIsZero(Compare):
+    pass
+
+
+@pytest.mark.parametrize(
+    "kwargs",
+    [
+        dict(bo=0.1, ro=0.1, f=0.1, theta=0.5),
+        dict(
+            bo=0.5,
+            ro=0.5,
+            f=0.1,
+            theta=0.5,
+            max_error_cpp_lin=1e-5,
+            max_error_cpp_num=5e-3,
+            max_error_cpp_brute=5e-3,
+        ),
+    ],
+)
+class TestBoEqualsRo(Compare):
+    pass
+
+
+@pytest.mark.parametrize(
+    "kwargs",
+    [dict(bo=0.9, ro=0.1, f=0.1, theta=0.5, max_error_cpp_lin=1e-6),],
+)
+class TestBoEqualsOneMinusRo(Compare):
+    pass
+
+
+@pytest.mark.parametrize(
+    "kwargs", [dict(bo=2.0, ro=0.1, f=0.1, theta=0.5),],
+)
+class TestNoOccultation(Compare):
+    pass
+
+
+@pytest.mark.parametrize(
+    "kwargs", [dict(bo=0.3, ro=2.0, f=0.1, theta=0.5),],
+)
+class TestCompleteOccultation(Compare):
+    pass
+
+
+@pytest.mark.parametrize(
+    "kwargs", [dict(bo=1.00058, ro=0.1, f=0.1, theta=-1.26629)],
+)
+class TestKsqLessThanHalf(Compare):
+    """
+    These configurations correspond to
+
+        1 - bo ** 2 - ro ** 2 < 0
+    
+    which is a special case for the J integral.
+    """
+
+    pass
+
+
+@pytest.mark.parametrize(
+    "kwargs",
+    [
+        dict(bo=1.0998151795596858, ro=0.1, f=0.1, theta=-1.4797462082413333),
+        dict(bo=1.0476234129628497, ro=0.05, f=0.1, theta=-1.28036820589747),
+        dict(
+            bo=1.004987562112089,
+            ro=0.1,
+            f=0.1,
+            theta=1.4711276743037345,
+            max_error_cpp_py=1e-6,
+            max_error_cpp_lin=1e-6,
+        ),
+        dict(bo=1.0956051937836104, ro=0.1, f=0.1, theta=1.5616688245611983),
+        dict(
+            bo=1.4251233204133293,
+            ro=0.5,
+            f=0.1,
+            theta=1.3587016514315267,
+            max_error_cpp_py=1e-6,
+            max_error_cpp_lin=1e-6,
+        ),
+        dict(
+            bo=0.9568856885688621,
+            ro=0.1,
+            f=0.25,
+            theta=-1.5707962222892058,
+            max_error_cpp_py=1e-5,
+            max_error_cpp_lin=1e-5,
+        ),
+        dict(
+            bo=1.0961596159615965,
+            ro=0.1,
+            f=0.1,
+            theta=0.5 * np.pi,
+            max_error_cpp_py=1e-3,  # TODO: Can we improve on this?
+            max_error_cpp_lin=1e-3,
+        ),
+    ],
+)
+class TestEdgeCases(Compare):
+    pass
+
+
+@pytest.mark.parametrize(
+    "kwargs", [],
+)
+class TestDebug(Compare):
+    def test_debug(self, kwargs):
+        # All the solvers
+        sTcpp = oblate.CppSolver(2).get_sT(**kwargs)
+        sTpy = oblate.PythonSolver(2).get_sT(**kwargs)
+        sTlin = oblate.NumericalSolver(2, linear=True).get_sT(**kwargs)
+        sTnum = oblate.NumericalSolver(2, linear=False).get_sT(**kwargs)
+        sTbrute = oblate.BruteSolver(2).get_sT(**kwargs)
+
+        # NaN indices
+        icpp = np.isnan(sTcpp)
+        ipy = np.isnan(sTpy)
+        inum = np.isnan(sTnum)
+        ilin = np.isnan(sTlin)
+        sTcpp[icpp] = -1.0
+        sTpy[ipy] = -1.0
+        sTnum[inum] = -1.0
+        sTlin[ilin] = -1.0
+
+        # Plot comparison
+        fig, ax = plt.subplots(3, figsize=(10, 8))
+
+        ax[0].plot(sTcpp, label="c++")
+        ax[0].plot(np.arange(len(sTcpp))[icpp], sTcpp[icpp], "ro")
+        ax[0].plot(sTpy, label="py")
+        ax[0].plot(np.arange(len(sTpy))[ipy], sTpy[ipy], "ro")
+        ax[0].legend()
+
+        ax[1].plot(sTcpp, label="c++")
+        ax[1].plot(np.arange(len(sTcpp))[icpp], sTcpp[icpp], "ro")
+        ax[1].plot(sTlin, label="lin")
+        ax[1].plot(np.arange(len(sTlin))[ilin], sTlin[ilin], "ro")
+        ax[1].legend()
+
+        ax[2].plot(sTnum, label="num")
+        ax[2].plot(np.arange(len(sTnum))[inum], sTlin[inum], "ro")
+        ax[2].plot(sTbrute, label="brute")
+        ax[2].legend()
+
+        # View geometry
+        oblate.draw(**kwargs)
+        plt.show()
+
+
+@pytest.mark.parametrize(
+    "kwargs", [dict(bo=0.95, ro=0.25, f=0.2, theta=0.5)],
+)
+class TestCppExact:
+    def test_debug(self, kwargs):
+        # All the solvers
+        lmax = 4
+        sTcpp = oblate.CppSolver(lmax).get_sT(**kwargs)
+        sTpy = oblate.PythonSolver(lmax).get_sT(**kwargs)
+        sTlin = oblate.NumericalSolver(lmax, linear=True).get_sT(**kwargs)
+        sTnum = oblate.NumericalSolver(lmax, linear=False).get_sT(**kwargs)
+        sTbrute = oblate.BruteSolver(lmax).get_sT(**kwargs)
+
+        # NaN indices
+        icpp = np.isnan(sTcpp)
+        ipy = np.isnan(sTpy)
+        inum = np.isnan(sTnum)
+        ilin = np.isnan(sTlin)
+        sTcpp[icpp] = -1.0
+        sTpy[ipy] = -1.0
+        sTnum[inum] = -1.0
+        sTlin[ilin] = -1.0
+
+        # Cases
+        cases = oblate.PythonSolver(lmax).cases
+        ibad = np.abs(sTcpp - sTnum) > 1e-3
+
+        # Plot comparison
+        fig, ax = plt.subplots(3, figsize=(10, 8))
+
+        ax[0].plot(sTcpp, label="c++")
+        ax[0].plot(np.arange(len(sTcpp))[icpp], sTcpp[icpp], "ro")
+        ax[0].plot(sTnum, label="num")
+        ax[0].plot(np.arange(len(sTnum))[inum], sTnum[inum], "ro")
+        ax[0].legend()
+        diff = np.log10(np.abs(sTcpp - sTnum))
+        diff[diff < -15] = -15
+        ax[1].plot(diff, "k-", label="diff")
+
+        ax[2].axhline(1, lw=1, alpha=0.3, color="k")
+        ax[2].axhline(2, lw=1, alpha=0.3, color="k")
+        ax[2].axhline(3, lw=1, alpha=0.3, color="k")
+        ax[2].axhline(4, lw=1, alpha=0.3, color="k")
+        ax[2].axhline(5, lw=1, alpha=0.3, color="k")
+
+        ax[2].set_yticks([1, 2, 3, 4, 5])
+
+        ax[2].plot(cases, "ro", alpha=0.25)
+        ax[2].plot(np.arange(len(cases))[ibad], cases[ibad], "ro")
+
+        # View geometry
+        oblate.draw(**kwargs)
+        plt.show()
diff --git a/starry/_core/ops/lib/include/oblate/tests/test_speed.py b/starry/_core/ops/lib/include/oblate/tests/test_speed.py
new file mode 100644
index 000000000..15fa646ab
--- /dev/null
+++ b/starry/_core/ops/lib/include/oblate/tests/test_speed.py
@@ -0,0 +1,53 @@
+import oblate
+import numpy as np
+import time
+import os
+import pytest
+
+
+@pytest.mark.skipif(
+    os.getenv("CI", "false") == "false", reason="Only run this on CI."
+)
+def test_speed(solvers, nruns=1000, niter=30):
+    """
+    Compare the speed of the C++ and Python implementations.
+
+    """
+    # Run the timing tests a bunch of times
+    cpptime = np.zeros(niter)
+    cppxtime = np.zeros(niter)
+    pytime = np.zeros(niter)
+    for k in range(niter):
+
+        # Time the C++ solver (linear)
+        tstart = time.time()
+        res = solvers["cpp"].get_sT(nruns=nruns)
+        cpptime[k] = (time.time() - tstart) / nruns
+
+        # Time the C++ solver (exact)
+        tstart = time.time()
+        res = solvers["cppx"].get_sT(nruns=nruns)
+        cppxtime[k] = (time.time() - tstart) / nruns
+
+        # Time the Python solver
+        tstart = time.time()
+        res = solvers["py"].get_sT(nruns=nruns)
+        pytime[k] = (time.time() - tstart) / nruns
+
+    # Log
+    print("")
+    print(
+        "C++ (linear):    {:.3} +/- {:.3} ms".format(
+            np.median(cppxtime) * 1e3, np.std(cppxtime) * 1e3
+        )
+    )
+    print(
+        "C++ (exact):     {:.3} +/- {:.3} ms".format(
+            np.median(cpptime) * 1e3, np.std(cpptime) * 1e3
+        )
+    )
+    print(
+        "Python (linear): {:.3} +/- {:.3} ms".format(
+            np.median(pytime) * 1e3, np.std(pytime) * 1e3
+        )
+    )
diff --git a/starry/extensions/tests/greedy/conftest.py b/starry/extensions/tests/greedy/conftest.py
index 473d611f0..09136cc7b 100644
--- a/starry/extensions/tests/greedy/conftest.py
+++ b/starry/extensions/tests/greedy/conftest.py
@@ -1,4 +1,9 @@
 # -*- coding: utf-8 -*-
+import pytest
 import starry
 
-starry.config.lazy = False
+
+@pytest.fixture(scope="module", autouse=True)
+def setup():
+    starry.config.lazy = False
+    yield
diff --git a/starry/extensions/tests/lazy/conftest.py b/starry/extensions/tests/lazy/conftest.py
index f513c9f55..3967c244b 100644
--- a/starry/extensions/tests/lazy/conftest.py
+++ b/starry/extensions/tests/lazy/conftest.py
@@ -1,4 +1,9 @@
 # -*- coding: utf-8 -*-
+import pytest
 import starry
 
-starry.config.lazy = True
+
+@pytest.fixture(scope="module", autouse=True)
+def setup():
+    starry.config.lazy = True
+    yield
diff --git a/tests/greedy/conftest.py b/tests/greedy/conftest.py
index 505d8dee9..3c70c7e21 100644
--- a/tests/greedy/conftest.py
+++ b/tests/greedy/conftest.py
@@ -1,6 +1,11 @@
 # -*- coding: utf-8 -*-
+import pytest
 import starry
 import matplotlib
 
-matplotlib.use("Agg")
-starry.config.lazy = False
+
+@pytest.fixture(scope="module", autouse=True)
+def setup():
+    matplotlib.use("Agg")
+    starry.config.lazy = False
+    yield
diff --git a/tests/greedy/test_compile.py b/tests/greedy/test_compile.py
index 507b2c8db..6bfc071c0 100644
--- a/tests/greedy/test_compile.py
+++ b/tests/greedy/test_compile.py
@@ -2,7 +2,7 @@
 """Test theano compiling."""
 import starry
 import pytest
-from starry.compat import theano
+from starry.compat import change_flags
 import numpy as np
 
 
@@ -14,7 +14,6 @@ def test_test_value_error():
     test values are missing.
 
     """
-    theano.config.compute_test_value = "raise"
-    map = starry.Map()
-    flux = map.flux()
-    assert theano.config.compute_test_value == "raise"
+    with change_flags(compute_test_value="raise"):
+        map = starry.Map()
+        flux = map.flux()
diff --git a/tests/greedy/test_cpp.py b/tests/greedy/test_cpp.py
new file mode 100644
index 000000000..6a72dc678
--- /dev/null
+++ b/tests/greedy/test_cpp.py
@@ -0,0 +1,14 @@
+# -*- coding: utf-8 -*-
+"""
+Tests for C++ functions. TODO.
+
+"""
+from starry import _c_ops as Ops
+import numpy as np
+import pytest
+
+
+# Only run tests if `STARRY_UNIT_TESTS=1` on compile
+cpp = pytest.mark.skipif(
+    not Ops.STARRY_UNIT_TESTS, reason="c++ unit tests not found"
+)
diff --git a/tests/greedy/test_intensity.py b/tests/greedy/test_intensity.py
index c6fac323f..b16bb14b1 100644
--- a/tests/greedy/test_intensity.py
+++ b/tests/greedy/test_intensity.py
@@ -5,7 +5,11 @@
 import pytest
 
 # Setup
-map = starry.Map(ydeg=1)
+@pytest.fixture(autouse=True)
+def map():
+    return starry.Map(ydeg=1)
+
+
 params = [
     [[0, 0, 0], 0, 0, 1.0 / np.pi],
     [[1, 0, 0], 0, np.linspace(-180, 180, 10), 1.0 / np.pi],
@@ -21,12 +25,12 @@
 
 
 @pytest.mark.parametrize("y,lat,lon,I", params)
-def test_intensity(y, lat, lon, I):
+def test_intensity(y, lat, lon, I, map):
     map[1:, :] = y
     assert np.allclose(map.intensity(lat=lat, lon=lon), I)
 
 
-def test_design_matrix():
+def test_design_matrix(map):
     map[1:, :] = [0.1, 0.2, 0.3]
     lat = [0, 30]
     lon = [45, 90]
@@ -34,7 +38,7 @@ def test_design_matrix():
     assert np.allclose(map.intensity(lat=lat, lon=lon), P.dot(map.y))
 
 
-def test_limb_darkened():
+def test_limb_darkened(map):
     map = starry.Map(udeg=2)
     map[1] = 0.5
     map[2] = 0.25
diff --git a/tests/greedy/test_misc.py b/tests/greedy/test_misc.py
index 5b37e332a..c2e545fb7 100644
--- a/tests/greedy/test_misc.py
+++ b/tests/greedy/test_misc.py
@@ -7,6 +7,7 @@
 import io
 import logging
 import warnings
+import pytest
 
 
 def test_quiet():
@@ -35,6 +36,7 @@ def test_quiet():
         logger.removeHandler(ch)
 
 
+@pytest.mark.skip(reason="Kwargs checks not yet implemented.")
 def test_check_kwargs_map(caplog):
     """Test that we capture bad keyword arguments."""
     starry.config.quiet = False
@@ -49,6 +51,7 @@ def test_check_kwargs_map(caplog):
     )
 
 
+@pytest.mark.skip(reason="Kwargs checks not yet implemented.")
 def test_check_kwargs_map_flux(caplog):
     """Test that we capture bad keyword arguments."""
     starry.config.quiet = False
diff --git a/tests/greedy/test_numerical.py b/tests/greedy/test_numerical.py
index f45126b68..56940e02e 100644
--- a/tests/greedy/test_numerical.py
+++ b/tests/greedy/test_numerical.py
@@ -8,11 +8,10 @@
 import numpy as np
 
 res = 1000
-map_tmp = starry.Map(1)
 
 
 def flux(map, xs=0, ys=0, zs=0, xo=0, yo=0, zo=1, ro=0, **kwargs):
-    (lat, lon), (x, y, z) = map_tmp.ops.compute_ortho_grid(res)
+    (lat, lon), (x, y, z) = starry.Map(1).ops.compute_ortho_grid(res)
     occulted = (xo - x) ** 2 + (yo - y) ** 2 < ro ** 2
     occulted = occulted.reshape(res, res)
     if map._ops_class_ == starry._core.OpsReflected:
diff --git a/tests/greedy/test_oblate.py b/tests/greedy/test_oblate.py
new file mode 100644
index 000000000..dfc8693ad
--- /dev/null
+++ b/tests/greedy/test_oblate.py
@@ -0,0 +1,283 @@
+import starry
+import numpy as np
+import pytest
+from scipy.integrate import dblquad
+
+
+class StarrySolver:
+    """
+    Wrapper for the `sT` solution vector computed by `starry`.
+
+    """
+
+    def __init__(self, lmax=5):
+        self.map = starry.Map(lmax, oblate=True, gdeg=0)
+
+    def get_sT(self, bo=0.58, ro=0.4, f=0.2, theta=0.5, **kwargs):
+        return self.map.ops.sT(
+            np.array(f, dtype="float64"),
+            np.array([theta], dtype="float64"),
+            np.array([bo], dtype="float64"),
+            np.array(ro, dtype="float64"),
+        ).flatten()
+
+
+class BruteSolver:
+    """
+    Brute-force grid integration for computing the `sT` solution vector.
+
+    """
+
+    def __init__(self, lmax=5):
+        self.lmax = lmax
+
+        # Integrate the Green's basis over
+        # the entire unit disk
+        self.sT0 = np.zeros((self.lmax + 1) ** 2)
+        for n in range((self.lmax + 1) ** 2):
+            self.sT0[n], _ = dblquad(
+                lambda x, y: self.g(n, x, y),
+                -1,
+                1,
+                lambda x: -np.sqrt(1 - x ** 2),
+                lambda x: np.sqrt(1 - x ** 2),
+            )
+
+    def g(self, n, x, y, z=None):
+        """
+        Return the nth term of the Green's basis.
+
+        """
+        if z is None:
+            z2 = 1 - x ** 2 - y ** 2
+            z = np.sqrt(np.abs(z2))
+            on_star = z2 > 0
+            if hasattr(z, "__len__"):
+                z[~on_star] = np.nan
+            else:
+                if not on_star:
+                    z2 = np.nan
+        l = int(np.floor(np.sqrt(n)))
+        m = n - l * l - l
+        mu = l - m
+        nu = l + m
+        if nu % 2 == 0:
+            I = [mu // 2]
+            J = [nu // 2]
+            K = [0]
+            C = [(mu + 2) // 2]
+        elif (l == 1) and (m == 0):
+            I = [0]
+            J = [0]
+            K = [1]
+            C = [1]
+        elif (mu == 1) and (l % 2 == 0):
+            I = [l - 2]
+            J = [1]
+            K = [1]
+            C = [3]
+        elif mu == 1:
+            I = [l - 3, l - 1, l - 3]
+            J = [0, 0, 2]
+            K = [1, 1, 1]
+            C = [-1, 1, 4]
+        else:
+            I = [(mu - 5) // 2, (mu - 5) // 2, (mu - 1) // 2]
+            J = [(nu - 1) // 2, (nu + 3) // 2, (nu - 1) // 2]
+            K = [1, 1, 1]
+            C = [(mu - 3) // 2, -(mu - 3) // 2, -(mu + 3) // 2]
+        res = z * 0
+        for i, j, k, c in zip(I, J, K, C):
+            if c != 0:
+                res += c * x ** i * y ** j * z ** k
+        return res
+
+    def get_sT(self, bo=0.58, ro=0.4, f=0.2, theta=0.5, res=999, **kwargs):
+
+        # Grid up the occultor
+        xo = bo * np.sin(theta)
+        yo = bo * np.cos(theta)
+        x = np.linspace(xo - ro, xo + ro, res)
+        y = np.linspace(yo - ro, yo + ro, res)
+        x, y = np.meshgrid(x, y)
+        b = 1 - f
+        on_star = 1 - x ** 2 - (y / b) ** 2 > 0
+
+        # Compute the Green's basis in (x, y') where
+        #
+        #     y' = y / b
+        #
+        # is the transformed y coordinate
+        g = np.array(
+            [self.g(n, x, y / b) for n in range((self.lmax + 1) ** 2)]
+        )
+
+        # Compute the masked pixels
+        under_occultor = (x - xo) ** 2 + (y - yo) ** 2 <= ro ** 2
+        inds = on_star & under_occultor
+        g /= np.count_nonzero(under_occultor)
+
+        # The solution vector for the complement of the visible region
+        sTbar = np.pi * ro ** 2 * np.sum(g[:, inds], axis=1)
+
+        # The solution vector over the entire star
+        sT0 = b * self.sT0
+
+        return sT0 - sTbar
+
+
+@pytest.fixture(scope="module")
+def solvers(lmax=5):
+    return StarrySolver(lmax), BruteSolver(lmax)
+
+
+class Compare:
+    """
+    Base class for comparing the `starry` solver to the brute solver.
+
+    """
+
+    def test_compare(self, solvers, kwargs):
+        atol = kwargs.get("atol", 2.0e-6)
+        sT_starry = solvers[0].get_sT(**kwargs)
+        sT_brute = solvers[1].get_sT(**kwargs)
+        diff = np.abs(sT_starry - sT_brute)
+        print(
+            "Med: {:.3e} | Max: {:.3e} | Tol: {:.3e}".format(
+                np.median(diff), np.max(diff), atol
+            )
+        )
+        if kwargs.get("complete", False):
+            assert np.allclose(sT_starry, 0.0, atol)
+        else:
+            assert np.allclose(sT_starry, sT_brute, atol=atol)
+
+
+@pytest.mark.parametrize("kwargs", [dict(bo=0.5, ro=0.1, f=0.1, theta=0.5)])
+class TestOccultorDoesntTouchLimb(Compare):
+    pass
+
+
+@pytest.mark.parametrize("kwargs", [dict(bo=0.85, ro=0.1, f=0.1, theta=0.5)])
+class TestKsqGreaterThanOne(Compare):
+    pass
+
+
+@pytest.mark.parametrize(
+    "kwargs",
+    [
+        dict(bo=0.91, ro=0.1, f=0.1, theta=0.8),
+        dict(bo=0.85, ro=0.1, f=0.1, theta=0.1),
+    ],
+)
+class TestKsqLessThanOne(Compare):
+    pass
+
+
+@pytest.mark.parametrize(
+    "kwargs",
+    [
+        dict(bo=0.95, ro=0.1, f=0.0, theta=0.8),
+        dict(bo=0.95, ro=0.1, f=1e-15, theta=0.8),
+        dict(bo=0.95, ro=0.1, f=1e-12, theta=0.8),
+    ],
+)
+class TestNearlySpherical(Compare):
+    pass
+
+
+@pytest.mark.parametrize(
+    "kwargs",
+    [
+        dict(bo=0.0, ro=0.1, f=0.2, theta=0.8),
+        dict(bo=1e-15, ro=0.1, f=0.1, theta=0.8),
+        dict(bo=1e-12, ro=0.1, f=0.1, theta=0.8),
+    ],
+)
+class TestSmallImpactParameter(Compare):
+    pass
+
+
+@pytest.mark.parametrize(
+    "kwargs",
+    [
+        dict(bo=0.5, ro=1e-6, f=0.1, theta=0.5),
+        dict(bo=0.5, ro=1e-4, f=0.1, theta=0.5),
+    ],
+)
+class TestSmallRadius(Compare):
+    pass
+
+
+@pytest.mark.parametrize(
+    "kwargs",
+    [
+        dict(bo=0.85, ro=0.1, f=0.1, theta=0.0),
+        dict(bo=0.85, ro=0.1, f=0.1, theta=0.5 * np.pi),
+    ],
+)
+class TestSinOrCosThetaIsZero(Compare):
+    pass
+
+
+@pytest.mark.parametrize(
+    "kwargs",
+    [
+        dict(bo=0.1, ro=0.1, f=0.1, theta=0.5),
+        dict(bo=0.5, ro=0.5, f=0.1, theta=0.5, atol=1e-4),
+    ],
+)
+class TestBoEqualsRo(Compare):
+    pass
+
+
+@pytest.mark.parametrize("kwargs", [dict(bo=0.9, ro=0.1, f=0.1, theta=0.5)])
+class TestBoEqualsOneMinusRo(Compare):
+    pass
+
+
+@pytest.mark.parametrize("kwargs", [dict(bo=2.0, ro=0.1, f=0.1, theta=0.5)])
+class TestNoOccultation(Compare):
+    pass
+
+
+@pytest.mark.parametrize(
+    "kwargs", [dict(bo=0.3, ro=2.0, f=0.1, theta=0.5, complete=True)]
+)
+class TestCompleteOccultation(Compare):
+    pass
+
+
+@pytest.mark.parametrize(
+    "kwargs", [dict(bo=1.00058, ro=0.1, f=0.1, theta=-1.26629)]
+)
+class TestKsqLessThanHalf(Compare):
+    pass
+
+
+@pytest.mark.parametrize(
+    "kwargs",
+    [
+        dict(bo=1.0998151795596858, ro=0.1, f=0.1, theta=-1.4797462082413333),
+        dict(bo=1.0476234129628497, ro=0.05, f=0.1, theta=-1.28036820589747),
+        dict(
+            bo=1.004987562112089,
+            ro=0.1,
+            f=0.1,
+            theta=1.4711276743037345,
+            atol=1e-5,
+        ),
+        dict(bo=1.0956051937836104, ro=0.1, f=0.1, theta=1.5616688245611983),
+        dict(
+            bo=1.4251233204133293,
+            ro=0.5,
+            f=0.1,
+            theta=1.3587016514315267,
+            atol=1e-4,
+        ),
+        dict(bo=0.9568856885688621, ro=0.1, f=0.25, theta=-1.5707962222892058),
+        dict(bo=1.0961596159615965, ro=0.1, f=0.1, theta=0.5 * np.pi),
+    ],
+)
+class TestEdgeCases(Compare):
+    pass
diff --git a/tests/greedy/test_oblate_system.py b/tests/greedy/test_oblate_system.py
new file mode 100644
index 000000000..21373676a
--- /dev/null
+++ b/tests/greedy/test_oblate_system.py
@@ -0,0 +1,34 @@
+import starry
+import numpy as np
+import matplotlib.pyplot as plt
+import pytest
+
+
+@pytest.mark.parametrize("ydeg,nw", [[0, None], [0, 10], [1, None], [1, 10]])
+def test_system(ydeg, nw):
+
+    # Oblate map
+    map = starry.Map(udeg=2, ydeg=ydeg, oblate=True, nw=nw)
+    map[1] = 0.5
+    map[2] = 0.25
+    map.omega = 0.5
+    map.beta = 1.23
+    map.tpole = 8000
+    map.f = 1 - 2 / (map.omega ** 2 + 2)
+    map.obl = 30
+
+    # Compute system flux
+    star = starry.Primary(map, r=1.5)
+    planet = starry.Secondary(starry.Map(amp=0, nw=nw), porb=1.0, r=0.1, m=0)
+    sys = starry.System(star, planet)
+    t = np.linspace(-0.1, 0.1, 1000)
+    flux_sys = sys.flux(t, integrated=True)
+
+    # Compute map flux manually
+    x, y, z = sys.position(t)
+    xo = x[1] / star._r
+    yo = y[1] / star._r
+    flux_map = map.flux(xo=xo, yo=yo, ro=planet._r / star._r, integrated=True)
+
+    # Check that they agree
+    assert np.allclose(flux_map, flux_sys)
diff --git a/tests/greedy/test_reflected_normalization.py b/tests/greedy/test_reflected_normalization.py
index f22fb6c3a..797a9d6d3 100644
--- a/tests/greedy/test_reflected_normalization.py
+++ b/tests/greedy/test_reflected_normalization.py
@@ -9,33 +9,35 @@
 
 import numpy as np
 import starry
+import pytest
 
 
-starry.config.lazy = False
-map = starry.Map(1, reflected=True)
+@pytest.fixture(autouse=True)
+def map():
+    return starry.Map(1, reflected=True)
 
 
-def test_amp():
+def test_amp(map):
     assert map.amp == 1.0
 
 
-def test_flux():
+def test_flux(map):
     assert np.allclose(map.flux(zs=1, xs=0, ys=0), 2.0 / 3.0)
 
 
-def test_design_matrix():
+def test_design_matrix(map):
     assert np.allclose(
         map.amp * map.design_matrix(zs=1, xs=0, ys=0)[0, 0], 2.0 / 3.0
     )
 
 
-def test_intensity():
+def test_intensity(map):
     assert np.allclose(
         map.intensity(lat=0, lon=0, zs=1, xs=0, ys=0), 1.0 / np.pi
     )
 
 
-def test_albedo():
+def test_albedo(map):
     assert np.allclose(
         map.amp
         * map.intensity(lat=0, lon=0, zs=1, xs=0, ys=0, illuminate=False),
@@ -43,23 +45,23 @@ def test_albedo():
     )
 
 
-def test_albedo_design_matrix():
+def test_albedo_design_matrix(map):
     assert np.allclose(
         map.amp * map.intensity_design_matrix(lat=0, lon=0)[0, 0], 1.0
     )
 
 
-def test_render():
+def test_render(map):
     assert np.allclose(np.nanmax(map.render(zs=1, xs=0, ys=0)), 1.0 / np.pi)
 
 
-def test_render_albedo():
+def test_render_albedo(map):
     assert np.allclose(
         np.nanmax(map.render(zs=1, xs=0, ys=0, illuminate=False)), 1.0
     )
 
 
-def test_numerical_flux(res=300):
+def test_numerical_flux(map, res=300):
     res = 300
     atol = 1e-4
     rtol = 1e-4
@@ -67,7 +69,7 @@ def test_numerical_flux(res=300):
     assert np.allclose(num_flux, 2.0 / 3.0, atol=atol, rtol=rtol)
 
 
-def test_solve():
+def test_solve(map):
     # Generate a light curve with two points: one at
     # full phase, one at new phase.
     np.random.seed(0)
@@ -92,7 +94,7 @@ def test_solve():
     assert np.allclose(map.amp, 1.0)
 
 
-def test_one_over_r_squared(n_tests=10, plot=False):
+def test_one_over_r_squared(map, n_tests=10, plot=False):
     """Test that the flux decreases as 1/r^2."""
     flux0 = map.flux()
     zs = np.linspace(1, 10, 100)
diff --git a/tests/greedy/test_solve_greedy.py b/tests/greedy/test_solve_greedy.py
index 429491838..9fda97562 100644
--- a/tests/greedy/test_solve_greedy.py
+++ b/tests/greedy/test_solve_greedy.py
@@ -10,36 +10,46 @@
 import pytest
 import itertools
 
+
+@pytest.fixture(autouse=True)
+def data():
+
+    # Instantiate a dipole map
+    map = starry.Map(ydeg=1, reflected=True)
+    amp_true = 0.75
+    inc_true = 60
+    y_true = np.array([1, 0.1, 0.2, 0.3])
+    map.amp = amp_true
+    map[1, :] = y_true[1:]
+    map.inc = inc_true
+
+    # Generate a synthetic light curve with just a little noise
+    theta = np.linspace(0, 360, 100)
+    phi = 3.5 * theta
+    xs = np.cos(phi * np.pi / 180)
+    ys = 0.1 * np.cos(phi * np.pi / 180)
+    zs = np.sin(phi * np.pi / 180)
+    kwargs = dict(theta=theta, xs=xs, ys=ys, zs=zs)
+    flux = map.flux(**kwargs)
+    sigma = 1e-5
+    np.random.seed(1)
+    flux += np.random.randn(len(theta)) * sigma
+
+    return (map, kwargs, amp_true, inc_true, y_true, sigma, flux)
+
+
 # Parameter combinations we'll test
 vals = ["scalar", "vector", "matrix", "cholesky"]
 woodbury = [False, True]
 solve_inputs = itertools.product(vals, vals)
 lnlike_inputs = itertools.product(vals, vals, woodbury)
 
-# Instantiate a dipole map
-map = starry.Map(ydeg=1, reflected=True)
-amp_true = 0.75
-inc_true = 60
-y_true = np.array([1, 0.1, 0.2, 0.3])
-map.amp = amp_true
-map[1, :] = y_true[1:]
-map.inc = inc_true
-
-# Generate a synthetic light curve with just a little noise
-theta = np.linspace(0, 360, 100)
-phi = 3.5 * theta
-xs = np.cos(phi * np.pi / 180)
-ys = 0.1 * np.cos(phi * np.pi / 180)
-zs = np.sin(phi * np.pi / 180)
-kwargs = dict(theta=theta, xs=xs, ys=ys, zs=zs)
-flux = map.flux(**kwargs)
-sigma = 1e-5
-np.random.seed(1)
-flux += np.random.randn(len(theta)) * sigma
-
 
 @pytest.mark.parametrize("L,C", solve_inputs)
-def test_solve(L, C):
+def test_solve(L, C, data):
+
+    map, kwargs, amp_true, inc_true, y_true, sigma, flux = data
+
     # Place a generous prior on the map coefficients
     if L == "scalar":
         map.set_prior(L=1)
@@ -76,8 +86,11 @@ def test_solve(L, C):
 
 
 @pytest.mark.parametrize("L,C,woodbury", lnlike_inputs)
-def test_lnlike(L, C, woodbury):
+def test_lnlike(L, C, woodbury, data):
     """Test the log marginal likelihood method."""
+
+    map, kwargs, amp_true, inc_true, y_true, sigma, flux = data
+
     # Place a generous prior on the map coefficients
     if L == "scalar":
         map.set_prior(L=1)
diff --git a/tests/greedy/test_system_rv_greedy.py b/tests/greedy/test_system_rv_greedy.py
index bc32f5280..cb29b41fd 100644
--- a/tests/greedy/test_system_rv_greedy.py
+++ b/tests/greedy/test_system_rv_greedy.py
@@ -12,17 +12,23 @@ def test_compare_to_map_rv():
     """Ensure we get the same result by calling `map.rv()` and `sys.rv()`.
     """
     # Define the map
-    map = starry.Map(ydeg=1, udeg=2, rv=True, amp=1, veq=1, alpha=0)
+    map = starry.Map(ydeg=1, udeg=2, rv=True, amp=1, veq=1, alpha=0, lazy=True)
     map[1, 0] = 0.5
 
     # Define the star
     A = starry.Primary(
-        map, r=1.0, m=1.0, prot=0, length_unit=u.Rsun, mass_unit=u.Msun
+        map,
+        r=1.0,
+        m=1.0,
+        prot=0,
+        length_unit=u.Rsun,
+        mass_unit=u.Msun,
+        lazy=True,
     )
 
     # Define the planet
     b = starry.Secondary(
-        starry.Map(rv=True, amp=1, veq=0),
+        starry.Map(rv=True, amp=1, veq=0, lazy=True),
         r=0.1,
         porb=1.0,
         m=0.01,
@@ -52,8 +58,8 @@ def test_compare_to_map_rv():
     ro = b.r / A.r
 
     # Compare
-    rv1 = map.rv(xo=xo, yo=yo, ro=ro)
-    rv2 = sys.rv(time, keplerian=False)
+    rv1 = map.rv(xo=xo, yo=yo, ro=ro).eval()
+    rv2 = sys.rv(time, keplerian=False).eval()
     assert np.allclose(rv1, rv2)
 
 
diff --git a/tests/lazy/conftest.py b/tests/lazy/conftest.py
index 194ef96ae..b6b054d3e 100644
--- a/tests/lazy/conftest.py
+++ b/tests/lazy/conftest.py
@@ -1,6 +1,11 @@
 # -*- coding: utf-8 -*-
+import pytest
 import starry
 import matplotlib
 
-matplotlib.use("Agg")
-starry.config.lazy = True
+
+@pytest.fixture(scope="module", autouse=True)
+def setup():
+    matplotlib.use("Agg")
+    starry.config.lazy = True
+    yield
diff --git a/tests/lazy/test_grid.py b/tests/lazy/test_grid.py
index 9f885cc39..77476920e 100644
--- a/tests/lazy/test_grid.py
+++ b/tests/lazy/test_grid.py
@@ -13,11 +13,9 @@
 import numpy as np
 
 
-map = starry.Map(1)
-
-
 @pytest.mark.parametrize("compile", [True, False])
 def test_ortho_grid(compile):
+    map = starry.Map(1)
     for res in np.arange(30, 101):
         if compile:
             (lat, lon), (x, y, z) = map.ops.compute_ortho_grid(res)
diff --git a/tests/lazy/test_op_gradients.py b/tests/lazy/test_op_gradients.py
index dd4619891..2455e1c23 100644
--- a/tests/lazy/test_op_gradients.py
+++ b/tests/lazy/test_op_gradients.py
@@ -478,3 +478,19 @@ def test_sT_reflected(abs_tol=1e-5, rel_tol=1e-5, eps=1e-7):
             n_tests=1,
             rng=np.random,
         )
+def test_sT_oblate(abs_tol=1e-5, rel_tol=1e-5, eps=1e-7):
+    with change_flags(compute_test_value="off"):
+        map = starry.Map(gdeg=4, oblate=True)
+        f = 0.2
+        theta = np.array([0.5, 0.5])
+        bo = np.array([0.75,0.75])
+        ro = 0.5
+        theano.gradient.verify_grad(
+            map.ops.sT,
+            (f, theta, bo, ro),
+            abs_tol=abs_tol,
+            rel_tol=rel_tol,
+            eps=eps,
+            n_tests=1,
+            rng=np.random,
+        )
diff --git a/tests/lazy/test_show_lazy.py b/tests/lazy/test_show_lazy.py
index ecd4f17dc..e3c621e56 100644
--- a/tests/lazy/test_show_lazy.py
+++ b/tests/lazy/test_show_lazy.py
@@ -12,10 +12,12 @@
 import pytest
 
 
-if starry.compat.USE_AESARA:
-    theano_config = dict(aesara_config=dict(compute_test_value="ignore"))
-else:
-    theano_config = dict(theano_config=dict(compute_test_value="ignore"))
+@pytest.fixture(autouse=True)
+def theano_config():
+    if starry.compat.USE_AESARA:
+        return dict(aesara_config=dict(compute_test_value="ignore"))
+    else:
+        return dict(theano_config=dict(compute_test_value="ignore"))
 
 
 # TODO: Several tests segfault on CI. Investigate.
@@ -88,7 +90,7 @@ def test_system_rv_show():
         os.remove("tmp.mp4")
 
 
-def test_show_pymc3():
+def test_show_pymc3(theano_config):
     with pm.Model(**theano_config) as model:
         map = starry.Map()
         ncoeff = map.Ny
@@ -99,7 +101,7 @@ def test_show_pymc3():
         os.remove("tmp.pdf")
 
 
-def test_show_reflected_pymc3():
+def test_show_reflected_pymc3(theano_config):
     with pm.Model(**theano_config) as model:
         map = starry.Map(ydeg=1, udeg=1, reflected=True)
         ncoeff = map.Ny
@@ -111,7 +113,7 @@ def test_show_reflected_pymc3():
 
 
 @pytest.mark.skipif(STARRY_ON_CI, reason="Segfaults on CI")
-def test_show_rv_pymc3():
+def test_show_rv_pymc3(theano_config):
     with pm.Model(**theano_config) as model:
         map = starry.Map(ydeg=1, udeg=1, rv=True)
         ncoeff = map.Ny
@@ -122,7 +124,7 @@ def test_show_rv_pymc3():
         os.remove("tmp.pdf")
 
 
-def test_show_ld_pymc3():
+def test_show_ld_pymc3(theano_config):
     with pm.Model(**theano_config) as model:
         map = starry.Map(udeg=2)
         map[1:] = pm.MvNormal("u", [0.5, 0.25], np.eye(2), shape=(2,))
@@ -130,7 +132,7 @@ def test_show_ld_pymc3():
         os.remove("tmp.pdf")
 
 
-def test_system_show_pymc3():
+def test_system_show_pymc3(theano_config):
     with pm.Model(**theano_config) as model:
         pri = starry.Primary(starry.Map())
         sec = starry.Secondary(starry.Map(), porb=1.0)
@@ -143,7 +145,7 @@ def test_system_show_pymc3():
 
 
 @pytest.mark.skipif(STARRY_ON_CI, reason="Segfaults on CI")
-def test_system_rv_show_pymc3():
+def test_system_rv_show_pymc3(theano_config):
     with pm.Model(**theano_config) as model:
         pri = starry.Primary(starry.Map(rv=True))
         sec = starry.Secondary(starry.Map(rv=True), porb=1.0)
diff --git a/tests/lazy/test_solve_lazy.py b/tests/lazy/test_solve_lazy.py
index 7b1ad46ac..82a268dc3 100644
--- a/tests/lazy/test_solve_lazy.py
+++ b/tests/lazy/test_solve_lazy.py
@@ -10,36 +10,46 @@
 import pytest
 import itertools
 
+
+@pytest.fixture(autouse=True)
+def data():
+
+    # Instantiate a dipole map
+    map = starry.Map(ydeg=1, reflected=True)
+    amp_true = 0.75
+    inc_true = 60
+    y_true = np.array([1, 0.1, 0.2, 0.3])
+    map.amp = amp_true
+    map[1, :] = y_true[1:]
+    map.inc = inc_true
+
+    # Generate a synthetic light curve with just a little noise
+    theta = np.linspace(0, 360, 100)
+    phi = 3.5 * theta
+    xs = np.cos(phi * np.pi / 180)
+    ys = 0.1 * np.cos(phi * np.pi / 180)
+    zs = np.sin(phi * np.pi / 180)
+    kwargs = dict(theta=theta, xs=xs, ys=ys, zs=zs)
+    flux = map.flux(**kwargs).eval()
+    sigma = 1e-5
+    np.random.seed(1)
+    flux += np.random.randn(len(theta)) * sigma
+
+    return (map, kwargs, amp_true, inc_true, y_true, sigma, flux)
+
+
 # Parameter combinations we'll test
 vals = ["scalar", "vector", "matrix", "cholesky"]
 woodbury = [False, True]
 solve_inputs = itertools.product(vals, vals)
 lnlike_inputs = itertools.product(vals, vals, woodbury)
 
-# Instantiate a dipole map
-map = starry.Map(ydeg=1, reflected=True)
-amp_true = 0.75
-inc_true = 60
-y_true = np.array([1, 0.1, 0.2, 0.3])
-map.amp = amp_true
-map[1, :] = y_true[1:]
-map.inc = inc_true
-
-# Generate a synthetic light curve with just a little noise
-theta = np.linspace(0, 360, 100)
-phi = 3.5 * theta
-xs = np.cos(phi * np.pi / 180)
-ys = 0.1 * np.cos(phi * np.pi / 180)
-zs = np.sin(phi * np.pi / 180)
-kwargs = dict(theta=theta, xs=xs, ys=ys, zs=zs)
-flux = map.flux(**kwargs).eval()
-sigma = 1e-5
-np.random.seed(1)
-flux += np.random.randn(len(theta)) * sigma
-
 
 @pytest.mark.parametrize("L,C", solve_inputs)
-def test_solve(L, C):
+def test_solve(L, C, data):
+
+    map, kwargs, amp_true, inc_true, y_true, sigma, flux = data
+
     # Place a generous prior on the map coefficients
     if L == "scalar":
         map.set_prior(L=1)
@@ -78,8 +88,11 @@ def test_solve(L, C):
 
 
 @pytest.mark.parametrize("L,C,woodbury", lnlike_inputs)
-def test_lnlike(L, C, woodbury):
+def test_lnlike(L, C, woodbury, data):
     """Test the log marginal likelihood method."""
+
+    map, kwargs, amp_true, inc_true, y_true, sigma, flux = data
+
     # Place a generous prior on the map coefficients
     if L == "scalar":
         map.set_prior(L=1)