From 1edcae4cb9d9d30897eb8f367b8ec12994323826 Mon Sep 17 00:00:00 2001
From: David Hendriks <davidhendriks93@gmail.com>
Date: Wed, 11 Sep 2019 16:59:35 +0100
Subject: [PATCH] fixed functionality for fetching default args via binaryc
 interface. making some evolution plots now

---
 binaryc_python_utils/functions.py             |  25 ++--
 .../single_star_full-checkpoint.ipynb         | 133 ++++++++++++++++++
 david_calculations/single_star_full.ipynb     | 133 ++++++++++++++++++
 david_results/s.png                           | Bin 0 -> 17847 bytes
 tests_david/testing_automatic_log_readout.py  |  17 +--
 5 files changed, 291 insertions(+), 17 deletions(-)
 create mode 100644 david_calculations/.ipynb_checkpoints/single_star_full-checkpoint.ipynb
 create mode 100644 david_calculations/single_star_full.ipynb
 create mode 100644 david_results/s.png

diff --git a/binaryc_python_utils/functions.py b/binaryc_python_utils/functions.py
index ace2c62cb..eb22fda38 100644
--- a/binaryc_python_utils/functions.py
+++ b/binaryc_python_utils/functions.py
@@ -7,7 +7,6 @@ def create_arg_string(arg_dict):
     """
     Function that creates the arg string
     """
-
     arg_string = '' 
     for key in arg_dict.keys():
         arg_string += "{key} {value} ".format(key=key, value=arg_dict[key])
@@ -27,17 +26,26 @@ def get_defaults():
             key, value = default.split(' = ')
 
             # Filter out NULLS (not compiled anyway)
-            if not value=='NULL':
-                default_dict[key] = value
+            if not value in ['NULL', 'Function']:
+                if not value=='':
+                    default_dict[key] = value
     return default_dict
 
+def get_arg_keys():
+    """
+    Function that return the list of possible keys to give in the arg string
+    """
+
+    return get_defaults().keys()
+
 def run_system(**kwargs):
     """
     Wrapper to run a system with settings 
     """
 
     # Load default args
-    physics_args = get_defaults()
+    args = get_defaults()
+    # args = {}
 
     # For example
     # physics_args['M_1'] = 20
@@ -46,11 +54,12 @@ def run_system(**kwargs):
 
     # Use kwarg value to override defaults and add new args
     for key in kwargs.keys():
-        physics_args[key] = kwargs[key]
+        args[key] = kwargs[key]
 
     # Construct arguments string and final execution string
-    arg_string = create_arg_string(physics_args)
+    arg_string = create_arg_string(args)
     arg_string = f'binary_c {arg_string}' 
+    # print(arg_string)
 
     # Run it and get output
     buffer = ""
@@ -58,7 +67,6 @@ def run_system(**kwargs):
 
     return output
 
