From 36edbd384801eddb86df6d787475aaa7f6ff5b17 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 15 May 2023 15:45:50 -0400 Subject: [PATCH 1/3] Added a Swiftest vs. Swifter plotting script --- .../1pl_1tp_encounter/.gitignore | 2 + .../swiftest_vs_swifter.ipynb | 159 ++++++++++++++++++ .../1pl_1tp_encounter/swiftest_vs_swifter.py | 9 + 3 files changed, 170 insertions(+) create mode 100644 examples/rmvs_swifter_comparison/1pl_1tp_encounter/swiftest_vs_swifter.ipynb create mode 100644 examples/rmvs_swifter_comparison/1pl_1tp_encounter/swiftest_vs_swifter.py diff --git a/examples/rmvs_swifter_comparison/1pl_1tp_encounter/.gitignore b/examples/rmvs_swifter_comparison/1pl_1tp_encounter/.gitignore index 80c7e7bd6..89c2a1c6c 100644 --- a/examples/rmvs_swifter_comparison/1pl_1tp_encounter/.gitignore +++ b/examples/rmvs_swifter_comparison/1pl_1tp_encounter/.gitignore @@ -1,3 +1,5 @@ * !.gitignore !init_cond.py +!swiftest_vs_swifter.py +!swiftest_vs_swifter.ipynb diff --git a/examples/rmvs_swifter_comparison/1pl_1tp_encounter/swiftest_vs_swifter.ipynb b/examples/rmvs_swifter_comparison/1pl_1tp_encounter/swiftest_vs_swifter.ipynb new file mode 100644 index 000000000..5f5f7b2ba --- /dev/null +++ b/examples/rmvs_swifter_comparison/1pl_1tp_encounter/swiftest_vs_swifter.ipynb @@ -0,0 +1,159 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import swiftest" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Reading Swifter file /home/daminton/git_debug/swiftest/examples/rmvs_swifter_comparison/1pl_1tp_encounter/swifter_sim/param.in\n", + "Reading in time 1.355e-01\n", + "Creating Dataset\n", + "Successfully converted 199 output frames.\n", + "Swifter simulation data stored as xarray DataSet .data\n" + ] + } + ], + "source": [ + "swiftersim = swiftest.Simulation(simdir=\"swifter_sim\", read_data=True, codename=\"Swifter\")" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Reading Swiftest file /home/daminton/git_debug/swiftest/examples/rmvs_swifter_comparison/1pl_1tp_encounter/swiftest_sim/param.in\n", + "\n", + "Creating Dataset from NetCDF file\n", + "Successfully converted 1 output frames.\n", + "\n", + "Creating Dataset from NetCDF file\n", + "Successfully converted 221 output frames.\n", + "Swiftest simulation data stored as xarray DataSet .data\n", + "Reading initial conditions file as .init_cond\n", + "Finished reading Swiftest dataset files.\n" + ] + } + ], + "source": [ + "swiftestsim = swiftest.Simulation(simdir=\"swiftest_sim\",read_data=True)\n", + "swiftestsim.data = swiftestsim.data.swap_dims({\"name\" : \"id\"})" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "swiftdiff = swiftestsim.data - swiftersim.data" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAABBsAAAEiCAYAAAC4D9oGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAxNElEQVR4nO3de1RVdf7/8dfhLsrFS4IUInZRHLsJI2KR1jioZWXaN9MirbT8zhQptUyzlZcmHKvVmJqWDnmZb6nzzZyxGXNkKqxGROWLZsnPzPCWEGkF3gLEz++PhjMeQQTc5xw2PB9rnbU6n/PZe7/3aa838HJfHMYYIwAAAAAAAIv4eLsAAAAAAADQvBA2AAAAAAAASxE2AAAAAAAASxE2AAAAAAAASxE2AAAAAAAASxE2AAAAAAAASxE2AAAAAAAASxE2AAAAAAAASxE2AAAAAAAASxE2AAAAAAAASxE2AAAAAAAASxE2AAAAAAAASxE2oMV55513dPXVV6tVq1Zq3769BgwYoBMnTkiSxowZo6FDh2rGjBnq2LGjQkND9eijj6qiosK5/Pr163XjjTcqPDxc7du315AhQ7R3716XbRw6dEj33nuv2rVrp9atWyshIUG5ubnOz9977z3Fx8crKChIXbt21YwZM3T69Gm37K8xRgMGDNCgQYNkjJEk/fjjj+rcubOmTp3qlm0CsF5L61379u2Tj4+Ptm3b5jI+b948xcTEOPsZgKarpfUtSerSpYscDkeNF9ASETagRSkqKtLIkSP10EMPqaCgQNnZ2Ro2bJjLL60ffPCBCgoK9NFHH2nFihVas2aNZsyY4fz8xIkTSk9P19atW/XBBx/Ix8dHd911l86cOSNJOn78uPr166fDhw9r7dq12rFjhyZNmuT8/B//+Ifuv/9+paWladeuXXrjjTe0dOlSvfDCC+et+6233lKbNm3qfL311lu1LutwOLRs2TJt2bJFc+fOlSSNHz9eERERmj59+sV+pQA8oCX2ri5dumjAgAFasmSJy/iSJUs0ZswYfnkHmriW2LckaevWrSoqKlJRUZEOHTqkPn36KDk5+WK/TsCeDNCC5OXlGUlm3759tX4+evRo065dO3PixAnn2MKFC02bNm1MVVVVrcuUlJQYSWbnzp3GGGPeeOMNExISYo4ePVrr/OTkZJORkeEy9qc//cl06tTpvHWXlZWZPXv21PkqKyurc9///Oc/m8DAQDNlyhQTHBxsdu/eXed8AE1HS+1dq1atMm3btjU//fSTMcaY7du3G4fDYQoLC8+7DICmoaX2rbOlpaWZmJgYU1JSUq/5QHND2IAW5fTp0+ZXv/qVCQkJMXfffbdZtGiR+f77752fjx492tx8880uy2zfvt3lh+VXX31lRo4caWJjY01ISIhp3bq1kWT+/ve/G2OM+e///m9z0003nbeG4OBgExQUZFq3bu18BQUFGUkuP3DdYeTIkUaSWbhwoVu3A8BaLbV3lZeXm0suucSsWLHCGPPzL+633HKLW7YFwFottW9Ve+ONN0zr1q3N9u3b3bodoCnjMgq0KL6+vsrKytL777+vHj16aN68eerWrZsKCwsvuGz1Kbu33367jh49qsWLFys3N9d5XWD1NYatWrWqcz1nzpzRjBkztH37dudr586d2rNnj4KCgmpd5mJP6ZOkkydPKi8vT76+vtqzZ88F9xdA09FSe1dAQIBSU1O1ZMkSVVRU6O2339ZDDz10wX0G4H0ttW9JUnZ2th5//HEtX75c11577QX3F2iu/LxdAOBpDodDN9xwg2644QY999xziomJ0Zo1a5Seni5J2rFjh06dOuX8AbZ582a1adNGl112mY4ePaqCggK98cYbzuvvPv30U5f1X3PNNfrjH/+o77//Xu3ataux/V69emn37t264oor6l3zHXfcocTExDrnRERE1Pn5k08+KR8fH73//vu69dZbddttt+mWW26pdw0AvKul9q6xY8eqZ8+eWrBggSorKzVs2LB6bx+Ad7XEvvXVV19p+PDheuaZZ+hXaPEIG9Ci5Obm6oMPPlBKSoo6duyo3Nxcfffdd4qLi3POqaio0MMPP6xnn31W+/fv17Rp0/TYY4/Jx8dHbdu2Vfv27bVo0SJ16tRJBw4c0OTJk122MXLkSGVkZGjo0KGaNWuWOnXqpPz8fEVFRSkpKUnPPfechgwZoujoaP3Xf/2XfHx89Nlnn2nnzp363e9+V2vdISEhCgkJafR+//3vf9ebb76pnJwc9erVS5MnT9bo0aP12WefqW3bto1eLwDPaKm9S5Li4uLUp08fPf3003rooYcu+C+ZAJqGlti3Tp06pdtvv13XXXedHnnkERUXFzs/i4yMbNQ6AVvz9nUcgCft2rXLDBw40FxyySUmMDDQXHXVVWbevHnOz0ePHm3uvPNO89xzz5n27dubNm3amLFjxzpvTmaMMVlZWSYuLs4EBgaaa665xmRnZxtJZs2aNc45+/btM8OHDzehoaEmODjYJCQkmNzcXOfn69evN3379jWtWrUyoaGhpnfv3mbRokVu2eeSkhITERHhcoOkyspK07t3b3PPPfe4ZZsArNUSe9fZMjMzjSSzZcsWt28LgDVaYt8qLCw0kmp9AS2RwxgeVA1UGzNmjH788Uf95S9/8XYpAFBvzb13vfDCC1q5cqV27tzp7VIAWKS59y0AEjeIBAAATdLx48e1detWzZs3T2lpad4uBwAANABhAwAAaJIee+wx3XjjjerXrx9PoQAAwGa4jAIAAAAAAFiKMxsAAAAAAIClCBsAAAAAAIClCBsAAAAAAICl/LxdQHNw5swZHT58WCEhIXI4HN4uB0AzZIzRsWPHFBUVJR+fi8+J6VsA3I2+BcCOrO5dLRlhgwUOHz6s6Ohob5cBoAU4ePCgLrvssoteD30LgKfQtwDYkVW9qyUjbLBASEiIpJ8PyNDQUC9XA6A5KisrU3R0tLPfXCz6FgB3o28BsCOre1dLRthggepT+UJDQ/nhB8CtrDp1mL4FwFPoWwDsiMu1Lh4XoQAAAAAAAEsRNgAAAAAAAEsRNgAAAAAAAEsRNgAAAAAAAEsRNgAAAAAAAEsRNgAAAAAAAEsRNgAAAAAAAEsRNgAAAAAAAEsRNgAAAAAAAEsRNgAAAAAAAEsRNgAAAAAAAEsRNgAAAAAAAEsRNgAAAAAAAEsRNgAAAAAAAEsRNgAAAAAAAEsRNgAAAAAAAEsRNgAAAAAAAEsRNgAAAAAAAEsRNgAAAAAAAEsRNgAAAAAAAEsRNgAAAAAAAEvZLmxYsGCBYmNjFRQUpPj4eH3yySd1zt+4caPi4+MVFBSkrl276vXXXz/v3JUrV8rhcGjo0KEWVw0AAAAAQMthq7Bh1apVmjBhgqZOnar8/HwlJydr8ODBOnDgQK3zCwsLdeuttyo5OVn5+fl65plnlJaWptWrV9eYu3//fj311FNKTk52924AAAAAANCs2SpseOWVV/Twww9r7NixiouL05w5cxQdHa2FCxfWOv/1119X586dNWfOHMXFxWns2LF66KGH9PLLL7vMq6qq0n333acZM2aoa9euntgVAAAAAACaLduEDRUVFcrLy1NKSorLeEpKijZt2lTrMjk5OTXmDxw4UNu2bVNlZaVzbObMmbrkkkv08MMPW184AAAAAAAtjJ+3C6ivI0eOqKqqShERES7jERERKi4urnWZ4uLiWuefPn1aR44cUadOnfSvf/1LmZmZ2r59e71rKS8vV3l5ufN9WVlZ/XcEALyAvgXAbuhbAGBvtjmzoZrD4XB5b4ypMXah+dXjx44d0/3336/FixerQ4cO9a5h1qxZCgsLc76io6MbsAcA4Hn0LQB2Q98CAHuzTdjQoUMH+fr61jiLoaSkpMbZC9UiIyNrne/n56f27dtr79692rdvn26//Xb5+fnJz89Py5cv19q1a+Xn56e9e/fWut4pU6aotLTU+Tp48KA1OwkAbkLfAmA39C0AsDfbXEYREBCg+Ph4ZWVl6a677nKOZ2Vl6c4776x1maSkJL333nsuYxs2bFBCQoL8/f3VvXt37dy50+XzZ599VseOHdOrr7563gQ9MDBQgYGBF7lHAOA59C0AdkPfAgB7s03YIEnp6elKTU1VQkKCkpKStGjRIh04cEDjx4+X9HMC/s0332j58uWSpPHjx2v+/PlKT0/XuHHjlJOTo8zMTK1YsUKSFBQUpJ49e7psIzw8XJJqjAMAAAAAgPqxVdgwYsQIHT16VDNnzlRRUZF69uypdevWKSYmRpJUVFSkAwcOOOfHxsZq3bp1mjhxol577TVFRUVp7ty5Gj58uLd2AQAAAACAZs9hqu+YiEYrKytTWFiYSktLFRoa6u1yADRDVvcZ+hYAd6NvAbAjeo11bHODSAAAAAAAYA+EDQAAAAAAwFKEDQAAAAAAwFKEDQAAAAAAwFKEDQAAAAAAwFKEDQAAAAAAwFKEDQAAAAAAwFKEDQAAAAAAwFKEDQAAAAAAwFKEDQAAAAAAwFKEDQAAAAAAwFKEDQAAAAAAwFKEDQAAAAAAwFKEDQAAAAAAwFKEDQAAAAAAwFKEDQAAAAAAwFKEDQAAAAAAwFKEDQAAAAAAwFKEDQAAAAAAwFKEDQAAAAAAwFKEDQAAAAAAwFKEDQAAAAAAwFKEDQAAAAAAwFKEDQAAAAAAwFKEDQAAAAAAwFKEDQAAAAAAwFKEDQAAAAAAwFK2CxsWLFig2NhYBQUFKT4+Xp988kmd8zdu3Kj4+HgFBQWpa9euev31110+X7x4sZKTk9W2bVu1bdtWAwYM0JYtW9y5CwAAAAAANGu2ChtWrVqlCRMmaOrUqcrPz1dycrIGDx6sAwcO1Dq/sLBQt956q5KTk5Wfn69nnnlGaWlpWr16tXNOdna2Ro4cqY8++kg5OTnq3LmzUlJS9M0333hqtwAAAAAAaFYcxhjj7SLqKzExUb169dLChQudY3FxcRo6dKhmzZpVY/7TTz+ttWvXqqCgwDk2fvx47dixQzk5ObVuo6qqSm3bttX8+fP1wAMP1KuusrIyhYWFqbS0VKGhoQ3cKwC4MKv7DH0LgLvRtwDYEb3GOn7eLqC+KioqlJeXp8mTJ7uMp6SkaNOmTbUuk5OTo5SUFJexgQMHKjMzU5WVlfL396+xzMmTJ1VZWal27dqdt5by8nKVl5c735eVlTVkVwDA4+hbAOyGvgUA9mabyyiOHDmiqqoqRUREuIxHRESouLi41mWKi4trnX/69GkdOXKk1mUmT56sSy+9VAMGDDhvLbNmzVJYWJjzFR0d3cC9AQDPom8BsBv6FgDYm23ChmoOh8PlvTGmxtiF5tc2LkkvvviiVqxYoXfffVdBQUHnXeeUKVNUWlrqfB08eLAhuwAAHkffAmA39C0AsDfbXEbRoUMH+fr61jiLoaSkpMbZC9UiIyNrne/n56f27du7jL/88svKyMjQP//5T11zzTV11hIYGKjAwMBG7AUAeAd9C4Dd0LcAwN5sc2ZDQECA4uPjlZWV5TKelZWlvn371rpMUlJSjfkbNmxQQkKCy/0aXnrpJT3//PNav369EhISrC8eAAAAAIAWxDZhgySlp6frj3/8o958800VFBRo4sSJOnDggMaPHy/p59Ptzn6CxPjx47V//36lp6eroKBAb775pjIzM/XUU08557z44ot69tln9eabb6pLly4qLi5WcXGxjh8/7vH9AwAAAACgObDNZRSSNGLECB09elQzZ85UUVGRevbsqXXr1ikmJkaSVFRUpAMHDjjnx8bGat26dZo4caJee+01RUVFae7cuRo+fLhzzoIFC1RRUaG7777bZVvTpk3T9OnTPbJfAAAAAAA0Jw5TfcdENBrPYgXgbjyvHoDd0LcA2BG9xjq2uowCAAAAAAA0fYQNAAAAAADAUoQNAAAAAADAUoQNAAAAAADAUoQNAAAAAADAUoQNAAAAAADAUoQNAAAAAADAUoQNAAAAAADAUoQNAAAAAADAUn7eLgAAAAAAgObGGKPTp0+rqqrK26W4ha+vr/z8/ORwOGr9nLABAAAAAAALVVRUqKioSCdPnvR2KW4VHBysTp06KSAgoMZnhA0AAAAAAFjkzJkzKiwslK+vr6KiohQQEHDef/23K2OMKioq9N1336mwsFBXXnmlfHxc79JA2AAAAAAAgEUqKip05swZRUdHKzg42NvluE2rVq3k7++v/fv3q6KiQkFBQS6fc4NIAAAAAAAsdu6/9DdHde1j8997AAAAAADgUYQNAAAAAADYSP/+/TVhwoTzft6lSxfNmTPHY/XUhns2AAAAAABgI++++678/f29XUadCBsAAAAAALCRdu3aebuEC+IyCgAAAAAAbOTsyyhKSkp0++23q1WrVoqNjdVbb73l3eL+jTMbAAAAAACwqTFjxujgwYP68MMPFRAQoLS0NJWUlHi7LMIGAAAAAADs6Msvv9T777+vzZs3KzExUZKUmZmpuLg4L1fGZRQAAAAAANhSQUGB/Pz8lJCQ4Bzr3r27wsPDvVfUvxE2AAAAAABgQ8YYSZLD4fByJTURNgAAAAAAYENxcXE6ffq0tm3b5hzbvXu3fvzxR+8V9W+EDQAAAAAA2FC3bt00aNAgjRs3Trm5ucrLy9PYsWPVqlUrb5dG2AAAAAAAgF0tWbJE0dHR6tevn4YNG6ZHHnlEHTt29HZZPI0CAAAAAAA7yc7Odv53ZGSk/va3v7l8npqa6uGKarqoMxsqKip06NAhHThwwOXlTgsWLFBsbKyCgoIUHx+vTz75pM75GzduVHx8vIKCgtS1a1e9/vrrNeasXr1aPXr0UGBgoHr06KE1a9a4q3wAAAAAAJq9RoUNe/bsUXJyslq1aqWYmBjFxsYqNjZWXbp0UWxsrNU1Oq1atUoTJkzQ1KlTlZ+fr+TkZA0ePPi8AUdhYaFuvfVWJScnKz8/X88884zS0tK0evVq55ycnByNGDFCqamp2rFjh1JTU3XPPfcoNzfXbfsBAAAAAEBz5jDVz8pogBtuuEF+fn6aPHmyOnXqVOMxG9dee61lBZ4tMTFRvXr10sKFC51jcXFxGjp0qGbNmlVj/tNPP621a9eqoKDAOTZ+/Hjt2LFDOTk5kqQRI0aorKxM77//vnPOoEGD1LZtW61YsaJedZWVlSksLEylpaUKDQ294Hxz5oxOnTxWr3UDaKb8g9UqwK/ejylqaJ+xen30LQDyD5YcDrXy961X76JvAWgSvPA7108//aTCwkLnGfnNWV372qh7Nmzfvl15eXnq3r27JQXWR0VFhfLy8jR58mSX8ZSUFG3atKnWZXJycpSSkuIyNnDgQGVmZqqyslL+/v7KycnRxIkTa8yZM2fOeWspLy9XeXm5831ZWVmD9uXUyWMKfrlzg5YB0LzE/fSm8mbeqeAAz9w6h74F4GLF/fSmTilIu2YO9Ejvom8BsIKnf+fCfzTqMooePXroyJEjVtdSpyNHjqiqqkoREREu4xERESouLq51meLi4lrnnz592ln/+eacb52SNGvWLIWFhTlf0dHRjdklAPAY+hYAu6FvAYC91TveOTtNnj17tiZNmqSMjAxdffXV8vf3d5lrxaly53Pu6S/GmDpPialt/rnjDV3nlClTlJ6e7nxfVlbWoB+ArYJDdPIp995IE0DTlucfrFb+vh7bHn0LwMXKO+syCk+gbwGwgqd/58J/1DtsCA8Pd/kD3BijX/3qVy5zqv9Ir6qqsq7Cf+vQoYN8fX1rnHFQUlJS48yEapGRkbXO9/PzU/v27eucc751SlJgYKACAwMbsxuSJIePj4LbhDV6eQBoKPoWALuhbwGAvdU7bPjoo48k/Xz9XEZGhkaOHOnRezYEBAQoPj5eWVlZuuuuu5zjWVlZuvPOO2tdJikpSe+9957L2IYNG5SQkOA8GyMpKUlZWVku923YsGGD+vbt64a9AAAAAACg+at32NCvXz/nf99333265ZZbdOWVV7qlqPNJT09XamqqEhISlJSUpEWLFunAgQMaP368pJ9Pt/vmm2+0fPlyST8/eWL+/PlKT0/XuHHjlJOTo8zMTJenTDzxxBO66aabNHv2bN15553661//qn/+85/69NNPPbpvAAAAAAA0F426QeQDDzygzMxMq2u5oBEjRmjOnDmaOXOmrrvuOn388cdat26dYmJiJElFRUU6cOA/1+bFxsZq3bp1ys7O1nXXXafnn39ec+fO1fDhw51z+vbtq5UrV2rJkiW65pprtHTpUq1atUqJiYke3z8AAAAAAJoDh6m+Y2IDPP7441q+fLmuuOIKJSQkqHXr1i6fv/LKK5YVaAdWP0caAM7l7efVA0BD0bcA2JEVveann35SYWGhYmNjFRQUZHGF7vXxxx/rpZdeUl5enoqKirRmzRoNHTr0vPPr2tdGPWz0888/V69evSRJX375pctndT3FAQAAAAAANE0nTpzQtddeqwcffNDlioDGaFTYUH2zSAAAAAAA0DwMHjxYgwcPtmRdjQobAAAAAABA/RhjdKqyyuPbbeXv67WrDwgbAAAAAABwo1OVVerx3D88vt1dMwcqOMA7f/Y36mkUAAAAAAAA58OZDQAAAAAAuFErf1/tmjnQK9v1FsIGAAAAAADcyOFweO1yBm9pWXsLAAAAAABqdfz4cX311VfO94WFhdq+fbvatWunzp07N2hdhA0AAAAAAEDbtm3TzTff7Hyfnp4uSRo9erSWLl3aoHURNgAAAAAAAPXv31/GGEvWxdMoAAAAAACApQgbAAAAAACApQgbAAAAAACApQgbAAAAAACApQgbAAAAAACApQgbAAAAAACApQgbAAAAAACApQgbAAAAAACApQgbAAAAAACApQgbAAAAAACApQgbAAAAAABo4WbNmqVf/vKXCgkJUceOHTV06FDt3r270esjbAAAAAAAoIXbuHGjfvvb32rz5s3KysrS6dOnlZKSohMnTjRqfX4W1wcAAAAAAGxm/fr1Lu+XLFmijh07Ki8vTzfddFOD10fYAAAAAACAOxkjVZ70/Hb9gyWHo1GLlpaWSpLatWvXqOUJGwAAAAAAcKfKk1JGlOe3+8xhKaB1gxczxig9PV033nijevbs2ahNEzYAAAAAAACnxx57TJ999pk+/fTTRq/DNmHDDz/8oLS0NK1du1aSdMcdd2jevHkKDw8/7zLGGM2YMUOLFi3SDz/8oMTERL322mv6xS9+IUn6/vvvNW3aNG3YsEEHDx5Uhw4dNHToUD3//PMKCwvzxG4BAAAAAJo7/+CfzzLwxnYb6PHHH9fatWv18ccf67LLLmv0pm0TNowaNUqHDh1y3rTikUceUWpqqt57773zLvPiiy/qlVde0dKlS3XVVVfpd7/7nX79619r9+7dCgkJ0eHDh3X48GG9/PLL6tGjh/bv36/x48fr8OHDeueddzy1awAAAACA5szhaNTlDJ5kjNHjjz+uNWvWKDs7W7GxsRe1PluEDQUFBVq/fr02b96sxMRESdLixYuVlJSk3bt3q1u3bjWWMcZozpw5mjp1qoYNGyZJWrZsmSIiIvT222/r0UcfVc+ePbV69WrnMpdffrleeOEF3X///Tp9+rT8/Gzx9QAAAAAAcFF++9vf6u2339Zf//pXhYSEqLi4WJIUFhamVq1aNXh9PlYX6A45OTkKCwtzBg2S1KdPH4WFhWnTpk21LlNYWKji4mKlpKQ4xwIDA9WvX7/zLiP9fMfN0NBQggYAAAAAQIuxcOFClZaWqn///urUqZPztWrVqkatzxZ/URcXF6tjx441xjt27OhMW2pbRpIiIiJcxiMiIrR///5alzl69Kief/55Pfroo3XWU15ervLycuf7srKyOucDgLfRtwDYDX0LADzLGGPp+rx6ZsP06dPlcDjqfG3btk2S5Kjl2aDGmFrHz3bu5+dbpqysTLfddpt69OihadOm1bnOWbNmKSwszPmKjo6+0K4CgFfRtwDYDX0LAOzNq2HDY489poKCgjpfPXv2VGRkpL799tsay3/33Xc1zlyoFhkZKUk1znwoKSmpscyxY8c0aNAgtWnTRmvWrJG/v3+ddU+ZMkWlpaXO18GDBxuy2wDgcfQtAHZD3wIAe/PqZRQdOnRQhw4dLjgvKSlJpaWl2rJli3r37i1Jys3NVWlpqfr27VvrMrGxsYqMjFRWVpauv/56SVJFRYU2btyo2bNnO+eVlZVp4MCBCgwM1Nq1axUUFHTBegIDAxUYGFifXQSAJoG+BcBu6FsAYG+2uEFkXFycBg0apHHjxmnz5s3avHmzxo0bpyFDhrg8iaJ79+5as2aNpJ8vn5gwYYIyMjK0Zs0aff755xozZoyCg4M1atQoST+f0ZCSkqITJ04oMzNTZWVlKi4uVnFxsaqqqryyrwAAAAAA2J0tbhApSW+99ZbS0tKcT5e44447NH/+fJc5u3fvVmlpqfP9pEmTdOrUKf3mN7/RDz/8oMTERG3YsEEhISGSpLy8POXm5kqSrrjiCpd1FRYWqkuXLm7cIwAAAAAAmifbhA3t2rXT//zP/9Q559y7ZzocDk2fPl3Tp0+vdX7//v0tv+MmAAAAAAAt4W/NuvbRFpdRAAAAAABgB9UPHDh58qSXK3G/6n2s7SELtjmzAQAAAACAps7X11fh4eEqKSmRJAUHB8vhcHi5KmsZY3Ty5EmVlJQoPDxcvr6+NeYQNgAAAAAAYKHIyEhJcgYOzVV4eLhzX89F2AAAAAAAgIUcDoc6deqkjh07qrKy0tvluIW/v3+tZzRUI2wAAAAAAMANfH196/yDvDnjBpEAAAAAAMBShA0AAAAAAMBShA0AAAAAAMBShA0AAAAAAMBShA0AAAAAAMBShA0AAAAAAMBShA0AAAAAAMBShA0AAAAAAMBShA0AAAAAAMBShA0AAAAAAMBShA0AAAAAAMBShA0AAAAAAMBShA0AAAAAAMBShA0AAAAAAMBShA0AAAAAAMBShA0AAAAAAMBShA0AAAAAAMBShA0AAAAAAMBShA0AAAAAAMBShA0AAAAAAMBShA0AAAAAAMBShA0AAAAAAMBStgkbfvjhB6WmpiosLExhYWFKTU3Vjz/+WOcyxhhNnz5dUVFRatWqlfr3768vvvjivHMHDx4sh8Ohv/zlL9bvAAAAAAAALYRtwoZRo0Zp+/btWr9+vdavX6/t27crNTW1zmVefPFFvfLKK5o/f762bt2qyMhI/frXv9axY8dqzJ0zZ44cDoe7ygcAAAAAoMXw83YB9VFQUKD169dr8+bNSkxMlCQtXrxYSUlJ2r17t7p161ZjGWOM5syZo6lTp2rYsGGSpGXLlikiIkJvv/22Hn30UefcHTt26JVXXtHWrVvVqVMnz+wUAAAAAADNlC3ObMjJyVFYWJgzaJCkPn36KCwsTJs2bap1mcLCQhUXFyslJcU5FhgYqH79+rksc/LkSY0cOVLz589XZGRkveopLy9XWVmZywsAmjL6FgC7oW8BgL3ZImwoLi5Wx44da4x37NhRxcXF511GkiIiIlzGIyIiXJaZOHGi+vbtqzvvvLPe9cyaNct574iwsDBFR0fXe1kA8Ab6FgC7oW8BgL15NWyYPn26HA5Hna9t27ZJUq33UzDGXPA+C+d+fvYya9eu1Ycffqg5c+Y0qO4pU6aotLTU+Tp48GCDlgcAT6NvAbAb+hYA2JtX79nw2GOP6d57761zTpcuXfTZZ5/p22+/rfHZd999V+PMhWrVl0QUFxe73IehpKTEucyHH36ovXv3Kjw83GXZ4cOHKzk5WdnZ2bWuOzAwUIGBgXXWDQBNCX0LgN3QtwDA3rwaNnTo0EEdOnS44LykpCSVlpZqy5Yt6t27tyQpNzdXpaWl6tu3b63LxMbGKjIyUllZWbr++uslSRUVFdq4caNmz54tSZo8ebLGjh3rstzVV1+tP/zhD7r99tsvZtcAAAAAAGixbPE0iri4OA0aNEjjxo3TG2+8IUl65JFHNGTIEJcnUXTv3l2zZs3SXXfdJYfDoQkTJigjI0NXXnmlrrzySmVkZCg4OFijRo2S9PPZD7XdFLJz586KjY31zM4BAAAAANDM2CJskKS33npLaWlpzqdL3HHHHZo/f77LnN27d6u0tNT5ftKkSTp16pR+85vf6IcfflBiYqI2bNigkJAQj9YOAAAAAEBL4jDGGG8XYXdlZWUKCwtTaWmpQkNDvV0OgGbI6j5D3wLgbvQtAHZEr7GOLR59CQAAAAAA7IOwAQAAAAAAWIqwAQAAAAAAWIqwAQAAAAAAWIqwAQAAAAAAWIqwAQAAAAAAWIqwAQAAAAAAWIqwAQAAAAAAWIqwAQAAAAAAWIqwAQAAAAAAWIqwAQAAAAAAWIqwAQAAAAAAWIqwAQAAAAAAWIqwAQAAAAAAWIqwAQAAAAAAWIqwAQAAAAAAWIqwAQAAAAAAWIqwAQAAAAAAWIqwAQAAAAAAWIqwAQAAAAAAWIqwAQAAAAAAWMrP2wU0B8YYSVJZWZmXKwHQXFX3l+p+c7HoWwDcjb4FwI6s7l0tGWGDBY4dOyZJio6O9nIlAJq7Y8eOKSwszJL1SPQtAO5H3wJgR1b1rpbMYYhsLtqZM2d0+PBhhYSEyOFw1GuZsrIyRUdH6+DBgwoNDXVzhdaidu+xc/3UfnGMMTp27JiioqLk43PxV8DRt+zFzvVTu/d4u3761sWzc/3U7j12rr8p1G5172rJOLPBAj4+PrrssssatWxoaKjtmkA1avceO9dP7Y1nZbpO37InO9dP7d7jzfrpW9awc/3U7j12rt/btXNGgzWIagAAAAAAgKUIGwAAAAAAgKUIG7wkMDBQ06ZNU2BgoLdLaTBq9x4710/t9mfn78HOtUv2rp/avcfu9VvB7t+Bneundu+xc/12rh01cYNIAAAAAABgKc5sAAAAAAAAliJsAAAAAAAAliJsAAAAAAAAliJsaKQFCxYoNjZWQUFBio+P1yeffFLn/I0bNyo+Pl5BQUHq2rWrXn/99RpzVq9erR49eigwMFA9evTQmjVrLnq7nqh98eLFSk5OVtu2bdW2bVsNGDBAW7ZscZkzffp0ORwOl1dkZKTXa1+6dGmNuhwOh3766aeL2q6n6u/fv3+t9d92223OOd747ouKijRq1Ch169ZNPj4+mjBhQq3zmuIxX5/aPXnMW4m+9R+e/n9o595F33LlqWPeHfXbsXfRt/6DvuW92j3Ztxpaf1PrXfQtuDBosJUrVxp/f3+zePFis2vXLvPEE0+Y1q1bm/3799c6/+uvvzbBwcHmiSeeMLt27TKLFy82/v7+5p133nHO2bRpk/H19TUZGRmmoKDAZGRkGD8/P7N58+ZGb9dTtY8aNcq89tprJj8/3xQUFJgHH3zQhIWFmUOHDjnnTJs2zfziF78wRUVFzldJSUm963ZX7UuWLDGhoaEudRUVFV3Udj1Z/9GjR13q/vzzz42vr69ZsmSJc443vvvCwkKTlpZmli1bZq677jrzxBNP1JjTVI/5+tTuqWPeSvQt7/Qtd9Xvqd5F33LlqWPeXfXbrXfRt+hbTeW791Tfakz9Tal30bdwLsKGRujdu7cZP368y1j37t3N5MmTa50/adIk0717d5exRx991PTp08f5/p577jGDBg1ymTNw4EBz7733Nnq7nqr9XKdPnzYhISFm2bJlzrFp06aZa6+9tt511sYdtS9ZssSEhYVZul2r1tOY7/4Pf/iDCQkJMcePH3eOeeO7P1u/fv1q/eHRVI/5+tR+Lncd81aib3mnbxlj795F33LlqWP+YtfTXHoXfYu+daHteqr2c7mrbxlj795F38K5uIyigSoqKpSXl6eUlBSX8ZSUFG3atKnWZXJycmrMHzhwoLZt26bKyso651SvszHb9VTt5zp58qQqKyvVrl07l/E9e/YoKipKsbGxuvfee/X111/Xq2531378+HHFxMTosssu05AhQ5Sfn39R2/V0/WfLzMzUvffeq9atW7uMe/q7r4+mesw3hjuOeSvRt7zTt9xdv7t7F32rJk8c81au50Kacu+ib9G36rNdT9Z+Nnf0rcbWXx9N9bhvjKbct1ATYUMDHTlyRFVVVYqIiHAZj4iIUHFxca3LFBcX1zr/9OnTOnLkSJ1zqtfZmO16qvZzTZ48WZdeeqkGDBjgHEtMTNTy5cv1j3/8Q4sXL1ZxcbH69u2ro0ePerX27t27a+nSpVq7dq1WrFihoKAg3XDDDdqzZ0+jt+vJ+s+2ZcsWff755xo7dqzLuDe++/poqsd8Y7jjmLcSfcs7fcud9Xuid9G3avLEMW/lei6kKfcu+hZ9qz7b9VTtZ3NX32ps/fXRVI/7xmjKfQs1+Xm7ALtyOBwu740xNcYuNP/c8fqss6HbrW8tF1t7tRdffFErVqxQdna2goKCnOODBw92/vfVV1+tpKQkXX755Vq2bJnS09O9VnufPn3Up08f5+c33HCDevXqpXnz5mnu3LmN3q6n6j9bZmamevbsqd69e7uMe+u7t2qd3jjmG8Ldx7yV6Fve6VvuqN+TvYu+1fB1eutnRkPYpXfRt+hbTe27d3ffakz9Vq2T37lgNc5saKAOHTrI19e3RkJXUlJSI8mrFhkZWet8Pz8/tW/fvs451etszHY9VXu1l19+WRkZGdqwYYOuueaaOmtp3bq1rr76amea7e3aq/n4+OiXv/ylsy4rvndP1H/y5EmtXLmyRspeG0989/XRVI/5hnDnMW8l+pZ3+pYn6q/mjt5F36rJE8e8les5Hzv0LvoWfas+2/V07e7sW42tvz6a6nHfEHboW6iJsKGBAgICFB8fr6ysLJfxrKws9e3bt9ZlkpKSaszfsGGDEhIS5O/vX+ec6nU2Zrueql2SXnrpJT3//PNav369EhISLlhLeXm5CgoK1KlTJ6/XfjZjjLZv3+6sy4rv3RP1//nPf1Z5ebnuv//+C9biie++PprqMV9f7j7mrUTf8k7fcnf9Z3NH76Jv1eSJY97K9dTGLr2LvkXfqs92PV27O/tWY+uvj6Z63NeXXfoWamH9PSebv+rHumRmZppdu3aZCRMmmNatW5t9+/YZY4yZPHmySU1Ndc6vfqTOxIkTza5du0xmZmaNR+r861//Mr6+vub3v/+9KSgoML///e/P+0ia823XW7XPnj3bBAQEmHfeecflkTPHjh1zznnyySdNdna2+frrr83mzZvNkCFDTEhIiNdrnz59ulm/fr3Zu3evyc/PNw8++KDx8/Mzubm5ln7v7qq/2o033mhGjBhR63a98d0bY0x+fr7Jz8838fHxZtSoUSY/P9988cUXzs+b6jFfn9o9dcxbib7lnb7lrvo91bvoW97pW+6q3269i75F32oq3301d/etxtRvTNPpXfQtnIuwoZFee+01ExMTYwICAkyvXr3Mxo0bnZ+NHj3a9OvXz2V+dna2uf76601AQIDp0qWLWbhwYY11/u///q/p1q2b8ff3N927dzerV69u0Ha9VXtMTIyRVOM1bdo055wRI0aYTp06GX9/fxMVFWWGDRvm0ki8VfuECRNM586dTUBAgLnkkktMSkqK2bRpU4O26836jTFm9+7dRpLZsGFDrdv01ndf2zERExPjMqepHvMXqt2Tx7yV6Fv/4en/h3buXfStGJc5njrm3VG/HXsXfes/6Fveq90Yz/WtxtTflHoXfQtncxjz7zugAAAAAAAAWIB7NgAAAAAAAEsRNgAAAAAAAEsRNgAAAAAAAEsRNgAAAAAAAEsRNgAAAAAAAEsRNgAAAAAAAEsRNgAAAAAAAEsRNgAAAAAAAEsRNgAWy87OlsPh0I8//ujtUgCgXuhbAOyGvgU0fQ5jjPF2EYCd9e/fX9ddd53mzJkjSaqoqND333+viIgIORwO7xYHALWgbwGwG/oWYD9+3i4AaG4CAgIUGRnp7TIAoN7oWwDshr4FNH1cRgFchDFjxmjjxo169dVX5XA45HA4tHTpUpfT+pYuXarw8HD97W9/U7du3RQcHKy7775bJ06c0LJly9SlSxe1bdtWjz/+uKqqqpzrrqio0KRJk3TppZeqdevWSkxMVHZ2tnd2FECzQd8CYDf0LcCeOLMBuAivvvqqvvzyS/Xs2VMzZ86UJH3xxRc15p08eVJz587VypUrdezYMQ0bNkzDhg1TeHi41q1bp6+//lrDhw/XjTfeqBEjRkiSHnzwQe3bt08rV65UVFSU1qxZo0GDBmnnzp268sorPbqfAJoP+hYAu6FvAfZE2ABchLCwMAUEBCg4ONh5Kt//+3//r8a8yspKLVy4UJdffrkk6e6779af/vQnffvtt2rTpo169Oihm2++WR999JFGjBihvXv3asWKFTp06JCioqIkSU899ZTWr1+vJUuWKCMjw3M7CaBZoW8BsBv6FmBPhA2ABwQHBzt/8ElSRESEunTpojZt2riMlZSUSJL+7//+T8YYXXXVVS7rKS8vV/v27T1TNIAWjb4FwG7oW0DTQtgAeIC/v7/Le4fDUevYmTNnJElnzpyRr6+v8vLy5Ovr6zLv7B+YAOAu9C0AdkPfApoWwgbgIgUEBLjcaMgK119/vaqqqlRSUqLk5GRL1w0A9C0AdkPfAuyHp1EAF6lLly7Kzc3Vvn37dOTIEWdafjGuuuoq3XfffXrggQf07rvvqrCwUFu3btXs2bO1bt06C6oG0JLRtwDYDX0LsB/CBuAiPfXUU/L19VWPHj10ySWX6MCBA5asd8mSJXrggQf05JNPqlu3brrjjjuUm5ur6OhoS9YPoOWibwGwG/oWYD8OY4zxdhEAAAAAAKD54MwGAAAAAABgKcIGAAAAAABgKcIGAAAAAABgKcIGAAAAAABgKcIGAAAAAABgKcIGAAAAAABgKcIGAAAAAABgKcIGAAAAAABgKcIGAAAAAABgKcIGAAAAAABgKcIGAAAAAABgKcIGAAAAAABgqf8PkSwGwQz1DMEAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "swiftdiff['rh'].plot(x=\"time\",hue=\"id\",col=\"space\")" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAABBwAAAEiCAYAAABa08F/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA4oUlEQVR4nO3de3hU1bnH8d/kSgJhIEACKSHBHiFQUBEEglCwcgIqoqIVxEbAAvVYRMFLRVsuWol6rAevaBEBlYpWmxZvUY41tBYSLhKrkhMQuQoxgJJw0SQk6/xhM2XIhZnJntvO9/M888jsWXuvdybjC+vNWms7jDFGAAAAAAAAFooIdgAAAAAAAMB+KDgAAAAAAADLUXAAAAAAAACWo+AAAAAAAAAsR8EBAAAAAABYjoIDAAAAAACwHAUHAAAAAABgOQoOAAAAAADAchQcAAAAAACA5Sg4AAAAAAAAy1FwAAAAAAAAlqPgAAAAAAAALEfBAbb02muvqW/fvoqLi1OHDh00cuRIHT9+XJI0efJkXXnllVqwYIGSkpLUtm1b/eIXv1BVVZXr/Ly8PA0dOlTt2rVThw4dNGbMGO3YscOtj3379mnChAlKTExU69atNWDAABUWFrpef+ONN9S/f3+1atVKZ511lhYsWKCTJ0/65f0aYzRy5EiNHj1axhhJ0pEjR9StWzfde++9fukTgP+0tBy2a9cuRUREaNOmTW7Hn3jiCaWlpbnyGoDw0dLymCSlp6fL4XDUewAtGQUH2M6BAwd03XXX6cYbb1RxcbHy8/M1btw4t3+wvv/++youLtYHH3ygl19+Wbm5uVqwYIHr9ePHj2v27NnauHGj3n//fUVEROiqq65SbW2tJOnYsWMaPny49u/fr9WrV+vjjz/WXXfd5Xr93Xff1c9+9jPNnDlTW7du1bPPPqvly5frgQceaDTulStXqk2bNk0+Vq5c2eC5DodDK1as0IYNG/T4449Lkm666SYlJydr/vz5zf1IAQRQS8xh6enpGjlypJYtW+Z2fNmyZZo8eTL/YAfCTEvMY5K0ceNGHThwQAcOHNC+ffs0ePBgDRs2rLkfJxDeDDy2du1aM2bMGNOlSxcjyeTm5oZEf1u3bjWXX365adu2rWnTpo0ZNGiQ2b17t19jC2WbN282ksyuXbsafH3SpEkmMTHRHD9+3HVs8eLFpk2bNqampqbBc8rKyowk88knnxhjjHn22WdNQkKCOXz4cIPthw0bZhYuXOh27MUXXzRdunRpNO6Kigqzffv2Jh8VFRVNvvdXX33VxMbGmjlz5pj4+HhTUlLSZHsAoael5rBXXnnFtG/f3nz33XfGGGOKioqMw+EwO3fubPQcAKGppeaxU82cOdOkpaWZsrIyj9oDdhUVvFJH+Dl+/LjOPfdcTZkyRVdffXVI9Ldjxw4NHTpUP//5z7VgwQI5nU4VFxerVatWfo8vVJ177rm6+OKL1bdvX40aNUpZWVm65ppr1L59e7c28fHxrueZmZk6duyY9u7dq7S0NO3YsUO/+c1vVFBQoEOHDrmq5Xv27FGfPn1UVFSkfv36KTExscEYNm/erI0bN7pV0WtqavTdd9/pxIkTbn3XSUhIUEJCQrPe+09/+lPl5uYqJydHixcvVo8ePZp1PQCB11Jz2JVXXqkZM2YoNzdXEyZM0PPPP6+LLrpI6enpPl8TQHC01DxW5/e//72WLl2qf/zjH+rUqVOzrweEM5ZUeOGSSy7Rb3/7W40bN67B16uqqnTXXXfpBz/4gVq3bq1BgwYpPz/fb/1J0r333qtLL71UDz/8sPr166ezzjpLl112mZKSknzuN9xFRkZqzZo1euedd9S7d2898cQT6tmzp3bu3HnGc+um7V5++eU6fPiwlixZosLCQtd6wLq1hXFxcU1ep7a2VgsWLFBRUZHr8cknn2j79u2NFoOaO41Pkk6cOKHNmzcrMjJS27dvP+P7BRB6WmoOi4mJUXZ2tpYtW6aqqir94Q9/0I033njG9wwg9LTUPCZJ+fn5uuWWW/TCCy/o3HPPPeP7BeyOGQ4WmjJlinbt2qVVq1YpJSVFubm5Gj16tD755BOdffbZlvdXW1urt956S3fddZdGjRqlLVu2qHv37pozZ46uvPJKy/sLJw6HQxdeeKEuvPBCzZ07V2lpacrNzdXs2bMlSR9//LG+/fZb119WBQUFatOmjbp27arDhw+ruLhYzz77rGvd3Ycffuh2/XPOOUfPPfecvv766wYr6+eff75KSkr0H//xHx7HPHbsWA0aNKjJNsnJyU2+fvvttysiIkLvvPOOLr30Ul122WX6yU9+4nEMAEJDS81hU6dOVZ8+ffT000+rurq6yYI7gNDWEvPY559/rquvvlr33HMP+QuoE+w1HeFKp+2p8PnnnxuHw2G+/PJLt3YXX3yxmTNnjuX9GWPMgQMHjCQTHx9vHn30UbNlyxaTk5NjHA6Hyc/Pb3af4aqgoMA88MADZuPGjWb37t3m1VdfNTExMebtt982xny/brBNmzbmuuuuM5999pl5++23TXJysrn77ruNMcbU1NSYDh06mJ/97Gdm+/bt5v333zcXXHCB28+gsrLS9OjRwwwbNsx8+OGHZseOHea1114z69atM8YYk5eXZ6Kiosy8efPMp59+arZu3WpWrVpl7r33Xr+97zfffNPExMSYzZs3G2OM+fWvf226du1qvv76a7/1CcB6LTWH1RkyZIiJiYkxN910k9/7AuAfLTGPnThxwmRkZJif/OQnZv/+/ebAgQOuB9CSUXDw0ekFgFdffdVIMq1bt3Z7REVFmWuvvdYYY8zOnTuNpCYfv/zlLz3qzxhjvvzySyPJXHfddW7HL7/8cjNhwgRL32842bp1qxk1apTp1KmTiY2NNT169DBPPPGE6/VJkyaZK664wsydO9d06NDBtGnTxkydOtW1UZkxxqxZs8b06tXLxMbGmnPOOcfk5+fX+xns2rXLXH311aZt27YmPj7eDBgwwBQWFrpez8vLM0OGDDFxcXGmbdu2ZuDAgeb3v/+9X95zWVmZSU5Odtscqbq62gwcOND1/QMQHlpiDjvV0qVLjSSzYcMGv/cFwD9aYh5r6t/5QEvmMIabW/vC4XAoNzfXtXThlVde0fXXX6/PPvtMkZGRbm3btGmjzp07q7q6ut79g0/Xvn37Bqdqnd6f9P0attatW2vevHn69a9/7Tr+q1/9Sh9++KH+8Y9/+P4GbWzy5Mk6cuSI/vznPwc7FADwmt1z2AMPPKBVq1bpk08+CXYoAPzE7nkMwL+xh4NF+vXrp5qaGpWVlTV6v93o6GhlZGRY1mdMTIwuuOAClZSUuB3ftm2b0tLSLOsHAAB/O3bsmIqLi/XEE0/o/vvvD3Y4AADAAhQcvHDs2DF9/vnnruc7d+5UUVGREhMT1aNHD11//fW64YYb9Lvf/U79+vXToUOH9Ne//lV9+/bVpZdeaml/3bp1kyTdeeedGj9+vH784x/roosuUl5ent54441m3R0DAIBAmzFjhl5++WVdeeWV3J0CAACbYEmFF/Lz83XRRRfVOz5p0iQtX75c1dXV+u1vf6sXXnhBX375pTp06KDMzEwtWLBAffv2tby/Os8//7xycnK0b98+9ezZUwsWLNAVV1zhdX8AAAAAAFiFggMAAAAAALBcRLADAAAAAAAA9kPBAQAAAAAAWI5NI8+gtrZW+/fvV0JCghwOR7DDAWBzxhgdPXpUKSkpiohofk2YHAYg0MhjAMKZ1TmspaPgcAb79+9XampqsMMA0MLs3btXXbt2bfZ1yGEAgoU8BiCcWZXDWjoKDmeQkJAg6fsvXNu2bYMcDQC7q6ioUGpqqiv3NBc5DECgkccAhDOrc1hLR8HhDOqm7rVt25a/5AAEjFXThslhAIKFPAYgnLGEyxosSgEAAAAAAJaj4AAAAAAAACwXdgWHp59+Wt27d1erVq3Uv39//f3vf2+y/dq1a9W/f3+1atVKZ511lp555pkARQoAAAAAQMsVVgWHV155RbfddpvuvfdebdmyRcOGDdMll1yiPXv2NNh+586duvTSSzVs2DBt2bJF99xzj2bOnKnXX389wJEDAAAAANCyhFXB4dFHH9XPf/5zTZ06Vb169dKiRYuUmpqqxYsXN9j+mWeeUbdu3bRo0SL16tVLU6dO1Y033qhHHnkkwJEDAAAAANCyhM1dKqqqqrR582bdfffdbsezsrK0bt26Bs9Zv369srKy3I6NGjVKS5cuVXV1taKjoy2P09ScVOW29y2/rsf9t+8uk3iW5ydUn1DEvkKpttZ/QQEtWG3K+VJce8VFR4bFbsemtlbfnjga7DAAhJLoeMXFRIVFDpPIYwBOE2Y5zG7CpuBw6NAh1dTUKDk52e14cnKySktLGzyntLS0wfYnT57UoUOH1KVLl3rnVFZWqrKy0vW8oqLCqzi//e6E4l+51qtzrHTSRCiz8kkdVDuP2j8VvUiXRW7wb1BACzah6tcqqO2trfeNUnyM/1Nus3PYiaOKf6Sb1WEBCGO9vntem++7IiA5TCKPAbBWoHMY3IXdp356ZcoY02S1qqH2DR2vk5OTowULFjQnQn1Sm96M832X4diraEeNkhzf6KBp59E5XR2HJEm7apN1VHF+jA5omY6bVgHtr/k5DACCizwGAPbhMHUj8BBXVVWl+Ph4/fGPf9RVV13lOn7rrbeqqKhIa9eurXfOj3/8Y/Xr10+PPfaY61hubq6uvfZanThxosElFQ1V1VNTU1VeXq62bdueMU5jjL6trvH27Vmi1ZPnKKLiS3035X9V26WfZ+c8/xNFlH6s78a/qtofXuznCIGWy9MlFRUVFXI6nR7nnNM1O4cxFRnA6bycjkweAxBSApzD4C5sZjjExMSof//+WrNmjVvBYc2aNbriiisaPCczM1NvvPGG27H33ntPAwYMaHT/htjYWMXGxvocp8PhCN50Hcf3e4C2ioyQPI7h+70bWkVHenEOgFDV7BwWEaH4Nk4LIwIA75DHAMA+wuouFbNnz9Zzzz2n559/XsXFxZo1a5b27Nmjm266SZI0Z84c3XDDDa72N910k3bv3q3Zs2eruLhYzz//vJYuXao77rgjWG/Bv1xVOy8mrdQ1dYTVVwEAAAAAEOLC6lfa48eP1+HDh3XffffpwIED6tOnj95++22lpaVJkg4cOKA9e/a42nfv3l1vv/22Zs2apaeeekopKSl6/PHHdfXVVwfrLfhXXdHAeHHHibq2FBwAAAAAABYKq4KDJN188826+eabG3xt+fLl9Y4NHz5cH330kZ+jChEUHAAAAAAAIYJRpp1QcAAAAAAAhAhGmbbyrz0cvLrxiGsTB6uDAQAAAAC0YBQc7IQZDgAAAACAEMEo004oOAAAAAAAQgSjTDuh4AAAAAAACBGMMu3EUbeHgy8FB/ZwAAAAAABYh4KDnbiKBl5sGunaM5KCAwAAAADAOhQc7MS1pMKbggNLKgAAAAAA1mOUaSfs4QAAAAAACBGMMu2kOQUHsaQCAAAAAGAdCg62UrdppBdLKuo2cWCGAwAAAADAQowy7YQlFQAAAACAEMEo004oOAAAAAAAQgSjTDuh4AAAAAAACBGMMu3EUbeHgy8FBzaNBAAAAABYh4KDnbhmKXixaWRdU2Y4AAAAAAAsxCjTThw+3KWCGQ4AAAAAAD+g4GAn7OEAAAAAAAgRjDLtpDkFBzHDAQAAAABgHQoOtuLDkoq6TRyY4QAAAAAAsBCjTDthSQUAAAAAIEQwyrQTCg4AAAAAgBDBKNNOKDgAAAAAAEJE2Iwyv/nmG2VnZ8vpdMrpdCo7O1tHjhxptH11dbV+9atfqW/fvmrdurVSUlJ0ww03aP/+/YELOtCaVXBg00gAAAAAgHXCpuAwceJEFRUVKS8vT3l5eSoqKlJ2dnaj7U+cOKGPPvpIv/nNb/TRRx/pT3/6k7Zt26axY8cGMOoAcxUNPNw08tTNJZnhAAAAAACwUFSwA/BEcXGx8vLyVFBQoEGDBkmSlixZoszMTJWUlKhnz571znE6nVqzZo3bsSeeeEIDBw7Unj171K1bt4DEHlB1BQdPZzhQcAAAAAAA+ElYjDLXr18vp9PpKjZI0uDBg+V0OrVu3TqPr1NeXi6Hw6F27dr5IcoQ4FpS4ekMh1MKEyypAAAAAABYKCxmOJSWliopKane8aSkJJWWlnp0je+++0533323Jk6cqLZt2zbarrKyUpWVla7nFRUV3gccLN7u4eDWjoIDYAdhncMAQOQxALCToM5wmD9/vhwOR5OPTZs2SZIcDfwG3hjT4PHTVVdXa8KECaqtrdXTTz/dZNucnBzXxpROp1Opqam+vbmg8HZJxakzHMJisguAMwjvHAYA5DEAsJOgjjJnzJih4uLiJh99+vRR586d9dVXX9U7/+DBg0pOTm6yj+rqal177bXauXOn1qxZ0+TsBkmaM2eOysvLXY+9e/c26z0GlLdLKsQeDoDdhHUOAwCRxwDAToK6pKJjx47q2LHjGdtlZmaqvLxcGzZs0MCBAyVJhYWFKi8v15AhQxo9r67YsH37dn3wwQfq0KHDGfuKjY1VbGys528ilDRnSQUFB8AWwjqHAYDIYwBgJ2ExyuzVq5dGjx6tadOmqaCgQAUFBZo2bZrGjBnjdoeKjIwM5ebmSpJOnjypa665Rps2bdLKlStVU1Oj0tJSlZaWqqqqKlhvxb8oOAAAAAAAQkTYjDJXrlypvn37KisrS1lZWTrnnHP04osvurUpKSlReXm5JGnfvn1avXq19u3bp/POO09dunRxPby5s0VYaVbBgU0jAQAAAADWCYu7VEhSYmKiXnrppSbbmFP2LkhPT3d73iK4igae3haTPRwAAAAAAP7BKNNOHNylAgAAAAAQGhhl2onXSyqY4QAAAAAA8A9GmXbi7W0x2cMBAAAAAOAnFBxsxdclFRQbAAAAAADWouBgJ97OcKjbXJLlFAAAAAAAizHStBNfb4tJwQEAAAAAYDFGmnZCwQEAAAAAECIYadqJzwUH9nAAAAAAAFiLgoOduAoHnt6lgj0cAAAAAAD+wUjTTlhSAQAAAAAIEYw07cTh420xKTgAAAAAACzGSNNOvL0tpqsdezgAAAAAAKxFwcFWfJ3hQMEBAAAAAGAtCg524u0MB7FpJAAAAADAPxhp2gmbRgIAAAAAQgQjTTuh4AAAAAAACBGMNO3E54IDezgAAAAAAKxFwcFOXIUDL+9SwQwHAAAAAIDFGGnaCUsqAAAAAAAhgpGmnTh8vS0mXwMAAAAAgLUYadqJ1zMc6pZesIcDAAAAAMBaFBzsxFVw8HQPBzaNBAAAAAD4BwUHW6lbUuFhwUFsGgkAAAAA8I+wGWl+8803ys7OltPplNPpVHZ2to4cOeLx+b/4xS/kcDi0aNEiv8UYdGwaCQAAAAAIEWEz0pw4caKKioqUl5envLw8FRUVKTs726Nz//znP6uwsFApKSl+jjLIKDgAAAAAAEJEVLAD8ERxcbHy8vJUUFCgQYMGSZKWLFmizMxMlZSUqGfPno2e++WXX2rGjBl69913ddlllwUq5ODwueDAHg4AAAAAAGuFxa+2169fL6fT6So2SNLgwYPldDq1bt26Rs+rra1Vdna27rzzTv3oRz8KRKjBxW0xAQAAAAAhIixmOJSWliopKane8aSkJJWWljZ63kMPPaSoqCjNnDnT474qKytVWVnpel5RUeFdsMHkKhx4epcKNo0E7CascxgAiDwGAHYS1JHm/Pnz5XA4mnxs2rRJkuRoYNq/MabB45K0efNmPfbYY1q+fHmjbRqSk5Pj2pjS6XQqNTXVtzcXDOzhALR4YZ3DAEDkMQCwk6CONGfMmKHi4uImH3369FHnzp311Vdf1Tv/4MGDSk5ObvDaf//731VWVqZu3bopKipKUVFR2r17t26//Xalp6c3GtOcOXNUXl7ueuzdu9eqt+t/vi6pEHs4AHYR1jkMAEQeAwA7CeqSio4dO6pjx45nbJeZmany8nJt2LBBAwcOlCQVFhaqvLxcQ4YMafCc7OxsjRw50u3YqFGjlJ2drSlTpjTaV2xsrGJjY714FyHENcPB2yUVFBwAuwjrHAYAIo8BgJ2ExR4OvXr10ujRozVt2jQ9++yzkqTp06drzJgxbneoyMjIUE5Ojq666ip16NBBHTp0cLtOdHS0Onfu3ORdLcJb3QwHDwsOYg8HAAAAAIB/hEXBQZJWrlypmTNnKisrS5I0duxYPfnkk25tSkpKVF5eHozwQgN7OAAAAABAyKipqVF1dXWww/Cb6OhoRUZGNvp62BQcEhMT9dJLLzXZxpzhN/u7du2yMKIQRMEBAAAAAELCsWPHtG/fvjOOU8OZw+FQ165d1aZNmwZfD5uCAzzgc8GBPRwAAAAAwCo1NTXat2+f4uPj1alTJ6/unBgujDE6ePCg9u3bp7PPPrvBmQ4UHOzE17tUMMMBAAAAACxTXV0tY4w6deqkuLi4YIfjN506ddKuXbtUXV3dYMGBkaaduAoH3t6lgq8BAAAAAFjNjjMbTnWm98dI007YwwEAAAAAECIYadqJrwUH2bvqBgAAAAB2M2LECN12222Nvp6enq5FixYFLJ6GsIeDnbj2cPB0SQWbRgIAAABAOPrTn/6k6OjoYIfRJAoOtuJlwcF1GhNdAAAAACCcJCYmBjuEM2KkaSfs4QAAAAAALcKpSyrKysp0+eWXKy4uTt27d9fKlSuDG9y/MMPBTnwuOLCkAgAAAADC1eTJk7V371799a9/VUxMjGbOnKmysrJgh0XBwVaY4QAAAAAALcq2bdv0zjvvqKCgQIMGDZIkLV26VL169QpyZCypsBfXppEUHAAAAACgJSguLlZUVJQGDBjgOpaRkaF27doFL6h/YaRpJ67Cgad3qTCnnQcAAAAACCfmX+M6RwgulWekaScsqQAAAACAFqVXr146efKkNm3a5DpWUlKiI0eOBC+of2GkaSe+FhwUepUwAAAAAMCZ9ezZU6NHj9a0adNUWFiozZs3a+rUqYqLiwt2aBQcbMXnPRwoOAAAAABAuFq2bJlSU1M1fPhwjRs3TtOnT1dSUlKww/L+LhXHjx/Xgw8+qPfff19lZWWqrXUf3H7xxReWBQdv1RUcPG3PHg4AAAAAEI7y8/Ndf+7cubPefPNNt9ezs7MDHFF9Xhccpk6dqrVr1yo7O1tdunQJyY0pWiyvl1RQcAAAAAAA+IfXBYd33nlHb731li688EJ/xIPm8HnTSIpGAAAAAABref2r7fbt2ysxMdEfsaC5uEsFAAAAACBEeD3SvP/++zV37lydOHHCH/GgOSg4AAAAAABChEdLKvr16+e2V8Pnn3+u5ORkpaenKzo62q3tRx99ZG2E8JzrZ+ThrpHs4QAAAAAA8BOPCg5XXnmln8OAJZjhAAAAAAAIER4VHObNm+f685QpU3T99dfr4osv5g4VocbXgoP4OQIAAAAArOX1r7YPHz6sMWPGqGvXrrrjjjtUVFTkh7Dq++abb5SdnS2n0ymn06ns7GwdOXLkjOcVFxdr7NixcjqdSkhI0ODBg7Vnzx7/BxwMzHAAAAAAAIQIr0eaq1evVmlpqebNm6dNmzapf//+6t27txYuXKhdu3b5IcTvTZw4UUVFRcrLy1NeXp6KioqUnZ3d5Dk7duzQ0KFDlZGRofz8fH388cf6zW9+o1atWvktzuD610wF4+keDtwWEwAAAADgHx4tqThdu3btNH36dE2fPl379u3Tyy+/rOeff15z587VyZMnrY5RxcXFysvLU0FBgQYNGiRJWrJkiTIzM1VSUqKePXs2eN69996rSy+9VA8//LDr2FlnnWV5fCHD4WXBQWwaCQAAAADwj2aNNKurq7Vp0yYVFhZq165dSk5OtiouN+vXr5fT6XQVGyRp8ODBcjqdWrduXYPn1NbW6q233lKPHj00atQoJSUladCgQfrzn//slxhDgtdLKuoKDsxwAAAAAABIf/vb33T55ZcrJSVFDoejWWNonwoOH3zwgaZNm6bk5GRNmjRJCQkJeuONN7R3716fA2lKaWmpkpKS6h1PSkpSaWlpg+eUlZXp2LFjevDBBzV69Gi99957uuqqqzRu3DitXbu20b4qKytVUVHh9ggb7OEAtHhhncMAQOQxAAi248eP69xzz9WTTz7Z7Gt5vaSia9euOnz4sEaNGqVnn31Wl19+uc97IsyfP18LFixoss3GjRslqcE7YhhjGr1TRm3t94PpK664QrNmzZIknXfeeVq3bp2eeeYZDR8+vMHzcnJyzhhTyKLgALR4YZ3DAEDkMQAItksuuUSXXHKJJdfyuuAwd+5c/fSnP1X79u2b3fmMGTM0YcKEJtukp6frn//8p7766qt6rx08eLDRZRwdO3ZUVFSUevfu7Xa8V69e+vDDDxvtb86cOZo9e7breUVFhVJTU5uMMWS4ii+ebhrJHg6A3YR1DgMAkccA2JMxRt9W1wSl77joyEZ/Ue9vXhccpk+fblnnHTt2VMeOHc/YLjMzU+Xl5dqwYYMGDhwoSSosLFR5ebmGDBnS4DkxMTG64IILVFJS4nZ827ZtSktLa7Sv2NhYxcbGevEuQggzHIAWL6xzGACIPAbAnr6trlHvue8Gpe+t941SfIxP94totrAYafbq1UujR4/WtGnTVFBQoIKCAk2bNk1jxoxxu0NFRkaGcnNzXc/vvPNOvfLKK1qyZIk+//xzPfnkk3rjjTd08803B+Nt+J+vBQexaSQAAAAAwFrBKXP4YOXKlZo5c6aysrIkSWPHjq23iUVJSYnKy8tdz6+66io988wzysnJ0cyZM9WzZ0+9/vrrGjp0aEBjDxhmOAAAAABAyImLjtTW+0YFre9gCZuCQ2Jiol566aUm2xhTf++CG2+8UTfeeKO/wgox/5qp0MDn0CBXwYEZDgAAAADgLw6HI2jLGoKp5b1jO3N4WXAQm0YCAAAAAP7t2LFj+vzzz13Pd+7cqaKiIiUmJqpbt25eXYuCg514vaSiruDADAcAAAAAgLRp0yZddNFFrud1dw6aNGmSli9f7tW1KDjYCXs4AAAAAACaYcSIEQ1uV+ALRpp2QsEBAAAAABAiGGnaiWtphKebRrKHAwAAAADAPxhp2gkzHAAAAAAAIYKRpp34WnAQm0YCAAAAAKxFwcFOmOEAAAAAAAgRjDRt5ZSZCp7sKkrBAQAAAADgJ4w07eTUwoFHtzGp2zSSJRUAAAAAAGtRcLCTUwsHniyrcM1woOAAAAAAALAWBQc7cZvh4EnBgdtiAgAAAAD8g5GmnXhdcGAPBwAAAACAfzDStBO3pRGebBrJDAcAAAAAwPdycnJ0wQUXKCEhQUlJSbryyitVUlLi8/UYadoJMxwAAAAAAD5au3atfvnLX6qgoEBr1qzRyZMnlZWVpePHj/t0vSiL40Mw+VpwEJtGAgAAAEBLl5eX5/Z82bJlSkpK0ubNm/XjH//Y6+tRcLATZjgAAAAAQOgxRqo+EZy+o+N9vjNheXm5JCkxMdGn8yk42AkFBwAAAAAIPdUnpIUpwen7nv1STGuvTzPGaPbs2Ro6dKj69OnjU9cUHGzllKqV8WDTyLqNJX2sdgEAAAAA7GnGjBn65z//qQ8//NDna1BwsBO3GQ6e3KWiboYDBQcAAAAA8Jvo+O9nGgSrby/dcsstWr16tf72t7+pa9euPndNwcFOTi0ceLSkgttiAgAAAIDfORw+LWsINGOMbrnlFuXm5io/P1/du3dv1vUoONiJw6Hvl1UY9nAAAAAAAHjll7/8pf7whz/oL3/5ixISElRaWipJcjqdiouL8/p6jDTtpm6WAwUHAAAAAIAXFi9erPLyco0YMUJdunRxPV555RWfrhc2I81vvvlG2dnZcjqdcjqdys7O1pEjR5o859ixY5oxY4a6du2quLg49erVS4sXLw5MwMHiKh54socDSyoAAAAAAN8zxjT4mDx5sk/XC5uR5sSJE1VUVKS8vDzl5eWpqKhI2dnZTZ4za9Ys5eXl6aWXXlJxcbFmzZqlW265RX/5y18CFHUQ1BUPvJnhIDaNBAAAAABYKywKDsXFxcrLy9Nzzz2nzMxMZWZmasmSJXrzzTdVUlLS6Hnr16/XpEmTNGLECKWnp2v69Ok699xztWnTpgBGH2C+FByY4QAAAAAAsFhYjDTXr18vp9OpQYMGuY4NHjxYTqdT69ata/S8oUOHavXq1fryyy9ljNEHH3ygbdu2adSoUYEIOzgoOAAAAAAAQkBY3KWitLRUSUlJ9Y4nJSW5ds1syOOPP65p06apa9euioqKUkREhJ577jkNHTq00XMqKytVWVnpel5RUdG84AOubtNID/ZwEHs4AHYT/jkMQEtHHgMA+wjqSHP+/PlyOBxNPuqWPzgc9fcZMMY0eLzO448/roKCAq1evVqbN2/W7373O91888363//930bPycnJcW1M6XQ6lZqa2vw3Gkg+zXBgDwfALsI+hwFo8chjAGAfQZ3hMGPGDE2YMKHJNunp6frnP/+pr776qt5rBw8eVHJycoPnffvtt7rnnnuUm5uryy67TJJ0zjnnqKioSI888ohGjhzZ4Hlz5szR7NmzXc8rKirC6y86V8HBm7tUUHAA7CLscxiAFo88BsBOjEczz8PXmd5fUAsOHTt2VMeOHc/YLjMzU+Xl5dqwYYMGDhwoSSosLFR5ebmGDBnS4DnV1dWqrq5WRIT7JI7IyEjV1jb+2//Y2FjFxsZ68S5CTF3xgD0cgBYp7HMYgBaPPAbADiIjIyVJVVVViouLC3I0/lNVVSXp3+/3dGGxh0OvXr00evRoTZs2Tc8++6wkafr06RozZox69uzpapeRkaGcnBxdddVVatu2rYYPH64777xTcXFxSktL09q1a/XCCy/o0UcfDdZb8T8KDgAAAAAQVFFRUYqPj9fBgwcVHR1d7xfhdlBbW6uDBw8qPj5eUVENlxbCouAgSStXrtTMmTOVlZUlSRo7dqyefPJJtzYlJSUqLy93PV+1apXmzJmj66+/Xl9//bXS0tL0wAMP6Kabbgpo7AHlKh54saRCLKkAAAAAAKs4HA516dJFO3fu1O7du4Mdjt9ERESoW7duje6tGDYFh8TERL300ktNtjl9/Ujnzp21bNkyf4YVergtJgAAAAAEXUxMjM4++2zXsgM7iomJaXL2RtgUHOAhCg4AAAAAEBIiIiLUqlWrYIcRNIw07YaCAwAAAAAgBDDStJ26TSM9uf1K3W0x+RoAAAAAAKzFSNNufJrhwKaRAAAAAABrUXCwG1fBwZO7VFBwAAAAAAD4BwUHu6krHng0w4ElFQAAAAAA/2CkaTdeFRzYNBIAAAAA4B+MNO3GVTzwZElFXRuWVAAAAAAArEXBwW64LSYAAAAAIAQw0rQbCg4AAAAAgBDASNNuKDgAAAAAAEIAI03bqds00oM9HMRdKgAAAAAA/sFI0258muHAppEAAAAAAGtRcLAbllQAAAAAAEIAI027cXixpKKuDTMcAAAAAAAWo+BgN8xwAAAAAACEAEaaduOareDFDAcxwwEAAAAAYC0KDnbDDAcAAAAAQAhgpGk3FBwAAAAAACGAkabdUHAAAAAAAIQARpq2U3eXCgoOAAAAAIDgYaRpN64ZDh5sGiluiwkAAAAA8A8KDnbDkgoAAAAAQAgIm5HmAw88oCFDhig+Pl7t2rXz6BxjjObPn6+UlBTFxcVpxIgR+uyzz/wbaLA5vFlSUTfDIWy+BgAAAACAMBE2I82qqir99Kc/1X/91395fM7DDz+sRx99VE8++aQ2btyozp076z//8z919OhRP0YaZN4sqXDNcGBJBQAAAADAWmFTcFiwYIFmzZqlvn37etTeGKNFixbp3nvv1bhx49SnTx+tWLFCJ06c0B/+8Ac/RxtEruKBJwWHujYUHAAAAAAA1gqbgoO3du7cqdLSUmVlZbmOxcbGavjw4Vq3bl0QI/Mz9nAAAAAAAISAqGAH4C+lpaWSpOTkZLfjycnJ2r17d6PnVVZWqrKy0vW8oqLCPwH6CwUHoEUL+xwGoMUjjwGAfQR1pDl//nw5HI4mH5s2bWpWH47T9icwxtQ7dqqcnBw5nU7XIzU1tVn9BxwFB6BFC/scBqDFI48BgH0EdaQ5Y8YMFRcXN/no06ePT9fu3LmzpH/PdKhTVlZWb9bDqebMmaPy8nLXY+/evT71Hzze3KWCggNgN+GfwwC0dOQxALCPoC6p6Nixozp27OiXa3fv3l2dO3fWmjVr1K9fP0nf3+li7dq1euihhxo9LzY2VrGxsX6JKSC8uUtF3caS3KUCsI2wz2EAWjzyGADYR9j8anvPnj0qKirSnj17VFNTo6KiIhUVFenYsWOuNhkZGcrNzZX0/VKK2267TQsXLlRubq4+/fRTTZ48WfHx8Zo4cWKw3ob/saQCAAAAABACwmbTyLlz52rFihWu53WzFj744AONGDFCklRSUqLy8nJXm7vuukvffvutbr75Zn3zzTcaNGiQ3nvvPSUkJAQ09oByeLik4tQZEBQcAAAAAAAWC5uCw/Lly7V8+fIm25jTlhE4HA7Nnz9f8+fP919gocbTGQ6nvk7BAQAAAABgMUaaduPpfgwe7fEAAAAAAIBvKDjYDTMcAAAAAAAhgJGm3VBwAAAAAACEAEaadkPBAQAAAAAQAhhp2o6nd6mg4AAAAAAA8B9GmnbjmuFwpk0hT70tpocbTQIAAAAA4CEKDnbDkgoAAAAAQAhgpGk3FBwAAAAAACGAkabdODzdw+HUJRV8DQAAAAAA1mKkaTeu/RjOsIeD2x4P7OEAAAAAALAWBQe78XTTSLclFRQcAAAAAADWigp2ALBYXcHhoxelnX9rvF1NVd0JFBwAAAAAAJaj4GA3bVO+/2/5nu8fnrYHAAAAAMBCFBzsZshMqfM5UtVxz9p3HeDfeAAAAAAALRIFB7uJipV6jAp2FAAAAACAFo5NIwEAAAAAgOUoOAAAAAAAAMtRcAAAAAAAAJaj4AAAAAAAACxHwQEAAAAAAFiOggMAAAAAALAct8U8A2OMJKmioiLIkQBoCepyTV3uaS5yGIBAI48BCGdW57CWjoLDGRw9elSSlJqaGuRIALQkR48eldPptOQ6EjkMQOCRxwCEM6tyWEvnMJRumlRbW6v9+/crISFBDofDo3MqKiqUmpqqvXv3qm3btn6O0HrEH1zhHH84xy6FRvzGGB09elQpKSmKiGj+qjdyGPEHGvEHVyjETx5rvnCOP5xjl4g/2EIhfqtzWEvHDIcziIiIUNeuXX06t23btmH5P3od4g+ucI4/nGOXgh+/ldV0chjxBwvxB1ew4yePWSOc4w/n2CXiD7Zgx8/MButQsgEAAAAAAJaj4AAAAAAAACxHwcEPYmNjNW/ePMXGxgY7FJ8Qf3CFc/zhHLsU/vFbJdw/B+IPLuIPrnCP3yrh/jmEc/zhHLtE/MEW7vGjPjaNBAAAAAAAlmOGAwAAAAAAsBwFBwAAAAAAYDkKDgAAAAAAwHIUHBrw9NNPq3v37mrVqpX69++vv//97022X7t2rfr3769WrVrprLPO0jPPPFOvzeuvv67evXsrNjZWvXv3Vm5ubrP7DVT8S5Ys0bBhw9S+fXu1b99eI0eO1IYNG9zazJ8/Xw6Hw+3RuXPnkIh/+fLl9WJzOBz67rvvmtVvoOIfMWJEg/FfdtllrjZWff7exH7gwAFNnDhRPXv2VEREhG677bYG24Xqd9+T+AP93bcKOYwcFkqffyBzmLfxk8e+Rx7zvd9AxR/InyM5jBxGDoOlDNysWrXKREdHmyVLlpitW7eaW2+91bRu3drs3r27wfZffPGFiY+PN7feeqvZunWrWbJkiYmOjjavvfaaq826detMZGSkWbhwoSkuLjYLFy40UVFRpqCgwOd+Axn/xIkTzVNPPWW2bNliiouLzZQpU4zT6TT79u1ztZk3b5750Y9+ZA4cOOB6lJWVeRW7v+JftmyZadu2rVtsBw4caFa/gYz/8OHDbnF/+umnJjIy0ixbtszVxorP39vYd+7caWbOnGlWrFhhzjvvPHPrrbfWaxPK331P4g/kd98q5DByWKh9/oHKYb7ETx77HnnMt34DGX+gfo7kMHIYOQxWo+BwmoEDB5qbbrrJ7VhGRoa5++67G2x/1113mYyMDLdjv/jFL8zgwYNdz6+99lozevRotzajRo0yEyZM8LnfQMZ/upMnT5qEhASzYsUK17F58+aZc88916tYG+KP+JctW2acTqel/Vp1HV8+///5n/8xCQkJ5tixY65jVnz+zfkMhg8f3uBfEqH83T9VY/Gfzp/ffauQw8hhnvRr1XVCKYcZQx4jj/0becw35DBymLf9noochoawpOIUVVVV2rx5s7KystyOZ2Vlad26dQ2es379+nrtR40apU2bNqm6urrJNnXX9KXfQMZ/uhMnTqi6ulqJiYlux7dv366UlBR1795dEyZM0BdffOFx7P6O/9ixY0pLS1PXrl01ZswYbdmypVn9Bjr+Uy1dulQTJkxQ69at3Y435/O36jM4XSh/933hr+++VchhDcd/OnJY4OM/lT9ymK/xeyKUv/++II+5twmFn2Oo5DFyWOPxn4oc1jByGBpDweEUhw4dUk1NjZKTk92OJycnq7S0tMFzSktLG2x/8uRJHTp0qMk2ddf0pd9Axn+6u+++Wz/4wQ80cuRI17FBgwbphRde0LvvvqslS5aotLRUQ4YM0eHDh4Mef0ZGhpYvX67Vq1fr5ZdfVqtWrXThhRdq+/btPvcbyPhPtWHDBn366aeaOnWq2/Hmfv5WfQanC+Xvvi/89d23Cjms4fhPRw4LbPyn8lcO8zV+T4Ty998X5LGGrxlO/x9J/vk5ksMajv9U5LDGkcPQmKhgBxCKHA6H23NjTL1jZ2p/+nFPrultv97E09z46zz88MN6+eWXlZ+fr1atWrmOX3LJJa4/9+3bV5mZmfrhD3+oFStWaPbs2UGNf/DgwRo8eLDr9QsvvFDnn3++nnjiCT3++OM+9xuo+E+1dOlS9enTRwMHDnQ7btXnb9Vn4O01g/XZeyMQ332rkMPIYZ70G6j4T+XvHOZL/FZdkzxmLfJY8PIYOYwc1px+yWE4HTMcTtGxY0dFRkbWq8KVlZXVq9bV6dy5c4Pto6Ki1KFDhybb1F3Tl34DGX+dRx55RAsXLtR7772nc845p8lYWrdurb59+7qq16EQf52IiAhdcMEFrtjC5fM/ceKEVq1aVa+q3hBvP3+rPoPThfJ33xv+/u5bhRzWcPx1yGHBjd+fOczX+D0Ryt9/b5DHQvfnGCp5jBzWdPzksKaRw9AYCg6niImJUf/+/bVmzRq342vWrNGQIUMaPCczM7Ne+/fee08DBgxQdHR0k23qrulLv4GMX5L++7//W/fff7/y8vI0YMCAM8ZSWVmp4uJidenSJSTiP5UxRkVFRa7YwuHzl6RXX31VlZWV+tnPfnbGWLz9/K36DE4Xyt99TwXiu28VcljD8UvksFCI3585zNf4PRHK339PkcdC++cYKnmMHNZ0/OSwppHD0Cjr96EMb3W3c1m6dKnZunWrue2220zr1q3Nrl27jDHG3H333SY7O9vVvu52OrNmzTJbt241S5curXc7nX/84x8mMjLSPPjgg6a4uNg8+OCDjd6OprF+gxn/Qw89ZGJiYsxrr73mdruZo0ePutrcfvvtJj8/33zxxRemoKDAjBkzxiQkJIRE/PPnzzd5eXlmx44dZsuWLWbKlCkmKirKFBYWetxvMOOvM3ToUDN+/PgG+7Xi8/c2dmOM2bJli9myZYvp37+/mThxotmyZYv57LPPXK+H8nffk/gD+d23CjmMHBZqn38df+cwX+I3hjxmDHnM036DGX+gfo7kMHIYOQxWo+DQgKeeesqkpaWZmJgYc/7555u1a9e6Xps0aZIZPny4W/v8/HzTr18/ExMTY9LT083ixYvrXfOPf/yj6dmzp4mOjjYZGRnm9ddf96rfYMaflpZmJNV7zJs3z9Vm/PjxpkuXLiY6OtqkpKSYcePGuSWLYMZ/2223mW7dupmYmBjTqVMnk5WVZdatW+dVv8GM3xhjSkpKjCTz3nvvNdinVZ+/t7E39L1IS0tzaxPK3/0zxR/o775VyGHksFD6/I0JXA7zJX7yGHnM036DGX8gf47kMHIYOQxWchjzr51RAAAAAAAALMIeDgAAAAAAwHIUHAAAAAAAgOUoOAAAAAAAAMtRcAAAAAAAAJaj4AAAAAAAACxHwQEAAAAAAFiOggMAAAAAALAcBQcAAAAAAGA5Cg6Al/Lz8+VwOHTkyJFghwIAPiGPAQhn5DAgfDiMMSbYQQChbMSIETrvvPO0aNEiSVJVVZW+/vprJScny+FwBDc4APAAeQxAOCOHAeErKtgBAOEmJiZGnTt3DnYYAOAz8hiAcEYOA8IHSyqAJkyePFlr167VY489JofDIYfDoeXLl7tN41u+fLnatWunN998Uz179lR8fLyuueYaHT9+XCtWrFB6errat2+vW265RTU1Na5rV1VV6a677tIPfvADtW7dWoMGDVJ+fn5w3igA2yKPAQhn5DAgvDHDAWjCY489pm3btqlPnz667777JEmfffZZvXYnTpzQ448/rlWrVuno0aMaN26cxo0bp3bt2untt9/WF198oauvvlpDhw7V+PHjJUlTpkzRrl27tGrVKqWkpCg3N1ejR4/WJ598orPPPjug7xOAfZHHAIQzchgQ3ig4AE1wOp2KiYlRfHy8a+re//3f/9VrV11drcWLF+uHP/yhJOmaa67Riy++qK+++kpt2rRR7969ddFFF+mDDz7Q+PHjtWPHDr388svat2+fUlJSJEl33HGH8vLytGzZMi1cuDBwbxKArZHHAIQzchgQ3ig4ABaIj493/QUnScnJyUpPT1ebNm3cjpWVlUmSPvroIxlj1KNHD7frVFZWqkOHDoEJGgBOQR4DEM7IYUBoouAAWCA6OtrtucPhaPBYbW2tJKm2tlaRkZHavHmzIiMj3dqd+hcjAAQKeQxAOCOHAaGJggNwBjExMW4bDFmhX79+qqmpUVlZmYYNG2bptQHgdOQxAOGMHAaEL+5SAZxBenq6CgsLtWvXLh06dMhVGW+OHj166Prrr9cNN9ygP/3pT9q5c6c2btyohx56SG+//bYFUQPAv5HHAIQzchgQvig4AGdwxx13KDIyUr1791anTp20Z88eS667bNky3XDDDbr99tvVs2dPjR07VoWFhUpNTbXk+gBQhzwGIJyRw4Dw5TDGmGAHAQAAAAAA7IUZDgAAAAAAwHIUHAAAAAAAgOUoOAAAAAAAAMtRcAAAAAAAAJaj4AAAAAAAACxHwQEAAAAAAFiOggMAAAAAALAcBQcAAAAAAGA5Cg4AAAAAAMByFBwAAAAAAIDlKDgAAAAAAADLUXAAAAAAAACW+38PTmcaR2/46AAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "swiftdiff['vh'].plot(x=\"time\",hue=\"id\",col=\"space\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python (My debug_env Kernel)", + "language": "python", + "name": "debug_env" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/rmvs_swifter_comparison/1pl_1tp_encounter/swiftest_vs_swifter.py b/examples/rmvs_swifter_comparison/1pl_1tp_encounter/swiftest_vs_swifter.py new file mode 100644 index 000000000..dce7c1358 --- /dev/null +++ b/examples/rmvs_swifter_comparison/1pl_1tp_encounter/swiftest_vs_swifter.py @@ -0,0 +1,9 @@ +import swiftest +import numpy as np + +swiftersim = swiftest.Simulation(simdir="swifter_sim", read_data=True, codename="Swifter") +swiftestsim = swiftest.Simulation(simdir="swiftest_sim",read_data=True) +swiftestsim.data = swiftestsim.data.swap_dims({"name" : "id"}) +swiftdiff = swiftestsim.data - swiftersim.data +swiftdiff['rh'].plot(x="time",hue="id",col="space") +swiftdiff['vh'].plot(x="time",hue="id",col="space") \ No newline at end of file From fed8c213ea075e32a1948064f44352bb3ea311db Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 15 May 2023 16:31:19 -0400 Subject: [PATCH 2/3] More change to clean up Swifter file saving --- python/swiftest/swiftest/simulation_class.py | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index 2ff9c8454..804eee651 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -1441,13 +1441,19 @@ def ascii_file_input_error_msg(codename): if "IN_FORM" in self.param: init_cond_format = self.param['IN_FORM'] else: - init_cond_format = "EL" + if self.codename.title() == "Swiftest": + init_cond_format = "EL" + else: + init_cond_format = "XV" if init_cond_file_type is None: if "IN_TYPE" in self.param: init_cond_file_type = self.param['IN_TYPE'] else: - init_cond_file_type = "NETCDF_DOUBLE" + if self.codename.title() == "Swiftest": + init_cond_file_type = "NETCDF_DOUBLE" + else: + init_cond_file_type = "ASCII" if self.codename.title() == "Swiftest": init_cond_keys = ["CB", "PL", "TP"] @@ -1664,10 +1670,11 @@ def set_output_files(self, else: self.param['BIN_OUT'] = output_file_name - if output_format != "XV" and self.codename != "Swiftest": - warnings.warn(f"{output_format} is not compatible with {self.codename}. Setting to XV",stacklevel=2) - output_format = "XV" - self.param["OUT_FORM"] = output_format + if output_format is not None: + if output_format != "XV" and self.codename != "Swiftest": + warnings.warn(f"{output_format} is not compatible with {self.codename}. Setting to XV",stacklevel=2) + output_format = "XV" + self.param["OUT_FORM"] = output_format if self.restart: self.param["OUT_STAT"] = "APPEND" From ca9eaa4cab1e9b03c6a64ae5b0350c3962d090c1 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 15 May 2023 17:45:10 -0400 Subject: [PATCH 3/3] Put back some of the original operations from the Swifter code in the drift functions --- src/swiftest/swiftest_drift.f90 | 84 ++++++++++++++++----------------- 1 file changed, 42 insertions(+), 42 deletions(-) diff --git a/src/swiftest/swiftest_drift.f90 b/src/swiftest/swiftest_drift.f90 index 4c7302908..5c0427c62 100644 --- a/src/swiftest/swiftest_drift.f90 +++ b/src/swiftest/swiftest_drift.f90 @@ -251,22 +251,22 @@ pure subroutine swiftest_drift_kepmd(dm, es, ec, x, s, c) ! executable code fac1 = 1.0_DP / (1.0_DP - ec) q = fac1 * dm - fac2 = es**2 * fac1 - ec / 3.0_DP + fac2 = es*es*fac1 - ec / 3.0_DP x = q * (1.0_DP - 0.5_DP * fac1 * q * (es - q * fac2)) - y = x**2 + y = x*x s = x * (a0 - y * (a1 - y * (a2 - y * (a3 - y * (a4 - y))))) / a0 - c = sqrt(1.0_DP - s**2) + c = sqrt(1.0_DP - s*s) f = x - ec * s + es * (1.0_DP - c) - dm fp = 1.0_DP - ec * c + es * s fpp = ec * s + es * c fppp = ec * c - es * s dx = -f / fp - dx = -f / (fp + dx * fpp * 0.5_DP) - dx = -f / (fp + dx * fpp * 0.5_DP + dx**2* fppp * SIXTH) + dx = -f / (fp + dx * fpp/2.0_DP) + dx = -f / (fp + dx * fpp/2.0_DP + dx*dx * fppp * SIXTH) x = x + dx - y = x**2 + y = x*x s = x * (a0 - y * (a1 - y * (a2 - y * (a3 - y * (a4 - y))))) / a0 - c = sqrt(1.0_DP - s**2) + c = sqrt(1.0_DP - s*s) return end subroutine swiftest_drift_kepmd @@ -352,18 +352,18 @@ pure subroutine swiftest_drift_kepu_guess(dt, r0, mu, alpha, u, s) if (alpha > 0.0_DP) then if (dt / r0 <= thresh) then - s = dt / r0 - (dt**2 * u) / (2 * r0**3) + s = dt/r0 - (dt*dt*u)/(2.0_DP*r0*r0*r0) else - a = mu / alpha - en = sqrt(mu / a**3) - ec = 1.0_DP - r0 / a - es = u / (en * a**2) - e = sqrt(ec**2 + es**2) - y = en * dt - es + a = mu/alpha + en = sqrt(mu/(a*a*a)) + ec = 1.0_DP - r0/a + es = u/(en*a*a) + e = sqrt(ec*ec + es*es) + y = en*dt - es call swiftest_orbel_scget(y, sy, cy) - sigma = sign(1.0_DP, es * cy + ec * sy) - x = y + sigma * danbyk * e - s = x / sqrt(alpha) + sigma = sign(1.0_DP, es*cy + ec*sy) + x = y + sigma*DANBYK*e + s = x/sqrt(alpha) end if else call swiftest_drift_kepu_p3solve(dt, r0, mu, alpha, u, s, iflag) @@ -408,16 +408,16 @@ pure subroutine swiftest_drift_kepu_lag(s, dt, r0, mu, alpha, u, fp, c1, c2, c3, do nc = 0, ncmax x = s * s * alpha call swiftest_drift_kepu_stumpff(x, c0, c1, c2, c3) - c1 = c1 * s - c2 = c2 * s**2 - c3 = c3 * s**3 + c1 = c1*s + c2 = c2*s*s + c3 = c3*s*s*s f = r0 * c1 + u * c2 + mu * c3 - dt fp = r0 * c0 + u * c1 + mu * c2 fpp = (-r0 * alpha + mu) * c1 + u * c0 - ds = -ln * f / (fp + sign(1.0_DP, fp) * sqrt(abs((ln - 1)**2 * fp**2 - (ln - 1) * ln * f * fpp))) + ds = -ln*f/(fp + sign(1.0_DP, fp)*sqrt(abs((ln - 1.0_DP)*(ln - 1.0_DP)*fp*fp - (ln - 1.0_DP)*ln*f*fpp))) s = s + ds fdt = f / dt - if (fdt**2 < DANBYB**2) then + if (fdt*fdt < DANBYB*DANBYB) then iflag = 0 return end if @@ -454,23 +454,23 @@ pure subroutine swiftest_drift_kepu_new(s, dt, r0, mu, alpha, u, fp, c1, c2, c3, real(DP) :: x, c0, ds, f, fpp, fppp, fdt do nc = 0, 6 - x = s**2 * alpha + x = s*s*alpha call swiftest_drift_kepu_stumpff(x, c0, c1, c2, c3) - c1 = c1 * s - c2 = c2 * s**2 - c3 = c3 * s**3 - f = r0 * c1 + u * c2 + mu * c3 - dt - fp = r0 * c0 + u * c1 + mu * c2 - fpp = (-r0 * alpha + mu) * c1 + u * c0 - fppp = (-r0 * alpha + mu) * c0 - u * alpha * c1 - ds = -f / fp - ds = -f / (fp + ds * fpp * 0.5_DP) - ds = -f / (fp + ds * fpp * 0.5_DP + ds**2 * fppp * SIXTH) + c1 = c1*s + c2 = c2*s*s + c3 = c3*s*s*s + f = r0*c1 + u*c2 + mu*c3 - dt + fp = r0*c0 + u*c1 + mu*c2 + fpp = (-r0*alpha + mu)*c1 + u*c0 + fppp = (-r0*alpha + mu)*c0 - u*alpha*c1 + ds = -f/fp + ds = -f/(fp + ds*fpp/2.0_DP) + ds = -f/(fp + ds*fpp/2.0_DP + ds*ds*fppp/6.0_DP) s = s + ds - fdt = f / dt - if (fdt**2 < DANBYB**2) then - iflag = 0 - return + fdt = f/dt + if (fdt*fdt < DANBYB*DANBYB) then + iflag = 0 + return end if end do iflag = 1 @@ -503,9 +503,9 @@ pure subroutine swiftest_drift_kepu_p3solve(dt, r0, mu, alpha, u, s, iflag) a2 = 0.5_DP * u / denom a1 = r0 / denom a0 = -dt / denom - q = (a1 - a2**2 * THIRD) * THIRD - r = (a1 * a2 - 3 * a0) * SIXTH - (a2 * THIRD)**3 - sq2 = q**3 + r**2 + q = (a1 - a2*a2 * THIRD) * THIRD + r = (a1*a2 - 3 * a0) * SIXTH - (a2*a2*a2)/27.0_DP + sq2 = q*q*q + r*r if (sq2 >= 0.0_DP) then sq = sqrt(sq2) if ((r + sq) <= 0.0_DP) then @@ -565,9 +565,9 @@ pure subroutine swiftest_drift_kepu_stumpff(x, c0, c1, c2, c3) if (n /= 0) then do i = n, 1, -1 c3 = (c2 + c0 * c3) / 4.0_DP - c2 = c1**2 / 2.0_DP + c2 = c1*c1 / 2.0_DP c1 = c0 * c1 - c0 = 2 * c0**2 - 1.0_DP + c0 = 2 * c0*c0 - 1.0_DP x = x * 4 end do end if