From 050386cd7d9577695dec86d89023ec96aef44f11 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 12 Aug 2021 22:11:44 -0400 Subject: [PATCH] By carefully reconstructing the plplenc_list object in rearray, and keeping track of recursion levels, I've enabled the ability to resolve fragments at the moment they are identified, and not wait until the recursion steps are done. --- examples/symba_mars_disk/mars.in | 24 +- examples/symba_mars_disk/param.in | 4 +- .../1pl_1pl_encounter/param.swifter.in | 1 + .../swiftest_vs_swifter.ipynb | 254 +++++++++--------- src/modules/swiftest_classes.f90 | 2 + src/modules/symba_classes.f90 | 2 +- src/setup/setup.f90 | 6 + src/symba/symba_encounter_check.f90 | 4 + src/symba/symba_util.f90 | 39 ++- src/util/util_append.f90 | 2 + src/util/util_copy.f90 | 2 + src/util/util_spill.f90 | 2 + 12 files changed, 197 insertions(+), 145 deletions(-) diff --git a/examples/symba_mars_disk/mars.in b/examples/symba_mars_disk/mars.in index 8760a493d..447d1e308 100644 --- a/examples/symba_mars_disk/mars.in +++ b/examples/symba_mars_disk/mars.in @@ -1,4 +1,16 @@ 1500 ! Mars System in SI units +727 1.71022032e+06 2.13948145e+04 ! particle number mass Rhill +1.24108926e+04 !particle radius in m +-8.12608230e+06 -4.37306608e+06 -9.62736144e+03 ! x y z +9.87984575e+02 -1.88769371e+03 1.06882012e+01 ! vx vy vz +0.4 0.4 0.4 ! Ip +0.0 0.0 0.0 ! rot +231 6.25152932e+05 1.58916481e+04 ! particle number mass Rhill +8.87389776e+03 !particle radius in m +-8.21586374e+06 -4.28792953e+06 2.41010139e+04 ! x y z +1.01581225e+03 -1.90933511e+03 -2.60449634e+00 ! vx vy vz +0.4 0.4 0.4 ! Ip +0.0 0.0 0.0 ! rot 2 9.90685589e+04 8.35558297e+03 ! particle number mass Rhill 7.07643092e+03 !particle radius in m -2.35807426e+06 8.60445552e+06 1.25224401e+04 ! x y z @@ -1373,12 +1385,6 @@ -1.24261883e+03 1.71209694e+03 -6.95777672e+00 ! vx vy vz 0.4 0.4 0.4 ! Ip 0.0 0.0 0.0 ! rot -231 6.25152932e+05 1.58916481e+04 ! particle number mass Rhill -8.87389776e+03 !particle radius in m --8.21586374e+06 -4.28792953e+06 2.41010139e+04 ! x y z -1.01581225e+03 -1.90933511e+03 -2.60449634e+00 ! vx vy vz -0.4 0.4 0.4 ! Ip -0.0 0.0 0.0 ! rot 232 1.02687634e+05 1.54462392e+04 ! particle number mass Rhill 4.85987440e+03 !particle radius in m 1.62736579e+07 2.82256969e+06 3.66384128e+04 ! x y z @@ -4349,12 +4355,6 @@ -1.87648989e+03 1.17295601e+03 -5.20044045e-01 ! vx vy vz 0.4 0.4 0.4 ! Ip 0.0 0.0 0.0 ! rot -727 1.71022032e+06 2.13948145e+04 ! particle number mass Rhill -1.24108926e+04 !particle radius in m --8.12608230e+06 -4.37306608e+06 -9.62736144e+03 ! x y z -9.87984575e+02 -1.88769371e+03 1.06882012e+01 ! vx vy vz -0.4 0.4 0.4 ! Ip -0.0 0.0 0.0 ! rot 728 7.56089690e+05 4.18623276e+04 ! particle number mass Rhill 9.45460600e+03 !particle radius in m -2.34951111e+05 2.32053308e+07 -5.22036528e+04 ! x y z diff --git a/examples/symba_mars_disk/param.in b/examples/symba_mars_disk/param.in index e2be7eece..ef35236ba 100644 --- a/examples/symba_mars_disk/param.in +++ b/examples/symba_mars_disk/param.in @@ -22,9 +22,9 @@ CHK_QMIN_RANGE 3389500.0 338950000000.0 EXTRA_FORCE no BIG_DISCARD no RHILL_PRESENT yes -GMTINY 1000.0 +GMTINY 1000.0 ENERGY yes -FRAGMENTATION yes +FRAGMENTATION no ROTATION yes MU2KG 1.0 DU2M 1.0 diff --git a/examples/symba_swifter_comparison/1pl_1pl_encounter/param.swifter.in b/examples/symba_swifter_comparison/1pl_1pl_encounter/param.swifter.in index 853815639..a67348c0e 100644 --- a/examples/symba_swifter_comparison/1pl_1pl_encounter/param.swifter.in +++ b/examples/symba_swifter_comparison/1pl_1pl_encounter/param.swifter.in @@ -24,3 +24,4 @@ ENC_OUT enc.swifter.dat EXTRA_FORCE no BIG_DISCARD no RHILL_PRESENT yes +ENERGY yes diff --git a/examples/symba_swifter_comparison/1pl_1pl_encounter/swiftest_vs_swifter.ipynb b/examples/symba_swifter_comparison/1pl_1pl_encounter/swiftest_vs_swifter.ipynb index ec0e145ef..79543cb09 100644 --- a/examples/symba_swifter_comparison/1pl_1pl_encounter/swiftest_vs_swifter.ipynb +++ b/examples/symba_swifter_comparison/1pl_1pl_encounter/swiftest_vs_swifter.ipynb @@ -81,8 +81,8 @@ { "data": { "text/plain": [ - "[,\n", - " ]" + "[,\n", + " ]" ] }, "execution_count": 6, @@ -91,7 +91,7 @@ }, { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAERCAYAAABl3+CQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAgUklEQVR4nO3de3Rd5X3m8e9jybZsSeZiXXyRjWVJNjaXOEQ1AbIIBMiA2+IGJhk8mQk0Sb3IrdPJpFN3WNO0zUpCFsmapI3bjpM0A2kSr4TWJQnGYEJS0qSUGMLFYIyNbbBsgWWDMTa+yfrNH+fICHPOsSydc/Y+0vNZS8tn7/3ufX4SHD3a+333uxURmJmZ5TMm6QLMzCzdHBRmZlaQg8LMzApyUJiZWUEOCjMzK8hBYWZmBY3YoJD095J2SVpfpON9SdL67Nd/KsYxzcwqwYgNCuD/AVcX40CSfhu4AFgAXAj8saRJxTi2mVnajdigiIgHgZcHrpPUJmmNpEck/ULS2YM83HzgXyKiNyIOAI9TpBAyM0u7ERsUeawAPhUR7wA+A/zNIPd7HLhG0kRJDcDlwIwS1WhmlirVSRdQLpLqgIuBH0rqXz0+u+064C9z7LYjIv5DRNwn6beAXwE9wL8BvaWv2swseRrJcz1JmgX8JCLOzfYpbIyIqUU47veAf4iI1cM9lplZ2o2aS08RsQ/YKun9AMp422D2lVQlaXL29fnA+cB9JSvWzCxFRuwZhaTvA5cBDcBLwGeBB4C/BaYCY4GVEZHrktOJx6oBHs0u7gNujojHil+1mVn6jNigMDOz4hg1l57MzGxoRuSop4aGhpg1a1bSZZiZVYxHHnlkd0Q05to2IoNi1qxZrFu3LukyzMwqhqTn823zpSczMyvIQWFmZgU5KMzMrCAHhZmZFeSgMDOzghwUZmZWkIPCzMwKSjQoJF0taaOkzZKW5dguSX+V3f6EpAuSqNPM7GR27z/MD9dtT7qMkkgsKCRVAcuBa8g8QW6JpPknNLsG6Mh+LSUzoZ+ZWer8yZ1P8Md3PsG23QeSLqXokrwzeyGwOSK2AEhaCSwGnh7QZjFwR2RmLnxI0umSpkZEdykKemz7Xvo8SaKZnaIN3fv46TO7AHjmxX3MaqhNuKLiSjIopgMDz9O6gAsH0WY6UJKgeOIbNzOHbaU4tJmNYG3AV+rb+cz+JTzz4mtcfe6wn4+WKkkGhXKsO/HP+cG0yTSUlpK5PMXMmTOHVNAV85qpfWXPkPY1s9Ht7TOn8dcbJrLxxdeSLqXokgyKLmDGgOUWYOcQ2gAQESuAFQCdnZ1Dun40fcnXhrKbmRkAZ7/yyIgMiiRHPf0a6JDUKmkccAPwoxPa/Aj4UHb00zuBV0vVP2FmNlxzp9Szbc8BDh09lnQpRZVYUEREL/BJ4F5gA/CDiHhK0s2Sbs42Ww1sATYD3wA+nkixZmaDcPaUevoCNr20P+lSiirR51FExGoyYTBw3d8NeB3AJ8pdl5nZUMydUg9kRj6d13JawtUUj+/MNjMrkrMm11IzdsyI66dwUJiZFUnVGNHRVM8zDgozM8tn7hQHhZmZFXD2lHp27z/Mnv2Hky6laBwUZmZFdPaUSQAjqp/CQWFmVkRvjHxyUJiZWQ6N9eOZXDvOZxRmZpbf3Cn1PPOSg8LMzPKYO6WeZ198jb6+kfHYAgeFmVmRnT2lnoNHj/HCy68nXUpROCjMzIpsbnbk00jp0HZQmJkV2ZzmOqSRM0TWQWFmVmQTx1Vz1pkT2fjSvqRLKQoHhZlZCYykqTwcFGZmJTB3yiS27R4ZDzFyUJiZlUD/Q4w276r8hxg5KMzMSmAkTeXhoDAzK4FZk2sZXz2GZ7orv0PbQWFmVgJVY0RHcx0bR8BUHg4KM7MSOXvKJF96GipJZ0paK2lT9t8z8rTbJulJSY9JWlfuOs3MhuPsKfX0vHaYlw8cSbqUYUnqjGIZ8NOI6AB+ml3O5/KIWBARneUpzcysON7o0K7sfoqkgmIxcHv29e3A7yVUh5lZyfQHRaVP5ZFUUDRHRDdA9t+mPO0CuE/SI5KWFjqgpKWS1kla19PTU+RyzcxOXWPdeM4cAQ8xqi7VgSXdD0zJsemWUzjMJRGxU1ITsFbSMxHxYK6GEbECWAHQ2dk5MiaBN7OKJok5I2DkU8mCIiKuzLdN0kuSpkZEt6SpwK48x9iZ/XeXpFXAQiBnUJiZpdGc5npWPbqDiEBS0uUMSVKXnn4E3Jh9fSNw14kNJNVKqu9/DbwXWF+2Cs3MiqCjuZ7XDvfS/eqhpEsZsqSC4lbgKkmbgKuyy0iaJml1tk0z8K+SHgceBu6OiDWJVGtmNkRzmzMd2s9W8OWnkl16KiQi9gBX5Fi/E1iUfb0FeFuZSzMzK6o5zXVAJigum5tv3E66+c5sM7MSOn3iOBrrx/PsS5U7i6yDwsysxOY017Gpgi89OSjMzEpsTnM9z760n76+yhy576AwMyuxOc31HDx6jB17DyZdypA4KMzMSmxgh3YlclCYmZVYx/EhspXZoe2gMDMrsUk1Y5l6Wo3PKMzMLL+O5noHhZmZ5TenqY7Nu/ZzrAJHPjkozMzKYM6Ueg739vHCy68nXcopc1CYmZXBnAqe88lBYWZWBh1NmSGylXiHtoPCzKwMasdXM/30CRU5RNZBYWZWJnOnVObIJweFmVmZdDTXsaXnAL3H+pIu5ZQ4KMzMymROUz1HjvWxbU9ljXxyUJiZlUmljnxyUJiZlUl7Ux2Sg8LMzPKYMK6KmWdOZFOFjXxyUJiZlVFHU+WNfEokKCS9X9JTkvokdRZod7WkjZI2S1pWzhrNzEphTnMdW3cf4Ehv5Yx8SuqMYj1wHfBgvgaSqoDlwDXAfGCJpPnlKc/MrDTmTqmnty/YuvtA0qUMWiJBEREbImLjSZotBDZHxJaIOAKsBBaXvjozs9LpaKq8kU9p7qOYDmwfsNyVXZeTpKWS1kla19PTU/LizMyGYnZjLWMqbORTdakOLOl+YEqOTbdExF2DOUSOdXknco+IFcAKgM7Ozsqb8N3MRoWasVXMaqh1UABExJXDPEQXMGPAcguwc5jHNDNL3JwKG/mU5ktPvwY6JLVKGgfcAPwo4ZrMzIZtTnMd2/Yc4NDRY0mXMihJDY99n6Qu4CLgbkn3ZtdPk7QaICJ6gU8C9wIbgB9ExFNJ1GtmVkwdzfX0BTzXUxk33pXs0lMhEbEKWJVj/U5g0YDl1cDqMpZmZlZyHc2Zhxht3rWfc6adlnA1J5fmS09mZiNSa0Nm5NNzuyrjjMJBYWZWZuOrqzhrci2bHBRmZpZPe1Odg8LMzPLraKpj2+4DHK2Ap905KMzMEtDeVEdvX/D8nvTP+eSgMDNLQP+cT5XwbAoHhZlZAtqaaoHMENm0c1CYmSVg4rhqpp8+oSI6tB0UZmYJ6WiujJFPDgozs4R0NNWxpWc/x/rSPeG1g8LMLCHtTXUc7u2j65XXky6lIAeFmVlC2itk5JODwswsIe1N2ckBUz6LrIPCzCwhp00YS1P9eJ9RmJlZfh3NdWzele6n3TkozMwS1NFUz+Zd+4lI78gnB4WZWYLam+o4cOQY3a8eSrqUvBwUZmYJ6u/QTvONdw4KM7MEdTS98VjUtHJQmJklaHLdeM6sHZfqDu1EgkLS+yU9JalPUmeBdtskPSnpMUnrylmjmVm5tDfWpXqIbFJnFOuB64AHB9H28ohYEBF5A8XMrJK1ZycHTOvIp0SCIiI2RMTGJN7bzCxtOprqePXgUXbvP5J0KTmlvY8igPskPSJpaaGGkpZKWidpXU9PT5nKMzMbvjdGPqWzn6JkQSHpfknrc3wtPoXDXBIRFwDXAJ+QdGm+hhGxIiI6I6KzsbFx2PWbmZVL/2NRn0vpyKfqUh04Iq4swjF2Zv/dJWkVsJDB9WuYmVWM5knjqR9fndp7KVJ76UlSraT6/tfAe8l0gpuZjSiSaGuqS+29FEkNj32fpC7gIuBuSfdm10+TtDrbrBn4V0mPAw8Dd0fEmiTqNTMrtY6m9D4WtWSXngqJiFXAqhzrdwKLsq+3AG8rc2lmZonoaK7jh4908errRzlt4tiky3mT1F56MjMbTd54iFH6Rj6dNCgkNeVYN7c05ZiZjU4dKX4s6mDOKH4h6QP9C5L+BzkuG5mZ2dBNP30CNWPHpLKfYjB9FJcBKyS9n0wH8wYyw1TNzKxIxowRbY3p7NA+6RlFRHQDa8iMUJoF3BER6ftOzMwqXFtjHVt60vfrdTB9FGuBC4FzyYxI+j+SvlzqwszMRpv2pjp27D3IwSPHki7lTQbTR3EP8L8iYm9ErAcuBl4tbVlmZqNPW2MdEbBld7rOKgYTFPXAvZJ+IekTwOSI+FyJ6zIzG3XammoBeK7nQMKVvNlg+ij+IiLOAT4BTAP+RdL9Ja/MzGyUmTW5ljFK3+SAp3LD3S7gRWAP8JZ7K8zMbHhqxlYx48yJPJeyDu3BdGZ/TNLPgZ8CDcAfRMT5pS7MzGw0amtM3+SAg7mP4izgjyLisRLXYmY26rU11vLLzbs51hdUjVHS5QCD66NY5pAwMyuPtsY6Dvf2sXPvwaRLOc6TApqZpUjb8ckB03P5yUFhZpYi7Y2ZoEjTyCcHhZlZipxRO44za8elauSTg8LMLGXaGmt5bld6brpzUJiZpUxbY53PKMzMLL+2xjr2HDjCKweOJF0KkFBQSLpN0jOSnpC0StLpedpdLWmjpM2SlpW5TDOzRPQ/FjUtZxVJnVGsBc7N3uH9LPCnJzaQVAUsB64B5gNLJM0va5VmZgloa3RQEBH3RURvdvEhoCVHs4XA5ojYEhFHgJXA4nLVaGaWlOlnTGBc9ZjUzCKbhj6KD5N55sWJpgPbByx3ZdflJGmppHWS1vX09BS5RDOz8qkaI2Y31KbmXoqSBYWk+yWtz/G1eECbW4Be4Lu5DpFjXeR7v4hYERGdEdHZ2Ng4/G/AzCxBbU11qbk7ezCTAg5JRFxZaLukG4HfAa6IiFwB0AXMGLDcAuwsXoVmZunV1ljHPU92c+joMWrGViVaS1Kjnq4G/gS4NiJez9Ps10CHpFZJ44AbgB+Vq0YzsyS1NdbSF/D8nny/IssnqT6Kr5N5xOpaSY9J+jsASdMkrQbIdnZ/ErgX2AD8ICKeSqheM7OyStPIp5JdeiokItrzrN8JLBqwvBpYXa66zMzSYnZj5vnZaXiIURpGPZmZ2Qkmjqtm+ukTUnFG4aAwM0uptqZ0zPnkoDAzS6n+WWT7+vLeGVAWDgozs5Rqa6zj4NFjdO87lGgdDgozs5RqS8nT7hwUZmYplZZZZB0UZmYp1VA3jkk11Q4KMzPLTVJmzidfejIzs3wyj0VNdrpxB4WZWYq1N9XR89phXj14NLEaHBRmZimWhjmfHBRmZinWP+fT1gQvPzkozMxSbMYZE6kaI7bs9hmFmZnlMK56DDPPnMgWn1GYmVk+sxtqHRRmZpbf7MZatu45wLGEJgd0UJiZpdzsxjqO9Paxc+/BRN7fQWFmlnKzGzIjn7bsTubyk4PCzCzlWrNDZLckdC+Fg8LMLOUa68ZTP746sQ7t6iTeVNJtwO8CR4DngN+PiL052m0DXgOOAb0R0VnGMs3MUkESsxtrE7uXIqkzirXAuRFxPvAs8KcF2l4eEQscEmY2ms1urEvsjCKRoIiI+yKiN7v4ENCSRB1mZpVidkMt3a8e4vUjvSdvXGRp6KP4MHBPnm0B3CfpEUlLCx1E0lJJ6ySt6+npKXqRZmZJmp2dHHBrAiOfShYUku6XtD7H1+IBbW4BeoHv5jnMJRFxAXAN8AlJl+Z7v4hYERGdEdHZ2NhY1O/FzCxprf1DZBO4/FSyzuyIuLLQdkk3Ar8DXBEROW83jIid2X93SVoFLAQeLHatZmZpl2RQJHLpSdLVwJ8A10bE63na1Eqq738NvBdYX74qzczSY8K4KqafPiGRkU9J9VF8HagH1kp6TNLfAUiaJml1tk0z8K+SHgceBu6OiDXJlGtmlrzZjclMDpjIfRQR0Z5n/U5gUfb1FuBt5azLzCzNZjfU8o+P7iAikFS2903DqCczMxuE2Y117D/cS89rh8v6vg4KM7MK0d+h/VyZLz85KMzMKkT/87PL3aGdSB9FEo4ePUpXVxeHDh1KupSyq6mpoaWlhbFjxyZdipkNw7TTJlAzdkzZO7RHTVB0dXVRX1/PrFmzytoJlLSIYM+ePXR1ddHa2pp0OWY2DGPGiFmTa8s+3fioufR06NAhJk+ePKpCAjKzTk6ePHlUnkmZjURtjXVln8Zj1AQFMOpCot9o/b7NRqLZjbVsf+UgR3r7yvaeoyoozMwqXWtDLcf6ghdeLt9ZhYOiRC6++OKc62+66SbuvPPOMldjZiNF/yyy5Rwi66AokV/96ldJl2BmI9DxIbJlDIpRM+qp3Orq6ti/fz8Rwac+9SkeeOABWltbyTNRrpnZoEyqGUtD3Xi2lvFeCp9RlNiqVavYuHEjTz75JN/4xjd8pmFmw1buyQEdFCX24IMPsmTJEqqqqpg2bRrvec97ki7JzCpcW2MtW8o4RNZBUQYenmpmxdTaUMvLB46w9/UjZXk/B0WJXXrppaxcuZJjx47R3d3Nz372s6RLMrMKN7uhvCOf3JldYu973/t44IEHOO+885gzZw7vfve7ky7JzCrcGyOf9vOOs84o+fs5KEpk//7MiARJfP3rX0+4GjMbSWacOZHqMSrbVB6+9GRmVmHGVo1h5uSJZRv55KAwM6tAsxvqyvZcCgeFmVkFmt1Yy7Y9r3Osr/Q38SYSFJI+J+kJSY9Juk/StDztrpa0UdJmScvKXaeZWVrNbqjlSG8fO145WPL3SuqM4raIOD8iFgA/Af7sxAaSqoDlwDXAfGCJpPllrdLMLKWOTw5YhstPiQRFROwbsFgL5Dp3WghsjogtEXEEWAksLkd9ZmZp1z9EdmsZOrQT66OQ9HlJ24EPkuOMApgObB+w3JVdl+94SyWtk7Sup6enuMUWwfbt27n88suZN28e55xzDl/72tfe0iYi+MM//EPa29s5//zzefTRRxOo1MwqweTacUyqqS5Lh3bJgkLS/ZLW5/haDBARt0TEDOC7wCdzHSLHury9NhGxIiI6I6KzsbGxON9EEVVXV/OVr3yFDRs28NBDD7F8+XKefvrpN7W555572LRpE5s2bWLFihV87GMfS6haM0s7ScxurCvLENmS3XAXEVcOsun3gLuBz56wvguYMWC5BdhZhNL4ix8/xdM795284SmYP20Sn/3dc/Junzp1KlOnTgWgvr6eefPmsWPHDubPf6Pb5a677uJDH/oQknjnO9/J3r176e7uPr6fmdlAsxtq+dVze0r+PkmNeuoYsHgt8EyOZr8GOiS1ShoH3AD8qBz1ldq2bdv4zW9+w4UXXvim9Tt27GDGjDeysaWlhR07dpS7PDOrELMba3lx3yEOHO4t6fskNYXHrZLmAn3A88DNANlhst+MiEUR0Svpk8C9QBXw9xHxVDHevNBf/qW2f/9+rr/+er761a8yadKkN23L9VAjzzxrZvn0j3zauvsA504/rWTvk0hQRMT1edbvBBYNWF4NrC5XXaV29OhRrr/+ej74wQ9y3XXXvWV7S0sL27e/0X/f1dXFtGk5bzExM6O1ITvyqcRB4TuzyyQi+MhHPsK8efP49Kc/nbPNtddeyx133EFE8NBDD3Haaae5f8LM8po1ORMU20o8OaBnjy2TX/7yl3znO9/hvPPOY8GCBQB84Qtf4IUXXgDg5ptvZtGiRaxevZr29nYmTpzIt7/97QQrNrO0mzCuimmn1ZR8FlkHRZm8613vytkHMZAkli9fXqaKzGwkaC3DY1F96cnMrILNmlzLlp79J/1DdDgcFGZmFay1oZZ9h3p55fWjJXsPB4WZWQU7PudTCafycFCYmVWw1obMvRSlnMrDQWFmVsFazphA9RixbY+DwszMchhbNYaZZ04s6RBZB0UZffjDH6apqYlzzz33+LqXX36Zq666io6ODq666ipeeeWV49u++MUv0t7ezty5c7n33ntzHrPQ/mY2OsxqqPWlp5HipptuYs2aNW9ad+utt3LFFVewadMmrrjiCm699VYAnn76aVauXMlTTz3FmjVr+PjHP86xY8fecsx8+5vZ6NHaUMu2PQfoK9Hzs0fnDXf3LIMXnyzuMaecB9cU/iV96aWXsm3btjetu+uuu/j5z38OwI033shll13Gl770Je666y5uuOEGxo8fT2trK+3t7Tz88MNcdNFFg9rfzEaP1oZaDh3t48V9h5h2+oSiH99nFAl76aWXjs/nNHXqVHbt2gUMfsrxfPub2egxe8DkgKUwOs8oTvKXfxp4ynEzG6zWxjeC4pL2hqIf32cUCWtubqa7uxuA7u5umpqagMFPOZ5vfzMbPZrra5gwtqpkZxQOioRde+213H777QDcfvvtLF68+Pj6lStXcvjwYbZu3cqmTZtYuHDhoPc3s9FjzBgxq6HWQTESLFmyhIsuuoiNGzfS0tLCt771LZYtW8batWvp6Ohg7dq1LFu2DIBzzjmHD3zgA8yfP5+rr76a5cuXU1VVBcBHP/pR1q1bB5B3fzMbXVobSncvhUo542BSOjs7o/8Xab8NGzYwb968hCpK3mj//s1Guu8//AKPb9/LF687b0j9mZIeiYjOXNtGZ2e2mdkIs2ThTJYsnFmSYycSFJI+BywG+oBdwE3Z52Wf2G4b8BpwDOjNl3ZmZlY6SfVR3BYR50fEAuAnwJ8VaHt5RCwoRkiMxMtsgzFav28zK45EgiIi9g1YrAVK/puspqaGPXv2jLpfmhHBnj17qKmpSboUM6tQifVRSPo88CHgVeDyPM0CuE9SAP83IlYM9f1aWlro6uqip6dnqIeoWDU1NbS0tCRdhplVqJKNepJ0PzAlx6ZbIuKuAe3+FKiJiM/mOMa0iNgpqQlYC3wqIh7M835LgaUAM2fOfMfzzz9fjG/DzGxUKDTqKfHhsZLOAu6OiHNP0u7Pgf0R8eWTHTPX8FgzM8uvUFAk0kchqWPA4rXAMzna1Eqq738NvBdYX54KzcysX1J9FLdKmktmeOzzwM2QudQEfDMiFgHNwKrsjSPVwPciYk2e45mZWYkkfumpFCT1kAmgoWgAdhexnGJLe32Q/hpd3/Clvca01wfpq/GsiGjMtWFEBsVwSFqX5hv70l4fpL9G1zd8aa8x7fVBZdTYz5MCmplZQQ4KMzMryEHxVkO+qa9M0l4fpL9G1zd8aa8x7fVBZdQIuI/CzMxOwmcUZmZWkIPCzMwKGjVBIelqSRslbZb0lueFKuOvstufkHTBYPdNukZJMyT9TNIGSU9J+m9pqm/A9ipJv5H0k1LUN9waJZ0u6U5Jz2R/lhelrL7/nv3vu17S9yUVfUrgQdR3tqR/k3RY0mdOZd+ka0zR5yTvzzC7veSfk1MWESP+C6gCngNmA+OAx4H5J7RZBNwDCHgn8O+D3TcFNU4FLsi+rgeeLXaNw6lvwPZPA98DfpK2/87ZbbcDH82+Hgecnpb6gOnAVmBCdvkHZB74Ve76moDfAj4PfOZU9k1BjWn5nOSsr1yfk6F8jZYzioXA5ojYEhFHgJVknrA30GLgjsh4CDhd0tRB7ptojRHRHRGPAkTEa8AGMr9YUlEfgKQW4LeBbxa5rqLUKGkScCnwLYCIOBIRe9NSX3ZbNTBBUjUwEXjLUyFLXV9E7IqIXwNHT3XfpGtMy+ekwM+wXJ+TUzZagmI6sH3Achdv/R8kX5vB7Jt0jcdJmgW8Hfj3lNX3VeB/kpnfq1SGU+NsoAf4dva0/5vKTEaZivoiYgfwZeAFoBt4NSLuS6C+Uux7KoryPgl/Tgr5KqX/nJyy0RIUyrHuxHHB+doMZt9iGE6NmY1SHfCPwB/Fm58iWAxDrk/S7wC7IuKRItd0ouH8DKuBC4C/jYi3AweAYl9nH87P8Awyf5m2AtOAWkn/JYH6SrHvqRj2+6Tgc5J7x/J9Tk7ZaAmKLmDGgOUW3nranq/NYPZNukYkjSXzP/93I+KfUlbfJcC1kraRORV/j6R/SFmNXUBXRPT/hXknmeBIS31XAlsjoicijgL/BFycQH2l2PdUDOt9UvI5yadcn5NTl3QnSTm+yPy1uIXMX2P9HUznnNDmt3lzJ+LDg903BTUKuAP4ahp/hie0uYzSdWYPq0bgF8Dc7Os/B25LS33AhcBTZPomRKbj/VPlrm9A2z/nzR3FqfmcFKgxFZ+TfPWdsK1kn5MhfV9JF1C2bzQzmuRZMiMSbsmuuxm4ecD/RMuz258EOgvtm6YagXeROb19Angs+7UoLfWdcIySfgCG+d95AbAu+3P8Z+CMlNX3F2Qe8rUe+A4wPoH6ppD5q3kfsDf7elLKPic5a0zR5yTvz7Bcn5NT/fIUHmZmVtBo6aMwM7MhclCYmVlBDgozMyvIQWFmZgU5KMzMrCAHhVkB2RllPz5geZqkO0v0Xr8n6c9O0ubLkt5Tivc3y8fDY80KyM4J9JOIOLcM7/Ur4NqI2F2gzVnANyLivaWux6yfzyjMCrsVaJP0mKTbJM2StB5A0k2S/lnSjyVtlfRJSZ/OTir4kKQzs+3aJK2R9IikX0g6+8Q3kTQHOBwRuyXVZ483NrttkqRtksZGxPPAZElTyvgzsFHOQWFW2DLguYhYEBF/nGP7ucB/JjO99OeB1yMzqeC/AR/KtllBZrqNdwCfAf4mx3EuAQZOgf1zMtN5ANwA/GNk5ngi2+6SYX5fZoNWnXQBZhXuZ9lf7K9JehX4cXb9k8D52ZlKLwZ+KB2fWHR8juNMJTPNeb9vkplu+p+B3wf+YMC2XWRmkDUrCweF2fAcHvC6b8ByH5nP1xhgb0QsOMlxDgKn9S9ExC+zl7neDVRFxPoBbWuy7c3KwpeezAp7jcxjM4ckMs872Crp/XD8mdhvy9F0A9B+wro7gO8D3z5h/RwyEwOalYWDwqyAiNgD/FLSekm3DfEwHwQ+IulxMlOF53pE6IPA2zXg+hTwXeAMMmEBHH+eQjuZWW7NysLDY81SQtLXgB9HxP3Z5f8ILI6I/zqgzfuACyLifydUpo1C7qMwS48vkHlAEZL+GriGzLMNBqoGvlLmumyU8xmFmZkV5D4KMzMryEFhZmYFOSjMzKwgB4WZmRXkoDAzs4L+P9l9DczSV9flAAAAAElFTkSuQmCC\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAERCAYAAABy/XBZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAXEklEQVR4nO3df5TddX3n8ef7JoGUEkChgcQhJhp+JBAadeSXFpAfrdJuWGFlYW0lK8oJWKhl2W62nG2Le0QsehZ2TduN/FjisuYIAlk9ECoFCotGGhAEEjEqESZE8gOjppQfkvf+ce8Mk2HmZiZz73zvfL/Pxzlzztzvr/uaO3Pvez7fz/f7+URmIkmqnlrRASRJxbAASFJFWQAkqaIsAJJUURYASaooC4AkVdS4KwARcUNEbIqIJ1t0vJURsS0ivjlgeUTEZyPihxGxNiIuacXzSVKnGHcFAPhfwAdbeLyrgT8aZPlC4GDg8MycAyxv4XNKUuHGXQHIzAeAF/svi4h3Nv6TfyQiHoyIw0dwvH8AfjXIqguBz2TmjsZ2m0aTW5I6zbgrAENYClycme8BLgP+pgXHfCfwbyNidUTcFRGHtOCYktQxJhYdYLQiYm/geOCWiOhdvGdj3ZnAZwbZbUNm/t4uDr0n8HJmdjeOcwPwO61JLUnFG/cFgHorZltmzh+4IjNvA27bzeP2AF9vfH87cONuHkeSOtK4PwWUmb8EnomIj0Df1Tu/3YJD3wGc3Pj+ROCHLTimJHWMGG+jgUbEV4GTgAOAF4C/BO4F/haYBkwClmfmYKd+Bjveg8DhwN7AVuD8zLw7IvYDbgZmANuBRZn5eEt/GEkq0LgrAJKk1hj3p4AkSbtnXHUCH3DAATlz5syiY0jSuPLII49syczfGrh8XBWAmTNnsnr16qJjSNK4EhE/HWy5p4AkqaIsAJJUURYASaooC4AkVZQFQJIqqtACEBEfjIinI+JHEbG4yCySVDWFFYCImAAsAT4EzAXOjYi5ReWRpKop8j6Ao4EfZeZPACJiOXAGsKbVT/TszZcQLzzR6sNKqoCX3jqXwxYuKTpGWxRZAN4GPNfvcQ9wzMCNIuIC4AKAGTNm7NYTbdj2ErHtX3ZrX0nVtubFLRz40qvst9ceRUdpuSILQAyy7E0j02XmUuozftHd3b1bI9cd96nrdmc3SRW37Dvr+cyKp1iwo5yDZhbZCdxDfdL1Xl3A8wVlkaQ36Z1lcEdJR00usgD8E3BIRMyKiD2Ac4D/W2AeSdpJrfc8RTk//4s7BZSZv46IPwbuBiYAN2TmU0XlkaSBan0tgIKDtEmho4Fm5p3AnUVmkKSh9DYAPAUkSRVTsw9Akqqp8flPST//LQCSNJTeFoAFQJIqptb4hPQUkCRVTGAfgCRVUm8fQFkvA7UASNIQ3ugDKGcFsABI0hD6CkDBOdrFAiBJQ6j1nQIqZwmwAEjSEPr6AHYUm6NdLACSNARHA5WkivJGMEmqqN4+gCxpN7AFQJKGUPbhoC0AkjQUrwKSpGryRjBJqqiaQ0FIUjV5FZAkVVTYByBJ1eRw0JJUUTWnhJSkaqrV7AOQpEpyNFBJqiz7ACSpkuwDkKSKqjkctCRVkzeCSVJFeSOYJFVUOBaQJFWTo4FKUkU5IYwkVZRTQkpSRYUtAEmqpui7EaycFcACIElD8EYwSaqovsHgdhSbo10KKQARcXVE/CAivh8Rt0fEfkXkkKRm+i4DLThHuxTVAvgWcGRmHgX8EPjPBeWQpCF5J3AbZObfZ+avGw9XAV1F5JCkZsIbwdru48BdQ62MiAsiYnVErN68efMYxpJUdbWSDwUxsV0Hjoh7gIMGWXV5Zq5obHM58Gvg5qGOk5lLgaUA3d3dJf01SOpEZb8KqG0FIDNPbbY+Is4D/gA4JcvavpI0rkXJJ4RpWwFoJiI+CPwn4MTMfKmIDJK0Kw4G1x5fAqYA34qIxyLi7wrKIUlDajQA7ANopcycXcTzStJIlL0PoBOuApKkjuSUkJJUUdH4hLQFIEkVYwtAkirqjU7gclYAC4AkDcEpISWposIpISWpmuwDkKSK6hsOuqTngCwAkjQE+wAkqaJqTggjSdUUTgkpSdVVC0cDlaRKighPAUlSFdXCTmBJqiRbAJJUUbWgtL3AFgBJaqJmC0CSqimwD0CSKskWgCRVVISDwUlSJdVq4Y1gklRF9VNARadoDwuAJDVR7wQuZwWwAEhSE2ELQJKqycHgJKmiahFeBSRJVVQfDK6cFcACIElN2AcgSRUV9gFIUjXVIso6GKgFQJKasQ9AkirKO4ElqapsAUhSNdXvA7AASFLl1BwOuj0i4rKIyIg4oMgckjQUJ4Rpg4g4GDgNeLaoDJI0HHYCt95/A/4MSnuJraQSsA+gxSJiAbAhMx8fxrYXRMTqiFi9efPmMUgnSW+o1crbApjYrgNHxD3AQYOsuhz4c+B3h3OczFwKLAXo7u4u6a9BUqcqcwugbQUgM08dbHlEzANmAY9HBEAX8GhEHJ2ZP2tXHknaHWUeDK5tBWAomfkEMLX3cUSsB7ozc8tYZ5GkXXFKSEmqqDLfBzDmLYCBMnNm0RkkaSjeByBJFVXpKSEjYuogyw5rTxxJ6ixR8cHgHoyIs3sfRMR/AG5vXyRJ6hxR8T6Ak4ClEfER4EBgLXB0O0NJUqeoRfDr3FF0jLbYZQsgMzcCK4HjgJnAsszc3uZcktQRytwJvMsWQER8C9gIHEn9pq0bIuKBzLys3eEkqWgR5R2wbDh9AHcBf56Z2zLzSeB44BftjSVJnaHqU0JOAe6OiAcj4lPA/pn5X9ucS5I6Qr0TuJwVYDh9AFdk5hHAp4DpwD82BnqTpNIrcx/ASG4E2wT8DNhKv7F8JKnMyjwUxHBuBLswIu4H/gE4APhkZh7V7mCS1AmqPhro24FPZ+Zjbc4iSR2nVuI+gF0WgMxcPBZBJKkTBfYBSFIllXlKSAuAJDURJZ4S0gIgSU1UejhoSaoyp4SUpIqqhX0AklRJ3gksSRUV9gFIUjWV+UYwC4AkNRH2AUhSNdkHIEkVVebB4CwAktRELaCsk0JaACSpiapPCSlJlVXvBC5nBbAASFITtQh2lLQJYAGQpCaiylNCSlKV1SJK2gVsAZCkpmr2AUhSNYU3gklSNTkUhCRVVM0pISWpmmpeBSRJ1eRgcG0QERdHxNMR8VRE/HVROSSpmfqcwEWnaI+JRTxpRHwAOAM4KjNfiYipReSQpF2JCKA+KUzv92VRVAvgQuCqzHwFIDM3FZRDkpqq9RWAgoO0QVEF4FDgdyLiuxHxjxHx3qE2jIgLImJ1RKzevHnzGEaUpN7hoMt5M1jbTgFFxD3AQYOsurzxvG8BjgXeC3wtIt6Rg1xrlZlLgaUA3d3d5fsNSOpotUYFKGM/QNsKQGaeOtS6iLgQuK3xgf9wROwADgD8F19SRypjC6CoU0B3ACcDRMShwB7AloKySNKQytwHUMhVQMANwA0R8STwKnDeYKd/JKlovX0AWcIxQQspAJn5KvCHRTy3JI1EbwugjH0A3gksSU1Eia8CsgBIUhN9N4LtKDhIG1gAJKmJMt8HYAGQpCb6rgIqOEc7WAAkqQlbAJJUVX1XAVkAJKlS+u4DKN/nvwVAkpqp2QKQpGqyBSBJFRW2ACSpmnrnACvh578FQJKasQ9Akiqq1viUdDA4SaqYWr9J4cvGAiBJTUSJh4MuakKYlnnttdfo6enh5ZdfLjpKISZPnkxXVxeTJk0qOopUSm90ApevAoz7AtDT08OUKVOYOXNmX6Wuisxk69at9PT0MGvWrKLjSKXkhDAd7OWXX2b//fev3Ic/1Jum+++/f2VbP9JYKPOUkOO+AACV/PDvVeWfXRoLfX0ATggjSdXicNDayfHHHz/o8oULF3LrrbeOcRpJ7dQ3JWT5Pv8tALvj29/+dtERJI2RMrcAxv1VQEXYe++92b59O5nJxRdfzL333susWbNKeZmYVHVOCalB3X777Tz99NM88cQTfPnLX7ZlIJVQlLgFYAEYhQceeIBzzz2XCRMmMH36dE4++eSiI0lqsXAoCA3FyzClcnujD6DYHO1gARiFE044geXLl/P666+zceNG7rvvvqIjSWqxvjuBS1gB7AQehQ9/+MPce++9zJs3j0MPPZQTTzyx6EiSWiz67gQuHwvAbti+fTtQP/3zpS99qeA0ktrJCWEkqaKcElKSKqpWswUgSZXkVUCSVFHeByBJFVVzMDhJqqbeTmD7AFokIuZHxKqIeCwiVkfE0UXkaJXnnnuOD3zgA8yZM4cjjjiCa6+99k3bZCaXXHIJs2fP5qijjuLRRx8tIKmkkXJKyNb7a+CKzJwP/EXj8bg1ceJEvvjFL7J27VpWrVrFkiVLWLNmzU7b3HXXXaxbt45169axdOlSLrzwwoLSShqJMg8GV9SNYAns0/h+X+D5Vhz0im88xZrnf9mKQ/WZO30f/vJfHdF0m2nTpjFt2jQApkyZwpw5c9iwYQNz587t22bFihV87GMfIyI49thj2bZtGxs3buzbT1JnKnMfQFEF4NPA3RHxBeqtkMGn2AIi4gLgAoAZM2aMSbjRWL9+Pd/73vc45phjdlq+YcMGDj744L7HXV1dbNiwwQIgdbha4zxJGa8CalsBiIh7gIMGWXU5cArwp5n59Yg4G7geOHWw42TmUmApQHd3d9PfwK7+U2+37du3c9ZZZ3HNNdewzz777LRusD8eRxKVOl9Q3j6AthWAzBz0Ax0gIpYBf9J4eAtwXbtyjJXXXnuNs846i49+9KOceeaZb1rf1dXFc8891/e4p6eH6dOnj2VESbuhzFNCFtUJ/DzQO3TmycC6gnK0RGZy/vnnM2fOHC699NJBt1mwYAHLli0jM1m1ahX77ruvp3+kcSBKPBhcUX0AnwSujYiJwMs0zvGPVw899BBf+cpXmDdvHvPnzwfgyiuv5NlnnwVg0aJFnH766dx5553Mnj2bvfbaixtvvLHAxJKGq1biM7WFFIDM/H/Ae4p47nZ4//vfv8sOoohgyZIlY5RIUqs4HLQkVVTffQA7is3RDhYASWrCFoAkVVSZp4S0AEhSEzWHg5akanIwOEmqqDIPBmcBaIGPf/zjTJ06lSOPPLJv2Ysvvshpp53GIYccwmmnncbPf/7zvnWf+9znmD17Nocddhh33333oMdstr+ksRNOCalmFi5cyMqVK3dadtVVV3HKKaewbt06TjnlFK666ioA1qxZw/Lly3nqqadYuXIlF110Ea+//vqbjjnU/pLGVq2vF7h8FaCoO4Hb467F8LMnWnvMg+bBh5p/+J5wwgmsX79+p2UrVqzg/vvvB+C8887jpJNO4vOf/zwrVqzgnHPOYc8992TWrFnMnj2bhx9+mOOOO25Y+0saW/YBaMReeOGFvrF+pk2bxqZNm4Chh4Ue7v6SxlaZp4QsVwtgF/+pdwKHhZbGF1sAGrEDDzyQjRs3ArBx40amTp0KDH9Y6KH2lzS2osQTwlgA2mTBggXcdNNNANx0002cccYZfcuXL1/OK6+8wjPPPMO6des4+uijh72/pLHllJBq6txzz+X+++9ny5YtdHV1ccUVV7B48WLOPvtsrr/+embMmMEtt9wCwBFHHMHZZ5/N3LlzmThxIkuWLGHChAkAfOITn2DRokV0d3cPub+ksdU7HPT/fODHfG31c803bqMrz5zHe2e+taXHjPHUrOnu7s7Vq1fvtGzt2rXMmTOnoESdwddAap/M5PMrn+bZF/+50BwXnTSbI9+2727tGxGPZGb3wOW2ACSpiYhg8YcOLzpGW9gHIEkVVYoCMJ5OY7ValX92SaMz7gvA5MmT2bp1ayU/CDOTrVu3Mnny5KKjSBqHxn0fQFdXFz09PWzevLnoKIWYPHkyXV1dRceQNA6N+wIwadIkZs2aVXQMSRp3xv0pIEnS7rEASFJFWQAkqaLG1Z3AEbEZ+Olu7n4AsKWFcdqh0zOab/Q6PWOn54POz9iJ+d6emb81cOG4KgCjERGrB7sVupN0ekbzjV6nZ+z0fND5GTs9X3+eApKkirIASFJFVakALC06wDB0ekbzjV6nZ+z0fND5GTs9X5/K9AFIknZWpRaAJKkfC4AkVVQpCkBEfDAino6IH0XE4kHWR0T898b670fEu4e7b5H5IuLgiLgvItZGxFMR8SedlK/f+gkR8b2I+GY78o02Y0TsFxG3RsQPGq/lcR2W708bv98nI+KrEdGW4V2HkfHwiPhORLwSEZeNZN8i843V+2Q0Gfutb/t7ZUQyc1x/AROAHwPvAPYAHgfmDtjmdOAuIIBjge8Od9+C800D3t34fgrww07K12/9pcD/Ab7Zab/jxrqbgE80vt8D2K9T8gFvA54BfqPx+GvAwoJew6nAe4HPApeNZN+C87X9fTLajGP1XhnpVxlaAEcDP8rMn2Tmq8By4IwB25wBLMu6VcB+ETFtmPsWli8zN2bmowCZ+StgLfUPjI7IBxARXcDvA9e1OFdLMkbEPsAJwPUAmflqZm7rlHyNdROB34iIicBewPMtzjesjJm5KTP/CXhtpPsWmW+M3iejyghj9l4ZkTIUgLcBz/V73MObf/lDbTOcfYvM1yciZgLvAr7bYfmuAf4M2NHiXMN9/l1t8w5gM3Bjo+l9XUT8Zqfky8wNwBeAZ4GNwC8y8+9bnG+4Gdux73C15Dna+D6B0We8hva/V0akDAUgBlk28NrWobYZzr6jNZp89ZURewNfBz6dmb9sYbZdPnezbSLiD4BNmflIizMNNJrXcCLwbuBvM/NdwD8DrT6HPZrX8C3U/4ucBUwHfjMi/rDF+YZ8/jHYd7hG/Rxtfp/AKDKO4XtlRMpQAHqAg/s97uLNTeihthnOvkXmIyImUf+jvjkzb2txttHmex+wICLWU28OnxwR/7vDMvYAPZnZ+x/hrdQLQqfkOxV4JjM3Z+ZrwG3A8S3ON9yM7dh3uEb1HGPwPoHRZRyr98rIFN0JMdov6v/h/YT6f1C9HTNHDNjm99m5A+7h4e5bcL4AlgHXdOLrN2Cbk2hfJ/CoMgIPAoc1vv8r4OpOyQccAzxF/dx/UO+wvriI17Dftn/Fzp2sHfE+aZKv7e+T0WYcsK5t75UR/0xFB2jRL+Z06j3/PwYubyxbBCzq9weypLH+CaC72b6dkg94P/Um5veBxxpfp3dKvgHHaOsf9Sh/x/OB1Y3X8Q7gLR2W7wrgB8CTwFeAPQt6DQ+i/l/uL4Ftje/36aD3yaD5xup9MtrXcKzeKyP5cigISaqoMvQBSJJ2gwVAkirKAiBJFWUBkKSKsgBIUkVZAFRZjVFCL+r3eHpE3Nqm5/rXEfEXu9jmCxFxcjueXxqMl4GqshrjxnwzM48cg+f6NrAgM7c02ebtwJcz83fbnUcCWwCqtquAd0bEYxFxdUTMjIgnASJiYUTcERHfiIhnIuKPI+LSxoByqyLirY3t3hkRKyPikYh4MCIOH/gkEXEo8EpmbomIKY3jTWqs2yci1kfEpMz8KbB/RBw0hq+BKswCoCpbDPw4M+dn5n8cZP2RwL+jPgzwZ4GXsj6g3HeAjzW2WUp96Ib3AJcBfzPIcd4H9B+u+H7qQ0MAnAN8PevjANHY7n2j/LmkYZlYdACpg93X+MD+VUT8AvhGY/kTwFGN0SePB26J6Bsocs9BjjON+pDUva6jPizwHcC/Bz7Zb90m6qOCSm1nAZCG9kq/73f0e7yD+nunBmzLzPm7OM6/APv2PsjMhxqnm04EJmTmk/22ndzYXmo7TwGpyn5FfQrB3ZL1MeefiYiPQN+8v789yKZrgdkDli0DvgrcOGD5odQHhZPazgKgysrMrcBDjcnYr97Nw3wUOD8iHqc+rPNgUyU+ALwr+p0nAm4G3kK9CAB9Y9rPpj5yqdR2XgYqjYGIuBb4Rmbe03j8b4AzMvOP+m3zYeqTm/+XgmKqYuwDkMbGldQnfyEi/gfwIepjy/c3EfjiGOdShdkCkKSKsg9AkirKAiBJFWUBkKSKsgBIUkVZACSpov4/+QGZUlrJI0YAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] @@ -466,134 +466,134 @@ " fill: currentColor;\n", "}\n", "
<xarray.DataArray 'vx' (time (y): 221)>\n",
-       "array([ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "array([ 0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
        "...\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
-       "        0.00000000e+00,  3.98117095e-11, -9.46069889e-11, -2.30539143e-10,\n",
-       "       -3.67965214e-10, -5.06869213e-10, -6.47232490e-10, -7.89038168e-10,\n",
-       "       -9.32266708e-10, -1.07690123e-09, -1.22292310e-09, -1.37031364e-09,\n",
-       "       -1.51905422e-09, -1.66912617e-09, -1.82051085e-09, -1.97318961e-09,\n",
-       "       -2.12714202e-09, -2.28234942e-09, -2.43879317e-09, -2.59645283e-09,\n",
-       "       -2.75530976e-09, -2.91534441e-09, -3.07653725e-09, -3.23886606e-09,\n",
-       "       -3.40231310e-09])\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,\n",
+       "        0.0000000e+00, -8.8817842e-16, -8.8817842e-16, -8.8817842e-16,\n",
+       "       -8.8817842e-16, -8.8817842e-16, -8.8817842e-16, -8.8817842e-16,\n",
+       "       -8.8817842e-16, -8.8817842e-16, -8.8817842e-16, -8.8817842e-16,\n",
+       "       -8.8817842e-16, -8.8817842e-16, -8.8817842e-16, -8.8817842e-16,\n",
+       "       -8.8817842e-16, -8.8817842e-16, -8.8817842e-16, -8.8817842e-16,\n",
+       "       -8.8817842e-16, -8.8817842e-16, -8.8817842e-16, -8.8817842e-16,\n",
+       "       -8.8817842e-16])\n",
        "Coordinates:\n",
        "    id        float64 2.0\n",
-       "  * time (y)  (time (y)) float64 0.0 0.0006845 0.001369 ... 0.1492 0.1499 0.1506
    • id
      ()
      float64
      2.0
      array(2.)
    • time (y)
      (time (y))
      float64
      0.0 0.0006845 ... 0.1499 0.1506
      array([0.      , 0.000684, 0.001369, ..., 0.149213, 0.149897, 0.150582])
    • " ], "text/plain": [ "\n", - "array([ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + "array([ 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", "...\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", - " 0.00000000e+00, 3.98117095e-11, -9.46069889e-11, -2.30539143e-10,\n", - " -3.67965214e-10, -5.06869213e-10, -6.47232490e-10, -7.89038168e-10,\n", - " -9.32266708e-10, -1.07690123e-09, -1.22292310e-09, -1.37031364e-09,\n", - " -1.51905422e-09, -1.66912617e-09, -1.82051085e-09, -1.97318961e-09,\n", - " -2.12714202e-09, -2.28234942e-09, -2.43879317e-09, -2.59645283e-09,\n", - " -2.75530976e-09, -2.91534441e-09, -3.07653725e-09, -3.23886606e-09,\n", - " -3.40231310e-09])\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,\n", + " 0.0000000e+00, -8.8817842e-16, -8.8817842e-16, -8.8817842e-16,\n", + " -8.8817842e-16, -8.8817842e-16, -8.8817842e-16, -8.8817842e-16,\n", + " -8.8817842e-16, -8.8817842e-16, -8.8817842e-16, -8.8817842e-16,\n", + " -8.8817842e-16, -8.8817842e-16, -8.8817842e-16, -8.8817842e-16,\n", + " -8.8817842e-16, -8.8817842e-16, -8.8817842e-16, -8.8817842e-16,\n", + " -8.8817842e-16, -8.8817842e-16, -8.8817842e-16, -8.8817842e-16,\n", + " -8.8817842e-16])\n", "Coordinates:\n", " id float64 2.0\n", " * time (y) (time (y)) float64 0.0 0.0006845 0.001369 ... 0.1492 0.1499 0.1506" @@ -987,7 +987,7 @@ " nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan])\n", "Coordinates:\n", " id float64 100.0\n", - " * time (y) (time (y)) float64 0.0 0.0006845 0.001369 ... 0.1492 0.1499 0.1506
    • id
      ()
      float64
      100.0
      array(100.)
    • time (y)
      (time (y))
      float64
      0.0 0.0006845 ... 0.1499 0.1506
      array([0.      , 0.000684, 0.001369, ..., 0.149213, 0.149897, 0.150582])
    • " ], "text/plain": [ "\n", @@ -1441,7 +1441,7 @@ " nan])\n", "Coordinates:\n", " id int64 100\n", - " * time (time) float64 0.0 0.0006845 0.001369 ... 0.1492 0.1499 0.1506
    • id
      ()
      int64
      100
      array(100)
    • time
      (time)
      float64
      0.0 0.0006845 ... 0.1499 0.1506
      array([0.      , 0.000684, 0.001369, ..., 0.149213, 0.149897, 0.150582])
    • " ], "text/plain": [ "\n", diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index bebf2acfa..ee2521916 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -317,6 +317,8 @@ module swiftest_classes integer(I4B), dimension(:), allocatable :: status !! status of the interaction integer(I4B), dimension(:), allocatable :: index1 !! position of the first body in the encounter integer(I4B), dimension(:), allocatable :: index2 !! position of the second body in the encounter + integer(I4B), dimension(:), allocatable :: id1 !! id of the first body in the encounter + integer(I4B), dimension(:), allocatable :: id2 !! id of the second body in the encounter real(DP), dimension(:,:), allocatable :: x1 !! the position of body 1 in the encounter real(DP), dimension(:,:), allocatable :: x2 !! the position of body 2 in the encounter real(DP), dimension(:,:), allocatable :: v1 !! the velocity of body 1 in the encounter diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 8516e3861..fb42b7434 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -711,7 +711,7 @@ end subroutine symba_util_spill_pl module subroutine symba_util_spill_encounter(self, discards, lspill_list, ldestructive) use swiftest_classes, only : swiftest_encounter implicit none - class(symba_encounter), intent(inout) :: self !! SyMBA pl-tp encounter list + class(symba_encounter), intent(inout) :: self !! SyMBA pl-tp encounter list class(swiftest_encounter), intent(inout) :: discards !! Discarded object logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 726da07fc..ea15bf1fe 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -88,6 +88,8 @@ module subroutine setup_encounter(self, n) if (allocated(self%status)) deallocate(self%status) if (allocated(self%index1)) deallocate(self%index1) if (allocated(self%index2)) deallocate(self%index2) + if (allocated(self%id1)) deallocate(self%id1) + if (allocated(self%id2)) deallocate(self%id2) if (allocated(self%x1)) deallocate(self%x1) if (allocated(self%x2)) deallocate(self%x2) if (allocated(self%v1)) deallocate(self%v1) @@ -98,6 +100,8 @@ module subroutine setup_encounter(self, n) allocate(self%status(n)) allocate(self%index1(n)) allocate(self%index2(n)) + allocate(self%id1(n)) + allocate(self%id2(n)) allocate(self%x1(NDIM,n)) allocate(self%x2(NDIM,n)) allocate(self%v1(NDIM,n)) @@ -108,6 +112,8 @@ module subroutine setup_encounter(self, n) self%status(:) = INACTIVE self%index1(:) = 0 self%index2(:) = 0 + self%id1(:) = 0 + self%id2(:) = 0 self%x1(:,:) = 0.0_DP self%x2(:,:) = 0.0_DP self%v1(:,:) = 0.0_DP diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index 326f5d257..eb230b7e0 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -43,6 +43,8 @@ module function symba_encounter_check_pl(self, system, dt, irec) result(lany_enc plplenc_list%lvdotr(1:nenc) = pack(loc_lvdotr(:), lencounter(:)) plplenc_list%index1(1:nenc) = pack(pl%k_plpl(1,:), lencounter(:)) plplenc_list%index2(1:nenc) = pack(pl%k_plpl(2,:), lencounter(:)) + plplenc_list%id1(1:nenc) = pl%id(plplenc_list%index1(1:nenc)) + plplenc_list%id2(1:nenc) = pl%id(plplenc_list%index2(1:nenc)) do k = 1, nenc plplenc_list%status(k) = ACTIVE plplenc_list%level(k) = irec @@ -178,6 +180,8 @@ module function symba_encounter_check_tp(self, system, dt, irec) result(lany_enc pltpenc_list%lvdotr(1:nenc) = pack(loc_lvdotr(:,:), lencounter(:,:)) pltpenc_list%index1(1:nenc) = pack(spread([(i, i = 1, npl)], dim=1, ncopies=ntp), lencounter(:,:)) pltpenc_list%index2(1:nenc) = pack(spread([(i, i = 1, ntp)], dim=2, ncopies=npl), lencounter(:,:)) + pltpenc_list%id1(1:nenc) = pl%id(pltpenc_list%index1(1:nenc)) + pltpenc_list%id2(1:nenc) = tp%id(pltpenc_list%index2(1:nenc)) select type(pl) class is (symba_pl) pl%lencounter(:) = .false. diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 8340f6e14..2232b9599 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -395,12 +395,13 @@ module subroutine symba_util_rearray_pl(self, system, param) class(symba_parameters), intent(in) :: param !! Current run configuration parameters ! Internals class(symba_pl), allocatable :: tmp !! The discarded body list. - integer(I4B) :: i + integer(I4B) :: i, j, k logical, dimension(:), allocatable :: lmask class(symba_plplenc), allocatable :: plplenc_old logical :: lencounter associate(pl => self, pl_adds => system%pl_adds) + allocate(tmp, mold=pl) ! Remove the discards and destroy the list, as the system already tracks pl_discards elsewhere allocate(lmask, source=pl%ldiscard(:)) @@ -424,6 +425,7 @@ module subroutine symba_util_rearray_pl(self, system, param) allocate(lmask(pl%nbody)) lmask(:) = pl%status(1:pl%nbody) == NEW_PARTICLE + call symba_io_dump_particle_info(system, param, plidx=pack([(i, i=1, pl%nbody)], lmask)) where(pl%status(:) /= INACTIVE) pl%status(:) = ACTIVE @@ -431,8 +433,39 @@ module subroutine symba_util_rearray_pl(self, system, param) pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY pl%nplm = count(pl%lmtiny(:)) - ! Reindex + ! Reindex the bodies and calculate the level 0 encounter list call pl%eucl_index() + lencounter = pl%encounter_check(system, param%dt, 0) + select type(tp => system%tp) + class is (symba_tp) + lencounter = tp%encounter_check(system, param%dt, 0) + end select + + associate(idnew1 => system%plplenc_list%id1, idnew2 => system%plplenc_list%id2, idold1 => plplenc_old%id1, idold2 => plplenc_old%id2) + do k = 1, system%plplenc_list%nenc + if ((idnew1(k) == idold1(k)) .and. (idnew2(k) == idold2(k))) then + ! This is an encounter we already know about, so save the old information + system%plplenc_list%lvdotr(k) = plplenc_old%lvdotr(k) + system%plplenc_list%status(k) = plplenc_old%status(k) + system%plplenc_list%x1(:,k) = plplenc_old%x1(:,k) + system%plplenc_list%x2(:,k) = plplenc_old%x2(:,k) + system%plplenc_list%v1(:,k) = plplenc_old%v1(:,k) + system%plplenc_list%v2(:,k) = plplenc_old%v2(:,k) + system%plplenc_list%t(k) = plplenc_old%t(k) + system%plplenc_list%level(k) = plplenc_old%level(k) + else if (((idnew1(k) == idold2(k)) .and. (idnew2(k) == idold1(k)))) then + ! This is an encounter we already know about, but with the order reversed, so save the old information + system%plplenc_list%lvdotr(k) = plplenc_old%lvdotr(k) + system%plplenc_list%status(k) = plplenc_old%status(k) + system%plplenc_list%x1(:,k) = plplenc_old%x2(:,k) + system%plplenc_list%x2(:,k) = plplenc_old%x1(:,k) + system%plplenc_list%v1(:,k) = plplenc_old%v2(:,k) + system%plplenc_list%v2(:,k) = plplenc_old%v1(:,k) + system%plplenc_list%t(k) = plplenc_old%t(k) + system%plplenc_list%level(k) = plplenc_old%level(k) + end if + end do + end associate end if @@ -956,4 +989,4 @@ module subroutine symba_util_spill_tp(self, discards, lspill_list, ldestructive) return end subroutine symba_util_spill_tp -end submodule s_symba_util \ No newline at end of file +end submodule s_symba_util diff --git a/src/util/util_append.f90 b/src/util/util_append.f90 index 6f35de54e..88755870f 100644 --- a/src/util/util_append.f90 +++ b/src/util/util_append.f90 @@ -200,6 +200,8 @@ module subroutine util_append_encounter(self, source, lsource_mask) call util_append(self%status, source%status, nold, nsrc, lsource_mask) call util_append(self%index1, source%index1, nold, nsrc, lsource_mask) call util_append(self%index2, source%index2, nold, nsrc, lsource_mask) + call util_append(self%id1, source%id1, nold, nsrc, lsource_mask) + call util_append(self%id2, source%id2, nold, nsrc, lsource_mask) call util_append(self%x1, source%x1, nold, nsrc, lsource_mask) call util_append(self%x2, source%x2, nold, nsrc, lsource_mask) call util_append(self%v1, source%v1, nold, nsrc, lsource_mask) diff --git a/src/util/util_copy.f90 b/src/util/util_copy.f90 index 87634a419..c8407416d 100644 --- a/src/util/util_copy.f90 +++ b/src/util/util_copy.f90 @@ -17,6 +17,8 @@ module subroutine util_copy_encounter(self, source) self%status(1:n) = source%status(1:n) self%index1(1:n) = source%index1(1:n) self%index2(1:n) = source%index2(1:n) + self%id1(1:n) = source%id1(1:n) + self%id2(1:n) = source%id2(1:n) self%x1(:,1:n) = source%x1(:,1:n) self%x2(:,1:n) = source%x2(:,1:n) self%v1(:,1:n) = source%v1(:,1:n) diff --git a/src/util/util_spill.f90 b/src/util/util_spill.f90 index a26fbfad7..954406d95 100644 --- a/src/util/util_spill.f90 +++ b/src/util/util_spill.f90 @@ -260,6 +260,8 @@ module subroutine util_spill_encounter(self, discards, lspill_list, ldestructive call util_spill(keeps%status, discards%status, lspill_list, ldestructive) call util_spill(keeps%index1, discards%index1, lspill_list, ldestructive) call util_spill(keeps%index2, discards%index2, lspill_list, ldestructive) + call util_spill(keeps%id1, discards%id1, lspill_list, ldestructive) + call util_spill(keeps%id2, discards%id2, lspill_list, ldestructive) call util_spill(keeps%x1, discards%x1, lspill_list, ldestructive) call util_spill(keeps%x2, discards%x2, lspill_list, ldestructive) call util_spill(keeps%v1, discards%v1, lspill_list, ldestructive)