<jupyter_start><jupyter_text>Calculate the rotation distribution for hot stars<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
from plotstuff import colours
cols = colours()
%matplotlib inline
plotpar = {'axes.labelsize': 20,
           'text.fontsize': 20,
           'legend.fontsize': 15,
           'xtick.labelsize': 20,
           'ytick.labelsize': 20,
           'text.usetex': True}
plt.rcParams.update(plotpar)
KID, Teff, logg, Mass, Prot, Prot_err, Rper, LPH, w, DC, Flag = \
    np.genfromtxt("Table_1_Periodic.txt", delimiter=",", skip_header=1).T
m = Teff > 6250
Prot, Rper, Teff = Prot[m], Rper[m], Teff[m]
plt.scatter(Prot, np.log(Rper), c=Teff)
plt.colorbar()
plt.hist(Prot, 50)
N, P_bins = np.histogram(Prot, 50)
m = N == max(N)
ind = int(np.arange(len(P_bins))[m][0] + 1)
plt.axvline((P_bins[m] + P_bins[ind])/2, color="r")
print((P_bins[m] + P_bins[ind])/2)<jupyter_output>[ 2.10053]<jupyter_text>Fit a Gaussian<jupyter_code>def Gaussian(par, x):
    A, mu, sig = par
    return A * np.exp(-.5*(x-mu)**2/sig**2)
def chi2(par, x, y):
    return sum((y - Gaussian(par, x))**2) 
import scipy.optimize as sco
par_init = 300, 2.10053, 5.
x, y = P_bins[1:], N
result1 = sco.minimize(chi2, par_init, args=(x, y))
A, mu, sig = result1.x
print(A, mu, sig)
plt.hist(Prot, 50)
xs = np.linspace(0, 70, 1000)
ys = Gaussian(result1.x, xs)
plt.plot(xs, ys, "r")<jupyter_output><empty_output><jupyter_text>Fit two Gaussians<jupyter_code>def Double_Gaussian(par, x):
    A1, A2, mu1, mu2, sig1, sig2 = par
    return A1 * np.exp(-.5*(x-mu1)**2/sig1**2) + A2 * np.exp(-.5*(x-mu2)**2/sig2**2)
def Double_chi2(par, x, y):
    return sum((y - Double_Gaussian(par, x))**2) 
double_par_init = A, mu, sig, 12, 5, 3
result2 = sco.minimize(Double_chi2, double_par_init, args=(x, y))
A1, A2, mu1, mu2, sig1, sig2 = result2.x
print(result2.x)
print(mu1, mu2)
print(sig1, sig2)
plt.hist(Prot, 50, color="w", histtype="stepfilled", 
         label="$P_{\mathrm{rot}}~(T_{\mathrm{eff}} > 6250)$")  # ,~\mathrm{McQuillan~et~al.~(2013)}$")
ys = Double_Gaussian(result2.x, xs)
ys1 = Gaussian([A1, mu1, sig1], xs)
ys2 = Gaussian([A2, mu2, sig2], xs)
plt.plot(xs, ys, color=cols.blue, lw=2, label="$G1 + G2$")
plt.plot(xs, ys1, color=cols.orange, lw=2, label="$G1:\mu={0:.1f}, \sigma={1:.1f}$".format(mu1, sig1))
plt.plot(xs, ys2, color=cols.pink, lw=2, label="$G2:\mu={0:.1f}, \sigma={1:.1f}$".format(mu2, sig2))
plt.xlim(0, 30)
plt.legend()
plt.xlabel("$P_{\mathrm{rot}}~\mathrm{(Days)}$")
plt.ylabel("$\mathrm{Number~of~stars}$")
plt.subplots_adjust(bottom=.25, left=.25)
plt.savefig("hot_star_hist.pdf")
print(chi2(result1.x, x, y)/(len(x)-3-1), Double_chi2(result2.x, x, y)/(len(x)-6-1))<jupyter_output>107.320483706 18.0529947306
{'fraction_non_alphanumeric': 0.2109689213893967, 'fraction_numerical': 0.0643510054844607, 'mean_word_length': 4.846153846153846, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '12434', 'license': 'MIT', 'n_tokens_mistral': 1339, 'n_tokens_neox': 1200, 'n_words': 239, 'path': 'code/hot_star_rotation.ipynb', 'repo_name': 'RuthAngus/LSST-max', 'size': '2735'}