{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Day 9, 2d & 3d movies!" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# import our usual things\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import ipywidgets\n", "import numpy as np" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "planets = pd.read_csv('https://jnaiman.github.io/csci-p-14110_su2020/lesson08/planets_2020.06.22_10.10.17.csv', \n", " sep=\",\", comment=\"#\")" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
rowidpl_hostnamepl_letterpl_namepl_discmethodpl_controvflagpl_pnumpl_orbperpl_orbpererr1pl_orbpererr2...st_bmyst_bmyerrst_bmylimst_m1st_m1errst_m1limst_c1st_c1errst_c1limst_colorn
0111 Comb11 Com bRadial Velocity01326.0300000.320000-0.320000...NaNNaNNaNNaNNaNNaNNaNNaNNaN7.0
1211 UMib11 UMi bRadial Velocity01516.2199703.200000-3.200000...NaNNaNNaNNaNNaNNaNNaNNaNNaN5.0
2314 Andb14 And bRadial Velocity01185.8400000.230000-0.230000...NaNNaNNaNNaNNaNNaNNaNNaNNaN7.0
3414 Herb14 Her bRadial Velocity011773.4000202.500000-2.500000...0.5370.0010.00.3660.0020.00.4380.0060.09.0
4516 Cyg Bb16 Cyg B bRadial Velocity01798.5000001.000000-1.000000...0.4180.0030.00.2220.0030.00.3510.0030.017.0
..................................................................
41594160tau Gembtau Gem bRadial Velocity01305.5000000.100000-0.100000...NaNNaNNaNNaNNaNNaNNaNNaNNaN7.0
41604161ups Andbups And bRadial Velocity034.6170330.000023-0.000023...NaNNaNNaNNaNNaNNaNNaNNaNNaN8.0
41614162ups Andcups And cRadial Velocity03241.2580000.064000-0.064000...NaNNaNNaNNaNNaNNaNNaNNaNNaN8.0
41624163ups Anddups And dRadial Velocity031276.4600000.570000-0.570000...NaNNaNNaNNaNNaNNaNNaNNaNNaN8.0
41634164xi Aqlbxi Aql bRadial Velocity01136.7500000.250000-0.250000...NaNNaNNaNNaNNaNNaNNaNNaNNaN9.0
\n", "

4164 rows × 356 columns

\n", "
" ], "text/plain": [ " rowid pl_hostname pl_letter pl_name pl_discmethod \\\n", "0 1 11 Com b 11 Com b Radial Velocity \n", "1 2 11 UMi b 11 UMi b Radial Velocity \n", "2 3 14 And b 14 And b Radial Velocity \n", "3 4 14 Her b 14 Her b Radial Velocity \n", "4 5 16 Cyg B b 16 Cyg B b Radial Velocity \n", "... ... ... ... ... ... \n", "4159 4160 tau Gem b tau Gem b Radial Velocity \n", "4160 4161 ups And b ups And b Radial Velocity \n", "4161 4162 ups And c ups And c Radial Velocity \n", "4162 4163 ups And d ups And d Radial Velocity \n", "4163 4164 xi Aql b xi Aql b Radial Velocity \n", "\n", " pl_controvflag pl_pnum pl_orbper pl_orbpererr1 pl_orbpererr2 ... \\\n", "0 0 1 326.030000 0.320000 -0.320000 ... \n", "1 0 1 516.219970 3.200000 -3.200000 ... \n", "2 0 1 185.840000 0.230000 -0.230000 ... \n", "3 0 1 1773.400020 2.500000 -2.500000 ... \n", "4 0 1 798.500000 1.000000 -1.000000 ... \n", "... ... ... ... ... ... ... \n", "4159 0 1 305.500000 0.100000 -0.100000 ... \n", "4160 0 3 4.617033 0.000023 -0.000023 ... \n", "4161 0 3 241.258000 0.064000 -0.064000 ... \n", "4162 0 3 1276.460000 0.570000 -0.570000 ... \n", "4163 0 1 136.750000 0.250000 -0.250000 ... \n", "\n", " st_bmy st_bmyerr st_bmylim st_m1 st_m1err st_m1lim st_c1 \\\n", "0 NaN NaN NaN NaN NaN NaN NaN \n", "1 NaN NaN NaN NaN NaN NaN NaN \n", "2 NaN NaN NaN NaN NaN NaN NaN \n", "3 0.537 0.001 0.0 0.366 0.002 0.0 0.438 \n", "4 0.418 0.003 0.0 0.222 0.003 0.0 0.351 \n", "... ... ... ... ... ... ... ... \n", "4159 NaN NaN NaN NaN NaN NaN NaN \n", "4160 NaN NaN NaN NaN NaN NaN NaN \n", "4161 NaN NaN NaN NaN NaN NaN NaN \n", "4162 NaN NaN NaN NaN NaN NaN NaN \n", "4163 NaN NaN NaN NaN NaN NaN NaN \n", "\n", " st_c1err st_c1lim st_colorn \n", "0 NaN NaN 7.0 \n", "1 NaN NaN 5.0 \n", "2 NaN NaN 7.0 \n", "3 0.006 0.0 9.0 \n", "4 0.003 0.0 17.0 \n", "... ... ... ... \n", "4159 NaN NaN 7.0 \n", "4160 NaN NaN 8.0 \n", "4161 NaN NaN 8.0 \n", "4162 NaN NaN 8.0 \n", "4163 NaN NaN 9.0 \n", "\n", "[4164 rows x 356 columns]" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "planets" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/opt/anaconda3/envs/myNewEnv2/lib/python3.7/site-packages/numpy/lib/histograms.py:839: RuntimeWarning: invalid value encountered in greater_equal\n", " keep = (tmp_a >= first_edge)\n", "/opt/anaconda3/envs/myNewEnv2/lib/python3.7/site-packages/numpy/lib/histograms.py:840: RuntimeWarning: invalid value encountered in less_equal\n", " keep &= (tmp_a <= last_edge)\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAQ+UlEQVR4nO3df6xfd13H8eeLlk1+6Tp22zRtZ4upwGbcgGudogaY2MIMnQlLigoNaVKNk0BiIh1/aIxpUv4xaHSSZiA1Ik3lh6tM0aY40QArdzAYbam7bKO9aV0vQ0QgGWl5+8c9xO/ae3dPe+/33t4Pz0dyc875nM/5nvent3nd00/POTdVhSSpLc9a7AIkSfPPcJekBhnuktQgw12SGmS4S1KDli92AQDXXXddrV+/frHLkKQl5cEHH/x6VY1Mt++KCPf169czNja22GVI0pKS5Gsz7XNaRpIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGnRFPKE6V+t33bco5318z22Lcl5Jmo1X7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDZo13JO8OMlDA1/fSvKOJNcmOZTkkW65YuCYu5KMJzmRZPNwhyBJutCs4V5VJ6rq5qq6GXgF8F3gY8Au4HBVbQQOd9skuQHYBtwIbAHuTrJsSPVLkqZxqdMytwJfraqvAVuBfV37PuD2bn0rsL+qnqqqx4BxYNN8FCtJ6udSw30b8KFufVVVnQHoliu79jXAqYFjJrq2p0myM8lYkrHJyclLLEOS9Ex6h3uSq4A3AH83W9dp2uqihqq9VTVaVaMjI9P+fldJ0mW6lCv31wGfr6onuu0nkqwG6JZnu/YJYN3AcWuB03MtVJLU36WE+5v4/ykZgIPA9m59O3DvQPu2JFcn2QBsBI7MtVBJUn+9XhyW5LnAa4HfGmjeAxxIsgM4CdwBUFVHkxwAjgHngDur6vy8Vi1Jeka9wr2qvgu88IK2J5m6e2a6/ruB3XOuTpJ0WXxCVZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWpQr3BPck2SDyf5SpLjSX4uybVJDiV5pFuuGOh/V5LxJCeSbB5e+ZKk6fS9cv9T4BNV9RLgJuA4sAs4XFUbgcPdNkluALYBNwJbgLuTLJvvwiVJM5s13JP8KPBLwPsAqup7VfVNYCuwr+u2D7i9W98K7K+qp6rqMWAc2DTfhUuSZtbnyv1FwCTwV0m+kOSeJM8DVlXVGYBuubLrvwY4NXD8RNf2NEl2JhlLMjY5OTmnQUiSnq5PuC8HXg78ZVW9DPgO3RTMDDJNW13UULW3qkaranRkZKRXsZKkfvqE+wQwUVUPdNsfZirsn0iyGqBbnh3ov27g+LXA6fkpV5LUx6zhXlX/BZxK8uKu6VbgGHAQ2N61bQfu7dYPAtuSXJ1kA7ARODKvVUuSntHynv3eBnwwyVXAo8BbmfrBcCDJDuAkcAdAVR1NcoCpHwDngDur6vy8Vy5JmlGvcK+qh4DRaXbdOkP/3cDuOdQlSZoDn1CVpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGtQr3JM8nuThJA8lGevark1yKMkj3XLFQP+7kownOZFk87CKlyRN71Ku3F9dVTdX1Q9+UfYu4HBVbQQOd9skuQHYBtwIbAHuTrJsHmuWJM1iLtMyW4F93fo+4PaB9v1V9VRVPQaMA5vmcB5J0iXqG+4F/EuSB5Ps7NpWVdUZgG65smtfA5waOHaia3uaJDuTjCUZm5ycvLzqJUnTWt6z3yur6nSSlcChJF95hr6Zpq0uaqjaC+wFGB0dvWi/JOny9bpyr6rT3fIs8DGmplmeSLIaoFue7bpPAOsGDl8LnJ6vgiVJs5s13JM8L8kLfrAO/ArwZeAgsL3rth24t1s/CGxLcnWSDcBG4Mh8Fy5JmlmfaZlVwMeS/KD/31bVJ5J8DjiQZAdwErgDoKqOJjkAHAPOAXdW1fmhVC9Jmtas4V5VjwI3TdP+JHDrDMfsBnbPuTpJ0mXxCVVJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQb3DPcmyJF9I8vFu+9okh5I80i1XDPS9K8l4khNJNg+jcEnSzC7lyv3twPGB7V3A4araCBzutklyA7ANuBHYAtydZNn8lCtJ6qNXuCdZC9wG3DPQvBXY163vA24faN9fVU9V1WPAOLBpfsqVJPXR98r9PcDvA98faFtVVWcAuuXKrn0NcGqg30TXJklaILOGe5JfBc5W1YM9PzPTtNU0n7szyViSscnJyZ4fLUnqo8+V+yuBNyR5HNgPvCbJ3wBPJFkN0C3Pdv0ngHUDx68FTl/4oVW1t6pGq2p0ZGRkDkOQJF1o1nCvqruqam1VrWfqP0o/WVW/CRwEtnfdtgP3dusHgW1Jrk6yAdgIHJn3yiVJM1o+h2P3AAeS7ABOAncAVNXRJAeAY8A54M6qOj/nSiVJvV1SuFfV/cD93fqTwK0z9NsN7J5jbZKky+QTqpLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJatCs4Z7kR5IcSfLFJEeT/FHXfm2SQ0ke6ZYrBo65K8l4khNJNg9zAJKki/W5cn8KeE1V3QTcDGxJcguwCzhcVRuBw902SW4AtgE3AluAu5MsG0bxkqTpzRruNeXb3eazu68CtgL7uvZ9wO3d+lZgf1U9VVWPAePApnmtWpL0jHrNuSdZluQh4CxwqKoeAFZV1RmAbrmy674GODVw+ETXduFn7kwylmRscnJyLmOQJF2gV7hX1fmquhlYC2xK8lPP0D3TfcQ0n7m3qkaranRkZKRftZKkXi7pbpmq+iZwP1Nz6U8kWQ3QLc923SaAdQOHrQVOz7lSSVJvfe6WGUlyTbf+HOCXga8AB4HtXbftwL3d+kFgW5Krk2wANgJH5rtwSdLMlvfosxrY193x8izgQFV9PMlngANJdgAngTsAqupokgPAMeAccGdVnR9O+ZKk6aTqounwBTc6OlpjY2OXffz6XffNYzVXvsf33LbYJUi6AiR5sKpGp9vnE6qS1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSg2YN9yTrkvxrkuNJjiZ5e9d+bZJDSR7plisGjrkryXiSE0k2D3MAkqSL9blyPwf8XlW9FLgFuDPJDcAu4HBVbQQOd9t0+7YBNwJbgLuTLBtG8ZKk6c0a7lV1pqo+363/L3AcWANsBfZ13fYBt3frW4H9VfVUVT0GjAOb5rtwSdLMLmnOPcl64GXAA8CqqjoDUz8AgJVdtzXAqYHDJrq2Cz9rZ5KxJGOTk5OXXrkkaUa9wz3J84GPAO+oqm89U9dp2uqihqq9VTVaVaMjIyN9y5Ak9dAr3JM8m6lg/2BVfbRrfiLJ6m7/auBs1z4BrBs4fC1wen7KlST10edumQDvA45X1Z8M7DoIbO/WtwP3DrRvS3J1kg3ARuDI/JUsSZrN8h59Xgm8GXg4yUNd27uAPcCBJDuAk8AdAFV1NMkB4BhTd9rcWVXn571ySdKMZg33qvoPpp9HB7h1hmN2A7vnUJckaQ76XLnrCrN+132Ldu7H99y2aOeW1J+vH5CkBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1KBZwz3J+5OcTfLlgbZrkxxK8ki3XDGw764k40lOJNk8rMIlSTPrc+X+AWDLBW27gMNVtRE43G2T5AZgG3Bjd8zdSZbNW7WSpF5mDfeq+hTwjQuatwL7uvV9wO0D7fur6qmqegwYBzbNU62SpJ4ud859VVWdAeiWK7v2NcCpgX4TXZskaQEtn+fPyzRtNW3HZCewE+D666+f5zI0LOt33bco5318z22Lcl5pqbrcK/cnkqwG6JZnu/YJYN1Av7XA6ek+oKr2VtVoVY2OjIxcZhmSpOlcbrgfBLZ369uBewfatyW5OskGYCNwZG4lSpIu1azTMkk+BLwKuC7JBPCHwB7gQJIdwEngDoCqOprkAHAMOAfcWVXnh1S7JGkGs4Z7Vb1phl23ztB/N7B7LkVJkubGJ1QlqUGGuyQ1yHCXpAYZ7pLUIMNdkho030+oSkOxWE/Ggk/Hamnyyl2SGmS4S1KDDHdJapBz7tIsfBOmliKv3CWpQYa7JDXIcJekBhnuktQgw12SGuTdMtIVyqdyNRdeuUtSg7xyl3QR7+1f+gx3SVcMf6jMH6dlJKlBQwv3JFuSnEgynmTXsM4jSbrYUKZlkiwD/gJ4LTABfC7Jwao6NozzSdJctHhn0rCu3DcB41X1aFV9D9gPbB3SuSRJFxjWf6iuAU4NbE8APzvYIclOYGe3+e0kJ+ZwvuuAr8/h+KXO8Tt+x79E5d1zOvzHZ9oxrHDPNG31tI2qvcDeeTlZMlZVo/PxWUuR43f8jv+Hd/wzGda0zASwbmB7LXB6SOeSJF1gWOH+OWBjkg1JrgK2AQeHdC5J0gWGMi1TVeeS/C7wz8Ay4P1VdXQY5+rMy/TOEub4f7g5fl0kVTV7L0nSkuITqpLUIMNdkhq0ZMJ9ttcZZMqfdfu/lOTli1HnsPQY/2904/5Skk8nuWkx6hyWvq+zSPIzSc4neeNC1jdsfcaf5FVJHkpyNMm/LXSNw9Tj7/+PJfmHJF/sxv/WxajzilJVV/wXU/8p+1XgRcBVwBeBGy7o83rgn5i6x/4W4IHFrnuBx//zwIpu/XU/bOMf6PdJ4B+BNy523Qv8/b8GOAZc322vXOy6F3j87wLe3a2PAN8Arlrs2hfza6lcufd5ncFW4K9rymeBa5KsXuhCh2TW8VfVp6vqv7vNzzL1bEEr+r7O4m3AR4CzC1ncAugz/l8HPlpVJwGqqqU/gz7jL+AFSQI8n6lwP7ewZV5Zlkq4T/c6gzWX0WeputSx7WDqXzGtmHX8SdYAvwa8dwHrWih9vv8/CaxIcn+SB5O8ZcGqG74+4/9z4KVMPSz5MPD2qvr+wpR3ZVoqv6xj1tcZ9OyzVPUeW5JXMxXuvzDUihZWn/G/B3hnVZ2funhrSp/xLwdeAdwKPAf4TJLPVtV/Dru4BdBn/JuBh4DXAD8BHEry71X1rWEXd6VaKuHe53UGLb/yoNfYkvw0cA/wuqp6coFqWwh9xj8K7O+C/Trg9UnOVdXfL0yJQ9X37//Xq+o7wHeSfAq4CWgh3PuM/63AnpqadB9P8hjwEuDIwpR45Vkq0zJ9XmdwEHhLd9fMLcD/VNWZhS50SGYdf5LrgY8Cb27kam3QrOOvqg1Vtb6q1gMfBn6nkWCHfn//7wV+McnyJM9l6i2sxxe4zmHpM/6TTP2rhSSrgBcDjy5olVeYJXHlXjO8ziDJb3f738vUHRKvB8aB7zL1k7wJPcf/B8ALgbu7q9dz1cib8nqOv1l9xl9Vx5N8AvgS8H3gnqr68uJVPX96fv//GPhAkoeZmsZ5Z1Ut2dcAzwdfPyBJDVoq0zKSpEtguEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QG/R8fe7GnCvOm5gAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.hist(planets['pl_orbeccen'])\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "from hermite_library import read_hermite_solution_from_file" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "planet_file = 'data/Kepler-11-savedSim.txt'\n", "#planet_file = 'Kepler-11-savedSim.txt' # if all .txt files are in the same directory as my notebook file\n", "\n", "t_h, E_h, r_h, v_h = read_hermite_solution_from_file(planet_file)\n", "#time, energy, postion, velocity" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(7, 3, 8800)" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r_h.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The shape of r_h (position) array is: `r_h[# of planets+host star, x/y/z - 3D position, time steps]`" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(3, 8800)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r_h[0,:,:].shape # all positions, all times of the first planet" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(8800,)" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# x position, all times for the 2nd body in this system\n", "r_h[1,0,:].shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First thing we are going to do is plot the orbits statically (not a movie).\n", "\n", "2D plot of x/y positions over ALL timesteps for ALL of the planets." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "7" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r_h.shape[0]" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAASIAAAEKCAYAAAC2QVjtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2dd3hURduH70lP6L1D6AkQahBBlK4INiwoooCigKJYUfS1oCiifvIiCiJWioDltVAFRJqA9BIgobfQCZCQhPT5/jhh92yS7bs5u2Hu69orp8zMmQPZX6Y8RUgpUSgUCiMJMLoDCoVCoYRIoVAYjhIihUJhOEqIFAqF4SghUigUhqOESKFQGI6hQiSE6C2E2CeEOCiEGGOjXHshRK4Q4v7i7J9CoSgeDBMiIUQgMAW4HWgGDBBCNLNS7kNgafH2UKFQFBdGjohuAA5KKQ9LKbOAecDdRZR7FvgfcK44O6dQKIqPIAOfXQs4oTtPBDroCwghagH9gO5Ae1uNCSGGAcMASpUq1S4qKsqjnVUoFPbZunXrBSllFWfrGSlEoohrBf1NJgGvSilzhSiquK6ilNOB6QCxsbFyy5YtHumkQqFwHCHEMVfqGSlEiUAd3Xlt4FSBMrHAvHwRqgz0EULkSCl/L54uKhSK4sBIIdoMNBZC1AdOAg8BD+sLSCnrXzsWQnwPLFQipFCUPAwTIilljhDiGbTdsEDgWynlHiHEiPz704zqm0KhKF6MHBEhpVwMLC5wrUgBklIOKY4+KRSK4kdZVisUCsNRQqRQKAxHCZFCoTAcJUQKhcJwlBApFArDUUKkUCgMRwmRQqEwHCVECoXCcJQQKRQKw1FCpFAoDEcJkUKhMBwlRAqFwnCUECkUCsNRQqRQKAxHCZFCoTAcJUQKhcJwlBApFArDUUKkUCgMx6dTTgsh7hZC7BJC7BBCbBFCdDainwqFwrsYFrNal3K6F1pqoc1CiPlSyr26YiuA+VJKKYRoCfwEqMyJCkUJw6dTTkspU6WU15IulqJwAkaFQlEC8OmU0wBCiH7AB0BVoG/xdE3hT0gpST53lrOHD3Dm0AEunkokNSmJ1EtJpCdftlk3JDyciHLlKV2hEhVr1qZq/QZUjWxI5br1CA4NK6Y3UPh6ymmklL8BvwkhbgHGAT2LbEyIYcAwgLp163qwmwpfIenkCRLWrebQ5n85f/yoR9rMunqVrKtXuXzmNInxu62WCwmPoEHb9kR37kpk67YEBAR65PkKDV9POW1CSrlGCNFQCFFZSnmhiPvTgekAsbGxagrn51y9ksK2xX+wbcl8sq5etVu+TKUqVGvQiGoNGlGpdh3KVKpC6YqViChbjoDAokVDSknW1XTSLl8i9WISSYnHOXf0COeOHuLckUMWZbOuppOwbjUJ61ZbXG8Y24Eb+z1I9UZNXH9ZBcK8BFPMDxYiCNgP9EBLOb0ZeFhKuUdXphFwKH+xui2wAKgt7XQ6NjZWbtmyxXudV3ic7IwMtiz8jfU//2C1TES58jTv0oNG7W+keqMmxToqSb2YxKGtm9i7diWn9u21Wq5+63Z0efQJKtWuY7VMSUYIsVVKGet0PaOECEAI0QeYhDnl9Pv6lNNCiFeBQUA2cBUYLaX8x167Soj8g5Tz51jx3TQOb91U5P3Wt91Bu773UL5a9WLumWPkZGWxd+3f/Pu/H7mSdL7Q/ZDwCG4d/ixNbuyMEEWtRJQ8/FKIvIUSIt8lPSWZpV9M4vC2zYXuNe/Sg079B1K2clUDeuY+uTk57F65jL++nlroXmBQEPe8+jaRLdsY0LPiQwmRDiVEvoXMy+PfX38sctrV5dGhtOl9J4FBRi5XeofTB/ex7MvPuFBgYb1q/Yb0e/VtSleoaEzHvIgSIh1KiHyD1ItJ/PL+myQlHre43umBgXTo19/qInJJJDFhD398/B4ZqVcsrvcdNZqom7oY1CvPo4RIhxIiY0mM382PYy09dmo3a8FdL75OeJmyBvXKN5BSsm3xfFbN/Mrietvb76LroCcQAf7t/qmESIcSImM4tHUjv380zuJatyHDadP7jutmsdYZLp5KZN5br3D1SorpWuMbOtH3udEEBgUb2DPXUUKkQwlR8XJk+xZ+nTDW4tqAcR9Ts0m0MR3yM7IyrvLrB29zMsFsFtC04830HTXa70ZISoh0KCEqHs4fP8rM0c9YXBv88edUrhtpTIf8nNycHBZ//gn7N6w1XevQrz+dHxpkYK+cQwmRDiVE3iUzPY1pwweRk5Vpuvboh5OpGtnAwF6VHHJzsvnlvTctXE76jXmbBm3aG9grx1BCpEMJkfdYNfNrti763XR+zytv0bDdDQb2qOSSmZ7Gt88Pt3DcfeqrH4goW87AXtlGCZEOJUSe5+zhg8x+7XnTefu77+eWh4cY16HriLNHDjF7zHOm8za976T7Y8MN7JF1lBDpUELkOaSU/DphLEd3bAUgMDiYp6b/QGhEhPeemSeRQECA2mnT8+//5rHup9mm8yenfOtzVuhKiHQoIfIMSYnH+f6lp03n7qxTXD6bzo4VJ9iz5qSnukdQSACtetQhpmttSpUL9Vi7vkxmejqfP9bfdN7ujn50fXSogT2yRAmRDiVE7rNq5ldsXfQHAOFlyjJ82gyHbVuuXMzgr+/2cuqA7aBk3qJh2yp0HRhFWCn/tMVxhB1LF7Hi2y9M56Nm/EJwmPGB3JQQ6VBC5Dq5OdlMGtjPdN5n1Gii7bgg5ObksWJGPAc2n7VZLjAogJbdatPs5pqUr+r61E5KycVTacStPunQCKtNr7p07NcQUcKmehmpqUwZ+pDpvP/bH1CnWYyBPVJCZIESItc4f+wIM1951nT+9NdzrLpkZGfm8vvEbZw7dqXI+2Urh3Hr0BZUq198Lh1SSo7tTmLpV7vJycorskyjdlXpNbR5iVp/WvjpR+xbvwaAlj170+vJZ+zU8B5KiHQoIXKeuL+XsezLyYDmF/bg2xOKLLf1z6P8+/vhIu/d90o7qjfwra3lIzvPs/iLuCLv9RrajCbtfTPWkbMc27WDX95/A4DQUqUY+c08Q9xqlBDpUELkHH9O/S97Vq8A4LYRz9GiWy+L+1kZOXz1/JpC9QKDAnhkXEdKV/CPheKkk6nMG1c4CFv5ahEMeOsGAgL9y52iIFdTrzB16ADT+TPf/UhoRKli7YMSIh1KiBxDSsm04Y+aDOYKWkenXLjKrDc2FKp313OtqRPt37F0EjacZsWM+ELXn5h4M6ER/rvILfPymDzofnKyswB4fNKXVKhRq9ier4RIhxIi++Tl5fLfAeY0ck9/M5fw0mUASL2UwYzX1luUDy8TzODxNxEY7N+jhoJkpGXzzUtrC11/ctIthIT5b7A2/brRgHH/R80mxZOXVAmRDiVEtsnJzubTR8w7Yy/M+YOAwECyM3OZ/pxllooGbarQe1iLEh/GIy9PMuftf0k+b84YIgIET33e1W9329b/PIcNv8wB4J5X3qRhu0JpAz2OXwqREKI38Cla8PyvpZQTCtwfCLyaf5oKPCWl3GmvXSVE1snOyGDy4PtN5y/OW4AQgj8mbScx4ZLpet3mFbnz2dZGdNFQZJ7ku1f/4eqVbNO11r3qctN9jQzslevoNyEcMcVwF78TIiFEIFo6oV5oOc42AwOklHt1ZToB8VLKS0KI24GxUkq7sq6EqGiyszKZ/Oh9AJSpXIVhU77jRMJF5k/aYSoTFBzAsE+7+O0owFPkZOXy5SjL0eHAd26kfDXvubZ4i0NbN/H7R+8CcOeLr9Gkw01ee5Y/ClFHNGG5Lf/8NQAp5QdWylcAdksp7a68KSEqTG5ODpMG3gNA+Wo1ePzT6Ux9aqVFmUHjO1GmovHWub7E2aMp/DLB/LtUtV4ZHnjN98NxFEQvRt6cprkqREauPNYCTujOE/OvWWMosMTaTSHEMCHEFiHElvPnC+eYup6ReXkmEQorVZq+z/+fhQh1uKs+I6d1VyJUBNUiyzJyWncatqkCwLljV5gy4m9SL2UY3DPnaNjuBu566XUAfv9oHCf2Fm1bZRRGClFRY/8ih2dCiG5oQvRqUfdBSzktpYyVUsZWqVLFQ10sGUwccJfpOLLda/zvo62m82GTuxDbp74R3fIreg+PYfAH5inNjNfWs2nhEQN75DyNb+hEn1GjAfjpnde4dMZqhvdix+enZkKIlsBvwO1Syv2OtK2mZmZmv/Y8Zw8fBCC0/Aum3a/mt9Si68NNi70/mYePcObdd0n/91+X2yhze2+qv/EGQZUqebBnjjNv3CaSTqaazkdO625IP1zl319/ZN2PswAY+c08wkqX9ljb/rhGFIS2WN0DOIm2WP2wlHKPrkxd4G9gkJRyfZENFYESIo2V309n25L5AISWH4X2Tw4PvXkDlWp57pfPGrmXL3Ok/4NkHz9uv7CblO7aldqTP0WEhHj9WQBHdl1g8dRdpvPhk7sQFOI/edoWTPyA/RvXAfDC3D8ICPBM3/1OiACEEH2ASWjb999KKd8XQowAkFJOE0J8DdwHHMuvkuPISyohggOb1jP/k/EAhJYbjgjQTP2fmtrNqw6fyQsWcmr0aLvlImJjqTzyaSI6dHAoU4XMyiJl+XKSpk0j88BBu+Uj580lvLV3zQ+yrubw1Qtm15dHxt1IuSr+s6v2xbBHSE++TFBwCM/N/tUjbfqlEHmL612IUi8m8eVTgwEIKfMQAUE1CS0VxBOf3OKV52XEx3Ok371W79d4/z3K33efV54tpeT8xP+S9NVXVss0XvePV6dxU0b8bTq+89lW1G1uzJTRWaSUTHzoTgBa39aXHo8/5XabSoh0XM9CJPPyTIvTgWE3EhzeiehONeg+yPM5xk6+9DIpixYVul6qc2fqTP/SkJxceenpHOzWndzk5EL3qr3+GhUHeSc1z6w31pNyQdtJ6z4omuhONbzyHE9z9UoKU594GIB7X3uH+q3budWeEiId17MQTR/5GFcunAcRQVj5Edz8YGNadqvjsfallOy/sSN5Bb7oQVWq0GjVSoQP5bPPvXKF/e0LZxgp2+d2ak2c6PHnLftmjyk43C0PNSGma22PP8MbnNizi5/e1bb2n/3+J0LCXZ9e+qMdkcLD7Fm9QhMhtHWhno8186gIHb77HhKim1mIUK3PJhOdEE/jtWt8SoQAAsuUITohnuiEeCqPHGm6nrJ4CfFR0Zz9+GOPPu/Woc1p2V0TnzXz9hO3KtGj7XuLOs1bmkK/fDakv53S3kEJUQkhIy2VP6f+F4CQskO4fXgMTTt4JujXpblziY+KJnPfPtO1BosWEp0QT9levWzU9B2qPPsM0Qnx1J46xXTt4jffEh8VTfr27R57zs39m9D2tnqAJkaHtp3zWNve5LYR5nRFa+fOKPbnKyEqIUx5XItdHBjahl5Db6JhW/fTzOSlpREfFc2Zd941XYv8+SeiE+IJbdjQ7faNoEz37kQnxFN9nPmdjg14mPioaDy1TNGxX0Na3KI5Cfw5fTdnj6R4pF1v89RXPwCw6fefST53plifrYSoBPD39z+Yjjv1f5yoG91fKD391tvsa2ee6lcc+jjRCfGEx3guOHtOXg5/HPyDe+ffS8yMmEKfO367gxl7ZpCVm+WxZ16jwgMPaIIaZY7TkxDdjOSFhRffXaHLw02p2bg8AL98uIX0FM+/g6eJKFuOnk9oU9ivn32iWJ+tFqv9nLTLqUwbro2GGt34Gne/4J5ntZSShOhmFtei4ve6HY9o/MbxzE2Y61YbejrX6syUHlMIEO7/Lc1NTWN/rOX6anRC4eiNrqDf2n/6i25+EdfpkwfvAKBT/4F0vG+AndKWqF0zHdeLEGl2IPcBWQSGtOT5WePdai9j/36O3GWO2lh97FgqPPSgS20dvnyYu/+4235BDzGp2yR61O3hVhuJo57jyrJlpvOmW7cQUMr9mM96MfIHd5D0lGS+eHIg4HzcayVEOq4XIZo8dAbZqT8D5gBnrnJ+8mQuTDUn7IvatdNpd4ms3CzazbZuhxIUEMTCfgupVdr1GMoXrl5g4KKBnEqz7rD51/1/Ua1UNZfazz55koM9eprO637/HaVuvNGltq6Rl5vHFyNXAVCzcXn6vdTWrfaKg1WzvmHrwt8IDg1j1MxfHK6nhEjH9SBEW/88yqrvtPxVD7z5IXVbNHe5rUN9+pJ12JwiyNlpyabTmxi6rOi0xyv7r6RyeGWX+2aP9Ox0OswpOrbOs22eZVjLYU63WXB6WnnkSKo8616usMtn0/nhbc3R954X2lCraQW32isOrk3RHpnwKdXqO7Y5oYRIR0kXovSULL4aNZmcq6sAeOnHhS63FR9ltrgu378/Nd59x+G6W85s4bGljxW6vvjexdQp4zn7JUe5nHGZm3+8udD1F9u9yGMtCvfTHgmtWiMzMwHNqbbOtC/s1LDN+l8Psn2Z5gDsbZ8/T3Bkx1Z+/eBtwPHfMWXQeB3x7ei1JhEa8eUsl9vRi1CNCR84LELp2enEzIgpJEK7Bu0ibnCcISIEUD6sPHGD44gbHEdYoDnI28StE4mZEcPh5KITQ1ojaucOytx2GwCpq1Zx+F7r/nSO0Olec9zrL55eaaOkb6B394hf693+KiHyM5Z9s4ec9OUAVGvQiFLlXRvi60Wo3qyZlL/nHofqxcyIKTQVuvbl96Udoc2PbCZucBxBwpwS6O7f7yZmhnPmB7U/nUTlpzVn0My98RwbNNitfj01tZvpOGHDabfaKg6GTv4agMWff+LV5ygh8iMy0rLZv+kMuVm7AXjwnQ9dasdChObMIaK9/RjMSVeTCn2Jdzy6g7jBvhVytCDbB21n16BdFtdiZsSw+cxmh9uoMmoUlfPXiNI3beLUG2+43J+AAEGvodr604oZ8cg8314aKV+tOiLfRGL3yuVee44SIj/im5fMU7KaTZsRHOJ8qucTw0eYjut8NZ2Itm3s1nllzSt0/amr6fyz7p8RNziOQA8F0/I2QgjiBsfx850/m649vvRxp0ZHVUaOpMLDmpd68i//I3nBApf706S92fVmqh9M0Z74TBsVLZ32qdeeoYTITzi0/RxSSnIzNb+o+19/106NwiQvWEDqai1FTtUxr1L65sILuwWJmRHDkiPmnAW7Bu2ia52uTj/bF4iqGFVoBOeMGFV/601Cm2rhdU+NfoWsxJMu92X4ZHN+sZQLV22UNJ6yVczuQvHrVtso6TpKiPyEP7/cTW6mthNYsWZtgsOcy7iRc+kSp0a/AkBYTAyVhgyxW0f/Je1Vr5fPrQO5StzgOJ6MedJ0HjMjhuzcbBs1zDT443fT8aGePV32TwsKCaRBay3Jw6w3NrjURnHy+KfTAVg82bMRC66hhMgPWPOjljMg56qWo92VtaEDHTuZjuv//JPNslJKCxH6+tavmdjV8/F7jGRU21EsuMc8vWo7uy2pWak2apjR21kVdIdxhttHmP+Nj8ZdcLmd4qBC9Zqm46STJ2yUdA1DhUgI0VsIsU8IcVAIMaaI+1FCiA1CiEwhxMtG9NEXiFuZSF7OWdN5RNlyTtXXL07bM1bMk3m0nNnSdP7PQ//QoYb3c6YbQWS5SLY/ag4B0nFuRy5mXHSobtPt20zHJ/NHmq7QZUATABZN2WWnpPHc8bz2Ff3h9Rc93rZhQpSfcnoKcDvQDBgghCj45+UiMAr4v2Luns+wYqYmHFlXNA/7+5xcG7qyapXpuOFf9nc9Ws1sZTrePHAz5UKdEz1/IyggyGJXrcuPXRyapgWEh1Pj/fcBSFmwgJxLl1x6fosu5iiOh3f4dmLQph07A5CdcZXcHMemso5i5IjoBuCglPKwlDILmAdYeElKKc9JKTcDnn1rPyJh/WmkzDOdR7Zyzk8pcYRmA1OqU0dCatsOXaqfjm17ZBthQV7O/JqTCWPLOf5J8U5CQCGEhRi1nd3WobWf8veZDRz1U19n6faoFopkyTTfNoUAaNJBi+7wzzzXDWmLwp9STtukJKac3rxIyySak6EtZjZq75zzpX5KVvfbb22W1YvQmgfXEBwY7NSzHEYvLO85GbxtYrRlfQ+6JwkhLKZp+umpLaLi95qOz7zvWvSDZjeZ1198PYjarSNGAbBlgWfSD13DL1JOO0JJTDm9aYEmRLkZGwG4dfgoh+vmXDSvdTT6e4XNsi+tesl0PKfPHCqEedghM36BWTxs8dR6GJts/rwYDxVspMN+p7zW5m8jrJdxgqCAIFb2N9v1OLK1L4Sger5rzKVZs1zeRWvbWwsv+8uHvu0jqQ8JknoxyWPtGilEiYDeKak24DvJuA3mWkpjKTNN18LLlHW4/oFO5gBpwTVrWi2XeCWRZce0GDy96vUipornIjCyfbYmFD8+Uvje25ctRWdsMlQrEEGgbE14bkfhcgWnjDvnas/5yv1YP5XDK/NyrHlf5Kd9tncYASr0Nwecd3UXreM9Zu/2nOxcl9ooLro88jgAf34xyWNtGilEm4HGQoj6QogQ4CFgvoH98SnmjdsEQOmy2l/Itrff5XDdq7vM6x36qUNR3P7r7aZjj23R5+VqwvDHSMvrb5wzi4k79khvnDW3o+fkVu255/e73jYwuLnZn2zcv+McGuU0WmUeSeVluRcWdu67m9yq723a9tWWco/t8lzSgSD7RbyDlDJHCPEMsBRzyuk9BVJOVwe2AGWBPCHE80AzKaVvT6TdRP+Ln3Rcm5bd9NCjDtc/2l+LqhjatKlNA0T91KOgP5bLfFgfrhbYAi8oGMDxpHRu+bho94aBHeryfj8HR2bX2tZP+6a0t/pcR4kbHGf692k5s6Vdn7rg6ma3jX0tW7kUavaJiTfz9YtrSTnv25bWATrXnozUVMJKl3a7TcOECEBKuRhYXODaNN3xGbQp23XF1iVHASx2y0LCwh2qm3nQnBdebwVckJQss5Z/0fMLz1hMF1wDemEvlNP2H75ac5j3Fzv25fxh43F+2Hjc4trRCX3tPDtZW7x+p7xlf9wQo8X3LqbPr30AWHF8hd1QtE22bDHFvpZ5eU5nug2NMG8QJJ1MpVIt97/g3qLdHf3YuvA31s79nl5Puhc0DpRltU+ycb62SN2802UA6rZwbAcH4PAdWi7zwPLlbZa7aa55Dalzrc7OdrEwBUVobDKUq0XkmEVEjlnksAhZ41o7f+62ETpDCO25YQXEKPOKS8/Ux1V6fuXzdssHljYv5CY0cy1iZpMOWojba1NzX6XTA5oD8K6//vRIe4aOiBSF0U/L9qzSvMW7DnrSWnHLutlmc6vGG9ZbLbfu5DrT8Y5HdzjbxcLoRajzC9BzLJFjik7LM+nB1tzTxjErDSkl9V+zGDAzYrZm0fz3S11oUMXKiGHMMTi/D6bkp5v+oDa8fBBKO7+bqp+iPbviWT7r8ZnN8o3+XsHB7q4H8e8xuBn7N561X9BgHB2hO4oaEfkYO/4ym1ZlpqcBUKWejS1sHQkx5pGTranWiL/M291uh/LQi9DwtXwR9GghEWpUtTRHJ/Tl6IS+DosQaO9wrd7ed2+zuNf9k9VWxQ6AKk3hLZ218/81ghzXFpFvqqWNHlclrrJbVr9DeXHmTKefpQ8fm3w+3en6xUnpipUAOHPogNttKSHyMdb/T1vjuWWA65lUGy5bavXe9nPmnQ63F6j1IjRqB5GfnuDDPxNMlypEBHN0Ql/+erFLEZWdIyIkiKMT+rLvvd4W1yPHLCI7N6/oSgEBmpnANd5zzb5sWk/TsiVvrnvTbvlKT2oj2LPjP3DpedUbaGYa//toq0v1i4uO92vTsw2/zHG7LSVEPorM0dZUGrS1Hz0RtJxk1wipW9dquUFLBpmO3Vqg1ovQwF+I/MjSTODohL5sf+tW19u3QmhQIEcn9LXY/W/8nyVsOmLFWVUIeFPn2W7PqNIKzSpp9kG/H7S+AXCNqi+ZnUJdMXDs+7Tm73f1im97NjXvok1BD29zPNqlNZQQ+RAZaeZfvO1LtBAVsXc6FrBdnxjRGjl5OaZjt0ZDe34zH4eWJfIbyymP3R0uD3Dkg77sGmsWuv5fbmBFvJW1lcBgeFhnmOiCGP14x4+m4/gkxxfeT456zulnhZX2knuNhwkM8twSs1UhEkIsEELM133+EEJ8I4QowkxW4Qn++UmbawsBF04cA6B2dAun2mhiI41Sm1nmsLBujYZ+HmI6jEw2T1ta1i5XLCJ0jbJhwRbPGzpjC8lXrYwimliuMZFnZTrnAP0X9rdb5ppn/pXl7sV5PrnfNa/+4iY7M8Ot+rZGRP8HfKL7TAQWAv2EEBPceqqiSPZtPAPAbcPM4uOIYOSlmxc19VvI1ph1uxue07rRRGSGeW2gSbXSzH/GA2YALqAXo1bvLLNeUG9T9K7z/nTrBqyzXygfvWe+K9Rtri0EL/t6j1vteJsG7bSdyQMbre/SOoJVIZJSri7i8xvQH+htrZ7Cfeo1d+5LcuwR+1bXV7LMtjStq7Z2uk8FSavcyuJ82QvuL0i7w4H3za4qNnfTRurscy4ctF6uCMqGmH39ftnveBrmqzucN5HoOlCLjZ2e4p67iLdp0VVLz71nzd9uteP0GpGU0rc98koAJxO0v4I1mkQ5VD5jr7ZQXH3s21bLdJrrerwcE7rRUPPEV03HxTkds0ZwYAD3tTUb4U9ZaUVkqjQ1H3/erugyDvDOBscz4h59aIDT7Zep6OVYUB4isqUWH+t4nHv2aLbWiCoW8WkohHgH8O3xoh9y9Yr5L9+xOG2LvX4r574oFR56yG6ZL3q6lzYZ4Jsc8+jjyAd93G7PU3zS3zxK+3jpPusFX9Nl38jNsV6uCFb1X+Vw2QYLXU855C84m8TBGrZGRFvRHE636o7nARHAUx55usJE3Grzl+PEHm1Hq05z+46fzm4Pu+zOMdOcCXZcjnkq6GtZPfSjs44fWInDFKqzyB5Xyan2K4Wby+dJ2wveoY0a2bzvKMk+7gTrCWytEdWXUjbI/3ntuL2UcjRQ8v9lipmdKzSL6uibapgsVWs0bmqrCgBXbeySXePEFQ9kXThc2FPeF6ZktjidbGMn54bhbrf/4SbHs6m4GjANIG51ost1/QWH14iERnchxNdoQc0UHiTrqjZF0IcNDQyyb1A0IqQAACAASURBVE9y8uXRdsv0X2B/u9lR6mfM9lhb3kIvkOlZVqZefT4yHzs5PbvGnATHLYrT1ju/qxTdqQYA+zaccbpucRIQqNkTZWW4Pj6xK0RCiA5CiE+BY2iBy9YCjq2iKpymar0yTpXPOasZ8ZV/4H6rZVKztWiPDcu56DZy9B/Tocz/lfH10dA1mr1l3d3FhJPTs8ndJjvdj3MfO5+IpkEbzSVFb+jqi1Sqo1nyX0x0feRta7H6fSHEAWA8EAe0Ac5LKWdIKf3DysoPCQh0zdi90tChdst83uNzl9rme/8QneKiW91uDpcNj9U2HDITEuyULEytph6OHe4lKtfR4m1fSDxup6R1bP3WDwPOAl8As6WUSbgR3F7hXUIiI+2WqV3m+okx59CI7dWjXu9HxcGD7ReyQnCIm5ERiolKtbS4TUleEqLqwPvAXcBBIcQsIFwIoWIYXac0yZgBwNg7XU+zbARZOVZ2t8K9P+Io1bGj159hNKXKa/+OV1Ncj4Zpa9csV0q5REo5CGgE/AGsB04KIdz3+8ehlNNCCDE5//4uIYRz2QUVHiULbfF8yE2OxUfyFZq8scRrbdvbDQsoZd/lxt8JzX/HjLQ0l9twaEFCSpkhpfxFSnkf0Bgt4L1bOJhy+vb85zVGmyq6b42nUHiQ5EzbowBfs7PyBqERml1WZnqqy2244uKRIqWc4fITzdhNOZ1/PlNq/AuUF0LU8MCzfRbphle4v5M0N4Gkuc4v6noEF+18MnLd8zovCYRGRACQle7F7Xsv4kjKaYfTUpeUlNO5uZor3zXbjOuJnAtXkRmu2fS4jYtCFBRw/f0/FSQgPy6RO2mFfD3ltMNpqUtKyunAQG2nJM9FIzt/purI1lR+zLn4Sx7DydQ/1whSezcEhYQA5oiNLrXhSCEhRCcgUl9eSul8ZHBLHEk5fd2lpXY2F1ZJQgT433pK6RDbowB3XDv8hbyca6N4180NHLGsnoUWJK0z0D7/E+vyE804knJ6PjAof/fsRiBZSmkjsZXCmwSgrV/tSrxsp6Rvsek/rv+ltoe9qZnMKPlrSLk5muW3Iy5J1nBkRBSLlubZo9LuSMpptCywfYCDQDrwmCf7oHCO3aFDaZb5HXd9vs5vXDwAqpaxEqoiy/vpejIPHfb6M4wmI1XbLQst5foakSNCtBvNuNHjIxEHUk5LYKSnn1sSyU1OJrCc7aDwFzMuUjGsosvPiBCZLtctbv7zm+1c9QCM9/4GbNraNS7X9ZdpXUaaFv0zvIxzfpJ6HFmQqAzsFUIs1QfTd/mJCpvkWcvRZYeLM+xbVHyw0bU8W9zxX9fqGcgPG113N7DH3qS99gvlc/lnLaRsmV69nH5O6iVN+MPL+HZWj4wrmhCFlfauEI0F7kFzftUH01d4gUtnnJwu5G+dXphq39bzz6Mu5imPfbzQJZtxoX2In0c44GLhpM/ZkD+HOFw2+5S2t1LhEeeT35w6oK3F1Wxc3um6xcnVKykAhHlzRGQliP5ql5+osEnCBvMM2JGhec0P7I9yfrvrN7tlHOVo2MMea8tb6EWyfaSVqegOnZeSkz5nV3M0w7321R1LfgkQ0d75/Z1rqYRqNfFtL/yU8+cIL1OW4JBQl9uwFQbkn/yfV4QQKbrPFSFEistPVBRJTBfNTnPXqkQq1NCCo13LbWaLsnfYXzRuVMEzIUsLUv81/xgVFcnv7kc71qeitocrZhmJ8f4hRMnnz1KuajW32rDl9No5/2cZKWVZ3aeMlLKstXoK12jRRQvRkZcjicwPmn9ku/0wsM76Mu2+sNv5zoFFTrDeAVpKHl9cS9WPhqzu7OXpEtGMcS6Ylz5bbkhgiM2yucmue6MDXLmobf1XqBHhVjveJuX8WcpW8ZIQKYoX/S9b/db5QrTDvhDpSfv3X7tlBixyPrVNQaaFTDId+9Ja0b4zV+wXAnhXN10Lc+5vqiNZXq9x9sOP7BdyAF92nM3LzSXl/DnvjYgUxYv+l61WM83NIXGvc6OX40Osm1ktuMcDqW3eNhsyvlPJnEq52/+tcr9tD3DbJPNWudXRUKZOrB78welnHLikJTYY0WqE3bLJv/4KQO0vpjr9nFxrMZR8jEunT5Gbk0Ol2nXdakcJkQ9y6bRz2T0brbSfZTOyXKTp2OWsHjqxHJz2nen4yIU01h+64FqbHkI/Mjuoy/paiA90USqj73DqGfr0QSNbO27eVqab46Flr3Fo+zkA6kT79vrQ+eNHAKhSz70YVY64eDwjhPDtf40SQvUG2jRh+bfO5a8MrmE2zHMkjEifX91IiqhbK9LvoD381Ub2njJmD0MvQm/d0Ywga3G/3zdnSOGVI04/p82sNg6XzTzsnkX1tj+1jYrWvdwbaXib88eOEBAYSMVadewXtoEjI6LqwGYhxE/5ERV9d8Lq5/R8rDkAyeeuUr6aJi5nDzuXn93W9Gztg2td75ye5v1Mh0d7m6ePfSav5aM/izeekF6EAgMEj3e28pc55RRk6yIIRjhvYX5tRPRpt0/tlj3cR5salrr5ZqefA5B0UutrnWjXLeGLg/NHD1OxZm2Cgt0zunTEjugNtAiJ3wBDgANCiPFCCBdz0yisUa5KuOm4Zc/eAOxYtthacQtqf/4ZAOmbNlktUz7MbBh35293utJFjQe+Nx+vGs+RNzuYTqeuOlQsC9g7T1y2eE7p0CAOjbcy0pMSJkabz8c6v5u14JB5ja173e4O16s9xfnMKTLPvB3py3/3ZV4epw4kOJQI1B6OhoqVwJn8Tw5QAfhFCOGZbQFFIapEal/u3SuXOVS+TM+epuO8dOvW2e92eheAoylHXe8cWHyZxccNOTrecl0mcswiDp5zPXSoLSLHLOLuKevMXbmzGbvfuc16hXd0lslvuZYJ6/V/XgegVuki4/JZkL5tm+k4IMT2Fn9R7NuoJVT0dYvqCyeOkZmWRq2o5m635cga0SghxFbgI2AdECOlfApoB9zndg8UFkTlZ/dcNHW/y23sa9vO6r1+jc3TqieWPuHyMwCLXTTercDR92+1uN1z4moixyziUppzi+/WiByzqNBo69D4PraD+Y/VOQI/vdGlAGiLDpuf+ed99t1kjj08EIBKT7r27/v3LG16e8tDTVyqX1wkJmhrmbWj3Q9m56jT671SytuklD9LKbMBpJR5gHPbDgq7dH3YPMwVQvvvuXzWsZTDjdf9Y78QMDpWS1O98cxG9zy8hYDXdXHqxlXm6PP12Pdeb4tibcYtL1JEHGHhrlNW6x6d0JdAa8HUcrIsRej+76CqawmKx6zVEsxUL1Xdbtm8LLPoVnnxRZeed21qVqmW62E1ioPEvbspXakyZatUdbstu2FApJRv2bgX73YPFBYEBpn/NrTp059ti+ax/ucf6PPMS3brBlUyp04+NeY1ak4o2g9tUPNBfLzlYwBazmxJ3GAHQmZYI6QUvHEe3ssPzzvtJkJLV+PohP38eziJh6ZbGlkWJShhwQE0qVaGXYmOrd3YjYX07xfwpy471ZDFEHmTQ20XZPASc4LE5fcvt1FS41AP8xTZlfWdxISLgO973Ofm5HBs13Yad+jkkXUsZUfkg0S2rAzAnnXalzt+7UqH69aapIXsSP79d5vlfuhjNuY7mXrS2S5aEhRiuQCcehbGluPG+hU5OqEvRz6wbS6QkZ1nV4S+GNiWoxP62hehseUsRejVYy6LUE5eDtvOaes9L8e+bLe8zMsjJz9xQ5Mtm1165oLPdgLQe1iMS/WLi1P79pKZnkaDto47/tpCRf72QXoPa8G0Z1YhhPmvYm5OtkOhOMv27s1JXgDg7McfU2306CLLtazS0vy8//V2b1R0jbHJltOh/EViMTbZQkC2H79Ev6nrbTbVsEopVrzU1fFnT+0E5wrYX7mwO6ZHbzc0uLn91NFH7jZnwwp0IaNFXp4kL1eblvn6QvWhbZsJDAqiXkxrj7SnhMgH0U/PqjVoy9nD29g8/1duvPdBh+pXH/cuZ958i4vffGtViADiBscRM0P7yxszI8ZzYnRqB0zvortWznwPaFO3gudCzX7cGNLOWV7r/CL0fNutZp9Z8YzpeOPDG+2Wlzk5ZB7QbL4ab7AtstbY9udRAGpH+bb9sJSSw1s3UrtZDCHhnnHINWRqJoSoKIRYLoQ4kP+zyH95IcS3QohzQggXXcb9l9uHawJxOUmLY7Pux1kO163wwAOm4yMP2HbSnNJjiul4yREPpWau2VoTnYLe6WPLaZ9VE9xr//Aqc1sFRWhsstsidDzlOKsTtZBbA6MHEhFs/8uW0MI8lQqq4JqQbJyvWXvfOtT97XBvcvbwQS6dPkWTG12b8haFUSOiMcAKKeWE/Jz3Y4BXiyj3PfA54G7qIr+jQRttfUgEmIf4OVlZphxS9qg35weOPTyQjLg48rKyrNqz3FL7FtPxK2te4cYaN1IhzEN/kd/MT3Q5oS5k6KZJqz7QPnpePgCli9h9yUqDyW20dSdbuDkNu0ZuXi59fzOP1sbcMMZGaY2sY+a4UU137nDpuacPmk0hwss4b3tUnMT/s4rAoCCadOjssTaNEqK7ga75xzOAVRQhRFLKNUKIyOLqlK9Rt3klju9JIjCkOblZe1g18yt6PuGYs2VE27am430tWxGdYH2DUz9Fu+XHW9g5aCcBwoOD5TH58aPTL8JHVmx+/q+x8+0Omg8Nutgv5wStZ5nXPBydqh66TTNXKNO7NwGhrkUp/PX/tEXxfi857s9mBHl5uexbv4b6bWLdyuxaEKN2zapdy0+W/9NtQ4SSknJazx3PaAvKQRGaS8HO5c5NnaJ2m79IF2fantrpv3StZrbyTgaJiIr5C9r5nyJiYdukYkPL+h4WoWtiDLDlEcdiQSW+8ILpuPYk15IMpCSZc8bXbOzb60PHdm4n7fIlojt39Wi7XhsRCSH+QnOYLch/vPE8KeV0YDpAbGysD8YOdB4hBBHlQkjXzTounDhG5Tr1HKsfFETFIUO4+P33nB0/ngoPPYiwMbXbPHAz7X/QtmNbzmzJrkG7vOvrdMd/fSZDiF6EFtyzgNBA+yOb7NOnubJEs7RusND1eE+z/rMBgB6Do+2UNJ4dyxYRUa48DWM72C/sBF4bEUkpe0opWxTx+QM4K4SoAZD/85zt1q5fBo3vBEBwKc1J9YfXnbPWrTbGPONNaNnKZtmwoDALD/2WM1v6TW4td9CL0Jw+cyxiN1lDSsnBbtpINbxVK0IbuRYX/PJZs19gVEfv51lzh+RzZzi8fQste9zmVlbXojBqajYfuGaYMRj4w6B++DyBgQFUqFGKwBBtDSUnK5NsJ9MYR+0129fs79jJZtnyYeVZ/aA5SUvLmS25kuVgCFY/I0/mWYjQd7d9R0wVxwwJE6KbmY4jf5znch9+eFuzPL/1Cd/eKQNtaUAIQcueNgLPuYhRQjQB6CWEOAD0yj9HCFFTCGGKeyGEmAtsAJoKIRKFEEMN6a3BDHjrBgACQ7WF1AWTnNv+FgEB1JulbTzmXrpkd72oYlhFNg80WwZ3mtuJv4795dQzfZ0Dlw7QaqZ5hLj0vqXEVncs5c+Zd981HTfdttXlPhzbk2Q6bhzrXsxnb5N1NZ24FUtpFHsjZSpV9nj7hgiRlDJJStlDStk4/+fF/OunpJR9dOUGSClrSCmDpZS1pZTfGNFfoxFC0Kp7HYLCuwJadg9HIjHqiWjfnrJ9tL9kZ8ePJyPBdgCzsKAwiwXsF1a9YDF68GdeWPkC986/13S+eeBmapauaaOGmeRFi7g0Zy4Adb/7loAI1wz6pJQszHfn6P+6Z9wkvMnO5UvISEul/d3eCbihfM38hM79GyNEAAFBWkjOpV9OdrqNWhMnmo6P3NOPnAv240wX3MKOmRFD0tUkK6V9m+zcbGJmxPDXcfPoLm5wHGFBYQ7Vz9i7l1MvaT5nlZ58glIdHcgia4XVc7UwL+FlgqlS1/UMqcVBdlYmWxb+Rr2WbajRyP0gaEWhhMiPePCN9gSX1uIJ7Vn1l9OjIsDCnuhA55vJTU2zUVojbnAc4zuPN513/akrrWd6xseouBi5YiRtZ5ttq+5vcr9TLi2Zh49w5F5tNBDepg1VX7IfDcEaqZcy2LNGczR+ZJzrYlZc7P57GenJl+nQz/FUSs6ihMiPqFy7DAEBQaZR0Z9fTLJTo2j0YrQ/NpbcFPtB7+9seCe7Bu0ynefKXGJmxPDpNvvxm41k+bHlxMyIYU2iOdXQ1ke28nZHx91AMg8d4nAfbcVAhIYSOXeOnRq2mfGa5ovW+YHGhIT5trtnVsZVNv72E7WimnskAJo1lBD5GU9N7WYaFe1d8zc5Wa5FP7QQoxs6kH36tN06QgjiBsdZBI//Ou5rYmbEMGXHFBs1i58/j/xJzIwYXlxlNnd4ts2zxA2Os5uhVc/VnTs53FeL/yeCg4ly0YXjGv/8dMB03KqHe5kvioMtC34j7fIlbhk4xKs2ZUqI/AwhBPe/2oHAEG27d8boZ11uSy9GB7t1J229Y17j3et2J25wHH0bmH2ypu2cRsyMGNrMamOo7dHAxQOJmRHD6DXmqANlQ8oSNziOYS2HOdXWpZ9/5uiDDwEQULYsUXG77NSwzYXEK+z8W8sp98RE17J7FCdply+xZcGvNO7QiZpNvGtsqYTID6neoBw1o7X1istnTnLpjP3RjDX0YnT88aGcfsvxKcuEmycQNziOW+uZY1Xn5OXQcmZLYmbE8N6/77ncL2f4bvd3xMyIIWZGDLvOm8WieqnqxA2OY92AdTZqF03iqOc486YWnLTUzTfTdJP9UCC2yM3J48f3NJOIXkObERrh2xEYAdb9NJvcnGxuHmA/FpO7iJJoORsbGyu3bHEub7w/Mvmxz8lO11wMXvpxoVtt7WsXS16aeeHalpOsNZYfW24xFSrIxK4T6VWvl0v907Pj3A4eXfKo1fsvx77sUCCzopBSWhgrVn35JSo94WaSAWDKCC0bb83G5en3Uls7pY3n1P545r45mnZ976broCcdrieE2CqldMwgS19PCZH/kpOdy6ePaFEBqzW6i0fed27qUZDzU6dyYfJnpvNGK/+2yCLrDMOWDWPD6Q0Ol+9Wp5spamRWbhZn0s6w5MgSMnIdsyKvEl6Fvx74y62oARkJCRy5x5zlJPLnnwmPcX+BdunXuzm4RfNievqLbj6dqwy0eNSzX3uejLRUHvtkqlPBz5QQ6bhehAjgyM4j/DpeWyfqNfwzWnZ3Lwd59unTJh8qgPDYdkTOnu1Wm5vPbObxpU562jvApK6T6FGvh0faOtT3DrIOHTKdR+2OQwS5v6O1a+UJ1v6oLVA/8d9bCA337V0ygM0LfmXN7G+566XXaXyDbZeggrgqRL7/r6KwSf1W9anRpCOn929g+ZfPU7nODLfiHQfXqEF0QjzxUdri5NUtW4mPiqbBgvmENnYhZhDQvnr7QjY7lzMuM2HzBIucYda4ofoNvHnjmw45ozpL+pYtHHvEPM0r168fNT8Yb6OG4yQmXDSJ0IC3OviFCF08lcj6n36gQbsbaNS++Gyc1IiohPDJg9oWc3DEbQyaMJTy1dyPJZy6bh0nhlquj0TF7UK4mefcF8hNTWV/rKVrReN1/1ikZHKHc8dS+PkD7Xew9/AWNGzjfu4vb5Obk8O8t0Zz+ewZBn/8OaUrOv9v4eqISO2alRCGTPwCgOz0pcx+c6VFsC1XKX3TTUQnxBNQtqzpWkJMS+Kjov02PEheRgbxUdEWIlR51LNEJ8R7TIQunkoziVDnBxr7hQgBbPztJ84cOkCvJ0e6JELuoISohFCpVh3a330/AJnJXzDrPxtIueC+GAE03bSRpju2W1xLiG5GfFQ0eU6GJDGK7FOniI+KZl9rcyjWoGrViIrfS5Wnn/bYcy6eTmPuu9pW/w131vcLo0XQdsn+/XUezW7uRpMbPReL2lGUEJUgbnl4iOk4K3URs97YwIXEVI+0HRAWRnRCPI1Wr7K4vq91G+KjoklZtswjz/E05ydPJj4qmoPddYvawcFE7dpJ49WrPLqDdeZwMnPf0USoda+6tO/r3sZBcZGefJkF/51A2SpV6fbYcEP6oNaIShjZGRlMHqyNjIJL3UFgSBP6PN2S+i09G0OmqDWWa9T/43fCmnrHS9sRrqxYQeLIZ4q856ndsIIc3XWBRVM1Y8oOdzUgtk+kx5/hDfJyc/nf+Dc5tS+BAe/9H1UjG7jVntq+13E9CxHA6QP7mPOG5h0eUnYoAYHliO0bSYc73fsls0bis6O4srzovPBl77yTmh+M98qX/xq5qakcGzSIzL1FG2FWeelFKj/puFGes2xZfJSN8w8D0O2RKJp1diy2kS+wds73bPrjF2576nladO3pdntKiHRc70IEsOmPX1g753sAQss/hxCBlK4YyuDxnkuKV5C8tDT23dABcnNtlit3zz1UefEFgqs6t4grpSRz/wHOvvce6Ztt55YPj21HvZkzEQHeXX3430dbOHNYi15wx7OtqNe8eBd53WHP6hX8OfW/tOzRm17Dih5BOosSIh1KiDRmjXmOc0c0I73Q8i+Y1kNGTOlKYKB3v6BSSs5PnEjSV1979Tl6an44gXK6/PPeJCsjh6+eN4cWefS9jpStHF4sz/YEJ/bs4pf336J2dDPufe0djwXD9yshEkJUBH4EIoGjQH8p5aUCZeqgZXitDuQB06WUDgW/UUJkxmRfFBZBYPgI0/WHx3agQvVSxdqXzMOHOfXyaDL27nW7rVJdbqHGu+MIrlb8W+PH9yaxYPJO0/nwz7oQFBxY7P1wlaTEE8x962VKV6jEQ+9+RFgpzyVK9Dch+gi4qEs5XUFK+WqBMjWAGlLKbUKIMsBW4B4ppd3fYiVEZvLycvnvAG2UUKdZS86fNq8DtOlVl073uZYG53pESsnCz3dxPD/ofVTH6vQY3MxOLd8i5fw55o19ldzsbB5+7xPKVfVs0H5/M2i8Gy3VNPk/7ylYQEp5Wkq5Lf/4ChAP1Cq2HpYQAgICeXbGzwCc2LuLyGZbTYn8ti8/zpQRf5N2OdPILvoF546lMPWplSYR6vdSW78TodSLSfw87j9kpadz72vveFyE3MEvUk4LISKBNoDVoDAlMeW0pwgJC2fEl1oKoYR1qzmybS5DPzEH5vp+zDpWzIz3W2tpb5KXJ/nlwy0mS+mIsiGM+LyrW/58RpCekszP771BWvJl7n3tHarVb2h0lyzw2tTMTsrpGVLK8rqyl6SURSb9FkKUBlYD70spf3Xk2WpqVjRXLl5g+lNDAGjRrRe3jXiOnStO8M/P5vClfUe2JDLG83mr/JEdfx1n3S8HTef++m+Teuki/3v/TS6fOc29r79DnWbeSwvlb2tE+4CuUsrT+WtBq6SUhSzghBDBwEJgqZRyYsH71lBCZJ2UC+f5auRjADSM7cA9o98kOyuXma+vJyM121RuwNsdqFijeBezfYXEfZf4479ml5baURW4c1RrAgJ8O45QUSSfO8sv771B2uVL3D36DerFeDf7ir8J0cdAkm6xuqKU8pUCZQTa+tFFKeXzzrSvhMg2V5IuMP3pIQCUr16DoZ9+BWgxla+FM72GEbtrRnFy/yV+n2jpUzfkw5soVS7UoB65R9LJE/zy/ptkZ1zl3jFjvR53GvxPiCoBPwF1gePAA1LKi0KImsDXUso+QojOwFogDm37HuB1KeXiIhvVoYTIPhlpqUx5/CHT+Ytz55uM/w5vP8+SLy3jB907uh01GpYr1j4WF/s3n2H5N5absfe/Gku1+mWt1PB9EuN388cn4wkICOD+/4yjSr3i8XvzKyHyNkqIHCMnO5tPHzGHRn366zmElzF/+Q5tO8ef03db1Gnbux433tUA4YfTFD052bmsnrOPhA1nLK7f90o7qjfwb8Hds3oFy778jHJVq9FvzNtUqF58LidKiHQoIXIcKSVTn3iYjNQrADz83ifUaGy5XHfqwCV++8RyyhIUHECfp1pSp1nFYuuru0gp2b/xDH99b+mTFloqiPtfifVIMDkjkXl5rPvpBzb+9iN1W7TkzhdeJ6y054wVHUEJkQ4lRM6zfPrn7FqhZQS5+eEh3JAf20hPZno2S7/ew4m9Fy2uh5cNoevDTanfqrLPBYbPy81jz9pTrJm3v9C9qI7V6TKgKUEh/mMVbY2rV1JY8vknHNmxlZjut9Jj6NMEetHR2BpKiHQoIXKNfRv+YeGkCQCElirF01/NISCw6C/pqYOX+XP6bq6mFM4026hdVdrcWpcqdcsUuzBJKTkRf5Eti49y+mByofsVqkdw+4iYErUAf/rAPhZMmkD65Ut0HTyMVr1uN+wPghIiHUqIXCflwjm+GmnOuDH448+pXDfSZp1TBy6x9qcDXDhRdBC2sFLBNGxXlYZtqlCzSXmPOdzmZOeSmHCJA1vOsn/jWavl6jSryC0PNvH7qVdBZF4e25bMZ80P31O6YiXufGEM1Ru6luDAUygh0qGEyD2klMwe8zznjmqe+zHdb+XW4aMcqpudlcvetafYseI4qRcdcx0JjQgiomwIEeVCCCsVQkCgIC9XkpebR3ZmLmnJWaScv0puTp79xoBKtUoT2yeSBm2q+KXtjyMknzvL0i8mcWJvHA1jO9D7qReKfT2oKJQQ6VBC5Bni161m8eSPTeePfzrdpR2Y5PPpHNp+nsPbz3P2SIonu0jNxuVp1K4qjWKrEl46xKNt+yJSSvas+ouVM6YjJXQb8iQtuvbymbU5JUQ6lBB5Dn3oWYCmnW6h76jRHvvFl3mSjLRs0lOySE/JIiMtm7xcSUCgIDAogKCQAEqVC6VMpTBCwnw/L5g3uXDiGCu/n87x3Tup3awFvZ96waccV0EJkQVKiDzPtWh+17j75Tdo1P5GA3t0/ZCZnsYXwx4hN1tzwen++Aha9+rj9eiTrqCESIcSIu+Qk5XFs/mNqgAACi9JREFUjNEjuXzmtOnasKnfU6aS/zmC+gMyL4+Fn37E/n//MV17ZMKnPuc5r0cJkQ4lRN7l/PGjzBxtjnFcpnIVBn/8OaERJWdL3Gj0MccBYu+8ly6PPG69go+ghEiHEqLiYefyxfz19VTTeY0mUTzw5vsEh/ink6gvsH3pQv7+dprpvGbTZvR/632PxZT2NkqIdCghKl7+mTeTjb/9ZDqvWLM2D737kYXfmsI2BUdAwaFhPPH5N0SU9S+/NyVEOpQQFT8yL4/lX08hbsVS07WAwCAGfTSZSrXrGtgz3yU7K5Mln3/CgY3rTdfCSpdh8MefF3vueU+hhEiHEiLjkFKy6fef+WfeTIvrNz34KB369fcZexcjOX/sCHPeeJmcLLPBZ7lq1Xn4vU/8bgRUECVEOpQQ+QYHNq1n/ifjLa5FlCvP3S+/Qc0mUQb1yhgy09NYNv1z9m9Ya3G9Va/b6TZkuCEOqt5ACZEOJUS+RXpKMvM/Gc/JhD0W16tGNqTPsy+V2Klb1tV01v04m21L5he6d/foN2kU28GAXnkXJUQ6lBD5Lke2b+H3j8eRVyAtdUh4OD2HPk1U565+PX27eCqRNT98z6Et/xa6F3vnvXR+6FG/2QFzBSVEOpQQ+T5SSvau+dvCWltPtQaNufHeB2nY7gaftCC+RvK5s2xe8Cs7ly0q8n6rW/ty84BB142NlV8JkYMpp8OANUAoEAT8IqV825H2lRD5H2cPH2TFd9M4vT/Bapnozl1p0a0XdZrFGCJOl86cYv+/64j7eynJZ88UWSYoNJSbBwymVa8+JWbdxxn8TYgcSTktgFJSytT8tEL/AM9JKQuPeQughMi/kVJyZPsW/v11HqcP7LNbvlqDxlRv1ISq9epTuW49KtepR0i4c7GH8nJzSb2YxNmjhzhzcD9nDu7j9IH9ZGdm2KwXFBpKbN97aNe3n0+E4TAaV4XIKMm+G+iafzwDWAVYCJHUFPJapK3g/E/Jm0cqCiGEoEHb9jRo29507WrqFfatW8PuVX9x9vABi/JnDx8odM0blKpQkcY3dKT5LT2o1rCxX69l+RpGCZFFymkhRJEpp4UQgcBWoBEwRUppNeW0omQTXroMrW/rS+vb+lpcz83J5vSBfZw7cojzx49y/thRLhw/Qm5OjtPPKF2pMpXr1KNGo6bUaNyU6g0bK+vwYsJrQmQn5bRDSClzgdZCiPLAb0KIFlLK3UWVFUIMA4YB1K1bMreDFYUJDAqmdnQLake3MLorCjfwmhBJKXtauyeEOCuEqKFLOX3OTluXhRCrgN5AkUIkpZwOTAdtjcjljisUimLHqH3R+cDg/OPBwB8FCwghquSPhBBChAM9AetbKgqFwm8xSogmAL2EEAeAXvnnCCFqCiGupZSuAawUQuwCNgPLpZQLDemtQqHwKoYsVkspk4AeRVw/BfTJP94FtCnmrikUCgPwXZNVhUJx3aCESKFQGI4SIoVCYThKiBQKheEoIVIoFIajhEihUBiOEiKFQmE4SogUCoXhKCFSKBSGo4RIoVAYjhIihUJhOEqIFAqF4SghUigUhqOESKFQGI4SIoVCYThKiBQKheEoIVIoFIajhEihUBiOEiKFQmE4hgiREKKiEGK5EOJA/s8KNsoGCiG2CyFU4HyFooRi1IhoDLBCStkYWJF/bo3ngPhi6ZVCoTAEo4TobrSc9+T/vKeoQkKI2kBf4Oti6pdCoTAAQ9IJAdWklKcB8rO9VrVSbhLwClDGXoP6lNNAphCiyIywJYDKwAWjO+FF1Pv5N01dqeQ1IRJC/AVUL+LWfxysfwdwTkq5VQjR1V55fcppIcQWKWWsE931G0ryu4F6P39HCLHFlXpeEyIpZU9r94QQZ4UQNfJHQzWAc0UUuwm4SwjRBwgDygohZkspH/FSlxUKhUEYtUY0Hy3nPfk//yhYQEr5mpSytpQyEngI+FuJkEJRMjFKiCYAvYQQB4Be+ecIIWoKIRZ7oP3pHmjDVynJ7wbq/fwdl95PSCk93RGFQqFwCmVZrVAoDEcJkUKhMBy/F6KS7i7iyPsJIeoIIVYKIeKFEHuEEM8Z0VdnEEL0FkLsE0IcFEIUsqwXGpPz7+8SQrQ1op+u4sD7Dcx/r11CiPVCiFZG9NMV7L2brlx7IUSuEOJ+u41KKf36A3wEjMk/HgN8aKPsi8AcYKHR/fbk+wE1gLb5x2WA/UAzo/tu450CgUNAAyAE2Fmwv0AfYAkggBuBjUb328Pv1wmokH98u7+8nyPvpiv3N7AYuN9eu34/IqLku4vYfT8p5Wkp5bb84ytovnm1iq2HznMDcFBKeVhKmQXMQ3tPPXcDM6XGv0D5fJszf8Du+0kp10spL+Wf/gvULuY+uooj/3cAzwL/o2gbwUKUBCGycBcB7LmL5BVXxzyEo+8HgBAiEmgDbPR6z1ynFnBCd55IYeF0pIyv4mzfh6KN/vwBu+8mhKgF9AOmOdqoUb5mTlHc7iLFjbvvp2unNNpfoeellCme6JuXEEVcK2hH4kgZX8XhvgshuqEJUWev9shzOPJuk4BXpZS5QhRVvDB+IUSyhLuLeOD9EEIEo4nQD1LKX73UVU+RCNTRndcGTrlQxldxqO9CiJZoSwW3SymTiqlv7uLIu8UC8/JFqDLQRwiRI6X83WqrRi9+eWDx7GMsF3M/slO+K/61WG33/dD+Ss0EJhndXwffKQg4DNTHvODZvECZvlguVm8yut8efr+6wEGgk9H99fS7FSj/PQ4sVhv+Yh74h6mEFlztQP7PivnXawKLiyjvb0Jk9/3QhvUS2AXsyP/0Mbrvdt6rD9ru3iHgP/nXRgAj8o8FMCX/fhwQa3SfPfx+XwOXdP9fW4zus6ferUBZh4RIuXgoFArDKQm7ZgqFws9RQqRQKAxHCZFCoTAcJUQKhcJwlBApFArDUUKkKFaEEOtdqBMkhLgghPigwPWjQojKuvOu/hRZQWFGCZGiWJFSdnKh2q3APqC/cNRnQOFXKCFSeIT82DO7hBBhQohS+XGRWhRRLjX/Z1chxCohxC9CiAQhxA82RGYA8ClwHM3KWlHC8AtfM4XvI6XcLISYD7wHhAOzpZT2kly2AZqj+SqtQ/MJ/EdfQAgRDvQAhgPl0URpg2d7rzAaNSJSeJJ30bKyxKIFdLPHJillopQyD83NIbKIMncAK6WU6WhOvf2EEIH594pyC1CuAn6IEiKFJ6kIlEaLEhnmQPlM3XEuRY/QBwA9hRBHga1ovnfd8u8lAfrQuRUp2emcSyxKiBSeZDrwJvAD8KG7jQkhyqI59NaVUkZKLdnmSDRxAlgFPJpfNhB4BFjp7nMVxY8SIoVHEEIMAnKklHPQEma2F0J0d7PZe9Ey/OpHTn+gxZYKBcYBjYQQO4HtaGE1Zrv5TIUBKO97hUJhOGpEpFAoDEcJkUKhMBwlRAqFwnCUECkUCsNRQqRQKAxHCZFCoTAcJUQKhcJw/h+vRwqh4eobMgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(1,1,figsize=(4,4))\n", "\n", "# want to loop over all bodies in this simulation and plot their TOTAL orbits -- all times\n", "# x/y positions only\n", "# r_h[# bodies, x/y/z 3D positions, all time steps]\n", "for i in range(r_h.shape[0]): # looping over all bodies\n", " # now we want to plot (for given body) the x/y positions\n", " # over ALL times\n", " # r_h[ith planet, x position, all times], r_h[ith planet, y position, all times]\n", " ax.plot(r_h[i,0,:], r_h[i,1,:])\n", " \n", "ax.set_xlabel('x in AU')\n", "ax.set_ylabel('y in AU')\n", "ax.set_xlim(-0.4, 0.4)\n", "ax.set_ylim(-0.4, 0.4)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAScAAAEICAYAAAAdoDKiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3hUVfrA8e+bkNA7oRMSaoAgIFGkSYLEtWPvgl0U9+euDSy7uqurYlnbIuiia9e1odhWQRIUAQWUEiD0UEOvoaSe3x8z3Jkkk8kkmZk7M3k/z8OTc++cO+cdiS937j33vGKMQSmlQk2U3QEopZQnmpyUUiFJk5NSKiRpclJKhSRNTkqpkKTJSSkVkmxNTiJyloisFpF1IjLRw+siIi85X18mIidX4dh7RcSISKtAfw6llP/VsWtgEYkGJgPpwFZgoYjMMMasdOt2NtDd+WcQMAUYVNmxItLJ+dpmX2Jp1aqVSUhI8MvnUkpVzeLFi/cYY+LK7rctOQGnAuuMMRsARORDYDTgnpxGA28bx0zRBSLSTETaAQmVHPs8cD/whS+BJCQksGjRopp/IqVUlYnIJk/77fxa1wHY4ra91bnPlz4VHisiFwDbjDFL/R2wUip47DxzEg/7yj5LU1Efj/tFpAHwEHBmpYOL3ArcChAfH19Zd6VUkNl55rQV6OS23RHY7mOfivZ3BRKBpSKS49z/m4i0LTu4MeY1Y0yKMSYlLq7c112llM3sTE4Lge4ikigiscCVwIwyfWYAY5x37U4DDhpjcis61hiz3BjT2hiTYIxJwJHETjbG7Ajap1JK+YVtX+uMMUUicifwHRANvGGMWSEi45yvTwW+Ac4B1gFHgRu8HWvDx1BKBYjokimQkpJi9G6dUvYQkcXGmJSy+3WGuFIqJGlyUsqP1p91Nnv/86bdYUQETU5K+cmmMWMpyMlh16RJdocSETQ5KeUHpqSEo7/+CkDr+++3OZrIoMlJKT/I7t3Hare88QYbI4kcmpyUqqG8OXOsdg+96+s3mpyUqqEtt42z2tGNGtoYSWTR5KRUDaxy+zrXK3uVjZFEHk1OSlVTyZEjUFICQMcpr9gcTeTR5KRUNa0e6JrU3DgtzcZIIpMmJ6WqwX2iZdKKLPsCiWCanJSqhhMTLeunDESio22OJjJpclKqilYl9bLaCe++a2MkkU2Tk1JVULhtm9VOnOHTEvWqmjQ5KVUF684YZbXr9ehhYySRT5OTUj7a5vbMXNKqlV56Kn/Q5KSUD4wxHJrxJQAtrr8eEU81NpQ/aXJSygfZvXpb7TYTJ9gYSe0RceXIReQxZ98lIvK9iLQP1udRkenob79b7e4/z7UxktrFtuTkVlL8bKA3cJWI9C7Tzb0c+a04ypFXduwzxpiTjDH9ga+Avwb6s6jItunqq612nZYtbYykdrHzzMkqR26MKQBOlBR3Z5UjN8YsAE6UI6/wWGPMIbfjG1K+UKdSPluX7qrPqg/2BpedFX89lRQf5EOfisqRW8eKyD+AMcBBQB96UtViCgoo3OL4NWv3xBM2R1P72Hnm5Pdy5FbDmIeMMZ2A94A7PQ4ucquILBKRRbt37/YxZFWbZJ/Uz2o3u/giGyOpnSKtHHlZ7wOXeBpcy5Erbw5+4Zr93XPZUhsjqb0iqhw5gIh0dzv+AiA70B9ERZ7tExw3gOu0a0dUbKzN0dROkViO/CkR6QmUAJuAcShVBe4P9nbPmG1jJLWbnRfEMcZ8gyMBue+b6tY2wHhfj3Xu9/g1TilfFO3bZ7U7v/uOjZEonSGulJu1Q4Za7QYpKV56qkDT5KSU085nnrHaSStXeOmpgkGTk1I4Huzd9/obADQ552wkSv/XsJv+DShF6Qd7O/zznzZGok7Q5KRqvfy1a61215nf2xiJcqfJSdV6G86/wGrHdurkpacKJk1OqlbbdN0Yq62rW4YWTU6q1jLFxRxduBCA1vfeo6tbhhhNTqrWyu6TbLVb3nyzjZEoTzQ5qVrpcEaG1e6xaKGNkaiKaHJStdLW2++w2tGNGtkYiaqIJidV67g/2KurW4YuTU6qVinOy7PaHV95xcZIVGU0OalaZU3KKVa78cjQXMG5uLiEnOV77A7DdpqcVK2xd9o0q52UtdzGSCq2fe0Bpo7P5OvJyziw86jd4djK1vWclAqmXc8+B0D9lIFIndD71Z/1n5Ws/mWHtd2sTQMbo7Ff6P0NKRUA7hfBE95918ZIyjPG8MrtrqkNfU7vQOrVPW2MKDRoclIRr2DrVqudOOMLLz2DL2//cd56YJ61femEFNokNrExotARieXInxGRbGf/6SLSLFifR4Wm9aPSrXa9Hj1sjKS0rB+3lUpM4/6VqonJTSSWI58JJBtjTgLWAA8E+KOoELbt7rutdig92PvWAz8z5/3VADSNq8/4qSOJrqP3p9zZ+bXOKikOICInSoq7/wZZ5ciBBSJyohx5QkXHGmPcF+RZAFwa8E+iQpIpKeHQN98C0OKmG0Piwd6iwmJe/eMca3vkmCR6DWlvY0ShKyLLkbu5EfhvjSNVYSm7dx+r3ea++2yMxCF3/UE+e2axtX39pKE0bFrXxohCm53JKWDlyAFE5CGgCEdJ8vKDi9yK46si8fHxlcWqwsyJpVAAus/72cZIHGa/vYpV83Kt7TumpIXEmVwoszM51aQceay3Y0VkLHAecIbzK2E5xpjXgNcAUlJSPPZR4ct9Ebk6LVrYFkfZaQK9h7Un7dok2+IJJ3YmJ6ukOLANR0nxq8v0mQHc6bymNAhnOXIR2V3RsSJyFjABGGGMqd1TbGuptamux1LsfLD3yIF83pzoOmu75P6BtO3S1LZ4wk0kliP/F1AXmOk8bV5gjNGS5LVEyfHjFO1wzLJuP+kp2+JYOXc7Ge9mW9vjXk4lOkbvxlWFVPCtp1ZJSUkxixYtsjsM5QehsBzKOw/P49Ce4wA0alGXsU8MreSI2k1EFhtjypVX1hniKmIc+OQTq91z2dKgj19cWMLUP2Za26nX9KTP8A5BjyNSaHJSESP34b8AENM5nqjY2KCOvWPDQT592jVNYOyTQ2nUXKcJ1IQmJxUR3L/Odfvuu6COnfFeNit/ct1o1mkC/qHJSYW9ot27rXbnD94P2rhlpwkkDW7LGWPLPoGlqkuTkwp7a4efbrUbDBgQlDGPHMznzQmuaQIX3Xsy7bvpM+b+pMlJhbWca6+12kkrV3jp6T+r5uUy+23XncDbXh5BnZjooIxdm2hyUmHLFBdzbJHjInSjkSORqMDPI3rvkQXW8rkNmsZyw6RhAR+zttLkpMKWe8XeTq9MDuhYxUUlTL0z09oecVUPkkd0DOiYtZ0mJxWW8ua6rvd0+2FWQMfatekQHz/pmqQ75okhNG5RL6BjKk1OKkxtuflmqx3TIXATHX/8YDXL52yztnWaQPBoclJhx31OUyBXt5w8brbV7nFqG9Jv7OOlt/I3TU4qrJQcdS00EffnPwfkLObooQL+c/9ca/uiewbQvntzv4+jvNPkpMLK6pMHWu1Wt93q//dfkMusN92mCbw0gjqxOk3ADpqcVNjY+583rXaPAKwi8cHff2Hf9iMA1GsYw03PDff7GMp3mpxU2Ng1aZLVjm7U0G/vW1xcwtTxmdb28Cu6c1Jap4oPUEGhyUmFhUCt07R782E+esK13vh1/xhMk5b1/fb+qvo0OamQV7hzp9XuOOUVv73vTx+tYdlsVzVgnSYQWjQ5qZC3bkSq1W6cllZxxypwnybQLaU1f7g52UtvZYdILEd+mYisEJESESm39KcKL6Uq9mYtr/H7HTtcUCoxjf5Tf01MIcq2Mye3kuLpOEpALRSRGcYY91l17uXIB+EoRz6okmOzgIuBV4P2YVRAGGOsir31evdG6tTs13XNwh3MfN3163XrSyOI0WkCISsSy5Gvcu4L2gdRgZHdy7VwW+Jnn9bovf77j1/ZsyUPgNh60dzywogavZ8KvEgvR67C1LElS6x2ly9nVPt9SopLmOI2TWDopd3oP0orPIeDiC1HXungWo48pOVceZXVrtu9e7XeY/eWw3z0D7dpAo8PpkkrnSYQLiKyHLkvtBx56Foz1LWAW3Uf7P3503UsmbnZ2tZpAuHHzrt1VjlyEYnFUVK87Pn7DGCM867daTjLkft4rApDJQUFFO/dC0DzMddVK6FMHjfbSkxd+scxfurIar3P1KVT6ftWXxbvXFx551rOlJSwf8d2Vs+fy08fvMVnTz5CwfFjNXrPiCtHLiIXAS8DccDXIrLEGPOH4H46VV2rT+pntds++GCVjj2eV8jr9/5kbV9wV3869WpR5RjunXMv3+W4yktd/7/rWT625tMYIkVxUSF7t25hV84GduWsZ9fGDezetIGCY45kFBUdTcsOnTh64ACxbav/NVrLkaPlyEPFgc8/J3fiAwB0nz+POs19X6Zk7aKdfD/NVeDglhdOJ7Ze1f7tPW/6eWw6tKnc/tqcmAqOH2N3zkZHEsrZwK6NG9i7dRPFRUUAxNStR1znRFondiGucxfaJHalZcd46lShqKmWI1ch70RiAqqUmD5+ahG7cg4BEB0TxbiXU6s0bt+3+nrcX9uS0tGDB9i1cT07czawO2cDu3I2sH/HdnCewNRv3ITWiV05+ZzRtE7oQuvErjRr246oqMDMFdPkpEJCdR7sLTtNYPDFXTn5zM6uDo82dfx8aCfElF/zu7YmJWMMh3bvZNfGDW5nROvJ27/P6tMkrg2tE7rQa3gqrRO60jqhC41atAzqTQVNTsp2Rfv3W+12Tz7p0zF7t+Xx4WO/WtvXPnYaTeMaODZ+extm/NHVeekHkHKDtVmbklJJcTF7t21h18b11jWi3TkbyT/qWLdKoqJo0b4jnZL7Oc6GnImoXqNGNkeuyUmFgLWDh1jtZhddWGn/+dPX89t3rmtD1jSB4iJ4rGXpzrGNrcQU6UmpMP84uzfllLpQvWdLDsWFhQDUia1LXHwCSUNPt5JQy/jOxMTWtTlyzzQ5KVvtdDtT6rlsaaX93R/aTejbknPHO+/unfgK5+7RgxhjOCkCk9Kxw4fKfS3bn7sdY0oAqNewEa0Tu9D/D+fRxnl9qHm7DkRFh8+zhBUmJxE5TOlZ1wbYA2QAE4wxewMcm4pwxhj2vfU2ANHNmxPl5Q7P8SOFvH6Pa5rAeX/sR+c+LWHZx/DZzaU73zoH066fx6TUtWlXPr/wc/98gCAwxrBv21b2524rdcfs8N7dVp/GLeNondiFnkOGO86IErvQuGVc2E86rTA5GWMal90nIs2B64GpwGWBC0vVBu4P9vaYP6/Cfut/38X/Xs2ytm954XRiY/B4tlT0170MeGdAuf1nJ5zN0yOerlnAAVZcVMTG3xeRlTmL9YsWlHtdJIrm7TvQIam3dX0oLiGRBk08nDVGgCp9rTPG7AeeF5HrAhSPqiXy162z2gkfflBhv8+eXUzuuoMAiMAdU0Z6TErHHsrl1PcHQZnEdHu/27mj/x1+itp/8o8eYfW8n8jKmEnuutVe+7bt2p20628jLj6BmHq1p9Jwla85iUhMdY5Tyt2G88632vX79y/3ekmJYcodGdb2aRd2YWD80nKJaf81H3P6vHvg/dKLUjw+9HFGdxvt56ir59CeXayY8wNZGbM4tHun176develT+ooegwaWqsSkSferjld7GF3c+AK4JOARaQi3sZLLrXaSStXlHt97/Y8Pvy7a5rANY+eQrOpHcHtm86WOnU4p1N7mHdPqWOnnTmNQe3sWT3HGMOujetZnjGTFRkzKSos8Nq/5+DhJKelE9+3X8AmMoYzb2dA55fZNsBe4EVjzNeBC0lFMlNUxPEVjoTU+KyzkKjSz57/MmMDi77JsbbvaHMxMtV1X2Zp3Viubd+23PtOv2A63Zp3C0zQHpQUF7Np2e9kZcxkzS8/e+0bW78+fVJHkZyaTlznxLC/UB0s3i6I31DRayJyijFmYUWvK1WR7GTXHbSOLzxf6jX3aQLxnYs5P991hjWzQX3ubhNX7v1+uOwHWjdoHYBIXQqOH2PN/LlkZc5kW7b3JVyatW1Hcmo6vU8fSeOWrQIaV6Tz+dqRiPTGsTTJVcBBQIsHqCo5nOG6htQtw5WI8o8WMu1u1zSBc5s9TkK+Y5mSt5o05tmW5Z+zm3/VfBrF+n8W8+F9e1g5ZzYr5sxif673JcI6JPWmT+ooep42jNj6DfweS23nNTmJSGccyegqoAjoDKQYY3ICH5qKNFtvd901i2nXDoANS3bz7VTXZMhbWl9NbNQxnmjRnA+alpvNwuJrFxMb7fsT797s3rSRrIyZZGXOouDYUa99u586hOS0dBL6nRxWExnDmbcL4vOApsCHwKXGmLUislETk6oO9wd7T6xu+fnzv7Ft9QFr//i2F/FQqxbMaFz+69vSMUuJkuqtjVhSUszmrGVkZcxk9bwfvfatExNLn7R0klNH0aZLN70+ZCNvZ067cSx/2wbHwm1rqeI63UoBFOcdsdqt77sPDEy+3fW17tRGH/BcjzlMrVt+LfdlY5ZVKUEU5h9n7a/zycqYyZYVy7z2bRLXmuTUdPqMOIMmcYG9bqWqztsF8dEi0hS4BPibiHQDmonIqcaYXys6Tqmy1qS4Lk9GnXsFr7jNX/q43+NMbbAbKP3wqS/PvR05sJ+VP2WQlTGTfdu2eO3btlsPklPT6TlkOPUa2v/Evaqc12tOxpiDwBvAGyLSGsccpxdEpJMxppO3Y5UC2PPav6324Sc/Z/ajv1jbU0/7E0jpk/GKktLerVvIynRcHzp++JDXMbsMPJXktHS6DEghuk5MDaJXdqrWMr0i0tkYU3490zCly/QGzolrTbNTJ1v7tjZZw1d9JpfqdyIpmZIStq7KIitjJit/ysAbiYoiOS2d5NR02nXvqdeHwpRfl+n1V2ISkbOAF3EUKZhmjHmqzOvifP0cHAUOrjfG/ObtWBFpAfwXR1XgHOBy5zOBKshWJfWiKLoePw5/ztr3bc9/s6mF4yHe6GL4sOe/yMqcxXNXnOf1vRq1aElyWjp9Tj+DZm3bBTRuFRpse0ZORKKByUA6jvp0C0VkhjHGfZbb2UB3559BwBRgUCXHTgR+MMY8JSITndsTgvW5lEPhtm3sadGHZSc5pg+YkmPMa/s0/VbUJe2wayndr797ptyxrRO6kpw2iqShI6jfuEnQYlahxc4HeE8F1hljNgCIyIfAaMA9OY0G3jaO754LRKSZiLTDcVZU0bGjgVTn8W8BmWhyossDX1MSpHutzQv38sje7WxJaEDB/n9a+wceLD9vKad+PCsbJ7GxQQIl4jZ/KBPI/Klc/1LHPnWunyJWoajS5CQiccAtOBKC1d8Yc2MNx+4AuN9i2Yrj7KiyPh0qObaNs/Amxphc54X8cmpbOfKAJCZjaJe/g16Hs+mTl23tFqJISbyP+vuFrAOuBLOyUU9WNU5ie912jvVPaij34DHaNdXy4pHKlzOnL4CfgFlAsR/H9vTbWfZ/oYr6+HKsV7WtHPnfR/fhr1+UXwHAV1GmmMSjOfQ+nE3Csc1e+xqpx6L9q2DTXF7uf3u1x/QmNjpKE1OE8yU5NTDGBOJr0VbAfTpCR6Dsw0wV9Yn1cuxOEWnnPGtqB+zya9RhaszgBMYMTvCp7/G8PLJ/nkNW5ix2bljrte++xgWs65hHQ0lhYK6rOEGPzDvpnb2SyssVKOWZL8npKxE5xxjzjZ/HXgh0F5FEYBuOh4qvLtNnBnCn85rSIOCgM+ns9nLsDGAs8JTz5xd+jjuiHNy1g6zMH8jKnEne3j1e+8b37U9yWjoNesVz3pcXWPtvm/88guPRkqYH1jFwyfN0+fqrgMatIp8vyeku4EERyQcKcXylMsaYGt1GMcYUicidwHc4pgO8YYxZISLjnK9PBb7BMY1gHY6pBDd4O9b51k8BH4nITcBmdK1zwLEQ2o51axwTGTNmUVJc5LV/r+FpJKeOolPvvtaaS9n7srnsy8scfxtATHFdbvrVtS53ctZrtN7jqKBSt2vXwHwQVWtUaxJmpIm0SZjFRUVsXLKYrIyZHhfKd1evYSPnQmijaBWf4LHPgtwF3PL9LaX2ddqfxLnZrutJF8T/Tt7b0wDHg706IVL5qsqTMEUkyRiTLSIne3r9xGRIZa/8o0dZPe9HsjJnkrvW+0L5Ldp3JDktnV7D02jUvEWl7/31hq+Z+NPEcvsf/+06tua7fpfumDyC7D7jAWjz0EOamJRfePtadzeOW+3PeXjNACMDEpGq0KE9u1k5x3F96OAu7wvld+ydTHJqerUWyn8j6w2eX/x8uf2LN2zh3zs/Y6tz++SzOjP4wq6llkNpcd21VRpLqYp4W5XgVufPtOCFo8C1UP6J60NFBfle+/ccPJzk1FHEn9S/RgvlP77gcf67+r/l9i/ZuJm8orb8e89n1r4r/3IqLTs04thy14O6PRbMr/bYSpWlJZ5sVlJSzKalvi2UH1OvPsmpo+iTOoo2if674Hzr97cyP7d8Ylm2cTMC/JZ3IfPzxlr7b5+cSlS04yJ5zmWXA1CnbVuimzXzW0xKaXIKooLjx1iz4GeyMmayLdv7hMhmbdqRnBbYhfLP/ORMco/kltt/IikBTNn1KSUljkTUJrEJl05wXWvafKPrIYHumd5XEFCqqjQ5BUje/n2s/HE2WRkz2Z+7zWvf9j17k5w6ip6Dg7NQft+3+nrcv3yja+Z3QUk9/r3LVYn3D7ck022g60mgkmPHODLPcbbVceqUAEWqajNfnq37AXjOfRKmiLx24pqUgt2bcxwL5WfM9Gmh/D6po0jsPzDoC+V7SkpRCEs3ll4BZ8voZcx4db21fdOzw6nXqPSibasHuG7iNk5N9W+gSuHbmVMiMMFZq+5vzn21vizUvI/fZ/4n71f4enRMDMmpzoXyu3a39fa6p6TUrWkXpi/JLL1z0O18k3M1G90S0/ip5W/KHvrue6vdc9lSv8WplDtfktMB4AzgJRH5EtB7xcC21a6VXRq3irMWym/auo2NUbkYYzjp7ZPK7U/vnM4/M1/HMXnerf8jB3jl9gzA8QjLgPR4hlxSvoKuMYZtd90FQKMzziAq1j9lmpQqy5fkJMaYIuAOEbkemAuUr3JYy1z28ON2h+BRRUnphuQbuLuoEXx7X+kX7lvPwaMNefd21wXtKx4+hVYdy6+9BLDmlFOtdqfJ//JP0Ep54EtymnqiYYx5U0SWA+MDF5KqjuKSYvq/07/c/ocHPcwV3S+Bx1qWfmHAtTB6MktmbebnT1wllNynCZRVtG8fJXl5ACR+Pt1/wSvlQaXJyRjzapntxUBNF5pTfpJfnE/Ku+UvAb6Y9iIj40fCo01xrqnn8uhBAF67aw6F+Y4luuLiG3P5g6d4HWvtkKFWu15SUs0CV6oSOpUgTB3MP8iwD4eV2//uOe/SL66fY2PPutIv3rMaGrelML+Y1+6aY+1Ov6k3PU5p63W8vdOmWe2kldVftE4pX2lyCjM7juwg/ZP0cvu/vPBLEpomlN7ZynlBO+k8uPI9ALZm7+OLF5ZYXW58Zhj1G3u/qG1KStj1rOMRy1Z33G4toaJUIGlyChOr963m0i8vLbc/4/IMWtX3MoPc+RUO4H+vZbH+N9fCoHdMSfNpikN27z5WO+7//s/HiJWqGU1OIc7TWkoAv1z9Cw1ifJtNboxxThNw6DeyE8Mu7+7TsfnrXXOeus2Z46WnUv6lySlEzVg/g4fmPlRu/2/X/UZMlO8ltg/tPcY7D7ke6r38wVOIi/c8TcCTDee6il3GtPFYyEapgNDkFII8zeheNmZZlWeZL529hbkfuQoUjJucSnQF0wQ82T7BVdciadVKLz2V8j9brmyKSAsRmSkia50/PU7qFJGzRGS1iKxzVu/1eryItBSRDBHJE5GImCG4fOxylo9dXuXENO2eH63E1LJDQ8ZPHVmlxGQKCzn4xQwA2j89SVe3VEFnyxriIvI0sM+tZHjzsuWnnCXH1+BWchy4yhizsqLjRaQhMABIBpKNMXf6Ek8krSFeWFDMa//nujY06obe9BzkfZqAJ+6rW/bKXuWX2JTypKI1xO26JzwaR6lwnD89lTezypUbYwqAEyXHKzzeGHPEGDMXOB6owEPZttX7SyWmG54eVq3EdGSBqyhCz8WRkbRV+LHrmpMvJcNrXHK8Nvn+9RWsXehaV9zXaQKebL7+BgDq9e1LVMOGfolPqaoKWHISkVmAp3+2y9+CquAtPOzz23dQEbkV53Md8fHx/nrboCs7TaBvakdOv7JHtd9vw/muYpmJH39Uo9iUqomAJSdjzKiKXhMRX0qGeytXXuOS48aY14DXwHHNqarHh4LD+47z9oPzrO3LHkihdefq1zotzssjf63jInrnd9+pcXxK1YRd15xOlAyHikuGW+XKRSQWR8nxGVU4PqItz9xaKjGN+1dqjRITwJoU14O/DVJq/XqCymZ2XXPyWDJcRNoD04wx51S35LiI5ABNgFgRuRA40xgTUZN03rh/LscOFQDQrE0DrvnbaTV+zwOffmq1k7KWe+mpVHDYkpyMMXtxrK5Zdv924By37W+Abzz083i887UEvwUaYooKinnV7W7cyDG96DWkXY3f1xhD7kMPA9Ds8suROjo3V9lPfwvDxPZ1B5j+rKsC/PWThtKwaV2/vHd2r95Wu93f/+alp1LBo8kpDPzw5kqyF+ywtmsyTaCswlxX3bqu3/3PL++plD9ocgphZacJ9BnentRr/LsC5bo0V3WV2M6d/freStWEJqcQdexwAW/cN9favmTCQNomNvXrGDufecZq64O9KtRocgpBm7L28tW/XPXgxr2cSnSMf2d9mOJi9r3+BgBtHpioD/aqkKPJKcT88PYqsuc5rgP1H9WJoZf6tihcVWX3SbbaLcaO9dJTKXtocgoRxUUlTL0z09q+6J4BtO8emPKAx7JcBQq6z5/npadS9tHkFAL2bT/CB3//xdq++fnTqVs/cH81OZc61iKvExdHnea1vj6qClGanGzmvlplu25NufjegQEdb/PNrvXIu//0Y0DHUqomNDnZxBjDh4/9yr7tRwBIuy6J3kPbB3TMktSLFOwAAAwDSURBVPx8jsx13AHsOOWVgI6lVE1pcrLBsbwC3rjXNU3g2sdOo2mcb5VUamJ1P1e58sZpaQEfT6ma0OQUZJtX7OXLlx3TBCRKGPfyCKKqsLZ3dR2eNctq91y21EtPpUKDJqcgynhnFSt/dkwTqErtuJoyxrD1zj8C0Cg1lahY7xV+lQoFmpyCoOw0gQv/PIAOPYN3l2ztaYOtdqepU4I2rlI1ockpwPbvOML7j7pNE/jncOo28L0oZk0V7d9P8UFHSfLE6Z8FbVylakqTUwAtz9zKjx+uAaBtl6Zccn9gpwl4snbwEKtdr1cvLz2VCi2anALAGMNHTyxkz5Y8AFKv6Umf4R2CHsde57NzAEkrV3jpqVTo0eTkZ8ePFPL6PT9Z29f87TSatQn8NIGyjDHscq460HLcbUiUXcvFK1U9kVaOPF1EFovIcufPkZ7eN1C2rNpXKjGNm5xqS2KC0qtbtv7Tn2yJQamasOuf04nAD8aY7sAPzu1SnOXIJwNnA72Bq0SkdyXH7wHON8b0xVGVJWj1jea8v5oZLy4BoO+IDoyfOpLoIMxf8iR/w0ar3W1Opi0xKFVTdn2tGw2kOttvAZnAhDJ9rHLkACJyohz5yoqON8b87nb8CqCeiNQ1xuT7/RM4FReXMHV8prU9+k/96ZjUIlDD+WTDOVaNCGLatLExEqWqL5LLkV8C/B7IxHRg51Hee2SBtX3Tc8Op1zB40wQ82f7Ag1ZbV7dU4Swiy5GLSB9gEnCmlz41Kkee9eM25ry/GoDWnRtz6cQU21eTNIWFHJw+HYD2k56yPR6laiLiypGLSEdgOjDGGLPeS3zVLkf+yaRF7Nx4CIARV/UgeUTHqhweMNl9T7LaTUePtjESpWouosqRi0gz4GvgAWPMz/4O+viRQiaPm20lpqsfHRQyienIL79a7Z6LF9kYiVL+YVdyegpIF5G1QLpzGxFpLyLfABhjioAT5chXAR+VKUde7nhn/27AX0RkifOPp+tRVVZYUFxumkDztg398dZ+sdm5Dni9vn2Jahg6cSlVXRFVjtwY8zjwuF+Ddco/UgQEpnZcTW24wPUVLvHjj2yMRCn/0RniPmrUvC7jpwZ1TqdPivPyyF/jeH6v83vv2hyNUv6jzzSEuTUpp1jtBgOD/2CxUoGiySmMHfjkE6udtCLLxkiU8j9NTmHKGEPuw38BoNlVVyLR0TZHpJR/aXIKU+4P9rZ75BEbI1EqMDQ5haHC7dutdteZ39sYiVKBo8kpDK0b6ZpFEdupk5eeSoUvTU5hZudTk6y2PtirIpkmpzBiSkrY9+abALT561/0wV4V0TQ5hZHs3n2sdourr7YxEqUCT5NTmDi2bJnV7rHwVy89lYoMmpzCRM7lVwAQ06kT0Y0b2xyNUoGnySkMbBp7vdXuplMHVC2hySnEleTnc/QXR8XgTtOm2RyNUsGjySnEre7X32o3GjbUxkiUCi5NTiHs0LffWu2k5cu89FQq8mhyCmHb/nw3AE3OPReJsbeqi1LBpskpRLkXK+jw3LM2RqKUPSKtHPmpbmuHLxWRi4L1mfypaO9eTGEhAF2+/srmaJSyR6SVI88CUowx/YGzgFdFJOyWIl47dJjVrtu1q42RKGUfu5LTaBxlxHH+vNBDH6scuTGmADhRjrzC440xR51VWwDq4WMRzlCyZ8oUq60P9qrazK7kVKqcOOBrOfIOlR0vIoNEZAWwHBjnlqxCnjGG3S++BEDc3Xfrg72qVou4cuTGmF+APiLSC3hLRL41xhz3EF+NypEHgvvqlq1uvcXGSJSyX8SVI3cbf5WIHAGSgXIlcGtSjjwQ8tets9rd58+zMRKlQkOklSNPPHEBXEQ6Az2BnEB8AH/bcN75AEQ1bEid5h5vXipVq0RaOfJhwFIRWQJMB+4wxuwJ0meqtq1//rPV7rm43EmeUrWSGGP7NxrbpaSkmEWL7EkKpqiI7OS+AHR46UWanHmmLXEoZRcRWWyMSSm7X2eI2+xEYgI0MSnlRpOTjfJ+/tlq91y21MZIlAo9mpxstOWmmwFoOGQwUbGxNkejVGjR5GSTdWf+wWrHv/GGjZEoFZo0OdmgOC+Pws2bAUj4+GObo1EqNGlyssGalFOsdv2+yTZGolTo0uQUZPs//K/V1gd7laqYJqcg2/HoowC0uOEGfbBXKS80OQXRqqReVrvNhPttjESp0KfJyQbdMjPsDkGpkBd2q0SGs5jO8TQ58w/EtPW0koxSyp0mpyDq9t13doegVNjQr3VKqZCkyUkpFZI0OSmlQpImJ6VUSNLkpJQKSZqclFIhSZOTUiokaXJSSoUkLXAAiMhuYFMA3roVYGf1Fx1fx7e7+pAvMXQ2xsSV3anJKYBEZJGnqhI6vo5fG8avaQz6tU4pFZI0OSmlQpImp8B6TcfX8Wvx+FCDGPSak1IqJOmZk1IqJGly8hMRaSEiM0VkrfNncw99OolIhoisEpEVInJXsGNw9ntDRHaJSJafxj1LRFaLyDoRmejhdRGRl5yvLxORk/0xbhXGTxKR+SKSLyL3+nNsH8e/xvm5l4nIPBHpF+TxRzvHXiIii0RkWDDHd+t3iogUi8ilPr2xMUb/+OEP8DQw0dmeCEzy0KcdcLKz3RhYA/QOZgzO104HTgay/DBmNLAe6ALEAkvLfibgHOBbQIDTgF/8+Jl9Gb81cArwD+BeP/+9+zL+EKC5s322DZ+/Ea5LOCcB2cEc363fbOAb4FJf3lvPnPxnNPCWs/0WcGHZDsaYXGPMb872YWAV0CGYMTjH/hHY56cxTwXWGWM2GGMKgA+dcZSN623jsABoJiLtgjW+MWaXMWYhUOinMas6/jxjzH7n5gKgY5DHzzPODAE0BPx5odmXv3+APwKfArt8fWNNTv7TxhiTC44khONf6wqJSAIwAPjFrhj8pAOwxW17K+UTri99Ajl+IFV1/JtwnEUGdXwRuUhEsoGvgRuDOb6IdAAuAqZW5Y11DfEqEJFZgKfqBA9V8X0a4fhX5E/GmEN2xOBHnorvlf2X2Zc+gRw/kHweX0TScCQnf17z8Wl8Y8x0YLqInA48BowK4vgvABOMMcVVqdWoyakKjDEV/oWKyE4RaWeMyXV+ZfF4+ioiMTgS03vGmM/siMHPtgKd3LY7Atur0SeQ4weST+OLyEnANOBsY8zeYI9/gjHmRxHpKiKtjDH+eO7Ol/FTgA+diakVcI6IFBljPvf2xvq1zn9mAGOd7bHAF2U7iONv53VglTHmn3bEEAALge4ikigiscCVzjjKxjXGedfuNODgia+fQRo/kCodX0Tigc+A64wxa2wYv5vzdw/nndJYwF8JstLxjTGJxpgEY0wC8AlwR2WJ6cSB+sc/dy1aAj8Aa50/Wzj3twe+cbaH4TjlXQYscf45J5gxOLc/AHJxXCDeCtxUw3HPwXHncT3wkHPfOGCcsy3AZOfry4EUP/+3r2z8ts7PeQg44Gw3CeL404D9bn/ni4L8+ScAK5xjzweGBXP8Mn3fxMe7dTpDXCkVkvRrnVIqJGlyUkqFJE1OSqmQpMlJKRWSNDkppUKSJicV8kRkXjWOqSMie0TkyTL7c0Skldt2qoh85Y84lX9pclIhzxgzpBqHnQmsBi4/MQFRhRdNTso2zvV9lolIPRFp6FzjKtlDvzznz1QRyRSRT0QkW0Te85J4rgJeBDbjWKZFhRl9tk7ZxhizUERmAI8D9YF3jTGVLYA3AOiD4/mtn4GhwFz3DiJSHzgDuA1ohiNRzfdv9CrQ9MxJ2e3vQDqOh0Of9qH/r8aYrcaYEhyPYyR46HMekGGMOYrjIeuLRCTa+ZqnRyL0MYkQpMlJ2a0FjpUaGwP1fOif79YuxvPZ/1XAKBHJARbjeOYwzfnaXsB9+eIW2F8VV3mgyUnZ7TXgL8B7wKSavpmINMHxgHW8cT0JPx5HwgLIBK5z9o0GrgUyajqu8j9NTso2IjIGKDLGvA88BZwiIiNr+LYXA7ONMe5nWF8AF4hIXRwLrXUTkaXA78A64N0ajqkCQFclUEqFJD1zUkqFJE1OSqmQpMlJKRWSNDkppUKSJielVEjS5KSUCkmanJRSIUmTk1IqJP0/5Io8ljefSoEAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# X VS. Z plot\n", "\n", "fig, ax = plt.subplots(1,1,figsize=(4,4))\n", "\n", "# want to loop over all bodies in this simulation and plot their TOTAL orbits -- all times\n", "# x/y positions only\n", "# r_h[# bodies, x/y/z 3D positions, all time steps]\n", "for i in range(r_h.shape[0]): # looping over all bodies\n", " # now we want to plot (for given body) the x/z positions\n", " # over ALL times\n", " # r_h[ith planet, x position, all times], r_h[ith planet, z position, all times]\n", " ax.plot(r_h[i,0,:], r_h[i,2,:])\n", " \n", "ax.set_xlabel('x in AU')\n", "ax.set_ylabel('z in AU')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAScAAAEICAYAAAAdoDKiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2dd3hUVfrHP296AgQSCL2EntBBkKYIqAjY1rqKP0RFUUFXXd1VV8XCuta1o9hYBUXArthQmoB0pCf0FggkpADp7fz+uJPJTDKZTCYzcyfJ+TzPPHPvuefc+04y873nnPue9xWlFBqNRuNvBJhtgEaj0ThCi5NGo/FLtDhpNBq/RIuTRqPxS7Q4aTQav0SLk0aj8UtMFScRGSsiu0Vkn4g84uC4iMgbluPbRGRANdo+JCJKRJp5+3NoNBrPE2TWhUUkEJgJXAwkARtE5Dul1C6bauOArpbXYOAdYHBVbUWkneXYEVdsadasmYqNjfXI59JoNNVj06ZNp5RSMeXLTRMn4Fxgn1LqAICIzAeuBGzF6UpgjjI8RdeKSBMRaQXEVtH2VeCfwLeuGBIbG8vGjRtr/ok0Gk21EZHDjsrNHNa1AY7a7CdZylypU2lbEbkCOKaU2uppgzUaje8ws+ckDsrKr6WprI7DchGJAB4DxlR5cZEpwBSA9u3bV1Vdo9H4GDN7TklAO5v9tsBxF+tUVt4Z6AhsFZFDlvLNItKy/MWVUu8ppQYqpQbGxFQY7mo0GpMxU5w2AF1FpKOIhAA3AN+Vq/MdcLPlqd0Q4LRSKrmytkqp7Uqp5kqpWKVULIaIDVBKnfDZp9JoNB7BtGGdUqpIRO4BfgECgdlKqZ0icpfl+CzgR2A8sA/IAW511taEj6HRaLyE6JApMHDgQKWf1mk05iAim5RSA8uXaw9xjUbjl5j5tE6jcZmSnByyVq/m7E8/k71+PcWnTlVat8Hw4UTfeisNhg1FAvT9t7aixUnjd5Tk5XFq1izSZr3rVvvs1avJXr26QnnLGc/Q5NprEXHkiaLxN/ScE3rOyR8oSDrGgSuuQOXkOK0X3Lo1EcOGEt6nDyHtOyBBgZTk5lKYnEzutm1kLVlKcUZGldcL6diRjt9+Q0BIiKc+gsZNKptz0uKEFiezUEVFHL17KtkrVzo83uKJx2lyzTUEhIXV6DqFx49z5NbbKDjscJUEMff9jWZ3312ja2jcR4uTE7Q4+ZaitDT2Dj+vQnlYnz60/+B9AiMjvXZtpRRZy5eTdPfUCsckLIzumzfpeSofU5k46Tknjc8oSk1l7/kjKpR3XvwLIT5aQiQiNBo1ivjEBJRSpL7+unVuS+XlkdijJwBxCbv03JTJ6J4TuufkbVRBAYl9+tqVhfaIp+PChUiQf9wf83bv4eCVV9qVRQwcSIdP5ppkUf1B+zlpTOHEs/+xE6awnj2J27WTTl995TfCBBDWvRvxiQl0/u1Xa1nOxo0kxMWTtXKViZbVX/zn26GpUxQmJ7Nv1Gi7srhtWxE/fzoW0rYt8YkJZK1cxdE77gCwvuuhnm/RPSeNx0me/qSdMMV+vpD4xAS/FyZbGp5/HvGJCYT3Lev1Jcb3IMuB/5TGO+g5J/Sck6dQRUUk9upt3Q/r0YOOX31pokWeoSg9nb3Dhlv3A2Oa0a0S9wdN9dFzThqvkr93r50wxc7/rE4IE0BQdDTxiQmEdOgAQHHqKRLi4lHFxSZbVrfR4qSpMenz5nHg8ius+3E7dxDer5+JFnmHzr/8TPuP/mfdT+zZi8Jjx0y0qG6jxUlTI47edTcnn5kBQHj//sbcUmCgyVZ5jwZDhtB9y5/W/X0XXsTZ5cvNM6gOo8VJ4zZ7hg4jy/LDbDH9CWI/m2euQT4iICyMuISyJEFJd93NKTcXKWsqR4uTxi0S4uKtC2w7fPoJ0RMmmGyRbxER4hMTwOKrlfraayQ/+ZS5RtUxtDhpqk1CXLx1u/PPPxFxzjkmWmMu8Tu2Ez7ASESduWABx/75T5MtqjvUuXTkIjLDUneLiCwWkda++jz1AVth6vL7CkJ0pmRi531Ko7FjATjz3fckT3/SZIvqBqaJk01K8XFAD+BGEelRrpptOvIpGOnIq2r7klKqj1KqH7AImO7tz1JfSLBxFeiyYgXBzZubaI1/0fa1V2k0zhCozIULOfXueyZbVPsxs+dkTUeulCoASlOK22JNR66UWguUpiOvtK1S6oxN+wZUTNSpcYNDE26CoiIAOv/2K8EttDCVp+2rr1qHeKmvvsrZpctMtqh2U+fSkQOIyLMichS4Cd1zqjGpb75F7ubNAMR+8QUhbduabJH/EjvvU+synaSpUyk4erSKFprKMFOcPJ6O3Lqh1GNKqXbAp8A9Di8uMkVENorIxtTUVBdNrn9kr1nDqZkzAWj90ouE9+ppskX+T9y2rdbt/RePQRUWmmhN7aWupSMvzzzgGkcX1+nIq6Y4M5Mjt94GQJPrrqXx5ZebbFHtwdYPKrF3HxMtqb3UqXTkACLS1ab9FUCitz9IXWXPkKHW7VYzZphoSe1DRIjbvs26b/uUU+MapomTUqoIY8j1C5AALCxNR16akhwjHfkBjHTk7wNTnbW1tHleRHaIyDZgDHCfrz5TXcL2xxSfmGCiJbUXCQ6m06Lvrfvai7x66JAp6JAp5Un7cDYpL70EQLeNGwhs2NBki2o3J56ZQcY8Y2lP1zV/EBQVZbJF/oUOmaJxiaKMDKswtXnjdS1MHqDl9Ces23uHDjPRktqFFieNHaU/noDGjYkcM8Zka+oOthPkuwfU3+U+1UGLk8ZK0r1/s253W7vGREvqHiJC7OefA1CSk0Negp7HqwotThoACk+mcPZXI/NIpx9/1IH8vUB4715IaCgAB6+62mRr/B8tThoA9l1wAWCkbgrt1NFka+oucVu3WLf36Pknp2hx0pD61kzrdscvvzDRkvpBaWz14owMik6dMtka/0WLUz1HFRdz6q23AHR2Wx8R1qMs+Mbe88430RL/RotTPSexZy/rdsTACq4mGi8Rt2undTv9k09NtMR/0eJUjyk8XrYcsbvNYlWN95GAAJrecTsAJ//9b5Ot8U+0ONVj9o2+EIAG559PQC3KxltXaP7gg9btA1f+xURL/BMtTvWU7DVlfkzt39dRG82igyVjTf7u3aiSEpOt8S+0ONVTSkOhNH/4YZMtqd9E9O9v3U7soWNl2aLFqR5yetEP1u2mt95iniEawEgUUYoqKDDREv8iyGwDNL7n+EMPAdDassC3NlBSosg9W0B+dhEFeZZXbjFpx7KQACE0PIiQ8EBCwoIIbxRMVMsGhDUMrhWe7raJIhL79NUhaixocapnnPn5Z+t248svM9ESg5ISRcqhMxzbk0Hy/tMk7ztNQW6RT67deUBz+l3UjhYdI00XsW5r11iD+6miIiRI/zT1X6Cecez+BwBo/cLzPr922rEsdv5+jD0bT5KfXX0BCm8UTFiDYILDgggJCyQkPIjjezPJy3IvRvf+zSns35xSobxN9yguuLEbUS0buHVedwhs0sS6ndirt+49ocWpXpFjyaAC0PjK8lm4PEtRQTE7fj/Gn78eIee083mUpm0a0LZ7NK26NKZl58Y0aBzqNbuUUpxKyuLAn6kkrk0mKz2/Qp1juzOY99Q6u7KB42MZdGksAYHem6bt8vsK9o24wGqn2b05s9GRMKk/kTBLQ+82mzqVmL/d69FzqxLFnvUnWPvtAbIyKv7gASIiQ+g9sg3dzm1JZLNwj17fEyilSNqdwYpPd3M6Nddp3Qtu7EavCzyfIqv0fxTSpTOdFy3y+Pn9kcoiYZoqTiIyFngdCAQ+UEo9X+64WI6PB3KAW5RSm521FZGXgMuBAmA/cKtSKtOZHfVBnIoyMqyB5Dw1ZCjML2b9ooNs+fWIw+NdBzbnnHGxNG1Te6NpKqVIXJPM0jmV58mI7dOMS27vSVBIYI2vl7tlC4duuBGoP7Hb/U6cLCnF9wAXY6R62gDcqJTaZVNnPHAvhjgNBl5XSg121lZExgBLlVJFIvICgFLKqTNPfRCn0jtyaNeudPq+fJIb1ykqLGbNV/vZtiypwrHWXZsw9KrOtOzU2O3z+ztKKXasOMbv8/c4PN55QHPG3N6TgAD3h2Sl/6sW058gesIEt89TW6hMnMycc7KmFAcQkdKU4rts6ljTkQNrRaQ0HXlsZW2VUott2q8FrvX6J/FzbG9AHb/52q32lfUeepzXmqFXdSasQXCNbKwtiAi9R7al90hjSJd2LIv5M9Zbj+/fnMI7U41J9gsmdKfXiPJJrKsm+tZbSf/f/zj5zIx6IU6VYaY4OUopPtiFOpWlIy/fFuA2YEGNLa3lnJxRtrBUAl0feuTnFPLTu9s5ttt+VNx9cEtG3NCNkHD9PKVpm4ZMmzUagPTkbD57umwifcW83ayYt5uwBsHc/Nwwgl0c9rV4+J+k/+9/gDEcr6/ZWsz8dnktHTmAiDwGFGGkJK94cZEpwBSA9u3bV2VrraY0LVFpkLOqOHHwNF++sMmurHHzcC67py9Nmkd43L66QnSrBlahOvBnKj+9ux2AvOxC3vub4QV+wxPnVmsObu/QYfVm7qk8ZopTTdKRhzhrKyKTgMuAC1Ulk2pKqfeA98CYc3LvI/g/hcnJ1m3bIGeOOLIzje/ftA+dcs64Dgy+vBNSgzmU+kin/jFMmzWagrwiPn5kNQV5xQDWIeBl9/alQ8+mlbbvtn4de851NBioP5gpTtaU4sAxjJTi5QfY3wH3WOaUBmNJRy4iqZW1tTzFexi4QCmV45uP4r/sG2XcySMGV/5FP7AllZ9mbbcr+8sD/WnTvX4OJzxJSFgQd7xm+C79/tlutq84BsAiy03g8nv70t6BSAVGRlq3M+YvIOqGv/rAWv/CbFeC8cBrGO4As5VSz5amIldKzbK4ErwFjMVwJbhVKbWxsraW8n1AKJBmucxapdRdOKEuP60rffITt2N7hSURKYfP8Plz9p/7+n8NIqZ9I5/ZVx/ZvjypwtO+vz5+Ls3a2g/3kp98iswFxpRpXR7a+Z0rgT9RV8Upa8UKjt5p6LLtlzv7dD4fPbzaru6N0wcT3dp3yzU0sG1ZEisX2IvUHa+NICTMuIkopUiMN4bicQm76qzHuD+6Emi8TKkwxdx/H2B82X/9cCd7N5atJ7vivn60i482xb76Tp9Rbekzqi2/L9jDdovf2Pv3/077nk25/N6+dmJ07P4HaPv6a2aZagq650Td7TlZh3QJuziakM73b5RNdg+/tgv9LqrbTylrE0op3n/gdwotE+cAVz04gIgdy0l+9FGg7g7tdM+pnnFmseGLqhA+nb7WulYsslkYE54cQmCwjjPoT4gIU167gLPpecz51x8AfP3fzQSHRlGaPKq+LQbWPSfqZs8pIS6ezMhObB5QFkT/2kcG0iI20kkrjb+w5pv9bP75sHV/4KYXaDuyL21eetFEq7xDZT0nffuso2zvMdkqTM1jI5n69igtTLWIoX/pzO2vjrDubzznYVYfaGGiRb5H95yoWz2nwvxi3ruvLCb1lff3o22cnvCuzSybm8Cu1WXOtHe/PapGC4v9Dd1zqgecOHjaTpguWP2QFqY6wKiJ8Qza8B/r/jtTl5F5su77F2txqiPsXHnMuh6uecomRi+fRtzvS0y2SuMpmrVtyKjl91j3P31yLXs3njTRIu+jxakOsHRuAss/3Q3AxbfF02vXbIB6u5q9LhK7YD6CYvTyaVZn2cUf7GTJnLrpXgBanGo982esI8EyH3H9Y4No9JPO3lsXCQgLs27f8OgABl3WEYDEP5Lt4knVJbQ41WL+989VpB3LBuC2l88jpl0jMj+bD0BUPQ5SVtdJeuABzr2sI1fc1w8wAt598PffTbbK82hxqoUopXjnnmXknDGymkx54wLCG4bY1Wnx+GNmmKbxIs2mTgUg6zdjLrFdfDTXPWo85MrPKWLmXUtNs80baHGqhbx3/++UFBkuIHe9OdIaYVEVly19kAD9r61rNLtnWoWy5h0iuenpIdb99+9fUaFObUV/g2sZC55dT1G+IUJ3zxxptwzl5AsvmGWWxgfY3nBUSYl1u0mLCP5vhiFQBXnFfP7cBp/b5g20ONUifn5vO6eOZgHGUK58gseMOXMBaHTxRT63TeNbUl58yW6/cUyEdYiXcvgsSz7e5ahZrUKLUy1h65Kj7N+cCsDkl893Giy/9Yt1b/2VxqDBBcaSlvSPPqpwrHmHSC67ty8AiWtOcODPVF+a5nG0ONUCTh48w6rP9wJw09NDCGvoPA1TQLj/ZdPVeIY2/33F6fEOPZsycHwsAD+9u52z6Xk+sMo7aHHyc/Jzi/jiBWPd35jJPWnSwnH2k4KkY740S2MSgQ2rjlY6+IpONGhsPL2d868/qK3rZ00VJxEZKyK7RWSfiDzi4LiIyBuW49tEZEBVbUXkOhHZKSIlIlJhMWFt44MHDP+VuCEt6Tqo8lXpSZbHzJr6Q/Hp05Ueu+WF86zbb9+9zBfmeBzTxMmSUnwmMA7oAdwoIuVzF40DulpeU4B3XGi7A7gaqPVeabZB8C+8xXlap/w9Rt2md9zhVZs0/sPxRx51evyuN0datzf9fMi7xngBM3tO1nTkSqkCoDSluC3WdORKqbVAaTryStsqpRKUUrt99zG8Q9qxLLYvN+JK3/7K+VXULqPZvfdUXUlTq2k0ZgwAWcuc94gCgwO4dGofANZ+c4CCvCKv2+ZJzBSnylKNu1LHlba1mtL1UmNu70lohPMJcFsCQkKqrqSp1bR8+imX68b2aUZwqPFk9/37a9dgwkxx8mo68iovLjJFRDaKyMbUVP965LpqofFkDoGuA6uOfmjrkKep+1Q32sQdr5VF1CztjdcGzBSnmqQjd6WtU5RS7ymlBiqlBsbExFSnqVfJyy5k61KjU3jHKyOqqG1w5ocfvGmSppYjIoy/uzdgzGPWlqd3dS4deW3nk+lrABg4PpaQcNf+PcmPP+FNk9xCKUVhXi4FubkUFuRTVFBAUX4+RZZta1lBPkX5+RQWFFBcWIgqKaakuJji4mLyCnI5cvoQhzIPkpZzioASQRQEqLKOc8dk9xKBBoWGEhoeQUhEA0LDwwkJjyAkPIKwho0Ij4wkvGEjGjSJokGTaBpERdEgKpqwBg39LvtJQVISIW3bVlmvY9+yG/D8Geu5cXrl6en9BdPESSlVJCL3AL9QllJ8p206cuBHYDywD0s6cmdtAUTkKuBNIAb4QUS2KKUu8e2nc4/MlBzys41Jy8FXdHK5ncrPB6DpXXd61B6lFPk52eSePUPumTPknj1teT9DbtZZCnKyyc/JoSA3h/xy2wU5uShV/eGmQlESACWiUGK8RwiESbixH2CU15SifEMUszMzan4yGxo2bUZUy9Y0adnKeG/VmqZt2tGkRSsCAiv36neXlBdeoO2bb7pUd/LL5/PhQytJP55NUWExQcGet8eT6AQH+E+Cg9KQFxdP7kG3QS1dbleaPLPb+nUERladYaW4qJDsjAzOpqeRnZFGVnoaWRnpZGWkk52RRs7p04YAnT1DiU2kA1sCAoMIbdDA0vuIIDSiAaERRu/Duh3RgJCwcILDwvgzbStz935KcaCiKFBRHGB5D1QUBRjvrghPWGAY3aK7ERsZS0x4DEUlRWQVZpFVmMWp3FMcPXOUlNwU5ydREFQsBBcFEFJkvAcXBhBSFEBoYQChBQGEFgYSF9qR2IDW5GRmkJ2ZQX5OdpV/2+rSuHkLYjp0okWnLrTq0p2WXboSGlF1b/DAVVeTn2BEwaxOsk3bsCrTZo2uvsFeQCfV9HNsA9ZXR5hsKRWmwrw8Tqee5HTKSc5Y3o3tFM6mnyL3TEXnvcCgIBpENaVBVBRNWraiVdfuhEc2JrxRJOGNIoko3Y6MJDyyMcGhYU6HOG9sfoP3t5db49fccd2wwDBeHvEiI9uN9NmwSSnFymMrmbakYhiSUjaSAWy27vdp1oc54+YQGOC4x1FSUkxWWhoZycfJPHmcjBPJZCQfIz3pKJknkx22Kf3f7NuwplI7olq1oU1cT9r36kOH3v2IaNyE5g89yNHJt7v2YW24e+ZI3pm23LC3uKTC4nF/Qvec8I+e06x7llNcVOJyr6mkuJjMkydIP3aUxL8/QFZYMCWDB3E65SQ5pzPt6gaFhBIZ05zGzVvQqGkzGkY3pWFUUxpGRdMwuikNoqIJbxRZI2HIKshi6GdDq6z37kXvMqzNMNdPrBT8/Aism+WeYRMWQrfqjerXJq/ljsVVO7P+cs0vtG7Y2j27MJ6yZp5MJuXQAU7s38uJfXs4vieRkmLX/ZFCCovoNmYcnQYMIrbPAIJtwvlWRmnvqWFUKJOeG+62/Z6isp6TFifMF6f8nEI++PtKoGJXWynF2VOppBw6QMqhA6QePkj68SQyTyTbfYlDC4to3rc/jVu0onHzFjavlkQ0buKVHklhcSEDPhlQ6fHGoY359dpfCQ9ycSFy+kF4o5+HrKuCJ9Ig0LWBQ35xPsM/G05+cX6ldVbfuJrIEM8nLVUlJaQlHSEpYSdHdmzl8PY/KcjNrbJdUGgoPc4bRfyIUbTpFm8XC6ogt4j3Lcui/GFop8XJCWaL0w9vb+PQtlP0H9OOHsMbkbwngZMH95Fy6CCphw6Ql23EcEKEqJatadq2HdGt2xLdph2RIaFkTL6T4JKSas091ITZO2bz6qZXHR77YMwHDG7l4pMgpeDpJtU3oFl3OP/v0HUMhEZCSRFkHISE72HZs9U/X9xlcMOnLldfn7yeyYsnOzx2QdsLeOvCt6pvgxskxMWTExxE4MvPs3/TOg5t3Vxlm25Dz2fA2Mv5+tVkRIShV3VmwCUdfGBt5WhxcoJZ4lRUWMiJ/Xv48vlFlBQdJyz8FDmW+aCgkFBi2scSE9uR5rGdaB7bmWbtOlTotqe+8San3n4bqN7EqDsM/2w4ZwrOVCifM24O/Zv3d/1E/24BRVWE8qhGz8Zljm+B9y5wXufSV2CQY+FxxLbUbdz0400Oj229eSsB4r05ndIHIY7+72fTT5G4+ncSVi4j9fDByk8iEYyaNJE+F40lKNj1lQiepNriJCJnsfe6VsApYBnwsFIqzRuGmoGvxKmkpJjkPbtJ3rebw9u3kJSwgyKLG4AENCF+eH9ad4+nVdc4mrXr4NKj573nj6DI4uHuLXHq/XHvCmWDWg5i9iWzXT9JVb2kJzPB1z5EZ47DK/GVH3+q8lX/jrj151vZeLLi92jbzdu8MqwuFacuv68guHklTxtsKMjLZdfvy/jzp+9IP+7YUzwypgUjbrqFbkPO89nDCY/0nEQkCrgFGKaUus5z5pmLN8WpIC+XVfPn8OdP39uVR7duS/ve/di/JZSCvBiuvH8Y7XpUP3V46Rc0fMAAYue5PjRxhQeXP8jiw4vtyl4f9Tqj21djnsKZKFXzx+9VctLhxY6Oj1XTzuVHl3Pv0nvtymLCY1h6vWezo5T+71s8+gjRkyZVu/3bU5dRXJBCcPAmstMdh/UdfNVfGXz19QSHhNbIVmd4dFgnIpuVUpXPhNYyPC1OxUVFrJz3EZt++KbCsS6DhjDqljuJbGZ47JY+OXF3YrL0C9ryyelE3Xijmxbbk3Q2iXFfjbMrm33JbAa1HFS9Ez3VuGLZX96Bfn7uzO/I7vMegIueqtZp/jj+B3f+au8Y++UVX9Itqpv7ttlgvTH17UvsgvnVbp95ModPn1wLlH3/khJ2sGLuh5zYv7dC/fNunMS5V1zj8cw+HhMnEQkGNiml+njKOLPxlDgl79vNvMcerFDe9+JxnD/hlgrOdY6+HNXF2rVfvozglu75R9kydN5QsgqzrPt39rmTe/pXMwxL/ll4rtySir9+AvGX19g+n+JIpNzo7b244UXm7pprV7Z90nZ3rbKyZ9hwitPTAfeH9KU3x6nvjKowjMs5c5oVcz9k1+/2Pb6Q8HAmPPsKTdu0wxO4M+d0tYPiKOCvwCql1DMescwPqKk47V3/B9/99z92ZV0HD+OSu+4nNMJxWF2AxR/uZO+Gkwy7pgv9L27v1rVLxSkuYVeN5giUUvSZY3+/cesH9PXdsHWefZk/Dd/cobxI3b8DmlT/h1l+7m7LxC2VOnS6wvHHHuP0l18BNReni27tQffBld/csjMz+PHNlziyY5td+fVPPke7HhXnJKuDO+L0v3JFCkgDliul6tQyeHfFKeXQAeY+/De7suueeJb2vfq61L70i3HHayMICXPvyZSzJzaukpmXyfkLygLaPT3saa7u6ujeVAXlf8SPp0CQ9+YqfMqWefDN3WX7l/wHhlbuXV4ZH2z/gNc3v27dX3XDKhqHOuihucDpH37g+IMPAe7//1cu3MO2pcbkuKu99z1rV/H9q8/bld3+5oc0bl51eB9HeHrOaZBSqm5k7qP64qSU4otnn+DI9i3Wslteeafa3dyazjdBzcXpWNYxxn451rq/bsI6IoIr7+1VSnlhqu29JUeUFMMzNg8txjwLw6ofefRswVmGfVbmJf/rtb/SskH1h+SFKSnsG2G4Rrj7/y8qLObde40swdX9Hh7fk8hnTzxk3e93yaVceNvdTlo4pjJxcnlmS0R6iMgzIrIXSyzv+khRQQGv3HC5VZiuevhJHlywqNrCVFjgeEGtL8kuzLYTpu2TtrsnTB+US+JZF4UJICDQ/rMtfgzS9lf7NI1CGrHt5rLh0cVfXOzU+7wygjwQh6wmkQlad4vjwQWLuHiK8WRyyy8/8N+/XuaxeFFOxUlEOojIIyKyFZgLTAUudqRy9QFVUsLrE8uGO/d/+jWdBlTzCZaFpARjIrNtXPWiGnqSIfOGWLfdnqDNSYckm050XRUmW2w/45vuPbQWEbu/+cBPqv+T8rQfUkmJe6LS58JLmDa77GnhKzd45sFHpeIkIn9gxFMKBq5VSp0DnFVKHfLIlWshHz1opF8KCY/gwQWLCAxy36O2NBtrp37mROG0nZy1vYtXG1vfoPogTKVMt4kD5eipnotsmVg2NeDI2dWX1CRDcFiDhvz9s++s+0tm13xw5aznlAo0AlpgBG6DasbprksUFuRbvWrv+d+CGuGPSfgAACAASURBVJ/v8E7Dwb5D76Y1Pld1+T2pLND94msWu38H/vzWsu0n6syCAdcICIBxL5XtZ7v3+QMDAu087ZOzHIdW8SbNOzQCYOfKmiVmlYAA7v1oIWAM8WpKpeKklLoS6I0R0OZpETkIRInIuTW+ai0k4XcjDU/Xc4d5pDude7YQgEZRVYe48DS2MYxaNWzl/ol2flW27el1cLWBwVPKtl9yPXJpeWydW8d8OaYmFrlF+57GDTIpseZRQUPCIwgKNZ7QnjlVRdC/KnA656SUOq2Umq2Uuhgjhvd04DUROeqsXV2kNDJAg6jqLzFxhgT4dj3ZlpSyYUSNHAFtYw7Vp+FceaZ55qG17dA6t6jqkCiepF28Z7/TEZHGcqXCvCoWd1eBy0/rlFIpSqk3lVLDgPOqbFDHiBtuZELZ8ssiky2pGRN/muiZE83w/XDUL4mxWYqy/Qu3T2PbGz/3U98OTprHNvLYuZRSnEk9CUB0DT3I3Voko5Q6XKOrWhCRsSKyW0T2icgjDo6LiLxhOb5NRAZU1VZEokXkVxHZa3n3yOOwyGZlq77Xfb3QE6c0lY3/Z37M9DrHl66HWnHErIvcjPZZQwKDPLdW7rPp/wCgcYuWNZ7+MC2AsIgEAjOBcUAP4EYR6VGu2jigq+U1BYt/VRVtHwGWKKW6Akss+x7hjpmG0/yq+XNY+2X1F1r6E6GBHvLcnlD7hdpfGN7GnJC5nnJJWPDUIyTvSQTglv/W/GmdmbOY5wL7lFIHACy56a4EbGM3XAnMUYZX11oRaSIirYBYJ22vBEZa2n8MLAce9oTBkc1iuOGZl5g//R+sXvgJ67/9gmmzP6uRS4Ev8UpgwdYV/XyWJJzkka+2U1yiKCwuobhEUVSiKLa8PElbAihAkRMaSIkyzl/27tFLAfDydX259pyq88TVJ7Iy0nn3rput+3fOmuORwHVV9pxEJEZE/iUi74nI7NJXja8MbQDbifUkS5krdZy1baGUSgawvDuMwuVuOvI23eO5/c0PACjMz+O1m65i5byPXG5f59jwQYWiyR9vJPVsPunZBZzNKyKnoJiCohKPC1MgMJ+GfEUjsvKN6+QXlVBY7B1hAnjo863eObGJlBS7l86+pLiYhc/8y06Y/jbnCxp66KGRKz2nb4GVwG+AJ9dcOOpLlv9KVVbHlbZOUUq9B7wHxtq66rRt3Lwlf5//PYtef5E9a1ay/tsvWP/tFwwYfyUjJ072eLwbT+GVyIYrnodRj9oVPX1FT578bmelTQIDhCDLKyBACAwQAsXYDhCs24E2xwIDhAARAgJsjovwTp4iKCKY21p3JDAASx3781mvI2Jz3LBDxMH1LeWlfy0REAQRGB1XdcTJ2kZWRvWWzhTk5fLNC89wdFfZ096RN9/OOZf+xaN2uSJOEUopjwyLypEE2E7ntwWOu1gnxEnbkyLSypK2vBVQM2eLShARLr//YfKn3MNHD00jK+0Um3/8ls0/fkvTtu258h+PE9XS/bRBtZlJw2KZNCzWbDN8z9Xv16h5Zl5m1ZW8QPpx15KFHt6+hS/+/bhdWZdBQ7n8748QUIPQL5XhijgtEpHxSqkfPXztDUBXEekIHANuAMqHSPwOuMcypzQYOG0RnVQnbb8DJgHPW96/9bDddoRGNODOtz+iIC+Xb19+liPbt5CWdITZ9xkOevHnj2LULVMIb+j4cW1xcQmBJiQ2/OXQL1wSW4Ms7U9mloXfzTsDYZ5Pi1QrOF7mN0af62t0KtuwNb5k/xZjWiOsQcV5olNHDvHzO69z8oB9ZExfjBJcEaf7gH+JSD5QiDGkUkqpGn0blVJFInIP8AvG9MFspdROEbnLcnwWxtq+8cA+IAe41Vlby6mfBxaKyGTgCOCTWOchYeFc9/i/UUqxc8USfnnnNQASVi4jYaXhXd6qWxwX3z6NmA4dadqmIWnHskg9cpaWHd1fm+UuD614qGbiZDtEfL5d/XXErCqbixt8Mv4Tj5/TGYlrjCUzccNaoZRi/6b1LJn9Dllpp+zqBYWGct3j/6Z1NydJITxIleKklPKch1bFc/+IIUC2ZbNsthXgMKKXo7aW8jTgQs9a6joiQq+RF9Fr5EUUFRSw5ot5rP/WcM5L3pPInH+WBb4PDO3L/k1RtOxYjbRKNWTbzdusES+Ts5Jrtnzl1p/hf5aQK4k/Qtx4D1hYi1gyo2z74UM1OtV9S++zbveNcS1YoacoKUylKG81axceYK0Dz5CxUx+gx4jRPsvGUkql4iQicUqpRFvHR1uUUlVn8KvnBIWEcP6EWzh/wi0opdj1+1J++/Btazqo4vytrP18K2s/N+p3G3Ie/S+5jDbxPb32RbA975gvx9RsCUsHm/Tj8280J72TWZQUw8qXy/bD3ff1LS4pZulRS+DBfq5H1yxxY3lISXExe9auYuOirzl5YJ/DOjHtYxk5aQrte5mbJsBZmN73lFJTRGSZg8NKKWV+HmMPYUZSzZwzZ3jv3tcpzltXaZ2AwCB6X3gJPUeMpmWXbg4Fy51ImLbxwhsFN+KPCX9U03q7k9mnfqovwzvbMCk1/My2oVKqc7PIXr+eIzcbKaEc/f+LCgvZv3EdO1f8xsE/nXy/JYygsGHc88G9pvjsVRYJs9Kek1JqiuV9lDcNq69EREYSHD6c4PDhTJs1moK8XBJWLuPPnxeRlnQEgJLiIrYu/oGti+3DT8R06Ejc8AvoOrgs1KtSyuXelojwxJAnmLF2BmcLzzJ311wm9nBzzZ0ITP4NPrREw3yqcd0XKFthur9mWVRsham6S4qyV60GoAQ4ums7+zasZd+GNZxJdf6Aun2vvgy87Cpi+w7gTFo+nzyxBsDvnIl1OnLMS0deGkN80nPDaRhVcTnJ2fRT7FqxlF0rl5F+rOpAEEHBIXTsP5CO/QcS228AjaKbOa1v+8N47vznuKzTZdX8BDaseAmW/bts//FUCApx/3z+SGEePGsTxP/Kt6G/41TkrmD793915Ktc1OGiSuuWlBSTcvAAR3Zs5ciOrRze9qdL1+jQpz89R15El0FDHCbGnP/v9aQlZdGkRQQ3PT3EwRm8j0cTHNQ1zBKn9d8fYMMPh4gf3orRE117AqJKSji2exeJq39n7/o/yDntum9Mw6hoWsf1pE33HrTpHk9Mh470/aSf9fjtvW/nvgH3OTlDFWz4AH6wyds36nG44B/un8+f+OUxWPNW2f41H0Lva90+na0w3dbrNu7ufgcnDuzl5P69HN+byPHdCeRlnXX5fC27dKPLwCF0GTSE6DbtXO5Fl94gJz47lMim4dX7EB5Ci5MTzBKngrwi3r/fiEpZ06Sa4X370vz9dzm8dTMH/tzIoa2bq/XltiWqdVtiOnQkpn0sMR1iiW7TjsYxLQgIdMHR7kwyvBJnX1abJ8odpVN/aC80rNpTPC87i/RjSaQfO0rasaOkH0/i1NHDnEk56ZYpzdp1oF2vPrTv1Y/2vfqwv5/xrCpqwgRaTn/CrXN6IgNQTan2nJNNwyXAf22dMEsnyz1sY73DNlddSYkiwI3Ac0GtWlGUnEzu1q1ERDYm/vxRxJ9fcZrQiLOTwrHduziWuJNjibusc1vlyTieRMbxJPasWemyHQ2aRNGoaTMaNY0hot3zNNg8kwZBhTQIKiDikbY0CCok4skkgkL8f6inlKIgN5e8f8eSWxxMdlE02UUhnC0K5WzPWzn76mucOZXK2bRTFBcWeuSawaFhtOzclRadu9KqSzdad4unYbRrMbMaX+FeQoFdq8ovyPAvXHHC7Ag8bMlV97SlrF5mX/EG7XtEc2RXOlt+PcKASzpUu33U9deR+vobVdYTERo3b0Hj5i3o4UC8AH488COPrHiYxtnBRJ0JJq64LUOCepF+PInTVdztszMzyM7M4MT+Uk/i2IqVJrqRqNNUHAR9W7HE5dZBoaE0bdOe6DZt+eTkV5xuWEhmw0LONihiw8SNHgtbE9bXPb+oZZ8Y4U0uurV8pCL/wBVxysRwanxDRL4H/s+7JtUvRk+K56OHV7Pm6/1uiVPja65xSZxcYXyn8YztOJa+c/qS2aiQg+zhJ/Ywc+JMRrQdUWm7kpJisjMyOJt2irNpp8jOzCDndIblPZPsXUvJLgohpyiYEvNCiFWL4IAiwgKKCA8qpGGP0UQ0jqJR06Y0ahZDo6Yxll5iM0LCnM/TXPT5RZzM+c1IFQI0Dm3MthtWedRWd3zibKdznKUhNxNXxEmUUkXAVBG5BVgFmJdsrY7RoHHZ3TMvu9Dh+iZnBDf37Cr5AAlg+6TtzN01lxc3vAiUJUSYO24u/Zr3q9gmIND6Y3XMk8bba70h0/FQksteg4G3Oj7mLTbPhe8qydgbGgmPuh8qf/Ivk1l/Yr1dWU1Sj3uapR+7n77eV1Q5IS4idyql3rXZPweYppS6zdvG+QqzJsRLKX1q17JTY6755znVbl/TlOTOOGfuORSUFNiVPTzoYf6vRw060Kf2wltVzAy07g93LPPcRLpS8NGlcHi183pTlhvXdoMSVULfORWHWC+NeImxHcc6aFEzavJ/L50Iv/HJwUS3auBRu6qL2xPitsJk2d8E1Blh8gcGXdaRDT8c4sSB09VypvQFmyZuAuwffb+w4QVe2PACAJsnbiY4oJrOe826ljlqHlgOc66sWOf4nxWfknmL6+dCjyvcbr4ueR23L769Qvk/B/3TfefWKig8ccLttkd3pVu3zRYmZ9TDZGP+h4jQslMkJw6cYc3X+xl2dRe3zlOcmUlgE+/8oEuXVUz6aRKbU8qWVQ6YazzObtuwLT9e/WP1hbXTSHuP8iUz7NeseYMBN8Plb9SoV3bkzBEu/fpSh8c+HvsxA1q4l6bcVU7NfNvttt+9YYR5GXaNe98zX6H9nDB/WAeQn1vEBw8YPk9T3xlVrR95afc+evJttPiHb5weE9ISuH5R5fGLfr32V1o28NBEq1Kw6X+w6IHqtRvzLAy+y2MJP23n4Ryx9eatBIhvJvxL/+eBMc3ottJ1l4/jezP4+r+Gd7mZvk22aCdMJ/iDOAEs/M8GUo+cpffItoy4oVvVDSzsPmcgJdlGNENvzDtVRVU/WoAfr/6Rdo1qlsfMlyilmLNrDi9vdN6LW3LdEppH+D50b6k4tZ35Fo0udD1CUOlc05C/dOKcsbHeMK3aaHFygr+IU2F+Me/dtwKAu94cSWCwa3fhM4sXc+xvxrITM8TJlg0nNnDbL65NSX4w5gMGtxrsZYuqprC4kOfXP8/CPVWnuQqUQNbdtM5zqbXcpFSc4hJ2udzLTlyTzBLLUzp/6TVBDSbENb4jODSQzgOas39zCgueXc+Ep1xbiNno4ou9bJnrDGo5yC7sx5ydc3hp40sO6zqaRC4lLjqOCXETuCT2EiKCI9y2RynFkbNHWLh7IXN2zXHrHD9f8zNtGpZPDOQfuCpMSimrMI25vac3TfIYuueE//ScwPgSvX23EULrxumDiW7t2tMUb7oTeJLC4kIeW/0YPx38yWxTHPLW6Le4oJ3nQ+96EqUUifGGV7er/+/5M9aTdiwL8K9eE/hZz0lEooEFGGscDgHXK6UyHNQbC7yOESf8A6XU887ai0hT4AtgEPCRUqoSDzv/RUQYe2cvfn53B589s67ak+MlBQUE+PH6teDAYF4c8SIvjqg4R7U/cz8vbXiJ1cer8EWqIff2v5eJPSYSHmTOKvyakrlgQfXqp+RYhenWF8/zhklewaxhXWnK8OdF5BHLvl36KZuU4xdjpIjaICLfKaV2OWmfBzwB9LK8aiWd+zdHxHhI9euHOxlzu+sfJfWVV2nxiDcyeXmfzk06M+viWVVXrOeceMpY4iphYS7V/3T6WgA69GpKRKT/3rjKY9ZCpysxUoVjeXeUjc+arlwpVQCUphyvtL1SKlsptQpDpGo1U143hhZ7N6aQerTq0CdhvQwBS//oI2+apfEj2s+uOvH2vKfLwkBfdo9vEyfUFLPEyZWU4TVOOV6bCQoJZNxdhlf2wmc3UFzkPGV025kzfWGWxo+IGOB8mc3+zSlkJBsuJre8MNwXJnkUr4mTiPwmIjscvBysVXB8CgdlHpu9F5EpIrJRRDampqZ66rQepVO/GFp2MhaKzrpnudO6wS3qnD5rakBediE/v7cDMHyabBeY1xa8Jk5KqYuUUr0cvL7FkjIcwEnKcGfpyl1pX5V97ymlBiqlBsbExFS3uc+4+h9lyyB++98ul9oUZ7mWXlpT+8j4/PMq6yil+PDBMq9xf3G2rC5mDetKU4ZD5SnDrenKRSQEI+X4d9VoXycQEe541YiltHvdCXauPFZlm2P33+9tszQmceKJ6QAEta48GWqpKwr4n9tAdTBLnJ4HLhaRvRhP40pdBFqLyI9gpBwHSlOOJwALy6Ucr9Deco5DwCvALSKSJCL+GeavGoSEB3H9vwYBsPzT3RxNSHdYr+mddwKQvcqzwcw0/kfs/PkOy22H/3e+6d/+WlWhnTDxLydMZxzYkspPswzv6xumn0vT1g3tjrvjnKepXThztv3m1c0c221k46ks3Zg/UpkTZu2ImaoBjAnyoVd3BmD+M+tJP24/t2TrrKmKi31qm8b75G7bVumx79/YYhWmqx4aUGuEyRlanGoZA8Z04JyxRqzxz55ZV0GgSkl+7HFfmqXxAYeu/6vD8h/e3sYRSwC5y+/tS+suPgrS52W0ONVChvyls51A7d1YlhklerIREeD0N9+YYpvG+3T+5Wfr9vwZ6zm07RQAl07rQ/uerqWTqg1ocaqlDPlLZ2L7GAkFFn+wkzXf7Aeg+UMPmWmWxgeEdDBuTDPvWmpdM3fptD7E9naefr62ocWpFnPp1D70v7g9AJt/PszX/91sN+9UmFJt9y+Nn5I+9xPrtipR1qBxAH99/Nw6J0ygn9YBtedpXWXsWX+CX2eXOWiOWj7N6l6vn9rVDUqf0gWNHMdiLrOW3/riebVqMa8j9NO6Oky3c1ty/WODrPvLRs4kPyTSRIs03uBkzAA7Ybpr5shaL0zO0OJUR4hp14i73hxp3V897DmSWp+PKnG+YFjj/xQcPszqITPY2XMyAGENg5k2azSBgXX751u3P109IzA4gGmzRhPawAjTtafbDbw9dTl66F57yT6dz/vP7Sc/LBowFvFOfvl8k63yDVqc6iC3/3cEw/uXZel9++5lnDhw2kkLjT/y+4I9fPRwWVTQif/qVWsX8bqDnhCn9k+IV8aOHr1ZMeJ1u7Kpb49CAvwno7CmIlkZeXz86B/W/bDcUwxb92SdfbihJ8TrIYElRYxePo1ORTutZW9PXcauVcedtNKYhVKKL17YaCdMAze9wLB1T5polXlocarDtH3biI4Zu+ptuxXqyz5JZOZdS8k8mWOWaZpybF+exNt3L+PkwTMAtOgYybRZo4k8ewSArmv+cNa8TqLz1tVhGo0ui+UTUFTAtFmj2bcphV/eNyIkfvqkEfj+zjcuICgk0BQb6zunks6y4N8b7MpKIwrkbN5sLQuKivK1aaaj55you3NOUOa8B/YOmYs/3MneDSft6lYny7CmZqQfz+azZ9bZlY2f2oeOfco8vUv/dyEdO9L5px99ap8v0enInVCXxanw5En2XTASqOgtXlKimDVtGeW/Ane+eQFBwbon5Q3Sjmcx/5n1dmW9R7VlxF+7VahrTTm+cwcSWHf/H36VVFPjO4JbtLBun1m8mMgxY6z7AQHC1HdGU1RYzLv3rrCWl27f9MwQmjR3PxW4pozda5P57SP7m0P88FaMnhjvsH7SfWWhluuyMDlD95yo2z0ngKT7H+Dsz0aYDWePo4sLS3j3b8sr9KRGTYyjx/DW3jSxTlJSXMIPM8tiLZXSY3grRlUiSqWU9pravP46kZeMcVq3tuNXwzovpiMvjSceAhQA/1BKLS1/3vLUdXGCsi979y1/ElBFplilFIs/3Mm+jRWjGtz8n2E0inYt02x95fCONBa9tbVC+ZjJPek6qIWDFvbkJSZy8C9XAfVj4ba/idOLQLpNOvEopZSjdOR7sElHDtyolNpVWXsR6Q+cVEodF5FewC9KqTZUQX0SJ6jeF9726Z4trbo05vJ7+xEcWj+HHOVJO5bF/BnrHR6b+O+hRDYLd/lcpf+rwJhmdFu5soratR9/E6fdwEilVLIl79xypVT3cnWGAk8ppS6x7D8KoJR6zsX2ApwCWiul8p3ZUx/EqTgzkz1DhgLu3Y2Li0v47rUtHN+bWeFY4+bhXPvPgYQ1DK6xnbWJ5H2ZfPXyZofHhl7dmQFjOlT7nKqggMQ+RtrwuF07kYC6//TU3ybE7dKJi4ir6cgHV6P9NcCfVQlTfSGwSVlc6UMTbiJ23qfVax8YwFUPGgk+T6fmMu+ptZQUGze20ym5fPhQ2R1+2DVd6Hdhuzq3TKYwv5ilcxMcDncBOg9ozsW39SAwyH1BKRUmoF4IkzO8Jk4i8hvQ0sGhx1w9hYMyl7p5ItITeAGodCZRRKYAUwDat2/vokm1m/ZzPubIzZPI3ez4bu8qjWPCuXvmKAByzhTw+XMbyMoouwf88eU+/vhyn3W/z+i2DL6iEyFhtevh8Jm0XJZ+nMCxPRV7i6V0H9ySUTfHeSR8ie0opj56hJfHa98WpdRFlR0TkZMi0spmWOZWOnJH7UWkLfA1cLNSar8T+94D3gNjWOfq56rNNDj3XOv2yeeeo8Wjj9b4nBGRIUx6bjhg/Lh2rjzOinm77epsW5rEtqVJdmXdh7Rk0KWxNI4x31VBlShOHDzDmq/2kby/6ugN1/zzHFp2auxxO/aeP8K6XR89wstj1pzTS0CazYR2tFLqn+XqBGFMiF8IHMOYEJ+glNpZWXsRaQKsAJ5RSn3pqj31Yc6plMwvvyL5MaPz6u0nQapEkbAmmWVzE6vVLrZ3U9p0j6JNtyiiWkXU2CG0IK+IkwfOcGRXGkd2pVeaTssRjaLDuOSOXrTo6P3IoqUT4R2//Zaw7hWdMusq/jYh3hRYCLQHjgDXKaXSRaQ1hsvAeEu98cBrGK4Es5VSz1bR/nHgUWCvzeXGKKWcRvqvT+IEZT+CprdP9nm2lqKCYhL+SGbtN/spyPOvxJ9N2zZk+DVdaBsXZZcowhfsHzuOgkOHgPrhPmCLX4mTv1HfxCljwUJOPGmE4fCXH4JSijOnctm97iQH/kwh7ZjrvRtXCQgU2vdsSvse0XQ5pznhjfwj/rZtGvnYL74gvFdPky3yLVqcnFDfxAnKek+NLrmEtq+/ZrI19Rt3fdDqCjrYnMaOtrPeAeDsL7+YbEn9RhWUhVPu/OtiEy3xP7Q41VMajRxp3ba9c2t8i61fU0i7dk5q1j+0ONVjbO/UJdmen+PROCf/wEHrdvc/a+Z7VhfR4lSPsb1T7z6nwpBf42UOjB8PgISGEhDu+tq7+oIWp3pO3K6y5AenF/1goiX1ixP/+Y91O27rFhMt8V+0ONVzJCCAxldfDcBxH/s81VdUcTEZc+YC0PKp+plZxRW0OGlo/Z9nrdt6ctz7JPbsZd2OuuEGEy3xb7Q4aQDosmK5dTtn85/mGVLHyVi40LrdfYv+OztDi5MGMGKNh8TGAnB4wgRzjamjqOJiTkw3hnFNrru2yoik9R0tThornX/+ybqth3eex3Y412rGDBMtqR1ocdLY0fWP1dbttI8+Ms+QOsaxB8seNnTfVjG+uKYiWpw0dgRFRxNlGdalPP8CxZmVB1rTuEb+gQOc+cFw02j5zNMEhPjHgmN/R4uTpgItpz9h3S6NO65xD6UUB8Zfat2Puv56E62pXWhx0jgkLmGXdVvPP7lPaSgUqJ8RB2qCFieNQ0SErqtXWfcPWBw1Na5jK+rabaD6aHHSVEpQ06a0eeW/AOTvSiDtw9kmW1R7OP5YWR6PDp/N024DbqDFSeOUyPHjaXjRhQCkvPQS2escJ47UlJH59Tec/vIrAJreeScR/fubbFHtxBRxEpFoEflVRPZa3h2mmhCRsSKyW0T2WRIZOG0vIueKyBbLa6uIXOWrz1SXaffWW9btI5Mm2YX60NiTvX49yZasNiGdOtH8gftNtqj2YlbP6RFgiVKqK7DEsm+HJR35TGAc0AO4UUR6VNF+BzBQKdUPGAu8a8nioqkhtpO5B8aPpzA52URr/JO8hASO3DzJut/5Rx3loSaYJU5XAh9btj8G/uKgzrnAPqXUAaVUATDf0q7S9kqpHKVUkaU8DBeTcGpcw1ag9o0aTeHJkyZa41/kJSRw8Kqyhwb6yVzNMUuc7NKJA66mI29TVXsRGSwiO4HtwF02YqXxAHYCdcFIPcQDcjZu1MLkBbwmTiLym4jscPC6surWxikclFXZE1JKrVNK9QQGAY+KiMPHJCIyRUQ2isjG1NRUF03SQMUhXvYf9Td19unvF3H4/yZa97UweQ6viZNS6iKlVC8Hr2+xpBMHqEk6cmftlVIJQDbQq/wxy/H3lFIDlVIDY2Ji3P2Y9RbbH+GR2yZz6p13TLTGHJKfmM7xf/zDuq+FybOYNaz7DiidOZwEfOugzgagq4h0FJEQ4AZLu0rbW+oGWbY7AN2BQ974ABr7H2Pq62+wd8QFJlrjWxLi4sn8/HMAgtu318LkBcwSp+eBi0VkL3CxZR8RaS0iPwJY5oruAX4BEoCFSqmdztoD5wFbRWQL8DUwVSl1ykefqV4Sn5hA5OWXA1CUkkJCXDx1OVGrKiqy8/xuescddFmsc/95A53xl/qZ8dfTZK1cxdE77rDux375BeE961Za7azVqzk6+XbrfuzCBYT36WOiRXUDnY7cCVqcPENJbi67+w+w7ge3bk2XpUtMtMhzlF/8HLd9GxIcbJI1dQudjlzjdQLCw4lPTEAiIgAoPH6chLh48vfvN9ky98ndssVOmEK6dDY+oxYmr6O9pzUeJ27zJnK3bOHQDTcCcODSy4zyhF2IOPIQ5spOpwAAB1BJREFU8T9UURGJvXrblXVa9D2hXbqYZFH9Q/ecNF4hvF8/o4cRGmotS4zvwcG//tVEq1xj7+jRdsIU3K4d8YkJWph8jO45abxK3NYtFGVksHfoMADytm4jIS6eiEGD6DB3jsnWlaGU4sDYcRQcPmxX3m3jRgIbNjDJqvqNFieN1wmKiiI+MYGczZs5POEmAHI2bLDO5XTfvIkAyzyVryk+fZo9g4dUKO/41ZeE9ejhoIXGV2hx0viMiAEDiE9MIHvdeo5MKlu9v3vAOQBEXnoprV9+yevzUqqkhKNT7iR71aoKx7R7gP+gXQnQrgRmUZSezt5hwys93unHHwjt1Mkj18rduZND11xb6fFu69cRGBnpkWtpqkdlrgS656QxjaDoaOuyj4zPP+fEE9PtjttmLSkl8orLaXLttYT36lVhKFiSk0POn39y+quvramYnNHmlf8SOX58DT6BxpvonhO65+RvVNXLqQmdvv+O0K5dvXJujXvonpOm1hDes6fdQlpVXMzp777n5LPPUpKV5dI5AqOiaPH4Y0SOGaMdJmspWpw0fo8EBtLkqr/Q5CpHAVM1dRXthKnRaPwSLU4ajcYv0eKk0Wj8Ei1OGo3GL9HipNFo/BItThqNxi/R4qTRaPwSLU4ajcYv0ctXABFJBQ5XWRGaAf6QzUXbYY+2w57aZkcHpVSF5JFanKqBiGx0tAZI26Ht0HZ43g49rNNoNH6JFieNRuOXaHGqHu+ZbYAFbYc92g576oQdes5Jo9H4JbrnpNFo/BItTk4QkWgR+VVE9lreoxzUCROR9SKyVUR2isjTJtnRTkSWiUiCxY77zLDDUm+2iKSIyA4PX3+siOwWkX0i8oiD4yIib1iObxORAY7O4wM74kRkjYjki8hD3rDBRTtusvwdtonIHyLS1yQ7rrTYsEVENorIeS6dWCmlX5W8gBeBRyzbjwAvOKgjQEPLdjCwDhhigh2tgAGW7UbAHqCHr+2wHBsBDAB2ePDagcB+oBMQAmwt//mA8cBPlv/JEGCdF74TrtjRHBgEPAs85KXvpit2DAOiLNvjTPx7NKRsCqkPkOjKuXXPyTlXAh9btj8GKoRiVAalsWODLS9PT+S5YkeyUmqzZfsskAC08bUdluv/DqR7+NrnAvuUUgeUUgXAfIs95e2bY/mfrAWaiEgrX9uhlEpRSm0ACj187era8YdSKsOyuxZoa5IdWcqiTEADXPx9aHFyTgulVDIYP36MO2IFRCRQRLYAKcCvSql1ZthhY08s0B+jF2eaHR6mDXDUZj+JiuLrSh1f2OELqmvHZIxepSl2iMhVIpII/ADc5sqJ630McRH5DWjp4NBjrp5DKVUM9BORJsDXItJLKVWt+RZP2GE5T0PgS+B+pdSZ6rT1pB1ewFGmzfJ3YFfq+MIOX+CyHSIyCkOcXJvr8YIdSqmvMX4bI4AZwEVVnbjei5NSqtI/koicFJFWSqlky/AgpYpzZYrIcmAsUC1x8oQdIhKMIUyfKqW+qs71PWmHl0gC2tnstwWOu1HHF3b4ApfsEJE+wAfAOKVUmll2lKKU+l1EOotIM6WU03V3eljnnO+A0rzZk4Bvy1cQkRhLjwkRCce4IySaYIcAHwIJSqlXPHx9l+3wIhuAriLSUURCgBss9pS372bLU7shwOnSYaiP7fAFVdohIu2Br4CJSqk9JtrRxfL9xPIENQSoWii98SShrryApsASYK/lPdpS3hr4UZU9ffgT2IbRW5pukh3nYXSntwFbLK/xvrbDsv8ZkIwxIZwETPbQ9cdjPIXcDzxmKbsLuMuyLcBMy/HtwEAvfS+qsqOl5XOfATIt25Em2PEBkGHzfdho0t/jYWCnxYY1wHmunFd7iGs0Gr9ED+s0Go1fosVJo9H4JVqcNBqNX6LFSaPR+CVanDQajV+ixUnj94jIH260CRKRUyLyXLnyQyLSzGZ/pIgs8oSdGs+ixUnj9yilhrnRbAywG7i+1AFQU7vQ4qQxDRGZYRt3SkSeFZG/OaiXZXkfKSLLReQLEUkUkU+dCM+NwOvAEYzwKZpahhYnjZl8iGU5jIgEYCx9+LSKNv2B+4EeGDGEhpevYFlGdCGwCMNb/UbPmazxFVqcNKahlDoEpIlIf4xh2J+q6sWp65VSSUqpEozlELEO6lwGLFNK5WAshL5KRAJLL+vIFHfs13iXeh+VQGM6HwC3YKxHm+1C/Xyb7WIcf4dvBIaLyCHLflNgFPAbxoLTKMoy0UbjH9lxNeXQPSeN2XyNEWJmEPBLTU8mIpEYi6DbK6VilVKxwDTKhnbLgYmWuoHA/wHLanpdjefR4qQxFWWEdl0GLFRG0L6acjWwVCll28P6FrhCREIxAp11EZGtGNEk9gGfeOC6Gg+joxJoTMUyEb4ZuE4ptddsezT+g+45aUxDRHpg9FyWaGHSlEf3nDQajV+ie04ajcYv0eKk0Wj8Ei1OGo3GL9HipNFo/BItThqNxi/R4qTRaPyS/wf1QPXWRKUczwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Y vs. Z.\n", "fig, ax = plt.subplots(1,1,figsize=(4,4))\n", "\n", "# want to loop over all bodies in this simulation and plot their TOTAL orbits -- all times\n", "# x/y positions only\n", "# r_h[# bodies, x/y/z 3D positions, all time steps]\n", "for i in range(r_h.shape[0]): # looping over all bodies\n", " # now we want to plot (for given body) the x/z positions\n", " # over ALL times\n", " # r_h[ith planet, y position, all times], r_h[ith planet, z position, all times]\n", " ax.plot(r_h[i,1,:], r_h[i,2,:])\n", " \n", "ax.set_xlabel('y in AU')\n", "ax.set_ylabel('z in AU')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Attempting! to make and save a 2D movie" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "from animations_library import plot_animations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note: to use this animation library we have to use the `matplotlib.animation` sub-library." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "from matplotlib import animation" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "animation" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "# create a figure objection (but don't plot anything on it)\n", "#fig, ax = plt.subplots(1,1, figsize=(5,5))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before creating our animation on this figure object, we might want to make it a little more efficient by downsampling in time." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "stepSize = 50 # I'm going to create a image every 50th timestep out of the 8800" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "# downsampled position vector\n", "r = r_h[:,:,0:-1:stepSize] # technically this goes up until the 2nd to last timestep \n", "# can also do: r_h[:,:, 0::stepSize]" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "((7, 3, 8800), (7, 3, 176))" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r_h.shape, r.shape" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(8800,)" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "t_h.shape" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0.00000000e+00, 1.13649278e+03, 2.27298557e+03, ...,\n", " 9.99772701e+06, 9.99886351e+06, 1.00000000e+07])" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "t_h" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "# downsample in time the time variable\n", "t = t_h[0:-1:stepSize]" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "((8800,), (176,))" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "t_h.shape, t.shape" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAASIAAAEKCAYAAAC2QVjtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2dd1yV1R/H34fpxL0XblBx4jZXWo7STCtX2S/LbO/SSjNHmS0blplZWjlyz4Zp7j1SNEDNjVsURUDgcn5/PHAHcLmb514679eLF884zznnUe7nnvEdQkqJQqFQ6Imf3h1QKBQKJUQKhUJ3lBApFArdUUKkUCh0RwmRQqHQHSVECoVCd3QVIiFEDyFErBDimBBiVB7lWgohDEKIAfnZP4VCkT/oJkRCCH9gGtATaAAMEkI0sFLuA+D3/O2hQqHIL/QcEbUCjkkpj0spU4H5QN9cyj0HLAYu5WfnFApF/hGgY9tVgDNm52eB1uYFhBBVgH5AV6BlXpUJIUYAIwCKFi3aIiwszK2dVSgUttm7d+8VKWU5R5/TU4hELtey+5tMBd6QUhqEyK242YNSzgBmAERGRso9e/a4pZMKhcJ+hBCnnHlOTyE6C1QzO68KnMtWJhKYnylCZYFeQoh0KeWy/OmiQqHID/QUot1AXSFETSAOGAgMNi8gpayZdSyE+AFYpURIoSh46CZEUsp0IcSzaLth/sAsKeVhIcTIzPvT9eqbQqHIX/QcESGlXAOsyXYtVwGSUj6aH31SKBT5j7KsVigUuqOESKFQ6I4SIoVCoTtKiBQKhe4oIVIoFLqjhEihUOiOEiKFQqE7SogUCoXuKCFSKBS6o4RIoVDojhIihUKhO0qIFAqF7ighUigUuqOESKFQ6I4SIoVCoTtKiBQKhe4oIVIoFLqjhEihUOiOV6ecFkL0FUIcFEL8LYTYI4TooEc/FQqFZ9EtZrVZyunuaKmFdgshVkgp/zErtg5YIaWUQojGwC+AypyoUBQwvDrltJQyUUqZlXSxKDkTMCoUigKAV6ecBhBC9APeB8oDvfOnawpfQmZkcOPKJS6dPM6lk8e5di6OhMsXuXH5Esk3biBlhtVnAwsVpniZshQvU5bSlatStnoNylYLpXzN2gQEBubjW/y38faU00gplwJLhRAdgQlAt1wrE2IEMAKgevXqbuymwlu4cfkS0Vs2cHTXdi4eP+qWOtNSkomPO0N83BlOHdxvtVxQ4SLUa9OesPadqNYwAj8/f7e0r9Dw9pTTRqSUm4QQtYUQZaWUV3K5PwOYARAZGammcD5Oyq1EDv75G3tWLSX5RoLN8sVKl6F8aC3K16xN6cpVKVG+AiHlKlC4eHH8A3If2ciMDG4nJ5F49Qo3rl4mPu4sV86c4vKpE1w68a9F2dTkJA79tZZDf621uF6vTQda9R1AhVp1nH9ZBcK0BJPPDQsRABwB7kRLOb0bGCylPGxWpg7wb+ZidXNgJVBV2uh0ZGSk3LNnj+c6r3A76Wlp7P9tJVvmzSbDYMi1TFDhIjTo2JV6bdpTpX4D/Pzzb1RyM/4KR3dsJXrrRi4cO2K1XN1W7egw6BFKV66ab33zJoQQe6WUkQ4/p5cQAQghegFTMaWcnmSecloI8QbwCJAGJAOvSSm32KpXCZFvkHgtnnXffcWx3Ttyvd+wczcie99H2eqh+dsxO0m7ncKhDX+ya/kiEq/mGKQTVLgwPZ99lTqROZY+Cyw+KUSeQgmR95KceJM/pn+Wq/jUa92eDoMeoVSlKjr0zHXS09I4uHYNf83+Nsc9P/8A7h81jhqNm+rQs/xDCZEZSoi8Cykle1ctZeNPs3Lc6zDwESLv7Wd1HceXOX80ll+nfcy185ZLn5XrhdPnlTcpWrKUTj3zHEqIzFBC5B0kXotnyfvvcPnUCYvrbfoPom3/gfm6xqM3cTH/sOzDCaQk3rS4fu9Lo6jXpuA4DCghMkMJkb6cPxbL3LdesbhWsXZd+r42hmKlSuvUK+9ASsmeVUvZlG102OKefnQa8j+En2+7fyohMkMJkT6cPLifxZPGWFzr/MjjNO/VFyFyMxv7b3M17gxz33qF1OQk47X6be+g13Ov+uxoUQmRGUqI8pfThw6wcMJbFtceGDOJ6o2a6NQj3yI1OYmlH4znbPQh47XwDp3p+czLPjdCUkJkhhKi/OHq2TP88MpTFteGvj9VGfc5iSE9jZWffsC/e0w7im0HDKLdA0N07JVjKCEyQwmRZ0m7ncKsF0aQeC3eeG3IpE+oWKeejr0qOKSnpbHkvbGc+SfKeO3+0e9Ss2kLHXtlH0qIzFBC5Dk2z/2BXcsXGc/7vPwmdVu307FHBZeUxERmPjec20m3jNee+vZnioSU0LFXeaOEyAwlRO4n/txZvn9ppPG8cbcedHv8GbUInQ9cPH6Mn0a/aDxvcU8/Oj88XMceWUcJkRlKiNyHlJIVH0+ysIR+9vsFBBcpqmOv/ptsXzSPbQt/Np4/8eUsQsqV17FHOVFCZIYSIvdw/eIFvnv+ceN5r+dfI7x9J6fqSk5M5dDGOA6sO8PtpHS39K9YqWCa3FmNhndUITDYN7e7HeV20i2+/N9DxvNW9z3AHYOG6dgjS5QQmaGEyHW2LfyZ7YvmARAQGMQzs+YTEBRk17Mpt9LYODeWY3svebKLVmnctSqt+9QiqJCeUW48y4G1a/hz5lfG8+dmLySoUGEde6ShhMgMJUTOY0hP54thAzCka6OW7k88S+NuPfJ8JiNDsn3pv/y99rTN+sPbVaJBh8qUr1EcP3/nbGQyDBnExV7n8JZz/LvPttjd8VA9IjpXKXDrWUk3Evj6CdPW/oNj36Naw8Y69kgJkQVKiJzj2oVzzHphhPF85Dc/WnXMzDBk8NuMQ5w4kDP8BYB/gB93Pd6Qmk3K5psAyAxJzI7zrJ8TY7VM4y5V6fBg3QIlSqs//5CYrRsBaN6zD10eHWHjCc+hhMgMJUSOE7t9M6umfgBofmGDJ32S64f10MazbJyXe2CwXk83pmbjsh7tpyNIKTm8+Rwb58bmev+eZ5tQo1GZfO6VZzB3rwkqXJhnZy3QxSpbCZEZSogcY+2MLzm47jcAOj/yBC16WyRTIT3NwJw3t5F8My3Hs4PGtqZ0Zd/YQTt/7DpLPtqX43qV+qXo80JT/Px8e5SUfaqmx+6mEiIzlBDZz5zXnuXy6ZNATuvo5JupzHotZ0DMHk82onYz79o2dpT9a0+zbfGxHNdHfN6JwCDf3YGTGRl8PuwB0lNvAzD8s28pWbFSvrWvhMgMJUS2kRkZfDKoj/H8qRk/UaRESUDbap/1ak4BeuLTjgQVLlg7UYnXUpg9eluO609+3okAHxakFR+/x9Fd2nsNHP8hVeqH50u7SojMUEKUNxkGA58ONk2/XvhpKQGBgRjSM5j+7AaLsuVDQxjwRosCtbibG4b0DGa+tIn0NFMOtOKlC/HwpLY+++7bFs5l+6K5APR74x1qNW/p8TZ9UoiEED2Az9CC58+UUk7Odn8I8EbmaSLwlJTygK16lRBZJz0tjc+G9jOevzxvBcLPjz9mHuLoHtNWeJkqRXno7VY++yF0FoMhU4zNPhat+9YismeoXl1yiai//uCP6Z8D+RMN0ueESAjhj5ZOqDtajrPdwCAp5T9mZdoB0VLKa0KInsA4KaXNlAhKiHLHkJ7G1CGaCAUEBvH8j4u5fPomC9+3/Lca+UVn/AN9Kw6Ou7mdlMbMlzdbXBv2fjuKlSqkU4+c5+ju7az4aBIAvV94nbB2HT3Wli8KUVs0Ybk783w0gJTyfSvlSwGHpJQ2UzwoIcqJ+XSsUNFiPP3dPL57ZbOFu8WQd9tQskIRvbrolcQducayT0wZYKs3LMO9z/lewLd/9+5i2ZTxANz3+lhqt2jlkXacFSI9v/aqAGfMzs9mXrPGcOBXazeFECOEEHuEEHsuX77spi4WDKSUFmtCgyZ+x1dP/WUUoYhOVXhmelclQrlQpV4pnpnelYq1tNAbpw9fZdrI9aTcymnK4M3UbtGKPq9qUTSXTRlPXMw/Np7IX/QUotwWH3IdngkhuqAJ0Ru53Qct5bSUMlJKGVmuXDk3dbFg8PkjA4zHTXtPYd74ncbz4R/fQcdB9fXolk/R//UWDBnfxnj+3SubOfjXWR175Dh1W7al5zMvAzD/ndeJPxenc49MeP3UTAjRGFgK9JRSWs/1a4aampn45d3Rxkh/wSVfRAjtu6dmk7L0eir//ZLSr17l4gcfcGPFSqfrqDDmbUoNGqRbPOfvXtlsHBEFBPnx5OeddemHs+xYsoCtC34E4JlZ8ylUtJjb6vbFNaIAtMXqO4E4tMXqwVLKw2ZlqgPrgUeklDmNPayghEhj6y8/sWPxfACCSz6HEFoSw/6vtzBONTyJTEvj8mefcXXmdx5vq9SQIVR4czQin7JfxO68wJ/fm6Y3vrbAv/LTyRzZodmKZe2cugOfEyIAIUQvYCra9v0sKeUkIcRIACnldCHETKA/cCrzkXR7XlIJkWVmjeASTyD8igPw5BedCAj03Ic1+fBhTvYfYLOcf7mylOjTh5C77qJQRESeHwRpMJB2/jw3f/+dS59OhXTb8Yzq7dmDfzHPujdktzwf9n57ipUK9mib7mTaYwNJuZVIQFAwL/y42C11+qQQeYr/uhDdun6N6U8+DEBgsX74B9YE4JnpXT3SniExkSOR1o3lirRtQ9UvvvSIMBgSEjj35lskrltntUy9PbvxL+a+6Yc5Ukq+euov43l+jTbdgZSSTwbeC2ihf7s/8azLdSohMuO/LETmf1z+wS0ILNKJ0MZl6f20+9eDLr4/mfjZs3Nc9wsJoe7mTfgF5//o4PaxYxy/595c71X7dgbF7rjDI+1+8/wG0lM1q+weIxpRu7lv+OKlJCYybfhAAAa8NZEajZu6VJ8SIjP+y0I0b8xrnDsSDfhRqNSLRPYKpXWfWm5tI+7V17ixalWO63U2biCwQgW3tuUKqadP8+9dd+e4Xmny+5S87z63t7fs033ExV4HoNv/GlC/dUW3t+EJzKfxrkZ69EU7IoWbid2+OVOEtMXpOx6q61YRuvzVV0SHhVuIUMV33yU8JprwmGivEiGAoOrVjX0r1tU0LT0/ajTRYeHEz5nj1vbue6k59dto4vPn9/9wdM9Ft9bvKao3akLDTncC8MWwB3TpgxKiAkJaSooxsFlQ8YfpMrQBjbtUc0vdN379leiwcK58/oXxWtUvvyA8JppSDz3oljY8TbWvpmn9HTrUeO3ie+8THRbO7ePH3dZOt0cb0KBDZQD+mHmY04evuq1uT9Lj6ZeMx1vmu1eg7UEJUQHh82HaTpV/UEPuGNiWhnfY9ISxiUxLIzosnLiXXjZeqzB6FOEx0RTv1s3l+vWg4ttvER4TTdmnTamyj/fqTXSE+9bQugwNo3Yzzah25RcHiD93y8YT3sHTMzVP/Z1Lf+HGlfxNfKCEqACwbeEy43Hz3sNp2q26y3UmrFhBjNmHM6R3b8Jjoik9zH2pa6SUXE66TPTVaH765yeeW/ccPRb3oOfinjz959MsPLKQuMQ4PLGOWe755wmPiTZdyBTd5EOHrT/kAD2ejKBURc1lZt74ndxOdk8KJU9SuHgIXR/Tkmh++8xj+dq2Wqz2cZJvJvPV49q8vmLY8wx59y6X64wOswyiFRZ1EBEY6FKdF25doPui7i7VYc6UjlPoEdrDLWFKcttpsxApF5g2cr3x+Omvu/hEWJWPH7oHgDb9B9H+wSE2Sluids3M+C8J0SeDBiMzbuAXWI+XfvrEpbpSz57l324msSh+991U/WyqU3XdTL1Ju3ntXOqPI+wdupcgf/vyrlnjSPsOGK6a1nTCDh5A2JnLLS/MxchTtlzuJDnxJl8NHwTAcz/8QlBh+52hlRCZ8V8Roi9HzON2gpaC+KV5K/BzwUw/fu5cLo6fYDyvs+5PAqs4ts4kpaTPsj6cvHHSapmnmjzFiMYjCPBzPOSsIcPAin9XMHbb2DzL7X94v1P1A6TExnKir2lrv85f6wms5FrMZ0NaBtOf26DV16I8dz/RyKX68oP1P3zD/l9XElykKM9+v8Du55QQmfFfEKJDm+L4fdqTAPR5+R3qtnY+DOi5UaNJWGZaZ3J0WnIr7RZt5rbJ9d6YNmN4sL7ndtaWHl1qVZhW3reS0BKhDtcpMzKIadDQeF5txjcU6+haMLHLp2/yy3u7ARjwRiQVaoa4VF9+kDVFG/r+VCrUqmPXM0qIzCjoQpSanM7XT39FetIfALyyIKdxob2cevgRknZrH5AiLVtS40f7t24v3rpIt0U5d8/eav0WA8MGOt0nZ1l3eh0v/vVijuu/9f+NKsUc30U0Xysr/dhjVHj9NZf6t35ONNHbzgO+sV5knivN3r8xJURmFHQh+vLJddy+/ikAT3w5i5ByzrkTHO3chfQLFwAoOWggld55x67n0jPSafZjsxzX1z2wjvJF9HdtSMtIo/mPzXNc3zJwCyWCHfMDi23ZioybNwEo2qkj1b/5xqW++dp6Udao6J4X36B+W9vuMcqy+j/C5gVHSE/RIqKUKF/BaRGKDgs3ilC5l16yW4Q+3ftpDhHaOXgnUcOivEKEAAL9AokaFsX+h/dbXO8wvwMRsyMcqqv+7l2U6KulXbq1cRPn3nzLpb6N/KKz8fj4394fSfSxqZrwZhnLegolRD5E2m0DB9afwZCiRVgc+v5nTtVztFNn43GlSRMp+6TtXOnpGelEzI5g1qFZxmvrH1hP1LAoigR6Z4jZAL8AooZFsemhTRbXI2ZHcC3lmt31VP7gA0oNHgxAwpIlXF+6zMYT1vEP9OOOh+oC8Ov0KI/YSLmTUpWqGEO0ZGUD9gRKiHyIGS9sJD1lKwClK1elkBOhLS59/DHpFzUfqMofTqFk//42n1l9fLXFKKhDlQ5EDYuiXBHfCMlbqlApooZF8Xyz543XOi7oyMQdE+2uo+LYMQSHhQFwfvRoUqKdtzMyd735/o2tTteTXzz+hRbYbu2MLz3WhhIiH+FMdDwAhpRdAAye9LHDddw+epSr384EoNTQoZS4N/dwGea0/KklozaPMp5vGbiFr7t97XDb3sATjZ/gwCOmtHgLYhc4NFWrtWyp8fhEv/tJv2b/qCo7j3+q7cIl30jlVsJtp+vJD0LKmr5wsqI6uhslRD7Cis/+xnBbC00aXKQowUUcCzImDQaO32tKMV3xbdtrHRGzI0gxpBjPo4ZFObzY6234CT+ihkVZXIuYHWH3FMnctOFo23ZOT62CCwcYt/B/8IFR0f8+nQ5oIWY9gRIiH2D36hMApCVpc/Shkx1fG4ppaDKis8dOyHyk8Fyz53J8eH2dqGFRvN36beN54zmNMWQY7Ho2LNoUqzomvIHTfej/egvj8YXjCU7Xkx+UrlzVeHzl9Em316+rEAkhegghYoUQx4QQo3K5HyaE2C6EuC2EeFWPPnoDu1aeIMNgcj0oWcGxgFvm9jBhUQdtljcXoe/v/p4RjW0vZvsiD4U9xJp+a4znTX9sSobMsPmcEIJ6O3cYz+Pn/OhU+0IIWvYOBWDxlL1O1ZGfZOVFmzfWNXuq3NBNiDJTTk8DegINgEFCiOxfL/HA88BH+dw9r2Hb4mMApN5cBEDv5x37I0g+YFoTCV240KbzqrkIzb9nPpEVHTYJ8SmqhVRj+6DtxvMmc5rYNd3yL1GC0v/7HwAX33sPmWFbwHKj1b2mwHVnY+KdqiO/qNuyLQCpyckY7Ehg4Ah6johaAceklMellKnAfKCveQEp5SUp5W7At9JqupH9a09rHwypxbQJa9/JoedPPqRZOAdWqULhiLx9nMxFaG6vuTQs0zCP0m5ASrhwCH7sB+NKWP/5ugOc3AJOfthtUSyoGNsGmbJVNZ7T2C4xqvDG68Zjc5cQR2l3v+Y+sXzq307XkV/Uba05Mrs7eJovpZzOk4KYcvrwZi0Tp+G2NmyvEubYH/uZJ0caj+us+zPPsuYi9E33b4go55jhn93E7TMJzLslYXp7+Hd93s9cjIIfesP4UqZn49w7lSkeVJyV95mSPjaeY1+gtPoHTOIR/+NPTrXd7C5T/KhrF7w7iNrdI18AYM/KJW6t1ydSTttDQUw5veHnWADSkzWDvHteeD2v4hZIg4HEjRsBzWkzL5YcNf1Rvdj8RdpVdnP4jrQUk4B828U9dX7b1VTntZNuqTK0RCiTOkwynvdd1jeP0hp+wcGUekRL3XRx0iSnd9HCMmNdzx2300ZJfTHfrb113Xnzhew4FyvBPZwFzIMqVwXO6dQXryPhchIAUprm4sVKl7H7efNdsrw8xzNkBu9sM7l3DI8Y7kg38yblBkzOI2728/uhtAPB/VMSYHIu0Sc/a5JZ399QuqZjfcxGn9p9+OHwDxy9dpTjCcc5fv04tUrm3ceKb77JtcwF6yMtIqm/z/HRWtdh4cTs0FxuMjIkfn7e6xDbYdAwtsybze/TP+P+UePcUqeeI6LdQF0hRE0hRBAwEFihY3+8ip/GaLsyRYtpi82OrA3dPn7CeFx/b97Ov03mNDEeu22LXkptpJKbCI2Og3EJ2o8jIgRQqITp2Xeu57z/eVOtXRfXkpb0MY0Q+y63PSoCqLVmNQAZSUlOtWnuib/0o31O1ZFftLz3fgBO7HefY7luQiS1r/pngd+BaOAXKeVhIcTIrLTTQoiKQoizwMvA20KIs0II7w/k4kbiz2rTqy6P2r+FfrxXL+3A3x+/otYNH81dHHYOdtOU4FK0tvaTnSwBCdbcUjIyJClpBjbEXuKtpVHsPXWNaX8dY+GeMySl2rEjI4Spzp5TLO+NLwUHf3HpNQ4+YjJzsMf6OriWSVSPdnRsQyGLRye3B7zfpsjP35SyPOmGe/qq59QMKeUaYE22a9PNji+gTdn+U2QtUpuvNxQJsc+iOT3etAUcdsj6CCfNkMaCWC3yXqeqndzjuDpvEMSusbw2TvtDzciQNJ+4lutJuW+A/rzztPH4tUWWtk7Ln2lPk2q5iFsWrZ/Ufj6qB4mZucSWPKH9jHPugyKEYEL7CYzZqsXjiU+Jp3Sh0nk+U3PFck706Uv6JecyYBQtacqMe/1iEiUreKczMUDLPv3ZvWIx2375iW6PP+Nyfcqy2gvJWqSOuEPzQSpbPdTuZ4+2a288zivwVvOfTPF6vrzTDc6M40pYitCja2BcAov3niV01GpqvbnGqgjZou+0rYSOWs1D32zPe7T06hF460LOfjm5gHxfHVPI2E4LbI9yCtWrZzw+0X+AU21Wb6CJ3c/v7LBRUl9a93sIgANrf3VLfUqIvAzzUdDRndpaRdf/PelwPfX3W19nuJRk+sZeO2Ctw3XnYFy20do715l5phKho1bzysIDuT7SsV459o/pzsnJvS1+jkzsyaPtQnN9ZueJeBqM/Z1BM/L4kAYWzjkKerek02JkPmU13120Ro2ftEXrlMPOpSXq9ZT78qt5kuAi7h2t6To1U+TkyE7TN3rCJe24arh9wdbjXjNt7/sVtp6//M6FdxqPKxZ1MT97NhFKevMqDUavybXosUk9CfDP+7svKMCPcX0aMq6PZjMlpaRmtvq2H79K6KjVvHZ3fZ7pYiWW8rgE+DgcbmZuxL5bEsZeAwcTDJhPWd/Z9g73170/7/KRJkv0hJUr7YpwYI5/oKl/txJuU7REcB6l9aV42XLcvHKZc0eiqVwv3PYDeaBGRF7Gnz9oDqkte9cwXrM3tvGNlZpBXpVPracVSjOYpkf7hrq4O5NNhHY+cpwGY3+3uFazbFHjaMeWCOWGEML4/Lh7LT2APvw9ltBRq63b7rwSDXXN8ryNL+Vw+2C5cL3u9Dqb5Yv36AHAudfst/syp0R57Utk+af7bZTUl7YDtJRDO5e6tjEASoi8lpLlM4OX1bfPuzvtnMkEK6RnT6vlzNeGAv1dSJq4z9LRs7n/Ih7KNmU6Obk3f73a2fk2svFo+5qaoGWzsak5eg2p6Va27IcsBH+zUUX2aaQdmH8R5BacPzt5fRHYQ98XtSB01y44ZwqQX4R30IxTj+/b7XJdSoi8iPRUUxiKA2u16UiLXn2sFbfgWNc7bZYx9yz/9X4XFhnTU2HFs8bTu0usJP5WqvH8oweacHJyb+frt8Gx93qx603L96339q+kGayI0Zhsu1gXHLeXMh8VJacn51nWXLgSli93uK3ipQs5/IweBLiY/dccq0IkhFgphFhh9rNcCPGdEGKo21pXWLBrpckQ8XSU5sNUq0Vrh+qos2mj1XsvrH/BeFy1uAtWERNNLjTDKy0h9uJN4/mBsXcxoIXnLS7KhxTixPu9LK7VfetXMjKsTNPMF7Cnd3C4PXNxafVzK5vli3XRRgvn3sgR3cYhLp++abuQF5Cakrc42yKvEdFHwMdmP58Aq4B+QgjPhGn7j7N/rWZL03GgaRvY0W+dwPLWM2lsOLsBgH51+jneuSyiTfmtUss0YN0JUwTHmAk9KFHEfd+StshaPzKn1pu5L5QD8JzZmti4POySrDCl4xTbhTKp+uUXDtdvTvlQzW73zx/+sVFSX2pHal+Ux3Ztt1Eyb6wKkZRyYy4/S4EHgR4utarIk/C2ju1kXZs/32aZ9AyT/c277d51uE9GFgwxHtaLM0U43DemO4UC/XN7wuNkHxm9tdTK1KtMbbMTqU0xHaBnTdPaW0x8TJ5lhZn1cUpsrEPtAHR9WAvUH3/Ou73xwzt0BiB6q/WRuD04vEYkpbQvnqbCaRIuawvPJcpXsKv8hXGasBSKsO6K8Ow605qO0xlGj5v+2AalmmJeP9OlNqWLBjlXpxsQQrB9tClZobmVdg7MfdQmOh+l4YGVD9hd9kTf+2wXykaZKo5naNGDms00c4WTf7sWliWvNaLSufzUFkK8CzhnraWwSprZQnWWM2HWf7K9VP9+ltV7W89pAdpbV3JszcmCOaaF8+0ZpthIr90d5nydbqJSCUu7qdBRq3MvmF2EHTR0/LDTh3aXrfK5c3nnfImgQtbt1RwhrxHRXmBP5u+s4/lAEeApt7SuMBK7w2TIeNxJIfK3I8+Z06mAbpliZo9Ke9x4nH1apCfZ14vMdy/EdmwAACAASURBVPIsMB8V/ebYYnKPUPtXJYp37+5Q3dZIvunYFNIXyWuNqKaUslbm76zjllLK1wDXlsgVOTi4XgtWWb1haeJitAVKe6xVU0/nMQ3JxDw7RaCfk4vJH5q8y+cbTNMgp6d5HuK7YSbxbj7BivuKeZ93Ts+9jB0cu3Ysz/vu+reJ3nbeLfV4M3avEQmNrkKImWhBzRRuJMt4LbxdZTIM2sJyoaK2RziXPvnUZpnFRxe71jkreNJWyFnuDLdcV7NqdX3Pp+aFnGqr3wr7dx/TMrPrOkKNRlogvNidF2yU1JeslNSubOHbFCIhRGshxGfAKbTAZZsB/RcFCihV6jm2rXzzNy3XmV+I9TBNE3ZMcKlP5oHG7rr9gWt15QPli5ssqX89ZOVDHPmY6fis+wJ8WePSB/Zv/WdRp4VmiuHtO2flQ7XdSFfyneW1WD1JCHEUeA+IApoBl6WUs6WU7gtWq7CgUDHnpk6lh9q2M3U6P1nUQuPhEalFXZzxcAtrpXVnx2iT1fXTP9vhT/ddN4fqn3nXTLvLBlbR8kHcWJOHfZMVakTYHxpYT8pU1f4m4uOcnyjlNSIaAVwEvgZ+klJexYXg9gr7cHZdoWT/vL3CAYaGO2kUvzSngN3V0EWvfQ/i6XjPrSratqzOouQDzsUlAihcTD+TCEcoWaESANcvOj+FzEuIKgKTgD7AMSHEj0BhIYQKHeKFZH3z5kWpQs55nxdYqrVx6jFHviwKN2liu5CPk5XUIfHaVRslrZPXrplBSvmrlPIRoA6wHNgGxAkh5jrdohl2pJwWQojPM+8fFEI0z60ehcIp+ts/xXKW4DpW4iUVIApnhjFOdiF+tV27ZlLKFCnlIillf6AuWsB7l7Az5XTPzPbqok0VnTSCUfzXGNzalHZo81ErCTdL5pHqyE34lXA87IivkRWt8baTGUzAORePG1LK2U63aMJmyunM8zlSYwdQUghRyQ1tK5zgUEaoR+u/Ojeaq3Oj3VJX7wjTn0lSqn5eSSKg4K9k+AdoGyyGdOczw3t7ymm701IXxJTT3sYN6dmsEunxKWQkOv/HbE4Zs4XeDnXK2n4gzTkbGHNn4twQDoam9UUCgrR/66DCzv99eHvKabvTUhfElNPeRl0/z9qxln+6KeVGuCd4/PHLJtubdTFW0vuYJ2IMdM5nSuT6J/rfIisNdZYnvjPYNW4UQrQDQs3LSynnON2qhj0pp1Vaai+inLjh0fqFG7fdRy02RVRsZi0nWqLj1s6OYtWyuwBhSNdGhf7+zoeBsSlEmdv2tYG/gazJtgRcFSJjymkgDi3l9OBsZVYAzwoh5gOtgQQpZcF3vFG4zI0U05SpSkkro52NrluJ+/vl/eFLPZa3P1pBIDVZW6QOciHFkD0jokiggXSztEsp04UQWSmn/YFZWSmnM+9PR8sC2ws4BiQB/3NnHwoS0mCwCMaVG8npyRQOcE/YBl/CqoHj3u893vbNvzZ4vA29MQqRh9eIDqEZN7odKeUaKWU9KWVtKeWkzGvTs9JOZ+6WPZN5P0JK6XmnIJ1xVu9vHz1qs8y3B791qm7aPZ/jUswFz07TvJnbhtt2l01cvx6AIi1bOtxOVowqLwtwkIPb+SREZYF/hBC/mwfTd7pFRZ6kpji31XzpY9spbL6NclKIOr1hPKwmtHWVHlM3O1dXPvDMXJN/2W8v3mH7gfscCwUyccdEu8sm/60lQSjWtauNkjmJj9MW3MtU9e5ojSk3tQD/9kSLsIY9QjQOuA/N+dU8mL7CA1w87ph1ala2iFubrQvDqFauZZIg2PQHtjn4JdfqygdWHzQtI4ZVtBKV4Mwu03HTQQ7Vv+zYMof7VLR9O4efuXJW+4CXrVbc4Wfzk8R4zbWjWOnSTtdhU4isBNF3LVK2IgeFQzRbjNhdJsfBtNsp1oobqTTJ9rfzfXUcj5lsD/dN2+qRel3hyEU70+9853r0xPn32E5akIUzrh7njmqRJCvWtB7ixRu4GX+FIiVKGg0bnSGvMCBbMn/fFELcMPu5KYSH93H/gzTurNlpHtl5kUp16gNw/ugRm88F2PEtVDSwqPE4LcNJg8EXTZkxmgutX3+fuW6ttG7c9ekm43HMBCthXc3X4drbztxqjfDS9ud7d8aw8fTheACqNXB+pJEfJF69YnR8dZa8nF47ZP4uLqUMMfspLqX0bon2QcLbVzYehzbVYv2cPOBYZgSZYSXTqRlf7v/SsY5lUdLku7UkeJzx2GqQeh34dK2lcFtNb/SumV1Rt3EOtXH0mmlTwE/kLS5Z60POknJL+9Lw9syvCZcuElLWNSPigm9/7iMULWGKKlizmSZEJxxM0XJ94SKr90oGax++WYesZ/qwSfNhxsNymEZDh+Kc97p2F4YMyWfrTCKxf4ydUy8Ht6TuX2E77lMWZ0ZqOSaKdrRjwTwPvC0uuDmG9DSuXzxPmarVbRfOAyVEXkjZ6lqgentDbwaFhgJw4Z13rJb5Y8AfrnYL7jWlx9ld6Gnj8T1fbLGe6jmfqG2W4bVWuaKUspZnbZyZN/ybztvGzrrbtqAbrmtiXXnSJIfrzxoNeTvXL5wnw2CgdBXXIhkoIfJCzvzj2NpL6KKFNsuYGzJ+/beT0VSEgC6m7K4L628wHueZ6tnDZJ8ernu5U+4FU7KN3IIcs3u5cMu0kdCyov12QQFO+D4e26v5x9WNtJ5C3Bu4Gqf5pJfxtBAJIZ4VQqjQfvlAncw/uj++M+WvtGfdxzyfWdolKw6eQGQFLdXOVwe+craL0Ok142HLUzPwx2T3FDpqdb77VmUXoZgJPaxPZSabTR/GxjvcVvdF9u+0ubo+9PefWpqoRp2rulSPp7l86gRC+FG6imv9tGdEVBHYLYT4JTOiovdOWH2cTgO13TJDWgYVamnbvacPHczrkRwc62hlNIDldOJWmguZIcx20P4t9LDFrZqj15CQ7PlphZQyhwj99uId1heozadkfaeBDR+x3NrLYtugbTbLnxyo2SaF9LnXoXaySLikhSWpVNu7A6udPxpL2RqhBAa7tqBujx3R22gREr8DHgWOCiHeE0LUdqllRQ7MM3g07qZtPe/7dbldz9bdvMlmGfPvkDZznYvXDGg7aC2fMJ6eLGTpq9zk3T8YMnOH8/Xb4Oedp6g52nIqeGDsXdaNF395xPK8meNJBO5bbrLFKh5kv4Fh5fffd7it9DTTKNObv/dlRgbnj8ZSuW59l+uyN1SsBC5k/qQDpYBFQgjHkzUp7CKoSCMAju/bbVd583WItLg4q+UW3LPAeJwhbU/7rNL7I4vTk4UGE+Rv+nPaeuwqoaNWc8CNtka3bqcTOmo1by09ZHE9ZkIPShSxYky3eyb8Yybm7zjXn+MJxwF4tumzNste+dbkSmPLETk3DqzT1l18YX0oNTmJSnVdT3NozxrR80KIvcAUYCsQIaV8CmgB9He5BwoLWvYOBWDdbNvGjNY4dqf1PF0NypjCgjeZ42KGiXGWi79HAgfy58sdLa71nbaV0FGrmbDqH6ebOXj2OqGjVtPwHctQ6cPa1uDk5N7Wp2M7vobVr5jOx8Y75UE6evNo4/GTTZ60Wf5ypt9ftZnOBeffsUwTvbb3e3fg/bgYbS2zcj3XhcieMCBlgfullKfML0opM4QQ97jcA4UFLXqFsnv1SQCCChcmNTmZy6dOUK5GTZvP1tuzhyOR2oK0lNLqsH5Ozzk88qs2XUlJT6FQgAvz+3euWxgI1vmqKifev55j6vTdlhN8t+UEAGWKBrH+lc5WRzGp6Rk8N28fvx+2Hrhs7UsdqVshjynSlNqQdMV0/sYph9eFslh1fBUAQ8KH2CybHm9aBC/Wob3DbZmvRXm7IePJA/sIKVeekhUr2y5sA5tCJKUcm8c990Q6VxjxN5veNOs5lJ1LvmXDnJk8MMa2LYp/MZMrx7Gud1L3r/W5lmtWvpnxuOXPLYkaFpVrObsQIocYiXdLcnLsCTIKlcp1W//qrVSajHfOrmnJ0+1oXj2PTVwpLS2nAcZeAydjR0fMjjAe2+M8fLSdJj5Ztl2OcuJvTTxDynq3CBnS0zl96ABh7Tq5ZR1L2RF5IQ07an5nf6/XhOX0oQN2P1tz2VIA0s/nbay372FTqIxfYn9xtIuWCJFjmsaUmviNL8nJyb05Obk33z/qeDyeLOpXKG6sJ08RWj/JrSJ05Jppejyv9zyb5aXBtMhcc8lip9r89RvtS+Guxxs59Xx+cf5YLKnJydRo0sx2YTso+LlOfJBOg+pxeFMcwsyXKe12il1bpIXCTPP1o506U3fjhlzLBfoFUq14Nc7cPMOEHRPoX7e/zbCnNhmXAPvmwIrnzK6VgA4v0aXbOE5O7g1AmiGDIxdvsmjvWb7fetKiiic71aJNzTKEVwqhYgk7RwWXY2FatjTQVSLhiXXOvwvQf4VpCbRRWdvCENPQVMbPibCphjTT5kGFUO925zy+dxd+/v5Ub+SeTLZKiLwQ86Fu+VqtuXR8J1t/+ZnODw+36/may5dzom9f0i9ezHOtaM39a4xTj6Y/NnVtipZF80eg8UCYaGZNvOVT7adiBIzcQqC/Hw0rl6Bh5RK8c29D59s6tBgWPZbz+iuxUNy1oKLmU7KDj9i25Uq/akq3XP+Ac8aMmxdqvnL1WlVw6vn8QkpJ7PYt1Iho6lIwNHN0mZoJIUoLIdYKIY5m/s51vC2EmCWEuCSEOJTb/YJM35e0IW/CVc0Bdu+qpXY/W6h+PeNxbOO8v7FW9VtlPH5146uOdNE6AUHa6OjVbOFrL0RpI6RxJWDuQMt0PvYgJWyZaqojuwgNnKe166II/XHStH71YccP7VoDOdq+AwCFGjbELzjYRuncObxJM7voPMT1XShPcvH4MW5cvki9tq4585qj14hoFLBOSjk5M+f9KOCNXMr9AHyJ6xlDfI6q9TVtFn6m6UlqcpLdcYFrrVnN8V69kWlpGG7exL947jtMNUJq0LFqRzad3cTvJ3/nwXoP0qpSq1zLOkyx8powXD8NUyMs7x35FcZbWe/xCwAbiQstGDgPwno5308zEm4n8MpG05Z/j5pWYhqZkbTbZOtlj99fbsQduWY8Dgx2cYrsYWK3b8bPP4A6kS4YxWZDr8XqvkBW2urZaKFocyCl3AQ47hRUQMhatPYP1kZHf35nv7NqcK1axuMjLfMWlml3TjMeD/9juIVzp1soWV0TpLHXoOPrtsvbI0Ilqmne8+MS3CZC6RnpdJjfwXhu71T11MOaKUTZZ55xegdp2Sf7Aej/egunns8vMgwGYrZsILRJMwoVc18sbb2EqEJWfrLM3y6bkBbElNOdB2um8wGFtSFw9Oa/HHo+7B+T8+zpx5/Io6Tlh677ou4kpiY61JZd+PlB17c08RiXAKPOQJ8v7Hu2/Qsw6rRmKjAuAV465LD3fF5kyAya/WjaAbJnXQjg9PDHjcflnrNtdZ0bV+NM/9YVa3m3b9nx/XtIvBZPo653ubVej03NhBB/knsaorc80Z6UcgYwAyAyMrLApNcsHxrCpZOmyLxnow9RNdy+rV3h50eVLz4n7rnnubVlC2kXLhBY0fr6SdSwKOMibdt5bdk5eCdFAj2Y775QiLa43fwR22U9iCHDQNMfmxrP9z+8366RTdLevdzaqsXtrv37b063P3+CFsi/99PuSbftSaLW/UbRUqWp1cx5c4zc8NiISErZTUrZKJef5cBFIUQlgMzf1mNX/McZkDlUDyr2AAALxjmWkSOkuyl0xbHOXWyWN7cvaj23NQm39Y++6EnSMtIsRGjH4B0E+Nn+fpapqZwaojnPlho8mKAaNZxq//JpU7D/0MZlnaojv7h59Qon9u+lUedu+Ae4dwyj19RsBZAVd3QYYJ+L+X8Q4SeoGlYKv0BT4KnkRDszVWQSFm3y84oOyzvge6BfIHuHmkLUdpjfgb8vuRZbx1u5nnKd5j82N57/OeBPi0QDeRFjthtZcewYp/vwy3vaQvc9z7rHHseT7P9tJQARbp6WgX5CNBnoLoQ4CnTPPEcIUVkIYfQJEELMA7YD9YUQZ4UQ9hnSFDD6PK99YwcUagvAoolv51U8B0IIaq4waf3FD/IOmhDkH8TfD5vE5+FfH2bknyMdatPbWXhkIXcsMG0/bx+0nQpF7bPfMV8XMl+Hc5QTB0xrmTUauZYFw9PcTkriwNpfqdemPSXKuz/xsy5CJKW8KqW8U0pZN/N3fOb1c1LKXmblBkkpK0kpA6WUVaWU3+nRX70RfoJ299fBv5C2XXrpxL8Y0h0LPlaoXj1KPawFMYv//nsSVq7Ms7y/n7/FAvbWuK1EzI7I9wiMniBidgTjt483nh985CDFguzbAbq+eIlxXajWyhVOpQkCzShwzdfav+/AMW4yl/AgB9f9RmpyEi37eCbghvI18xGa3VUdIQR+gXUBWPnpZIfrqPjWmwRU0L71z732Okn79tl4QlvAfqj+Q8bzxnMas/uCfTGSvI2TCSctLKaLBBQhaliU3VvuN9au5fxb2l5L2eefI7huXaf7snqatitXskIRylTx7pTS6amp7Fu9jOqNGhsjh7obJUQ+xJB32xBYVPPX+nfPTodHRYCF79mpwUNIPXXKeuFM3m7zNn/0N1kbP/b7Y0TMjsCQYcjjKe9BSknE7AjuXWYK2/pjzx/ZOWSn3XXcXP8Xcc89D0BIr56Ue/ppG09Y5/rFJE4d0lxCBo71gdHQn7+SeC2eVn0f9FgbSoh8iJIVilCiXBH8AjUXjuUfOZ6mBiA8xhS95d+7e+QZcD+LSsUq5bCtafpjU1r86N0GeK9seIXGcyy3xfc/vJ+m5ZtaeSIniRs3cjZTeIp370aVTz5xqU8/v6OF0e32aLhF2BdvJDUlmR1Lf6F6o8bUaGz/v5mjePe/giIHQye0JbCotox2Yv8ebic5FwTfXIyOdexESkyMzWeEEEQNi+LX+381XkvNSCVidoTXrR8NXDWQiNkR/HHKNJKb0X0GUcOi7NqezyJx82bOPKkt1Bft1JGqX9hpgGmFdXNM/+7121Ryqa78YN/q5STfSKD9Q5619VJC5GMIIRj0Thv8g7VIjDOfy9tiOi/MxejEff24amdo06rFqxI1LIrPu3xucb3xnMZEzI4gPkUfr5zk9GSjKB6+atrNejzicaKGRdG2cluH6rv6ww+ceWIEAEXatqH6N9+41L8LxxOI2abFiXpiakcbpfUnKeE6u1cuoXZkG7eEg80LFQbEBylTuRh129xPzMY9pCTe4PLpU5Sr7pxBXXhMtNG26NJHH3P1u1nU2247XQ5Al+pdiBoWxaazm3hm3TPG650WmFIabX5oMyULlcztcbeQZkij26JuuYrfm63fZFDYIKfqjW3dhowEzZiz+N13U/Wzqa7187aBxVM0+6xeTzcmqJD3f/Q2/fwD6ampdBzyqMfbEt40nHYXkZGRcs+ePXp3w+N8/th00m5pYTxeWbDKRum8OffGGyQsX2E8D4v+x2EHzrM3z9JzSU+r99tXbs9HnT6ye6s8N1LSU5iwYwIr/l1htczsHrNpXqG51ft5kXHrFrEtIo3nlT+YTIm+fZ2qy5xpI7WwvXValOfuJ7w7+iLAuSPRzBvzGi37DqDj4Eftfk4IsVdKGWm7ZLbnlBD5LgZDBlMH9wEgpGJnnvjMtXhCt7Zt4/RjJpvRWqtXEVzbufR1a0+t5eUNLzv0TPGg4pQKLsW129e4meqY9fjzzZ7n8YjHXYqfnLh5C2eeME11627dQkAZ1w0NV007wKkobZfsmeldXa7P02RkGPh59Msk3Uzgf598TVChwrYfykQJkRn/FSECOHf0PPPe1j487QdNoc19DWw8kTeGGzc40qq1xTXztSRnSExNZMiaIcbcYO6gWGAxfuv/GyWC3eOtfrRTZ9IvmrKGODMizI2//zzN1kXHAG1dyBemZLtXLmHTT7O458VR1G/bwfYDZjgrRN7/r6LIk8p1KxHa9E5O/r2OrfNep2LtHwiNcN550j8kxGLdCDT/NFemKMWCirH8Pkt3wvSMdK4kX2HRkUV8czD3ReAAEcDAsIH0qd2HuqXqOrTbZS/xP//MxQkTjedlHh9O+VfdE6ny2N5LRhEaPK61T4jQ1bNn2LrgR+q0bEu9No6nQ3IWNSIqIHz8kJZiLqBwJwaPH0m56vanRbZG6qlT/Hu3ZYTCenv2WKQt8lVuHz/O8V69La7V27kD/xLuGWHFxV5j2adasLM+LzSlWnhpt9TrSTIMBuaNfY3rFy/w6EfTKFoyj4wpVnB2RKS27wsII776AYD05I0smLjBItiWswTVqEF4TDQlHzRZ1B6JjCQ6LJyMlBSX69eDjKQkosPCLUSo2oxvCI+JdpsIXTl70yhC3Yc38AkRAti1fBEXjh3hzsdGOiVCrqCEqIBQvExZ2g7QtqpvJ3zDvPE7uXLWPVEWK41/l7DDlvkLYps2IzosnNv//uuWNjxNSrQ23YxtbrIEL3HffYTHRFOso/tsei6fvsmCiZovXvsBdajX0v2e6p7gbPQhtv3yM/XbdaS+G4Pi24sSogJEuwdMKZFTE39hwcRdXDjunsBmwt+f8Jho6u2y9M863vseosPCuTrre7e0406klJx/912iw8I50e9+43W/kBDCog5SefL7bm3vbEy8Mb5Qq3tr0rRbdbfW7ymSEq6z+rMplKxYke5PPOuWRXpHUUJUwHjxZy3tkEyPw3D7HxZP2culUzdsPGU/WYvZYf8cxi/ElATw0pQpRIeFEx0WzvUl9qc+8gSJGzcSHRZOTHgDrs+bb7xeasgQwmOiqb9rJyIw0K1t/rP1HMunajGc2t1fh5a9a7q1fk8hMzJY8+XHpCQmcs+Lowh2IjGkO1CL1QWQSyeP8+Mbmqd4UMgw/PzLcNfwhtRt6ZnEfZc//4IrX32V672AcuWotWa11XRG7iBpzx5ODX3Y6v2aS5dQKDzvyJSusH3pMfb9fhqAHiMaUbu5y7kg8o1Nc39g9/JFdB/xLI3vtJ06yRbKjsiM/7oQAfz9+2rWzdLSDwWXfA4hAmnQoTJdhnrOZ0impnJ6+OMWeb6sUWHsGEo9+CDCgdjH6fHxXPt5LlemTwdD3iFISj/2GOVffcXpwGX2IKVkwcTdxo2BAW9EUqGmd6eKNufQhj/5/eupNOneizuHP+WWKZkSIjOUEGksmjSGUwe13Zvgki8Z/9BGftEZ/0DPz8pvbd/O6f/lkhLaQ9RZ9yeBVarkS1vJianMenWL8fzhSW0JKWO/BbLenP3nEAsnvk3VBo24f9Q4twXD9ykhEkKUBhYAocBJ4EEp5bVsZaqhZXitCGQAM6SUn9lTvxIiE1n2RQCFSplcLh56uxVlq+ZvZEBDQgLX5s3n8lTXHEhBG/GUe+Zp/Irmv03T0T0X+WOmybt/5Jed8Q/wneXWq2dPM/+dNygcUoLBEz9yW/568D0hmgLEm6WcLiWlfCNbmUpAJSnlPiFEcWAvcJ+U8p9cqrRACZEJKSWfDNQiE5YLrUVIxUeJi70OaJlks5I4KmwjMySLpuw15plr2r067ft7JnSqp0i4dIH5Y18nIyODQeM/pGRF98ZE8jWDRpspp6WU56WU+zKPbwLRQP6MuwsQQgie/3ExAJdPHic9aSndh2v+aIc3xTFt5HquX0zSs4s+wYUTCXz19F9GEXrwzZY+J0KJ1+JZNHEM6ampDHh7ottFyBX0GhFdl1KWNDu/JqW0asophAgFNgGNpJS57kULIUYAIwCqV6/e4pQdsZj/S6QkJjJt+EAA6rVuz91Pvca3L20y3q/esDT3PNtEFxsSbyY91cBPY7ZzKyEVgBLlCjNoXGuvD/GanaQbCSyc8BYJFy8w4O2JHgt05nVTMxspp2fbK0RCiGLARmCSlHKJPW2rqVnuJCVc5+sRWnbSWi1a0e/1scTuOM+fP5i86+9+ohF1WvjO9rMn2fvbSXYsM0UM8BWfsewkxl9l4cS3uXH5Eve9PoYaEZ6LPe11QpRno0LEAp2llOcz14I2SClzLFYIIQKBVcDvUkq7I5YrIbJO4rV4vhmpxR+uXL8Bg8ZPwWDIYMGEXVy7YJqi9X+jBRVrusf3ytc4GXXFmO4HIKxtRbo+Eu6To8Ubly+xcMJb3Eq4Tr83xlKtQYTth1zA14ToQ+Cq2WJ1aSnl69nKCLT1o3gp5YuO1K+EKG+Sb97gq8cHA1C0ZCmenD4HIQTXLtxi7jhLF44HRkdSvobv2Ma4wvH9l/n1G1NSSeEneHRye4qEBOnYK+eJP3eWRRPHkJqSxP2j3vV43GnwPSEqA/wCVAdOAw9IKeOFEJWBmVLKXkKIDsBmIApt+x7gTSnlmlwrNUMJkW1uJyXx5f9MXvUv/rzMaEtyNvYayzO9x7O49/kmVG/g3WmRnSV62znWz7HMYqKHeYM7OfNPFCs+moTw96f/m+OpUNO5SJuO4lNC5GmUENlHhsHAp4NNwc5GfvOjRfiHMzHxrMj0n8qiQYfK3PFQXQIC/fOtn54g5VYa62ZHc/LgFYvrg8a2pnRl34639M+m9fw+/XNKVqjI/aPHeSRXvTWUEJmhhMh+pJTMfO5xblzWwqQ+OPY9qjW0TEh46dQNFr6f89+z99ONCW3sfDTI/EZmSA5timPT/CMW14uUCOL+V1tQopzvWEbnRkaGge0L57JjyQKqNWxMn5ffpFCx/B3VKSEyQwmR42yYM5O9q5cB0Lrfg3QYmDOhXlqqgU3zjxhzc2Xh5ye4Y2A9GrSvhJ+XbWunpxk4uP4s25fmjJvUvEcNWt9b0+v67AxJNxJY88VHnDq4n4adu9H9iWfwD3BvhAF7UEJkhhIi5zi+bzdLP3gXACH8eOGnJVZ9kC6evMHvMw5xMz5npMbqDcvQomcNKtUuke87TTJDcvLQVXavOsHl0zkzgZQPDeHuJxr6lF+Yfl3xGgAADBZJREFULc4fi2XlJ5NJSrhG1/+NJOLOu3Xb4VNCZIYSIue5df0a0580hdQYMukTKtapl+czV+MS2br4GGf+sZ7htW7LCtRrVYEq9UsRGOSe9aXbSWmcOnSV2J0XOH04j7Yjy9PmvtqElC044gNaLKG9a5azee5sipUuQ5+XR1Ohlr7W3kqIzFBC5BpSSuaNfY3zR7SdpLqt2nHvy6Pt+pY1pGdwZNcF9v56ioTLyXa15+cvKFYqmCIhQfj5+yEEIARCQGpyOsk303IdeVmjfGgILe6uQc0mZRF+vmf7Yw8Jly7y29efcvafQ9SObM3dT71I4WKei/lkL0qIzFBC5B7Mp2oAj378FWWqOh7+NOVWGsf2XuLYnovEHbnuvg4KCG1UhnqtK1KzcVkC3DTS8maklBzasJYNs78FoMuwETTs3M1rjC2VEJmhhMh9pKem8sWjD5JhSAegRuNm9B/9rlsDjqWmpJOUkErSzVSkQSKlRErtQxdUOIDCxYIoEhJEYHDBF5q8iD93lnXffc3pQweo2qARPZ56iRLlPRN101mUEJmhhMj9xG7fwqqpk43n97z4hi7ZHv6LpN1OYefShexesZjA4GDaD3yYpt17eTT6pLMoITJDCZFnMKSn89PoF7ly+qTx2vDPvvWqcBIFCSklR3ZsZfPc70m4dJHwO7rQaehj+Z5zzBGUEJmhhMizXI07ww8vP2U8L1w8hMemzsh347mCzM6lv7Bl/hwAylStzp2PjcxhaOqNKCEyQwlR/hC9+S/WfPmx8bxM1eoMmvAhwUV820VCT/7+Yw3rvjNlRGn34BBa93sQPz/fWB9TQmSGEqL8ZesvP7FjsSl/WMkKlRg4fopXTyG8CSklu5YtNI6AQBtlPvrJ1xQJ8a1QLEqIzFBClP9IKdkybza7li+yuD70/am6G9l5K2mpt/nty084snOr8VqhosUY9tE0ipX2zUgHSojMUEKkH1JK9qxayqafZllcbztgEG37D/LKnZ785sKxI8wb+7rRJAKgbLUaPDD2PZ8bAWVHCZEZSoi8gxN/72XJ++9YXAsIDKL/W+OpGt5Ip17pQ8qtRP6c+RWx2zZZXG98Zw+6PjbSbXnF9EYJkRlKiLyL5MSbrPr0fU4fOmhxvVSlyvR+/vUCO3VLTU5i87zZ/P376hz3+o9+l9CmLXTolWdRQmSGEiLv5fShgyx5fyyG9HSL6/6BgXQb/jQNOnX1mR2i3Lhy5hQb5sw0Ztg1p2Wf/rR/aKgu4TnyCyVEZigh8g2O7trG6s8/xJCWluNe+Zq1adV3AHVbtcPP33uFKf7cWXavWMKhv/7I9X6Tu3pzx6BhBBcpks890wefEiI7U04XQstlFgwEAIuklJYLDlZQQuR7XD51gg1zZnL60AGrZeq1vYOGnbpSI6Jpvo8qpJTEx50lZttGDm9cx80rl3MtFxAcTOeHHyfizrt8emTnLL4mRPaknBZAUSllYmZaoS3AC1LKHbbqV0Lk20gpOfH3HvasXMqZwwdtli9TtTqV64VRPrQ25WvWolSlKhQqVtwhj3RDehrXL5znyplTXD51gnOx0cTFRlvsbOVGcJGitOw7gCbde7o1h7yv4qwQ6bVU3xfonHk8G9gAWAiR1BQyMfM0MPOn4M0jFTkQQlCrWUtqNWtpvHY76RYxWzcRvWUDcTGHLcpfPXuaq2dPe7xfxUqXIax9J8I7dKZ8aC2Pt/dfQi8hqiClPA9ajnshRK6pRYUQ/sBeoA4wTUq5M7dyioJPcJGiNOnekybde1pcz8gwcPnkCc4fO8Klk/9y6cS/XL9wnttJtxxuo1SlKpSpWp1yNUKpXC+cSnXrK3eVfMJjQmQj5bRdSCkNQFMhRElgqRCikZTykJX2RgAjAKpXdzx4l8I38fPzp0KtOgXWBOC/gseESErZzdo9IcRFIUQls5TTl2zUdV0IsQHoAeQqRFLKGcAM0NaInO64QqHId/Syt18BDMs8HgYsz15ACFEucySEEKIw0A2IyV5OoVD4PnoJ0WSguxDiKNA98xwhRGUhRFZK6UrAX0KIg8BuYK2UcpUuvVUoFB5Fl8VqKeVV4M5crp8DemUeHwSa5XPXFAqFDihXaIVCoTtKiBQKhe4oIVIoFLqjhEihUOiOEiKFQqE7SogUCoXuKCFSKBS6o4RIoVDojhIihUKhO0qIFAqF7ighUigUuqOESKFQ6I4SIoVCoTtKiBQKhe4oIVIoFLqjhEihUOiOEiKFQqE7SogUCoXuKCFSKBS6o4sQCSFKCyHWCiGOZv4ulUdZfyHEfiGECpyvUBRQ9BoRjQLWSSnrAusyz63xAhCdL71SKBS6oJcQ9UXLeU/m7/tyKySEqAr0BmbmU78UCoUO6JJOCKggpTwPkJnttbyVclOB14Hitio0TzkN3BZC5JoRtgBQFriidyc8iHo/36a+Mw95TIiEEH8CFXO59Zadz98DXJJS7hVCdLZV3jzltBBij5Qy0oHu+gwF+d1AvZ+vI4TY48xzHhMiKWU3a/eEEBeFEJUyR0OVgEu5FGsP9BFC9AIKASFCiJ+klEM91GWFQqETeq0RrUDLeU/m7+XZC0gpR0spq0opQ4GBwHolQgpFwUQvIZoMdBdCHAW6Z54jhKgshFjjhvpnuKEOb6Ugvxuo9/N1nHo/IaV0d0cUCoXCIZRltUKh0B0lRAqFQnd8XogKuruIPe8nhKgmhPhLCBEthDgshHhBj746ghCihxAiVghxTAiRw7JeaHyeef+gEKK5Hv10Fjveb0jmex0UQmwTQjTRo5/OYOvdzMq1FEIYhBADbFYqpfTpH2AKMCrzeBTwQR5lXwbmAqv07rc73w+oBDTPPC4OHAEa6N33PN7JH/gXqAUEAQey9xfoBfwKCKANsFPvfrv5/doBpTKPe/rK+9nzbmbl1gNrgAG26vX5EREF313E5vtJKc9LKfdlHt9E882rkm89dJxWwDEp5XEpZSowH+09zekLzJEaO4CSmTZnvoDN95NSbpNSXss83QFUzec+Oos9/3cAzwGLyd1GMAcFQYgs3EUAW+4iGfnVMTdh7/sBIIQIBZoBOz3eM+epApwxOz9LTuG0p4y34mjfh6ON/nwBm+8mhKgC9AOm21upXr5mDpHf7iL5javvZ1ZPMbRvoRellDfc0TcPIXK5lt2OxJ4y3ordfRdCdEETog4e7ZH7sOfdpgJvSCkNQuRWPCc+IUSygLuLuOH9EEIEoonQz1LKJR7qqrs4C1QzO68KnHOijLdiV9+FEI3Rlgp6Simv5lPfXMWed4sE5meKUFmglxAiXUq5zGqtei9+uWHx7EMsF3On2CjfGd9arLb5fmjfUnOAqXr31853CgCOAzUxLXg2zFamN5aL1bv07reb3686cAxop3d/3f1u2cr/gB2L1bq/mBv+YcqgBVc7mvm7dOb1ysCaXMr7mhDZfD+0Yb0EDgJ/Z/700rvvNt6rF9ru3r/AW5nXRgIjM48FMC3zfhQQqXef3fx+M4FrZv9fe/Tus7veLVtZu4RIuXgoFArdKQi7ZgqFwsdRQqRQKHRHCZFCodAdJUQKhUJ3lBApFArdUUKkyFeEENuceCZACHFFCPF+tusnhRBlzc47+1JkBYUJJUSKfEVK2c6Jx+4CYoEHhb0+AwqfQgmRwi1kxp45KIQoJIQomhkXqVEu5RIzf3cWQmwQQiwSQsQIIX7OQ2QGAZ8Bp9GsrBUFDJ/wNVN4P1LK3UKIFcBEoDDwk5TSVpLLZkBDNF+lrWg+gVvMCwghCgN3Ak8CJdFEabt7e6/QGzUiUriT8WhZWSLRArrZYpeU8qyUMgPNzSE0lzL3AH9JKZPQnHr7CSH8M+/l5hagXAV8ECVECndSGiiGFiWykB3lb5sdG8h9hD4I6CaEOAnsRfO965J57ypgHjq3NAU7nXOBRQmRwp3MAMYAPwMfuFqZECIEzaG3upQyVGrJNp9BEyeADcDDmWX9gaHAX662q8h/lBAp3IIQ4hEgXUo5Fy1hZkshRFcXq70fLcOv+chpOVpsqWBgAlBHCHEA2I8WVuMnF9tU6IDyvlcoFLqjRkQKhUJ3lBApFArdUUKkUCh0RwmRQqHQHSVECoVCd5QQKRQK3VFCpFAodOf/Vcq1aLt4iN0AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# use the downsampled versions!!\n", "fig, ax = plt.subplots(1,1,figsize=(4,4))\n", "\n", "# want to loop over all bodies in this simulation and plot their TOTAL orbits -- all times\n", "# x/y positions only\n", "# r_h[# bodies, x/y/z 3D positions, all time steps]\n", "for i in range(r.shape[0]): # looping over all bodies\n", " # now we want to plot (for given body) the x/y positions\n", " # over ALL times\n", " # r_h[ith planet, x position, all times], r_h[ith planet, y position, all times]\n", " ax.plot(r[i,0,:], r[i,1,:])\n", " \n", "ax.set_xlabel('x in AU')\n", "ax.set_ylabel('y in AU')\n", "ax.set_xlim(-0.4, 0.4)\n", "ax.set_ylim(-0.4, 0.4)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After that bit of downsampling tangent, let's actually make that animation!" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUQAAAEzCAYAAABJzXq/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAANbUlEQVR4nO3cf6jd9X3H8edrSQOt7ao0aenyg2UjrWZDR711UvbDTrYm7o9Q8A+1TCaFIGjpn8r+aAf+s/4xKMUfIUiQ/tP8U+nSkSpjo3Vg0+YGNBpFuYvM3KZgrKUDC5Poe3+c9+bp3Y33m5tzzvWG5wMu3O/3fO657w9Xnn7PufebVBWSJPittR5Akt4vDKIkNYMoSc0gSlIziJLUDKIktRWDmORQkteSPH+Bx5PkW0kWkpxM8pnJjylJ0zfkCvExYM97PL4X2NUf+4FHLn0sSZq9FYNYVU8Bb7zHkn3At2vkGHBlkk9OakBJmpVJvIe4FTgzdrzY5yRpXdk4gefIMueWvR8wyX5GL6u54oorrr/66qsn8O0l6V0nTpx4vaq2rOZrJxHERWD72PE24OxyC6vqIHAQYG5urubn5yfw7SXpXUn+c7VfO4mXzEeAO/u3zTcCv6qqn0/geSVppla8QkzyHeAmYHOSReDrwAcAquoAcBS4BVgAfg3cNa1hJWmaVgxiVd2+wuMF3DOxiSRpjXiniiQ1gyhJzSBKUjOIktQMoiQ1gyhJzSBKUjOIktQMoiQ1gyhJzSBKUjOIktQMoiQ1gyhJzSBKUjOIktQMoiQ1gyhJzSBKUjOIktQMoiQ1gyhJzSBKUjOIktQMoiQ1gyhJzSBKUjOIktQMoiQ1gyhJzSBKUjOIktQMoiQ1gyhJzSBKUjOIktQMoiQ1gyhJzSBKUjOIktQMoiQ1gyhJzSBKUjOIktQMoiQ1gyhJbVAQk+xJ8lKShST3L/P4R5N8P8mzSU4luWvyo0rSdK0YxCQbgIeAvcBu4PYku5csuwd4oaquA24C/jHJpgnPKklTNeQK8QZgoapOV9VbwGFg35I1BXwkSYAPA28A5yc6qSRN2ZAgbgXOjB0v9rlxDwLXAGeB54CvVtU7S58oyf4k80nmz507t8qRJWk6hgQxy5yrJcdfAJ4Bfgf4I+DBJL/9/76o6mBVzVXV3JYtWy56WEmapiFBXAS2jx1vY3QlOO4u4PEaWQBeAa6ezIiSNBtDgngc2JVkZ/+i5DbgyJI1rwI3AyT5BPBp4PQkB5Wkadu40oKqOp/kXuBJYANwqKpOJbm7Hz8APAA8luQ5Ri+x76uq16c4tyRN3IpBBKiqo8DRJecOjH1+FviryY4mSbPlnSqS1AyiJDWDKEnNIEpSM4iS1AyiJDWDKEnNIEpSM4iS1AyiJDWDKEnNIEpSM4iS1AyiJDWDKEnNIEpSM4iS1AyiJDWDKEnNIEpSM4iS1AyiJDWDKEnNIEpSM4iS1AyiJDWDKEnNIEpSM4iS1AyiJDWDKEnNIEpSM4iS1AyiJDWDKEnNIEpSM4iS1AyiJDWDKEnNIEpSM4iS1AyiJDWDKEnNIEpSM4iS1AYFMcmeJC8lWUhy/wXW3JTkmSSnkvxosmNK0vRtXGlBkg3AQ8BfAovA8SRHquqFsTVXAg8De6rq1SQfn9bAkjQtQ64QbwAWqup0Vb0FHAb2LVlzB/B4Vb0KUFWvTXZMSZq+IUHcCpwZO17sc+M+BVyV5IdJTiS5c1IDStKsrPiSGcgy52qZ57keuBn4IPDjJMeq6uXfeKJkP7AfYMeOHRc/rSRN0ZArxEVg+9jxNuDsMmueqKo3q+p14CnguqVPVFUHq2ququa2bNmy2pklaSqGBPE4sCvJziSbgNuAI0vW/BPwp0k2JvkQ8MfAi5MdVZKma8WXzFV1Psm9wJPABuBQVZ1Kcnc/fqCqXkzyBHASeAd4tKqen+bgkjRpqVr6duBszM3N1fz8/Jp8b0mXryQnqmpuNV/rnSqS1AyiJDWDKEnNIEpSM4iS1AyiJDWDKEnNIEpSM4iS1AyiJDWDKEnNIEpSM4iS1AyiJDWDKEnNIEpSM4iS1AyiJDWDKEnNIEpSM4iS1AyiJDWDKEnNIEpSM4iS1AyiJDWDKEnNIEpSM4iS1AyiJDWDKEnNIEpSM4iS1AyiJDWDKEnNIEpSM4iS1AyiJDWDKEnNIEpSM4iS1AyiJDWDKEnNIEpSM4iS1AYFMcmeJC8lWUhy/3us+2ySt5PcOrkRJWk2Vgxikg3AQ8BeYDdwe5LdF1j3DeDJSQ8pSbMw5ArxBmChqk5X1VvAYWDfMuu+AnwXeG2C80nSzAwJ4lbgzNjxYp/7P0m2Al8EDkxuNEmarSFBzDLnasnxN4H7qurt93yiZH+S+STz586dGzqjJM3ExgFrFoHtY8fbgLNL1swBh5MAbAZuSXK+qr43vqiqDgIHAebm5pZGVZLW1JAgHgd2JdkJ/Ay4DbhjfEFV7fzfz5M8Bvzz0hhK0vvdikGsqvNJ7mX02+MNwKGqOpXk7n7c9w0lXRaGXCFSVUeBo0vOLRvCqvrbSx9LkmbPO1UkqRlESWoGUZKaQZSkZhAlqRlESWoGUZKaQZSkZhAlqRlESWoGUZKaQZSkZhAlqRlESWoGUZKaQZSkZhAlqRlESWoGUZKaQZSkZhAlqRlESWoGUZKaQZSkZhAlqRlESWoGUZKaQZSkZhAlqRlESWoGUZKaQZSkZhAlqRlESWoGUZKaQZSkZhAlqRlESWoGUZKaQZSkZhAlqRlESWoGUZKaQZSkZhAlqRlESWqDgphkT5KXkiwkuX+Zx7+U5GR/PJ3kusmPKknTtWIQk2wAHgL2AruB25PsXrLsFeDPq+pa4AHg4KQHlaRpG3KFeAOwUFWnq+ot4DCwb3xBVT1dVb/sw2PAtsmOKUnTNySIW4EzY8eLfe5Cvgz8YLkHkuxPMp9k/ty5c8OnlKQZGBLELHOull2YfJ5REO9b7vGqOlhVc1U1t2XLluFTStIMbBywZhHYPna8DTi7dFGSa4FHgb1V9YvJjCdJszPkCvE4sCvJziSbgNuAI+MLkuwAHgf+pqpenvyYkjR9K14hVtX5JPcCTwIbgENVdSrJ3f34AeBrwMeAh5MAnK+quemNLUmTl6pl3w6curm5uZqfn1+T7y3p8pXkxGovyLxTRZKaQZSkZhAlqRlESWoGUZKaQZSkZhAlqRlESWoGUZKaQZSkZhAlqRlESWoGUZKaQZSkZhAlqRlESWoGUZKaQZSkZhAlqRlESWoGUZKaQZSkZhAlqRlESWoGUZKaQZSkZhAlqRlESWoGUZKaQZSkZhAlqRlESWoGUZKaQZSkZhAlqRlESWoGUZKaQZSkZhAlqRlESWoGUZKaQZSkZhAlqRlESWoGUZLaoCAm2ZPkpSQLSe5f5vEk+VY/fjLJZyY/qiRN14pBTLIBeAjYC+wGbk+ye8myvcCu/tgPPDLhOSVp6oZcId4ALFTV6ap6CzgM7FuyZh/w7Ro5BlyZ5JMTnlWSpmpIELcCZ8aOF/vcxa6RpPe1jQPWZJlztYo1JNnP6CU1wH8neX7A91+vNgOvr/UQU+T+1q/LeW8An17tFw4J4iKwfex4G3B2FWuoqoPAQYAk81U1d1HTriPub327nPd3Oe8NRvtb7dcOecl8HNiVZGeSTcBtwJEla44Ad/Zvm28EflVVP1/tUJK0Fla8Qqyq80nuBZ4ENgCHqupUkrv78QPAUeAWYAH4NXDX9EaWpOkY8pKZqjrKKHrj5w6MfV7APRf5vQ9e5Pr1xv2tb5fz/i7nvcEl7C+jlkmSvHVPktrUg3i53/Y3YH9f6n2dTPJ0kuvWYs7VWGlvY+s+m+TtJLfOcr5LNWR/SW5K8kySU0l+NOsZL8WA/zY/muT7SZ7t/a2b9/6THEry2oX+dG/VXamqqX0w+iXMfwC/B2wCngV2L1lzC/ADRn/LeCPwk2nOtAb7+xxwVX++d73sb8jextb9G6P3mG9d67kn/LO7EngB2NHHH1/ruSe8v78DvtGfbwHeADat9ewD9/dnwGeA5y/w+Kq6Mu0rxMv9tr8V91dVT1fVL/vwGKO/0VwPhvzsAL4CfBd4bZbDTcCQ/d0BPF5VrwJU1Xra45D9FfCRJAE+zCiI52c75upU1VOM5r2QVXVl2kG83G/7u9jZv8zo/1rrwYp7S7IV+CJwgPVnyM/uU8BVSX6Y5ESSO2c23aUbsr8HgWsY3UTxHPDVqnpnNuNN3aq6MujPbi7BxG77e58aPHuSzzMK4p9MdaLJGbK3bwL3VdXbo4uMdWXI/jYC1wM3Ax8EfpzkWFW9PO3hJmDI/r4APAP8BfD7wL8k+feq+q9pDzcDq+rKtIM4sdv+3qcGzZ7kWuBRYG9V/WJGs12qIXubAw53DDcDtyQ5X1Xfm82Il2Tof5uvV9WbwJtJngKuA9ZDEIfs7y7gH2r0pttCkleAq4GfzmbEqVpdV6b8xudG4DSwk3ff2P2DJWv+mt988/Ona/2G7YT3t4PRHTyfW+t5J723JesfY339UmXIz+4a4F977YeA54E/XOvZJ7i/R4C/788/AfwM2LzWs1/EHn+XC/9SZVVdmeoVYl3mt/0N3N/XgI8BD/eV1PlaBzfWD9zbujVkf1X1YpIngJPAO8CjVbUu/oWmgT+/B4DHkjzHKBz3VdW6+FdwknwHuAnYnGQR+DrwAbi0rniniiQ171SRpGYQJakZRElqBlGSmkGUpGYQJakZRElqBlGS2v8AhbvZaFRTi14AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# create a figure objection (but don't plot anything on it)\n", "fig, ax = plt.subplots(1,1, figsize=(5,5))" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [], "source": [ "# 2nd step is to use plot_animations to setup the animation stuff we need\n", "init, animate, nFrames = plot_animations(fig, ax, t, r) # t is need here because Jill wrote the library funny!" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": [ "# 3rd step is to use matplotlib to put animation files together\n", "anim = animation.FuncAnimation(fig, animate, init_func=init, \n", " frames=nFrames, interval=20, blit=True)" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [], "source": [ "anim.save?" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [], "source": [ "anim.save('anim_trial.mp4')" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [], "source": [ "# if the above didn't work, you can try this:\n", "Writer = animation.writers['ffmpeg']\n", "writer = Writer(fps=15, metadata=dict(artist='Me'), bitrate=1800)" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [], "source": [ "anim.save('anim_trial2.mp4', writer=writer)" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [], "source": [ "from IPython.display import Video" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Video(\"anim_trial.mp4\", width=400)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Stuff in 3D!" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [], "source": [ "from hermite_library import read_hermite_solution_from_file\n", "planet_file = 'data/Kepler-11-savedSim.txt'\n", "#planet_file = 'Kepler-11-savedSim.txt' # if all .txt files are in the same directory as my notebook file\n", "\n", "t_h, E_h, r_h, v_h = read_hermite_solution_from_file(planet_file)" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [], "source": [ "# import numpy as np\n", "# import ipyvolume as ipv\n", "# V = np.zeros((128,128,128)) # our 3d array\n", "# # outer box\n", "# V[30:-30,30:-30,30:-30] = 0.75\n", "# V[35:-35,35:-35,35:-35] = 0.0\n", "# # inner box\n", "# V[50:-50,50:-50,50:-50] = 0.25\n", "# V[55:-55,55:-55,55:-55] = 0.0\n", "# ipv.quickvolshow(V, level=[0.25, 0.75], opacity=0.03, level_width=0.1, data_min=0, data_max=1)" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [], "source": [ "#!conda install -c conda-forge ipyvolume --yes\n", "import ipyvolume\n", "import numpy as np" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "1c709f7fafbd4e07a13dae55cd59ce52", "version_major": 2, "version_minor": 0 }, "text/plain": [ "VBox(children=(Figure(camera=PerspectiveCamera(fov=45.0, position=(0.0, 0.0, 2.0), quaternion=(0.0, 0.0, 0.0, …" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "x, y, z = np.random.random((3,10))\n", "ipyvolume.quickscatter(x,y,z, size=1, marker=\"sphere\")" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(10,)" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x.shape" ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(7, 8800)" ] }, "execution_count": 40, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r_h[:,0,:].shape # [all the planets, x-values, all times]" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([-8.69460729e-02, -8.69427997e-02, -8.69329803e-02, ...,\n", " -2.30013480e-05, -2.30219366e-05, -2.30424531e-05])" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r_h[:,0,:].ravel() # take our x array for all of the planets and making it into a 1D array (ipyvolume)" ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(61600,)" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r_h[:,0,:].ravel().shape" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [], "source": [ "# if we wanna downsample\n", "stepSize = 50\n", "r = r_h[:,:,0:-1:stepSize] # technically this goes up until the 2nd to last timestep \n", "\n", "# do this for all of the positions x/y/z\n", "x = r[:,0,:].ravel()\n", "y = r[:,1,:].ravel()\n", "z = r[:,2,:].ravel()" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "955db879a5ab49c3a04ef52ec39603f3", "version_major": 2, "version_minor": 0 }, "text/plain": [ "VBox(children=(Figure(camera=PerspectiveCamera(fov=45.0, position=(0.0, 0.0, 2.0), quaternion=(0.0, 0.0, 0.0, …" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ipyvolume.quickscatter(x,y,z, size=1, marker=\"sphere\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Going to make a more complex figure and loop and plot each orbit!" ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(176,)" ] }, "execution_count": 45, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r[0,0,:].shape" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "13c37c09ace64a539c022a2a1d53f5ab", "version_major": 2, "version_minor": 0 }, "text/plain": [ "VBox(children=(Figure(camera=PerspectiveCamera(fov=45.0, position=(0.0, 0.0, 2.0), quaternion=(0.0, 0.0, 0.0, …" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ipyvolume.figure()\n", "colors = ['red', 'blue', 'green', 'orange', 'gray', 'yellow', 'magenta']\n", "#colors = [(1,0,0),(0,0,1), (0,1,0), (1,0.5,0), (0.5,0.5,0.5), (0.5,1,0), (1,0,1)]\n", "for i in range(r.shape[0]): # looping over planets\n", " ipyvolume.scatter(r[i,0,:], r[i,1,:], r[i,2,:], \n", " color=colors[i], marker='sphere')\n", " \n", "ipyvolume.show()" ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [], "source": [ "from flip_colors import flip_colors" ] }, { "cell_type": "code", "execution_count": 48, "metadata": {}, "outputs": [], "source": [ "color = [(1,0,0), (0,0,1), (0,1,0), (0,1,1), (1,1,0), (1,0,1), (0.5, 0.5, 0.5)]" ] }, { "cell_type": "code", "execution_count": 49, "metadata": {}, "outputs": [], "source": [ "colors = flip_colors(color,r)" ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "((176, 7, 3), (176, 3, 7))" ] }, "execution_count": 50, "metadata": {}, "output_type": "execute_result" } ], "source": [ "colors.shape, r.T.shape # so this is a way to do more data formatting specific to animations" ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "f8623a4aafe04cab90f6768840334f48", "version_major": 2, "version_minor": 0 }, "text/plain": [ "VBox(children=(Figure(animation=200.0, camera=PerspectiveCamera(fov=45.0, position=(0.0, 0.0, 2.0), quaternion…" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# all of this weird data manipulation is so we can do animations in 3D\n", "ipyvolume.figure()\n", "\n", "s = ipyvolume.scatter(r[:,0,:].T+0.5, r[:,1,:].T+0.5, r[:,2,:].T+0.5, marker='sphere', color=colors)\n", "\n", "ani = ipyvolume.animation_control(s, interval=200)\n", "\n", "ipyvolume.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "export to HTML" ] }, { "cell_type": "code", "execution_count": 52, "metadata": {}, "outputs": [], "source": [ "ipyvolume.figure()\n", "\n", "s = ipyvolume.scatter(r[:,0,:].T+0.5, r[:,1,:].T+0.5, r[:,2,:].T+0.5, marker='sphere', color=colors)\n", "\n", "ani = ipyvolume.animation_control(s, interval=200)" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [], "source": [ "import ipywidgets" ] }, { "cell_type": "code", "execution_count": 54, "metadata": {}, "outputs": [], "source": [ "myVBox = ipywidgets.VBox([ipyvolume.gcc()]) # grabbing the last plotted thing" ] }, { "cell_type": "code", "execution_count": 55, "metadata": {}, "outputs": [], "source": [ "ipyvolume.embed.layout = myVBox.children[0].layout\n", "ipyvolume.embed.layout.min_width='400px'" ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [], "source": [ "ipyvolume.embed.embed_html(\"myNewPage.html\", myVBox, offline=False, devmode=False)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.7.6" } }, "nbformat": 4, "nbformat_minor": 4 }