Skip to content
Snippets Groups Projects
Commit d1eb8dc3 authored by Folkert Nobels's avatar Folkert Nobels
Browse files

Update plotSolutions of star formation example and remove the two other examples

parent c9692302
No related branches found
No related tags found
1 merge request!705Star formation following Schaye08
Hide whitespace changes
Inline Side-by-side
Showing
with 58 additions and 1543 deletions
examples/IsolatedGalaxy_starformation/plotSolution.py
+ 58
− 3
View file @ d1eb8dc3
......@@ -39,6 +39,9 @@ k_in_cgs = 1.38064852e-16
mH_in_cgs = 1.6737236e-24
year_in_cgs = 3600.0 * 24 * 365.0
Msun_in_cgs = 1.98848e33
G_in_cgs = 6.67259e-8
pc_in_cgs = 3.08567758e18
Msun_p_pc2 = Msun_in_cgs / pc_in_cgs**2
# Gemoetry info
boxsize = f["/Header"].attrs["BoxSize"]
......@@ -49,6 +52,9 @@ unit_length_in_cgs = f["/Units"].attrs["Unit length in cgs (U_L)"]
unit_mass_in_cgs = f["/Units"].attrs["Unit mass in cgs (U_M)"]
unit_time_in_cgs = f["/Units"].attrs["Unit time in cgs (U_t)"]
# Calculate Gravitational constant in internal units
G = G_in_cgs * ( unit_length_in_cgs**3 / unit_mass_in_cgs / unit_time_in_cgs**2)**(-1)
# Read parameters of the model
KS_law_slope = float(f["/Parameters"].attrs["SchayeSF:SchmidtLawExponent"])
KS_law_norm = float(f["/Parameters"].attrs["SchayeSF:SchmidtLawCoeff_MSUNpYRpKPC2"])
......@@ -60,7 +66,18 @@ KS_thresh_max_norm = float(f["/Parameters"].attrs["SchayeSF:thresh_max_norm_HpCM
KS_gamma_effective = float(
f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_gamma_effective"]
)
KS_high_den_thresh = float(f["/Parameters"].attrs["SchayeSF:SchmidtLawHighDens_thresh_HpCM3"])
KS_law_slope_high_den = float(f["/Parameters"].attrs["SchayeSF:SchmidtLawHighDensExponent"])
EAGLE_Z = float(f["/Parameters"].attrs["EAGLEChemistry:init_abundance_metal"])
EAGLEfloor_rho_norm = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_density_threshold_H_p_cm3"])
EAGLEfloor_T = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Jeans_temperature_norm_K"])
EAGLEfloor_cool_rho = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_density_threshold_H_p_cm3"])
EAGLEfloor_cool_temperature_norm_K = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_temperature_norm_K"])
EAGLEfloor_cool_gamma_effective = float(f["/Parameters"].attrs["EAGLEEntropyFloor:Cool_gamma_effective"])
KS_law_norm_cgs = KS_law_norm * Msun_in_cgs / ( 1e6 * pc_in_cgs**2 * year_in_cgs )
gamma = 5./3.
KS_press_norm = k_in_cgs * EAGLEfloor_T * EAGLEfloor_rho_norm
# Read gas properties
gas_pos = f["/PartType0/Coordinates"][:, :]
......@@ -71,6 +88,11 @@ gas_SFR = f["/PartType0/SFR"][:]
gas_XH = f["/PartType0/ElementAbundance"][:, 0]
gas_Z = f["/PartType0/Metallicity"][:]
gas_hsml = f["/PartType0/SmoothingLength"][:]
gas_sSFR = f["/PartType0/sSFR"][:]
# Read the Star properties
stars_pos = f["/PartType4/Coordinates"][:, :]
stars_BirthDensity = f["/PartType4/BirthDensity"][:]
# Centre the box
gas_pos[:, 0] -= centre[0]
......@@ -90,12 +112,15 @@ gas_nH = gas_rho * unit_mass_in_cgs / unit_length_in_cgs ** 3
gas_nH /= mH_in_cgs
gas_nH *= gas_XH
stars_BirthDensity = gas_rho * unit_mass_in_cgs / unit_length_in_cgs ** 3
stars_BirthDensity /= mH_in_cgs
stars_BirthDensity *= gas_XH
# Equations of state
eos_cool_rho = np.logspace(-5, 5, 1000)
eos_cool_T = eos_cool_rho ** 0.0 * 8000.0
eos_cool_T = EAGLEfloor_cool_temperature_norm_K * eos_cool_rho **( EAGLEfloor_cool_gamma_effective - 1.0 )
eos_Jeans_rho = np.logspace(-1, 5, 1000)
eos_Jeans_T = (eos_Jeans_rho / 10 ** (-1)) ** (1.0 / 3.0) * 8000.0
eos_Jeans_T = EAGLEfloor_T * (eos_Jeans_rho / EAGLEfloor_rho_norm) ** (KS_gamma_effective - 1.0 )
# Plot the phase space diagram
figure()
subplot(111, xscale="log", yscale="log")
......@@ -135,6 +160,36 @@ savefig("rho_SFR.png", dpi=200)
########################################################################3
# Histogram of the birth density
figure()
subplot(111, xscale="linear", yscale="linear")
hist(np.log10(stars_BirthDensity),density=True,bins=100)
xlabel("${\\rm Birth Density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0)
ylabel("Probability", labelpad=-7)
savefig("BirthDensity.png", dpi=200)
# Plot of the specific star formation rate in the galaxy
rhos = 10**np.linspace(-1,np.log10(KS_high_den_thresh),100)
rhoshigh = 10**np.linspace(np.log10(KS_high_den_thresh),5,100)
P_effective = KS_press_norm * ( rhos / EAGLEfloor_rho_norm)**(KS_gamma_effective)
P_norm_high = KS_press_norm * (KS_high_den_thresh / EAGLEfloor_rho_norm)**(KS_gamma_effective)
sSFR = KS_law_norm_cgs * (Msun_p_pc2)**(-KS_law_slope) * (gamma/G_in_cgs * KS_gas_fraction *P_effective)**((KS_law_slope-1.)/2.)
KS_law_norm_high_den_cgs = KS_law_norm_cgs * (Msun_p_pc2)**(-KS_law_slope) * (gamma/G_in_cgs * KS_gas_fraction * P_norm_high)**((KS_law_slope-1.)/2.)
sSFR_high_den = KS_law_norm_high_den_cgs * ((rhoshigh/KS_high_den_thresh)**KS_gamma_effective)**((KS_law_slope_high_den-1)/2.)
figure()
subplot(111)
hist2d(np.log10(gas_nH), np.log10(gas_sSFR), bins=50,range=[[-1.5,5],[-.5,2.5]])
plot(np.log10(rhos),np.log10(sSFR)+np.log10(year_in_cgs)+9.,'k--',label='sSFR low density EAGLE')
plot(np.log10(rhoshigh),np.log10(sSFR_high_den)+np.log10(year_in_cgs)+9.,'k--',label='sSFR high density EAGLE')
xlabel("${\\rm Density}~n_{\\rm H}~[{\\rm cm^{-3}}]$", labelpad=0)
ylabel("${\\rm sSFR}~[{\\rm Gyr^{-1}}]$", labelpad=-7)
savefig("density-sSFR.png", dpi=200)
########################################################################3
# Select gas in a pillow box around the galaxy
mask = (
(gas_pos[:, 0] > -15)
......
examples/IsolatedGalaxy_starformation_25p/KS_law.py deleted 100644 → 0
+ 0
− 134
View file @ c9692302
#!/usr/bin/env python
import numpy as np
import h5py
import matplotlib.pyplot as plt
# for the moment only use a single snapshot
i = 29
# Load the snapshot
#f = h5py.File('output_%04dcopy.hdf5'%i, 'r')
f = h5py.File('output_%04d.hdf5'%i, 'r')
# Load the stars
starparticle = f['PartType4']
# Load the properties of the stars
coordinates = starparticle['Coordinates'][:,:]
newstarflag = starparticle['NewStarFlag'][:]
boxsize = f['Header'].attrs['BoxSize']
masses_star = starparticle['Masses'][:]
# Load the positions of the star particles and only take
# the new star particles
x = coordinates[:,0] - boxsize[0]/2.
y = coordinates[:,1] - boxsize[1]/2.
r = (x**2 + y**2)**.5
birthdensity = starparticle['BirthDensity'][:]
birthdensity=birthdensity[newstarflag==1]
r_new = r[newstarflag==1]
masses_new = masses_star[newstarflag==1]
# Load the gas particles
gasparticle = f['PartType0']
# Load the properties of the gas particles
coordinates_gas = gasparticle['Coordinates'][:,:]
masses_gas = gasparticle['Masses'][:]
# Load the positions of the gas particles
x_gas = coordinates_gas[:,0] - boxsize[0]/2.
y_gas = coordinates_gas[:,0] - boxsize[1]/2.
r_gas = (x_gas**2 + y_gas**2)**.5
r_array = np.linspace(0,15,100)
area = np.zeros(len(r_array)-1)
SigmaSFR = np.zeros(len(r_array)-1)
Sigmagas = np.zeros(len(r_array)-1)
massSFR_array = np.zeros(len(r_array)-1)
massGAS_array = np.zeros(len(r_array)-1)
for i in range(1,len(r_array)):
area_loop = np.pi*(r_array[i]**2 - r_array[i-1]**2)
mask = (r_new > r_array[i-1]) & (r_new < r_array[i])
mask2 = (r_gas > r_array[i-1]) & (r_gas < r_array[i])
massSFR = np.sum(masses_new[mask])
massSFR_array[i-1] = massSFR
massGAS = np.sum(masses_gas[mask2])
massGAS_array[i-1] = massGAS
#print('%d area=%1.4f, M_SFR=%1.4f, M_gas=%1.4f'%(i,area_loop,massSFR,massGAS))
area[i-1] = area_loop
SigmaSFR[i-1] = massSFR/area_loop
Sigmagas[i-1] = massGAS/area_loop
# Put SF law in correct units
SigmaSFR *= 1e10/0.29e8
Sigmagas *= 1e10/(1e3**2)
SigmaSFR[SigmaSFR==0] = 1e-6
print(Sigmagas)
Sgas_array =np.logspace(-1,3,40)
Delta_Sgas = np.log10(Sgas_array[1]) - np.log10(Sgas_array[0])
SigmaSFR_avg = np.zeros(len(Sgas_array))
Sigmagas_avg = np.zeros(len(Sgas_array))
for i in range(0,len(Sgas_array)):
totgasmass = 0.0
totnewstarmass = 0.0
area_loop = 0.0
for j in range(0,len(Sigmagas)):
#print('%d %1.4e %1.4e %1.4e %1.4e'%(j,Sigmagas[j],Sgas_array[i],10**(-Delta_Sgas),10**(Delta_Sgas)))
if (Sigmagas[j]>Sgas_array[i]*10**(-Delta_Sgas)) and (Sigmagas[j]<Sgas_array[i]*10**(Delta_Sgas)):
totgasmass += massGAS_array[j]
totnewstarmass += massSFR_array[j]
area_loop += area[j]
if area_loop==0.0:
SigmaSFR_avg[i] = 0.0
Sigmagas_avg[i] = 0.0
else:
SigmaSFR_avg[i] = totnewstarmass/area_loop
Sigmagas_avg[i] = totgasmass/area_loop
print(SigmaSFR_avg)
print(Sigmagas_avg)
# Put SF law in correct units
SigmaSFR_avg *= 1e10/0.29e8
Sigmagas_avg *= 1e10/(1e6)
SigmaSFR_avg[SigmaSFR_avg==0]=1e-6
def KS_law(Sgas,n=1.4,A=1.515):
return A*1e-4*(Sgas)**n
Sigmagaslaw = np.logspace(-1,3,100)
SigmaSFRlaw = KS_law(Sigmagaslaw)
Z = 0
ncrit = 10 # g/cm^3
ncrit = 10 # M_sun /pc^2
# plot SF law
plt.loglog(Sigmagas,SigmaSFR,'o',label='Data')
#plt.loglog(Sigmagas_avg,SigmaSFR_avg,'o',label='binned Data')
plt.loglog(Sigmagaslaw, SigmaSFRlaw,'r-.',label='KS law')
#plt.loglog(np.logspace(-1,3,40),1e-5*np.ones(40),'o')
plt.axvline(x=10,linestyle='--', color='r',label='$\Sigma_c = 10$')
plt.xlabel('$\Sigma_{gas}$')
plt.ylabel('$\Sigma_{SFR}$')
plt.legend()
plt.show()
'''
plt.hist(birthdensity,bins=100,range=[0,2.5])
plt.show()
plt.hist(x_gas,bins=100)
plt.show()
plt.hist(r,bins=100)
plt.show()
plt.hist(r_gas,bins=100)
plt.show()
'''
examples/IsolatedGalaxy_starformation_25p/StarFormationHistory.py deleted 100644 → 0
+ 0
− 36
View file @ c9692302
#!/usr/bin/env python3
import numpy as np
import h5py
import matplotlib.pyplot as plt
snaps = 30
SFH = np.zeros(snaps)
for i in range(0,snaps):
f = h5py.File('output_%04d.hdf5'%i, 'r')
header = f["Header"]
time = header.attrs['Time'][0]
print(time)
particles = f["PartType4"]
newstar = particles["NewStarFlag"][:]
SFH[i] = np.sum(newstar)
plt.plot(SFH)
plt.xlabel('Snapshot number')
plt.ylabel('Total number of stars formed')
plt.show()
f = h5py.File('output_%04d.hdf5'%i, 'r')
particles = f["PartType4"]
newstar = particles["NewStarFlag"][:]
Coordinates = particles["Coordinates"][:,:]
x = Coordinates[:,0][newstar==1]
y = Coordinates[:,1][newstar==1]
print(len(x), len(y))
plt.plot(x,y,'.')
plt.show()
examples/IsolatedGalaxy_starformation_25p/gasdensity.py deleted 100644 → 0
+ 0
− 563
View file @ c9692302
"""
Makes a movie using sphviewer and ffmpeg.
Edit low_frac and up_frac to focus on a certain view of the box.
The colour map can also be changed via colour_map.
Usage: python3 makeMovie.py CoolingHalo_
Written by James Willis (james.s.willis@durham.ac.uk)
"""
import glob
import h5py as h5
import numpy as np
import matplotlib.pyplot as plt
from multiprocessing import Pool
from tqdm import tqdm
from pylab import sys
from sphviewer.tools import QuickView
from matplotlib.animation import FuncAnimation
from matplotlib.colors import LogNorm
from scipy.optimize import curve_fit
def gen_data(filename,pixels,case):
"""
Data generation function. Reads a snapshot, filters the particle data and
creates an image using QuickView from sphviewer.
"""
# Read the data
with h5.File(filename, "r") as f:
box_size = f["/Header"].attrs["BoxSize"][0]
coordinates = f["/PartType0/Coordinates"][:,:]
masses = f["/PartType0/Masses"][:]
smoothing_lengths = f["/PartType0/SmoothingLength"][:]
coordinates_star= f["/PartType4/Coordinates"][:,:]
star_flag = f["/PartType4/NewStarFlag"][:]
SFR = f["/PartType0/SFR"][:]
# Filter data
low_frac=0.48
up_frac=0.52
absmaxz = 2 #kpc
absmaxxy = 10 #kpc
part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) <
absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) <
absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) <
absmaxz)
sfr_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) <
absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) <
absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) <
absmaxz) & (SFR > 0 )
masses = masses[part_mask]
SFR = SFR[sfr_mask]
if case ==0:
coordinates = coordinates[part_mask]
smoothing_lengths = smoothing_lengths[part_mask]
masvar = masses
elif case==1:
coordinates = coordinates[part_mask]
smoothing_lengths = smoothing_lengths[part_mask]
masvar = masses**2
elif case==2:
coordinates = coordinates[sfr_mask]
smoothing_lengths = smoothing_lengths[sfr_mask]
masvar = SFR
# Use sphviewer QuickView to generate an image
qv = QuickView(
coordinates,
mass=masvar,
r="infinity",
xsize=pixels,
ysize=pixels,
p=0, # Viewing angle theta
roll=0, # Viewing angle phi
plot=False,
logscale=False,
hsml=smoothing_lengths
)
image_data = qv.get_image()
extent = qv.get_extent()
return image_data, extent
def getnewstars(filename):
# Read the data
with h5.File(filename, "r") as f:
box_size = f["/Header"].attrs["BoxSize"][0]
coordinates= f["/PartType4/Coordinates"][:,:]
star_flag = f["/PartType4/NewStarFlag"][:]
absmaxz = 2 #kpc
absmaxxy = 10 #kpc
part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) <
absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) <
absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) <
absmaxz) & (star_flag == 1)
return coordinates[:,0][part_mask]-box_size/2., coordinates[:,1][part_mask]-box_size/2.
def getdensfrtemp(filename):
# Read the data
with h5.File(filename, "r") as f:
box_size = f["/Header"].attrs["BoxSize"][0]
coordinates= f["/PartType0/Coordinates"][:,:]
SFR = f["/PartType0/SFR"][:]
density = f["/PartType0/Density"][:]
temperature = f["/PartType0/Temperature"][:]
absmaxz = 2 #kpc
absmaxxy = 10 #kpc
part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) <
absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) <
absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) <
absmaxz)
SFR = SFR[part_mask]
density = density[part_mask]
temperature = temperature[part_mask]
return density, SFR, temperature
def getbirthdensity(filename):
# Read the data
with h5.File(filename, "r") as f:
box_size = f["/Header"].attrs["BoxSize"][0]
coordinates= f["/PartType4/Coordinates"][:,:]
birthdensity = f["/PartType4/BirthDensity"][:]
absmaxz = 2 #kpc
absmaxxy = 10 #kpc
part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) <
absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) <
absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) <
absmaxz) & (birthdensity>0)
birthdensity = birthdensity[part_mask]
return birthdensity
def getdensityrate(filename):
# Read the data
with h5.File(filename, "r") as f:
box_size = f["/Header"].attrs["BoxSize"][0]
coordinates= f["/PartType0/Coordinates"][:,:]
sfrrate = f["/PartType0/sSFR"][:]
density = f["/PartType0/Density"][:]
absmaxz = 2 #kpc
absmaxxy = 10 #kpc
part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) <
absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) <
absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) <
absmaxz) & (sfrrate>0)
sfrrate = sfrrate[part_mask]
density = density[part_mask]
return sfrrate, density
def getgasmass(filename):
# Read the data
with h5.File(filename, "r") as f:
mass = f["/PartType0/Masses"][:][0]
return mass*1e10
def getSFH(filename):
# Read the data
with h5.File(filename, "r") as f:
box_size = f["/Header"].attrs["BoxSize"][0]
coordinates= f["/PartType4/Coordinates"][:,:]
mass = f["/PartType4/Masses"][:]
flag = f["/PartType4/NewStarFlag"][:]
birth_time = f["/PartType4/Birth_time"][:]
absmaxz = 2 #kpc
absmaxxy = 10 #kpc
part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) <
absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) <
absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) <
absmaxz) & (flag==1)
birth_time = birth_time[part_mask]
mass = mass[part_mask]
histogram = np.histogram(birth_time,bins=200,weights=mass*2e4,range=[0,.1])
values = histogram[0]
xvalues =( histogram[1][:-1] +histogram[1][1:])/2.
return xvalues, values
def KS_law(Sigmagas,n=1.4,A=1.515e-4):
return A * Sigmagas**n
def specificSFR(rho,fg=.25,gamma_eos=4/3):
return .363 * (10 * fg * (rho/0.1)**gamma_eos )**.2
def specificSFRhighden(rho,gamma_eos=4/3):
return specificSFR(1e3)* (( rho/1e3 )**gamma_eos)**.5
def linearfunc(x,a,b):
return a*x + b
def makefullplot(filename,binsize=200):
# Get the Surface density from the snapshot
image_data,extent = gen_data(filename,binsize,case=0)
image_data2 = image_data[:,:]*1e10/(1e2**2)
surface_density = image_data[:,:]*1e10/(1e2**2)
if binsize==20:
image_data2 /= 1e2
surface_density /= 1e2
minvalue = np.min(image_data2[image_data2!=0])
image_data2[image_data2==0] = minvalue*.9
# Get the SFR surface density from the snapshot
image_data,extent = gen_data(filename,binsize,case=2)
image_data3 = image_data[:,:] * 1e10/1e9/(.1**2)
SFR_smoothed = image_data3/image_data2
SFR_surface_density = image_data[:,:] * 1e10/1e9/(.1**2)/surface_density
if binsize==20:
image_data3 /= 1e2
SFR_surface_density /= 1e2
SFR_smoothed /= 1e2
plot_SFR_smoothed = np.ones(np.shape(image_data3))*np.min(SFR_smoothed[SFR_smoothed!=0])*.9
plot_SFR_smoothed[image_data3!=0] = SFR_smoothed[image_data3!=0]
# Define the font used in the plots
font = {'color':'white', 'size':12}
# Make the first plot of the Gas density
fig = plt.figure(1,figsize=(16.5,16.5))
plt.subplot(3,3,1)
plt.imshow(np.log10(image_data2), extent=np.array([-10,10,-10,10]), cmap='viridis') #,vmin=-3.5, vmax=0.25)
plt.title('$\log \Sigma_{gas}$ [ $M_\odot \\rm pc^{-2}$]')
plt.plot([-7.5,-2.5],[-9.0,-9.0],linewidth=3,color='white')
plt.text(-7.5,-8.5,'5 kpc',fontdict=font)
plt.tick_params(
axis='x', # changes apply to the x-axis
which='both', # both major and minor ticks are affected
bottom=False, # ticks along the bottom edge are off
top=False, # ticks along the top edge are off
labelbottom=False) # labels along the bottom edge are off
plt.tick_params(
axis='y', # changes apply to the x-axis
which='both', # both major and minor ticks are affected
left=False, # ticks along the bottom edge are off
right=False, # ticks along the top edge are off
labelleft=False) # labels along the bottom edge are off
cbar = plt.colorbar()
plt.clim(-1.5,3.0)
#cbar.set_label('Surface density ($M_\odot pc^{-2}$)', rotation=270)
# Make the second plot of the SFR density
plt.subplot(3,3,2)
plt.imshow(np.log10(plot_SFR_smoothed), extent=np.array([-10,10,-10,10]), cmap='viridis') #,vmin=-3.5, vmax=0.25)
plt.title('$\log \Sigma_{\star}$ [ $M_\odot \\rm yr^{-1} kpc^{-2}$]')
plt.plot([-7.5,-2.5],[-9.0,-9.0],linewidth=3,color='white')
plt.text(-7.5,-8.5,'5 kpc',fontdict=font)
plt.tick_params(
axis='x', # changes apply to the x-axis
which='both', # both major and minor ticks are affected
bottom=False, # ticks along the bottom edge are off
top=False, # ticks along the top edge are off
labelbottom=False) # labels along the bottom edge are off
plt.tick_params(
axis='y', # changes apply to the x-axis
which='both', # both major and minor ticks are affected
left=False, # ticks along the bottom edge are off
right=False, # ticks along the top edge are off
labelleft=False) # labels along the bottom edge are off
cbar = plt.colorbar()
plt.clim(-5,-1)
# Make the third plot of the stars that are formed
plt.subplot(3,3,3)
xnewstar, ynewstar = getnewstars(filename)
newstarmatrix = np.histogram2d(xnewstar,ynewstar,bins=binsize)[0]
newstarmatrix[newstarmatrix==0] = 1e-1
plt.imshow(np.log10(newstarmatrix), extent=np.array([-10,10,-10,10]), cmap='viridis') #,vmin=-3.5, vmax=0.25)
plt.title('New stars')
plt.plot([-7.5,-2.5],[-9.0,-9.0],linewidth=3,color='white')
plt.text(-7.5,-8.5,'5 kpc',fontdict=font)
plt.tick_params(
axis='x', # changes apply to the x-axis
which='both', # both major and minor ticks are affected
bottom=False, # ticks along the bottom edge are off
top=False, # ticks along the top edge are off
labelbottom=False) # labels along the bottom edge are off
plt.tick_params(
axis='y', # changes apply to the x-axis
which='both', # both major and minor ticks are affected
left=False, # ticks along the bottom edge are off
right=False, # ticks along the top edge are off
labelleft=False) # labels along the bottom edge are off
cbar = plt.colorbar()
plt.clim(-1,3)
# Make the fourth plot of the KS law
plt.subplot(3,3,4)
sigma_gas = image_data2.flatten()
sigma_star = SFR_smoothed.flatten()
sigma_gas=sigma_gas[sigma_star>0]
sigma_star=sigma_star[sigma_star>0]
sigma_gas2 = surface_density.flatten()
sigma_star2 = SFR_surface_density.flatten()
#KShistogram = np.histogram2d(np.log10(sigma_gas),np.log10(sigma_star),bins=50,range=[[-1,4],[-5,2]])
#plt.imshow(np.transpose(KShistogram[0]), extent=np.array([-1,4,-5,2]))
plt.title('KS law')
fitresult = curve_fit(linearfunc,np.log10(sigma_gas),np.log10(sigma_star),p0=(1.4,np.log10(1.515e-4)),bounds=((0,-6),(3,-2)))
#print(fitresult[0], fitresult[1])
sigma_gas_range = 10**np.linspace(-1,2,100)
sigma_star_range = KS_law(sigma_gas_range)
sigma_star_one_range = KS_law(sigma_gas_range,n=0.7)
if binsize==200:
plt.hist2d(np.log10(sigma_gas),np.log10(sigma_star),range=[[-1,2],[-5,-2]],bins=50)
else:
plt.hist2d(np.log10(sigma_gas),np.log10(sigma_star),range=[[-1,2],[-5,-2]],bins=25)
#plt.hist2d(np.log10(sigma_gas2),np.log10(sigma_star2),range=[[-1,2],[-5,-2]],bins=50)
plt.plot(np.log10(sigma_gas_range), np.log10(sigma_star_range),'k--',label='EAGLE SF')
plt.plot(np.log10(sigma_gas_range), np.log10(sigma_star_one_range),'k-.',label='n=.7')
plt.plot(np.log10(sigma_gas_range), np.log10(KS_law(sigma_gas_range,n=0.6)),'k:',label='n=.6')
#plt.plot(np.log10(sigma_gas_range), linearfunc(np.log10(sigma_gas_range), *fitresult[0]),'r--',label='fit')
plt.xlabel('$\log \Sigma_{gas}$ [ $M_\odot \\rm pc^{-2}$]')
plt.ylabel('$\log \Sigma_{\star}$ [ $M_\odot \\rm yr^{-1} kpc^{-2}$]')
cbar = plt.colorbar()
if binsize==200:
plt.clim(0,700)
elif binsize==20:
plt.clim(0,40)
#plt.scatter(np.log10(sigma_gas),np.log10(sigma_star))
plt.legend()
# make a plot of the gas density vs SFR
plt.subplot(3,3,5)
critden = .1*((0.0129)/.002)**(-.64)
mu = 1.3
density, SFR, temperature = getdensfrtemp(filename)
density *= 1e10/(1e3)**3*40.4759/mu
SFR *= 1e10/1e9
SFR_plot = SFR[:]
plt.hist2d(np.log10(density[SFR_plot>0]),np.log10(SFR_plot[SFR_plot>0]), bins=50,range=[[np.log10(critden),5],[-4.4,-2.5]])
#plt.loglog(density*1e10/(1e3)**3*40.4759/mu,SFR_plot*1e10/1e9,'.')
xx = 10**np.linspace(-1,2,100)
def sfr_model(density,critden=.1,gammaeff=4./3.,mgas=1.4995e4):
return 6.9e-11 * mgas * (density/critden)**(gammaeff/5.)
rhos = 10**np.linspace(-1,4,100)
rhoshigh = 10**np.linspace(2,5,100)
massparticle =getgasmass(filename)
#plt.plot(np.log10(xx),np.log10(1e-4*xx**(4./15.)),'k--',label='Slope of EAGLE SF')
#plt.plot(np.log10(xx),np.log10(sfr_model(xx)),'r--',label='Slope of EAGLE SF')
plt.plot(np.log10(rhos),np.log10(specificSFR(rhos)*massparticle/1e9),'k--',label='SFR low density EAGLE')
plt.xlabel('$\log$ density [$\\rm cm^{-3}$]')
plt.ylabel('$\log$ Star Formation rate (SFR) [$ M_\odot \\rm yr^{-1} $]')
plt.xlim(np.log10(critden),5)
cbar = plt.colorbar()
plt.legend()
SFR_tracer = np.log10(SFR_plot[SFR_plot>0])
plt.subplot(3,3,6)
plt.scatter(np.log10(density[SFR_plot>0]),np.log10(temperature[SFR_plot>0]),s=1,c=SFR_tracer,cmap='viridis')
plt.axvline(x=np.log10(critden),linestyle='--',color='gray',linewidth=3,label='Threshold density')
plt.axvline(x=0,linestyle='--',color='k',linewidth=3,label='End of Wiersma Tables')
cbar = plt.colorbar()
plt.clim(-4.4,-2.5)
plt.scatter(np.log10(density[SFR_plot<=0]),np.log10(temperature[SFR_plot<=0]),s=1,c='gray')
plt.xlabel('$\log$ density [$\\rm cm^{-3}$]')
plt.ylabel('$\log$ temperature [$\\rm K$]')
plt.ylim(2,4)
plt.xlim(-3.5,5)
plt.legend()
plt.subplot(3,3,7)
birthdensity = getbirthdensity(filename)*404.759
plt.hist(np.log10(birthdensity),bins=50,density=True,range=(np.log10(critden),5))
plt.axvline(x=np.log10(critden),linestyle='--',color='gray',linewidth=3,label='Threshold density')
plt.xlabel('$\log$ density [$\\rm cm^{-3}$]')
plt.xlim(np.log10(critden),5)
plt.legend()
plt.subplot(3,3,8)
sfrrate, gasdensity = getdensityrate(filename)
gasdensity *= 1e10/(1e3)**3*40.4759/mu
plt.hist2d(np.log10(gasdensity),np.log10(sfrrate), bins=50,range=[[-1.5,5],[-.5,2.5]])
cbar = plt.colorbar()
rhos = 10**np.linspace(-1,4,100)
rhoshigh = 10**np.linspace(2,5,100)
plt.plot(np.log10(rhos),np.log10(specificSFR(rhos)),'k--',label='sSFR low density EAGLE')
plt.plot(np.log10(rhoshigh),np.log10(specificSFRhighden(rhoshigh)),'k--',label='sSFR high density EAGLE')
plt.xlabel('density [\\rm cm^{-3}]')
plt.ylabel('sSFR ($\\rm Gyr^{-1}$)')
plt.legend()
plt.subplot(3,3,9)
timearray, SFH = getSFH(filename)
plt.plot(timearray*1e3,SFH, label='spart birth_time')
plt.xlabel('Time (Myr)')
plt.ylabel('SFH ($\\rm M_\odot \\rm yr^{-1}$)')
plt.legend()
if binsize==200:
plt.savefig('./ksplot/'+filename+'_0.1kpc.png')
elif binsize==20:
plt.savefig('./ksplot_1kpc/'+filename+'_0.1kpc.png')
plt.close()
def getmassandsfr(filename):
# Read the data
with h5.File(filename, "r") as f:
box_size = f["/Header"].attrs["BoxSize"][0]
coordinates= f["/PartType0/Coordinates"][:,:]
SFR = f["/PartType0/SFR"][:]
coordinates_star = f["/PartType4/Coordinates"][:,:]
masses_star = f["/PartType4/Masses"][:]
absmaxz = 2 #kpc
absmaxxy = 10 #kpc
part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) <
absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) <
absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) <
absmaxz) & (SFR>0)
mask = ((coordinates_star[:,0]-box_size/2.) > -absmaxxy) & ((coordinates_star[:,0]-box_size/2.) <
absmaxxy) & ((coordinates_star[:,1]-box_size/2.) > -absmaxxy) & ((coordinates_star[:,1]-box_size/2.) <
absmaxxy) & ((coordinates_star[:,2]-box_size/2.) > -absmaxz) & ((coordinates_star[:,2]-box_size/2.) <
absmaxz)
SFR_snap = np.sum(SFR[part_mask])
mass_star = np.sum(masses_star[mask])
return SFR_snap, mass_star
def getsfrsnapwide():
time = np.arange(1,101,1)
SFR_sparticles = np.zeros(100)
SFR_gparticles = np.zeros(100)
new_sparticles = np.zeros(100)
previous_mass = 0
previous_numb = 0
for i in range(1,43):
print(i)
# Read the data
filename= 'output_%04d.hdf5'%i
with h5.File(filename, "r") as f:
box_size = f["/Header"].attrs["BoxSize"][0]
coordinates= f["/PartType0/Coordinates"][:,:]
SFR = f["/PartType0/SFR"][:]
coordinates_star = f["/PartType4/Coordinates"][:,:]
masses_star = f["/PartType4/Masses"][:]
newstar = f["/PartType4/NewStarFlag"][:]
absmaxz = 2 #kpc
absmaxxy = 10 #kpc
part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) <
absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) <
absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) <
absmaxz) & (SFR>0)
mask = ((coordinates_star[:,0]-box_size/2.) > -absmaxxy) & ((coordinates_star[:,0]-box_size/2.) <
absmaxxy) & ((coordinates_star[:,1]-box_size/2.) > -absmaxxy) & ((coordinates_star[:,1]-box_size/2.) <
absmaxxy) & ((coordinates_star[:,2]-box_size/2.) > -absmaxz) & ((coordinates_star[:,2]-box_size/2.) <
absmaxz) & (newstar==1)
SFR = SFR[part_mask]
masses_star = masses_star[mask]*1e4
newstar = newstar[mask]
total_SFR = np.sum(SFR)
total_mass = np.sum(masses_star)
new_spart = np.sum(newstar)
new_sparticles[i] = new_spart - previous_numb
SFR_sparticles[i] = total_mass - previous_mass
SFR_gparticles[i] = total_SFR*10
previous_numb = new_spart
previous_mass = total_mass
factor = masses_star[0]
print(new_sparticles)
print(len(time),len(SFR_sparticles),len(SFR_gparticles))
return time[:-1], SFR_sparticles[1:], SFR_gparticles[1:], np.sqrt(new_sparticles[1:])*factor
if __name__ == "__main__":
# Some global variables to define the whole run
input_file_name = sys.argv[1]
pixels = 1000
num_procs = 4 # No. of processors to use
colour_map = "inferno"
# Find the list of files and sort them in numerical order
file_list = sorted(glob.glob(input_file_name + "*.hdf5"), key=lambda x:
int(x[len(input_file_name):len(input_file_name)+4]))
# get the image data
#for i in [1,5,10,33,66,90]:
# print('output_%04d.hdf5'%i)
# makefullplot('output_%04d.hdf5'%i)
# makefullplot('output_%04d.hdf5'%i,binsize=20)
time, SFR1, SFR2, SFR_error = getsfrsnapwide()
time2, SFR3 = getSFH('output_%04d.hdf5'%42)
plt.plot(time,SFR1,label='Using created star particles')
print(SFR_error)
plt.fill_between(time, SFR1+SFR_error,SFR1-SFR_error,label='Uncertainty',alpha=.3)
plt.plot(time,SFR2,label='Using SFR of gas particles')
plt.plot(time2*1e3,SFR3,label='Using birth_time of star particles')
plt.xlabel('Time (Myr)')
plt.ylabel('SFH ($\\rm M_\odot \\rm yr^{-1}$)')
plt.legend()
plt.show()
'''
snapshots = 52
mass = np.zeros(snapshots)
sfr = np.zeros(snapshots)
for i in range(1,snapshots):
if i!=6:
print(i)
sfr[i], mass[i] = getmassandsfr('output_%04d.hdf5'%i)
plt.plot(np.log10(mass*1e10))
plt.xlabel('Snapshot number (Myr)')
plt.ylabel('Mass')
plt.savefig('Mass increase')
plt.close()
plt.plot(np.log10(sfr*10))
plt.xlabel('Snapshot number (Myr)')
plt.ylabel('SFR ($M_\odot / \\rm yr$)')
plt.savefig('SFR')
plt.close()
'''
for i in [1,5,10,33,66,90]:
if i!=6:
print('output_%04d.hdf5'%i)
makefullplot('output_%04d.hdf5'%i)
makefullplot('output_%04d.hdf5'%i,binsize=20)
examples/IsolatedGalaxy_starformation_25p/getIC.sh deleted 100755 → 0
+ 0
− 2
View file @ c9692302
#!/bin/bash
wget https://home.strw.leidenuniv.nl/~nobels/data/starforming_25p.hdf5
examples/IsolatedGalaxy_starformation_25p/isolated_galaxy.yml deleted 100644 → 0
+ 0
− 115
View file @ c9692302
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.9891E43 # 10^10 solar masses
UnitLength_in_cgs: 3.08567758E21 # 1 kpc
UnitVelocity_in_cgs: 1E5 # km/s
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
# Parameters for the self-gravity scheme
Gravity:
mesh_side_length: 32 # Number of cells along each axis for the periodic gravity mesh.
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.7 # Opening angle (Multipole acceptance criterion).
comoving_softening: 0.01 # Comoving softening length (in internal units).
max_physical_softening: 0.01 # Physical softening length (in internal units).
# Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.)
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 2. # The end time of the simulation (in internal units).
dt_min: 1e-9 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
basename: output # Common part of the name of output files
time_first: 0. # (Optional) Time of the first output if non-cosmological time-integration (in internal units)
delta_time: 0.001 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1e-2 # Time between statistics output
time_first: 0. # (Optional) Time of the first stats output if non-cosmological time-integration (in internal units)
# Parameters related to the initial conditions
InitialConditions:
file_name: starforming_25p.hdf5 # The file to read
periodic: 0 # Are we running with periodic ICs?
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
h_tolerance: 1e-4 # (Optional) Relative accuracy of the Netwon-Raphson scheme for the smoothing lengths.
h_max: 10. # (Optional) Maximal allowed smoothing length in internal units. Defaults to FLT_MAX if unspecified.
max_volume_change: 1.4 # (Optional) Maximal allowed change of kernel volume over one time-step.
max_ghost_iterations: 30 # (Optional) Maximal number of iterations allowed to converge towards the smoothing length.
minimal_temperature: 100 # (Optional) Minimal temperature (in internal units) allowed for the gas particles. Value is ignored if set to 0.
H_mass_fraction: 0.755 # (Optional) Hydrogen mass fraction used for initial conversion from temp to internal energy. Default value is derived from the physical constants.
H_ionization_temperature: 1e4 # (Optional) Temperature of the transition from neutral to ionized Hydrogen for primoridal gas.
EAGLECooling:
dir_name: ./coolingtables/ # Location of the Wiersma+08 cooling tables
H_reion_z: 11.5 # Redshift of Hydrogen re-ionization
He_reion_z_centre: 3.5 # Redshift of the centre of the Helium re-ionization Gaussian
He_reion_z_sigma: 0.5 # Spread in redshift of the Helium re-ionization Gaussian
He_reion_eV_p_H: 2.0 # Energy inject by Helium re-ionization in electron-volt per Hydrogen atom
# Primordial abundances
EAGLEChemistry:
init_abundance_metal: 0.0129 # Inital fraction of particle mass in *all* metals
init_abundance_Hydrogen: 0.7065 # Inital fraction of particle mass in Hydrogen
init_abundance_Helium: 0.2806 # Inital fraction of particle mass in Helium
init_abundance_Carbon: 0.00207 # Inital fraction of particle mass in Carbon
init_abundance_Nitrogen: 0.000836 # Inital fraction of particle mass in Nitrogen
init_abundance_Oxygen: 0.00549 # Inital fraction of particle mass in Oxygen
init_abundance_Neon: 0.00141 # Inital fraction of particle mass in Neon
init_abundance_Magnesium: 0.000591 # Inital fraction of particle mass in Magnesium
init_abundance_Silicon: 0.000683 # Inital fraction of particle mass in Silicon
init_abundance_Iron: 0.0011 # Inital fraction of particle mass in Iron
# Hernquist potential parameters
HernquistPotential:
useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions
position: [0.,0.,0.] # Location of centre of isothermal potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units)
idealizeddisk: 1 # Run with an idealized galaxy disk
M200: 30.0 # M200 of the galaxy disk
h: 0.704 # reduced Hubble constant (value does not specify the used units!)
concentration: 7.1 # concentration of the Halo
diskfraction: 0.040 # Disk mass fraction
bulgefraction: 0.010 # Bulge mass fraction
timestep_mult: 0.01 # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration
epsilon: 0.1 # Softening size (internal units)
# Schaye and Dalla Vecchia 2008 star formation
SchayeSF:
thresh_MinOverDens: 57.7 # The critical density contrast to form stars
thresh_temp: 1e5 # The critical temperature to form stars
fg: 0.25 # The mass fraction in gas
SchmidtLawCoeff_MSUNpYRpKPC2: 1.515e-4 # The normalization of the Kennicutt-Schmidt law
SchmidtLawExponent: 1.4 # The power law of the Kennicutt-Schmidt law
SchmidtLawHighDensExponent: 2.0 # The high density exponent for the Kennicutt-Schmidt law
SchmidtLawHighDens_thresh_HpCM3: 1e3
Schaye2004: 1 # whether to use the metallicity dependent critical star formation of Schaye (2004) (1) or not (0).
thresh_norm_HpCM3: 0.1 # Critical sf normalization to use (is not a normalization when Schaye2004=0, than it is the value.
thresh_max_norm_HpCM3: 10.0 # Maximum norm of the critical sf density
MetDep_Z0: 0.002 # Scale metallicity to use for the equation (not used for Schaye2004=0)
MetDep_SFthresh_Slope: -0.64 # Scaling of the critical density with the metallicity (not used for Schaye2004=0)
thresh_MaxPhysDensOn: 0 # Default is 0.
thresh_MaxOverDens_HpCM3: 1e5 # Density at which the SF law changes
EOS_Jeans_GammaEffective: 1.33333 # The polytropic index
EOS_Jeans_TemperatureNorm_K: 1e5 # No idea how this works
EOS_JEANS_DensityNorm_HpCM3: 0.1 # No idea what the value is.
# Parameters for the EAGLE "equation of state"
EAGLEEntropyFloor:
Jeans_density_threshold_H_p_cm3: 0.1 # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
Jeans_over_density_threshold: 10. # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in.
Jeans_temperature_norm_K: 8000 # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold expressed in Kelvin.
Jeans_gamma_effective: 1.3333333 # Slope the of the EAGLE Jeans limiter entropy floor
Cool_density_threshold_H_p_cm3: 1e-5 # Physical density above which the EAGLE Cool limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
Cool_over_density_threshold: 10. # Overdensity above which the EAGLE Cool limiter entropy floor can kick in.
Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin.
Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor
examples/IsolatedGalaxy_starformation_25p/run.sh deleted 100755 → 0
+ 0
− 3
View file @ c9692302
#!/bin/bash
../swift --threads=32 --external-gravity --self-gravity --stars --star-formation --cooling --temperature --hydro isolated_galaxy.yml 2>&1 | tee output.log
examples/IsolatedGalaxy_starformation_25p_intermediate/gasdensity.py deleted 100644 → 0
+ 0
− 563
View file @ c9692302
"""
Makes a movie using sphviewer and ffmpeg.
Edit low_frac and up_frac to focus on a certain view of the box.
The colour map can also be changed via colour_map.
Usage: python3 makeMovie.py CoolingHalo_
Written by James Willis (james.s.willis@durham.ac.uk)
"""
import glob
import h5py as h5
import numpy as np
import matplotlib.pyplot as plt
from multiprocessing import Pool
from tqdm import tqdm
from pylab import sys
from sphviewer.tools import QuickView
from matplotlib.animation import FuncAnimation
from matplotlib.colors import LogNorm
from scipy.optimize import curve_fit
def gen_data(filename,pixels,case):
"""
Data generation function. Reads a snapshot, filters the particle data and
creates an image using QuickView from sphviewer.
"""
# Read the data
with h5.File(filename, "r") as f:
box_size = f["/Header"].attrs["BoxSize"][0]
coordinates = f["/PartType0/Coordinates"][:,:]
masses = f["/PartType0/Masses"][:]
smoothing_lengths = f["/PartType0/SmoothingLength"][:]
coordinates_star= f["/PartType4/Coordinates"][:,:]
star_flag = f["/PartType4/NewStarFlag"][:]
SFR = f["/PartType0/SFR"][:]
# Filter data
low_frac=0.48
up_frac=0.52
absmaxz = 2 #kpc
absmaxxy = 10 #kpc
part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) <
absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) <
absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) <
absmaxz)
sfr_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) <
absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) <
absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) <
absmaxz) & (SFR > 0 )
masses = masses[part_mask]
SFR = SFR[sfr_mask]
if case ==0:
coordinates = coordinates[part_mask]
smoothing_lengths = smoothing_lengths[part_mask]
masvar = masses
elif case==1:
coordinates = coordinates[part_mask]
smoothing_lengths = smoothing_lengths[part_mask]
masvar = masses**2
elif case==2:
coordinates = coordinates[sfr_mask]
smoothing_lengths = smoothing_lengths[sfr_mask]
masvar = SFR
# Use sphviewer QuickView to generate an image
qv = QuickView(
coordinates,
mass=masvar,
r="infinity",
xsize=pixels,
ysize=pixels,
p=0, # Viewing angle theta
roll=0, # Viewing angle phi
plot=False,
logscale=False,
hsml=smoothing_lengths
)
image_data = qv.get_image()
extent = qv.get_extent()
return image_data, extent
def getnewstars(filename):
# Read the data
with h5.File(filename, "r") as f:
box_size = f["/Header"].attrs["BoxSize"][0]
coordinates= f["/PartType4/Coordinates"][:,:]
star_flag = f["/PartType4/NewStarFlag"][:]
absmaxz = 2 #kpc
absmaxxy = 10 #kpc
part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) <
absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) <
absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) <
absmaxz) & (star_flag == 1)
return coordinates[:,0][part_mask]-box_size/2., coordinates[:,1][part_mask]-box_size/2.
def getdensfrtemp(filename):
# Read the data
with h5.File(filename, "r") as f:
box_size = f["/Header"].attrs["BoxSize"][0]
coordinates= f["/PartType0/Coordinates"][:,:]
SFR = f["/PartType0/SFR"][:]
density = f["/PartType0/Density"][:]
temperature = f["/PartType0/Temperature"][:]
absmaxz = 2 #kpc
absmaxxy = 10 #kpc
part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) <
absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) <
absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) <
absmaxz)
SFR = SFR[part_mask]
density = density[part_mask]
temperature = temperature[part_mask]
return density, SFR, temperature
def getbirthdensity(filename):
# Read the data
with h5.File(filename, "r") as f:
box_size = f["/Header"].attrs["BoxSize"][0]
coordinates= f["/PartType4/Coordinates"][:,:]
birthdensity = f["/PartType4/BirthDensity"][:]
absmaxz = 2 #kpc
absmaxxy = 10 #kpc
part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) <
absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) <
absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) <
absmaxz) & (birthdensity>0)
birthdensity = birthdensity[part_mask]
return birthdensity
def getdensityrate(filename):
# Read the data
with h5.File(filename, "r") as f:
box_size = f["/Header"].attrs["BoxSize"][0]
coordinates= f["/PartType0/Coordinates"][:,:]
sfrrate = f["/PartType0/sSFR"][:]
density = f["/PartType0/Density"][:]
absmaxz = 2 #kpc
absmaxxy = 10 #kpc
part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) <
absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) <
absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) <
absmaxz) & (sfrrate>0)
sfrrate = sfrrate[part_mask]
density = density[part_mask]
return sfrrate, density
def getgasmass(filename):
# Read the data
with h5.File(filename, "r") as f:
mass = f["/PartType0/Masses"][:][0]
return mass*1e10
def getSFH(filename):
# Read the data
with h5.File(filename, "r") as f:
box_size = f["/Header"].attrs["BoxSize"][0]
coordinates= f["/PartType4/Coordinates"][:,:]
mass = f["/PartType4/Masses"][:]
flag = f["/PartType4/NewStarFlag"][:]
birth_time = f["/PartType4/Birth_time"][:]
absmaxz = 2 #kpc
absmaxxy = 10 #kpc
part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) <
absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) <
absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) <
absmaxz) & (flag==1)
birth_time = birth_time[part_mask]
mass = mass[part_mask]
histogram = np.histogram(birth_time,bins=200,weights=mass*2e4,range=[0,.1])
values = histogram[0]
xvalues =( histogram[1][:-1] +histogram[1][1:])/2.
return xvalues, values
def KS_law(Sigmagas,n=1.4,A=1.515e-4):
return A * Sigmagas**n
def specificSFR(rho,fg=.25,gamma_eos=4/3):
return .363 * (10 * fg * (rho/0.1)**gamma_eos )**.2
def specificSFRhighden(rho,gamma_eos=4/3):
return specificSFR(1e3)* (( rho/1e3 )**gamma_eos)**.5
def linearfunc(x,a,b):
return a*x + b
def makefullplot(filename,binsize=200):
# Get the Surface density from the snapshot
image_data,extent = gen_data(filename,binsize,case=0)
image_data2 = image_data[:,:]*1e10/(1e2**2)
surface_density = image_data[:,:]*1e10/(1e2**2)
if binsize==20:
image_data2 /= 1e2
surface_density /= 1e2
minvalue = np.min(image_data2[image_data2!=0])
image_data2[image_data2==0] = minvalue*.9
# Get the SFR surface density from the snapshot
image_data,extent = gen_data(filename,binsize,case=2)
image_data3 = image_data[:,:] * 1e10/1e9/(.1**2)
SFR_smoothed = image_data3/image_data2
SFR_surface_density = image_data[:,:] * 1e10/1e9/(.1**2)/surface_density
if binsize==20:
image_data3 /= 1e2
SFR_surface_density /= 1e2
SFR_smoothed /= 1e2
plot_SFR_smoothed = np.ones(np.shape(image_data3))*np.min(SFR_smoothed[SFR_smoothed!=0])*.9
plot_SFR_smoothed[image_data3!=0] = SFR_smoothed[image_data3!=0]
# Define the font used in the plots
font = {'color':'white', 'size':12}
# Make the first plot of the Gas density
fig = plt.figure(1,figsize=(16.5,16.5))
plt.subplot(3,3,1)
plt.imshow(np.log10(image_data2), extent=np.array([-10,10,-10,10]), cmap='viridis') #,vmin=-3.5, vmax=0.25)
plt.title('$\log \Sigma_{gas}$ [ $M_\odot \\rm pc^{-2}$]')
plt.plot([-7.5,-2.5],[-9.0,-9.0],linewidth=3,color='white')
plt.text(-7.5,-8.5,'5 kpc',fontdict=font)
plt.tick_params(
axis='x', # changes apply to the x-axis
which='both', # both major and minor ticks are affected
bottom=False, # ticks along the bottom edge are off
top=False, # ticks along the top edge are off
labelbottom=False) # labels along the bottom edge are off
plt.tick_params(
axis='y', # changes apply to the x-axis
which='both', # both major and minor ticks are affected
left=False, # ticks along the bottom edge are off
right=False, # ticks along the top edge are off
labelleft=False) # labels along the bottom edge are off
cbar = plt.colorbar()
plt.clim(-1.5,3.0)
#cbar.set_label('Surface density ($M_\odot pc^{-2}$)', rotation=270)
# Make the second plot of the SFR density
plt.subplot(3,3,2)
plt.imshow(np.log10(plot_SFR_smoothed), extent=np.array([-10,10,-10,10]), cmap='viridis') #,vmin=-3.5, vmax=0.25)
plt.title('$\log \Sigma_{\star}$ [ $M_\odot \\rm yr^{-1} kpc^{-2}$]')
plt.plot([-7.5,-2.5],[-9.0,-9.0],linewidth=3,color='white')
plt.text(-7.5,-8.5,'5 kpc',fontdict=font)
plt.tick_params(
axis='x', # changes apply to the x-axis
which='both', # both major and minor ticks are affected
bottom=False, # ticks along the bottom edge are off
top=False, # ticks along the top edge are off
labelbottom=False) # labels along the bottom edge are off
plt.tick_params(
axis='y', # changes apply to the x-axis
which='both', # both major and minor ticks are affected
left=False, # ticks along the bottom edge are off
right=False, # ticks along the top edge are off
labelleft=False) # labels along the bottom edge are off
cbar = plt.colorbar()
plt.clim(-5,-1)
# Make the third plot of the stars that are formed
plt.subplot(3,3,3)
xnewstar, ynewstar = getnewstars(filename)
newstarmatrix = np.histogram2d(xnewstar,ynewstar,bins=binsize)[0]
newstarmatrix[newstarmatrix==0] = 1e-1
plt.imshow(np.log10(newstarmatrix), extent=np.array([-10,10,-10,10]), cmap='viridis') #,vmin=-3.5, vmax=0.25)
plt.title('New stars')
plt.plot([-7.5,-2.5],[-9.0,-9.0],linewidth=3,color='white')
plt.text(-7.5,-8.5,'5 kpc',fontdict=font)
plt.tick_params(
axis='x', # changes apply to the x-axis
which='both', # both major and minor ticks are affected
bottom=False, # ticks along the bottom edge are off
top=False, # ticks along the top edge are off
labelbottom=False) # labels along the bottom edge are off
plt.tick_params(
axis='y', # changes apply to the x-axis
which='both', # both major and minor ticks are affected
left=False, # ticks along the bottom edge are off
right=False, # ticks along the top edge are off
labelleft=False) # labels along the bottom edge are off
cbar = plt.colorbar()
plt.clim(-1,3)
# Make the fourth plot of the KS law
plt.subplot(3,3,4)
sigma_gas = image_data2.flatten()
sigma_star = SFR_smoothed.flatten()
sigma_gas=sigma_gas[sigma_star>0]
sigma_star=sigma_star[sigma_star>0]
sigma_gas2 = surface_density.flatten()
sigma_star2 = SFR_surface_density.flatten()
#KShistogram = np.histogram2d(np.log10(sigma_gas),np.log10(sigma_star),bins=50,range=[[-1,4],[-5,2]])
#plt.imshow(np.transpose(KShistogram[0]), extent=np.array([-1,4,-5,2]))
plt.title('KS law')
fitresult = curve_fit(linearfunc,np.log10(sigma_gas),np.log10(sigma_star),p0=(1.4,np.log10(1.515e-4)),bounds=((0,-6),(3,-2)))
#print(fitresult[0], fitresult[1])
sigma_gas_range = 10**np.linspace(-1,2,100)
sigma_star_range = KS_law(sigma_gas_range)
sigma_star_one_range = KS_law(sigma_gas_range,n=0.7)
if binsize==200:
plt.hist2d(np.log10(sigma_gas),np.log10(sigma_star),range=[[-1,2],[-5,-2]],bins=50)
else:
plt.hist2d(np.log10(sigma_gas),np.log10(sigma_star),range=[[-1,2],[-5,-2]],bins=25)
#plt.hist2d(np.log10(sigma_gas2),np.log10(sigma_star2),range=[[-1,2],[-5,-2]],bins=50)
plt.plot(np.log10(sigma_gas_range), np.log10(sigma_star_range),'k--',label='EAGLE SF')
plt.plot(np.log10(sigma_gas_range), np.log10(sigma_star_one_range),'k-.',label='n=.7')
plt.plot(np.log10(sigma_gas_range), np.log10(KS_law(sigma_gas_range,n=0.6)),'k:',label='n=.6')
#plt.plot(np.log10(sigma_gas_range), linearfunc(np.log10(sigma_gas_range), *fitresult[0]),'r--',label='fit')
plt.xlabel('$\log \Sigma_{gas}$ [ $M_\odot \\rm pc^{-2}$]')
plt.ylabel('$\log \Sigma_{\star}$ [ $M_\odot \\rm yr^{-1} kpc^{-2}$]')
cbar = plt.colorbar()
if binsize==200:
plt.clim(0,700)
elif binsize==20:
plt.clim(0,40)
#plt.scatter(np.log10(sigma_gas),np.log10(sigma_star))
plt.legend()
# make a plot of the gas density vs SFR
plt.subplot(3,3,5)
critden = .1*((0.0129)/.002)**(-.64)
mu = 1.3
density, SFR, temperature = getdensfrtemp(filename)
density *= 1e10/(1e3)**3*40.4759/mu
SFR *= 1e10/1e9
SFR_plot = SFR[:]
plt.hist2d(np.log10(density[SFR_plot>0]),np.log10(SFR_plot[SFR_plot>0]), bins=50,range=[[np.log10(critden),5],[-4.4,-2.5]])
#plt.loglog(density*1e10/(1e3)**3*40.4759/mu,SFR_plot*1e10/1e9,'.')
xx = 10**np.linspace(-1,2,100)
def sfr_model(density,critden=.1,gammaeff=4./3.,mgas=1.4995e4):
return 6.9e-11 * mgas * (density/critden)**(gammaeff/5.)
rhos = 10**np.linspace(-1,4,100)
rhoshigh = 10**np.linspace(2,5,100)
massparticle =getgasmass(filename)
#plt.plot(np.log10(xx),np.log10(1e-4*xx**(4./15.)),'k--',label='Slope of EAGLE SF')
#plt.plot(np.log10(xx),np.log10(sfr_model(xx)),'r--',label='Slope of EAGLE SF')
plt.plot(np.log10(rhos),np.log10(specificSFR(rhos)*massparticle/1e9),'k--',label='SFR low density EAGLE')
plt.xlabel('$\log$ density [$\\rm cm^{-3}$]')
plt.ylabel('$\log$ Star Formation rate (SFR) [$ M_\odot \\rm yr^{-1} $]')
plt.xlim(np.log10(critden),5)
cbar = plt.colorbar()
plt.legend()
SFR_tracer = np.log10(SFR_plot[SFR_plot>0])
plt.subplot(3,3,6)
plt.scatter(np.log10(density[SFR_plot>0]),np.log10(temperature[SFR_plot>0]),s=1,c=SFR_tracer,cmap='viridis')
plt.axvline(x=np.log10(critden),linestyle='--',color='gray',linewidth=3,label='Threshold density')
plt.axvline(x=0,linestyle='--',color='k',linewidth=3,label='End of Wiersma Tables')
cbar = plt.colorbar()
plt.clim(-4.4,-2.5)
plt.scatter(np.log10(density[SFR_plot<=0]),np.log10(temperature[SFR_plot<=0]),s=1,c='gray')
plt.xlabel('$\log$ density [$\\rm cm^{-3}$]')
plt.ylabel('$\log$ temperature [$\\rm K$]')
plt.ylim(3.5,5.5)
plt.xlim(-3.5,5)
plt.legend()
plt.subplot(3,3,7)
birthdensity = getbirthdensity(filename)*404.759
plt.hist(np.log10(birthdensity),bins=50,density=True,range=(np.log10(critden),5))
plt.axvline(x=np.log10(critden),linestyle='--',color='gray',linewidth=3,label='Threshold density')
plt.xlabel('$\log$ density [$\\rm cm^{-3}$]')
plt.xlim(np.log10(critden),5)
plt.legend()
plt.subplot(3,3,8)
sfrrate, gasdensity = getdensityrate(filename)
gasdensity *= 1e10/(1e3)**3*40.4759/mu
plt.hist2d(np.log10(gasdensity),np.log10(sfrrate), bins=50,range=[[-1.5,5],[-.5,2.5]])
cbar = plt.colorbar()
rhos = 10**np.linspace(-1,4,100)
rhoshigh = 10**np.linspace(2,5,100)
plt.plot(np.log10(rhos),np.log10(specificSFR(rhos)),'k--',label='sSFR low density EAGLE')
plt.plot(np.log10(rhoshigh),np.log10(specificSFRhighden(rhoshigh)),'k--',label='sSFR high density EAGLE')
plt.xlabel('density [\\rm cm^{-3}]')
plt.ylabel('sSFR ($\\rm Gyr^{-1}$)')
plt.legend()
plt.subplot(3,3,9)
timearray, SFH = getSFH(filename)
plt.plot(timearray*1e3,SFH, label='spart birth_time')
plt.xlabel('Time (Myr)')
plt.ylabel('SFH ($\\rm M_\odot \\rm yr^{-1}$)')
plt.legend()
if binsize==200:
plt.savefig('./ksplot/'+filename+'_0.1kpc.png')
elif binsize==20:
plt.savefig('./ksplot_1kpc/'+filename+'_0.1kpc.png')
plt.close()
def getmassandsfr(filename):
# Read the data
with h5.File(filename, "r") as f:
box_size = f["/Header"].attrs["BoxSize"][0]
coordinates= f["/PartType0/Coordinates"][:,:]
SFR = f["/PartType0/SFR"][:]
coordinates_star = f["/PartType4/Coordinates"][:,:]
masses_star = f["/PartType4/Masses"][:]
absmaxz = 2 #kpc
absmaxxy = 10 #kpc
part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) <
absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) <
absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) <
absmaxz) & (SFR>0)
mask = ((coordinates_star[:,0]-box_size/2.) > -absmaxxy) & ((coordinates_star[:,0]-box_size/2.) <
absmaxxy) & ((coordinates_star[:,1]-box_size/2.) > -absmaxxy) & ((coordinates_star[:,1]-box_size/2.) <
absmaxxy) & ((coordinates_star[:,2]-box_size/2.) > -absmaxz) & ((coordinates_star[:,2]-box_size/2.) <
absmaxz)
SFR_snap = np.sum(SFR[part_mask])
mass_star = np.sum(masses_star[mask])
return SFR_snap, mass_star
def getsfrsnapwide():
time = np.arange(1,101,1)
SFR_sparticles = np.zeros(100)
SFR_gparticles = np.zeros(100)
new_sparticles = np.zeros(100)
previous_mass = 0
previous_numb = 0
for i in range(1,100):
print(i)
# Read the data
filename= 'output_%04d.hdf5'%i
with h5.File(filename, "r") as f:
box_size = f["/Header"].attrs["BoxSize"][0]
coordinates= f["/PartType0/Coordinates"][:,:]
SFR = f["/PartType0/SFR"][:]
coordinates_star = f["/PartType4/Coordinates"][:,:]
masses_star = f["/PartType4/Masses"][:]
newstar = f["/PartType4/NewStarFlag"][:]
absmaxz = 2 #kpc
absmaxxy = 10 #kpc
part_mask = ((coordinates[:,0]-box_size/2.) > -absmaxxy) & ((coordinates[:,0]-box_size/2.) <
absmaxxy) & ((coordinates[:,1]-box_size/2.) > -absmaxxy) & ((coordinates[:,1]-box_size/2.) <
absmaxxy) & ((coordinates[:,2]-box_size/2.) > -absmaxz) & ((coordinates[:,2]-box_size/2.) <
absmaxz) & (SFR>0)
mask = ((coordinates_star[:,0]-box_size/2.) > -absmaxxy) & ((coordinates_star[:,0]-box_size/2.) <
absmaxxy) & ((coordinates_star[:,1]-box_size/2.) > -absmaxxy) & ((coordinates_star[:,1]-box_size/2.) <
absmaxxy) & ((coordinates_star[:,2]-box_size/2.) > -absmaxz) & ((coordinates_star[:,2]-box_size/2.) <
absmaxz) & (newstar==1)
SFR = SFR[part_mask]
masses_star = masses_star[mask]*1e4
newstar = newstar[mask]
total_SFR = np.sum(SFR)
total_mass = np.sum(masses_star)
new_spart = np.sum(newstar)
new_sparticles[i] = new_spart - previous_numb
SFR_sparticles[i] = total_mass - previous_mass
SFR_gparticles[i] = total_SFR*10
previous_numb = new_spart
previous_mass = total_mass
factor = masses_star[0]
print(new_sparticles)
print(len(time),len(SFR_sparticles),len(SFR_gparticles))
return time[:-1], SFR_sparticles[1:], SFR_gparticles[1:], np.sqrt(new_sparticles[1:])*factor
if __name__ == "__main__":
# Some global variables to define the whole run
input_file_name = sys.argv[1]
pixels = 1000
num_procs = 4 # No. of processors to use
colour_map = "inferno"
# Find the list of files and sort them in numerical order
file_list = sorted(glob.glob(input_file_name + "*.hdf5"), key=lambda x:
int(x[len(input_file_name):len(input_file_name)+4]))
# get the image data
#for i in [1,5,10,33,66,90]:
# print('output_%04d.hdf5'%i)
# makefullplot('output_%04d.hdf5'%i)
# makefullplot('output_%04d.hdf5'%i,binsize=20)
'''
snapshots = 52
mass = np.zeros(snapshots)
sfr = np.zeros(snapshots)
for i in range(1,snapshots):
if i!=6:
print(i)
sfr[i], mass[i] = getmassandsfr('output_%04d.hdf5'%i)
plt.plot(np.log10(mass*1e10))
plt.xlabel('Snapshot number (Myr)')
plt.ylabel('Mass')
plt.savefig('Mass increase')
plt.close()
plt.plot(np.log10(sfr*10))
plt.xlabel('Snapshot number (Myr)')
plt.ylabel('SFR ($M_\odot / \\rm yr$)')
plt.savefig('SFR')
plt.close()
'''
for i in [1,5,10,33,66,90]:
if i!=6:
print('output_%04d.hdf5'%i)
makefullplot('output_%04d.hdf5'%i)
makefullplot('output_%04d.hdf5'%i,binsize=20)
time, SFR1, SFR2, SFR_error = getsfrsnapwide()
time2, SFR3 = getSFH('output_%04d.hdf5'%99)
plt.plot(time,SFR1,label='Using created star particles')
print(SFR_error)
plt.fill_between(time, SFR1+SFR_error,SFR1-SFR_error,label='Uncertainty',alpha=.3)
plt.plot(time,SFR2,label='Using SFR of gas particles')
plt.plot(time2*1e3,SFR3,label='Using birth_time of star particles')
plt.xlabel('Time (Myr)')
plt.ylabel('SFH ($\\rm M_\odot \\rm yr^{-1}$)')
plt.legend()
plt.show()
examples/IsolatedGalaxy_starformation_25p_intermediate/getIC.sh deleted 100755 → 0
+ 0
− 2
View file @ c9692302
#!/bin/bash
wget https://home.strw.leidenuniv.nl/~nobels/data/starforming_intermediate.hdf5
examples/IsolatedGalaxy_starformation_25p_intermediate/isolated_galaxy.yml deleted 100644 → 0
+ 0
− 114
View file @ c9692302
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.9891E43 # 10^10 solar masses
UnitLength_in_cgs: 3.08567758E21 # 1 kpc
UnitVelocity_in_cgs: 1E5 # km/s
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
# Parameters for the self-gravity scheme
Gravity:
mesh_side_length: 32 # Number of cells along each axis for the periodic gravity mesh.
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.7 # Opening angle (Multipole acceptance criterion).
comoving_softening: 0.01 # Comoving softening length (in internal units).
max_physical_softening: 0.01 # Physical softening length (in internal units).
# Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.)
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 0.1 # The end time of the simulation (in internal units).
dt_min: 1e-9 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
basename: output # Common part of the name of output files
time_first: 0. # (Optional) Time of the first output if non-cosmological time-integration (in internal units)
delta_time: 0.001 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1e-2 # Time between statistics output
time_first: 0. # (Optional) Time of the first stats output if non-cosmological time-integration (in internal units)
# Parameters related to the initial conditions
InitialConditions:
file_name: starforming_intermediate.hdf5 # The file to read
periodic: 0 # Are we running with periodic ICs?
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
h_tolerance: 1e-4 # (Optional) Relative accuracy of the Netwon-Raphson scheme for the smoothing lengths.
h_max: 10. # (Optional) Maximal allowed smoothing length in internal units. Defaults to FLT_MAX if unspecified.
max_volume_change: 1.4 # (Optional) Maximal allowed change of kernel volume over one time-step.
max_ghost_iterations: 30 # (Optional) Maximal number of iterations allowed to converge towards the smoothing length.
minimal_temperature: 100 # (Optional) Minimal temperature (in internal units) allowed for the gas particles. Value is ignored if set to 0.
H_ionization_temperature: 1e4 # (Optional) Temperature of the transition from neutral to ionized Hydrogen for primoridal gas.
EAGLECooling:
dir_name: ../coolingtables/ # Location of the Wiersma+08 cooling tables
H_reion_z: 11.5 # Redshift of Hydrogen re-ionization
He_reion_z_centre: 3.5 # Redshift of the centre of the Helium re-ionization Gaussian
He_reion_z_sigma: 0.5 # Spread in redshift of the Helium re-ionization Gaussian
He_reion_eV_p_H: 2.0 # Energy inject by Helium re-ionization in electron-volt per Hydrogen atom
# Primordial abundances
EAGLEChemistry:
init_abundance_metal: 0.0129 # Inital fraction of particle mass in *all* metals
init_abundance_Hydrogen: 0.7065 # Inital fraction of particle mass in Hydrogen
init_abundance_Helium: 0.2806 # Inital fraction of particle mass in Helium
init_abundance_Carbon: 0.00207 # Inital fraction of particle mass in Carbon
init_abundance_Nitrogen: 0.000836 # Inital fraction of particle mass in Nitrogen
init_abundance_Oxygen: 0.00549 # Inital fraction of particle mass in Oxygen
init_abundance_Neon: 0.00141 # Inital fraction of particle mass in Neon
init_abundance_Magnesium: 0.000591 # Inital fraction of particle mass in Magnesium
init_abundance_Silicon: 0.000683 # Inital fraction of particle mass in Silicon
init_abundance_Iron: 0.0011 # Inital fraction of particle mass in Iron
# Hernquist potential parameters
HernquistPotential:
useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions
position: [0.,0.,0.] # Location of centre of isothermal potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units)
idealizeddisk: 1 # Run with an idealized galaxy disk
M200: 30.0 # M200 of the galaxy disk
h: 0.704 # reduced Hubble constant (value does not specify the used units!)
concentration: 7.1 # concentration of the Halo
diskfraction: 0.040 # Disk mass fraction
bulgefraction: 0.010 # Bulge mass fraction
timestep_mult: 0.01 # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration
epsilon: 0.01 # Softening size (internal units)
# Schaye and Dalla Vecchia 2008 star formation
SchayeSF:
thresh_MinOverDens: 57.7 # The critical density contrast to form stars
thresh_temp: 1e5 # The critical temperature to form stars
fg: 0.25 # The mass fraction in gas
SchmidtLawCoeff_MSUNpYRpKPC2: 1.515e-4 # The normalization of the Kennicutt-Schmidt law
SchmidtLawExponent: 1.4 # The power law of the Kennicutt-Schmidt law
SchmidtLawHighDensExponent: 2.0 # The high density exponent for the Kennicutt-Schmidt law
SchmidtLawHighDens_thresh_HpCM3: 1e3
Schaye2004: 1 # whether to use the metallicity dependent critical star formation of Schaye (2004) (1) or not (0).
thresh_norm_HpCM3: 0.1 # Critical sf normalization to use (is not a normalization when Schaye2004=0, than it is the value.
thresh_max_norm_HpCM3: 10.0 # Maximum norm of the critical sf density
MetDep_Z0: 0.002 # Scale metallicity to use for the equation (not used for Schaye2004=0)
MetDep_SFthresh_Slope: -0.64 # Scaling of the critical density with the metallicity (not used for Schaye2004=0)
thresh_MaxPhysDensOn: 0 # Default is 0.
thresh_MaxOverDens_HpCM3: 1e5 # Density at which the SF law changes
EOS_Jeans_GammaEffective: 1.3333333 # The polytropic index
EOS_Jeans_TemperatureNorm_K: 1e5 # No idea how this works
EOS_JEANS_DensityNorm_HpCM3: 0.1 # No idea what the value is.
# Parameters for the EAGLE "equation of state"
EAGLEEntropyFloor:
Jeans_density_threshold_H_p_cm3: 0.1 # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
Jeans_over_density_threshold: 10. # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in.
Jeans_temperature_norm_K: 8000 # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold expressed in Kelvin.
Jeans_gamma_effective: 1.3333333 # Slope the of the EAGLE Jeans limiter entropy floor
Cool_density_threshold_H_p_cm3: 1e-5 # Physical density above which the EAGLE Cool limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
Cool_over_density_threshold: 10. # Overdensity above which the EAGLE Cool limiter entropy floor can kick in.
Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin.
Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor
examples/IsolatedGalaxy_starformation_25p_intermediate/run.sh deleted 100755 → 0
+ 0
− 8
View file @ c9692302
#!/bin/bash
if [ ! -e starforming_intermediate.hdf5 ]
then
./getIC.sh
fi
../swift --threads=32 --external-gravity --self-gravity --stars --star-formation --cooling --temperature --hydro isolated_galaxy.yml 2>&1 | tee output.log
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment