diff --git a/doc/spectrum.png b/docs/spectrum.png
similarity index 100%
rename from doc/spectrum.png
rename to docs/spectrum.png
diff --git a/doc/u.png b/docs/u.png
similarity index 100%
rename from doc/u.png
rename to docs/u.png
diff --git a/pyTurbulence/poissonSolver.py b/pyTurbulence/poissonSolver.py
index f249b06d59ebc7d8c01be06538c07d4797c8b39b..e16dbff150a20ac3e3f3dade4c5f644f0433caca 100644
--- a/pyTurbulence/poissonSolver.py
+++ b/pyTurbulence/poissonSolver.py
@@ -48,7 +48,8 @@ def poisson_2d(f, nx, ny, hx, hy)->np.ndarray:
     res[:, ny-1] = res[:, 0]
     return res
 
-def poisson_3d(f, nx, ny, nz, hx, hy, hz)->np.ndarray:
+def poisson_3d(f : np.ndarray, nx : int, ny : int, nz : int, 
+               hx : float, hy : float, hz : float)->np.ndarray:
     """
     Solve a 3D Poisson equation with periodic boundary conditions.
     
@@ -102,7 +103,8 @@ def poisson_3d(f, nx, ny, nz, hx, hy, hz)->np.ndarray:
     res[:, :, nz-1] = res[:, :, 0]
     return res
 
-def dilatational_velocities(theta, nx, ny, nz, hx, hy, hz)->np.ndarray:
+def dilatational_velocities(theta : np.ndarray, nx : int, ny : int, nz : int, 
+                            hx : float, hy : float, hz : float)->np.ndarray:
     """
     Compute the dilatational velocity field using Fourier space projection.
     
diff --git a/pyTurbulence/syntheticTurbulence.py b/pyTurbulence/syntheticTurbulence.py
index e2b7dfef178fba5467eaf58ecce11c4981be269b..ee0962739212f2063e5ed352280a9176cb5a579a 100644
--- a/pyTurbulence/syntheticTurbulence.py
+++ b/pyTurbulence/syntheticTurbulence.py
@@ -280,17 +280,12 @@ def compute_dilatational_velocities(u: np.ndarray,v: np.ndarray,w: np.ndarray,p:
     hessian_rhs = np.array([np.gradient(grad, hx, hy, hz, edge_order=2) for grad in grad_rhs])
     rhs_ij      = np.einsum('ij...->...', hessian_rhs) 
     rhs_ij      *= 2.0
-
-    print(f"mean rhs_ij = {np.mean(rhs_ij)}")
-
     ddt_p = poisson_3d(rhs_ij, nx, ny, nz, hx, hy, hz)
 
     # Compute dilatation from Eq.7
     # -gamma*dilatation = ddt_p+v_k*dpdxk
     dilatation = -overGamma*(np.sum(velocity * gradPressure, axis=0) + ddt_p)
 
-    print(f"mean dilatation = {np.mean(dilatation)}")
-
     # Compute dilatational velocity from Eq.12
     vd_x, vd_y, vd_z = dilatational_velocities(dilatation, nx, ny, nz, hx, hy, hz)
     return vd_x, vd_y, vd_z