From 82b597c01f766f6e38cd40af15395e19c00b5789 Mon Sep 17 00:00:00 2001 From: David Minton Date: Sat, 22 May 2021 21:21:21 -0400 Subject: [PATCH] Improved energy and momentum bookkeeping --- .../collision_visualization.ipynb | 170 ++++++++---------- examples/symba_energy_momentum/escape.in | 18 ++ .../symba_energy_momentum/param.escape.in | 33 ++++ .../symba_energy_momentum/param.merger.in | 3 +- examples/symba_energy_momentum/sun.in | 4 +- src/io/io_conservation_report.f90 | 7 +- src/main/swiftest_symba.f90 | 6 +- src/modules/module_interfaces.f90 | 4 +- src/modules/swiftest_data_structures.f90 | 1 + src/symba/symba_discard_conserve_mtm.f90 | 61 ++++--- 10 files changed, 176 insertions(+), 131 deletions(-) create mode 100644 examples/symba_energy_momentum/escape.in create mode 100644 examples/symba_energy_momentum/param.escape.in diff --git a/examples/symba_energy_momentum/collision_visualization.ipynb b/examples/symba_energy_momentum/collision_visualization.ipynb index e3e542221..e5ec77c73 100644 --- a/examples/symba_energy_momentum/collision_visualization.ipynb +++ b/examples/symba_energy_momentum/collision_visualization.ipynb @@ -25,17 +25,17 @@ "name": "stdout", "output_type": "stream", "text": [ - "Reading Swiftest file param.sun.in\n", - "Reading in time 3.810e-03\n", + "Reading Swiftest file param.escape.in\n", + "Reading in time 1.000e-04\n", "Creating Dataset\n", "\n", "Adding particle info Dataset\n", - "Successfully converted 382 output frames.\n" + "Successfully converted 11 output frames.\n" ] } ], "source": [ - "config_file_name = 'param.sun.in'\n", + "config_file_name = 'param.escape.in'\n", "config = swio.read_swiftest_config(config_file_name)\n", "ds = swio.swiftest2xr(config)" ] @@ -57,7 +57,7 @@ { "data": { "text/plain": [ - "[]" + "[]" ] }, "execution_count": 4, @@ -66,7 +66,7 @@ }, { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlsElEQVR4nO3dd3yV9d3/8dcng4QdICFAEiBsAjLDHoKroFK0jjqxtUqR4qh6e2u923u0/d1dtkrrqKtWRRFnQUVUlKGywgp7rxBGIBI2Wd/fHzntHeNhhOTKdc7J+/l4nAc55/pe53p7qbxzbXPOISIiUlGU3wFERCQ0qSBERCQoFYSIiASlghARkaBUECIiEpQKQkREglJBiJyFma0xsxGnmTbCzHJqNpFIzVBBiJyFc66bc25OTS3PzDqZ2T/MLM/M8s1slpl1rqnli/yTCkIk9CQA04HOQDKwGPiHn4GkdlJBiJyFmW03s0sCP9c1s5fM7GszWwv0q+7lOecWO+decM7lO+eKgD8Bnc2sWXUvS+RMYvwOIBJm/hNoH3jVB2aeabCZZQOtTzP5NefcxHNY5nBgr3PuYGWCilSVCkKkcq4HJjrn8oF8M5sM/OJ0g51zPaqyMDNLBZ4E7q/K94icD+1iEqmcVsCucu93eLUgM0sCPgaecs697tVyRE5HBSFSOXuAtHLvT7f7CPjXKbJHT/N65gzzNaGsHKY7535dPdFFKke7mEQqZxrwiJktouwYxN1nGuyc61bZBZhZI2AW8KVz7uHzSilSDbQFIVI5/03ZbqVtlP2G/4oHy7iasrOjflhhi+OMWysi1c30wCAREQlGWxAiIhKUCkJERIJSQYiISFAqCBERCSqiTnNNTEx0bdu29TuGiEjYWLp06QHnXFKwaRFVEG3btiUrK8vvGCIiYcPMTns3AO1iEhGRoFQQIiISlApCRESCUkGIiEhQKggREQnK04Iws1FmtsHMNpvZt+5KaWZdzGyBmZ0yswcrM6+IiHjLs4Iws2jKnoQ1GsgAbjSzjArD8oF7gD+cx7wiIuIhL7cg+gObnXNbnXOFwFRgbPkBzrn9zrklQFFl561Ok2dvYvrKXA4ePeXVIkREwo6XF8ql8M1HM+YAA6p7XjMbD4wHaN268rfLP1lUwktfbSf/WCEA3VMaMbRDEsM6JtK3TRPiY6Mr/Z0iIpHAy4KwIJ+d68Mnznle59yzwLMAmZmZlX64RXxsNEsevYRVuwv4YlMe8zYd4Pn5W3lm7hbiY6Pon96MYR0SGdYpkc7JDTELFk1EJPJ4WRA5fPPZvalAbg3MW2nRUUavtAR6pSUw6aKOHDtVzKJtB5m38QBfbD7Arz9cBx9CUsM4hnZIZFjHRIZ2SKR5o3ivIomI+M7LglgCdDSzdGA3cANwUw3MW2X142K4qEsyF3VJBmBPwQnmbzrA/E0HmLsxj3eX7wagS4uGXNg5iZGdm9O3TRNio3XWsIhEDk8fOWpmlwOPA9HAi865X5vZBADn3DNm1gLIAhoBpcBRIMM5dzjYvGdbXmZmpvP6Zn2lpY61ew4zf9MB5m3MY8n2fIpLHQ3jYhjWKZERnZozonOSti5EJCyY2VLnXGbQaZH0TOqaKIiKjpws4svNB5mzYT+fb9jPvsNlZ0J1a9WIkZ2bM7JLEr3SmhAdpWMXIhJ6VBA1xDnH+r1H+HzDfuasz2Ppzq8pKXUk1ItlWMckRnZO4sJOSTRrEOdbRhGR8lQQPik4XsT8zXnM2VD2OnD0FGbQt3UTLslI5pKuyXRo3sDvmCJSi6kgQkBpqWNN7mE+XbeP2ev3sXr3YQDaJdb/V1n0baNdUSJSs1QQISj30Almr9vHx2v3sXDrQYpKHE3qxXJRl2QuzWjOsI5J1I+LqAf+iUgIUkGEuCMni5i7MY9P1+7js/X7OXyymDoxUQxp34xLMpK5tGuyzooSEU+oIMJIUUkpS7bn8+na/Xyybi+78k/867jFqO4tGNW9BalN6vkdU0QihAoiTDnn2LDvCLNW72Pm6j2s33sEgB6pjRnVvQWju7ckPbG+zylFJJypICLE9gPHmLl6Lx+t3sPKnAKg7GruUd1bcPkFLenYvIHuFSUilaKCiEC7D53go0BZZO34GuegXVJ9Rge2LLq1aqSyEJGzUkFEuP2HTzJr7T4+Wr2HhVvzKSl1tG5ajzE9WzKmZyvdhVZETksFUYvkHyvkk7V7eT97D19tOUhJqaNj8waM6dmKK3u0pF2SLswTkf+jgqilDhw9xczVe5mxMpcl2/NxruyBSGN6tOKKHi11NpSIqCCk7JblH2TvYUb2HlbuOgRA3zZNGNOjJZf3aEnzhrrOQqQ2UkHIN+w8eJwZ2bnMWJnL+r1HMIOB6c24qncrRl/QkkbxsX5HFJEaooKQ09q07wgzsvcwY2Uu2w4cIy4miksykvle7xSGd0rSQ5BEIpwKQs7KOcfKnALeXZbDjOw95B8rpGn9Oozp0ZKr+6TSM7WxzoQSiUAqCKmUopJS5m4oe7TqJ+v2UVhcSruk+lzdK4WreqeQ1lQHt0UihQpCzlvBiSJmrtrDO8t3s3hbPgD92zbl6j4pXNFDxytEwp0KQqrFrvzjTF+ZyzvLctiSV3a84vILWnJd31QGtmtGlJ5lIRJ2VBBSrZxzZOcU8ObSXfxjRS5HThaT2qQu1/VN49rMVFIS6vodUUTOkQpCPHOyqIRZa/YyLWsXX24+iBkM7ZDIdZlpXJaRTHxstN8RReQMVBBSI3blH+etpTm8tTSH3YdO0Cg+hqt6p3Bd3zS6p+jmgSKhSAUhNaq01PHVloO8uXQXM1fvpbC4lC4tGnJj/9Zc1TuFxnV1YFskVKggxDcFx4uYnp3LtCW7WLW7gPjYKMb0aMVNA1rTKy1BWxUiPlNBSEhYlVPAa4t38I8VuRwvLKFry0bcNKA1V/VqRUOdLiviCxWEhJQjJ4v4x4pcpizaybo9h6lXJ5qxvVpxU/82XJDa2O94IrWKCkJCknOOFbsO8dqinczIzuVkUSkXpDTm5gGtGdOzFfXjYvyOKBLxVBAS8gpOFPHe8t1MWbSDjfuO0iAuhmv6pHDroLZ0aK6HHIl4RQUhYcM5x9IdXzNl0U4+yN5DYUkpwzomctugtozs0pxoXa0tUq1UEBKWDhw9xdTFO3l14U72Hj5JWtO63DqwDddnppFQr47f8UQiggpCwlpRSSmfrN3HS19tZ/G2fOJjo7iqVwrjBrUlo1Ujv+OJhDUVhESMdXsO8/KC7by7fDcni0rp37Yp4wa34TvdWujhRiLnQQUhEefQ8ULezMrh5YXb2ZV/gpaN4/nB4Lbc0L+1rtQWqYQzFYSnv3KZ2Sgz22Bmm83s4SDTzcwmB6Znm1mfctN+amZrzGy1mb1uZvFeZpXwklCvDncOb8ecB0fy/LhM0hPr878z1zPof2fzX9PXsPPgcb8jioQ9z7YgzCwa2AhcCuQAS4AbnXNry425HLgbuBwYADzhnBtgZinAF0CGc+6EmU0DPnTOvXSmZWoLonZbk1vAC19sY8bKXEpKHZdltOCOYen0bdNEt/QQOY0zbUF4eSVSf2Czc25rIMRUYCywttyYscDLrqylFppZgpm1LJetrpkVAfWAXA+zSgTo1qoxf7y+Fw99pwsvL9jOlEU7+WjNXnqlJXDHsHRGdWtBjI5TiJwzL/9vSQF2lXufE/jsrGOcc7uBPwA7gT1AgXPu42ALMbPxZpZlZll5eXnVFl7CV4vG8Tw0qgsLHrmIX47txqHjhUx6bTkX/n4Oz8/fyuGTRX5HFAkLXhZEsG36ivuzgo4xsyaUbV2kA62A+mZ2S7CFOOeedc5lOucyk5KSqhRYIku9OjHcOqgtnz0wgufGZZLapC6/+mAdQ/73M34zcz37j5z0O6JISPNyF1MOkFbufSrf3k10ujGXANucc3kAZvYOMBh41bO0ErGiooxLM5K5NCOZVTkF/HXeFp6dt4UXv9zGtX1TGT+sHW0T6/sdUyTkeLkFsQToaGbpZlYHuAGYXmHMdGBc4GymgZTtStpD2a6lgWZWz8qOLl4MrPMwq9QSF6Q25i839eGzB0Zwbd9U3srK4aLH5vCT15axeneB3/FEQoqn10EEzlJ6HIgGXnTO/drMJgA4554J/OX/F2AUcBz4oXMuKzDvfwPfB4qB5cAdzrlTZ1qezmKSytp/+CQvfrmdKQt3cORUMcM6JnLXiPYMatdMZz5JraAL5UTO4vDJIl5duIMXv9jOgaOn6JmWwF0XtuOyjBZE6QaBEsFUECLn6GRRCW8vy+Gvc7eyM/847ZLqM2lkB77bs5VOkZWIpIIQqaTiklJmrt7Lk59vZv3eI7RpVo+fjOzA1b1TdM8niSgqCJHzVFrq+GTdPv782SZW7z5MapO6TBzRgWv7plInRkUh4U8FIVJFzjk+37CfJ2ZvZuWuQ7RqHM9dI9pzXWYa8bHRfscTOW8qCJFq4pxj3qYDTJ69iaU7via5URwTLmzPjf1bqygkLKkgRKqZc44FWw7yxOxNLNqWT2KDOO4a0Z6bB6goJLyoIEQ8tHDrQSbP3sRXWw7SolE8ky7qwPWZaTpGIWFBBSFSA77acoDHPt7I0h1fk9qkLvde3JGre6fo9FgJab49MEikNhncPpG3JgzipR/2o0m9OvzbW9lc9vg8ZqzMpbQ0cn4Rk9pDBSFSjcyMEZ2bM33SEP56a19io6K4+/XlXD55Ph+v2UskbbFL5FNBiHjAzPhOtxbMvHcYk2/sTWFxKeNfWcrYJ79k7sY8FYWEBRWEiIeioozv9mzFxz8dzu+u7cHBo4Xc9uJibn5+Edk5h/yOJ3JGOkgtUoMKi0t5bdEOJn+2mfxjhVzZoyUPXtZZz6MQ3+gsJpEQc+RkEc/N28pz87dRVFLKTQNac8/FHUlsEOd3NKllVBAiIWr/kZNMnr2J1xfvIj4mijuHt+OOYe1oEOflwx5F/o8KQiTEbc07yh8+3sCHq/aS2KAO91zckRv7t9adY8Vzug5CJMS1S2rAUzf35d2Jg2mX1IBf/GMN3/nTPD5du09nPIlvVBAiIaR36ya8MX4gL9yWCQZ3vJzFLS8sYm3uYb+jSS2kghAJMWbGxV2TmXXfcP77u91Yk3uYK/48n4ffzmb/kZN+x5NaRAUhEqJio6O4bXBb5j44ktuHpPPW0hxG/n4OT36+mZNFJX7Hk1pABSES4hrXi+XnV2bwyf0XMrhDIr+ftYGLH5vL9JW5Oj4hnlJBiISJ9MT6PDcuk9fuHEDjurHc8/pyrnn6K12RLZ5RQYiEmcHtE5lx91B+d00PduafYOyTX/LIO6vIP1bodzSJMCoIkTAUHWVc3y+Nzx68kB8NSWda1i5G/mEOryzcQYluLS7VRAUhEsYaxcfyH1dmMPPeYXRr1Yifv7eaMX/+gqU78v2OJhFABSESATolN2TKHQN48qY+fH28kGueXsD901botFipEhWESIQwM67o0ZLZD1zIT0a25/2Ve7joD3N5fv5WikpK/Y4nYUgFIRJh6tWJ4d++04VZPx1OZtsm/OqDdVwxeT5Z27XbSSpHBSESodIT6/O3H/TjuXGZHDtVwrXPLODht7M5dFxnO8m5UUGIRDAz49KMZD65fzjjh7fjzaU5XPzYXN5ZlqOL7OSsVBAitUC9OjH87PKuzJg0lLSm9bh/2kpufn4RW/OO+h1NQpgKQqQWyWjViHfuGsyvrurOqt0FjHp8Po9/upFTxbq3k3ybCkKklomKMm4Z2IbZD1zIqO4tePzTTYx+fD4Lthz0O5qEGE8LwsxGmdkGM9tsZg8HmW5mNjkwPdvM+pSblmBmb5nZejNbZ2aDvMwqUts0bxjP5Bt78/Lt/Skuddz43EJ+9u4qjpws8juahAjPCsLMooEngdFABnCjmWVUGDYa6Bh4jQeeLjftCeAj51wXoCewzqusIrXZ8E5JzLpvOHcOS2fq4p1c9qd5fL5+v9+xJAR4uQXRH9jsnNvqnCsEpgJjK4wZC7zsyiwEEsyspZk1AoYDLwA45wqdc4c8zCpSq9WtE82jV2TwzsQhNIyP4YcvLeG+qct1A8BazsuCSAF2lXufE/jsXMa0A/KAv5nZcjN73szqB1uImY03sywzy8rLy6u+9CK1UK+0BGbcPZR7Lu7I+9l7uPSPc3k/W8+dqK28LAgL8lnF/8pONyYG6AM87ZzrDRwDvnUMA8A596xzLtM5l5mUlFSVvCICxMVEc/+lnZhx91BaJdRl0mvL+fErS9l3WPd1qm28LIgcIK3c+1Qg9xzH5AA5zrlFgc/foqwwRKSGdG3ZiHcnDuaR0V2YuzGPS/44l7eX6gK72sTLglgCdDSzdDOrA9wATK8wZjowLnA200CgwDm3xzm3F9hlZp0D4y4G1nqYVUSCiImO4scXtmfmvcPo0qIhD7y5kh+/spQDR0/5HU1qgGcF4ZwrBiYBsyg7A2mac26NmU0wswmBYR8CW4HNwHPAxHJfcTcwxcyygV7A//Mqq4icWbukBkwdP4ifXd6FORvyuOxP8/ho9R6/Y4nHLJI2FzMzM11WVpbfMUQi2sZ9R7h/2gpW7z7M1b1T+K8x3WhcL9bvWHKezGypcy4z2DRdSS0ildIpuSHvThzCvRd3ZPrKXL7z+DzmbtQZhJFIBSEilRYbHcVPL+3EuxMH0yA+htteXMyj767ieGGx39GkGqkgROS89UhN4P27h3LH0HReW7yTKyd/QXbOIb9jSTVRQYhIlcTHRvMfV2bw2h0DOVFUwvee+oqn52yhtDRyjm/WVioIEakWg9o3Y+a9w7g0I5nffrSem59fxJ6CE37HkipQQYhItUmoV4enbu7Db6+5gBW7DjHq8fk6HTaMqSBEpFqZGd/v15oP7hlKm2b1mPDqMh5+O1sHsMOQCkJEPNEuqQFvTRjMXSPa80bWLq6c/AWrdxf4HUsq4awFEbgNRtrZxomIVFQnJop/H9WFKXcM4Hhh2QHsVxbu0P2cwsRZC8KV/Zt8z/soIhKpBrdP5MN7hzGofTN+/t5qJr2+XE+uCwPnuotpoZn18zSJiES0pvXr8Lcf9OOhUZ35aPVexvxZu5xC3bkWxEhggZltCTw7elXgJnoiIucsKsqYOKIDr985kJNFpXzvae1yCmUx5zhutKcpRKRW6Z/elA/uGcr901by8/dWs3DrQX7zvQtoGK+b/oWScyoI59wOr4OISO3SrEEcf/tBP56Zt4XHPt7Imt0FPHVzXzJaNfI7mgToNFcR8U35XU4nikr43tNf8t7y3X7HkgAVhIj4rn96U96/exg9UhO4740V/Nf0NRSVlPodq9ZTQYhISEhqGMeUOwZw+5B0XvpqOzc/t4j9R076HatWU0GISMiIjY7iF2MyeOKGXmTvPsSVk79g6Y6v/Y5Va6kgRCTkjO2VwrsThxAfG80Nzy7QqbA+UUGISEjq2rIRMyYNZWiHRH7+3moeeiubU8UlfseqVVQQIhKyGteL5YXb+nH3RR14c2kONz+3iANHT/kdq9ZQQYhISIuKMh64rDN/vrE3q3YXMPYvX7I297DfsWoFFYSIhIUxPVvx5oRBFJeWcu0zXzFrzV6/I0U8FYSIhI0eqQlMnzSUjs0b8ONXlvLk55t18NpDKggRCSvJjeJ548eDGNurFb+ftYF7p67gZJEOXnvhXG/WJyISMuJjo3n8+73olNyQ38/awI784zw/LpOkhnF+R4so2oIQkbBkZvxkZAeeuaUvG/Ye5ntPf8nm/Uf9jhVRVBAiEtZGdW/B1PGDOFFYwjVPf8WirQf9jhQxVBAiEvZ6pSXwzl1DaNagDre+sJh/rNAdYauDCkJEIkLrZvV4567B9GqdwL1TV+gMp2qgghCRiJFQrw6v/Kg/3+1ZdobTz95dRbFuG37edBaTiESUuJiyM5xSm9TlqTlbyDtyir/c1If42Gi/o4UdbUGISMSJijIeGtWF/xnbjdnr9zPuhcUUnCjyO1bY8bQgzGyUmW0ws81m9nCQ6WZmkwPTs82sT4Xp0Wa23Mze9zKniESmcYPa8sQNvVm+62u+/9cF7D+sBxBVhmcFYWbRwJPAaCADuNHMMioMGw10DLzGA09XmH4vsM6rjCIS+b7bsxUv3NaPnfnHufaZBew4eMzvSGHDyy2I/sBm59xW51whMBUYW2HMWOBlV2YhkGBmLQHMLBW4Anjew4wiUgsM75TElDsGcORkEdc8vYA1uQV+RwoLXhZECrCr3PucwGfnOuZx4CHgjKcgmNl4M8sys6y8vLwqBRaRyNW7dRPenDCIOtHGDX9dqAvqzoGXBWFBPqt4UnLQMWZ2JbDfObf0bAtxzj3rnMt0zmUmJSWdT04RqSU6NG/IW3cNpnmjOG7722LmbdQvlWfiZUHkAGnl3qcCuec4ZgjwXTPbTtmuqYvM7FXvoopIbdEqoS5v/HgQ6YkNuOPvWXy6dp/fkUKWlwWxBOhoZulmVge4AZheYcx0YFzgbKaBQIFzbo9z7hHnXKpzrm1gvs+cc7d4mFVEapHEBnG8fucAurZsyIRXl/JB9h6/I4UkzwrCOVcMTAJmUXYm0jTn3Bozm2BmEwLDPgS2ApuB54CJXuURESkvoV4dXr1jAL3SErj79WW8uzzH70ghxyLpXiWZmZkuKyvL7xgiEkaOnSrmzpezWLD1IP/v6gu4sX9rvyPVKDNb6pzLDDZNV1KLSK1WPy6GF3/Qjws7JfHIO6t4deEOvyOFDBWEiNR68bHR/PXWvlzcpTn/8d5qpi7e6XekkKCCEBGh7CZ/T93ShxGdk3jk3VVMy9p19pkinApCRCQgLiaaZ27py9AOifz729m8s6x2H7hWQYiIlBMfG81z4zIZ1K4ZD765slY/nU4FISJSQXxsNC/c1o/+6U356RsrmLGy4jW+tYMKQkQkiLp1ykois01ZSXy+fr/fkWqcCkJE5DTqx8Xwwg8y6RK44nrxtny/I9UoFYSIyBk0jI/l7z/sT0qTuvzopSWs3l17bhWughAROYtmDeJ49UcDaBgfw20vLmZr3lG/I9UIFYSIyDlolVCXV+4YAMCtLyxmT8EJnxN5TwUhInKO2ic14O+39+fwiSJueX4R+ccK/Y7kKRWEiEgldE9pzPO3ZZLz9Ql+9PclnCwq8TuSZ1QQIiKVNKBdM564oRcrdh3i3qnLKSmNnLtil6eCEBE5D6O6t+TnV2Qwa80+fvXBWr/jeCLG7wAiIuHq9qHp5Hx9ghe/3EZqk3r8aGi635GqlQpCRKQKHr2iK7mHTvCrD9bSqnE8oy9o6XekaqNdTCIiVRAdZTx+Qy96pSVw3xsrWL7za78jVRsVhIhIFf3z5n7NG8Ux/pWl7C046XekaqGCEBGpBk3r1+GF2/px/FQx41/JiojTX1UQIiLVpFNyQ564oTerdhfwb29l41x4n/6qghARqUaXZCTz4GWdmbEyl6fmbPE7TpXoLCYRkWo2cUR7Nu47wh8+3kCn5IZcmpHsd6Tzoi0IEZFqZmb89poe9EhpzH1Tl7N5/xG/I50XFYSIiAfiY6P5662Z1K0TzYRXl3HsVLHfkSpNBSEi4pEWjeOZfENvtuYd5ZF3VoXdQWsVhIiIhwZ3SOSByzozfWUuLy/Y4XecSlFBiIh47K4L23Nxl+b86oO1LAujK61VECIiHouKMv54fS9aNI7nJ1OW8XWYPGhIBSEiUgMa14vl6Zv7cuDoKR56OzwuolNBiIjUkO4pjfn3UV34ZO0+piza6Xecs1JBiIjUoNuHpDO8UxK/fH8tG/eF9vURKggRkRoUFWU8dl1PGsbHcM/ry0P6pn6eFoSZjTKzDWa22cweDjLdzGxyYHq2mfUJfJ5mZp+b2TozW2Nm93qZU0SkJiU1jOMP1/Vk/d4j/Gbmer/jnJZnBWFm0cCTwGggA7jRzDIqDBsNdAy8xgNPBz4vBh5wznUFBgI/CTKviEjYGtG5ObcPSeelr7Yzd2Oe33GC8nILoj+w2Tm31TlXCEwFxlYYMxZ42ZVZCCSYWUvn3B7n3DIA59wRYB2Q4mFWEZEa99CoznRo3oCH387m8Mkiv+N8i5cFkQLsKvc+h2//JX/WMWbWFugNLAq2EDMbb2ZZZpaVlxeaLSwiEkx8bDSPXdeT/UdO8av31/od51u8LAgL8lnFE3/POMbMGgBvA/c55w4HW4hz7lnnXKZzLjMpKem8w4qI+KFnWgI/Ht6OaVk5fLZ+n99xvsHLgsgB0sq9TwVyz3WMmcVSVg5TnHPveJhTRMRX917SkU7JDXj47VUUHA+dXU1eFsQSoKOZpZtZHeAGYHqFMdOBcYGzmQYCBc65PWZmwAvAOufcHz3MKCLiu7iYaB67rhcHjxXym49C56wmzwrCOVcMTAJmUXaQeZpzbo2ZTTCzCYFhHwJbgc3Ac8DEwOdDgFuBi8xsReB1uVdZRUT8dkFqY24f0pbXF+8ka3u+33EAsHC4H8i5yszMdFlZWX7HEBE5L8dOFXPZn+ZRPy6a9+8eRp0Y769lNrOlzrnMYNN0JbWISIioHxfD/4ztxsZ9R3lu/la/46ggRERCycVdkxndvQWTZ29ix8FjvmZRQYiIhJj/HNON2Ogofvn+Ol9zqCBEREJMi8bx/GRkBz5dt48vNh3wLYcKQkQkBP1wSFvSmtbll++vpbik1JcMKggRkRAUHxvNz0Z3ZcO+I7yRtevsM3hABSEiEqJGdW9B//SmPPbxRl9u5qeCEBEJUWbGL67M4OvjhTz52eYaX74KQkQkhHVPacz3eqfy0lfb2Xf4ZI0uWwUhIhLi7rukIyWljr/U8FaECkJEJMSlNa3H9/ulMXXJTnblH6+x5aogRETCwKSLOmBmTJ69qcaWqYIQEQkDLRvX5daBbXh7WU6N3YJDBSEiEibGD29HTFQUz8/fViPLU0GIiISJ5EbxXN07hWlZuzhw9JTny1NBiIiEkfEXtqOwpJSXv9ru+bJUECIiYaR9UgMu7ZrM3xfs4NipYk+XpYIQEQkzE0a0p+BEEdM8vkeTCkJEJMz0ad2EXmkJvLpwB14+NloFISIShm4Z2IYtecdYuDXfs2WoIEREwtCVPVrSuG4sUxbt8GwZKggRkTAUHxvNtX1TmbVmr2envKogRETC1Pf7pVFU4pi+IteT71dBiIiEqU7JDeme0oi3l+V48v0qCBGRMDZuYFt6t06gsLj6n1sdU+3fKCIiNeb6fmlc3y/Nk+/WFoSIiASlghARkaBUECIiEpQKQkREglJBiIhIUCoIEREJSgUhIiJBqSBERCQo8/Je4jXNzPKA8721YSJwoBrjVDflqxrlqxrlq5pQztfGOZcUbEJEFURVmFmWcy7T7xyno3xVo3xVo3xVE+r5Tke7mEREJCgVhIiIBKWC+D/P+h3gLJSvapSvapSvakI9X1A6BiEiIkFpC0JERIJSQYiISFARUxBmNsrMNpjZZjN7OMh0M7PJgenZZtbnbPOaWVMz+8TMNgX+bFJu2iOB8RvM7DuhlM/M2prZCTNbEXg941O+68xsjZmVmllmhe8LhfUXNF8Irb/fm9n6wPh3zSwhxNZf0HwhtP5+GRi7wsw+NrNWIbb+guY7n/XnGedc2L+AaGAL0A6oA6wEMiqMuRyYCRgwEFh0tnmB3wEPB35+GPht4OeMwLg4ID0wf3QI5WsLrA6B9dcV6AzMATLLfVeorL/T5QuV9XcZEBP4+bch+N/f6fKFyvprVG7+e4BnQmz9nS5fpdafl69I2YLoD2x2zm11zhUCU4GxFcaMBV52ZRYCCWbW8izzjgX+Hvj578BV5T6f6pw75ZzbBmwOfE+o5KssT/I559Y55zYEWV5IrL8z5Kssr/J97JwrDsy/EEgt912hsP5Ol6+yvMp3uNz89QFX7rtCYf2dLl/IiJSCSAF2lXufE/jsXMacad5k59wegMCfzSuxPD/zAaSb2XIzm2tmw86Qzct8VVmen/kg9Nbf7ZT9hnquy/MzH4TI+jOzX5vZLuBm4BeVWJ6f+aBy688zkVIQFuSzim18ujHnMu/5LK+y46sz3x6gtXOuN3A/8JqZNQqhfFp/lchnZo8CxcCUSizPz3whs/6cc48659IC2SZVYnl+5qvs+vNMpBREDpBW7n0qkHuOY840777AZiKBP/dXYnm+5QtsOh8M/LyUsn2gnXzIV5Xl+ZYvlNafmd0GXAnc7Jz7518sIbP+guULpfVXzmvANZVYnm/5zmP9eedcD1aE8guIAbZSdsDpnweCulUYcwXfPIi0+GzzAr/nmweBfxf4uRvfPMi1lTMf5KrpfEn/zEPZwbHdQNOazldu3jl88yBwSKy/M+QLifUHjALWAkkVvisk1t8Z8oXK+utYbv67gbdCbP2dLl+l1p+XL9//cq+2f5Cyswg2Uta2jwY+mwBMCPxswJOB6av45l8I35o38HkzYDawKfBn03LTHg2M3wCMDqV8lP0msibwH+MyYIxP+a6m7DeoU8A+YFaIrb+g+UJo/W2mbP/1isDrmRBbf0HzhdD6extYDWQDM4CUEFt/QfOdz/rz6qVbbYiISFCRcgxCRESqmQpCRESCUkGIiEhQKggREQlKBSEiIkGpIETOk5klmNnEwM+tzOwtvzOJVCed5ipynsysLfC+c66731lEvBDjdwCRMPYboL2ZraDsYsWuzrnuZvYDyu6sGw10Bx6j7CraWym7KO9y51y+mbWn7OKqJOA4cKdzbn1N/0OInI52MYmcv4eBLc65XsC/VZjWHbiJsts9/xo47spuvrYAGBcY8yxwt3OuL/Ag8FRNhBY5V9qCEPHG5865I8ARMyug7FYKUHYbhh5m1gAYDLxp9q8bfsbVfEyR01NBiHjjVLmfS8u9L6Xs/7so4FBg60MkJGkXk8j5OwI0PJ8ZXdnTxLaZ2XXwr2ca96zOcCJVpYIQOU+u7J79X5rZaspuvV5ZNwM/MrOVlN29c2x15hOpKp3mKiIiQWkLQkREglJBiIhIUCoIEREJSgUhIiJBqSBERCQoFYSIiASlghARkaD+Pz7lU8R6mOlmAAAAAElFTkSuQmCC\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZoAAAEWCAYAAABfdFHAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAqe0lEQVR4nO3dd5xU5dn/8c9F720pgvTeUVhA1GCJsSAqiCaW2KMm4ZeoeaIssWFsYEwsUcNDQEOKGsOCgkgQiYoiliXKFlg60mHpZYFt1++POeZZCWUX9uzM7Hzfr9e+ZvY+5565blj2y5lz5hpzd0RERMJSKdoFiIhIxaagERGRUCloREQkVAoaEREJlYJGRERCpaAREZFQKWhEQmBmWWZ27lG2nWtm68u3IpHoUdCIhMDde7j7B+X1fGbW2czeMrMcM9thZrPNrEt5Pb/IsShoRCqGBsB0oAvQDPgceCuaBYl8Q0EjEgIzW2NmFwT3a5rZn8xsp5ktBvqX9fO5++fuPsndd7h7PvAM0MXMksr6uURKq0q0CxBJAA8DHYKv2sCsY+1sZulA66NsftXdf1qC5xwMbHb37aUpVCQMChqR8H0f+Km77wB2mNnzwENH29nde5/Mk5lZS+BF4Bcn8zgiZUUvnYmErwWwrtj3X4f1RGbWBHgXeMndXwvreURKQ0EjEr5NQKti3x/tZTHgP5dG7zvK1/hjzGtIJGSmu/vjZVO6yMnTS2ci4XsDGG1mnxE5R/OzY+3s7j1K+wRmVg+YDcx395QTqlIkJDqiEQnfI0ReLltN5IjjLyE8x3AiV7PdctgR0DGPnkTKg+mDz0REJEw6ohERkVApaEREJFQKGhERCZWCRkREQpXQlzc3btzY27ZtG+0yRETiysKFC7e5e5OS7p/QQdO2bVvS0tKiXYaISFwxs1J1t9BLZyIiEioFjYiIhEpBIyIioVLQiIhIqBQ0IiISqlCDxszuMrPMoO353cFYHzNbYGYZZjYj6DpbornB+NXBWJGZJR82Z7SZrTCzpWZ2UZhrExGRkgktaMysJ3A7MADoAww1s07ARCDF3XsB04B7SzEXIBO4Eph32JzuwDVAD+Bi4CUzqxzC0kREpBTCPKLpBnzq7rnuXgB8SKSVeRf+LyTmACNKMRd3X+LuS48w5wrgdXc/5O6rgRVEgkpERALuzt+/WMt7i7eU23OGGTSZwGAzSzKzWsAQIp8ymAlcHuxzNd/+5MHjzT2WU/n2x+WuD8a+xczuMLM0M0vLyckp1YJEROLZ2u25XD/xM0alZvDmVxvK7XlD6wzg7kvMbByRo5Z9wCKgALgVeN7MHgKmA3mlmHssdqQyjvDYE4AJAMnJyfowHhGp8AqLnD99soanZy+lciXj8eE9ubZ/+X0mXqgtaNx9EjAJwMyeANa7ezZwYTDWGbi0pHOP83Tr+fZRT0tg48nULyIS75Zt2ct9U9L5at0uzu/alMeH96R5/ZrlWkOoQWNmTd19a/BxslcCg4qNVQIeAMaXdO5xnm468KqZ/Q5oAXQCPi+zxYiIxJG8giL+8MFKXnh/OXVrVOW5a07j8j4tMDvSiz/hCrupZqqZJQH5wEh33xlctjwy2D4VeAXAzFoAE919yNHmBvsNB34PNAFmmtlX7n6Ru2eZ2RvAYiIvs41098KQ1yciEnMWrdvFqNR0sjfv5fI+LXj4su4k1aketXrMPXFPUyQnJ7u6N4tIRXEgr5Bn3lvGxI9W0bRuDR4b1pMLujcr8+cxs4Xunnz8PSMS+mMCREQqigUrtzN6ajprtudy3cDWpFzSlXo1qka7LEBBIyIS1/YczGfsrGxe/WwtbZJq8ertAzmzQ+Nol/UtChoRkTg1d8kW7p+Wyda9B7ljcHvuuaAzNavFXkMUBY2ISJzZvu8Qj8xYzPRFG+nSrC7jb+jHaa0aRLuso1LQiIjECXdn+qKNPDJjMXsP5nPPBZ35ybkdqFYlthvxK2hEROLApt0HeGBaJnOzt9KnVQOeGtGbLqfUjXZZJaKgERGJYUVFzutfrOPJd5aQX1TEA5d245az2lG5Uvm/8fJEKWhERGLUmm37SZmazqerdnBmhyTGXtmb1km1ol1WqSloRERiTEFhES/PX81v311GtcqVGHtlL37Qv1VU2seUBQWNiEgMyd68h1FT0lm0fjcXdGvGY8N6ckr9GtEu66QoaEREYsChgkJefH8lL72/gvo1q/L7a09naO/mcXsUU5yCRkQkyr5cu5NRqeks27KP4aefyoNDu9OodrVol1VmFDQiIlGSm1fAb99dxsvzV3NKvRq8cnN/zuvaNNpllTkFjYhIFHyyYhspUzNYuyOXH57RmlEXd6VujDTBLGsKGhGRcrT7QD5PvrOE179YR9ukWvz9jjMY2D4p2mWFSkEjIlJO3s3azANvZrJt3yHuPCfSBLNG1dhrglnWFDQiIiHbtu8QY6Zn8Xb6JrqeUpeJNyXTu2WDaJdVbhQ0IiIhcXfe/GoDj8xYTO6hQn55YWfuPKcDVSvHdhPMsqagEREJwYZdB7h/WgYfLM2hb+sGPHVVbzo2jY8mmGVNQSMiUoaKipy/fb6Wse8socjh4cu6c+OgtnHVBLOsKWhERMrIqpx9pKRm8PmaHZzdsTFPXtmLVo3irwlmWVPQiIicpILCIiZ+vJpn5iyjepVKPHVVb67u17JCtI8pCwoaEZGTsHjjHu5LXUTmhj1c1KMZj17Rk6b14rsJZllT0IiInICD+YW88K8VjP9wJQ1qVeMP1/flkl7No11WTFLQiIiU0sKvd3DflHRW5uxnRN+WPDi0Gw1qVZwmmGVNQSMiUkL7DxXwm9lLmbxgDS3q12TyrQM4p3OTaJcV8xQ0IiIlMG9ZDqOnZrBx9wFuPKMN917clTrV9Su0JPSnJCJyDLtz83l05mKmLFxP+ya1eePOQfRv2yjaZcUVBY2IyFH8M3MTD76VxY79efz03A78/LudEqIJZllT0IiIHGbr3oM8/FYWszI30715PV65uT89T60f7bLiloJGRCTg7qT+ewOPvr2YA/mF3HdxF27/TvuEa4JZ1kL90zOzu8ws08yyzOzuYKyPmS0wswwzm2Fm9Uo6NxhvZGZzzGx5cNswGK9qZpODx11iZqPDXJuIVCzrduRy48uf88t/LKJzszrMuus7/PTcjgqZMhDan6CZ9QRuBwYAfYChZtYJmAikuHsvYBpwbynmAqQAc929EzA3+B7gaqB68Lj9gDvNrG1IyxORCqKoyPnT/NVc9Ow8/v31Tn59RQ/+fscgOjSpE+3SKowwo7ob8Km757p7AfAhMBzoAswL9pkDjCjFXIArgMnB/cnAsOC+A7XNrApQE8gD9pTpikSkQlmxdR/f/98FjJmxmOS2jZh9z2BuHNSWSgncaTkMYQZNJjDYzJLMrBYwBGgVjF8e7HN1MFbSuQDN3H0TQHDbNBifAuwHNgFrgafdfcfhD2xmd5hZmpml5eTklMU6RSTO5BcW8eL7Kxjy3Ecs37qPp6/uw+Rb+tOyoTothyG0iwHcfYmZjSNy1LIPWAQUALcCz5vZQ8B0IkceJZ17LAOAQqAF0BD4yMzec/dVhz32BGACQHJysp/4CkUkHmVu2M19U9JZvGkPl/ZqzpjLe9CkbvVol1WhhXrVmbtPAiYBmNkTwHp3zwYuDMY6A5eWdG6waYuZNXf3TWbWHNgajF8H/NPd84GtZjYfSAZWISIJ72B+Ic/NXc6EeatoVLsa43/Yj4t7nhLtshJC2FedNQ1uWwNXAq8VG6sEPACML+ncYNN04Kbg/k3AW8H9tcD5FlEbOAPILus1iUj8+WLNDoY89xF/+GAlI/qeynv3nKOQKUdhv48m1cySgHxgpLvvDC5bHhlsnwq8AmBmLYCJ7j7kaHOD8bHAG2Z2G5FwuToYfzF4rEzAgFfcPT3k9YlIDNt3qICn/pnNnxd8TcuGNfnrbQM5u1PjaJeVcMw9cU9TJCcne1paWrTLEJEQfLB0K/dPy2Tj7gPcfGZbfnlhF2qrCWaZMLOF7p5c0v31py4iFcrO/Xk8OnMxU/+9gY5N6zDlx2fSr03DaJeV0BQ0IlIhuDuzMjfz0FuZ7MrN52fnd+T/nd+R6lXUBDPaFDQiEve27jnIg29lMjtrC71Orc+fbx1I9xZH7G4lUaCgEZG45e78I209j81czKGCIlIu6cqPzm5HFfUniykKGhGJS+t25DJ6agYfr9jGgLaNGDuiF+3VnywmKWhEJK4UFjmTP1nDb2YvpZLBo8N6cv2A1upPFsMUNCISN5Zv2cuo1HT+vXYX53ZpwuPDe3Fqg5rRLkuOQ0EjIjEvv7CI8R+s5Pf/WkHt6pV59genccVpLTDTUUw8UNCISEzLWL+be6csInvzXob2jjTBbFxHTTDjiYJGRGLSwfxCnnlvGX+ct4rGdaoz4YZ+XNhD/cnikYJGRGLOp6u2M3pqBqu37eea/q0YPaQb9WtWjXZZcoIUNCISM/YezGfsrGz+9tlaWjWqyd9+NJCzOqoJZrxT0IhITHg/eyu/mpbB5j0Hue3sdvzPhZ2pVU2/oioC/S2KSFTt2J/Hr2dk8eZXG+nUtA6pPzmTvq3VBLMiUdCISFS4O2+nb2LM9Cx2H8jn59/txMjzOqgJZgWkoBGRcrdlz0Hun5bJe0u20Ltlff76o4F0a64mmBWVgkZEyo278/cv1vH4O0vIKyji/iHduOWstmqCWcEpaESkXHy9fT8pqRksWLWdge0aMW5Eb9o2rh3tsqQcKGhEJFSFRc4r81fz9LtLqVKpEk8M78U1/VupCWYCUdCISGiWbt7LfanpLFq3i+92bcpjw3vSvL6aYCYaBY2IlLm8giJe+mAFL76/gro1qvLcNadxeR81wUxUChoRKVOL1u3ivinpLN2ylytOa8FDQ7uTpCaYCU1BIyJl4kBeIb+bs5RJH6+mad0aTLwxmQu6N4t2WRIDFDQictI+WbmNlNQM1u7I5bqBrUm5pCv1aqgJpkQoaETkhO05mM+T72Tz2udraZNUi9duP4NBHZKiXZbEGAWNiJyQ9xZv4f43M8jZe4g7Brfnngs6U7Oa2sfIf1PQiEipbN93iEdmLGb6oo10PaUuE25Ipk+rBtEuS2KYgkZESsTdmb5oI2OmZ7HvUAH3XNCZn5zbgWpV1D5Gjk1BIyLHtWn3AR6Ylsnc7K2c1qoBT13Vm87N6ka7LIkTChoROaqiIue1L9by5DvZFBY5Dw7tzs1ntqWy2sdIKYR6zGtmd5lZppllmdndwVgfM1tgZhlmNsPMjtgb/Ehzg/FGZjbHzJYHtw2LbesdPHZW8Pg1wlyfSEW2ett+rv3jp9w/LZM+reoz++7B3HZ2O4WMlFpoQWNmPYHbgQFAH2ComXUCJgIp7t4LmAbcW4q5ACnAXHfvBMwNvsfMqgB/BX7s7j2Ac4H8sNYnUlEVFBYxYd5KLn52Hos37WHciF789baBtE6qFe3SJE6FeUTTDfjU3XPdvQD4EBgOdAHmBfvMAUaUYi7AFcDk4P5kYFhw/0Ig3d0XAbj7dncvLNsliVRs2Zv3cOUfPuGJd7IZ3LkJ7/3iHH7Qv7V6lMlJCTNoMoHBZpZkZrWAIUCrYPzyYJ+rg7GSzgVo5u6bAILbpsF4Z8DNbLaZ/dvM7jtSUWZ2h5mlmVlaTk5OGSxTJP4dKijkd3OWMfT5j9mw8wAvXHc6E27oR7N6evVZTl5oFwO4+xIzG0fkqGUfsAgoAG4Fnjezh4DpQF4p5h5LFeBsoD+QC8w1s4XuPvewx54ATABITk72E1+hSMXw77U7GTUlneVb93Hl6afy4NDuNKxdLdplSQUS6lVn7j4JmARgZk8A6909m8jLXJhZZ+DSks4NNm0xs+buvsnMmgNbg/H1wIfuvi2Y8w7Ql8h5HBE5TG5eAU/PXsYrn6ymeb0avHJLf87r0vT4E0VKKeyrzpoGt62BK4HXio1VAh4Axpd0brBpOnBTcP8m4K3g/mygt5nVCi4MOAdYXNZrEqkI5q/YxkXPzuPl+av54cA2zL5nsEJGQhP2+2hSzSyJyNVfI919Z3DZ8shg+1TgFQAzawFMdPchR5sbjI8F3jCz24C1RM7zEDz274AvAAfecfeZIa9PJK7sPpDPEzOX8Pe0dbRrXJu/33EGA9urCaaEy9wT9zRFcnKyp6WlRbsMkXIxO2szD76Zyfb9edz+nfbcfUEnalRVE0wpveD8d3JJ91dnAJEKLmfvIcZMz2Jmxia6Na/HpJv606tl/WiXJQlEQSNSQbk7077cwK/fXkzuoUJ+eWFn7jynA1UrqwmmlC8FjUgFtGHXAe6flsEHS3Po2zrSBLNjUzXBlOhQ0IhUIEVFzt8++5qxs7JxYMxl3blhkJpgSnQpaEQqiJU5+xidmsHna3bwnU6NeWJ4L1o1Un8yiT4FjUicKygsYsJHq3j2veXUqFKJ31zVm6v6tVR/MokZChqROJa1cTejUtPJ3LCHi3o049EretJU/ckkxihoROLQwfxCfv+v5Yz/cBUNa1XjD9f35ZJezaNdlsgRKWhE4kzamh2MSk1nZc5+RvRtyYNDu9GglppgSuxS0IjEif2HCvjN7KVMXrCGFvVrMvnWAZzTuUm0yxI5LgWNSByYtyyH0VMz2Lj7ADee0YZ7L+5Kner65yvxQT+pIjFsV24ej81cwpSF62nfpDZv3DmI/m0bRbsskVJR0IjEqFkZm3jwrSx25ubx03M78PPvqgmmxKfjBo1FLsZv6e7ryqEekYS3de9BHn4ri1mZm+nevB5/uqU/PU9VE0yJX8cNGnd3M3sT6Bd+OSKJy92ZsnA9j769mIMFRdx7URfuGNxeTTAl7pX0pbNPzay/u38RajUiCWrdjlx+NS2Dj5ZvI7lNQ8aO6E3HpnWiXZZImShp0JwH3GlmXwP7ASNysNM7tMpEEkBRkfPnBWt4avZSDPj1FT344cA2VFITTKlASho0l4RahUgCWrF1L6NSM1j49U4Gd27CE8N70rKhmmBKxVOioHH3r8MuRCRR5BcWMWHeKp57bzk1q1Xmt1f34cq+p6oJplRYurxZpBxlbtjNfVPSWbxpD5f2as6Yy3vQpG71aJclEioFjUg5OJhfyHNzlzNh3ioa1a7G+B/24+Kep0S7LJFyoaARCdnnq3eQkprOqm37+X5yS+4f0p36tapGuyyRcqOgEQnJvkMFjJuVzV8+/ZqWDWvy19sGcnanxtEuS6TcKWhEQvD+0q3cPzWDTXsOcstZbfnlhV2orSaYkqD0ky9Shnbuz+PRtxcz9csNdGxahyk/PpN+bRpGuyyRqFLQiJQBd+edjM08PD2TXbn5/Pz8jow8vyPVq6gJpoiCRuQkbd1zkAfezOTdxVvodWp9/nzrQLq3qBftskRihoJG5AS5O/9IW8+jMxeTV1DE6Eu6ctvZ7aiiJpgi36KgETkBa7fnMnpaOvNXbGdAu0aMG9Gbdo1rR7sskZikoBEphcIi50+frOHp2UupXMl4bFhPrhvQWk0wRY5BQSNSQsu37OW+1HS+XLuL87o04fHhvWjRoGa0yxKJeaG+mGxmd5lZppllmdndwVgfM1tgZhlmNsPMjnjW9Ehzg/FGZjbHzJYHtw0Pm9fazPaZ2S/DXJskjryCIp6fu5xLn/+YNdv28+wPTuPlm/srZERKKLSgMbOewO3AAKAPMNTMOgETgRR37wVMA+4txVyAFGCuu3cC5gbfF/cMMKvsVySJKH39Li5/4WN+N2cZF/U8hTm/OIdhp6vTskhphHlE0w341N1z3b0A+BAYDnQB5gX7zAFGlGIuwBXA5OD+ZGDYN5PMbBiwCsgq05VIwjmYX8iT7yxh2Ivz2Zmbxx9vTOb3155O4zrqtCxSWmEGTSYw2MySzKwWMARoFYxfHuxzdTBW0rkAzdx9E0Bw2xTAzGoDo4BHjlWUmd1hZmlmlpaTk3NSC5SK6dNV27n42Xn877xVfD+5Fe/ecw7f694s2mWJxK3QLgZw9yVmNo7IUcs+YBFQANwKPG9mDwHTgbxSzD2WR4Bn3H3fsV7WcPcJwASA5ORkL+26pOLaezCfsbOy+dtna2ndqBav/mggZ3ZUE0yRkxXqVWfuPgmYBGBmTwDr3T0buDAY6wxcWtK5waYtZtbc3TeZWXNgazA+ELjKzJ4CGgBFZnbQ3V8IZXFSofwrewv3T8tky56D/OjsdvzPhV2oWU3tY0TKQqhBY2ZN3X2rmbUGrgQGFRurBDwAjC/p3GDTdOAmYGxw+xaAu3+n2NwxwD6FjBzPjv15/HpGFm9+tZHOzerw0vVncnprNcEUKUthv48m1cySgHxgpLvvDC5bHhlsnwq8AmBmLYCJ7j7kaHOD8bHAG2Z2G7CWyHkekVJxd2akb2LM9Cz2Hsznru92YuR5HalWRe1jRMqauSfuaYrk5GRPS0uLdhlSzjbvjjTBfG/JFvq0rM+4q3rT9RQ1wRQpKTNb6O7JJd1fnQEkYbg7r3+xjidmLiG/qIgHLu3GLWe1o7Lax4iESkEjCeHr7ftJSc1gwartDGqfxNgRvWiTpCaYIuVBQSMVWmGR88r81Tz97lKqVqrEk1f24pr+rfTOfpFypKCRCmvp5kgTzEXrdnFBt6Y8NqwXp9SvEe2yRBKOgkYqnLyCIl76YAUvvr+CujWq8vy1p3NZ7+Y6ihGJEgWNVChfrdvFqCnpLN2yl2GnteChy3rQqHa1aJclktAUNFIhHMgr5LfvLuXl+atpVq8GL9+czPld1Z9MJBYoaCTufbJyGympGazdkcv1A1uTcklX6taoGu2yRCSgoJG4tedgPk++s4TXPl9H26RavH7HGZzRPinaZYnIYRQ0EpfeW7yF+9/MIGfvIe4c3J67L+isJpgiMUpBI3Fl+75DjJmxmBmLNtL1lLr88cZkerdsEO2yROQYFDQSF9yd6Ys2MmZ6FvsPFfI/3+vMned0UBNMkTigoJGYt3HXAR54M5N/ZW/l9NYNeGpEbzo1qxvtskSkhBQ0ErOKipxXP1/L2FnZFBY5Dw3tzk1ntlUTTJE4o6CRmLR6235SUtP5bPUOzuqYxJPDe9M6qVa0yxKRE6CgkZhSUFjEpI9X87s5y6hWpRJPjejN1ckt1T5GJI4paCRmLNm0h1Gp6aSv3833ujfjsWE9aVZPTTBF4p2CRqLuUEEhL/5rBS99sJIGtary4nV9GdLrFB3FiFQQChqJqoVf72RUajortu7jyr6n8uCl3WmoJpgiFYqCRqIiN6+A38xeyp8+WUPzejV45Zb+nNelabTLEpEQKGik3H28fBspU9NZv/MANw5qw30Xd6VOdf0oilRU+tct5Wb3gXwen7mYN9LW065xbd64cxAD2jWKdlkiEjIFjZSL2VmbefDNTLbvz+Mn53bgru92okZVNcEUSQQKGglVzt5DjJmexcyMTXRrXo9JN/WnV8v60S5LRMqRgkZC4e5M+3IDv357MbmHCrn3oi7cMbg9VSurCaZIolHQSJnbsOsAv5qawYfLcujXpiHjRvSmY9M60S5LRKJEQSNlpqjI+etnXzNuVjYOjLmsOzcOakslNcEUSWgKGikTK3P2kZKazhdrdvKdTo15YngvWjVSE0wRUdDISSooLGLCR6t49r3l1KhSid9c1Zur+qkJpoj8HwWNnLCsjbsZlZpO5oY9XNzjFH49rAdN66oJpoh8W6iXAJnZXWaWaWZZZnZ3MNbHzBaYWYaZzTCzeiWdG4w3MrM5ZrY8uG0YjH/PzBYGj7vQzM4Pc22J7GB+Ib+Znc3lL8xn8+5D/OH6voy/oZ9CRkSOKLSgMbOewO3AAKAPMNTMOgETgRR37wVMA+4txVyAFGCuu3cC5gbfA2wDLgse9ybgL2GtLZEt/HoHlz7/ES++v5Lhp5/Ke78YzCW9mke7LBGJYWEe0XQDPnX3XHcvAD4EhgNdgHnBPnOAEaWYC3AFMDm4PxkYBuDuX7r7xmA8C6hhZtXLdkmJa/+hAsZMz+Kq8Qs4mF/En28dwNNX96FBLXVaFpFjCzNoMoHBZpZkZrWAIUCrYPzyYJ+rg7GSzgVo5u6bAILbI7X8HQF86e6HDt9gZneYWZqZpeXk5JzE8hLHvGU5XPjMPCYvWMNNg9ry7j2DGdy5SbTLEpE4EdrFAO6+xMzGETlq2QcsAgqAW4HnzewhYDqQV4q5x2VmPYBxwIVHqWsCMAEgOTnZS7mshLIrN4/HZi5hysL1tG9Sm3/cOYjktmqCKSKlE+pVZ+4+CZgEYGZPAOvdPZsgBMysM3BpSecGm7aYWXN332RmzYGt38wxs5ZEzvvc6O4rw1lVYpiVsYkH38piZ24eI8/rwM/OVxNMETkxoQaNmTV1961m1hq4EhhUbKwS8AAwvqRzg03TiZzsHxvcvhXs3wCYCYx29/lhrqsi27r3IA+/lcWszM30aFGPybf2p0cLNcEUkRMX9vtoUs0sCcgHRrr7zuCy5ZHB9qnAKwBm1gKY6O5DjjY3GB8LvGFmtwFriZznAfh/QEfgQTN7MBi70N3/c8QjR+fuTFm4nsdmLuFAfiGjLu7K7d9pRxU1wRSRk2TuiXuaIjk52dPS0qJdRtSt25HLr6Zl8NHybfRv25CxI3rToYmaYIrIkZnZQndPLun+6gyQwIqKnD8vWMNTs5diwKNX9OD6gW3UBFNEypSCJkGt2LqXUakZLPx6J+d0bsLjw3vSsqGaYIpI2VPQJJj8wiImzFvFc+8tp1b1yvzu+30YfvqpaoIpIqFR0CSQzA27uXdKOks27eHS3s0Zc1kPmtRV8wQRCZeCJgEczC/k2feW88ePVtGodjX+94Z+XNTjlGiXJSIJQkFTwX2+egcpqems2rafHyS34ldDulG/VtVolyUiCURBU0HtO1TAuFnZ/OXTr2nVqCZ/vW0gZ3dqHO2yRCQBKWgqoPeXbuX+qRls2nOQW89qxy8v6kytavqrFpHo0G+fCmTn/jwefXsxU7/cQKemdUj9yZn0bd0w2mWJSIJT0FQA7s7MjE08/FYWuw/k8/PzOzLy/I5Ur6ImmCISfQqaOLdlz0EefDOTdxdvodep9fnrjwbSrfkRPx1bRCQqFDRxyt15I20dj81cQl5BEb8a0pVbz1ITTBGJPQqaOLR2ey6jp6Uzf8V2BrZrxLgRvWnbuHa0yxIROSIFTRwpLHL+9Mkanp69lMqVjMeH9+Ta/q3VBFNEYpqCJk4s27KX+6ak89W6XZzftSmPD+9J8/o1o12WiMhxKWhiXF5BEeM/XMnv/7WcOtWr8Nw1p3F5nxZqgikicUNBE8MWrdvFqNR0sjfv5fI+LXj4su4k1VETTBGJLwqaGHQgr5Bn3lvGxI9W0bRuDSbemMwF3ZtFuywRkROioIkxC1ZuZ/TUdNZsz+XaAa0ZPaQr9WqoCaaIxC8FTYzYczCfsbOyefWztbRJqsWrtw/kzA5qgiki8U9BEwP+lb2FX03NZOveg9z+nXb84ntdqFlN7WNEpGJQ0ETR9n2H+PXbi3nrq410aVaX8Tf047RWDaJdlohImVLQRIG7M33RRh6ZsZi9B/O5+4JO/PTcjlSrovYxIlLxKGjK2abdB3hgWiZzs7fSp1UDnhrRmy6n1I12WSIioVHQlJOiIuf1L9bx5DtLyC8q4oFLu3HLWe2orPYxIlLBKWjKwZpt+0mZms6nq3YwqH0SY0f0ok2SmmCKSGJQ0ISosMh5+ePV/HbOUqpWqsTYK3vxg/6t1D5GRBKKgiYkSzfv5b4pi1i0fjcXdGvKY8N6cUr9GtEuS0Sk3CloytihgkJeen8lL32wgno1qvL7a09naO/mOooRkYSloClDX67dyajUdJZt2cew01rw0GU9aFS7WrTLEhGJKgVNGcjNK+C37y7j5fmrOaVeDV6+OZnzu6oJpogIQKjvEDSzu8ws08yyzOzuYKyPmS0wswwzm2Fm9Uo6NxhvZGZzzGx5cNuw2LbRZrbCzJaa2UVhru0bn6zYxsXPfsSkj1dz3YDWvHvPYIWMiEgxoQWNmfUEbgcGAH2AoWbWCZgIpLh7L2AacG8p5gKkAHPdvRMwN/geM+sOXAP0AC4GXjKz0BqG7T6QT0pqOtdN/IxKBq/fcQaPD+9FXXVaFhH5ljCPaLoBn7p7rrsXAB8Cw4EuwLxgnznAiFLMBbgCmBzcnwwMKzb+ursfcvfVwAoiQVXm0tfv4sJnPuSNtHXceU57/nn3YM5onxTGU4mIxL0wgyYTGGxmSWZWCxgCtArGLw/2uToYK+lcgGbuvgkguG0ajJ8KrCv2GOuDsW8xszvMLM3M0nJyck5oYa0b1aJzs7q8OfIsRl/SjRpV1WlZRORoQgsad18CjCNy1PJPYBFQANwKjDSzhUBdIK8Uc4/lSNcP+xEee4K7J7t7cpMmTUq+oGIa1KrGX24bSO+WDU5ovohIIgn1YgB3n+Tufd19MLADWO7u2e5+obv3A14DVpZ0brBpi5k1Bwhutwbj6/n20VFLYGPZr0pEREoj7KvOmga3rYErgdeKjVUCHgDGl3RusGk6cFNw/ybgrWLj15hZdTNrB3QCPi/rNYmISOmE/T6aVDNLAvKBke6+M7hseWSwfSrwCoCZtQAmuvuQo80NxscCb5jZbcBaIud5cPcsM3sDWEzkZbaR7l4Y8vpEROQ4zP2/TmMkjOTkZE9LS4t2GSIiccXMFrp7ckn310c6iohIqBQ0IiISKgWNiIiESkEjIiKhSuiLAcwsB/j6JB6iMbCtjMqJB4m2XtCaE4XWXDpt3L3E73hP6KA5WWaWVporL+Jdoq0XtOZEoTWHSy+diYhIqBQ0IiISKgXNyZkQ7QLKWaKtF7TmRKE1h0jnaEREJFQ6ohERkVApaEREJFQJFTRmdrGZLTWzFWaWcoTtZmbPB9vTzazv8eaaWSMzm2Nmy4PbhsW2jQ72X2pmFxUb72dmGcG2583sSB/aVmHWbGa1zGymmWWbWZaZjQ1rvbGy5sOeb7qZZYax1uPVXWx7ef1sVzOzCWa2LPj7PtJHtVe0NV8b/HtON7N/mlnjirBmi3zC8ftmts/MXjjseUr3O8zdE+ILqEzkQ9baA9WIfGpn98P2GQLMIvJpnWcAnx1vLvAUkBLcTwHGBfe7B/tVB9oF8ysH2z4HBgXPMwu4pCKvGagFnBfsUw34qKKvudhzXQm8CmQmyM/2I8Bjwf1KQOOKvGYiH7Wy9Zt1BvPHVJA11wbOBn4MvHDY85Tqd1giHdEMAFa4+yp3zwNeB644bJ8rgD97xKdAA4t8iuex5l4BTA7uTwaGFRt/3d0PuftqYAUwIHi8eu6+wCN/Y38uNqesxcSa3T3X3d8HCB7r30Q+ATUMMbFmADOrA/wCeCyEdRYXM2sm8lHtTwK4e5G7h/Vu+1hZswVftYP/1dcjvE/2Ldc1u/t+d/8YOFj8CU7kd1giBc2pwLpi368Pxkqyz7HmNnP3TQDBbdMSPNb649RRVmJlzf9hZg2Ay4C5pVtKicXSmh8FfgvknshCSiEm1hz83QI8amb/NrN/mFmzE1rR8cXEmt09H/gJkEEkYLoDk05sScdV3ms+Vh2l+h2WSEFzpNcQD7+2+2j7lGRuSZ/vRB7rRMXKmiMbzaoQ+Uju59191XEe60TFxJrN7DSgo7tPO878shATaybyMlJLYL679wUWAE8f57FOVEys2cyqEgma04EWQDow+jiPdaLKe80nU8e3JFLQrAdaFfu+Jf99iHu0fY41d0twKPnNIeXWEjxWyyOMhyFW1vyNCcByd3+2tAsphVhZ8yCgn5mtAT4GOpvZBye0ouOLlTVvJ3L09k24/gPoSzhiZc2nAbj7yuBlpDeAM09oRcdX3ms+Vh2l+x12oiem4u2LyP+2VhE5kffNybAeh+1zKd8+kfb58eYCv+HbJ9KeCu734NsnD1fxfydMvwge/5sTaUMSYM2PAalApUT5ey72fG0J92KAmFkzkdf+zw/u3wz8oyKvmchRzCagSbDfo8BvK8Kaiz3mzfz3xQCl+h0W2j/4WPwickXGMiJXX9wfjP0Y+HFw34AXg+0ZQPKx5gbjSUTONywPbhsV23Z/sP9Sil2VASQDmcG2Fwg6NFTUNRP5H48DS4Cvgq8fVeQ1H1ZPW0IMmlhaM9AGmEfkJaS5QOsEWPOPg5/tdGAGkFSB1rwG2AHsI3Ik882VaqX6HaYWNCIiEqpEOkcjIiJRoKAREZFQKWhERCRUChoREQmVgkZEREKloBEpR2bWwMx+GtxvYWZTol2TSNh0ebNIOTKztsDb7t4z2rWIlJcq0S5AJMGMBTqY2VdE3iDXzd17mtnNRDrgVgZ6EmnGWQ24AThE5J3XO8ysA5E35DUh0u7ldnfPLu9FiJSGXjoTKV8pwEp3Pw2497BtPYHriLR0fxzIdffTiTSnvDHYZwLwM3fvB/wSeKk8ihY5GTqiEYkd77v7XmCvme0m0s4EIq1Eegefb3Mm8I9iH2hYvfzLFCkdBY1I7DhU7H5Rse+LiPxbrQTsCo6GROKGXjoTKV97gbonMtHd9wCrzexq+M/nw/cpy+JEwqCgESlH7r4dmG9mmUTas5fW9cBtZrYIyOK/P8pXJObo8mYREQmVjmhERCRUChoREQmVgkZEREKloBERkVApaEREJFQKGhERCZWCRkREQvX/AZuuXuwzC9wiAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] @@ -139,157 +139,141 @@ " \n", " 0\n", " 0.00000\n", - " -3.427842e-05\n", - " 0.000000\n", + " 0.000505\n", " 0.0\n", " 0.0\n", - " 2.000000e-08\n", + " 0.0\n", + " 0.0001\n", " 39.478418\n", - " -0.000034\n", + " 0.000505\n", " \n", " \n", " 1\n", " 0.00001\n", - " -3.427842e-05\n", - " 0.000000\n", + " 0.000505\n", + " 0.0\n", " 0.0\n", " 0.0\n", - " 2.000000e-08\n", + " 0.0001\n", " 39.478418\n", - " -0.000034\n", + " 0.000505\n", " \n", " \n", " 2\n", " 0.00002\n", - " -3.427842e-05\n", - " 0.000000\n", + " 0.000505\n", " 0.0\n", " 0.0\n", - " 2.000000e-08\n", + " 0.0\n", + " 0.0001\n", " 39.478418\n", - " -0.000034\n", + " 0.000505\n", " \n", " \n", " 3\n", " 0.00003\n", - " -3.427842e-05\n", - " 0.000000\n", + " 0.000505\n", + " 0.0\n", " 0.0\n", " 0.0\n", - " 2.000000e-08\n", + " 0.0001\n", " 39.478418\n", - " -0.000034\n", + " 0.000505\n", " \n", " \n", " 4\n", " 0.00004\n", - " -3.427842e-05\n", - " 0.000000\n", + " 0.000505\n", " 0.0\n", " 0.0\n", - " 2.000000e-08\n", + " 0.0\n", + " 0.0001\n", " 39.478418\n", - " -0.000034\n", + " 0.000505\n", " \n", " \n", - " ...\n", - " ...\n", - " ...\n", - " ...\n", - " ...\n", - " ...\n", - " ...\n", - " ...\n", - " ...\n", + " 5\n", + " 0.00005\n", + " 0.000505\n", + " 0.0\n", + " 0.0\n", + " 0.0\n", + " 0.0001\n", + " 39.478418\n", + " 0.000505\n", " \n", " \n", - " 377\n", - " 0.00377\n", - " -3.427842e-05\n", - " 0.000000\n", + " 6\n", + " 0.00006\n", + " 0.000505\n", + " 0.0\n", " 0.0\n", " 0.0\n", - " 2.000000e-08\n", + " 0.0001\n", " 39.478418\n", - " -0.000034\n", + " 0.000505\n", " \n", " \n", - " 378\n", - " 0.00378\n", - " -3.427842e-05\n", - " 0.000000\n", + " 7\n", + " 0.00007\n", + " 0.000505\n", " 0.0\n", " 0.0\n", - " 2.000000e-08\n", + " 0.0\n", + " 0.0001\n", " 39.478418\n", - " -0.000034\n", + " 0.000505\n", " \n", " \n", - " 379\n", - " 0.00379\n", - " -3.427842e-05\n", - " 0.000000\n", + " 8\n", + " 0.00008\n", + " 0.000505\n", + " 0.0\n", " 0.0\n", " 0.0\n", - " 2.000000e-08\n", + " 0.0001\n", " 39.478418\n", - " -0.000034\n", + " 0.000505\n", " \n", " \n", - " 380\n", - " 0.00380\n", - " -3.427842e-05\n", - " 0.000000\n", + " 9\n", + " 0.00009\n", + " 0.000505\n", " 0.0\n", " 0.0\n", - " 2.000000e-08\n", + " 0.0\n", + " 0.0001\n", " 39.478418\n", - " -0.000034\n", + " 0.000505\n", " \n", " \n", - " 381\n", - " 0.00381\n", - " 5.066057e-13\n", - " 0.000034\n", + " 10\n", + " 0.00010\n", + " 0.000505\n", + " 0.0\n", " 0.0\n", " 0.0\n", - " 2.000000e-08\n", + " 0.0001\n", " 39.478418\n", - " -0.000034\n", + " 0.000505\n", " \n", " \n", "\n", - "

382 rows × 8 columns

\n", "" ], "text/plain": [ - " t Eorbit Ecollisions Lx Ly Lz Mtot \\\n", - "0 0.00000 -3.427842e-05 0.000000 0.0 0.0 2.000000e-08 39.478418 \n", - "1 0.00001 -3.427842e-05 0.000000 0.0 0.0 2.000000e-08 39.478418 \n", - "2 0.00002 -3.427842e-05 0.000000 0.0 0.0 2.000000e-08 39.478418 \n", - "3 0.00003 -3.427842e-05 0.000000 0.0 0.0 2.000000e-08 39.478418 \n", - "4 0.00004 -3.427842e-05 0.000000 0.0 0.0 2.000000e-08 39.478418 \n", - ".. ... ... ... ... ... ... ... \n", - "377 0.00377 -3.427842e-05 0.000000 0.0 0.0 2.000000e-08 39.478418 \n", - "378 0.00378 -3.427842e-05 0.000000 0.0 0.0 2.000000e-08 39.478418 \n", - "379 0.00379 -3.427842e-05 0.000000 0.0 0.0 2.000000e-08 39.478418 \n", - "380 0.00380 -3.427842e-05 0.000000 0.0 0.0 2.000000e-08 39.478418 \n", - "381 0.00381 5.066057e-13 0.000034 0.0 0.0 2.000000e-08 39.478418 \n", - "\n", - " Etot \n", - "0 -0.000034 \n", - "1 -0.000034 \n", - "2 -0.000034 \n", - "3 -0.000034 \n", - "4 -0.000034 \n", - ".. ... \n", - "377 -0.000034 \n", - "378 -0.000034 \n", - "379 -0.000034 \n", - "380 -0.000034 \n", - "381 -0.000034 \n", - "\n", - "[382 rows x 8 columns]" + " t Eorbit Ecollisions Lx Ly Lz Mtot Etot\n", + "0 0.00000 0.000505 0.0 0.0 0.0 0.0001 39.478418 0.000505\n", + "1 0.00001 0.000505 0.0 0.0 0.0 0.0001 39.478418 0.000505\n", + "2 0.00002 0.000505 0.0 0.0 0.0 0.0001 39.478418 0.000505\n", + "3 0.00003 0.000505 0.0 0.0 0.0 0.0001 39.478418 0.000505\n", + "4 0.00004 0.000505 0.0 0.0 0.0 0.0001 39.478418 0.000505\n", + "5 0.00005 0.000505 0.0 0.0 0.0 0.0001 39.478418 0.000505\n", + "6 0.00006 0.000505 0.0 0.0 0.0 0.0001 39.478418 0.000505\n", + "7 0.00007 0.000505 0.0 0.0 0.0 0.0001 39.478418 0.000505\n", + "8 0.00008 0.000505 0.0 0.0 0.0 0.0001 39.478418 0.000505\n", + "9 0.00009 0.000505 0.0 0.0 0.0 0.0001 39.478418 0.000505\n", + "10 0.00010 0.000505 0.0 0.0 0.0 0.0001 39.478418 0.000505" ] }, "execution_count": 7, diff --git a/examples/symba_energy_momentum/escape.in b/examples/symba_energy_momentum/escape.in new file mode 100644 index 000000000..41a44b9a6 --- /dev/null +++ b/examples/symba_energy_momentum/escape.in @@ -0,0 +1,18 @@ +3 +1 39.47841760435743 +0.0 0.0 0.0 +0.0 0.0 0.0 +0.4 0.4 0.4 !Ip +0.0 0.0 0.0 !rot +2 1e-07 0.0009 +7e-06 +99.9 0.0 0.0 +100.00 10.00 0.0 +0.4 0.4 0.4 !Ip +0.0 0.0 0.0 !rot +3 1e-08 0.0004 +3.25e-06 +1.0 4.20E-05 0.0 +0.00 -6.28 0.0 +0.4 0.4 0.4 !Ip +0.0 0.0 0.0 !rot diff --git a/examples/symba_energy_momentum/param.escape.in b/examples/symba_energy_momentum/param.escape.in new file mode 100644 index 000000000..ee8628cca --- /dev/null +++ b/examples/symba_energy_momentum/param.escape.in @@ -0,0 +1,33 @@ +T0 0.0e0 +TSTOP 0.00100 ! simulation length in seconds = 100 years +DT 0.00010 ! stepsize in seconds +PL_IN escape.in +TP_IN tp.in +IN_TYPE ASCII +ISTEP_OUT 1 ! output cadence every year +BIN_OUT bin.escape.dat +PARTICLE_FILE particle.escape.dat +OUT_TYPE REAL8 ! double precision real output +OUT_FORM XV ! osculating element output +OUT_STAT REPLACE +ISTEP_DUMP 1 ! system dump cadence +J2 0.0 ! no J2 term +J4 0.0 ! no J4 term +CHK_CLOSE yes ! check for planetary close encounters +CHK_RMIN 0.005 +CHK_RMAX 1e2 +CHK_EJECT -1.0 ! ignore this check +CHK_QMIN -1.0 ! ignore this check +!CHK_QMIN_COORD HELIO ! commented out here +!CHK_QMIN_RANGE 1.0 1000.0 ! commented out here +ENC_OUT enc.escape.dat +EXTRA_FORCE no ! no extra user-defined forces +BIG_DISCARD no ! output all planets if anything discarded +RHILL_PRESENT yes ! Hill's sphere radii in input file +MTINY 1.0e-16 +FRAGMENTATION yes +MU2KG 1.98908e30 +TU2S 3.1556925e7 +DU2M 1.49598e11 +ENERGY yes +ROTATION yes diff --git a/examples/symba_energy_momentum/param.merger.in b/examples/symba_energy_momentum/param.merger.in index 9a51b03c7..7c52e6ba0 100644 --- a/examples/symba_energy_momentum/param.merger.in +++ b/examples/symba_energy_momentum/param.merger.in @@ -1,8 +1,6 @@ ! ! Parameter file for Mars in SI units ! -NPLMAX 3 ! not used -NTPMAX 0 ! not used T0 0.0e0 TSTOP 0.10 ! simulation length in seconds = 100 years DT 0.01 ! stepsize in seconds @@ -11,6 +9,7 @@ TP_IN tp.in IN_TYPE ASCII ISTEP_OUT 1 ! output cadence every year BIN_OUT bin.merger.dat +PARTICLE_FILE particle.merger.dat OUT_TYPE REAL8 ! double precision real output OUT_FORM XV ! osculating element output OUT_STAT REPLACE diff --git a/examples/symba_energy_momentum/sun.in b/examples/symba_energy_momentum/sun.in index 8359affdb..7fe65bfa5 100644 --- a/examples/symba_energy_momentum/sun.in +++ b/examples/symba_energy_momentum/sun.in @@ -1,4 +1,4 @@ -2 +3 1 39.47841760435743 0.0 0.0 0.0 0.0 0.0 0.0 @@ -10,7 +10,7 @@ -10.00 2.00 0.0 0.4 0.4 0.4 !Ip 0.0 0.0 0.0 !rot -3 1e-08 0.0004 +3 1e-05 0.0004 3.25e-06 1.0 4.20E-05 0.0 0.00 -6.28 0.0 diff --git a/src/io/io_conservation_report.f90 b/src/io/io_conservation_report.f90 index 159cce80b..40c8580bd 100644 --- a/src/io/io_conservation_report.f90 +++ b/src/io/io_conservation_report.f90 @@ -32,7 +32,8 @@ module subroutine io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param, "; DM/M0 = ", ES12.5)' associate(Ecollisions => symba_plA%helio%swiftest%Ecollisions, Lescape => symba_plA%helio%swiftest%Lescape, Mescape => symba_plA%helio%swiftest%Mescape, & - mass => symba_plA%helio%swiftest%mass, dMcb => symba_plA%helio%swiftest%dMcb, Mcb_initial => symba_plA%helio%swiftest%Mcb_initial) + Eescape => symba_plA%helio%swiftest%Eescape, mass => symba_plA%helio%swiftest%mass, dMcb => symba_plA%helio%swiftest%dMcb, & + Mcb_initial => symba_plA%helio%swiftest%Mcb_initial) if (lfirst) then if (param%out_stat == "OLD") then open(unit = egyiu, file = energy_file, form = "formatted", status = "old", action = "write", position = "append") @@ -58,8 +59,8 @@ module subroutine io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param, Lmag_now = norm2(Ltot_now) Lerror = norm2(Ltot_now - Ltot_orig) / Lmag_orig Eorbit_error = (Eorbit - Eorbit_orig) / abs(Eorbit_orig) - Ecoll_error = -Ecollisions / abs(Eorbit_orig) - Etotal_error = (Eorbit - (Eorbit_orig + Ecollisions)) / abs(Eorbit_orig) + Ecoll_error = Ecollisions / abs(Eorbit_orig) + Etotal_error = (Eorbit - Ecollisions - Eorbit_orig - Eescape) / abs(Eorbit_orig) Merror = (Mtot_now - Mtot_orig) / Mtot_orig write(*, egytermfmt) Lerror, Ecoll_error, Etotal_error, Merror if (Ecoll_error > 0.0_DP) then diff --git a/src/main/swiftest_symba.f90 b/src/main/swiftest_symba.f90 index 77c1f5609..b39890096 100644 --- a/src/main/swiftest_symba.f90 +++ b/src/main/swiftest_symba.f90 @@ -21,7 +21,7 @@ program swiftest_symba logical :: lfrag_add, ldiscard_pl, ldiscard_tp integer(I4B) :: nplm, ntp, ntp0, nsppl, nsptp, iout, idump, iloop, i integer(I4B) :: nplplenc, npltpenc, nmergeadd, nmergesub - real(DP) :: t, tfrac, tbase, msys + real(DP) :: t, tfrac, tbase, msys, Eorbit_orig real(DP) :: Ecollision, Eorbit_before, Eorbit_after, ke_orbit_before, ke_spin_before, ke_orbit_after, ke_spin_after, pe_before, pe_after real(DP), dimension(NDIM) :: Ltot character(STRMAX) :: inparfile @@ -153,6 +153,7 @@ program swiftest_symba if (param%lenergy) then call coord_h2b(npl, symba_plA%helio%swiftest, msys) call io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param, lterminal=.false.) + call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_before, ke_spin_before, pe_before, Eorbit_orig, Ltot) end if write(*, *) " *************** Main Loop *************** " @@ -175,7 +176,6 @@ program swiftest_symba call symba_collision(t, symba_plA, nplplenc, plplenc_list, lfrag_add, mergeadd_list, nmergeadd, param) ldiscard_pl = ldiscard_pl .or. lfrag_add - ! These next two blocks should be streamlined/improved but right now we treat discards separately from collisions so it has to be this way if (ldiscard_pl .or. ldiscard_tp) then if (param%lenergy) call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_before, ke_spin_before, pe_before, Eorbit_before, Ltot) call symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA, nmergeadd, mergeadd_list, discard_plA, & @@ -190,7 +190,7 @@ program swiftest_symba if (param%lenergy) then call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_after, ke_spin_after, pe_after, Eorbit_after, Ltot) - Ecollision = Eorbit_before - Eorbit_after ! Energy change resulting in this collisional event Total running energy offset from collision in this step + Ecollision = Eorbit_after - Eorbit_before ! Energy change resulting in this collisional event Total running energy offset from collision in this step symba_plA%helio%swiftest%Ecollisions = symba_plA%helio%swiftest%Ecollisions + Ecollision end if !if (ntp > 0) call util_dist_index_pltp(nplm, ntp, symba_plA, symba_tpA) diff --git a/src/modules/module_interfaces.f90 b/src/modules/module_interfaces.f90 index a5918240f..1481104cb 100644 --- a/src/modules/module_interfaces.f90 +++ b/src/modules/module_interfaces.f90 @@ -877,14 +877,14 @@ END SUBROUTINE symba_discard_sun_pl END INTERFACE INTERFACE - subroutine symba_discard_conserve_mtm(param, swiftest_plA, ipl, lescape) + subroutine symba_discard_conserve_mtm(param, swiftest_plA, ipl, lescape_body) use swiftest_globals use swiftest_data_structures implicit none type(user_input_parameters), intent(inout) :: param type(swiftest_pl), intent(inout) :: swiftest_plA integer(I4B), intent(in) :: ipl - logical, intent(in) :: lescape + logical, intent(in) :: lescape_body end subroutine END INTERFACE diff --git a/src/modules/swiftest_data_structures.f90 b/src/modules/swiftest_data_structures.f90 index 46b30ae6a..5aebc5a87 100644 --- a/src/modules/swiftest_data_structures.f90 +++ b/src/modules/swiftest_data_structures.f90 @@ -52,6 +52,7 @@ module swiftest_data_structures real(DP) :: Rcb_initial !! Initial radius of the central body real(DP) :: dRcb = 0.0_DP!! Change in the radius of the central body real(DP) :: Ecollisions = 0.0_DP !! Energy lost from system due to collisions + real(DP) :: Eescape = 0.0_DP !! Energy gained from system due to escaped bodies integer(I4B), dimension(:,:), allocatable :: k_plpl integer(I8B) :: num_plpl_comparisons contains diff --git a/src/symba/symba_discard_conserve_mtm.f90 b/src/symba/symba_discard_conserve_mtm.f90 index a18e54e71..df478d598 100644 --- a/src/symba/symba_discard_conserve_mtm.f90 +++ b/src/symba/symba_discard_conserve_mtm.f90 @@ -1,4 +1,4 @@ -subroutine symba_discard_conserve_mtm(param, swiftest_plA, ipl, lescape) +subroutine symba_discard_conserve_mtm(param, swiftest_plA, ipl, lescape_body) !! author: David A. Minton !! !! Conserves system momentum when a body is lost from the system or collides with central body @@ -10,26 +10,19 @@ subroutine symba_discard_conserve_mtm(param, swiftest_plA, ipl, lescape) type(swiftest_pl), intent(inout) :: swiftest_plA type(user_input_parameters), intent(inout) :: param integer(I4B), intent(in) :: ipl - logical, intent(in) :: lescape + logical, intent(in) :: lescape_body ! Internals real(DP), dimension(NDIM) :: Lpl, Lcb, xcom, vcom real(DP) :: pe, ke_orbit, ke_spin + integer(I4B) :: i associate(npl => swiftest_plA%nbody, xb => swiftest_plA%xb, vb => swiftest_plA%vb, & rot => swiftest_plA%rot, Ip => swiftest_plA%Ip, radius => swiftest_plA%radius, mass => swiftest_plA%mass, & Lcb_initial => swiftest_plA%Lcb_initial, dLcb => swiftest_plA%dLcb, Ecollisions => swiftest_plA%Ecollisions, & Rcb_initial => swiftest_plA%Rcb_initial, dRcb => swiftest_plA%dRcb, & - Mcb_initial => swiftest_plA%Mcb_initial, dMcb => swiftest_plA%dMcb, Mescape => swiftest_plA%Mescape) - - xcom(:) = (mass(ipl) * xb(:, ipl) + mass(1) * xb(:,1)) / (mass(1) + mass(ipl)) - vcom(:) = (mass(ipl) * vb(:, ipl) + mass(1) * vb(:,1)) / (mass(1) + mass(ipl)) - - call util_crossproduct(xb(:,ipl) - xcom(:), vb(:,ipl) - vcom(:), Lpl) - Lpl(:) = mass(ipl) * (Lpl(:) + radius(ipl)**2 * Ip(3,ipl) * rot(:, ipl)) - - call util_crossproduct(xb(:, 1) - xcom(:), vb(:, 1) - vcom(:), Lcb) - Lcb(:) = mass(1) * Lcb(:) + Mcb_initial => swiftest_plA%Mcb_initial, dMcb => swiftest_plA%dMcb, & + Mescape => swiftest_plA%Mescape, Lescape => swiftest_plA%Lescape, Eescape => swiftest_plA%Eescape) ! Add the potential and kinetic energy of the lost body to the records pe = -mass(1) * mass(ipl) / norm2(xb(:, ipl) - xb(:, 1)) @@ -37,36 +30,52 @@ subroutine symba_discard_conserve_mtm(param, swiftest_plA, ipl, lescape) ke_spin = 0.5_DP * mass(ipl) * radius(ipl)**2 * Ip(3, ipl) * dot_product(rot(:, ipl), rot(:, ipl)) ! Add the pre-collision ke of the central body to the records - ke_spin = 0.5_DP * mass(1) * radius(1)**2 * Ip(3, 1) * dot_product(rot(:, 1), rot(:, 1)) + ke_spin ! Add planet mass to central body accumulator - if (lescape) then + if (lescape_body) then Mescape = Mescape + mass(ipl) - ke_orbit = 0.0_DP + call util_crossproduct(xb(:,ipl), vb(:,ipl), Lpl) + Lpl(:) = mass(ipl) * (Lpl(:) + radius(ipl)**2 * Ip(3,ipl) * rot(:, ipl)) + Lescape(:) = Lescape(:) + Lpl(:) + do i = 2, npl + if (i == ipl) cycle + pe = pe - mass(i) * mass(ipl) / norm2(xb(:, ipl) - xb(:, i)) + end do else + xcom(:) = (mass(ipl) * xb(:, ipl) + mass(1) * xb(:,1)) / (mass(1) + mass(ipl)) + vcom(:) = (mass(ipl) * vb(:, ipl) + mass(1) * vb(:,1)) / (mass(1) + mass(ipl)) + + call util_crossproduct(xb(:,ipl) - xcom(:), vb(:,ipl) - vcom(:), Lpl) + Lpl(:) = mass(ipl) * (Lpl(:) + radius(ipl)**2 * Ip(3,ipl) * rot(:, ipl)) + + call util_crossproduct(xb(:, 1) - xcom(:), vb(:, 1) - vcom(:), Lcb) + Lcb(:) = mass(1) * Lcb(:) ke_orbit = 0.5_DP * mass(1) * dot_product(vb(:, 1), vb(:, 1)) + ke_orbit + ke_spin = ke_spin + 0.5_DP * mass(1) * radius(1)**2 * Ip(3, 1) * dot_product(rot(:, 1), rot(:, 1)) ! Update mass of central body to be consistent with its total mass dMcb = dMcb + mass(ipl) dRcb = dRcb + 1.0_DP / 3.0_DP * (radius(ipl) / radius(1))**3 - 2.0_DP / 9.0_DP * (radius(ipl) / radius(1))**6 mass(1) = Mcb_initial + dMcb radius(1) = Rcb_initial + dRcb param%rmin = radius(1) - ! Update position and velocity of central body + ! Add planet angular momentum to central body accumulator + dLcb(:) = Lpl(:) + Lcb(:) + dLcb(:) + ! Update rotation of central body to by consistent with its angular momentum + rot(:,1) = (Lcb_initial(:) + dLcb(:)) / (Ip(3, 1) * mass(1) * radius(1)**2) + ke_spin = ke_spin - 0.5_DP * mass(1) * radius(1)**2 * Ip(3, 1) * dot_product(rot(:, 1), rot(:, 1)) xb(:, 1) = xcom(:) vb(:, 1) = vcom(:) ke_orbit = ke_orbit - 0.5_DP * mass(1) * dot_product(vb(:, 1), vb(:, 1) ) end if - - ! Add planet angular momentum to central body accumulator - dLcb(:) = Lpl(:) + Lcb(:) + dLcb(:) - - ! Update rotation of central body to by consistent with its angular momentum - rot(:,1) = (Lcb_initial(:) + dLcb(:)) / (Ip(3, 1) * mass(1) * radius(1)**2) - - ! Remove the new kinetic energy of the central body from the records - ke_spin = ke_spin - 0.5_DP * mass(1) * radius(1)**2 * Ip(3, 1) * dot_product(rot(:, 1), rot(:, 1)) + ! Update position and velocity of central body ! We must do this for proper book-keeping, since we can no longer track this body's contribution to energy directly - Ecollisions = Ecollisions - 2 * (ke_orbit + ke_spin + pe) + if (lescape_body) then + Ecollisions = Ecollisions + ke_orbit + ke_spin + pe + Eescape = Eescape - (ke_orbit + ke_spin + pe) + else + Ecollisions = Ecollisions + pe + Eescape = Eescape - pe + end if ! Update the heliocentric coordinates of everything else call coord_b2h(npl, swiftest_plA)