<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Apsidal Motion Age for HD 144548
Here, I am attempting to derive an age for the triple eclipsing hierarchical triple HD 144548 (Upper Scoripus member) based on the observed orbital precession (apsidal motion) of the inner binary system's orbit about the tertiary companion (star A). A value for the orbital precession is reported by Alonso et al. (2015, arXiv
Step1:  What complicates the issue is that the interior structure constants for the B components also vary as a function of age, so we need to compute a mean mass track using the $c_2$ coefficients and the individual $k_2$ values.
Step2:  Create tracks with equally spaced time steps.
Step3:  Now, compute the $c_2$ coefficients for each age.
Step4:  With that, we have an estimate for the mean B component properties as a function of age. One complicating factor is the "radius" of the average B component. If we are modeling the potential created by the Ba/Bb components as that of a single star, we need to assume that the A component never enters into any region of the combined potential that is dominated by either component.
Unfortunately, it is very likely that the ratio of the Ba/Bb binary "radius" to the semi-major axis of the A/B orbit is going to be a dominant contributor to the apsidal motion.
Attempt 1
Step5:  We are now in a position to compute the classical apsidal motion rate from the combined A/B tracks.
Step6:  One can see from this that the general relativistic component is a very small contribution to the total apsidal motion of the system. 
Let's look at the evolution of the apsidal motion for the A/B binary system.
Step7:  How sensitive is this to the properties of the A component, which are fairly uncertain?
<ASSISTANT_TASK:>
Python Code:
def c2(masses, radii, e, a, rotation=None):
    f = (1.0 - e**2)**-2
    g = (8.0 + 12.0*e**2 + e**4)*f**(5.0/2.0) / 8.0
    if rotation == None:
        omega_ratio_sq = 0.0
    elif rotation == 'synchronized':
        omega_ratio_sq = (1.0 + e)/(1.0 - e)**3 
    else:
        omega_ratio_sq = 0.0
    
    c2_0 = (omega_ratio_sq*(1.0 + masses[1]/masses[0])*f + 15.0*g*masses[1]/masses[0])*(radii[0]/a)**5
    c2_1 = (omega_ratio_sq*(1.0 + masses[0]/masses[1])*f + 15.0*g*masses[0]/masses[1])*(radii[1]/a)**5
    
    return c2_0, c2_1
# parameters for the orbit of Ba/Bb
e = 0.0015
a = 7.249776
masses = [0.984, 0.944]
# c2_B = c2(masses, radii, e, a)
Explanation: Apsidal Motion Age for HD 144548
Here, I am attempting to derive an age for the triple eclipsing hierarchical triple HD 144548 (Upper Scoripus member) based on the observed orbital precession (apsidal motion) of the inner binary system's orbit about the tertiary companion (star A). A value for the orbital precession is reported by Alonso et al. (2015, arXiv: 1510.03773) as $\dot{\omega} = 0.0235 \pm 0.002\ {\rm deg\, cycle}^{-1}$, obtained from photo-dynamical modeling of a Kepler/K2 lightcurve. The technique of determining an age from apsidal motion observed in young binary systems is detailed by Feiden & Dotter (2013, ApJ, 765, 86). Their technique relies heavily on the analytical framework for the classical theory of orbital precession due to tidal and rotational distortions of graviational potentials by Kopal (1978, ASSL, 68) and the inclusion of general relativistic orbital precession by Giménez (1985, ApJ, 297, 405).  
In brief, the technique outlined by Feiden & Dotter (2013) relies on the fact that young stars are contracting quasi-hydrostatically as they approach the main-sequence. As they contract, the mean density of the increases (assuming the star has constant mass), thereby altering the distribution of mass with respect to the mean density. This alters the interior structure parameter, which is related to the deviation from sphericity of the star and its resulting gravitational potential. A non-symmetric potential induces a precession of the point of periastron in a companion star's orbit, provided the orbit is eccentric. 
Since the internal density structure of a young star is changing as it contracts, the inferred interior structure parameter, and thus the induced pertrubation on the star's gravitational potential, also changes. Therefore, the rate at which the precision of a binary companions point of periastron occurs changes as a function of time. By measuring the rate of precession, one can then estimate the age of the system by inferring the required density distribution to induce that rate of precession, subject to the constraint that the orbital and stellar fundamental properties must be well determined - hence the reason why Feiden & Dotter (2013) focused exclusively on eclipsing binary systems.
While a rate of orbital precession was measured by Alonso et al. (2015) for HD 144548, and the properties of all three stars were determined with reasonable precision, there is a fundamental difficulty: it's a triple system. The method outlined by Feiden & Dotter (2013) was intended for binary systems, with no discussion of the influence of a tertiary companion.
Fortunately, the measured orbtial precision is for the orbit of the inner binary (Ba/Bb) about the tertiary star (A). Below, I focus on modeling the inner binary as a single object orbiting the tertiary star with a mass equal to the sum of the component masses (thus more massive than component A). 
The first big hurdle is to figure out how to treat the Ba/Bb component as a single star. For an initial attempt, we can assume that the B component is a "single" star with a mass equal to the total mass of the binary system with an interior structure constant equal to the weighted mean of the two individual interior structure constants.
To compute the mean interior structure constants, we first need to compute the individual weights $c_{2, i}$. For $e = 0$, we have $f(e) = g(e) = 1$.
End of explanation
import numpy as np
trk_Ba = np.genfromtxt('/Users/grefe950/evolve/dmestar/trk/gs98/p000/a0/amlt1884/m0980_GS98_p000_p0_y28_mlt1.884.trk')
trk_Bb = np.genfromtxt('/Users/grefe950/evolve/dmestar/trk/gs98/p000/a0/amlt1884/m0940_GS98_p000_p0_y28_mlt1.884.trk')
Explanation: What complicates the issue is that the interior structure constants for the B components also vary as a function of age, so we need to compute a mean mass track using the $c_2$ coefficients and the individual $k_2$ values.
End of explanation
from scipy.interpolate import interp1d
log10_age = np.arange(6.0, 8.0, 0.01)  # log10(age/yr)
ages = 10**log10_age
icurve = interp1d(trk_Ba[:,0], trk_Ba, kind='linear', axis=0)
new_trk_Ba = icurve(ages)
icurve = interp1d(trk_Bb[:,0], trk_Bb, kind='linear', axis=0)
new_trk_Bb = icurve(ages)
Explanation: Create tracks with equally spaced time steps.
End of explanation
mean_trk_B = np.empty((len(ages), 3))
for i, age in enumerate(ages):
    c2s = c2(masses, [10**new_trk_Ba[i, 4], 10**new_trk_Bb[i, 4]], e, a, 
             rotation='synchronized')
    avg_k2 = (c2s[0]*new_trk_Ba[i, 10] + c2s[1]*new_trk_Bb[i, 10])/(sum(c2s))
    
    mean_trk_B[i] = np.array([age, 10**new_trk_Ba[i, 4] + 10**new_trk_Bb[i, 4], avg_k2])
