%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

"""
Explanation: Preparation for AGB Group Meeting
Sofia requested that each participant bring a plot (of something!) for discussion at the AGB group meeting on 27 Aug 2015. A few ideas about which plot(s) to bring: (1) demonstration that more recent solar abundance distributions (Asplund et al. 2009; Caffau et al. 2011) are in better agreement with NIR broadband photometric colors and should be implemented in AGB atmosphere models. Lower oxygen abundance is much better for M dwarfs owing to decreased H20, whereas M giant wind driving mechanisms may be inhibited due to decreased oxygen available to capture into grains? (2) Magnetic fields have been tentatively detected at the surface of AGB stars with typical magnetic field strength estimates in the 1 &ndash; 5 G range (Ramstedt, priv. comm.). What magnetic field strengths are required to interact and inhibit convection in AGB star envelopes (neglecting large-scale dynamics)?
End of explanation
"""


Pleiades_NIR = np.genfromtxt('../Projects/pleiades_colors/data/Stauffer_Pleiades_nir.dat', 
                             usecols=(2, 3, 5, 6, 8, 9))  # J, errJ, H, errH, K, errK

dist_mod = 5.62 # distance modulus for the Pleiades

"""
Explanation: Solar Abundances and NIR Broadband Photometric Colors
We'll begin by loading necessary data, here NIR photometric data for the Pleiades (and M67?).
End of explanation
"""


iso_gs98  = np.genfromtxt('../Projects/pleiades_colors/data/dmestar_00120.0myr_z+0.00_a+0.00_phx.iso')
iso_gas07 = np.genfromtxt('../Projects/pleiades_colors/data/dmestar_00120.0myr_z+0.00_a+0.00_marcs.iso')

"""
Explanation: Now we'll need to load appropriate stellar evolution isochrones with bolometric corrections derived from MARCS model atmospheres. MARCS uses Grevesse et al. (2007), but differences between colors derived from GAS07, AGSS09, and CIFIST11 are negligible.
End of explanation
"""


fig, ax = plt.subplots(1, 3, figsize=(15., 8.), sharey=True)

ax[0].set_ylabel('$M_K$', fontsize=22.)
for axis in ax:
    axis.grid(True)
    axis.set_ylim(9., 3.)
    axis.tick_params(which='major', axis='both', length=15., labelsize=16.)

# (J-K)
ax[0].set_xlim(0.2, 1.2)
ax[0].set_xlabel('$(J - K)$', fontsize=22.)
ax[0].plot(Pleiades_NIR[:,0] - Pleiades_NIR[:,4] - 0.02, Pleiades_NIR[:,4] - 0.01 - dist_mod, 'o', c='#555555', 
        markersize=4.0, alpha=0.1)
ax[0].plot(iso_gs98[:,10] - iso_gs98[:,12], iso_gs98[:,12], dashes=(20., 5.), lw=3, c='#B22222')
ax[0].plot(iso_gas07[:,11] - iso_gas07[:,13], iso_gas07[:,13], lw=3, c='#0094b2')

# (J-H)
ax[1].set_xlim(0.0, 1.0)
ax[1].set_xlabel('$(J - H)$', fontsize=22.)
ax[1].plot(Pleiades_NIR[:,0] - Pleiades_NIR[:,2] - 0.01, Pleiades_NIR[:,4] - 0.01 - dist_mod, 'o', c='#555555', 
        markersize=4.0, alpha=0.1)
ax[1].plot(iso_gs98[:,10] - iso_gs98[:,11], iso_gs98[:,12], dashes=(20., 5.), lw=3, c='#B22222')
ax[1].plot(iso_gas07[:,11] - iso_gas07[:,12], iso_gas07[:,13], lw=3, c='#0094b2')

# (H-K)
ax[2].set_xlim(0.0, 0.5)
ax[2].set_xlabel('$(H - K)$', fontsize=22.)
ax[2].plot(Pleiades_NIR[:,2] - Pleiades_NIR[:,4] - 0.01, Pleiades_NIR[:,4] - 0.01 - dist_mod, 'o', c='#555555', 
        markersize=4.0, alpha=0.1)
ax[2].plot(iso_gs98[:,11] - iso_gs98[:,12], iso_gs98[:,12], dashes=(20., 5.), lw=3, c='#B22222')
ax[2].plot(iso_gas07[:,12] - iso_gas07[:,13], iso_gas07[:,13], lw=3, c='#0094b2')

