Skip to content

Commit

Permalink
Updating several plots
Browse files Browse the repository at this point in the history
Nothing major though
  • Loading branch information
cdplagos committed Mar 31, 2022
1 parent 84d9133 commit c82c76d
Show file tree
Hide file tree
Showing 11 changed files with 783 additions and 69 deletions.
197 changes: 197 additions & 0 deletions standard_plots/age_size_disk.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,197 @@
#
# ICRAR - International Centre for Radio Astronomy Research
# (c) UWA - The University of Western Australia, 2018
# Copyright by UWA (in the framework of the ICRAR)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
#
"""Size plots"""

import functools

import numpy as np

import os
import common
import utilities_statistics as us

# Initialize arguments
zlist = [0]

##################################
#Constants
RExp = 1.67
MpcToKpc = 1e3
G = 4.299e-9 #Gravity constant in units of (km/s)^2 * Mpc/Msun

mlow = 6.5
mupp = 12.5
dm = 0.2
mbins = np.arange(mlow,mupp,dm)
xmf = mbins + dm/2.0

def prepare_data(hdf5_data, age_bins, sfh, delta_t, LBT, index, disk_size, disk_size_age):

(h0, _, mdisk, mbulge, rdisk, typeg, age) = hdf5_data
(bulge_diskins_hist, bulge_mergers_hist, disk_hist) = sfh

bin_it = functools.partial(us.wmedians, xbins=xmf)

BTthresh = 0.5

mstar = mdisk + mbulge
ind = np.where(mstar > 0)
mdisk = mdisk[ind]
mbulge = mbulge[ind]
rdisk = rdisk[ind]
typeg = typeg[ind]
age = age[ind]
mstar = mstar[ind]

age_disk = np.zeros(shape = len(mdisk))
for j in range(0,len(mdisk)):
sfrtimesage = np.sum(disk_hist[j,:] * delta_t[:] * LBT[:])
massformed = np.sum(disk_hist[j,:] * delta_t[:])
age_disk[j] = sfrtimesage / massformed

ind = np.where((mdisk > 0) & (mdisk / mstar >= BTthresh))
disk_size[index,:] = bin_it(x=np.log10(mdisk[ind]) - np.log10(float(h0)),
y=np.log10(rdisk[ind]*MpcToKpc) - np.log10(float(h0)))

BtoT = mbulge / mstar
ind = np.where((mdisk > 0))
print("#log10(mdisk/Msun) log10(rdisk/kpc) age_disk age_gal B/T")
for a,b,c,d,e in zip(np.log10(mdisk[ind]) - np.log10(float(h0)), np.log10(rdisk[ind]*MpcToKpc) - np.log10(float(h0)), age_disk[ind], 13.7969 - age[ind], BtoT[ind]):
print(a,b,c,d,e)

age_disk = 13.7969 - age
age_meds = np.zeros(shape = len(age_bins)-1)
for j in range(0,len(age_bins)-1):
ind = np.where((mdisk > 0) & (age_disk >= age_bins[j]) & (age_disk < age_bins[j+1]) & (mdisk / mstar >= BTthresh))
disk_size_age[index,j,:] = bin_it(x=np.log10(mdisk[ind]) - np.log10(float(h0)),
y=np.log10(rdisk[ind]*MpcToKpc) - np.log10(float(h0)))
age_meds[j] = np.mean(age_disk[ind])

ind = np.where((mdisk >= 5e7) & (mdisk / mstar >= BTthresh))
return(np.log10(mdisk[ind]) - np.log10(float(h0)), np.log10(rdisk[ind]*MpcToKpc) - np.log10(float(h0)), age_disk[ind], age_meds)

def plot_age_disk(plt, outdir, obsdir, mdisk_z0, rdisk_z0, age_z0, disk_size, disk_size_age, age_bins, age_meds):

fig = plt.figure(figsize=(5,4.5))
xtit = "$\\rm log_{10} (\\rm M_{\\rm disk}/M_{\odot})$"
ytit = "$\\rm log_{10} (\\rm r_{\\star,disk}/kpc)$"
xmin, xmax, ymin, ymax = 8, 12, -0.5, 2

# LTG ##################################
ax = fig.add_subplot(111)
plt.subplots_adjust(bottom=0.15, left=0.15)

common.prepare_ax(ax, xmin, xmax, ymin, ymax, xtit, ytit, locators=(0.1, 1, 0.1, 1))

