Example use case: Solar system using the API functionality
We use the API interface to construct a model of the Solar system.
First we must construct the argument string that we pass to binary_c
[1]:
import os
from binarycpython.utils.functions import temp_dir
from binarycpython import _binary_c_bindings
TMP_DIR = temp_dir("notebooks", "notebook_solar_system")
M_1 = 1.0 # Msun
metallicity = 0.02
max_evolution_time = 15000 # Myr. You need to include this argument.
api_log_filename_prefix = TMP_DIR,
orbiting_objects = """\
orbiting_object name=Mercury,M=1MMercury,orbital_separation=1AMercury,orbits=star1,orbital_eccentricity=0.206 \
orbiting_object name=Venus,M=1MVenus,orbital_separation=1AVenus,orbits=star1,orbital_eccentricity=0.007 \
orbiting_object name=Earth,M=1MEarth,orbital_separation=1AEarth,orbits=star1,orbital_eccentricity=0.017 \
orbiting_object name=Mars,M=1MMars,orbital_separation=1AMars,orbits=star1,orbital_eccentricity=0.093 \
orbiting_object name=Jupiter,M=1MJupiter,orbital_separation=1AJupiter,orbits=star1,orbital_eccentricity=0.048 \
orbiting_object name=Saturn,M=1MSaturn,orbital_separation=1ASaturn,orbits=star1,orbital_eccentricity=0.056 \
orbiting_object name=Uranus,M=1MUranus,orbital_separation=1AUranus,orbits=star1,orbital_eccentricity=0.047 \
orbiting_object name=Neptune,M=1MNeptune,orbital_separation=1ANeptune,orbits=star1,orbital_eccentricity=0.009 \
orbiting_object name=Pluto,M=1MPluto,orbital_separation=1APluto,orbital_eccentricity=0.2444,orbits=star1,orbital_eccentricity=0.244 \
orbiting_objects_log 1
"""
# Here we set up the argument string that is passed to the bindings
argstring = """
binary_c M_1 {M_1} metallicity {metallicity} max_evolution_time {max_evolution_time} api_log_filename_prefix {api_log_filename_prefix} {orbiting_objects}
""".format(
M_1=M_1,
metallicity=metallicity,
max_evolution_time=max_evolution_time,
api_log_filename_prefix=TMP_DIR,
orbiting_objects=orbiting_objects
).strip()
output = _binary_c_bindings.run_system(argstring)
#print(output)
There is a lot of data here if you uncomment the print statement!
Let’s split it into a dict of lists of data, one list for each planet, and let’s select only objects that are still orbiting their central star.
[13]:
import re
import pandas as pd
import math
#pd.set_option("display.max_rows", None, "display.max_columns", None)
data = {}
for line in output.splitlines():
#print (line)
match = re.search('Object (\d+) (\S+)',line)
if match:
number = match.group(1)
name = match.group(2)
if not name in data:
data[name] = []
x = line.split()
if x[4] == 'CS1':
x.pop(0) # remove first element of the list "Object" - this is superfluous
x.pop(0) # remove second element of the list "index" - this is superfluous
x.pop(0) # remove third element of the list "name" - this is superfluous (it's the dict key)
x.pop(1) # remove (originally) fourth element "CS1" - this is superfluous (we select this already)
data[name].append(x)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
/tmp/ipykernel_216542/1241251901.py in <module>
14 x = line.split()
15 if x[4] == 'CS1':
---> 16 pop(x[0:2])
17 x.pop(0) # remove first element of the list "Object" - this is superfluous
18 x.pop(0) # remove second element of the list "index" - this is superfluous
NameError: name 'pop' is not defined
Now convert this data to Pandas dataframes
[11]:
dataframes = {}
for planet in data:
dataframes[planet] = pd.DataFrame(data[planet],
dtype=float, # required! argh!
columns=['time',
'mass',
'angular momentum',
'separation',
'period',
'eccentricity',
'temperature',
'angular frequency',
'spin of central object',
'in disc'],
)
#print (dataframes['Earth'])
We now make a plot of the separation (distance from the object to the Sun) as a function of time.
[12]:
import seaborn as sns
import matplotlib.pyplot as plt
# set up seaborn for use in the notebook
sns.set(rc={'figure.figsize':(20,10)})
sns.set_context("notebook",
font_scale=1.5,
rc={"lines.linewidth":2.5})
x = 'time'
y = 'separation'
for planet in dataframes:
#print (dataframes[planet])
df = dataframes[planet]
p = sns.lineplot(data=df,
sort=False,
x=x,
y=y,
estimator=None,
ci=None,
legend='full'
)
xx = df.head(1).loc[0,x]
yy = df.head(1).loc[0,y]
p.text(xx,yy,planet)
p.set_xlabel("time/Myr")
p.set_ylabel("separation/au")
#p.set(xlim=(0,5)) # might be necessary?
p.set(yscale="log")
[12]:
[None]

The inner objects are swallowed by the Sun when it becomes a red giant. Earth survives, although mass loss from the red-giant Sun and tides mess with its orbit somewhat. Jupiter is pushed out beyond the orbits of Saturn and Uranus, and this simple model assumes they are ejected in the interaction because Jupiter is (far) more massive. There are options to detect when its orbit is too close to Neptune, and hence possibly eject Neptune, but I’ll let you explore these.
We now construct a plot of the temperature of each planet vs time.
[5]:
import seaborn as sns
import matplotlib.pyplot as plt
# set up seaborn for use in the notebook
sns.set(rc={'figure.figsize':(20,10)})
sns.set_context("notebook",
font_scale=1.5,
rc={"lines.linewidth":2.5})
x = 'time'
y = 'temperature'
for planet in dataframes:
df = dataframes[planet]
p = sns.lineplot(data=df,
sort=False,
x=x,
y=y,
estimator=None,
ci=None,
legend='full'
)
xx = df.head(1).loc[0,x]
yy = df.head(1).loc[0,y]
p.text(xx,yy,planet)
p.set_xlabel("time/Myr")
p.set_ylabel("Temperature/K")
p.set(ylim=(30,3000)) # might be necessary?
p.set(yscale="log")
[5]:
[None]

It gets a little toasty on Earth in the not too distant future!