<jupyter_start><jupyter_text>Simulating Kepler exoplanet orbitsUpload some actual data from the Kepler telescope confirmed planets database to simulate the stability of planetary systems. Start by importing our usual things:<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline<jupyter_output><empty_output><jupyter_text>If we recall we have a lot of columns in this dataset for each planet. We can check them out inline, or just simply open up the<jupyter_code>!head -n 40 downloads/kepler101data.txt<jupyter_output># columns are:
# (1) Row ID from table
# (2) System Name
# (3) planet letter
# (4) Number of total planets in system
# (5) Orbital period in days
# (6) Upper error in Orbital Period (P_orb +/- (upper error)/(lower error)) in days
# (7) Lower error in orbital period in days
# (8) Semi-major axis  in AU
# (9) Error in semi major axis (a +/- ea)
# (10) Eccentricity
# (11) Upper error in eccentricity
# (12) Lower error in eccentricity
# (13) Inclination in degrees
# (14) Upper error in inclination
# (15) Lower error in inclination
# (16) Planet's mass in Jupiter mass (or M*sin(i) if given)
# (15) Upper error in planet's mass in jupiter masses
# (16) Lower error in planet's mass in jupieter masses
# (17) Mass type ('Mass' or 'Msini' - tells you what is the "mass" you actually measure)
# (18) Star's mass in solar masses
# (19) Error in star's mass
# (20) Star's radius in solar radii
# (21) Upper error in star's radius
# (22) Lower error in star's radius
# (23) Transi[...]<jupyter_text>Okay dokey, read in the file! Because there are so many parameters, let's use another library to do this for us.Make sure you've downloaded the kepler 10 and 101 data sets!<jupyter_code>from convert_kepler_data import read_kepler_data
kepler_data = read_kepler_data('/Users/jillnaiman1/Downloads/kepler101data.txt')
kepler_data<jupyter_output><empty_output><jupyter_text>The above are all the different parameters for each planet, as listed by the "head" of the text file.We can also look at individual parameters that might be interesting:<jupyter_code>kepler_data['ecc']<jupyter_output><empty_output><jupyter_text>So, the first planet has a small eccentricity and the second has zero, so a circular orbit. Let's use a library to convert this to units that the Hermite solver can use.  Note: this uses a few assumptions to convert the Kepler data into initial conditions for our simulations.Please check out the library if you want to see these assumptions if you are curious!<jupyter_code>from convert_kepler_data import convert_kepler_data
star_mass, \
planet_masses, \
planet_initial_position, \
planet_initial_velocity, ecc = convert_kepler_data(kepler_data)<jupyter_output><empty_output><jupyter_text>Let's see how these things look:<jupyter_code>star_mass # Msun
planet_masses # in Jupiter masses
planet_initial_position # in AU
planet_initial_velocity # in km/s<jupyter_output><empty_output><jupyter_text>ExerciseUse the Hermite solver to model this system & plot the orbits of the planets. How large/small does your timestep have to be to be to get a stable orbit?Bonus: what about for the 5 planet system?Once you have something you like, make a note of the time parameters you used.  We'll eventually save your favorite model to a file.  Possible ans Now import the hermite stuffs to solve this:<jupyter_code>from hermite_library import do_hermite
# h is for hermite!
r_h, v_h, t_h, e_h = do_hermite(star_mass, 
                                planet_masses, 
                                planet_initial_position, 
                                planet_initial_velocity, 
                               tfinal=1e7, Nsteps=8800)
# let's plot!
fig, ax = plt.subplots(1, 2, figsize = (10*2, 10))
fig.suptitle('Coordinates Plot')

ax[0].set_xlabel('x in AU')
ax[0].set_ylabel('y in AU')

for i in range(len(planet_masses)):
    ax[0].plot(r_h[i,0,:], r_h[i,1,:], lw=3)

ax[1].set_xlabel('Time in years')
ax[1].set_ylabel('Energy')
# re-norm energy
ax[1].plot(t_h, e_h)

plt.show()<jupyter_output><empty_output><jupyter_text>Saving data to a file, and reading it back in!We'll also want to save our models and then be able to read them back in so you don't have to keep running models once you find one you like.Per usual, we'll use a function from a library!<jupyter_code>from hermite_library import save_hermite_solution_to_file<jupyter_output><empty_output><jupyter_text>Now we'll simply call this and save our data!<jupyter_code># you can save it to your local directory, or another place
save_hermite_solution_to_file('myPlanetSystem.txt', 
                              t_h, e_h, r_h, v_h)<jupyter_output><empty_output><jupyter_text>Let's quickly take a look at this:<jupyter_code>!head myPlanetSystem.txt<jupyter_output>0.000000000000000000e+00,-6.799198885028552943e-05,-4.331852307244812156e-02,-1.107107656726798300e-07,0.000000000000000000e+00,6.702900216666182742e-02,1.365104365687684444e-02,0.000000000000000000e+00,5.078102049476259360e-06,-1.107107656726798300e-07,0.000000000000000000e+00,1.201416740318806473e-06,-8.523511935193556210e-01,0.000000000000000000e+00,-1.481390926396770891e-01,7.274383576207541502e-01,0.000000000000000000e+00,1.201416740318806473e-06,1.047159119364062479e-04,0.000000000000000000e+00
1.136492783270826521e+03,-6.799198824681156647e-05,-4.330247022065980556e-02,-1.077698343988626843e-03,0.000000000000000000e+00,6.683539105535182823e-02,1.456950090025424180e-02,0.000000000000000000e+00,5.077589147273007199e-06,2.168645587296400320e-08,0.000000000000000000e+00,2.539002113577943157e-02,-8.520353387632382258e-01,0.000000000000000000e+00,-1.581064578856623548e-01,7.253369660846603617e-01,0.000000000000000000e+00,-2.012641502274654666e-06,1.046919656279508951e-04,0.0000000000[...]<jupyter_text>We want to be able to read this simulation data back in to use in plots, and visualizations next week.  So let's practice doing that with yet another function.<jupyter_code>from hermite_library import read_hermite_solution_from_file
t_h2, e_h2, r_h2, v_h2 = read_hermite_solution_from_file('myPlanetSystem.txt')
# let's plot!
fig, ax = plt.subplots(1, 2, figsize = (10*2, 10))
fig.suptitle('Coordinates Plot')

ax[0].set_xlabel('x in AU')
ax[0].set_ylabel('y in AU')

for i in range(len(planet_masses)):
    ax[0].plot(r_h2[i,0,:], r_h2[i,1,:], lw=3)

ax[1].set_xlabel('Time in years')
ax[1].set_ylabel('Energy')
# re-norm energy
ax[1].plot(t_h2, e_h2)

plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0965316144874156, 'fraction_numerical': 0.14303253529772866, 'mean_word_length': 5.254318618042227, 'pattern_counts': {'":': 0, '<': 43, '<?xml version=': 0, '>': 43, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 13, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '83745', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2929, 'n_tokens_neox': 2304, 'n_words': 643, 'path': 'lesson05/.ipynb_checkpoints/_Part1_programming_readWriteData_keplersData_lesson05-checkpoint.ipynb', 'repo_name': 'jnaiman/csci-p-14110_su2019', 'size': '6516'}