im = ax.hexbin(mdisk_z0, rdisk_z0, age_z0, xscale='linear', yscale='linear', gridsize=(20,20), cmap='magma', mincnt=30)
#cbar_ax = fig.add_axes([0.86, 0.15, 0.025, 0.7])
cbar = fig.colorbar(im)#, cax=cbar_ax)
cbar.ax.set_ylabel('mean stellar age')

#Predicted size-mass for disks in disk=dominated galaxies
ind = np.where(disk_size[0,0,:] != 0)
xplot = xmf[ind]
yplot = disk_size[0,0,ind]
errdn = disk_size[0,1,ind]
errup = disk_size[0,2,ind]
ax.errorbar(xplot,yplot[0],yerr=[errdn[0],errup[0]], ls='None', mfc='None', ecolor = 'k', mec='k',marker='o',label="Shark")

def plot_robotham_size_relation_use_bins(ax, age_meds):
for age in age_meds:
y1 = 0.35 * (xmf - 10) - 0.06 * (age) + 1.01
ax.plot(xmf, y1, linestyle='dotted', color='k')


def plot_lange_relation(ax):
#Lange et al. (2016)
a = 5.56
aerr = 1.745
b = 0.274
ind = np.where(xmf < 10.3)
rL16 = np.log10(a*pow((pow(10.0,xmf)/1e10),b))
ax.plot(xmf[ind],rL16[ind],'b', linestyle='solid', label ='L16 disks')

m,r = common.load_observation(obsdir, 'SizesAndAM/rdisk_L16.dat', [0,1])
ax.plot(m[0:36], r[0:36], linestyle='dotted',color='b',label="50th, 68th, 90th")
ax.plot(m[38:83], r[38:83], linestyle='dotted',color='b')
ax.plot(m[85:128], r[85:129], linestyle='dotted',color='b')

plot_robotham_size_relation_use_bins(ax, np.array([3, 6, 9, 12]))
plot_lange_relation(ax)

common.prepare_legend(ax, ['b','b','k','r'], loc=2)
common.savefig(outdir, fig, 'age_disks.pdf')


fig = plt.figure(figsize=(5,4.5))
xtit = "$\\rm log_{10} (\\rm M_{\\rm disk}/M_{\odot})$"
ytit = "$\\rm log_{10} (\\rm r_{\\star,disk}/kpc)$"
xmin, xmax, ymin, ymax = 8, 12, -0.5, 2

# LTG ##################################
ax = fig.add_subplot(111)
plt.subplots_adjust(bottom=0.15, left=0.15)

common.prepare_ax(ax, xmin, xmax, ymin, ymax, xtit, ytit, locators=(0.1, 1, 0.1, 1))

for j in range(0,len(age_bins)-1):
ind = np.where(disk_size_age[0,j,0,:] != 0)
xplot = xmf[ind]
yplot = disk_size_age[0,j,0,ind]
errdn = disk_size_age[0,j,1,ind]
errup = disk_size_age[0,j,2,ind]
ax.errorbar(xplot,yplot[0],yerr=[errdn[0],errup[0]], ls='None', mfc='None', marker='o', label='age=[%s,%s]' % (str(age_bins[j]), str(age_bins[j+1])))

plot_robotham_size_relation_use_bins(ax, age_meds)
common.prepare_legend(ax, ['k','k','k','k','k','k','k'], loc=2)
common.savefig(outdir, fig, 'age_bins_disks.pdf')


def main(modeldir, outdir, redshift_table, subvols, obsdir):

plt = common.load_matplotlib()
fields = {'galaxies': ('mstars_disk', 'mstars_bulge', 'rstar_disk', 'type',
'mean_stellar_age')}

sfh_fields = {'bulges_diskins': ('star_formation_rate_histories'),
'bulges_mergers': ('star_formation_rate_histories'),
'disks': ('star_formation_rate_histories')}

# Loop over redshift and subvolumes

age_bins=[0,3,6,10,14]
disk_size = np.zeros(shape = (len(zlist), 3, len(xmf)))
disk_size_age = np.zeros(shape = (len(zlist), len(age_bins)-1, 3, len(xmf)))

for index, snapshot in enumerate(redshift_table[zlist]):
hdf5_data = common.read_data(modeldir, snapshot, fields, subvols)
sfh, delta_t, LBT = common.read_sfh(modeldir, 199, sfh_fields, subvols)

(mdisk, rdisk, age, age_meds) = prepare_data(hdf5_data, age_bins, sfh, delta_t, LBT, index, disk_size, disk_size_age)
if(index == 0):
mdisk_z0 = mdisk
rdisk_z0 = rdisk
age_z0 = age
age_meds_z0 = age_meds

print("Will produce plots")
plot_age_disk(plt, outdir, obsdir, mdisk_z0, rdisk_z0, age_z0, disk_size, disk_size_age, age_bins, age_meds)

if __name__ == '__main__':
main(*common.parse_args())
23 changes: 23 additions & 0 deletions standard_plots/common.py
Original file line number Diff line number Diff line change
Expand Up @@ -274,6 +274,29 @@ def read_data(model_dir, snapshot, fields, subvolumes, include_h0_volh=True):

return list(data.values())

def read_data_ext(model_dir, snapshot, fields, subvolumes, dustmodel):
"""Read the galaxies.hdf5 file for the given model/snapshot/subvolume"""

data = collections.OrderedDict()
for idx, subv in enumerate(subvolumes):

fname = os.path.join(model_dir, str(snapshot), str(subv), 'extinction-' + dustmodel + '.hdf5')
logger.info('Reading galaxies data from %s', fname)
with h5py.File(fname, 'r') as f:
for gname, dsnames in fields.items():
group = f[gname]
for dsname in dsnames:
full_name = '%s/%s' % (gname, dsname)
l = data.get(full_name, None)
if l is None:
l = group[dsname][()]
else:
l = np.concatenate([l, group[dsname][()]])
data[full_name] = l

return list(data.values())


def read_sfh(model_dir, snapshot, fields, subvolumes, include_h0_volh=True):
"""Read the galaxies.hdf5 file for the given model/snapshot/subvolume"""

Expand Down
2 changes: 1 addition & 1 deletion standard_plots/extinction_calculation.py
Original file line number Diff line number Diff line change
Expand Up @@ -254,7 +254,7 @@ def main(model_dir, output_dir, redshift_table, subvols, obs_dir):
#zlist = np.arange(2,10,0.25)
#zlist = (0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75, 0.85, 0.95, 0.1, 0.2, 0.3, 0.4, 0.6, 0.7, 0.8, 0.9, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 0, 0.25, 0.5, 1, 2, 3, 4, 6, 8, 9, 10)

zlist = [0.194739, 0.254144, 0.359789, 0.450678, 0.8, 0.849027, 0.9, 1.20911, 1.28174, 1.39519, 1.59696, 2.00392, 2.47464723643932, 2.76734390952347, 3.01916, 3.21899984389701, 3.50099697082904, 3.7248038025221, 3.95972, 4.465197621546, 4.73693842543988] #[5.02220991014863, 5.52950356184419, 5.96593, 6.55269895697227, 7.05756323172746, 7.45816170313544, 8.02352, 8.94312532315157, 9.95655]
zlist = [9.95655] #0.194739, 0.254144, 0.359789, 0.450678, 0.8, 0.849027, 0.9, 1.20911, 1.28174, 1.39519, 1.59696, 2.00392, 2.47464723643932, 2.76734390952347, 3.01916, 3.21899984389701, 3.50099697082904, 3.7248038025221, 3.95972, 4.465197621546, 4.73693842543988] #[5.02220991014863, 5.52950356184419, 5.96593, 6.55269895697227, 7.05756323172746, 7.45816170313544, 8.02352, 8.94312532315157, 9.95655]
#[0.016306640039433, 0.066839636933135, 0.084236502339783, 0.119886040396529, 0.138147164704691, 0.175568857770275, 0.214221447279112, 0.23402097095238, 0.274594901875312, 0.316503156974571]

plt = common.load_matplotlib()
Expand Down
39 changes: 24 additions & 15 deletions standard_plots/global_quantities.py
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@

def prepare_data(hdf5_data, redshifts):