"""
Explanation: We've loaded isochrones with both Grevesse & Sauval (1998) and Grevesse et al. (2007) solar abundance distributions. The former is similar to Anders & Grevesse (1989) and Grevesse & Noels (1993) distributions in that the solar surface $(Z/X) = 0.023$ &ndash; $0.025$. On the other hand, the Grevesse et al. (2007) distribution has a characteristic $(Z/X) = 0.014$.
Starting with CMDs against $M_K$.
End of explanation
"""


fig, ax = plt.subplots(1, 3, figsize=(15., 8.), sharey=True)

ax[0].set_ylabel('$M_J$', fontsize=22.)
for axis in ax:
    axis.grid(True)
    axis.set_ylim(10., 3.)
    axis.tick_params(which='major', axis='both', length=15., labelsize=16.)

# (J-K)
ax[0].set_xlim(0.2, 1.2)
ax[0].set_xlabel('$(J - K)$', fontsize=22.)
ax[0].plot(Pleiades_NIR[:,0] - Pleiades_NIR[:,4] - 0.02, Pleiades_NIR[:,0] - 0.01 - dist_mod, 'o', c='#555555', 
        markersize=4.0, alpha=0.1)
ax[0].plot(iso_gs98[:,10] - iso_gs98[:,12], iso_gs98[:,10], dashes=(20., 5.), lw=3, c='#B22222')
ax[0].plot(iso_gas07[:,11] - iso_gas07[:,13], iso_gas07[:,11], lw=3, c='#0094b2')

# (J-H)
ax[1].set_xlim(0.0, 1.0)
ax[1].set_xlabel('$(J - H)$', fontsize=22.)
ax[1].plot(Pleiades_NIR[:,0] - Pleiades_NIR[:,2] - 0.01, Pleiades_NIR[:,0] - 0.01 - dist_mod, 'o', c='#555555', 
        markersize=4.0, alpha=0.1)
ax[1].plot(iso_gs98[:,10] - iso_gs98[:,11], iso_gs98[:,10], dashes=(20., 5.), lw=3, c='#B22222')
ax[1].plot(iso_gas07[:,11] - iso_gas07[:,12], iso_gas07[:,11], lw=3, c='#0094b2')

# (H-K)
ax[2].set_xlim(0.0, 0.5)
ax[2].set_xlabel('$(H - K)$', fontsize=22.)
ax[2].plot(Pleiades_NIR[:,2] - Pleiades_NIR[:,4] - 0.01, Pleiades_NIR[:,0] - 0.01 - dist_mod, 'o', c='#555555', 
        markersize=4.0, alpha=0.1)
ax[2].plot(iso_gs98[:,11] - iso_gs98[:,12], iso_gs98[:,10], dashes=(20., 5.), lw=3, c='#B22222')
ax[2].plot(iso_gas07[:,12] - iso_gas07[:,13], iso_gas07[:,11], lw=3, c='#0094b2')

"""
Explanation: Now CMDs against $M_J$
End of explanation
"""


fig, ax = plt.subplots(1, 3, figsize=(15., 8.), sharey=True)

ax[0].set_ylabel('$M_H$', fontsize=22.)
for axis in ax:
    axis.grid(True)
    axis.set_ylim(10., 3.)
    axis.tick_params(which='major', axis='both', length=15., labelsize=16.)

# (J-K)
ax[0].set_xlim(0.2, 1.2)
ax[0].set_xlabel('$(J - K)$', fontsize=22.)
ax[0].plot(Pleiades_NIR[:,0] - Pleiades_NIR[:,4] - 0.02, Pleiades_NIR[:,2] - 0.01 - dist_mod, 'o', c='#555555', 
        markersize=4.0, alpha=0.1)
ax[0].plot(iso_gs98[:,10] - iso_gs98[:,12], iso_gs98[:,11], dashes=(20., 5.), lw=3, c='#B22222')
ax[0].plot(iso_gas07[:,11] - iso_gas07[:,13], iso_gas07[:,12], lw=3, c='#0094b2')

# (J-H)
ax[1].set_xlim(0.0, 1.0)
ax[1].set_xlabel('$(J - H)$', fontsize=22.)
ax[1].plot(Pleiades_NIR[:,0] - Pleiades_NIR[:,2] - 0.01, Pleiades_NIR[:,2] - 0.01 - dist_mod, 'o', c='#555555', 
        markersize=4.0, alpha=0.1)
