Solutions¶
Mass histograms
Below is an example solution for stars with $M < M_{\odot}$. First I created a new dataframe with just the low mass stars. I then created the histograms and the caption.
low_mass_df = planet_df.where(planet_df['star_mass'] <=1.0).dropna()
low_mass_df
fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(1,1,1)
ax.hist(low_mass_df['star_mass'].where(low_mass_df['detection_type']=='Primary Transit'), bins=15, label='Primary Transit', alpha=0.6, color='DeepPink', edgecolor='DeepPink')
ax.hist(low_mass_df['star_mass'].where(low_mass_df['detection_type']=='Radial Velocity'), bins=15, label='Radial Velocity', alpha=0.7, color='DodgerBlue', edgecolor='DodgerBlue')
ax.set_xlabel('Host star mass ($M_\odot$)')
ax.set_ylabel('Count')
ax.legend()
ax.set_title('Number of exoplanets detected using transit and radial velocity\n techniques for host star masses below 1 $M_{\odot}$')
plt.show()
Number of exoplanets detected as a function of host star mass and detection technique for low mass ($M < 1 M_\odot$) stars. Pink bars denote exoplanets detected using primary transit observations. Blue bars denote those detected using radial velocity measurements. The fraction of stars detected using primary transit is higher for low mass stars.
Clean your data
kepler_df = planet_df.where((planet_df['orbital_period'] <= 10000) & (planet_df['star_mass'] > 0.7)).dropna()
kepler_df
Check that we've actually got what we want by looking at the data:
kepler_df['orbital_period'].max()
kepler_df['star_mass'].min()
kepler_df['orbital_period'].hist()
kepler_df['star_mass'].hist()
(Yes, I am aware the plots aren't pretty. These aren't supposed to be. Get in the habit of making quick and dirty plots to look at your data)
Estimate $G \pm \sigma_G$
We're using curve_fit, so we need to import it. Remember to add the import statement to the top of your notebook if you don't already have it loaded.
First things first, define the function that we want to fit. We're estimating $G$ by rearranging the equation for Kepler's 3rd law:
$$ T^2 = \dfrac{4 \pi^2}{G} \dfrac{a^3}{M_*} $$Our function is then
def keplers_law(a_3_m, g):
return (4.0 * np.pi**2 * a_3_m / g)
Now we use curve_fit to find $G$ and $\sigma_G$
Note: The first (few...) times I tried this it didn't work. curve_fit got a very strange answer. By providing a starting guess with the p0 parameter the fit worked properly. The problem was the very large numbers ($10^6$ and $10^{17}$). As we can see from the units in the $T^2$ equation, $G$ will be approximately $\dfrac{10^6}{10^{17}} = 10^{-11}$ m$^{3}$ kg$^{-1}$ s$^{-2}$, so p0=1e-11 is a reasonable starting guess.
popt, pcov = curve_fit(keplers_law, (kepler_df['a_3_m']), (kepler_df['T_sq']), p0=1e-11)
g = popt[0]
g_unc = np.sqrt(pcov[0][0])
print("G = {0:.3e} +/- {1:.3e}".format(g, g_unc))
fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(1,1,1)
ax.plot(kepler_df['a_3_m']/1e6, kepler_df['T_sq'] / 1e17, 'k.',ls='None', label='Data')
ax.plot(kepler_df['a_3_m']/1e6, keplers_law(kepler_df['a_3_m'], g)/1e17, 'r--', label='Best fit')
ax.set_ylabel('$T^2$ ($10^{17}$ s$^2$)')
ax.set_xlabel('$a^3 / M_*$ ($10^6$ m$^3$ kg$^{-1}$)')
title_string = (r'$G = ({0:.3f} \pm {1:.3f}) \times 10^{{-11}}$ m$^{{3}}$ kg$^{{-1}}$ s$^{{-2}}$'.format((g / 1e-11), (g_unc / 1e-11)))
ax.set_title(title_string)
ax.legend()
plt.show()