Explanation: Now, compute the $c_2$ coefficients for each age.
End of explanation
e2 = 0.2652
a2 = 66.2319
masses_2 = [1.44, 1.928]
trk_A = np.genfromtxt('/Users/grefe950/evolve/dmestar/trk/gs98/p000/a0/amlt1884/m1450_GS98_p000_p0_y28_mlt1.884.trk', 
                      usecols=(0,1,2,3,4,5,6,7,8,9,10))
icurve = interp1d(trk_A[:,0], trk_A, kind='linear', axis=0)
new_trk_A = icurve(ages)
Explanation: With that, we have an estimate for the mean B component properties as a function of age. One complicating factor is the "radius" of the average B component. If we are modeling the potential created by the Ba/Bb components as that of a single star, we need to assume that the A component never enters into any region of the combined potential that is dominated by either component.
Unfortunately, it is very likely that the ratio of the Ba/Bb binary "radius" to the semi-major axis of the A/B orbit is going to be a dominant contributor to the apsidal motion.
Attempt 1: Semi-major axis + radius of B component
Let's define orbtial properties of the (A, B) system.
End of explanation
cl_apsidal_motion_rate = np.empty((len(ages), 2))
for i, age in enumerate(ages):
    c2_AB = c2(masses_2, [10**new_trk_A[i, 4], a + 0.5*mean_trk_B[i, 1]], e2, a2)
    cl_apsidal_motion_rate[i] = np.array([age, 360.0*(c2_AB[0]*new_trk_A[i, 10] + c2_AB[1]*mean_trk_B[i, 2])])
GR_apsidal_motion_rate = 5.45e-4*(sum(masses)/33.945)**(2./3.) / (1.0 - e2**2)  # Giménez (1985)
GR_apsidal_motion_rate
Explanation: We are now in a position to compute the classical apsidal motion rate from the combined A/B tracks.
End of explanation
%matplotlib inline
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 1, figsize=(8., 8.), sharex=True)
ax.grid(True)
ax.tick_params(which='major', axis='both', length=15., labelsize=18.)
ax.set_xlabel('Age (yr)', fontsize=20., family='serif')
ax.set_ylabel('Apsidal Motion Rate (deg / cycle)', fontsize=20., family='serif')
ax.plot([1.0e6, 1.0e8], [0.0215, 0.0215], '--', lw=1, c='#555555')
ax.plot([1.0e6, 1.0e8], [0.0255, 0.0255], '--', lw=1, c='#555555')
ax.semilogx(cl_apsidal_motion_rate[:, 0], cl_apsidal_motion_rate[:, 1], '-', lw=2, c='#b22222')
Explanation: One can see from this that the general relativistic component is a very small contribution to the total apsidal motion of the system. 
Let's look at the evolution of the apsidal motion for the A/B binary system.
End of explanation
icurve = interp1d(cl_apsidal_motion_rate[:,1], cl_apsidal_motion_rate[:,0], kind='linear')
print icurve(0.0235)/1.0e6, icurve(0.0255)/1.0e6, icurve(0.0215)/1.0e6
Explanation: How sensitive is this to the properties of the A component, which are fairly uncertain?
End of explanation
<END_TASK>