From fcfddcdebb89b443fb33e2d19658e1ba41ce830d Mon Sep 17 00:00:00 2001
From: Tom Sandnes <tdsandnes@gmail.com>
Date: Wed, 26 Mar 2025 10:22:45 +0000
Subject: [PATCH] Tweak plotting of Planetary  KH, RT, and 3D square test
 examples

---
 .../KelvinHelmholtz_EarthLike_3D/plotSnapshots.py |  9 +++++++++
 .../KelvinHelmholtz_IdealGas_3D/plotSnapshots.py  |  8 ++++++++
 .../plotSnapshots.py                              |  8 ++++++++
 .../plotSnapshots.py                              |  8 ++++++++
 .../plotSnapshots.py                              | 13 +++++++++++--
 .../RayleighTaylor_EarthLike_3D/plotSnapshots.py  | 11 ++++++++++-
 .../RayleighTaylor_IdealGas_3D/plotSnapshots.py   |  8 ++++++++
 .../plotSnapshots.py                              | 15 ++++++++++++---
 examples/Planetary/SquareTest_3D/plotSolution.py  | 10 ++++++++++
 9 files changed, 84 insertions(+), 6 deletions(-)

diff --git a/examples/Planetary/KelvinHelmholtz_EarthLike_3D/plotSnapshots.py b/examples/Planetary/KelvinHelmholtz_EarthLike_3D/plotSnapshots.py
index 6e33363ee3..391ab35958 100644
--- a/examples/Planetary/KelvinHelmholtz_EarthLike_3D/plotSnapshots.py
+++ b/examples/Planetary/KelvinHelmholtz_EarthLike_3D/plotSnapshots.py
@@ -78,6 +78,15 @@ def plot_kh(ax, snap, mat_id1, mat_id2, cmap1, cmap2, norm1, norm2):
         A1_m = f["/PartType0/Masses"][:] * m
         A1_mat_id = f["/PartType0/MaterialIDs"][:]
 
+    # Sort arrays based on z position
+    sort_indices = np.argsort(A1_z)
+    A1_x = A1_x[sort_indices]
+    A1_y = A1_y[sort_indices]
+    A1_z = A1_z[sort_indices]
+    A1_rho = A1_rho[sort_indices]
+    A1_m = A1_m[sort_indices]
+    A1_mat_id = A1_mat_id[sort_indices]
+
     # Mask to select slice
     slice_thickness = 0.2 * (np.max(A1_z) - np.min(A1_z))
     slice_pos_z = 0.5 * (np.max(A1_z) + np.min(A1_z))
diff --git a/examples/Planetary/KelvinHelmholtz_IdealGas_3D/plotSnapshots.py b/examples/Planetary/KelvinHelmholtz_IdealGas_3D/plotSnapshots.py
index 95b26e55d6..85df472ec9 100644
--- a/examples/Planetary/KelvinHelmholtz_IdealGas_3D/plotSnapshots.py
+++ b/examples/Planetary/KelvinHelmholtz_IdealGas_3D/plotSnapshots.py
@@ -70,6 +70,14 @@ def plot_kh(ax, snap, cmap, norm):
         A1_rho = f["/PartType0/Densities"][:]
         A1_m = f["/PartType0/Masses"][:]
 
+    # Sort arrays based on z position
+    sort_indices = np.argsort(A1_z)
+    A1_x = A1_x[sort_indices]
+    A1_y = A1_y[sort_indices]
+    A1_z = A1_z[sort_indices]
+    A1_rho = A1_rho[sort_indices]
+    A1_m = A1_m[sort_indices]
+
     # Mask to select slice
     slice_thickness = 0.2 * (np.max(A1_z) - np.min(A1_z))
     slice_pos_z = 0.5 * (np.max(A1_z) + np.min(A1_z))
diff --git a/examples/Planetary/KelvinHelmholtz_IdealGas_LargeDensityContrast_3D/plotSnapshots.py b/examples/Planetary/KelvinHelmholtz_IdealGas_LargeDensityContrast_3D/plotSnapshots.py
index ad1ec66201..2f5d1404d7 100644
--- a/examples/Planetary/KelvinHelmholtz_IdealGas_LargeDensityContrast_3D/plotSnapshots.py
+++ b/examples/Planetary/KelvinHelmholtz_IdealGas_LargeDensityContrast_3D/plotSnapshots.py
@@ -70,6 +70,14 @@ def plot_kh(ax, snap, cmap, norm):
         A1_rho = f["/PartType0/Densities"][:]
         A1_m = f["/PartType0/Masses"][:]
 
+    # Sort arrays based on z position
+    sort_indices = np.argsort(A1_z)
+    A1_x = A1_x[sort_indices]
+    A1_y = A1_y[sort_indices]
+    A1_z = A1_z[sort_indices]
+    A1_rho = A1_rho[sort_indices]
+    A1_m = A1_m[sort_indices]
+
     # Mask to select slice
     slice_thickness = 0.2 * (np.max(A1_z) - np.min(A1_z))
     slice_pos_z = 0.5 * (np.max(A1_z) + np.min(A1_z))
diff --git a/examples/Planetary/KelvinHelmholtz_IdealGas_Smooth_3D/plotSnapshots.py b/examples/Planetary/KelvinHelmholtz_IdealGas_Smooth_3D/plotSnapshots.py
index d6f024d4e2..0c4394ecef 100644
--- a/examples/Planetary/KelvinHelmholtz_IdealGas_Smooth_3D/plotSnapshots.py
+++ b/examples/Planetary/KelvinHelmholtz_IdealGas_Smooth_3D/plotSnapshots.py
@@ -70,6 +70,14 @@ def plot_kh(ax, snap, cmap, norm):
         A1_rho = f["/PartType0/Densities"][:]
         A1_m = f["/PartType0/Masses"][:]
 
+    # Sort arrays based on z position
+    sort_indices = np.argsort(A1_z)
+    A1_x = A1_x[sort_indices]
+    A1_y = A1_y[sort_indices]
+    A1_z = A1_z[sort_indices]
+    A1_rho = A1_rho[sort_indices]
+    A1_m = A1_m[sort_indices]
+
     # Mask to select slice
     slice_thickness = 0.2 * (np.max(A1_z) - np.min(A1_z))
     slice_pos_z = 0.5 * (np.max(A1_z) + np.min(A1_z))
diff --git a/examples/Planetary/KelvinHelmholtz_JupiterLike_3D/plotSnapshots.py b/examples/Planetary/KelvinHelmholtz_JupiterLike_3D/plotSnapshots.py
index cc7dbe496d..18df751c10 100644
--- a/examples/Planetary/KelvinHelmholtz_JupiterLike_3D/plotSnapshots.py
+++ b/examples/Planetary/KelvinHelmholtz_JupiterLike_3D/plotSnapshots.py
@@ -78,6 +78,15 @@ def plot_kh(ax, snap, mat_id1, mat_id2, cmap1, cmap2, norm1, norm2):
         A1_m = f["/PartType0/Masses"][:] * m
         A1_mat_id = f["/PartType0/MaterialIDs"][:]
 
+    # Sort arrays based on z position
+    sort_indices = np.argsort(A1_z)
+    A1_x = A1_x[sort_indices]
+    A1_y = A1_y[sort_indices]
+    A1_z = A1_z[sort_indices]
+    A1_rho = A1_rho[sort_indices]
+    A1_m = A1_m[sort_indices]
+    A1_mat_id = A1_mat_id[sort_indices]
+
     # Mask to select slice
     slice_thickness = 0.2 * (np.max(A1_z) - np.min(A1_z))
     slice_pos_z = 0.5 * (np.max(A1_z) + np.min(A1_z))