ax[1].plot(iso_gs98[:,10] - iso_gs98[:,11], iso_gs98[:,11], dashes=(20., 5.), lw=3, c='#B22222')
ax[1].plot(iso_gas07[:,11] - iso_gas07[:,12], iso_gas07[:,12], lw=3, c='#0094b2')

# (H-K)
ax[2].set_xlim(0.0, 0.5)
ax[2].set_xlabel('$(H - K)$', fontsize=22.)
ax[2].plot(Pleiades_NIR[:,2] - Pleiades_NIR[:,4] - 0.01, Pleiades_NIR[:,2] - 0.01 - dist_mod, 'o', c='#555555', 
        markersize=4.0, alpha=0.1)
ax[2].plot(iso_gs98[:,11] - iso_gs98[:,12], iso_gs98[:,11], dashes=(20., 5.), lw=3, c='#B22222')
ax[2].plot(iso_gas07[:,12] - iso_gas07[:,13], iso_gas07[:,12], lw=3, c='#0094b2')

"""
Explanation: And finally, against $M_H$
End of explanation
"""


# log(g) = -0.5; [m/H] = 0.0
Teff = np.array([2500., 2700., 2800.])
Ttau = np.array([2850., 3000., 3200.])
Peff = np.array([5.0e1, 1.3e2, 1.7e2])
Ptau = np.array([3.8e2, 4.7e2, 5.5e2])
Rtau = np.array([1.41e-9, 2.33e-9, 2.61e-9])

gammas = (np.log(Ptau) - np.log(1.0e13))/np.log(Rtau)

"""
Explanation: Spot Temperatures on AGB Stars
Using fairly simple energy equipartition arguments with a polytropic equation of state, I was able to estimate spot temperature contrasts for low mass stars with convective envelopes (G &ndash; M) with reasonble accuracy. Can we do this for AGB stars?
We'll use MARCS model hydrostatic model atmospheres to derive approximations for atmospheric properties where $\tau \approx 1$.
End of explanation
"""


fig, ax = plt.subplots(1, 1, figsize=(8,4))

ax.set_xlabel('Effective Temperature (K)', fontsize=20.)
ax.set_ylabel('Spot Temperature (K)', fontsize=20.)
ax.tick_params(which='major', axis='both', length=10., labelsize=16.)
ax.grid(True)

ax.plot(Teff, Ttau*0.4**((gammas - 1)/gammas), '--o', c="#b22222", markersize=10.)

"""
Explanation: Energy equipartition dictates that the temperature ratio of a spot to its surroundings is roughly equal to $0.4^{(\gamma - 1)/\gamma}$. This leads to the following spot temperatures,
End of explanation
"""


fig, ax = plt.subplots(1, 1, figsize=(8,4))

ax.set_xlabel('Effective Temperature (K)', fontsize=20.)
ax.set_ylabel('Equipartition Strength (G)', fontsize=20.)
ax.tick_params(which='major', axis='both', length=10., labelsize=16.)
ax.grid(True)

ax.plot(Teff, np.sqrt(8.0*np.pi*Ptau), '--o', c="#b22222", markersize=10.)

"""
Explanation: Temperature contrasts are on the order of 86 &ndash; 87%, meaning spots are only slightly cooler than the surrounding photosphere.
Estimating magnetic field strengths required to generate spots is more difficult to assess. One could compare the local Alfven velocity to the convective velocity. Alternatively, following an analogy with the Sun, umbral magnetic field strengths often reach equipartition with the gas pressure outside of a spot. Therefore,
End of explanation
"""


# arbitrary velocity values; rho ~ 10**-9
velocities = np.arange(1.0, 5.0, 0.5) # log-scale

# B-fields to get V_A = 0.1 V_C
B_010 = np.sqrt(4.0e-9*np.pi*(0.1*10**velocities)**2)
B_100 = np.sqrt(4.0e-9*np.pi*(10**velocities)**2)

print B_100

"""
Explanation: These values of the magnetic field strength are two orders of magnitude larger than observed magnetic field strengths. Are they dynamically important? At the risk of over-relying on simple theories, I'll adopt the convective velocities from mixing length theory.
End of explanation
"""


