<jupyter_start><jupyter_text>Lets see if we can do something about the bug with the Battistini and Bitsch model application to the simulation...Loading in the data first, scraped from the BB paper:<jupyter_code>BBmodel = ascii.read("../sav/Combineddata.csv",data_start=2)
BBmodel.rename_column('\ufeffFe_H', 'Fe_H')

Comp = Table(BBmodel)
Comp<jupyter_output><empty_output><jupyter_text>OK - these are the data from the appropriate figure. We only really care about the $\mathrm{[Fe/H]}$ vs water mass fractionNow we'll do the polynomial fit...<jupyter_code># Now interpolate between points

model_poly = models.Polynomial1D(degree=3)
fitter_poly = fitting.LinearLSQFitter()
best_fit_poly = fitter_poly(model_poly, BBmodel['Fe_H'],BBmodel['H20'])

print(best_fit_poly)

plt.errorbar(BBmodel['Fe_H'],BBmodel['H20'], fmt='k.')
plt.plot(BBmodel['Fe_H'],best_fit_poly(BBmodel['Fe_H']),color='r', linewidth=3)

newxs = np.linspace(-3,1.,100)

plt.plot(newxs, best_fit_poly(newxs))
plt.ylim(-0.1,1.1)<jupyter_output>Model: Polynomial1D
Inputs: ('x',)
Outputs: ('y',)
Model set size: 1
Degree: 3
Parameters:
            c0                  c1                  c2                  c3        
    ------------------ ------------------- ------------------- -------------------
    0.3480026126904762 -0.5659026289285716 -0.3847110642857163 0.12633408333333815<jupyter_text>This doesnt look ideal at the low $\mathrm{[Fe/H]}$ end, so we can try some other options here... what about a spline?<jupyter_code>s = UnivariateSpline(BBmodel['Fe_H'],BBmodel['H20'], k=4)


plt.errorbar(BBmodel['Fe_H'],BBmodel['H20'], fmt='k.')
plt.plot(BBmodel['Fe_H'],s(BBmodel['Fe_H']),color='r', linewidth=3)

newxs = np.linspace(-3,1.,100)

plt.plot(newxs, s(newxs))
plt.ylim(-0.1,1.1)<jupyter_output><empty_output><jupyter_text>Again, not ideal since the mass fraction can get to above 1 (and below 0....). What actually happens in the data we have for stars?<jupyter_code>#we'll grab the APOGEE data using jobovy's apogee code (this requires some set-up, ask Ted for details...)
allstar = apread.allStar(main=True, exclude_star_bad=True, exclude_star_warn=True)<jupyter_output>/Users/tedmackereth/anaconda3/lib/python3.7/site-packages/apogee-1.-py3.7.egg/apogee/tools/read.py:303: RuntimeWarning: Extinction-corrected J,H,K not added because esutil is not installed
  warnings.warn("Extinction-corrected J,H,K not added because esutil is not installed",RuntimeWarning)<jupyter_text>in the BB paper, they looked at C/O abundances in GALAH, to see what happens with O as the $\mathrm{[Fe/H]}$ increases (more O is locked into CO and CO2 at increasing Fe). We can use this idea to figure out what might be the best course of action at the low Fe end as well?<jupyter_code>fig = plt.figure()
fig.set_size_inches(4,4)

#we'll clean up the APOGEE data a bit to only look at stars where the abundances are probably good (red giants)
bad = (allstar['LOGG_ERR'] < 0.1) & (allstar['LOGG'] < 4)&  (allstar['LOGG'] > 1) & (allstar['C_FE'] != -9999.9902) & (allstar['O_FE'] != -9999.9902)



plt.scatter(allstar['FE_H'][bad], allstar['C_FE'][bad]-allstar['O_FE'][bad], s=0.1, lw=0., alpha=0.8, rasterized=True)
plt.xlim(-2.,0.7)
plt.ylim(-1,0.5)