@@ -132,13 +141,13 @@ def plot_kh(ax, snap, mat_id1, mat_id2, cmap1, cmap2, norm1, norm2):
 if __name__ == "__main__":
 
     # Set colormap
-    cmap1 = plt.get_cmap('Purples')
+    cmap1 = plt.get_cmap('winter_r')
     mat_id1 = 304
     rho_min1 = 6800
     rho_max1 = 8700
     norm1 = mpl.colors.Normalize(vmin=rho_min1,vmax=rho_max1)
 
-    cmap2 = plt.get_cmap('Oranges_r')
+    cmap2 = plt.get_cmap('autumn')
     mat_id2 = 307
     rho_min2 = 3450
     rho_max2 = 4050
diff --git a/examples/Planetary/RayleighTaylor_EarthLike_3D/plotSnapshots.py b/examples/Planetary/RayleighTaylor_EarthLike_3D/plotSnapshots.py
index 692c22c881..edb29f98e1 100644
--- a/examples/Planetary/RayleighTaylor_EarthLike_3D/plotSnapshots.py
+++ b/examples/Planetary/RayleighTaylor_EarthLike_3D/plotSnapshots.py
@@ -81,6 +81,15 @@ def plot_kh(ax, snap, mat_id1, mat_id2, cmap1, cmap2, norm1, norm2):
         A1_m = f["/PartType0/Masses"][:] * m
         A1_mat_id = f["/PartType0/MaterialIDs"][:]
 
+    # Sort arrays based on z position
+    sort_indices = np.argsort(A1_z)
+    A1_x = A1_x[sort_indices]
+    A1_y = A1_y[sort_indices]
+    A1_z = A1_z[sort_indices]
+    A1_rho = A1_rho[sort_indices]
+    A1_m = A1_m[sort_indices]
+    A1_mat_id = A1_mat_id[sort_indices]
+
     # Mask to select slice
     slice_thickness = 0.1 * (np.max(A1_z) - np.min(A1_z))
     slice_pos_z = 0.5 * (np.max(A1_z) + np.min(A1_z))
@@ -160,7 +169,7 @@ if __name__ == "__main__":
         time = times[i]
 
         plot_kh(ax, snap, mat_id1, mat_id2, cmap1, cmap2, norm1, norm2)
-        ax.text(0.5,-0.05,r"$t =\;$"+time+"$\,$s", horizontalalignment='center',size=18, transform=ax.transAxes)
+        ax.text(0.5,-0.05,r"$t =\;$"+time+r"$\,$s", horizontalalignment='center',size=18, transform=ax.transAxes)
 
     # Colour bar
     sm1 = plt.cm.ScalarMappable(cmap=cmap1, norm=norm1)
diff --git a/examples/Planetary/RayleighTaylor_IdealGas_3D/plotSnapshots.py b/examples/Planetary/RayleighTaylor_IdealGas_3D/plotSnapshots.py
index ff2f5c2d62..99787dea94 100644
--- a/examples/Planetary/RayleighTaylor_IdealGas_3D/plotSnapshots.py
+++ b/examples/Planetary/RayleighTaylor_IdealGas_3D/plotSnapshots.py
@@ -73,6 +73,14 @@ def plot_kh(ax, snap, cmap, norm):
         A1_rho = f["/PartType0/Densities"][:]
         A1_m = f["/PartType0/Masses"][:]
 
+    # Sort arrays based on z position
+    sort_indices = np.argsort(A1_z)
+    A1_x = A1_x[sort_indices]
+    A1_y = A1_y[sort_indices]
+    A1_z = A1_z[sort_indices]
+    A1_rho = A1_rho[sort_indices]
+    A1_m = A1_m[sort_indices]
+
     # Mask to select slice
     slice_thickness = 0.1 * (np.max(A1_z) - np.min(A1_z))
     slice_pos_z = 0.5 * (np.max(A1_z) + np.min(A1_z))
