From 41b4eb781f62bce7747fc613271bfd51b95d7424 Mon Sep 17 00:00:00 2001 From: David Minton Date: Fri, 21 May 2021 17:04:14 -0400 Subject: [PATCH] Fixed energy balancing due to collisions with the central body --- .../collision_visualization.ipynb | 244 ++++++++++++++++-- src/io/io_conservation_report.f90 | 3 +- src/symba/symba_discard_conserve_mtm.f90 | 48 ++-- 3 files changed, 254 insertions(+), 41 deletions(-) diff --git a/examples/symba_energy_momentum/collision_visualization.ipynb b/examples/symba_energy_momentum/collision_visualization.ipynb index 80d8991e7..e3e542221 100644 --- a/examples/symba_energy_momentum/collision_visualization.ipynb +++ b/examples/symba_energy_momentum/collision_visualization.ipynb @@ -26,11 +26,11 @@ "output_type": "stream", "text": [ "Reading Swiftest file param.sun.in\n", - "Reading in time 5.500e-03\n", + "Reading in time 3.810e-03\n", "Creating Dataset\n", "\n", "Adding particle info Dataset\n", - "Successfully converted 551 output frames.\n" + "Successfully converted 382 output frames.\n" ] } ], @@ -46,16 +46,39 @@ "metadata": {}, "outputs": [], "source": [ - "egy = pd.read_csv(\"energy.out\")" + "ds['r'] = np.sqrt(ds['px']**2+ds['py']**2+ds['pz']**2)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + }, + { + "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", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], "source": [ - "egy['Etot'] = egy['Eorbit'] + egy['Ecollisions']" + "ds['r'].sel(id=2).plot(x='time')" ] }, { @@ -64,39 +87,218 @@ "metadata": {}, "outputs": [], "source": [ - "ds['r'] = np.sqrt(ds['px']**2+ds['py']**2+ds['pz']**2)" + "egy = pd.read_csv(\"energy.out\")" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, + "outputs": [], + "source": [ + "egy['Etot'] = egy['Eorbit'] - egy['Ecollisions']" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "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", + "
tEorbitEcollisionsLxLyLzMtotEtot
00.00000-3.427842e-050.0000000.00.02.000000e-0839.478418-0.000034
10.00001-3.427842e-050.0000000.00.02.000000e-0839.478418-0.000034
20.00002-3.427842e-050.0000000.00.02.000000e-0839.478418-0.000034
30.00003-3.427842e-050.0000000.00.02.000000e-0839.478418-0.000034
40.00004-3.427842e-050.0000000.00.02.000000e-0839.478418-0.000034
...........................
3770.00377-3.427842e-050.0000000.00.02.000000e-0839.478418-0.000034
3780.00378-3.427842e-050.0000000.00.02.000000e-0839.478418-0.000034
3790.00379-3.427842e-050.0000000.00.02.000000e-0839.478418-0.000034
3800.00380-3.427842e-050.0000000.00.02.000000e-0839.478418-0.000034
3810.003815.066057e-130.0000340.00.02.000000e-0839.478418-0.000034
\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]" ] }, - "execution_count": 6, + "execution_count": 7, "metadata": {}, "output_type": "execute_result" - }, - { - "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", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" } ], "source": [ - "ds['r'].sel(id=2).plot(x='time')" + "egy" ] }, { diff --git a/src/io/io_conservation_report.f90 b/src/io/io_conservation_report.f90 index 0463cb193..159cce80b 100644 --- a/src/io/io_conservation_report.f90 +++ b/src/io/io_conservation_report.f90 @@ -59,13 +59,12 @@ module subroutine io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param, 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) + Etotal_error = (Eorbit - (Eorbit_orig + Ecollisions)) / 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 write(*,*) 'Something has gone wrong! Collisional energy should not be positive!' end if - end if ke_orb_last = ke_orbit ke_spin_last = ke_spin diff --git a/src/symba/symba_discard_conserve_mtm.f90 b/src/symba/symba_discard_conserve_mtm.f90 index 929bd3580..a18e54e71 100644 --- a/src/symba/symba_discard_conserve_mtm.f90 +++ b/src/symba/symba_discard_conserve_mtm.f90 @@ -13,11 +13,14 @@ subroutine symba_discard_conserve_mtm(param, swiftest_plA, ipl, lescape) logical, intent(in) :: lescape ! Internals real(DP), dimension(NDIM) :: Lpl, Lcb, xcom, vcom - reaL(DP) :: pe, ke + real(DP) :: pe, ke_orbit, ke_spin + 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, & - statipl => swiftest_plA%status(ipl), xbpl => swiftest_plA%xb(:,ipl), xbcb => swiftest_plA%xb(:,1)) + 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)) @@ -28,33 +31,42 @@ subroutine symba_discard_conserve_mtm(param, swiftest_plA, ipl, lescape) call util_crossproduct(xb(:, 1) - xcom(:), vb(:, 1) - vcom(:), Lcb) Lcb(:) = mass(1) * Lcb(:) + ! Add the potential and kinetic energy of the lost body to the records + pe = -mass(1) * mass(ipl) / norm2(xb(:, ipl) - xb(:, 1)) + ke_orbit = 0.5_DP * mass(ipl) * dot_product(vb(:, ipl), vb(:, ipl)) + 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 - ! Add the potential and kinetic energy of the lost body to the records - pe = -mass(1) * mass(ipl) / norm2(xb(:, ipl) - xb(:, 1)) - ke = 0.5_DP * mass(ipl) * dot_product(vb(:, ipl), vb(:, ipl)) + Mescape = Mescape + mass(ipl) + ke_orbit = 0.0_DP else - ! Add the potential energy of the lost body to the records - pe = -mass(1) * mass(ipl) / norm2(xb(:, ipl) - xb(:, 1)) - ke = 0.0_DP - swiftest_plA%dMcb = swiftest_plA%dMcb + mass(ipl) - swiftest_plA%dRcb = swiftest_plA%dRcb + 1.0_DP / 3.0_DP * (radius(ipl) / radius(1))**3 - 2.0_DP / 9.0_DP * (radius(ipl) / radius(1))**6 + ke_orbit = 0.5_DP * mass(1) * dot_product(vb(:, 1), vb(:, 1)) + ke_orbit ! Update mass of central body to be consistent with its total mass - mass(1) = swiftest_plA%Mcb_initial + swiftest_plA%dMcb - radius(1) = swiftest_plA%Rcb_initial + swiftest_plA%dRcb + 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 + xb(:, 1) = xcom(:) + vb(:, 1) = vcom(:) + ke_orbit = ke_orbit - 0.5_DP * mass(1) * dot_product(vb(:, 1), vb(:, 1) ) end if - swiftest_plA%Ecollisions = swiftest_plA%Ecollisions - (ke + pe) ! Add planet angular momentum to central body accumulator - swiftest_plA%dLcb(:) = Lpl(:) + Lcb(:) + swiftest_plA%dLcb(:) + dLcb(:) = Lpl(:) + Lcb(:) + dLcb(:) ! Update rotation of central body to by consistent with its angular momentum - swiftest_plA%rot(:,1) = (swiftest_plA%Lcb_initial(:) + swiftest_plA%dLcb(:)) / (Ip(3, 1) * mass(1) * radius(1)**2) + rot(:,1) = (Lcb_initial(:) + dLcb(:)) / (Ip(3, 1) * mass(1) * radius(1)**2) - ! Update position and velocity of central body - xb(:, 1) = xcom(:) - vb(:, 1) = vcom(:) + ! 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)) + + ! 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) ! Update the heliocentric coordinates of everything else call coord_b2h(npl, swiftest_plA)