Coupling sims work

3 files changed, 348 insertions, 81 deletions

diff --git a/Untitled.ipynb b/Untitled.ipynb
new file mode 100644
index 0000000..984cb8b
--- /dev/null
+++ b/Untitled.ipynb
@@ -0,0 +1,72 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "df53270d-9a10-4794-9560-4d168db3670b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from matplotlib import pyplot as plt\n",
+    "%matplotlib inline"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "id": "4e568696-b208-4ad0-b886-ed25addb7bec",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "Text(0.5, 0, 'd [mm]')"
+      ]
+     },
+     "execution_count": 8,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAAGzCAYAAAArAc0KAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABGf0lEQVR4nO3dd3hUZd7/8c9MykDaJIFUUiBUIaEjRVBQiaKCoKu4WLA8igo8ssLqqltwHwV117Y/7KvYwd0VXAu6oEIQqdKboQUSICEhCZNCMmnn90dg1iglA0nOZPJ+Xde5kpxzz5nvnOvofLjPfe5jMQzDEAAAgBezml0AAABAYyPwAAAAr0fgAQAAXo/AAwAAvB6BBwAAeD0CDwAA8HoEHgAA4PUIPAAAwOsReAAAgNcj8AAAAK/na+abv/LKK3rllVe0f/9+SVKPHj30xz/+UaNGjZIkGYahxx9/XK+//roKCws1cOBAvfTSS+rRo4drH06nUzNmzNC8efNUVlamyy67TC+//LLi4uLqXUdNTY0OHz6s4OBgWSyWBv2MAACgcRiGoeLiYsXGxspqPUsfjmGiTz/91Pjiiy+M9PR0Iz093Xj00UcNPz8/Y9u2bYZhGMZTTz1lBAcHGx9//LGxdetWY/z48UZMTIxRVFTk2se9995rtGvXzliyZImxYcMGY8SIEUavXr2MqqqqeteRlZVlSGJhYWFhYWFphktWVtZZv+sthuFZDw8NDw/XX/7yF915552KjY3VtGnT9PDDD0uq7c2JiorS008/rUmTJsnhcCgiIkLvvfeexo8fL0k6fPiw4uPjtWjRIl1xxRX1ek+Hw6HQ0FBlZWUpJCSk0T4bAABoOEVFRYqPj9exY8dkt9vP2NbUS1o/VV1drX/+858qLS3V4MGDlZGRoZycHKWmprra2Gw2XXLJJVq5cqUmTZqk9evXq7Kysk6b2NhYJScna+XKlacNPE6nU06n0/V3cXGxJCkkJITAAwBAM1Of4SimD1reunWrgoKCZLPZdO+992rhwoXq3r27cnJyJElRUVF12kdFRbm25eTkyN/fX2FhYadtcyqzZ8+W3W53LfHx8Q38qQAAgCcxPfB07dpVmzZt0urVq3Xfffdp4sSJ2rFjh2v7z1ObYRhnTXJna/PII4/I4XC4lqysrPP7EAAAwKOZHnj8/f3VqVMn9e/fX7Nnz1avXr304osvKjo6WpJ+0VOTm5vr6vWJjo5WRUWFCgsLT9vmVGw2m+vyFZexAADwfqYHnp8zDENOp1MdOnRQdHS0lixZ4tpWUVGhtLQ0DRkyRJLUr18/+fn51WmTnZ2tbdu2udoAAACYOmj50Ucf1ahRoxQfH6/i4mLNnz9fy5Yt01dffSWLxaJp06Zp1qxZ6ty5szp37qxZs2YpICBAEyZMkCTZ7Xbdddddmj59utq0aaPw8HDNmDFDKSkpuvzyy838aAAAwIOYGniOHDmiW2+9VdnZ2bLb7erZs6e++uorjRw5UpL00EMPqaysTPfff79r4sHFixcrODjYtY/nn39evr6+uvHGG10TD7799tvy8fEx62MBAAAP43Hz8JihqKhIdrtdDoeD8TwAADQT7nx/e9wYHgAAgIZG4AEAAF6PwAMAALwegQcAAHg9Ag8AAPB6BJ5GtutIsbIdZWaXAQBAi0bgaUSzF+1U6vPLNff7/WaXAgBAi0bgaUT9Emuf4v7x+oOqqKoxuRoAAFouAk8jurRbpCKDbcovrdDXO4+YXQ4AAC0WgacR+fpYdUP/OEnSvLWZJlcDAEDLReBpZOP7J0iSVuw5qqyC4yZXAwBAy0TgaWQJbQI0tFNbGYb0jx+yzC4HAIAWicDTBG66MF5SbeCpqmbwMgAATY3A0wRGdo9SeKC/jhQ5lbYrz+xyAABocQg8TcDm66Pr+rSTJM1by2UtAACaGoGniZy8rLU0PVdHispNrgYAgJaFwNNEOkUGa0D7MFXXGPong5cBAGhSBJ4mdNOA2lvUP/ohSzU1hsnVAADQchB4mtBVKTEKbuWrrIIyfb/3qNnlAADQYhB4mlBrfx+NOzF4ef46LmsBANBUCDxNbPyA2sHLi7fnKL/EaXI1AAC0DASeJtYj1q6ecXZVVhtasOGQ2eUAANAiEHhMcHLw8rx1mTIMBi8DANDYCDwmGNM7VgH+PtqXV6p1+wvNLgcAAK9H4DFBkM1Xo3vGSpLmr800uRoAALwfgcckJ2de/mJrthxllSZXAwCAdyPwmKR3fKi6RgXLWVWjf29i8DIAAI2JwGMSi8Xi6uWZtzaLwcsAADQiAo+JxvVpJ39fq3ZmF2nLQYfZ5QAA4LUIPCYKDfDXVcnRkqT56xi8DABAYyHwmOymC2vn5Pl002GVOqtMrgYAAO9E4DHZwA7hSmobqNKKan2+5bDZ5QAA4JUIPCazWCyu52vNW8sDRQEAaAwEHg9wXd84+Vot2pR1TD/mFJldDgAAXofA4wEigm0a2T1KkjSfXh4AABocgcdDnBy8vGDDQZVXVptcDQAA3oXA4yGGdWqrdqGtVVRepS+3ZZtdDgAAXoXA4yGs1v8OXuayFgAADYvA40Fu6B8nq0Vak1GgfXklZpcDAIDXIPB4kBh7aw3vGilJ+mgdvTwAADQUAo+HuenEZa1/rT+oiqoak6sBAMA7EHg8zKXdIhUZbFN+aYW+3nnE7HIAAPAKBB4P4+tj1Q394yRJ89byQFEAABoCgccDje9fOyfPij1HlVVw3ORqAABo/gg8HiihTYCGdmorw5D++QODlwEAOF8EHg9104W1g5f/8cNBVVUzeBkAgPNB4PFQI7tHKSzATzlF5UrblWd2OQAANGsEHg9l8/XR9X1PDl7mshYAAOeDwOPBTl7WWpqeqyNF5SZXAwBA80Xg8WCdIoM1oH2YqmsM/Wv9QbPLAQCg2SLweLibBtTeoj5/XaZqagyTqwEAoHki8Hi4q1JiFNzKV1kFZVq5N9/scgAAaJYIPB6utb+PxvZuJ0mat46ZlwEAOBcEnmbg5ODlxdtzlF/iNLkaAACaHwJPM9Aj1q6ecXZVVhuav45b1AEAcJepgWf27NkaMGCAgoODFRkZqbFjxyo9Pb1Om9tvv10Wi6XOMmjQoDptnE6npk6dqrZt2yowMFBjxozRwYPedVfT7UPaS5Lmfp+h8spqc4sBAKCZMTXwpKWlafLkyVq9erWWLFmiqqoqpaamqrS0tE67K6+8UtnZ2a5l0aJFdbZPmzZNCxcu1Pz587VixQqVlJTommuuUXW19wSD0b1i1S60tY6WVPB8LQAA3GQxDMNj7nXOy8tTZGSk0tLSdPHFF0uq7eE5duyYPvnkk1O+xuFwKCIiQu+9957Gjx8vSTp8+LDi4+O1aNEiXXHFFWd936KiItntdjkcDoWEhDTY52lo767arz/+e7viwlpr2Yzh8vXhiiQAoOVy5/vbo74xHQ6HJCk8PLzO+mXLlikyMlJdunTR3XffrdzcXNe29evXq7KyUqmpqa51sbGxSk5O1sqVK0/5Pk6nU0VFRXWW5uDG/vFqE+ivg4Vl+mzLYbPLAQCg2fCYwGMYhh588EENHTpUycnJrvWjRo3SBx98oG+//VbPPvus1q1bp0svvVROZ+3dSjk5OfL391dYWFid/UVFRSknJ+eU7zV79mzZ7XbXEh8f33gfrAG18vPRnUM7SJJeWbaXiQgBAKgnjwk8U6ZM0ZYtWzRv3rw668ePH6+rr75aycnJGj16tL788kvt2rVLX3zxxRn3ZxiGLBbLKbc98sgjcjgcriUrq/mMibllUKKCbL7adaRE3/6Ye/YXAAAAzwg8U6dO1aeffqqlS5cqLi7ujG1jYmKUmJio3bt3S5Kio6NVUVGhwsLCOu1yc3MVFRV1yn3YbDaFhITUWZoLe2s/3TIoUZL08rI98qAhWAAAeCxTA49hGJoyZYoWLFigb7/9Vh06dDjra/Lz85WVlaWYmBhJUr9+/eTn56clS5a42mRnZ2vbtm0aMmRIo9VupjuHtpe/r1UbMo9pTUaB2eUAAODxTA08kydP1vvvv68PP/xQwcHBysnJUU5OjsrKyiRJJSUlmjFjhlatWqX9+/dr2bJlGj16tNq2batx48ZJkux2u+666y5Nnz5d33zzjTZu3KhbbrlFKSkpuvzyy838eI0mMriVbuhX2xP28rK9JlcDAIDnMzXwvPLKK3I4HBo+fLhiYmJcy0cffSRJ8vHx0datW3XttdeqS5cumjhxorp06aJVq1YpODjYtZ/nn39eY8eO1Y033qiLLrpIAQEB+uyzz+Tj42PWR2t0ky7uKKtFWr4rT9sOOcwuBwAAj+ZR8/CYpbnMw/NzD8zfqH9vOqyre8bopQl9zS4HAIAm1Wzn4YF77hveUZL05dZsZRwtPUtrAABaLgJPM9YtOkSXdYtUjSG9lsZYHgAATofA08yd7OX5eMNB5TjKTa4GAADPROBp5vq3D9eF7cNVWW3ozRX7zC4HAACPRODxAveNqO3l+WBNpo4drzC5GgAAPA+BxwsM7xKhC2JCdLyiWu+sPGB2OQAAeBwCjxewWCyusTxvr8zQ8YoqkysCAMCzEHi8xFXJ0UpsE6DC45Wat7b5PAwVAICmQODxEr4+Vk26uLaX5+/f7VNFVY3JFQEA4DkIPF7k+n7tFBlsU7ajXJ9sOmR2OQAAeAwCjxex+frof4bVPnH+1bS9qq5p8U8NAQBAEoHH60wYmKiQVr7al1eqxdtzzC4HAACPQODxMkE2X00c0l6S9PKyveLZsAAAEHi80u1D2quVn1VbDzn0/Z58s8sBAMB0BB4v1CbIppsGJEiSXl62x+RqAAAwH4HHS919cZJ8rRat3JuvTVnHzC4HAABTEXi8VLvQ1rq2dztJ0stL6eUBALRsBB4vdt/wJFks0uIdR7T7SLHZ5QAAYBoCjxfrFBms1O5RkqRX0/aZXA0AAOYh8Hi5+4Z3kiT9e9MhHTpWZnI1AACYg8Dj5XrHh2pIxzaqqjH0xnJ6eQAALROBpwW4/0Qvz/x1mcovcZpcDQAATY/A0wJc1KmNesbZVV5Zo7nf7ze7HAAAmhyBpwWwWCy6f3hHSdI7q/aruLzS5IoAAGhaBJ4WIrV7tJIiAlVcXqUP12SaXQ4AAE2KwNNCWK0W3XtJbS/P31dkqLyy2uSKAABoOgSeFmRs73aKsbdSXrFTH284aHY5AAA0GQJPC+Lva9Xdw5IkSa+l7VNVdY3JFQEA0DQIPC3MTRfGKyzAT5kFx+nlAQC0GASeFibA31eTR9TOy/Ps4l06XlFlckUAADQ+Ak8LdOvgRMWHt1ZusVN//y7D7HIAAGh0BJ4WyObro99e0U2S9FraXuUVM/syAMC7EXhaqGtSYtQzzq7Simq9+M0us8sBAKBREXhaKKvVokevukCSNG9tlvbklphcEQAAjYfA04INSmqjyy+IVHWNoae/+tHscgAAaDQEnhbud6O6ycdq0ZIdR7Q2o8DscgAAaBQEnhauU2Swxg+IlyQ9uWinDMMwuSIAABoegQeadnlnBfj7aHPWMX2xNdvscgAAaHAEHigyuJXuubj2kRPPfJUuZxUPFgUAeBcCDyRJdw9LUkSwTZkFx/X+6kyzywEAoEEReCBJCrT56sGRXSRJ/+/b3XKUVZpcEQAADYfAA5cb+sWpc2SQjh2v1MvL9phdDgAADYbAAxdfH6t+N6r2kRNzv9+vg4XHTa4IAICGQeBBHZd2i9SgpHBVVNXo2cU8cgIA4B0IPKjDYrHosau6S5IWbjykbYccJlcEAMD5I/DgF1Li7Lq2d6wkaRaTEQIAvACBB6c0I7Wr/H2sWrk3X8t25ZldDgAA54XAg1OKDw/QxCGJkqSnFv2o6hp6eQAAzReBB6c1ZURn2Vv7Kf1Isf61PsvscgAAOGcEHpyWPcBPUy/tJEl6dvEuHa+oMrkiAADODYEHZ3Tr4ETFhbVWbrFTf/8uw+xyAAA4JwQenJHN10e/vaKrJOm1tL3KK3aaXBEAAO4j8OCsRveMVc84u0orqvXiN0xGCABofgg8OCur1aJHr7pAkjRvbZb25JaYXBEAAO4xNfDMnj1bAwYMUHBwsCIjIzV27Filp6fXaWMYhmbOnKnY2Fi1bt1aw4cP1/bt2+u0cTqdmjp1qtq2bavAwECNGTNGBw8ebMqP4vUGJbXR5RdEqrrG0NNf/Wh2OQAAuMXUwJOWlqbJkydr9erVWrJkiaqqqpSamqrS0lJXm2eeeUbPPfec5syZo3Xr1ik6OlojR45UcXGxq820adO0cOFCzZ8/XytWrFBJSYmuueYaVVdXm/GxvNbDV3aT1SIt2XFEazMKzC4HAIB6sxge9NyAvLw8RUZGKi0tTRdffLEMw1BsbKymTZumhx9+WFJtb05UVJSefvppTZo0SQ6HQxEREXrvvfc0fvx4SdLhw4cVHx+vRYsW6Yorrjjr+xYVFclut8vhcCgkJKRRP2Nz98iCrZq3NlO940O18P4hslgsZpcEAGih3Pn+9qgxPA5H7YMqw8PDJUkZGRnKyclRamqqq43NZtMll1yilStXSpLWr1+vysrKOm1iY2OVnJzsavNzTqdTRUVFdRbUz29GdlaAv482ZR3TF1uzzS4HAIB68ZjAYxiGHnzwQQ0dOlTJycmSpJycHElSVFRUnbZRUVGubTk5OfL391dYWNhp2/zc7NmzZbfbXUt8fHxDfxyvFRncSvdcnCRJeuardDmruGwIAPB8HhN4pkyZoi1btmjevHm/2PbzyyaGYZz1UsqZ2jzyyCNyOByuJSuLxya44+5hSYoItimz4LjeX51pdjkAAJyVRwSeqVOn6tNPP9XSpUsVFxfnWh8dHS1Jv+ipyc3NdfX6REdHq6KiQoWFhadt83M2m00hISF1FtRfoM1Xv7m8iyTp/327W46ySpMrAgDgzEwNPIZhaMqUKVqwYIG+/fZbdejQoc72Dh06KDo6WkuWLHGtq6ioUFpamoYMGSJJ6tevn/z8/Oq0yc7O1rZt21xt0PBu7B+nTpFBOna8Ui8v22N2OQAAnJGpgWfy5Ml6//339eGHHyo4OFg5OTnKyclRWVmZpNpLWdOmTdOsWbO0cOFCbdu2TbfffrsCAgI0YcIESZLdbtddd92l6dOn65tvvtHGjRt1yy23KCUlRZdffrmZH8+r+fpY9ciobpKkuSv2a/eR4rO8AgAA8/ia+eavvPKKJGn48OF11s+dO1e33367JOmhhx5SWVmZ7r//fhUWFmrgwIFavHixgoODXe2ff/55+fr66sYbb1RZWZkuu+wyvf322/Lx8Wmqj9IiXdotUiO6Rmhpep4e+niL/nXvEPlYuU0dAOB5PGoeHrMwD8+5O3ysTKnPL1eJs0p/uKa77hra4ewvAgCgATTbeXjQ/MSGttbvTlza+ut/0pWZf9zkigAA+CUCD87bhAsTNCgpXGWV1frdgi2i0xAA4GkIPDhvVqtFT1/fU638rFq5N1/z1zGvEQDAsxB40CAS2wRqRmpXSdKsL3Yq21FmckUAAPwXgQcN5o6LOqh3fKiKnVX6/cJtXNoCAHgMAg8ajI/Vomd+1VN+PhZ982OuPt182OySAACQROBBA+sSFaypl3aWJM38dLuOljhNrggAAAIPGsF9wzvqgpgQFR6v1J8+3W52OQAAEHjQ8Px8rPrLr3rKx2rRF1uy9Z/tOWd/EQAAjYjAg0aR3M6uey5OkiT9/pNtchznieoAAPMQeNBoHriss5IiApVX7NQTX+wwuxwAQAtG4EGjaeXno2eu7ymLRfrn+oNavivP7JIAAC0UgQeNqn/7cE0c3F6S9MiCrSp1VplbEACgRSLwoNH99oquigtrrUPHyvTMVz+aXQ4AoAUi8KDRBdp89dR1PSVJ76w6oLUZBSZXBABoaXzr0+jBBx90e8e///3vFR4e7vbr4J2Gdm6r8f3j9dEPWXr44y368oFhauXnY3ZZAIAWwmLU44FHVqtVgwcPlr+/f712umLFCqWnpyspKem8C2wKRUVFstvtcjgcCgkJMbscr+Uoq1Tq82k6UuTUpEuS9MioC8wuCQDQjLnz/V2vHh5JWrhwoSIjI+vVNjg4uL67RQtib+2nJ8am6O53f9Aby/fp6pQY9YwLNbssAEALUK8xPHPnzpXdbq/3Tl977TVFRUWdc1HwXiO7R2l0r1jVGNJD/9qiiqoas0sCALQA9bqk5e24pNW08kucGvn8chWUVug3l3fRA5d3NrskAEAz5M73N3dpocm1CbJp5pgekqQ5S3crPafY5IoAAN6u3mN4wsLCZLFYztquoIBbjnF2o3vG6NNNh/X1ziN66F+b9fF9Q+TrQ/4GADSOegeeF154wfW7YRi677779Oc//7neA5mBn7JYLHpyXLLWZORr80GH3vo+Q/dc3NHssgAAXuqcx/AEBwdr8+bNzebW8zNhDI95PlqXqYc/3iqbr1VfTbtYHdoGml0SAKCZYAwPmo0b+8draKe2clbV6OGPt6impsWPoQcANAICD0xlsVg0+7oUtfbz0dqMAn2w5oDZJQEAvBCBB6aLDw/QQ1d2lSQ98cVO/ZhTZHJFAABvU+9Byz9/nlZFRYWefPLJX0xI+NxzzzVMZWhRJg5ur2XpeUrblafJH2zQp1OGKtBW79MTAIAzqveg5REjRpx9ZxaLvv322/MuqqkxaNkz5Jc4dfXfViinqFxje8fq+fG96zUVAgCgZXLn+5uZlkXg8STr9hfoptdXq7rG0FPXpeimCxPMLgkA4KG4SwvN1oD24ZqRWjue50+fbtfObMbzAADOn9uDJKqrq/X222/rm2++UW5urmpq6j78sTle0oJnmXRxktZm5Gtp+onxPFOHKojxPACA8+B2D88DDzygBx54QNXV1UpOTlavXr3qLMD5slotevbG3oqxt9K+o6V6dMFWceUVAHA+3B7D07ZtW7377ru66qqrGqumJscYHs+0/kCBxr+2WlU1hp4cl6ybByaaXRIAwIM06hgef39/derU6ZyLA+qrX2K4a36exz/boW2HHCZXBABortwOPNOnT9eLL77IJQY0if8ZmqTLukWqoqpGUz7coOLySrNLAgA0Q25f0ho3bpyWLl2q8PBw9ejRQ35+fnW2L1iwoEELbApc0vJsx45X6Oq/rdChY2W6umeM5vy6D/PzAADc+v52+9aX0NBQjRs37pyLA9wVGuCv/zehj258dZW+2JKtQUltdOsgxvMAAOqPiQdFD09z8ffv9umJL3bK38eqBfcPUXI7+9lfBADwWkw8CK9019AOGtk9ShXVNbr/gw0qYjwPAKCe6hV4+vbtq8LCwnrvdOjQoTp06NA5FwWcisVi0V9/1UvtQlsrs+C4fvfxFgbPAwDqpV5jeDZt2qTNmzcrPDy8XjvdtGmTnE7neRUGnIo9wE8v3dxXN7y6Uou25ujdVQc0cUh7s8sCAHi4eg9avuyyy+r9r2nuoEFj6h0fqkdGXaA/f75DT3yxQ30SQtUzLtTssgAAHqxegScjI8PtHcfFxbn9GqC+7riovdZk5Os/249o8ocb9PnUYbK39jv7CwEALVK9Ak9iIrcAw7NYLBY986te2pH9nbIKyvTQvzbr1Vv60bsIADgl7tJCs2Vv7aeXJvSVv49V/9l+RHO/3292SQAAD0XgQbPWMy5Uj119gSRp9pc7tSnrmLkFAQA8EoEHzd5tgxN1VUq0KqsNTf5ggxzHmZ8HAFAXgQfNnsVi0VPX91RimwAdOlamGf/azPw8AIA6zivwlJSUqKioqM4CmCGk1X/H8yzZcURvrnD/zkIAgPdyO/BkZGTo6quvVmBgoOx2u8LCwhQWFqbQ0FCFhYU1Ro1AvSS3s+sP19SO53nqyx+1/kD9ZwcHAHg3t5+WfvPNN0uS3nrrLUVFRXEbMDzKLYMStSajQJ9vydY97/6ghfdfpIQ2AWaXBQAwmduBZ8uWLVq/fr26du3aGPUA58Visejp63sq42ipth8u0u1vr9WC+4YoNMDf7NIAACZy+5LWgAEDlJWV1Ri1AA0i0Oart24foBh7K+3LK9Wk99bLWVVtdlkAABNZDDdvZ9m7d6/uvfde3XLLLUpOTpafX93p/Hv27NmgBTaFoqIi2e12ORwOhYSEmF0OGsiPOUX61SurVOKs0rg+7fTcjb24BAsAXsSd72+3e3jy8vK0d+9e3XHHHRowYIB69+6tPn36uH66Y/ny5Ro9erRiY2NlsVj0ySef1Nl+++23y2Kx1FkGDRpUp43T6dTUqVPVtm1bBQYGasyYMTp48KC7HwteqFt0iF65pa98rRYt3HhIz3+92+ySAAAmcTvw3HnnnerTp49WrVqlffv2KSMjo85Pd5SWlqpXr16aM2fOadtceeWVys7Odi2LFi2qs33atGlauHCh5s+frxUrVqikpETXXHONqqu5hAFpWOcIPTkuWZL0t292658/cDkWAFoitwctHzhwQJ9++qk6dep03m8+atQojRo16oxtbDaboqOjT7nN4XDozTff1HvvvafLL79ckvT+++8rPj5eX3/9ta644orzrhHN3/gBCcosOK6Xlu7VIwu2Kja0tS7q1NbssgAATcjtHp5LL71UmzdvboxaTmnZsmWKjIxUly5ddPfddys3N9e1bf369aqsrFRqaqprXWxsrJKTk7Vy5crT7tPpdDJhYgszfWRXjekVq6oaQ/e+t17pOcVmlwQAaEJu9/CMHj1av/nNb7R161alpKT8YtDymDFjGqy4UaNG6YYbblBiYqIyMjL0hz/8QZdeeqnWr18vm82mnJwc+fv7/2LCw6ioKOXk5Jx2v7Nnz9bjjz/eYHXC81mtFv3lhp7KcZRr7f4C3fn2Oi28f4giQ1qZXRoAoAm4fZeW1Xr6TiGLxXLOY2csFosWLlyosWPHnrZNdna2EhMTNX/+fF133XX68MMPdccdd8jpdNZpN3LkSHXs2FGvvvrqKffjdDrrvKaoqEjx8fHcpdUCFJZW6PpXVmrf0VKltLPro0mDFODvdu4HAHiARr1Lq6am5rRLYw8UjomJUWJionbvrr3bJjo6WhUVFSosrPsIgdzcXEVFRZ12PzabTSEhIXUWtAxhgf6ae8cAhQf6a+shh/533kZV1/CgUQDwdo32tPSUlJQGn6AwPz9fWVlZiomJkST169dPfn5+WrJkiatNdna2tm3bpiFDhjToe8N7JLYJ1Bu39ZfN16qvd+bqz59t5+nqAODlGi3w7N+/X5WVlWdsU1JSok2bNmnTpk2Sah9MumnTJmVmZqqkpEQzZszQqlWrtH//fi1btkyjR49W27ZtNW7cOEmS3W7XXXfdpenTp+ubb77Rxo0bdcsttyglJcV11xZwKv0Sw/T8+N6SpHdWHdBb3+83tR4AQONqtMBTHz/88IP69OnjmrDwwQcfVJ8+ffTHP/5RPj4+2rp1q6699lp16dJFEydOVJcuXbRq1SoFBwe79vH8889r7NixuvHGG3XRRRcpICBAn332mXx8fMz6WGgmrkqJ0aNXdZMkPfHFDn217fQD3QEAzZvbg5brKzg4WJs3b1ZSUlJj7L5B8WiJlsswDP3h39v0/upMtfKzav49g9U7PtTssgAA9dCog5YBb2KxWDRzdA+N6Bqh8soa/c8765RVcNzssgAADYzAgxbP18eqORP6qkdsiI6WVOj2uWvlOH7m8WcAgOaFwANICrT56q3bByjG3kp780o16f0f5KzieWwA4C0aLPBkZWXpzjvvdP392muvnXEuHMDTRIW00tw7BijI5qvV+wr0u4+3crs6AHiJBgs8BQUFeuedd1x/T5gwQYGBgQ21e6BJdIsO0cs395WP1aKFGw/p+a93m10SAKABcEkL+JmLu0ToybHJkqS/fbNb/1p/0OSKAADni8ADnMJNFyZo8oiOkqTffbxF3+w8YnJFAIDzQeABTmP6yK66tnesqmoM3fv+ev1nOxMTAkBzVe/HRF933XVn3H7s2LHzrQXwKFarRX+9oZeqawx9viVbkz/YoL/9uo+uSokxuzQAgJvqHXjsdvtZt992223nXRDgSfx8rHphfG/5+Vi1cOMhTZ23UZXVNbq2dzuzSwMAuKHegWfu3LmNWQfgsXx9rPrrDb3kY7XoX+sP6jcfbVJ1jaHr+saZXRoAoJ4YwwPUg4/Vomeu76lfXxivGkOa/s/N+se6LLPLAgDUE4EHqCer1aInx6bo1kGJMgzpoY+36IM1B8wuCwBQDwQewA1Wq0V/vraH7riovSTpsYXb9O6q/abWBAA4OwIP4CaLxaI/XtNd91ycJEn647+36+/f7TO5KgDAmRB4gHNgsVj0yKhurskJn/hip15N22tyVQCA0yHwAOfIYrFoRmpXPXBZZ0nSU1/+qDnf8uwtAPBEBB7gPFgsFv1mZBfNSO0iSfrr4l16fskunrIOAB6GwAM0gCmXdtbvRnWTJL34zW79dXE6oQcAPAiBB2gg917SUb+/+gJJ0ktL92r2lz8SegDAQxB4gAb0P8OS9PiYHpKk15fv058/30HoAQAPQOABGtjEIe315LhkSdLc7/frj//erpoaQg8AmInAAzSCmwcm6pnre8pikd5bfUCPfbKV0AMAJiLwAI3kxgHxevaGXrJapHlrs/TQx1tUTegBAFMQeIBGdF3fOD0/vrfrSevT/7FJldU1ZpcFAC0OgQdoZNf2bqe/3dRHvlaLPtl0WLfPXStHWaXZZQFAi0LgAZrA1T1j9Ppt/RTg76Pv9+TrV6+sVFbBcbPLAoAWg8ADNJFLu0XpH5MGKyrEpt25JRr38kptyjpmdlkA0CIQeIAmlNzOrk8mX6QLYkJ0tMSpm15fpa+25ZhdFgB4PQIP0MRi7K31z3sHa0TXCJVX1ui+D9brjeX7mKAQABoRgQcwQZDNV2/c1l+3DkqUYUhPLtqp33+yTVXcwQUAjYLAA5jE18eqP1/bQ7+/+gJZLNIHazJ11zs/qLicO7gAoKEReAATWSwW/c+wJL16Sz+18rMqbVeebnh1lQ4fKzO7NADwKgQewANc0SNa/5g0WG2DbPoxp1jjXv5e2w45zC4LALwGgQfwED3jQvXJ5CHqEhWkI0VO3fjaKn2z84jZZQGAVyDwAB4kLixA/7pviIZ1bqvjFdW6+90f9Pb3GWaXBQDNHoEH8DAhrfz01u0DdNOAeNUY0szPdmjmp9t58CgAnAcCD+CB/Hysmn1din43qpsk6e2V+zXpvR9U6qwyuTIAaJ4IPICHslgsuveSjnppQl/5+1r19c5cjX99lY4UlZtdGgA0OwQewMNd3TNG8+4epDaB/tp2qEhjX/peO7OLzC4LAJoVAg/QDPRLDNPC+y9Sx4hAZTvKdcOrq7Q0PdfssgCg2SDwAM1EQpsALbjvIg1OaqMSZ5XufHudnl2czuMoAKAeCDxAM2IP8NM7d16oCQMTZBjS//t2j27++xrG9QDAWRB4gGbG39eqWeNS9OJNvRXo76M1GQW66sXvtHxXntmlAYDHIvAAzdS1vdvps6lDdUFMiPJLKzRx7lr99T9c4gKAUyHwAM1YUkSQFt4/RDefuMQ1Z+keTXhjjXIcXOICgJ8i8ADNXCs/Hz05LkVzJvRRkM1Xa/cX6Kq/fadl3MUFAC4EHsBLXNMzVp9PHaoesSEqKK3Q7XPX6emvfuQSFwCIwAN4lfZtA/XxfUN066BESdIry/bq12+sVrajzOTKAMBcBB7Ay7Ty89H/jU3WSxP6Ksjmq3X7C3XVi98xUSGAFo3AA3ipq3vG6POpQ5XcLkSFxyt1x9x1eurLH1XJJS4ALRCBB/BiJy9xTRxce4nr1bS9uun11Tp8jEtcAFoWAg/g5Wy+Pnr82mS9fHNfBdt8tf5Aoa7623f69scjZpcGAE2GwAO0EFelxOjz/x2qlHZ2HTteqTvf/kGzF+3kEheAFoHAA7QgiW0C9a/7Buv2Ie0lSa8t36fxr63SIS5xAfBypgae5cuXa/To0YqNjZXFYtEnn3xSZ7thGJo5c6ZiY2PVunVrDR8+XNu3b6/Txul0aurUqWrbtq0CAwM1ZswYHTx4sAk/BdC82Hx9NHNMD716S18Ft/LVhsxjGvXCcn2+5bDZpQFAozE18JSWlqpXr16aM2fOKbc/88wzeu655zRnzhytW7dO0dHRGjlypIqLi11tpk2bpoULF2r+/PlasWKFSkpKdM0116i6urqpPgbQLF2ZHKMvpg5Tr/hQFZVXacqHGzXjn5tV4qwyuzQAaHAWwzAMs4uQJIvFooULF2rs2LGSant3YmNjNW3aND388MOSantzoqKi9PTTT2vSpElyOByKiIjQe++9p/Hjx0uSDh8+rPj4eC1atEhXXHFFvd67qKhIdrtdDodDISEhjfL5AE9VWV2jv32zWy8t3aMaQ0psE6AXxvdWn4Qws0sDgDNy5/vbY8fwZGRkKCcnR6mpqa51NptNl1xyiVauXClJWr9+vSorK+u0iY2NVXJysqvNqTidThUVFdVZgJbKz8eq6aldNf+ewWoX2loH8o/rV6+u0pxvd6u6xiP+PQQA581jA09OTo4kKSoqqs76qKgo17acnBz5+/srLCzstG1OZfbs2bLb7a4lPj6+gasHmp8LO4Rr0QPDNLpXrKprDP118S79+vXVDGgG4BU8NvCcZLFY6vxtGMYv1v3c2do88sgjcjgcriUrK6tBagWaO3trP/3tpt569oZeCvT30dr9BbryheX6bDMDmgE0bx4beKKjoyXpFz01ubm5rl6f6OhoVVRUqLCw8LRtTsVmsykkJKTOAqCWxWLR9f3itOiBYeqTEKri8ipNnbdRD/5jEwOaATRbHht4OnTooOjoaC1ZssS1rqKiQmlpaRoyZIgkqV+/fvLz86vTJjs7W9u2bXO1AXBuEtsE6h+TBut/L+0kq0VasOGQrnrxO23ILDz7iwHAw/ia+eYlJSXas2eP6++MjAxt2rRJ4eHhSkhI0LRp0zRr1ix17txZnTt31qxZsxQQEKAJEyZIkux2u+666y5Nnz5dbdq0UXh4uGbMmKGUlBRdfvnlZn0swGv4+Vj1YGpXDesSoWnzNymz4LhueHWVHrissyaP6CQf65kvLwOApzD1tvRly5ZpxIgRv1g/ceJEvf322zIMQ48//rhee+01FRYWauDAgXrppZeUnJzsalteXq7f/va3+vDDD1VWVqbLLrtML7/8slsDkbktHTg7R1ml/vDJNn16YjzPgPZhen58b8WFBZhcGYCWyp3vb4+Zh8dMBB6g/hZuPKg/fLJdJc4qBbfy1ZPjUjSmV6zZZQFogbxiHh4Anmlcnzgt+t9h6ntiQPP/ztuoBz/apOLySrNLA4DTIvAAcFtCmwD9Y9JgPXBZ59oBzRsP6aq/faf1BxjQDMAzEXgAnBNfH6t+M7KL/jFpsOLCWiuroEw3vrZK//f5DjnK6O0B4FkIPADOS//2tTM0j+1dO0PzmysydOlfl2ne2kweTQHAYzBoWQxaBhpK2q48/fmz7dqbVypJ6hEbopljemhA+3CTKwPgjbhLy00EHqDhVFbX6N1VB/TC17tUXF47M/PoXrF6ZFQ3xYa2Nrk6AN6EwOMmAg/Q8PJLnPrr4l2avy5ThiG18rPqvks6adIlSWrl52N2eQC8AIHHTQQeoPFsO+TQnz/bobX7CyRJ7UJb67GrL9Co5OizPggYAM6EwOMmAg/QuAzD0OdbsjV70U4ddpRLkgYlhetPo3voghj+mwNwbgg8biLwAE2jrKJar6bt1atpe+WsqpHVIv36wgRNT+2q8EB/s8sD0MwQeNxE4AGa1sHC45r95Y/6Yku2JCmkla9+M7KLbhmUKD8fZssAUD8EHjcReABzrN6Xr8c/26Gd2UWSpM6RQfrj6O4a1jnC5MoANAcEHjcReADzVNcYmr8uU3/9T7oKj9fO0Dyye5R+f/UFSmwTaHJ1ADwZgcdNBB7AfI7jlXrhm116d9UBVdcY8vOxaEyvdrrn4iR1jQ42uzwAHojA4yYCD+A5dh8p1p8/36Hvdh91rbukS4QmXZykwR3bcCs7ABcCj5sIPIDn2ZhZqDe+26evtuXo5CO5esSG6J6Lk3RVSgyDmwEQeNxF4AE814H8Ur25IkP/+CFL5ZU1kmonL7zjova66cIEBdl8Ta4QgFkIPG4i8ACer7C0Qu+vPqB3Vu3X0ZIKSVJwK1/dPDBRd1zUXlEhrUyuEEBTI/C4icADNB/lldVauPGQ3vhun/adeCq7n49F1/Zup7uHMcAZaEkIPG4i8ADNT02NoW9+zNUby/e5ntMlMcAZaEkIPG4i8ADNGwOcgZaJwOMmAg/gHc40wPnGAfEKaeVncoUAGhKBx00EHsC7nGqAc5DNVzf0j9MdQzoooU2AyRUCaAgEHjcReADvVF5ZrQUbDumt7zO0J7dEkmS11D664q6hSRrQPoxxPkAzRuBxE4EH8G41NYaW787Tmysy6szgnNLOrjuHttfVKbHy92WcD9DcEHjcROABWo5dR4o19/sMLdhwSM6q2nE+USE23Ta4vSZcmKCwQH+TKwRQXwQeNxF4gJYnv8SpD9dk6t3VB5RX7JQktfKz6rq+cbrzovbqFMl8PoCnI/C4icADtFzOqmp9sSVbb67I0PbDRa71l3SJ0F1DO2hY57aM8wE8FIHHTQQeAIZhaE1Ggd5ckaGvdx7Ryf8zdokK0p0XddDYPu3Uys/H3CIB1EHgcROBB8BP7T9aqrdX7tc/f8hSaUW1JCk80F83D0zQrYMSFclzuwCPQOBxE4EHwKk4yir1j3VZenvlfh06ViZJ8rVaNLJ7lG4emKghHdvIauVyF2AWAo+bCDwAzqSqukaLdxzRWysy9MOBQtf69m0C9OsLE/SrfnFqE2QzsUKgZSLwuInAA6C+fswp0odrMrVgwyGVOKskSf4+Vo1KidbNAxOZzBBoQgQeNxF4ALir1FmlzzYf1gdrMrX1kMO1vnNkkG4emKBxfeNkb82zu4DGROBxE4EHwPnYcvCYPlyTqX9vOqyyytpBzq38rBrdM1Y3D0pUrzg7vT5AIyDwuInAA6AhFJVX6pONh/TB6kylHyl2re8RG6KbBybq2t6xCrT5mlgh4F0IPG4i8ABoSIZhaP2BQn2wJlNfbM1WxYlHWATZfDW2T6wmXJio7rH8vwY4XwQeNxF4ADSWwtIKfbzhoD5Yk6mMo6Wu9X0SQvWrfnEalRyjcJ7fBZwTAo+bCDwAGpthGFq1N18frMnUf7bnqKqm9n+9PlaLLurUVqN7xii1RzQDnQE3EHjcROAB0JRyi8u1cMMhfbblsLYd+u/zu/x9rLqka4Su6Rmjyy+IYrwPcBYEHjcReACYZV9eiT7fkq3PNh/W7twS1/pWflZddkGURveM0fCukTzHCzgFAo+bCDwAPEF6TrE+23xYn285rP35x13rg2y+Su0epdG9YnVRp7by97WaWCXgOQg8biLwAPAkhmFo26EifbblsD7ffFiHHeWubfbWfhqVHK3RvWI1KKmNfHiWF1owAo+bCDwAPFVNjaGNWYX6bHO2Pt+SraMlTte2tkE2XZVSG376JoQRftDiEHjcROAB0BxU1xhak5GvzzZn68tt2Tp2vNK1rW2Qvy7rFqWR3aM0tHNbxvygRSDwuInAA6C5qayu0Yo9R/X55mwt3pGj4vIq17ZWflYN6xyhkd2jdFm3SJ7kDq9F4HETgQdAc1ZRVaO1GQX6eucRLdlxRIeOlbm2WSxSv4Qwjexe2/uTFBFkYqVAwyLwuInAA8BbGIahHdlFWrLjiL7eeaTOPD+S1DEiUJd3j1Jq9yj1jmfcD5o3Ao+bCDwAvNXhY2Wunp/V+/JVWf3f/+Uz7gfNHYHHTQQeAC1BUXml0tLztGTHES1Nzz31uJ8LojS8a4QiQ1qZWClQPwQeNxF4ALQ0Zxr3I0k9YkM0vGuEhneNVJ/4UPn6MNkhPA+Bx00EHgAt2U/H/Sz9MVebDzrqbA9p5athXSI0omukLukSoYhg7vqCZyDwuInAAwD/dbTEqeW78rQsPU/Ld+fVme9HklLa2U/0/kQw8BmmIvC4icADAKdWXWNoU9YxLUvP1bL0PG09VLf3JzTAT8M6R2hE1whd3CVCbZnzB02IwOMmAg8A1E9esVNpu/K0LD1X3+0+KkfZf3t/LJaTvT+RGt41Qr3iQun9QaPyqsAzc+ZMPf7443XWRUVFKScnR1LttefHH39cr7/+ugoLCzVw4EC99NJL6tGjR73fg8ADAO6rqq450fuTp6Xpudp+uO6cPyGtfDWgfbgGJoVrYIc26hEbwuBnNCh3vr99m6im89KjRw99/fXXrr99fP47V8Qzzzyj5557Tm+//ba6dOmiJ554QiNHjlR6erqCg4PNKBcAWgRfH6v6tw9X//bhmnFFV+UWlZ/o/akd+1NUXqVvfszVNz/mSpIC/X3Ur324BnaoXXrGhcrflwCEptEsAo+vr6+io6N/sd4wDL3wwgt67LHHdN1110mS3nnnHUVFRenDDz/UpEmTmrpUAGixIkNa6Yb+8bqhf7yqqmu0I7tIa/YVaE1GvtZmFKiovErLd+Vp+a48SbVz//RNCNOFHWp7gPokhDL5IRpNswg8u3fvVmxsrGw2mwYOHKhZs2YpKSlJGRkZysnJUWpqqqutzWbTJZdcopUrV5428DidTjmdTtffRUVFp2wHADg3vj5W9YwLVc+4UN19cZJqagz9mFOsNRn5WrOvQGv3F6igtEIr9+Zr5d58Sbvl72NVr3i7BnZoo4FJ4eqbEKZAW7P4mkIz4PFn0sCBA/Xuu++qS5cuOnLkiJ544gkNGTJE27dvd43jiYqKqvOaqKgoHThw4LT7nD179i/GBQEAGo/ValH32BB1jw3RHRd1kGEY2pNbotUZBVqbUaA1+/KVW+zUuv2FWre/UHOWSr5Wi5Lb2TWwQ+1ls74JoTz5HefM4wct/1xpaak6duyohx56SIMGDdJFF12kw4cPKyYmxtXm7rvvVlZWlr766qtT7uNUPTzx8fEMWgYAkxiGof35x7VmX+3lrzUZBb+Y/VmS2rcJUN+EMPVNDFPfhDB1jQ7mTrAWzOsGLf9UYGCgUlJStHv3bo0dO1aSlJOTUyfw5Obm/qLX56dsNptsNv6VAACewmKxqEPbQHVoG6ibLkyQJGUVHD8RfvK1IfOY9uSWaH/+ce3PP64FGw9Jqh0I3TshtDYEJYSpT0KoQgP8zfwo8FDNLvA4nU7t3LlTw4YNU4cOHRQdHa0lS5aoT58+kqSKigqlpaXp6aefNrlSAMD5iA8PUHx4gK7vFydJOna8QhuzjmnjgUJtyDymjZmFKq2o1vd78vX9nnzX6zpGBKpvQpj6Jdb2BHWKCJKVXqAWz+MDz4wZMzR69GglJCQoNzdXTzzxhIqKijRx4kRZLBZNmzZNs2bNUufOndW5c2fNmjVLAQEBmjBhgtmlAwAaUGiAv0Z0jdSIrpGSameB3nWkWBsyC7X+QKE2Zh5TxtFS7c2rXf65/qAkKbiVr/okhKnviZ6gXnGhsgf4mflRYAKPDzwHDx7Ur3/9ax09elQREREaNGiQVq9ercTEREnSQw89pLKyMt1///2uiQcXL17MHDwA4OV8rBZdEBOiC2JCdPPA2u+EgtIKbTwRgDZkFmpzlkPFP7sdXpIS2wQopZ1dPePsSmkXquR2IQpuRQjyZs1u0HJjYKZlAPBOVdU1+jGnthdow4FCbcw6pgP5x0/ZNikiUCnt7CeCUKh6xIZwW7yH86pHSzQFAg8AtBzHjldo26EibTl0TFsPOrTloOOUd4RZLFKniCClxNnVs51dKXGh6h4Totb+TI7oKQg8biLwAEDLll/i1NZDDm07VBuAth5yKNtR/ot2PlaLOkcGnbgUZlfX6BB1jQ6WvTWXw8xA4HETgQcA8HO5xeX/DUAHHdp80KGjJc5Tto2xt1KXqGB1iw5W1xNLx4ggHpXRyAg8biLwAADOxjAMHSlyasvBY67eoF1HSk55OUyq7Q1q3yZA3U70AnWNDlbXqGAlhAdwm3wDIfC4icADADhXReWV2pVTrB9zirXrSO3P9JxiOcoqT9m+tZ+PukQFqWt08IleoRB1iQ5SRJBNFgtByB0EHjcReAAADelkb1D6kWKl5xS5QtDu3BJVVNWc8jXBNl8lRQapY9tAJUUEqmNEkJIigpTYJoBLY6dB4HETgQcA0BSqqmt0oOC40k/2COUUK/1IsQ7kl6rmNN/GVosUFxbwkxAUqKS2QeoYGdjie4UIPG4i8AAAzOSsqtaB/OPam1uifUdLtTe3RHuPlmpfXomKy6tO+7pgm2/dIHTiZ/s2gS2iV4jA4yYCDwDAExmGobwSp/bllWpvXon25dWGoL15pTpYePy0vUKSFGtvpQ4RtQ9kbd8m0BWE4sMD5OdjbboP0YgIPG4i8AAAmpvyytpeoX15P+kVOvH7mXqFfKwWJYQHqH2bAHVoG1QbitoEqkNEoGJCWjWrO8gIPG4i8AAAvIVhGCoorVDG0VLXsj+/VPvyan+WV5560LQk2Xytat/mRK9Q20AlhAcoPry14sICFBvaSjZfz7pM5s73Nw8JAQDAi1gsFrUJsqlNkE3924fX2VZTY+hIcbky8kqVkV9a+/No7e+Z+cflrKqpvbPsSPEp9itFBtsUHxaguLDaEBQX1lrx4bU/Y+yt5e/ruZfK6OERPTwAAFRV1+jQsTLtO1obhPbnlyqr4LgOFpbpYGGZyiqrz/h6i0WKDmn1k0B0IhSFt1Z8WICi7a0afOwQl7TcROABAOD0Tl4mO1hYpqzCkyGo9ufJUOQ8zfxCJ/36wgTNvi6lQevikhYAAGgwP71M1is+9BfbDcPQ0ZKK/4agwv/2DJ1cFxfWuukL/wkCDwAAOC8Wi0URwTZFBNvUJyHsF9tragxV1py5B6ixEXgAAECjslotslnNvcPLc4dTAwAANBACDwAA8HoEHgAA4PUIPAAAwOsReAAAgNcj8AAAAK9H4AEAAF6PwAMAALwegQcAAHg9Ag8AAPB6BB4AAOD1CDwAAMDrEXgAAIDX42npkgzDkCQVFRWZXAkAAKivk9/bJ7/Hz4TAI6m4uFiSFB8fb3IlAADAXcXFxbLb7WdsYzHqE4u8XE1NjQ4fPqzg4GBZLBYVFRUpPj5eWVlZCgkJMbu8FoPjbg6Ouzk47ubguDe9xjzmhmGouLhYsbGxslrPPEqHHh5JVqtVcXFxv1gfEhLCfxAm4Libg+NuDo67OTjuTa+xjvnZenZOYtAyAADwegQeAADg9Qg8p2Cz2fSnP/1JNpvN7FJaFI67OTju5uC4m4Pj3vQ85ZgzaBkAAHg9engAAIDXI/AAAACvR+ABAABej8ADAAC8HoHnZ15++WV16NBBrVq1Ur9+/fTdd9+ZXZJXmzlzpiwWS50lOjra7LK8zvLlyzV69GjFxsbKYrHok08+qbPdMAzNnDlTsbGxat26tYYPH67t27ebU6wXOdtxv/32239x/g8aNMicYr3I7NmzNWDAAAUHBysyMlJjx45Venp6nTac8w2vPsfdzHOewPMTH330kaZNm6bHHntMGzdu1LBhwzRq1ChlZmaaXZpX69Gjh7Kzs13L1q1bzS7J65SWlqpXr16aM2fOKbc/88wzeu655zRnzhytW7dO0dHRGjlypOs5czg3ZzvuknTllVfWOf8XLVrUhBV6p7S0NE2ePFmrV6/WkiVLVFVVpdTUVJWWlrracM43vPocd8nEc96Ay4UXXmjce++9ddZ169bN+N3vfmdSRd7vT3/6k9GrVy+zy2hRJBkLFy50/V1TU2NER0cbTz31lGtdeXm5YbfbjVdffdWECr3Tz4+7YRjGxIkTjWuvvdaUelqS3NxcQ5KRlpZmGAbnfFP5+XE3DHPPeXp4TqioqND69euVmppaZ31qaqpWrlxpUlUtw+7duxUbG6sOHTropptu0r59+8wuqUXJyMhQTk5OnXPfZrPpkksu4dxvAsuWLVNkZKS6dOmiu+++W7m5uWaX5HUcDockKTw8XBLnfFP5+XE/yaxznsBzwtGjR1VdXa2oqKg666OiopSTk2NSVd5v4MCBevfdd/Wf//xHb7zxhnJycjRkyBDl5+ebXVqLcfL85txveqNGjdIHH3ygb7/9Vs8++6zWrVunSy+9VE6n0+zSvIZhGHrwwQc1dOhQJScnS+KcbwqnOu6Suec8T0v/GYvFUudvwzB+sQ4NZ9SoUa7fU1JSNHjwYHXs2FHvvPOOHnzwQRMra3k495ve+PHjXb8nJyerf//+SkxM1BdffKHrrrvOxMq8x5QpU7RlyxatWLHiF9s45xvP6Y67mec8PTwntG3bVj4+Pr9I97m5ub/4VwAaT2BgoFJSUrR7926zS2kxTt4Vx7lvvpiYGCUmJnL+N5CpU6fq008/1dKlSxUXF+dazznfuE533E+lKc95As8J/v7+6tevn5YsWVJn/ZIlSzRkyBCTqmp5nE6ndu7cqZiYGLNLaTE6dOig6OjoOud+RUWF0tLSOPebWH5+vrKysjj/z5NhGJoyZYoWLFigb7/9Vh06dKiznXO+cZztuJ9KU57zXNL6iQcffFC33nqr+vfvr8GDB+v1119XZmam7r33XrNL81ozZszQ6NGjlZCQoNzcXD3xxBMqKirSxIkTzS7Nq5SUlGjPnj2uvzMyMrRp0yaFh4crISFB06ZN06xZs9S5c2d17txZs2bNUkBAgCZMmGBi1c3fmY57eHi4Zs6cqeuvv14xMTHav3+/Hn30UbVt21bjxo0zsermb/Lkyfrwww/173//W8HBwa6eHLvdrtatW8tisXDON4KzHfeSkhJzz3lT7g3zYC+99JKRmJho+Pv7G3379q1zOx0a3vjx442YmBjDz8/PiI2NNa677jpj+/btZpfldZYuXWpI+sUyceJEwzBqb9P905/+ZERHRxs2m824+OKLja1bt5pbtBc403E/fvy4kZqaakRERBh+fn5GQkKCMXHiRCMzM9Psspu9Ux1zScbcuXNdbTjnG97ZjrvZ57zlRJEAAABeizE8AADA6xF4AACA1yPwAAAAr0fgAQAAXo/AAwAAvB6BBwAAeD0CDwAA8HoEHgAA4PUIPAA8zvDhwzVt2rTTbp85c6YsFossFoteeOGFJqvrp9q3b++q4dixY6bUAKD+CDwAmqUePXooOztb99xzjynvv27dOn388cemvDcA9/HwUADNkq+vr6Kjo017/4iICIWHh5v2/gDcQw8PAFOVlpbqtttuU1BQkGJiYvTss8+e874sFotee+01XXPNNQoICNAFF1ygVatWac+ePRo+fLgCAwM1ePBg7d271/WamTNnqnfv3nrrrbeUkJCgoKAg3XfffaqurtYzzzyj6OhoRUZG6sknn2yIjwvAJAQeAKb67W9/q6VLl2rhwoVavHixli1bpvXr15/z/v7v//5Pt912mzZt2qRu3bppwoQJmjRpkh555BH98MMPkqQpU6bUec3evXv15Zdf6quvvtK8efP01ltv6eqrr9bBgweVlpamp59+Wr///e+1evXq8/qsAMzDJS0ApikpKdGbb76pd999VyNHjpQkvfPOO4qLizvnfd5xxx268cYbJUkPP/ywBg8erD/84Q+64oorJEkPPPCA7rjjjjqvqamp0VtvvaXg4GB1795dI0aMUHp6uhYtWiSr1aquXbvq6aef1rJlyzRo0KBzrg2AeQg8AEyzd+9eVVRUaPDgwa514eHh6tq16znvs2fPnq7fo6KiJEkpKSl11pWXl6uoqEghISGSau+4Cg4OrtPGx8dHVqu1zrrc3NxzrguAubikBcA0hmE0+D79/Pxcv1ssltOuq6mpOeVrTrY51bqfvgZA80LgAWCaTp06yc/Pr87YmMLCQu3atcvEqgB4Iy5pATBNUFCQ7rrrLv32t79VmzZtFBUVpccee6zOpSQAaAgEHgCm+stf/qKSkhKNGTNGwcHBmj59uhwOh9llAfAyFqMxLqIDQCOaOXOmPvnkE23atMnUOpYtW6YRI0aosLBQoaGhptYC4MzoNwbQLG3dulVBQUF6+eWXTXn/Hj16aNSoUaa8NwD30cMDoNkpKChQQUGBpNpHPNjt9iav4cCBA6qsrJQkJSUlMe4I8HAEHgAA4PX4JwkAAPB6BB4AAOD1CDwAAMDrEXgAAIDXI/AAAACvR+ABAABej8ADAAC8HoEHAAB4vf8PH8tzdEbcpx4AAAAASUVORK5CYII=",
+      "text/plain": [
+       "<Figure size 640x480 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "x = list(range(1, 26))\n",
+    "l = [289.4, 253.1, 222.8, 196.8, 174.9, 155.9, 139.4, 125.4, 113.1, 102.4, 92.2, 83.9, 77.4, 71.1, 65.3, 60.4, 56.0, 51.9, 48.4, 45.2, 42.3, 39.8, 37.6, 35.5, 33.7]\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(x, l)\n",
+    "ax.set_ylabel('L_m [nH]')\n",
+    "ax.set_xlabel('d [mm]')"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "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.11.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/coil_parasitics.py b/coil_parasitics.py
index 08a3603..0eddf6f 100644
--- a/coil_parasitics.py
+++ b/coil_parasitics.py
@@ -108,10 +108,15 @@ def elmer_solver(cwd, stdout_log=None, stderr_log=None):
     return result
 
 
-@click.command()
+@click.group()
+def cli():
+    pass
+
+
+@cli.command()
 @click.option('-d', '--sim-dir', type=click.Path(dir_okay=True, file_okay=False, path_type=Path))
 @click.argument('mesh_file', type=click.Path(dir_okay=False, path_type=Path))
-def run_capacitance_simulation(mesh_file, sim_dir):
+def self_capacitance(mesh_file, sim_dir):
     physical = dict(enumerate_mesh_bodies(mesh_file))
     if sim_dir is not None:
         sim_dir = Path(sim_dir)
@@ -168,10 +173,10 @@ def run_capacitance_simulation(mesh_file, sim_dir):
         
         capacitance_matrix = np.loadtxt(tmpdir / 'capacitance.txt')
 
-@click.command()
+@cli.command()
 @click.option('-d', '--sim-dir', type=click.Path(dir_okay=True, file_okay=False, path_type=Path))
 @click.argument('mesh_file', type=click.Path(dir_okay=False, path_type=Path))
-def run_inductance_simulation(mesh_file, sim_dir):
+def inductance(mesh_file, sim_dir):
     physical = dict(enumerate_mesh_bodies(mesh_file))
 
     if sim_dir is not None:
@@ -264,16 +269,149 @@ def run_inductance_simulation(mesh_file, sim_dir):
         elif P is None or R is None or U_mag is None:
             raise click.ClickException(f'Error during solver execution. Electrical parameters could not be calculated. See log files for details:\n{solver_stdout.absolute()}\n{solver_stderr.absolute()}')
 
-        U = math.sqrt(P*R)
+        V = math.sqrt(P*R)
         I = math.sqrt(P/R)
         L = 2*U_mag / (I**2)
 
-        assert math.isclose(U, 1.0, abs_tol=1e-3)
+        assert math.isclose(V, 1.0, abs_tol=1e-3)
 
+        print(f'Total magnetic field energy: {format_si(U_mag, "J")}')
+        print(f'Reference coil current: {format_si(I, "Ω")}')
         print(f'Coil resistance calculated by solver: {format_si(R, "Ω")}')
         print(f'Inductance calucated from field: {format_si(L, "H")}')
 
 
+@cli.command()
+@click.option('-r', '--reference-field', type=float, required=True)
+@click.option('-d', '--sim-dir', type=click.Path(dir_okay=True, file_okay=False, path_type=Path))
+@click.argument('mesh_file', type=click.Path(dir_okay=False, path_type=Path))
+def mutual_inductance(mesh_file, sim_dir, reference_field):
+    physical = dict(enumerate_mesh_bodies(mesh_file))
+
+    if sim_dir is not None:
+        sim_dir = Path(sim_dir)
+        sim_dir.mkdir(exist_ok=True)
+
+    sim = elmer.load_simulation('3D_steady', 'coil_mag_sim.yml')
+    mesh_dir = '.'
+    mesh_fn = 'mesh'
+    sim.header['Mesh DB'] = f'"{mesh_dir}" "{mesh_fn}"'
+    sim.constants.update({
+        'Permittivity of Vacuum': str(constants.epsilon_0),
+        'Gravity(4)': f'0 -1 0 {constants.g}',
+        'Boltzmann Constant': str(constants.Boltzmann),
+        'Unit Charge': str(constants.elementary_charge)})
+
+    air = elmer.load_material('air', sim, 'coil_mag_materials.yml')
+    ro4003c = elmer.load_material('ro4003c', sim, 'coil_mag_materials.yml')
+    copper = elmer.load_material('copper', sim, 'coil_mag_materials.yml')
+
+    solver_current = elmer.load_solver('Static_Current_Conduction', sim, 'coil_mag_solvers.yml')
+    solver_magdyn = elmer.load_solver('Magneto_Dynamics', sim, 'coil_mag_solvers.yml')
+    solver_magdyn_calc = elmer.load_solver('Magneto_Dynamics_Calculations', sim, 'coil_mag_solvers.yml')
+
+    copper_eqn = elmer.Equation(sim, 'copperEqn', [solver_current, solver_magdyn, solver_magdyn_calc])
+    air_eqn = elmer.Equation(sim, 'airEqn', [solver_magdyn, solver_magdyn_calc])
+
+    bdy_trace1 = elmer.Body(sim, 'trace1', [physical['trace1'][1]])
+    bdy_trace1.material = copper
+    bdy_trace1.equation = copper_eqn
+
+    bdy_trace2 = elmer.Body(sim, 'trace2', [physical['trace2'][1]])
+    bdy_trace2.material = copper
+    bdy_trace2.equation = copper_eqn
+
+    bdy_sub1 = elmer.Body(sim, 'substrate1', [physical['substrate1'][1]])
+    bdy_sub1.material = ro4003c
+    bdy_sub1.equation = air_eqn
+
+    bdy_sub2 = elmer.Body(sim, 'substrate2', [physical['substrate2'][1]])
+    bdy_sub2.material = ro4003c
+    bdy_sub2.equation = air_eqn
+
+
+    bdy_ab = elmer.Body(sim, 'airbox', [physical['airbox'][1]])
+    bdy_ab.material = air
+    bdy_ab.equation = air_eqn
+
+    bdy_if_top1 = elmer.Body(sim, 'interface_top1', [physical['interface_top1'][1]])
+    bdy_if_top1.material = copper
+    bdy_if_top1.equation = copper_eqn
+
+    bdy_if_bottom1 = elmer.Body(sim, 'interface_bottom1', [physical['interface_bottom1'][1]])
+    bdy_if_bottom1.material = copper
+    bdy_if_bottom1.equation = copper_eqn
+
+    bdy_if_top2 = elmer.Body(sim, 'interface_top2', [physical['interface_top2'][1]])
+    bdy_if_top2.material = copper
+    bdy_if_top2.equation = copper_eqn
+
+    bdy_if_bottom2 = elmer.Body(sim, 'interface_bottom2', [physical['interface_bottom2'][1]])
+    bdy_if_bottom2.material = copper
+    bdy_if_bottom2.equation = copper_eqn
+
+    potential_force = elmer.BodyForce(sim, 'electric_potential', {'Electric Potential': 'Equals "Potential"'})
+    bdy_trace1.body_force = potential_force
+    bdy_trace2.body_force = potential_force
+
+    # boundaries
+    boundary_airbox = elmer.Boundary(sim, 'FarField', [physical['airbox_surface'][1]])
+    boundary_airbox.data['Electric Infinity BC'] = 'True'
+
+    boundary_vplus1 = elmer.Boundary(sim, 'Vplus1', [physical['interface_top1'][1]])
+    boundary_vplus1.data['Potential'] = 1.0
+    boundary_vplus1.data['Save Scalars'] = True
+
+    boundary_vminus1 = elmer.Boundary(sim, 'Vminus1', [physical['interface_bottom1'][1]])
+    boundary_vminus1.data['Potential'] = 0.0
+
+    boundary_vplus2 = elmer.Boundary(sim, 'Vplus2', [physical['interface_top2'][1]])
+    boundary_vplus2.data['Potential'] = 1.0
+    boundary_vplus2.data['Save Scalars'] = True
+
+    boundary_vminus2 = elmer.Boundary(sim, 'Vminus2', [physical['interface_bottom2'][1]])
+    boundary_vminus2.data['Potential'] = 0.0
+
+    with tempfile.TemporaryDirectory() as tmpdir:
+        tmpdir = sim_dir if sim_dir else Path(tmpdir)
+
+        sim.write_startinfo(tmpdir)
+        sim.write_sif(tmpdir)
+        # Convert mesh from gmsh to elemer formats. Also scale it from 1 unit = 1 mm to 1 unit = 1 m (SI units)
+        elmer_grid(mesh_file.absolute(), 'mesh', cwd=tmpdir, scale=[1e-3, 1e-3, 1e-3],
+                   stdout_log=(tmpdir / 'ElmerGrid_stdout.log'),
+                   stderr_log=(tmpdir / 'ElmerGrid_stderr.log'))
+        solver_stdout, solver_stderr = (tmpdir / 'ElmerSolver_stdout.log'), (tmpdir / 'ElmerSolver_stderr.log')
+        res = elmer_solver(tmpdir,
+                   stdout_log=solver_stdout,
+                   stderr_log=solver_stderr)
+        
+        P, R, U_mag = None, None, None
+        solver_error = False
+        for l in res.stdout.splitlines():
+            if (m := re.fullmatch(r'StatCurrentSolve:\s*Total Heating Power\s*:\s*([0-9.+-Ee]+)\s*', l)):
+                P = float(m.group(1))
+            elif (m := re.fullmatch(r'StatCurrentSolve:\s*Effective Resistance\s*:\s*([0-9.+-Ee]+)\s*', l)):
+                R = float(m.group(1))
+            elif (m := re.fullmatch(r'MagnetoDynamicsCalcFields:\s*ElectroMagnetic Field Energy\s*:\s*([0-9.+-Ee]+)\s*', l)):
+                U_mag = float(m.group(1))
+            elif re.fullmatch(r'IterSolve: Linear iteration did not converge to tolerance', l):
+                solver_error = True
+
+        if solver_error:
+            raise click.ClickException(f'Error: One of the solvers did not converge. See log files for details:\n{solver_stdout.absolute()}\n{solver_stderr.absolute()}')
+        elif P is None or R is None or U_mag is None:
+            raise click.ClickException(f'Error during solver execution. Electrical parameters could not be calculated. See log files for details:\n{solver_stdout.absolute()}\n{solver_stderr.absolute()}')
+
+        V = math.sqrt(P*R)
+        I = math.sqrt(P/R)
+        Lm = (U_mag - 2*reference_field) / ((I/2)**2)
+
+        assert math.isclose(V, 1.0, abs_tol=1e-3)
+
+        print(f'Mutual inductance calucated from field: {format_si(Lm, "H")}')
+
+
 if __name__ == '__main__':
-    run_inductance_simulation()
+    cli()
 
diff --git a/twisted_coil_gen_twolayer.py b/twisted_coil_gen_twolayer.py
index 5f7e685..d57bae7 100644
--- a/twisted_coil_gen_twolayer.py
+++ b/twisted_coil_gen_twolayer.py
@@ -58,8 +58,6 @@ def traces_to_gmsh(traces, mesh_out, bbox, model_name='gerbonara_board', log=Tru
     occ = gmsh.model.occ
     eps = 1e-6
 
-    board_thickness -= 2*copper_thickness
-
     gmsh.initialize()
     gmsh.model.add('gerbonara_board')
     if log:
@@ -144,19 +142,32 @@ def traces_to_gmsh(traces, mesh_out, bbox, model_name='gerbonara_board', log=Tru
     print('Writing')
     gmsh.write(str(mesh_out))
 
-
-def traces_to_gmsh_mag(traces, mesh_out, bbox, model_name='gerbonara_board', log=True, copper_thickness=0.035, board_thickness=0.8, air_box_margin_h=30.0, air_box_margin_v=80.0):
+@contextmanager
+def model_delta():
     import gmsh
-    occ = gmsh.model.occ
-    eps = 1e-6
+    gmsh.model.occ.synchronize()
+    entities = {i: set() for i in range(4)}
+    for dim, tag in gmsh.model.getEntities():
+        entities[dim].add(tag)
 
-    board_thickness -= 2*copper_thickness
+    yield
 
-    gmsh.initialize()
-    gmsh.model.add('gerbonara_board')
-    if log:
-        gmsh.logger.start()
+    gmsh.model.occ.synchronize()
+    new_entities = {i: set() for i in range(4)}
+    for dim, tag in gmsh.model.getEntities():
+        new_entities[dim].add(tag)
+
+    for i, dimtype in enumerate(['points', 'lines', 'surfaces', 'volumes']):
+        delta = entities[i] - new_entities[i]
+        print(f'Removed {dimtype} [{len(delta)}]: {", ".join(map(str, delta))[:180]}')
 
+        delta = new_entities[i] - entities[i]
+        print(f'New {dimtype} [{len(delta)}]: {", ".join(map(str, delta))[:180]}')
+
+
+def _gmsh_coil_inductance_geometry(traces, mesh_out, bbox, copper_thickness, board_thickness, air_box_margin_h):
+    import gmsh
+    occ = gmsh.model.occ
     trace_tags = []
     trace_ends = set()
     render_cache = {}
@@ -218,35 +229,8 @@ def traces_to_gmsh_mag(traces, mesh_out, bbox, model_name='gerbonara_board', log
     (_dim, toplevel_tag), = tags
 
     (x1, y1), (x2, y2) = bbox
-    #print('first disk', first_disk)
-    #print('bbox', [occ.getBoundingBox(2, tag) for tag in first_disk[1]])
-    #print('last disk', last_disk)
-    #print('bbox', [occ.getBoundingBox(2, tag) for tag in last_disk[1]])
 
     first_geom = traces[0][0]
-    #contact_tag_top = occ.addCylinder(first_geom.start.x, first_geom.start.y, copper_thickness, 0, 0, copper_thickness, first_geom.width/2)
-    #contact_tag_bottom = occ.addCylinder(first_geom.start.x, first_geom.start.y, -board_thickness-copper_thickness, 0, 0, -copper_thickness, first_geom.width/2)
-
-    @contextmanager
-    def model_delta():
-        occ.synchronize()
-        entities = {i: set() for i in range(4)}
-        for dim, tag in gmsh.model.getEntities():
-            entities[dim].add(tag)
-
-        yield
-
-        occ.synchronize()
-        new_entities = {i: set() for i in range(4)}
-        for dim, tag in gmsh.model.getEntities():
-            new_entities[dim].add(tag)
-
-        for i, dimtype in enumerate(['points', 'lines', 'surfaces', 'volumes']):
-            delta = entities[i] - new_entities[i]
-            print(f'Removed {dimtype} [{len(delta)}]: {", ".join(map(str, delta))[:180]}')
-
-            delta = new_entities[i] - entities[i]
-            print(f'New {dimtype} [{len(delta)}]: {", ".join(map(str, delta))[:180]}')
 
     with model_delta():
         print('Fragmenting disks')
@@ -254,46 +238,35 @@ def traces_to_gmsh_mag(traces, mesh_out, bbox, model_name='gerbonara_board', log
         interface_tag_bottom = occ.addDisk(first_geom.start.x, first_geom.start.y, -board_thickness, first_geom.width/2, first_geom.width/2)
         occ.fragment([(3, toplevel_tag)], [(2, interface_tag_top), (2, interface_tag_bottom)], removeObject=True, removeTool=True)
 
-    #occ.synchronize()
-    #_, toplevel_adjacent = gmsh.model.getAdjacencies(3, toplevel_tag)
+    substrate = occ.addBox(x1, y1, -board_thickness, x2-x1, y2-y1, board_thickness)
 
-    #x0, y0, w = traces[0][0].start.x, traces[0][0].start.y, traces[0][0].width
-    #print(x0, y0, w)
-    #in_bbox = occ.getEntitiesInBoundingBox(x0-w/2-eps, y0-w/2-eps, -board_thickness-eps, x0+w/2+eps, y0+w/2+eps, eps, dim=1)
-    #print('in bbox', in_bbox)
-    #for dim, tag in in_bbox:
-    #    print(tag, 'adjacent', gmsh.model.getAdjacencies(dim, tag))
+    print('cut')
+    with model_delta():
+        print(occ.cut([(3, substrate)], [(3, toplevel_tag)], removeObject=True, removeTool=False))
 
-    #print('fragment', occ.fragment([(2, tag) for tag in toplevel_adjacent], [(2, interface_tag_top), (2, interface_tag_bottom)]))
+    return toplevel_tag, interface_tag_top, interface_tag_bottom, substrate
 
-    substrate = occ.addBox(x1, y1, -board_thickness, x2-x1, y2-y1, board_thickness)
 
+def traces_to_gmsh_mag(traces, mesh_out, bbox, model_name='gerbonara_board', log=True, copper_thickness=0.035, board_thickness=0.8, air_box_margin_h=30.0, air_box_margin_v=80.0):
+    import gmsh
+    occ = gmsh.model.occ
+    eps = 1e-6
+
+    gmsh.initialize()
+    gmsh.model.add('gerbonara_board')
+    if log:
+        gmsh.logger.start()
+
+    toplevel_tag, interface_tag_top, interface_tag_bottom, substrate = _gmsh_coil_inductance_geometry(traces, mesh_out, bbox, copper_thickness, board_thickness, air_box_margin_h)
+
+    (x1, y1), (x2, y2) = bbox
     x1, y1 = x1-air_box_margin_h, y1-air_box_margin_h
     x2, y2 = x2+air_box_margin_h, y2+air_box_margin_h
     w, d = x2-x1, y2-y1
-    z0 = -board_thickness-air_box_margin_v
-    ab_h = board_thickness + 2*air_box_margin_v
+    z0 = -2*copper_thickness-board_thickness-air_box_margin_v
+    ab_h = 2*copper_thickness + board_thickness + 2*air_box_margin_v
     airbox = occ.addBox(x1, y1, z0, w, d, ab_h)
 
-    occ.synchronize()
-
-    #trace_surface = gmsh.model.getBoundary([(3, toplevel_tag)], oriented=False)
-    #print('Fragmenting trace surface')
-    #with model_delta():
-    #    print(occ.fragment(trace_surface, [(2, interface_tag_top), (2, interface_tag_bottom)], removeObject=True, removeTool=False))
-
-    #print('Fragmenting trace')
-    #with model_delta():
-    #    print(occ.fragment([(3, toplevel_tag)], [(3, contact_tag_top), (3, contact_tag_bottom)], removeObject=True, removeTool=False))
-
-    print('cut')
-    with model_delta():
-        print(occ.cut([(3, substrate)], [(3, toplevel_tag)], removeObject=True, removeTool=False))
-
-    #print('Fragmenting substrate')
-    #with model_delta():
-    #    print(occ.fragment([(3, substrate)], [(3, toplevel_tag), (3, contact_tag_top), (3, contact_tag_bottom)], removeObject=True, removeTool=False))
-
     print('cut')
     with model_delta():
         print(occ.cut([(3, airbox)], [(3, toplevel_tag), (3, substrate)], removeObject=True, removeTool=False))
@@ -302,12 +275,10 @@ def traces_to_gmsh_mag(traces, mesh_out, bbox, model_name='gerbonara_board', log
     with model_delta():
         print(occ.fragment([(3, airbox)], [(3, toplevel_tag), (3, substrate)], removeObject=True, removeTool=False))
 
-    #occ.fragment([(3, substrate)], [(2, interface_tag_top), (2, interface_tag_bottom)])
-    #occ.fragment([(3, airbox)], [(3, substrate), (3, toplevel_tag)])
-
     print('Synchronizing')
     occ.synchronize()
 
+    first_geom = traces[0][0]
     pcx, pcy = first_geom.start.x, first_geom.start.y
     pcr = first_geom.width/2
     (_dim, plane_top), = gmsh.model.getEntitiesInBoundingBox(pcx-pcr-eps, pcy-pcr-eps, -eps, pcx+pcr+eps, pcy+pcr+eps, eps, 2)
@@ -341,7 +312,85 @@ def traces_to_gmsh_mag(traces, mesh_out, bbox, model_name='gerbonara_board', log
     gmsh.option.setNumber('Mesh.MeshSizeMin', 0.08)
     gmsh.option.setNumber('General.NumThreads', multiprocessing.cpu_count())
 
-    gmsh.write('/tmp/test.msh')
+    print('Meshing')
+    gmsh.model.mesh.generate(dim=3)
+    print('Writing to', str(mesh_out))
+    gmsh.write(str(mesh_out))
+
+
+def traces_to_gmsh_mag_mutual(traces, mesh_out, bbox, model_name='gerbonara_board', log=True, copper_thickness=0.035, board_thickness=0.8, air_box_margin_h=30.0, air_box_margin_v=80.0, mutual_offset=(0, 0, 5), mutual_rotation=(0, 0, 0)):
+    import gmsh
+    occ = gmsh.model.occ
+    eps = 1e-6
+
+    gmsh.initialize()
+    gmsh.model.add('gerbonara_board')
+    if log:
+        gmsh.logger.start()
+
+    m_dx, m_dy, m_dz = mutual_offset
+    m_dz += 2*copper_thickness + board_thickness
+
+    toplevel_tag1, interface_tag_top1, interface_tag_bottom1, substrate1 = _gmsh_coil_inductance_geometry(traces, mesh_out, bbox, copper_thickness, board_thickness, air_box_margin_h)
+
+    occ.translate([(3, toplevel_tag1), (2, interface_tag_top1), (2, interface_tag_bottom1), (3, substrate1)], m_dx, m_dy, m_dz)
+
+    toplevel_tag2, interface_tag_top2, interface_tag_bottom2, substrate2 = _gmsh_coil_inductance_geometry(traces, mesh_out, bbox, copper_thickness, board_thickness, air_box_margin_h)
+
+    (x1, y1), (x2, y2) = bbox
+    x1, y1 = x1-air_box_margin_h, y1-air_box_margin_h
+    x2, y2 = x2+air_box_margin_h, y2+air_box_margin_h
+    w, d = x2-x1, y2-y1
+    z0 = -2*copper_thickness-board_thickness-air_box_margin_v
+    ab_h = 4*copper_thickness + 2*board_thickness + 2*air_box_margin_v + m_dz
+    airbox = occ.addBox(x1, y1, z0, w, d, ab_h)
+
+    print('cut')
+    with model_delta():
+        print(occ.cut([(3, airbox)], [(3, toplevel_tag1), (3, toplevel_tag2), (3, substrate1), (3, substrate2)], removeObject=True, removeTool=False))
+
+    print(f'Fragmenting airbox ({airbox}) with {toplevel_tag1=} {substrate1=} {toplevel_tag2=} {substrate2=}')
+    with model_delta():
+        print(occ.fragment([(3, airbox)], [(3, toplevel_tag1), (3, toplevel_tag2), (3, substrate1), (3, substrate2)], removeObject=True, removeTool=False))
+
+    print('Synchronizing')
+    occ.synchronize()
+
+    first_geom = traces[0][0]
+    pcx, pcy = first_geom.start.x, first_geom.start.y
+    pcr = first_geom.width/2
+    (_dim, plane_top1), = gmsh.model.getEntitiesInBoundingBox(pcx-pcr-eps, pcy-pcr-eps, m_dz-eps, pcx+pcr+eps, pcy+pcr+eps, m_dz+eps, 2)
+    (_dim, plane_bottom1), = gmsh.model.getEntitiesInBoundingBox(pcx-pcr-eps, pcy-pcr-eps, m_dz-board_thickness-eps, pcx+pcr+eps, pcy+pcr+eps, m_dz-board_thickness+eps, 2)
+    (_dim, plane_top2), = gmsh.model.getEntitiesInBoundingBox(pcx-pcr-eps, pcy-pcr-eps, -eps, pcx+pcr+eps, pcy+pcr+eps, eps, 2)
+    (_dim, plane_bottom2), = gmsh.model.getEntitiesInBoundingBox(pcx-pcr-eps, pcy-pcr-eps, -board_thickness-eps, pcx+pcr+eps, pcy+pcr+eps, -board_thickness+eps, 2)
+
+    substrate1_physical = gmsh.model.add_physical_group(3, [substrate1], name='substrate1')
+    trace1_physical = gmsh.model.add_physical_group(3, [toplevel_tag1], name='trace1')
+    substrate2_physical = gmsh.model.add_physical_group(3, [substrate2], name='substrate2')
+    trace2_physical = gmsh.model.add_physical_group(3, [toplevel_tag2], name='trace2')
+    airbox_physical = gmsh.model.add_physical_group(3, [airbox], name='airbox')
+
+    interface_top1_physical = gmsh.model.add_physical_group(2, [plane_top1], name='interface_top1')
+    interface_bottom1_physical = gmsh.model.add_physical_group(2, [plane_bottom1], name='interface_bottom1')
+    interface_top2_physical = gmsh.model.add_physical_group(2, [plane_top2], name='interface_top2')
+    interface_bottom2_physical = gmsh.model.add_physical_group(2, [plane_bottom2], name='interface_bottom2')
+
+    airbox_adjacent = set(gmsh.model.getAdjacencies(3, airbox)[1])
+    in_bbox = {tag for _dim, tag in gmsh.model.getEntitiesInBoundingBox(x1+eps, y1+eps, z0+eps, x2-eps, y2-eps, z0+ab_h-eps, dim=2)}
+    airbox_physical_surface = gmsh.model.add_physical_group(2, list(airbox_adjacent - in_bbox), name='airbox_surface')
+    
+    points_airbox_adjacent = {tag for _dim, tag in gmsh.model.getBoundary([(3, airbox)], recursive=True, oriented=False)}
+    print(f'{points_airbox_adjacent=}')
+    points_inside = {tag for _dim, tag in gmsh.model.getEntitiesInBoundingBox(x1+eps, y1+eps, z0+eps, x1+w-eps, y1+d-eps, z0+ab_h-eps, dim=0)}
+    #gmsh.model.mesh.setSize([(0, tag) for tag in points_airbox_adjacent - points_inside], 300e-3)
+
+    gmsh.option.setNumber('Mesh.MeshSizeFromCurvature', 32)
+    gmsh.option.setNumber('Mesh.Smoothing', 10)
+    gmsh.option.setNumber('Mesh.Algorithm3D', 10)
+    gmsh.option.setNumber('Mesh.MeshSizeMax', 10)
+    gmsh.option.setNumber('Mesh.MeshSizeMin', 0.08)
+    gmsh.option.setNumber('General.NumThreads', multiprocessing.cpu_count())
+
     print('Meshing')
     gmsh.model.mesh.generate(dim=3)
     print('Writing to', str(mesh_out))
@@ -475,6 +524,10 @@ def print_valid_twists(ctx, param, value):
 @click.option('--arc-tolerance', type=float, default=0.02)
 @click.option('--mesh-out', type=click.Path(writable=True, dir_okay=False, path_type=Path))
 @click.option('--mag-mesh-out', type=click.Path(writable=True, dir_okay=False, path_type=Path))
+@click.option('--mag-mesh-mutual-out', type=click.Path(writable=True, dir_okay=False, path_type=Path))
+@click.option('--mutual-offset-x', type=float, default=0)
+@click.option('--mutual-offset-y', type=float, default=0)
+@click.option('--mutual-offset-z', type=float, default=5)
 @click.option('--magneticalc-out', type=click.Path(writable=True, dir_okay=False, path_type=Path))
 @click.option('--clipboard/--no-clipboard', help='Use clipboard integration (requires wl-clipboard)')
 @click.option('--counter-clockwise/--clockwise', help='Direction of generated spiral. Default: clockwise when wound from the inside.')
@@ -482,7 +535,7 @@ def print_valid_twists(ctx, param, value):
 def generate(outfile, turns, outer_diameter, inner_diameter, via_diameter, via_drill, via_offset, trace_width, clearance,
              footprint_name, layer_pair, twists, clipboard, counter_clockwise, keepout_zone, keepout_margin,
              arc_tolerance, pcb, mesh_out, magneticalc_out, circle_segments, mag_mesh_out, copper_thickness,
-             board_thickness):
+             board_thickness, mag_mesh_mutual_out, mutual_offset_x, mutual_offset_y, mutual_offset_z):
     if 'WAYLAND_DISPLAY' in os.environ:
         copy, paste, cliputil = ['wl-copy'], ['wl-paste'], 'xclip'
     else:
@@ -827,6 +880,10 @@ def generate(outfile, turns, outer_diameter, inner_diameter, via_diameter, via_d
         if mag_mesh_out:
             traces_to_gmsh_mag(traces, mag_mesh_out, ((-r, -r), (r, r)), copper_thickness=copper_thickness, board_thickness=board_thickness)
 
+        if mag_mesh_mutual_out:
+            m_dx, m_dy, m_dz = mutual_offset_x, mutual_offset_y, mutual_offset_z
+            traces_to_gmsh_mag_mutual(traces, mag_mesh_mutual_out, ((-r, -r), (r, r)), copper_thickness=copper_thickness, board_thickness=board_thickness, mutual_offset=(m_dx, m_dy, m_dz))
+
         if magneticalc_out:
             traces_to_magneticalc(traces, magneticalc_out)