Skip to content

Commit

Permalink
Merge branch 'master' of https://github.com/philrosenfield/tpagb-cali…
Browse files Browse the repository at this point in the history 
…bration
  • Loading branch information
@philrosenfield
philrosenfield committed Feb 13, 2017
2 parents a145ce2 + dcd14d4 commit 14ff3e9
Show file tree
Hide file tree
Showing 2 changed files with 76 additions and 68 deletions.
131 changes: 69 additions & 62 deletions consistency/sim_as_phot.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,86 +24,93 @@ def cull_data(simgalaxname):
mass = sgal.data['m_ini'][good]
tot_mass = np.sum(mass)
tp_mass = sgal.data['m_ini'][inds]
import pdb; pdb.set_trace()
return sgal, tot_mass, tp_mass, mass, good


def csfr_masshist():
def csfr_masshist(csfr=True):
# use padua sfh from Dan Weisz and PARSEC sfh from me (third attempt)
pd_loc = os.path.join(snap_src, 'varysfh/extpagb/sim_as_phot/padua/')
pc_loc = os.path.join(snap_src, 'varysfh/extpagb/sim_as_phot/parsec/')
# trilegal simulations from input sfh
sims = ['out_eso540-030_f606w_f814w.mcmc.zc_caf09_v1.2s_m36_s12d_ns_nas_bestsfr.dat',
'out_ugc5139_f555w_f814w.mcmc.zc_caf09_v1.2s_m36_s12d_ns_nas_bestsfr.dat']
sim = 'out_ugc5139_f555w_f814w.mcmc.zc_caf09_v1.2s_m36_s12d_ns_nas_bestsfr.dat'
# zc merged sfh solutions
sfhs= ['eso540-030_f606w_f814w.mcmc.zc.dat',
'ugc5139_f555w_f814w.mcmc.zc.dat']
sfh = 'ugc5139_f555w_f814w.mcmc.zc.dat'
# calc sfh output for Av, dmod, etc (not really needed for csfr plots)
metas = ['eso540-030_f606w_f814w.sfh',
'ugc5139_f555w_f814w.sfh']
meta = 'ugc5139_f555w_f814w.sfh'

dm = 0.1 # could be input... mass step
bins = np.arange(0.8, 20 + dm, dm)
for i, sim in enumerate(sims):
if i == 0:
continue
# make tp-agb and full mass histograms
sgalpd, tot_masspd, tp_masspd, masspd, goodpd = \
cull_data(os.path.join(pd_loc, 'caf09_v1.2s_m36_s12d_ns_nas', sim))
hpd = np.histogram(tp_masspd, bins=bins)[0] / tot_masspd
hfpd = np.histogram(masspd, bins=bins)[0] / tot_masspd

sgalpc, tot_masspc, tp_masspc, masspc, goodpc = \
cull_data(os.path.join(pc_loc, 'caf09_v1.2s_m36_s12d_ns_nas', sim))
hpc = np.histogram(tp_masspc, bins=bins)[0] / tot_masspc
hfpc = np.histogram(masspc, bins=bins)[0] / tot_masspc

# compare tp-agb in both sfh solutions
a, p = ks_2samp(hpd, hpc)
target = sgalpd.name.split('_')[1]
print('TPAGB {} KS: {:.4f} {:.4f}'.format(target, a, p))

# load sfh for csfr plots
sfhpd = SFH(os.path.join(pd_loc, sfhs[i]),
meta_file=os.path.join(pd_loc, metas[i]))
sfhpc = SFH(os.path.join(pc_loc, sfhs[i]),
meta_file=os.path.join(pc_loc, metas[i]))

fig, (ax1, ax2) = plt.subplots(nrows=2, figsize=(16, 8))
fig.subplots_adjust(hspace=.4, top=0.97, bottom=0.11)
bins = np.arange(0.55, 20 + dm, dm)

