%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
pd.options.display.max_columns = 150

%config InlineBackend.figure_format = 'retina'

import astropy
from astropy.table import Table
from astropy.io import ascii
import numpy as np

df_abc = pd.read_csv('../data/Rebull2016/Rebull_Stauffer_merge.csv')

df_abc.head()

"""
Explanation: ApJdataFrames Rebull 2016 a & b Extra analysis
Titles:
a.: Rotation in the Pleiades with K2: I. Data and First Results
b.: ROTATION IN THE PLEIADES WITH K2. II. MULTIPERIOD STARS    
Authors: L. M. Rebull, J. R. Stauffer, J. Bouvier, A. M. Cody, L. A. Hillenbrand, et al.
Data is from this paper:
http://iopscience.iop.org/article/10.3847/0004-6256/152/5/114/meta
End of explanation
"""


tab1 = pd.read_fwf('../data/Fang2016/Table_1+4_online.dat', na_values=['-99.00000000', '-9999.0', '99.000, 99.0'])
df = tab1.rename(columns={'#             Object_name':'Object_name'})

df.head()

df.Object_name.values[0:30]

"""
Explanation: Import the Fang et al. 2016 data
End of explanation
"""


ra1 = df.RAJ2000
dec1 = df.DEJ2000

ra2 = df_abc.RAdeg
dec2 = df_abc.DEdeg

from astropy.coordinates import SkyCoord
from astropy import units as u
c = SkyCoord(ra=ra1*u.degree, dec=dec1*u.degree)  
catalog = SkyCoord(ra=ra2*u.degree, dec=dec2*u.degree)  
idx, d2d, d3d = c.match_to_catalog_sky(catalog)  

plt.figure(figsize=(10,10))
plt.plot(ra1, dec1, '.', alpha=0.3)
plt.plot(ra2, dec2, '.', alpha=0.3)

plt.hist(d2d.to(u.arcsecond)/u.arcsecond, bins=np.arange(0,3, 0.1));
plt.yscale('log')
plt.axvline(x=0.375,color='r', linestyle='dashed')

"""
Explanation: Ugh, the naming convention is non-standard in a way that is likely more work than it's worth to try to match to other catalogs.  Whyyyyyyyyy.
Try full-blown coordinate matching
From astropy: http://docs.astropy.org/en/stable/coordinates/matchsep.html
End of explanation
"""


boolean_matches = d2d.to(u.arcsecond).value < 0.375

"""
Explanation: Ok, we'll accept all matches with better than 0.375 arcsecond separation.
End of explanation
"""


boolean_matches.sum()

"""
Explanation: How many matches are there?
End of explanation
"""


df['EPIC'] = ''

matched_idx = idx[boolean_matches]

matched_idx

df.shape, df_abc.shape

idx.shape

df['EPIC'][boolean_matches] = df_abc['EPIC'].iloc[matched_idx].values

fang_K2 = pd.merge(df_abc, df, how='left', on='EPIC')

fang_K2.columns

fang_K2[['Name_adopt', 'Object_name']][fang_K2.Object_name.notnull()].tail(10)

fang_K2.to_csv('../data/Fang2016/Rebull_Fang_merge.csv', index=False)

fang_K2

"""
Explanation: 120 matches--- not bad. Only keep the subset of fang sources that also have K2
End of explanation
"""


plt.figure(figsize=(7, 7))
plt.plot(fang_K2.Teff, fang_K2.Tspec, 'o')
plt.xlabel(r'$T_{\mathrm{eff}}$ Stauffer et al. 2016')
plt.ylabel(r'$T_{\mathrm{spec}}$ Fang et al. 2016')
plt.plot([3000, 6300], [3000, 6300], 'k--')
plt.ylim(3000, 6300)
plt.xlim(3000, 6300);

"""
Explanation: Great correspondence!  Looks like there are 120 targets in both categories.
Let's spot-check if they use similar temperatures:
End of explanation
"""


delta_Tspec = fang_K2.Teff - fang_K2.Tspec
delta_Tspec = delta_Tspec.dropna()
RMS_Tspec = np.sqrt((delta_Tspec**2.0).sum()/len(delta_Tspec))
print('{:0.0f}'.format(RMS_Tspec))

"""
Explanation: What's the scatter?
End of explanation
"""


fang_K2['flux_amp'] = 1.0 - 10**(fang_K2.Ampl/-2.5)

plt.hist(fang_K2.flux_amp, bins=np.arange(0, 0.15, 0.005));

plt.hist(fang_K2.fs1.dropna(), bins=np.arange(0, 0.8, 0.03));

sns.set_context('talk')

plt.figure(figsize=(7, 7))
plt.plot(fang_K2.fs1, fang_K2.flux_amp, '.')
plt.plot([0,1], [0,1], 'k--')
plt.xlim(0.0,0.2)
plt.ylim(0,0.2)
plt.xlabel('LAMOST-measured $f_{spot}$ \n Fang et al. 2016')
plt.ylabel('K2-measured spot amplitude $A \in [0,1)$ \n Rebull et al. 2016')
plt.xticks(np.arange(0, 0.21, 0.05))
plt.yticks(np.arange(0, 0.21, 0.05))
plt.savefig('K2_LAMOST_starspots_data.png', bbox_inches='tight', dpi=300);

