From c20c3680fd07fde016a148bae680526f1c3e8db1 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 28 Sep 2021 09:56:08 -0400 Subject: [PATCH] Fixed bug that was causing test particle encounters to be computed incorrectly in SyMBA. Fixed another initial conditions file --- .../swiftest_vs_swifter.ipynb | 12 +- .../1pl_1tp_encounter/cb.swiftest.in | Bin 80 -> 59 bytes .../1pl_1tp_encounter/init_cond.py | 108 ++++++++---------- .../1pl_1tp_encounter/param.swiftest.in | 12 +- .../1pl_1tp_encounter/pl.swifter.in | 4 +- .../1pl_1tp_encounter/pl.swiftest.in | Bin 176 -> 167 bytes .../swiftest_vs_swifter.ipynb | 32 +++--- .../1pl_1tp_encounter/tp.swifter.in | 2 +- .../1pl_1tp_encounter/tp.swiftest.in | Bin 128 -> 54 bytes src/symba/symba_encounter_check.f90 | 20 ++-- src/symba/symba_kick.f90 | 18 +-- 11 files changed, 103 insertions(+), 105 deletions(-) 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 090b756b6..eb88c82e3 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 @@ -75,23 +75,23 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "[,\n", - " ]" + "[,\n", + " ]" ] }, - "execution_count": 6, + "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAfMklEQVR4nO3de3Qc5Znn8e8jybZ8vwp8kW0JYoNtLo6jYCDEXDIQQzLxEGZ2MAwshAzDbODkcrIbZndz25nsSZZwDslCwjjcBjILkyEkC6wDmcQJzkBIbAeCbcAgbMcWNliW46t8k/rZP6pkt1qtVkvdpSp3/T7ntEtV71tVT0vu9+l6q+otc3dERCS9quIOQERE4qVEICKSckoEIiIpp0QgIpJySgQiIilXE3cA/TVp0iRvaGiIOwwRkRPKmjVrdrp7Xb6yEy4RNDQ0sHr16rjDEBE5oZjZH3orU9eQiEjKKRGIiKScEoGISMqdcOcIRETicvToUVpaWjh06FDcofSqtraW+vp6hgwZUvQ6SgQiIkVqaWlh9OjRNDQ0YGZxh9ODu9PW1kZLSwuNjY1Fr6euIRGRIh06dIiJEycmMgkAmBkTJ07s9xGLEoGISD8kNQl0GUh8SgQiIoXsb4UV/wDvvhp3JJFRIhARKeTADlh5B+x8oyybO//88/Muv+GGG3j88cfLso/+UiIQESkk0xlMq8pzbc0LL7xQlu2Uk64aEhEpJNMRTKuqIVP65kaNGsX+/ftxd2677TZWrFhBY2MjcT4tUkcEIiKFeNj6W3VZN/ujH/2IDRs2sHbtWr73ve/FeqSgRCAiUsixrqHyJoKVK1eydOlSqqurmTp1KpdccklZt98fSgQiIoVkdw2VWVIuRVUiEBEpxMMjgjJ3DS1atIjHHnuMzs5Otm/fzi9+8Yuybr8/dLJYRKSQiLqGrrzySlasWMGZZ57J7NmzufDCC8u6/f5QIhARKaTMl4/u378fCLqF7r777rJss1TqGhIRKSSirqEkUSIQESnk2BFB5TaXlfvORETKwcvbNZREkSUCM3vAzHaY2bpeyq81s1fC1wtmdnZUsYiIDFjX5aPqGhqQh4DFBco3ARe6+1nA3wPLIoxFRGRgIrpqKEkiO9Zx95Vm1lCgPPt+6heB+qhiEREZsG5DTJRhsKEESso5gpuAn/RWaGY3m9lqM1vd2to6iGGJSOpFeGfxQHziE5/gpJNO4owzzijbNmNPBGZ2MUEi+EJvddx9mbs3uXtTXV3d4AUnIpKwrqEbbriBZ555pqzbjDURmNlZwH3AEndvizMWEZG8EnYfwaJFi5gwYUJZtxnb9VBmNgN4ArjO3cvz6B8RkXLr5c7irz61nle37S3rruZOHcOX/3ReWbdZjMgSgZk9ClwETDKzFuDLwBAAd78X+BIwEfhOOAJfh7s3RRWPiMiAJKxrKApRXjW0tI/yTwKfjGr/IiJlcaxrqHtPehzf3KMS+8liEZFEK/Ogc0mkRCAiUkjCLh9dunQp5513Hhs2bKC+vp7777+/5G1WbooTESmHhF019Oijj5Z9mzoiEBEpJBPeTZyQI4IoKBGIiBSSsCOCKCgRiIgUkukATM8jEBFJrUxnRXcLgRKBiEhh3lnRl46CEoGISGGZzoo+PwBKBCIihSWoa2jr1q1cfPHFzJkzh3nz5vGtb32rLNut7OMdEZFSeWeP4SXiUlNTw5133smCBQvYt28f73vf+7j00kuZO3duSdtNxrsTEUmqTEdizhFMmTKFBQsWADB69GjmzJnD22+/XfJ2k/HuRESSqreuoZ/cDu+sLe++Jp8Jl3+9qKqbN2/mpZdeYuHChSXvVkcEIiKFePJOFu/fv5+rrrqKu+66izFjxpS8PR0RiIgUksnk7xoq8pt7uR09epSrrrqKa6+9lo9//ONl2aaOCERECsl0JOauYnfnpptuYs6cOXzuc58r23aT8e5ERJIqQV1Dzz//PI888ggrVqxg/vz5zJ8/n+XLl5e8XXUNiYgUkqD7CC644ALcvezb1RGBiEghCbp8NCpKBCIihXgmMV1DUVEiEBEpJNPZ7WRxFF0z5TSQ+CJLBGb2gJntMLN1vZSbmX3bzJrN7BUzWxBVLCIiA5Y1+mhtbS1tbW2JTQbuTltbG7W1tf1aL8qOr4eAu4GHeym/HJgVvhYC3w2nIiLJkek41jVUX19PS0sLra2tMQfVu9raWurr6/u1TmSJwN1XmllDgSpLgIc9SK0vmtk4M5vi7tujiklEpN+yrhoaMmQIjY2NMQdUfnGeI5gGbM2abwmX9WBmN5vZajNbneRMLCIVSCeLI2V5luXteHP3Ze7e5O5NdXV1EYclIpIl05GY+wiiEmciaAGmZ83XA9tiikVEJL8E3VAWlTgTwZPA9eHVQ+cCe3R+QEQSJ0FDTEQlspPFZvYocBEwycxagC8DQwDc/V5gOXAF0Ay0AzdGFYuIyIBlKv/h9VFeNbS0j3IHPhXV/kVEykJdQyIiKZegZxZHpbLfnYhIqXREICKSchp9VEQk5VJw1ZASgYhIIZmMuoZERFLNdY5ARCTdskYfrVRKBCIiheiqIRGRlPPKv7NYiUBEpJCMrhoSEUm3nGcWV6LKfnciIqXSfQQiIimnO4tFRFLMPXhUpa4aEhFJKc8EU3UNiYikVKYzmOqIQEQkpTIdwVSJQEQkpTw8IlDXkIhISqlrSEQk5Y4lAl0+KiKSTse6hiq7qYz03ZnZYjPbYGbNZnZ7nvKxZvaUmf3ezNab2Y1RxiMi0i/qGiqNmVUD9wCXA3OBpWY2N6fap4BX3f1s4CLgTjMbGlVMIiL94uoaKtU5QLO7b3T3I8BjwJKcOg6MNjMDRgG7gI4IYxIRKV7X5aO6amjApgFbs+ZbwmXZ7gbmANuAtcCn3btu5TvOzG42s9Vmtrq1tTWqeEVEulPXUMkszzLPmf8w8DIwFZgP3G1mY3qs5L7M3Zvcvamurq7ccYqI5KchJkrWAkzPmq8n+Oaf7UbgCQ80A5uA0yOMSUSkeLqzuGSrgFlm1hieAL4aeDKnzhbgQwBmdjJwGrAxwphERIqXkq6hyE6Fu3uHmd0KPAtUAw+4+3ozuyUsvxf4e+AhM1tL0JX0BXffGVVMIiL9kpIhJiK9JsrdlwPLc5bdm/XzNuCyKGMQERkw3VksIpJyKekaUiIQEemNhpgQEUk5HRGIiKTcsctHdY5ARCSdUnLVkBKBiEhvMuGdxeoaEhFJKd1ZLCKScuoaEhFJOV01JCKScjoiEBFJOQ0xISKScuoaCpjZTTnz1Wb25ehCEhFJCA0xccyHzGy5mU0xszOAF4HREcclIhK/lNxZ3Oe7c/drzOwvCZ4p3A4sdffnI49MRCRu6hoKmNks4NPAD4HNwHVmNiLiuERE4qerho55Cviiu/8NcCHwJsFjKEVEKtuxISZS3jUEnOPuewHc3YE7zSz32cMiIpXn2DmCyj5ZXMw5gr1mdj7QkFP/zaiCEhFJhJR0DfWZCMzsEeBU4GUg/K3gwMPRhSUikgApOVlcTNdQEzA37BYSEUmPlFw+WkzH1zpg8kA2bmaLzWyDmTWb2e291LnIzF42s/Vm9txA9iMiEgkPTxantWvIzJ4i6AIaDbxqZr8FDneVu/vHCm3YzKqBe4BLgRZglZk96e6vZtUZB3wHWOzuW8zspBLei4hIealriG+G0/cB9wNbAOvHts8Bmt19I4CZPQYsAV7NqnMN8IS7bwFw9x392L6ISLS8MxhewvrT9J14eu0acvfn3P05giOCvwO+BpwBvB4u78s0YGvWfEu4LNtsYLyZ/dLM1pjZ9fk2ZGY3m9lqM1vd2tpaxK5FRMog01Hx3UJQxDkCd/+qu88DPgVMBZ4zs58Vse18KTT3hHMNwRHHR4APA180s9l5Yljm7k3u3lRXV1fErkVEyiDTWfHdQlDcVUNddgDvAG1AMX35LcD0rPl6YFueOjvd/QBwwMxWAmcDb/QjLhGRaHhGRwQAZva3ZvZL4OfAJOCv3f2sIra9CphlZo1mNhS4Gsi9I/n/Ah80s5pw/KKFwGv9eQMiIpE5sh+GVv7QasUcEcwEPuPuL/dnw+7eYWa3As8C1cAD7r7ezG4Jy+9199fM7BngFSAD3Ofu6/r1DkREotLeBiMmxh1F5IoZYiLv9f/FcPflwPKcZffmzN8B3DHQfYiIRKb9j6lIBJU9kpKISCna22D4+LijiJwSgYhIbw7u0hGBiEhquYfnCCbEHUnklAhERPI5vDe4oUxHBCIiKdW+K5gO1xGBiEg6dSUCHRGIiKTUwa5EoCMCEZF0am8LpjoiEBFJqXYdEYiIpFt7W/AsgmFj444kckoEIiL5HNwVXDFUVfnNZOW/QxGRgUjJzWSgRCAikl97OoaXACUCEZH82nel4mYyUCIQEclPXUMiIinWcQQOtMLoyXFHMiiUCEREcv1xM3gnTJwVdySDQolARCRXW3MwnfieeOMYJEoEIiK52t4MphNPiTeOQaJEICKSq60ZRkxKxWMqQYlARKSntrdgUjrOD0DEicDMFpvZBjNrNrPbC9R7v5l1mtmfRxmPiEhR2pph4qlxRzFoIksEZlYN3ANcDswFlprZ3F7qfQN4NqpYRESKdmgv7H83NSeKIdojgnOAZnff6O5HgMeAJXnq3Qb8ENgRYSwiIsVJ2RVDEG0imAZszZpvCZcdY2bTgCuBeyOMQ0SkeDteC6aTTos3jkEUZSKwPMs8Z/4u4Avu3llwQ2Y3m9lqM1vd2tparvhERHra+iLUjk3VEUFNhNtuAaZnzdcD23LqNAGPmRnAJOAKM+tw9x9nV3L3ZcAygKamptxkIiJSPlt+A9MXpuI5BF2iTASrgFlm1gi8DVwNXJNdwd0bu342s4eAp3OTgIjIoGnfBTs3wFn/Ie5IBlVkicDdO8zsVoKrgaqBB9x9vZndEpbrvICIJMvW3wbTGefGG8cgi/KIAHdfDizPWZY3Abj7DVHGIiLSp60vQlUNTF0QdySDKj2dYCIifdm0EqbMh6Ej4o5kUCkRiIgA7N0Gb6+B0y6PO5JBp0QgIgLw+v8LpnM+Fm8cMVAiEBEBeO0pmDQb6mbHHcmgUyIQEdnfCpv/HU7/aNyRxEKJQETkpUeCR1OefXXckcRCiUBE0i3TCWsehIYPQl16xhfKpkQgIunW/DPYvQWaPhF3JLFRIhCR9HKHX90Jo6em9vwAKBGISJq99XPY+htY9HmoGRp3NLFRIhCRdMpkYMU/wNjp8N7r4o4mVpGONSQikli/+yfY9hJcuSzVRwOgIwIRSaP9rfCzr8DMC1I35HQ+SgQiki7u8OStcPQgfOROsHwPU0wXdQ2JSLr8dhm88Qws/gacdHrc0SSCjghEJD02/Qqe/a8wezEs/Ju4o0kMJQIRSYfWDfCD62HCKfDxZeoSyqJEICKVb9cmePjPgqePXfMvUDs27ogSRecIRKSyvbMWvn8VdB6B//h0cEQg3eiIQEQq1+bn4cErgiOBG5+ByWfEHVEiKRGISOVxhzUPwSNXwujJcNNPdYVQAeoaEpHKcng/PP1ZWPsDOOUiuOoBGDkx7qgSLdIjAjNbbGYbzKzZzG7PU36tmb0Svl4ws7OjjEdEKtzmf4d/XATrHoeL/zv81RNKAkWI7IjAzKqBe4BLgRZglZk96e6vZlXbBFzo7n80s8uBZcDCqGISkQp18I/wb18Oxg8aNxOufxIaPxh3VCeMKLuGzgGa3X0jgJk9BiwBjiUCd38hq/6LQH2E8YhIpek4AqsfgOe+AYd2w/m3wUV/B0NHxh3ZCSXKRDAN2Jo130Lhb/s3AT/JV2BmNwM3A8yYMaNc8YnIiarzKKz7YZAAdm2ExkVw2ddgyllxR3ZCijIR5Lttz/NWNLuYIBFckK/c3ZcRdBvR1NSUdxsikgJHD8JL34fnvw17tsBJ8+Dax+E9f6I7hUsQZSJoAaZnzdcD23IrmdlZwH3A5e7eFmE8InKi2tMCv3sEVt8PB1qh/hy44n/BrA9Dla6CL1WUiWAVMMvMGoG3gauBa7IrmNkM4AngOnd/I8JYROREk+mEN/8N1jwIb/4UPAPvuRQu+AzM/ICOAMooskTg7h1mdivwLFANPODu683slrD8XuBLwETgOxb8UTvcvSmqmEQk4dzh7TWw9nFY/yPY/w6MOhku+CwsuB7GN8QdYUUy9xOry72pqclXr14ddxgiUi6ZTmhZBRt+EjT+u/8A1cNg1qVw1l/CaZdD9ZC4ozzhmdma3r5o685iERl8h/bCWyuCB8S8+VNobwvGA2pcBBd+AeZ8VCOEDiIlAhGJXseRoMtn08rgtfU3kDkKteNg1mVw2mI49UMwfFzckaaSEoGIlN/RQ7D9ZdjyImz+Ffzh13D0AGAw5Ww492+Dp4RNXwjVaobipr+AiJRu77bgW/7WVcF0+++Db/wAk06D+dfAKRcGV/uMmBBvrNKDEoGIFM89aPTfeSVo7Ltee98OymtqYeoCOO9Twbf9+vfDqLp4Y5Y+KRGISH6dR4PhG3a8CtuzGv72nWEFg0mzYOb5QeM/fSFMPhNqhsYatvSfEoFI2nU1+K2vw47XofW1YNrWfLx7p6oGTpoT9OtPOTt4nTwPho2KN3YpCyUCkTRwDy7RbHsLdr0VTNuaYecbsPPN4w0+BuNnQt3pMPvDQeNfd3owrRkW61uQ6CgRiFQK92Bc/l0buzf4u96Cto1weM/xulYN42bApNnBjVt1c4JHOU6arSGcU0iJQOREcnA37N4S3H27e8vx1x/D+SP7siobjJ0OE0+Fs/4CJpwa/Dzh1CAJqC9fQkoEIknReRT2vRNclbP37ePT7IY++1s9wNBRwRO5xs2AhguC6YTGoLEf3wBDamN5K3JiUSIQGQwdh2Hf9rBxz2ro97QcX7b/XXo8smPIyKBxHz8TZp4X/HzsNROGj9conFIyJQKRUhw5EDTg+94Npl2vfe8GI2fufzf4ln+gtee6w8bAmKnB6+R5MGZaOD/t+PLasWroJXJKBCK5Ojvg4C7YvyNszHcEjXm++W598qGqGhh5Eow6KWjUpy6AsfXHG/cx02D0FKgdM/jvTSQPJQKpfJlOaN8V3Ah1oBUO7Axeeed3Blfe5Huq6tDRMPrkYHz8KWcH0675rtfoyTB8gp6aJScUJQI5sbjD0fagYT+4q/u0W+PeFkzbdwZleR+XbcG4NyMmwci64Fr5kXXhfLhs9OTgm/2ok3VZpVQsJQKJT2cHHNqdv1HPnh7MqdN5uPdtDh9/vDGvOw1GfiCncZ90fH7EBKiqHqx3K5JYSgRSuqMHg8b60O5wuifr53Car4E/tKf3bVbVBI368AlBgz2+Aaa99/h8b1M9yUqk35QIJOhuOXIgaLQP7clp1PNMc+sU+oYOQd/6iOxGvTFPIz6++/yw0bpaRmSQKBFUgq5+80N74fDecLqn+3yhb+uH9mSNNZOPBVe41I4LniBVOy4YjqB2XHB5Y9eyHtPxwSWSevCISKLpExo396Br5VgDvjdsrPfkLMud5pRnOgrvx6p7NtrjZvTRkI8LyoeN1VUwIhUs0kRgZouBbwHVwH3u/vWccgvLrwDagRvc/XdRxlRWvTXiuQ13j2X9bMSx4Jt17ZhwOja4Hn3Y6VnLssq65rN/HjpKXS0ikldkicDMqoF7gEuBFmCVmT3p7q9mVbscmBW+FgLfDafRK7YR7+uber8b8THBzUR1p/VswLMb7uzp0FH6Ri4ikYnyiOAcoNndNwKY2WPAEiA7ESwBHnZ3B140s3FmNsXdt5c7mHUrf0ztc/+Dkd7OCD/ACG9nCIUb8QzGQRtOu42k3UbQXjWSAzaSA9bAwaqRHKgdEZaNpL0qfNkIDobzB2wkh2x490b8CNgRg64bUnO+pAezDuwB9vT4Em9ZK/Qo60fdXJZTwbqV5YtxYOvmrp1d3nO7uWsWWLePuhTcT+/x9x1j4XW777f3+HPX7c977zkf4X66lRX+T9Wv/7s561WZYWbhz8fnu36uMrLKrVsdI1ynqvs6wfKsdaqCaU1VFUNrqhhWc3w6fGg1k0YNY8KIoVRVBdHdveJNfv76joLvuS+e73aWfvj4gmlcf15DaRvJI8pEMA3YmjXfQs9v+/nqTAO6JQIzuxm4GWDGjBkDCqZq2AgODJ3ITptBe9VIDoYN98GqkRy0kRyoGkW7HV/eXhU04m5VeB9/vexiz7lxaYgXXzd3N567zI/9U9S6x8u817K+1s0t7Htd772sHzHm6td+eqzrvZblLujPuoViyi3vqwEo2376WJcB/t5yy3v+bQv/P6FgjN2XZDx4jxkPytwh48H8YKupMuZNG8vieZP55k/f4MxpYxk3orRLlPtKoIXU1kRz30uUiSDfu839UxZTB3dfBiwDaGpqGtB/h7kLL4OFlw1kVRFJCA8TQpAYgiThWfP5kkehdTrd6cw4h49mONLZyeGjGQ53Zmg/3MnO/Yd5Z+8h/mXVVr7xzOvMnz6Of73lPIZUV143bZSJoAWYnjVfD2wbQB0RESD4Nl1tUJ33O2Q0rjt3JstWbuSTH2ysyCQAEOW7WgXMMrNGMxsKXA08mVPnSeB6C5wL7Ini/ICIyEBNHTecr3xsHvXjR8QdSmQiOyJw9w4zuxV4luDy0Qfcfb2Z3RKW3wssJ7h0tJng8tEbo4pHRETyi/Q+AndfTtDYZy+7N+tnBz4VZQwiIlJYZXZ4iYhI0ZQIRERSTolARCTllAhERFJOiUBEJOWUCEREUs76GuMlacysFfjDAFefBOwsYzhRSHqMSY8Pkh9j0uOD5Meo+PpvprvX5Ss44RJBKcxstbs3xR1HIUmPMenxQfJjTHp8kPwYFV95qWtIRCTllAhERFIubYlgWdwBFCHpMSY9Pkh+jEmPD5Ifo+Iro1SdIxARkZ7SdkQgIiI5lAhERFKuYhKBmS02sw1m1mxmt+cpNzP7dlj+ipktKHbdOOMzs+lm9gsze83M1pvZp6OIr5QYs8qrzewlM3s6afGZ2Tgze9zMXg9/l+clMMbPhn/jdWb2qJnVxhDf6Wb2azM7bGaf78+6ccaXsM9Jr7/DsDzSz8mABM/3PLFfBA++eQs4BRgK/B6Ym1PnCuAnBM9JPhf4TbHrxhzfFGBB+PNo4I1yx1dqjFnlnwP+D/B00uID/gn4ZPjzUGBckmIEpgGbgOHh/A+AG2KI7yTg/cDXgM/3Z92Y40vS5yRvjIPxORnoq1KOCM4Bmt19o7sfAR4DluTUWQI87IEXgXFmNqXIdWOLz923u/vvANx9H/AaQaNRbqX8DjGzeuAjwH0RxFZSfGY2BlgE3A/g7kfcfXeSYgzLaoDhZlYDjKD8z+/uMz533+Huq4Cj/V03zviS9Dkp8DscjM/JgFRKIpgGbM2ab6Hnf4Le6hSzbpzxHWNmDcB7gd+UOb6i9t9HnbuA/wJkIoit1PhOAVqBB8ND8vvMbGSSYnT3t4FvAluA7QTP7/5pDPFFsW6xyrKPBHxOCrmLaD8nA1IpicDyLMu9Lra3OsWsW6pS4gsKzUYBPwQ+4+57yxhbUfsvVMfMPgrscPc15Q+r8L6LrFMDLAC+6+7vBQ4AUfRxl/I7HE/wzbIRmAqMNLO/iiG+KNYtVsn7SMjnJP+Kg/M5GZBKSQQtwPSs+Xp6Hlb3VqeYdeOMDzMbQvCf+5/d/Ykyx1aOGD8AfMzMNhMcKl9iZt9PUHwtQIu7d31DfJwgMZRbKTH+CbDJ3Vvd/SjwBHB+DPFFsW6xStpHgj4nvRmMz8nAxH2Sohwvgm98Gwm+TXWdwJmXU+cjdD9J99ti1405PgMeBu5K6u8wp85FRHOyuKT4gF8Bp4U/fwW4I0kxAguB9QTnBozg5PZtgx1fVt2v0P1kbCI+JwXiS8znpLcYc8oi+ZwM+H3FHUAZ/0BXEFwp8Bbw38JltwC3ZP1HuScsXws0FVo3KfEBFxAcer4CvBy+rkhSjDnbiOw/eIl/4/nA6vD3+GNgfAJj/CrwOrAOeAQYFkN8kwm+9e4Fdoc/j0nQ5yRvfAn7nPT6OxyMz8lAXhpiQkQk5SrlHIGIiAyQEoGISMopEYiIpJwSgYhIyikRiIiknBKBpFo4Kul/ypqfamaPR7SvPzOzL/VR55tmdkkU+xfpjS4flVQLx6V52t3PGIR9vQB8zN13FqgzE/ieu18WdTwiXXREIGn3deBUM3vZzO4wswYzWwdgZjeY2Y/N7Ckz22Rmt5rZ58KB6140swlhvVPN7BkzW2NmvzKz03N3YmazgcPuvtPMRofbGxKWjTGzzWY2xN3/AEw0s8mD+DuQlFMikLS7HXjL3ee7+3/OU34GcA3B8MNfA9o9GLju18D1YZ1lBMNBvA/4PPCdPNv5AJA9TPIvCYabALga+KEHYwwR1vtAie9LpGg1cQcgknC/CBvufWa2B3gqXL4WOCsc7fJ84F/Njg1MOSzPdqYQDIXd5T6C4Yh/DNwI/HVW2Q6CEUhFBoUSgUhhh7N+zmTNZwg+P1XAbnef38d2DgJju2bc/fmwG+pCoNrd12XVrQ3riwwKdQ1J2u0jeLThgHgw5v0mM/sLOPZM4rPzVH0NeE/OsoeBR4EHc5bPJhh4TmRQKBFIqrl7G/B8+MD4Owa4mWuBm8zs9wRDSed7hONK4L2W1X8E/DMwniAZAMfG1H8PwUipIoNCl4+KDBIz+xbwlLv/LJz/c2CJu1+XVedKgoewfzGmMCWFdI5AZPD8T4IH0GBm/xu4nGBs+2w1wJ2DHJeknI4IRERSTucIRERSTolARCTllAhERFJOiUBEJOWUCEREUu7/A7uC4qt29TFlAAAAAElFTkSuQmCC\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAUoUlEQVR4nO3dfZBldX3n8fcnM0NGwxCew0BDZoJDMgNmcewAAgURJQVoGJWkltEoKBuWMlhmWdadrJU1+SO7VBmz6EqkBjEFxGXKxYeANUpUMKQgKIMgD44jE8DQzBhwquRhCQL63T/uhWranpnbv+7b97b9flV1zT3n/M45n2k48+lzzu1zU1VIkjRVvzDoAJKkuckCkSQ1sUAkSU0sEElSEwtEktRk4aADzKb999+/li1bNugYkjSn3HnnnT+sqgMmzp9XBbJs2TI2bdo06BiSNKck+f5k872EJUlqYoFIkppYIJKkJvPqHogkDcLzzz/P2NgYzz777KCj7NLixYsZGRlh0aJFPY23QCSpz8bGxliyZAnLli0jyaDjTKqq2LFjB2NjYyxfvryndbyEJUl99uyzz7LffvsNbXkAJGG//fab0lmSBSJJs2CYy+NFU81ogUiSmlggkjRHHH/88ZPOP/fcc7nuuutmOY0FIklzxm233TboCC/ju7AkaY7Yc889efrpp6kq3ve+93HTTTexfPlyBvXJsp6BSNIc8/nPf54tW7Zw7733csUVVwzszMQCkaQ55pZbbmHt2rUsWLCAgw8+mFNOOWUgOSwQSZqDhuFtwRaIJM0xJ510Ehs2bOAnP/kJ27dv5+abbx5IDm+iS9Ic89a3vpWbbrqJV7/61RxxxBGcfPLJA8lhgUjSHPH0008DnctXH//4xwecxktYkqRGFogkqYkFIklqYoFIkppYIJKkJhaIJKmJBSJJ88R73vMeDjzwQI466qgZ2Z4FIknzxLnnnsuXv/zlGdveQAskyWlJtiTZmmTdJMuT5GPd5fckWT1h+YIkdyX54uyllqS56aSTTmLfffedse0N7DfRkywALgNOBcaAO5JcX1XfGTfsdGBF9+tY4BPdP1/0fmAzsNeshJakafrzG+7nO9uenNFtrjp4Lz70u0fO6DZ7McgzkGOArVX1YFU9B2wA1kwYswa4ujpuB/ZOshQgyQjwJuCTsxlaktQxyGdhHQI8Mm56jJefXexszCHAduBS4APAkl3tJMn5wPkAhx122LQCS9J0DeJMoV8GeQYy2cPsJ34u46RjkrwZeKyq7tzdTqpqfVWNVtXoAQcc0JJTkjSJQRbIGHDouOkRYFuPY04AzkzyMJ1LX6ck+dv+RZWkuW/t2rW87nWvY8uWLYyMjHDllVdOa3uDvIR1B7AiyXLgUeBs4O0TxlwPXJhkA53LW09U1XbgT7pfJPlt4OKq+oNZyi1Jc9K11147o9sbWIFU1QtJLgRuBBYAn6qq+5Nc0F1+ObAROAPYCjwDvHtQeSVJLzfQD5Sqqo10SmL8vMvHvS7gj3azja8DX+9DPEnSLvib6JKkJhaIJKmJBSJJamKBSJKaWCCSNA888sgjvP71r2flypUceeSRfPSjH532Ngf6LixJ0uxYuHAhH/nIR1i9ejVPPfUUr33tazn11FNZtWpV8zY9A5GkeWDp0qWsXt35RIwlS5awcuVKHn300Wlt0zMQSZpNX1oHP7h3Zrd50Kvh9Et6Hv7www9z1113ceyxE59fOzWegUjSPPL0009z1llncemll7LXXtP7KCXPQCRpNk3hTGGmPf/885x11lm84x3v4G1ve9u0t+cZiCTNA1XFeeedx8qVK7noootmZJsWiCTNA7feeivXXHMNN910E0cffTRHH300Gzdu3P2Ku+AlLEmaB0488UQ6z6edOZ6BSJKaWCCSpCYWiCTNgpm+fNQPU81ogUhSny1evJgdO3YMdYlUFTt27GDx4sU9r+NNdEnqs5GREcbGxnj88ccHHWWXFi9ezMjISM/jLRBJ6rNFixaxfPnyQceYcV7CkiQ1sUAkSU0sEElSEwtEktTEApEkNbFAJElNLBBJUhMLRJLUxAKRJDWxQCRJTSwQSVKTgRZIktOSbEmyNcm6SZYnyce6y+9Jsro7/9AkNyfZnOT+JO+f/fSSNL8NrECSLAAuA04HVgFrk6yaMOx0YEX363zgE935LwD/uapWAscBfzTJupKkPhrkGcgxwNaqerCqngM2AGsmjFkDXF0dtwN7J1laVdur6lsAVfUUsBk4ZDbDS9J8N8gCOQR4ZNz0GD9bArsdk2QZ8BrgGzMfUZK0M4MskEwyb+LHde1yTJI9gc8Cf1xVT066k+T8JJuSbBr2D3ORpLlkkAUyBhw6bnoE2NbrmCSL6JTHp6vqczvbSVWtr6rRqho94IADZiS4JGmwBXIHsCLJ8iR7AGcD108Ycz3wru67sY4Dnqiq7UkCXAlsrqq/mt3YkiQY4EfaVtULSS4EbgQWAJ+qqvuTXNBdfjmwETgD2Ao8A7y7u/oJwDuBe5Pc3Z3336pq4yz+FSRpXkvVxNsOP79GR0dr06ZNg44hSXNKkjuranTifH8TXZLUxAKRJDWxQCRJTSwQSVITC0SS1MQCkSQ1sUAkSU0sEElSEwtEktTEApEkNbFAJElNLBBJUhMLRJLUxAKRJDWxQCRJTSwQSVITC0SS1MQCkSQ1sUAkSU0sEElSEwtEktTEApEkNbFAJElNLBBJUhMLRJLUxAKRJDWxQCRJTSwQSVKT3RZIkvMmTC9I8qH+RZIkzQW9nIG8IcnGJEuTHAXcDizpcy5J0pBbuLsBVfX2JP8euBd4BlhbVbf2PZkkaaj1cglrBfB+4LPAw8A7k7xyJnae5LQkW5JsTbJukuVJ8rHu8nuSrO51XUlSf/VyCesG4E+r6j8CJwMPAHdMd8dJFgCXAacDq4C1SVZNGHY6sKL7dT7wiSmsK0nqo91ewgKOqaonAaqqgI8kuX4G9n0MsLWqHgRIsgFYA3xn3Jg1wNXd/d6eZO8kS4FlPaw7Y27/6z9kyY8292PTkn7OfX/h4Vz1yxcMOgarDt6LD/3ukTO6zV7ugTyZ5Hg6/2iPH//ANPd9CPDIuOkx4NgexhzS47oAJDmfztkLhx122PQSS5JestsCSXINcDhwN/CT7uwCrp7mvjPJvOpxTC/rdmZWrQfWA4yOjk46ZneOe+8VLatJEkcCZww6RJ/0cglrFFjVvYw0k8aAQ8dNjwDbehyzRw/rSpL6qJeb6PcBB/Vh33cAK5IsT7IHcDYw8d7K9cC7uu/GOg54oqq297iuJKmPdnoGkuQGOpeFlgDfSfJN4Md0Lh9VVZ05nR1X1QtJLgRuBBYAn6qq+5Nc0F1+ObCRztnfVjq/g/LuXa07nTySpKnJzq5MJTm5+/JIYOI/zqmqr/cxV1+Mjo7Wpk2bBh1DkuaUJHdW1ejE+Tu9hFVV/1BV/wC8FzgOuAX4JvB7wP/sV1BJ0tzQyz2QY+ncsL6Nzr2HbcAJ/QwlSRp+vRTI88C/Aa8AFgMPVdVP+5pKkjT0eimQO+gUyG8BJ9J5bMh1fU0lSRp6vfweyHlV9eKd5x8Aa5K8s4+ZJElzwG7PQMaVx/h51/QnjiRprvAjbSVJTSwQSVITC0SS1MQCkSQ1sUAkSU0sEElSEwtEktTEApEkNbFAJElNLBBJUhMLRJLUxAKRJDWxQCRJTSwQSVITC0SS1MQCkSQ1sUAkSU0sEElSEwtEktTEApEkNbFAJElNLBBJUhMLRJLUxAKRJDUZSIEk2TfJV5I80P1zn52MOy3JliRbk6wbN//DSb6b5J4kn0+y96yFlyQBgzsDWQd8rapWAF/rTr9MkgXAZcDpwCpgbZJV3cVfAY6qqt8Evgf8yaykliS9ZFAFsga4qvv6KuAtk4w5BthaVQ9W1XPAhu56VNXfV9UL3XG3AyP9jStJmmhQBfIrVbUdoPvngZOMOQR4ZNz0WHfeRO8BvjTjCSVJu7SwXxtO8lXgoEkWfbDXTUwyrybs44PAC8Cnd5HjfOB8gMMOO6zHXUuSdqdvBVJVb9zZsiT/mmRpVW1PshR4bJJhY8Ch46ZHgG3jtnEO8GbgDVVV7ERVrQfWA4yOju50nCRpagZ1Cet64Jzu63OAv5tkzB3AiiTLk+wBnN1djySnAf8VOLOqnpmFvJKkCQZVIJcApyZ5ADi1O02Sg5NsBOjeJL8QuBHYDHymqu7vrv9xYAnwlSR3J7l8tv8CkjTf9e0S1q5U1Q7gDZPM3wacMW56I7BxknGv6mtASdJu+ZvokqQmFogkqYkFIklqYoFIkppYIJKkJhaIJKmJBSJJamKBSJKaWCCSpCYWiCSpiQUiSWpigUiSmlggkqQmFogkqYkFIklqYoFIkppYIJKkJhaIJKmJBSJJamKBSJKaWCCSpCYWiCSpiQUiSWpigUiSmlggkqQmFogkqYkFIklqYoFIkppYIJKkJhaIJKmJBSJJajKQAkmyb5KvJHmg++c+Oxl3WpItSbYmWTfJ8ouTVJL9+59akjTeoM5A1gFfq6oVwNe60y+TZAFwGXA6sApYm2TVuOWHAqcC/zIriSVJLzOoAlkDXNV9fRXwlknGHANsraoHq+o5YEN3vRf9L+ADQPUxpyRpJwZVIL9SVdsBun8eOMmYQ4BHxk2PdeeR5Ezg0ar69u52lOT8JJuSbHr88cenn1ySBMDCfm04yVeBgyZZ9MFeNzHJvEryyu42fqeXjVTVemA9wOjoqGcrkjRD+lYgVfXGnS1L8q9JllbV9iRLgccmGTYGHDpuegTYBhwOLAe+neTF+d9KckxV/WDG/gKSpF0a1CWs64Fzuq/PAf5ukjF3ACuSLE+yB3A2cH1V3VtVB1bVsqpaRqdoVlsekjS7BlUglwCnJnmAzjupLgFIcnCSjQBV9QJwIXAjsBn4TFXdP6C8kqQJ+nYJa1eqagfwhknmbwPOGDe9Edi4m20tm+l8kqTd8zfRJUlNLBBJUhMLRJLUxAKRJDWxQCRJTSwQSVITC0SS1MQCkSQ1sUAkSU0sEElSEwtEktTEApEkNbFAJElNLBBJUhMLRJLUxAKRJDWxQCRJTSwQSVITC0SS1MQCkSQ1sUAkSU0sEElSEwtEktTEApEkNUlVDTrDrEnyOPD9xtX3B344g3H6Ydgzmm/6hj3jsOeD4c84jPl+taoOmDhzXhXIdCTZVFWjg86xK8Oe0XzTN+wZhz0fDH/GYc83npewJElNLBBJUhMLpHfrBx2gB8Oe0XzTN+wZhz0fDH/GYc/3Eu+BSJKaeAYiSWpigUiSmlggQJLTkmxJsjXJukmWJ8nHusvvSbK613UHmS/JoUluTrI5yf1J3j9M+cYtX5DkriRf7Ee+6WZMsneS65J8t/u9fN2Q5ftP3f++9yW5Nsnimc7XY8bfSPJPSX6c5OKprDvIfLN1nEwn47jlfT9WpqSq5vUXsAD4Z+DXgD2AbwOrJow5A/gSEOA44Bu9rjvgfEuB1d3XS4DvDVO+ccsvAv4P8MVh+2/cXXYV8B+6r/cA9h6WfMAhwEPAK7rTnwHOHdD38EDgt4C/AC6eyroDztf342S6GWfrWJnql2cgcAywtaoerKrngA3Amglj1gBXV8ftwN5Jlva47sDyVdX2qvoWQFU9BWym8w/OUOQDSDICvAn45AznmpGMSfYCTgKuBKiq56rqR8OSr7tsIfCKJAuBVwLbZjhfTxmr6rGqugN4fqrrDjLfLB0n08oIs3asTIkF0vkf5ZFx02P87P88OxvTy7qDzPeSJMuA1wDfGLJ8lwIfAH46w7l63f/uxvwa8DjwN91LB59M8kvDkq+qHgX+EvgXYDvwRFX9/Qzn6zVjP9bt1Yzso4/HCUw/46X0/1iZEgukc0lgoonvbd7ZmF7Wna7p5OssTPYEPgv8cVU9OYPZdrvvXY1J8mbgsaq6c4YzTTSd7+FCYDXwiap6DfD/gJm+hj+d7+E+dH6KXQ4cDPxSkj+Y4Xw73f8srNurae+jz8cJTCPjLB4rU2KBdH4KOHTc9Ag/ewlgZ2N6WXeQ+UiyiM5B8emq+twMZ5tuvhOAM5M8TOd0/pQkfztkGceAsap68SfS6+gUyrDkeyPwUFU9XlXPA58Djp/hfL1m7Me6vZrWPmbhOIHpZZytY2VqBn0TZtBfdH7CfJDOT3Av3tg6csKYN/HyG5jf7HXdAecLcDVw6TB+/yaM+W36dxN9WhmBfwR+vfv6z4APD0s+4Fjgfjr3PkLnhv/7BvE9HDf2z3j5TeqhOE52ka/vx8l0M05Y1rdjZcp/p0EHGIYvOu9w+R6dd0h8sDvvAuCC7usAl3WX3wuM7mrdYckHnEjnFPke4O7u1xnDkm/CNvp6UEzzv/HRwKbu9/ELwD5Dlu/Pge8C9wHXAL84oO/hQXR+yn4S+FH39V5DdJxMmm+2jpPpfg9n61iZypePMpEkNfEeiCSpiQUiSWpigUiSmlggkqQmFogkqYkFIjXqPqX3veOmD05yXZ/29ZYk/303Y/4yySn92L80Gd/GKzXqPjfpi1V11Czs6zbgzKr64S7G/CpwRVX9Tr/zSOAZiDQdlwCHJ7k7yYeTLEtyH0CSc5N8IckNSR5KcmGSi7oPZLw9yb7dcYcn+XKSO5P8Y5LfmLiTJEcAP66qHyZZ0t3eou6yvZI8nGRRVX0f2C/JQbP4PdA8ZoFI7dYB/1xVR1fVf5lk+VHA2+k8xvsvgGeq80DGfwLe1R2zns6jR14LXAz89STbOQEY/7jxr9N5tAnA2cBnq/McLLrjTpjm30vqycJBB5B+jt3c/Qf/qSRPADd0598L/Gb36a/HA/83eelBrb84yXaW0nmk/Is+Seex3l8A3g384bhlj9F5Kq/UdxaI1D8/Hvf6p+Omf0rn2PsF4EdVdfRutvNvwC+/OFFVt3Yvl50MLKiq+8aNXdwdL/Wdl7Ckdk/R+QjUJtX5zImHkvw+vPS55/9ukqGbgVdNmHc1cC3wNxPmH0HnoYpS31kgUqOq2gHcmuS+JB9u3Mw7gPOSfJvOY9kn+6jXW4DXZNx1LuDTwD50SgR46TMtXkXnycFS3/k2XmkOSPJR4Iaq+mp3+veANVX1znFj3gqsrqo/HVBMzTPeA5Hmhv9B58OjSPK/gdPpfLbEeAuBj8xyLs1jnoFIkpp4D0SS1MQCkSQ1sUAkSU0sEElSEwtEktTk/wMOj27WwTLGZgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] @@ -103,7 +103,7 @@ } ], "source": [ - "swiftdiff['vhx'].plot.line(x=\"time (y)\")" + "swiftdiff['xhx'].plot.line(x=\"time (y)\")" ] }, { diff --git a/examples/symba_swifter_comparison/1pl_1tp_encounter/cb.swiftest.in b/examples/symba_swifter_comparison/1pl_1tp_encounter/cb.swiftest.in index d0ae0ed15fe3ea8dd15557055a926fce3c60b59c..06036db504f8af9b71dc87951a1e51600548a923 100644 GIT binary patch literal 59 zcmW;Bu@QhU3+AsI-}OJ>6US3*597mhVB-S-U7iOf diff --git a/examples/symba_swifter_comparison/1pl_1tp_encounter/init_cond.py b/examples/symba_swifter_comparison/1pl_1tp_encounter/init_cond.py index b516bda9c..f9974a19a 100755 --- a/examples/symba_swifter_comparison/1pl_1tp_encounter/init_cond.py +++ b/examples/symba_swifter_comparison/1pl_1tp_encounter/init_cond.py @@ -18,7 +18,7 @@ swiftest_pl = "pl.swiftest.in" swiftest_tp = "tp.swiftest.in" swiftest_cb = "cb.swiftest.in" -swiftest_bin = "bin.swiftest.dat" +swiftest_bin = "bin.swiftest.nc" swiftest_enc = "enc.swiftest.dat" swiftest_dis = "discard.swiftest.dat" @@ -43,11 +43,16 @@ rmax = 1000.0 npl = 1 -plid = 2 -tpid = 100 +ntp = 1 +plid = 1 +tpid = 2 + +cbname = "Sun" +plname = "Planet" +tpname = "TestParticle" radius = np.double(4.25875607065041e-05) -mass = np.double(0.00012002693582795244940133) +Gmass = np.double(0.00012002693582795244940133) apl = np.longdouble(1.0) atp = np.longdouble(1.01) vpl = np.longdouble(2 * np.pi) @@ -62,23 +67,21 @@ rhill = np.double(apl * 0.0100447248332378922085) #Make Swifter files -plfile = open(swifter_pl, 'w') -print(npl+1, f'! Planet input file generated using init_cond.py',file=plfile) -print(1,GMSun,file=plfile) -print('0.0 0.0 0.0',file=plfile) -print('0.0 0.0 0.0',file=plfile) -print(plid,"{:.23g}".format(mass),rhill, file=plfile) -print(radius, file=plfile) -print(*p_pl, file=plfile) -print(*v_pl, file=plfile) -plfile.close() - -tpfile = open(swifter_tp, 'w') -print(1,file=tpfile) -print(tpid, file=tpfile) -print(*p_tp, file=tpfile) -print(*v_tp, file=tpfile) -tpfile.close() +with open(swifter_pl, 'w') as plfile: + print(npl+1, f'! Planet input file generated using init_cond.py',file=plfile) + print(0,GMSun,file=plfile) + print('0.0 0.0 0.0',file=plfile) + print('0.0 0.0 0.0',file=plfile) + print(plid,"{:.23g}".format(Gmass),rhill, file=plfile) + print(radius, file=plfile) + print(*p_pl, file=plfile) + print(*v_pl, file=plfile) + +with open(swifter_tp, 'w') as tpfile: + print(1,file=tpfile) + print(tpid, file=tpfile) + print(*p_tp, file=tpfile) + print(*v_tp, file=tpfile) sys.stdout = open(swifter_input, "w") print(f'! Swifter input file generated using init_cond.py') @@ -110,41 +113,25 @@ sys.stdout = sys.__stdout__ #Now make Swiftest files -cbfile = FortranFile(swiftest_cb, 'w') -Msun = np.double(1.0) -cbfile.write_record(0) -cbfile.write_record(np.double(GMSun)) -cbfile.write_record(np.double(rmin)) -cbfile.write_record(np.double(J2)) -cbfile.write_record(np.double(J4)) -cbfile.close() - -plfile = FortranFile(swiftest_pl, 'w') -plfile.write_record(npl) - -plfile.write_record(plid) -plfile.write_record(p_pl[0]) -plfile.write_record(p_pl[1]) -plfile.write_record(p_pl[2]) -plfile.write_record(v_pl[0]) -plfile.write_record(v_pl[1]) -plfile.write_record(v_pl[2]) -plfile.write_record(mass) -plfile.write_record(rhill) -plfile.write_record(radius) -plfile.close() -tpfile = FortranFile(swiftest_tp, 'w') -ntp = 1 -tpfile.write_record(ntp) -tpfile.write_record(tpid) -tpfile.write_record(p_tp[0]) -tpfile.write_record(p_tp[1]) -tpfile.write_record(p_tp[2]) -tpfile.write_record(v_tp[0]) -tpfile.write_record(v_tp[1]) -tpfile.write_record(v_tp[2]) - -tpfile.close() +with open(swiftest_cb, 'w') as cbfile: + print(cbname,file=cbfile) + print(GMSun, file=cbfile) + print(rmin, file=cbfile) + print(J2, file=cbfile) + print(J4, file=cbfile) + +with open(swiftest_pl, 'w') as plfile: + print(npl, f'! Planet input file generated using init_cond.py',file=plfile) + print(plname,"{:.23g}".format(Gmass),rhill, file=plfile) + print(radius, file=plfile) + print(*p_pl, file=plfile) + print(*v_pl, file=plfile) + +with open(swiftest_tp, 'w') as tpfile: + print(ntp,file=tpfile) + print(tpname, file=tpfile) + print(*p_tp, file=tpfile) + print(*v_tp, file=tpfile) sys.stdout = open(swiftest_input, "w") print(f'! Swiftest input file generated using init_cond.py') @@ -154,12 +141,13 @@ print(f'CB_IN {swiftest_cb}') print(f'PL_IN {swiftest_pl}') print(f'TP_IN {swiftest_tp}') -print(f'IN_TYPE REAL8') +print(f'IN_TYPE ASCII') +print(f'IN_FORM XV') print(f'ISTEP_OUT {iout:d}') -print(f'ISTEP_DUMP {iout:d}') +print(f'ISTEP_DUMP {10*iout:d}') print(f'BIN_OUT {swiftest_bin}') -print(f'OUT_TYPE REAL8') -print(f'OUT_FORM XV') +print(f'OUT_TYPE NETCDF_DOUBLE') +print(f'OUT_FORM XVEL') print(f'OUT_STAT REPLACE') print(f'CHK_CLOSE yes') print(f'CHK_RMIN {rmin}') diff --git a/examples/symba_swifter_comparison/1pl_1tp_encounter/param.swiftest.in b/examples/symba_swifter_comparison/1pl_1tp_encounter/param.swiftest.in index 9fb0bf743..ed04ce1b3 100644 --- a/examples/symba_swifter_comparison/1pl_1tp_encounter/param.swiftest.in +++ b/examples/symba_swifter_comparison/1pl_1tp_encounter/param.swiftest.in @@ -5,12 +5,13 @@ DT 0.0006844626967830253 CB_IN cb.swiftest.in PL_IN pl.swiftest.in TP_IN tp.swiftest.in -IN_TYPE REAL8 +IN_TYPE ASCII +IN_FORM XV ISTEP_OUT 1 -ISTEP_DUMP 1 -BIN_OUT bin.swiftest.dat -OUT_TYPE REAL8 -OUT_FORM XV +ISTEP_DUMP 10 +BIN_OUT bin.swiftest.nc +OUT_TYPE NETCDF_DOUBLE +OUT_FORM XVEL OUT_STAT REPLACE CHK_CLOSE yes CHK_RMIN 0.004650467260962157 @@ -30,3 +31,4 @@ DU2M 149597870700.0 TU2S 31557600.0 RHILL_PRESENT yes GMTINY 1e-12 +INTERACTION_LOOPS TRIANGULAR diff --git a/examples/symba_swifter_comparison/1pl_1tp_encounter/pl.swifter.in b/examples/symba_swifter_comparison/1pl_1tp_encounter/pl.swifter.in index 17d461561..fb6b10800 100644 --- a/examples/symba_swifter_comparison/1pl_1tp_encounter/pl.swifter.in +++ b/examples/symba_swifter_comparison/1pl_1tp_encounter/pl.swifter.in @@ -1,8 +1,8 @@ 2 ! Planet input file generated using init_cond.py -1 39.476926408897625196 +0 39.476926408897625196 0.0 0.0 0.0 0.0 0.0 0.0 -2 0.00012002693582795244940133 0.010044724833237892 +1 0.00012002693582795244940133 0.010044724833237892 4.25875607065041e-05 1.0 0.0 0.0 0.0 6.283185307179586 0.0 diff --git a/examples/symba_swifter_comparison/1pl_1tp_encounter/pl.swiftest.in b/examples/symba_swifter_comparison/1pl_1tp_encounter/pl.swiftest.in index c94c6ae61581655ba2acb43d632ec99b4c8d1cc3..3bc644cdb8e9f9226a6c2b67b786315583c61672 100644 GIT binary patch literal 167 zcmXYp%MHRX6a??C!W4kD@8?e!bbt^sk)_B);2h|l1aVl+VMa^-J_!7DiVUnUPXfAtc3 T)sRUoaK4cFghIv;HZUImPoWF< diff --git a/examples/symba_swifter_comparison/1pl_1tp_encounter/swiftest_vs_swifter.ipynb b/examples/symba_swifter_comparison/1pl_1tp_encounter/swiftest_vs_swifter.ipynb index 78c6cf2b0..56c39bb4c 100644 --- a/examples/symba_swifter_comparison/1pl_1tp_encounter/swiftest_vs_swifter.ipynb +++ b/examples/symba_swifter_comparison/1pl_1tp_encounter/swiftest_vs_swifter.ipynb @@ -43,8 +43,8 @@ "output_type": "stream", "text": [ "Reading Swiftest file param.swiftest.in\n", - "Reading in time 1.506e-01\n", - "Creating Dataset\n", + "\n", + "Creating Dataset from NetCDF file\n", "Successfully converted 221 output frames.\n", "Swiftest simulation data stored as xarray DataSet .ds\n" ] @@ -81,8 +81,8 @@ { "data": { "text/plain": [ - "[,\n", - " ]" + "[,\n", + " ]" ] }, "execution_count": 6, @@ -91,7 +91,7 @@ }, { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEGCAYAAABGnrPVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAWK0lEQVR4nO3dfZBddZ3n8fd3kkCGITwT6NDB9JjAJAEWY2+IaKGCmQrRSVRmLDLOEHyiIuLDMqybGWt3xtoaTZXjLrpmpIJIJY47KRflQSvARMDFwgkSHgRCjMkAkg4tiVGQrMuj3/3j3mRvOjfdN31/9yHk/aq61fec8/2d8+2bPvn0Oef2uZGZSJJUyu91ugFJ0muLwSJJKspgkSQVZbBIkooyWCRJRY3tdAPtdMIJJ+SUKVM63YYkHVTuv//+X2bmiY3WH1LBMmXKFNavX9/pNiTpoBIRPz+Qek+FSZKKMlgkSUUZLJKkog6payySNJyXX36ZgYEBXnjhhU630hHjx4+nt7eXcePGNbUeg0WSqgYGBpgwYQJTpkwhIjrdTltlJjt37mRgYIC+vr6m1uWpMEmqeuGFFzj++OMPuVABiAiOP/74IkdrBosk1TgUQ2W3Ut+7wSJJKspgkaQOOvfcc+vOv/TSS7nhhhva3E0ZBoskddCPfvSjTrdQnO8Kk6QOOvLII9m1axeZycc//nHuvPNO+vr6OJg/3dcjFknqAjfeeCObNm3ikUce4dprrz2oj2QMFknqAnfffTeLFi1izJgxTJo0ifPPP7/TLY2awSJJXeK18lZng0WSusB5553H6tWrefXVVxkcHOSuu+7qdEuj5sV7SeoC73nPe7jzzjs588wzOe2003jrW9/a6ZZGzWCRpA7atWsXUDkN9pWvfKXD3ZThqTBJUlEGiySpKINFklSUwSJJKspgkSQVZbBIkooyWCSpi2zdupW3v/3tTJ8+nZkzZ/KlL31pn5rM5BOf+ARTp07lrLPO4oEHHuhAp/vn37FIUhcZO3YsX/ziF5k1axbPP/88b3zjG5k7dy4zZszYU3PrrbeyefNmNm/ezL333stHP/pR7r333g52vbeOHrFExLyI2BQRWyJiaZ3lERFfri5/OCJmDVk+JiIejIjvta9rSWqdnp4eZs2q/Fc3YcIEpk+fzrZt2/aqufnmm7nkkkuICObMmcOzzz7L4OBgJ9qtq2NHLBExBlgOzAUGgPsi4pbMfKym7EJgWvVxDvDV6tfdPglsBI5qS9OSDhmf/e4GHnv6N0XXOWPSUfztn8xsuP7JJ5/kwQcf5Jxzztlr/rZt25g8efKe6d7eXrZt20ZPT0+xXpvRySOW2cCWzHw8M18CVgMLh9QsBFZlxTrgmIjoAYiIXuCdwNfa2bQktcOuXbu46KKLuPrqqznqqL1/d673IWDddGfkTl5jOQXYWjM9wN5HI/urOQUYBK4GPg1MGG4jEXEZcBnAqaee2lTDkg4dB3JkUdrLL7/MRRddxPvf/37e+9737rO8t7eXrVv//3+NAwMDTJo0qZ0tDquTRyz14nVoDNetiYh3Adsz8/6RNpKZKzKzPzP7TzzxxNH0KUltk5l86EMfYvr06Vx55ZV1axYsWMCqVavITNatW8fRRx/dNafBoLNHLAPA5JrpXuDpBmv+FFgQEfOB8cBREfFPmfkXLexXklrunnvu4Rvf+AZnnnkmZ599NgCf+9zneOqppwBYsmQJ8+fPZ82aNUydOpUjjjiC66+/voMd76uTwXIfMC0i+oBtwMXAnw+puQW4IiJWUzlN9lxmDgJ/XX0QEW8DrjJUJL0WvOUtb6l7DaVWRLB8+fI2dXTgOhYsmflKRFwB3A6MAb6emRsiYkl1+TXAGmA+sAX4LfCBTvUrSWpMR/9AMjPXUAmP2nnX1DxP4GMjrOMHwA9a0J4kaRS8pYskqSiDRZJUlMEiSSrKYJEkFWWwSFIX+eAHP8jEiRM544wz9sz71a9+xdy5c5k2bRpz587l17/+9Z5ln//855k6dSqnn346t99+e911Dje+FQwWSeoil156Kbfddtte85YtW8YFF1zA5s2bueCCC1i2bBkAjz32GKtXr2bDhg3cdtttXH755bz66qv7rHN/41vFYJGkLnLeeedx3HHH7TXv5ptvZvHixQAsXryYm266ac/8iy++mMMPP5y+vj6mTp3Kj3/8433Wub/xreIHfUlSPbcuhV88UnadJ58JFx740cIzzzyz515gPT09bN++HajcPn/OnDl76nbfPr/R8a3iEYskHaS69fb5HrFIUj2jOLJolZNOOonBwUF6enoYHBxk4sSJQOO3z9/f+FbxiEWSutyCBQtYuXIlACtXrmThwoV75q9evZoXX3yRJ554gs2bNzN79uyGx7eKwSJJXWTRokW86U1vYtOmTfT29nLdddexdOlS1q5dy7Rp01i7di1Lly4FYObMmbzvfe9jxowZzJs3j+XLlzNmzBgAPvzhD7N+/XqA/Y5vlRjp9syvJf39/bn7hZakoTZu3Mj06dM73UZH1XsNIuL+zOxvdB0esUiSijJYJElFGSySVONQujwwVKnv3WCRpKrx48ezc+fOQzJcMpOdO3cyfvz4ptfl37FIUlVvby8DAwPs2LGj0610xPjx4+nt7W16PQaLJFWNGzeOvr6+Trdx0PNUmCSpKINFklSUwSJJKspgkSQVZbBIkooyWCRJRRkskqSiDBZJUlEGiySpKINFklSUwSJJKqqjwRIR8yJiU0RsiYh9PiszKr5cXf5wRMyqzp8cEXdFxMaI2BARn2x/95KkejoWLBExBlgOXAjMABZFxIwhZRcC06qPy4CvVue/AvxVZk4H5gAfqzNWktQBnTximQ1syczHM/MlYDWwcEjNQmBVVqwDjomInswczMwHADLzeWAjcEo7m5ck1dfJYDkF2FozPcC+4TBiTURMAd4A3Fu+RUnSgepksESdeUM/tm3Ymog4Evg28KnM/E3djURcFhHrI2L9ofrhPZLUTp0MlgFgcs10L/B0ozURMY5KqHwzM7+zv41k5orM7M/M/hNPPLFI45Kk/etksNwHTIuIvog4DLgYuGVIzS3AJdV3h80BnsvMwYgI4DpgY2b+t/a2LUkaTsc+mjgzX4mIK4DbgTHA1zNzQ0QsqS6/BlgDzAe2AL8FPlAd/mbgL4FHIuKh6ry/ycw1bfwWJEl1RObQyxqvXf39/bl+/fpOtyFJB5WIuD8z+xut9y/vJUlFGSySpKIMFklSUQaLJKkog0WSVJTBIkkqymCRJBVlsEiSijJYJElFGSySpKIMFklSUQaLJKkog0WSVJTBIkkqymCRJBVlsEiSijJYJElFGSySpKIMFklSUQaLJKkog0WSVJTBIkkqymCRJBVlsEiSijJYJElFGSySpKIMFklSUQaLJKkog0WSVJTBIkkqymCRJBU1YrBExMQ6804vsfGImBcRmyJiS0QsrbM8IuLL1eUPR8SsRsdKkjqjkSOWH0bE+3ZPRMRfATc2u+GIGAMsBy4EZgCLImLGkLILgWnVx2XAVw9grCSpA8Y2UPM2YEVE/BlwErARmF1g27OBLZn5OEBErAYWAo/V1CwEVmVmAusi4piI6AGmNDC2mHX/+BEmPLuxFauWpJb6+djXs/LoJcyYdBR/+ycz27LNEY9YMnMQuA14E5X/0Fdl5q4C2z4F2FozPVCd10hNI2MBiIjLImJ9RKzfsWNH001LkoY34hFLRKwFBoEzgF7g6xFxd2Ze1eS2o868bLCmkbGVmZkrgBUA/f39dWtGMufya0czTJI6biYwv83bbOQay63A32Tms5n5KHAu8FyBbQ8Ak2ume4GnG6xpZKwkqQMaCZYJwO0R8cOI+BhwfGb+1wLbvg+YFhF9EXEYcDFwy5CaW4BLqu8OmwM8Vz0118hYSVIHNHKN5bOZORP4GDAJ+N8R8f1mN5yZrwBXALdTeUPAtzJzQ0QsiYgl1bI1wOPAFuBa4PLhxjbbkySpeY28K2y37cAvgJ3APn/bMhqZuYZKeNTOu6bmeVIJtIbGSpI6r5E/kPxoRPwAuAM4AfhIZp7V6sYkSQenRo5YXgd8KjMfanEvkqTXgBGDJTO9XYokqWHehFKSVJTBIkkqymCRJBVlsEiSijJYJElFGSySpKIMFklSUQaLJKkog0WSVJTBIkkqymCRJBVlsEiSijJYJElFGSySpKIMFklSUQaLJKkog0WSVJTBIkkqymCRJBVlsEiSijJYJElFGSySpKIMFklSUQaLJKkog0WSVJTBIkkqymCRJBXVkWCJiOMiYm1EbK5+PXY/dfMiYlNEbImIpTXzvxARP42IhyPixog4pm3NS5KG1akjlqXAHZk5DbijOr2XiBgDLAcuBGYAiyJiRnXxWuCMzDwL+Bnw123pWpI0ok4Fy0JgZfX5SuDddWpmA1sy8/HMfAlYXR1HZv5LZr5SrVsH9La2XUlSozoVLCdl5iBA9evEOjWnAFtrpgeq84b6IHBr8Q4lSaMytlUrjojvAyfXWfSZRldRZ14O2cZngFeAbw7Tx2XAZQCnnnpqg5uWJI1Wy4IlM9+xv2UR8UxE9GTmYET0ANvrlA0Ak2ume4Gna9axGHgXcEFmJvuRmSuAFQD9/f37rZMkldGpU2G3AIurzxcDN9epuQ+YFhF9EXEYcHF1HBExD/hPwILM/G0b+pUkNahTwbIMmBsRm4G51WkiYlJErAGoXpy/Argd2Ah8KzM3VMd/BZgArI2IhyLimnZ/A5Kk+lp2Kmw4mbkTuKDO/KeB+TXTa4A1deqmtrRBSdKo+Zf3kqSiDBZJUlEGiySpKINFklSUwSJJKspgkSQVZbBIkooyWCRJRRkskqSiDBZJUlEGiySpKINFklSUwSJJKspgkSQVZbBIkooyWCRJRRkskqSiDBZJUlEGiySpKINFklSUwSJJKspgkSQVZbBIkooyWCRJRRkskqSiDBZJUlEGiySpKINFklSUwSJJKspgkSQVZbBIkorqSLBExHERsTYiNle/HrufunkRsSkitkTE0jrLr4qIjIgTWt+1JKkRnTpiWQrckZnTgDuq03uJiDHAcuBCYAawKCJm1CyfDMwFnmpLx5KkhnQqWBYCK6vPVwLvrlMzG9iSmY9n5kvA6uq43f478GkgW9inJOkAdSpYTsrMQYDq14l1ak4BttZMD1TnERELgG2Z+ZORNhQRl0XE+ohYv2PHjuY7lyQNa2yrVhwR3wdOrrPoM42uos68jIgjquv440ZWkpkrgBUA/f39Ht1IUou1LFgy8x37WxYRz0RET2YORkQPsL1O2QAwuWa6F3gaeD3QB/wkInbPfyAiZmfmL4p9A5KkUenUqbBbgMXV54uBm+vU3AdMi4i+iDgMuBi4JTMfycyJmTklM6dQCaBZhookdYdOBcsyYG5EbKbyzq5lABExKSLWAGTmK8AVwO3ARuBbmbmhQ/1KkhrUslNhw8nMncAFdeY/DcyvmV4DrBlhXVNK9ydJGj3/8l6SVJTBIkkqymCRJBVlsEiSijJYJElFGSySpKIMFklSUQaLJKkog0WSVJTBIkkqymCRJBVlsEiSijJYJElFGSySpKIMFklSUQaLJKkog0WSVJTBIkkqymCRJBVlsEiSijJYJElFGSySpKIMFklSUQaLJKmoyMxO99A2EbED+Pkoh58A/LJgO+1gz+1hz+1hz+1Rr+fXZeaJja7gkAqWZkTE+szs73QfB8Ke28Oe28Oe26NEz54KkyQVZbBIkooyWBq3otMNjII9t4c9t4c9t0fTPXuNRZJUlEcskqSiDBZJUlEGCxAR8yJiU0RsiYildZZHRHy5uvzhiJjV6Nhu6zkiJkfEXRGxMSI2RMQnu7nfmuVjIuLBiPheO/pttueIOCYiboiIn1Zf6zcdBD3/h+rPxKMR8c8RMb5Lev6jiPjXiHgxIq46kLHd1nOn9r9meq5Z3vg+mJmH9AMYA/wb8IfAYcBPgBlDauYDtwIBzAHubXRsF/bcA8yqPp8A/KzVPTfTb83yK4H/CXyv238uqstWAh+uPj8MOKabewZOAZ4Afr86/S3g0i7peSLw74G/B646kLFd2HPb979me65Z3vA+6BELzAa2ZObjmfkSsBpYOKRmIbAqK9YBx0RET4Nju6rnzBzMzAcAMvN5YCOV/1S6sl+AiOgF3gl8rcV9Fuk5Io4CzgOuA8jMlzLz2W7uubpsLPD7ETEWOAJ4uht6zsztmXkf8PKBju22nju0/zXVMxz4PmiwVP5Rt9ZMD7DvP/T+ahoZ2wrN9LxHREwB3gDcW77FA+tlhJqrgU8Dv2tRf/U00/MfAjuA66unDr4WEX/QymZH6GfEmszcBvwD8BQwCDyXmf/Swl6H7acNY5tRZLtt3P+g+Z6v5gD2QYOlckpgqKHvwd5fTSNjW6GZnisLI44Evg18KjN/U7C3ekbdb0S8C9iemfeXb2tYzbzGY4FZwFcz8w3A/wHacf6/mdf5WCq/wfYBk4A/iIi/KNxfPc3sQ928/w2/gvbuf9BEz6PZBw2WSnJPrpnuZd9TAPuraWRsKzTTMxExjsoP9Tcz8zst7HPEXhqoeTOwICKepHL4fn5E/FPrWh2xn0ZqBoCBzNz9m+gNVIKm1Zrp+R3AE5m5IzNfBr4DnNvCXkfqp9Vjm9HUdjuw/0FzPR/4Ptjqi0bd/qDy2+XjVH5T231Ra+aQmney9wXPHzc6tgt7DmAVcPXB8BoPqXkb7bt431TPwA+B06vP/w74Qjf3DJwDbKBybSWovPng493Qc03t37H3hfCu3f+G6bnt+1+zPQ9Z1tA+2LZvrJsfVN4p8zMq75r4THXeEmBJzQ/D8uryR4D+4cZ2c8/AW6gcAj8MPFR9zO/Wfoeso6Ef6m7oGTgbWF99nW8Cjj0Iev4s8FPgUeAbwOFd0vPJVH7j/g3wbPX5Ufsb2809d2r/a/Z1rllHQ/ugt3SRJBXlNRZJUlEGiySpKINFklSUwSJJKspgkSQVZbBIo1S9g/HlNdOTIuKGFm3r3RHxX0ao+YeIOL8V25cOhG83lkapeq+n72XmGW3Y1o+ABZn5y2FqXgdcm5l/3Op+pOF4xCKN3jLg9RHxUER8ISKmRMSjABFxaUTcFBHfjYgnIuKKiLiyelPKdRFxXLXu9RFxW0TcHxE/jIg/GrqRiDgNeDEzfxkRE6rrG1dddlREPBkR4zLz58DxEXFyG18DaR8GizR6S4F/y8yzM/M/1ll+BvDnVG5Z/vfAb7NyU8p/BS6p1qygcuuUNwJXAf9YZz1vBmpvtf4DKrdmAbgY+HZW7u9Fte7NTX5fUlPGdroB6TXsrmoQPB8RzwHfrc5/BDireofbc4H/FbHn5rOH11lPD5Xb8O/2NSq3ML8J+ADwkZpl26ncnVjqGINFap0Xa57/rmb6d1T2vd8Dns3Ms0dYz/8Fjt49kZn3VE+7vRUYk5mP1tSOr9ZLHeOpMGn0nqfy8bKjkpXP4XgiIv4M9nwe/b+rU7oRmDpk3irgn4Hrh8w/jcpNJKWOMVikUcrMncA9EfFoRHxhlKt5P/ChiPgJldvW1/to3buBN0TN+TLgm8CxVMIF2PM5H1Op3FVZ6hjfbiwdBCLiS8B3M/P71ek/BRZm5l/W1LwHmJWZ/7lDbUqA11ikg8XnqHwYFxHxP4ALqXy+Rq2xwBfb3Je0D49YJElFeY1FklSUwSJJKspgkSQVZbBIkooyWCRJRf0/eWhzqnV1OZoAAAAASUVORK5CYII=\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEGCAYAAABGnrPVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAUoElEQVR4nO3dfZBddZ3n8fdnkjAZxzA8DwkNJqNxJwFdjL2AaIGiTAE6RGW3lugoKCvLjli6rutm1pp1/GO2rHKcRVdGKz5MAeOacvFhwIoyCjhO4eAQBHkwRjKI0hAFsyUPyyKg3/3jXqymvUnf9P3de7vJ+1XVlXvO+Z1zPmn68Mk55/a5qSokSWrlN8YdQJL09GKxSJKaslgkSU1ZLJKkpiwWSVJTi8cdYJQOOeSQWrly5bhjSNKCcuONN/60qg7td/w+VSwrV65k69at444hSQtKkh/uzXgvhUmSmrJYJElNWSySpKb2qXsskjQOjz/+OFNTUzz66KPjjrJHS5cuZWJigiVLlgy0HYtFkoZsamqKZcuWsXLlSpKMO05PVcWuXbuYmppi1apVA23LS2GSNGSPPvooBx988LwtFYAkHHzwwU3OqiwWSRqB+VwqT2qV0WKRJDVlsUjSAnHiiSf2nH/uuedy+eWXjzjN7lkskrRAfPOb3xx3hL74rjBJWiCe+cxn8vDDD1NVvO1tb+Oaa65h1apVzLdPAvaMRZIWmC984Qts376dW2+9lY9//OPz7kzGYpGkBeYb3/gGGzZsYNGiRaxYsYJTTjll3JGewmKRpAVoPr992WKRpAXmpJNOYvPmzfziF79g586dXHvtteOO9BTevJekBeY1r3kN11xzDc973vN47nOfy8knnzzuSE9hsUjSAvHwww8DnctgH/nIR8acZve8FCZJaspikSQ1ZbFIkpqyWCRJTVkskqSmLBZJUlMWiyTtI9785jdz2GGHccwxxwx1PxaLJO0jzj33XL7yla8MfT9jLZYkpyXZnmRHko09lifJh7vLb0mybsbyRUluSvKl0aWWpIXppJNO4qCDDhr6fsb2m/dJFgEXA6cCU8ANSa6oqu9OG3Y6sLr7dTzw0e6fT3o7sA3YfyShJWlA77vydr5774NNt7l2xf689w+PbrrNQYzzjOU4YEdV3VlVjwGbgfUzxqwHLq2O64EDkiwHSDIBvBL4xChDS5L2bJzPCjsCuHva9BRPPRvZ3ZgjgJ3ARcC7gWV72kmS84HzAY466qiBAkvSoObTmcWwjPOMpdeHCcz8fM2eY5K8Crivqm6cbSdVtamqJqtq8tBDD51LTknSXhhnsUwBR06bngDu7XPMi4Ezk9xF5xLaKUn+ZnhRJWnh27BhAy960YvYvn07ExMTfPKTnxzKfsZ5KewGYHWSVcA9wNnA62aMuQK4MMlmOpfJHqiqncCfdL9I8lLgXVX1RyPKLUkL0mc+85mR7GdsxVJVTyS5ELgKWAR8qqpuT3JBd/nHgC3AGcAO4BHgTePKK0nqz1g/6KuqttApj+nzPjbtdQFvnWUbXwe+PoR4kqQ58DfvJUlNWSySpKYsFklSUxaLJKkpi0WS9gF33303L3vZy1izZg1HH300H/rQh4a2r7G+K0ySNBqLFy/mgx/8IOvWreOhhx7ihS98Iaeeeipr165tvi/PWCRpH7B8+XLWret88siyZctYs2YN99xzz1D25RmLJI3SlzfCj29tu83Dnwenv7/v4XfddRc33XQTxx8/87m/bXjGIkn7kIcffpizzjqLiy66iP33H85HWXnGIkmjtBdnFq09/vjjnHXWWbz+9a/nta997dD24xmLJO0DqorzzjuPNWvW8M53vnOo+7JYJGkfcN1113HZZZdxzTXXcOyxx3LssceyZcuW2VecAy+FSdI+4CUveQmd5/oOn2cskqSmLBZJUlMWiySNwKguQw2iVUaLRZKGbOnSpezatWtel0tVsWvXLpYuXTrwtrx5L0lDNjExwdTUFPfff/+4o+zR0qVLmZiYGHg7FoskDdmSJUtYtWrVuGOMjJfCJElNWSySpKYsFklSUxaLJKkpi0WS1JTFIklqymKRJDVlsUiSmrJYJElNWSySpKYsFklSU2MtliSnJdmeZEeSjT2WJ8mHu8tvSbKuO//IJNcm2Zbk9iRvH316SVIvYyuWJIuAi4HTgbXAhiRrZww7HVjd/Tof+Gh3/hPAf6qqNcAJwFt7rCtJGoNxnrEcB+yoqjur6jFgM7B+xpj1wKXVcT1wQJLlVbWzqr4NUFUPAduAI0YZXpLU2ziL5Qjg7mnTU/x6Ocw6JslK4AXAt9pHlCTtrXEWS3rMm/nxansck+SZwOeAd1TVgz13kpyfZGuSrfP9Q3Yk6elgnMUyBRw5bXoCuLffMUmW0CmVT1fV53e3k6raVFWTVTV56KGHNgkuSdq9cRbLDcDqJKuS7AecDVwxY8wVwBu77w47AXigqnYmCfBJYFtV/eVoY0uS9mRsH01cVU8kuRC4ClgEfKqqbk9yQXf5x4AtwBnADuAR4E3d1V8MvAG4NcnN3Xn/taq2jPCvIEnqIVUzb2s8fU1OTtbWrVvHHUOSFpQkN1bVZL/j/c17SVJTFoskqSmLRZLUlMUiSWrKYpEkNWWxSJKaslgkSU1ZLJKkpiwWSVJTFoskqSmLRZLUlMUiSWrKYpEkNWWxSJKaslgkSU1ZLJKkpiwWSVJTFoskqSmLRZLUlMUiSWrKYpEkNWWxSJKaslgkSU1ZLJKkpiwWSVJTFoskqSmLRZLUlMUiSWpq1mJJct6M6UVJ3ju8SJKkhayfM5aXJ9mSZHmSY4DrgWVDziVJWqAWzzagql6X5N8CtwKPABuq6rqhJ5MkLUj9XApbDbwd+BxwF/CGJM9osfMkpyXZnmRHko09lifJh7vLb0myrt91JUnj0c+lsCuBP62qfw+cDNwB3DDojpMsAi4GTgfWAhuSrJ0x7HRgdffrfOCje7GuJGkMZr0UBhxXVQ8CVFUBH0xyRYN9HwfsqKo7AZJsBtYD3502Zj1waXe/1yc5IMlyYGUf6zZz/V+9hWU/2zaMTUvSUP1w8bO55HcuYO2K/XnvHx49kn32c4/lwSQn0vmf+fTxdwy47yOAu6dNTwHH9zHmiD7XBSDJ+XTOdjjqqKMGSyxJmtWsxZLkMuDZwM3AL7qzC7h0wH2nx7zqc0w/63ZmVm0CNgFMTk72HDObE/7443NZTZLG7mjgjBHvs59LYZPA2u7lqJamgCOnTU8A9/Y5Zr8+1pUkjUE/N+9vAw4fwr5vAFYnWZVkP+BsYOa9myuAN3bfHXYC8EBV7exzXUnSGOz2jCXJlXQuLy0Dvpvkn4CfP7m8qs4cZMdV9USSC4GrgEXAp6rq9iQXdJd/DNhC5yxuB53foXnTntYdJI8kqY3s7gpXkpO7L18I3Af8iGn3Nqrq74eerrHJycnaunXruGNI0oKS5Maqmux3/G7PWJ4sjiQvBc4D/g+wGbi8qn4yWExJ0tPVrPdYqup9VXU08FZgBfD3Sb429GSSpAVpbx6bfx/wY2AXcNhw4kiSFrp+nhX2H5J8HbgaOAR4S1U9f9jBJEkLUz+/x/Is4B1VdfOQs0iSngb6eaSLTw6WJPXNjyaWJDVlsUiSmrJYJElNWSySpKYsFklSUxaLJKkpi0WS1JTFIklqymKRJDVlsUiSmrJYJElNWSySpKYsFklSUxaLJKkpi0WS1JTFIklqymKRJDVlsUiSmrJYJElNWSySpKYsFklSUxaLJKkpi0WS1JTFIklqaizFkuSgJF9Nckf3zwN3M+60JNuT7Eiycdr8DyT5XpJbknwhyQEjCy9J2qNxnbFsBK6uqtXA1d3pp0iyCLgYOB1YC2xIsra7+KvAMVX1fOD7wJ+MJLUkaVbjKpb1wCXd15cAr+4x5jhgR1XdWVWPAZu761FVf1dVT3THXQ9MDDeuJKlf4yqW362qnQDdPw/rMeYI4O5p01PdeTO9Gfhy84SSpDlZPKwNJ/kacHiPRe/pdxM95tWMfbwHeAL49B5ynA+cD3DUUUf1uWtJ0lwNrViq6hW7W5bkJ0mWV9XOJMuB+3oMmwKOnDY9Adw7bRvnAK8CXl5VxW5U1SZgE8Dk5ORux0mS2hjXpbArgHO6r88B/rbHmBuA1UlWJdkPOLu7HklOA/4LcGZVPTKCvJKkPo2rWN4PnJrkDuDU7jRJViTZAtC9OX8hcBWwDfhsVd3eXf8jwDLgq0luTvKxUf8FJEm9De1S2J5U1S7g5T3m3wucMW16C7Clx7jnDDWgJGnO/M17SVJTFoskqSmLRZLUlMUiSWrKYpEkNWWxSJKaslgkSU1ZLJKkpiwWSVJTFoskqSmLRZLUlMUiSWrKYpEkNWWxSJKaslgkSU1ZLJKkpiwWSVJTFoskqSmLRZLUlMUiSWrKYpEkNWWxSJKaslgkSU1ZLJKkpiwWSVJTFoskqSmLRZLUlMUiSWrKYpEkNWWxSJKaslgkSU2NpViSHJTkq0nu6P554G7GnZZke5IdSTb2WP6uJJXkkOGnliT1Y1xnLBuBq6tqNXB1d/opkiwCLgZOB9YCG5Ksnbb8SOBU4EcjSSxJ6su4imU9cEn39SXAq3uMOQ7YUVV3VtVjwObuek/6H8C7gRpiTknSXhpXsfxuVe0E6P55WI8xRwB3T5ue6s4jyZnAPVX1ndl2lOT8JFuTbL3//vsHTy5J2qPFw9pwkq8Bh/dY9J5+N9FjXiV5Rncbf9DPRqpqE7AJYHJy0rMbSRqyoRVLVb1id8uS/CTJ8qramWQ5cF+PYVPAkdOmJ4B7gWcDq4DvJHly/reTHFdVP272F5Akzcm4LoVdAZzTfX0O8Lc9xtwArE6yKsl+wNnAFVV1a1UdVlUrq2olnQJaZ6lI0vwwrmJ5P3BqkjvovLPr/QBJViTZAlBVTwAXAlcB24DPVtXtY8orSerT0C6F7UlV7QJe3mP+vcAZ06a3AFtm2dbK1vkkSXPnb95LkpqyWCRJTVkskqSmLBZJUlMWiySpKYtFktSUxSJJaspikSQ1ZbFIkpqyWCRJTVkskqSmLBZJUlMWiySpKYtFktSUxSJJaspikSQ1ZbFIkpqyWCRJTVkskqSmLBZJUlMWiySpKYtFktSUxSJJaspikSQ1laoad4aRSXI/8MM5rn4I8NOGcUbBzKNh5tEw82j0yvysqjq03w3sU8UyiCRbq2py3Dn2hplHw8yjYebRaJHZS2GSpKYsFklSUxZL/zaNO8AcmHk0zDwaZh6NgTN7j0WS1JRnLJKkpiwWSVJTFguQ5LQk25PsSLKxx/Ik+XB3+S1J1vW77nzLnOTIJNcm2Zbk9iRvn895py1flOSmJF8aRd5BMyc5IMnlSb7X/V6/aAFk/o/dn4nbknwmydJ5kvn3k/xjkp8nedferDvfMo/r+Bsk87Tl/R+DVbVPfwGLgH8Gfg/YD/gOsHbGmDOALwMBTgC+1e+68zDzcmBd9/Uy4PvDzjxI3mnL3wn8L+BL8/3norvsEuDfdV/vBxwwnzMDRwA/AH6rO/1Z4Nx5kvkw4F8Bfw68a2/WnYeZR378DZp52vK+j0HPWOA4YEdV3VlVjwGbgfUzxqwHLq2O64EDkizvc915lbmqdlbVtwGq6iFgG53/qczLvABJJoBXAp8Ycs4mmZPsD5wEfBKgqh6rqp/N58zdZYuB30qyGHgGcO98yFxV91XVDcDje7vufMs8puNvoMyw98egxdL5j3r3tOkpfv0/9O7G9LPuMAyS+VeSrAReAHyrfcS9yzLLmIuAdwO/HFK+XgbJ/HvA/cBfdy8dfCLJbw8z7Cx5Zh1TVfcAfwH8CNgJPFBVfzfErHvMM4J1B9FkvyM8/mDwzBexF8egxdK5JDDTzPdg725MP+sOwyCZOwuTZwKfA95RVQ82zNbLnPMmeRVwX1Xd2D7WHg3yPV4MrAM+WlUvAP4vMIrr/4N8nw+k8y/YVcAK4LeT/FHjfL0McgzN5+NvzxsY7fEHA2SeyzFosXSa+8hp0xP8+iWA3Y3pZ91hGCQzSZbQ+aH+dFV9fog5Z83Sx5gXA2cmuYvO6fspSf5meFFnzdPPmClgqqqe/Jfo5XSKZtgGyfwK4AdVdX9VPQ58HjhxiFlnyzPsdQcx0H7HcPzBYJn3/hgc9k2j+f5F51+Xd9L5l9qTN7WOnjHmlTz1huc/9bvuPMwc4FLgooXwPZ4x5qWM7ub9QJmBfwD+Rff1nwEfmM+ZgeOB2+ncWwmdNx+8bT5knjb2z3jqjfB5e/ztIfPIj79BM89Y1tcxOLK/2Hz+ovNOme/TedfEe7rzLgAumPbDcHF3+a3A5J7Wnc+ZgZfQOQW+Bbi5+3XGfM07Yxt9/VDPh8zAscDW7vf5i8CBCyDz+4DvAbcBlwG/OU8yH07nX9wPAj/rvt5/d+vO58zjOv4G/T5P20Zfx6CPdJEkNeU9FklSUxaLJKkpi0WS1JTFIklqymKRJDVlsUhz1H2C8R9Pm16R5PIh7evVSf7bLGP+Iskpw9i/tDd8u7E0R91nPX2pqo4Zwb6+CZxZVT/dw5hnAR+vqj8Ydh5pTzxjkebu/cCzk9yc5ANJVia5DSDJuUm+mOTKJD9IcmGSd3YfSnl9koO6456d5CtJbkzyD0l+f+ZOkjwX+HlV/TTJsu72lnSX7Z/kriRLquqHwMFJDh/h90D6NRaLNHcbgX+uqmOr6j/3WH4M8Do6jyz/c+CR6jyU8h+BN3bHbKLz6JQXAu8C/qrHdl4MTH/U+tfpPJoF4Gzgc9V5vhfdcS8e8O8lDWTxuANIT2PXdovgoSQPAFd2598KPL/7hNsTgf+d/Orhs7/ZYzvL6TyG/0mfoPMI8y8CbwLeMm3ZfXSeTiyNjcUiDc/Pp73+5bTpX9I59n4D+FlVHTvLdv4f8DtPTlTVdd3LbicDi6rqtmljl3bHS2PjpTBp7h6i8/Gyc1Kdz+H4QZJ/A7/6PPp/2WPoNuA5M+ZdCnwG+OsZ859L5yGS0thYLNIcVdUu4LoktyX5wBw383rgvCTfofPY+l4frfsN4AWZdr0M+DRwIJ1yAX71OR/PofNUZWlsfLuxtAAk+RBwZVV9rTv9r4H1VfWGaWNeA6yrqj8dU0wJ8B6LtFD8dzofxkWS/wmcTufzNaZbDHxwxLmkX+MZiySpKe+xSJKaslgkSU1ZLJKkpiwWSVJTFoskqan/D1xMo0J6AqrAAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] @@ -103,12 +103,12 @@ } ], "source": [ - "swiftdiff['vx'].plot.line(x=\"time (y)\")" + "swiftdiff['vhx'].plot.line(x=\"time (y)\")" ] }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 8, "metadata": {}, "outputs": [ { @@ -465,7 +465,7 @@ " stroke: currentColor;\n", " fill: currentColor;\n", "}\n", - "
<xarray.DataArray 'vx' (time (y): 199)>\n",
+       "
<xarray.DataArray 'vhx' (time (y): 199)>\n",
        "array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
        "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
        "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