def running_percentile(x, y, bins):
    '''quick function to get the running median/percentiles'''
    bin_inds = np.digitize(x, bins)
    values = np.ones((len(bins),3))*np.nan
    for i in np.unique(bin_inds):
        if i == 0 or i == 15:
            continue
        in_bin = bin_inds == i
        if sum(in_bin) < 10:
            continue
        values[i] = np.percentile(y[in_bin], [16,50,84])
    bin_centers = (bins[1:]+bins[:-1])/2.
    return values, bin_centers
        
bins = np.linspace(-1.5,0.5,15)
medians, bin_centers = running_percentile(allstar['FE_H'][bad], allstar['C_FE'][bad]-allstar['O_FE'][bad], bins)  

plt.plot(bins-((bins[1]-bins[0])/2.), medians[:,1], c='Black')
plt.fill_between(bins-((bins[1]-bins[0])/2.), medians[:,0], medians[:,2], color='Black', alpha=0.3)

plt.axvline(-0.4, color='Black', linestyle='dashed')
plt.axvline(0.4, color='Black', linestyle='dashed')
plt.text(-1.5,0.3, 'BB20 [Fe/H] limits')

plt.xlabel(r'$\mathrm{[Fe/H]}$')
plt.ylabel(r'$\mathrm{[C/O]}$')

plt.savefig('../plots/CO_FEH_APOGEEDR16.pdf')<jupyter_output><empty_output><jupyter_text>the APOGEE behaviour matches GALAH quite well in the BB limits, which is reassuring, but the trend changes quite significantly at low $\mathrm{[Fe/H]}$. Since the C/O drops right off, we might assume that the Oxygen available in forming the ISO's at low metallicities is much higher? Since the accreted dwarfs probably dominate in terms of mass at this metallicity regime, maybe this means that the ISO's accreted in dwarfs would be genuinely disentangle-able from the MW ones...Just to illustrate that there is a lot of accreted debris down there at low metallicity - you can see that this stands out in the Tinsley diagram as a larger scatter in $\mathrm{[Mg/Fe]}$ (alpha elements) at low $\mathrm{[Fe/H]}$.<jupyter_code>plt.scatter(allstar['FE_H'], allstar['MG_FE'], s=0.1, lw=0., color='Black')
plt.xlim(-2,0.7)
plt.ylim(-0.2,0.5)<jupyter_output><empty_output><jupyter_text>In that case, maybe the polynomial is ok. I think the best we can do is to just set the function to be the upper and lower limit:<jupyter_code># Now interpolate between points

def piecewise_poly(x):
    '''this function allows us to set the extrapolation to the limits of the data in x'''
    model_poly = models.Polynomial1D(degree=3)
    fitter_poly = fitting.LinearLSQFitter()
    best_fit_poly = fitter_poly(model_poly, BBmodel['Fe_H'],BBmodel['H20'])
    minx, maxx = np.min(BBmodel['Fe_H']), np.max(BBmodel['Fe_H'])
    minxy = best_fit_poly(minx)
    maxxy = best_fit_poly(maxx)
    if not hasattr(x, '__iter__'):
        if x < minx:
            return minxy
        elif x > maxx:
            return maxxy
        else:
            return best_fit_poly(x)
    else:
        out = np.zeros(len(x))
        out[x < minx] = minxy
        out[x > maxx] = maxxy
        out[(x >= minx) & (x <= maxx)] = best_fit_poly(x[(x >= minx) & (x <= maxx)])
        return out

plt.errorbar(BBmodel['Fe_H'],BBmodel['H20'], fmt='k.')
plt.plot(BBmodel['Fe_H'],piecewise_poly(BBmodel['Fe_H']),color='r', linewidth=3)

newxs = np.linspace(-3,1.,100)

plt.plot(newxs, piecewise_poly(newxs))
files = ['GalaxyA_FOF507.dat', 'EAGLE_MW_L0025N0376_REFERENCE_ApogeeRun_30kpc_working.dat', 'GalaxyC_FOF526.dat']

sims = [ascii.read('../sav/%s' % file) for file in files]
fig, ax = plt.subplots(1,3, sharex=True, sharey=True)
fig.set_size_inches(10,2)
for i in range(len(sims)):
    ax[i].scatter(sims[i]['fe_h'], sims[i]['mg_h']-sims[i]['fe_h'], s=0.01, color='Black')
    ax[i].set_xlabel(r'$\mathrm{[Fe/H]}$')

