diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index e197d7b..75e1cd8 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -67,7 +67,7 @@ jobs:
           pip install .
       - name: Test with pytest
         run: |
-          pytest --full-trace --cov-report xml --cov enterprise_gwecc
+          pytest --full-trace --cov-report xml --cov test/test_enterprise_gwecc.py
       # - name: Codecov
       #   uses: codecov/codecov-action@v3
     
@@ -97,7 +97,6 @@ jobs:
       - name: Run examples
         run: |
           cd examples
-          julia benchmark.jl
           julia lbar_gammabar_vs_e.jl
           julia residuals_1psr.jl
           julia residuals_components.jl
diff --git a/Project.toml b/Project.toml
index e7d630e..066735b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,14 +6,20 @@ version = "0.1.1"
 [deps]
 Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
-FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
 JLD = "4138dd39-2aa7-5051-a626-17a0bb65d9c8"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-NumericalIntegration = "e7bfaba1-d571-5449-8927-abc22e82249b"
 PhysicalConstants = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab"
+ThreadsX = "ac1d9e8a-700a-412c-b207-f0111f4b6c0d"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
+UnitfulAstro = "6112ee07-acf9-5e0f-b108-d242c714bf9f"
+
+[extras]
+FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
+NumericalIntegration = "e7bfaba1-d571-5449-8927-abc22e82249b"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
-Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
-UnitfulAstro = "6112ee07-acf9-5e0f-b108-d242c714bf9f"
+
+[targets]
+test = ["FiniteDifferences", "NumericalIntegration", "Statistics", "Test", "UnPack"]
diff --git a/src/antennapattern.jl b/src/antennapattern.jl
index eacf061..3abe25a 100644
--- a/src/antennapattern.jl
+++ b/src/antennapattern.jl
@@ -41,7 +41,7 @@ function gw_polarization_tensors(pos::SkyLocation)
     e32c = cos(λ) * cos(β)
     e13c = -sin(λ) * cos(β)
     e23c = cos(λ) * cos(β)
-    e33c = 0
+    e33c = 0.0
 
     ep = [e11p e12p e13p; e21p e22p e23p; e31p e32p e33p]
     ec = [e11c e12c e13c; e21c e22c e23c; e31c e32c e33c]
diff --git a/src/residual.jl b/src/residual.jl
index 67303c6..0555f39 100644
--- a/src/residual.jl
+++ b/src/residual.jl
@@ -1,5 +1,7 @@
 export waveform, residual, residuals, residual_and_waveform, residuals_and_waveform
 
+import ThreadsX
+
 "PTA waveform."
 function waveform(
     mass::Mass,
@@ -99,13 +101,12 @@ function residuals(
     tref::Time,
     tEs::Vector{Time},
 )
-    dts = [tE - tref for tE in tEs]
-
     coeffs = EvolvCoeffs(mass, n_init, e_init)
     ap = AntennaPattern(psrpos, gwpos)
     Δp = pulsar_term_delay(ap, dp)
 
-    ss = [residual(mass, coeffs, l0p, proj, ap, terms, Δp, dt) for dt in dts]
+    resid_fn = tE -> residual(mass, coeffs, l0p, proj, ap, terms, Δp, tE - tref)
+    ss = ThreadsX.map(resid_fn, tEs)
 
     return ss
 end
diff --git a/test/test_enterprise_gwecc.py b/test/test_enterprise_gwecc.py
index c4607ef..26df25b 100644
--- a/test/test_enterprise_gwecc.py
+++ b/test/test_enterprise_gwecc.py
@@ -25,9 +25,12 @@
 )
 
 testdatadir = pathlib.Path(__file__).resolve().parent / "testdata"
-par = str(testdatadir / "J1909-3744_NANOGrav_12yv4.gls.par")
-tim = str(testdatadir / "J1909-3744_NANOGrav_12yv4.tim")
-psr = Pulsar(par, tim)
+
+@pytest.fixture
+def psr():
+    par = str(testdatadir / "J1909-3744_NANOGrav_12yv4.gls.par")
+    tim = str(testdatadir / "J1909-3744_NANOGrav_12yv4.tim")
+    return Pulsar(par, tim, timing_package="pint")
 
 
 year = 365.25 * 24 * 3600
@@ -136,7 +139,7 @@ def test_eccentric_pta_signal_target(psrTerm, spline):
     "psrTerm, tie_psrTerm, spline",
     it.product([True, False], [True, False], [True, False]),
 )
-def test_gwecc_block(psrTerm, tie_psrTerm, spline):
+def test_gwecc_block(psr, psrTerm, tie_psrTerm, spline):
     tref = max(psr.toas)
 
     tm = MarginalizingTimingModel()
@@ -180,7 +183,7 @@ def test_gwecc_block(psrTerm, tie_psrTerm, spline):
 
 
 @pytest.mark.parametrize("psrTerm, spline", it.product([True, False], [True, False]))
-def test_gwecc_1psr_block(psrTerm, spline):
+def test_gwecc_1psr_block(psr, psrTerm, spline):
     tref = max(psr.toas)
 
     tm = MarginalizingTimingModel()
@@ -218,7 +221,7 @@ def test_gwecc_1psr_block(psrTerm, spline):
     "psrTerm, tie_psrTerm, spline",
     it.product([True, False], [True, False], [True, False]),
 )
-def test_gwecc_target_block(psrTerm, tie_psrTerm, spline):
+def test_gwecc_target_block(psr, psrTerm, tie_psrTerm, spline):
     tref = max(psr.toas)
 
     tm = MarginalizingTimingModel()
@@ -266,7 +269,7 @@ def test_gwecc_target_block(psrTerm, tie_psrTerm, spline):
     assert lnprior_fn(x0) == -np.inf
 
 
-def test_psrdist_prior():
+def test_psrdist_prior(psr):
     parfiles = sorted(glob.glob(f"{testdatadir}/*.par"))
     timfiles = sorted(glob.glob(f"{testdatadir}/*.tim"))