Merge pull request #2 from mreinacampos/master

Added an input parameter to change the slope of the turbulent power spectrum

MakeCloud.py

@@ -13,8 +13,8 @@
    --density_exponent=<f>   Power law exponent of the density profile [default: 0.0]
    --spin=<f>           Spin parameter: fraction of binding energy in solid-body rotation [default: 0.0]
    --omega_exponent=<f>  Powerlaw exponent of rotational frequency as a function of cylindrical radius [default: 0.0]
-   --turb_type=<s>      Type of initial turbulent velocity (and possibly density field): 'gaussian' or 'full' [default: gaussian]
-   --turb_sol=<f>       Fraction of turbulence in solenoidal modes, used when turb_type is 'gaussian' [default: 0.5]
+   --turb_slope=<f>      Slope of the turbulent power spectra [default: 2.0]
+   --turb_sol=<f>       Fraction of turbulence in solenoidal modes [default: 0.5]
    --alpha_turb=<f>     Turbulent virial parameter (BM92 convention: 2Eturb/|Egrav|) [default: 2.]
    --bturb=<f>          Magnetic energy as a fraction of the binding energy [default: 0.1]
    --bfixed=<f>         Magnetic field in magnitude in code units, used instead of bturb if not set to zero [default: 0]
@@ -84,7 +84,7 @@ def get_glass_coords(N_gas, glass_path):
     return x
 
 
-def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42):
+def TurbField(res=256, minmode=2, maxmode=64, slope=2.0, sol_weight=1.0, seed=42):
     freqs = fftpack.fftfreq(res)
     freq3d = np.array(np.meshgrid(freqs, freqs, freqs, indexing="ij"))
     intfreq = np.around(freq3d * res)
@@ -98,7 +98,7 @@ def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42):
         rand_phase = fftpack.fftn(
             np.random.normal(size=kSqr.shape)
         )  # fourier transform of white noise
-        vk = rand_phase * (float(minmode) / res) ** 2 / (kSqr + 1e-300)
+        vk = rand_phase * (float(minmode) / res) ** 2 / (np.power(kSqr, slope/2.0) + 1e-300)
         vk[intkSqr == 0] = 0.0
         vk[intkSqr < minmode**2] *= (
             intkSqr[intkSqr < minmode**2] ** 2 / minmode**4
@@ -149,6 +149,7 @@ def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42):
 seed = int(float(arguments["--turb_seed"]) + 0.5)
 tmax = int(float(arguments["--tmax"]))
 nsnap = int(float(arguments["--nsnap"]))
+turb_slope = float(arguments["--turb_slope"])
 turb_sol = float(arguments["--turb_sol"])
 magnetic_field = float(arguments["--bturb"])
 bfixed = float(arguments["--bfixed"])
@@ -215,7 +216,7 @@ def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42):
     "M%3.2g_" % (M_gas)
     + ("Mstar%g_" % (M_star) if M_star > 0 else "")
     + ("rho_exp%g_" % (-density_exponent) if density_exponent < 0 else "")
-    + "R%g_Z%g_S%g_A%g_B%g_I%g_Res%d_n%d_sol%g"
+    + "R%g_Z%g_S%g_A%g_B%g_I%g_Res%d_n%d_beta%g_sol%g"
     % (
         R,
         metallicity,
@@ -225,6 +226,7 @@ def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42):
         ISRF,
         res_effective,
         minmode,
+        turb_slope,
         turb_sol,
     )
     + ("_%d" % seed)
@@ -403,9 +405,9 @@ def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42):
 
 if not os.path.exists(turb_path):
     os.makedirs(turb_path)
-fname = turb_path + "/vturb%d_sol%g_seed%d.npy" % (minmode, turb_sol, seed)
+fname = turb_path + "/vturb%d_beta%g_sol%g_seed%d.npy" % (minmode, turb_slope, turb_sol, seed)
 if not os.path.isfile(fname):
-    vt = TurbField(minmode=minmode, sol_weight=turb_sol, seed=seed)
+    vt = TurbField(minmode=minmode, slope = turb_slope, sol_weight=turb_sol, seed=seed)
     nmin, nmax = vt.shape[-1] // 4, 3 * vt.shape[-1] // 4
     vt = vt[
         :, nmin:nmax, nmin:nmax, nmin:nmax