<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Phoenix BT-Settl Bolometric Corrections
Figuring out the best method of handling Phoenix bolometric correction files.
Step1:  Change to directory containing bolometric correction files.
Step2:  Load a bolometric correction table, say for the Cousins AB photometric system.
Step3:  Now, the structure of the file is quite irregular. The grid is not rectangular, which is not an immediate problem. The table is strucutred such that column 0 contains Teff in increasing order, followed by logg in column 1 in increasing order. However, metallicities in column 2 appear to be in decreasing order, which may be a problem for simple interpolation routines. Alpha abundances follow and are in increasing order, but since this is a "standard" grid, whereby alpha enrichment is a function of metallicity, we can ignore it for the moment.
Let's take a first swing at the problem by using the LinearND Interpolator from SciPy.
Step4:  The surface compiled, but that is not a guarantee that the interpolation will work successfully. Some tests are required to confirm this is the case. Let's try a few Teffs at logg = 5 with solar metallicity.
Step5:  This agrees with data in the bolometric correciton table.
Teff      logg   [Fe/H]  [a/Fe]         B           V           R           I
1500.00    5.00    0.00    0.00       -15.557     -16.084     -11.560      -9.291
Now, let's raise the temperature.
Step6:  Again, we have a good match to tabulated values,
Teff      logg   [Fe/H]  [a/Fe]         B           V           R           I
3000.00    5.00    0.00    0.00        -6.603      -5.641      -4.566      -3.273
However, since we are using a tabulated metallicity, the interpolation may proceed without too much trouble. If we select a metallicity between grid points, how do we fare?
Step7:  This appears consistent. What about progressing to lower metallicity values?
Step8:  For reference, at [Fe/H] = $-0.5$ dex, we have
Teff      logg   [Fe/H]  [a/Fe]         B           V           R           I
3000.00    5.00   -0.50    0.20        -6.533      -5.496      -4.424      -3.154
The interpolation routine has seemingly handled the non-monotonic nature of the metallicity column, as all interpolate values lie between values at the two respective nodes.
Now let's import an isochrone and calcuate colors for stellar models for comparison against MARCS bolometric corrections.
Step9:  Make sure there are magnitudes and colors associated with this isochrone.
Step10:  A standard isochrone would only have 6 columns, so 11 indicates this isochrone does have photometric magnitudes computed, likely BV(Ic) (JK)2MASS.
Step11:  For each Teff and logg combination we now have BCs for BV(RI)c from BT-Settl models. Now we need to convert the bolometric corrections to absolute magnitudes.
Step12:  Let's try something different
Step13:  Create an interpolation surface from the magnitude table.
Step14:  Compute magnitudes for a Dartmouth isochrone.
Step15:  Convert surface magnitudes to absolute magnitudes using the distance modulus and the radius of the star.
Step16:  Now compare against MARCS values.
Step17:  Load an isochrone from the Lyon-Phoenix series.
Step18:  Export a new isochrone with colors from AGSS09 (PHX)
Step19:  Separate Test Case
These are clearly not correct and are between 1 and 2 magnitudes off from expected values. Need to reproduce the Phoenix group's results, first.
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import scipy.interpolate as scint
Explanation: Phoenix BT-Settl Bolometric Corrections
Figuring out the best method of handling Phoenix bolometric correction files.
End of explanation
cd /Users/grefe950/Projects/starspot/starspot/color/tab/phx/
Explanation: Change to directory containing bolometric correction files.
End of explanation
bc_table = np.genfromtxt('colmag.BT-Settl.server.JOHNSON.AB.bolcor', comments='!')
Explanation: Load a bolometric correction table, say for the Cousins AB photometric system.
End of explanation
test_surface = scint.LinearNDInterpolator(bc_table[:, :3], bc_table[:, 4:])
Explanation: Now, the structure of the file is quite irregular. The grid is not rectangular, which is not an immediate problem. The table is strucutred such that column 0 contains Teff in increasing order, followed by logg in column 1 in increasing order. However, metallicities in column 2 appear to be in decreasing order, which may be a problem for simple interpolation routines. Alpha abundances follow and are in increasing order, but since this is a "standard" grid, whereby alpha enrichment is a function of metallicity, we can ignore it for the moment.
Let's take a first swing at the problem by using the LinearND Interpolator from SciPy.
End of explanation
test_surface(np.array([1500., 5.0, 0.0]))
Explanation: The surface compiled, but that is not a guarantee that the interpolation will work successfully. Some tests are required to confirm this is the case. Let's try a few Teffs at logg = 5 with solar metallicity.
End of explanation
test_surface(np.array([3000., 5.0, 0.0]))
Explanation: This agrees with data in the bolometric correciton table.
Teff      logg   [Fe/H]  [a/Fe]         B           V           R           I
1500.00    5.00    0.00    0.00       -15.557     -16.084     -11.560      -9.291
Now, let's raise the temperature.
End of explanation
test_surface(np.array([3000., 5.0, 0.1]))
Explanation: Again, we have a good match to tabulated values,
Teff      logg   [Fe/H]  [a/Fe]         B           V           R           I
3000.00    5.00    0.00    0.00        -6.603      -5.641      -4.566      -3.273
However, since we are using a tabulated metallicity, the interpolation may proceed without too much trouble. If we select a metallicity between grid points, how do we fare?
End of explanation
test_surface(np.array([3000., 5.0, -0.2]))
Explanation: This appears consistent. What about progressing to lower metallicity values?
End of explanation
iso = np.genfromtxt('/Users/grefe950/evolve/dmestar/iso/dmestar_00120.0myr_z+0.00_a+0.00_marcs.iso')
Explanation: For reference, at [Fe/H] = $-0.5$ dex, we have
Teff      logg   [Fe/H]  [a/Fe]         B           V           R           I
3000.00    5.00   -0.50    0.20        -6.533      -5.496      -4.424      -3.154
The interpolation routine has seemingly handled the non-monotonic nature of the metallicity column, as all interpolate values lie between values at the two respective nodes.
Now let's import an isochrone and calcuate colors for stellar models for comparison against MARCS bolometric corrections.
End of explanation
iso.shape
Explanation: Make sure there are magnitudes and colors associated with this isochrone.
End of explanation
test_bcs = test_surface(10**iso[:,1], iso[:, 2], 0.0)
test_bcs.shape
Explanation: A standard isochrone would only have 6 columns, so 11 indicates this isochrone does have photometric magnitudes computed, likely BV(Ic) (JK)2MASS.
End of explanation
bol_mags = 4.74 - 2.5*iso[:, 3]
for i in range(test_bcs.shape[1]):
    bcs = -1.0*np.log10(10**iso[:, 1]/5777.) + test_bcs[:, i] - 5.0*iso[:, 4]
    if i == 0:
        test_mags = bol_mags - bcs
    else:
        test_mags = np.column_stack((test_mags, bol_mags - bcs))
