From 9b3c6ca251379b1a4c933ee163844f68435de403 Mon Sep 17 00:00:00 2001 From: loikki Date: Fri, 31 Jan 2020 09:50:51 +0100 Subject: [PATCH] GEAR: 1D phase plot done --- examples/GEAR/ZoomIn/cleanupSwift.py | 2 +- examples/GEAR/ZoomIn/halo_distribution.py | 13 +++++++ examples/GEAR/ZoomIn/make_image.py | 46 ++++++++++++++++------ examples/GEAR/ZoomIn/phase_plot.py | 47 ++++++++++++++++++----- examples/GEAR/ZoomIn/projection_plot.py | 4 +- 5 files changed, 88 insertions(+), 24 deletions(-) create mode 100644 examples/GEAR/ZoomIn/halo_distribution.py diff --git a/examples/GEAR/ZoomIn/cleanupSwift.py b/examples/GEAR/ZoomIn/cleanupSwift.py index 5889bca9e..75f7174c4 100644 --- a/examples/GEAR/ZoomIn/cleanupSwift.py +++ b/examples/GEAR/ZoomIn/cleanupSwift.py @@ -46,7 +46,7 @@ for i in range(NPartType): ("Velocities", "Velocities", a**0.5), ("Density", "Densities", 1. / h**2), ("Entropies", "Entropies", 1.), - ("InternalEnergy", "InternalEnergies", 1.), + ("InternalEnergy", "InternalEnergies", 1. / a**2), ("SmoothingLength", "SmoothingLengths", h) ] diff --git a/examples/GEAR/ZoomIn/halo_distribution.py b/examples/GEAR/ZoomIn/halo_distribution.py new file mode 100644 index 000000000..5a2ec7cce --- /dev/null +++ b/examples/GEAR/ZoomIn/halo_distribution.py @@ -0,0 +1,13 @@ +#!/usr/bin/env + +import yt +from yt.extensions.astro_analysis.halo_analysis.api import HaloCatalog + + +def doPlot(f, name, i): + hc = HaloCatalog(data_ds=f, finder_method="fof", + finder_kwargs={"padding": 0.02}) + hc.create() + print(hc) + print(dir(hc)) + exit() diff --git a/examples/GEAR/ZoomIn/make_image.py b/examples/GEAR/ZoomIn/make_image.py index e9fbc68a5..d1900eb18 100644 --- a/examples/GEAR/ZoomIn/make_image.py +++ b/examples/GEAR/ZoomIn/make_image.py @@ -7,6 +7,7 @@ import sys import matplotlib import projection_plot import phase_plot +import halo_distribution import add_fields import matplotlib.pyplot as plt from mpl_toolkits.axes_grid1 import AxesGrid @@ -24,8 +25,10 @@ gear = "gear/snapshot_%04i.hdf5" % snap do_plot = { "projection_density": False, "projection_temperature": False, - "projection_mass": True, - "phase": False + "projection_mass": False, + "halo_distribution": True, + "phase_1d": False, + "phase_2d": False } # Generate the figures @@ -35,7 +38,7 @@ figures = { "projection_density": plt.figure(figsize=figsize), "projection_temperature": plt.figure(figsize=figsize), "projection_mass": plt.figure(figsize=figsize), - "phase": plt.figure(figsize=figsize) + "phase_2d": plt.figure(figsize=figsize) } # define global variables @@ -58,13 +61,19 @@ axes = { figures["projection_mass"], fig_size, subgrid, add_all=True, share_all=True, cbar_mode="single", cbar_size="2%", cbar_pad=0.02), - "phase": AxesGrid( - figures["phase"], fig_size, subgrid, axes_pad=0.05, + "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) } +# Data +data = { + "phase_1d": ([], []) +} + + def savePlot(): if do_plot["projection_density"]: projection_plot.savePlot(figures["projection_density"], @@ -77,13 +86,16 @@ def savePlot(): projection_plot.savePlot(figures["projection_mass"], "mass") - if do_plot["phase"]: - phase_plot.savePlot(figures["phase"]) + if do_plot["phase_1d"]: + phase_plot.save1DPlot(data["phase_1d"]) + + if do_plot["phase_2d"]: + phase_plot.save2DPlot(figures["phase_2d"]) def doPlot(filename, i, name): f = yt.load(filename) - if (do_plot["projection_temperature"] or do_plot["phase"]): + if (do_plot["projection_temperature"] or do_plot["phase_2d"]): add_fields.addTemperature(f) global scale_factor @@ -112,10 +124,20 @@ def doPlot(filename, i, name): f, name, i, figures["projection_mass"], axes["projection_mass"]) - # Phase plot - if do_plot["phase"]: - phase_plot.doPlot(f, name, i, figures["phase"], - axes["phase"]) + # 1D Phase plot + 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"]: + phase_plot.do2DPlot(f, name, i, figures["phase_2d"], + axes["phase_2d"]) + + # halo distribution + if do_plot["halo_distribution"]: + halo_distribution.doPlot(f, name, i) doPlot(swift, 0, "SWIFT") diff --git a/examples/GEAR/ZoomIn/phase_plot.py b/examples/GEAR/ZoomIn/phase_plot.py index 7a6917d31..63f78e11e 100644 --- a/examples/GEAR/ZoomIn/phase_plot.py +++ b/examples/GEAR/ZoomIn/phase_plot.py @@ -1,25 +1,45 @@ #!/usr/bin/env python3 import yt -from yt.units import Msun, amu, cm, K +from yt.units import Msun, amu, cm, K, kpc from matplotlib.offsetbox import AnchoredText 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_mass = (1e3 * Msun, 1e7 * Msun) -def savePlot(fig): - fig.savefig("phase.png", bbox_inches="tight", + +def save2DPlot(fig): + fig.savefig("phase_2d.png", bbox_inches="tight", pad_inches=0.03, dpi=300) -def doPlot(f, name, i, fig, axes): - # limits in g / cm3 - # u_rho = g / cm**3 - # limits = [1e-33 * u_rho, 1e-24 * u_rho] +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]): + rho = p.x.in_units("g/cm**3").d + mass = p["Masses"].in_units("Msun").d + 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{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("phase_1d.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"), ("PartType0", "mass"), @@ -54,6 +74,15 @@ def doPlot(f, name, i, fig, axes): # p.hide_axes() p._setup_plots() - at = AnchoredText(name, loc=2, prop=dict(size=6), frameon=True) plot.axes.add_artist(at) + + +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", "density"), ("PartType0", "Masses"), + weight_field=None, n_bins=50, + accumulation=False) + + return p diff --git a/examples/GEAR/ZoomIn/projection_plot.py b/examples/GEAR/ZoomIn/projection_plot.py index e4b35f3ff..cbc4faedf 100644 --- a/examples/GEAR/ZoomIn/projection_plot.py +++ b/examples/GEAR/ZoomIn/projection_plot.py @@ -11,7 +11,7 @@ center = (1441.9954833984375000, 1666.1169433593750000, 1891.6248779296875000) small_width = 200 * kpc -large_width = 2500 * kpc +large_width = 1000 * kpc limits_density = (1e-10 * g / cm**2, 1e-2 * g / cm**2) limits_temperature = (10 * K, 1e5 * K) limits_mass = None @@ -205,7 +205,7 @@ def doMassPlot(f, name, i, fig, axes): field = ("all", "Masses") # compute the projection - p = yt.ParticleProjectionPlot(f, direction, field, center="center", + p = yt.ParticleProjectionPlot(f, direction, field, center=center, width=large_width) # # Compute the limits -- 2.26.2