plt.xlim(-2.5,1.)
plt.ylim(-0.3,0.5)
ax[0].set_ylabel(r'$\mathrm{[Mg/Fe]}$')

fig, ax = plt.subplots(1,3, sharex=True, sharey=True)
fig.set_size_inches(6,2)
for i in range(len(sims)):
    ax[i].scatter(sims[i]['x_g'], sims[i]['z_g'], s=0.1, lw=0., color='Black')
    ax[i].set_xlabel(r'$x\ \mathrm{[kpc]}$')
ax[0].set_ylabel(r'$z\ \mathrm{[kpc]}$')<jupyter_output><empty_output><jupyter_text>lets re-make Chris' plot with the fixed limits (note I also do a thing to 'spread out' the probability density where we extrapolate... not sure if this is useful though...)<jupyter_code>colors = ['#4477AA', '#BB5566', '#DDAA33']
labels = ['late', 'bi-modal', 'early']

frac_water_rich = []
mean_age = []
for i in range(len(sims)):
    water_mass_frac = piecewise_poly(sims[i]['fe_h'])
    bins = np.linspace(np.min(water_mass_frac),np.max(water_mass_frac),20)
    hist, bins = np.histogram(water_mass_frac, bins=bins, density=True)
    old_end = np.copy(bins[-2])
    #figure out the correction to get the final bin probability density right...
    end_correct = (bins[1]-bins[0])/(1-old_end)
    bins[-2] = 1.
    hist[-1] *= end_correct
    plt.step(bins[:-1], hist, color=colors[i], lw=2., label=labels[i])
    frac_water_rich.append(sum(water_mass_frac > 0.4)/len(water_mass_frac))
    mean_age.append(np.mean(sims[i]['age']))
    
plt.xlim(0.,1.)
plt.ylim(0.,3)
plt.legend()
plt.xlabel(r'$\mathrm{H_2O\ mass\ fraction}$')
plt.ylabel(r'$p(\mathrm{H_2O\ mass\ fraction})_i$')

#mark out the regions where we extrapolate
plt.gca().axvspan(bins[-3], 1., alpha=0.1, color='Black')
plt.gca().axvspan(0, bins[0], alpha=0.1, color='Black')
#just for fun, how does it look with the spline?
frac_water_rich = []
mean_age = []
for i in range(len(sims)):
    water_mass_frac = s(sims[i]['fe_h'])
    bins = np.linspace(-0.1,1.1,50)
    plt.hist(water_mass_frac, histtype='step', density=True, bins=bins, lw=2., color=colors[i])
    frac_water_rich.append(sum(water_mass_frac > 0.4)/len(water_mass_frac))
    mean_age.append(np.mean(sims[i]['age']))
    
#mark out the regions where we extrapolate
plt.gca().axvspan(np.max(BBmodel['H20']), 1., alpha=0.1, color='Black')
plt.gca().axvspan(0, np.min(BBmodel['H20']), alpha=0.1, color='Black')

plt.xlim(0.,1.)

plt.xlabel(r'$\mathrm{H_2O\ mass\ fraction}$')
plt.ylabel(r'$p(\mathrm{H_2O\ mass\ fraction})_i$')<jupyter_output><empty_output><jupyter_text>If the behaviour is anything like the spline fit, then (i.e. the water mass fraction goes much higher at low Fe) then we actually still see a peak at ~0.5 for the EAGLE galaxies...<jupyter_code>plt.plot(mean_age, frac_water_rich)
#plot the MDF of each galaxy...
for i in range(len(sims)):
    plt.hist(sims[i]['fe_h'], range=[-2,1.], histtype='step', density=True, bins=30)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1732538686741949, 'fraction_numerical': 0.03398159765788373, 'mean_word_length': 4.45014245014245, 'pattern_counts': {'":': 0, '<': 42, '<?xml version=': 0, '>': 42, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '490863', 'license': 'MIT', 'n_tokens_mistral': 3863, 'n_tokens_neox': 3607, 'n_words': 838, 'path': 'py/iso_abundances.ipynb', 'repo_name': 'jmackereth/ISOgal', 'size': '9564'}