Commit 9b3c6ca2 authored by Loic Hausammann's avatar Loic Hausammann

GEAR: 1D phase plot done

parent 57fc59ac
......@@ -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)
]
......
#!/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()
......@@ -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")
......
#!/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
......@@ -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
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment