GEAR: update plots

examples/GEAR/ZoomIn/halo_distribution.py

 #!/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

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")

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)
 
 

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

examples/GEAR/ZoomIn/velocity_plot.py

+#!/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