From 4ea1b26293914aff6f146bded586ce2ee8aa68cb Mon Sep 17 00:00:00 2001 From: jaseg Date: Mon, 9 Oct 2023 19:56:42 +0200 Subject: Coupling sims work --- Untitled.ipynb | 72 +++++++++++++++ coil_parasitics.py | 152 ++++++++++++++++++++++++++++++-- twisted_coil_gen_twolayer.py | 205 +++++++++++++++++++++++++++---------------- 3 files changed, 348 insertions(+), 81 deletions(-) create mode 100644 Untitled.ipynb 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": [ + "
" + ] + }, + "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) -- cgit