diff --git a/examples/Planetary/RayleighTaylor_JupiterLike_3D/plotSnapshots.py b/examples/Planetary/RayleighTaylor_JupiterLike_3D/plotSnapshots.py
index 23b48fbfff..92ec546b39 100644
--- a/examples/Planetary/RayleighTaylor_JupiterLike_3D/plotSnapshots.py
+++ b/examples/Planetary/RayleighTaylor_JupiterLike_3D/plotSnapshots.py
@@ -81,6 +81,15 @@ def plot_kh(ax, snap, mat_id1, mat_id2, cmap1, cmap2, norm1, norm2):
         A1_m = f["/PartType0/Masses"][:] * m
         A1_mat_id = f["/PartType0/MaterialIDs"][:]
 
+    # Sort arrays based on z position
+    sort_indices = np.argsort(A1_z)
+    A1_x = A1_x[sort_indices]
+    A1_y = A1_y[sort_indices]
+    A1_z = A1_z[sort_indices]
+    A1_rho = A1_rho[sort_indices]
+    A1_m = A1_m[sort_indices]
+    A1_mat_id = A1_mat_id[sort_indices]
+
     # Mask to select slice
     slice_thickness = 0.1 * (np.max(A1_z) - np.min(A1_z))
     slice_pos_z = 0.5 * (np.max(A1_z) + np.min(A1_z))
@@ -135,13 +144,13 @@ def plot_kh(ax, snap, mat_id1, mat_id2, cmap1, cmap2, norm1, norm2):
 if __name__ == "__main__":
 
     # Set colormap
-    cmap1 = plt.get_cmap('Purples')
+    cmap1 = plt.get_cmap('winter_r')
     mat_id1 = 304
     rho_min1 = 6800
     rho_max1 = 8700
     norm1 = mpl.colors.Normalize(vmin=rho_min1,vmax=rho_max1)
 
-    cmap2 = plt.get_cmap('Oranges_r')
+    cmap2 = plt.get_cmap('autumn')
     mat_id2 = 307
     rho_min2 = 3450
     rho_max2 = 4050
@@ -160,7 +169,7 @@ if __name__ == "__main__":
         time = times[i]
 
         plot_kh(ax, snap, mat_id1, mat_id2, cmap1, cmap2, norm1, norm2)
-        ax.text(0.5,-0.05,r"$t =\;$"+time+"$\,$s", horizontalalignment='center',size=18, transform=ax.transAxes)
+        ax.text(0.5,-0.05,r"$t =\;$"+time+r"$\,$s", horizontalalignment='center',size=18, transform=ax.transAxes)
 
     # Colour bar
     sm1 = plt.cm.ScalarMappable(cmap=cmap1, norm=norm1)
diff --git a/examples/Planetary/SquareTest_3D/plotSolution.py b/examples/Planetary/SquareTest_3D/plotSolution.py
index a758c5b2f7..f1010413b7 100644
--- a/examples/Planetary/SquareTest_3D/plotSolution.py
+++ b/examples/Planetary/SquareTest_3D/plotSolution.py
@@ -119,6 +119,16 @@ if __name__ == "__main__":
         A1_P = f["/PartType0/Pressures"][:]
         A1_m = f["/PartType0/Masses"][:]
 
+    # Sort arrays based on z position
+    sort_indices = np.argsort(A1_z)
+    A1_x = A1_x[sort_indices]
+    A1_y = A1_y[sort_indices]
+    A1_z = A1_z[sort_indices]
+    A1_rho = A1_rho[sort_indices]
+    A1_u = A1_u[sort_indices]
+    A1_P = A1_P[sort_indices]
+    A1_m = A1_m[sort_indices]
+
     # Mask to select slice
     slice_thickness = 0.1
     slice_pos_z = 0.5 * (np.max(A1_z) + np.min(A1_z))
-- 
GitLab