-
 def parse_output(output, selected_header):
     """
     Function that parses output of binaryc when it is construction like this:
@@ -98,5 +106,4 @@ def parse_output(output, selected_header):
         for key in keys:
             final_values_dict[key].append(value_dict[key])
 
-    return final_values_dict
-
+    return final_values_dict
\ No newline at end of file
diff --git a/david_calculations/.ipynb_checkpoints/single_star_full-checkpoint.ipynb b/david_calculations/.ipynb_checkpoints/single_star_full-checkpoint.ipynb
new file mode 100644
index 000000000..30c50c50e
--- /dev/null
+++ b/david_calculations/.ipynb_checkpoints/single_star_full-checkpoint.ipynb
@@ -0,0 +1,133 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import os, sys, time\n",
+    "\n",
+    "import matplotlib.pyplot as plt\n",
+    "from collections import defaultdict\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "\n",
+    "# sys.path.append('../')\n",
+    "import binary_c\n",
+    "\n",
+    "\n",
+    "from binaryc_python_utils.functions import create_arg_string, parse_output, run_system\n",
+    "\n",
+    "result_dir = '../david_results/'"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Execute command and parse args"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Took 0.20s to run single system\n",
+      "The following keys are present in the results:\n",
+      "dict_keys(['t', 'mass', 'zams_mass', 'stellar_type', 'prev_stellar_type', 'metallicity', 'probability', 'dM_in_timestep', 'mdot', 'core_mass', 'core_radius', 'luminosity', 'radius', 'dt', 'dtm', 'Teff', 'omega', 'vwind', 'drdt', 'tm', 'tn', 'tkh', 'angular_momentum', 'he_core_mass', 'CO_core_mass', 'GB_core_mass', 'v_eq', 'v_eq_ratio'])\n"
+     ]
+    }
+   ],
+   "source": [
+    "start = time.time()\n",
+    "output = run_system(M_1=10, M_2=20, separation=0, orbital_period=100000000000)\n",
+    "result = parse_output(output, 'DAVID_SINGLE_ANALYSIS')\n",
+    "stop = time.time()\n",
+    "print(\"Took {:.2f}s to run single system\".format(stop-start))\n",
+    "print(\"The following keys are present in the results:\\n{}\".format(result.keys()))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Cast data into pandas framework and do analysis"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#### Now do whatever you want with it: \n",
+    "\n",
+    "# Cast the data into a dataframe. \n",
+    "df = pd.DataFrame.from_dict(result, dtype=np.float64)\n",
+    "\n",
+    "# Get last change moment\n",
+    "last_st = df['stellar_type'].unique()[-1]\n",
+    "last_stellar_type_change_time_1 = df[df.stellar_type==last_st]['t'].iloc[0]\n",
+    "\n",
+    "# slice to get that last time\n",
+    "sliced_df = df[df.t < last_stellar_type_change_time_1] # Cut off late parts of evolution"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZgd5X3g+++v6qy9qxdtLbV2gcHIYGQwxI4J18Z4ggzBG8STZGzGGu4NXjJ37sR2PLYz9/EymXFu4oHYlwAx9tgQOw4YJsIwjjHYMQaJVRKLEEJLtyRaUku9d5/tN39Une7Tre7TZ+3T5/Tv8zz1VJ23qt56j86j+vW71FuiqhhjjDGzcSpdAGOMMQubBQpjjDFZWaAwxhiTlQUKY4wxWVmgMMYYk1Wg0gUoh/b2dl27dm2li2GMMVXj6aefPqmqHTPtq8lAsXbtWnbt2lXpYhhjTNUQkUOz7bOmJ2OMMVlZoDDGGJNVTQUKEdkmIrf39/dXuijGGFMzaipQqOqDqrq9ubm50kUxxpiaUVOBwhhjTOlZoDDGGJNVTQUK66MwxpjSq6nnKFT1QeDBrVu3fqKgDPb/M3TvgmAEAhEIhCEQ9dZBfx2ITC4zHSdS2i9ljDEVVlOBomg//RycfKW4PNxwRgDJCCQTgSYj4OQakKYcO0uerv2UxpjysLvLFApv2gbXfRsS45AYm7rEp31OjEN89Oxj4xn7E6NTjxsbmOX80eKKLm5uwScY9T/762DdtPVsaRmfAxGrORmziFigmM4JQLjBW+aTKiRjpQtI088fG4BEr3/eKMRH/GPHCivvTMHkrOCTLehM2xeq95Zgeh21YGTMAmGBYqEQ8f/qD0NkHp8DSaW8IBPPXEamrTMDy/Rjph0fG4bhk2fnkRzPs2CSETzqINTgf04HlYaMANMwmZ4ONNMXC0DGFMwCxWLnOJM303JKJSdrMNODSGwE4sNekMlc4iMQG/I/+9tjAzB43E/3A1M+zXbiQKgRwplLw+R2aHp6k5/eMO2cRi+oG7MI1FSgEJFtwLaNGzdWuihmOsctX5NeKjlZmzkr2GR+9oPO+ODUZWwABo5OTUPnvq4bygggTVMDTrjRqxlGmiHSMrmOtmSkN1uwMVWhpgJF0cNjTXVy3MmbcymkUl7gSQeN2LTAMj4E4wMZ+4cmt4d64dRr3v7RM5CKZ79WIDoZNKYHkUjL2fuiS6CuDaKtXnObMfOgpgKFMSXhOBm1nxWF56PqNbWN9XtBY6w/YznjL/1T9w/1wslXJ/dpavb8A1Goa/WCRl1rxnbbtO0lk9vhRuujMXmzQGFMuYhMjvBqXJ7/+ap+v0w6kJyB0dMwcgpG+mC0z1unt4/v9rdPM2vTmRP0aiX17d7SsAzql0JDh79eCvUdk2s3WNQ/gakNFiiMWahEJpvUmlflfl4q6QWXkT4vqEwElGnbQ73QvdNbx0dmzivaOi14LIXGZdDUCU0rvaVxpff8jqlZFiiMqTWOO9kURY4DO8aHYLgXhk74a39Jbw+fgJ5n/KAyfPb5dW1ewEgHj8xA0tTpBTrrU6laFiiMMZN9Mq3r5z52bAAGj3kjxSaWHj+tB3p2eTWW6eqXwpI1sGQttKzxtlv8z02dNg3NAma/jDEmP5Emb+k4Z/Zj4mMZwaQHzhyC04e89ZEnYc8/giYnjxfXq3UsWQvtm6B9s7du2+QFEaemJrquOgs+UIjIm4BPA+3AP6vqtypcJGPMXIIRaF3nLTNJxr0AcvrgZAA5fRD6XocXfugNL57Iqw7aNvrBww8gy86H1g1WC5knFflXFpG7gGuAXlV9c0b61cBfAy5wh6p+XVVfAm4WEQf4LmCBwphq5wa92sOStWfvU/WHCe/zl1e9dfdTsOfHTIzoCkRg6Ztg2Zth+QXesuz8+Z0CJ63vdW89W2CscpUKx98BbsW78QMgIi5wG/AeoBvYKSIPqOqLIvJ+4P8EvleBshpj5pOIN7KqcRmse+fUfbEROLUf3tjrDQd+Yze8/E/wbMatoWUNdL4VVr0NOrfCireUf1TWNy/01l+uzZemVSRQqOrjIrJ2WvIlwH5VPQAgIvcC1wIvquoDwAMi8k/AD2bKU0S2A9sBurq6ylRyY0xFhepgxRZv4UYvTdXrDzm+B46/4C3du2Dvfd5+J+DVNjq3wupLYM1vQXNnxb5CNVpIDXydwJGMz93ApSJyBXA9EAZ2zHayqt4uIseAbaFQ6OJyFtQYs4CITA7F3XzVZPrgcS9g9Ozy1s/9AHb+rbevdT2sfae/vAOaingCfxFYSIFiRqr6C+AXOR5rcz0ZYzyNy+FN13gLeA8ivrEHDv4KXv8l7L0fnrnb29e2ETZcCZve6wUOe4BwioUUKHqA1RmfV/lpObPZY40xs3Jcr79ixVvgsj/2AsfxF/zA8Tg88z146nZvDq3174JNV8Hm9879VHxshgcQK+DZw6d59vAZPvr2LsIBt6R5L6RAsRPYJCLr8ALEDcDv55OB1SiMMTlzXFh5kbdc/knv3SgHfwX7HoZXH4Z9P4V/wuvbOP86OO86aFl9dj7P3zPvRZ/JY/tO8Fc/e5U/uGxNyfOu1PDYe4ArgHYR6Qa+pKp3isgtwMN4w2PvUtW9eeZrNQpjTGGCUdj0Hm/R/+oNyX1lh9dE9cgXvGXV27yAcf51Xk0jlYLffNs7f9mbs+dfZqeHYzRHgwTd0j+cWKlRTzfOkr6DLB3WOeRrNQpjTPFEvCfPO86Bd/wJ9B3wAsbe++CRP/OWVZd4o69Oveqdk21K+HlwajhGa32oLHkvpKYnY4xZmFrXwzv/vbeces0LGHvvh513eJMhrrzQCyYVdHokxpK68kwLX1OBwpqejDFl17YBfvs/eMup17ynzP/XF73O8Qo63j/GOctL9JbHaWpqpi1VfVBVtzc3V+ARfmPM4tO2AVq6vEkNK9j0pKocPTPGyuZoWfKvqUBhjDEVIc7U2XDn2emROKPxJCtbLFDMSUS2icjt/f21Od+KMWaBcipbozh6ZhTAAkUurOnJGFMR4nhzTlVIjx8oVi2xQGGMMQuTSEU7s3tOW40iZ9b0ZIypCHE5PTLGV3e8VJHLHz0zSiTolG14bE0FCmt6MsZUhDgkEkluf/wAvzkww/vCy+xo/ygrW6KISFnyr6lAYYwxFSEODl5n9qfvfXaiz2C+9JwZo7NMzU5ggcIYY4rnuDgo79rcwUgsyR/d9RS9A2Pzdvme06MWKIwxZkETBwflzZ1N3PGHWzl6ZpRrb/sXnnq9r+yXHosnOTk0XraObKixQGGd2caYivCbnlwRLl3fxo9uvgzXET78/z/Bzd97msf2nWAsXp5RUcf7vZpLOQNFTc31ZLPHGmMqQfEDheP97X3+ymYe+ZPf5lu/eI3vPnGIn+49TtAVOluitNaHCDgOdWGXr/zeBUU3GaX7Q8rZ9FRTgaJYn4/EeGTkGQI/eDuCN3ogvZ5cycTIgunHZI44mJ425ViZ/ZiZ8pp+7Gx5z3q96WXLckxmXnOVbcZ/m3zKNi1vEcERZ6Icjjg4OJPpGcfMmp75OeOY6enTr5OZL8KUY9LnTllnljF9nH9Oegk4AW8t3tp1XFxxJ9fTts86x3HmPM4Vt2wjXUzu/rDvV8RXLuEdGW00daEA//dV5/B/XbGR3xw4xc6DfRzuG6F/NM5YPMkvXjnBL17p5aOXFveiIQsU82y3q4yjfGjj701JV7wnLlV1ynbmvpnSph87U9pMeWfL76zrTqzyLNtMZc3Ia9bvPcN1czlmpuuq6lmfU6S8taYm0zQ1e3rmZ/zP04/LSFd0yjEpUqBMHFON0gEjHVAccQg6QQJOgKATnHnbnX1/5jHZ8gm6QSJuhJAbIuJGCAfCs34OOLV9q3ku3geREO9yzm7Nj4ZcfufcpfzOuUsn0lIp5fwvPcyBE8W/RrXn9CgisKw5XHRes6ntXy9PovBedwl/esmfVroopkJyCTITgYaz09NLUpPekkqS0hQJTZBKeemJVGIyTVMkU8nJ4zPP8Y/LTE9vZ+6f6bhEKkFCE8STceKpOIlUgnhqcjuWjDEcG55x3/Tjpv8BU4iABAgHwoTdySUSiExsRwNR6oJ11AfrqQ/UT2zXBeuoC/jpwXrqAnWT+/z0hVSjyvXlco4jrGuv57UTQ0Vf80jfCMubIiV/T3YmCxTGZBARXCnff7hqlEwlZw0i48lxxpJjxJIxxhJjOX2eWBLevvHkOAMjAwzHhxmJjzCSGGE0kdtzCK64NIWaaAo30RxqpjHcSHOomaZQE83hZprDk9utkVbaom20RdqIBCJl+bdyZ6hRzGZ9Rz3Pd58p+pqH+kboaq0rOp9saipQ2IuLjCk91/GatOZTMpVkNDHqBY/ECCPxkYntiXVsmIHYgLeMD9Af66d/rJ/DA4fpH+9nMDY4a22oMdjoBQ0/cLRH22mPtrOsfhkr6lewsmElS+uWEnTymxIj4OReu1nfXs+O3ceIJVKEAoUPQD10aoT/I6NZqxxqKlDYqCdjaoPruDSEGmgINRScR0pTDMYGJ4JI31gfJ0dPcmr0lLce89b7Tu/jiaNPMBgfnHK+Iw5L65aysn4ly+uX09nQyerG1axpWsOapjW0RlrPavbKI06wtr2elMLhvhE2Li3sew6NJzg5NE5Xm9UojDEmb444E81Pq1k95/FjiTGODx/n6PBRjg0d4+jwUe/z0FGeP/E8jxx8hIQmJo5vDDbS1dRFR7RjIi2fisHa9noADp4cLjhQHD41AsAaCxTGGFN+kUCEtc1rWdu8dsb9iVSCY0PHODhwkEMDhzg0cIiDAwc5PnJ84phAHjWKdW1+oDhV+Minw33euWta6wvOIxcWKIwxJgcBJ8DqptWsblrNO3nnlH1//sMP8PDQy7iS+wixJfUhmqPBogLFQb9GUe6mp5qawsMYYyrBxSEp+dUowGt+ev1k4YFif+8QHY1hmqPleQ9FmgUKY4wpkiMOCSTn5yjSNnY08OobhT9L8crxQc5Z1ljw+bmyQGGMMUWarFHk93Dim1Y00js4zqmh8byvmUwpr/YOcs5yCxSIyHUi8rci8vciclWly2OMMdM5OCRFcCS/aWDOXd4EeDWDfB3uG2EsnqrdGoWI3CUivSKyZ1r61SLyiojsF5HPAqjq/ar6CeBm4COVKK8xxmTj+LdSR/KbSvzcFd5N/qUCAkU6uNRyjeI7wNWZCSLiArcB7wPOA24UkfMyDvmCv98YYxYUR7xbqZJfoGhvCNPeEGZvT/7v0NnT04/rCJtrtUahqo8D01/9dAmwX1UPqGoMuBe4Vjz/BXhIVZ+ZLU8R2S4iu0Rk14kTJ8pXeGOMmSZdoxDieZ97UVcLzx7Jf86nZ4+c5k0rGomGyj+9ykLqo+gEjmR87vbTPgm8G/igiNw828mqeruqblXVrR0dHbMdZowxJefg3awlzxoFwMVrlvD6yeG8OrSTKeX5I/1ctHpJ3tcrxEIKFDNS1W+q6sWqerOqfjvbsfYqVGNMJQjppqfEHEee7eI13s3+mcO51yr2vTHI0HiCC1e35H29QiykQNEDUyZkWeWn5UxVH1TV7c3NzSUtmDHGZDPZ9JR/jeKCzmZCAYcnXjuV8zm/fNVrXr98Y1ve1yvEQgoUO4FNIrJORELADcAD+WRgNQpjTCU4RdQoIkGXy9a38egrvTmf8/i+k2xe1sCK5vK9/jRTpYbH3gM8AZwjIt0icpOqJoBbgIeBl4AfqurefPK1GoUxphKK6aMAuPLcpbx+cjinN971j8Z56vU+3rV5/vpiKzXq6UZVXaGqQVVdpap3+uk7VHWzqm5Q1a/km6/VKIwxlSB4kzylCqhRAFx1/jIcgfufnbu1fcfuY8SSKa7ZsrKgaxViITU9Fc1qFMaYSii2RrGiOcq7Nnfww11HiCdnf7pbVbn3qcNsXNrAllXzd5+rqUBhNQpjTGUU9sBdpj+8bC1vDIxzz1OHZz3m0Vd6eb67n4//1rqz3q5XTjUVKKxGYYypBEf9QKGFB4orzungsvVt/LeHX+HADH0VJwbH+cJ9e1jXXs8HL15V8HUKYS8uMsaYIonf9FTIqKeJPET4iw9u4f23PcYH7vwRH7y4i4tXnIuIS/fpUe761eucHolz7/aLCeXzztUSqKlAISLbgG0bN26sdFGMMYuIlKDpCeDI6HM0bPw6Z2J93HsU7jkSIX7mbcT63sl5S1dx6+9fxFvm6SG7TNb0ZIwxRZKJpqfCaxR7Tu7hkz//JB31bfzXS7/If1pzI5d3XESk/de0nfuXfPQ9R7ioa36m7JiupmoUxhhTGV7TU6rAPoqUpvjPT/xnlkSWcOeGf82SH98C8WE+DHSv/22+2tLE1576GpFAhOs3XV/CcuempmoUxhhTCcXM9QTws0M/46W+l/j0+R9jyQOfgbb18Af3w+/8GauOv8g3n36IyxrX89Unv8qBMwdKWfSc1FSgsOGxxpjKKK6P4kf7fkRnQye/2/0KxAbhg38HG34H3vUf4VPPElhzOV996TeEJMA3nv5GKQuek5oKFNZHYYypCPWanpIFND0dGzrGk8ee5Nr123Ce+z6c+7vQvmnygEgTfOhu2p0wH0/V8Xj347x46sVSlTwncwYKEblMRG4TkRdE5ISIHBaRHSLyxyJid2RjjKHwzuyfH/k5ivK7kZUwchLe/IGzD6pvg8tv4cOvP0PECfGjfT8qtsB5yRooROQh4N/iTdR3NbAC7zWlXwAiwE9E5P3lLqQxxixkqVT6Vjr79Buz+Zeef6GrsYuuw7vADcHGd8984IUfpUmFq6KdPPT6Q8SSscILnKe5ahR/oKo3qeoDqnpUVROqOqSqz6jqN1T1CuDX81DOnFgfhTGmEtSfFFDQvM4bT46z8/hOLl95ORx5EjovhvAs78Bu7oSuy7iq7w2G48PsOr6r2GLnLGugUNWTc2WQyzHzxfoojDGVkCI971J+geL53ucZS47xW8svgWPPw6qt2U/YcCWXHnuFiBvmF92/KKishZir6WlQRAYylsHM9XwV0hhjFrLUxAN3+TU9PX/ieQAuSgUhGYNVl2Q/YcOVRFS5tGENvz46f405WR+4U9VZ6kDGGGPSUn5FIt9A8cLJF1jbtJbm3pe9hLlqFMsvADfMW5Mujw0d4tToKdqi5X8das7DY0XkLSJyi79sKWehjDGmmkzUKPJoelJVdp/YzQXtF8CJlyHSAo0rsp/kBmH5BVw04L1f+7kTzxVc5nzkFChE5NPA94Gl/vJ9EflkOQtmjDHVIl2jSKVyf47i+PBxTo2d4oKOC+DkPug4B3J5x0TnWznv6EsEnSDP9S6gQAHcBFyqql9U1S8Cbwc+Ub5iGWNM9UilvBt8PjWK1wdeB2BTyyYvULRvzu3E5VsIx4fZ3LiGl/tezrushcg1UAhMeTY96actKDY81hhTCan0A3ep3Psojg8fB2CFG4XhE16NIhd+QNkcamXf6X35FbRAuQaKvwOeFJEvi8ifA78B7ixfsQpjw2ONMZWQUvHXuTc9jcRHAGgY8AIG7bkGCm96j80aoG+sj5Oj5X9CIadpxlX1L0XkF8A7/KSPqeqzZSuVMcZUkaQKiDddeK7iqTgAwTPdXkLbhtxOrGuFunY2jXmBZt/pfbRH2/Mqb75y7czeAOxV1W8Cu4F3isj8v2bJGGMWoHSNIp8+iolAMejXKJo6c79g+yY2978BwKunX839vALl2vT0YyApIhuBbwOrgR+UrVTGGFNFkv7wWDT/QBHoPwL1SyEYyf2CbRtY0neIlnALhwYO5VPUguQaKFLqTYt4PXCrqv4/eBMEGmPMopeaiA95ND0l4wSdIDLQA82r8rtgcxcMHaerYRWHBw/nd24Bcg0UcRG5EfhD4H/6acHyFMkYY6pLYuLJ7PxqFEEnCGeOQMvq/C7oB5bVkTaODBzJ79wC5BooPgZcBnxFVV8XkXXA98pXrEkisl5E7hSRf5iP6xljTK76hmMMjMUnphkvKFD0d0NzYYGiy63n2PCxsk85nuuopxeBT2V8fl1EflXoRUXkLuAaoFdV35yRfjXw13hvKr9DVb+uqgeAmyxQGGMq6cxIjD09A+zu6Wd3zxle6O6n+/Qoq1ujXBZNQiS/uZ68QOFCYrTgQLEaF0XpHupmffP6/PLIQ9ZAISIu8GGgE/ipqu4RkWuAzwNR4KICr/sd4Fbgu9OudRvwHqAb2CkiD/hByhhj5s2poXF29/Sz9+gAe3r62XO0nyN9oxP7u1rreMvqFq7ZspJvP/Yaze4wNOU56ikZJ5h+bjnvPgo/UIx7NYmewZ7KBQq8h+pWA08B3xSRo8BW4LOqen+hF1XVx0Vk7bTkS4D9fg0CEbkXuBawQGGMKYvh8QSvnRhi3xtDvNo7yKtvDPHSsQGO9Y9NHLOmrY4tnS38/iVrOH9lE1tWNdNSF5rYf2YkxgvP+FNp5FujSMeV5jyGxgIEwtCwjGWjZwDoHenN7/w8zRUotgJbVDUlIhHgOLBBVU+VoSydQGavTDdwqYi0AV8BLhKRz6nq12Y6WUS2A9sBurq6ylA8Y0w1UlV6B8c5eHKYQ30jvNY7xL43Bnm1d4ju05O1hJDrsL6jnkvWtXJBZzPnr2zmvJVNNEezj9v57PvO5cZnfuZdK8/nKILp4xuW5//FmlfRPnACQXhj5I38z8/DXIEipn6jm6qOiciBMgWJWfnXuzmH424XkWPAtlAodHH5S2aMWShiiRTH+kc53DfCwVMjHD417K9HONQ3zFh88i/9UMBhQ0cDb+1awke2rmbTskY2LWtgTWsdATfnNy9MaKkL8Vc3vJUPPvN3+Xdmp8fV1hfwZHXDcoKnD9LW3lbxGsW5IvKCvy3ABv+zAKqqpXwvRQ9eM1faKj8tZ6r6IPDg1q1bbWZbY2pEKqWcHB7n2Jkxjp4ZpefMKMf6ve2j/vrk0PiUZ93CAYc1bXWsaavnnZvaWdNez5rWOta01dHZEi0oIGSzekkDUEBntiahrs17z0S+GpZC91Ms7bqY4yPH8z8/D3MFijeV9epT7QQ2+UNve4AbgN/PJwMR2QZs27hxYxmKZ4wppVgixcmhcU4MjtM76K297TFve2ic3gEvLZacegOOBl1WtkRY2RLl3HOWssLf7mqtY21bPUsbwzjO/E1wLU7+r0JNJBMEkgnvqexCNCyD4ZMsiy6lezivv6nzNlegOKxz1KVEROY6ZoZz7gGuANpFpBv4kqreKSK3AA/jDY+9S1X35pOv1SiMqZyxeJLTIzH6hmOcHo7TNxLj9LD/eWRynQ4MZ0biM+bTWh+ioyFMR2OYS9fV09EUprMlyormqBccmqO01AWRXF7yM0+8QZv59VGkSCHJODTkOeIpraEDUJaFmni695nC8sjRXIHiURH5MfATVZ14TlxEQngzyf4R8CjecNecqeqNs6TvAHbkk1cmq1EYUxxVZTiWZGA0Tv9onIHROANjicnPY3HOjMSn3PhPD8fpG44xGp99iu2WuiCtdSGW1IdY397Apeva6GgMs7TRCwjppb0hTLDEzULzQST/uZ5UFScZh6ZlhV20wTtvmRNhIDbAaGKUaCBaWF5zmCtQXA18HLjHbxI6A0Tw/uJ/BPirhTTduNUozGIXS6QYHk8w5C/D4wkGxxIMjM1+4x8YTWRsxzPmLZpZYzhAa0OIJXUhljZGOGdZE631QVrqQrTWe+mt9SFa64MsqQvRHA2WvE9goRHyf2d2SlO4yZjX11AIP1As9W/jvSO9rGlaU1hec8gaKFR1DPgb4G9EJAi0A6OqeqYspTFmERpPJBkaSzA8npxyg59+sx8eTzAcm9z29ienHBtLzN1GHg26NEUDNEWCNEWDtDeE2NBRT1M06KcFaJ7Y9tbNUS+9IRyo+Zt+ISb7KPJohU8lcTRVRKDwzmuIe898DMeHC8snBzlN4QGgqnHgWNlKUgLW9GTKRVUZi6cYiSUYiSUZjXs36NFYkpFYkpF4kpHxqftGYklvf8a+kXiS0ZgXFEZi3s09nszt5lIXcqkPB2gMB6gPB6gPu3S2RGkIuzREAtP2TW57N3tv3RgJEA64Zf7XWoQk/87sVHLc2yi0M9s/LzI+BMB4Or8yyDlQVANrelqcVJVYMsVYLMVYwrs5T6zjKcbiScbi3g18LJ7y194y4t/oR/0A4C2Z2/6+eDKf5mdcR6gLuf4SIBr0tpujQVY0RSZu+vXhAI2RAPUhl4ZIkIawl96QXvwAUB8K4M7jKB6Tn3QfRT5NT5pK4KAQLfAdcKE6CNYRjXk1idHE6BwnFK6mAoVZOJIpJZZIMZ5IMp5ITb1R+zfysYkb+uT+8Vlu6JmfveOmfs5v3J1HBOqCLtGQ99d5+mZeFwrQ3hCmLuTtqwu51GdsR/0AUB8KTGynz0vvD7nOghqVY8qrkOGxqklEATc057GzirYS8gPFeKLCNQoRqcfrm0iJyGbgXOAhvzlqwbCmJ0/6L+zxRIrxeMrbjns3bC9tcjvzZj57esZ2ln2xdP6JZM7NKdOJQCTg3WwjAYdIyJ38HHRYUhckHPRu6pGg46+9Jb0dDTlEAu5Z52YeGwk6djM3JZMeHks+NQpVb0rAQh62S6tbQmBsAICkzj7qrFi51igex3tP9hK80U47gY8AHy1XwQqxUJqeEumbtH/TnLiBTrnhJjNu4lNvxOPxJOMT6WfflCdvyJP5eDfqyTyK5QhEgi7hgEMo4BAOeNvh4OR2Qzgw677p6Zk37Ck39Yy1d77dvE31mXiOIo+qbUqTfqAookZR10ZwdACCkNBE4fnMIddAIao6IiI3AX+jqn8hIs+VrVQVoiiDY3H+y09fnvyredrNOjbDDXr6vuRc4wtzEA44/k12hhtvwKEpGsxIn/lmHc7cF5zcDs2Z7tjIFmPyIAV0ZquqN6i2mBpFtBV38AgEIZmqfI1CROQyvBrETX7aghs6UWzTUzKlDI4luOOXB6bebKfcrL0mjJZocOKmHHKdKTfc9A3bS3dnvJFn3pQnj/GuYU0ixlSZQjqzNeXVKJximp5acUfPQEP9gmh6+gzwOeA+Vd0rIuvxnsheUErR9BQOOLz6lX9Vwk6NruYAABInSURBVFIZY2qeOIhqXk1PqilEtejO7MBoP1Bf+RqFqj4GPJbx+QAZr0Y1xphFTSQ9pXbOp6TSNYqiOrNbcfGau+Kp8o0tynXU06PM0J2vqleWvETGGFNtxPFfappHjYJSjHpqw/UvuRCanv5DxnYE+ABQvi52Y4ypJn6gyLePwuvMLq7pyfWvuRCanp6elvQvIvJUGcpTFHuOwhhTUXnOHlt0Z3a0hcA81ChyGgMpIq0ZS7uIvBdoLlupCqSqD6rq9ubmBVc0Y0wt80cp5tWZTbozu4hAEW6aGH6aSFX+OYqn8RrfBK/J6XUmh8kaY8wi53dm53HGZGd2EU1PkSZcPzhVvI9CVdeVrQTGGFMDRPPro6AUD9yFmyZu4hWrUYjIlar6cxG5fqb9qvqP5SmWMcZUEREEJZ86RSrdSOMU8exyMIo4AVykok1P7wJ+DmybYZ8CFiiMMYb8n6NQVRwpcqocEb+fQio315Oqfslff6xsJSghG/VkjKkImb2PYiQ+wpPHnuTl0y9zfPg4KU0RdIKcTo2XZqqeiNehnUoVPxnobOZqevr32far6l+WtjjFWSizxxpjFqfMPorB2CC3v3A7f//K30+8VKg92o4rLvFUnFFNsrYU/c/hJgIMVHT22EZ/fQ7wNuAB//M2YME9R2GMMZUxtelp76m9fObRz9A70sv71r2P6zdez5aOLUQCkYkzUj/5JM6Jh4u/dKSZgA5Uro9CVf8cQEQeB96qqoP+5y8D/1S2UhljTDWZaHpSnut9jpt/djPNoWa+977vsaVjy4ynOKlEcUNj08JNuOO6IJ6jWAbEMj7H/DRjjDF+jeJYYohP/fxTtEfbueOqO1hev3z2U5Kx4obGpoUbCYxp5Z+jAL4LPCUi9/mfrwPuLk+RjDGmyog3OPafhw/REGzgv1/537MHCYBUvLjpO9IiTbinU5WvUajqV0Tkp8A7/KSPqeqzZSuVMcZUmWHHG+r6+Us/z7rmHJ5RTsZL1vQU1CTJSgcK8CYGFJEjeLPHIiJdqnq4bCUzxpiqITQmU0gwwjXrr8ntlFI1PfnTeCQS48XnNYtcJwV8v4i8ijfH02P++qGylWrqtetF5G4R+VsR+eh8XNMYY/Iiwk96jvLz1R/O/dmIZLxEfRRNBBQSybHi85pFro8F/r/A24F9/rxP7wZ+U+hFReQuEekVkT3T0q8WkVdEZL+IfNZPvh74B1X9BPD+Qq9pjDHlI3QkU4TzmY6jVE1PkSZclESi8oEirqqnAEdEHFV9FNhaxHW/A1ydmSAiLnAb8D7gPOBGETkPWAUc8Q8rX7f+ZEnKfwljTI3KY1LAVBycnFv/ZxesI6AsiD6KMyLSADwOfF9EeoHhQi+qqo+LyNppyZcA+/33cSMi9wLXAt14weI5sgQ2EdkObAfo6uoqtGjGGJO/dHNTPvOMJ2OlqVG4QQIoiWT53pmda43iWmAE+BPgp8BrzDxRYDE6maw5gBcgOvEmHvyAiHwLeHC2k1X1dlXdqqpbOzo6Slw0Y4zJJt0SkUekKFUfhRvy+igqOIUHAKqarj2kgLtFxAFuBL5froJNu3ZOkxLapIDGmIooZHK/UgUKJ4iLMlbGpqesNQoRaRKRz4nIrSJylXhuAQ4AHy5xWXqA1RmfV/lpxhhTHfKYZrx0TU8hv4+icu/M/h7ehIC7gX8LPAp8CLhOVa8tcVl2AptEZJ2IhIAbmJyEMCf2zmxjTGVUsukp6D1HUcGmp/WqegGAiNwBHAO6VLWocVgicg9wBdAuIt3Al1T1Tr+28jDgAnep6t5irmOMMfNiojM731FPJeqjoLw1irkCxUQ3uqomRaS72CDh53XjLOk7gB2F5mt9FMaYipBCahQlanpyAgTKXKOYq+npLSIy4C+DwJb0togMlK1UBbKmJ2NM1UgmStSZ7eACCa3QG+5UtYi3fs8/q1EYYyoq787sUgSKdI2icp3ZVcVqFMaYypntrdkzUPX6KErR9CSu10dRxhpFTQUKY4ypGJHcaxTpZx5K0ZntuH6NwgJFTkRkm4jc3t/fX+miGGMWnTweukv6LwwtRdOTuLha3j6KmgoU1vRkjKmsHGsUE4GiFKOeXG+uJ+ujMMaYBS6fpqek3/RUkhqFQ1MqxagmiJdpYsCaChTW9GSMqZw8OrNL2fTkuLQmvdpE31hf8fnNdImy5Foh1vRkjKmYfCYGLGnTU4DWpNc/YYHCGGMWukqMehKXtjLXKErweiVjjDGzNj2pwv6fwXM/gJ5dMHpm8s12bgluwY5b9hqFBQpjjCmFmTqzx/rhvpvhlR1QvxTWvgMalkJ8BBIx6Lq8BNd1yt5HUVOBoiRTeNgrs40xBZlWoxgfgrvfD2/sgau+Apdsh0AJ+iTOuqxQjxDC4dToqdLnT431UVhntjGmYjI7s1XhJ38Mx3fDR74Pl99SniAxcWmXVjfMqTELFMYYs7Clm5723gcv3g9X/hmcc3X5r+u4tDphG/VkjDELm1+jSIzDI/8JVlwIl396ni7tskSC9I+X5xkyCxTGGFMK6c7sp78DA93w7i+XZlRTLpwAUXEYTYyWJ/uy5GqMMYuOgCbh17dC12Ww/or5u7TjEMECRU5sCg9jTEXt/xn0H4ZL/11+T2oXS1wiOIwlin5T9YxqKlDYqCdjTMWIwKn90LAMzr1mfq/tuEQRxpIWKIwxZuEaH/DWF3yoNJP95UNcIgijiVE0n9ex5sgChTHGlNJ5183/NR2XKA4pTRFPlX6qcQsUxhhTSp0Xz/81xeFDwQ52XL+DgFP6kVY1NYWHMcZU1OarwanA399OgBZ1aGlcXZbsLVAYY0wpfPH0/I50yuS4k1OXl4EFCmOMKYVK1CTSxPWe4SiTBd9HISLrReROEfmHSpfFGGMWJMeFVKp82ZctZ0BE7hKRXhHZMy39ahF5RUT2i8hns+WhqgdU9aZyltMYY6qaOGWtUZS76ek7wK3Ad9MJIuICtwHvAbqBnSLyAOACX5t2/sdVtbfMZTTGmOrmBKq3j0JVHxeRtdOSLwH2q+oBABG5F7hWVb8GFPw4o4hsB7YDdHV1FZqNMcZUnzIHikr0UXQCRzI+d/tpMxKRNhH5NnCRiHxutuNU9XZV3aqqWzs6OkpXWmOMWejcICSrtEZRCqp6Crg5l2NL8SpUexOqMabqOAEo04SAUJkaRQ+Q+VTIKj+taDYpoDFmUXKDkIyVLftKBIqdwCYRWSciIeAG4IFSZGzTjBtjFiWnvE1P5R4eew/wBHCOiHSLyE2qmgBuAR4GXgJ+qKp7S3E9q1EYYxYlNwBlmAwwrdyjnm6cJX0HsKPU1ytFH4UxxlQdJwjJ8gWKBf9kdj6sRmGMWZTcUFlrFDUVKKyPwhizKLmB6u2jmG9WozDGLEpO0GoUxhhjsnCtjyJn1vRkjFmUglGIj5Qt+5oKFNb0ZIxZlEKN3gN3ifI8dFdTgcIYYxalcIO3jg2VJXsLFMYYU+1CfqAYHyxL9jUVKKyPwhizKFmNInfWR2GMWZRCjd563AKFMcaYmUzUKKzpyRhjzEwm+iisRjEn66MwxixK1keRO+ujMMYsStZHYYwxJivrozDGGJNVIOxNDGg1CmOMMbMKN1gfhTHGmCxCjVajMMYYk4XVKHJjw2ONMYtWqMHmespFscNjtcTlMcaYeROMQGK8LFnXVKAwxphFyw1DYqwsWVugMMaYWhAIey8vKgMLFMYYUwsCYWt6MsYYk4VrNQpjjDHZBEJlq1EEypJrCYnIdcDvAk3Anar6SIWLZIwxC48bhmQVNj2JyF0i0isie6alXy0ir4jIfhH5bLY8VPV+Vf0EcDPwkXKW1xhjqlYV1yi+A9wKfDedICIucBvwHqAb2CkiDwAu8LVp539cVXv97S/45xljjJnOLV9ndlkDhao+LiJrpyVfAuxX1QMAInIvcK2qfg24ZnoeIiLA14GHVPWZcpbXGGOqViAMKKSS4LglzboSndmdwJGMz91+2mw+Cbwb+KCI3DzbQSKyXUR2iciuEydOlKakxhhTLTrOhfN/zwsUJbbgO7NV9ZvAN3M47nYROQZsC4VCF5e/ZMYYs4Cc935vKYNK1Ch6gNUZn1f5aUWzV6EaY0zpVSJQ7AQ2icg6EQkBNwAPlCJjmz3WGGNKr9zDY+8BngDOEZFuEblJVRPALcDDwEvAD1V1bymuZzUKY4wpvXKPerpxlvQdwI5SX09EtgHbNm7cWOqsjTFm0aqpKTysRmGMMaVXU4HCGGNM6dVUoLDObGOMKb2aChTW9GSMMaUnqrX3pmgROQEcKvD0duBkCYuzkNTyd4Pa/n61/N2gtr9ftXy3NaraMdOOmgwUxRCRXaq6tdLlKIda/m5Q29+vlr8b1Pb3q4XvVlNNT8YYY0rPAoUxxpisLFCc7fZKF6CMavm7QW1/v1r+blDb36/qv5v1URhjjMnKahTGGGOyskBhjDEmKwsUPhG5WkReEZH9IvLZSpen1ETkoIjsFpHnRGRXpctTLBG5S0R6RWRPRlqriPwvEXnVXy+pZBkLNct3+7KI9Pi/33Mi8q8qWcZCichqEXlURF4Ukb0i8mk/vep/uyzfrep/O+ujAETEBfYB78F7NetO4EZVfbGiBSshETkIbFXVanjwZ04i8tvAEPBdVX2zn/YXQJ+qft0P9ktU9U8rWc5CzPLdvgwMqep/q2TZiiUiK4AVqvqMiDQCTwPXAf+GKv/tsny3D1Plv53VKDyXAPtV9YCqxoB7gWsrXCaThao+DvRNS74WuNvfvhvvP2nVmeW71QRVPaaqz/jbg3jvpOmkBn67LN+t6lmg8HQCRzI+d1MjP3AGBR4RkadFZHulC1Mmy1T1mL99HFhWycKUwS0i8oLfNFV1TTPTicha4CLgSWrst5v23aDKfzsLFIvHO1T1rcD7gD/2mzdqlnptqrXUrvotYANwIXAM+EZli1McEWkAfgx8RlUHMvdV+283w3er+t/OAoWnB1id8XmVn1YzVLXHX/cC9+E1t9WaN/x24nR7cW+Fy1MyqvqGqiZVNQX8LVX8+4lIEO9G+n1V/Uc/uSZ+u5m+Wy38dhYoPDuBTSKyTkRCwA3AAxUuU8mISL3fuYaI1ANXAXuyn1WVHgD+yN/+I+AnFSxLSaVvor7fo0p/PxER4E7gJVX9y4xdVf/bzfbdauG3s1FPPn/I2l8BLnCXqn6lwkUqGRFZj1eLAO896T+o9u8nIvcAV+BN4fwG8CXgfuCHQBfeNPMfVtWq6xSe5btdgdd0ocBB4N9ltOlXDRF5B/BLYDeQ8pM/j9eWX9W/XZbvdiNV/ttZoDDGGJOVNT0ZY4zJygKFMcaYrCxQGGOMycoChTHGmKwsUBhjjMnKAoVZ9ESkLWNmz+PTZvr8dZmueZGI3Olv/xsRURF5d8b+6/y0D+aZ770isqnU5TWLmwUKs+ip6ilVvVBVLwS+Dfx/6c+qenmZLvt54JsZn3fjPeiZdiPwfD4Z+rMgfwv4j0WXzpgMFiiMyUJEhvz1FSLymIj8REQOiMjXReSjIvKU/56PDf5xHSLyYxHZ6S+/NUOejcAWVc0MBL8ELhGRoD9X0EbgOf/4K0Xk/ozz3yMi96XLJyLfEJHngcv8fN4tIoHy/IuYxcgChTG5ewtwM/Am4A+Azap6CXAH8En/mL/Gq5G8DfiAv2+6rZw9jYMCPwPeizflduYUMo8C54pIh//5Y8Bd/nY98KSqvkVVf+XPJ7TfL6sxJWGBwpjc7fTfOTAOvAY84qfvBtb62+8GbhWR5/Bu9k1+DSHTCuDEDPnfi9f8dANwTzrRn031e8C/FpEWvJrDQ/7uJN4kdJl6gZV5fztjZmHVU2NyN56xncr4nGLy/5IDvF1Vx7LkMwpEpieq6lMicgEwoqr7vDnmJvwd8CAwBvxIVRN++piqJqdlFfGvYUxJWI3CmNJ6hMlmKETkwhmOeQmvD2Imn8Xr6J5CVY8CR4Ev4AWNbDZThTOUmoXLAoUxpfUpYKv/NrMX8fo0plDVl4Hm9NTv0/Y9pKqPzpL394EjqvrSbBcXkWXAqKoeL6z4xpzNZo81pgJE5E+AQVWdqbN7tnNuBZ5V1TvnyHcg2zHG5MtqFMZUxreY2ueRlYg8DWwB/scch54B7i6iXMacxWoUxhhjsrIahTHGmKwsUBhjjMnKAoUxxpisLFAYY4zJygKFMcaYrP437BUn7ZQo6w4AAAAASUVORK5CYII=\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "plt.plot(sliced_df['t'], sliced_df['radius'], label='radius')\n",
+    "plt.plot(sliced_df['t'], sliced_df['omega'], label='Omega')\n",
+    "plt.plot(sliced_df['t'], sliced_df['v_eq'], label='Equatorial velocity')\n",
+    "plt.xlabel('Time (Myr)')\n",
+    "plt.ylabel('Radius (Rsol)')\n",
+    "plt.yscale('log')\n",
+    "plt.savefig(os.path.join(result_dir, 's'))\n",
+    "plt.show()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.6.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/david_calculations/single_star_full.ipynb b/david_calculations/single_star_full.ipynb
new file mode 100644
index 000000000..30c50c50e
--- /dev/null
+++ b/david_calculations/single_star_full.ipynb
@@ -0,0 +1,133 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import os, sys, time\n",
+    "\n",
+    "import matplotlib.pyplot as plt\n",
+    "from collections import defaultdict\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "\n",
+    "# sys.path.append('../')\n",
+    "import binary_c\n",
+    "\n",
+    "\n",
+    "from binaryc_python_utils.functions import create_arg_string, parse_output, run_system\n",
+    "\n",
+    "result_dir = '../david_results/'"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Execute command and parse args"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Took 0.20s to run single system\n",
+      "The following keys are present in the results:\n",
+      "dict_keys(['t', 'mass', 'zams_mass', 'stellar_type', 'prev_stellar_type', 'metallicity', 'probability', 'dM_in_timestep', 'mdot', 'core_mass', 'core_radius', 'luminosity', 'radius', 'dt', 'dtm', 'Teff', 'omega', 'vwind', 'drdt', 'tm', 'tn', 'tkh', 'angular_momentum', 'he_core_mass', 'CO_core_mass', 'GB_core_mass', 'v_eq', 'v_eq_ratio'])\n"
+     ]
+    }
+   ],
+   "source": [
+    "start = time.time()\n",
+    "output = run_system(M_1=10, M_2=20, separation=0, orbital_period=100000000000)\n",
+    "result = parse_output(output, 'DAVID_SINGLE_ANALYSIS')\n",
+    "stop = time.time()\n",
+    "print(\"Took {:.2f}s to run single system\".format(stop-start))\n",
+    "print(\"The following keys are present in the results:\\n{}\".format(result.keys()))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Cast data into pandas framework and do analysis"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#### Now do whatever you want with it: \n",
+    "\n",
+    "# Cast the data into a dataframe. \n",
+    "df = pd.DataFrame.from_dict(result, dtype=np.float64)\n",
+    "\n",
+    "# Get last change moment\n",
+    "last_st = df['stellar_type'].unique()[-1]\n",
+    "last_stellar_type_change_time_1 = df[df.stellar_type==last_st]['t'].iloc[0]\n",
+    "\n",
+    "# slice to get that last time\n",
+    "sliced_df = df[df.t < last_stellar_type_change_time_1] # Cut off late parts of evolution"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZgd5X3g+++v6qy9qxdtLbV2gcHIYGQwxI4J18Z4ggzBG8STZGzGGu4NXjJ37sR2PLYz9/EymXFu4oHYlwAx9tgQOw4YJsIwjjHYMQaJVRKLEEJLtyRaUku9d5/tN39Une7Tre7TZ+3T5/Tv8zz1VJ23qt56j86j+vW71FuiqhhjjDGzcSpdAGOMMQubBQpjjDFZWaAwxhiTlQUKY4wxWVmgMMYYk1Wg0gUoh/b2dl27dm2li2GMMVXj6aefPqmqHTPtq8lAsXbtWnbt2lXpYhhjTNUQkUOz7bOmJ2OMMVlZoDDGGJNVTQUKEdkmIrf39/dXuijGGFMzaipQqOqDqrq9ubm50kUxxpiaUVOBwhhjTOlZoDDGGJNVTQUK66MwxpjSq6nnKFT1QeDBrVu3fqKgDPb/M3TvgmAEAhEIhCEQ9dZBfx2ITC4zHSdS2i9ljDEVVlOBomg//RycfKW4PNxwRgDJCCQTgSYj4OQakKYcO0uerv2UxpjysLvLFApv2gbXfRsS45AYm7rEp31OjEN89Oxj4xn7E6NTjxsbmOX80eKKLm5uwScY9T/762DdtPVsaRmfAxGrORmziFigmM4JQLjBW+aTKiRjpQtI088fG4BEr3/eKMRH/GPHCivvTMHkrOCTLehM2xeq95Zgeh21YGTMAmGBYqEQ8f/qD0NkHp8DSaW8IBPPXEamrTMDy/Rjph0fG4bhk2fnkRzPs2CSETzqINTgf04HlYaMANMwmZ4ONNMXC0DGFMwCxWLnOJM303JKJSdrMNODSGwE4sNekMlc4iMQG/I/+9tjAzB43E/3A1M+zXbiQKgRwplLw+R2aHp6k5/eMO2cRi+oG7MI1FSgEJFtwLaNGzdWuihmOsctX5NeKjlZmzkr2GR+9oPO+ODUZWwABo5OTUPnvq4bygggTVMDTrjRqxlGmiHSMrmOtmSkN1uwMVWhpgJF0cNjTXVy3MmbcymkUl7gSQeN2LTAMj4E4wMZ+4cmt4d64dRr3v7RM5CKZ79WIDoZNKYHkUjL2fuiS6CuDaKtXnObMfOgpgKFMSXhOBm1nxWF56PqNbWN9XtBY6w/YznjL/1T9w/1wslXJ/dpavb8A1Goa/WCRl1rxnbbtO0lk9vhRuujMXmzQGFMuYhMjvBqXJ7/+ap+v0w6kJyB0dMwcgpG+mC0z1unt4/v9rdPM2vTmRP0aiX17d7SsAzql0JDh79eCvUdk2s3WNQ/gakNFiiMWahEJpvUmlflfl4q6QWXkT4vqEwElGnbQ73QvdNbx0dmzivaOi14LIXGZdDUCU0rvaVxpff8jqlZFiiMqTWOO9kURY4DO8aHYLgXhk74a39Jbw+fgJ5n/KAyfPb5dW1ewEgHj8xA0tTpBTrrU6laFiiMMZN9Mq3r5z52bAAGj3kjxSaWHj+tB3p2eTWW6eqXwpI1sGQttKzxtlv8z02dNg3NAma/jDEmP5Emb+k4Z/Zj4mMZwaQHzhyC04e89ZEnYc8/giYnjxfXq3UsWQvtm6B9s7du2+QFEaemJrquOgs+UIjIm4BPA+3AP6vqtypcJGPMXIIRaF3nLTNJxr0AcvrgZAA5fRD6XocXfugNL57Iqw7aNvrBww8gy86H1g1WC5knFflXFpG7gGuAXlV9c0b61cBfAy5wh6p+XVVfAm4WEQf4LmCBwphq5wa92sOStWfvU/WHCe/zl1e9dfdTsOfHTIzoCkRg6Ztg2Zth+QXesuz8+Z0CJ63vdW89W2CscpUKx98BbsW78QMgIi5wG/AeoBvYKSIPqOqLIvJ+4P8EvleBshpj5pOIN7KqcRmse+fUfbEROLUf3tjrDQd+Yze8/E/wbMatoWUNdL4VVr0NOrfCireUf1TWNy/01l+uzZemVSRQqOrjIrJ2WvIlwH5VPQAgIvcC1wIvquoDwAMi8k/AD2bKU0S2A9sBurq6ylRyY0xFhepgxRZv4UYvTdXrDzm+B46/4C3du2Dvfd5+J+DVNjq3wupLYM1vQXNnxb5CNVpIDXydwJGMz93ApSJyBXA9EAZ2zHayqt4uIseAbaFQ6OJyFtQYs4CITA7F3XzVZPrgcS9g9Ozy1s/9AHb+rbevdT2sfae/vAOaingCfxFYSIFiRqr6C+AXOR5rcz0ZYzyNy+FN13gLeA8ivrEHDv4KXv8l7L0fnrnb29e2ETZcCZve6wUOe4BwioUUKHqA1RmfV/lpObPZY40xs3Jcr79ixVvgsj/2AsfxF/zA8Tg88z146nZvDq3174JNV8Hm9879VHxshgcQK+DZw6d59vAZPvr2LsIBt6R5L6RAsRPYJCLr8ALEDcDv55OB1SiMMTlzXFh5kbdc/knv3SgHfwX7HoZXH4Z9P4V/wuvbOP86OO86aFl9dj7P3zPvRZ/JY/tO8Fc/e5U/uGxNyfOu1PDYe4ArgHYR6Qa+pKp3isgtwMN4w2PvUtW9eeZrNQpjTGGCUdj0Hm/R/+oNyX1lh9dE9cgXvGXV27yAcf51Xk0jlYLffNs7f9mbs+dfZqeHYzRHgwTd0j+cWKlRTzfOkr6DLB3WOeRrNQpjTPFEvCfPO86Bd/wJ9B3wAsbe++CRP/OWVZd4o69Oveqdk21K+HlwajhGa32oLHkvpKYnY4xZmFrXwzv/vbeces0LGHvvh513eJMhrrzQCyYVdHokxpK68kwLX1OBwpqejDFl17YBfvs/eMup17ynzP/XF73O8Qo63j/GOctL9JbHaWpqpi1VfVBVtzc3V+ARfmPM4tO2AVq6vEkNK9j0pKocPTPGyuZoWfKvqUBhjDEVIc7U2XDn2emROKPxJCtbLFDMSUS2icjt/f21Od+KMWaBcipbozh6ZhTAAkUurOnJGFMR4nhzTlVIjx8oVi2xQGGMMQuTSEU7s3tOW40iZ9b0ZIypCHE5PTLGV3e8VJHLHz0zSiTolG14bE0FCmt6MsZUhDgkEkluf/wAvzkww/vCy+xo/ygrW6KISFnyr6lAYYwxFSEODl5n9qfvfXaiz2C+9JwZo7NMzU5ggcIYY4rnuDgo79rcwUgsyR/d9RS9A2Pzdvme06MWKIwxZkETBwflzZ1N3PGHWzl6ZpRrb/sXnnq9r+yXHosnOTk0XraObKixQGGd2caYivCbnlwRLl3fxo9uvgzXET78/z/Bzd97msf2nWAsXp5RUcf7vZpLOQNFTc31ZLPHGmMqQfEDheP97X3+ymYe+ZPf5lu/eI3vPnGIn+49TtAVOluitNaHCDgOdWGXr/zeBUU3GaX7Q8rZ9FRTgaJYn4/EeGTkGQI/eDuCN3ogvZ5cycTIgunHZI44mJ425ViZ/ZiZ8pp+7Gx5z3q96WXLckxmXnOVbcZ/m3zKNi1vEcERZ6Icjjg4OJPpGcfMmp75OeOY6enTr5OZL8KUY9LnTllnljF9nH9Oegk4AW8t3tp1XFxxJ9fTts86x3HmPM4Vt2wjXUzu/rDvV8RXLuEdGW00daEA//dV5/B/XbGR3xw4xc6DfRzuG6F/NM5YPMkvXjnBL17p5aOXFveiIQsU82y3q4yjfGjj701JV7wnLlV1ynbmvpnSph87U9pMeWfL76zrTqzyLNtMZc3Ia9bvPcN1czlmpuuq6lmfU6S8taYm0zQ1e3rmZ/zP04/LSFd0yjEpUqBMHFON0gEjHVAccQg6QQJOgKATnHnbnX1/5jHZ8gm6QSJuhJAbIuJGCAfCs34OOLV9q3ku3geREO9yzm7Nj4ZcfufcpfzOuUsn0lIp5fwvPcyBE8W/RrXn9CgisKw5XHRes6ntXy9PovBedwl/esmfVroopkJyCTITgYaz09NLUpPekkqS0hQJTZBKeemJVGIyTVMkU8nJ4zPP8Y/LTE9vZ+6f6bhEKkFCE8STceKpOIlUgnhqcjuWjDEcG55x3/Tjpv8BU4iABAgHwoTdySUSiExsRwNR6oJ11AfrqQ/UT2zXBeuoC/jpwXrqAnWT+/z0hVSjyvXlco4jrGuv57UTQ0Vf80jfCMubIiV/T3YmCxTGZBARXCnff7hqlEwlZw0i48lxxpJjxJIxxhJjOX2eWBLevvHkOAMjAwzHhxmJjzCSGGE0kdtzCK64NIWaaAo30RxqpjHcSHOomaZQE83hZprDk9utkVbaom20RdqIBCJl+bdyZ6hRzGZ9Rz3Pd58p+pqH+kboaq0rOp9saipQ2IuLjCk91/GatOZTMpVkNDHqBY/ECCPxkYntiXVsmIHYgLeMD9Af66d/rJ/DA4fpH+9nMDY4a22oMdjoBQ0/cLRH22mPtrOsfhkr6lewsmElS+uWEnTymxIj4OReu1nfXs+O3ceIJVKEAoUPQD10aoT/I6NZqxxqKlDYqCdjaoPruDSEGmgINRScR0pTDMYGJ4JI31gfJ0dPcmr0lLce89b7Tu/jiaNPMBgfnHK+Iw5L65aysn4ly+uX09nQyerG1axpWsOapjW0RlrPavbKI06wtr2elMLhvhE2Li3sew6NJzg5NE5Xm9UojDEmb444E81Pq1k95/FjiTGODx/n6PBRjg0d4+jwUe/z0FGeP/E8jxx8hIQmJo5vDDbS1dRFR7RjIi2fisHa9noADp4cLjhQHD41AsAaCxTGGFN+kUCEtc1rWdu8dsb9iVSCY0PHODhwkEMDhzg0cIiDAwc5PnJ84phAHjWKdW1+oDhV+Minw33euWta6wvOIxcWKIwxJgcBJ8DqptWsblrNO3nnlH1//sMP8PDQy7iS+wixJfUhmqPBogLFQb9GUe6mp5qawsMYYyrBxSEp+dUowGt+ev1k4YFif+8QHY1hmqPleQ9FmgUKY4wpkiMOCSTn5yjSNnY08OobhT9L8crxQc5Z1ljw+bmyQGGMMUWarFHk93Dim1Y00js4zqmh8byvmUwpr/YOcs5yCxSIyHUi8rci8vciclWly2OMMdM5OCRFcCS/aWDOXd4EeDWDfB3uG2EsnqrdGoWI3CUivSKyZ1r61SLyiojsF5HPAqjq/ar6CeBm4COVKK8xxmTj+LdSR/KbSvzcFd5N/qUCAkU6uNRyjeI7wNWZCSLiArcB7wPOA24UkfMyDvmCv98YYxYUR7xbqZJfoGhvCNPeEGZvT/7v0NnT04/rCJtrtUahqo8D01/9dAmwX1UPqGoMuBe4Vjz/BXhIVZ+ZLU8R2S4iu0Rk14kTJ8pXeGOMmSZdoxDieZ97UVcLzx7Jf86nZ4+c5k0rGomGyj+9ykLqo+gEjmR87vbTPgm8G/igiNw828mqeruqblXVrR0dHbMdZowxJefg3awlzxoFwMVrlvD6yeG8OrSTKeX5I/1ctHpJ3tcrxEIKFDNS1W+q6sWqerOqfjvbsfYqVGNMJQjppqfEHEee7eI13s3+mcO51yr2vTHI0HiCC1e35H29QiykQNEDUyZkWeWn5UxVH1TV7c3NzSUtmDHGZDPZ9JR/jeKCzmZCAYcnXjuV8zm/fNVrXr98Y1ve1yvEQgoUO4FNIrJORELADcAD+WRgNQpjTCU4RdQoIkGXy9a38egrvTmf8/i+k2xe1sCK5vK9/jRTpYbH3gM8AZwjIt0icpOqJoBbgIeBl4AfqurefPK1GoUxphKK6aMAuPLcpbx+cjinN971j8Z56vU+3rV5/vpiKzXq6UZVXaGqQVVdpap3+uk7VHWzqm5Q1a/km6/VKIwxlSB4kzylCqhRAFx1/jIcgfufnbu1fcfuY8SSKa7ZsrKgaxViITU9Fc1qFMaYSii2RrGiOcq7Nnfww11HiCdnf7pbVbn3qcNsXNrAllXzd5+rqUBhNQpjTGUU9sBdpj+8bC1vDIxzz1OHZz3m0Vd6eb67n4//1rqz3q5XTjUVKKxGYYypBEf9QKGFB4orzungsvVt/LeHX+HADH0VJwbH+cJ9e1jXXs8HL15V8HUKYS8uMsaYIonf9FTIqKeJPET4iw9u4f23PcYH7vwRH7y4i4tXnIuIS/fpUe761eucHolz7/aLCeXzztUSqKlAISLbgG0bN26sdFGMMYuIlKDpCeDI6HM0bPw6Z2J93HsU7jkSIX7mbcT63sl5S1dx6+9fxFvm6SG7TNb0ZIwxRZKJpqfCaxR7Tu7hkz//JB31bfzXS7/If1pzI5d3XESk/de0nfuXfPQ9R7ioa36m7JiupmoUxhhTGV7TU6rAPoqUpvjPT/xnlkSWcOeGf82SH98C8WE+DHSv/22+2tLE1576GpFAhOs3XV/CcuempmoUxhhTCcXM9QTws0M/46W+l/j0+R9jyQOfgbb18Af3w+/8GauOv8g3n36IyxrX89Unv8qBMwdKWfSc1FSgsOGxxpjKKK6P4kf7fkRnQye/2/0KxAbhg38HG34H3vUf4VPPElhzOV996TeEJMA3nv5GKQuek5oKFNZHYYypCPWanpIFND0dGzrGk8ee5Nr123Ce+z6c+7vQvmnygEgTfOhu2p0wH0/V8Xj347x46sVSlTwncwYKEblMRG4TkRdE5ISIHBaRHSLyxyJid2RjjKHwzuyfH/k5ivK7kZUwchLe/IGzD6pvg8tv4cOvP0PECfGjfT8qtsB5yRooROQh4N/iTdR3NbAC7zWlXwAiwE9E5P3lLqQxxixkqVT6Vjr79Buz+Zeef6GrsYuuw7vADcHGd8984IUfpUmFq6KdPPT6Q8SSscILnKe5ahR/oKo3qeoDqnpUVROqOqSqz6jqN1T1CuDX81DOnFgfhTGmEtSfFFDQvM4bT46z8/hOLl95ORx5EjovhvAs78Bu7oSuy7iq7w2G48PsOr6r2GLnLGugUNWTc2WQyzHzxfoojDGVkCI971J+geL53ucZS47xW8svgWPPw6qt2U/YcCWXHnuFiBvmF92/KKishZir6WlQRAYylsHM9XwV0hhjFrLUxAN3+TU9PX/ieQAuSgUhGYNVl2Q/YcOVRFS5tGENvz46f405WR+4U9VZ6kDGGGPSUn5FIt9A8cLJF1jbtJbm3pe9hLlqFMsvADfMW5Mujw0d4tToKdqi5X8das7DY0XkLSJyi79sKWehjDGmmkzUKPJoelJVdp/YzQXtF8CJlyHSAo0rsp/kBmH5BVw04L1f+7kTzxVc5nzkFChE5NPA94Gl/vJ9EflkOQtmjDHVIl2jSKVyf47i+PBxTo2d4oKOC+DkPug4B3J5x0TnWznv6EsEnSDP9S6gQAHcBFyqql9U1S8Cbwc+Ub5iGWNM9UilvBt8PjWK1wdeB2BTyyYvULRvzu3E5VsIx4fZ3LiGl/tezrushcg1UAhMeTY96actKDY81hhTCan0A3ep3Psojg8fB2CFG4XhE16NIhd+QNkcamXf6X35FbRAuQaKvwOeFJEvi8ifA78B7ixfsQpjw2ONMZWQUvHXuTc9jcRHAGgY8AIG7bkGCm96j80aoG+sj5Oj5X9CIadpxlX1L0XkF8A7/KSPqeqzZSuVMcZUkaQKiDddeK7iqTgAwTPdXkLbhtxOrGuFunY2jXmBZt/pfbRH2/Mqb75y7czeAOxV1W8Cu4F3isj8v2bJGGMWoHSNIp8+iolAMejXKJo6c79g+yY2978BwKunX839vALl2vT0YyApIhuBbwOrgR+UrVTGGFNFkv7wWDT/QBHoPwL1SyEYyf2CbRtY0neIlnALhwYO5VPUguQaKFLqTYt4PXCrqv4/eBMEGmPMopeaiA95ND0l4wSdIDLQA82r8rtgcxcMHaerYRWHBw/nd24Bcg0UcRG5EfhD4H/6acHyFMkYY6pLYuLJ7PxqFEEnCGeOQMvq/C7oB5bVkTaODBzJ79wC5BooPgZcBnxFVV8XkXXA98pXrEkisl5E7hSRf5iP6xljTK76hmMMjMUnphkvKFD0d0NzYYGiy63n2PCxsk85nuuopxeBT2V8fl1EflXoRUXkLuAaoFdV35yRfjXw13hvKr9DVb+uqgeAmyxQGGMq6cxIjD09A+zu6Wd3zxle6O6n+/Qoq1ujXBZNQiS/uZ68QOFCYrTgQLEaF0XpHupmffP6/PLIQ9ZAISIu8GGgE/ipqu4RkWuAzwNR4KICr/sd4Fbgu9OudRvwHqAb2CkiD/hByhhj5s2poXF29/Sz9+gAe3r62XO0nyN9oxP7u1rreMvqFq7ZspJvP/Yaze4wNOU56ikZJ5h+bjnvPgo/UIx7NYmewZ7KBQq8h+pWA08B3xSRo8BW4LOqen+hF1XVx0Vk7bTkS4D9fg0CEbkXuBawQGGMKYvh8QSvnRhi3xtDvNo7yKtvDPHSsQGO9Y9NHLOmrY4tnS38/iVrOH9lE1tWNdNSF5rYf2YkxgvP+FNp5FujSMeV5jyGxgIEwtCwjGWjZwDoHenN7/w8zRUotgJbVDUlIhHgOLBBVU+VoSydQGavTDdwqYi0AV8BLhKRz6nq12Y6WUS2A9sBurq6ylA8Y0w1UlV6B8c5eHKYQ30jvNY7xL43Bnm1d4ju05O1hJDrsL6jnkvWtXJBZzPnr2zmvJVNNEezj9v57PvO5cZnfuZdK8/nKILp4xuW5//FmlfRPnACQXhj5I38z8/DXIEipn6jm6qOiciBMgWJWfnXuzmH424XkWPAtlAodHH5S2aMWShiiRTH+kc53DfCwVMjHD417K9HONQ3zFh88i/9UMBhQ0cDb+1awke2rmbTskY2LWtgTWsdATfnNy9MaKkL8Vc3vJUPPvN3+Xdmp8fV1hfwZHXDcoKnD9LW3lbxGsW5IvKCvy3ABv+zAKqqpXwvRQ9eM1faKj8tZ6r6IPDg1q1bbWZbY2pEKqWcHB7n2Jkxjp4ZpefMKMf6ve2j/vrk0PiUZ93CAYc1bXWsaavnnZvaWdNez5rWOta01dHZEi0oIGSzekkDUEBntiahrs17z0S+GpZC91Ms7bqY4yPH8z8/D3MFijeV9epT7QQ2+UNve4AbgN/PJwMR2QZs27hxYxmKZ4wppVgixcmhcU4MjtM76K297TFve2ic3gEvLZacegOOBl1WtkRY2RLl3HOWssLf7mqtY21bPUsbwzjO/E1wLU7+r0JNJBMEkgnvqexCNCyD4ZMsiy6lezivv6nzNlegOKxz1KVEROY6ZoZz7gGuANpFpBv4kqreKSK3AA/jDY+9S1X35pOv1SiMqZyxeJLTIzH6hmOcHo7TNxLj9LD/eWRynQ4MZ0biM+bTWh+ioyFMR2OYS9fV09EUprMlyormqBccmqO01AWRXF7yM0+8QZv59VGkSCHJODTkOeIpraEDUJaFmni695nC8sjRXIHiURH5MfATVZ14TlxEQngzyf4R8CjecNecqeqNs6TvAHbkk1cmq1EYUxxVZTiWZGA0Tv9onIHROANjicnPY3HOjMSn3PhPD8fpG44xGp99iu2WuiCtdSGW1IdY397Apeva6GgMs7TRCwjppb0hTLDEzULzQST/uZ5UFScZh6ZlhV20wTtvmRNhIDbAaGKUaCBaWF5zmCtQXA18HLjHbxI6A0Tw/uJ/BPirhTTduNUozGIXS6QYHk8w5C/D4wkGxxIMjM1+4x8YTWRsxzPmLZpZYzhAa0OIJXUhljZGOGdZE631QVrqQrTWe+mt9SFa64MsqQvRHA2WvE9goRHyf2d2SlO4yZjX11AIP1As9W/jvSO9rGlaU1hec8gaKFR1DPgb4G9EJAi0A6OqeqYspTFmERpPJBkaSzA8npxyg59+sx8eTzAcm9z29ienHBtLzN1GHg26NEUDNEWCNEWDtDeE2NBRT1M06KcFaJ7Y9tbNUS+9IRyo+Zt+ISb7KPJohU8lcTRVRKDwzmuIe898DMeHC8snBzlN4QGgqnHgWNlKUgLW9GTKRVUZi6cYiSUYiSUZjXs36NFYkpFYkpF4kpHxqftGYklvf8a+kXiS0ZgXFEZi3s09nszt5lIXcqkPB2gMB6gPB6gPu3S2RGkIuzREAtP2TW57N3tv3RgJEA64Zf7XWoQk/87sVHLc2yi0M9s/LzI+BMB4Or8yyDlQVANrelqcVJVYMsVYLMVYwrs5T6zjKcbiScbi3g18LJ7y194y4t/oR/0A4C2Z2/6+eDKf5mdcR6gLuf4SIBr0tpujQVY0RSZu+vXhAI2RAPUhl4ZIkIawl96QXvwAUB8K4M7jKB6Tn3QfRT5NT5pK4KAQLfAdcKE6CNYRjXk1idHE6BwnFK6mAoVZOJIpJZZIMZ5IMp5ITb1R+zfysYkb+uT+8Vlu6JmfveOmfs5v3J1HBOqCLtGQ99d5+mZeFwrQ3hCmLuTtqwu51GdsR/0AUB8KTGynz0vvD7nOghqVY8qrkOGxqklEATc057GzirYS8gPFeKLCNQoRqcfrm0iJyGbgXOAhvzlqwbCmJ0/6L+zxRIrxeMrbjns3bC9tcjvzZj57esZ2ln2xdP6JZM7NKdOJQCTg3WwjAYdIyJ38HHRYUhckHPRu6pGg46+9Jb0dDTlEAu5Z52YeGwk6djM3JZMeHks+NQpVb0rAQh62S6tbQmBsAICkzj7qrFi51igex3tP9hK80U47gY8AHy1XwQqxUJqeEumbtH/TnLiBTrnhJjNu4lNvxOPxJOMT6WfflCdvyJP5eDfqyTyK5QhEgi7hgEMo4BAOeNvh4OR2Qzgw677p6Zk37Ck39Yy1d77dvE31mXiOIo+qbUqTfqAookZR10ZwdACCkNBE4fnMIddAIao6IiI3AX+jqn8hIs+VrVQVoiiDY3H+y09fnvyredrNOjbDDXr6vuRc4wtzEA44/k12hhtvwKEpGsxIn/lmHc7cF5zcDs2Z7tjIFmPyIAV0ZquqN6i2mBpFtBV38AgEIZmqfI1CROQyvBrETX7aghs6UWzTUzKlDI4luOOXB6bebKfcrL0mjJZocOKmHHKdKTfc9A3bS3dnvJFn3pQnj/GuYU0ixlSZQjqzNeXVKJximp5acUfPQEP9gmh6+gzwOeA+Vd0rIuvxnsheUErR9BQOOLz6lX9Vwk6NruYAABInSURBVFIZY2qeOIhqXk1PqilEtejO7MBoP1Bf+RqFqj4GPJbx+QAZr0Y1xphFTSQ9pXbOp6TSNYqiOrNbcfGau+Kp8o0tynXU06PM0J2vqleWvETGGFNtxPFfappHjYJSjHpqw/UvuRCanv5DxnYE+ABQvi52Y4ypJn6gyLePwuvMLq7pyfWvuRCanp6elvQvIvJUGcpTFHuOwhhTUXnOHlt0Z3a0hcA81ChyGgMpIq0ZS7uIvBdoLlupCqSqD6rq9ubmBVc0Y0wt80cp5tWZTbozu4hAEW6aGH6aSFX+OYqn8RrfBK/J6XUmh8kaY8wi53dm53HGZGd2EU1PkSZcPzhVvI9CVdeVrQTGGFMDRPPro6AUD9yFmyZu4hWrUYjIlar6cxG5fqb9qvqP5SmWMcZUEREEJZ86RSrdSOMU8exyMIo4AVykok1P7wJ+DmybYZ8CFiiMMYb8n6NQVRwpcqocEb+fQio315Oqfslff6xsJSghG/VkjKkImb2PYiQ+wpPHnuTl0y9zfPg4KU0RdIKcTo2XZqqeiNehnUoVPxnobOZqevr32far6l+WtjjFWSizxxpjFqfMPorB2CC3v3A7f//K30+8VKg92o4rLvFUnFFNsrYU/c/hJgIMVHT22EZ/fQ7wNuAB//M2YME9R2GMMZUxtelp76m9fObRz9A70sv71r2P6zdez5aOLUQCkYkzUj/5JM6Jh4u/dKSZgA5Uro9CVf8cQEQeB96qqoP+5y8D/1S2UhljTDWZaHpSnut9jpt/djPNoWa+977vsaVjy4ynOKlEcUNj08JNuOO6IJ6jWAbEMj7H/DRjjDF+jeJYYohP/fxTtEfbueOqO1hev3z2U5Kx4obGpoUbCYxp5Z+jAL4LPCUi9/mfrwPuLk+RjDGmyog3OPafhw/REGzgv1/537MHCYBUvLjpO9IiTbinU5WvUajqV0Tkp8A7/KSPqeqzZSuVMcZUmWHHG+r6+Us/z7rmHJ5RTsZL1vQU1CTJSgcK8CYGFJEjeLPHIiJdqnq4bCUzxpiqITQmU0gwwjXrr8ntlFI1PfnTeCQS48XnNYtcJwV8v4i8ijfH02P++qGylWrqtetF5G4R+VsR+eh8XNMYY/Iiwk96jvLz1R/O/dmIZLxEfRRNBBQSybHi85pFro8F/r/A24F9/rxP7wZ+U+hFReQuEekVkT3T0q8WkVdEZL+IfNZPvh74B1X9BPD+Qq9pjDHlI3QkU4TzmY6jVE1PkSZclESi8oEirqqnAEdEHFV9FNhaxHW/A1ydmSAiLnAb8D7gPOBGETkPWAUc8Q8rX7f+ZEnKfwljTI3KY1LAVBycnFv/ZxesI6AsiD6KMyLSADwOfF9EeoHhQi+qqo+LyNppyZcA+/33cSMi9wLXAt14weI5sgQ2EdkObAfo6uoqtGjGGJO/dHNTPvOMJ2OlqVG4QQIoiWT53pmda43iWmAE+BPgp8BrzDxRYDE6maw5gBcgOvEmHvyAiHwLeHC2k1X1dlXdqqpbOzo6Slw0Y4zJJt0SkUekKFUfhRvy+igqOIUHAKqarj2kgLtFxAFuBL5froJNu3ZOkxLapIDGmIooZHK/UgUKJ4iLMlbGpqesNQoRaRKRz4nIrSJylXhuAQ4AHy5xWXqA1RmfV/lpxhhTHfKYZrx0TU8hv4+icu/M/h7ehIC7gX8LPAp8CLhOVa8tcVl2AptEZJ2IhIAbmJyEMCf2zmxjTGVUsukp6D1HUcGmp/WqegGAiNwBHAO6VLWocVgicg9wBdAuIt3Al1T1Tr+28jDgAnep6t5irmOMMfNiojM731FPJeqjoLw1irkCxUQ3uqomRaS72CDh53XjLOk7gB2F5mt9FMaYipBCahQlanpyAgTKXKOYq+npLSIy4C+DwJb0togMlK1UBbKmJ2NM1UgmStSZ7eACCa3QG+5UtYi3fs8/q1EYYyoq787sUgSKdI2icp3ZVcVqFMaYypntrdkzUPX6KErR9CSu10dRxhpFTQUKY4ypGJHcaxTpZx5K0ZntuH6NwgJFTkRkm4jc3t/fX+miGGMWnTweukv6LwwtRdOTuLha3j6KmgoU1vRkjKmsHGsUE4GiFKOeXG+uJ+ujMMaYBS6fpqek3/RUkhqFQ1MqxagmiJdpYsCaChTW9GSMqZw8OrNL2fTkuLQmvdpE31hf8fnNdImy5Foh1vRkjKmYfCYGLGnTU4DWpNc/YYHCGGMWukqMehKXtjLXKErweiVjjDGzNj2pwv6fwXM/gJ5dMHpm8s12bgluwY5b9hqFBQpjjCmFmTqzx/rhvpvhlR1QvxTWvgMalkJ8BBIx6Lq8BNd1yt5HUVOBoiRTeNgrs40xBZlWoxgfgrvfD2/sgau+Apdsh0AJ+iTOuqxQjxDC4dToqdLnT431UVhntjGmYjI7s1XhJ38Mx3fDR74Pl99SniAxcWmXVjfMqTELFMYYs7Clm5723gcv3g9X/hmcc3X5r+u4tDphG/VkjDELm1+jSIzDI/8JVlwIl396ni7tskSC9I+X5xkyCxTGGFMK6c7sp78DA93w7i+XZlRTLpwAUXEYTYyWJ/uy5GqMMYuOgCbh17dC12Ww/or5u7TjEMECRU5sCg9jTEXt/xn0H4ZL/11+T2oXS1wiOIwlin5T9YxqKlDYqCdjTMWIwKn90LAMzr1mfq/tuEQRxpIWKIwxZuEaH/DWF3yoNJP95UNcIgijiVE0n9ex5sgChTHGlNJ5183/NR2XKA4pTRFPlX6qcQsUxhhTSp0Xz/81xeFDwQ52XL+DgFP6kVY1NYWHMcZU1OarwanA399OgBZ1aGlcXZbsLVAYY0wpfPH0/I50yuS4k1OXl4EFCmOMKYVK1CTSxPWe4SiTBd9HISLrReROEfmHSpfFGGMWJMeFVKp82ZctZ0BE7hKRXhHZMy39ahF5RUT2i8hns+WhqgdU9aZyltMYY6qaOGWtUZS76ek7wK3Ad9MJIuICtwHvAbqBnSLyAOACX5t2/sdVtbfMZTTGmOrmBKq3j0JVHxeRtdOSLwH2q+oBABG5F7hWVb8GFPw4o4hsB7YDdHV1FZqNMcZUnzIHikr0UXQCRzI+d/tpMxKRNhH5NnCRiHxutuNU9XZV3aqqWzs6OkpXWmOMWejcICSrtEZRCqp6Crg5l2NL8SpUexOqMabqOAEo04SAUJkaRQ+Q+VTIKj+taDYpoDFmUXKDkIyVLftKBIqdwCYRWSciIeAG4IFSZGzTjBtjFiWnvE1P5R4eew/wBHCOiHSLyE2qmgBuAR4GXgJ+qKp7S3E9q1EYYxYlNwBlmAwwrdyjnm6cJX0HsKPU1ytFH4UxxlQdJwjJ8gWKBf9kdj6sRmGMWZTcUFlrFDUVKKyPwhizKLmB6u2jmG9WozDGLEpO0GoUxhhjsnCtjyJn1vRkjFmUglGIj5Qt+5oKFNb0ZIxZlEKN3gN3ifI8dFdTgcIYYxalcIO3jg2VJXsLFMYYU+1CfqAYHyxL9jUVKKyPwhizKFmNInfWR2GMWZRCjd563AKFMcaYmUzUKKzpyRhjzEwm+iisRjEn66MwxixK1keRO+ujMMYsStZHYYwxJivrozDGGJNVIOxNDGg1CmOMMbMKN1gfhTHGmCxCjVajMMYYk4XVKHJjw2ONMYtWqMHmespFscNjtcTlMcaYeROMQGK8LFnXVKAwxphFyw1DYqwsWVugMMaYWhAIey8vKgMLFMYYUwsCYWt6MsYYk4VrNQpjjDHZBEJlq1EEypJrCYnIdcDvAk3Anar6SIWLZIwxC48bhmQVNj2JyF0i0isie6alXy0ir4jIfhH5bLY8VPV+Vf0EcDPwkXKW1xhjqlYV1yi+A9wKfDedICIucBvwHqAb2CkiDwAu8LVp539cVXv97S/45xljjJnOLV9ndlkDhao+LiJrpyVfAuxX1QMAInIvcK2qfg24ZnoeIiLA14GHVPWZcpbXGGOqViAMKKSS4LglzboSndmdwJGMz91+2mw+Cbwb+KCI3DzbQSKyXUR2iciuEydOlKakxhhTLTrOhfN/zwsUJbbgO7NV9ZvAN3M47nYROQZsC4VCF5e/ZMYYs4Cc935vKYNK1Ch6gNUZn1f5aUWzV6EaY0zpVSJQ7AQ2icg6EQkBNwAPlCJjmz3WGGNKr9zDY+8BngDOEZFuEblJVRPALcDDwEvAD1V1bymuZzUKY4wpvXKPerpxlvQdwI5SX09EtgHbNm7cWOqsjTFm0aqpKTysRmGMMaVXU4HCGGNM6dVUoLDObGOMKb2aChTW9GSMMaUnqrX3pmgROQEcKvD0duBkCYuzkNTyd4Pa/n61/N2gtr9ftXy3NaraMdOOmgwUxRCRXaq6tdLlKIda/m5Q29+vlr8b1Pb3q4XvVlNNT8YYY0rPAoUxxpisLFCc7fZKF6CMavm7QW1/v1r+blDb36/qv5v1URhjjMnKahTGGGOyskBhjDEmKwsUPhG5WkReEZH9IvLZSpen1ETkoIjsFpHnRGRXpctTLBG5S0R6RWRPRlqriPwvEXnVXy+pZBkLNct3+7KI9Pi/33Mi8q8qWcZCichqEXlURF4Ukb0i8mk/vep/uyzfrep/O+ujAETEBfYB78F7NetO4EZVfbGiBSshETkIbFXVanjwZ04i8tvAEPBdVX2zn/YXQJ+qft0P9ktU9U8rWc5CzPLdvgwMqep/q2TZiiUiK4AVqvqMiDQCTwPXAf+GKv/tsny3D1Plv53VKDyXAPtV9YCqxoB7gWsrXCaThao+DvRNS74WuNvfvhvvP2nVmeW71QRVPaaqz/jbg3jvpOmkBn67LN+t6lmg8HQCRzI+d1MjP3AGBR4RkadFZHulC1Mmy1T1mL99HFhWycKUwS0i8oLfNFV1TTPTicha4CLgSWrst5v23aDKfzsLFIvHO1T1rcD7gD/2mzdqlnptqrXUrvotYANwIXAM+EZli1McEWkAfgx8RlUHMvdV+283w3er+t/OAoWnB1id8XmVn1YzVLXHX/cC9+E1t9WaN/x24nR7cW+Fy1MyqvqGqiZVNQX8LVX8+4lIEO9G+n1V/Uc/uSZ+u5m+Wy38dhYoPDuBTSKyTkRCwA3AAxUuU8mISL3fuYaI1ANXAXuyn1WVHgD+yN/+I+AnFSxLSaVvor7fo0p/PxER4E7gJVX9y4xdVf/bzfbdauG3s1FPPn/I2l8BLnCXqn6lwkUqGRFZj1eLAO896T+o9u8nIvcAV+BN4fwG8CXgfuCHQBfeNPMfVtWq6xSe5btdgdd0ocBB4N9ltOlXDRF5B/BLYDeQ8pM/j9eWX9W/XZbvdiNV/ttZoDDGGJOVNT0ZY4zJygKFMcaYrCxQGGOMycoChTHGmKwsUBhjjMnKAoVZ9ESkLWNmz+PTZvr8dZmueZGI3Olv/xsRURF5d8b+6/y0D+aZ770isqnU5TWLmwUKs+ip6ilVvVBVLwS+Dfx/6c+qenmZLvt54JsZn3fjPeiZdiPwfD4Z+rMgfwv4j0WXzpgMFiiMyUJEhvz1FSLymIj8REQOiMjXReSjIvKU/56PDf5xHSLyYxHZ6S+/NUOejcAWVc0MBL8ELhGRoD9X0EbgOf/4K0Xk/ozz3yMi96XLJyLfEJHngcv8fN4tIoHy/IuYxcgChTG5ewtwM/Am4A+Azap6CXAH8En/mL/Gq5G8DfiAv2+6rZw9jYMCPwPeizflduYUMo8C54pIh//5Y8Bd/nY98KSqvkVVf+XPJ7TfL6sxJWGBwpjc7fTfOTAOvAY84qfvBtb62+8GbhWR5/Bu9k1+DSHTCuDEDPnfi9f8dANwTzrRn031e8C/FpEWvJrDQ/7uJN4kdJl6gZV5fztjZmHVU2NyN56xncr4nGLy/5IDvF1Vx7LkMwpEpieq6lMicgEwoqr7vDnmJvwd8CAwBvxIVRN++piqJqdlFfGvYUxJWI3CmNJ6hMlmKETkwhmOeQmvD2Imn8Xr6J5CVY8CR4Ev4AWNbDZThTOUmoXLAoUxpfUpYKv/NrMX8fo0plDVl4Hm9NTv0/Y9pKqPzpL394EjqvrSbBcXkWXAqKoeL6z4xpzNZo81pgJE5E+AQVWdqbN7tnNuBZ5V1TvnyHcg2zHG5MtqFMZUxreY2ueRlYg8DWwB/scch54B7i6iXMacxWoUxhhjsrIahTHGmKwsUBhjjMnKAoUxxpisLFAYY4zJygKFMcaYrP437BUn7ZQo6w4AAAAASUVORK5CYII=\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "plt.plot(sliced_df['t'], sliced_df['radius'], label='radius')\n",
+    "plt.plot(sliced_df['t'], sliced_df['omega'], label='Omega')\n",
+    "plt.plot(sliced_df['t'], sliced_df['v_eq'], label='Equatorial velocity')\n",
+    "plt.xlabel('Time (Myr)')\n",
+    "plt.ylabel('Radius (Rsol)')\n",
+    "plt.yscale('log')\n",
+    "plt.savefig(os.path.join(result_dir, 's'))\n",
+    "plt.show()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.6.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/david_results/s.png b/david_results/s.png
new file mode 100644
index 0000000000000000000000000000000000000000..99edac26ddbed3c3483a289314d03d838c750ffc
GIT binary patch
literal 17847
zcmdtKby(DIyFEI9Al=ew07@g>9ZGk1Nq47Gk}869D=>6-i<E%m&?(*B@O$R_{@%U!
z+56h({C~WzfpKQ|#1r?u*1Fd{AC(nlu+YiSArJ_boUEiO1cFcw9yF*Z;FSpmx_98$
zBUf=bbyV=r57j&Z{2$FpR>u_r!N7+<5Ju+%Zo!MZZc^HAYL1p}o+d695PK6hCp$+s
zyLYBE9u_XH?;IUov#_(UGtpSPxjFH(vi|4$ERHT#thkoMq!0)VL{3sn-79Tx-rG<8
ztmWYdmi!`FUK90EFQOc(ZKtBLt=J!;M;2o<?-`z2m!fxT+83QUQ*R%LeO~a<BX0gi
zE~OHsTvXYhEG1Tv)64Q6pC<VE&ris`8NXcu0@IRdgG0>e;^i&ZPfpq6_%?p)w9<EY
zZ`=qio(M%kp|8YIX+ttct{NfmC-l$TLvY0}U!bwkhTsOsA<Bc-jJ`k!#8H>O{r};u
z4`mpCRD4c1#avt}%c9`-_=1%+Zgf=Dh&WdPk-Pm%eEf8kRouVcve=!fXk`&kl_SW(
z6-R9&bmA&!)~hYI7$9QNuV3290b{c%ilBv*D(ufcJv}WcBXiTp07j8txg9wPS3F{g
z27xXFx06*a2#q%6@a^cwRn@ZuCmZ5TEMrt82?i2*W6pO6F-H7-f#B_gAdK>@vHaKm
z*KWfH7E@;;<!0WqlM~^*olK|bCZ3@X&J>o#g(+(h=8}A#_c>l30ph55G>9ziGlQzR
zEsHHrirr%j&x{Kv7D|RL_{+y@nhh1EyVnw>Gxco6YU*-7y}`(y`UQ^CZtdY`Ia$jd
z#+B0g*T5!KJc30Qv6$~#dDBy&XU|BZ!@-^~!XB;Pa<YY!FnyMLTvXxe#f_W%?6weN
zRfn5gJr}W;No<haauwZx{NUs7<4mEr_RkU4F(PxgkhWve1MQkxR3MYq!JgB;t6Hv^
z29?D>c8`aP-S@vKj0Ce&h<DQd-1Un4gE+*ur@yaz@I}jY-~~x2sw}>7hjhMcND%4f
zf^$#BY>lo-RmYS*iKZ==%v&r@^><@f6Ekk}-USS9zXc5U$Dgxe{6a;IX})wBl;K+-
z_ou@5zG@szdnNe|*9n16Le3qJf4-YcS*BcvGTXCn;7L%dsGA~>p`p=t1+BCZpY%yB
zn-lvIE`FqgQ5riILIrl?$S}5^rV|UyoutetT(Rd%h{f|QxS@f6$VCrOksK?yt;N<c
zxva`3w0PIa12sq7=C?<q26h<5YHC#P|NIE{c-Nks!y1<9I#Z(<?Qt*ui0`UR6qQ~A
zm89`QL%Ko{6gsiK{i<5Ds40u<!(@i1yitAqZ$1NG?0#2PmO8yMzvw^KY@SH2Od)*z
zyD1JP4I|^NC-^y0pCJUDvxIA=!vnpB`;KMM8TXFHw}@rux7{?wFF(w<l+?X*NBepY
z;o8W`u3Wn#8;oV!@nE0o^c9z+2<ZVsGe|Bj_&|y}f(CZ*rVTifOlpDGMQzQbZB0W8
z2lJfG0J(xu<?Y<dlC7%me_nTUqSs~Mj_jDWOP!AP6OopMW}hAt3WLQgNy8o@Ztl1@
zp^3$%%Urt6q@L;CY{X&YCab4u<a~90t<3y82cHJ2l>1JmU8>tQJAcggKrmDkgb_T`
zuulPhN{n`PbZu&AikQOsU+XFj{VKbiSJrt{8}?)Nr4RFX<JpI#fIrL)TW*qR<0~Pv
zw@=$+sf7AnhlOO`y?e<s=Y({z7}!=8^<2`4BTP)*Xm0{vV<P1hBRYAg9<4q@kgK<o
zAG!M7lpFa6v$i4QD1{GHNHiWnM1C#@>DFPxh<3}Jq0JWyzCXLWoi7hfV8<nTb;nB|
za7?SF>o+ody-r7DNIA?=i!@6*9qL{O-ur<t*vKIjd`LA@G#s%Gs~GRR8k=#iKgq8)
z7WkI0ka=CX<<6Hojc1KS$=c(#3)#{k{p&|$j6uX?-1p_h?@W=!n%|aH6<$ZC!uOYE
z_yh!fzki$miYCRR<o`OD%5&28aML#Dy^-GfsuH(TFWTgf>!<U^uTUG<e7{T9Dalv<
z1&r3<t^)|PHX{dS^R9&>%M<04>_>~t>9K-q-h8RSGt8*1TeHV6Gm)o4&<J`525c+s
zF81g1zJLFml7ilz#^*X!{f<B)9NPk%^>-5chg65gkF$<#^N(K;iMydojW~!z>>~Gx
z*k31qC+b(U4)!=lz4anFZVT&jG}7;xlfbCFRQZ&IhGr(zl-^;GDTwwEfPL@uM!Y!Q
zt<csuG&J=2ix&;+1Ox<E*E?m?6=shQT5nj+&d%ghxl6kDs;a8`hK5kG^78V!1{gy&
z&3~~mr}<bx>h$#u)v&u4#g0sv#{&c(8MQhcg&l4_m!&8Q)iI^}88_}@=QHL>UUM<y
zif^2r_ZC_EEuCUo@6V3}dIqeFj8VnKOc2sHZz7CM_px8D)_jZqP85kP-a*AMm)jCp
zCwB827e#Kh`YgUj;&$_VvA%?tKc#v{C$N*!;8i<5K7Q};u;(B#A>pS_pOz+t&Do4-
zA7xEUUK^IQq`ivDHs<+zM}F+lnLj|WlR~@oN;Vx8iB4h;_o6J9@4z)4vcn0E!twfb
zUtb?B56`!dkVlhcM#HI&ty}*(mSogZ-;PgPGZ@uBdYV}&%F@gC(Dt6B=OIrA^?l;#
zIJ(T%{JLRtz`y8UGRZwP6?jeY;YM;qNboQm@nl-bXMe5%B%Ti-DS-v!y<M2EU2ToW
zs8wp-|MlgR9XCimvOYeIdN&DVNw_(KzvoRhG53&Xglze!J1XLqmVBnFbbV=T9>{6@
zxiK46eC3C}+YB4{AqFpab~raShs_$@rzF(Ht*kFQU`worgXOw2Sz_yPRsC+PbV?<c
zOBIPg+>7ZI)i0($hLb)emL9uCPHwMn4_mqICWA6<O}G!G_m=Y^Hw*Rkag$Wv$AbHd
zYVefbXP+Y*6r+uh=b_Mr9HL7w<jrTbIOu1)!$c}ocI+?je02xy+J{XxnWobHodu{3
zk>VjIuK%1XUGec<>4>3PM%Y+WHJd;sTCBs!Bq@hPQQ@?|Hd4`N$cZidb8Zaj*k8&f
z40%fbbQbT{PjWqMTG~IyqkdUfAdxTX>FV@IyFKO;AIP7TRaGBj&sTESLveHX_*=Mj
zdy|l38b@wQQ8c|;V;5BDFt3AdY#cS6(!S%0yK!0tU6rF%9Qt7$v6~1B1X#t>NrHSF
z#=?f*lYHVHEsvgIQjA!NLZ=WRj+&|LCx9P5U#jYsofzWKp@xJ}4+lq?ht62Vc~Gh-
z1uQ7tPmsw0PWY6N@B;iv#LyaQIoE@5@b-s{h|mpN1lb6*jBq~`{UA83=%lAe=jGmK
zG@jqT)zmz1eH{yl2h%E!dciAQ1~b`sXt|&a#}KG{vKd{3LiZ1xpp2#a8G%bVKS;~u
zbw+;trbuA#Fm%J_BPVGDC#k?*QtLk-TUbIZ&%i@5G~&cAnjd0*xOvf}_|?6C8DkiR
ztg(lI$kv`m5Dq{6{N|acuKx3bg$IJ1uIS(XTZ~jbI}n+uSNzqk;QnjG-0v_D52Y<h
ze`}#G6U34wiIs(_+YF;F+?*+S-^Wu-`X%=?-JpMr<7Vd+jx<!E6UXGo<PS}LM$jv3
z9hhI$c~TR-1uL1qIig1(Hr_i^@JD+fCjM||gn^(;@xW)PKcJua@=d&jy-xm7<2RNj
zyLst2E=PxSqPQXk=KZFldtLqv+qS`k!V(KxOC#olMIFP9CdH9#hM6X8N`@>6n0bS3
zi}a2_zQNAO?J-@C(>~1E59TxE?WT~3q%he~jX5-Y8hO9+HeV)xlWh}IL1fH8zo(()
zQ0eh2m#a2Ip^qgkykFJ53W9rnU@hN&yqk6@h$vJ1>YS4o>7PGrRxw>zTIT}x%|la_
zf1s1rgkFQP?PoIW@|&F2q%wu<Xq7z$=*_E{YS-lL%VGAhZT@fElm!I}leP{wZ!acR
z_y;L3zU-yU8oAE$4d;Jf<yu%w=$~8}{?Xu0u^XG&Kq*jQ+!H$=-gSK0`MqP(W^>Dp
zgs^8LCQ$(PdOMoLT^f`B19u6JsWK+X*{=bSCbt&pfwJm2Tz^!$7CB#@B-D}6M;5pu
z9;Z*xf1AkL7u4@ok+`CI8qUWFl<mXJy(f2kd38qbW2vM{Sx!4#aY$r?xELihbX-@-
zjEE(83!u=@)MrWQ9qX84O*W4B0fWq*Lp;+94I&c()A!`8k5h`r{@#sN@=aEm^;p)K
zsDVURZBS!YqgPRMPM4{p(Zk_8c2-s$b248#!#dmAu@E$9(`8y`T%lJ!Q=*V>PgK0!
z?dhnhj5-mHZ-e?meZvc8<d=;3!*W9bg~0MdvXdF_@-L^~D6cJZr@1H{FnjFA6G`0-
z%JK~g3>o8dY;T5N*H}H84djiQuyll>9_J6(LS%O>&%LPcWfA4wyd!h6JY!$eHxH5e
z4F%a!xC|mMeb0o=ZWq=~s9BZdb}Wdd*OXp)>rrCWqCm-@3^fUc`aP4P=_RrJw(}M>
z$B7}E;`+fvr1)~`vV<eqBa)$#cu&Zr++?e5h51ic+<4S@Ix#87eiN;rYG}N>u|6#n
z>Z+JhKOUa=GWu+8)hzSX7-F$e@uMW8_>s*}r#P@Q$9CAZ{q>4eG!k@enYkv(Y?!#z
z^`M2d>O9G175od#T%}vkuvaXc`b$5r2-G>%owDx{iJHE5uO?wx`h=X(a5HK8@^o`2
z$700gDtS>ihI+Uby&NYs-kd{ald)b)(5z!@SNsOU#j~N|md*INO9=2T0IG*dkva7#
zsPAN6$Z2S7+qvr(v@M}mM~0!l+Gb=ZW6UoujUJs*?Y9}=^}c`JwBzBLeSa+C)rDzt
zihZEq=hb{pg{PExDEwTCZiE#|&7RG!fH<^d{X+R`tx=tVM~J4wJFa*^5sTw~=+^8!
zg&T|ti3O0w4h0p$PE3<hhz&;6ptz0!sYQ(iA&gK>-QZ_3^mXY)?wZWYGDgQr1#{Ww
zeJXvYQQI<pl7ha&4@(T9vq((bf+4-iNZuBEJH5`0#bIJC&t_|I0-u9Ohw3BU@_Avo
zWs7nl)onD@Ab(Q#t%^Gyasu*SgWgrQwRgv_j;GX`!>N_-y_}7zTZr-)(EEfs`|Nu@
zsc|>Gxo;o~$XFOy^}l3<lUeHzY!W1eBZYa=`5d&);_<`k8f%M-Eq!;sT3Mb|Pgomd
zcpo1Rd+V3J-mnGJ`^8MCEo{9E(ShM5AbEC?Ik-^Ixesg71W!fZ6eoAO){N`TE?nhS
z5WLWU;2Q(ZCusJEkF0tQv;5u}Z%U(=vE*p8z9*6qfrKIi)(&tix5Hd-%f-S#S-=^i
z@b%_1BVrj;T9#EzHcuZ~M9h6k<n{{@L*e%_dbfvWPn<tzyJ1waoGx`i^y-~xot>S{
z&CNM=THSXvQLxA)b#=)YHH#s}lf^pCr=xF+b!#F6Z&#^X8A@@*^*=A*8Nb>`*0vdT
ziuMW3l}ezZya->3%-Rh`6~-(_#vm$p-<hNvu;V6=B4%Ch!V&q~`}tY1UR~@nnb>+^
zqgp#s@cRC44mz!5k)XXe=eGOl0qL>|)^AI#zL%w@%B*PRid-h0j~YD>7R;yj=bPl-
zzJ12UmGJ2k?)vZKt{)1i${GnGx63&4>HG;JxiU8eKgCf;Dx(nA1FWKH{d0I-e%)<^
zrZW}hEZA5%Wir&ok)k$^bCq{?ihud?MgQdE$B%?WM8};NtmXEz+zSCW_L?QS6D#L{
zU`JC5BvzR9sH&=tl1dSYqe3u;vhxTpsCD@`UWQdAoTW7=Z%lgJC}lo2R$nB6ANY1f
zxcvLuZ<TIaqnfaR{{E+oj2NO}U+~4<UbUDp$$6FqBSh~P2}QUZ6uLgIKawxIyV{7r
ze*FbM#9%W!YT|!F#ag>*aZSxc%yj@})ls^2sGZu-oSiadJEMKh(_Ie2*8TH`UO#<_
zlFm~9OVn*!6(KP>Sx!+AyP=_hlasSN;LgYLchd5t@DHU-MGp_|k&(SX$vRx|>IAVe
zcQY;beI8n}Gk$5h?n8G6?;S>gw36c0@Rc{EZ?+b$`w{xB5V$?}^-CMp$pS8y9|Lr#
z)8fO0@*y~ZQH#T(pD+xQw39XP>h0)QD)qPG7jV#Y`MzT!W}ym?iw<+QZzKmjz;4yO
zJiX@)LcN2rQd}%Cs-Fjz-F&0_AWJ5<J=xd~g%hu~W8NpUA*{Rt&R#B9m`m3k&bQt>
z7k<}wIu316$D*CU>n<H1$q)`sSzj7)Vq5W3l)_*%T0b#oeWWViKNxt%#$~A2-6_*9
zhw?=`nv*)0f34ZoMtza@&R9{~A@R0gw-_lCZV;3R6$QVH>}2hstJpEuDT>qbK2XMU
zClaWOl{?+k4t~0XC~pHkyZEP7agPdxq}|>OpMCquz17VT5v*crxA=PZ5h<EV^c(b{
zo3=EdgZu+Lkp=iMp);kBs3?w#CHlVRma?ESr9`f#p?;+}S7>qA;=LVBB!2T?IBpOe
zOtl(oJUQLNWD7D{WDELAwLEun%~O0s_wDwWERtp1rMUjf8a2gQgOW?Vg*qY+CR!^!
zn^8Zj^omnYBD3C9r9~7cEjlfKg26<qvIfuNlNR-ej$cX~ZN?{kF}GVyD7!bl#w~1Q
z16R(*I97_5_z#Yj0I8rO0w6RRV~Dzp!L5b$yUH*E;xRGwEK{h{3OUnRZF`V}X_<oQ
zA-4HXAl`&8{Rs;53QzLb`_`>LFTt3g()x9*7YEnmqB%<pp$IcxM)&*@?A(*7>+&RT
zH(t36k`pAy>tvKM5*>jZrn0vk8Ft(y>4T+!%35AnT|^+$Qj+n&cPJ9Q$@jjY;il&B
zy06US)NfSrK9#t2`ZsG%So!s#jqt$3?}6_ZO8;@CnoaM;&YH>Z=<|h_?e=u2WYc3R
z54h^Lg_D{wf7~GloG|zQAr<sEG{VH(8O&_GCpfG1Va*JFc&8S$H&YYU9Yy@w@6yTV
zd{=I<)jt`mHP+Tk({C8)Thgo{hsfmUKI(Q<C{sTOQX_EUHny8vvNo6V+gdJ7v@z$7
zIcQtLk7S+kq$Z^7dZ(AZf4M2sM6%l+FKoU&9*|Gw2n94=RbBn(=Vvn7Abz>KyZ=t(
ziwFGLQ(6a5fPu=wEUlvgf}X3!lU@F*s8m%)@ryg6O%b%GzB~qfxQHUYk+%KEayUC(
zA_14-vr10boFyJBUaK)Yn;?M2Tj1g0MN<i-fafQNT{vpGx-%xogyN`y##M~EG3_5J
za6O;wiA!1J@SzQ5h8@LmwwXTcWb1L%O^lzI+vA15#d%%4QSOm2n^~ia)Wq2=Xpmpk
z%LQ^0tG0%+vJmviUDM&BWVm#EK&!UipNH#pEp<YY^S1#W*SJ1X!{@hqdkhs1wp{H`
zpFR5+8R^oJb~1EPDqRZSbR9q5dw3@LBK&~9Y8u*I*%sEyWt~I%>`!xz{^1n+l5fOA
zAh$FAy>Pp(y?AkP@ouyWT+gkV)^K~}z(bo+DBLwYQ)PbMAA&SM>s+ikW+QiN%_jAF
z<i|!bboEyI?paHhZ(VZ$#C1n{{ljCNO4kiV2%UUV`{8nzK%O!Xg&plYnCs{4MQ>~W
z>OsP!o_90|k~AImlW$^uVoR)8v!cN~&zv6~|5~_D%nPZ$5-tZaORV?Fy7T$2&QONX
zK(qIW*K#N^f2AEb7>}LDL@wP=iSOY6EOKNmI2&mlaJ%)(k_=|8n1zsbP6c|jHA!KS
z*p~*E7`CL8wcV#29y&qp55q;A2DxB|J4O$UUFh^7xc8l;#`_*Cyy`9ywDaM{Os;}R
z#XaqlJdxBfVm0i7IY!ElwNq`4q0k)Ij=_mfi}2=w7~ov=sZg75ZF3Y>7k{tfF}H<L
z-?-8|N6f+3T`U+zg@kQA;|<v71u*ggg90+6kJNOt-_)=TElLvd<%U8l-;L1`G3$<h
zwm69pRFndnTeB%L$34&Sgj)kkLtXckT!d_p*`jD1g4TZfGT{>(T0G2Z?Fy)-3+rDS
zS`zaa3S%6hQ(OXi#lKOkt&pV&EhCHZ0u@S4Hl@B*b$F<%D`F_9cf~5s$wmxH`a{CN
zhw_qOrZzhV91^pb>G9H3gP9nqvAu)1p!3S3t*tF4gT}8WU11eq(uRN%_1OZ%fHspF
z+r998eeH@VW(bMKlUKu{an-esb7nMm{fZ{zq<kLxNqQDog!b6wdX6nmS`3GeR(l$N
zrUMj{jq_Q@gWgJ@WiL9puQ)yHWa5?ynh@Bpi`kz&M>;17S{4_+%T;~eu4|x@-(}2p
zxZiY8Y2Jr_fBdaYSzR5B4EzF?<#>Fwu!`5yYH8>S_Z~-2*x0`OL6rvOojBh^lubqY
zZOvKfFwbuoc<1iWdSwY=B6=JC1Vp{tk0;WWp7n^ulX%8jWa8VZ?7H?PnG5nN?kdXX
zl6I~^E!W#cyCBr`PzR=c{YvQK>e_HJBp~qbmGwlOCl%L#tK3D|4VQHB97)=?rl!>c
zxydF!3))_CjPjoQLY#u%jbT{iq1jgW#f#9)`_p_%)q=O_g5F<5?zWWl>z#t@76Y2p
zVqb>f?w9j5jx=>WF`oFY>$NpcGCw?rbeer?|3~_;x2%_*i+`eI^`=SF)|fjr2H>ru
z->-+AfZ_;Q@I4R;IF8RO2jdcl)9kb?zOu5?+uLhh=?7F8Oft@h>S|8SVr>Lb8=F@}
zT4m>#+F)ZaV`&P6qYs^2lqnvC*Ws|AHKgme+!81zYT7c^v$uZ$QC;JsFRTsX#vh@W
z<rw2JytCCmwr@`j(T%y0iX!@)nK@kVysF*k_VVcHXm`3Q{0XBbo9!eMxCG-4#6%{Y
z&fBXKgMWj<)i@D4cJ|M|f6D`~((xx-LP$6kSKJFaWV&jtf*d}qerT|ajGfugjw7OJ
zTJer{V$jZ9+Z6L^jSHTXVRTJA_&G%=ni`l#S`9|+o6k*OMMkj6C@a5Ii1-GoD>&iU
zl;ZhgDq31H2MaBHV6j$IRP6B|{tp2Ys_|jtUyu@BdZK)zp@BXO`{G&~nYQw#^`YX+
z0XieEq3YUp+Q^{Y^JlavQX&c#UZ$K~_&IXBL}57$&F2<WN_O`{Cx#IkFZSJ+oFWy)
z`Hs+9d@z61|CqYlyVEa1Ai~Fgub5xMnIcc8s-`weX~_JNiRt;vEmJOjFa+=rQ(4JD
z6?p4>ey}(^rX(_G7AG{h3rC~oR#r&L0E*ID5Q<;2Mhy-Du_TSDa9!K^J7allW=T$C
z$GUErl{o|Hk&3+Hft0j3Cn?&<PnlF@U6S>!*rJ=GpX&#Z3maq2rYDdOgNS%kf@6&Q
zc&`HUzIh8;*$y?fwPb>*-}-Rx3v+D)bOnf?nAE})8X78up+M1f_o6j3Gt;TD(J0!n
zu&@}(SE2^!wDWwg<1SB?_Ps`W_phyEGR1@mCunhsng!S21PxapW4w2*QD|)HkqYC;
z(t1=EqhKmphKh0?`MkvH(R`=G9RCl){KM%3y5pY${JRdPt4<T`%CtH8dOeS5DtWT-
zt)f>huj2O|4>okM%5!DnJf!&`va+%?A|i`iM#{>f;^O4emJX7?#0j4|6++*<(?F6U
zw&8PaO?<tvQ+<I!oIk8N$umvYZ%zXUR~La<lCWS#K!v5MPsnQ;mM_*-BeoVh#N@GJ
z`MC=~hI<mRxppr`IeR)9uA-|IIpNTq<+{Vy9Z)jydT+MQ5hw!mR`$g$|C<!c=|g)P
z2Zg4URXEjFW_;=Qc>~?%)mu(0&9Bdi&B2%{>acOFEa@wFCQKK7xyeR{y)9?T^ZL$o
z@H7gGo^Mh<db(vs-Pd{;$@b=Y$S?v$<BY9+1rni#Ii(@+z2^HV_C?mqQKcF(<9IrS
zl}SBE3qXKQ_k5kxO2fiXS}#{H41M>w=oC_ZzDLHG<3F5haD@xx-@ku<PDgh<uBKFB
zJ@QnuRKM@*OSST8J0|`5Pp7ADb?1x#8~7<{XxK1w_ZN}#YT0Ku{|I)Sz)CB6m(|=D
zJo-w8S7Xu@+|el1M{BtCm%KEk^Jh9A=`apK_Av4>g`YhkBlNnyN)c$!e%YP&Q~z05
z+SiWY)@a)QArtC^lY&c$Fgbfv#aQ|>5n^ix!cOyXKjZIofp2qj`lhC)1w}<&mR0QD
z$L0Yisf%8x+6f?laW(nPRKEkF{@~snpB$Dv>VZx9OkqKlvWo4pcu8m^KO?aXb6U}4
zXD{iS4O=FK^&Qs2mwa1<n^aV9hI^HJxC&%pZs3WJG_P<J2Od?<0mt$<QRK7Enf@zf
z<8&56`H_8UHhqk8Ln3>ULgSs@5u2c7c{(buB?S0`i3kG>TS#aq9s$7^>6&ZmczcL-
zzfF~zhU52~$!$~`=IB@36v67YElgQtfreJ^gd@^Se@izeOKe(l<5QbG6%vlzU#4}#
zz{zf!KB6iu850r=VTi{|9|{!jV|K_$tg8JH(o?f1miUb$rJYbg-yaLoJ<&0duuxLu
zFP-~Y7^uqXIw82&M5*c9K8y+Ono#m$^$DFqXq~EoV=i<-T%#rC)rQsC%*Mo-=Bu5W
z3#jx#3Gye!-~zU1Ds<0OOi+Nh`cq08*MbAzo*UnSqh2@;du;L}2k90SmAG3H-iKD}
z=R1`D_5ut_Kv4ueGkbTT!BVG(w^_^1l3KPnwZ(P?r+y*<3`>ngvIIGs@<X;lTv*&T
z<aX*j+6H6rrxYhM9E=mC!z?VGVQQ;?-W6W$snX@fMkpgKXQRM;4;ToRW4C1EgAZ|$
z@aUM6(oFJ@KXX&#*&G&GHkJ;hD^hq+v{}ihWyw7C$u45)u;x^?DmwY+l!)>vxKd^w
zO>bHz6id3L)4SJwV1^zx<kc1>VR2GLcjMes1r^#gXFE0a^CshrA*5$X5|~L6;OD7i
z2w6lNy8-6Y_~3|7(ql1*uX6e}d<yT}*fS>Ezm_{>P;&|r4N-Ly+U5+sLYOSx-nQxN
z?EJIkMj_x4ebDxhQK*)0%IR|;jE+Z<mOaX=Y(p8W>e=T9E~Dj{@pds7f!%Q(`|*Z#
z)>7E>B#JTnl0MEKu$M_OKAqy#LG0lqY~h<gAK!fP<c(KE$C*yjYSp3im<g><C%FR8
zLbwVdeO3c3bajP`EuD7!<U@*C_76-oeU7u*>q!_64ZuQZ1a#$a7?m#2^S3V^6!HF2
zKU_Yxg_Z909ychY;Obe`6>(1Oj_d7>>)DU%ktyiyacWbZ#6p|`Je;q(-SrFuuB?XI
zkDjZGVkENVQpKyp^%b}(gp;00hLV1a!hFFb8P}JEEVKDyCn=i0r}mm|piA6M_mWqf
zaXc%Yj%yiD>Zly$hog)XO|iguMpps(QX@ORb~RtSx<<Dhca_hX4SbsfJlrh4cKO>I
z%akh9=tGbLg=*WyLK>daY%_cG9REOA#!N!ld4XWk@Z2daQzr#2_t#B@HxW$?sT+DW
z3T3wR9)^+c7MinsYBWDsuBF4`n2mHD*$M)Qp_^##A+dc~V+)i!+}N~*EfJ-P2z5~7
zQp1RAk2yR0eQYRYOHTbhBL<NJ@>HoolXjCw$<DMsx8BlS`EG1mL@m|*QDu2O35!Bk
z&N~t=tk3JU&Kpl1T?vupV>+LX7lg{JmPfFM>l3ljhDq;mRVfsR9cV9{R=xcePU5zS
zR6|~kyO(8vWGm1Fy<?wlg^ZzS!HSHl9xr$vCxB7N@ITkd2)KS{Ih>*3<HJv{lwm%a
zCr7En9amk4R?*yupk9A);H#(K^MoR%byI>YR7$cSc<OMy>*-(g>_2nIbhEtYj(=W`
z{(K=xyk<y{LuAZvAb`kJERf}&{1&Si8D`)~2_EUT!a(WT$P8?iK%b!6vlp)xkcI}_
zwAyYGQMu?0zK|@#m`h~~#ob%+V3q7yRiAOk`dy5iWT>%52I`jzDdW3FcbNg#g<2Tn
z8kqY8Qj*r~^A<F~L=w4T$to!mpJ{5XgeAZ8zD#RxG-T?FRfhL3E(0yC9)W}A4;?f{
zhX4xxb-?=VP?uBi%4->vVl=*^RlG3!75$DQaIXf&;V{QrMl1`kasqdf<E-x$)Ps>f
z=$w-_B}Z(;tb8^PW%A4Mu?wtv-oXkWFN0%4#s(aD)MduVl?7#0sLfd8!p;nyaU`M%
z`{o{ZF?8bO?vl37#u!NK&{sfa0;pQc^SU&zXz?+SXU+HBbG7y~ncHW^LQrwZLKhsm
zo^qA1udT0pnDYVK%fZW64nnHk4t4KltMk$D+6hAq?b~g>k!cfvA$SQh3KRmg!l5d2
z(Pe<ykX2oahH@<w*FkNr2?w)Dqo}Her(HakZ7S+JLC<&Lo88x?ZTCI^-?&U$Mdm!R
zgF9-P=BPe)O?=EvSBUswjX4qOg0-uvBy>gyDQms92`82Y<B;6M%ro?vcV^-vQnjt|
z6QZq9ufn$bvofvQ+w5?vV_lPGF(07d_^rK(Qv>?EYu+v@(3i)jxBq0XN@|GRG>v@l
z0LtfQTjOug;&CdKJD`^y4)QFrnf<~P#fU8RqE0#$2k>w7JFGo6k#g<3qmila{kbL5
zkGlE&1OORQky9UOt0@T{J<1Zou>;bFq_i~R;oo>AmoWt%ce^dWtaPMq4BGce9z~zI
z?W<<Jj@qM&g_~1ARmAo2@LGy6>=Q+L^tAD>(wPRQU5T99+g|`4AaM|xMny(Np`kTg
zYKhCrqIP$8Q>rq6s-?^jon+B`|EiJT(BlSbHUWF?wt-v9&$OSvl`Oul7MQX#+zKtw
z(kKkBF~7gPiYDiY1%Mo>+<8@czj3D&n=z@ZZ|n`AH!dhsuYEI?Kg{8HIh}r(d%c^H
z$WzW`eX!>s9z9N#{}f@e8P2oF9h*M`AqVk~&x=?xE@&4z(A4gw#j=$vD$X}`_eorO
z61p(jNqP77Laie^v!K8N^#&LbBve)Lfl32O6ZTJO^}q3u?Evd#QJo*<crapcX}JQY
zt~Iy_WpaV;r!}U{q0k>y*yIx)$A#ZP2uQMWa?9XKp8<jfvVw>Cv_7<`s9~qHsUG+`
zbecTu^C{F&Y3D4Q=FrZ6DrK~1*-G5CuW!07$2go8Xusc;Ugv&16?nf9cnk<OyVH^=
z;EB_<_T*Yd{xv<14gYtweQDlvagDg`5}CE-6hD>6@%%Cg=-QvDG<^;uUxancKd!S!
z(hBm<CT~IKuGx!oL3?D6&6YbNYiGyg(L~uaKa?jIKLxC1;3Vb72zqr6_Lcv(hTI-X
z#8U8eHoEWdy8J~&Mn(pn3WO}!j*##VP@m3|d=Vl)5i|3(4j?x!wkPjxIAghukfJ_*
z%!WcEL1toeTo42-5)}vM;<BB9PM(I9wXCmlzR5E_F;N=ONQmqDZxbNx0@SIhtlSRd
z{)$<F4OUF4u#+d^^#k*d-}4>J+Bni<=aS~0+`pa(%UHUUXdC1zS4ibH;ZD7pUA+Gl
z!?Pg!aDSIdxl|{Fj)7so#T@kWrzCjRu$HLjwUO>w|8RekEbQ+SHEU&M_1a@s`&TUG
zTRFLpS_5E7l(?SY7&&s$zC!nry5Yz9kr4qT96nn;{`A<U{YDHhIN;uscpMGGpOhX|
z%yJ|-X=)Pb=;%N=zJA5+l$k3Co+g7JyX}>+lM>ofcW^%SFXAeXj6{crd^mZPq4p1f
zLW}hqzKoBnFZk>*ml4s?(gxq&`oa9f#1PQX&>Z&61`?UgEiC$sTR};F!+{`YFvsBG
zxU&~c70LEjy61viox^<ndWx;r)!H|mRzE(V=zv`z9ZRt_(P?gJ35<af4h}5fP{l@V
z89a-AV<rB*;8>`O4WnMCgXKCii+-EN%GG7ATR%s7`{L=ZfLek8cNJ}p<>TZlW&Sef
zsM%?E-59dyiKYCSnD`|rsaH)&s5?h0O1r`3h1w7K4&W35g23q5*!)b{KN8aG{mHyt
zBI|)WQ+s14PKh^I#dx_#rUxV~hDrdaH+5kc|2OIWAIbl#o<ancVR7&l0#`w4FrX)o
zaj2fZdW8<r&?qb^c>$jNLm?w0Be7ZY1e!IF;ye68araLQ2yzUqAn$jw6H8E%2G@=D
z3jaA=bWZ(+vREqesiBRy+i4#V@B`%k5g&WzZ`zxuLJC|@#W8c}esVAN?y~+zU1Ow`
zB@*{3!&L#XDeh!Dkmtb;w-O|M)Bm;wH~k`5C#_2j8Q~p>ov3BU3apYXX3t}dedx5?
zC))GSR8W^+?|T^4>a<M72wzICp#IBa;c_F)cEA}DDj;25R;li75qJ+bmeDr<f?o$L
zWNM!MtGc&AgVy;5&R@8P!YgVvows-|kSIyQQwpQ=cxDN$Jg6wR4=%kPc)3s~c_K1N
zTG6O#<e?D#o@@AM2nuqr&bo;2M6I)s27&BK*`g8yh7e&xA85FuQnAKAW>|V%8*#yV
zG&|j5<Zlg+-INd}us^rZUg@$@>M3~^24C^~<C72=bX(!Q#F_G<KBd5{GnY8NjJ<k4
z@(v-99$~!=ZZ`@U)+dwCZ+owRBid$bvhRZ$vsrJp1K3~DE=CRnqh65$<a)Isx}RY-
zQ%=;9FWS>;D9H_%ti{#&v^-g?e@T3W+fK^fjfY&WVK}i6TehP9P~3b{GDcvgrlQjD
z2bZD_?;qP~9jN1sf9)zIHtBVfC0&&3nKbhlzqcu=Nq3c}j{oN%pl~yLhd-{>{!_$b
zA=}-jbJqj1UF@FjsZad{VjH}&kM&U;6=7zR%oipKgl*WLA%W5yS7hz-@pi;g;vCEu
zFG3DUrGhc8W=^Zy-1aNTC)(e?0B@lH6DH88dq!Lc76?!lruxDs&0J1}|GDyIDk~?<
z<o%kSyz1nkpuQs}wg%uE!tE5IueCjpuynhxaQ9hZ*1Wgv3}kQKI2;mb#d$Qx#Wr%L
z!$gFRj$TlN5!;3@rA8P35rm;=>4WIbh3_HK?^<x;181&!G~Fbp?Bx}b!P9d*G`5i~
z@Fx+`(a-t#NTpEHF<(KUs+yY8ATqps`LdNI_G114d0fk95AUVJ$QG@v^pewYD&5a`
zdZkK0&;yR5jNJN<{`%Iej!#VBk&`nnKM5g}qUGl1cCuYGnWy<t<``SDDV6{^X$-|6
z#wRA0k(QQ@Ujq@|;!l<saN<>%c1JSq!)v4yMOvQb>{9e!4$%#dn+|AHk#<PnUO~XY
z6%wa`wTsp8cQ|zbh|n_Hbq<s=u@uX|tGBkE>gaIH?|W`%MW`0RI0G~U5Px1nKoW8|
zi(>)xr<}*?PIDMUnEY#G|8F+HIH7ZrgO&?JE5f&^qD?Hh&{zBhvkw`m>}J4yaha8$
zjLSU2*E*y@9^mF)VBqJG_+RXOKozrdz5`kq5OszBuNF2>Dn*FE1pL!tA5hqVTatr=
z18!9_>xq3KB$Q5{;ni8QFj=BUD(J-qq@s0@K;4;*sw}bICNd#Z6NLdWb@oh0GUR+J
zIw$8TysQ9jysx*{cw;C%-Ob}E!`A6c(<;c8OFS*O|B|{Up1^>gkWN5=Jm||8{9)|8
z_n03@C|q1zya4>I2LkX*R@Qvqi+x~|bEscN!KUQ@*&dA40HDFu>6j9zRIh5aT_@<>
z1fC4?By(7x$BEqY0);juB?X(GfVO_>WW*A|EXi?TR1VnUK&=0}vC*`SfsPI&E^C-@
zRb?eSnCsQrncST3?ah2M*WU258StdwQ(K<!$Bv0x`M9-JQVJ_G5-9}`6avpgz>8L)
zTSEBciSyrHyj>3gu;MUD*dRdy`}JLCk5%<RYN-T49|#;(mV>6C7N3#+3D-W`^^rF2
zM7!T)NYa0=AD8K?eX78f(hCg_NBWeSI`Dc}u!}y^uQOXByngFP>fH9aJsU~60885J
zo(CBj8Dwwf01y^MadAYj(w5F(Mi=*TfJRU>_FEp=^YN?XkGHQ6r4PhX3NSD;zW|zN
zeuhsZ(2Ofw*53S(PkujWX$D-gdHHqsTCTR{oIVeui7BcwMFj;kf`VzAo7Mn^UV*Au
zGOzPfP$8lNs{xcS{<H-Kf|-Zyo&%5M?tIfXBTxjeTWBVYlDa|UbKfQeRsK2OgO+g^
ztjcjwczrl?F*%>^{l_X4OWVVLuIiEY+lN)9RPLB}V?Vslw8!$5tgcVDfCFH%=;leV
z?av6jXfm$Ijt`IHzkdB$12e?1mY}k*vT_Bm8ekSukpRJjgoJPQfndFJxUS1fN0W2_
z46_NH-UHUkw%avUk{IcyPoILKF@NF-S(1<!pgB<f3^7WYG_G<WF2={lub!U=dacpb
zzRbV$xu?avUfPVDnTGq3Kha9~8DK6;96M$qQC!nE`{FR&qxfK-dhvjzZ!w%9WXPo5
zgI4co4UcK4Z-#7*#nk_06pYBq!x>rlP7wZZtVxeeTBd-ytd&RJ@&LCSN-BSGIK5YD
z907FTO9%UPSj}YJ+NLBt2S<qFLSs7Q{1Y4pMn1xGsdd!*_lTB5X?*|@Z!Wfp;0Rx4
z**?Oey;Y5P)=mc0a#U?76!5A#xX2H-2d~F2iH3%T;Qi%lD@zJ)21d;I!iNGBI#|0d
zjXy3DEI%_Kvd?@4>}?d$GT2F7+g%JGd%$RMxtSYRP+CeyPaiHr<sAZa%>}PjEO7~m
zD`1=6Txb!XS4=ZLTIv4!^(!cKcyE3GMgn#_{KbW4Z*?uL*I?}7<Qxix0;icN5J=Mm
zJh7v?J1Z>)lK5PI#0BC%&8oLC{uW0>PykOEuy#q?)|-6+z!JvHb=1`f{cp~@ZE9RV
z@npgO!l3nPogA*l|F(R6b8)bdVq5Br^%%@_HwJ5S`$CJa_j<CG*HZ9fh}U7qV{h<)
z!!S8{dAGeqbB+p=E|l}V*&gulzFZDRD;IMLDI*oH3D<ptB_E-9yqtTuTDLH7bzs8g
z9@>+tl*0M3w^tTkc8!g-1gj5%NzNS&ly?;kjqYzo0e$v$^XoG|T3e}t44OPD09u4N
zwB317dH=-|wXFq)oQGbp^I{eAbIZ#=t*n^(T<q^T^MzWAt5oZ1;~D@1EP<HS=AUeZ
zj4aHVh*LO5G-kskBiq(39oy1B;ZMgdnuwftpr@k?Lwbx|eU$|Crtf)qAuYZarkrHJ
zRpPPpnOVGA&N>gMp_`n01uwm~!^j<$Mr5ehz&^JE1P{bPs(?!qA|~yiX*&n^n}Pxc
z_@r)#JisqL3<?DmuL!(n(oSa+TA)~Va9Ct?7kRb=YpQK<T1EuEBd|+q9T$fH)#5)U
zlqyI6_{rAz-tJ#e+I#M!gjDI<NBmMMQ{^T|54Rg4q-@3zaIHVVSo>Wab03(2+J7-9
z)q+!e%E%Z3<{*gI`Lxv)sMPj@8ZuwNZT*h!E+_&6W`xrfhP4X6?!%=HN*9Nxy`j>N
zG_aVKsPb3-Zao9qRzIh(isH$j!D7ysoYBu%+}^m+g2F;tCMHZN745bFA&^B~LsRP+
zK0RPvj_*_f!hd0R;zfLa%KPM<xLt7c8lzx^=Yb(04DX@PpFe*d_A{1UHp)It3Xv9F
zi)&?(%7=GaME#1R9tLs6-vAS@Sp4Y+D6wg~0elZ;Jt0Fl8~rE&D(G#_3y+DC49Hms
zoD?G9b{WY(=rMn&tlyP9&h3MJ-*N&?LGjsnA-G}#oAq%cs;jduZeM~h*c?^z$W2AZ
zcFv+H_(sokD>`p@)SFz5bj$!dTLl$3Eq!`EzGPthjypeyxdXngelC{3DO|R^U>7wF
zJZ%=G&S{x_zd9m{01IJ|)B}2W0~g^}`b1WvA)rV0aF>-q0k%2rcuehPT-dw{_K~YX
znNC$$bTsD5>Z-r?Qv_?I>-Uo(B#%JclY9n;C}6PR-UeWjasfFJ2`)_>gJNa<<UnsP
zAd)yiI#rhJmcK~=D}g)t;fw2gZEYezC>(%xLhi3&ZC^fr7L%4n85X`_f_Qs-!;6?#
zCmZqL)qsG2cwH6O8;3rqM$zT`LBK}Jzw`gnLEwYH>aQ3xkn(PO0&gF4r|wnmr7}{t
z@07I-fwe)(Zu-b`!P_c6vj<?H!3-f^KoF^bN$v05#iA9kzLnaGfSwAV_|1o9HXow5
zJZpcShZ07X1tQ1lh)0{AkOt-5o5pW%pOFc;{EaA!%G!40cL#Z+6n1^KGchxx%f-&F
zBIL1~WK%agG3(Hz&;rt@x`xIC050Np14j2NGr~eb27`lxAaZH|^BDZSloT>9@FFua
zGyKJ_t}cU}gXt=n!QNhF(BCuO-~YxP?C5PEYK((w3aHe%!%7WG-EYoz@)N&))puTu
zfjI*tS3Cm}B<vPsqypgfR}%sk=mUno2n!1fM9OV&hh?qjGqwunW2#r!)WXHsBy3BS
z$DH@x`kwRd+W`ILf{)GYdAI|FxeXK?a4i!Jt{YGQ30rqJ7N6W>J<2=ahHew6Ky}@E
zR}H+b=eICMLm2$FhnpULU#u<qx@Uf9!db^zy?5$%Pq7sjbXtP^7if|@p^MXgGAuG2
zhD8nmrGI)5|EH#==yv>$x`}^h2&KY(7Y+?;L>~KdebJ<^H(_Pe|0pOo9^i853kwVV
zKtusRCkQUz+m^K?5*zt_27&8Kax$8zy*)DsE%vW6a~S44X6;WmN05Smk<jSKDw_cv
z%sG7^(6JkUyJEV=_S=qb^5!EHnPCZAr}6e%x_^w6HwSG>Qc}oKM9ikgfBRZlzMwuu
z`GC+gzn|hMj1K45cC)oHASMGO0B}jpY@D2P*st$`(-)v0#%rgf4xXc`tcImk2tir3
zv}u23D?ihtNI#<yoD+x`{D4+~1%dUU!?D%!HWO6Wo%E4*z6R=G5VQORCeCLpEU{hT
zI8BL;^Nk6Qsr_(M2Y(4&%ys3)3D2wuTMnjl7wc5<dhS!i2_8p0(hjmN`7xf1f<q;5
zXh;Q8V^-k9JuqAjo$0u^Z~%Nl1p#RC84nNfAFBA#5<MP}+ZF>KgpB;p5^TwHM<TjS
zK^=x8;MNO-@pf>INGQ~9syvtnWa7xRhS63%4nPfnSwW;C@FCzHL}8?vRWmp^a16ZB
z1Hupkz{Hi6Ujpu~8)|I}{b*@v;WHRc<r6VArQJU;z*-b4B`vM0<L!JJetGEyB83ic
zQO(wR3l|h0ihgom+58I81DX|p?Is|nlv)794E$UtM?KW>NlCJnmUOQzf45s!9Zi|V
zaZTSRXc;<lc{AzNMuQJ|N>@lhTs_L28jpUF>r|K`gTnco+n7Sqpe2a22iK;pXOnu~
zcV{p;C8Z5$QK{)tdvNI<aC~+HkVV32tsMWO9rXHy;ZP3&FBZ`EnBjbBuHG3k2fD%l
z1ni%{%j46U@q7yuq+j4~fPxmWqYZC(ny$2v&k!VSnvWf^{@YD#WQT7IScpTd9|Vs2
z;flUbONL6WbOo3&x*mOSgBSCj)Fvh-a1%*+IotYRstgdYqInmm^9W>mdO;i*8)Nl|
zy*wUZod5z))A~17qf!>5*5>ZLJ&&gAo5csiyUW$MSvNr-0l8ldAm0KU*~IURdckDa
ztnCkUCB54@qkF8&^c)FY`q*dRJk#DWf{&v77n%j{qSDiU2qiW&o~^ZqTTWJc;&>Wx
zx?c(i%({=X<mLUm2_+J*wgJ6dA>rXCiw}2;tt@X2*;8nLEOw`82c4uA6}|n&VzD*)
z9b}V#bxMi)GNj_}u{*Cwf*z{P9n;HnuoJHvQ9BrGMQ+zWmAxH#_dn^*MSuldr$EIP
zCja>C_3QYwGzuV<zO%MAzPmmHXYbPL#!HvS`qCp5G0=Bl==C=L(G;j~HvkC@U@K}_
zd3p1*?FkZ0l>b}dCIO9gaNG?h+VP^oob$SL*qlQ}Q<KifTufF&qrbnOEgA7ag*s=|
z{SATX>diZ{)PFfrccPz0usjFE1VG*`Yl1v^Is(8RT3NzSApv6_b4&h`X~V1l9+Hzg
za%Tq(y(~-8x!>R|JgC@k><k*0;m%OdhWpYc$OCz2B=<4kBvkf#s;c4us<-H6u{{4Y
ztTZ$jh<to?HP*1Uw$?yZFf)N{7gW<TRbd9AD`*EEF+se#{X2K7MU!bnlbOnjtPk%A
z{*MmdtxUh;K6-zHh*E1XG@xfV9R2S}yz7ZkuwOv1pF`_Y7|_2!TlGJJ7rc|2rQ-kL
z<vw(GZhi$)Kd@`DgZf>Nm*GrjD8|RAsAtcfb=4F%H>YI+cm<EA@ZSJJ4~4%t)py-E
z2kH9Dr%!O=Ov>l-oZ9a&80YO~0}wvpwgceYfp4pZjSM8j4NutU6>|SU8RKWVpVGkn
z4Pe^e|1|mlvepAW?r*?ZCab6jlBy7h%3Uz2qe$R~N{fg<1r%*~xhovu^6nr|33M2G
zgI-Daeg<59bvW?8D(>xe>-x}MAARQYr%&GluJE(xyKN@m^+0AM5%S>#7#-w^uo!Zl
zO0eI8F6V390vub;G(~Uk8}SQ+xhj1x97ab+Wfc@&Jbxa-=ek}F{B%adZ>x<vL(tqn
zwW_ZY$kEPV+afrDA8@p7L|8JPB2Xs^cfLr+QDKRin?IKb!^{SI*EqnFIUV#epp{Dl
zjau%_8<#-<qfy311f^7x1id*x$&40kC|LIvE-0OxZjC(-0^~tWUw?NM-y7DSz(^|W
z#|xtU=`3CCujAw6skT57EdzHA(ffU}OzKAq2g9)E?eBX*QtJSz{u1zCZ;%+TZcm2&
z4f3$RGoVN2TsD#zx<1@pxUT#}!=~UB%@p?c0viLo?@<Va!tEk!(^i>4jfM7E1s^sw
zo!%nuUBI4z?i?_ogr^xn;QI&e;NJ@%8W7NG+nSXX1o{lA6MLJ|Hh{*|158}gExvVe
zZ~ywY_?+bw6rlbol>&XOSM`z$V?|o5r#?bJT?Z5!Ru&jY6cq~Gd<k&hCE$dR;Ik&%
zty{#k3z!D{bpKspfIHN8OuePt-MPT5m2SleUl+*X)QPf#q8HrJ0)hVk`DF`?Rouw{
zr7-~R0XH{l8VNu)U^nl@rSjRNg&m7t!RI*m_-02{edH%d*<}F5h$<*3XnbSVTi;pn
z8vIB=KMvMZ5y+)lz!OG>yTNTWL{S483!ojw$^`iM<A7!9z_AwaVfk<|Gb^+#hmns1
z5z2^o?Ac)YET!jMwxzjw0k{A)K%Yu8_*{R%x=(vSonAwflG`qc-}EOc9P@31(^Ipx
zH6;$QDhC7_e00kt{wFX1srmo&`wb8o7#Uq-@E4jq59o5k9{$pcY<)7@Uu??+cZ*PN
zIY>q*WddS#=`fGOTm&$k-JQ=m?0F8Q@%{b#*Q=&j+)WBJh(LhFje^%H2<$qiMgNN<
zP*VX7jBDSFJbKwd8W|e>{?_0G+=dwluBBJNaKjF+yzcxRXiNr$m8-IcYu2&^wTRL;
zi1HDdl9G~;a<d+6fL!>&Jjj6>511ncKo2cIW`vVeAi+R@H?6azZs85c&&kQjw;(;W
z+GT<X1(dHqPBeyEBm<O)%qp}2FJHMlS_PB@KS)XB?kfQN4U_OhmMr93TjYl$a?AEi
zI@%cE*l6%K_OrDnHM#=3R!=3hKu358IX0DGFF^OaKBo+h%||Eze3^mp0ME|Nhn*M@
z!ONu&x94?jYoI-L8I)O;Gko_O>_mLecN^>&Xr)vxr-uW6fs!mMxBbkxIZ2gVJlzq1
zziO(gi0wsXt^FVjBG_A(7nL?f0(t3pEc2c)T~S6?mkbm&b_=LK;2OMtl!WWh`}=pE
zQl>BvC`JLYrc=tmv>pruWCQ|7f-l}Cvbp?~1?yT2phkd|70%Vv)g?iYW1W2fknkV!
z2rIJ!P`q?4GRFy&vyLy?9@0-XhO2C0EI$;}(Sr=jHo%s(ugt8k=LK`SI;cjRGr)cc
zzBS`A{?8v}ShZrl>)?(XzgB$uyu${#@=j3ism}zlXnT8`MNb%%=3DR1^#Ha5OTg}U
z8i=8TAnvdy%flfs3Fgd0D_uNVi{YaAd}AX&u)hcj$kE=00A_z!2?rY~1PCCCI8F0|
z2)7KTYtW(1v>P2TzkuoE{{(HGU~N4aaN3<>C*9M3BnTRxns=`@GOfU^0_$KJ&<HDS
zCZGe8a7mV!Mc+wi<wphJ_8jj<^DHlf0HFgdI3Q8WM)k0=u+T9vMT4m58-Oo<i9iba
z6{z8*Kxxp5RFy8M3wQQP$B?Zx1{Qz6la`V)aHgb142ECtXRz7Y`nYN);Q=2YG3Opi
zPkrZ>q$Cn}jrp8NDS|}_-q}HW`vrXbfBUwPdlBr|G8=9o%3<&&D-by;Mafce<M;mu
Dm%8V|

literal 0
HcmV?d00001

diff --git a/tests_david/testing_automatic_log_readout.py b/tests_david/testing_automatic_log_readout.py
index 5cb370e2b..5e8f636c5 100644
--- a/tests_david/testing_automatic_log_readout.py
+++ b/tests_david/testing_automatic_log_readout.py
@@ -9,7 +9,6 @@ import pandas as pd
 import binary_c
 
 
-from binaryc_python_utils.defaults import physics_defaults
 from binaryc_python_utils.functions import create_arg_string, parse_output, run_system
 
 """
@@ -31,13 +30,15 @@ print("The following keys are present in the results:\n{}".format(result.keys())
 
 # Cast the data into a dataframe. 
 df = pd.DataFrame.from_dict(result, dtype=np.float64)
-print(df)
 
-
-# sliced_df = df[df.t < 1000] # Cut off late parts of evolution
+sliced_df = df[df.t < 1000] # Cut off late parts of evolution
 # print(sliced_df["t"])
 
-# plt.plot(sliced_df['t'], sliced_df['radius'])
-# plt.xlabel('Time (Myr)')
-# plt.ylabel('Radius (Rsol)')
-# plt.show()
\ No newline at end of file
+
+
+
+
+plt.plot(sliced_df['omega'], sliced_df['radius'])
+plt.xlabel('Time (Myr)')
+plt.ylabel('omega (Rsol)')
+plt.show()
\ No newline at end of file
-- 
GitLab