diff --git a/TurTLE/DNS.py b/TurTLE/DNS.py
index 1c83e8a4bd6bb84302892fd36092cb6b60b313c7..3aea17babb57b56e02f9bdb5a3635caff925cbf5 100644
--- a/TurTLE/DNS.py
+++ b/TurTLE/DNS.py
@@ -725,9 +725,11 @@ class DNS(_code):
         if dns_type == 'kraichnan_field':
             pars['output_velocity'] = int(1)
             pars['field_random_seed'] = int(1)
-            pars['spectrum_slope'] = float(-5./3)
-            pars['spectrum_k_cutoff'] = float(16)
-            pars['spectrum_coefficient'] = float(0.1)
+            pars['spectrum_dissipation'] = float(0.027)
+            pars['spectrum_Lint'] = float(1.3)
+            pars['spectrum_etaK'] = float(0.3)
+            pars['spectrum_large_scale_const'] = float(6.78)
+            pars['spectrum_small_scale_const'] = float(0.40)
         if dns_type == 'NSVE_Stokes_particles':
             pars['initial_field_amplitude'] = float(0.0)
             pars['initial_particle_vel'] = float(0.05)
diff --git a/TurTLE/TEST.py b/TurTLE/TEST.py
index 110e51e90cbadde1bd6ac7eca0a4564b95c8a310..27cf3ad07eb721a402c2ab9dc435381d83457f38 100644
--- a/TurTLE/TEST.py
+++ b/TurTLE/TEST.py
@@ -133,9 +133,11 @@ class TEST(_code):
         if dns_type == 'Gauss_field_test':
             pars['histogram_bins'] = int(129)
             pars['max_velocity_estimate'] = float(8)
-            pars['spectrum_slope'] = float(-5./3)
-            pars['spectrum_k_cutoff'] = float(16)
-            pars['spectrum_coefficient'] = float(1)
+            pars['spectrum_dissipation'] = float(0.027)
+            pars['spectrum_Lint'] = float(1.3)
+            pars['spectrum_etaK'] = float(0.3)
+            pars['spectrum_large_scale_const'] = float(6.78)
+            pars['spectrum_small_scale_const'] = float(0.40)
         return pars
     def get_kspace(self):
         kspace = {}
diff --git a/TurTLE/test/test_Gaussian_field.py b/TurTLE/test/test_Gaussian_field.py
index 00b1692e084cbb346a3a19d3f046cdba5cb1d01b..6238aadcde14fb26844d3a6f76f2d5b097beac07 100644
--- a/TurTLE/test/test_Gaussian_field.py
+++ b/TurTLE/test/test_Gaussian_field.py
@@ -16,14 +16,36 @@ 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 = 1024
-    slope = -5./3.
-    k_cutoff = 64.
-    func = lambda k, k_c=k_cutoff, s=slope : k**s*np.exp(-k/k_c)
-    total_energy = quad(func, 0.6, k_cutoff*8)[0]
-    coeff = 1./total_energy
+    # from a 1024^3 simulation
+    dissipation = 0.37127850066981344
+    Lint = 0.9029294339535114
+    etaK = 0.003730966576882254
+
+    # matching the simulation's energy and dissipation
+    c_L = 6.031276834972215
+    c_eta = 0.40331391388620164
+
     bin_no = 100
     rseed = int(time.time())
     simname = 'Gaussianity_test'
@@ -43,25 +65,40 @@ def main():
              '--wd', './',
              '--histogram_bins', str(bin_no),
              '--max_velocity_estimate', '8.',
-             '--spectrum_slope', str(slope),
-             '--spectrum_k_cutoff', str(k_cutoff),
-             '--spectrum_coefficient', str(coeff),
+             '--spectrum_dissipation', str(dissipation),
+             '--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:])
     plot_stuff(simname, total_energy = total_energy)
     return 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 f_L(x, *, c_L):
+    return (x/(x**2 + c_L)**(1./2.))**(11./3.)
+
+def f_eta(x, *, beta, c_eta):
+    return np.exp(-beta*((x**4 + c_eta**4)**(1./4.) - c_eta))
+
 def plot_stuff(simname, total_energy = 1.):
     df = h5py.File(simname + '.h5', 'r')
-    for kk in ['spectrum_slope',
-               'spectrum_k_cutoff',
-               'spectrum_coefficient',
+    for kk in ['spectrum_dissipation',
+               'spectrum_Lint',
+               'spectrum_etaK',
+               'spectrum_large_scale_const',
+               'spectrum_small_scale_const',
                'field_random_seed',
                'histogram_bins']:
         print(kk, df['parameters/' + kk][...])
-    coeff = df['parameters/spectrum_coefficient'][()]
-    k_cutoff = df['parameters/spectrum_k_cutoff'][()]
-    slope = df['parameters/spectrum_slope'][()]
+    dissipation = df['parameters/spectrum_dissipation'][()]
+    Lint = df['parameters/spectrum_Lint'][()]
+    etaK = df['parameters/spectrum_etak'][()]
+    c_L = df['parameters/spectrum_large_scale_const'][()]
+    c_eta = df['parameters/spectrum_small_scale_const'][()]
     bin_no = df['parameters/histogram_bins'][()]
 
     f = plt.figure()
@@ -73,7 +110,10 @@ def plot_stuff(simname, total_energy = 1.):
     energy = (phi_ij[..., 0, 0] + phi_ij[..., 1, 1] + phi_ij[..., 2, 2])/2
 
     a.scatter(kk[1:-2], energy[1:-2])
-    a.plot(kk[1:-2], coeff*kk[1:-2]**slope*np.exp(-kk[1:-2]/k_cutoff), ls='--', c='C0')
+    a.plot(kk[1:-2],
+            E(kk[1:-2],
+                epsilon=dissipation, Lint=Lint, etaK=etaK, c_eta=c_eta, c_L=c_L),
+            ls='--', c='C0')
     a.set_xscale('log')
     a.set_yscale('log')
     # test isotropy
@@ -85,7 +125,6 @@ def plot_stuff(simname, total_energy = 1.):
     hist_vel = df['statistics/histograms/velocity'][0, :, :3]
     f_vel = hist_vel / np.sum(hist_vel, axis=0, keepdims=True).astype(float) / velbinsize
 
-    print('Energy analytically: {}'.format(total_energy))
     print(np.sum(energy*np.arange(len(energy))**2))
     print('Energy sum: {}'.format(np.sum(energy*df['kspace/dk'][()])))
     print('Moment sum: {}'.format(df['statistics/moments/velocity'][0,2,3]/2))
@@ -108,9 +147,9 @@ def plot_stuff(simname, total_energy = 1.):
             df['statistics/moments/velocity_gradient'][0, 2].mean()))
 
     print('----------- k2-premultiplied spectrum -----------')
-    k2func = lambda k, k_c=k_cutoff, s=slope : k**(2+s)*np.exp(-k/k_c)
-    k2sum_analytic = quad(k2func, 0, k_cutoff*20)[0]
-    print('Analytically: {}'.format(k2sum_analytic))
+    #k2func = lambda k, k_c=k_cutoff, s=slope : k**(2+s)*np.exp(-k/k_c)
+    #k2sum_analytic = quad(k2func, 0, k_cutoff*20)[0]
+    #print('Analytically: {}'.format(k2sum_analytic))
     k2spec_trace = (
             df['statistics/spectra/k*velocity_k*velocity'][..., 0, 0]
             + df['statistics/spectra/k*velocity_k*velocity'][..., 1, 1]
diff --git a/cpp/full_code/Gauss_field_test.cpp b/cpp/full_code/Gauss_field_test.cpp
index 17f2c2023dd858cc657328cdc56f515b06452a71..8f4c907bb8255ff1b0c5838b9a06f70e2e9f3e16 100644
--- a/cpp/full_code/Gauss_field_test.cpp
+++ b/cpp/full_code/Gauss_field_test.cpp
@@ -82,9 +82,11 @@ int Gauss_field_test<rnumber>::read_parameters()
             H5F_ACC_RDONLY,
             H5P_DEFAULT);
     this->max_velocity_estimate = hdf5_tools::read_value<double>(parameter_file, "/parameters/max_velocity_estimate");
-    this->spectrum_slope = hdf5_tools::read_value<double>(parameter_file, "/parameters/spectrum_slope");
-    this->spectrum_k_cutoff = hdf5_tools::read_value<double>(parameter_file, "/parameters/spectrum_k_cutoff");
-    this->spectrum_coefficient = hdf5_tools::read_value<double>(parameter_file, "/parameters/spectrum_coefficient");
+    this->spectrum_dissipation = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_dissipation");
+    this->spectrum_Lint = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_Lint");
+    this->spectrum_etaK = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_etaK");
+    this->spectrum_large_scale_const = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_large_scale_const");
+    this->spectrum_small_scale_const = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_small_scale_const");
     this->random_seed = hdf5_tools::read_value<int>(parameter_file, "/parameters/field_random_seed");
     H5Fclose(parameter_file);
     return EXIT_SUCCESS;
