# HG changeset patch
# User calpe  <miguel.calpe.linares@gmail.com>
# Date 1538738931 -7200
#      Fri Oct 05 13:28:51 2018 +0200
# Node ID 4cd6c9b9050a5d3ece2020e21f159c1c298a3a11
# Parent  14d6d55c07bd845dad26f17fd9de3763e9b00f7d
Changes scripts/ns2d.strat

diff --git a/scripts/ns2d.strat/compute_anisotropy.py b/scripts/ns2d.strat/compute_anisotropy.py
--- a/scripts/ns2d.strat/compute_anisotropy.py
+++ b/scripts/ns2d.strat/compute_anisotropy.py
@@ -11,6 +11,8 @@
 import h5py
 import numpy as np
 
+from fluidsim import load_params_simul
+
 def _compute_array_times_from_path(path_simulation):
     """
     Compute array with times from path simulation.
@@ -55,28 +57,35 @@
 
     print("Compute anisotropy nx = {} and gamma {}..".format(res_out, gamma_out))
 
-    # Compute index start average time imin.
-    times = _compute_array_times_from_path(path_simulation)
-    dt_state_phys = np.median(np.diff(times))
+    # Load data energy spectra file.
+    with h5py.File(path_simulation + "/spectra1D.h5", "r") as f:
+        kx = f["kxE"].value
+        times = f["times"].value
+        spectrumkx_EK_ux = f["spectrum1Dkx_EK_ux"].value
+        spectrumkx_EK = f["spectrum1Dkx_EK"].value
+        spectrumkx_EK_uy = f["spectrum1Dkx_EK_uy"].value
+        dt_state_phys = np.median(np.diff(times))
 
+    # Compute start time average itmin
+    dt = np.median(np.diff(times))
     if not tmin:
         nb_files = 10
-        tmin = np.max(times) - (nb_files * dt_state_phys)
-    imin = np.argmin(abs(times - tmin))
+        tmin = np.max(times) - (nb_files * dt)
+    itmin = np.argmin(abs(times - tmin))
 
-    # Compute anisotropy
-    anisotropies = []
-    for path in glob(path_simulation + "/state_phys_t*")[imin:]:
-        with h5py.File(path, "r") as f:
-            ux = f["state_phys"]["ux"].value
-            uz = f["state_phys"]["uy"].value
-        anisotropies.append(np.mean(ux**2 / (ux**2 + uz**2)))
+    # Compute delta_kx
+    params = load_params_simul(path_simulation)
+    delta_kx = 2 * np.pi / params.oper.Lx
 
-    return np.mean(anisotropies)
+    # Compute spatial averaged energy from spectra
+    EK_ux = np.sum(np.mean(spectrumkx_EK_ux[itmin:,:], axis=0) * delta_kx)
+    EK = np.sum(np.mean(spectrumkx_EK[itmin:,:], axis=0) * delta_kx)
+
+    return EK_ux / EK
 
 if __name__ == "__main__":
-   path_simulation = ("/fsnet/project/meige/2015/15DELDUCA/DataSim/" +
-                      "sim1920_no_shear_modes/NS2D.strat_1920x480_S2pix1.571_F07_gamma1_2018-08-14_10-01-22")
+    path_simulation = ("/fsnet/project/meige/2015/15DELDUCA/DataSim/" +
+                    "sim1920_no_shear_modes/NS2D.strat_1920x480_S2pix1.571_F07_gamma05_2018-08-14_10-01-22")
 
-   anisotropy = compute_anisotropy(path_simulation)
-   print("anisotropy = {}".format(anisotropy))
+    anisotropy = compute_anisotropy(path_simulation)
+    print("anisotropy = {}".format(anisotropy))
diff --git a/scripts/ns2d.strat/plot_reynolds_froude.py b/scripts/ns2d.strat/plot_reynolds_froude.py
--- a/scripts/ns2d.strat/plot_reynolds_froude.py
+++ b/scripts/ns2d.strat/plot_reynolds_froude.py
@@ -17,10 +17,19 @@
 from compute_reynolds_froude import compute_buoyancy_reynolds
 from fluiddyn.output.rcparams import set_rcparams
 
+def _get_resolution_from_dir(path_simulation):
+    return path_simulation.split("NS2D.strat_")[1].split("x")[0]
 
 # Create path simulations
 path_root = "/fsnet/project/meige/2015/15DELDUCA/DataSim"
-directories = ["sim1920_no_shear_modes", "sim1920_modif_res_no_shear_modes"]
+directories = ["sim960_no_shear_modes",
+               "sim1920_no_shear_modes",
+               "sim1920_modif_res_no_shear_modes",
+               "sim3840_modif_res_no_shear_modes"]
+
+# directories = ["sim960_no_shear_modes",
+#                "sim1920_no_shear_modes"]
+
 paths_simulations = []
 for directory in directories:
     paths_simulations += sorted(glob(os.path.join(path_root, directory, "NS2D*")))
@@ -29,6 +38,7 @@
 reynoldsb = []
 anisotropies = []
 dissipations = []
+markers = []
 
 set_rcparams(fontsize=14, for_article=True)
 
@@ -39,23 +49,57 @@
 ax.set_yscale("log")
 fig.text(0.8, 4e-7, r"$\frac{D(k_{fx})}{D(k_x)}$", fontsize=16)
 
+
 for path in paths_simulations:
     F_h, Re_8, R_b = compute_buoyancy_reynolds(path)
     anisotropy = compute_anisotropy(path)
     dissipation = compute_ratio_dissipation(path)
+    res = _get_resolution_from_dir(path)
 
     froudes.append(F_h)
     reynoldsb.append(R_b)
     anisotropies.append(anisotropy)
     dissipations.append(dissipation)
 
+    if res == "960":
+        markers.append("o")
+    elif res == "1920":
+        markers.append("s")
+    elif res == "3840":
+        markers.append("^")
+
     print("F_h", F_h)
     print("Re_8", Re_8)
     print("R_b", R_b)
 
-areas = 500 * np.asarray(anisotropies)**2
-scatter = ax.scatter(froudes, reynoldsb, s=areas, c=dissipations, alpha=0.7)
-ax.scatter(0.7, 1e-4, s=500 * np.asarray(0.5)**2, c="red")
-ax.text(0.64, 1e-5, "isotropy", fontsize=14, color="r")
+for _f, _r, _a, _d, _m in zip(froudes, reynoldsb, anisotropies, dissipations, markers):
+    scatter = ax.scatter(_f, _r, s=250 * (_a**2), c=_d, vmin=0, vmax=0.3, marker=_m)
+# plt.show()
+
+# areas = 250 * np.asarray(anisotropies)**2
+# scatter = ax.scatter(froudes, reynoldsb, s=areas, c=dissipations, alpha=0.7, vmin=0, vmax=0.3)
+ax.scatter(max(froudes), 1e2 * min(reynoldsb), s=250 * np.asarray(0.5)**2, c="red")
+ax.scatter(max(froudes), min(reynoldsb), s=250 * np.asarray(1.0)**2, c="red")
+ax.text(0.12, 90 * min(reynoldsb), "anisotropy=0", fontsize=12, color="r")
+ax.text(0.12, min(reynoldsb), "anisotropy=1", fontsize=12, color="r")
 fig.colorbar(scatter)
+
+# Legend
+import matplotlib.lines as mlines
+import matplotlib.pyplot as plt
+
+blue_star = mlines.Line2D([], [], color='red', marker='o', linestyle='None',
+                          markersize=8, label=r'$n_x = 960$')
+red_square = mlines.Line2D([], [], color='red', marker='s', linestyle='None',
+                          markersize=8, label=r'$n_x = 1920$')
+purple_triangle = mlines.Line2D([], [], color='red', marker='^', linestyle='None',
+                          markersize=8, label=r'$n_x = 3840$')
+
+ax.legend(handles=[blue_star, red_square, purple_triangle],
+          loc="upper center",
+          bbox_to_anchor=(0.5,1.1),
+          borderaxespad=0.,
+          ncol=len(markers),
+          handletextpad=0.1,
+          fontsize=12)
 plt.show()