diff --git a/TurTLE/test/test_Gaussian_field.py b/TurTLE/test/test_Gaussian_field.py
index 3c8646ea1f6f2eb0d649f5fab7bd878c4be59a6d..41f3cf31efa8783a561cbcf9623fbb41883f971b 100644
--- a/TurTLE/test/test_Gaussian_field.py
+++ b/TurTLE/test/test_Gaussian_field.py
@@ -4,6 +4,7 @@ import numpy as np
 from scipy import trapz
 from scipy.stats import norm
 from scipy.integrate import quad
+from scipy.optimize import root
 import h5py
 import sys, os
 import time
@@ -16,35 +17,15 @@ try:
 except:
     plt = None
 
-def E(k, *, epsilon, C=1.5, Lint, etaK, c_L=6.78, c_eta=0.4, beta=5.2):
-    return C*epsilon**(2./3.)*k**(-5./3.)*f_L(k*Lint, c_L=c_L)*f_eta(k*etaK, beta=beta, c_eta=c_eta)
-
-def calculate_constants(*, epsilon, Lint, etaK):
-    Ee_optimize1 = lambda c_L : energy - quad(
-            lambda k : E(k, epsilon=epsilon, Lint=Lint, etaK=etaK, c_L=c_L),
-            0, np.inf)[0]
-    sol = root_scalar(Ee_optimize1, x0 = 6.78, bracket=(6., 10.))
-    assert sol.converged
-    c_L = sol.root
-    print('New c_L: {}'.format(c_L))
-    Ee_optimize2 = lambda c_eta : epsilon - quad(
-            lambda k : 2*nu*k**2*E(k, epsilon=epsilon, Lint=Lint, etaK=etaK, c_L=c_L, c_eta=c_eta),
-            0, np.inf)[0]
-    sol = root_scalar(Ee_optimize2, x0 = 0.4, bracket=(0.1, 1.))
-    assert sol.converged
-    c_eta = sol.root
-
 def main():
     # size of grid
     n = 256
     # from a 1024^3 simulation
     dissipation = 0.37127850066981344
-    Lint = 0.9029294339535114
+    Lint = 0.9029294339535114 # probably the wrong (inconsistent with Pope) definition of Lint
     etaK = 0.003730966576882254
 
-    # matching the simulation's energy and dissipation
-    c_L = 6.031276834972215
-    c_eta = 0.40331391388620164
+    c_L, c_eta = calculate_constants(epsilon=dissipation, Lint=Lint, etaK=etaK)
 
     bin_no = 100
     rseed = int(time.time())
@@ -75,6 +56,70 @@ def main():
     plot_stuff(simname)
     return None
 
+def main_convergence_test():
+    n_arr = [64, 128, 256, 512]
+    diss = 
+    energy = 
+    etaK = 
+    Lint = energy**(3./2.)/diss
+
+    c_L, c_eta = calculate_constants(epsilon=diss, Lint=Lint, etaK=etaK)
+
+    bin_no = 100
+
+    if not os.path.exists(simname+'.h5'):
+        for n in n_arr:
+            for run in range(3):
+                rseed = int(time.time())
+                simname = f'conv_test_n{n}_run{run}'
+
+                c = TEST()
+                opt = c.launch(
+                    ['Gauss_field_test',
+                     '--nx', str(n),
+                     '--ny', str(n),
+                     '--nz', str(n),
+                     '--simname', simname,
+                     '--np', '4',
+                     '--environment', 'short',
+                     '--minutes', '60',
+                     '--ntpp', '1',
+                     '--wd', './',
+                     '--histogram_bins', str(bin_no),
+                     '--max_velocity_estimate', '8.',
+                     '--spectrum_dissipation', str(diss),
+                     '--spectrum_Lint', str(Lint),
+                     '--spectrum_etaK', str(etaK),
+                     '--spectrum_large_scale_const', str(c_L),
+                     '--spectrum_small_scale_const', str(c_eta),
+                     '--field_random_seed', str(rseed)] +
+                     sys.argv[1:])
+    else:
+        plt.figure()
+        # TODO
+        for n in n_arr:
+            for run in range(3):
+                df = h5py.File(simname + '.h5', 'r')
+        
+
+def calculate_constants(*, epsilon, Lint, etaK):
+    sol = root(optimize_func,
+            x0 = (6.78, 0.40),
+            args = (epsilon, Lint, etaK))
+    assert sol.success
+    return sol.x
+
+def optimize_func(c_L, c_eta, 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),
+        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),
+        0, np.inf)[0]
+    return (energy_model - energy, diss_model - diss)
+
 def E(k, *, epsilon, C=1.5, Lint, etaK, c_L=6.78, c_eta=0.4, beta=5.2):
     return C*epsilon**(2./3.)*k**(-5./3.)*f_L(k*Lint, c_L=c_L)*f_eta(k*etaK, beta=beta, c_eta=c_eta)
 
@@ -160,5 +205,5 @@ def plot_stuff(simname):
     return None
 
 if __name__ == '__main__':
-    main()
+    main_convergence_test()