@@ -483,8 +483,8 @@
        "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
        "        0.,  0.,  0., nan])\n",
        "Coordinates:\n",
-       "    id        float64 100.0\n",
-       "  * time (y)  (time (y)) float64 0.0 0.0006845 0.001369 ... 0.1342 0.1348 0.1355
    • " ], "text/plain": [ - "\n", + "\n", "array([ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", @@ -548,17 +548,17 @@ " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", " 0., 0., 0., nan])\n", "Coordinates:\n", - " id float64 100.0\n", + " id int64 2\n", " * time (y) (time (y)) float64 0.0 0.0006845 0.001369 ... 0.1342 0.1348 0.1355" ] }, - "execution_count": 7, + "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "swiftdiff['vx'].sel(id=100)" + "swiftdiff['vhx'].sel(id=2)" ] }, { diff --git a/examples/symba_swifter_comparison/1pl_1tp_encounter/tp.swifter.in b/examples/symba_swifter_comparison/1pl_1tp_encounter/tp.swifter.in index 9c026369e..ea062b86c 100644 --- a/examples/symba_swifter_comparison/1pl_1tp_encounter/tp.swifter.in +++ b/examples/symba_swifter_comparison/1pl_1tp_encounter/tp.swifter.in @@ -1,4 +1,4 @@ 1 -100 +2 1.01 0.0 0.0 0.0 6.252003053624663 0.0 diff --git a/examples/symba_swifter_comparison/1pl_1tp_encounter/tp.swiftest.in b/examples/symba_swifter_comparison/1pl_1tp_encounter/tp.swiftest.in index e1506974ae338f098f33af16f09e91ed31946bbc..15ed8c26290abce29999a075e27b5cec1080501e 100644 GIT binary patch literal 54 zcmXry3P~+42}mp|$xO~kIX>i_@% diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index 91525b821..6ef91c960 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -30,13 +30,17 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l associate(pl => self, plplenc_list => system%plplenc_list) if (param%ladaptive_interactions) then - if (lfirst) then - write(itimer%loopname, *) "symba_encounter_check_pl" - write(itimer%looptype, *) "ENCOUNTERS" - call itimer%time_this_loop(param, pl, pl%nplplm) - lfirst = .false. + if (self%nplplm > 0) then + if (lfirst) then + write(itimer%loopname, *) "symba_encounter_check_pl" + write(itimer%looptype, *) "ENCOUNTERS" + call itimer%time_this_loop(param, pl, pl%nplplm) + lfirst = .false. + else + if (itimer%check(param, pl%nplplm)) call itimer%time_this_loop(param, pl, pl%nplplm) + end if else - if (itimer%check(param, pl%nplplm)) call itimer%time_this_loop(param, pl, pl%nplplm) + param%lflatten_encounters = .false. end if end if @@ -96,7 +100,7 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l end do end if - if (param%ladaptive_interactions) then + if (param%ladaptive_interactions .and. self%nplplm > 0) then if (itimer%is_on) call itimer%adapt(param, pl, pl%nplplm) end if @@ -174,7 +178,7 @@ module function symba_encounter_check(self, param, system, dt, irec) result(lany j = self%index2(k) xr(:) = tp%xh(:,j) - pl%xh(:,i) vr(:) = tp%vb(:,j) - pl%vb(:,i) - call encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%rhill(i), dt, lencounter(lidx), self%lvdotr(k)) + call encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%renc(i), dt, lencounter(lidx), self%lvdotr(k)) if (lencounter(lidx)) then rlim2 = (pl%radius(i))**2 rji2 = dot_product(xr(:), xr(:))! Check to see if these are physically overlapping bodies first, which we should ignore diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index 529f00215..36043a4d2 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -18,13 +18,17 @@ module subroutine symba_kick_getacch_int_pl(self, param) logical, save :: lfirst = .true. if (param%ladaptive_interactions) then - if (lfirst) then - write(itimer%loopname, *) "symba_kick_getacch_int_pl" - write(itimer%looptype, *) "INTERACTION" - call itimer%time_this_loop(param, self, self%nplplm) - lfirst = .false. + if (self%nplplm > 0) then + if (lfirst) then + write(itimer%loopname, *) "symba_kick_getacch_int_pl" + write(itimer%looptype, *) "INTERACTION" + call itimer%time_this_loop(param, self, self%nplplm) + lfirst = .false. + else + if (itimer%check(param, self%nplplm)) call itimer%time_this_loop(param, self, self%nplplm) + end if else - if (itimer%check(param, self%nplplm)) call itimer%time_this_loop(param, self, self%nplplm) + param%lflatten_interactions = .false. end if end if @@ -34,7 +38,7 @@ module subroutine symba_kick_getacch_int_pl(self, param) call kick_getacch_int_all_triangular_pl(self%nbody, self%nplm, self%xh, self%Gmass, self%radius, self%ah) end if - if (param%ladaptive_interactions) then + if (param%ladaptive_interactions .and. self%nplplm > 0) then if (itimer%is_on) call itimer%adapt(param, self, self%nplplm) end if