(h0, volh, _, mHI, mH2, mcold, mcold_metals, mhot, meje, mstar,
(h0, volh, _, mHI, mH2, mcold, mcold_metals, mhot, meje, mlost, mcreated, mstar,
mstar_burst_mergers, mstar_burst_diskins, mBH, sfrdisk, sfrburst,
mDM, mcold_halo, number_major_mergers, number_minor_mergers,
number_disk_instabil, max_smbh) = hdf5_data
Expand All @@ -65,7 +65,9 @@ def prepare_data(hdf5_data, redshifts):
#Add up cold halo component to hot gas.
mhot = mhot + mcold_halo

massbar = mcold+mhot+meje+mstar+mBH
print(mlost)
massbar = mcold+mhot+meje+mstar+mBH+mlost

sfr = sfrall / volh / GyrToYr
sfrd = sfrdisk / volh / GyrToYr
sfrb = sfrburst / volh / GyrToYr
Expand Down Expand Up @@ -130,6 +132,8 @@ def prepare_data(hdf5_data, redshifts):
mbar_dm_plot = np.zeros(shape = len(mcold))
mHI_dm_plot = np.zeros(shape = len(mcold))
mH2_dm_plot = np.zeros(shape = len(mcold))
mlost_dm_plot = np.zeros(shape = len(mcold))
mcreated_dm_plot = np.zeros(shape = len(mcold))

ind = np.where(mDM > 0.0)
mcold_dm_plot[ind] = np.log10(mcold[ind]/(mDM[ind]+massbar[ind]))
Expand All @@ -139,12 +143,14 @@ def prepare_data(hdf5_data, redshifts):
mbar_dm_plot[ind] = np.log10(massbar[ind]/(mDM[ind]+massbar[ind]))
mHI_dm_plot[ind] = np.log10(mHI[ind]/(mDM[ind]+massbar[ind]))
mH2_dm_plot[ind] = np.log10(mH2[ind]/(mDM[ind]+massbar[ind]))

mlost_dm_plot[ind] = np.log10(mlost[ind]/(mDM[ind]+massbar[ind]))
mcreated_dm_plot[ind] = np.log10(mcreated[ind]/(mDM[ind]+massbar[ind]))

return (mstar_plot, mcold_plot, mhot_plot, meje_plot,
mstar_dm_plot, mcold_dm_plot, mhot_dm_plot, meje_dm_plot, mbar_dm_plot,
sfr, sfrd, sfrb, mstarden, mstarbden_mergers, mstarbden_diskins, sfre, sfreH2, mhrat,
mHI_plot, mH2_plot, mH2den, mdustden, omegaHI, mdustden_mol, mcoldden, mhotden, mejeden,
history_interactions, mDMden)
history_interactions, mDMden, mlost_dm_plot, mcreated_dm_plot)

def plot_mass_densities(plt, outdir, obsdir, h0, redshifts, mstar, mcold, mhot, meje, mstarden, mcoldden, mhotden, mejeden):