# make tp-agb and full mass histograms
sgalpd, tot_masspd, tp_masspd, masspd, goodpd = \
cull_data(os.path.join(pd_loc, 'caf09_v1.2s_m36_s12d_ns_nas', 'lower_mass', sim))
hpd = np.histogram(tp_masspd, bins=bins)[0] / tot_masspd
hfpd = np.histogram(masspd, bins=bins)[0] / tot_masspd

sgalpc, tot_masspc, tp_masspc, masspc, goodpc = \
cull_data(os.path.join(pc_loc, 'caf09_v1.2s_m36_s12d_ns_nas', 'lower_mass', sim))
hpc = np.histogram(tp_masspc, bins=bins)[0] / tot_masspc
hfpc = np.histogram(masspc, bins=bins)[0] / tot_masspc

# compare tp-agb in both sfh solutions
a, p = ks_2samp(hpd, hpc)
target = sgalpd.name.split('_')[1]
print('TPAGB {} KS: {:.4f} {:.4f}'.format(target, a, p))

# load sfh for csfr plots
sfhpd = SFH(os.path.join(pd_loc, sfh),
meta_file=os.path.join(pd_loc, meta))
sfhpc = SFH(os.path.join(pc_loc, meta),
meta_file=os.path.join(pc_loc, meta))

fig, (ax1, ax2) = plt.subplots(nrows=2, figsize=(16, 8))
fig.subplots_adjust(hspace=.4, top=0.97, bottom=0.11)

if not csfr:
# Padua
ax1 = sfhpd.age_plot(ax=ax1, plt_kw={'color': 'k'})
# Parsec
ax1 = sfhpc.age_plot(ax=ax1, plt_kw={'color': '#30a2da'})
ax1.set_ylim(0, ax1.get_ylim()[1])
ax1.set_xlim(14.1, 0)
else:
# Padua
ax1 = sfhpd.plot_csfr(ax=ax1, plt_kw={'color': 'k'},
fill_between_kw={'color': 'gray', 'alpha': 0.3})
# Parsec
ax1 = sfhpc.plot_csfr(ax=ax1, plt_kw={'color': '#30a2da'},
fill_between_kw={'color': '#30a2da', 'alpha': 0.3})

# tpagb
tpd, = ax2.plot(bins[:-1], hpd, linestyle='steps-mid', label=r'$\rm{Padua}$',
color='k')
tpc, = ax2.plot(bins[:-1], hpc, linestyle='steps-mid', label=r'$\rm{PARSEC}$',
color='#30a2da')
# full
pd, = ax2.plot(bins[:-1], hfpd, linestyle='steps-mid', color='k', alpha=0.3)
pc, = ax2.plot(bins[:-1], hfpc, linestyle='steps-mid', color='#30a2da',
alpha=0.3)

ax2.set_yscale('log')
ax2.set_xlim(0.8, 12.)
ax2.set_ylim(1e-7, 1)
ax2.set_xlabel(r'$\rm{Mass\ (M_\odot)}$')
ax2.set_ylabel(r'$\rm{\#/Mass\ bin\ (%.1f\ M_\odot)}$' % dm)
ax1.set_xlabel(r'$\rm{Time\ (Gyr)}$')
ax1.set_ylabel(r'$\rm{Cumulative\ SF}$')
ax2.legend(((tpd, pd), (tpc, pc)), (r'$\rm{Padua}$', r'$\rm{PARSEC}$'), loc='best')
for ax in [ax1, ax2]:
ax.grid()
ax.tick_params(direction='in', which='both')
lab = r'$\rm{{{}}}$'.format(target.upper().replace('-1', '').replace('-','\!-\!'))
ax1.text(0.05, 0.05, lab, ha='right', fontsize=16, **emboss())
outfile = '{}_csfr_mass{}'.format(target, EXT)
import pdb; pdb.set_trace()
plt.savefig(outfile)
print('wrote {}'.format(outfile))
# tpagb
tpd, = ax2.plot(bins[:-1], hpd, linestyle='steps-mid', label=r'$\rm{Padua}$',
color='k')
tpc, = ax2.plot(bins[:-1], hpc, linestyle='steps-mid', label=r'$\rm{PARSEC}$',
color='#30a2da')
# full
pd, = ax2.plot(bins[:-1], hfpd, linestyle='steps-mid', color='k', alpha=0.3)
pc, = ax2.plot(bins[:-1], hfpc, linestyle='steps-mid', color='#30a2da',
alpha=0.3)