plt.figure(figsize=(7, 7))
plt.plot(fang_K2.fs1, fang_K2.flux_amp, '.')
plt.plot([0,1], [0,1], 'k--')
plt.xlim(0.0,0.5)
plt.ylim(0,0.5)
plt.xlabel('LAMOST-measured $f_{spot}$ \n Fang et al. 2016')
plt.ylabel('K2-measured spot amplitude $A \in [0,1)$ \n Rebull et al. 2016')
plt.xticks(np.arange(0, 0.51, 0.1))
plt.yticks(np.arange(0, 0.51, 0.1))
plt.savefig('K2_LAMOST_starspots_wide.png', bbox_inches='tight', dpi=300);

"""
Explanation: The authors disagree on temperature by about $\delta T \sim$ 100 K RMS.
Let's make the figure we really want to make:  K2 Amplitude versus spectroscopically measured filling factor of starspots $f_{spot}$.  We expect that the plot will be a little noisy due to differences in temperature assumptions and such, but it is absolutely fundamental.  
First we need to convert the amplitude in magnitudes to a faction $\in [0,1]$.  The $\Delta V$ in Stauffer et al. 2016 has negative values, so I'm not sure what it is!  The Ampl from Rebull et al. 2016 is: 

Amplitude, in mag, of the 10th to the 90th percentile
End of explanation
"""


fang_K2.columns

fang_K2.beat.value_counts()

"""
Explanation: Awesome!  The location of points indicate that starspots have a large longitudinally-symmetric component that evades detection in K2 amplitudes.
What effects can cause / mimic this behavior?
- Unresolved binarity could cause an errant TiO measurement, biasing the Fang et al. measurement.
- Increased rotation (Rossby number) could make stronger or weaker dipolar magnetic fields
- EW H$\alpha$ could be correlated, from an activity sense.
-
End of explanation
"""


plt.figure(figsize=(7, 7))
cross_tab = 'resc'
c1 = fang_K2[cross_tab] == 'yes'
c2 = fang_K2[cross_tab] == 'no'
plt.plot(fang_K2.fs1[c1], fang_K2.flux_amp[c1], 'r.', label='{} = yes'.format(cross_tab))
plt.plot(fang_K2.fs1[c2], fang_K2.flux_amp[c2], 'b.', label='{} = no'.format(cross_tab))

plt.legend(loc='best')

plt.plot([0,1], [0,1], 'k--')
plt.xlim(-0.01,0.8)
plt.ylim(0,0.15)
plt.xlabel('LAMOST-measured $f_{spot}$ \n Fang et al. 2016')
plt.ylabel('K2-measured spot amplitude $A \in [0,1)$ \n Rebull et al. 2016')
#plt.xticks(np.arange(0, 0.51, 0.1))
#plt.yticks(np.arange(0, 0.51, 0.1))
plt.savefig('K2_LAMOST_starspots_crosstab.png', bbox_inches='tight', dpi=300);

"""
Explanation: Crosstabs with discreate variables: Legend
End of explanation
"""


plt.figure(figsize=(7, 7))
cross_tab = 'Mass'
cm = plt.cm.get_cmap('Blues')
sc = plt.scatter(fang_K2.fs1, fang_K2.flux_amp, c=fang_K2[cross_tab], cmap=cm)
cb = plt.colorbar(sc)
#cb.set_label(r'$T_{spot}$ (K)')

plt.plot([0,1], [0,1], 'k--')
plt.xlim(-0.01,0.8)
plt.ylim(0,0.15)
plt.xlabel('LAMOST-measured $f_{spot}$ \n Fang et al. 2016')
plt.ylabel('K2-measured spot amplitude $A \in [0,1)$ \n Rebull et al. 2016')
#plt.xticks(np.arange(0, 0.51, 0.1))
#plt.yticks(np.arange(0, 0.51, 0.1))
plt.savefig('K2_LAMOST_starspots_cb.png', bbox_inches='tight', dpi=300);

"""
Explanation: Crosstabs with continuous variables: Colorbar
End of explanation
"""


import astropy.units as u

sini = fang_K2.vsini * u.km/u.s * fang_K2.Per1* u.day /(2.0*np.pi *u.solRad)

vec = sini.values.to(u.dimensionless_unscaled).value

sns.distplot(vec[vec == vec], bins=np.arange(0,2, 0.1), kde=False)
plt.axvline(1.0, color='k', linestyle='dashed')

inclination = np.arcsin(vec)*180.0/np.pi

sns.distplot(inclination[inclination == inclination], bins=np.arange(0,90.0, 5), kde=False);

fang_K2['sini'] = vec

plt.figure(figsize=(7, 7))
cross_tab = 'sini'
cm = plt.cm.get_cmap('hot')
sc = plt.scatter(fang_K2.fs1, fang_K2.flux_amp, c=fang_K2[cross_tab], cmap=cm)
cb = plt.colorbar(sc)
cb.set_label(r'$\sin{i}$')

plt.plot([0,1], [0,1], 'k--')
plt.xlim(-0.01,0.7)
plt.ylim(0,0.15)
plt.xlabel('LAMOST-measured $f_{spot}$ \n Fang et al. 2016')
plt.ylabel('K2-measured spot amplitude $A \in [0,1)$ \n Rebull et al. 2016')
#plt.xticks(np.arange(0, 0.51, 0.1))
#plt.yticks(np.arange(0, 0.51, 0.1))
plt.savefig('K2_LAMOST_starspots_cb.png', bbox_inches='tight', dpi=300);

"""
Explanation: What about inclination?
$$ V = \frac{d}{t} = \frac{2 \pi R}{P} $$
$$ V \sin{i} = \frac{2 \pi R}{P} \sin{i}$$
$$ V \sin{i} \cdot \frac{P}{2 \pi R} = \sin{i}$$
$$ \arcsin{\lgroup V \sin{i} \cdot \frac{P}{2 \pi R} \rgroup} = i$$
End of explanation
"""


