Skip to content

Commit

Permalink
Adding new data and small changes to python scripts
Browse files Browse the repository at this point in the history
  • Loading branch information
cdplagos committed Aug 14, 2023
1 parent 6c773e2 commit 94a8a3b
Show file tree
Hide file tree
Showing 7 changed files with 6,912 additions and 71 deletions.
5 changes: 5 additions & 0 deletions data/SMHM/Taylor20.dat
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
#log10(stellar mass/Msun) Mstar_low Mstar_high Mhalo / 1e12Msun Mhalo_low Mhalo_high
10.3436 10.3028 10.3894 0.284179 0.150448 0.534925
10.4374 10.3916 10.4870 0.284179 0.158806 0.509851
10.5389 10.4859 10.5908 0.743881 0.534925 1.04060
10.6421 10.5886 10.7000 0.229851 0.121194 0.459701
6,665 changes: 6,665 additions & 0 deletions data/SizeMass/ProFuse_cat_BestModel_v1.0.csv

Large diffs are not rendered by default.

193 changes: 152 additions & 41 deletions standard_plots/angularmomentum.py

Large diffs are not rendered by default.

29 changes: 19 additions & 10 deletions standard_plots/global_quantities.py
Original file line number Diff line number Diff line change
Expand Up @@ -391,6 +391,7 @@ def plot_cosmic_sfr(plt, outdir, obsdir, redshifts, h0, sfr, sfrd, sfrb, history
#D'Silva+23 (Chabrier IMF)
sfrD23, sfrD23errup, sfrD23errdn, zD23, zD23errup, zD23errdn = common.load_observation(obsdir, 'Global/DSilva23_sfr.dat', [0,1,2,3,4,5])


#Adams23 (Chabrier IMF)
zA23, sfrA23, sfrA23errdn, sfrA23errup = common.load_observation(obsdir, 'Global/Adams23_CSFRDCompilation.dat', [0,1,2,3])
sfrA23errdn = sfrA23 - sfrA23errdn #make them relative errors
Expand Down Expand Up @@ -531,20 +532,28 @@ def plot_cosmic_sfr(plt, outdir, obsdir, redshifts, h0, sfr, sfrd, sfrb, history
sfrM14errdn = abs(sfrM14errdn)
hobs = 0.7
sfrM14 = sfrM14 + np.log10(pow(hobs/h0, 2.0)) + np.log10(0.63)
ax.errorbar((reddnM14 + redupM14) / 2.0, sfrM14, xerr=abs(redupM14-reddnM14)/2.0, yerr=[sfrM14errdn, sfrM14errup], ls='None', mfc='None', ecolor = 'grey', mec='grey',marker='D', markersize=1.5, label='Madau+14')
#ax.errorbar((reddnM14 + redupM14) / 2.0, sfrM14, xerr=abs(redupM14-reddnM14)/2.0, yerr=[sfrM14errdn, sfrM14errup], ls='None', mfc='None', ecolor = 'grey', mec='grey',marker='D', markersize=1.5, label='Madau+14')

ax.errorbar(zD23, yobsD23, xerr=[zD23errdn, zD23errup], yerr=[sfrD23errdn, sfrD23errup], ls='None', mfc='None', ecolor = 'navy', mec='navy',marker='o', label ='D\'Silva+23')
ax.errorbar(zA23, yobsA23, yerr=[sfrA23errdn, sfrA23errup], ls='None', mfc='None', ecolor = 'darkgreen', mec='darkgreen',marker='s', label = 'Adams+23')
#Driver (Chabrier IMF)
redD17d, redD17u, sfrD17, err1, err2, err3, err4 = common.load_observation(obsdir, 'Global/Driver18_sfr.dat', [0,1,2,3,4,5,6])
hobs = 0.7
xobsD17 = (redD17d+redD17u)/2.0
yobsD17 = sfrD17 + np.log10(hobs/h0)
errD17 = yobsD17*0. - 999.
errD17 = np.sqrt(pow(err1,2.0)+pow(err2,2.0)+pow(err3,2.0)+pow(err4,2.0))
ax.errorbar(xobsD17, yobsD17, yerr=[errD17,errD17], ls='None', mfc='None', ecolor = 'darkorange', mec='darkorange',marker='o', label = 'Driver+18')


#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 (v2.0)')

ind = np.where(sfrd > 0)
ax.plot(redshifts[ind], np.log10(sfrd[ind]*pow(h0,2.0)), 'b', linestyle='solid', linewidth=1, label ='disks (v2.0)')
ax.plot(redshifts[ind], np.log10(sfrd[ind]*pow(h0,2.0)), 'b', linestyle='solid', linewidth=1, label ='disks')
ind = np.where(sfrb > 0)
ax.plot(redshifts[ind], np.log10(sfrb[ind]*pow(h0,2.0)),'r', linestyle='solid', linewidth=1, label ='bursts (v2.0)')
ax.plot(redshifts[ind], np.log10(sfrb[ind]*pow(h0,2.0)),'r', linestyle='solid', linewidth=1, label ='bursts')
#print("#will print SFR density")
#print("#redshift LBG/Gyr SFRD_tot[Msun/yr/Mpc^3] SFRD_disk[Msun/yr/Mpc^3] SFRD_bulges[Msun/yr/Mpc^3]")
#for a,b,c,d,e in zip(redshifts, us.look_back_time(redshifts), np.log10(sfr*pow(h0,2.0)), np.log10(sfrd*pow(h0,2.0)), np.log10(sfrb*pow(h0,2.0))):
Expand All @@ -557,7 +566,7 @@ def plot_cosmic_sfr(plt, outdir, obsdir, redshifts, h0, sfr, sfrd, sfrb, history
ax.plot(zin, sfr_l18d, 'b', linewidth=1,linestyle='dashed') #, label ='disks (L18)')
ax.plot(zin, sfr_l18b, 'r', linewidth=1,linestyle='dashed') #, label ='bursts (L18)')

common.prepare_legend(ax, ['k','b','r','k','grey','navy', 'darkgreen'], loc=3) #bbox_to_anchor=(0.52, 0.47))
common.prepare_legend(ax, ['k','b','r','k','navy', 'darkgreen','DarkOrange'], loc=3) #bbox_to_anchor=(0.52, 0.47))

common.savefig(outdir, fig, "cosmic_sfr_compL18.pdf")

Expand Down Expand Up @@ -664,10 +673,10 @@ def plot_stellar_mass_cosmic_density(plt, outdir, obsdir, redshifts, h0, mstarde


ind = np.where(mstardisk > 0)
ax.plot(redshifts[ind],np.log10(mstardisk[ind]*pow(h0,2.0)), 'b', linestyle='solid', label='disks (v2.0)')
ax.plot(redshifts[ind],np.log10(mstardisk[ind]*pow(h0,2.0)), 'b', linestyle='solid', label='disks')

ind = np.where(mstarbden_mergers + mstarbden_diskins > 0)
ax.plot(redshifts[ind],np.log10((mstarbden_mergers[ind]+mstarbden_diskins[ind])*pow(h0,2.0)), 'r', linestyle='solid', label='bursts (v2.0)')
ax.plot(redshifts[ind],np.log10((mstarbden_mergers[ind]+mstarbden_diskins[ind])*pow(h0,2.0)), 'r', linestyle='solid', label='bulges')


zin, sd_l18, sd_l18d, sd_l18b = common.load_observation(obsdir, 'Models/SharkVariations/Global_SMD_Lagos18.dat', [0, 2, 3, 4])
Expand All @@ -692,14 +701,14 @@ def plot_stellar_mass_cosmic_density(plt, outdir, obsdir, redshifts, h0, mstarde

err = yobs*0. - 999.
err = np.sqrt(pow(err1,2.0)+pow(err2,2.0)+pow(err3,2.0)+pow(err4,2.0))
ax.errorbar(xobs, yobs, yerr=[err,err], ls='None', mfc='None', ecolor = 'navy', mec='navy',marker='o', label="Driver+18")
ax.errorbar(xobs, yobs, yerr=[err,err], ls='None', mfc='None', ecolor = 'darkorange', mec='darkorange',marker='o', label="Driver+18")

zdnW22, zupW22, rhoW22, errrhoW22u, errrhoW22dn = common.load_observation(obsdir, 'Global/Weaver22_SMD.dat', [0,1,2,3,4])
rhoW22_cosmocorr = rhoW22 * hobs/h0 * 1e7
zW22 = (zupW22 + zdnW22)/2.0
errrhoW22u = np.log10(rhoW22 + errrhoW22u) - np.log10(rhoW22)
errrhoW22dn = np.log10(rhoW22) - np.log10(rhoW22 - errrhoW22dn)
ax.errorbar(zW22, np.log10(rhoW22_cosmocorr), xerr=[zW22 - zdnW22, zupW22 - zW22], yerr=[errrhoW22dn, errrhoW22u], ls='None', ecolor = 'grey', mec='grey', marker='*', label='Weaver+22')
ax.errorbar(zW22, np.log10(rhoW22_cosmocorr), xerr=[zW22 - zdnW22, zupW22 - zW22], yerr=[errrhoW22dn, errrhoW22u], ls='None', ecolor = 'navy', mec='navy', marker='*', label='Weaver+22')

zdnS23, zmidS23, zupS23, rhoS23, rhoS23_dn_s, rhoS23_up_s, rhoS23_dn_l, rhoS23_up_l = common.load_observation(obsdir, 'Global/Santini23_JWST_SMD.dat', [0, 1, 2, 3, 4, 5, 6, 7])
rhoS23_cosmocorr = rhoS23 * hobs/h0
Expand All @@ -708,9 +717,9 @@ def plot_stellar_mass_cosmic_density(plt, outdir, obsdir, redshifts, h0, mstarde
rhoS23_errdn2 = np.log10(rhoS23) - np.log10(rhoS23_dn_l)
rhoS23_errup2 = np.log10(rhoS23_up_l) - np.log10(rhoS23)

ax.errorbar(zmidS23, np.log10(rhoS23_cosmocorr), xerr=[zmidS23 - zdnS23, zupS23 - zmidS23], yerr=[rhoS23_errdn2, rhoS23_errup2], ls='None', ecolor = 'grey', mec='grey', marker='s', label='Santini+23')
ax.errorbar(zmidS23, np.log10(rhoS23_cosmocorr), xerr=[zmidS23 - zdnS23, zupS23 - zmidS23], yerr=[rhoS23_errdn2, rhoS23_errup2], ls='None', ecolor = 'darkgreen', mec='darkgreen', marker='s', label='Santini+23')

common.prepare_legend(ax, ['k','b','r','k','grey', 'navy', 'grey', 'grey'], loc=3)
common.prepare_legend(ax, ['k','b','r','k','grey', 'darkorange', 'navy', 'darkgreen'], loc=3)

common.savefig(outdir, fig, "cosmic_smd_compL18.pdf")

Expand Down
50 changes: 43 additions & 7 deletions standard_plots/smf.py
Original file line number Diff line number Diff line change
Expand Up @@ -432,7 +432,7 @@ def plot_stellarmf_z_molcomp(plt, outdir, obsdir, h0, plotz, hist_smf):

common.savefig(outdir, fig, 'stellarmf_z_resolutioncomparison.pdf')

def plot_HImf_z0(plt, outdir, obsdir, h0, plotz_HImf, hist_HImf, hist_HImf_cen, hist_HImf_sat):
def plot_HImf_z0(plt, outdir, obsdir, h0, plotz_HImf, hist_HImf, hist_HImf_cen, hist_HImf_sat, hist_HImf2):

fig = plt.figure(figsize=(5,4.5))
xtit = "$\\rm log_{10} (\\rm M_{\\rm HI}/M_{\odot})$"
Expand Down Expand Up @@ -470,6 +470,9 @@ def plot_HImf_z0(plt, outdir, obsdir, h0, plotz_HImf, hist_HImf, hist_HImf_cen,

# Predicted HIMF
if plotz_HImf[0]:
#y = hist_HImf2[0,:]
#ind = np.where(y < 0.)
#ax.plot(xmf[ind],y[ind],'r', linestyle='solid', label ='Shark v2.0')
y = hist_HImf[0,:]
ind = np.where(y < 0.)
ax.plot(xmf[ind],y[ind],'r', linestyle='solid', label ='Shark v2.0')
Expand Down Expand Up @@ -1350,9 +1353,27 @@ def plot_katsianis_data(ax, i):
plt.tight_layout()
common.savefig(outdir, fig, 'SSFR_functions_z0.pdf')

def calculate_HImass_above_threshold(MHI, rgas, thresh):

rlow = np.log10(0.01*rgas) #from 0.1*rgas
rupp = np.log10(10*rgas) #up to 10 times r50 gas
dr = 0.2
rbins = np.arange(rlow,rupp,dr)
xr = rbins + dr/2.0
re = rgas/1.67 #for an exponential disk

sigma_HI = MHI/(2.0 * np.pi * re**2.0) * np.exp(-10**xr/re)
ind = np.where(sigma_HI > thresh)
if(len(sigma_HI[ind] > 0)):
rin = 10**xr[ind]
MHIout = MHI * (1 - ( 1 + max(rin)/re) * np.exp(-max(rin)/re))
else:
MHIout = 0
return MHIout


def prepare_data(hdf5_data, index, hist_smf, hist_smf_offset, hist_smf_cen, hist_smf_sat,
hist_smf_30kpc, hist_HImf, hist_HImf_cen, hist_HImf_sat, hist_H2mf,
hist_smf_30kpc, hist_HImf, hist_HImf2, hist_HImf_cen, hist_HImf_sat, hist_H2mf,
hist_H2mf_cen, hist_H2mf_sat, mainseq, mainseqsf, sfe, mainseq_cen,
mainseqsf_cen, sfe_cen, mainseq_sat, mainseqsf_sat, sfe_sat, mzr,
fmzr, mzr_cen, mzr_sat, plotz, plotz_HImf, passive_fractions, hist_ssfr,
Expand All @@ -1362,8 +1383,16 @@ def prepare_data(hdf5_data, index, hist_smf, hist_smf_offset, hist_smf_cen, hist
(h0, volh, sfr_disk, sfr_burst, mdisk, mbulge, rstar_disk, mBH, mHI, mH2,
mgas_disk, mHI_bulge, mH2_bulge, mgas_bulge, mgas_metals_disk, mgas_metals_bulge,
mstars_metals_disk, mstars_metals_bulge, typeg, mvir_hosthalo, rstar_bulge,
mbulge_mergers, mbulge_diskins, mbulge_mergers_assembly, mbulge_diskins_assembly) = hdf5_data
mbulge_mergers, mbulge_diskins, mbulge_mergers_assembly, mbulge_diskins_assembly,
sAM_atomic_disk, vmax, rgas_disk) = hdf5_data

MHI_abovethresh = np.zeros(shape = len(mdisk))
thresh = 1 * 1e6 #Msun/kpc^2
for j in range(0, len(mdisk)):
if(mHI[j] > 0):
MHI_abovethresh[j] = calculate_HImass_above_threshold(mHI[j]/h0, rgas_disk[j]/h0*1e3, thresh)

print(mHI/h0, MHI_abovethresh)
zstar = (mstars_metals_disk + mstars_metals_bulge) / (mdisk + mbulge)
rcomb = (rstar_disk * mdisk + rstar_bulge * mbulge) / (mdisk + mbulge) / h0 * 1e3
mgas = mgas_disk+mgas_bulge
Expand Down Expand Up @@ -1450,9 +1479,12 @@ def prepare_data(hdf5_data, index, hist_smf, hist_smf_offset, hist_smf_cen, hist
#gas mass functions and scaling relations
ind = np.where((mHI+mHI_bulge) > 0)
mass_atom[ind] = np.log10(mHI[ind]+mHI_bulge[ind]) - np.log10(float(h0)) + np.log10(XH)

H_HI, _ = np.histogram(mass_atom,bins=np.append(mbins,mupp))
hist_HImf[index,:] = hist_HImf[index,:] + H_HI
ind = np.where(MHI_abovethresh+mHI_bulge > 0)
mass_atom[ind] = np.log10(MHI_abovethresh[ind]+mHI_bulge[ind]) - np.log10(float(h0)) + np.log10(XH)
H_HI, _ = np.histogram(mass_atom,bins=np.append(mbins,mupp))
hist_HImf2[index,:] = hist_HImf2[index,:] + H_HI

ind = np.where(((mHI+mHI_bulge) > 0) & (typeg == 0))
mass_atom_cen[ind] = np.log10(mHI[ind]+mHI_bulge[ind]) - np.log10(float(h0)) + np.log10(XH)
Expand Down Expand Up @@ -1554,6 +1586,7 @@ def prepare_data(hdf5_data, index, hist_smf, hist_smf_offset, hist_smf_cen, hist

hist_HImf_err[index,:] = (hist_HImf[index,:] - np.sqrt(hist_HImf[index,:]))/vol/dm
hist_HImf[index,:] = hist_HImf[index,:]/vol/dm
hist_HImf2[index,:] = hist_HImf2[index,:]/vol/dm
hist_HImf_cen[index,:] = hist_HImf_cen[index,:]/vol/dm
hist_HImf_sat[index,:] = hist_HImf_sat[index,:]/vol/dm

Expand Down Expand Up @@ -1617,6 +1650,7 @@ def main(modeldir, outdir, redshift_table, subvols, obsdir):

plotz = np.empty(shape=(len(zlist)), dtype=np.bool_)
hist_HImf = np.zeros(shape = (len(zlist), len(mbins)))
hist_HImf2 = np.zeros(shape = (len(zlist), len(mbins)))
hist_HImf_cen = np.zeros(shape = (len(zlist), len(mbins)))
hist_HImf_sat = np.zeros(shape = (len(zlist), len(mbins)))
plotz_HImf = np.empty(shape=(len(zlist)), dtype=np.bool_)
Expand All @@ -1633,12 +1667,13 @@ def main(modeldir, outdir, redshift_table, subvols, obsdir):
'mgas_metals_disk', 'mgas_metals_bulge',
'mstars_metals_disk', 'mstars_metals_bulge', 'type',
'mvir_hosthalo', 'rstar_bulge', 'mstars_burst_mergers',
'mstars_burst_diskinstabilities', 'mstars_bulge_mergers_assembly', 'mstars_bulge_diskins_assembly')}
'mstars_burst_diskinstabilities', 'mstars_bulge_mergers_assembly', 'mstars_bulge_diskins_assembly',
'specific_angular_momentum_disk_gas_atom', 'vmax_subhalo', 'rgas_disk')}

for index, snapshot in enumerate(redshift_table[zlist]):
hdf5_data = common.read_data(modeldir, snapshot, fields, subvols)
mass = prepare_data(hdf5_data, index, hist_smf, hist_smf_offset, hist_smf_cen,
hist_smf_sat, hist_smf_30kpc, hist_HImf, hist_HImf_cen, hist_HImf_sat,
hist_smf_sat, hist_smf_30kpc, hist_HImf, hist_HImf2, hist_HImf_cen, hist_HImf_sat,
hist_H2mf, hist_H2mf_cen, hist_H2mf_sat, mainseq, mainseqsf,
sfe, mainseq_cen, mainseqsf_cen, sfe_cen, mainseq_sat,
mainseqsf_sat, sfe_sat, mzr, fmzr, mzr_cen, mzr_sat, plotz,
Expand Down Expand Up @@ -1682,6 +1717,7 @@ def take_log(array):
take_log(hist_smf_sat)
take_log(hist_smf_offset)
take_log(hist_HImf)
take_log(hist_HImf2)
take_log(hist_H2mf)
take_log(hist_HImf_cen)
take_log(hist_H2mf_cen)
Expand All @@ -1692,7 +1728,7 @@ def take_log(array):
plot_stellarmf_z(plt, outdir, obsdir, h0, plotz, hist_smf, hist_smf_cen, hist_smf_sat, hist_smf_offset, hist_smf_30kpc)
plot_stellarmf_galcomponents(plt, outdir, obsdir, h0, plotz, hist_smf, hist_smf_comp)
plot_stellarmf_z_molcomp(plt, outdir, obsdir, h0, plotz, hist_smf)
plot_HImf_z0(plt, outdir, obsdir, h0, plotz_HImf, hist_HImf, hist_HImf_cen, hist_HImf_sat)
plot_HImf_z0(plt, outdir, obsdir, h0, plotz_HImf, hist_HImf, hist_HImf_cen, hist_HImf_sat, hist_HImf2)
plot_H2mf_z0(plt, outdir, obsdir, h0, plotz_HImf, hist_H2mf, hist_H2mf_cen, hist_H2mf_sat)
plot_SSFR_Mstars(plt, outdir, mainseq, mainseq_cen, mainseq_sat)
plot_mzr(plt, outdir, obsdir, h0, mzr, mzr_cen, mzr_sat)
Expand Down
17 changes: 13 additions & 4 deletions standard_plots/smf_all_z.py
Original file line number Diff line number Diff line change
Expand Up @@ -127,6 +127,7 @@ def add_weaver22_data(zappend, file_read='0.2z0.5', label = True):

# Driver al. (2022, z=0). Chabrier IMF
z0obs = []
z0obsPSO = []
lm, p, dp = common.load_observation(obsdir, 'mf/SMF/GAMAIV_Driver22.dat', [0,1,2])
hobs = 0.7
xobs = lm + 2.0 * np.log10(hobs/h0)
Expand All @@ -144,7 +145,8 @@ def add_weaver22_data(zappend, file_read='0.2z0.5', label = True):
lm, p, dpdn, dpup = common.load_observation(obsdir, 'mf/SMF/SMF_Li2009.dat', [0,1,2,3])
xobs = lm - 2.0 * np.log10(hobs) + 2.0 * np.log10(hobs/h0)
yobs = p + 3.0 * np.log10(hobs) - 3.0 * np.log10(hobs/h0)
z0obs.append((observation("Li&White+2009", xobs, yobs, abs(dpdn), dpup, err_absolute=False), 'd'))
#z0obs.append((observation("Li&White+2009", xobs, yobs, abs(dpdn), dpup, err_absolute=False), 'd'))
z0obsPSO.append((observation("Li+2009 (PSO)", xobs, yobs, abs(dpdn), dpup, err_absolute=False), 'd'))


# Moustakas (Chabrier IMF), ['Moustakas+2013, several redshifts']
Expand Down Expand Up @@ -226,7 +228,7 @@ def add_weaver22_data(zappend, file_read='0.2z0.5', label = True):
z7obs = []
add_weaver22_data(z7obs,file_read='6.5z7.5')

return (z0obs, z05obs, z1obs, z2obs, z3obs, z4obs, z5obs, z6obs, z7obs)
return (z0obs, z05obs, z1obs, z2obs, z3obs, z4obs, z5obs, z6obs, z7obs, z0obsPSO)

def plot_SMHM_z(plt, outdir, zlist, halo_mass_rel):

Expand Down Expand Up @@ -343,7 +345,7 @@ def plot_berhoozi13(ax, z, labels, label):

def plot_stellarmf_z(plt, outdir, obsdir, h0, hist_smf, hist_smf_err, hist_smf_30kpc, hist_smf_pass, hist_smf_pass_err):

(z0obs, z05obs, z1obs, z2obs, z3obs, z4obs, z5obs, z6obs, z7obs) = load_smf_observations(obsdir, h0)
(z0obs, z05obs, z1obs, z2obs, z3obs, z4obs, z5obs, z6obs, z7obs, z0obsPSO) = load_smf_observations(obsdir, h0)

PlotLagos18 = True
def plot_lagos18_smf(ax, z):
Expand Down Expand Up @@ -391,6 +393,10 @@ def plot_lagos18_smf(ax, z):
common.prepare_ax(ax, xmin, xmax, ymin, ymax, xtit, ytitle, locators=(0.1, 1, 0.1))
ax.text(xleg, yleg, 'z=%s' % str(z))

if( idx == 0):
for obs, marker in z0obsPSO:
common.errorbars(ax, obs.x, obs.y, obs.yerrdn, obs.yerrup, 'darkgreen',
marker, err_absolute=obs.err_absolute, label=obs.label, markersize=6)
# Observations
for obs, marker in obs_and_markers:
common.errorbars(ax, obs.x, obs.y, obs.yerrdn, obs.yerrup, 'grey',
Expand Down Expand Up @@ -423,7 +429,10 @@ def plot_lagos18_smf(ax, z):
# colors += ['r']
elif idx > 1:
colors = ['r']
colors += ['grey', 'grey','grey']
if idx == 0:
colors += ['darkgreen', 'grey', 'grey','grey']
else:
colors += ['grey', 'grey','grey']

if idx ==0 or idx == 3:
common.prepare_legend(ax, colors)
Expand Down
Loading

0 comments on commit 94a8a3b

Please sign in to comment.