ax2.set_yscale('log')
ax2.set_xlim(0.5, 12.)
ax2.xaxis.set_ticks(np.arange(0.5, 12.5, 0.5))
ax2.set_ylim(1e-7, 1)

ax2.set_xlabel(r'$\rm{Mass\ (M_\odot)}$')
ax2.set_ylabel(r'$\rm{\#/Mass\ bin\ (%.1f\ M_\odot)}$' % dm)
ax1.set_xlabel(r'$\rm{Time\ (Gyr)}$')
#ax1.set_ylabel(r'$\rm{Cumulative\ SF}$')
ax2.legend(((tpd, pd), (tpc, pc)), (r'$\rm{Padua}$', r'$\rm{PARSEC}$'), loc='best')
for ax in [ax1, ax2]:
ax.grid()
ax.tick_params(direction='in', which='both')
lab = r'$\rm{{{}}}$'.format(target.upper().replace('-1', '').replace('-','\!-\!'))
ax1.text(0.98, 0.05, lab, ha='right', fontsize=20, transform=ax1.transAxes, **emboss())
outfile = '{}_csfr_mass{}'.format(target, EXT)

plt.savefig(outfile)
print('wrote {}'.format(outfile))

if __name__ == "__main__":
csfr_masshist()
13 changes: 7 additions & 6 deletions plotting/plot_tracks.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -172,14 +172,14 @@ def vw93_plot(self, agescale=1e5, outfile=None, xlim=None, ylims=None,
Make a plot similar to Vassiliadis and Wood 1993. Instead of Vesc,
I plot C/O.
"""
ylims = ylims or [None] * 6
ycols = ['T_star', 'Tbot', 'L_star', 'CO', 'M_star', 'dMdt']
ylims = ylims or [None] * len(ycols)

if axs is None:
fig, axs = plt.subplots(nrows=6, sharex=True, figsize=(5.4, 10))
fig, axs = plt.subplots(nrows=len(ycols), sharex=True, figsize=(5.4, 10))
fig.subplots_adjust(hspace=0.05, right=0.97, top=0.97, bottom=0.07,
left=0.2)

ycols = ['T_star', 'L_star', 'period', 'CO', 'M_star', 'dMdt']

for i in range(len(axs)):
ycol = ycols[i]
ylim = ylims[i]
Expand Down Expand Up @@ -241,7 +241,8 @@ def str_agescale(scale=1.):
u = '10^%i\ ' % int(np.log10(scale))
return u

tdict = {'T_star': r'$log\ \rm{T}_{\rm{eff}}\ \rm{(K)}$',
tdict = {'Tbot': r'$log\ \rm{T}_{\rm{bce}}\ \rm{(K)}$',
'T_star': r'$log\ \rm{T}_{\rm{eff}}\ \rm{(K)}$',
'L_star': r'$log\ L\ (L_\odot)$',
'period': r'$\rm{P\ (days)}$',
'CO': r'$\rm{C/O}$',
Expand All @@ -262,7 +263,7 @@ def compare_vw93(agbs, outfile=None, xlim=None, ylims=None):
fmt = r'$\rm{Z}=%g$'
annotations = [fmt % agbs[0].Z, fmt % agbs[1].Z]
# sharex is off because I want one column pruned.
fig, axs = plt.subplots(nrows=6, ncols=2, figsize=(8.4, 10))
fig, axs = plt.subplots(nrows=6, ncols=2, figsize=(9, 10))
col1 = axs.T[0]
col2 = axs.T[1]

Expand Down

0 comments on commit 14ff3e9

Please sign in to comment.