diff --git a/README.md b/README.md
index 8040677f8..846d1a959 100644
--- a/README.md
+++ b/README.md
@@ -49,9 +49,12 @@ the figures in
 `SLEPLET` has been tested with
 [![Python](https://img.shields.io/pypi/pyversions/sleplet)](https://www.python.org).
 Windows is not currently supported as `SLEPLET` relies on
-[pyssht](https://pypi.org/project/pyssht) which does not work on Windows. This
-can hopefully be replaced with [s2fft](https://github.com/astro-informatics/s2fft)
-in the future when it is available on [PyPI](https://pypi.org).
+[pyssht](https://pypi.org/project/pyssht) and
+[pys2let](https://pypi.org/project/pys2let) which do not work on Windows.
+These can hopefully be replaced with
+[s2fft](https://github.com/astro-informatics/s2fft) and
+[s2wav](https://github.com/astro-informatics/s2wav) in the future when they
+are available on [PyPI](https://pypi.org).
 
 ## Example Usage
 
diff --git a/documentation/DOCUMENTATION.md b/documentation/DOCUMENTATION.md
index 032ce27d4..b345c2b71 100644
--- a/documentation/DOCUMENTATION.md
+++ b/documentation/DOCUMENTATION.md
@@ -53,13 +53,13 @@ done
 
 ```python
 import numpy as np
-import s2fft
+import pyssht as ssht
 import sleplet
 
 for ell in range(2, 0, -1):
     f = sleplet.functions.HarmonicGaussian(128, l_sigma=10**ell, m_sigma=10)
     flm = f.translate(alpha=0.75 * np.pi, beta=0.125 * np.pi)
-    f_sphere = s2fft.inverse(flm, f.L, method="jax", sampling="mwss")
+    f_sphere = ssht.inverse(flm, f.L, Method="MWSS")
     sleplet.plotting.PlotSphere(
         f_sphere,
         f.L,
@@ -75,12 +75,12 @@ sphere earth -L 128
 ```
 
 ```python
-import s2fft
+import pyssht as ssht
 import sleplet
 
 f = sleplet.functions.Earth(128)
 flm = sleplet.harmonic_methods.rotate_earth_to_south_america(f.coefficients, f.L)
-f_sphere = s2fft.inverse(flm, f.L, method="jax", sampling="mwss")
+f_sphere = ssht.inverse(flm, f.L, Method="MWSS")
 sleplet.plotting.PlotSphere(f_sphere, f.L, "fig_2").execute()
 ```
 
@@ -93,7 +93,7 @@ done
 ```
 
 ```python
-import s2fft
+import pyssht as ssht
 import sleplet
 
 for ell in range(2, 0, -1):
@@ -101,7 +101,7 @@ for ell in range(2, 0, -1):
     g = sleplet.functions.Earth(128)
     flm = f.convolve(f.coefficients, g.coefficients.conj())
     flm_rot = sleplet.harmonic_methods.rotate_earth_to_south_america(flm, f.L)
-    f_sphere = s2fft.inverse(flm_rot, f.L, method="jax", sampling="mwss")
+    f_sphere = ssht.inverse(flm_rot, f.L, Method="MWSS")
     sleplet.plotting.PlotSphere(f_sphere, f.L, f"fig_3_ell_{ell}").execute()
 ```
 
@@ -126,13 +126,13 @@ sphere slepian_south_america -L 128 -s 2 -u
 ```
 
 ```python
-import s2fft
+import pyssht as ssht
 import sleplet
 
 # a
 f = sleplet.functions.Earth(128, smoothing=2)
 flm = sleplet.harmonic_methods.rotate_earth_to_south_america(f.coefficients, f.L)
-f_sphere = s2fft.inverse(flm, f.L, method="jax", sampling="mwss")
+f_sphere = ssht.inverse(flm, f.L, Method="MWSS")
 sleplet.plotting.PlotSphere(f_sphere, f.L, "fig_3_a", normalise=False).execute()
 # b
 region = sleplet.slepian.Region(mask_name="south_america")
@@ -282,13 +282,13 @@ sphere slepian_africa -L 128 -s 2 -u
 ```
 
 ```python
-import s2fft
+import pyssht as ssht
 import sleplet
 
 # a
 f = sleplet.functions.Earth(128, smoothing=2)
 flm = sleplet.harmonic_methods.rotate_earth_to_africa(f.coefficients, f.L)
-f_sphere = s2fft.inverse(flm, f.L, method="jax", sampling="mwss")
+f_sphere = ssht.inverse(flm, f.L, Method="MWSS")
 sleplet.plotting.PlotSphere(f_sphere, f.L, "fig_9_a", normalise=False).execute()
 # b
 region = sleplet.slepian.Region(mask_name="africa")
@@ -639,13 +639,13 @@ done
 ```
 
 ```python
-import s2fft
+import pyssht as ssht
 import sleplet
 
 for ell in range(5):
     for m in range(ell + 1):
         f = sleplet.functions.SphericalHarmonic(128, ell=ell, m=m)
-        f_sphere = s2fft.inverse(f.coefficients, f.L, method="jax", sampling="mwss")
+        f_sphere = ssht.inverse(f.coefficients, f.L, Method="MWSS")
         sleplet.plotting.PlotSphere(
             f_sphere,
             f.L,
@@ -670,17 +670,17 @@ sphere elongated_gaussian -e -1 -1 -L 128 -m rotate -a 0.25 -b 0.25 -g 0.25
 
 ```python
 import numpy as np
-import s2fft
+import pyssht as ssht
 import sleplet
 
 # a
 f = sleplet.functions.ElongatedGaussian(128, p_sigma=0.1, t_sigma=0.1)
-f_sphere = s2fft.inverse(f.coefficients, f.L, method="jax", sampling="mwss")
+f_sphere = ssht.inverse(f.coefficients, f.L, Method="MWSS")
 sleplet.plotting.PlotSphere(f_sphere, f.L, "fig_2_2_a", annotations=[]).execute()
 # b-d
 for a, b, g in [(0, 0, 0.25), (0, 0.25, 0.25), (0.25, 0.25, 0.25)]:
     glm_rot = f.rotate(alpha=a * np.pi, beta=b * np.pi, gamma=g * np.pi)
-    g_sphere = s2fft.inverse(glm_rot, f.L, method="jax", sampling="mwss")
+    g_sphere = ssht.inverse(glm_rot, f.L, Method="MWSS")
     sleplet.plotting.PlotSphere(
         g_sphere,
         f.L,
@@ -713,12 +713,12 @@ done
 ```
 
 ```python
-import s2fft
+import pyssht as ssht
 import sleplet
 
 for j in [None, *list(range(4))]:
     f = sleplet.functions.AxisymmetricWavelets(128, B=3, j_min=2, j=j)
-    f_sphere = s2fft.inverse(f.coefficients, f.L, method="jax", sampling="mwss")
+    f_sphere = ssht.inverse(f.coefficients, f.L, Method="MWSS")
     sleplet.plotting.PlotSphere(
         f_sphere,
         f.L,
@@ -758,16 +758,16 @@ sphere gaussian -a 0.75 -b 0.125 -L 128 -m translate -o
 
 ```python
 import numpy as np
-import s2fft
+import pyssht as ssht
 import sleplet
 
 # a
 f = sleplet.functions.Gaussian(128)
-f_sphere = s2fft.inverse(f.coefficients, f.L, method="jax", sampling="mwss")
+f_sphere = ssht.inverse(f.coefficients, f.L, Method="MWSS")
 sleplet.plotting.PlotSphere(f_sphere, f.L, "fig_3_1_a", annotations=[]).execute()
 # b
 glm_trans = f.translate(alpha=0.75 * np.pi, beta=0.125 * np.pi)
-g_sphere = s2fft.inverse(glm_trans, f.L, method="jax", sampling="mwss")
+g_sphere = ssht.inverse(glm_trans, f.L, Method="MWSS")
 sleplet.plotting.PlotSphere(g_sphere, f.L, "fig_3_1_b", annotations=[]).execute()
 ```
 
@@ -782,16 +782,16 @@ sphere squashed_gaussian -a 0.75 -b 0.125 -L 128 -m translate -o
 
 ```python
 import numpy as np
-import s2fft
+import pyssht as ssht
 import sleplet
 
 # a
 f = sleplet.functions.SquashedGaussian(128)
-f_sphere = s2fft.inverse(f.coefficients, f.L, method="jax", sampling="mwss")
+f_sphere = ssht.inverse(f.coefficients, f.L, Method="MWSS")
 sleplet.plotting.PlotSphere(f_sphere, f.L, "fig_3_2_a", annotations=[]).execute()
 # b
 glm_trans = f.translate(alpha=0.75 * np.pi, beta=0.125 * np.pi)
-g_sphere = s2fft.inverse(glm_trans, f.L, method="jax", sampling="mwss")
+g_sphere = ssht.inverse(glm_trans, f.L, Method="MWSS")
 sleplet.plotting.PlotSphere(g_sphere, f.L, "fig_3_2_b", annotations=[]).execute()
 ```
 
@@ -806,16 +806,16 @@ sphere elongated_gaussian -a 0.75 -b 0.125 -L 128 -m translate -o
 
 ```python
 import numpy as np
-import s2fft
+import pyssht as ssht
 import sleplet
 
 # a
 f = sleplet.functions.ElongatedGaussian(128)
-f_sphere = s2fft.inverse(f.coefficients, f.L, method="jax", sampling="mwss")
+f_sphere = ssht.inverse(f.coefficients, f.L, Method="MWSS")
 sleplet.plotting.PlotSphere(f_sphere, f.L, "fig_3_3_a", annotations=[]).execute()
 # b
 glm_trans = f.translate(alpha=0.75 * np.pi, beta=0.125 * np.pi)
-g_sphere = s2fft.inverse(glm_trans, f.L, method="jax", sampling="mwss")
+g_sphere = ssht.inverse(glm_trans, f.L, Method="MWSS")
 sleplet.plotting.PlotSphere(g_sphere, f.L, "fig_3_3_b", annotations=[]).execute()
 ```
 
@@ -830,12 +830,12 @@ done
 ```
 
 ```python
-import s2fft
+import pyssht as ssht
 import sleplet
 
 for ell in range(2, 0, -1):
     f = sleplet.functions.HarmonicGaussian(128, l_sigma=10**ell, m_sigma=10)
-    f_sphere = s2fft.inverse(f.coefficients, f.L, method="jax", sampling="mwss")
+    f_sphere = ssht.inverse(f.coefficients, f.L, Method="MWSS")
     sleplet.plotting.PlotSphere(
         f_sphere,
         f.L,
diff --git a/examples/arbitrary/_slepian_wavelet_covariance.py b/examples/arbitrary/_slepian_wavelet_covariance.py
index 862ea4e58..cf5fef8e8 100644
--- a/examples/arbitrary/_slepian_wavelet_covariance.py
+++ b/examples/arbitrary/_slepian_wavelet_covariance.py
@@ -1,11 +1,11 @@
 import numpy as np
 import numpy.typing as npt
 
-import s2fft
+import pyssht as ssht
 
 import sleplet
 
-SAMPLING_SCHEME = "mwss"
+SAMPLING_SCHEME = "MWSS"
 
 
 def compute_slepian_wavelet_covariance(
@@ -29,7 +29,7 @@ def compute_slepian_wavelet_covariance(
         covariance.shape[0],
     )
     np.testing.assert_equal(
-        s2fft.samples.f_shape(slepian_wavelets.L, sampling=SAMPLING_SCHEME),
+        ssht.sample_shape(slepian_wavelets.L, Method=SAMPLING_SCHEME),
         covariance.shape[1:],
     )
     return covariance
diff --git a/examples/arbitrary/africa/slepian_wavelet_covariance_africa.py b/examples/arbitrary/africa/slepian_wavelet_covariance_africa.py
index 92f2910e2..a0f1b38e6 100644
--- a/examples/arbitrary/africa/slepian_wavelet_covariance_africa.py
+++ b/examples/arbitrary/africa/slepian_wavelet_covariance_africa.py
@@ -3,8 +3,6 @@
 
 import numpy as np
 
-import s2fft
-
 import sleplet
 
 sys.path.append(str(pathlib.Path(__file__).resolve().parents[1]))
@@ -40,11 +38,8 @@ def main() -> None:
     rng = np.random.default_rng(RANDOM_SEED)
 
     for i in range(RUNS):
-        # Generate random complex signal
-        f_p = s2fft.samples.flm_2d_to_1d(
-            s2fft.utils.signal_generator.generate_flm(rng, L),
-            L,
-        )
+        # Generate normally distributed random complex signal
+        f_p = sleplet.harmonic_methods.compute_random_signal(L, rng, var_signal=VAR_FP)
 
         # compute wavelet coefficients
         w_p = sleplet.wavelet_methods.slepian_wavelet_forward(
diff --git a/examples/arbitrary/africa/sparsity_wavelet_slepian_comparison_africa.py b/examples/arbitrary/africa/sparsity_wavelet_slepian_comparison_africa.py
index 5661aa5fd..4f2c35bed 100644
--- a/examples/arbitrary/africa/sparsity_wavelet_slepian_comparison_africa.py
+++ b/examples/arbitrary/africa/sparsity_wavelet_slepian_comparison_africa.py
@@ -43,11 +43,7 @@ def _plot_axisymmetric_coefficients(shannon: int) -> None:
     awc = sleplet.functions.AxisymmetricWaveletCoefficientsAfrica(L, B=B, j_min=J_MIN)
 
     # find sorted coefficients
-    w_lm = np.sort(np.sort(np.abs(awc.wavelet_coefficients), axis=1), axis=2)[
-        :,
-        ::-1,
-        ::-1,
-    ]
+    w_lm = np.sort(np.abs(awc.wavelet_coefficients), axis=1)[:, ::-1]
 
     # perform plot
     plt.plot(w_lm[0, :shannon], "--", label=r"$|W^{{\phi}}_{\ell m}|$")
diff --git a/examples/arbitrary/south_america/slepian_wavelet_covariance_south_america.py b/examples/arbitrary/south_america/slepian_wavelet_covariance_south_america.py
index 939072fdc..4ef6018da 100644
--- a/examples/arbitrary/south_america/slepian_wavelet_covariance_south_america.py
+++ b/examples/arbitrary/south_america/slepian_wavelet_covariance_south_america.py
@@ -3,8 +3,6 @@
 
 import numpy as np
 
-import s2fft
-
 import sleplet
 
 sys.path.append(str(pathlib.Path(__file__).resolve().parents[1]))
@@ -40,11 +38,8 @@ def main() -> None:
     rng = np.random.default_rng(RANDOM_SEED)
 
     for i in range(RUNS):
-        # Generate random complex signal
-        f_p = s2fft.samples.flm_2d_to_1d(
-            s2fft.utils.signal_generator.generate_flm(rng, L),
-            L,
-        )
+        # Generate normally distributed random complex signal
+        f_p = sleplet.harmonic_methods.compute_random_signal(L, rng, var_signal=VAR_FP)
 
         # compute wavelet coefficients
         w_p = sleplet.wavelet_methods.slepian_wavelet_forward(
diff --git a/examples/arbitrary/south_america/sparsity_wavelet_slepian_comparison_south_america.py b/examples/arbitrary/south_america/sparsity_wavelet_slepian_comparison_south_america.py
index d6131e988..e72baf12b 100644
--- a/examples/arbitrary/south_america/sparsity_wavelet_slepian_comparison_south_america.py
+++ b/examples/arbitrary/south_america/sparsity_wavelet_slepian_comparison_south_america.py
@@ -47,11 +47,7 @@ def _plot_axisymmetric_coefficients(shannon: int) -> None:
     )
 
     # find sorted coefficients
-    w_lm = np.sort(np.sort(np.abs(awc.wavelet_coefficients), axis=1), axis=2)[
-        :,
-        ::-1,
-        ::-1,
-    ]
+    w_lm = np.sort(np.abs(awc.wavelet_coefficients), axis=1)[:, ::-1]
 
     # perform plot
     plt.plot(w_lm[0, :shannon], "--", label=r"$|W^{{\phi}}_{\ell m}|$")
diff --git a/examples/misc/_denoising_axisym.py b/examples/misc/_denoising_axisym.py
index c1568c1bb..210f0e255 100644
--- a/examples/misc/_denoising_axisym.py
+++ b/examples/misc/_denoising_axisym.py
@@ -1,12 +1,11 @@
 import numpy as np
 import numpy.typing as npt
 
-import s2fft
+import pyssht as ssht
 
 import sleplet
 
-EXECUTION_MODE = "jax"
-SAMPLING_SCHEME = "mwss"
+SAMPLING_SCHEME = "MWSS"
 
 
 def denoising_axisym(  # noqa: PLR0913
@@ -61,5 +60,5 @@ def denoising_axisym(  # noqa: PLR0913
         else flm
     )
 
-    f = s2fft.inverse(flm, signal.L, method=EXECUTION_MODE, sampling=SAMPLING_SCHEME)
+    f = ssht.inverse(flm, signal.L, Method=SAMPLING_SCHEME)
     return f, noised_signal.snr, denoised_snr
diff --git a/examples/misc/translation_normalisation.py b/examples/misc/translation_normalisation.py
index 8eabed718..b43d76df6 100644
--- a/examples/misc/translation_normalisation.py
+++ b/examples/misc/translation_normalisation.py
@@ -3,7 +3,6 @@
 import seaborn as sns
 
 import pyssht as ssht
-import s2fft
 
 import sleplet
 
@@ -11,13 +10,13 @@
 
 ALPHA_DEFAULT = 0.75
 L = 128
-SAMPLING_SCHEME = "mwss"
+SAMPLING_SCHEME = "MWSS"
 
 
 def compute_translation_normalisation_theta() -> None:
     """Analysis of the translation norm for referee."""
     hg = sleplet.functions.HarmonicGaussian(L)
-    thetas = s2fft.samples.thetas(L, sampling=SAMPLING_SCHEME)
+    thetas, _ = ssht.sample_positions(L, Method=SAMPLING_SCHEME)
     norm = np.zeros(len(thetas))
     for i, theta in enumerate(thetas):
         print(f"compute norm {i+1}/{len(thetas)}")
diff --git a/examples/polar_cap/slepian_coefficients.py b/examples/polar_cap/slepian_coefficients.py
index 1e2b5dbe6..731866249 100644
--- a/examples/polar_cap/slepian_coefficients.py
+++ b/examples/polar_cap/slepian_coefficients.py
@@ -3,8 +3,6 @@
 import numpy.typing as npt
 import seaborn as sns
 
-import s2fft
-
 import sleplet
 
 sns.set(context="paper")
@@ -20,7 +18,7 @@ def _earth_region_harmonic_coefficients(
     """Harmonic coefficients of the Earth for the polar cap region."""
     region = sleplet.slepian.Region(theta_max=np.deg2rad(theta_max))
     earth = sleplet.functions.Earth(L, region=region)
-    coefficients = np.abs(s2fft.samples.flm_2d_to_1d(earth.coefficients, L))
+    coefficients = np.abs(earth.coefficients)
     coefficients[::-1].sort()
     return coefficients
 
diff --git a/examples/polar_cap/slepian_error.py b/examples/polar_cap/slepian_error.py
index ac97b90dd..9371ee401 100644
--- a/examples/polar_cap/slepian_error.py
+++ b/examples/polar_cap/slepian_error.py
@@ -3,15 +3,14 @@
 import numpy.typing as npt
 import seaborn as sns
 
-import s2fft
+import pyssht as ssht
 
 import sleplet
 
 sns.set(context="paper")
 
-EXECUTION_MODE = "jax"
 L = 16
-SAMPLING_SCHEME = "mwss"
+SAMPLING_SCHEME = "MWSS"
 THETA_MAX = 40
 
 
@@ -19,12 +18,7 @@ def main() -> None:
     """Creates a plot of Slepian coefficients against rank."""
     region = sleplet.slepian.Region(theta_max=np.deg2rad(THETA_MAX))
     earth = sleplet.functions.Earth(L, region=region)
-    field = s2fft.inverse(
-        earth.coefficients,
-        L,
-        method=EXECUTION_MODE,
-        sampling=SAMPLING_SCHEME,
-    )
+    field = ssht.inverse(earth.coefficients, L, Method=SAMPLING_SCHEME)
     integrate_region = _helper_region(L, region, field, earth.coefficients)
     integrate_sphere = _helper_sphere(L, region, field, earth.coefficients)
     N = sleplet.slepian.SlepianPolarCap(L, np.deg2rad(THETA_MAX)).N
diff --git a/examples/wavelets/axisymmetric_wavelet_covariance.py b/examples/wavelets/axisymmetric_wavelet_covariance.py
index d7ae13ccf..bbfa0a252 100644
--- a/examples/wavelets/axisymmetric_wavelet_covariance.py
+++ b/examples/wavelets/axisymmetric_wavelet_covariance.py
@@ -1,16 +1,15 @@
 import numpy as np
 import numpy.typing as npt
 
-import s2fft
+import pyssht as ssht
 
 import sleplet
 
 B = 3
-EXECUTION_MODE = "jax"
 J_MIN = 2
 L = 128
 RANDOM_SEED = 30
-SAMPLING_SCHEME = "mwss"
+SAMPLING_SCHEME = "MWSS"
 
 
 def _compute_wavelet_covariance(
@@ -19,7 +18,7 @@ def _compute_wavelet_covariance(
     var_signal: float,
 ) -> npt.NDArray[np.float_]:
     """Computes the theoretical covariance of the wavelet coefficients."""
-    covar_theory = (np.abs(wavelets) ** 2).sum(axis=(1, 2))
+    covar_theory = (np.abs(wavelets) ** 2).sum(axis=1)
     return covar_theory * var_signal
 
 
@@ -76,20 +75,15 @@ def axisymmetric_wavelet_covariance(
     for i in range(runs):
         print(f"start run: {i+1}/{runs}")
 
-        # Generate random complex signal
-        flm = s2fft.utils.signal_generator.generate_flm(rng, L)
+        # Generate normally distributed random complex signal
+        flm = sleplet.harmonic_methods.compute_random_signal(L, rng, var_signal=var_flm)
 
         # compute wavelet coefficients
         wlm = sleplet.wavelet_methods.axisymmetric_wavelet_forward(L, flm, aw.wavelets)
 
         # compute covariance from data
         for j, coefficient in enumerate(wlm):
-            f_wav_j = s2fft.inverse(
-                coefficient,
-                L,
-                method=EXECUTION_MODE,
-                sampling=SAMPLING_SCHEME,
-            )
+            f_wav_j = ssht.inverse(coefficient, L, Method=SAMPLING_SCHEME)
             covar_data[i, j] = (
                 f_wav_j.var() if _is_ergodic(j_min, j=j) else f_wav_j[0, 0]
             )
diff --git a/pyproject.toml b/pyproject.toml
index f694086c2..8dbc9e997 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -26,7 +26,6 @@ dependencies = [
     "cmocean>=3.0.3",
     "gmpy2>=2.1.5",
     "hypothesis>=6.70.2",
-    "jaxlib>=0.4.14",
     "libigl>=2.4.1",
     "matplotlib>=3.7.1",
     "multiprocess>=0.70.14",
@@ -36,8 +35,8 @@ dependencies = [
     "plotly>=5.14.0",
     "pooch>=1.7.0",
     "pydantic>=2.3.0",
+    "pys2let>=2.2.6",
     "pyssht>=1.5.2",
-    "s2wav@git+https://github.com/astro-informatics/s2wav.git@main#egg=s2wav",
     "scipy>=1.10.1",
     "seaborn>=0.12.2",
     "tomli>=2.0.1",
@@ -121,9 +120,8 @@ isort = {known-first-party = [
     "first-party",
     "local-folder",
 ], sections = {"astro-info" = [
+    "pys2let",
     "pyssht",
-    "s2fft",
-    "s2wav",
 ]}}
 per-file-ignores = {"examples/*" = [
     "D100",
diff --git a/src/sleplet/_data/create_earth_flm.py b/src/sleplet/_data/create_earth_flm.py
index 7f5377640..d09bed72f 100644
--- a/src/sleplet/_data/create_earth_flm.py
+++ b/src/sleplet/_data/create_earth_flm.py
@@ -3,7 +3,7 @@
 import numpy.typing as npt
 import scipy.io as sio
 
-import s2fft
+import pyssht as ssht
 
 import sleplet._data.setup_pooch
 import sleplet._smoothing
@@ -14,18 +14,16 @@ def create_flm(L: int, *, smoothing: int | None = None) -> npt.NDArray[np.comple
     # load in data
     flm = _load_flm()
 
-    # don't take the full L, and convert to 2D
-    flm = s2fft.samples.flm_1d_to_2d(flm[: L**2], L)
-
     # fill in negative m components so as to avoid confusion with zero values
     for ell in range(1, L):
         for m in range(1, ell + 1):
-            ind_pm = L - 1 + m
-            ind_nm = L - 1 - m
-            flm[ell, ind_nm] = (-1) ** m * flm[ell, ind_pm].conj()
+            ind_pm = ssht.elm2ind(ell, m)
+            ind_nm = ssht.elm2ind(ell, -m)
+            flm_pm = flm[ind_pm]
+            flm[ind_nm] = (-1) ** m * flm_pm.conj()
 
-    # invert dataset as Earth backwards
-    flm = flm.conj()
+    # don't take the full L, invert dataset as Earth backwards
+    flm = flm[: L**2].conj()
 
     if isinstance(smoothing, int):
         flm = sleplet._smoothing.apply_gaussian_smoothing(flm, L, smoothing)
diff --git a/src/sleplet/_data/create_wmap_flm.py b/src/sleplet/_data/create_wmap_flm.py
index 9659366f9..bac442774 100644
--- a/src/sleplet/_data/create_wmap_flm.py
+++ b/src/sleplet/_data/create_wmap_flm.py
@@ -3,7 +3,7 @@
 import numpy.typing as npt
 import scipy.io as sio
 
-import s2fft
+import pyssht as ssht
 
 import sleplet._data.setup_pooch
 import sleplet._vars
@@ -18,19 +18,20 @@ def create_flm(L: int) -> npt.NDArray[np.complex_]:
     rng = np.random.default_rng(sleplet._vars.RANDOM_SEED)
 
     # Simulate CMB in harmonic space.
-    flm = np.zeros(s2fft.samples.flm_shape(L), dtype=np.complex_)
+    flm = np.zeros(L**2, dtype=np.complex_)
     for ell in range(2, L):
         sigma = np.sqrt(2 * np.pi / (ell * (ell + 1)) * cl[ell - 2])
-        flm[ell, L - 1] = sigma * rng.standard_normal()
+        ind = ssht.elm2ind(ell, 0)
+        flm[ind] = sigma * rng.standard_normal()
         for m in range(1, ell + 1):
-            ind_pm = L - 1 + m
-            ind_nm = L - 1 - m
-            flm[ell, ind_pm] = (
+            ind_pm = ssht.elm2ind(ell, m)
+            ind_nm = ssht.elm2ind(ell, -m)
+            flm[ind_pm] = (
                 sigma
                 / np.sqrt(2)
                 * (rng.standard_normal() + 1j * rng.standard_normal())
             )
-            flm[ell, ind_nm] = (-1) ** m * flm[ell, ind_pm].conj()
+            flm[ind_nm] = (-1) ** m * flm[ind_pm].conj()
     return flm
 
 
diff --git a/src/sleplet/_integration_methods.py b/src/sleplet/_integration_methods.py
index 190156e11..fae0e8402 100644
--- a/src/sleplet/_integration_methods.py
+++ b/src/sleplet/_integration_methods.py
@@ -4,20 +4,17 @@
 import numpy as np
 import numpy.typing as npt
 
-import s2fft
+import pyssht as ssht
 
 import sleplet._vars
 
 
 def calc_integration_weight(L: int) -> npt.NDArray[np.float_]:
     """Computes the spherical Jacobian for the integration."""
-    thetas = np.tile(
-        s2fft.samples.thetas(L, sampling=sleplet._vars.SAMPLING_SCHEME),
-        (s2fft.samples.nphi_equiang(L, sampling=sleplet._vars.SAMPLING_SCHEME), 1),
-    ).T
-    phis = np.tile(
-        s2fft.samples.phis_equiang(L, sampling=sleplet._vars.SAMPLING_SCHEME),
-        (s2fft.samples.ntheta(L, sampling=sleplet._vars.SAMPLING_SCHEME), 1),
+    thetas, phis = ssht.sample_positions(
+        L,
+        Grid=True,
+        Method=sleplet._vars.SAMPLING_SCHEME,
     )
     delta_theta = np.ediff1d(thetas[:, 0]).mean()
     delta_phi = np.ediff1d(phis[0]).mean()
diff --git a/src/sleplet/_mask_methods.py b/src/sleplet/_mask_methods.py
index 9b908ad99..401933176 100644
--- a/src/sleplet/_mask_methods.py
+++ b/src/sleplet/_mask_methods.py
@@ -5,7 +5,7 @@
 import numpy.typing as npt
 import platformdirs
 
-import s2fft
+import pyssht as ssht
 
 import sleplet._data.create_earth_flm
 import sleplet._data.setup_pooch
@@ -32,13 +32,10 @@ def create_mask_region(
                                    phi_min or phi_max is provided
     * arbitrary - just checks the shape of the input mask.
     """
-    thetas = np.tile(
-        s2fft.samples.thetas(L, sampling=sleplet._vars.SAMPLING_SCHEME),
-        (s2fft.samples.nphi_equiang(L, sampling=sleplet._vars.SAMPLING_SCHEME), 1),
-    ).T
-    phis = np.tile(
-        s2fft.samples.phis_equiang(L, sampling=sleplet._vars.SAMPLING_SCHEME),
-        (s2fft.samples.ntheta(L, sampling=sleplet._vars.SAMPLING_SCHEME), 1),
+    thetas, phis = ssht.sample_positions(
+        L,
+        Grid=True,
+        Method=sleplet._vars.SAMPLING_SCHEME,
     )
 
     match region.region_type:
@@ -86,23 +83,21 @@ def ensure_masked_flm_bandlimited(
     spin: int,
 ) -> npt.NDArray[np.complex_]:
     """Ensures the coefficients is bandlimited for a given region."""
-    field = s2fft.inverse(
+    field = ssht.inverse(
         flm,
         L,
-        method=sleplet._vars.EXECUTION_MODE,
-        reality=reality,
-        sampling=sleplet._vars.SAMPLING_SCHEME,
-        spin=spin,
+        Method=sleplet._vars.SAMPLING_SCHEME,
+        Reality=reality,
+        Spin=spin,
     )
     mask = create_mask_region(L, region)
     field = np.where(mask, field, 0)
-    return s2fft.forward(
+    return ssht.forward(
         field,
         L,
-        method=sleplet._vars.EXECUTION_MODE,
-        reality=reality,
-        sampling=sleplet._vars.SAMPLING_SCHEME,
-        spin=spin,
+        Method=sleplet._vars.SAMPLING_SCHEME,
+        Reality=reality,
+        Spin=spin,
     )
 
 
@@ -160,17 +155,17 @@ def _create_africa_mask(
 ) -> npt.NDArray[np.float_]:
     """Creates the Africa region mask."""
     rot_flm = sleplet.harmonic_methods.rotate_earth_to_africa(earth_flm, L)
-    earth_f = s2fft.inverse(
+    earth_f = ssht.inverse(
         rot_flm,
         L,
-        method=sleplet._vars.EXECUTION_MODE,
-        reality=True,
-        sampling=sleplet._vars.SAMPLING_SCHEME,
+        Method=sleplet._vars.SAMPLING_SCHEME,
+        Reality=True,
+    )
+    thetas, _ = ssht.sample_positions(
+        L,
+        Grid=True,
+        Method=sleplet._vars.SAMPLING_SCHEME,
     )
-    thetas = np.tile(
-        s2fft.samples.thetas(L, sampling=sleplet._vars.SAMPLING_SCHEME),
-        (s2fft.samples.nphi_equiang(L, sampling=sleplet._vars.SAMPLING_SCHEME), 1),
-    ).T
     return (thetas <= _AFRICA_RANGE) & (earth_f >= 0)
 
 
@@ -180,17 +175,17 @@ def _create_south_america_mask(
 ) -> npt.NDArray[np.float_]:
     """Creates the Africa region mask."""
     rot_flm = sleplet.harmonic_methods.rotate_earth_to_south_america(earth_flm, L)
-    earth_f = s2fft.inverse(
+    earth_f = ssht.inverse(
         rot_flm,
         L,
-        method=sleplet._vars.EXECUTION_MODE,
-        reality=True,
-        sampling=sleplet._vars.SAMPLING_SCHEME,
+        Method=sleplet._vars.SAMPLING_SCHEME,
+        Reality=True,
+    )
+    thetas, _ = ssht.sample_positions(
+        L,
+        Grid=True,
+        Method=sleplet._vars.SAMPLING_SCHEME,
     )
-    thetas = np.tile(
-        s2fft.samples.thetas(L, sampling=sleplet._vars.SAMPLING_SCHEME),
-        (s2fft.samples.nphi_equiang(L, sampling=sleplet._vars.SAMPLING_SCHEME), 1),
-    ).T
     return (thetas <= _SOUTH_AMERICA_RANGE) & (earth_f >= 0)
 
 
diff --git a/src/sleplet/_smoothing.py b/src/sleplet/_smoothing.py
index 261ff293a..f37e2a1f8 100644
--- a/src/sleplet/_smoothing.py
+++ b/src/sleplet/_smoothing.py
@@ -4,7 +4,6 @@
 import numpy.typing as npt
 
 import pyssht as ssht
-import s2fft
 
 _logger = logging.getLogger(__name__)
 
@@ -24,7 +23,4 @@ def apply_gaussian_smoothing(
     sigma = np.pi / (smoothing_factor * L)
     fwhm = 2 * np.sqrt(np.log(2)) * sigma
     _logger.info(f"FWHM = {np.rad2deg(fwhm):.2f}degrees")
-    return s2fft.samples.flm_1d_to_2d(
-        ssht.gaussian_smoothing(s2fft.samples.flm_2d_to_1d(flm, L), L, sigma),
-        L,
-    )
+    return ssht.gaussian_smoothing(flm, L, sigma)
diff --git a/src/sleplet/_vars.py b/src/sleplet/_vars.py
index c6f31accf..72feb92c7 100644
--- a/src/sleplet/_vars.py
+++ b/src/sleplet/_vars.py
@@ -1,11 +1,10 @@
 import numpy as np
 
-EXECUTION_MODE = "jax"
 PHI_0 = np.pi
 PHI_MAX_DEFAULT = 2 * np.pi
 PHI_MIN_DEFAULT = 0.0
 RANDOM_SEED = 30
-SAMPLING_SCHEME = "mwss"
+SAMPLING_SCHEME = "MWSS"
 SPHERE_UNSEEN = -1.56e30
 THETA_0 = 0.0
 THETA_MAX_DEFAULT = np.pi
diff --git a/src/sleplet/functions/africa.py b/src/sleplet/functions/africa.py
index 55c5925b2..e5755d63f 100644
--- a/src/sleplet/functions/africa.py
+++ b/src/sleplet/functions/africa.py
@@ -3,7 +3,7 @@
 import numpy.typing as npt
 import pydantic.v1 as pydantic
 
-import s2fft
+import pyssht as ssht
 
 import sleplet._data.create_earth_flm
 import sleplet._data.setup_pooch
@@ -61,12 +61,11 @@ def _grid_fun(
             smoothing=self.smoothing,
         )
         rot_flm = sleplet.harmonic_methods.rotate_earth_to_africa(earth_flm, self.L)
-        earth_f = s2fft.inverse(
+        earth_f = ssht.inverse(
             rot_flm,
             self.L,
-            method=sleplet._vars.EXECUTION_MODE,
-            reality=self.reality,
-            sampling=sleplet._vars.SAMPLING_SCHEME,
+            Method=sleplet._vars.SAMPLING_SCHEME,
+            Reality=self.reality,
         )
         mask_name = f"{self.name}_L{self.L}.npy"
         mask_location = sleplet._data.setup_pooch.find_on_pooch_then_local(
diff --git a/src/sleplet/functions/axisymmetric_wavelet_coefficients_africa.py b/src/sleplet/functions/axisymmetric_wavelet_coefficients_africa.py
index 81c9b5cd8..4fdd0818f 100644
--- a/src/sleplet/functions/axisymmetric_wavelet_coefficients_africa.py
+++ b/src/sleplet/functions/axisymmetric_wavelet_coefficients_africa.py
@@ -5,7 +5,7 @@
 import numpy.typing as npt
 import pydantic.v1 as pydantic
 
-import s2wav
+import pys2let
 
 import sleplet._string_methods
 import sleplet._validation
@@ -81,10 +81,7 @@ def _create_wavelet_coefficients(
 
     @pydantic.validator("j")
     def _check_j(cls, v, values):  # noqa: N805
-        j_max = s2wav.utils.shapes.j_max(
-            values["L"],
-            values["B"],
-        )
+        j_max = pys2let.pys2let_j_max(values["B"], values["L"], values["j_min"])
         if v is not None and v < 0:
             raise ValueError("j should be positive")
         if v is not None and v > j_max - values["j_min"]:
diff --git a/src/sleplet/functions/axisymmetric_wavelet_coefficients_earth.py b/src/sleplet/functions/axisymmetric_wavelet_coefficients_earth.py
index f3657e9cb..f16d8fb04 100644
--- a/src/sleplet/functions/axisymmetric_wavelet_coefficients_earth.py
+++ b/src/sleplet/functions/axisymmetric_wavelet_coefficients_earth.py
@@ -5,7 +5,7 @@
 import numpy.typing as npt
 import pydantic.v1 as pydantic
 
-import s2wav
+import pys2let
 
 import sleplet._string_methods
 import sleplet._validation
@@ -78,10 +78,7 @@ def _create_wavelet_coefficients(
 
     @pydantic.validator("j")
     def _check_j(cls, v, values):  # noqa: N805
-        j_max = s2wav.utils.shapes.j_max(
-            values["L"],
-            values["B"],
-        )
+        j_max = pys2let.pys2let_j_max(values["B"], values["L"], values["j_min"])
         if v is not None and v < 0:
             raise ValueError("j should be positive")
         if v is not None and v > j_max - values["j_min"]:
diff --git a/src/sleplet/functions/axisymmetric_wavelet_coefficients_south_america.py b/src/sleplet/functions/axisymmetric_wavelet_coefficients_south_america.py
index 51bd11697..6b2e26710 100644
--- a/src/sleplet/functions/axisymmetric_wavelet_coefficients_south_america.py
+++ b/src/sleplet/functions/axisymmetric_wavelet_coefficients_south_america.py
@@ -5,7 +5,7 @@
 import numpy.typing as npt
 import pydantic.v1 as pydantic
 
-import s2wav
+import pys2let
 
 import sleplet._string_methods
 import sleplet._validation
@@ -80,10 +80,7 @@ def _create_wavelet_coefficients(
 
     @pydantic.validator("j")
     def _check_j(cls, v, values):  # noqa: N805
-        j_max = s2wav.utils.shapes.j_max(
-            values["L"],
-            values["B"],
-        )
+        j_max = pys2let.pys2let_j_max(values["B"], values["L"], values["j_min"])
         if v is not None and v < 0:
             raise ValueError("j should be positive")
         if v is not None and v > j_max - values["j_min"]:
diff --git a/src/sleplet/functions/axisymmetric_wavelets.py b/src/sleplet/functions/axisymmetric_wavelets.py
index 394462427..4009c2e8f 100644
--- a/src/sleplet/functions/axisymmetric_wavelets.py
+++ b/src/sleplet/functions/axisymmetric_wavelets.py
@@ -5,7 +5,7 @@
 import numpy.typing as npt
 import pydantic.v1 as pydantic
 
-import s2wav
+import pys2let
 
 import sleplet._string_methods
 import sleplet._validation
@@ -71,10 +71,7 @@ def _create_wavelets(self) -> npt.NDArray[np.complex_]:
 
     @pydantic.validator("j")
     def _check_j(cls, v, values):  # noqa: N805
-        j_max = s2wav.utils.shapes.j_max(
-            values["L"],
-            values["B"],
-        )
+        j_max = pys2let.pys2let_j_max(values["B"], values["L"], values["j_min"])
         if v is not None and v < 0:
             raise ValueError("j should be positive")
         if v is not None and v > j_max - values["j_min"]:
diff --git a/src/sleplet/functions/dirac_delta.py b/src/sleplet/functions/dirac_delta.py
index 9b45f05dc..e26e12064 100644
--- a/src/sleplet/functions/dirac_delta.py
+++ b/src/sleplet/functions/dirac_delta.py
@@ -3,7 +3,7 @@
 import numpy.typing as npt
 import pydantic.v1 as pydantic
 
-import s2fft
+import pyssht as ssht
 
 import sleplet._string_methods
 import sleplet._validation
@@ -18,9 +18,10 @@ def __post_init_post_parse__(self) -> None:
         super().__post_init_post_parse__()
 
     def _create_coefficients(self) -> npt.NDArray[np.complex_ | np.float_]:
-        flm = np.zeros(s2fft.samples.flm_shape(self.L), dtype=np.complex_)
+        flm = np.zeros(self.L**2, dtype=np.complex_)
         for ell in range(self.L):
-            flm[ell, self.L - 1] = np.sqrt((2 * ell + 1) / (4 * np.pi))
+            ind = ssht.elm2ind(ell, 0)
+            flm[ind] = np.sqrt((2 * ell + 1) / (4 * np.pi))
         return flm
 
     def _create_name(self) -> str:
diff --git a/src/sleplet/functions/flm.py b/src/sleplet/functions/flm.py
index d8ac4c178..3a3966fdf 100644
--- a/src/sleplet/functions/flm.py
+++ b/src/sleplet/functions/flm.py
@@ -6,7 +6,6 @@
 import pydantic.v1 as pydantic
 
 import pyssht as ssht
-import s2fft
 
 import sleplet._validation
 import sleplet.noise
@@ -27,26 +26,14 @@ def rotate(  # noqa: D102
         *,
         gamma: float = 0,
     ) -> npt.NDArray[np.complex_]:
-        return s2fft.samples.flm_1d_to_2d(
-            ssht.rotate_flms(
-                s2fft.samples.flm_2d_to_1d(self.coefficients, self.L),
-                alpha,
-                beta,
-                gamma,
-                self.L,
-            ),
-            self.L,
-        )
+        return ssht.rotate_flms(self.coefficients, alpha, beta, gamma, self.L)
 
     def _translation_helper(
         self,
         alpha: float,
         beta: float,
     ) -> npt.NDArray[np.complex_]:
-        return s2fft.samples.flm_1d_to_2d(
-            ssht.create_ylm(beta, alpha, self.L).conj().flatten(),
-            self.L,
-        )
+        return ssht.create_ylm(beta, alpha, self.L).conj().flatten()
 
     def _add_noise_to_signal(
         self,
diff --git a/src/sleplet/functions/gaussian.py b/src/sleplet/functions/gaussian.py
index f586a546b..3b04c4f7f 100644
--- a/src/sleplet/functions/gaussian.py
+++ b/src/sleplet/functions/gaussian.py
@@ -3,7 +3,7 @@
 import numpy.typing as npt
 import pydantic.v1 as pydantic
 
-import s2fft
+import pyssht as ssht
 
 import sleplet._string_methods
 import sleplet._validation
@@ -21,9 +21,10 @@ def __post_init_post_parse__(self) -> None:
         super().__post_init_post_parse__()
 
     def _create_coefficients(self) -> npt.NDArray[np.complex_ | np.float_]:
-        flm = np.zeros(s2fft.samples.flm_shape(self.L), dtype=np.complex_)
+        flm = np.zeros(self.L**2, dtype=np.complex_)
         for ell in range(self.L):
-            flm[ell, self.L - 1] = np.exp(-ell * (ell + 1) / (2 * self.sigma**2))
+            ind = ssht.elm2ind(ell, 0)
+            flm[ind] = np.exp(-ell * (ell + 1) / (2 * self.sigma**2))
         return flm
 
     def _create_name(self) -> str:
diff --git a/src/sleplet/functions/harmonic_gaussian.py b/src/sleplet/functions/harmonic_gaussian.py
index b36dcbfe1..e91fbc42a 100644
--- a/src/sleplet/functions/harmonic_gaussian.py
+++ b/src/sleplet/functions/harmonic_gaussian.py
@@ -3,7 +3,7 @@
 import numpy.typing as npt
 import pydantic.v1 as pydantic
 
-import s2fft
+import pyssht as ssht
 
 import sleplet._string_methods
 import sleplet._validation
@@ -27,14 +27,12 @@ def __post_init_post_parse__(self) -> None:
         super().__post_init_post_parse__()
 
     def _create_coefficients(self) -> npt.NDArray[np.complex_ | np.float_]:
-        flm = np.zeros(s2fft.samples.flm_shape(self.L), dtype=np.complex_)
+        flm = np.zeros(self.L**2, dtype=np.complex_)
         for ell in range(self.L):
             upsilon_l = np.exp(-((ell / self.l_sigma) ** 2) / 2)
             for m in range(-ell, ell + 1):
-                ind_pm = self.L - 1 + m
-                flm[ell, ind_pm] = upsilon_l * np.exp(
-                    -((m / self.m_sigma) ** 2) / 2,
-                )
+                ind = ssht.elm2ind(ell, m)
+                flm[ind] = upsilon_l * np.exp(-((m / self.m_sigma) ** 2) / 2)
         return flm
 
     def _create_name(self) -> str:
diff --git a/src/sleplet/functions/identity.py b/src/sleplet/functions/identity.py
index f89a42862..2fe60fca0 100644
--- a/src/sleplet/functions/identity.py
+++ b/src/sleplet/functions/identity.py
@@ -3,8 +3,6 @@
 import numpy.typing as npt
 import pydantic.v1 as pydantic
 
-import s2fft
-
 import sleplet._string_methods
 import sleplet._validation
 from sleplet.functions.flm import Flm
@@ -18,7 +16,7 @@ def __post_init_post_parse__(self) -> None:
         super().__post_init_post_parse__()
 
     def _create_coefficients(self) -> npt.NDArray[np.complex_ | np.float_]:
-        return np.ones(s2fft.samples.flm_shape(self.L), dtype=np.complex_)
+        return np.ones(self.L**2, dtype=np.complex_)
 
     def _create_name(self) -> str:
         return sleplet._string_methods._convert_camel_case_to_snake_case(
diff --git a/src/sleplet/functions/ridgelets.py b/src/sleplet/functions/ridgelets.py
index bb463bf2f..855283297 100644
--- a/src/sleplet/functions/ridgelets.py
+++ b/src/sleplet/functions/ridgelets.py
@@ -6,8 +6,8 @@
 import pydantic.v1 as pydantic
 import scipy.special
 
-import s2fft
-import s2wav
+import pys2let
+import pyssht as ssht
 
 import sleplet._string_methods
 import sleplet._validation
@@ -70,21 +70,16 @@ def _create_wavelets(self) -> npt.NDArray[np.complex_]:
         """Compute all wavelets."""
         ring_lm = self._compute_ring()
         kappas = sleplet.wavelet_methods.create_kappas(self.L, self.B, self.j_min)
-        wavelets = np.zeros(
-            (kappas.shape[0], *s2fft.samples.flm_shape(self.L)),
-            dtype=np.complex_,
-        )
-        ind_m0 = self.L - 1
+        wavelets = np.zeros((kappas.shape[0], self.L**2), dtype=np.complex_)
         for ell in range(self.L):
-            wavelets[0, ell, ind_m0] = kappas[0, ell] * ring_lm[ell, ind_m0]
-            wavelets[1:, ell, ind_m0] = (
-                kappas[1:, ell] * ring_lm[ell, ind_m0] / np.sqrt(2 * np.pi)
-            )
+            ind = ssht.elm2ind(ell, 0)
+            wavelets[0, ind] = kappas[0, ell] * ring_lm[ind]
+            wavelets[1:, ind] = kappas[1:, ell] * ring_lm[ind] / np.sqrt(2 * np.pi)
         return wavelets
 
     def _compute_ring(self) -> npt.NDArray[np.complex_]:
         """Compute ring in harmonic space."""
-        ring_lm = np.zeros(s2fft.samples.flm_shape(self.L), dtype=np.complex_)
+        ring_lm = np.zeros(self.L**2, dtype=np.complex_)
         for ell in range(abs(self.spin), self.L):
             logp2 = (
                 scipy.special.gammaln(ell + self.spin + 1)
@@ -93,7 +88,8 @@ def _compute_ring(self) -> npt.NDArray[np.complex_]:
                 - scipy.special.gammaln((ell - self.spin) / 2 + 1)
             )
             p0 = np.real((-1) ** ((ell + self.spin) / 2)) * np.exp(logp2)
-            ring_lm[ell, self.L - 1] = (
+            ind = ssht.elm2ind(ell, 0)
+            ring_lm[ind] = (
                 2
                 * np.pi
                 * np.sqrt((2 * ell + 1) / (4 * np.pi))
@@ -110,10 +106,7 @@ def _compute_ring(self) -> npt.NDArray[np.complex_]:
 
     @pydantic.validator("j")
     def _check_j(cls, v, values):  # noqa: N805
-        j_max = s2wav.utils.shapes.j_max(
-            values["L"],
-            values["B"],
-        )
+        j_max = pys2let.pys2let_j_max(values["B"], values["L"], values["j_min"])
         if v is not None and v < 0:
             raise ValueError("j should be positive")
         if v is not None and v > j_max - values["j_min"]:
diff --git a/src/sleplet/functions/slepian_dirac_delta.py b/src/sleplet/functions/slepian_dirac_delta.py
index 3a56f69d5..9a5bc7eb5 100644
--- a/src/sleplet/functions/slepian_dirac_delta.py
+++ b/src/sleplet/functions/slepian_dirac_delta.py
@@ -5,7 +5,7 @@
 import numpy.typing as npt
 import pydantic.v1 as pydantic
 
-import s2fft
+import pyssht as ssht
 
 import sleplet._string_methods
 import sleplet._validation
@@ -52,25 +52,15 @@ def _setup_args(self) -> None:
 
     def _compute_angles(self) -> None:
         """Computes alpha/beta if not provided."""
-        thetas = np.tile(
-            s2fft.samples.thetas(self.L, sampling=sleplet._vars.SAMPLING_SCHEME),
-            (
-                s2fft.samples.nphi_equiang(
-                    self.L,
-                    sampling=sleplet._vars.SAMPLING_SCHEME,
-                ),
-                1,
-            ),
-        ).T
-        phis = np.tile(
-            s2fft.samples.phis_equiang(self.L, sampling=sleplet._vars.SAMPLING_SCHEME),
-            (s2fft.samples.ntheta(self.L, sampling=sleplet._vars.SAMPLING_SCHEME), 1),
+        thetas, phis = ssht.sample_positions(
+            self.L,
+            Grid=True,
+            Method=sleplet._vars.SAMPLING_SCHEME,
         )
-        sp = s2fft.inverse(
+        sp = ssht.inverse(
             self.slepian.eigenvectors[0],
             self.L,
-            method=sleplet._vars.EXECUTION_MODE,
-            sampling=sleplet._vars.SAMPLING_SCHEME,
+            Method=sleplet._vars.SAMPLING_SCHEME,
         )
         idx = tuple(np.argwhere(sp == sp.max())[0])
         self._alpha = phis[idx]
diff --git a/src/sleplet/functions/slepian_wavelet_coefficients_africa.py b/src/sleplet/functions/slepian_wavelet_coefficients_africa.py
index 1b447fc23..85be5285a 100644
--- a/src/sleplet/functions/slepian_wavelet_coefficients_africa.py
+++ b/src/sleplet/functions/slepian_wavelet_coefficients_africa.py
@@ -5,7 +5,7 @@
 import numpy.typing as npt
 import pydantic.v1 as pydantic
 
-import s2wav
+import pys2let
 
 import sleplet._string_methods
 import sleplet._validation
@@ -91,10 +91,7 @@ def _create_wavelet_coefficients(
 
     @pydantic.validator("j")
     def _check_j(cls, v, values):  # noqa: N805
-        j_max = s2wav.utils.shapes.j_max(
-            values["L"] ** 2,
-            values["B"],
-        )
+        j_max = pys2let.pys2let_j_max(values["B"], values["L"] ** 2, values["j_min"])
         if v is not None and v < 0:
             raise ValueError("j should be positive")
         if v is not None and v > j_max - values["j_min"]:
diff --git a/src/sleplet/functions/slepian_wavelet_coefficients_south_america.py b/src/sleplet/functions/slepian_wavelet_coefficients_south_america.py
index d7358033b..ddcda1b39 100644
--- a/src/sleplet/functions/slepian_wavelet_coefficients_south_america.py
+++ b/src/sleplet/functions/slepian_wavelet_coefficients_south_america.py
@@ -5,7 +5,7 @@
 import numpy.typing as npt
 import pydantic.v1 as pydantic
 
-import s2wav
+import pys2let
 
 import sleplet._string_methods
 import sleplet._validation
@@ -91,10 +91,7 @@ def _create_wavelet_coefficients(
 
     @pydantic.validator("j")
     def _check_j(cls, v, values):  # noqa: N805
-        j_max = s2wav.utils.shapes.j_max(
-            values["L"] ** 2,
-            values["B"],
-        )
+        j_max = pys2let.pys2let_j_max(values["B"], values["L"] ** 2, values["j_min"])
         if v is not None and v < 0:
             raise ValueError("j should be positive")
         if v is not None and v > j_max - values["j_min"]:
diff --git a/src/sleplet/functions/slepian_wavelets.py b/src/sleplet/functions/slepian_wavelets.py
index 2a09712b4..7f4e64dcc 100644
--- a/src/sleplet/functions/slepian_wavelets.py
+++ b/src/sleplet/functions/slepian_wavelets.py
@@ -5,7 +5,7 @@
 import numpy.typing as npt
 import pydantic.v1 as pydantic
 
-import s2wav
+import pys2let
 
 import sleplet._string_methods
 import sleplet._validation
@@ -65,10 +65,7 @@ def _create_wavelets(self) -> npt.NDArray[np.float_]:
 
     @pydantic.validator("j")
     def _check_j(cls, v, values):  # noqa: N805
-        j_max = s2wav.utils.shapes.j_max(
-            values["L"] ** 2,
-            values["B"],
-        )
+        j_max = pys2let.pys2let_j_max(values["B"], values["L"] ** 2, values["j_min"])
         if v is not None and v < 0:
             raise ValueError("j should be positive")
         if v is not None and v > j_max - values["j_min"]:
diff --git a/src/sleplet/functions/south_america.py b/src/sleplet/functions/south_america.py
index c0af2b7d4..26b10f516 100644
--- a/src/sleplet/functions/south_america.py
+++ b/src/sleplet/functions/south_america.py
@@ -3,7 +3,7 @@
 import numpy.typing as npt
 import pydantic.v1 as pydantic
 
-import s2fft
+import pyssht as ssht
 
 import sleplet._data.create_earth_flm
 import sleplet._data.setup_pooch
@@ -61,12 +61,11 @@ def _grid_fun(
             earth_flm,
             self.L,
         )
-        earth_f = s2fft.inverse(
+        earth_f = ssht.inverse(
             rot_flm,
             self.L,
-            method=sleplet._vars.EXECUTION_MODE,
-            reality=self.reality,
-            sampling=sleplet._vars.SAMPLING_SCHEME,
+            Method=sleplet._vars.SAMPLING_SCHEME,
+            Reality=self.reality,
         )
         mask_name = f"{self.name}_L{self.L}.npy"
         mask_location = sleplet._data.setup_pooch.find_on_pooch_then_local(
diff --git a/src/sleplet/functions/spherical_harmonic.py b/src/sleplet/functions/spherical_harmonic.py
index c8569f7ca..5f6a72251 100644
--- a/src/sleplet/functions/spherical_harmonic.py
+++ b/src/sleplet/functions/spherical_harmonic.py
@@ -3,6 +3,8 @@
 import numpy.typing as npt
 import pydantic.v1 as pydantic
 
+import pyssht as ssht
+
 import sleplet._string_methods
 import sleplet._validation
 import sleplet.harmonic_methods
@@ -19,11 +21,8 @@ class SphericalHarmonic(Flm):
     r"""Order \(\leq |\ell|\)"""
 
     def _create_coefficients(self) -> npt.NDArray[np.complex_ | np.float_]:
-        return sleplet.harmonic_methods._create_spherical_harmonic(
-            self.L,
-            self.ell,
-            self.m,
-        )
+        ind = ssht.elm2ind(self.ell, self.m)
+        return sleplet.harmonic_methods._create_spherical_harmonic(self.L, ind)
 
     def _create_name(self) -> str:
         return (
diff --git a/src/sleplet/harmonic_methods.py b/src/sleplet/harmonic_methods.py
index 2ade4cfe0..06cf9afd4 100644
--- a/src/sleplet/harmonic_methods.py
+++ b/src/sleplet/harmonic_methods.py
@@ -6,7 +6,6 @@
 import numpy.typing as npt
 
 import pyssht as ssht
-import s2fft
 
 import sleplet._data.create_earth_flm
 import sleplet._integration_methods
@@ -21,10 +20,10 @@
 _SOUTH_AMERICA_GAMMA = np.deg2rad(63)
 
 
-def _create_spherical_harmonic(L: int, ell: int, m: int) -> npt.NDArray[np.complex_]:
+def _create_spherical_harmonic(L: int, ind: int) -> npt.NDArray[np.complex_]:
     """Create a spherical harmonic in harmonic space for the given index."""
-    flm = np.zeros(s2fft.samples.flm_shape(L), dtype=np.complex_)
-    flm[ell, L - 1 + m] = 1
+    flm = np.zeros(L**2, dtype=np.complex_)
+    flm[ind] = 1
     return flm
 
 
@@ -58,17 +57,13 @@ def invert_flm_boosted(
         The boosted field value.
     """
     boost = resolution**2 - L**2
-    flm = s2fft.samples.flm_1d_to_2d(
-        _boost_coefficient_resolution(s2fft.samples.flm_2d_to_1d(flm, L), boost),
-        resolution,
-    )
-    return s2fft.inverse(
+    flm = _boost_coefficient_resolution(flm, boost)
+    return ssht.inverse(
         flm,
         resolution,
-        method=sleplet._vars.EXECUTION_MODE,
-        reality=reality,
-        sampling=sleplet._vars.SAMPLING_SCHEME,
-        spin=spin,
+        Method=sleplet._vars.SAMPLING_SCHEME,
+        Reality=reality,
+        Spin=spin,
     )
 
 
@@ -86,22 +81,18 @@ def _ensure_f_bandlimited(
     If the function created is created in pixel space rather than harmonic
     space then need to transform it into harmonic space first before using it.
     """
-    thetas = np.tile(
-        s2fft.samples.thetas(L, sampling=sleplet._vars.SAMPLING_SCHEME),
-        (s2fft.samples.nphi_equiang(L, sampling=sleplet._vars.SAMPLING_SCHEME), 1),
-    ).T
-    phis = np.tile(
-        s2fft.samples.phis_equiang(L, sampling=sleplet._vars.SAMPLING_SCHEME),
-        (s2fft.samples.ntheta(L, sampling=sleplet._vars.SAMPLING_SCHEME), 1),
+    thetas, phis = ssht.sample_positions(
+        L,
+        Grid=True,
+        Method=sleplet._vars.SAMPLING_SCHEME,
     )
     f = grid_fun(thetas, phis)
-    return s2fft.forward(
+    return ssht.forward(
         f,
         L,
-        method=sleplet._vars.EXECUTION_MODE,
-        reality=reality,
-        sampling=sleplet._vars.SAMPLING_SCHEME,
-        spin=spin,
+        Method=sleplet._vars.SAMPLING_SCHEME,
+        Reality=reality,
+        Spin=spin,
     )
 
 
@@ -117,6 +108,29 @@ def _create_emm_vector(L: int) -> npt.NDArray[np.float_]:
     return emm
 
 
+def compute_random_signal(
+    L: int,
+    rng: np.random.Generator,
+    *,
+    var_signal: float,
+) -> npt.NDArray[np.complex_]:
+    """
+    Generates a normally distributed random signal of a
+    complex signal with mean `0` and variance `1`.
+
+    Args:
+        L: The spherical harmonic bandlimit.
+        rng: The random number generator object.
+        var_signal: The variance of the signal.
+
+    Returns:
+        The coefficients of a random signal on the sphere.
+    """
+    return np.sqrt(var_signal / 2) * (
+        rng.standard_normal(L**2) + 1j * rng.standard_normal(L**2)
+    )
+
+
 def mesh_forward(
     mesh: "sleplet.meshes.mesh.Mesh",
     u: npt.NDArray[np.complex_ | np.float_],
@@ -173,14 +187,11 @@ def rotate_earth_to_south_america(
     Returns:
         The spherical harmonic coefficients of the Earth centered on South America.
     """
-    return s2fft.samples.flm_1d_to_2d(
-        ssht.rotate_flms(
-            s2fft.samples.flm_2d_to_1d(earth_flm, L),
-            _SOUTH_AMERICA_ALPHA,
-            _SOUTH_AMERICA_BETA,
-            _SOUTH_AMERICA_GAMMA,
-            L,
-        ),
+    return ssht.rotate_flms(
+        earth_flm,
+        _SOUTH_AMERICA_ALPHA,
+        _SOUTH_AMERICA_BETA,
+        _SOUTH_AMERICA_GAMMA,
         L,
     )
 
@@ -199,13 +210,4 @@ def rotate_earth_to_africa(
     Returns:
         The spherical harmonic coefficients of the Earth centered on Africa.
     """
-    return s2fft.samples.flm_1d_to_2d(
-        ssht.rotate_flms(
-            s2fft.samples.flm_2d_to_1d(earth_flm, L),
-            _AFRICA_ALPHA,
-            _AFRICA_BETA,
-            _AFRICA_GAMMA,
-            L,
-        ),
-        L,
-    )
+    return ssht.rotate_flms(earth_flm, _AFRICA_ALPHA, _AFRICA_BETA, _AFRICA_GAMMA, L)
diff --git a/src/sleplet/meshes/mesh_slepian_wavelet_coefficients.py b/src/sleplet/meshes/mesh_slepian_wavelet_coefficients.py
index 56acac9fc..89053d18f 100644
--- a/src/sleplet/meshes/mesh_slepian_wavelet_coefficients.py
+++ b/src/sleplet/meshes/mesh_slepian_wavelet_coefficients.py
@@ -5,7 +5,7 @@
 import numpy.typing as npt
 import pydantic.v1 as pydantic
 
-import s2wav
+import pys2let
 
 import sleplet._string_methods
 import sleplet._validation
@@ -76,9 +76,10 @@ def _create_wavelet_coefficients(
 
     @pydantic.validator("j")
     def _check_j(cls, v, values):  # noqa: N805
-        j_max = s2wav.utils.shapes.j_max(
-            values["mesh"].mesh_eigenvalues.shape[0],
+        j_max = pys2let.pys2let_j_max(
             values["B"],
+            values["mesh"].mesh_eigenvalues.shape[0],
+            values["j_min"],
         )
         if v is not None and v < 0:
             raise ValueError("j should be positive")
diff --git a/src/sleplet/meshes/mesh_slepian_wavelets.py b/src/sleplet/meshes/mesh_slepian_wavelets.py
index 90be2f3d5..7ec36d52a 100644
--- a/src/sleplet/meshes/mesh_slepian_wavelets.py
+++ b/src/sleplet/meshes/mesh_slepian_wavelets.py
@@ -5,7 +5,7 @@
 import numpy.typing as npt
 import pydantic.v1 as pydantic
 
-import s2wav
+import pys2let
 
 import sleplet._string_methods
 import sleplet._validation
@@ -62,9 +62,10 @@ def _create_wavelets(self) -> npt.NDArray[np.float_]:
 
     @pydantic.validator("j")
     def _check_j(cls, v, values) -> int | None:  # noqa: N805
-        j_max = s2wav.utils.shapes.j_max(
-            values["mesh"].mesh_eigenvalues.shape[0],
+        j_max = pys2let.pys2let_j_max(
             values["B"],
+            values["mesh"].mesh_eigenvalues.shape[0],
+            values["j_min"],
         )
         if v is not None and v < 0:
             raise ValueError("j should be positive")
diff --git a/src/sleplet/noise.py b/src/sleplet/noise.py
index 64006dc83..7e3fbc9ae 100644
--- a/src/sleplet/noise.py
+++ b/src/sleplet/noise.py
@@ -4,7 +4,7 @@
 import numpy as np
 import numpy.typing as npt
 
-import s2fft
+import pyssht as ssht
 
 import sleplet._vars
 import sleplet.harmonic_methods
@@ -73,23 +73,24 @@ def _create_noise(
     rng = np.random.default_rng(sleplet._vars.RANDOM_SEED)
 
     # initialise
-    nlm = np.zeros(s2fft.samples.flm_shape(L), dtype=np.complex_)
+    nlm = np.zeros(L**2, dtype=np.complex_)
 
     # std dev of the noise
     sigma_noise = compute_sigma_noise(signal, snr_in)
 
     # compute noise
     for ell in range(L):
-        nlm[ell, L - 1] = sigma_noise * rng.standard_normal()
+        ind = ssht.elm2ind(ell, 0)
+        nlm[ind] = sigma_noise * rng.standard_normal()
         for m in range(1, ell + 1):
-            ind_pm = L - 1 + m
-            ind_nm = L - 1 - m
-            nlm[ell, ind_pm] = (
+            ind_pm = ssht.elm2ind(ell, m)
+            ind_nm = ssht.elm2ind(ell, -m)
+            nlm[ind_pm] = (
                 sigma_noise
                 / np.sqrt(2)
                 * (rng.standard_normal() + 1j * rng.standard_normal())
             )
-            nlm[ell, ind_nm] = (-1) ** m * nlm[ell, ind_pm].conj()
+            nlm[ind_nm] = (-1) ** m * nlm[ind_pm].conj()
     return nlm
 
 
@@ -100,11 +101,10 @@ def _create_slepian_noise(
     snr_in: float,
 ) -> npt.NDArray[np.complex_]:
     """Computes Gaussian white noise in Slepian space."""
-    flm = s2fft.forward(
+    flm = ssht.forward(
         sleplet.slepian_methods.slepian_inverse(slepian_signal, L, slepian),
         L,
-        method=sleplet._vars.EXECUTION_MODE,
-        sampling=sleplet._vars.SAMPLING_SCHEME,
+        Method=sleplet._vars.SAMPLING_SCHEME,
     )
     nlm = _create_noise(L, flm, snr_in)
     return sleplet.slepian_methods.slepian_forward(L, slepian, flm=nlm)
@@ -141,18 +141,12 @@ def harmonic_hard_thresholding(
     _logger.info("begin harmonic hard thresholding")
     for j, coefficient in enumerate(wav_coeffs[1:]):
         _logger.info(f"start Psi^{j + 1}/{len(wav_coeffs)-1}")
-        f = s2fft.inverse(
-            coefficient,
-            L,
-            method=sleplet._vars.EXECUTION_MODE,
-            sampling=sleplet._vars.SAMPLING_SCHEME,
-        )
+        f = ssht.inverse(coefficient, L, Method=sleplet._vars.SAMPLING_SCHEME)
         f_thresholded = _perform_hard_thresholding(f, sigma_j[j], n_sigma)
-        wav_coeffs[j + 1] = s2fft.forward(
+        wav_coeffs[j + 1] = ssht.forward(
             f_thresholded,
             L,
-            method=sleplet._vars.EXECUTION_MODE,
-            sampling=sleplet._vars.SAMPLING_SCHEME,
+            Method=sleplet._vars.SAMPLING_SCHEME,
         )
     return wav_coeffs
 
@@ -223,7 +217,7 @@ def _compute_sigma_j(
 ) -> npt.NDArray[np.float_]:
     """Compute sigma_j for wavelets used in denoising the signal."""
     sigma_noise = compute_sigma_noise(signal, snr_in)
-    wavelet_power = (np.abs(psi_j) ** 2).sum(axis=(1, 2))
+    wavelet_power = (np.abs(psi_j) ** 2).sum(axis=1)
     return sigma_noise * np.sqrt(wavelet_power)
 
 
diff --git a/src/sleplet/plot_methods.py b/src/sleplet/plot_methods.py
index 1dd8dba35..0045cb2ae 100644
--- a/src/sleplet/plot_methods.py
+++ b/src/sleplet/plot_methods.py
@@ -5,7 +5,7 @@
 import numpy as np
 import numpy.typing as npt
 
-import s2fft
+import pyssht as ssht
 
 import sleplet._mask_methods
 import sleplet._vars
@@ -67,8 +67,7 @@ def _calc_nearest_grid_point(
     values - the translation needs to be at the same position
     as the rotation such that the difference error is small.
     """
-    thetas = s2fft.samples.thetas(L, sampling=sleplet._vars.SAMPLING_SCHEME)
-    phis = s2fft.samples.phis_equiang(L, sampling=sleplet._vars.SAMPLING_SCHEME)
+    thetas, phis = ssht.sample_positions(L, Method=sleplet._vars.SAMPLING_SCHEME)
     pix_j = np.abs(phis - alpha_pi_fraction * np.pi).argmin()
     pix_i = np.abs(thetas - beta_pi_fraction * np.pi).argmin()
     alpha, beta = phis[pix_j], thetas[pix_i]
@@ -102,11 +101,10 @@ def find_max_amplitude(
             function.slepian,
         )
     else:
-        field = s2fft.inverse(
+        field = ssht.inverse(
             function.coefficients,
             function.L,
-            method=sleplet._vars.EXECUTION_MODE,
-            sampling=sleplet._vars.SAMPLING_SCHEME,
+            Method=sleplet._vars.SAMPLING_SCHEME,
         )
 
     # find resolution of final plot for boosting if necessary
@@ -154,7 +152,7 @@ def _set_outside_region_to_minimum(
     mask = sleplet._mask_methods.create_mask_region(L, region)
 
     # adapt for closed plot
-    n_phi = s2fft.samples.nphi_equiang(L, sampling=sleplet._vars.SAMPLING_SCHEME)
+    _, n_phi = ssht.sample_shape(L, Method=sleplet._vars.SAMPLING_SCHEME)
     closed_mask = np.insert(mask, n_phi, mask[:, 0], axis=1)
 
     # set values outside mask to negative infinity
@@ -188,13 +186,12 @@ def _boost_field(  # noqa: PLR0913
     """Inverts and then boosts the field before plotting."""
     if not upsample:
         return field
-    flm = s2fft.forward(
+    flm = ssht.forward(
         field,
         L,
-        method=sleplet._vars.EXECUTION_MODE,
-        reality=reality,
-        sampling=sleplet._vars.SAMPLING_SCHEME,
-        spin=spin,
+        Method=sleplet._vars.SAMPLING_SCHEME,
+        Reality=reality,
+        Spin=spin,
     )
     return sleplet.harmonic_methods.invert_flm_boosted(
         flm,
@@ -263,12 +260,11 @@ def _coefficients_to_field_sphere(
     return (
         sleplet.slepian_methods.slepian_inverse(coefficients, f.L, f.slepian)
         if hasattr(f, "slepian")
-        else s2fft.inverse(
+        else ssht.inverse(
             coefficients,
             f.L,
-            method=sleplet._vars.EXECUTION_MODE,
-            reality=f.reality,
-            sampling=sleplet._vars.SAMPLING_SCHEME,
-            spin=f.spin,
+            Method=sleplet._vars.SAMPLING_SCHEME,
+            Reality=f.reality,
+            Spin=f.spin,
         )
     )
diff --git a/src/sleplet/plotting/_create_plot_sphere.py b/src/sleplet/plotting/_create_plot_sphere.py
index 5e5fb890a..578448a9e 100644
--- a/src/sleplet/plotting/_create_plot_sphere.py
+++ b/src/sleplet/plotting/_create_plot_sphere.py
@@ -2,7 +2,6 @@
 import logging
 
 import cmocean
-import jax
 import numpy as np
 import numpy.typing as npt
 import plotly.graph_objs as go
@@ -10,7 +9,6 @@
 import pydantic.v1 as pydantic
 
 import pyssht as ssht
-import s2fft
 
 import sleplet._plotly_methods
 import sleplet._validation
@@ -25,7 +23,7 @@
 class PlotSphere:
     """Creates surface sphere plot via `plotly`."""
 
-    f: jax.Array | npt.NDArray[np.complex_ | np.float_]
+    f: npt.NDArray[np.complex_ | np.float_]
     """The field value sampled on the sphere."""
     L: int
     """The spherical harmonic bandlimit."""
@@ -84,12 +82,10 @@ def execute(self) -> None:
         camera = sleplet._plotly_methods.create_camera(-0.1, -0.1, 10, 7.88)
 
         # pick largest tick max value
-        tick_mark = float(
-            sleplet._plotly_methods.create_tick_mark(
-                vmin,
-                vmax,
-                amplitude=self.amplitude,
-            ),
+        tick_mark = sleplet._plotly_methods.create_tick_mark(
+            vmin,
+            vmax,
+            amplitude=self.amplitude,
         )
 
         data = [
@@ -126,7 +122,7 @@ def _setup_plot(  # noqa: PLR0913
         f: npt.NDArray[np.float_],
         resolution: int,
         *,
-        method: str = "mw",
+        method: str = "MW",
         close: bool = True,
         parametric: bool = False,
         parametric_scaling: list[float] | None = None,
@@ -143,14 +139,7 @@ def _setup_plot(  # noqa: PLR0913
         if parametric_scaling is None:
             parametric_scaling = [0.0, 0.5]
 
-        thetas = np.tile(
-            s2fft.samples.thetas(resolution, sampling=method),
-            (s2fft.samples.nphi_equiang(resolution, sampling=method), 1),
-        ).T
-        phis = np.tile(
-            s2fft.samples.phis_equiang(resolution, sampling=method),
-            (s2fft.samples.ntheta(resolution, sampling=method), 1),
-        )
+        thetas, phis = ssht.sample_positions(resolution, Grid=True, Method=method)
 
         if thetas.size != f.size:
             raise AttributeError("Bandlimit L deos not match that of f")
@@ -180,7 +169,7 @@ def _setup_plot(  # noqa: PLR0913
 
         # Close plot.
         if close:
-            n_phi = s2fft.samples.nphi_equiang(resolution, sampling=method)
+            _, n_phi = ssht.sample_shape(resolution, Method=method)
             f_plot = np.insert(f_plot, n_phi, f[:, 0], axis=1)
             if parametric:
                 f_normalised = np.insert(
diff --git a/src/sleplet/slepian/_slepian_decomposition.py b/src/sleplet/slepian/_slepian_decomposition.py
index 8196a2f70..223a8e72c 100644
--- a/src/sleplet/slepian/_slepian_decomposition.py
+++ b/src/sleplet/slepian/_slepian_decomposition.py
@@ -1,12 +1,11 @@
 import dataclasses
 import logging
 
-import jax
 import numpy as np
 import numpy.typing as npt
 import pydantic.v1 as pydantic
 
-import s2fft
+import pyssht as ssht
 
 import sleplet._integration_methods
 import sleplet._validation
@@ -21,8 +20,8 @@ class SlepianDecomposition:
     L: int
     slepian: SlepianFunctions
     _: dataclasses.KW_ONLY
-    f: jax.Array | npt.NDArray[np.complex_] | None = None
-    flm: jax.Array | npt.NDArray[np.complex_ | np.float_] | None = None
+    f: npt.NDArray[np.complex_] | None = None
+    flm: npt.NDArray[np.complex_ | np.float_] | None = None
     mask: npt.NDArray[np.float_] | None = None
 
     def __post_init_post_parse__(self) -> None:
@@ -56,11 +55,10 @@ def _integrate_region(self, rank: int) -> complex:
         \int\limits_{R} \dd{\Omega(\omega)}
         f(\omega) \overline{S_{p}(\omega)}.
         """
-        s_p = s2fft.inverse(
+        s_p = ssht.inverse(
             self.slepian.eigenvectors[rank],
             self.L,
-            method=sleplet._vars.EXECUTION_MODE,
-            sampling=sleplet._vars.SAMPLING_SCHEME,
+            Method=sleplet._vars.SAMPLING_SCHEME,
         )
         weight = sleplet._integration_methods.calc_integration_weight(self.L)
         integration = sleplet._integration_methods.integrate_region_sphere(
@@ -77,11 +75,10 @@ def _integrate_sphere(self, rank: int) -> complex:
         \int\limits_{S^{2}} \dd{\Omega(\omega)}
         f(\omega) \overline{S_{p}(\omega)}.
         """
-        s_p = s2fft.inverse(
+        s_p = ssht.inverse(
             self.slepian.eigenvectors[rank],
             self.L,
-            method=sleplet._vars.EXECUTION_MODE,
-            sampling=sleplet._vars.SAMPLING_SCHEME,
+            Method=sleplet._vars.SAMPLING_SCHEME,
         )
         weight = sleplet._integration_methods.calc_integration_weight(self.L)
         return sleplet._integration_methods.integrate_whole_sphere(
diff --git a/src/sleplet/slepian/slepian_arbitrary.py b/src/sleplet/slepian/slepian_arbitrary.py
index fa8a893ae..4b86c39b2 100644
--- a/src/sleplet/slepian/slepian_arbitrary.py
+++ b/src/sleplet/slepian/slepian_arbitrary.py
@@ -9,7 +9,7 @@
 import platformdirs
 import pydantic.v1 as pydantic
 
-import s2fft
+import pyssht as ssht
 
 import sleplet._array_methods
 import sleplet._data.setup_pooch
@@ -158,10 +158,10 @@ def _matrix_helper(
         integral = self._integral(i, i)
         D_r[i][i] = integral.real
         D_i[i][i] = integral.imag
-        _, m_i = s2fft.samples.ind2elm(i)
+        _, m_i = ssht.ind2elm(i)
 
         for j in range(i + 1, D_r.shape[0]):
-            ell_j, m_j = s2fft.samples.ind2elm(j)
+            ell_j, m_j = ssht.ind2elm(j)
             # if possible to use previous calculations
             if m_i == 0 and m_j != 0 and ell_j < self.L:
                 # if positive m then use conjugate relation
@@ -169,7 +169,7 @@ def _matrix_helper(
                     integral = self._integral(j, i)
                     D_r[j][i] = integral.real
                     D_i[j][i] = integral.imag
-                    k = s2fft.samples.elm2ind(ell_j, -m_j)
+                    k = ssht.elm2ind(ell_j, -m_j)
                     D_r[k][i] = (-1) ** m_j * D_r[j][i]
                     D_i[k][i] = (-1) ** (m_j + 1) * D_i[j][i]
             else:
@@ -180,16 +180,14 @@ def _matrix_helper(
     def _integral(self, i: int, j: int) -> complex:
         """Calculates the D integral between two spherical harmonics."""
         if i not in self._fields:
-            ell, m = s2fft.samples.ind2elm(i)
             self._fields[i] = sleplet.harmonic_methods.invert_flm_boosted(
-                sleplet.harmonic_methods._create_spherical_harmonic(self.L, ell, m),
+                sleplet.harmonic_methods._create_spherical_harmonic(self.L, i),
                 self.L,
                 self.resolution,
             )
         if j not in self._fields:
-            ell, m = s2fft.samples.ind2elm(j)
             self._fields[j] = sleplet.harmonic_methods.invert_flm_boosted(
-                sleplet.harmonic_methods._create_spherical_harmonic(self.L, ell, m),
+                sleplet.harmonic_methods._create_spherical_harmonic(self.L, j),
                 self.L,
                 self.resolution,
             )
diff --git a/src/sleplet/slepian/slepian_functions.py b/src/sleplet/slepian/slepian_functions.py
index 8f9c7f27d..bba5c6f29 100644
--- a/src/sleplet/slepian/slepian_functions.py
+++ b/src/sleplet/slepian/slepian_functions.py
@@ -6,8 +6,6 @@
 import numpy.typing as npt
 import pydantic.v1 as pydantic
 
-import s2fft
-
 import sleplet._validation
 
 _logger = logging.getLogger(__name__)
@@ -29,10 +27,7 @@ def __post_init_post_parse__(self) -> None:
         _logger.info(f"Shannon number N={self.N}")
         self.matrix_location = self._create_matrix_location()
         _logger.info("start solving eigenproblem")
-        self.eigenvalues, evecs = self._solve_eigenproblem()
-        self.eigenvectors = np.array(
-            [s2fft.samples.flm_1d_to_2d(e, self.L) for e in evecs],
-        )
+        self.eigenvalues, self.eigenvectors = self._solve_eigenproblem()
         _logger.info("finished solving eigenproblem")
 
     @abc.abstractmethod
diff --git a/src/sleplet/slepian_methods.py b/src/sleplet/slepian_methods.py
index 771069428..402c491e0 100644
--- a/src/sleplet/slepian_methods.py
+++ b/src/sleplet/slepian_methods.py
@@ -5,7 +5,6 @@
 import numpy.typing as npt
 
 import pyssht as ssht
-import s2fft
 
 import sleplet._vars
 import sleplet.harmonic_methods
@@ -137,16 +136,15 @@ def compute_s_p_omega(
     Returns:
         The complex \(S_{p}(\omega)\) values.
     """
-    n_theta, n_phi = s2fft.samples.f_shape(L, sampling=sleplet._vars.SAMPLING_SCHEME)
+    n_theta, n_phi = ssht.sample_shape(L, Method=sleplet._vars.SAMPLING_SCHEME)
     sp = np.zeros((slepian.N, n_theta, n_phi), dtype=np.complex_)
     for p in range(slepian.N):
         if p % L == 0:
             _logger.info(f"compute Sp(omega) p={p+1}/{slepian.N}")
-        sp[p] = s2fft.inverse(
+        sp[p] = ssht.inverse(
             slepian.eigenvectors[p],
             L,
-            method=sleplet._vars.EXECUTION_MODE,
-            sampling=sleplet._vars.SAMPLING_SCHEME,
+            Method=sleplet._vars.SAMPLING_SCHEME,
         )
     return sp
 
@@ -159,8 +157,8 @@ def _compute_s_p_omega_prime(
 ) -> npt.NDArray[np.complex_]:
     """Method to pick out the desired angle from Sp(omega)."""
     sp_omega = compute_s_p_omega(L, slepian)
-    p = ssht.theta_to_index(beta, L, Method=sleplet._vars.SAMPLING_SCHEME.upper())
-    q = ssht.phi_to_index(alpha, L, Method=sleplet._vars.SAMPLING_SCHEME.upper())
+    p = ssht.theta_to_index(beta, L, Method=sleplet._vars.SAMPLING_SCHEME)
+    q = ssht.phi_to_index(alpha, L, Method=sleplet._vars.SAMPLING_SCHEME)
     sp_omega_prime = sp_omega[:, p, q]
     # pad with zeros so it has the expected shape
     boost = L**2 - slepian.N
diff --git a/src/sleplet/wavelet_methods.py b/src/sleplet/wavelet_methods.py
index d8b7f0c0b..c16798484 100644
--- a/src/sleplet/wavelet_methods.py
+++ b/src/sleplet/wavelet_methods.py
@@ -2,8 +2,8 @@
 import numpy as np
 import numpy.typing as npt
 
-import s2fft
-import s2wav
+import pys2let
+import pyssht as ssht
 
 import sleplet._convolution_methods
 import sleplet.slepian_methods
@@ -79,10 +79,11 @@ def axisymmetric_wavelet_forward(
     """
     w = np.zeros(wavelets.shape, dtype=np.complex_)
     for ell in range(L):
-        wav_0 = np.sqrt((4 * np.pi) / (2 * ell + 1)) * wavelets[:, ell, L - 1].conj()
+        ind_m0 = ssht.elm2ind(ell, 0)
+        wav_0 = np.sqrt((4 * np.pi) / (2 * ell + 1)) * wavelets[:, ind_m0].conj()
         for m in range(-ell, ell + 1):
-            ind_pm = L - 1 + m
-            w[:, ell, ind_pm] = wav_0 * flm[ell, ind_pm]
+            ind = ssht.elm2ind(ell, m)
+            w[:, ind] = wav_0 * flm[ind]
     return w
 
 
@@ -102,12 +103,13 @@ def axisymmetric_wavelet_inverse(
     Returns:
         Spherical harmonic coefficients of the signal.
     """
-    flm = np.zeros(s2fft.samples.flm_shape(L), dtype=np.complex_)
+    flm = np.zeros(L**2, dtype=np.complex_)
     for ell in range(L):
-        wav_0 = np.sqrt((4 * np.pi) / (2 * ell + 1)) * wavelets[:, ell, L - 1]
+        ind_m0 = ssht.elm2ind(ell, 0)
+        wav_0 = np.sqrt((4 * np.pi) / (2 * ell + 1)) * wavelets[:, ind_m0]
         for m in range(-ell, ell + 1):
-            ind_pm = L - 1 + m
-            flm[ell, ind_pm] = (wav_coeffs[:, ell, ind_pm] * wav_0).sum()
+            ind = ssht.elm2ind(ell, m)
+            flm[ind] = (wav_coeffs[:, ind] * wav_0).sum()
     return flm
 
 
@@ -118,13 +120,11 @@ def _create_axisymmetric_wavelets(
 ) -> npt.NDArray[np.complex_]:
     """Computes the axisymmetric wavelets."""
     kappas = create_kappas(L, B, j_min)
-    wavelets = np.zeros(
-        (kappas.shape[0], *s2fft.samples.flm_shape(L)),
-        dtype=np.complex_,
-    )
+    wavelets = np.zeros((kappas.shape[0], L**2), dtype=np.complex_)
     for ell in range(L):
         factor = np.sqrt((2 * ell + 1) / (4 * np.pi))
-        wavelets[:, ell, L - 1] = factor * kappas[:, ell]
+        ind = ssht.elm2ind(ell, 0)
+        wavelets[:, ind] = factor * kappas[:, ell]
     return wavelets
 
 
@@ -140,15 +140,8 @@ def create_kappas(xlim: int, B: int, j_min: int) -> npt.NDArray[np.float_]:
     Returns:
         The Slepian wavelet generating functions.
     """
-    kappa, kappa0 = s2wav.filter_factory.filters.filters_axisym(
-        xlim,
-        J_min=j_min,
-        lam=B,
-    )
-    kappas = np.concatenate((kappa0[np.newaxis], kappa[j_min:]))
-    # this step is required when migrating from S2LET to S2WAV
-    kappas[kappas == np.inf] = 1
-    return kappas
+    kappa0, kappa = pys2let.axisym_wav_l(B, xlim, j_min)
+    return np.concatenate((kappa0[np.newaxis], kappa.T))
 
 
 def find_non_zero_wavelet_coefficients(
diff --git a/tests/conftest.py b/tests/conftest.py
index c0af9d68b..ee3a18f9b 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -2,8 +2,6 @@
 import numpy.typing as npt
 import pytest
 
-import s2fft
-
 import sleplet
 
 ARRAY_DIM = 10
@@ -107,14 +105,17 @@ def slepian_wavelets_lim_lat_lon(
 @pytest.fixture(scope="session")
 def random_flm() -> npt.NDArray[np.complex_]:
     """Creates random flm."""
-    return s2fft.utils.signal_generator.generate_flm(RNG, L)
+    return sleplet.harmonic_methods.compute_random_signal(L, RNG, var_signal=1)
 
 
 @pytest.fixture(scope="session")
 def random_nd_flm() -> npt.NDArray[np.complex_]:
     """Creates multiple random flm."""
     return np.array(
-        [s2fft.utils.signal_generator.generate_flm(RNG, L) for _ in range(ARRAY_DIM)],
+        [
+            sleplet.harmonic_methods.compute_random_signal(L, RNG, var_signal=1)
+            for _ in range(ARRAY_DIM)
+        ],
     )
 
 
diff --git a/tests/test_decomposition.py b/tests/test_decomposition.py
index 1437af31a..0573478c7 100644
--- a/tests/test_decomposition.py
+++ b/tests/test_decomposition.py
@@ -1,17 +1,16 @@
 import numpy as np
 
-import s2fft
+import pyssht as ssht
 
 import sleplet
 
 
 def test_decompose_all_polar(slepian_polar_cap, earth_polar_cap) -> None:
     """Tests that all three methods produce the same coefficients for polar cap."""
-    field = s2fft.inverse(
+    field = ssht.inverse(
         earth_polar_cap.coefficients,
         slepian_polar_cap.L,
-        method=sleplet._vars.EXECUTION_MODE,
-        sampling=sleplet._vars.SAMPLING_SCHEME,
+        Method=sleplet._vars.SAMPLING_SCHEME,
     )
     harmonic_sum_p = sleplet.slepian_methods.slepian_forward(
         slepian_polar_cap.L,
@@ -46,11 +45,10 @@ def test_decompose_all_lim_lat_lon(slepian_lim_lat_lon, earth_lim_lat_lon) -> No
     Tests that all three methods produce the same coefficients for
     limited latitude longitude region.
     """
-    field = s2fft.inverse(
+    field = ssht.inverse(
         earth_lim_lat_lon.coefficients,
         slepian_lim_lat_lon.L,
-        method=sleplet._vars.EXECUTION_MODE,
-        sampling=sleplet._vars.SAMPLING_SCHEME,
+        Method=sleplet._vars.SAMPLING_SCHEME,
     )
     harmonic_sum_p = sleplet.slepian_methods.slepian_forward(
         slepian_lim_lat_lon.L,
@@ -95,11 +93,10 @@ def test_equality_to_harmonic_transform_polar(
         slepian_polar_cap.L,
         slepian_polar_cap,
     )
-    f_harmonic = s2fft.inverse(
+    f_harmonic = ssht.inverse(
         earth_polar_cap.coefficients,
         slepian_polar_cap.L,
-        method=sleplet._vars.EXECUTION_MODE,
-        sampling=sleplet._vars.SAMPLING_SCHEME,
+        Method=sleplet._vars.SAMPLING_SCHEME,
     )
     mask = sleplet._mask_methods.create_mask_region(
         slepian_polar_cap.L,
@@ -123,11 +120,10 @@ def test_equality_to_harmonic_transform_lim_lat_lon(
         slepian_lim_lat_lon.L,
         slepian_lim_lat_lon,
     )
-    f_harmonic = s2fft.inverse(
+    f_harmonic = ssht.inverse(
         earth_lim_lat_lon.coefficients,
         slepian_lim_lat_lon.L,
-        method=sleplet._vars.EXECUTION_MODE,
-        sampling=sleplet._vars.SAMPLING_SCHEME,
+        Method=sleplet._vars.SAMPLING_SCHEME,
     )
     mask = sleplet._mask_methods.create_mask_region(
         slepian_lim_lat_lon.L,
diff --git a/tests/test_harmonic.py b/tests/test_harmonic.py
index a0ae75a78..b09cbe962 100644
--- a/tests/test_harmonic.py
+++ b/tests/test_harmonic.py
@@ -1,6 +1,6 @@
 import numpy as np
 
-import s2fft
+import pyssht as ssht
 
 import sleplet
 
@@ -11,21 +11,15 @@
 def test_harmonic_coefficients_padded(random_flm) -> None:
     """Tests that harmonic coefficients are zero padded for plotting."""
     boost = L_LARGE**2 - L_SMALL**2
-    flm_boosted = s2fft.samples.flm_1d_to_2d(
-        sleplet.harmonic_methods._boost_coefficient_resolution(
-            s2fft.samples.flm_2d_to_1d(random_flm, L_SMALL),
-            boost,
-        ),
-        L_LARGE,
+    flm_boosted = sleplet.harmonic_methods._boost_coefficient_resolution(
+        random_flm,
+        boost,
     )
-    np.testing.assert_equal(flm_boosted.shape, s2fft.samples.flm_shape(L_LARGE))
+    np.testing.assert_equal(len(flm_boosted), L_LARGE**2)
 
 
 def test_invert_flm_and_boost(random_flm) -> None:
     """Tests that the flm has been boosted and has right shape."""
-    n_theta, n_phi = s2fft.samples.f_shape(
-        L_LARGE,
-        sampling=sleplet._vars.SAMPLING_SCHEME,
-    )
+    n_theta, n_phi = ssht.sample_shape(L_LARGE, Method=sleplet._vars.SAMPLING_SCHEME)
     f = sleplet.harmonic_methods.invert_flm_boosted(random_flm, L_SMALL, L_LARGE)
     np.testing.assert_equal(f.shape, (n_theta, n_phi))
diff --git a/tests/test_translation.py b/tests/test_translation.py
index b48f509dd..b1349e71e 100644
--- a/tests/test_translation.py
+++ b/tests/test_translation.py
@@ -1,7 +1,7 @@
 import hypothesis
 import numpy as np
 
-import s2fft
+import pyssht as ssht
 
 import sleplet
 
@@ -56,13 +56,11 @@ def test_slepian_translation_changes_max_polar(slepian_dirac_delta_polar_cap) ->
         slepian_dirac_delta_polar_cap.slepian,
     )
     new_max = tuple(np.argwhere(field == field.max())[0])
-    thetas = np.tile(
-        s2fft.samples.thetas(
-            slepian_dirac_delta_polar_cap.L,
-            sampling=sleplet._vars.SAMPLING_SCHEME,
-        ),
-        (s2fft.samples.nphi_equiang(L, sampling=sleplet._vars.SAMPLING_SCHEME), 1),
-    ).T
+    thetas, _ = ssht.sample_positions(
+        slepian_dirac_delta_polar_cap.L,
+        Grid=True,
+        Method=sleplet._vars.SAMPLING_SCHEME,
+    )
     np.testing.assert_raises(
         AssertionError,
         np.testing.assert_equal,
@@ -91,13 +89,11 @@ def test_slepian_translation_changes_max_lim_lat_lon(
         slepian_dirac_delta_lim_lat_lon.slepian,
     )
     new_max = tuple(np.argwhere(field == field.max())[0])
-    thetas = np.tile(
-        s2fft.samples.thetas(
-            slepian_dirac_delta_lim_lat_lon.L,
-            sampling=sleplet._vars.SAMPLING_SCHEME,
-        ),
-        (s2fft.samples.nphi_equiang(L, sampling=sleplet._vars.SAMPLING_SCHEME), 1),
-    ).T
+    thetas, _ = ssht.sample_positions(
+        slepian_dirac_delta_lim_lat_lon.L,
+        Grid=True,
+        Method=sleplet._vars.SAMPLING_SCHEME,
+    )
     np.testing.assert_raises(
         AssertionError,
         np.testing.assert_equal,
diff --git a/tests/test_wavelets.py b/tests/test_wavelets.py
index 359259ffc..24c437f7d 100644
--- a/tests/test_wavelets.py
+++ b/tests/test_wavelets.py
@@ -1,6 +1,6 @@
 import numpy as np
 
-import s2wav
+import pys2let
 
 import sleplet
 
@@ -103,6 +103,6 @@ def test_only_wavelet_coefficients_within_shannon_returned() -> None:
 def test_create_kappas() -> None:
     """Checks that the method creates the scaling function and wavelets."""
     wavelets = sleplet.wavelet_methods.create_kappas(L_LARGE**2, B, J_MIN)
-    j_max = s2wav.utils.shapes.j_max(L_LARGE**2, B)
+    j_max = pys2let.pys2let_j_max(B, L_LARGE**2, J_MIN)
     np.testing.assert_equal(j_max - J_MIN + 2, wavelets.shape[0])
     np.testing.assert_equal(L_LARGE**2, wavelets.shape[1])