Expand Down Expand Up @@ -334,7 +340,7 @@ def plot_mass_densities(plt, outdir, obsdir, h0, redshifts, mstar, mcold, mhot,
common.savefig(outdir, fig, "global_hotgas.pdf")


def plot_baryon_fractions(plt, outdir, redshifts, mstar_dm, mcold_dm, mhot_dm, meje_dm, mbar_dm):
def plot_baryon_fractions(plt, outdir, redshifts, mstar_dm, mcold_dm, mhot_dm, meje_dm, mbar_dm, mlost_dm, mcreated_dm):

fig = plt.figure(figsize=(9.5,9.5))

Expand All @@ -348,14 +354,17 @@ def plot_baryon_fractions(plt, outdir, redshifts, mstar_dm, mcold_dm, mhot_dm, m
ax.plot(redshifts, mcold_dm, 'b', label='ISM gas')
ax.plot(redshifts, mhot_dm, 'r', label='halo gas')
ax.plot(redshifts, meje_dm, 'g', label='ejected gas')
ax.plot(redshifts, mlost_dm, 'orange', label='lost')
ax.plot(redshifts, mcreated_dm, 'orange', linestyle='dashed', label='created')

ax.plot(redshifts, mbar_dm, 'm', label='total baryons')

yplot = [0.1866920152, 0.1866920152]
xplot = [0, 10]

ax.plot(xplot, np.log10(yplot), 'k', linestyle='dashed', label ='Universal $f_{\\rm baryon}$')

common.prepare_legend(ax, ['k','b','r','g','m','k'])
common.prepare_legend(ax, ['k','b','r','g','orange','orange','m','k'])
common.savefig(outdir, fig, "baryon_frac.pdf")


Expand Down Expand Up @@ -406,9 +415,9 @@ def plot_cosmic_sfr(plt, outdir, obsdir, redshifts, h0, sfr, sfrd, sfrb, history
#note that only h^2 is needed because the volume provides h^3, and the SFR h^-1.
ind = np.where(sfr > 0)
ax.plot(redshifts[ind], np.log10(sfr[ind]*pow(h0,2.0)), 'k', linewidth=1, label ='total')
print("SFR density of the Universe")
for a,b in zip(redshifts[ind], np.log10(sfr[ind]*pow(h0,2.0))):
print(a,b)
#print("SFR density of the Universe")
#for a,b in zip(redshifts[ind], np.log10(sfr[ind]*pow(h0,2.0))):
# print(a,b)

ind = np.where(sfrd > 0)
ax.plot(redshifts[ind], np.log10(sfrd[ind]*pow(h0,2.0)), 'b', linestyle='dashed', linewidth=1, label ='quiescent')
Expand Down Expand Up @@ -666,9 +675,9 @@ def plot_omega_h2(plt, outdir, obsdir, redshifts, h0, mH2den):
ind = np.where(mH2den > 0)
ax.plot(us.look_back_time(redshifts[ind]), np.log10(mH2den[ind]*pow(h0,2.0)) + np.log10(XH), 'r')

print("H2 density of the Universe")
for a,b in zip(redshifts[ind], np.log10(mH2den[ind]*pow(h0,2.0))):
print(a,b)
#print("H2 density of the Universe")
#for a,b in zip(redshifts[ind], np.log10(mH2den[ind]*pow(h0,2.0))):
# print(a,b)

z, h2_modelvar = common.load_observation(obsdir, 'Models/SharkVariations/Global_OtherModels.dat', [0, 2])
h2_modelvar_burst3 = h2_modelvar[0:179]
Expand Down Expand Up @@ -832,7 +841,7 @@ def main(modeldir, outdir, redshift_table, subvols, obsdir):

plt = common.load_matplotlib()
fields = {'global': ('redshifts', 'm_hi', 'm_h2', 'mcold', 'mcold_metals',
'mhot_halo', 'mejected_halo', 'mstars', 'mstars_bursts_mergers', 'mstars_bursts_diskinstabilities',
'mhot_halo', 'mejected_halo', 'mbar_lost', 'mbar_created', 'mstars', 'mstars_bursts_mergers', 'mstars_bursts_diskinstabilities',
'm_bh', 'sfr_quiescent', 'sfr_burst', 'm_dm', 'mcold_halo', 'number_major_mergers',
'number_minor_mergers', 'number_disk_instabilities', 'smbh_maximum')}

Expand All @@ -859,10 +868,10 @@ def main(modeldir, outdir, redshift_table, subvols, obsdir):
mstar_dm_plot, mcold_dm_plot, mhot_dm_plot, meje_dm_plot, mbar_dm_plot,
sfr, sfrd, sfrb, mstarden, mstarbden_mergers, mstarbden_diskins, sfre, sfreH2, mhrat,
mHI_plot, mH2_plot, mH2den, mdustden, omegaHI, mdustden_mol, mcoldden, mhotden,
mejeden, history_interactions, mDMden) = prepare_data(hdf5_data, redshifts)
mejeden, history_interactions, mDMden, mlost_dm_plot, mcreated_dm_plot) = prepare_data(hdf5_data, redshifts)

plot_mass_densities(plt, outdir, obsdir, h0, redshifts, mstar_plot, mcold_plot, mhot_plot, meje_plot, mstarden, mcoldden, mhotden, mejeden)
plot_baryon_fractions(plt, outdir, redshifts, mstar_dm_plot, mcold_dm_plot, mhot_dm_plot, meje_dm_plot, mbar_dm_plot)
plot_baryon_fractions(plt, outdir, redshifts, mstar_dm_plot, mcold_dm_plot, mhot_dm_plot, meje_dm_plot, mbar_dm_plot, mlost_dm_plot, mcreated_dm_plot)
plot_cosmic_sfr(plt, outdir, obsdir, redshifts, h0, sfr, sfrd, sfrb, history_interactions, mDMden)
plot_stellar_mass_cosmic_density(plt, outdir, obsdir, redshifts, h0, mstarden, mstarbden_mergers, mstarbden_diskins)
plot_sft_efficiency(plt, outdir, redshifts, sfre, sfreH2, mhrat)
Expand Down
Loading

0 comments on commit c82c76d

Please sign in to comment.