@@ -97,15 +99,18 @@ int Gauss_field_test<rnumber>::do_work(void)
 
     /// initialize vector Gaussian random field
     make_gaussian_random_field(
-            this->kk,
-            this->vector_field,
-            this->random_seed,
-            this->spectrum_slope,
-            this->spectrum_k_cutoff,
-            this->spectrum_coefficient);
+        this->kk,
+        this->vector_field,
+        this->spectrum_dissipation,
+        this->spectrum_Lint,
+        this->spectrum_etaK,
+        this->spectrum_large_scale_const,
+        this->spectrum_small_scale_const,
+        3./2., //incompressibility projection factor
+        this->random_seed);
 
     /// impose divergence free condition while maintaining the energy of the field
-    this->kk->template project_divfree<rnumber>(this->vector_field->get_cdata(), true);
+    this->kk->template project_divfree<rnumber>(this->vector_field->get_cdata());
 
     /// initialize statistics file.
     hid_t stat_group, stat_file;
diff --git a/cpp/full_code/Gauss_field_test.hpp b/cpp/full_code/Gauss_field_test.hpp
index 855b19279c27da7a058b5c010ef2197308d317ef..b21df6699def682fb1ae3f62c6b14a37442cf3cb 100644
--- a/cpp/full_code/Gauss_field_test.hpp
+++ b/cpp/full_code/Gauss_field_test.hpp
@@ -46,9 +46,11 @@ class Gauss_field_test: public test
 
         /* parameters that are read in read_parameters */
         double max_velocity_estimate;
-        double spectrum_slope;
-        double spectrum_k_cutoff;
-        double spectrum_coefficient;
+        double spectrum_dissipation;
+        double spectrum_Lint;
+        double spectrum_etaK;
+        double spectrum_large_scale_const;
+        double spectrum_small_scale_const;
         int random_seed;
 
         /* other stuff */
diff --git a/cpp/full_code/kraichnan_field.cpp b/cpp/full_code/kraichnan_field.cpp
index 2c9aa904775db6e3685af471e540df894c89e578..bac9a4a5789ad362e187f5c1e8748d6040027ebd 100644
--- a/cpp/full_code/kraichnan_field.cpp
+++ b/cpp/full_code/kraichnan_field.cpp
@@ -127,8 +127,8 @@ int kraichnan_field<rnumber>::generate_random_velocity(void)
         this->spectrum_dissipation,
         this->spectrum_Lint,
         this->spectrum_etaK,
-        this->spectrum_c_L,
-        this->spectrum_c_eta,
+        this->spectrum_large_scale_const,
+        this->spectrum_small_scale_const,
         3./2., // incompressibility projection factor
         this->field_random_seed + 10*this->iteration); // not an ideal choice because resulting field sequence depends on MPI/OpenMP configuration. See note below
     // this->velocity is now in Fourier space
@@ -300,8 +300,8 @@ int kraichnan_field<rnumber>::read_parameters(void)
     this->spectrum_dissipation = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_dissipation");
     this->spectrum_Lint = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_Lint");
     this->spectrum_etaK = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_etaK");
-    this->spectrum_c_L = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_c_L");
-    this->spectrum_c_eta = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_c_eta");
+    this->spectrum_large_scale_const = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_large_scale_const");
+    this->spectrum_small_scale_const = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_small_scale_const");
     this->field_random_seed = hdf5_tools::read_value<int>(parameter_file, "parameters/field_random_seed");
     this->output_velocity = hdf5_tools::read_value<int>(parameter_file, "parameters/output_velocity");
     this->max_velocity_estimate = hdf5_tools::read_value<double>(parameter_file, "parameters/max_velocity_estimate");
diff --git a/cpp/full_code/kraichnan_field.hpp b/cpp/full_code/kraichnan_field.hpp
index d7fb215efcdfa57a6ac22c03717526839e1ab016..dc434c2d1f9d819e958960e22bb989b8de6063c1 100644
--- a/cpp/full_code/kraichnan_field.hpp
+++ b/cpp/full_code/kraichnan_field.hpp
@@ -60,9 +60,11 @@ class kraichnan_field: public direct_numerical_simulation
         /* other stuff */
         kspace<FFTW, SMOOTH> *kk;
         field<rnumber, FFTW, THREE> *velocity;
-        double spectrum_slope;
-        double spectrum_k_cutoff;
-        double spectrum_coefficient;
+        double spectrum_dissipation;
+        double spectrum_Lint;
+        double spectrum_etaK;
+        double spectrum_large_scale_const;
+        double spectrum_small_scale_const;
         double max_velocity_estimate;
         int field_random_seed;