diff --git a/examples/GEAR/ZoomIn/halo_distribution.py b/examples/GEAR/ZoomIn/halo_distribution.py
deleted file mode 100644
index 4824c20952cf1b602f784c11aca4fa6442798b81..0000000000000000000000000000000000000000
--- a/examples/GEAR/ZoomIn/halo_distribution.py
+++ /dev/null
@@ -1,44 +0,0 @@
-#!/usr/bin/env
-
-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")
- hc.create()
-
- 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 a87c523ba5ec35cab96aa2f227bf2ac02720442f..efee7dec18d54c7444d4abba54241fd3ee0536da 100644
--- a/examples/GEAR/ZoomIn/make_image.py
+++ b/examples/GEAR/ZoomIn/make_image.py
@@ -1,21 +1,13 @@
#!/usr/bin/env python3
import yt
-import sys
import projection_plot
-import velocity_plot
import phase_plot
-import halo_distribution
-import metal_plot
import add_fields
import star_plot
-import profile_plot
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import AxesGrid
from yt.units import kpc
-import numpy as np
-
-# Parameters
swift = "swift_final.hdf5"
gear = "h050_final.hdf5"
@@ -27,26 +19,14 @@ do_hydro = True
do_stars = True
do_feedback = True
do_plot = {
- # DMO
- "projection_mass": do_dmo,
- "velocity": do_dmo,
- "halo_distribution": True, # very slow
- "profile": do_dmo,
# hydro
"projection_density": do_hydro,
"projection_temperature": do_hydro,
- "phase_1d_density": do_hydro,
- "phase_1d_temperature": do_hydro,
"phase_2d": do_hydro,
- "metals": do_hydro,
# stars
"SFR": do_stars,
- "projection_mass_stars": do_stars,
# feedback
"abundances": do_feedback,
- "projection_metals": do_feedback,
- "iron_distribution": do_feedback,
- "mass_distribution": do_feedback
}
# Generate the figures
@@ -55,9 +35,6 @@ figsize = (100, 100)
figures = {
"projection_density": plt.figure(figsize=figsize),
"projection_temperature": plt.figure(figsize=figsize),
- "projection_metals": plt.figure(figsize=figsize),
- "projection_mass": plt.figure(figsize=figsize),
- "projection_mass_stars": plt.figure(figsize=figsize),
"phase_2d": plt.figure(figsize=figsize)
}
@@ -77,37 +54,17 @@ axes = {
figures["projection_temperature"], fig_size, subgrid,
add_all=True, share_all=True, cbar_mode="single",
cbar_size="2%", cbar_pad=0.02),
- "projection_mass": AxesGrid(
- figures["projection_mass"], fig_size, subgrid,
- add_all=True, share_all=True, cbar_mode="single",
- cbar_size="2%", cbar_pad=0.02),
- "projection_mass_stars": AxesGrid(
- figures["projection_mass_stars"], fig_size, subgrid,
- add_all=True, share_all=True, cbar_mode="single",
- cbar_size="2%", cbar_pad=0.02),
"phase_2d": AxesGrid(
figures["phase_2d"], fig_size, subgrid, axes_pad=0.05,
add_all=True, share_all=True, cbar_mode="single",
cbar_size="2%", cbar_pad=0.05, aspect=False),
- "projection_metals": AxesGrid(
- figures["projection_metals"], fig_size, subgrid,
- add_all=True, share_all=True, cbar_mode="single",
- cbar_size="2%", cbar_pad=0.02),
}
# Data
data = {
- "phase_1d_density": ([], []),
- "phase_1d_temperature": ([], []),
- "halo_distribution": ([], []),
- "velocity": ([], []),
- "metals": ([], []),
"SFR": ([], []),
- "profile": ([], []),
"abundances": ([], []),
- "iron_distribution": ([], []),
- "mass_distribution": ([], []),
}
names = []
@@ -121,62 +78,17 @@ def savePlot():
projection_plot.savePlot(figures["projection_temperature"],
"temperature")
- if do_plot["projection_metals"]:
- projection_plot.savePlot(figures["projection_metals"],
- "metals")
-
- if do_plot["projection_mass"]:
- projection_plot.savePlot(figures["projection_mass"],
- "mass")
-
- if do_plot["projection_mass_stars"]:
- projection_plot.savePlot(figures["projection_mass_stars"],
- "mass_stars")
-
- if do_plot["phase_1d_density"]:
- data["phase_1d_density"][1].extend(names)
- phase_plot.save1DPlotDensity(data["phase_1d_density"])
-
- if do_plot["phase_1d_temperature"]:
- data["phase_1d_temperature"][1].extend(names)
- phase_plot.save1DPlotTemperature(data["phase_1d_temperature"])
-
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"])
-
- if do_plot["metals"]:
- data["metals"][1].extend(names)
- metal_plot.save1DPlot(data["metals"])
-
if do_plot["SFR"]:
data["SFR"][1].extend(names)
star_plot.saveSFRPlot(data["SFR"])
- if do_plot["profile"]:
- data["profile"][1].extend(names)
- profile_plot.savePlot(data["profile"])
-
if do_plot["abundances"]:
data["abundances"][1].extend(names)
star_plot.saveAbundancesPlot(data["abundances"])
- if do_plot["iron_distribution"]:
- data["iron_distribution"][1].extend(names)
- star_plot.saveIronDistributionPlot(data["iron_distribution"])
-
- if do_plot["mass_distribution"]:
- data["mass_distribution"][1].extend(names)
- star_plot.saveMassDistributionPlot(data["mass_distribution"])
-
def doPlot(filename, i, name, center):
names.append(name)
@@ -215,72 +127,19 @@ def doPlot(filename, i, name, center):
f, name, i, figures["projection_temperature"],
axes["projection_temperature"])
- if do_plot["projection_metals"]:
- projection_plot.doMetalsPlot(
- f, name, i, figures["projection_metals"],
- axes["projection_metals"])
-
- # Do mass projection plot
- if do_plot["projection_mass"]:
- projection_plot.doMassPlot(
- f, name, i, figures["projection_mass"],
- axes["projection_mass"], "all")
-
- # Do stellar mass projection plot
- if do_plot["projection_mass_stars"]:
- projection_plot.doMassPlot(
- f, name, i, figures["projection_mass_stars"],
- axes["projection_mass_stars"], "stars")
-
- # 1D Phase plot
- if do_plot["phase_1d_density"]:
- p = phase_plot.do1DPlotDensity(f, name, i)
- data["phase_1d_density"][0].append(p)
-
- if do_plot["phase_1d_temperature"]:
- p = phase_plot.do1DPlotTemperature(f, name, i)
- data["phase_1d_temperature"][0].append(p)
-
# 2D Phase plot
if do_plot["phase_2d"]:
phase_plot.do2DPlot(f, name, i, figures["phase_2d"],
axes["phase_2d"])
- # halo distribution
- if do_plot["halo_distribution"]:
- 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)
-
- # metal plot
- if do_plot["metals"]:
- p = metal_plot.do1DPlot(f, name, i)
- data["metals"][0].append(p)
-
if do_plot["SFR"]:
p = star_plot.doSFRPlot(f, name, i)
data["SFR"][0].append(p)
- if do_plot["profile"]:
- p = profile_plot.doPlot(f, name, i)
- data["profile"][0].append(p)
-
if do_plot["abundances"]:
p = star_plot.doAbundancesPlot(f, name, i)
data["abundances"][0].append(p)
- if do_plot["iron_distribution"]:
- p = star_plot.doIronDistributionPlot(f, name, i)
- data["iron_distribution"][0].append(p)
-
- if do_plot["mass_distribution"]:
- p = star_plot.doMassDistributionPlot(f, name, i)
- data["mass_distribution"][0].append(p)
-
return center
diff --git a/examples/GEAR/ZoomIn/metal_plot.py b/examples/GEAR/ZoomIn/metal_plot.py
deleted file mode 100644
index 467bef8af8254f658ebb92c6e4188b840034bd9e..0000000000000000000000000000000000000000
--- a/examples/GEAR/ZoomIn/metal_plot.py
+++ /dev/null
@@ -1,38 +0,0 @@
-#!/usr/bin/env python3
-
-import yt
-import matplotlib.pyplot as plt
-
-
-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]):
- z = p.x.in_units("")
- m = p["Masses"].in_units("Msun").d
- plt.plot(z, m, linestyle="-", marker=markers[i],
- markeredgecolor='none', linewidth=1.2, alpha=0.8)
- plt.semilogx()
- plt.semilogy()
-
- plt.xlabel(r"Metal fraction", fontsize='large')
- plt.ylabel("$\mathrm{Mass\ (Msun)}$", 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("metals.png", bbox_inches='tight', pad_inches=0.03, dpi=300)
-
-
-def do1DPlot(f, name, i):
- sp = f.sphere(f.center, f.width)
- # Because ParticleProfilePlot doesn't exist, I will do the following trick.
- p = yt.create_profile(sp, ("PartType0", "Metallicity"),
- ("PartType0", "Masses"),
- weight_field=None, n_bins=50,
- accumulation=False)
-
- return p
diff --git a/examples/GEAR/ZoomIn/phase_plot.py b/examples/GEAR/ZoomIn/phase_plot.py
index c6c6e1d1a3765837e520de5c855e028bd47987e9..83fdc973669e87ec2f8534daff981cd9716a1b3d 100644
--- a/examples/GEAR/ZoomIn/phase_plot.py
+++ b/examples/GEAR/ZoomIn/phase_plot.py
@@ -17,56 +17,6 @@ def save2DPlot(fig):
pad_inches=0.03, dpi=300)
-def save1DPlotDensity(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]):
- rho = p.x.in_units("g/cm**3").d
- mass = p["Masses"].in_units("Msun").d
- mass[mass == 0] = np.nan
- plt.plot(rho, mass, linestyle="-", marker=markers[i],
- markeredgecolor='none', linewidth=1.2, alpha=0.8)
- plt.semilogx()
- plt.semilogy()
-
- plt.xlabel("$\mathrm{Density\ (g/cm^3)}$", fontsize='large')
- plt.ylabel(r"$\mathrm{Mass}\/\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("phase_1d_density.png", bbox_inches='tight',
- pad_inches=0.03, dpi=300)
-
-
-def save1DPlotTemperature(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]):
- rho = p.x.in_units("K").d
- mass = p["Masses"].in_units("Msun").d
- mass[mass == 0] = np.nan
- plt.plot(rho, mass, linestyle="-", marker=markers[i],
- markeredgecolor='none', linewidth=1.2, alpha=0.8)
- plt.semilogx()
- plt.semilogy()
-
- plt.xlabel("$\mathrm{Temperature\ K}$", fontsize='large')
- plt.ylabel(r"$\mathrm{Mass}\/\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("phase_1d_temperature.png", bbox_inches='tight',
- pad_inches=0.03, dpi=300)
-
-
def do2DPlot(f, name, i, fig, axes):
p = yt.PhasePlot(
f, ("PartType0", "density"), ("PartType0", "Temperature"),
@@ -118,22 +68,3 @@ def do2DPlot(f, name, i, fig, axes):
line = ln.Line2D(x_ticks, [y]*N, linestyle=":", linewidth=0.6,
color="k", alpha=0.7)
plot.axes.add_line(line)
-
-
-def do1DPlotDensity(f, name, i):
- sp = f.sphere(f.center, f.width)
- # Because ParticleProfilePlot doesn't exist, I will do the following trick.
- p = yt.create_profile(sp, ("PartType0", "density"),
- ("PartType0", "Masses"), weight_field=None,
- n_bins=40, accumulation=False)
- return p
-
-
-def do1DPlotTemperature(f, name, i):
- sp = f.sphere(f.center, f.width)
- # Because ParticleProfilePlot doesn't exist, I will do the following trick.
- p = yt.create_profile(sp, ("PartType0", "Temperature"),
- ("PartType0", "Masses"), weight_field=None,
- n_bins=40, accumulation=False)
-
- return p
diff --git a/examples/GEAR/ZoomIn/profile_plot.py b/examples/GEAR/ZoomIn/profile_plot.py
deleted file mode 100644
index dcabf608d1e60fc9a59bae494b4db14a26e3956d..0000000000000000000000000000000000000000
--- a/examples/GEAR/ZoomIn/profile_plot.py
+++ /dev/null
@@ -1,42 +0,0 @@
-#!/usr/bin/env python3
-
-import yt
-from yt.units import Msun, kpc
-import matplotlib.pyplot as plt
-
-limits_mass = (1e3 * Msun, 1e7 * Msun)
-
-
-def savePlot(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]):
- r = p.x.in_units("pc").d
- mass = p["Masses"].in_units("Msun").d
- plt.plot(r, mass, linestyle="-", marker=markers[i],
- markeredgecolor='none', linewidth=1.2, alpha=0.8)
- plt.semilogx()
- plt.semilogy()
-
- plt.xlabel(r"$\mathrm{Radius}\/\mathrm{(pc)}$", fontsize='large')
- plt.ylabel("$\mathrm{Integrated}\ \mathrm{mass}\ (Msun)}$",
- 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("integrated_mass.png", bbox_inches='tight',
- pad_inches=0.03, dpi=300)
-
-
-def doPlot(f, name, i):
- sp = f.sphere(f.center, f.width)
- # Because ParticleProfilePlot doesn't exist, I will do the following trick.
- p = yt.create_profile(sp, "particle_radius", "Masses",
- weight_field=None, n_bins=50,
- accumulation=True)
-
- return p
diff --git a/examples/GEAR/ZoomIn/projection_plot.py b/examples/GEAR/ZoomIn/projection_plot.py
index f77da7057b65b0ec7fecf2f0066a6eac77c3f720..d50e31b1421927145c13edaf77800117cfcc69d1 100644
--- a/examples/GEAR/ZoomIn/projection_plot.py
+++ b/examples/GEAR/ZoomIn/projection_plot.py
@@ -176,162 +176,3 @@ def doTemperaturePlot(f, name, i, fig, axes):
# Add code name
at = AnchoredText(name, loc=2, prop=dict(size=6), frameon=True)
plot.axes.add_artist(at)
-
-
-def doMassPlot(f, name, i, fig, axes, parttype):
- """
- Generate a mass projection (including dark matter) plot.
-
- Parameters
- ----------
-
- f: yt dataset
- The data set to use in the projection
-
- name: str
- The name to print on the plot
-
- i: int
- The index of the plot
-
- fig: matplotlib figure
- The figure to use
-
- axes: list
- The list of axes to use for the plot
-
- parttype: str
- The name of the particle type to use
- """
- if parttype == "stars":
- if name == "GEAR":
- parttype = "PartType1"
- else:
- parttype = "PartType4"
-
- direction = "x"
- field = (parttype, "Masses")
-
- # compute the projection
- p = yt.ParticleProjectionPlot(f, direction, field, center=f.center,
- width=f.width)
-
- # # Compute the limits
- p.set_unit("Masses", "Msun")
- data = p.to_fits_data()["Masses"].data
- global limits_mass
- if limits_mass is None:
- dmin = np.nanmin(data)
- if np.nanmin(data) == 0:
- dmin = 1e-6 * np.nanmax(data)
- else:
- dmin *= 0.5
- dmax = np.nanmax(data) * 2
- limits_mass = (dmin, dmax)
- else:
- if limits_mass[0] > np.nanmin(data):
- print("WARNING: data below min", np.nanmin(data))
- if limits_mass[1] < np.nanmax(data):
- print("WARNING: data above max", np.nanmax(data))
-
- # Adjust the plot
- p.set_cmap(field=field, cmap=cmap_mass)
- # p.set_log(field, True)
- p.set_zlim(field, limits_mass[0], limits_mass[1])
-
- # Draw it into the correct figure
- plot = p.plots[field]
-
- plot.figure = fig
- plot.axes = axes[i].axes
-
- # Add code name
- at = AnchoredText(name, loc=2, prop=dict(size=6), frameon=True)
- plot.axes.add_artist(at)
-
- # plot color bar
- if i != 0:
- plot.cax = axes.cbar_axes[0]
- p.hide_axes()
- p._setup_plots()
-
- # Add code name
- at = AnchoredText(name, loc=2, prop=dict(size=6), frameon=True)
- plot.axes.add_artist(at)
-
-
-def doMetalsPlot(f, name, i, fig, axes):
- """
- Generate a metal projection plot.
-
- Parameters
- ----------
-
- f: yt dataset
- The data set to use in the projection
-
- name: str
- The name to print on the plot
-
- i: int
- The index of the plot
-
- fig: matplotlib figure
- The figure to use
-
- axes: list
- The list of axes to use for the plot
- """
- direction = "x"
- field = ("PartType0", "Metallicity")
-
- # compute the projection
- p = yt.ProjectionPlot(f, direction, field, center=f.center,
- width=f.width, buff_size=(800, 800),
- weight_field=("PartType0", "Density"))
-
- # Compute the limits
- # p.set_unit("metallicity", "g / cm**2")
- data = p.to_fits_data()["Metallicity"].data
- global limits_metals
- if limits_metals is None:
- limits_metals = (data.min(), data.max())
- if limits_metals[0] > data.min():
- print("WARNING: data below min", data.min())
- if limits_metals[1] < data.max():
- print("WARNING: data above max", data.max())
-
- # Adjust the plot
- p.set_cmap(field=field, cmap=cmap_metals)
- p.set_log(field, True)
- if limits_metals[0] != limits_metals[1]:
- p.set_zlim(field, limits_metals[0], limits_metals[1])
-
- # Draw it into the correct figure
- plot = p.plots[field]
-
- plot.figure = fig
- plot.axes = axes[i].axes
-
- # Add code name
- at = AnchoredText(name, loc=2, prop=dict(size=6), frameon=True)
- plot.axes.add_artist(at)
-
- # plot color bar
- if i != 0:
- plot.cax = axes.cbar_axes[0]
- p.hide_axes()
- p._setup_plots()
-
- if i == 0:
- z = 1. / f.scale_factor - 1.
- text = "Redshift = %.2g" % z
- prop = {
- "color": "w"
- }
- at = AnchoredText(text, loc=3, prop=prop, frameon=False)
- plot.axes.add_artist(at)
-
- # Add code name
- at = AnchoredText(name, loc=2, prop=dict(size=6), frameon=True)
- plot.axes.add_artist(at)
diff --git a/examples/GEAR/ZoomIn/star_plot.py b/examples/GEAR/ZoomIn/star_plot.py
index a171642f759f4a7fa8634d651e5bd9b79f1ac5a3..855784bfa929126a368b23663e3d305c04001c4a 100644
--- a/examples/GEAR/ZoomIn/star_plot.py
+++ b/examples/GEAR/ZoomIn/star_plot.py
@@ -144,134 +144,3 @@ def doAbundancesPlot(f, name, i):
fe = getElement(f, sp, name, "Fe")
return mg, fe
-
-
-def saveIronDistributionPlot(profiles):
- plt.figure(figsize=(8, 8))
- plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
-
- fe_max = max(p.max() for p in profiles[0])
- fe_min = min(p.min() for p in profiles[0])
- fe_min = np.log10((fe_min + 1e-10) / fe_sol_abund)
- fe_max = np.log10(fe_max / fe_sol_abund)
-
- for i, p in enumerate(profiles[0]):
- fe = np.log10(p / fe_sol_abund)
-
- plt.hist(fe, bins=20, range=[fe_min, fe_max], histtype="step",
- density=True, alpha=0.8)
-
- plt.xlabel('[Fe / H]')
- plt.ylabel("Fraction of stars")
- plt.legend(profiles[1])
- plt.savefig("iron_distribution.png")
-
-
-def doIronDistributionPlot(f, name, i):
- sp = f.sphere(f.center, f.width)
- # Because ParticleProfilePlot doesn't exist, I will do the following trick.
- fe = getElement(f, sp, name, "Fe")
- return fe
-
-
-def saveMassDistributionPlot(profiles):
- # Do the time - stars
- plt.figure(figsize=(8, 8))
- plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
- markers = ["s", "o"]
- zorder = [2.5, 2.]
-
- for i, p in enumerate(profiles[0]):
- plt.plot(p["birth_time"], p["stars"], marker=markers[i],
- markeredgecolor="none", linestyle="", alpha=0.8,
- zorder=zorder[i])
- plt.semilogx()
- plt.semilogy()
-
- plt.xlabel("Star formation time (Gyr)", fontsize='large')
- plt.ylabel("Mass of the star $(M_{\odot})$", fontsize='large')
- plt.legend(profiles[1], 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("mass_lost_stars.png", bbox_inches='tight',
- pad_inches=0.03, dpi=300)
-
- # Do the stars
- plt.figure(figsize=(8, 8))
- plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
-
- mass_max = max(p["stars"].max() for p in profiles[0])
- mass_min = min(p["stars"].min() for p in profiles[0])
-
- for i, mass in enumerate(profiles[0]):
- plt.hist(mass["stars"], bins=100, range=[mass_min, mass_max],
- histtype="step", density=True, alpha=0.8)
- plt.semilogx()
- plt.semilogy()
-
- plt.xlabel("Mass of the stars $(M_{\odot})$", fontsize='large')
- plt.ylabel("Fraction of stars", fontsize='large')
- plt.legend(profiles[1], 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("mass_distribution_stars.png", bbox_inches='tight',
- pad_inches=0.03, dpi=300)
-
- # Do the gas
- plt.figure(figsize=(8, 8))
- plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=0.5)
-
- mass_max = max(p["gas"].max() for p in profiles[0])
- mass_min = min(p["gas"].min() for p in profiles[0])
-
- for i, mass in enumerate(profiles[0]):
- plt.hist(mass["gas"], bins=100, range=[mass_min, mass_max],
- histtype="step", density=True, alpha=0.8)
- plt.semilogx()
- plt.semilogy()
-
- plt.xlabel("Mass of the gas $(M_{\odot})$", fontsize='large')
- plt.ylabel("Fraction of gas", fontsize='large')
- plt.legend(profiles[1], 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("mass_distribution_gas.png", bbox_inches='tight',
- pad_inches=0.03, dpi=300)
-
-
-def doMassDistributionPlot(f, name, i):
- sp = f.sphere(f.center, f.width)
-
- gas = sp[("PartType0", "Masses")].in_units("Msun").d
-
- if name == "GEAR":
- stars = sp[("PartType1", "Masses")].in_units("Msun").d
- a = sp["PartType1", "StarFormationTime"]
- else:
- stars = sp[("PartType4", "Masses")].in_units("Msun").d
- a = sp["PartType4", "BirthScaleFactors"]
-
- cosmo = Cosmology(hubble_constant=f.hubble_constant,
- omega_matter=f.omega_matter,
- omega_lambda=f.omega_lambda)
- z = np.array(1. / a - 1.)
- with Pool() as p:
- formation_time = p.map(
- partial(cosmo.lookback_time, z_f=1e6), z)
- formation_time = f.arr(list(formation_time), "s").in_units("Gyr")
-
- d = {
- "gas": gas,
- "stars": stars,
- "birth_time": formation_time
- }
- return d
diff --git a/examples/GEAR/ZoomIn/velocity_plot.py b/examples/GEAR/ZoomIn/velocity_plot.py
deleted file mode 100644
index 5e1f5cf0fd8e3ccb22cd9e5999a83c75c5f878f6..0000000000000000000000000000000000000000
--- a/examples/GEAR/ZoomIn/velocity_plot.py
+++ /dev/null
@@ -1,52 +0,0 @@
-#!/usr/bin/env python3
-
-import yt
-import matplotlib.pyplot as plt
-
-DMO = True
-
-
-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')
- if DMO:
- plt.ylabel(r"$\mathrm{Mass}\/\mathrm{all}\/\mathrm{(M_{\odot})}$",
- fontsize='large')
- else:
- plt.ylabel(r"$\mathrm{Mass}\/\mathrm{Gas}\/\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):
- global DMO
- part_type = "all"
- if f.particle_type_counts["PartType0"] != 0:
- DMO = False
- part_type = "PartType0"
-
- sp = f.sphere("c", f.width)
-
- # Because ParticleProfilePlot doesn't exist, I will do the following trick.
- p = yt.create_profile(sp, (part_type, "particle_velocity_magnitude"),
- (part_type, "Masses"),
- weight_field=None, n_bins=50,
- accumulation=False)
-
- return p