diff --git a/TurTLE/test/test_Gaussian_field.py b/TurTLE/test/test_Gaussian_field.py
index e1d63989c6c1cfac70bbaad033a72bf82011df59..a47f20ff97c3659197561e42f05388340fa0ce88 100644
--- a/TurTLE/test/test_Gaussian_field.py
+++ b/TurTLE/test/test_Gaussian_field.py
@@ -68,12 +68,12 @@ def main_convergence_test():
 
     bin_no = 100
 
-    if not os.path.exists(simname+'.h5'):
+    if not os.path.exists('conv_test_n64_run0.h5'):
         for n, dk in zip(n_arr, dk_arr):
             for run in range(3):
-                rseed = int(time.time())
                 simname = f'conv_test_n{n}_run{run}'
 
+                rseed = int(time.time())
                 c = TEST()
                 opt = c.launch(
                     ['Gauss_field_test',
@@ -100,27 +100,40 @@ def main_convergence_test():
                      sys.argv[1:])
     else:
         plt.figure()
-
+        k = np.geomspace(1e-1, 4e1, 1000)
+        plt.plot(k, E(k, epsilon=diss, Lint=Lint, etaK=etaK, c_L=c_L, c_eta=c_eta))
         for n in n_arr:
-            for run in range(3):
-                df = h5py.File(simname + '.h5', 'r')
-        
+            run = 0
+            simname = f'conv_test_n{n}_run{run}'
+            df = h5py.File(simname + '.h5', 'r')
+
+            kk = df['kspace/kshell'][...]
+            phi_ij = df['statistics/spectra/velocity_velocity'][0] / df['kspace/dk'][()]
+            energy_spec = (phi_ij[..., 0, 0] + phi_ij[..., 1, 1] + phi_ij[..., 2, 2])/2
+            plt.scatter(kk[1:-2], energy_spec[1:-2], label=f'n={n}, r={run}')
+        plt.xscale('log')
+        plt.yscale('log')
+        plt.legend()
+        plt.savefig('spectra_convergence.pdf')
+                 
 
 def calculate_constants(*, epsilon, Lint, etaK):
     sol = root(optimize_func,
             x0 = (6.78, 0.40),
-            args = (epsilon, Lint, etaK))
+            args = (epsilon, Lint, etaK),
+            method='hybr',
+            options={'eps':0.001, 'factor':0.1})
     assert sol.success
     return sol.x
 
-def optimize_func(c_L, c_eta, diss, Lint, etaK):
+def optimize_func(constants, diss, Lint, etaK):
     energy = (Lint*diss)**(2./3.)
     energy_model = quad(
-        lambda k : E(k, epsilon=diss, Lint=Lint, etaK=etaK, c_L=c_L, c_eta=c_eta),
+        lambda k : E(k, epsilon=diss, Lint=Lint, etaK=etaK, c_L=constants[0], c_eta=constants[1]),
         0, np.inf)[0]
     nu = (etaK**4*diss)**(1./3.)
     diss_model = quad(
-        lambda k : 2*nu*k**2*E(k, epsilon=epsilon, Lint=Lint, etaK=etaK, c_L=c_L, c_eta=c_eta),
+        lambda k : 2*nu*k**2*E(k, epsilon=diss, Lint=Lint, etaK=etaK, c_L=constants[0], c_eta=constants[1]),
         0, np.inf)[0]
     return (energy_model - energy, diss_model - diss)