diff --git a/pyproject.toml b/pyproject.toml
index 60915ba..4e79885 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -8,7 +8,7 @@ build-backend = "poetry.core.masonry.api"
 
 [tool.poetry]
 name = "pyelq-sdk"
-version = "1.0.8b"
+version = "1.0.9"
 description = "Package for detection, localization and quantification code."
 authors = ["Bas van de Kerkhof", "Matthew Jones", "David Randell"]
 homepage = "https://sede-open.github.io/pyELQ/"
diff --git a/src/pyelq/meteorology.py b/src/pyelq/meteorology.py
index 87cc358..22ff420 100644
--- a/src/pyelq/meteorology.py
+++ b/src/pyelq/meteorology.py
@@ -16,7 +16,6 @@
 import plotly.express as px
 import plotly.graph_objects as go
 from pandas.arrays import DatetimeArray
-from scipy.stats import circstd
 
 from pyelq.coordinate_system import Coordinate
 from pyelq.sensor.sensor import SensorGroup
@@ -102,7 +101,12 @@ def calculate_uv_from_wind_speed_direction(self) -> None:
     def calculate_wind_turbulence_horizontal(self, window: str) -> None:
         """Calculate the horizontal wind turbulence values from the wind direction attribute.
 
-        Wind turbulence values are calculated as the circular standard deviation based on a rolling window.
+        Wind turbulence values are calculated as the circular standard deviation of wind direction
+        (https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.circstd.html).
+        The implementation here is equivalent to using the circstd function from scipy.stats as an apply
+        function on a rolling window. However, using the rolling mean on sin and cos speeds up
+        the calculation by a factor of 100.
+
         Outputted values are calculated at the center of the window and at least 3 observations are required in a
         window for the calculation. If the window contains less values the result will be np.nan.
         The result of the calculation will be stored as the wind_turbulence_horizontal attribute.
@@ -113,9 +117,9 @@ def calculate_wind_turbulence_horizontal(self, window: str) -> None:
 
         """
         data_series = pd.Series(data=self.wind_direction, index=self.time)
-        aggregated_data = data_series.rolling(window=window, center=True, min_periods=3).apply(
-            circstd, kwargs={"low": 0, "high": 360}
-        )
+        sin_rolling = (np.sin(data_series * np.pi / 180)).rolling(window=window, center=True, min_periods=3).mean()
+        cos_rolling = (np.cos(data_series * np.pi / 180)).rolling(window=window, center=True, min_periods=3).mean()
+        aggregated_data = np.sqrt(-2 * np.log((sin_rolling**2 + cos_rolling**2) ** 0.5)) * 180 / np.pi
         self.wind_turbulence_horizontal = aggregated_data.values
 
     def plot_polar_hist(self, nof_sectors: int = 16, nof_divisions: int = 5, template: object = None) -> go.Figure():
diff --git a/tests/test_meteorology.py b/tests/test_meteorology.py
index 2c1ebae..56543e9 100644
--- a/tests/test_meteorology.py
+++ b/tests/test_meteorology.py
@@ -201,7 +201,7 @@ def test_meteorology_group():
 def test_calculate_wind_turbulence_horizontal():
     """Checks that the wind turbulence values are calculated correctly.
 
-    To verify circstd, we define winds as draws from a normal distribution. We then check that the mean of the
+    To verify horizontal wind turbulence calculations, we define winds as draws from a normal distribution. We then check that the mean of the
     calculated turbulence values is within 3 standard deviations of the true value.
 
     """
@@ -219,10 +219,8 @@ def test_calculate_wind_turbulence_horizontal():
         pd.date_range(dt.datetime(2023, 1, 1), dt.datetime(2023, 1, 2), freq="5s"), dtype="datetime64[ns]"
     )
     met.wind_direction = np.random.normal(180, sigma, met.time.shape[0])
-
     met.calculate_wind_turbulence_horizontal(window="300s")
 
     tolerance = 3 * np.std(met.wind_turbulence_horizontal)
     mean_turbulence = np.mean(met.wind_turbulence_horizontal)
-
     assert (mean_turbulence - tolerance) < sigma < (mean_turbulence + tolerance)