# 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!