Commit 82d556e8 authored by Folkert Nobels's avatar Folkert Nobels
Browse files

Rename python file and format with black

parent e2342e4d
import matplotlib
matplotlib.use("Agg")
from pylab import *
from scipy import stats
......@@ -36,12 +37,12 @@ f = h5.File(filename, "r")
# Physical constants
k_in_cgs = 1.38064852e-16
mH_in_cgs = 1.6737236e-24
year_in_cgs = 3600. * 24 * 365.
year_in_cgs = 3600.0 * 24 * 365.0
Msun_in_cgs = 1.98848e33
# Gemoetry info
boxsize = f["/Header"].attrs["BoxSize"]
centre = boxsize / 2.
centre = boxsize / 2.0
# Read units
unit_length_in_cgs = f["/Units"].attrs["Unit length in cgs (U_L)"]
......@@ -56,31 +57,33 @@ KS_thresh_slope = float(f["/Parameters"].attrs["SchayeSF:MetDep_SFthresh_Slope"]
KS_thresh_norm = float(f["/Parameters"].attrs["SchayeSF:thresh_norm_HpCM3"])
KS_gas_fraction = float(f["/Parameters"].attrs["SchayeSF:fg"])
KS_thresh_max_norm = float(f["/Parameters"].attrs["SchayeSF:thresh_max_norm_HpCM3"])
KS_gamma_effective = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_gamma_effective"])
KS_gamma_effective = float(
f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_gamma_effective"]
)
EAGLE_Z = float(f["/Parameters"].attrs["EAGLEChemistry:init_abundance_metal"])
# Read gas properties
gas_pos = f["/PartType0/Coordinates"][:,:]
gas_pos = f["/PartType0/Coordinates"][:, :]
gas_mass = f["/PartType0/Masses"][:]
gas_rho = f["/PartType0/Density"][:]
gas_T = f["/PartType0/Temperature"][:]
gas_SFR = f["/PartType0/SFR"][:]
gas_XH = f["/PartType0/ElementAbundance"][:,0]
gas_XH = f["/PartType0/ElementAbundance"][:, 0]
gas_Z = f["/PartType0/Metallicity"][:]
gas_hsml = f["/PartType0/SmoothingLength"][:]
# Centre the box
gas_pos[:,0] -= centre[0]
gas_pos[:,1] -= centre[1]
gas_pos[:,2] -= centre[2]
gas_pos[:, 0] -= centre[0]
gas_pos[:, 1] -= centre[1]
gas_pos[:, 2] -= centre[2]
# Turn the mass into better units
gas_mass *= (unit_mass_in_cgs / Msun_in_cgs)
gas_mass *= unit_mass_in_cgs / Msun_in_cgs
# Turn the SFR into better units
gas_SFR = np.maximum(gas_SFR, np.zeros(np.size(gas_SFR)))
gas_SFR /= (unit_time_in_cgs / year_in_cgs)
gas_SFR *= (unit_mass_in_cgs / Msun_in_cgs)
gas_SFR /= unit_time_in_cgs / year_in_cgs
gas_SFR *= unit_mass_in_cgs / Msun_in_cgs
# Make it a Hydrogen number density
gas_nH = gas_rho * unit_mass_in_cgs / unit_length_in_cgs ** 3
......@@ -89,32 +92,32 @@ gas_nH *= gas_XH
# Equations of state
eos_cool_rho = np.logspace(-5, 5, 1000)
eos_cool_T = eos_cool_rho**0. * 8000.
eos_cool_T = eos_cool_rho ** 0.0 * 8000.0
eos_Jeans_rho = np.logspace(-1, 5, 1000)
eos_Jeans_T = (eos_Jeans_rho/ 10**(-1))**(1./3.) * 8000.
eos_Jeans_T = (eos_Jeans_rho / 10 ** (-1)) ** (1.0 / 3.0) * 8000.0
# Plot the phase space diagram
figure()
subplot(111, xscale="log", yscale="log")
plot(eos_cool_rho, eos_cool_T, 'k--', lw=0.6)
plot(eos_Jeans_rho, eos_Jeans_T, 'k--', lw=0.6)
plot(eos_cool_rho, eos_cool_T, "k--", lw=0.6)
plot(eos_Jeans_rho, eos_Jeans_T, "k--", lw=0.6)
scatter(gas_nH, gas_T, s=0.2)
xlabel("${\\rm Density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0)
ylabel("${\\rm Temperature}~T~[{\\rm K}]$", labelpad=2)
xlim(1e-4, 3e3)
ylim(500., 2e5)
ylim(500.0, 2e5)
savefig("rhoT.png", dpi=200)
# Plot the phase space diagram for SF gas
figure()
subplot(111, xscale="log", yscale="log")
plot(eos_cool_rho, eos_cool_T, 'k--', lw=1)
plot(eos_Jeans_rho, eos_Jeans_T, 'k--', lw=1)
scatter(gas_nH[gas_SFR > 0.], gas_T[gas_SFR > 0.], s=0.2)
plot(eos_cool_rho, eos_cool_T, "k--", lw=1)
plot(eos_Jeans_rho, eos_Jeans_T, "k--", lw=1)
scatter(gas_nH[gas_SFR > 0.0], gas_T[gas_SFR > 0.0], s=0.2)
xlabel("${\\rm Density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0)
ylabel("${\\rm Temperature}~T~[{\\rm K}]$", labelpad=2)
xlim(1e-4, 3e3)
ylim(500., 2e5)
ylim(500.0, 2e5)
savefig("rhoT_SF.png", dpi=200)
########################################################################3
......@@ -123,7 +126,7 @@ savefig("rhoT_SF.png", dpi=200)
figure()
subplot(111, xscale="log", yscale="log")
scatter(gas_nH, gas_SFR, s=0.2)
plot([1, 100], 2e-5 * np.array([1, 100])**0.266667, 'k--', lw=1)
plot([1, 100], 2e-5 * np.array([1, 100]) ** 0.266667, "k--", lw=1)
xlabel("${\\rm Density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0)
ylabel("${\\rm SFR}~[{\\rm M_\\odot~\\cdot~yr^{-1}}]$", labelpad=-7)
xlim(1e-4, 3e3)
......@@ -133,7 +136,14 @@ savefig("rho_SFR.png", dpi=200)
########################################################################3
# Select gas in a pillow box around the galaxy
mask = (gas_pos[:, 0] > -15) & (gas_pos[:, 0] < 15) & (gas_pos[:, 1] > -15) & (gas_pos[:, 1] < 15) & (gas_pos[:, 2] < 1.) & (gas_pos[:, 2] > -1.)
mask = (
(gas_pos[:, 0] > -15)
& (gas_pos[:, 0] < 15)
& (gas_pos[:, 1] > -15)
& (gas_pos[:, 1] < 15)
& (gas_pos[:, 2] < 1.0)
& (gas_pos[:, 2] > -1.0)
)
gas_pos = gas_pos[mask, :]
gas_SFR = gas_SFR[mask]
gas_nH = gas_nH[mask]
......@@ -163,7 +173,7 @@ ylim(-12, 12)
savefig("edge_on.png", dpi=200)
# Now a SF map
rcParams.update({ "figure.figsize": (4.15, 3.15)})
rcParams.update({"figure.figsize": (4.15, 3.15)})
figure()
subplot(111)
scatter(gas_pos[:, 0], gas_pos[:, 1], s=0.1, c=gas_SFR)
......@@ -182,38 +192,42 @@ savefig("SF_face_on.png", dpi=200)
x_edges = np.linspace(-15, 15, 31)
y_edges = np.linspace(-15, 15, 31)
map_mass,_,_,_ = stats.binned_statistic_2d(gas_pos[:,0], gas_pos[:,1], gas_mass, statistic='sum', bins=(x_edges, y_edges))
map_SFR,_,_,_ = stats.binned_statistic_2d(gas_pos[:,0], gas_pos[:,1], gas_SFR, statistic='sum', bins=(x_edges, y_edges))
map_mass, _, _, _ = stats.binned_statistic_2d(
gas_pos[:, 0], gas_pos[:, 1], gas_mass, statistic="sum", bins=(x_edges, y_edges)
)
map_SFR, _, _, _ = stats.binned_statistic_2d(
gas_pos[:, 0], gas_pos[:, 1], gas_SFR, statistic="sum", bins=(x_edges, y_edges)
)
qv = QuickView(
gas_pos,
mass=gas_mass,
r="infinity",
xsize=len(x_edges)-1,
ysize=len(y_edges)-1,
xsize=len(x_edges) - 1,
ysize=len(y_edges) - 1,
p=0, # Viewing angle theta
roll=0, # Viewing angle phi
plot=False,
logscale=False,
hsml=gas_hsml
hsml=gas_hsml,
)
map_mass2 = qv.get_image()
extent_mass = qv.get_extent()
gas_SFR[gas_SFR<=0] = 1e-10
gas_SFR[gas_SFR <= 0] = 1e-10
qv = QuickView(
gas_pos,
mass=gas_SFR,
r="infinity",
xsize=len(x_edges)-1,
ysize=len(y_edges)-1,
xsize=len(x_edges) - 1,
ysize=len(y_edges) - 1,
p=0, # Viewing angle theta
roll=0, # Viewing angle phi
plot=False,
logscale=False,
hsml=gas_hsml
hsml=gas_hsml,
)
map_SFR2 = qv.get_image()
......@@ -234,7 +248,7 @@ savefig("Map_mass.png", dpi=200)
# Mass map 2
figure()
subplot(111)
imshow(np.log10(map_mass2),extent=extent_mass)
imshow(np.log10(map_mass2), extent=extent_mass)
colorbar()
xlim(-12, 12)
ylim(-12, 12)
......@@ -245,7 +259,7 @@ savefig("Map_mass_SPHVIEWER.png", dpi=200)
# SF map
figure()
subplot(111)
pcolormesh(x_edges, y_edges, np.log10(map_SFR), vmax = -.5, vmin=-4.5)
pcolormesh(x_edges, y_edges, np.log10(map_SFR), vmax=-0.5, vmin=-4.5)
colorbar()
xlim(-12, 12)
ylim(-12, 12)
......@@ -254,11 +268,11 @@ ylabel("${\\rm Pos}~y~[{\\rm kpc}]$", labelpad=-3)
savefig("Map_SFR.png", dpi=200)
plot_map_SFR2 = np.zeros(np.shape(map_SFR2))
plot_map_SFR2[map_SFR2>1e-6] = map_SFR2[map_SFR2>1e-6]
plot_map_SFR2[map_SFR2 > 1e-6] = map_SFR2[map_SFR2 > 1e-6]
# SF map 2
figure()
subplot(111)
imshow(np.log10(plot_map_SFR2),extent=extent_SFR, vmax = -.5, vmin=-4.5)
imshow(np.log10(plot_map_SFR2), extent=extent_SFR, vmax=-0.5, vmin=-4.5)
colorbar()
xlim(-12, 12)
ylim(-12, 12)
......@@ -272,45 +286,89 @@ savefig("Map_SFR_SPHVIEWER.png", dpi=200)
map_SFR[map_SFR < 1e-6] = 1e-6
# Theoretical threshold (assumes all gas has the same Z)
KS_n_thresh = KS_thresh_norm * (gas_Z[0] / KS_thresh_Z0)**KS_thresh_slope
if np.isfinite(KS_n_thresh)==False:
KS_n_thresh = KS_thresh_norm * (gas_Z[0] / KS_thresh_Z0) ** KS_thresh_slope
if np.isfinite(KS_n_thresh) == False:
KS_n_thresh = KS_thresh_max_norm
KS_sigma_thresh = 29. * np.sqrt(KS_gas_fraction) * np.sqrt(KS_n_thresh)
KS_sigma_thresh = 29.0 * np.sqrt(KS_gas_fraction) * np.sqrt(KS_n_thresh)
# Theoretical KS law
KS_sigma_mass = np.logspace(-1, 3, 100)
KS_sigma_SFR = KS_law_norm * KS_sigma_mass**KS_law_slope
KS_sigma_SFR = KS_law_norm * KS_sigma_mass ** KS_law_slope
# KS relation
rcParams.update({ "figure.figsize": (3.15, 3.15), "figure.subplot.left": 0.18,})
rcParams.update({"figure.figsize": (3.15, 3.15), "figure.subplot.left": 0.18})
figure()
subplot(111, xscale="log", yscale="log")
plot(KS_sigma_mass, KS_sigma_SFR, 'k--', lw=0.6)
plot([KS_sigma_thresh, KS_sigma_thresh], [1e-8, 1e8], 'k--', lw=0.6)
text(KS_sigma_thresh*0.95, 2.2, "$\\Sigma_{\\rm c} = %.2f~{\\rm M_\\odot\\cdot pc^{-2}}$"%KS_sigma_thresh, va="top", ha="right", rotation=90, fontsize=7)
text(16, 10**(-3.5), "$n_{\\rm H,c} = %.3f~{\\rm cm^{-3}}$"%KS_n_thresh, fontsize=7)
text(16, 2e-6, "${\\rm K\\textendash S~law}$:\n$\\Sigma_{\\rm SFR} = A \\times \\Sigma_g^n$\n$n=%.1f$\n$A=%.3f\\times10^{-4}~{\\rm M_\\odot / yr^{1} / kpc^{2}}$\n$f_{\\rm g} = %.1f$\n$\gamma_{\\rm eos} = %.3f$\n$Z=%1.4f$"%(KS_law_slope, KS_law_norm*10**4, KS_gas_fraction, KS_gamma_effective, EAGLE_Z), fontsize=7)
plot(KS_sigma_mass, KS_sigma_SFR, "k--", lw=0.6)
plot([KS_sigma_thresh, KS_sigma_thresh], [1e-8, 1e8], "k--", lw=0.6)
text(
KS_sigma_thresh * 0.95,
2.2,
"$\\Sigma_{\\rm c} = %.2f~{\\rm M_\\odot\\cdot pc^{-2}}$" % KS_sigma_thresh,
va="top",
ha="right",
rotation=90,
fontsize=7,
)
text(16, 10 ** (-3.5), "$n_{\\rm H,c} = %.3f~{\\rm cm^{-3}}$" % KS_n_thresh, fontsize=7)
text(
16,
2e-6,
"${\\rm K\\textendash S~law}$:\n$\\Sigma_{\\rm SFR} = A \\times \\Sigma_g^n$\n$n=%.1f$\n$A=%.3f\\times10^{-4}~{\\rm M_\\odot / yr^{1} / kpc^{2}}$\n$f_{\\rm g} = %.1f$\n$\gamma_{\\rm eos} = %.3f$\n$Z=%1.4f$"
% (
KS_law_slope,
KS_law_norm * 10 ** 4,
KS_gas_fraction,
KS_gamma_effective,
EAGLE_Z,
),
fontsize=7,
)
scatter(map_mass.flatten() / 1e6, map_SFR.flatten(), s=0.4)
xlim(0.3, 900)
ylim(3e-7, 3)
xlabel("$\\Sigma_g~[{\\rm M_\\odot\\cdot pc^{-2}}]$", labelpad=0)
ylabel("$\\Sigma_{\\rm SFR}~[{\\rm M_\\odot \\cdot yr^{-1} \\cdot kpc^{-2}}]$", labelpad=0)
ylabel(
"$\\Sigma_{\\rm SFR}~[{\\rm M_\\odot \\cdot yr^{-1} \\cdot kpc^{-2}}]$", labelpad=0
)
savefig("KS_law.png", dpi=200)
close()
plot_map_SFR2[plot_map_SFR2<=0] = 1e-6
plot_map_SFR2[plot_map_SFR2 <= 0] = 1e-6
rcParams.update({ "figure.figsize": (3.15, 3.15), "figure.subplot.left": 0.18,})
rcParams.update({"figure.figsize": (3.15, 3.15), "figure.subplot.left": 0.18})
figure()
subplot(111, xscale="log", yscale="log")
plot(KS_sigma_mass, KS_sigma_SFR, 'k--', lw=0.6)
plot([KS_sigma_thresh, KS_sigma_thresh], [1e-8, 1e8], 'k--', lw=0.6)
text(KS_sigma_thresh*0.95, 2.2, "$\\Sigma_{\\rm c} = %.2f~{\\rm M_\\odot\\cdot pc^{-2}}$"%KS_sigma_thresh, va="top", ha="right", rotation=90, fontsize=7)
text(16, 10**(-3.5), "$n_{\\rm H,c} = %.3f~{\\rm cm^{-3}}$"%KS_n_thresh, fontsize=7)
text(16, 2e-6, "${\\rm K\\textendash S~law}$:\n$\\Sigma_{\\rm SFR} = A \\times \\Sigma_g^n$\n$n=%.1f$\n$A=%.3f\\times10^{-4}~{\\rm M_\\odot / yr^{1} / kpc^{2}}$\n$f_{\\rm g} = %.1f$\n$\gamma_{\\rm eos} = %.3f$\n$Z=%1.4f$"%(KS_law_slope, KS_law_norm*10**4, KS_gas_fraction, KS_gamma_effective, EAGLE_Z), fontsize=7)
plot(KS_sigma_mass, KS_sigma_SFR, "k--", lw=0.6)
plot([KS_sigma_thresh, KS_sigma_thresh], [1e-8, 1e8], "k--", lw=0.6)
text(
KS_sigma_thresh * 0.95,
2.2,
"$\\Sigma_{\\rm c} = %.2f~{\\rm M_\\odot\\cdot pc^{-2}}$" % KS_sigma_thresh,
va="top",
ha="right",
rotation=90,
fontsize=7,
)
text(16, 10 ** (-3.5), "$n_{\\rm H,c} = %.3f~{\\rm cm^{-3}}$" % KS_n_thresh, fontsize=7)
text(
16,
2e-6,
"${\\rm K\\textendash S~law}$:\n$\\Sigma_{\\rm SFR} = A \\times \\Sigma_g^n$\n$n=%.1f$\n$A=%.3f\\times10^{-4}~{\\rm M_\\odot / yr^{1} / kpc^{2}}$\n$f_{\\rm g} = %.1f$\n$\gamma_{\\rm eos} = %.3f$\n$Z=%1.4f$"
% (
KS_law_slope,
KS_law_norm * 10 ** 4,
KS_gas_fraction,
KS_gamma_effective,
EAGLE_Z,
),
fontsize=7,
)
scatter(map_mass2.flatten() / 1e6, plot_map_SFR2.flatten(), s=0.4)
xlim(0.3, 900)
ylim(3e-7, 3)
xlabel("$\\Sigma_g~[{\\rm M_\\odot\\cdot pc^{-2}}]$", labelpad=0)
ylabel("$\\Sigma_{\\rm SFR}~[{\\rm M_\\odot \\cdot yr^{-1} \\cdot kpc^{-2}}]$", labelpad=0)
ylabel(
"$\\Sigma_{\\rm SFR}~[{\\rm M_\\odot \\cdot yr^{-1} \\cdot kpc^{-2}}]$", labelpad=0
)
savefig("KS_law_SPHVIEWER.png", dpi=200)
Supports Markdown
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