diff --git a/.gitignore b/.gitignore
index 7060ca2..12a6e5d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -5,4 +5,5 @@ Manifest.toml
 /docs/build/
 
 dev
-test/output
\ No newline at end of file
+test/output
+.CondaPkg
diff --git a/CondaPkg.toml b/CondaPkg.toml
new file mode 100644
index 0000000..3ded1f4
--- /dev/null
+++ b/CondaPkg.toml
@@ -0,0 +1,4 @@
+[deps]
+python = "3.9"
+batman-package = ""
+numpy = ""
diff --git a/Project.toml b/Project.toml
index ea209f1..66a8b55 100644
--- a/Project.toml
+++ b/Project.toml
@@ -38,3 +38,19 @@ StatsFuns = "0.9"
 Unitful = "1.9"
 UnitfulAstro = "1"
 julia = "1.5"
+
+[extras]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"
+CondaPkg = "992eb4ea-22a4-4c89-a5bb-47a3300528ab"
+FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["BenchmarkTools", "ChainRulesTestUtils", "CondaPkg", "FiniteDifferences", "ForwardDiff", "HypothesisTests", "Pkg", "PythonCall", "QuadGK", "StableRNGs", "Test"]
diff --git a/test/Project.toml b/test/Project.toml
index aeea72d..1dad88b 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -4,12 +4,13 @@ Bijectors = "76274a88-744f-5084-9051-94815aaf08c4"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"
 Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d"
+CondaPkg = "992eb4ea-22a4-4c89-a5bb-47a3300528ab"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
-PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
+PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/orbits/keplerian.jl b/test/orbits/keplerian.jl
index 56df961..b8f96d0 100644
--- a/test/orbits/keplerian.jl
+++ b/test/orbits/keplerian.jl
@@ -6,11 +6,7 @@ using Transits.Orbits: KeplerianOrbit, flip,
                        stringify_units
 
 # Setup python env
-import Pkg
-ENV["PYTHON"] = ""
-Pkg.build("PyCall")
-using PyCall, Conda
-Conda.add(["batman-package"]; channel="conda-forge")
+using PythonCall, CondaPkg
 
 # Constants
 const G_nom = 2942.2062175044193 # Rsun^3/Msun/d^2
@@ -32,49 +28,46 @@ as_matrix(pos) = reinterpret(reshape, Float64, pos) |> permutedims
 
 @testset "KeplerianOrbit: sky coords" begin
     # Comparison coords from `batman`
-    py"""
+    @pyexec """
     import numpy as np
     from batman import _rsky
 
-    def sky_coords():
-        t = np.linspace(-100, 100, 1_000)
-
-        t0, period, a, e, omega, incl = (
-            x.flatten()
-            for x in np.meshgrid(
-                np.linspace(-5.0, 5.0, 2),
-                np.exp(np.linspace(np.log(5.0), np.log(50.0), 3)),
-                np.linspace(50.0, 100.0, 2),
-                np.linspace(0.0, 0.9, 5),
-                np.linspace(-np.pi, np.pi, 3),
-                np.arccos(np.linspace(0, 1, 5)[:-1]),
-            )
+    t = np.linspace(-100, 100, 1_000)
+
+    t0, period, a, e, omega, incl = (
+        x.flatten()
+        for x in np.meshgrid(
+            np.linspace(-5.0, 5.0, 2),
+            np.exp(np.linspace(np.log(5.0), np.log(50.0), 3)),
+            np.linspace(50.0, 100.0, 2),
+            np.linspace(0.0, 0.9, 5),
+            np.linspace(-np.pi, np.pi, 3),
+            np.arccos(np.linspace(0, 1, 5)[:-1]),
         )
+    )
 
-        r_batman = np.empty((len(t), len(t0)))
-
-        for i in range(len(t0)):
-            r_batman[:, i] = _rsky._rsky(
-                t, t0[i], period[i], a[i], incl[i], e[i], omega[i], 1, 1
-            )
-
-        m = r_batman < 100.0
+    r_batman = np.empty((len(t), len(t0)))
 
-        return {
-            "m_sum" : m.sum().item(), # Save native Int format
-            "r_batman" : r_batman,
-            "m" : m,
-            "t" : t,
-            "t0" : t0,
-            "period" : period,
-            "a" : a,
-            "e" : e,
-            "omega" : omega,
-            "incl" : incl,
-        }
+    for i in range(len(t0)):
+        r_batman[:, i] = _rsky._rsky(
+            t, t0[i], period[i], a[i], incl[i], e[i], omega[i], 1, 1
+        )
 
-    """
-    sky_coords = py"sky_coords"()
+    m = r_batman < 100.0
+
+    sky_coords = {
+        "m_sum" : m.sum().item(), # Save native Int format
+        "r_batman" : r_batman,
+        "m" : m,
+        "t" : t,
+        "t0" : t0,
+        "period" : period,
+        "a" : a,
+        "e" : e,
+        "omega" : omega,
+        "incl" : incl,
+    }
+    """ => sky_coords::Dict
 
     # Create comparison orbits from Transits.jl
     orbits = [
@@ -288,7 +281,8 @@ end
     )
 
     # Comparison coords from `batman`
-    py"""
+    @pyexec """
+    global np, _rsky, small_star
     import numpy as np
     from batman import _rsky
 
@@ -314,7 +308,8 @@ end
             "m": m,
         }
     """
-    small_star = py"small_star"(orbit.period, orbit.t0, orbit.aR_star, orbit.incl, orbit.ecc, orbit.omega)
+    small_star_py = @pyeval("small_star")(orbit.period, orbit.t0, orbit.aR_star, orbit.incl, orbit.ecc, orbit.omega)
+    small_star = pyconvert(Dict{String, Vector}, small_star_py)
 
     # Compare
     t = small_star["t"]