In [1]:
try:
reload(simple_blurring_churning)
reload(numpy)
except (ImportError, NameError):
import simple_blurring_churning
import numpy
from galpy.util import bovy_plot
from matplotlib import pyplot
from scipy import special
figsize(8,6)
bovy_plot.bovy_print(fig_width=7.,fig_height=5.,axes_labelsize=20,
text_fontsize=20,legend_fontsize=12,
xtick_labelsize=14,ytick_labelsize=14)
We obtain the distribution of guiding-center radii at a given R due to blurring using a Dehnen DF. The default parameters are a disk with a flat rotation curve, a scale length of 3 kpc, and a velocity dispersion that is constant with radius at 31.4 km/s. These parameters can be changed by specifying the hr=, sr=, and hs= keywords. For example, at R=8 kpc we obtain the distribution of guiding center radii and plot it:
In [2]:
Rgs= numpy.linspace(0.,18.,101)
In [3]:
pRgsb4= simple_blurring_churning.blurring_pRgR(Rgs,4.)
pRgsb8= simple_blurring_churning.blurring_pRgR(Rgs,8.)
pRgsb14= simple_blurring_churning.blurring_pRgR(Rgs,14.)
In [4]:
plot(Rgs,pRgsb4/numpy.sum(pRgsb4)/(Rgs[1]-Rgs[0]),color='k',lw=2.)
plot(Rgs,pRgsb8/numpy.sum(pRgsb8)/(Rgs[1]-Rgs[0]),color='c',lw=2.)
plot(Rgs,pRgsb14/numpy.sum(pRgsb14)/(Rgs[1]-Rgs[0]),color='y',lw=2.)
xlabel(r'$R_g$')
ylabel(r'$p(R_g|R)$')
Out[4]:
We construct the distribution of guiding-center radii including churning, by assuming that churning induces a Gaussian spreading of the initial guiding-center radii $R_{g,i}$:
$p(R_{g_f}|R_{g,i},\tau) = \mathcal{N}\left(R_{g,f}|R_{g,i},0.01+0.2\,\tau\,R_{g,i}\,e^{-(R_{g,i}-8\,\mathrm{kpc})^2/16}\right)\,,$
where $\tau$ is age and $\mathcal{N}(\cdot|m,V)$ is a Gaussian with mean $m$ and variance $V$. This spread increases as the square-root of time and in this model induces the largest spread around 8 kpc (the code has a parameter fmig=1 that multiplies the $0.2\tau\ldots$ term; setting fmig=0 turns migration off in the calculations). The distribution of final guiding-center radii as a function of initial guiding-center radii at three different final radii after 10 Gyr is shown next:
In [5]:
pRgfc4= simple_blurring_churning.churning_pRgfRgi(Rgs,4.,10.)
pRgfc8= simple_blurring_churning.churning_pRgfRgi(Rgs,8.,10.)
pRgfc14= simple_blurring_churning.churning_pRgfRgi(Rgs,14.,10.)
In [6]:
plot(Rgs,pRgfc4/numpy.sum(pRgfc4)/(Rgs[1]-Rgs[0]),color='k',lw=2.)
plot(Rgs,pRgfc8/numpy.sum(pRgfc8)/(Rgs[1]-Rgs[0]),color='c',lw=2.)
plot(Rgs,pRgfc14/numpy.sum(pRgfc14)/(Rgs[1]-Rgs[0]),color='y',lw=2.)
xlabel(r'$R_{g,f}$')
ylabel(r'$p(R_{g,f}|R_{g,i})$')
Out[6]:
This leads to the following distribution of initial guiding-center radii at different radii today after 1 Gyr and after 10 Gyr (we also display the blurring effect as the dashed line)
$p(R_{g,i}|R,\tau) = \int \mathrm{d}R_{g,f} p_{\mathrm{blurring}}(R_{g,f}|R)\,p(R_{g,i}|R_{g,f},\tau)\,,$
where $p_{\mathrm{bluring}}$ is the blurring distribution from above. We can write this as (because we assume that migration does not change the surface density of the disk $p(R_{g,i}|R_{g,f},\tau) = p(R_{g,f}|R_{g,i},\tau)$):
$p(R_{g,i}|R,\tau) = \int \mathrm{d}R_{g,f} p_{\mathrm{blurring}}(R_{g,i}|R)\,p(R_{g,f}|R_{g,i},\tau)\,,$
In [7]:
pRgRtau4_1= simple_blurring_churning.churning_pRgRtau(Rgs,4.,1.)
pRgRtau4_10= simple_blurring_churning.churning_pRgRtau(Rgs,4.,10.)
pRgRtau8_1= simple_blurring_churning.churning_pRgRtau(Rgs,8.,1.)
pRgRtau8_10= simple_blurring_churning.churning_pRgRtau(Rgs,8.,10.)
pRgRtau14_1= simple_blurring_churning.churning_pRgRtau(Rgs,14.,1.)
pRgRtau14_10= simple_blurring_churning.churning_pRgRtau(Rgs,14.,10.)
In [8]:
line1= plot(Rgs,pRgRtau4_1/numpy.sum(pRgRtau4_1)/(Rgs[1]-Rgs[0]),color='k',lw=2.)
line2= plot(Rgs,pRgRtau4_10/numpy.sum(pRgRtau4_10)/(Rgs[1]-Rgs[0]),ls='--',color='k',lw=2.)
plot(Rgs,pRgRtau8_1/numpy.sum(pRgRtau8_1)/(Rgs[1]-Rgs[0]),color='c',lw=2.)
plot(Rgs,pRgRtau8_10/numpy.sum(pRgRtau8_10)/(Rgs[1]-Rgs[0]),ls='--',color='c',lw=2.)
plot(Rgs,pRgRtau14_1/numpy.sum(pRgRtau14_1)/(Rgs[1]-Rgs[0]),color='y',lw=2.)
plot(Rgs,pRgRtau14_10/numpy.sum(pRgRtau14_10)/(Rgs[1]-Rgs[0]),ls='--',color='y',lw=2.)
line3= plot(Rgs,pRgsb4/numpy.sum(pRgsb4)/(Rgs[1]-Rgs[0]),ls=':',color='k',lw=2.)
plot(Rgs,pRgsb8/numpy.sum(pRgsb8)/(Rgs[1]-Rgs[0]),ls=':',color='c',lw=2.)
plot(Rgs,pRgsb14/numpy.sum(pRgsb14)/(Rgs[1]-Rgs[0]),ls=':',color='y',lw=2.)
xlabel(r'$R_g$')
ylabel(r'$p(R_g|R,\tau)$')
l1= pyplot.legend((line3[0],line1[0],line2[0]),
(r'$\mathrm{blurring}$',
r'$\mathrm{churning,\ age=1\, Gyr}$',
r'$\mathrm{churning,\ age=10\, Gyr}$'),
loc='upper right',#bbox_to_anchor=(.91,.375),
numpoints=8,
prop={'size':16},
frameon=False)
We see that at time increases, stars come from further away. We can also demonstrate the two-dimensional behavior of this function:
In [9]:
ages= numpy.linspace(1.,10.,101)
pRgage4= numpy.array([simple_blurring_churning.churning_pRgRtau(Rgs,4.,a) for a in ages])
pRgage8= numpy.array([simple_blurring_churning.churning_pRgRtau(Rgs,8.,a) for a in ages])
pRgage14= numpy.array([simple_blurring_churning.churning_pRgRtau(Rgs,14.,a) for a in ages])
In [10]:
figsize(16,6)
subplot(1,3,1)
imshow(pRgage4,origin='lower',cmap='jet',extent=[Rgs[0],Rgs[-1],ages[0],ages[-1]],aspect=(Rgs[-1]-Rgs[0])/(ages[-1]-ages[0]))
bovy_plot.bovy_text(r'$R=4\, \mathrm{kpc}$',top_right=True,size=18.,color='w')
xlabel(r'$R_g$')
ylabel(r'$\tau$')
subplot(1,3,2)
imshow(pRgage8,origin='lower',cmap='jet',extent=[Rgs[0],Rgs[-1],ages[0],ages[-1]],aspect=(Rgs[-1]-Rgs[0])/(ages[-1]-ages[0]))
bovy_plot.bovy_text(r'$R=8\, \mathrm{kpc}$',top_right=True,size=18.,color='w')
xlabel(r'$R_g$')
subplot(1,3,3)
imshow(pRgage14,origin='lower',cmap='jet',extent=[Rgs[0],Rgs[-1],ages[0],ages[-1]],aspect=(Rgs[-1]-Rgs[0])/(ages[-1]-ages[0]))
bovy_plot.bovy_text(r'$R=14\, \mathrm{kpc}$',top_left=True,size=18.,color='w')
xlabel(r'$R_g$')
tight_layout()
figsize(8,6)
Integrating this over time, we get the churning+blurring distribution of guiding-center radii at a given radius:
In [11]:
pRgRc4= simple_blurring_churning.churning_pRgR(Rgs,4.)
pRgRc8= simple_blurring_churning.churning_pRgR(Rgs,8.)
pRgRc14= simple_blurring_churning.churning_pRgR(Rgs,14.)
In [12]:
line1= plot(Rgs,pRgRc4/numpy.sum(pRgRc4)/(Rgs[1]-Rgs[0]),color='k',lw=2.)
plot(Rgs,pRgRc8/numpy.sum(pRgRc8)/(Rgs[1]-Rgs[0]),color='c',lw=2.)
plot(Rgs,pRgRc14/numpy.sum(pRgRc14)/(Rgs[1]-Rgs[0]),color='y',lw=2.)
line2= plot(Rgs,pRgsb4/numpy.sum(pRgsb4)/(Rgs[1]-Rgs[0]),ls='--',color='k',lw=2.)
plot(Rgs,pRgsb8/numpy.sum(pRgsb8)/(Rgs[1]-Rgs[0]),ls='--',color='c',lw=2.)
plot(Rgs,pRgsb14/numpy.sum(pRgsb14)/(Rgs[1]-Rgs[0]),ls='--',color='y',lw=2.)
xlabel(r'$R_g$')
ylabel(r'$p(R_g|R)$')
l1= pyplot.legend((line2[0],line1[0]),
(r'$\mathrm{blurring}$',
r'$\mathrm{churning,\ marginalized\ over\ age}$'),
loc='upper right',#bbox_to_anchor=(.91,.375),
numpoints=8,
prop={'size':16},
frameon=False)
To get the metallicity distribution function (MDF), we start from an assumed intrinsic MDF modeled as a skew-normal distribution with fixed shape but with a mean changing assuming a metallicity gradient. We then mix these initial MDFs using the effect from blurring and churning as modeled above.
We first plot the initial MDF:
In [13]:
fehs= numpy.linspace(-1.3,0.7,102)
In [14]:
pFehInit4= simple_blurring_churning.pFehRg(fehs,4.)
pFehInit8= simple_blurring_churning.pFehRg(fehs,8.)
pFehInit14= simple_blurring_churning.pFehRg(fehs,14.)
In [15]:
plot(fehs,pFehInit4/numpy.sum(pFehInit4)/(fehs[1]-fehs[0]),color='k',lw=2.)
plot(fehs,pFehInit8/numpy.sum(pFehInit8)/(fehs[1]-fehs[0]),color='c',lw=2.)
plot(fehs,pFehInit14/numpy.sum(pFehInit14)/(fehs[1]-fehs[0]),color='y',lw=2.)
xlabel(r'$[\mathrm{Fe/H}]$')
ylabel(r'$p([\mathrm{Fe/H}]|R)$')
Out[15]:
To apply the churning prescriptions above, we need the relation between age and metallicity at a given initial radius. We assume that the metallicity increases logarithmically with age over 10 Gyr, starting at -1.3:
In [16]:
age4= simple_blurring_churning.ageFehRg(fehs,4.)
age8= simple_blurring_churning.ageFehRg(fehs,8.)
age14= simple_blurring_churning.ageFehRg(fehs,14.)
In [17]:
plot(fehs,age4,color='k',lw=2.)
plot(fehs,age8,color='c',lw=2.)
plot(fehs,age14,color='y',lw=2.)
xlabel(r'$[\mathrm{Fe/H}]$')
ylabel(r'$\mathrm{age}([\mathrm{Fe/H}],R)$')
ylim(0.,10.)
xlim(-1.5,0.7)
Out[17]:
We obtain the MDF under the influence of blurring as a special case of the blurring+churning calculation below. The blurring MDF is given by
$p([\mathrm{Fe/H}]|R) = \int \mathrm{d}R_g\,p([\mathrm{Fe/H}]|R_g)\,p(R_g|R)\,.$
where $p(R_g|R)$ is the blurring distribution of guiding-center radii displayed in the first figure in this document.
In [18]:
pFehb4= simple_blurring_churning.blurring_pFehR(fehs,4.)
pFehb6= simple_blurring_churning.blurring_pFehR(fehs,6.)
pFehb8= simple_blurring_churning.blurring_pFehR(fehs,8.)
pFehb10= simple_blurring_churning.blurring_pFehR(fehs,10.)
pFehb12= simple_blurring_churning.blurring_pFehR(fehs,12.)
pFehb14= simple_blurring_churning.blurring_pFehR(fehs,14.)
#some more init
pFehInit6= simple_blurring_churning.pFehRg(fehs,6.)
pFehInit10= simple_blurring_churning.pFehRg(fehs,10.)
pFehInit12= simple_blurring_churning.pFehRg(fehs,12.)
In [19]:
saveToFile= False
bovy_plot.bovy_print(fig_width=7.,fig_height=5.,axes_labelsize=20,
text_fontsize=20,legend_fontsize=12,
xtick_labelsize=14,ytick_labelsize=14)
line1= plot(fehs,pFehb4/numpy.sum(pFehb4)/(fehs[1]-fehs[0]),color='k',zorder=1,lw=2.)
line7= plot(fehs,pFehInit4/numpy.sum(pFehInit4)/(fehs[1]-fehs[0]),ls='--',color='k',zorder=0)
line2= plot(fehs,pFehb6/numpy.sum(pFehb6)/(fehs[1]-fehs[0]),color='b',zorder=1,lw=2.)
plot(fehs,pFehInit6/numpy.sum(pFehInit6)/(fehs[1]-fehs[0]),ls='--',color='b',zorder=0)
line3= plot(fehs,pFehb8/numpy.sum(pFehb8)/(fehs[1]-fehs[0]),color='c',zorder=1,lw=2.)
plot(fehs,pFehInit8/numpy.sum(pFehInit8)/(fehs[1]-fehs[0]),ls='--',color='c',zorder=0)
line4= plot(fehs,pFehb10/numpy.sum(pFehb10)/(fehs[1]-fehs[0]),color='g',zorder=1,lw=2.)
plot(fehs,pFehInit10/numpy.sum(pFehInit10)/(fehs[1]-fehs[0]),ls='--',color='g',zorder=0)
line5= plot(fehs,pFehb12/numpy.sum(pFehb12)/(fehs[1]-fehs[0]),color='r',zorder=1,lw=2.)
plot(fehs,pFehInit12/numpy.sum(pFehInit12)/(fehs[1]-fehs[0]),ls='--',color='r',zorder=0)
line6= plot(fehs,pFehb14/numpy.sum(pFehb14)/(fehs[1]-fehs[0]),color='y',zorder=1,lw=2.)
plot(fehs,pFehInit14/numpy.sum(pFehInit14)/(fehs[1]-fehs[0]),ls='--',color='y',zorder=0)
xlim(-1.8,0.7)
ylim(0.,6.0)
xlabel(r'$[\mathrm{Fe/H}]$')
ylabel(r'$p([\mathrm{Fe/H}]|R)$')
bovy_plot._add_ticks()
#Add legend
l1= pyplot.legend((line1[0],line2[0],line3[0],line4[0],line5[0],line6[0]),
(r'$3 < R < 5$'+'\n'+r'$ \gamma_1 = %.2f$' % simple_blurring_churning.skewness(fehs,pFehb4),
r'$5 < R < 7$'+'\n'+r'$ \gamma_1 = %.2f$' % simple_blurring_churning.skewness(fehs,pFehb6),
r'$7 < R < 9$'+'\n'+r'$ \gamma_1 = %.2f$' % simple_blurring_churning.skewness(fehs,pFehb8),
r'$9 < R < 11$'+'\n'+r'$ \gamma_1 = %.2f$' % simple_blurring_churning.skewness(fehs,pFehb10),
r'$11 < R < 13$'+'\n'+r'$ \gamma_1 = %.2f$' % simple_blurring_churning.skewness(fehs,pFehb12),
r'$13 < R < 15$'+'\n'+r'$ \gamma_1 = %.2f$' % simple_blurring_churning.skewness(fehs,pFehb14)),
loc='upper left',#bbox_to_anchor=(.91,.375),
numpoints=8,
prop={'size':16},
frameon=False)
l2= pyplot.legend((line7[0],line1[0]),
(r'$\mathrm{circular\ orbits}$',
r'$\mathrm{blurring}$'),
loc='upper right',#bbox_to_anchor=(.91,.375),
numpoints=8,
prop={'size':16},
frameon=False)
pyplot.gca().add_artist(l1)
pyplot.gca().add_artist(l2)
if saveToFile:
bovy_plot.bovy_end_print('mdf_blurring_r.ps')
To calculate the churning MDF, we need to integrate over the current age distribution at the given radius. We can calculate the current and initial age distribution (without churning) as
$p(\tau|R) = \int \mathrm{d}R_g\,p(\tilde{[\mathrm{Fe/H}]}|R_g)\,p(R_g|R)/\left|\mathrm{d}\tau/\mathrm{d}[\mathrm{Fe/H}]\right|\,.$
where $\tilde{[\mathrm{Fe/H}]}$ is the metallicity corresponding to $R_g$ and $\tau$ (see the relations between age, metallicity, and $R_g$ above. $p(R_g|R)$ is the distribution of guiding-center radii due to churning, which we approximate by integrating $p(R_g|R,\tau)$ over $\tau$ with a constant $p(\tau)$ (this isn't quite correct).
In [20]:
ages= numpy.linspace(1.,10.,101)
In [21]:
pAgec4= simple_blurring_churning.churning_pAgeR(ages,4.)
pAgei4= simple_blurring_churning.pAgeRg(ages,4.)
pAgec8= simple_blurring_churning.churning_pAgeR(ages,8.)
pAgei8= simple_blurring_churning.pAgeRg(ages,8.)
pAgec14= simple_blurring_churning.churning_pAgeR(ages,14.)
pAgei14= simple_blurring_churning.pAgeRg(ages,14.)
In [22]:
plot(ages,pAgec4/numpy.sum(pAgec4)/(ages[1]-ages[0]),color='k',lw=2.)
plot(ages,pAgei4/numpy.sum(pAgei4)/(ages[1]-ages[0]),ls='--',color='k',lw=2.)
plot(ages,pAgec8/numpy.sum(pAgec8)/(ages[1]-ages[0]),color='c',lw=2.)
plot(ages,pAgei8/numpy.sum(pAgei8)/(ages[1]-ages[0]),ls='--',color='c',lw=2.)
plot(ages,pAgec14/numpy.sum(pAgec14)/(ages[1]-ages[0]),color='y',lw=2.)
plot(ages,pAgei14/numpy.sum(pAgei14)/(ages[1]-ages[0]),ls='--',color='y',lw=2.)
xlabel(r'$\mathrm{Age}$')
ylabel(r'$p(Age|R)$')
Out[22]:
As the final and initial distributions are quite similar, we will just use the initial distribution in the evaluation of the churning DF. We now obtain the MDF due to churning and blurring by doing the integral
$p([\mathrm{Fe/H}]|R) = \iint \mathrm{d}\tau \mathrm{d} R_g p([\mathrm{Fe/H}]|\tau,R_g)\,p(R_g|R,\tau)\,p(\tau|R)\,,$
because the first factor is a delta function in the simple chemical-evolution model, this can be simplified to
$p([\mathrm{Fe/H}]|R) = \int \mathrm{d} R_g p(R_g|R,\tilde{\tau})\,p(\tilde{\tau}|R)/\left|\mathrm{d}[\mathrm{Fe/H}]/\mathrm{d}\tau\right|\,,$
where $\tilde{\tau}$ is the age corresponding to a given $R_g$ and metallicity in the simple-chemical evolution model (see above). Evaluating this expression at different $R$ gives
In [23]:
tfehs= numpy.linspace(-1.3,0.7,101)
pFehc4= simple_blurring_churning.churning_pFehR(tfehs,4.)
pFehc6= simple_blurring_churning.churning_pFehR(tfehs,6.)
pFehc8= simple_blurring_churning.churning_pFehR(tfehs,8.)
pFehc10= simple_blurring_churning.churning_pFehR(tfehs,10.)
pFehc12= simple_blurring_churning.churning_pFehR(tfehs,12.)
pFehc14= simple_blurring_churning.churning_pFehR(tfehs,14.)
In [24]:
saveToFile= False
tfehs= numpy.linspace(-1.3,0.7,101)
normInit= 5./3.
bovy_plot.bovy_print(fig_width=7.,fig_height=5.,axes_labelsize=20,
text_fontsize=20,legend_fontsize=12,
xtick_labelsize=14,ytick_labelsize=14)
line1= plot(tfehs,pFehc4/numpy.nansum(pFehc4)/(tfehs[1]-tfehs[0]),color='k',zorder=1,lw=2.)
line7= plot(fehs,pFehInit4/numpy.sum(pFehInit4)/(fehs[1]-fehs[0])/normInit,ls='--',color='k',zorder=0)
line2= plot(tfehs,pFehc6/numpy.nansum(pFehc6)/(tfehs[1]-tfehs[0]),color='b',zorder=1,lw=2.)
plot(fehs,pFehInit6/numpy.sum(pFehInit6)/(fehs[1]-fehs[0])/normInit,ls='--',color='b',zorder=0)
line3= plot(tfehs,pFehc8/numpy.nansum(pFehc8)/(tfehs[1]-tfehs[0]),color='c',zorder=1,lw=2.)
plot(fehs,pFehInit8/numpy.sum(pFehInit8)/(fehs[1]-fehs[0])/normInit,ls='--',color='c',zorder=0)
line4= plot(tfehs,pFehc10/numpy.nansum(pFehc10)/(tfehs[1]-tfehs[0]),color='g',zorder=1,lw=2.)
plot(fehs,pFehInit10/numpy.sum(pFehInit10)/(fehs[1]-fehs[0])/normInit,ls='--',color='g',zorder=0)
line5= plot(tfehs,pFehc12/numpy.nansum(pFehc12)/(tfehs[1]-tfehs[0]),color='r',zorder=1,lw=2.)
plot(fehs,pFehInit12/numpy.sum(pFehInit12)/(fehs[1]-fehs[0])/normInit,ls='--',color='r',zorder=0)
line6= plot(tfehs,pFehc14/numpy.nansum(pFehc14)/(tfehs[1]-tfehs[0]),color='y',zorder=1,lw=2.)
plot(fehs,pFehInit14/numpy.sum(pFehInit14)/(fehs[1]-fehs[0])/normInit,ls='--',color='y',zorder=0)
xlim(-1.8,0.7)
ylim(0.,3.5)
xlabel(r'$[\mathrm{Fe/H}]$')
ylabel(r'$p([\mathrm{Fe/H}]|R)$')
bovy_plot._add_ticks()
#Add legend
l1= pyplot.legend((line1[0],line2[0],line3[0],line4[0],line5[0],line6[0]),
(r'$3 < R < 5$'+'\n'+r'$ \gamma_1 = %.2f$' % simple_blurring_churning.skewness(tfehs,pFehc4),
r'$5 < R < 7$'+'\n'+r'$ \gamma_1 = %.2f$' % simple_blurring_churning.skewness(tfehs,pFehc6),
r'$7 < R < 9$'+'\n'+r'$ \gamma_1 = %.2f$' % simple_blurring_churning.skewness(tfehs,pFehc8),
r'$9 < R < 11$'+'\n'+r'$ \gamma_1 = %.2f$' % simple_blurring_churning.skewness(tfehs,pFehc10),
r'$11 < R < 13$'+'\n'+r'$ \gamma_1 = %.2f$' % simple_blurring_churning.skewness(tfehs,pFehc12),
r'$13 < R < 15$'+'\n'+r'$ \gamma_1 = %.2f$' % simple_blurring_churning.skewness(tfehs,pFehc14)),
loc='upper left',#bbox_to_anchor=(.91,.375),
numpoints=8,
prop={'size':16},
frameon=False)
l2= pyplot.legend((line7[0],line1[0]),
(r'$\mathrm{circular\ orbits}$',
r'$\mathrm{churning}$'),
loc='upper right',#bbox_to_anchor=(.91,.375),
numpoints=8,
prop={'size':16},
frameon=False)
pyplot.gca().add_artist(l1)
pyplot.gca().add_artist(l2)
if saveToFile:
bovy_plot.bovy_end_print('mdf_churning_r.ps')
In [25]:
print simple_blurring_churning.skewness(tfehs,pFehc4), simple_blurring_churning.skewness(fehs,pFehInit4)
print simple_blurring_churning.skewness(tfehs,pFehc6), simple_blurring_churning.skewness(fehs,pFehInit6)
print simple_blurring_churning.skewness(tfehs,pFehc8), simple_blurring_churning.skewness(fehs,pFehInit8)
print simple_blurring_churning.skewness(tfehs,pFehc10), simple_blurring_churning.skewness(fehs,pFehInit10)
print simple_blurring_churning.skewness(tfehs,pFehc12), simple_blurring_churning.skewness(fehs,pFehInit12)
print simple_blurring_churning.skewness(tfehs,pFehc14), simple_blurring_churning.skewness(fehs,pFehInit14)
We see that the inner disk MDFs remain negatively skewed like the initial MDF, but the outer disk MDFs become positively skewed.
We can use the same formalism to compute the age-metallicity relation at different radii:
$p(\tau,[\mathrm{Fe/H}]|R) = \int \mathrm{d}R_g\,p([\mathrm{Fe/H}]|R_g,\tau)\,p(R_g|R,\tau)\,p(\tau|R)\,,$
which can be simplified to (because the first factor is a delta function)
$p(\tau,[\mathrm{Fe/H}]|R) = p(\tilde{R_g}|R,\tau)\,p(\tau|R)/\left|\mathrm{d}[\mathrm{Fe/H}]/\mathrm{d}R_g\right|\,,$
where $\tilde{R_g}$ is the guiding-center radius corresponding to $\tau$ and $[\mathrm{Fe/H}]$. This gives
In [26]:
tages= numpy.linspace(1.,9.9,101)
tfehs= numpy.linspace(-1.,0.7,101)
pAgeFehc4= numpy.array([simple_blurring_churning.churning_pFehAgeR(tfehs,a,4.) for a in tages])
pAgeFehc8= numpy.array([simple_blurring_churning.churning_pFehAgeR(tfehs,a,8.) for a in tages])
pAgeFehc14= numpy.array([simple_blurring_churning.churning_pFehAgeR(tfehs,a,14.) for a in tages])
In [27]:
figsize(16,6)
subplot(1,3,1)
imshow(pAgeFehc4.T,origin='lower',cmap='jet',extent=[tages[0],tages[-1],tfehs[0],tfehs[-1]],aspect=(tages[-1]-tages[0])/(tfehs[-1]-tfehs[0]))
plot(tages,numpy.nansum(pAgeFehc4*tfehs,axis=1)/numpy.nansum(pAgeFehc4,axis=1))
indx= (age4 > 1.)*(age4 < 10.)
plot(age4[indx],fehs[indx],color='b',ls='--')
xlabel(r'$\mathrm{Age}$')
ylabel(r'$[\mathrm{Fe/H}]$')
subplot(1,3,2)
imshow(pAgeFehc8.T,origin='lower',cmap='jet',extent=[tages[0],tages[-1],tfehs[0],tfehs[-1]],aspect=(tages[-1]-tages[0])/(tfehs[-1]-tfehs[0]))
plot(tages,numpy.nansum(pAgeFehc8*tfehs,axis=1)/numpy.nansum(pAgeFehc8,axis=1))
indx= (age8 > 1.)*(age8 < 10.)
plot(age8[indx],fehs[indx],color='b',ls='--')
xlabel(r'$\mathrm{Age}$')
subplot(1,3,3)
imshow(pAgeFehc14.T,origin='lower',cmap='jet',extent=[tages[0],tages[-1],tfehs[0],tfehs[-1]],aspect=(tages[-1]-tages[0])/(tfehs[-1]-tfehs[0]))
plot(tages,numpy.nansum(pAgeFehc14*tfehs,axis=1)/numpy.nansum(pAgeFehc14,axis=1))
indx= (age14 > 1.)*(age14 < 10.)
plot(age14[indx],fehs[indx],color='b',ls='--')
xlabel(r'$\mathrm{Age}$')
tight_layout()
figsize(8,6)
We can also just plot the conditional distribution:
In [28]:
saveToFile= False
tfehs= numpy.linspace(-1.,0.7,101)
figsize(16,6)
if saveToFile:
bovy_plot.bovy_print(fig_width=5.,fig_height=5.,axes_labelsize=18,
text_fontsize=18,legend_fontsize=12,
xtick_labelsize=14,ytick_labelsize=14)
figsize(11,4)
subplot(1,3,1)
imshow((pAgeFehc4/numpy.tile(numpy.nansum(pAgeFehc4,axis=1),(len(tfehs),1)).T).T,origin='lower',cmap='gist_yarg',extent=[tages[0],tages[-1],tfehs[0],tfehs[-1]],aspect=(tages[-1]-tages[0])/(tfehs[-1]-tfehs[0]))
bovy_plot.bovy_dens2d((pAgeFehc4/numpy.tile(numpy.nansum(pAgeFehc4,axis=1),(len(tfehs),1)).T).T,origin='lower',
justcontours=True,overplot=True,cntrmass=True,levels=special.erf(numpy.arange(1,4)/numpy.sqrt(2.)),
xrange=[tages[0],tages[-1]],yrange=[tfehs[0],tfehs[-1]],cntrcolors='0.35')
line1= plot(tages,numpy.nansum(pAgeFehc4*tfehs,axis=1)/numpy.nansum(pAgeFehc4,axis=1),color='r',lw=2.)
line2= plot(tages,simple_blurring_churning.fehAgeRg(tages,4.),color='b',ls='-',lw=2.)
l1= pyplot.legend((line1[0],line2[0]),
(r'$\mathrm{stellar\ [Fe/H](Age)}$',
r'$\mathrm{local\ ISM\ [Fe/H](Age)}$'),
loc='lower left',
numpoints=8,
prop={'size':16-3*saveToFile},
frameon=False)
for legobj in l1.legendHandles:
legobj.set_linewidth(2.0)
xlabel(r'$\mathrm{Age\,(Gyr)}$')
ylabel(r'$[\mathrm{Fe/H}]$')
pyplot.annotate(r'$R = 4\,\mathrm{kpc}$',(0.5,1.07+0.02*saveToFile),xycoords='axes fraction',horizontalalignment='center',
verticalalignment='top',size=20.-3*saveToFile)
bovy_plot._add_ticks()
subplot(1,3,2)
imshow((pAgeFehc8/numpy.tile(numpy.nansum(pAgeFehc8,axis=1),(len(tfehs),1)).T).T,origin='lower',cmap='gist_yarg',extent=[tages[0],tages[-1],tfehs[0],tfehs[-1]],aspect=(tages[-1]-tages[0])/(tfehs[-1]-tfehs[0]))
bovy_plot.bovy_dens2d((pAgeFehc8/numpy.tile(numpy.nansum(pAgeFehc8,axis=1),(len(tfehs),1)).T).T,origin='lower',
justcontours=True,overplot=True,cntrmass=True,levels=special.erf(numpy.arange(1,4)/numpy.sqrt(2.)),
xrange=[tages[0],tages[-1]],yrange=[tfehs[0],tfehs[-1]],cntrcolors='0.35')
plot(tages,numpy.nansum(pAgeFehc8*tfehs,axis=1)/numpy.nansum(pAgeFehc8,axis=1),color='r',lw=2.)
plot(tages,simple_blurring_churning.fehAgeRg(tages,8.),color='b',ls='-',lw=2.)
xlabel(r'$\mathrm{Age\,(Gyr)}$')
pyplot.annotate(r'$R = 8\,\mathrm{kpc}$',(0.5,1.07+0.02*saveToFile),xycoords='axes fraction',horizontalalignment='center',
verticalalignment='top',size=20.-3*saveToFile)
bovy_plot._add_ticks()
subplot(1,3,3)
imshow((pAgeFehc14/numpy.tile(numpy.nansum(pAgeFehc14,axis=1),(len(tfehs),1)).T).T,origin='lower',cmap='gist_yarg',extent=[tages[0],tages[-1],tfehs[0],tfehs[-1]],aspect=(tages[-1]-tages[0])/(tfehs[-1]-tfehs[0]))
bovy_plot.bovy_dens2d((pAgeFehc14/numpy.tile(numpy.nansum(pAgeFehc14,axis=1),(len(tfehs),1)).T).T,origin='lower',
justcontours=True,overplot=True,cntrmass=True,levels=special.erf(numpy.arange(1,4)/numpy.sqrt(2.)),
xrange=[tages[0],tages[-1]],yrange=[tfehs[0],tfehs[-1]],cntrcolors='0.35')
plot(tages,numpy.nansum(pAgeFehc14*tfehs,axis=1)/numpy.nansum(pAgeFehc14,axis=1),color='r',lw=2.)
plot(tages,simple_blurring_churning.fehAgeRg(tages,14.),color='b',ls='-',lw=2.)
xlabel(r'$\mathrm{Age\,(Gyr)}$')
pyplot.annotate(r'$R = 14\,\mathrm{kpc}$',(0.5,1.07+0.02*saveToFile),xycoords='axes fraction',horizontalalignment='center',
verticalalignment='top',size=20.-3*saveToFile)
bovy_plot._add_ticks()
if saveToFile:
bovy_plot.bovy_end_print('agemdf_churning_r.ps',orientation='landscape')
else:
tight_layout()
figsize(8,6)
In [28]:
In [28]: