diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 6b23607f2315757faece92424525b16a34f1904d..e271a11be3f601e01683549c85bc8eedbc87f799 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -35,7 +35,7 @@ stages: # List of stages for jobs, and their order of execution
     - deploy
 
 build:
-    stage: build  
+    stage: build
     <<: *build_test
     script:
         - pip install .
@@ -86,10 +86,10 @@ pages:
     - if: $CI_COMMIT_BRANCH == $CI_DEFAULT_BRANCH
 
 build-singularity:
-  image: 
+  image:
     name: quay.io/singularity/singularity:v3.7.0
     entrypoint: [""]
-  stage: deploy 
+  stage: deploy
   only:
     - master
   script:
@@ -98,4 +98,14 @@ build-singularity:
     - singularity build --remote library://nyhuma/ntrflows/ntr.sif:latest singularitybuild/ntrfc_container.def
   interruptible: true
 
-
+publish:
+  stage: deploy
+  <<: *build_test
+  script:
+    - pip install twine
+    - python setup.py sdist bdist_wheel
+    - twine upload --username $PYPI_USERNAME --password $PYPI_PASSWORD --verbose dist/*
+  only:
+      - tags
+      - master
+  interruptible: true
diff --git a/ntrfc/cascade_case/domain/domain.py b/ntrfc/cascade_case/domain/domain.py
index 3dc658d53c6962857412decc75d8431d754df74f..efd3437efc87a2bb07f1cceecd083183da4a5644 100644
--- a/ntrfc/cascade_case/domain/domain.py
+++ b/ntrfc/cascade_case/domain/domain.py
@@ -77,6 +77,6 @@ class CascadeDomain2D:
 
         plotter.add_axes()
         camera = plotter.camera
-        camera.roll+=90
+        camera.roll += 90
         plotter.show(cpos=(0, 0, 1))
         plotter.update()
diff --git a/ntrfc/cascade_case/solution/case_modules/postprocessing.py b/ntrfc/cascade_case/solution/case_modules/postprocessing.py
index 83cf62329054d1d38bc3195ad7b752f86752c4fd..8d786c4e1f5abf15544ccbfdb7218823c8739397 100644
--- a/ntrfc/cascade_case/solution/case_modules/postprocessing.py
+++ b/ntrfc/cascade_case/solution/case_modules/postprocessing.py
@@ -18,14 +18,16 @@ class PostProcessing:
         self.outlet["v"] = self.outlet["U"][::, 1]
         self.outlet["w"] = self.outlet["U"][::, 2]
         rho_1 = massflowave_plane(self.inlet, "rho")
-        mag_u_1 = vecAbs(np.array([massflowave_plane(self.inlet, "u"),massflowave_plane(self.inlet, "v"),massflowave_plane(self.inlet, "w")]))
+        mag_u_1 = vecAbs(np.array([massflowave_plane(self.inlet, "u"), massflowave_plane(self.inlet, "v"),
+                                   massflowave_plane(self.inlet, "w")]))
         U_1 = np.stack([massflowave_plane(self.inlet, "u"), massflowave_plane(self.inlet, "v"),
                         massflowave_plane(self.inlet, "w")])
         beta_1 = vecAngle(U_1, np.array([1, 0, 0]))
         rho_2 = massflowave_plane(self.outlet, "rho")
         U_2 = np.stack([massflowave_plane(self.outlet, "u"), massflowave_plane(self.outlet, "v"),
                         massflowave_plane(self.outlet, "w")])
-        mag_u_2 =  vecAbs(np.array([massflowave_plane(self.outlet, "u"),massflowave_plane(self.outlet, "v"),massflowave_plane(self.outlet, "w")]))
+        mag_u_2 = vecAbs(np.array([massflowave_plane(self.outlet, "u"), massflowave_plane(self.outlet, "v"),
+                                   massflowave_plane(self.outlet, "w")]))
         beta_2 = vecAngle(U_2, np.array([1, 0, 0]))
         self.avdr = avdr(rho_1, mag_u_1, beta_1, rho_2, mag_u_2, beta_2)
 
diff --git a/ntrfc/fluid/fluid.py b/ntrfc/fluid/fluid.py
index f7a9c398ef9ce638364c4115b3b825206ca9b3a3..d587746d6837f404ec809bfe11d44b372414c19b 100644
--- a/ntrfc/fluid/fluid.py
+++ b/ntrfc/fluid/fluid.py
@@ -1,6 +1,3 @@
-import numpy as np
-
-
 def sutherland_viscosity(temperature, sutherland_constant=1.458e-06, reference_temperature=110.4):
     """
     Calculate the dynamic viscosity of a gas using Sutherland's Law.
@@ -48,6 +45,5 @@ def mach_number(speed, specific_heat_ratio, specific_gas_constant, temperature):
     Returns:
     - the Mach number at the point
     """
-    mach_number = speed/ ((specific_heat_ratio * specific_gas_constant * temperature) ** 0.5)
+    mach_number = speed / ((specific_heat_ratio * specific_gas_constant * temperature) ** 0.5)
     return mach_number
-
diff --git a/ntrfc/fluid/isentropic.py b/ntrfc/fluid/isentropic.py
index be2e385264f9f90d459105af0d657c1031fa9fe9..523e6dce405b288035565d44722c8b72a3d188d3 100644
--- a/ntrfc/fluid/isentropic.py
+++ b/ntrfc/fluid/isentropic.py
@@ -1,6 +1,7 @@
-from ntrfc.fluid.fluid import mach_number, total_pressure
 import numpy as np
 
+from ntrfc.fluid.fluid import mach_number, total_pressure
+
 
 def p_t_is_from_mach(kappa, mach_number, static_pressure):
     # Calculates total pressure in isentropic flow
@@ -22,6 +23,7 @@ def temp_t_is(kappa, ma, T):
     T_t_is = T / (((1.0 + (kappa - 1.0) * 0.5 * ma ** 2.0)) ** (-1.0))
     return T_t_is
 
+
 def temp_is(kappa, ma, Tt):
     # Calculates static temperature in isentropic flow
     # Source: https://www.grc.nasa.gov/www/BGH/isentrop.html
@@ -38,7 +40,7 @@ def mach_is_x(kappa, p, pt):
     return y
 
 
-def isentropic_mach_number(isentropic_pressure, kappa, static_pressure,  mach, gas_constant, static_temperature):
+def isentropic_mach_number(isentropic_pressure, kappa, static_pressure, mach, gas_constant, static_temperature):
     """
     Calculates the isentropic Mach number.
 
@@ -62,7 +64,8 @@ def isentropic_mach_number(isentropic_pressure, kappa, static_pressure,  mach, g
     dynamic_pressure = total_pressure - isentropic_pressure
 
     # Calculate the isentropic Mach number
-    isentropic_mach_number = np.sqrt(2.0 / (kappa - 1.0) * (pow(1.0 + (dynamic_pressure / isentropic_pressure), (kappa - 1.0) / kappa) - 1.0))
+    isentropic_mach_number = np.sqrt(
+        2.0 / (kappa - 1.0) * (pow(1.0 + (dynamic_pressure / isentropic_pressure), (kappa - 1.0) / kappa) - 1.0))
 
     return isentropic_mach_number
 
diff --git a/ntrfc/fluid/turbulence.py b/ntrfc/fluid/turbulence.py
index 9716f7809443dbb39d592c5045462844e6d92b9d..ab99c5e82285f8d3eb234ef795f244cb228ae3e1 100644
--- a/ntrfc/fluid/turbulence.py
+++ b/ntrfc/fluid/turbulence.py
@@ -15,11 +15,11 @@ def calcTkeByTu(Tu, Uabs):
     return Tke
 
 
-def calcTke(u_2,v_2,w_2):
-    #u_2: u-fluktuationen quadriert und gemittelt
-    #v_2: v-fluktuationen quadriert und gemittelt
-    #w_2: w-fluktuationen quadriert und gemittelt
-    tke=( 0.5 * (u_2 + v_2 + w_2) )
+def calcTke(u_2, v_2, w_2):
+    # u_2: u-fluktuationen quadriert und gemittelt
+    # v_2: v-fluktuationen quadriert und gemittelt
+    # w_2: w-fluktuationen quadriert und gemittelt
+    tke = (0.5 * (u_2 + v_2 + w_2))
     return tke
 
 
@@ -51,5 +51,3 @@ def calcRey(u, v, w):
     mean_v_w_ = np.mean(v_w_)
 
     return [mean_u_u_, mean_v_v_, mean_w_w_, mean_u_v_, mean_u_w_, mean_v_w_]
-
-
diff --git a/ntrfc/meshquality/nondimensionals.py b/ntrfc/meshquality/nondimensionals.py
index 3b80fad8a4e947125dd9e8442891b58e58ae6197..8fb14273121d12b4fb1f312bc936857cc7f50d88 100644
--- a/ntrfc/meshquality/nondimensionals.py
+++ b/ntrfc/meshquality/nondimensionals.py
@@ -5,14 +5,12 @@ Created on Fri Apr 24 23:33:53 2020
 @author: malte
 """
 
-from itertools import combinations
-
 import numpy as np
 import pyvista as pv
 import vtk
 from scipy.spatial import KDTree
 
-from ntrfc.math.vectorcalc import vecAbs, vecAbs_list, vecProjection, unitvec
+from ntrfc.math.vectorcalc import vecAbs_list, unitvec
 
 
 def cellDirections(cellUMean, wallNorm):
@@ -63,10 +61,11 @@ def get_wall_shear_stress_velocity(mesh, dynamic_viscosity, density_fieldname, v
 
     grad_velocity = mesh[f"grad_{velocity_fieldname}"].reshape(mesh.number_of_cells, 3, 3)
     wall_normals = mesh["wallNormal"]
-    velocity_gradient_normal = vecAbs_list([np.dot(grad_velocity[i], wall_normals[i]) for i in range(mesh.number_of_cells)])
+    velocity_gradient_normal = vecAbs_list(
+        [np.dot(grad_velocity[i], wall_normals[i]) for i in range(mesh.number_of_cells)])
     fluid_density = mesh[density_fieldname]
-    tau_w = dynamic_viscosity/fluid_density*velocity_gradient_normal
-    return np.sqrt(tau_w/fluid_density)
+    tau_w = dynamic_viscosity / fluid_density * velocity_gradient_normal
+    return np.sqrt(tau_w / fluid_density)
 
 
 def constructWallMesh(surfaces):
@@ -121,7 +120,7 @@ def calc_dimensionless_gridspacing(volmesh, surfaces, use_velfield, use_rhofield
     print("finding walladjacent cells")
     merged_mesh = surfacenormals_surface + volmesh
     # Compute cell-based derivative of a vector field
-    derivatives = compute_scalar_gradient(merged_mesh,"U")
+    derivatives = compute_scalar_gradient(merged_mesh, "U")
 
     derivatives = derivatives.ctp()
     # Extract the volumetric mesh of the derivative-field
@@ -135,7 +134,7 @@ def calc_dimensionless_gridspacing(volmesh, surfaces, use_velfield, use_rhofield
     volmesh_walladjacent = volmesh.extract_cells(walladjacentids)
     volmesh_walladjacent["cellCenters"] = volmesh_walladjacent.cell_centers().points
 
-    facemesh=facemesh.extract_surface().compute_normals()
+    facemesh = facemesh.extract_surface().compute_normals()
     # Extract the cell centers of the wall-adjacent cells and the surface cells
     volmesh_walladjacent_centers = volmesh_walladjacent.cell_centers().points
     surfacemesh_surface_centers = facemesh.cell_centers().points
@@ -152,7 +151,6 @@ def calc_dimensionless_gridspacing(volmesh, surfaces, use_velfield, use_rhofield
     volmesh_walladjacent["wallNormal"] = surface_normals[indices]
     volmesh_walladjacent["wallUGradient"] = surface_gradients[indices]
 
-
     print("calculating cell spans from WallNormals...")
     spanS = cellSpans(volmesh_walladjacent, use_velfield)
     volmesh_walladjacent["xSpan"] = np.array([i[0] for i in spanS])  # calculate cell span in flow direction
@@ -167,8 +165,8 @@ def calc_dimensionless_gridspacing(volmesh, surfaces, use_velfield, use_rhofield
         [np.dot(grad_velocity[i], wall_normals[i]) for i in range(volmesh_walladjacent.number_of_cells)])
     fluid_density = volmesh_walladjacent[use_rhofield]
     u_tau = np.sqrt(mu_0 * velocity_gradient_normal / fluid_density)
-    #uTaus = volmesh_walladjacent["uTaus"]
-    d_v = mu_0/volmesh_walladjacent[use_rhofield] / u_tau
+    # uTaus = volmesh_walladjacent["uTaus"]
+    d_v = mu_0 / volmesh_walladjacent[use_rhofield] / u_tau
     volmesh_walladjacent["DeltaXPlus"] = volmesh_walladjacent[
                                              "xSpan"] / d_v  # calculate cell span in flow direction
     volmesh_walladjacent["DeltaYPlus"] = volmesh_walladjacent[
diff --git a/ntrfc/timeseries/stationarity.py b/ntrfc/timeseries/stationarity.py
index 90c65f93e564bd55396bfc7f347af2da5a8c49eb..22edaca1ea0a5127bbe09f9a0679f9bebe2ce9aa 100644
--- a/ntrfc/timeseries/stationarity.py
+++ b/ntrfc/timeseries/stationarity.py
@@ -1,15 +1,12 @@
 import numpy as np
+import pandas as pd
 from matplotlib import pyplot as plt
 from scipy import signal
 from sklearn.decomposition import PCA
-from statsmodels.tsa.arima.model import ARIMA
-import pandas as pd
 
 from ntrfc.math.methods import minmax_normalize
 
 
-
-
 def optimal_bin_width(sample1, sample2):
     """
     Compute the optimal bin width using cross-validation estimator for mean integrated squared error.
@@ -46,8 +43,6 @@ def optimal_bin_width(sample1, sample2):
     return bins
 
 
-
-
 def smd_probability_compare(sample1, sample2, verbose=False):
     """Compare the probability distribution of two signals using the Freedman-Diaconis rule
     to determine the number of bins.
@@ -114,17 +109,15 @@ def optimal_window_size(time_series, min_interval=0.05, max_interval=0.25, verbo
 
     # Get the length of the time series and define a range of allowed window sizes
     series_length = len(normalized_series)
-    allowed_window_sizes = np.array(range(int(series_length*min_interval), int(series_length *max_interval)+1))
-
+    allowed_window_sizes = np.array(range(int(series_length * min_interval), int(series_length * max_interval) + 1))
 
     # Iterate through the allowed window sizes and perform the Kolmogorov-Smirnov test and correlation coefficient calculation
     mean_scores = []
     var_scores = []
     for window_size in allowed_window_sizes:
-
         check_window = normalized_series[-window_size * 2:]
         check_window_df = pd.DataFrame(check_window)
-        rolling_check =  check_window_df.rolling(window_size)
+        rolling_check = check_window_df.rolling(window_size)
         mean_scores.append(np.std(rolling_check.mean()).values[0])
         var_scores.append(np.std(rolling_check.var()).values[0])
         # Compute the correlation coefficient
@@ -139,8 +132,8 @@ def optimal_window_size(time_series, min_interval=0.05, max_interval=0.25, verbo
     assert len(opt_window) == opt_window_size * 2
     probability_similiarity = smd_probability_compare(opt_window[:opt_window_size], opt_window[opt_window_size:])
 
-    if probability_similiarity< 0.99:
-        return False,False,False
+    if probability_similiarity < 0.99:
+        return False, False, False
 
     # Compute the period of the time series
     freqs, psd = signal.welch(time_series, fs=1, nperseg=series_length // 4)
@@ -148,7 +141,7 @@ def optimal_window_size(time_series, min_interval=0.05, max_interval=0.25, verbo
 
     if maxpsdid != 0:
         tperiod = freqs[np.argmax(psd)] ** -1
-        nperiods = (opt_window_size+(-opt_window_size%tperiod))//tperiod
+        nperiods = (opt_window_size + (-opt_window_size % tperiod)) // tperiod
     else:
         tperiod = np.inf
         nperiods = 0
@@ -163,10 +156,9 @@ def optimal_window_size(time_series, min_interval=0.05, max_interval=0.25, verbo
     return opt_window, opt_window_size, nperiods
 
 
-
-def estimate_stationarity(timeseries,  verbose=False):
+def estimate_stationarity(timeseries, verbose=False):
     sigma_threshold = 3
-    normalized_series = minmax_normalize(timeseries)#minmax_normalize(timeseries)
+    normalized_series = minmax_normalize(timeseries)  # minmax_normalize(timeseries)
     datalength = len(normalized_series)
     opt_window, opt_window_size, nperiods = optimal_window_size(normalized_series)
     if not opt_window_size:
@@ -187,7 +179,6 @@ def estimate_stationarity(timeseries,  verbose=False):
     mean_uncertainty = np.std(rolling_means)
     var_uncertainty = np.std(rolling_vars)
 
-
     checkseries = normalized_series[:-opt_window_size * 2]
 
     checkseries_reversed = pd.DataFrame(checkseries[::-1])
@@ -197,13 +188,13 @@ def estimate_stationarity(timeseries,  verbose=False):
     rolling_vars_reversed = rolling_win_reversed.var().values
 
     outer_mean_uncertainty = 0.0015
-    outer_var_uncertainty = outer_mean_uncertainty *0.25# variance can only be 25% of minmax normed series
+    outer_var_uncertainty = outer_mean_uncertainty * 0.25  # variance can only be 25% of minmax normed series
 
     rolling_means_errors_reversed = np.abs(rolling_means_reversed - reference_mean)
     rolling_vars_errors_reversed = np.abs(rolling_vars_reversed - reference_variance)
 
-    mean_limits = sigma_threshold * (mean_uncertainty)+outer_mean_uncertainty
-    var_limits = sigma_threshold * (var_uncertainty)+outer_var_uncertainty
+    mean_limits = sigma_threshold * (mean_uncertainty) + outer_mean_uncertainty
+    var_limits = sigma_threshold * (var_uncertainty) + outer_var_uncertainty
 
     rolling_means_errors_inliers_reversed = rolling_means_errors_reversed[opt_window_size:] <= mean_limits
     rolling_vars_errors_inliers_reversed = rolling_vars_errors_reversed[opt_window_size:] <= var_limits
@@ -217,28 +208,32 @@ def estimate_stationarity(timeseries,  verbose=False):
         answer_index = np.where(success_rate >= threshold)[0][0]
         return len(arr) - answer_index
 
-    mean_index = datalength - last_coherent_interval(rolling_means_errors_inliers_reversed)-3*opt_window_size
+    mean_index = datalength - last_coherent_interval(rolling_means_errors_inliers_reversed) - 3 * opt_window_size
 
-    variance_index = datalength - last_coherent_interval(rolling_vars_errors_inliers_reversed)-3*opt_window_size
+    variance_index = datalength - last_coherent_interval(rolling_vars_errors_inliers_reversed) - 3 * opt_window_size
 
     stationary_start_index = max(mean_index, variance_index)
 
     if verbose:
-
-        fig ,axs = plt.subplots(3,1,figsize=(24,20))
+        fig, axs = plt.subplots(3, 1, figsize=(24, 20))
         axs[0].plot(normalized_series, label="normalized series", color="blue")
-        axs[0].vlines(x=stationary_start_index, ymin=0, ymax=np.nanmax(normalized_series), label="stationary_start",color="k")
+        axs[0].vlines(x=stationary_start_index, ymin=0, ymax=np.nanmax(normalized_series), label="stationary_start",
+                      color="k")
 
         axs[1].plot(np.nan_to_num(rolling_means_errors_reversed[::-1], 0), label="mean error", color="red")
-        axs[1].hlines(y=(mean_limits), xmin=0, xmax=len(normalized_series), label="mean_limits",color="k")
-        axs[1].vlines(x=stationary_start_index, ymin=0, ymax=max(mean_limits, np.nanmax(rolling_means_errors_reversed)), label="stationary_start", color="k")
-        axs[1].vlines(x=mean_index, ymin=0, ymax=max(mean_limits, np.nanmax(rolling_means_errors_reversed)), label="mean_index", color="green")
+        axs[1].hlines(y=(mean_limits), xmin=0, xmax=len(normalized_series), label="mean_limits", color="k")
+        axs[1].vlines(x=stationary_start_index, ymin=0, ymax=max(mean_limits, np.nanmax(rolling_means_errors_reversed)),
+                      label="stationary_start", color="k")
+        axs[1].vlines(x=mean_index, ymin=0, ymax=max(mean_limits, np.nanmax(rolling_means_errors_reversed)),
+                      label="mean_index", color="green")
         axs[1].legend()
 
         axs[2].plot(np.nan_to_num(rolling_vars_errors_reversed[::-1], 0), label="variance error", color="red")
         axs[2].hlines(y=(var_limits), xmin=0, xmax=len(normalized_series), label="var_limits", color="k")
-        axs[2].vlines(x=stationary_start_index, ymin=0, ymax=max(var_limits, np.nanmax(rolling_vars_errors_reversed)), label="stationary_start", color="k")
-        axs[2].vlines(x=variance_index, ymin=0, ymax=max(var_limits, np.nanmax(rolling_vars_errors_reversed)), label="variance_index", color="k")
+        axs[2].vlines(x=stationary_start_index, ymin=0, ymax=max(var_limits, np.nanmax(rolling_vars_errors_reversed)),
+                      label="stationary_start", color="k")
+        axs[2].vlines(x=variance_index, ymin=0, ymax=max(var_limits, np.nanmax(rolling_vars_errors_reversed)),
+                      label="variance_index", color="k")
         axs[2].legend()
 
         plt.show()
@@ -246,7 +241,6 @@ def estimate_stationarity(timeseries,  verbose=False):
     return stationary_start_index
 
 
-
 def estimate_error_jacknife(timeseries, block_size=20, n_samples=4000):
     """
     Estimates the errors of the mean, variance, and autocorrelation of a given time series using jackknife resampling method.
diff --git a/ntrfc/turbo/bladeloading.py b/ntrfc/turbo/bladeloading.py
index 5f7c1b68571fe8f7f2aa44778c51bdfb8ade9bf4..9ba30ce754ca125bd12224ae80253053605a70ab 100644
--- a/ntrfc/turbo/bladeloading.py
+++ b/ntrfc/turbo/bladeloading.py
@@ -8,8 +8,9 @@ def calc_inflow_cp(px, pt1, p1):
     cp = (px - pt1) / (p1 - pt1)
     return cp
 
+
 # Totaldruckverlustbeiwert
-def calc_zeta(pt1,pt2x,p2):
+def calc_zeta(pt1, pt2x, p2):
     """
     Calculates the Total Pressure Loss Coefficient (Zeta) for a fluid system.
 
@@ -23,7 +24,7 @@ def calc_zeta(pt1,pt2x,p2):
     --------
     zeta : float : Total Pressure Loss Coefficient (Zeta) [dimensionless]
     """
-    zeta=(pt1-pt2x)/(pt1-p2)
+    zeta = (pt1 - pt2x) / (pt1 - p2)
     return zeta
 
 
@@ -42,5 +43,3 @@ def calc_cf(tau_w, u_inf, rho_inf):
     """
     cf = 2 * tau_w / (rho_inf * u_inf ** 2)
     return cf
-
-
diff --git a/requirements_dev.txt b/requirements_dev.txt
index 517cdd91259bc1633f4f0b42bb781619d1cdccd2..90fc3b80ed25a3a60f17a6eef1786c077c89ef16 100644
--- a/requirements_dev.txt
+++ b/requirements_dev.txt
@@ -2,3 +2,4 @@ pipreqs==0.4.11
 pytest-cov
 coverage
 profilehooks
+twine
diff --git a/setup.py b/setup.py
index 25ca34ee0f504361285915b0164ef722832d4dc8..8d15512da68d59cda0fbec3fdfcb2f8183164266 100644
--- a/setup.py
+++ b/setup.py
@@ -22,7 +22,7 @@ setup(
     author_email='nyhuis@tfd.uni-hannover.de',
     python_requires='>=3.10',
     classifiers=[
-        'Development Status :: 0.1.3',
+        'Development Status :: 3 - Alpha',
         'Intended Audience :: Developers',
         'License :: OSI Approved :: MIT License',
         'Natural Language :: English',
diff --git a/tests/cascadecase/solution/test_ansys.py b/tests/cascadecase/solution/test_ansys.py
index 766ce254774e6e6c2f9f85d281609fb5bc4214d9..16219dad3a84aeede3472a047913da8d893860fe 100644
--- a/tests/cascadecase/solution/test_ansys.py
+++ b/tests/cascadecase/solution/test_ansys.py
@@ -9,7 +9,7 @@ def test_solution(tmpdir):
     fake_inlet["v"] = np.array([0] * fake_inlet.number_of_cells)
     fake_inlet["w"] = np.array([0] * fake_inlet.number_of_cells)
     fake_inlet["rho"] = np.array([1] * fake_inlet.number_of_cells)
-    fake_inlet["U"] = np.stack([fake_inlet["u"],fake_inlet["v"],fake_inlet["w"]]).T
+    fake_inlet["U"] = np.stack([fake_inlet["u"], fake_inlet["v"], fake_inlet["w"]]).T
 
     fake_outlet = pv.Plane()
 
diff --git a/tests/cascadecase/solution/test_openfoam.py b/tests/cascadecase/solution/test_openfoam.py
index dec1ffef36a164546b073d05f4547bcec02b96f0..539dfcf32db234b2b0fb6bc2b8c871e312643a39 100644
--- a/tests/cascadecase/solution/test_openfoam.py
+++ b/tests/cascadecase/solution/test_openfoam.py
@@ -8,7 +8,7 @@ def test_solution():
     fake_inlet["v"] = np.array([0] * fake_inlet.number_of_cells)
     fake_inlet["w"] = np.array([0] * fake_inlet.number_of_cells)
     fake_inlet["rho"] = np.array([1] * fake_inlet.number_of_cells)
-    fake_inlet["U"] = np.stack([fake_inlet["u"],fake_inlet["v"],fake_inlet["w"]]).T
+    fake_inlet["U"] = np.stack([fake_inlet["u"], fake_inlet["v"], fake_inlet["w"]]).T
 
     fake_outlet = pv.Plane()
 
diff --git a/tests/cascadecase/solution/test_solution.py b/tests/cascadecase/solution/test_solution.py
index 043c2141c1944db3fc9a03eaba83521ddfd82366..d0da204a0063ab905f34010d4d46207690dbdbec 100644
--- a/tests/cascadecase/solution/test_solution.py
+++ b/tests/cascadecase/solution/test_solution.py
@@ -20,7 +20,7 @@ def test_solution(tmpdir):
     fake_inlet["v"] = np.array([0] * fake_inlet.number_of_cells)
     fake_inlet["w"] = np.array([0] * fake_inlet.number_of_cells)
     fake_inlet["rho"] = np.array([1] * fake_inlet.number_of_cells)
-    fake_inlet["U"] = np.stack([fake_inlet["u"],fake_inlet["v"],fake_inlet["w"]]).T
+    fake_inlet["U"] = np.stack([fake_inlet["u"], fake_inlet["v"], fake_inlet["w"]]).T
     fake_inlet.save(inletname)
     fake_outlet = pv.Plane(direction=(-1, 0, 0))
 
diff --git a/tests/fluid/test_ntrfc_fluid.py b/tests/fluid/test_ntrfc_fluid.py
index 15a3857319c17562fbff76e0eba62d43ed623166..121c51318edc97dbed1f2b6cc33c684452398bb6 100644
--- a/tests/fluid/test_ntrfc_fluid.py
+++ b/tests/fluid/test_ntrfc_fluid.py
@@ -53,5 +53,3 @@ def test_isentropic_total_temperature():
     expected_output = 315
     output = isentropic_total_temperature(1.4, 0.5, 300)
     assert abs(output - expected_output) < 1e-2
-
-
diff --git a/tests/fluid/test_ntrfc_isentropic.py b/tests/fluid/test_ntrfc_isentropic.py
index 7675116fc93b5dc1e280eca0e5b742c91c353e19..8bf62737ad36262ee5e6d3accd89677b22a91832 100644
--- a/tests/fluid/test_ntrfc_isentropic.py
+++ b/tests/fluid/test_ntrfc_isentropic.py
@@ -24,7 +24,7 @@ def test_p_is_from_mach():
     kappa = 1.4
     mach_number = 0.8
     total_pressure = 101325.0  # Pa
-    expected_static_pressure = 66471.39048022314 # Pa
+    expected_static_pressure = 66471.39048022314  # Pa
     calculated_static_pressure = p_is_from_mach(kappa, mach_number, total_pressure)
     assert np.isclose(calculated_static_pressure, expected_static_pressure)
 
@@ -56,14 +56,15 @@ def test_mach_is_x():
 def test_isentropic_mach_number():
     from ntrfc.fluid.isentropic import isentropic_mach_number
     # Standard conditions
-    isentropic_pressure = 101325 # Pa
+    isentropic_pressure = 101325  # Pa
     kappa = 1.4
-    static_pressure = 101325 # Pa
+    static_pressure = 101325  # Pa
     mach_number = 0.5
-    gas_constant = 287.058 # J/(kg*K)
-    static_temperature = 288.15 # K
+    gas_constant = 287.058  # J/(kg*K)
+    static_temperature = 288.15  # K
 
     expected_output = 0.0014693046301270448
 
     # Call function and check output
-    assert np.isclose(isentropic_mach_number(isentropic_pressure, kappa, static_pressure, mach_number, gas_constant, static_temperature), expected_output)
+    assert np.isclose(isentropic_mach_number(isentropic_pressure, kappa, static_pressure, mach_number, gas_constant,
+                                             static_temperature), expected_output)
diff --git a/tests/fluid/test_ntrfc_turbulence.py b/tests/fluid/test_ntrfc_turbulence.py
index 90233c7f1b649f6a387445b3aa24bfd03c8416b9..908c6f990e961c9889bb36828733da6edd36f9db 100644
--- a/tests/fluid/test_ntrfc_turbulence.py
+++ b/tests/fluid/test_ntrfc_turbulence.py
@@ -37,19 +37,19 @@ def test_calcTke():
 
 
 def test_calcFluc():
-    series = np.random.randn(1000) # calculated as 0.5 * (u_2 + v_2 + w_2)
-    assert np.isclose(np.mean(calcFluc(series)),0.0)
-    series = np.random.randn(1000)+1 # calculated as 0.5 * (u_2 + v_2 + w_2)
-    assert np.isclose(np.mean(calcFluc(series)),0.0)
+    series = np.random.randn(1000)  # calculated as 0.5 * (u_2 + v_2 + w_2)
+    assert np.isclose(np.mean(calcFluc(series)), 0.0)
+    series = np.random.randn(1000) + 1  # calculated as 0.5 * (u_2 + v_2 + w_2)
+    assert np.isclose(np.mean(calcFluc(series)), 0.0)
 
 
 def test_calcRey():
-    series = np.stack([np.random.randn(20000),np.random.randn(20000),np.random.randn(20000)])
-
-    uu,vv,ww,uv,uw,vw = calcRey(series[0],series[1],series[2])
-    assert np.isclose(uu,1.0,atol=5e-2)
-    assert np.isclose(vv,1.0,atol=5e-2)
-    assert np.isclose(ww,1.0,atol=5e-2)
-    assert np.isclose(uv,0.0,atol=5e-2)
-    assert np.isclose(uw,0.0,atol=5e-2)
-    assert np.isclose(vw,0.0,atol=5e-2)
+    series = np.stack([np.random.randn(20000), np.random.randn(20000), np.random.randn(20000)])
+
+    uu, vv, ww, uv, uw, vw = calcRey(series[0], series[1], series[2])
+    assert np.isclose(uu, 1.0, atol=5e-2)
+    assert np.isclose(vv, 1.0, atol=5e-2)
+    assert np.isclose(ww, 1.0, atol=5e-2)
+    assert np.isclose(uv, 0.0, atol=5e-2)
+    assert np.isclose(uw, 0.0, atol=5e-2)
+    assert np.isclose(vw, 0.0, atol=5e-2)
diff --git a/tests/timeseries/test_ntrfc_stationarity.py b/tests/timeseries/test_ntrfc_stationarity.py
index 1c97067624e1311f95b6d7f270c3f77507ed4d0f..ebd9e25cb757cd92935eff7e5bfb6f53d3a6f8f2 100644
--- a/tests/timeseries/test_ntrfc_stationarity.py
+++ b/tests/timeseries/test_ntrfc_stationarity.py
@@ -17,7 +17,7 @@ def test_optimal_timewindow(verbose=False):
     # we have four periods, at least one period should be captured
     # thats res // 4 as a return
     ysin = np.sin(x)
-    opt_window, opt_window_size, nperiods = optimal_window_size(ysin,min_interval, max_interval)
+    opt_window, opt_window_size, nperiods = optimal_window_size(ysin, min_interval, max_interval)
     assert opt_window_size == res // nper * nperiods
     assert nperiods == 2
     # tanh
@@ -76,7 +76,6 @@ def test_stationarity_uncertainties_stationarysine(verbose=False):
                               time - timesteps[stationary_timestep]), "computation failed"
 
 
-
 def test_stationarity_uncertainties_abatingsine(verbose=False):
     from ntrfc.timeseries.stationarity import estimate_stationarity
     from ntrfc.math.methods import reldiff
@@ -112,7 +111,7 @@ def test_stationarity_uncertainties_abatingsine(verbose=False):
             plt.figure()
             plt.plot(timesteps, values)
             plt.axvline(timesteps[stationary_timestep], color="green")
-            plt.axvline(well_computed_stationarity_limit, color="red",label="computed")
+            plt.axvline(well_computed_stationarity_limit, color="red", label="computed")
             plt.legend()
             plt.show()
         assert 0.05 > reldiff(stationary_time, well_computed_stationary_time), "computation failed"
diff --git a/tests/turbo/test_ntrfc_bladeloading.py b/tests/turbo/test_ntrfc_bladeloading.py
index 09e7effe05c46580cbfab95063b217d4bf26b14d..37bf90d9ff55e4a8b61a2c6d860484d40f33e460 100644
--- a/tests/turbo/test_ntrfc_bladeloading.py
+++ b/tests/turbo/test_ntrfc_bladeloading.py
@@ -39,7 +39,7 @@ def test_calc_cf():
     rho_inf = 1.225
     expected_cf = 0.0008163265306122448
     actual_cf = calc_cf(tau_w, u_inf, rho_inf)
-    assert np.isclose(actual_cf,expected_cf)
+    assert np.isclose(actual_cf, expected_cf)
 
     # Test case 2
     tau_w = 0.1
@@ -47,4 +47,4 @@ def test_calc_cf():
     rho_inf = 1.5
     expected_cf = 0.0003333333333333334
     actual_cf = calc_cf(tau_w, u_inf, rho_inf)
-    assert np.isclose(actual_cf,expected_cf)
+    assert np.isclose(actual_cf, expected_cf)