diff --git a/examples/GEAR/ZoomIn/halo_distribution.py b/examples/GEAR/ZoomIn/halo_distribution.py
index 5a2ec7cce358207ec556823b2e09176a2335a842..4824c20952cf1b602f784c11aca4fa6442798b81 100644
--- a/examples/GEAR/ZoomIn/halo_distribution.py
+++ b/examples/GEAR/ZoomIn/halo_distribution.py
@@ -1,13 +1,44 @@
#!/usr/bin/env
-import yt
from yt.extensions.astro_analysis.halo_analysis.api import HaloCatalog
+import matplotlib.pyplot as plt
+import numpy as np
+
+Nbins = 20
+
+
+def savePlot(data):
+ plt.figure()
+ plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
+ markers = ["s", "o"]
+ d_min = None
+ d_max = None
+ for d in data[0]:
+ if d_min is None or d.min() < d_min:
+ d_min = d.min()
+ if d_max is None or d.max() > d_max:
+ d_max = d.max()
+
+ for i, d in enumerate(data[0]):
+ plt.hist(d, Nbins, histtype="step", log=True, range=[d_min, d_max])
+
+ plt.xlabel("$\mathrm{Mass\ \mathrm{(M_{\odot})}}$", fontsize='large')
+ plt.ylabel(r"$\mathrm{Number}\ o\mathrm{f}\ \mathrm{Halos}$", fontsize='large')
+ plt.legend(data[1], loc=4, frameon=True, ncol=2, fancybox=True)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='small')
+ plt.grid(True)
+
+ plt.savefig("halo_distribution.png")
def doPlot(f, name, i):
- hc = HaloCatalog(data_ds=f, finder_method="fof",
- finder_kwargs={"padding": 0.02})
+ hc = HaloCatalog(data_ds=f, finder_method="fof")
hc.create()
- print(hc)
- print(dir(hc))
- exit()
+
+ masses = np.zeros(len(hc.catalog))
+ for i, halo in enumerate(hc.catalog):
+ masses[i] = halo["particle_mass"].in_units("Msun")
+
+ return masses
diff --git a/examples/GEAR/ZoomIn/make_image.py b/examples/GEAR/ZoomIn/make_image.py
index d1900eb18c1f646920468598b4b97c1f2fe19b4f..73fb7d21c8e84264b3e35434cb043081ff42beba 100644
--- a/examples/GEAR/ZoomIn/make_image.py
+++ b/examples/GEAR/ZoomIn/make_image.py
@@ -6,6 +6,7 @@ import unyt
import sys
import matplotlib
import projection_plot
+import velocity_plot
import phase_plot
import halo_distribution
import add_fields
@@ -26,9 +27,10 @@ do_plot = {
"projection_density": False,
"projection_temperature": False,
"projection_mass": False,
- "halo_distribution": True,
+ "halo_distribution": False,
"phase_1d": False,
- "phase_2d": False
+ "phase_2d": False,
+ "velocity": True
}
# Generate the figures
@@ -70,9 +72,13 @@ axes = {
# Data
data = {
- "phase_1d": ([], [])
+ "phase_1d": ([], []),
+ "halo_distribution": ([], []),
+ "velocity": ([], []),
}
+names = []
+
def savePlot():
if do_plot["projection_density"]:
@@ -87,13 +93,25 @@ def savePlot():
"mass")
if do_plot["phase_1d"]:
+ data["phase_1d"][1].extend(names)
phase_plot.save1DPlot(data["phase_1d"])
if do_plot["phase_2d"]:
phase_plot.save2DPlot(figures["phase_2d"])
+ # halo distribution
+ if do_plot["halo_distribution"]:
+ data["halo_distribution"][1].extend(names)
+ halo_distribution.savePlot(data["halo_distribution"])
+
+ if do_plot["velocity"]:
+ data["velocity"][1].extend(names)
+ velocity_plot.save1DPlot(data["velocity"])
+
def doPlot(filename, i, name):
+ names.append(name)
+
f = yt.load(filename)
if (do_plot["projection_temperature"] or do_plot["phase_2d"]):
add_fields.addTemperature(f)
@@ -128,7 +146,6 @@ def doPlot(filename, i, name):
if do_plot["phase_1d"]:
p = phase_plot.do1DPlot(f, name, i)
data["phase_1d"][0].append(p)
- data["phase_1d"][1].append(name)
# 2D Phase plot
if do_plot["phase_2d"]:
@@ -137,7 +154,13 @@ def doPlot(filename, i, name):
# halo distribution
if do_plot["halo_distribution"]:
- halo_distribution.doPlot(f, name, i)
+ m = halo_distribution.doPlot(f, name, i)
+ data["halo_distribution"][0].append(m)
+
+ # Velocity plot
+ if do_plot["velocity"]:
+ p = velocity_plot.do1DPlot(f, name, i)
+ data["velocity"][0].append(p)
doPlot(swift, 0, "SWIFT")
diff --git a/examples/GEAR/ZoomIn/phase_plot.py b/examples/GEAR/ZoomIn/phase_plot.py
index 63f78e11ebce320af6988383c1e648db9f1bae46..af52e922af4cd120fbafb908b64b1fc6b9416df8 100644
--- a/examples/GEAR/ZoomIn/phase_plot.py
+++ b/examples/GEAR/ZoomIn/phase_plot.py
@@ -7,7 +7,7 @@ import matplotlib.pyplot as plt
width = 5 * kpc
limits_temperature = (1e1 * K, 1e7 * K)
-limits_density = (1e-7 * amu / cm**3, 1e7 * amu / cm**3)
+limits_density = (1e-10 * amu / cm**3, 1e3 * amu / cm**3)
limits_mass = (1e3 * Msun, 1e7 * Msun)
diff --git a/examples/GEAR/ZoomIn/projection_plot.py b/examples/GEAR/ZoomIn/projection_plot.py
index cbc4faedfbbabcceb20f79876ba8474c942bc589..40e4605d5307b249877cfeab6b98824beb6d2611 100644
--- a/examples/GEAR/ZoomIn/projection_plot.py
+++ b/examples/GEAR/ZoomIn/projection_plot.py
@@ -10,8 +10,8 @@ cmap_mass = "inferno"
center = (1441.9954833984375000,
1666.1169433593750000,
1891.6248779296875000)
-small_width = 200 * kpc
-large_width = 1000 * kpc
+small_width = 250 * kpc
+large_width = 3000 * kpc
limits_density = (1e-10 * g / cm**2, 1e-2 * g / cm**2)
limits_temperature = (10 * K, 1e5 * K)
limits_mass = None
diff --git a/examples/GEAR/ZoomIn/velocity_plot.py b/examples/GEAR/ZoomIn/velocity_plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..cb2371dc483f0eb3117a01b3c0c535dc5e1786a4
--- /dev/null
+++ b/examples/GEAR/ZoomIn/velocity_plot.py
@@ -0,0 +1,41 @@
+#!/usr/bin/env python3
+
+import yt
+from yt.units import kpc
+import matplotlib.pyplot as plt
+
+width = 1500 * kpc
+
+
+def save1DPlot(profiles):
+ plt.figure(figsize=(8, 8))
+ plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
+ markers = ["s", "o"]
+ for i, p in enumerate(profiles[0]):
+ velocity = p.x.in_units("km/s").d
+ mass = p["Masses"].in_units("Msun").d
+ plt.plot(velocity, mass, linestyle="-", marker=markers[i],
+ markeredgecolor='none', linewidth=1.2, alpha=0.8)
+ plt.semilogx()
+ plt.semilogy()
+
+ plt.xlabel("$\mathrm{Velocity\ (km/s)}$", fontsize='large')
+ plt.ylabel(r"$\mathrm{Mass,}\/\mathrm{d}M\mathrm{/dlog}\/\mathrm{\rho}\/\mathrm{(M_{\odot})}$", fontsize='large')
+ plt.legend(profiles[1], loc=4, frameon=True, ncol=2, fancybox=True)
+ leg = plt.gca().get_legend()
+ ltext = leg.get_texts()
+ plt.setp(ltext, fontsize='small')
+ plt.grid(True)
+
+ plt.savefig("velocity.png", bbox_inches='tight', pad_inches=0.03, dpi=300)
+
+
+def do1DPlot(f, name, i):
+ sp = f.sphere("max", width)
+
+ # Because ParticleProfilePlot doesn't exist, I will do the following trick.
+ p = yt.create_profile(sp, ("PartType0", "velocity_magnitude"), ("PartType0", "Masses"),
+ weight_field=None, n_bins=50,
+ accumulation=False)
+
+ return p