iso[50, 0:4], iso[50, 6:], test_mags[50]
Explanation: For each Teff and logg combination we now have BCs for BV(RI)c from BT-Settl models. Now we need to convert the bolometric corrections to absolute magnitudes.
End of explanation
col_table = np.genfromtxt('colmag.BT-Settl.server.COUSINS.AB', comments='!')
Explanation: Let's try something different: using the color tables provided by the Phoenix group, from which the bolometric corrections are calculated.
End of explanation
col_surface = scint.LinearNDInterpolator(col_table[:, :3], col_table[:, 4:8])
Explanation: Create an interpolation surface from the magnitude table.
End of explanation
phx_mags = col_surface(10.0**iso[:, 1], iso[:, 2], 0.0)
Explanation: Compute magnitudes for a Dartmouth isochrone.
End of explanation
for i in range(phx_mags.shape[1]):
    phx_mags[:, i] = phx_mags[:, i] - 5.0*np.log10(10**iso[:, 4]*6.956e10/3.086e18) + 5.0
Explanation: Convert surface magnitudes to absolute magnitudes using the distance modulus and the radius of the star.
End of explanation
iso[40, :5], iso[40, 6:], phx_mags[40]
Explanation: Now compare against MARCS values.
End of explanation
phx_iso = np.genfromtxt('/Users/grefe950/Notebook/Projects/ngc2516_spots/data/phx_isochrone_120myr.txt')
fig, ax = plt.subplots(1, 2, figsize=(12., 8.), sharey=True)
ax[0].set_xlim(0.0, 2.0)
ax[1].set_xlim(0.0, 4.0)
ax[0].set_ylim(16, 2)
ax[0].plot(iso[:, 6] - iso[:, 7], iso[:, 7], lw=3, c="#b22222")
ax[0].plot(phx_mags[:, 0] - phx_mags[:, 1], phx_mags[:, 1], lw=3, c="#1e90ff")
ax[0].plot(phx_iso[:, 7] - phx_iso[:, 8], phx_iso[:, 8], dashes=(20., 5.), lw=3, c="#555555")
ax[1].plot(iso[:, 7] - iso[:, 8], iso[:, 7], lw=3, c="#b22222")
ax[1].plot(phx_mags[:, 1] - phx_mags[:, 3], phx_mags[:, 1], lw=3, c="#1e90ff")
ax[1].plot(phx_iso[:, 8] - phx_iso[:, 10], phx_iso[:, 8], dashes=(20., 5.), lw=3, c="#555555")
Explanation: Load an isochrone from the Lyon-Phoenix series.
End of explanation
new_isochrone = np.column_stack((iso[:, :6], phx_mags))
np.savetxt('/Users/grefe950/Notebook/Projects/pleiades_colors/data/dmestar_00120.0myr_z+0.00_a+0.00_mixed.iso', 
            new_isochrone, fmt='%16.8f')
Explanation: Export a new isochrone with colors from AGSS09 (PHX)
End of explanation
tmp = -10.*np.log10(3681./5777.) + test_surface(3681., 4.78, 0.0) #+ 5.0*np.log10(0.477)
tmp 
4.74 - 2.5*(-1.44) - tmp
Explanation: Separate Test Case
These are clearly not correct and are between 1 and 2 magnitudes off from expected values. Need to reproduce the Phoenix group's results, first.
End of explanation
<END_TASK>