More Gaussian field test stuff

TurTLE/test/test_Gaussian_field.py

 #! /usr/bin/env python
 
 import numpy as np
-from scipy.special import gamma
+from scipy import trapz
+from scipy.stats import norm
 import h5py
 import sys
+import time
 
 import TurTLE
 from TurTLE import TEST
@@ -16,11 +18,12 @@ except:
 def main():
     c = TEST()
     # size of grid
-    n = 64
-    slope = 1.
-    k_cutoff = 20.
-    coeff = np.sqrt(1./0.93058)
-    bin_no = 33
+    n = 128
+    slope = -5./3.
+    k_cutoff = 16.
+    coeff = 1.
+    bin_no = 129
+    rseed = int(time.time())
 
     c.launch(
             ['Gauss_field_test',
@@ -36,20 +39,23 @@ def main():
              '--spectrum_slope', str(slope),
              '--spectrum_k_cutoff', str(k_cutoff),
              '--spectrum_coefficient', str(coeff),
-             '--field_random_seed', '33'] +
+             '--field_random_seed', str(rseed)] +
              sys.argv[1:])
 
 
     df = h5py.File(c.simname + '.h5', 'r')
     f = plt.figure()
 
+    # EVERYTHING SHOULD BE READ FROM FILE BECAUSE WE CANT BE SURE THE CODE RAN AGAIN
+
     # test spectrum
     a = f.add_subplot(121)
     kk = df['kspace/kshell'][...]
-    phi_ij = df['statistics/spectra/velocity_velocity'][0]
+    print('dk: {}'.format(df['kspace/dk'][()]))
+    phi_ij = df['statistics/spectra/velocity_velocity'][0] / df['kspace/dk'][()]
     energy = (phi_ij[..., 0, 0] + phi_ij[..., 1, 1] + phi_ij[..., 2, 2])/2
     
-    a.plot(kk[1:-2], energy[1:-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.set_xscale('log')
     a.set_yscale('log')
@@ -58,8 +64,21 @@ def main():
 
     max_vel_estimate = df['parameters/max_velocity_estimate'][()]
     velbinsize = 2*max_vel_estimate/bin_no
-    np.linspace(-max_vel_estimate+velbinsize/2, max_vel_estimate-velbinsize/2, bin_no)
-    a.plot(df['statistics/histograms/velocity'][0, :, :3])
+    vel = np.linspace(-max_vel_estimate+velbinsize/2, max_vel_estimate-velbinsize/2, bin_no)
+    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 integral: {}'.format(trapz(energy[1:-2], kk[1:-2])))
+    print('Energy sum: {}'.format(np.sum(energy[1:-2])))
+    print('Moment sum: {}'.format(df['statistics/moments/velocity'][0,2,3]/2))
+    print('Velocity variances: {}'.format(trapz(vel[:,None]**2*f_vel, vel[:,None], axis=0)))
+
+    vel_variance = df['statistics/moments/velocity'][0,2,3]/3.
+    a.plot(vel[:,None]/np.sqrt(vel_variance), f_vel*np.sqrt(vel_variance))
+    a.plot(vel/np.sqrt(vel_variance), norm.pdf(vel/np.sqrt(vel_variance)), ls='--')
+    a.set_yscale('log')
+    a.set_xlim(-5,5)
+    a.set_ylim(1e-5,1)
     f.tight_layout()
     f.savefig('spectrum_isotropy_test.pdf')
     df.close()

cpp/full_code/Gauss_field_test.cpp

@@ -82,9 +82,10 @@ int make_gaussian_random_field(
                     }
                     case THREE:
                     for (int cc = 0; cc<3; cc++)
-                    {
-                        output_field->cval(cindex,cc,0) = coefficient / 2 / M_PI * rdist(rgen[omp_get_thread_num()]) * pow(k2, slope/4. - 1) * exp(- sqrt(k2)/k_cutoff/2.);
-                        output_field->cval(cindex,cc,1) = coefficient / 2 / M_PI * rdist(rgen[omp_get_thread_num()]) * pow(k2, slope/4. - 1) * exp(- sqrt(k2)/k_cutoff/2.);
+                    {   // In the slope, a factor of 1/2 compensates the k^2, a factor of 1/2 goes from the energy to the Fourier coefficient
+                        // and the 2*pi*k^2 divides out the surface of the shell
+                        output_field->cval(cindex,cc,0) = coefficient / 2 / M_PI * rdist(rgen[omp_get_thread_num()]) * pow(k2, (slope - 2)/4.) * exp(- sqrt(k2)/k_cutoff/2.);
+                        output_field->cval(cindex,cc,1) = coefficient / 2 / M_PI * rdist(rgen[omp_get_thread_num()]) * pow(k2, (slope - 2)/4.) * exp(- sqrt(k2)/k_cutoff/2.);
                     }  // TODO: ADJUST OTHER FIELD TYPES
                         break;
                     case THREExTHREE: