{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from pystencils.session import *\n", "sp.init_printing()\n", "frac = sp.Rational" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Tutorial 06: Phase-field simulation of dentritic solidification\n", "\n", "This is the second tutorial on phase field methods with pystencils. Make sure to read the previous tutorial first. \n", "\n", "In this tutorial we again implement a model described in **Programming Phase-Field Modelling** by S. Bulent Biner.\n", "This time we implement the model from chapter 4.7 that describes dentritic growth. So get ready for some beautiful snowflake pictures.\n", "\n", "We start again by adding all required arrays fields. This time we explicitly store the change of the phase variable φ in time, since the dynamics is calculated using an Allen-Cahn formulation where a term $\\partial_t \\phi$ occurs." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "dh = ps.create_data_handling(domain_size=(300, 300), periodicity=True, \n", " default_target='cpu')\n", "φ_field = dh.add_array('phi', latex_name='φ')\n", "φ_field_tmp = dh.add_array('phi_temp', latex_name='φ_temp')\n", "φ_delta_field = dh.add_array('phidelta', latex_name='φ_D')\n", "t_field = dh.add_array('T')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This model has a lot of parameters that are created here in a symbolic fashion. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAA4CAYAAAAhDYavAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAXsklEQVR4Ae2djdXcNBOFSU4KSEIFkA4CVAB0QEgFQAfkUEHORweQCvjpIEkFCekAqIAkHeS7j1fyK2vlteS11rJ3dI7Xtqyf0dW1ZjySvbc+fPjwkQVDwBAwBAwBQ8AQaBuBW7du/U8S3tX2qbZ/tD2RDn+vvYUFEVgK51tmYC3YK1aUIWAIGAKGgCFQAQGn9H+Rzsaw+kjnf2j3qc4/q1Dd1Ra5JM63t4CiGvx8C3KajIaAIWAIGALTCKDEtP3C2O72eGUsnEbgm+jyU50/FH54sywsh8BiODdvYIk8NPar5bCzkgwBQ8AQMATWQkBjOtNceGJ+0Pa1ju9rezlHHpWFgXFN+iE0pvzUYBg3B0bLc4xAiGkSZ/Hux+Nsw5imDSw1gKeaL4Yi25khYAgYAobAhhFYxEOAcSUMfpKR9mLDWGSLrnY+iNpK+wmvDzv7XQKBApzfiINM046Gpg0sSf29NtygFgwBQ8AQMAT2g8CkhyCjqc+U5ruMdHtNgifQFrnX790kzs7Y/UdGFnZKMjS7yN09nbCA708da/fhVrIFFlkFAWEOqZ47ElWpYyuFCgs8qc+ExaOtyHwNchpHD728dX5KfjxaeALu6R7z0zEnKaw8vyjB30r/81hCpUHxce9+7NL0C8TH8mwlXm1jfL6r9v/Qgsx7xToHZ6X5W33wWYq7d1ronBEZHkvgJyPXLLoiAo5UHwn/q3C9T0HJjSNMWJD7l47tjZ0pwC5w3Th6A/IO+Jn0ENy0cHikvsf7hfF0b3jl5szxozeodI6h9QfxWx/X1IbOcFQ7WjGu6L/dYV2AM8Y+GBz1R5NThK5hCG3hwggIexbu4Tk04zbA3g3KnZEVRNvhCggYR49B3yo/1ZcopheSf9QTddzaTpH9qjxJb5fTHw91vfucAfldWsa0k2tmEnU1FaW2saCftVidMudcm1+LdXFZ94p1Ic5/CvjkNGFzBpYaxpMGrs/+5rg4a660QmGPq55Fo1lTYaTX9s712e5REy6/qpFv1V6UgoUVEIBzqjaLo6TVZvxcoZ9yqlTfzPXEkO/Up3vwMn8ey6D7943i7rp648vNn0tuDCnGnlc6htvcCxiNs3SlK2PyTTiVfyrsDutSnJ2twlos+mMQmluD5YTk1V3/dILBxQ2FcmOKhr2FhREQ7uD8r7ZHwnh0atClY3EpNzWuekiVvXZCaTcdcnHadCMbFT4He+Nn3n28dhern/DEfK2xpvOUu/O3OscIGg1Kx5jDmpdZY47ykxeP2dF0zmiljVyQ7O8kCuP0IKgts9YnqzzGbmYrSryHg7pPnWwV6zk4Kw8zbmCJ7dKH5tZgSUDcbWxdkODcUBhY9rbEAZJav7jOfxf+o8YVFes6hm/n4VLf0C9HVjvp9hpov9rNm0vdoty9trPRdk1y1PjZPj91/3hPzFMd+/EDgyfHc45h9t7180mauno+VyIMN69T0Cdsmwtqw+ias7UbsyesZ+KM4f5t3A+34wjOBVY1ArqOSFV7FKe0KHC/FuuZzv3NeJTWIuYj4HBl4LJ1VxkwusHapgozsFoqiXE0H8kN8POlWoORhcHst68kNw9vU4EpqZNTYuIKHx/9S+kwrn7XxveK+HK8X6t0Mr/SJ4PyN6EXk8KtFAmm2gxrcUxdwPTzwMN45MFSAt6y6BWtzpmj/Vjbf+zDazpPhok8n+r6fZVz0lNCwUrDdKBNCSZRXjSSef3RRaOL1nSBwsQvjEVctQyInsveOP9Cca/ErZ+VjgH3sTa4TTxvwkzyUukIlMtbSf3bM12s/dRCYDccvXZ+6h47xxPDPf12jGQOW4y2L1WPn27EcHuia37dFt4GHAnZuk1pm9GLyN5CKMEaeZW+81yqXwbTaHFbNoq15yT89LzrjBgMmW7TBYjJPKI/x4MEsfw5SuoPf57a6/pkHqXx85V93amyLO6Ae00cXH990L7v99z6gry8lNBMX0qujrPaM6Cy/ejl0zFPGLSXwZWn5pDb7/x5zl75GagxsJpp+x5lEcaMKcUcDfIZP3fCUfUp3pKkDlI8yo11Sv39Ht4Piud+hUc8gE3qKZ9XaXepF9Uu9HkSK9/2sb3ylWCNYYXOZ+PbZaPjpa5vEmvJDR4dt8L29Q3VRW9dhnEQsldCZHSFjA5Yuj6ZR2lQcsmbJBTOjseJuBQ2rr/4oGjf77nHyusV3ygfcstaKp1kgsffUJ4CXOTFiL5tOvc3Qv/g4NJ2bQnTTh2rrObaPyXzFq+7fizmaIv9I5mMn8H9WMpHx4XkQ42uoZy58ZPjkeIxvrqHKO0n9RSyKWxeL6oNeH8xTOMNDMAkjh81Yn1/KU821kEeDNtRA0vXNou1ZPcP7gN76bYu+PCTDn7zJ3LTkQFl5F1f/tJ7HQDUUcjNI8Apg6lCyrewEgLCH0LTB9yAewn/iF98/d/z108R+vbRZsLTw67/xW2dOz3oM7G+g4ChZaECAjvkqPHzPJ6gO8YCHpk3Tr8M0rjxgDGBhfV+bMjRbZvXi8KDF8T40vhgExaMjU/jeHc+9cJBFtaDTpg+2TLW91PNCw0sPsx2M3d4ULzkiQkNKZOFKd4bTDl5MOaShpriLWQgoIGChZs/auMDmKwHYuAoCby58179XmpYlNRx0bRqi+dex61E2zCkUHI+nZeP9H6Nho87uXdlcM+wjstChID4CD/h5nO3L+UnJe6KowHvjJ8RXzJPk/pH/PLcej1Sjn8IYk1viZ4yvRgBWoh1lPvk6Zax9vwbvEDRGVgCzBMu1XoIHQdfWBzvz3PyIMjJxW6+MNsfI4DCUmz3X1watFFCYI7btiTwWulujKuo4XDrTRTHKYrNv7LdXRaWxMHpLj4YQLrrEz88KPAmzal7aKKI/V0WHnhFmcr5QRt9wUPZyxkt3StHjZ8zyKAs77Wd0j9Me6UCY+TPgYFLmpN6auKePpk3JcBUfS7PlvRiCdYjkByid4C1dzoNeNEZWGoiymFgeQVo+Iw+Kj738eE+ThOfk3bqRgnLs+M0Ag+CaObNu6fiIG70UIR+qIsMVEVem9EC27twpJjVZtoL1+M24w5naoGv8XI9G0el9cZaSR5l231gCiEMT3VSZIjunKPGz5Ad+cevlPTzOLkznHhYDMfELpl4xIMnXut4uUCsl+Jz04sx0DqfiXWipEHU1rFGfmaDsGv6cLs/Ohg8wWlvcKGUwsB5yjNAGm+k5eTB0osJTRkWMhBQR+IZCAcMvg9T4o3yBkFJnliyuJ/j66ucOyMJ2WJDqmuzcIvbzIDt48DVG02T8istnH+vzbyxx2gx6PgARoQw7hAz/nsuR42f++Mnuufoe0OOQnip+G8+/kaGNKzzxdPPh0bD+7NET3neuiqKdNycPFvRi7lYewxy9lvGGv37Om7kHRdBpw4GIxGSL1ZDxJQRdFQQ5RTmoVxPdLJbmIkAg4my8kT8ZUERDDhY3MV94AYt+s8rwJeKgxOt/JURSpwnVm806bQLxKe+q4Z35Wu1gXUaDMilgbY/LM205/TCPvYkeHySY8cIFrM4avw8QnM3/OSeVv+iiMOHoq7Bbix7oOt4T1kwzRQWbwsPxjid5+o204sdssc/uVgf5xyN2TrW8JHlIsMgoIhAQR+9Pqm4H7X1r8TqGOL2n1fQMQoLL0H/WqyOT+ahPlcnymzwqry/ZvubzwpMYeHwZnqw6HsmSs/rucWvvk/Jc43XhSPrjdT0/H67trTCB2VnHF2BI3vjp9rDQ9BZugMuUo6/D3U80G3EK+xaL7r29brbY1Fr7zA+sjO2jrXDMfl5kP7PnmX1MwDyGidPB31QPETkS+4EOgPXYBd0jafSl9q+U/yfh9jBF3KJGuQJ0nCTYKzFXgafxPYFCKgvwPPozyZTRSgtAwcGFl9v7/szldbiphEQnt7zxf3zZjrHdaUQPhigyXFgDAnj6Bgy5fF746fa0+kd3WvnfBE+/JI7oCb5qbpML5ZTbpBDGOKIYTkLMx4cYyu8Uv8N/mR6q1hLbsZ//rj8+NMWivQeJQyp7gONPi53Pyef8gw+AJlbl6VLe0mEJ8TFip7sQ6VhgCJtkUfBsB/FfjE86b899Yva0hmfpdxRvsUwLa17b+mXxFJlNcFPycHC9cmx7ty+VB2mFy/kdd0q1pIbI5xPTHS2VLi/rQtdUCTWZLE3Q9YbN1zRU7vLg8fFwgwEhB8LON9pQwnF4X4ckTjHGCMM1iYcoux3BgIex3jd0Yyi9pNF/OSJ9YHGlm5c4XyEs6lGG0dTqMyL2yM/vxMUP82DIz+X6cV8rM5NuUWsNZ7xAPmnZE/aQL2B5cDhw4BkKAn3Vbi/gSfzqXymp3CnpRYbT+a3BIeXCYQDiwLD6VyUGef+6+I6HA3eCMvut9GS7AJPLb4fPK5Xj4ozpJgafKVj3uriQYxpglzOeSxz01895mMA7JGfrk18lf0SD+qmF8fItXz8ZrAW93gIfCQuhm/zDxC5E54pIW9oMBiylidrYFO6UkMJA25UoFAeOz6JAPO9P6iv/tOeNXJ09tEaOsWlAkYuwRsGhzP7PRcBj+u55ewhP2szwWPw8VuNF7mc81jmpt8DZrXb4DGtXc9FyheX+Ess1lLx/2/V1vJStuowvXiBXt0Y1tgxx+uuApwGBhbxkDa4vvihyh8sbFu8ggsVqBsu2wgtFUllx38ZcFSEcMQlmXRLHiU+jvAvLeAFs7AMAhgC3uuyTIlnlNIAP89agKymG0fP6P9E1qb4mZBvVlRtfeWFql3PXvSix+uc/VawlpyTS6qODKxzgLmWvFJevCLce+F07t+07LxJ4bUxTCby8IE8pl5rPZXt6kl2DOMV4pvAdQf8pOuawHIFDtWs0jCtia6VbQhECNzRYNx98COKb/pUhsetUwKqTQwkTFGUBD41MekRUtlMeYTGFVOeH3ujStdxJfP5iVHXoa6fzKO8uL75k1w+lpk1VVvSUKXtPC0q+31JPsmTzRWVvUoflcg41fapNkT58QZme7AkJ+uTWDcXB7h7X9cfxxd0Dh9GeUV65dsDP2lKMUdL+n6qb1VWlTGkREZAOBWm2hDlbYWf2WNIJL+dGgKbQ6D/DtbmJF9BYA2OD1XtYw1soYHFK5r8vUrvbXKD6D3FJQ0YXZ/MozQM8M9UxkmFOgcGle1fcT5pBM0p+1rzuD7F63jW1JjKYTE408/FU+nKuwt+wiHj6LJ3Ugv8XLZFVpoh0D4C8VuE7Uu8roS8FvybF0GDFkYQi8vjtUwYVinvBIojK48zzpgqpPylQyevk2Xpsq+5vJgHl8ZiL/wEN+Po8uxZm5/Lt8hKNAQaRsAMrLLOiRefe+MHgyoMDGRj00UleTDmkoZaWNmMYy/vmIwzirz6LGDpcV0LjL3wE/w8lsbRZdjUAj+XaYmVYghsBAEzsDI7asKTlHoyvDtRdE4e1l+F/wI/UWT2ZRbjE6ZkPKSy31wEaqyXy6p7Z/ykzcbRrJ4vSrQaP4uktMSGwE4QMAMrvyPxPI0NUPFTdnyeqiVOE5+Th6f4GkYQ5RJSdR6u2G8pAvRTymguLWdu+j3xEwyMo3OZkM63Nj/TUlmsIbBjBKp/pkFP1g+FH4vAJ78ZsQGc/aDvRfUGF4NXGDh/E0YExyV5UNg1jKCzDAG3dou34AgodsIT9fFYmw8pVvx1Hp6Qg/DyF8l89nffHB60LubHpVsc11/CNS9rSZ5a/ESWszjqG8N+C2NQzXuqIX6G3WLHFRGoyaeKYu+u6OoGlhDjjbVmFW9BjzLgDwwpFqKLyCiklBH0OlV2YR7K9QovVdzcOF+mN46yy3E37uDtRsVhbP2lPX+B1L9NmV1o5YSSi3ZiAPYGluJ4W+8P7TH+S/+NIJbY4+inteLrJedz8dsTP8FrNkcTYDc9BomDjCs176kW+JnoFouqgcAF+FRD7F2Webtmq9TRfIAzZXzUrLZW2Qz4fqAK6+C/sPpPKajNKG7+/LHzJuicNwGfawuNs5N5gsKpzyuaIPq8Q8nmDd7PZpSEMcUfrfZB5fHZCtqLImsx0Cffu77p5JPM3nPVf3LjDME9L87uK3jDNkOW3fCTtguDczjaw6c+38IYVPueaoGffZ/YQXUEavOpegP2UsHtWg3RwMZNjaI4W+nUkrGkXK/0IkMJRcD3iv5WPH9SCbHx4vQGl84xrD7X1r8NmJFHybuAAfTcHS+9o1+QqzR8qwx/JTLhebnr+j1xedUolDVcTBkuqbhSYb9wGbxRUJr/7PQ75CeYzOVoh6fjIv3b+hhU+55anZ9dh9jPpRCozadLtWPz9dScImTq5YkGuX5aZvNoac2O2oCh5L0fXZOcwZRsnq6hdO8JB7wofTiVp08kA0jpauGHQcQX5UsDyuqh2nNXsqWMk7ulBdZOLzlp6+ADoEF/0Kfnhs5DoHrWVuR74id9Mpejvj+3MgbVvqda4afvF9vXRaA2n+pKv6PSq3iwpLxQ3EsorqaglgLFW1Vs8DhlXuTdcHlqYth5oVRPkUEkDD7Txj8AxMYVi8YJaxsZBylO/Lo24238We04d/0VNWF0YwysGtSWPfETLGdxlIzq482MQeq32vdUE/ykXyzUR+ACfKrfiJ3UsLiB5ZQX3o3mFe3MPmQqsNTzw1+oZOPhMGSqcQnlP9ZMvwifwfesIHkpg6dkDJbY8Dqr7CUzS07Ww2FYobgxiJ6eWz5lqgyM1FpTuaUi7oWftHsWR9Un9MemxyC1YZF7qkF+lvLZ0i+AwFJ8WkCUqypicQNL6P0kJVv8P2pbQV1tQzG/dQNXltgzDCUMuCUWX4/KJ5n8uqTHo4nyL+BpY2F/VZnzxUmnlHz8WTJvEz5QCv4P8l838KQz5MX6qd/BtHFe1uVTqW274CfIqC1zObqHMWipe6opfi7PeCsxE4Gl+JRZnSUDgUX/7FnKipsZJdZPhykObwFx4cJv6rawMgLqG26679U3s//02ZWBMiyeOl25+UwjvXMyfCL5Z3neVAaeK/6cGaPNwsIIlHJU6Tc/Brk2L3JPGT8XJuQGi1uSTxts/qoiL+bBUifilv9CiqY3rlZtmVWeg0D3WQX13axpQuXjFfhFFEGOsHPTSM6xP81mCgrezmq/k4e8GKoW6iCQzdE9jEEV7injZx1ebqLUCnzaRLtbEfLOgoJwI6PIYmXD4mcf/7cMsN1OHy6I5UWKUl+8UH+xNgzvYtEibeXDU/BAZfSeK8WxHgmDK3u9GekvEJgOxGN1T7LN8lSlZFR5nYGpazXXyqWqvpq4Qo5uegxa+p4yfl7NbZJs6NJ8SlZikScRWHSKMFWTOplpmBcaKG2KMAXQynHqH/+2VbbxoTwYzY/Vp4M1V4rH4GAtVlMGluPga8k1+ONsxWN4YRRmt11p++DKNW73iNQ5EM7FHA0l2UI/ScbF76kttDvsJzteDoEafFpOuusp6VIG1pFyux6I22+pG4h/jQ2mlORKi0HCuqNwUTfTbPe1PVQZza1FkszdW5+Srfc0KQ4PHNNPs/4qJ8iPF68pg1Jt2l0Q3jyoZXE0brzL2+wYJPkWv6eMnzELrue8Bp+uB71lW1rNwFIn8zo8T2W47Qko5Oehkuti7Wd1BILBeNKTo7S8tEC/pgIvMzRnYCGo5IaHoRcVpfY/yVs0NUpZBIcDSrufIj1csd8aCJRw1NevPJsYgxyXFr2njJ+eBde3r8Gn60NxmRZXM7CWEc9KuRQCNiDnIy2s8H490/aJDKzF1nTlS3CdKY2jef1u/MzDyVIZArURMAOrNsIbKV+DMh4dvFNfymiwN0FH+k04MR36r7ZHc71fI0Vb9AQCxtEJgHTZ+DmNkaUwBC6FgBlYl0J6A/VocMYzw7RZk9N8LUAojFi39VYY2dTgCh1iHD0NuvHzND521RC4JAK3L1mZ1dU2AjIaWCfHpxu6bw+1Le3lpRMuLJbnL1jMuLo8/F2NxtFx4I2f49jYFUNgDQTMwFoD9Ybr9MaDBmv/jaeGpb2caMKDRfK8cTj41MPlJLCaPALGUY/Ezd74eYOFHRkCrSBgBlYrPdGQHFJgvG33wA3aDUm2jijCgXVXeK2+XEcCqzVGwDh6g4jx8wYLOzIEWkLg/26vftLJmWzcAAAAAElFTkSuQmCC\n", "text/latex": [ "$\\displaystyle \\frac{{{φ}_{(0,0)}}^{4}}{4} - {{φ}_{(0,0)}}^{3} \\left(\\frac{1}{2} - \\frac{m}{3}\\right) + {{φ}_{(0,0)}}^{2} \\left(\\frac{1}{4} - \\frac{m}{2}\\right) + \\frac{ε^{2} \\left({\\partial_{0} {{φ}_{(0,0)}}}^{2} + {\\partial_{1} {{φ}_{(0,0)}}}^{2}\\right)}{2}$" ], "text/plain": [ " 4 2 ⎛ 2 2⎞\n", "φ_C 3 ⎛1 m⎞ 2 ⎛1 m⎞ ε ⋅⎝D(φ[0,0]) + D(φ[0,0]) ⎠\n", "──── - φ_C ⋅⎜─ - ─⎟ + φ_C ⋅⎜─ - ─⎟ + ────────────────────────────\n", " 4 ⎝2 3⎠ ⎝4 2⎠ 2 " ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ε, m, δ, j, θzero, α, γ, Teq, κ, τ = sp.symbols(\"ε m δ j θ_0 α γ T_eq κ τ\")\n", "εb = sp.Symbol(\"\\\\bar{\\\\epsilon}\")\n", "\n", "φ = φ_field.center\n", "φ_tmp = φ_field_tmp.center\n", "T = t_field.center\n", "\n", "def f(φ, m):\n", " return φ**4 / 4 - (frac(1, 2) - m/3) * φ**3 + (frac(1,4)-m/2)*φ**2\n", "\n", "free_energy_density = ε**2 / 2 * (ps.fd.Diff(φ,0)**2 + ps.fd.Diff(φ,1)**2 ) + f(φ, m)\n", "free_energy_density" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The free energy is again composed of a bulk and interface part." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAboAAAEWCAYAAAAQKVIQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXxV9Z3/8dcnG1mAhIQkQBKWALLIEjAsgmvVCi7FpVZstepoUavTbdqpM/XXOm1nusx0s7UqWkdstUrdoIp7VVRQCPsaxBAgCySQkAXI/v39kYsTY4CEe5Nzc+/7+Xjcx71n/+Tkwjvne875HnPOISIiEqoivC5ARESkOynoREQkpCnoREQkpCnoREQkpCnoREQkpCnoREQkpCnoRE6RmTkzG+X7/JiZ/bSTy6Wb2XIzqzGzX3VvlSIS5XUBIl4xs0IgHWgGGoEVwO3Oub3dvOkFwAGgv9ONrCLdTkd0Eu4ud871BQYD+4Hf98A2hwFbjxdyZtYr/wDtrXVL6FPQiQDOuTrgGWD8sXFm9raZ3dpm+CYze+9k6zKzfmb2lpndZ2bWbtpjwI3Av5pZrZldaGb3mtkzZvYXM6sGbjKzCDO728w+NrODZrbYzJLbrGemma0ws0NmtsHMzjtBPUPM7FkzKzezXWb2jTbT7vWt+3FfU+oWM8vtwrLt644zs0VmVmlm28zsX82syDf/98zs2Xa1/d7MfnuyfSriDwWdCGBm8cC1wAd+ricFeBN43zn3jfZHbc65m4AngF865/o6597wTZpHa9Am+aZ/A7gCOBcYAlQC9/u2kQG8BPwUSAa+CzxrZqkd1BMB/B3YAGQAFwDfMrOL28z2BeAp37aXAn/owrLt6/4RMBzIBi4Crm8z71+AOWaW5Ft/FK37/M/H2Z0iAaGgk3D3gpkdAqpp/Y/5v/1Y1xDgHeBvzrl7urjsSufcC865FufcUeA24AfOuSLnXD1wL/BFXzhcDyxzzi3zzf86kAdc0sF6pwGpzrkfO+canHMFwMPA/DbzvOdbVzOtoTO5C8u2r/tLwH855yqdc0XAfcdmdM6VAsuBa3yj5gAHnHNrurivRLpEbeoS7q5wzr1hZpG0Hp28Y2bjnXP7TmFdlwK1wIOnsGz7C2CGAc+bWUubcc20XjwzDLjGzC5vMy0aeKuD9Q4DhvjC/JhI4N02w21/1iNArC9QO7Ns+7qHtBvXfvoi4A5aA/N6dDQnPUBHdCKAc67ZOfccrWFylm/0YSC+zWyDTrKah4FXgGVmltDVEtoN7wXmOueS2rxinXPFvml/bjctwTn38w7WuxfY1W7efs65jo7+TmXZ9nWXAplthrPaTX8BmGRmE4DLaG3uFOlWCjoRwFrNAwYA23yj1wNXmVm87365WzqxqruAfOBFM4vzo6QHgf80s2G++lJ99UHrua7LzexiM4s0s1gzO8/MMjtYzyqg2sy+77tQJNLMJpjZtE7UcCrLLgb+zcwG+M4l3tV2YpuLfp4EVjnn9nSiDhG/KOgk3P3dzGppPUf3n8CNzrktvmm/ARpove1gEZ04+vBdfLKA1qOhJWYWe4p1/Y7WC0NeM7MaWi+SmeHbxl5am1n/HSj3bet7dPDv2Xfe7XIgB9hF6/17jwCJnfhZTmXZHwNFvvnfoDXU6tvNswiYiJotpYeY7lcVke5iZncA851z57YZNxTYDgxyzlV7VpyEDR3RiUjAmNlgM5vtuw9wDPAvwPNtpkcA3wGeUshJT9FVlyISSDHAQ8AI4BCt9+f9EcB3gc5+YDettxaI9Ag1XYqISEhT06WIiIS0Xtl0OXDgQDd8+HCvyxARkSCyZs2aA865z3SF1yuDbvjw4eTl5XldhoiIBBEz293ReDVdiohISFPQiYhISFPQiYhISFPQiYhISFPQiYhISFPQiYhISFPQiYhISAvLoKuua+QP//iIdXsqvS5FRES6Wa+8Ydxf0RER3PfmTqrrmpgydIDX5YiISDcKyyO6uJhIcrKS+KDgoNeliIhINwvLoAOYmZ3M5uIqqusavS5FRES6UfgG3cgUWhys3lXhdSkiItKNwjbopg4dQExkhJovRURCXNgGXWx0JDlDk/igQEd0IiKhLGyDDuDM7BS2lFRRdVTn6UREQlVYB93MbJ2nExEJdQEJOjObY2b5ZrbTzO7uYPpXzGyj77XCzCZ3dtnuNGVoEjFROk8nIhLK/A46M4sE7gfmAuOB68xsfLvZdgHnOucmAT8BFnZh2W4TGx3J1KFJrFTQiYiErEAc0U0HdjrnCpxzDcBTwLy2MzjnVjjnjvW39QGQ2dllu9vM7BS2llZTdUTn6UREQlEggi4D2NtmuMg37nhuAV7u6rJmtsDM8swsr7y83I9yP21mdgrOwapCnacTEQlFgQg662Cc63BGs/NpDbrvd3VZ59xC51yucy43NTX1lArtSE5WEn2iIlj5sZovRURCUSA6dS4CstoMZwIl7Wcys0nAI8Bc59zBrizbnVrP0w3QBSkiIiEqEEd0q4HRZjbCzGKA+cDStjOY2VDgOeAG59yOrizbE2Zmp7BtXzWHjjT09KZFRKSb+R10zrkm4C7gVWAbsNg5t8XMbjez232z/RBIAf5oZuvNLO9Ey/pbU1edOdJ3nk7304mIhJyAPI/OObcMWNZu3INtPt8K3NrZZXva5KzE1vN0BQf5/OmDvCxFREQCLKx7RjmmT1QkZwwboH4vRURCkILOZ2Z2Ctt1nk5EJOQo6HyOnafTUZ2ISGhR0PlMykwkNlr9XoqIhBoFnc//nadT0ImIhBIFXRtnZqewfV8NFYd1nk5EJFQo6NqYmZ0CwKpdOqoTEQkVCro2JmUmERcdqQtSRERCiIKujZioCHKHD1AHzyIiIURB187M7BTy99dwsLbe61JERCQAFHTtzMxOBtTvpYhIqFDQtXPsPN1K3WYgIhISFHTtREdGMCM7mXd2lONch8+AFRGRXkRB14ELxqWz++ARPi6v9boUERHxk4KuAxeMTQPgjW1lHlciIiL+UtB1YEhSHOMH9+eNrfu9LkVERPwUkKAzszlmlm9mO83s7g6mjzWzlWZWb2bfbTet0Mw2tX3yeDC4cFwaa/dUqjswEZFezu+gM7NI4H5gLjAeuM7MxrebrQL4BvA/x1nN+c65HOdcrr/1BMqF49NpcfDWdjVfioj0ZoE4opsO7HTOFTjnGoCngHltZ3DOlTnnVgONAdhej5gwJJG0fn14c7uaL0VEerNABF0GsLfNcJFvXGc54DUzW2NmCwJQT0BERBgXjEtj+Y4DNDS1eF2OiIicokAEnXUwris3oM12zk2ltenzTjM7p8ONmC0wszwzyysvLz+VOrvsgrHp1NY38aGeZiAi0msFIuiKgKw2w5lASWcXds6V+N7LgOdpbQrtaL6Fzrlc51xuamqqH+V23uxRA+kTFaGrL0VEerFABN1qYLSZjTCzGGA+sLQzC5pZgpn1O/YZ+DywOQA1BURcTCRnjx7IG9vK1EuKiEgv5XfQOeeagLuAV4FtwGLn3BYzu93Mbgcws0FmVgR8B7jHzIrMrD+QDrxnZhuAVcBLzrlX/K0pkC4Yl07xoaPk76/xuhQRkZDzxIe7Wbenslu3ERWIlTjnlgHL2o17sM3nfbQ2abZXDUwORA3d5VgvKW9uK2PsoP4eVyMiEjrqm5q5d+kWbj07mylDB3TbdtQzykmk9Y9lUmYib2zTeToRkUDK31dDY7NjYkZit25HQdcJF45LZ/3eQ5TX6GGsIiKBsrm4GkBBFwwuGJeGUy8pIiIBtam4isS4aDIHxHXrdhR0nTB+cH+GJMaq+VJEJIC2lFQxIaM/Zh3djh04CrpOMDM+Ny6Ndz86QF1js9fliIj0eg1NLWwvrWFCNzdbgoKu0y4cl87RxmZWFqiXFBERf+3YX0NDcwsThijogsbM7BTiYyLVS4qISABsKakCuv9CFFDQdVpsdGsvKf/Yrl5SRET8tam4in6xUQxLie/2bSnouuCCcemUVtWxpaTa61JERHq1TcXVnD6k+y9EAQVdl3xubBpmrb2kiIjIqWlsbmFbaXWPNFuCgq5LBvbtw5SsJN1mICLih51ltTQ0tfTIFZegoOuyi08fxKbiKgoPHPa6FBGRXmlTceuFKAq6IDUvJwMzeG5dsdeliIj0SluKq0iIiWRESkKPbE9B10WDEmM5a9RAnltbREuLrr4UEemqTcVVnD4kkYiI7r8QBRR0p+SqqRkUVR4lb3f3PkNJRCTUNDW3sLW0useaLUFBd0ouPn0Q8TGRPLe2yOtSRER6lYIDh6lrbGFiZs8931NBdwriY6KYO2EwL20sVd+XIiJdsKnIdyFKD3T9dUxAgs7M5phZvpntNLO7O5g+1sxWmlm9mX23K8sGq6unZlBT38Rr6hJMRKTTNhVXER8TSXZq3x7bpt9BZ2aRwP3AXGA8cJ2ZjW83WwXwDeB/TmHZoDQzO4UhibFqvhQR6YItJVWMH9yfyB66EAUCc0Q3HdjpnCtwzjUATwHz2s7gnCtzzq0GGru6bLCKiDCunJrB8h3llNXUeV2OiEjQa25xbCnp2QtRIDBBlwHsbTNc5BsX0GXNbIGZ5ZlZXnl5+SkVGmhXTsmkxcHS9SVelyIiEvR2HajlSENzrwy6jo4/O3uDWaeXdc4tdM7lOudyU1NTO11cdxqV1pfJWUk8u1Y3j4uInMzm4tYO8Xuqj8tjAhF0RUBWm+FMoLOHOP4sGxSunprBttJqtuqJBiIiJ7SpuIrY6AhGpvZMjyjHBCLoVgOjzWyEmcUA84GlPbBsULhs0hCiI00XpYiInMSm4irGDe5PVGTP3tnm99acc03AXcCrwDZgsXNui5ndbma3A5jZIDMrAr4D3GNmRWbW/3jL+ltTT0pOiOH8MWm8sL6EpuYWr8sREQlKLS2OrSXVPXr/3DFRgViJc24ZsKzduAfbfN5Ha7Nkp5btba6amslrW/fz7s4DnD8mzetyRESCTuHBw9TWN/X4+TlQzygB8bmxaSTFR/OcLkoREelQTz+apy0FXQDEREXwhclDeG3LPqrr2t8qKCIiW0qqiYmKYHR6z/WIcoyCLkCumppJfVMLyzaWel2KiEjQ2VRUxbhB/Yju4QtRQEEXMJMzE8lOTVDzpYhIO845NpdUedJsCQq6gDEzrp6ayarCCgoPHPa6HBGRoLGn4gg1dU0KulDwxTMyiY40Hn1/l9eliIgEjWMXonhxxSUo6AIqvX8sV+RksDhvLxWHG7wuR0QkKGwuriY60jgtvZ8n21fQBdiCc7Kpa2zh8ZWFXpciIhIUNhdXMWZQP2KivIkcBV2AjU7vx+fGpvH4yt0cbdDTx0UkvDnn2FRc5VmzJSjousWCc7KpONzAM+r/UkTCXFHlUaqONnp2IQoo6LrFjBHJTM5M5JF3C2hu6ewTi0REQs/aPZUATM5M8qwGBV03MDMWnDOS3QeP8NqWfV6XIyLimdWFFfTtE8W4wf09q0FB103mTBjE0OR4HlpegHM6qhOR8LR6VyVThw0gMqKj52z3DAVdN4mMMG49ewTr9x5idWGl1+WIiPS4Q0cayN9fw/ThAzytQ0HXja45I4sB8dEsXP6x16WIiPS4Nbtb/8jPHZ7saR0Kum4UFxPJDWcO541tZewsq/G6HBGRHrWqsILoSCMny7sLUSBAQWdmc8ws38x2mtndHUw3M7vPN32jmU1tM63QzDaZ2XozywtEPcHkxjOH0ScqgoeXq1swEQkveYWVTMxIJDY60tM6/A46M4sE7gfmAuOB68xsfLvZ5gKjfa8FwAPtpp/vnMtxzuX6W0+wSenbhy+ekcnz64opq67zuhwRkR5R19jMxqJDTBvhbbMlBOaIbjqw0zlX4JxrAJ4C5rWbZx7wuGv1AZBkZoMDsO1e4dazs2lsaeGxFYVelyIi0iM27D1EY7Nj2rDQCLoMYG+b4SLfuM7O44DXzGyNmS0IQD1BZ8TABC4eP4i/fLCb2vomr8sREel2qwsrAMj1+IpLCEzQdXRzRPsbx040z2zn3FRamzfvNLNzOtyI2QIzyzOzvPLy8lOv1iO3nZtNdV2TOnsWkbCwurCSMen9SIqP8bqUgARdEZDVZjgTKOnsPM65Y+9lwPO0NoV+hnNuoXMu1zmXm5qaGoCye9aUoQO4cFwa9/9jp87ViUhIa25xrN1dGRRHcxCYoFsNjDazEWYWA8wHlrabZynwVd/VlzOBKudcqZklmFk/ADNLAD4PbA5ATUHpnkvH09js+MUr+V6XIiLSbbaVVlNT38T0ILgQBQIQdM65JuAu4FVgG7DYObfFzG43s9t9sy0DCoCdwMPA133j04H3zGwDsAp4yTn3ir81BavhAxP4p7NG8OzaItbtUW8pIhKa8j45PxccQRcViJU455bRGmZtxz3Y5rMD7uxguQJgciBq6C3u+twonltbxL1Lt/D812cT4WH/byIi3WF1YSUZSXFkJMV5XQqgnlF6XN8+UXx/zlg2FFXxrJ5XJyIhxjnH6sKKoDk/Bwo6T1w5JYOcrCR+8Uo+NXWNXpcjIhIweyqOUFZTz7QgabYEBZ0nIiKM//jC6RyorecP/9jpdTkiIgFz7GktwXIhCijoPDM5K4lrzsjk0fd3UVBe63U5IiIBsXpXBYlx0YxK7et1KZ9Q0Hnoe3PG0Ccqkp+8uNXrUkREAmL17gqmDR8QVBfaKeg8lNYvlm9cMIq38st5a3uZ1+WIiPjlQG09BeWHg+a2gmMUdB67adYIsgcm8JMXt9LQ1OJ1OSIipyzPd34umC5EAQWd52KiIvh/l42n4MBh/vSenlknIr3X6sIK+kRFMDEj0etSPkVBFwTOH5vG58en8+vX81m/95DX5YiInJLVhRXkZCURExVc0RJc1YSxX35xEun9Y7nzibVUHm7wuhwRkS45XN/ElpLqoLqt4BgFXZBIio/hj1+ZSnlNPd9evJ6WlvZPOhIRCV7r9hyiucUF3YUooKALKpMyk/jh5eN5O7+c+9/SjeQi0nusLqwgwmDq0CSvS/kMBV2Q+cqMoVyRM4Rfv7GD9z464HU5IiKdsrqwgnGD+9MvNtrrUj5DQRdkzIz/umoio9P68s2n1rGvSg9pFZHg1tjcwro9h4LutoJjFHRBKD4mij9+5QyONjZz55NraWzW/XUiEry2lFRztLFZQSddMyqtL7+4ehJrdlfyi5e3e12OiMhxfVhwEIBpQfRonrYUdEHs8slDuGnWcB55bxcvbyr1uhwRkQ69s6OcsYP6kdY/1utSOhSQoDOzOWaWb2Y7zezuDqabmd3nm77RzKZ2dtlw9++XjCMnK4nvLN6gi1NEJOjU1DWyurCC88akeV3KcfkddGYWCdwPzAXGA9eZ2fh2s80FRvteC4AHurBsWIuJimDhV89gaHI8//TYal7fut/rkkREPvH+zoM0NjvOG5PqdSnHFRWAdUwHdjrnCgDM7ClgHtD22TPzgMedcw74wMySzGwwMLwTy4a9tH6xPH3bTG58dBW3/2UNv/7SZOblZHhdVkhraXHsq66jqPIo1Ucbqa5r9L03UX20kZq6Jmrrm4iIMPpERRAbHUGfqEj6RPneoyNISYhhaHI8w1ISSOvXJ6geWyISKO/sKKNfnyjOGBac5+cgMEGXAextM1wEzOjEPBmdXBYAM1tA69EgQ4cO9a/iXigpPoYnvjaTWx5bzbeeXs+Rhmaumx5++yHQDtbWk7+/hsIDRyg8eJjCA4cpPHiY3QePUH+cp0nEx0TSLzaKvn2icA7qm1qoa2ymvqmF+qZmGps/26tNTFQEWQPiGJaSwNDkeEam9WX68GRGp/VVAEqv5Zzjre3lnH3aQKIjg/eSj0AEXUf/Stv/Sz/ePJ1ZtnWkcwuBhQC5ublh2T9W3z5RPHbzdO54Yg3/9twmDtc3cevZ2V6X1Wscrm9iU3EVG4sOsWFvFRuKDlFUefST6TFREQxLjmf4wATOPS2V4QMTyBoQz4D4GPrFRtE/Lpp+sVEn/Qfd3OKob2qmrLqePRVH/u91sPV91a4KauubAEiKj2ba8GRmjEhm+ohkxg/uT1QQ/4ch0tb2fTXsq67jvNOC9/wcBCboioCsNsOZQEkn54npxLLSRlxMJAtvyOVbT6/jpy9to7a+iW9eMBozHRW0V1Zdx8qCg6z8+CBr91Sys6yWY12IZg6IY3JmEl89cxjjBycyIjWBwf1jA3J0FRlhxMdEMXxgFMMHJnxmunOOvRVH+XDXQVbtqmBVYcUn514TYiKZNiKZSycO5uIJg+gfhL1MiBzzVn7rA6PPDeLzcxCYoFsNjDazEUAxMB/4crt5lgJ3+c7BzQCqnHOlZlbeiWWlnZioCO6bP4X4mE389o2PqK1r4t8uGUdkmDeBVR5u4IOCg6z4+CArPj7Ax+WHAegf23r+4JKJg5mcmcSkzERS+vbxrE4zY2hKPENT4rkmt/XvvP3Vda2ht6uCd3aU871nNvKDFzZz4bg0vjA5g/PGpBIbHelZzSIdeTu/nNOH9Cc9SG8rOMbvoHPONZnZXcCrQCTwqHNui5nd7pv+ILAMuATYCRwBbj7Rsv7WFA6iIiP45dWT6Nsnikfe20Xe7kp++cVJnJbez+vSekxjcwtrdlfyzo5y3skvZ2tpNdB6Dm36iGS+lJvFrJEDGT+kf9D/EZDeP5bLJw/h8slDcM6xfu8hlqwv4cWNJSzbtI9+sVHMnTCIK3IymJmdovN64rmqo42s2V3JHeeO9LqUk7LWCyF7l9zcXJeXl+d1GUHBOceS9SX8x9+3UFvfxB3njeLO80fSJyo0//ovOXSUd3aU83Z+Ge/vPEhtfRNREcbUYQM4e9RAZo1KYVJmUlCfGO+KpuYWVnx8kCXrS3h1yz5q65sYldaXBedkc0VORtA94FLCx0sbS7nzybU8c/uZQfNoHjNb45zL/cx4BV1oOFhbz09e3MoL60t83YdN5IxhwfHl80dDUwt5hRW87Qu3HftrARiSGMu5Y9I497RUZo9KCcoe0wOtrrGZVzbv46HlBWwrrSa9fx9uOWsE100fGhY/vwSX7/5tA69v3c+aey4MmguoFHRh4q38Mn7w3CZKq+v46sxhfG/OWPr2CcSp2J5TfOgob+eX8XZ+OSt2HuBwQzPRkcb0Ecmcd1oa545JZXRa37C9AMc5x7sfHeDBdz5mxccH6RcbxfUzh3Hz7OGk9QvucyUSGlpaHDN+9iYzs1P4/XVTvC7nEwq6MFJb38T/vJrPopWFDO4fyx3njeQLORkkxgXnX/3HuhB6f+dB3v2o/JOjtoykOM4bk8p5Y9KYNTKFhF4W2D1hY9EhHlpewMubSomKiODLM4byzQtGMyAhxuvSJIRtLq7ist+/x6+umczVZ2R6Xc4nFHRhaM3uSu5duoVNxVX0iYpg7oRBXDttKDOzkz09GqprbGbtnkpW7Gy9OnJDURXNLY6YqAimDR/Aeaelcd6YVEaF8VFbVxUeOMxDyz/m6dV76RcbzTcvGM0NZw4LmXOVElx+/+ZH/Or1HeTdcyEDPbyCuT0FXZhyzrG5uJqn8/awZH0JNXVNDEuJ50u5WVw9NZNBid3b1HXsnrGNxYfYWFTFhr2HWL/3EPVNLUQYTMpMYvaoFGaNHMgZwwboEno/5e+r4acvbeXdjw6QnZrAPZeO4/wxafqDQQLq6gdW0NTcwpK7zvK6lE9R0Al1jc28vLmUp1fv5YOCik+CZtzgfoxJ78eYQf0ZO6jfKTV7Oec4UNtAadVRiiqPsqWkio1FVWwqruLQkUYAYiIjGDe4H1OHDWD2yIFMz07WDdHdwDnHW/ll/PTFbRQcOMzZowfy/y4bH1a3nkj3qTzcwBk/fZ27Pjea71x0mtflfIqCTj5l98HDPLOmiLzCSrbvq6bSF0YA6f37MHZQf7KS44iKiCDCjMgIiIgwIs2IjDCafZ0elx6qo6TqKKVVdTS06RsyMsIYk96PSZmJTMxMZHJmEqel99Pl8D2osbmFP6/czW/f2EFtfRM3zBzGv84Zq3Od4pcl64v55lPref7rs5gyNLg6clbQyXE55yivqWfbvhry91WzvbSG7ftqKK06SnOLo8W19t/Y7BwtvvcIM9L79WFIUhyDk+IYkhTLkMQ4BifGMiQpjlFpfdUMGSQqDzfwmzd28OcPdpORFMcvvziJWSMHel2W9FLffno97+woZ/UPLgy6jhiOF3T6004wM9L6x5LWP5ZzT+tcn3XOOZ336SUGJMTw43kTuHzyEL73tw18+eEPufHMYXx/7ljiY/RfgHReS4vjnR3lnDN6YNCF3ImoHUlOiUKu95k2PJmXv3kON88ezuMf7GbOb9/lw4KDXpclvcjG4ioqDjdw/tjgflpBewo6kTASFxPJjy4/nae+NhOAaxd+wL1Lt3CkocnjyqQ3eGt7GWZwzujgflpBewo6kTA0IzuFV751NjfNGs5jKwq55Hfvsrm4yuuyJMi9nV9GTlZSr+uQQEEnEqbiY6K49wun89evzaS+qYWrHljBEx/upjdeoCbd70BtPRuLqzh/TO9qtgQFnUjYO3NkCi/+81nMzE7hB89v5ttPr+dwvZoy5dOW7yjHORR0ItI7pfTtw2M3TeM7F53Gkg0lzLv/fT7aX+N1WRJE3txexsC+MZw+pL/XpXSZgk5EgNYOAb5xwWj+cssMDh1p4At/eJ/n1xV5XZYEgdr6Jt7ctp+LTx/UKx/661fQmVmymb1uZh/53ju8Td7M5phZvpntNLO724y/18yKzWy973WJP/WIiP9mjxrIS984m4kZiXz76Q3823ObqGts9ros8dBrW/ZR19jClVMyvC7llPh7RHc38KZzbjTwpm/4U8wsErgfmAuMB64zs/FtZvmNcy7H91rmZz0iEgDp/WN58mszuP3ckfx11R6+/PAHlNfUe12WeOT5dcVkDojjjGHB1eVXZ/kbdPOARb7Pi4ArOphnOrDTOVfgnGsAnvItJyJBLCoygrvnjuWBr0xla2k1V9z/Ptv3VXtdlvSwspo63t95gCtyMnptRxH+Bl26c64UwPfe0eU4GcDeNsNFvnHH3GVmG83s0eM1fYqId+ZOHMzfboC8yXgAABKESURBVJtFU0sLV/9xBW9u2+91SdKD/r6hlBYHV0wZ4nUpp+ykQWdmb5jZ5g5enT0q6+hPgGM36jwAjARygFLgVyeoY4GZ5ZlZXnl5eSc3LSKBMDEzkSV3nkV2al9ufTyPh5cX6H67MPHCumImZiQyKq33PubppEHnnLvQOTehg9cSYL+ZDQbwvZd1sIoiIKvNcCZQ4lv3fudcs3OuBXiY1mbO49Wx0DmX65zLTU3tXd3PiISCQYmxLL7tTOacPoj/XLaNu5/d9KlHM0no2VlWy6biKubl9N6jOfC/6XIpcKPv843Akg7mWQ2MNrMRZhYDzPctdywcj7kS2OxnPSLSjeJiIrn/y1O56/xRPJ23lxv+9CGVhxu8Lku6yZL1xUQYfGFyeAfdz4GLzOwj4CLfMGY2xMyWATjnmoC7gFeBbcBi59wW3/K/NLNNZrYROB/4tp/1iEg3i4gwvnvxGH57bQ7r9h7i6gdWsLfiiNdlSYA553h+XTGzRw0krX+s1+X4xa+HUTnnDgIXdDC+BLikzfAy4DO3DjjnbvBn+yLinSumZJAxII5bHlvNVQ+sYNHN0xnfC3vNkI6t2V1JUeVRvn3haV6X4jf1jCIip2za8GSeuWMWURHGtQ+tZOXHer5dqHh+XTGx0RFcPGGQ16X4TUEnIn45Lb0fz94xi/TEWG58dBUvbSz1uiTxU0NTCy9tKuXz4wfRt0/vfwq9gk5E/DYkKY5nbj+TiZmJ3PXXtSxaUeh1SeKHd3aUc+hIY6/t8qs9BZ2IBERSfAxP3DqDC8am86OlW/jvV7frXrte6oV1xSQnxHDW6IFelxIQCjoRCZjY6EgevH4q103P4v63PuZfn9lIU7PutetNqusaeWPbfi6fNJjoyNCIiN7f+CoiQSUqMoL/unIiqX37cN8/dlJb38Rv5+fQJyrS69KkE17ZvI/6phauCJFmS9ARnYh0AzPjO58fwz2XjuPlzfv42uNrONqgR/30Bi+sK2ZYSjw5WUlelxIwCjoR6Ta3np3Nz6+ayLsflXPjo6uormv0uiQ5gX1VdawsONirn1TQEQWdiHSr+dOHct/8KazdU8lXHv6QCnUZFrSWbijGOUKq2RIUdCLSAy6fPISFXz2D/P01XPvQSvZX13ldkrTjnOO5tcVMzkpixMAEr8sJKAWdiPSIz41N57Gbp1Fy6CjXPLhS/WMGmQ8KKti+r4Yv5WZ6XUrAKehEpMfMGjmQJ742k6qjjVzz4Ep2ltV6XZL4LFz+MSkJMVw9VUEnIuKXnKwknr5tJk0tLcxfuJLt+6q9Lins5e+r4a38cm6aNZzY6NC7DURBJyI9buyg/jy14EwiI4z5Cz9gc3GV1yWFtYXLC4iLjuT6mcO8LqVbKOhExBOj0vqy+LYzSYiJ4rqHP2DN7kqvSwpLpVVHWbK+mGunZTEgIcbrcrqFgk5EPDMsJYHFt59JSkIMX/3Th3xQoMf89LT/fb8QB9xy1givS+k2CjoR8VRGUhxP33Ymg5PiuOl/V7F8R7nXJYWN6rpGnvxwD5dOHExWcrzX5XQbv4LOzJLN7HUz+8j3PuA48z1qZmVmtvlUlheR0JbeP5anFsxkeEoCty7K481t+70uKSw8+eEeauubWHBOtteldCt/j+juBt50zo0G3vQNd+QxYI4fy4tIiBvYtw9PLZjJ2MH9uO3Pa1i2SQ9w7U71Tc08+t4uzho1kAkZiV6X0638Dbp5wCLf50XAFR3N5JxbDlSc6vIiEh6S4mP4y60zyMlK4q4n1/L8uiKvSwpZS9aXUFZTH/JHc+B/0KU750oBfO9p3bW8mS0wszwzyysvVxu+SKjqHxvNon+azszsFL6zeANPfrjH65JCTkuLY+HyAsYN7s/ZIfJw1RM5adCZ2RtmtrmD17yeKPAY59xC51yucy43NTW1JzctIj0soU8Uj940jfNOS+Xfn9/Eo+/t8rqkkPJWfhk7y2q57ZzskHpKwfGc9MGrzrkLjzfNzPab2WDnXKmZDQbKurh9f5cXkRAVGx3JQzfk8o2/ruPHL27laGMzd54/yuuyQsJD7xSQkRTHpZMGe11Kj/C36XIpcKPv843Akh5eXkRCWExUBH/48hTm5Qzhv1/N51ev5eOc87qsXm3tnkpWFVbwT2eNIDoyPO4w8/en/DlwkZl9BFzkG8bMhpjZsmMzmdlfgZXAGDMrMrNbTrS8iMgxUZER/PpLOVybm8Xv/7GT/3xpm8LODwvfKSAxLpr507K8LqXHnLTp8kSccweBCzoYXwJc0mb4uq4sLyLSVmSE8bOrJhIXE8kj7+3iSGMzP5k3gciI0D+/FEg79tfw6tZ9fP28kST08eu//14lfH5SEenVIiKMH10+nriYSB54+2Oqjzby6y/lEBMVHs1v/nLOcc8Lm0mMi+aWs0L/loK2FHQi0muYGd+fM5akuGh+9vJ2auqaePD6M4iLCb1HywTas2uLWbWrgp9fNZHkEO28+Xj0p5CI9Dq3nTuSn181kXc/Kuf6P31I1ZFGr0sKaoeONPBfy7YxdWgSX8oNn3NzxyjoRKRXmj99KPd/eSqbiqq4duFKymrqvC4paP3ilXyqjjby0ysmEhGG5zUVdCLSa82dOJhHb5rGnoojXPPgSvZWHPG6pKCzdk8lT63ew02zhjN+SH+vy/GEgk5EerWzRg/kiVtncOhII1c/sIL8fTVelxQ0mppbuOf5zaT3i+XbF53mdTmeUdCJSK83ZegAFt92JgDXPLiClR/rAa4Aj6/czdbSan54+Xj6htHtBO0p6EQkJIwZ1I9n75hFWv9Yvvrohzy7JryffLC/uo5fv76Dc09LZe6EQV6X4ykFnYiEjKzkeJ69YxbThifzL3/bwG9e3xG2vaj8+MWtNDa38ON5p4dFx80noqATkZCSGBfNYzdP54tnZPK7Nz/iXxZvoKGpxeuyetTyHeW8tLGUO88fxbCUBK/L8Vz4NtqKSMiKiYrgv784iWHJ8fzq9R2UVB3loetzSYyP9rq0blfX2MwPl2xmxMAEbjs3vHpAOR4d0YlISDIz/vmC0fz22hzW7j7ElQ+8z56DoX37gXOOn7y4lcKDR/jJvAn0iVKPMaCgE5EQd8WUDP58y3QO1jZw5R/fZ8XOA16X1G0WLi/giQ/3cNs52ZwVBk8O7ywFnYiEvBnZKTz39VkMSIjh+j99yP1v7aSlJbQuUvn7hhJ+9vJ2Lp00mO/PGet1OUFFQSciYWFkal+W3Dmbyya1PsT1lkWrOXSkweuyAmJ1YQX/sngD04YP4FfXTA7Lbr5OREEnImEjoU8Uv5ufw0+umMD7Ow9y6X3vsWHvIa/L8svH5bV87fE8MgfEsfCGXGKjdV6uPb+CzsySzex1M/vI9z7gOPM9amZlZra53fh7zazYzNb7Xpd0tLyISKCYGTfMHMbfbj/Wk8pK/ryysFfeb3egtp6b/3c1kWY8dvN0BoTZ43c6y98juruBN51zo4E3fcMdeQyYc5xpv3HO5fhey/ysR0SkUyZnJfHiP5/F7FEp/L8lW/jmU+upqes9j/s52tDMLYvyKKup45EbcxmaEu91SUHL36CbByzyfV4EXNHRTM655UCFn9sSEQmoAQkx/OnGaXzv4jG8uLGEC3/9Dq9s3ud1WSfV3OL45lPr2Fh0iN/Nn8KUoR02pomPv0GX7pwrBfC9p53COu4ys42+5k39tkSkR0VEGHeeP4rnvj6b5IQ+3P6XNXzt8TxKq456XVqHjjQ08b1nNvDa1v388LLxXHx6ePdj2RknDToze8PMNnfwmheA7T8AjARygFLgVyeoY4GZ5ZlZXnl5eQA2LSLyf3Kyklh612z+be5Y3v2onAt/9Q6Pvb+L5iC6DWFLSRWX//49nl9XzLcuHM3Ns0d4XVKvYP6cgDWzfOA851ypmQ0G3nbOjTnOvMOBF51zE05lelu5ubkuLy/vVMsWETmhvRVH+MELm1m+o5zJWUn87MqJnj60tKXF8b8rCvnFy9tJio/mN9fmMHuUbghvz8zWOOdy24/3t+lyKXCj7/ONwJIuFjW4zeCVwObjzSsi0lOykuNZdPM0fjc/h6KKI1z+h/f49+c3UVBe2+O1lNfUc/Njq/nJi1s557RUXvnWOQq5LvL3iC4FWAwMBfYA1zjnKsxsCPCIc+4S33x/Bc4DBgL7gR855/5kZn+mtdnSAYXAbcfO+Z2IjuhEpKccOtLA/7yWz+LVRTS2tHDRuHQWnJNN7vDkbt/22/llfPdvG6ipa+Key8Zz/YyhYf/InRM53hGdX0HnFQWdiPS0spo6Hl+xmz9/sJuqo41MHZrEgnOyuWj8ICID2BOJc47NxdU8uWoPf121h7GD+nHfdVM4Lb1fwLYRqhR0IiIBcKShib/lFfHIewXsrTjKsJR45k8byozsZCYMSSQm6tTOCBWU17J0QwlL15dQcOAw0ZHGV2YM4+65Y9XbSScp6EREAqi5xfHqln0sXF7Ael83Yn2iIsjJSmLa8GRyhw9g6rAB9I/97DPwmlsc9U3NHKxt4NUt+1iyvoRNxVWYwYwRyczLyWDuhEEkxaunk65Q0ImIdJOymjrWFFayurCSvN0VbCmpprnFYQZZA+I/Cbb6xhbqmpppbP70/7uTMhP5wuQhXDZpCIMSYz36KXq/4wWdnjAuIuKntH6xzJ04mLkTWy8kP1zfxIa9h1hdWMnO8lpiIiPoEx1BbFTkp94TYiKZPWog2al9Pf4JQpuCTkQkwBL6RDFr1EBm6TaAoKDH9IiISEhT0ImISEhT0ImISEhT0ImISEhT0ImISEhT0ImISEhT0ImISEhT0ImISEjrlV2AmVk5sLubNzMQONDN2+httE86pv3SMe2Xz9I+6Vig9ssw51xq+5G9Muh6gpnlddRnWjjTPumY9kvHtF8+S/ukY929X9R0KSIiIU1BJyIiIU1Bd3wLvS4gCGmfdEz7pWPaL5+lfdKxbt0vOkcnIiIhTUd0IiIS0hR0IiIS0hR0PmZ2jZltMbMWMzvuZa5mNsfM8s1sp5nd3ZM19jQzSzaz183sI9/7gOPMV2hmm8xsvZnl9XSdPeVkv3trdZ9v+kYzm+pFnT2pE/vkPDOr8n031pvZD72osyeZ2aNmVmZmm48zPey+J9Cp/dJt3xUF3f/ZDFwFLD/eDGYWCdwPzAXGA9eZ2fieKc8TdwNvOudGA2/6ho/nfOdcTqjeI9TJ3/1cYLTvtQB4oEeL7GFd+Pfwru+7keOc+3GPFumNx4A5J5geVt+TNh7jxPsFuum7oqDzcc5tc87ln2S26cBO51yBc64BeAqY1/3VeWYesMj3eRFwhYe1eK0zv/t5wOOu1QdAkpkN7ulCe1C4/XvoFOfccqDiBLOE2/cE6NR+6TYKuq7JAPa2GS7yjQtV6c65UgDfe9px5nPAa2a2xswW9Fh1Paszv/tw+3509uc908w2mNnLZnZ6z5QW1MLte9IV3fJdiQrUinoDM3sDGNTBpB8455Z0ZhUdjOvV92ecaJ90YTWznXMlZpYGvG5m231/vYWSzvzuQ+77cRKd+XnX0tr/YK2ZXQK8QGuTXTgLt+9JZ3XbdyWsgs45d6GfqygCstoMZwIlfq7TUyfaJ2a238wGO+dKfU0rZcdZR4nvvczMnqe1SSvUgq4zv/uQ+36cxEl/XudcdZvPy8zsj2Y20DkXzh0bh9v3pFO687uipsuuWQ2MNrMRZhYDzAeWelxTd1oK3Oj7fCPwmaNeM0sws37HPgOfp/XCnlDTmd/9UuCrvqvqZgJVx5p+Q9RJ94mZDTIz832eTuv/OQd7vNLgEm7fk07pzu9KWB3RnYiZXQn8HkgFXjKz9c65i81sCPCIc+4S51yTmd0FvApEAo8657Z4WHZ3+zmw2MxuAfYA1wC03SdAOvC87/sZBTzpnHvFo3q7zfF+92Z2u2/6g8Ay4BJgJ3AEuNmrentCJ/fJF4E7zKwJOArMdyHeHZOZ/RU4DxhoZkXAj4BoCM/vyTGd2C/d9l1RF2AiIhLS1HQpIiIhTUEnIiIhTUEnIiIhTUEnIiIhTUEnIiIhTUEnIiIhTUEnIiIhTUEn0guZ2S/M7JtmttnMVpvZOK9rEglWCjqRXsbMZtHa1dp6oBj4D+B3nhYlEsQUdCK9z3Tg77T2gt8IvAKc4WlFIkFMQSfS+3TUb19zj1ch0kso6ER6n3eBS/F1iAtc7RsnIh3Q0wtEehnn3FozewZ4GEgGEoHrva1KJHjp6QUivZSZnQd81zl3mde1iAQzNV2KiEhI0xGdiIiENB3RiYhISFPQiYhISFPQiYhISFPQiYhISFPQiYhISPv/XhEpxHIj8SUAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.figure(figsize=(7,4))\n", "plt.sympy_function(f(φ, m=1), x_values=(-1.05, 1.5) )\n", "plt.xlabel(\"φ\")\n", "plt.title(\"Bulk free energy\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compared to last tutorial, this bulk free energy has also two minima, but at different values. \n", "\n", "Another complexity of the model is its anisotropy. The gradient parameter $\\epsilon$ depends on the interface normal." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAAaCAYAAACpbxvnAAAACXBIWXMAAA7EAAAOxAGVKw4bAAANiklEQVR4Ae2cjXVUtxLHYx8KcEwFIR0YUkFIB4FU8EIH4VABBzogVJBHOgBXwIs7gA4g7oD3/wmNrNVK947W9+7eXVbnXOtrNPprZvSt9cmXL1++m8qdnJxciNe5eL6biuecfObEOyfvOWVy5H0jgaMOb2TxrYSk8z80fr0caq+HhvJeuqG6tp3Xwnw6FRBVwCTxbM8miVnw7pssprKBQ+Jz1OEhabOrLVfS/ZuREh4aWHjpRqrbanYV88lUOwoJ9x8152fxu95qszasbE68c/LesLnHYp0SOOqwU2AHRC7dv1BzPmgs+7PVLA8NZb10rXp2kV7DPMlEIcav1CAEO7hl20Wja3XOiXdO3rW2HNOml4BXh6L7XbWf6bsbUbxSH/g4PaLlcjxUGahdHyT1+0MLXw8NmvPSLUnLJeZbTxRieE8NRKjfDwl1KULoxSv6X4X9tb4fxtrXy3spMjniuJGAV4eiY9WZJgbFmTA4snghO9mLO7qbVm8WOmQZqG1/SCo/SpdPWtLx0FDWS9eqZxfpJeYpJgo6zNmQQHfR0FadEkAXXtGzavxF7XvU4mnpvbyt3CH7ksnDfRo4PTqMNvFI7fol153Suae7VPr3efohhg9dBmpfWABLlyct/XloKOula9UzdbrwjC5+S8zhMpuC+l61vhGgDKRvR2iWlN2FV4byp77RSSI2sIv3roQS9c2KaVanelhle2U3K5YO5h4d3he/ByVP2cmV0s7UbngcujtoGUiXHCF+pK+0FOmhoayXrlWPpd+m39IX9b3Rx0L5N330zaYrMd+BUol/y+PrcqqUWZcK92KrPSfeOXl3KWVZxI+XBWcYjVeH6i8cR7SOJBhgGEQP2n0jMmBcQ89DY6OHBlvw0kE7uZO+rsU0LNpk5yxkmhNgVnnCHHYUWUZv8KEKXEcQvWV3QT8n3jl570JWt6pTxogh8shhn1yXDtXGCzpdbKu1k8UT3zfhDlwG3L2u7RwLxXpoKOKlK9jvNJow38lhSOnsDp7pwz/X91aTQPOJmPJZOblfeUT+bH2Y3XDUQx1pxs4wfAoUX1+UQMPsFhzGqQCdmrrBCZ/fRDO2knPhjRjs+ID7iZWzaNVVcy7etYL7lib5MAnYYHhX4U+SUXrxpnxkZ7p4qLhNGKWux/igZx4SUNd/9cEHveN+0vc+rzek3u6PS4dqj+ECD7jOlYZd/xWrd/eJSB888bin9mxUNudTC4NZvDkam8R9IzKwo0TuYG3MKuXnoaGMl67kv8t4wpwmCimeDsirDS7p0qA8gpIO/HmEJmTTCRTgtxb/Ef8wMURju5R/pTTOA+mAYGBwTh1G6W/JU9pL+UwKvCxZGbyV3joKEHlyXrzwD/zE95m+3xUfmjCpwMs7gekJCAN3CmFQxheepz3lp6IVDgZ3BsWkI6Vxv/VGmMLW1mSlNGzqneJrunHywVDvi5aVDfLlYjxMSEoj/kH+30pLtqI0Lg/Bh51AQ95T0bQ6urKTG9WheFs/4TdD4MPB+6ny7K4OvMEpDZtes9eYnTwwi1/SqeLd+h4pc0/5k/zXBPFxy2AEU2o/AdEuTQY2tmEXpusVzIp4aCjjpSv57zKeMJ+CQgpCEBg5g7h3kqDouT5PB4SWCeCj+Kfdg+LUi6NT46BZ6/hKowO9EE46HVvBBwpbGUWDs1WrxWv+KF7xZTXMQGOOekxgllbzR3nXCnnSIqYwOUh+DJTvlYasduVoay5/sPAgIk/zYPPyYbB/ENse+CpskwODVnLCgO54tvpEH4sJ6rhMBMOBQR2KN/ZKW5+L91WFldkzv27lWAqbZJK09EqR0P/gmexX5bDBLn2PlRFe+t0j0Q1iqQLMEmN5rwzc7RDfJcrgOjYdu2g5Dw1lvXStenaRfoNZBgQAlEQAY6Wj8bGK5qK7+YmGlROdskkT+Z+JLvBv0Sr/ItL8WqOJeaw4FAz1EmByY+XFln0QQ1ZuEK94sc0MvBRmEKKelGZ5pS+apiwor4/dVM/HDspwwJvVdGqj4qO4RIMea3XC799GHjuDVE8rrLLoC9mHgVk+eFb0oHhTJsZXNB4+6PkfK2O+0qiT/82T8CpOnQmHwmZXKS2nz8OxbNM+lG/9pGoPykem/xY8sSF+jJow5mHlgS/YtaVDr69L354yosEOXfo1LKVPeX1N21NekoHCrnaIbpEyEC4bt1Z0kcvEQwO9ly7SztZvI38m8KYOrX055ju0QA5jZieRtr4iviZjxHloYGGrGAyn5YxmiCcGheMs+Zk+cCNUdhs8Y1074lBe7oZ4B7qi3axI2QWNlhNdkyaWtzP7HM9oWO06ExGyKXc11Ef78x2aojdO9SZ93qSGFWy4G1B+OMbJ88bCEc+l6MDDqt2ODDE+t9uAT9n+obqQl+04TC95Wqus0bbykRvHpGt0sT1nyq/KvMVQ6djxc8vfRN/eMuAWLUdQTJomH6va67tk4MUUK12qDM4dQvHQwMZLx4KiakOS6cb91tGOGknCfBpzMfAwIGJMfLVSlTQ6b2JWybckM8ofLaHiGw1YWo5BiQ4fhCmc/MSeH8QwQfD6hIFzyHnxGg+bQC0+5PfyHuKV54X2KqHUyVz15XXXwqwoWVGnO4qcKA4QeVIKK48diLmN+RiDmi9c/Jo2Pz69iHT/q9EXaU2ZZu1q8bGJcuwuq6jyu/KSeRN995T5SwDG+kmJMcQ7ZdCDaakysLHIxqaaXDw0lPPS1erYVVrCfBoRXG2I5FrljFmThToudNTxoEYkA+S4iXzoVi6poVe+dXYGFwzQOiXZTBp0TlbWZpwk15wLLwVjp6Beu6Cs8cvT3LzzQh1hBrHSnZUJW4jXdjG53Lk/Mn0hk/MMU75Q6OGTsegOsuP0XmZ7dNjaFbNYeRlt3QVScsrlVpbZRN+eMgx6a32srHwk3iODQUwLl4HZbq0NJiIPDbReOuO7BD9hPo1o2Oo8LpHFwbJMzuPvFakO/jlRDPMihm3vyiCvOB35KtL8LP+x0i5i3LzXCtAJbaXIS6RykCRu+Vau9Hvw2qprjKfV0cPbyvT4pjQrU8YtfW6fwbTUD09VceiAwc9WYKy+W4NhDx+x6XfRtnh15T1ia+owTgDYQj7ZBVCq540C7MirRwYDyHNZlWSlfst4SU+8pCnj0CB39LTm6FP6vujjXmvNbSiDEkMZX5QMikaDbeyExUMDWy9dAWGWaFX/lZoS5jtkygDeyTh4EYHB01mCMSl9rIMxwGNcHEVQpumUz7HRDyJ4LZ/zelYlPPfk8jAMLPJ5LUIeTw2NH41iRZgP2OwgmFDIg47Oy4WgDVCKVp0br0r33E9QWQ/vKrhGorWJtuaOuE2wefrcYfTDnRAPH0yHnLF/0seEjj7RCXb1RHT824AwkCqJRYG5UT4iZGFhd1Hs8uBDXRgwK3gci4afVFd4mvs1KewIWZBgl0ZnWUP+mA7hxVNtzoqxRwY9JofPZf1K87ogq4x4E333lPmsusrBOlSvNnCHAa+hPu2VQQ+mxcgg0wNBbLR11GikHhovL+M5iy/d0mfRvS2C+WkC7eOhSO3I9KZtIqBDb/ypEl45NF8F3Ib3HGW9eEXHINh8AVPD5uVdKzuUFrGsyFhpo68WWjxVloFu5bVQi3Yf09U2OkJ6SRTj6RXZUJtEO2rPokF+THpMRoOvqZQPluqrJ6Vf6ONHiCv9D3p9Xfr2loGvvsGXT8qnfUyyK7jyeKQZlIFoRtshmkXKgLbKsbMa7CceGi+vXL61sOraar/N29Y0hBrQWpqYMUulTlmjWVJaC6/SMdjQOeSzYsVSXIOLtU/0s8hCfLkETpOWwhjMYGc3TDVfZdmNDA4EtXL7kKZ2oUc6ODKyj3smV3tFN6kOI4bWRIEe1vKUNqhv5WOfK20aK2O6Ex2T22B/Hcs3XmO+B5NolioDcA0uxiL2QRpk5KVzyHNr/bbEPMW/GadjXqqRe/GvlbXVWsOrNDoeqx/uQewXtpw59xxb2KX7LLIQRgaPu/pwDHpd2L4WO/y/khM7AjrUipO8mv8uOies2Uee7w1Hm+JYihU89sVji/fCsXKcKzrsjtd71/KTU3pT34ZRxPxAFr7BDZXJaJgIWWTkR7mWbY84+A3VCs5E0BlwYlqUDGiicDOh8rJv5Ugzb76Hxssr57uE8FrbJIgpdhWcHVd/KDcF/6l5COsa3pgWVnIKj/7YsIWpxrtFe0y/ve3NIcNt6lB1YXMb9Z1NyqnM2g8Xcxkqf2Pbz/n0hFXnomQAdjkmr8ETBQ+Nl1ePvLZBW7bt1jsKMbRVCCtpLj8W7zRbsuKcBe+cvBcv2AMBuG0dqj7uKbqerKoMx2rhf6R5xR7LnKuu2sWll80sdEuSgbCwm+D3OM1XbB4aBOWlm0WoGzKtYT7dkNdKMQmUbfNzVcC2dvFuTrxz8l68YA8E4A50yCsyBqcex4BvL4tGy4k/iyOOUhY3SUTwi5CB5MQxIf8YdWiSGKWhTR5ese2L8VqYJ9lRWCtVCaucawm5ev5pdEvx58Q7J++lyO/QcWxTh7Gurh1Cj/zFn+Md/s3Nyl1ID4+5aZcgA2FgsTv4A00PDbLy0s0t1x7+LcyTThQ9gI60RwkcJXCUwFEC+yGB/wM4FTJYOLhLzQAAAABJRU5ErkJggg==\n", "text/latex": [ "$\\displaystyle \\bar{\\epsilon} \\left(δ \\cos{\\left(j \\left(- θ_{0} + \\operatorname{atan_{2}}{\\left({\\partial_{1} {{φ}_{(0,0)}}},{\\partial_{0} {{φ}_{(0,0)}}} \\right)}\\right) \\right)} + 1\\right)$" ], "text/plain": [ "\\bar{\\epsilon}⋅(δ⋅cos(j⋅(-θ₀ + atan2(D(φ[0,0]), D(φ[0,0])))) + 1)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def σ(θ):\n", " return 1 + δ * sp.cos(j * (θ - θzero))\n", "\n", "θ = sp.atan2(ps.fd.Diff(φ, 1), ps.fd.Diff(φ, 0))\n", "\n", "ε_val = εb * σ(θ)\n", "ε_val" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def m_func(T):\n", " return (α / sp.pi) * sp.atan(γ * (Teq - T))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, we just insert these parameters into the free energy formulation before doing the functional derivative, to make the dependence of $\\epsilon$ on $\\nabla \\phi$ explicit." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "scrolled": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1AAAACqCAYAAACj+q7TAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAgAElEQVR4Ae2dj7XctPa2IYsCwkkFl3QQQgUXOiDcCoAOYKWCLOgAUkGADkIqANJB+CoATgd87+Mj+SdrZFueY89oPK/W8tiStra23q2RtK0/fv/ff/99z84IlBB4//33v1P4Q10f6fpT17eqL7e62xkBI2AEjkbAbcvR0DmhETACRsAINIDAgwZksAgNIhAGOD/IYPpa12cS8UbXmwZFXV0klf2Jrk9XZ2yGqyEg/XyzGjMzOikC19y2nBToK87MbfhlKN/t+GXoyVKWEbABVcbFoe+993kGwgv5MSyYjdqto+NV4Z7LaPx1t4XcR8HeSlc/76MoV1eKq2xbrk7LZyqw2/AzAX9ctm7Hj8PNqRpAwAZUA0poWITUWIpL99KwhkU/WrSXSvnl0amd8CQIBAP3Tw2WvjpJhs5kbQTSduRa2pa1MTS/MgJuw8u4NBfqdrw5lVigBQi8rwq8gNyk14qABqq8NeaN/4eqM3HAsys4VMYfVKB3Kt/3uyrYjgsjnb1T8T7ea53cser6ol1D29IX1g+bIuA2fFN4N2PudnwzaM14QwRsQG0I7p5YhwaOPVG7NC5UPt6IMxgfGIhJOEbj77qi8RiXIbHUL4bBgyWA3wunb3W/CndOjJQ3e6EeC++vrwLsrJAB+7Ts1D/+p79kpM16VQb+d7ttW5oFfmeCJe3QoA2nmEmc2/ERvZ8TI+V91e34iEoc3DgCNqAaV1AL4qlx607jW3uQKr6fimcTe43GyhjC35OcA4NI4f9IN38r/HGqI4Uzi/WHwn9Mw/f8fE6MlHdn+Arv9/eMcalsoeycjNkbUAqLM8Uc/tJ8HQx152FahlJZHdYmAtJf8204yIV65nZ8ohqdEyPlfbXt+IRKHNU4Ag8al8/inRkBNWzsMVl9gCO+D8X32ZmLl2ZPOV+nAeH5Iw3ucuOJxh75S8bfHwrnyPdrcmfDSLoBa/ZCxRnBa8KdMn+Vll14xJmnQZ1tEZSt2pYWy7pHmS6oDQf+s7VRF6T7s2F05e34BVURi5oiYAMqRcPPAwTUQXKUd788Cr8ulgit4b5Yg8kaPFSmokEUyloyquIR56U4DCuW+l2FawQjDNl+FuYqgL8r5FvdWJLElbtSWE5zNr/qzZZty9nKdWUZN9+Go49G2qimq0YjGF1rO9503bBw4wjYgBrH5qpjQoPK0r3f9Pw5l555q33v2ZXAi6VurTgGc7d6C3Yw6FRYaRkU38XClWag3pb43JHv87cBjNhD83Sf6I6XSrj/qov9Hn09DP8tErX0/xoUQjLyEmaTtmWQkT2bIXBJbTggNNBGbaaLtRg3gNFVtuNr6c98To+A90CdHvPqHNVJMaV+b4OllCGDGPHmDXbRKZ49Pg/zSKWZ3GsSOlZmdHCPdP2lNP3BE4pnqdzHurhTtjj4ey26uPyIt4YYbFN8GIRxXC00P+liwBhnhj7R829pvvKPOuVF2qeiR65ZJ/ri/qdSQtGeTYcleU4VtgSj+8qkvNA7s4EHm8fvy/uS0gsH/q8sIf1FdXl0Cd+562SoG4vblkvSxRJZz62PXFbJM9f2XnwbTplDPTzYx5rjEWjdjpeAWTFM+nA7viKeZrU9Ah9sn4VzOAYBNSbfpYMg+TmlpjNIuKdxY/xn0nyk+BvxiQbMgI3CPxwEVHjEj46XN8ufKX1n+CnsB10/y9/td9K9m9FRGI0lb9APll5V8sH4+1i0vLXCiGIzc2eoKQz/O90ZSNYYoND/rWvWBd4M/jDaJp1oz6rDSeE2jFyC0UpiRN2hx9GXAivl1RybgDf/I/5//J9fjAnZQp3Uf3Jx2zJWnksPb0EfKYaSZ/dtOOUN/xm346nys+clGGVJj/VedTt+LGhOdz4EHpwva+c8hoAaLr631C/BkZ83fp3RFIwEltVBM+rm0ogPsz3PRMegc013I2bp22XkZAlgGlaTXy0fDCRmj/pZLj1HowkjrcaR120NoWgiz9L+p55F0M+l6rAvx5EPVRgdybuULOoOPV6do77r4jQ+ToTkhcL/U/2LOujxuPI62ePQykPD+qhte+egrOVz6jYcueP/w+34uBarMBpPvjjmqtvxxWg5wfkRUKfL+mBfjWCgGvFEFzMXvU7kZ1DEDEsaxheQOR2vD0ufa9KQXhezQ0Uex4aLJ2VgxozZKIwIZGUJRJ+P/JSJb7/0Yfmz4mv40AFybPiAj8LI85s8vOQX3awsMZ1oMQjnsG9Wh5IdnbPMa8nFcs8BvlN+8Z7FaCr90rhQJnQy+I8s5bMXeuHAEtNuCW4sk/zN1sko4zXdW9dHkG+XbTj1TG62jWpZR5LN7fiCPuma2haXtX6scl+sPqAlsWsKgeeSpl+Co7eUNJSlJWa3CucNUb9vSM+dq02jynMrWpbyrbK+O+T7RkIwFc93aDhemoEbM2jV7gg+ceq/Oo+MECxrHZjzxn8qTbM6DHJX7fWqBaRAV4NRIdnRQTdHp7zghPxvEZ//WVYMToFEB2n70GydzGS/Fm+T+jii7S3q6wg+p2zDkbmmjWpSRwjvdhwU7IzAeRH44LzZl3NX48vMxV+6WJLC4OjLmQGrSHbj8sMd4hK7fMBOhzM2cFyS5pX40JmUTptbCipv9ZgVKw7Q6VTH9Kg4ZoviMryj+SwVONBPYdmzlIzg+lDX3P6nS9ZhX95jHmoxoi6IPwMU/uedS/TPHgXiaQeYHcS91UX9OXhhQDgEcrkhcRe6398OG2FVc3jG1dbJRtXfqj6ObntVD5tvw6kLktPt+MyfohYj0bkdn8HS0ftF4EFrRdMf8gfJ9I7BlC42RjO4pVHfvQuN1lg5wSF3ceCYh0d/TRoGnfFY7pju2Hv6xjvyiMYc/qcqIzNSuFtdqQGIsRzdEj4xzX3uyDKHJfyRCze6bn4HOrwr4fG/NRiBNTOVL4LRhPHOyYmdE4bEs8SQJZ4Y1RyMAOYYUSUX61Gpvpfo9xJGveUgFu6pi/+57oAY18kUmvM/N66PJW0v9S7+9wD2Etpw5Kxpo+J/CPrcldoZ2qwpV5Nmzb54SpaauBqMKLPb8Ro0TbNLBD5otFRpQ8xAatESsEbLVCMWjfbYW3Q6qnSglHZcY7xr0sBzrvEf45+Hw+tJFhgHxuRB+VhehOP+tHs6/FnC5zD18pDflITZkDnXnSQoouLJhSHxpetwDoO5+BqMXooJ9fwLDSY7fjICYjr8zDxhGKQG0638Y/8NMCeeenNN7uCocuHJKWrgwRLaiMe118nW6kTL+ljS9l5iG05diG2N2/Hxf0YNRm7Hx/FzzBUg0JwBpU6fWafUsRxsqqFLaffwHAc9sSxx0PgwBoQ7/nSAmUYvScObsRpjLOU/9oyuvtMgrptF1PMjXeznYpkWjS0zCl350LPofuZSOPIyaI5ulo8I2buF0dO9KQt8yKsbPAZGzxX+ifKKnUEIPriB40PRHiwxVBhywZMrGodvFI7MY9+aumQdqljL3BEYobNn0svY//oLxX+ZSAH9GC1k1JdomOO/Cif8fhT2n+ri/xYd9ZTPCOR4XVWdjGA0fG9VH7Ntr+rWRbXh1IEj2iiStaojZFvdHYGR2/HVtWCGl4TABy0Lqz/0Q8nHYOq/Lcu5omwYM5S5d3RWwoHBesnIKQ4aF6aBbzS4+nyPeVC+8CkZKyzDivubetaiL9HGTfGluJxPiQZj6Jc+k4oHyfGrMKazZEZsMPBU3MFb/hmWF63DmbIVo4/AiDo+qHP818WHuk4cF3qMjiWmLOEbc+jt1VjknsOpuyrfoM4Wynt1dbKAQUtBzepD9Wl3bTiKP6KNalZHW1XkIzByO76VMsz3IhB40KqUGkhxhGpcX5sOploVeQ256Lx4g5w73jD3xoKwYZkOH4m9hVB+ZmNe66JBi24yTSTSnfwGg9kk7poef1Jh19gLZh3O1xoG/HE2Lx4YwQxlPF2qq9f4VaehG52BCnUemjUOQSHLPTrXyba0an1so4+12nCks47mdeR2fB4jU+wYgfc1CG+6eBogYQhwzPYag9umy4pwKi8na30cjaMosMIxKFkSh+Ntfb/UUXEMIDE2Oa2wn32ZSgMTnGjAl29Bzb3F7uj3+hMxFA4f3reM4mUdToAofDD0WRoJTp0T7v0MpeIxmKjf7E2jznO6191mKXlSJ1r2R7JkrX/BkMb7+Q4B18m2aoL1sb4+hGnXD67RhiOddTStI+HjdnwaIsfuHIFLMKCYIWGgxZ6J3jjYq17UKGEo8Z2hxWVV2s+XplMaPkLLuverd8KC/VivlmKYA2cd5ogc7xeWGEgciFCso4qPbcO1zFIfBabrZBk24dIduqH61RvwZcp1Q62PdfGM3ITrKm04/KyjiOr978Ky+XZcMp6lLbg/uuZwLgSaWsKnCsxG/n908SYpdzd5wB79oSPvZ5dqyxj+/IsGkSENM1B2dwhwcMHz+4JhHd4XwUF6DKfi7KjqL50yS1kX1fsB9yvxuE62pWjrYzN9rNKGI511tKqO3I6vCqeZtYBAUwaUGiz2Pvytq98DoWeW8+BnffO1OE6yY3C4xN0IP9ZtVznxZ/qdpU/eOxIQC/XvhbBZw6i0Dqtq4jhR0AP/A06aG7xUkZ+ZaWallx7yMZ7h/mNcJ9vSsfWxsj5WbsORzjq6p47cjt8TQCdvFoHmlvCFgdL/hBhHXz/SxUDpWzWM1caB6C/eCQemk99uVW7xZ6ngj6HDuXi81ixAwJ7vChVnPmrzsg5rkVpOJ2wxcmkX0pctyxldWQrXyaHCAx7ssT3pEr4ohfURkVj3HnC9dxuOVNbRurpJuQnbZtrxoOeztQUpLn6+DASKBpQqEpVoE4NFvJ+It5fcXEb9sJRGwAgYgd0i4EHTblXrghmBRQi4LVgEl4mFwMF3oFSJvpOB0y+LkT+e/tbNCKVxYwjOpOHIbZab3evt/ljeDjcCRsAIGAEjYASMgBEwAkbACGyFwGAGSoYNJ9j0y+XkZ//B42g0yc+ysv/JP3pkcE0a0TBti6G2ySyXeFc5ycE+oDdVxP9HxFHhq8+gSZa2z5P/v/L7yQhcDAL6rxaPPz+mAG4vjkHtuDRr6g0JpDuOzWc/be7oAzigqNQXcRrqaF+XM6rxu52vQck0RuAOgbXbAbi20hZYx5ePQG9AqVKxSRvjKJ194ohgjhDuZ4tCB/Chwop7DxQ/m0Y0dFovxWPVzuny1eESGAEjYASMwKkQUF/ko4tPBbbzMQINI+C2oGHlNCpaegrfc8n4KsoZjBwOcOBUvNRhOJXe5GHZYxjNpgnGF0v5oLUzAkbACBgBI2AEjIARMAJGwAhcBAKpAZUf7hCNm3ymCYNq7JtMS9JgrBUNsYtAzkIaASNgBIyAETACRsAIGAEjcHUIdAbUzExQPgMFSMw0TbmaNKw5/2yKieOMgBEwAkbACBgBI2AEjIARMAItIRBnoJg5Km2iRdZ8tin3l8qT0+R+0jCzNWeIlXg7zAgYASNgBIyAETACRsAIGAEjcBYE0mPM86V60aDKjRz8b0ekXZJmaingCHsHGwEjYASMwNYIhP2snFyHi0uzOaF1rO2/o/SvETACRsAIGIErQCAaUBgzA0OJgx7UiWIQlWaPfi9hszANfKPBVWLnMCNgBIyAETgxAsF4GpySqjCMqT90/0ztfH8q64lF2yK7PZVlC3zM0whcCwJuC65F0yuVMy7hw5CJbxlT1nyvqT9qXJ0nR77+gqEEkfycpPdaV2p8TaYhXXBTywYjje9GwAgYASNwWgQwlr5Ms1Sbz+ctaPf5VuBuHH1Z7M92UygXxAgYgcUIuC1YDNnVJ0i/A8X3mz7OOxMZR98o/FFA6qHiv46oKe6JnvkQLR+X/SUJH02T0GBo/ax0tvojKL4bASNgBM6MgNr1fyTC32qbH6eiKBzjiZdofFzdqwdScPxsBIyAETACV4VAakBh9PDl9d4QqkVCHevnS9MpzR9K83FtHqYzAkbACBiB7RGgbVYuvBwbfDA9MaB40ea9UNurwjkYASNgBIxAowjEJXzvqUP8XjL2s0u18mI8iXZRZxrSMANlZwSaRkB19Rtd/+j6d+JyXW5aixZuCQK82NLFy7VuqXaSFqMK59mnOxz8e6EIuF2/UMVZbCPQEALxEIko0ndqWL5Sx/ljDKi434i+ukMVf/ZLsRF5sbFWIYtJjMBqCKiuxlPI/ium7NnjYoaWD0BTj7vZ2iX1X2nsjMDFIaD/AnWe+v+96ntuWF1ceSzw9SLgdv16de+SG4E1EeiX8EWmaly6GaWtBoXiz1LBH90JR8R9bxEB1VPetn+ketoZSaHT/YH/RajDv+p50cxri+W0TEagBgHVefbIvlWd7w8VqklnGiPQEgJu11vShmUxApeNwIEBddnFsfRGYBsE1PG+1uDxM7jrmc30HJziN/HbwG2uDSGg+t4tUVV996qBhvRiUe6PgNv1+2NoDkbgWhHo90BdKwAutxGYQ0CdLLNRqbH0xMbTHGqO3wMCqvusGGCPrI2nPSjUZegRcLveQ+EHI2AEjkDABtQRoDnJ1SHA4PF1Umq+f8YeKDsjsFsEVMfjkeW98aQw6j57oeyMwKUj4Hb90jVo+Y3AGRGwAXVG8J11+wgEQ+krSfpTJu3TzG+vEdgNAqr3zLp+Uph5wqiyMwIXjYDb9YtWn4U3Ak0gkJ/C14RQFsIINIQAxhPfR0uX8PHMfih/BLohRVmUdRAIM0zs8/tFz/EkSmZcb3SxfJVPXtgZgUtGwO36JWvPshuBBhCwAdWAEixC0whgKH2bSWjDKQPE3l0hgPHEMr1u/1NWsupPVmTp7DUCLSHgdr0lbVgWI3CBCPgUvgtUmkU2AkbACBgBI2AEjIARMAJG4DwITO6BCks5NpFMvFljb2cEjIARMAJGwAgYASNgBIyAEbgYBEaX8LH2XWvd+6VL8rOc45Guv7incWOlnUnDaU434uPlUGMAOtwIGAEjYASMgBEwAkbACBiBphAozkDJsGENfPfxRKSVnw2XndEUNhD/FmiILrq5NOLzixI+Ex1r7e2MgBEwAkbACBgBI2AEjIARMALNI3CwB0oGDUvr/pfOMCnsncK+TmeLFPavwj5UWHo6WV/gmjSi4WSnl+LxrE/oByNgBIyAETACRsAIGAEjYASMQKMIlJbwPZesL6K8wchhlujvGBbuGE6f6mImaeBq02B8ibb7MKOefbrTAEV79oSA6jnHQfPCgP8Sdf3b/OVDDY3SndVdgoxnBciZGwEjYASMgBEwArtHoLSEj+98vE1KHpfY5TNNGFR8F6TklqR5JQYYYnZGYJcIBKPjB/2vmMXl+Fz+N2/SwtbQpPTneL4EGc+Bi/M0AkbACBgBI2AErguBgQGlAVI0fEoo5DNQ0PBGfcrVpOFtPINKOyOwVwQ+zwrGDO+T7P9WQ5OxObl3NRlVdsrvFycnV+FpMrR+T4NzK7lI36Vvhh2IV0t3kNABRsAIGIHGEBgYUJINA2psKV0+25T7S0XLaXI/aZjZmjPESrwdZgQuCYH05USczU3DKEvqH6M5d5nvLSODaxXiuWbjfALnubW5Qf7W7wagts/yrfTO4VNzrpZujo/jjYARMAJnReBBIfc4cItR0aDKjRz86VK/SM99SZqppYApTz8bgYtEQIbC48xYwIDA/X53e++9GppIe677ijK+VBm+PFc5nO/mCFi/m0PcVgahfftTRhQn9o66WrpRBo4wAkbACDSCQH6IBMbMwFBSg8dBDxhEpdmjfgCYlmdhGvhGgytl42cjsFcEOFDi4BCJrLA1NFmSk3sXy6i2hM8jvKKNOLm0znBzBGr0GwbZ9DN8VxDH/sCr6gP2iIF0+K3K9U7XT1P/71q6u6rRf0Zl1fqyR/wjXr4bASNwGgQGx5irUaGR+kMN3OM0e4Wzvpm36F8TLj97ITjqvDt+XH6W9TAweqawbmA0l0a0nRMdb6zg3X+0N0T5ZgR2h4DqO0bHw/hfKhWwhqaU7pRhx8ioNLQTfBJh9PMHsQyipY1hJuM/sU2Jca3eL1HmNbGs0a9oqP+9wSQ/fQ5Lv/hw+1Us6dwzBirbYKwwVr8W0K1eX/aM/xjeDjcCRmB9BAZL+OJARQ0MnVrvFP69PLxZ+i40Pp8pLP12E/RPdfWbwivSRP4f6+F19PhuBPaKgP47vCyYM55mac6NT005RmTkBcyPsZ0ZoYnBzEz/Wkkb05z7vqrMwrlvT89dsMr8J/Ub6g2nvPazTUG/vDyr2T9TKUa7ZFeAAZ81mVzGF7QzS7cFVlvwbLe2WTIjYAS2RGAwA0VGamB4g/SnOjYauEVOaT9fmk5pmPHCiLIzAheJAPVego+eJKn6/bVoGAzz4qGbaQ3+v+Xv9xHW0JwboPvIqLT/SP4vVebFbcu5y33q/IUVL6WYlelm/U+d/zH5zelX8axS+EJl+jDnrzg+zM4x/z/mcXvyXwMGKiOzzCxRnvyfz9FtgdUWPPdUP10WI2AE6hE4MKBIqkbmtRq/0QFhib3SMIh8q3T928USXRoW0tzsvdNMy+zn60NA9fyJSs1ytBdJ6RkYp0teZ2mStGd5rCnHmGBKW718b4zHNYULL97if6y28SIMqPvqV+kZdDPjeBHl3aIu7gUDlQND+SPpcnIMUUtXwnoLrLbgWZLdYUbACOwDgfwQiVgqlup9tdCwwRBaYjzxhpU38lfbYUawfd89Anw0l/o+WKakun+blLyGJiE/y+N9ZGQG7jYr81kK0Xqmant5GcUg9JJmY6r1q/LxsoAl38zAxlkKDGyuq3A7xwBj+IsKRVbRbYHVFjwrymsSI2AEdoRA0YBSp/arGpjPdfEWqcooEt3Szp43rD44YkeV6dqLov8LRtJzXdzZD8NMLnt+DpYsKW7gamgGCYIn5MlG62iMkTf5xoEpM8pRrr9CMk4/g6bftC8aBrUMgvm/IztpOCimX157rIzig4PPZFsS5Iz7J3i5MvkGu+N65p+lMose4ygaCujhL5WTPaadU3w38xS8n8qPIYXLdTrHB30y60leP+mCD/rFfaLrtzTfLvR+PzX6jTIhCzLdqHzU3Vch68n6EWgu+qbyXgMGLEt+qLKy3zO2SyW9TdJtgdUWPEsFc5gRMAJXgIAaOL5B48sYuA7cow6oqWBwyh6fT0/1f1JeDI7Jk72Hnf70zACNMF5+6Nb5edPb+RM6Dm75JtBEo2tQBxT/LtLf9y5e5IcRMMgj9Sue09liOSgDs+ADeoWxR5NBd3fP40/tlxyzMkeZRIvRM9AF6XX9HGniPdD1vGM4d7mlfHpdh/TUGxgN6kSIA1tkIg13BsEDHZT8gX5Uv4qP/w8OkRjwC2mRp6uPQQ7q8Si/Ao+DsuQ0x/olx4HMx/ASn2oMwELXonou+iYwkBzoDn1O4jZFp7jVsdqCZ6wH4t0E9lEe34dtjPEwHlvUgQ/0x7czAkbgHgjorSYDUgac7GnqZ3XuwbI2KUsC8wNfkAWHUYSD5hfJlb/dZ/b3D8mOvDe6nuo5f2PMAHotRx65DD1v5c3MCwPG6JCf79L1LtA8UlniQRzMkmN8pCeC9vRbP9TIXJABHKJuiEY/r8Urx564KVfLB8yfCqN+No+6oPzgzSC1XzmgsM54Ip5I+ZHtja5+FpLwETeqX/GhTsLrhXj3h6YkfGKdfSvaJwqPy7pjeEJ6+Ijc4tuvZpAf46Ob3eOexh2mvguZSfOR4lmifvR/W+mXYMB/oZdbaWfruWhawuA24EydmHJFui2w2oJnLFhj2EexfDcCRmBrBNQpxDe+vDFq5opy+b78zYH0yADtj4XX5NtC62FcD8KZwSH/HQwOBqFcB7Mna2Io/ui4y3OMr+IZjELTz1CltCGOgZceu5kRHrrZCt1XfaMa+BdnVEL+/UyHaBnYI0sflsg4mOEr0aVljLxFt/r/IZVPz0WZc1mCPOglzjBQZyjrAG/5makaxWsBH/TJSaf5rA959jM+gV8+Oxbrz0C2nFeStiiv8on/j4E+Ix/FM9v4T/QHfuA5OwMaePfyyd8Z4pGX/MzWHczwxfiQ12wa8UBPfT5p+ppnZNB1UKdjWsX1GOgZPVTXc9HCu5dNz7PlifnGe00a0VRhILrYNg3KEPOK9zE6ha+O1RY8KUfg2wz2EVvfh+1diod05vFR1h+k+Ph5vO7k2Aw61TzS/nogjdX1YqUGmcEPA1Ua5u7auj4onzi4HQyC03xFw+CRXr44kAlx3TIpPSM3hh+GBmm4igPiNI/a58C3ip9okWMweJYf+ZBpYOjLD/ZFA7FWtjXoJMOBzDlf0VAG8KWudIMu3aMe+0EY6eRGDSjFLeFDXgdL4RRGJoO6I/9g4C4/sybQFetPWj7RUK6ifgOPAyMulDPqNZdl1oASX7DrXgBEWeQflCHkQRmKxluIn01Del2ThliUoXRXWmSYxSDkA21VPYdOV1MYSJ6qejNGp/DVsdqIZ3PYl+qew653bGLdb6v7B2pY7IyAEbgfAgyuWEp3G6/7satKHZfDPZ6gjjTIN+ZYzsWAh5cpfLuFo7NZ38Uyqq8Ux0B2DcdyvJtKRuSZL5fqZFT4bcZjCd8s6areksx5BrxZZyDPARlRNz2NsB7Vk+KYsYruaD6RQekumR7rSnFngIj7/e42+VvUQ1KmMR7MluD6pYR33qrf56J6FSlDXtQTZEkddaZYj2vTCBd4sJQv1sOU/+TzQgwi/9p63iIGsR4f1PEMqAO6LbDagmcoR4vYZxDbawSMwFYIPNiKsfkagStC4O2pyxoGdOT7tJS3Bg3MyhDPQKzf/xJpFR8Hx93yH4XHgWxHorQMaDnJLw7ouvB7/CBHHDCNsgmDHWRj5qTk8sExNLN8S4zWCquQOWbFIL4/HTEEpviyDy3qBbxuYkLdU0N5CZ+ExeJHZtUwqpFlzkEzpQdmeUoOQ/37yjzy9MzSpP+9iGEdWPoAACAASURBVGUuL3UmxTLlsyQNxlrREEsZTjwvwaC2nreIQcS6VIYUnim6LbBam2eL2Kf4+tkIGIENEXiwIW+zNgLXggAb2L/ICxsG1nnwmn4OT+Ct+MD4kZ+BbxxY/lfPXyQD85j/Sz0wcI0zDs8L8jIgjvEx3bH335SwaOxlDOMAdSzfOOiKyXJ/DD/lfU7mKAsD+yfRE+4cKY4Dawbz8a09MzZxcK/HgVvCZ5Cw1hPqEB+27Y9Yn0lb1G8wjNBlagB2rJQHxjszt/0BEDN59NFKO4YNNKWBO/hOuZo06Kb0MoIju//VxTLGA3ckBnm9zv0c8tEMBlmhkaubjc/Cc+8B3RZYbcGzYexzjO03AkZgIwQ+2Iiv2RqBq0FAHTTfTXumiwEhA0kGuCzVqh18HoWV+LP87j9K/FJ3TkrjDSunj7EXpRuI687JZsR9qzty4RhMMrOQGinMOGFoEQcdA172VsQBvbz3chh0Nd+GYYDaLYfMcotyIF/q8EdjMQ0/5fOYzLkM6IGPlLMZP+rqhZ7/0oVBi97AnuWUX4vuZy55KTtGcXSzfESIYc0So864C3zIi0Ersz44jOZPlNfgFEOFYZBTfyNdRzzzM6Vf+HDSIHvyqHMYAxhNfEh3kLfCah3liHUiTwP/DscQcWB85Ankr0kDT+rbwKkMtyobskzV71oMYpnyfPDn9bwZDAaA3J3aOLZkMyWlHpfotsBqbZ6tYp/i62cjYAQ2RGAVA0qdB4177OBpWHAM0PIG/y7Gv0ZgZwiorjMwTA2Sk5SQwZsymhyEioZBGQOIojuF7OShdgJZmYWawokB/0E85QyD1NJguDQIK5Z1o8CizHleQQ8lXWFoHxjboi/RYlyhz1JczqdEQ5v8Sy5b9AtjysJeqK6+BD+GzmRbrvhR/QZ5H4sXBhRGHcbjGsY59Sl1S4yPmG5Jmr+VqFT/0EksX+Q7uNdiILql9bwZDJIC8x9nueOcK9JtgdUWPFW4FrGfw9zxRsAIrITAvQ0odYoPJctLNVB9Z60wjCm+McNm6YPB0Eqym40RMAKXhcBPEpfZmr5NUBvxRP5uxknPvHzhYoam5AinnenSi54BOd+4ygcypbSrhS2UebV8t2YUykXb/SJgS5YYUn3bTsCEO9BvSis9YbiNGm8prZ7pV6YcxsyAhnoguTGISkZO0chemAa+0eAqycZs3mT5KjGorefNYSD80Qn/aZYOj7oaupWx6mRZkWdz2I+C7QgjYAQ2QeDBClzpcL9M+aiRYnkGgxqWn9gZASNgBECAgWG/X0uDKIwl9o0814Uj/ke1H287X/ajcGZY3ikdy+Bod3hBUzu4z7gd510q83G5nC3VG+XM4Jd2O14cYU5bXuMG+q1JkNOAry74oF+eWcaYnkAYk2DIUH9yR9q+TijtwMiWH57dh4uThJNpEjryKxpQ4onhwFLMe7sF9bwpDELB2Qta81Kjlm4SzwVYTfJJIyt5toh9Wgw/GwEjsDEC76uxuFcW6jj4DgtLPAabhBVOB0znxXKQYqdzr4yd2AgYgYtDILQLr9QmdG/qg/83FYT2g+/kHHOc9UlxuESZTwVQrt8t81VeLAXk2P2BgadwDK5HIe/BXi7FYSBiKH4Z6yB0U2mIx4kGQ4tvQfUzqF3EXRwfzj553ZVMzWAAFkGeZ8Ki+BIkwQu5Z+kifYv31rBvESPLZAT2jMAaBhRvkOmUPlSj2XdkalyiAUUHN9mY7hlgl80IGIH/Q0DtAm/q36hNYAO53c4QOKV+lReGEss/J5fNlSBWWo75X5ROaTDwm6q3LWEgWZhd5oXp5KmKtXQlvbUU1hL2LeFiWYzAtSBwbwNqDCg1LrxhYsnDwLAao3e4ETAC14GA2gZmpll+N3qwxXUgsc9SnlK/yuu16hH76qpdkO+t0lWvjAhpbpTm5LNMcwVrAQPJQF/PKZKTuqilmytzK/EtYN8KFpbDCFwbAg+2KLAaFU5yokE99gOJW4hlnkbACDSAgAZZvPlnHwrthN3OEDixftkP1++rq4QSQ2iJ8cSsKQZ/c8ZTKG8LGDDr1O89m9BDLd0Ei6aiWsC+KUAsjBG4FgQ2mYFSh8bsE2/4ahrUa8Ha5TQCRsAIGIGVEVB/w4zmohmlJSKIP0sFOdykX6K+JP0paI3BKVAu52Hsy7g41AjsHYHVDSg1Jmy07T4EuXfwXD4jYASMgBEwAkbACBgBI2AErguBVZfwhTd1Np6uqw65tEbACBgBI2AEjIARMAJG4GoQWM2ACtPY/RfsQVBhfHODvVB2RsAIGAEjYASMgBEwAkbACBiBi0dgFQNKRhLHmPMV9vxULdam2xkBI2AEjIARMAJGwAgYASNgBHaBwL33QIUZptdCI/2mBqcW3eh6IqNq8IHdXaDmQhgBI2AEjIARMAJGwAgYASNwlQisYUDFD+mWAOQjhzagSsg4zAgYASNgBIyAETACRsAIGIGLQ+DeBtTFldgCGwEjYASMgBEwAkbACBgBI2AEjkTgg9p0Wqr3nWhZmsehEHyE8NvSdzFq6ZT+LK51+c4CijM1AkbACBgBI2AEjIARMAJGoAqBqkMkgtHxA4dE6PpMnNnf9CbPoZYuT3cq/5ryidcTXZ+eSnbnc1oErN/T4t1CbtI5H0yddDU0kwwcaQSMgBEwAkbACFw8AlUGlEqZn6b3QmEYEMxGpa6WLk1zyudV5GNwLaGfy5j89ZTCO6/TIGD9ngbnBnN5K93/PCNXDc0MC0cbASNgBIyAETACl4xArQFFGVNj6TYUOg2LOKRhU3SR/tT3NeR7KaG/PLXgzu9kCFi/J4O6nYzCC5E/ZUR9NSZVDc1YWocbASNgBIyAETAC+0DgqEMkNMBgJoc3tR9qQBGNpANEaukOEp4o4Bj5lOYHifdO5f7+RGI6mxMiUKtf0THIZk/goyAeS1zZG3g1bq8YqFzvpMSPZ9q2WZq0IqyN1dr8Uln9bASMgBEwAkbACEwjcKwBxeCBAeOkEREGIrN00yJuF7tUPtEze0XZJw3HKLHoMTSZzfjP1GAs0p/7fmnyro1XrX5Fx4EqvcEkP4YULxS+k56vYlnnnjFQ2dgL9Vi6zD8MruA7V0OT0K5aX/aMfcTMdyNgBIyAETACLSOw2IAKnffDqcEFBa6lOxc4x8i3NI3omaX4TFg9O1c5l+S7trzi9+klGRQ1+g0YPVO5OEyldwpnX9wbhX/YB+70Ye8YqHzdixLp8v0xFdbQkHZtrNbmN1Y+hxsBI2AEjIARMALjCDwgSp3y57p+GLti8tB51xhP3fImDUBG3+BGnue415ajIBvlel0ILwap/D/qugjjiQKsKa8wZlbmYsoeFFij349F+zTQ9zdh91aeh6Fu9eE7fdg1BtIlSzHZC5UfOtOrs4YmEK+N1dr8+jL5wQgYASNgBIyAEahDoHoGSoMJjuxmNuVbWAf/32Hg2OdWS9cnOPHDsfIp3aLleycuVnPZCS+MEfaRNGlE54CtoV/xYHnnr5dS5hyDNfx7wUDlYK/jR9LlYKYxxaiGJqXPn9fGam1+ubz2GwEjYASMgBEwAncIdDNQc2CoY2Z5Euv4f9Mzs1W8mcWQGmyar6VTurO4e8qHAXmrAdXooRlnKVSDmYb6wQD0ktwi/VKXdH0VyhrLiZHNdRVu5xhgDB/MNGaKraHpkqyN1dr8snLZawSMgBEwAkbACEwg8EEap06ZZVfPdXHnY7mvZTD8qPubEDb4RkrBmKilE7v/cyFfDLRonJA/ef8SqRLZ/gphjwJNv2mfQYXCGAhj2CE/fP4nPix7wR0l313S9+AxMBhD+OAW5GT2BceM3egb7DuS8/4ulVf0GM/RSEAHf6mM/WEiiu9mnkKpPpU/GlK5Puf4oEsO4CCvn3TBB93iPtH1W5pvF3q/n1r9RrmQB7luVEbq7quQ/WwdCXQXe1N5rwGDuCSTJcuxXcp1NkuzNlZr88sLZL8RMAJGwAgYASNQgYAGB+x9gZLB6T/cY9gp7sqPATL5fh7z0zMDNMJYQqNb5+dtb+dP6NiP9E2giUZXV56EhiPHB2HH+JUPeWEETPISDaezRUwpw1fRz12OE74YcHf3NO4cz5JjUt5UJtFi9Az0QHpdP6d0PAe6nncar7ilfHo9B97UGTIZ1IcQB7bIRBruDIAndRbSzepXvOJ/5EnOM+SHTF19DDypx7N1JtAelCXP41i/ZDiQ9x68qjEAC13VdV20zWAAZrrQ5yh2czSKr8Yq1IHJ+rI2v7QOiHcz2Kdy+Xm+7TJGxsh1wHXAdeD0deADdZzsZ2JAygCS08X6GR3iTuCY1fpT+fazTfIjDw6jCAfNL6LJ3+6zjPAPyY/MN7qe6jl/Y8wgeg0H/zz/AV/lzewLA8bokP/v6Anxj1SOuI+M5ZAYH2c5bGFO3ih3dgeHqBei0M1r8cpxJ27K1fIB86fCqJ/Jox4oP3gzQGWGtHMK64wn4gmQH9mYdYwzkASPuUn9ihd1En4vxJ+Zh9zFOvtWtE8UGfd+xfCcvvcjt3h2dYJA+TE8utk97mkc8SU3k+Yjxd+Iz73+2+KxBAP+C73sSjtZ1xXfGga3AWfqxZgbpVmI1Wx9WZtfWqAGsU/F87MRMAJGwAgYgfYQ0KAKoRgY8oCxwSCUazBzAt3al/JgIN7lO8Zb8QwuoOlnqFLaEMfgS4/d7AgP3YyF7qu9VQ28izMqUR7R9LMdemZwjyxpGLM3gxm+nCbySu/w0PXHwmv0zXnkDd/k+UDeGJfflQ6dxNkF6gzlHGAtP2Wdw6uGD7r8oyADefazPcTL5bNjse4MZMt5JWlH5RXv+B/pMUv5KJ7Zxn+yMDCdnAENfHv55O+M8MhHfmbrDmb4Yjz3mjSiQU99Pmn62melr8ZAtNV1PfDtZZP/7BhIhtg2Df6vKVZTNKFMKGeV+rI2v1iOwLcp7KNsvq/f5xpTY3rtdYA2Wdfq46lrx/Uay58uNWOgSsXqrlOAobziAHcwEE7zFg0DSAYixYFMiOuWSekZ2TH++HOQhmt0UJzmM/cceFbzEj1y9INnPSMb8gwMG/kHyxfn5NgqXnIM5C3lIxrKALbUlW7QpXvUYT8II63cqAGluCV8yOtgGZzCyGRQb+QfDNrlZ8YEumLdScsoGso1qt/A58CQC2WNus3lmTSgxBPsOuM/yiL/oAyBP2UoDsRD/Gwa0uuaNMSiDGN3pUeOWQxCXtDO1nVodDWHgWSarTtTNIqrwirFWmlG68va/EK9aRL7FBM/e8DvOuA64DrgOtBiHXigjhnH4IpldLfx6kK3/4lL4h5PZBVpkHHMsaSLAQ8G4be6OD6bNV4so+KkNAYm93UsxbtZwIQ80yVTnXwKu814LOWbJV/Nm8tbYswMBAN5DseIeunphPOojhTHjFV0R/OJDEp3yfRYV4o5A0Tc73e3yd9RPSTlGuPDjAmuX0545539fS6KV5Eq5EM9QZbUUWeKdbg2jXCBB0v5Yj1M+c8+L8Qg5lFT11vFINblg3qegFWkWYhVwq78uDa/JJdWsU9E9KMRMAJGwAgYgfYQeBBEensO0cKgjryflvLXwIFle8QzEOv3wERaxccBMgNyBm1xINuRKC0DWvZWxQFdF37kDzLEAdMkizDgQTZmT3KXD46Jr+KbM1rLPyNvmg2D+HSvGnEptuxBizoBrxsIgkuN5CV8Yvpj7syqYVAjy5yDZk4PzPSUHIb695X5pOmZoUn/exHLXF7qTIplymNJGoy1oiGWMpx5XoJBTV1vFYOId6kMEaI5miVYRZ5T97X5tYr9FAaOMwJGwAgYASNwdgQeBAnYwP5FLk0YWOfBa/s5QIE34wPjR34Gv3Fw+V89f5EMzqMML/XAwDXOOjwvyMygOMbHdMfcf1OioqFXYBYHqaV846ArJsv9MfyU9yl5UzkY2D9JA/TMkeI4cGYwH9/YM1sTB/d6HLglfAYJaz2h/vBR2/6I9Zm0o/oNhhG6TI3Ajp3ywXhn9rY/BGImn5huDBviS4N28J1yNWnQzcGLCJjyv9H1ry6WMh64IzHI6/bAr7yawiArNLJ1M/JZeOot0hyJVcp38Lw2P5g3jv2g/PYYASNgBIyAEWgNgQ8QSB30r+pQn+liMMhAkgEuS7VqB58iP84pD5bf/UepX+rOaWm8ZeUEMvajdINx3TnZjLhvdUc2HANKZhdSI4UZJwwt4qBjwMv+ijiol/dohzHHIBNcogxjzBikdksiE4IoA7KlDn80FNPwUz6X5C3ljw6+EwYcSBD19ELPf+nCmEVnHTa6fy26n7kUTtkxiKOb5SNCjGqWGHXGXeBDXgxamfHBYTB/orwGpxgqDGMcPUW6jnjmZ06/8OK0QfbkUecwBjCa/s7zV1iNoxyxTuT08E7r2MDwyImDvyYNPKlvB05luFXZkGeqjtdiEMuV54U/retNYZCBQh0dW7IZSadoarGKvObua/NrGfs5LBxvBIyAETACRuCsCHQGFBJoAMWgMDVGTiYYgzdlNhgE55mLhkEZg4ii21p++GuAiZzMQs3hxKB/QEMZwwC1NBieG6gVy7xi4IG8Jd5BByU9YWgfGNuiL9FS19BlKS7nU6JhAJ4vI+zFFcaUhb1QXV0JfoycdODe08cHxU/qN8j8WPwwoDDsMCDva5xTn1JXa3gcm+ZvJSzVv46fyhjLl/Lvn2sxEN2Sut4UBn1h7/7nLHmccrQFRZparKaYp3Fr8wu8W8U+LbqfjYARMAJGwAg0h0BvQDUnWZsC/SSxmK0ZGEcaVD9RWDfjpGfe7HIxS5M7wjAKuvRhMM73rfKBTJ5uVf8CeVfNd2tmoVzMdL0I2JIlhlTJECMud0X9pkTSFcbbqAGX0ur5YeZPvRgzg3jqgeTGiCoZOUUje2Ea+EYjLZUlfWZGb7J8lRjU1PUmMZAO0Av/aZYOF10NDQkrsYp5DOpDDEzvK/JrEvu0rH42AkbACBgBI9AqAg9aFaxRuRgU5nu1MJbYN/I8yAzNjxrovA3+/qYwZljeafDFMjgG+pxmVzu47/nc50H5Vst7n3zOlPaN8mXgy7LBeHGE+W2lPAf6rUw3IANjXfBCxzyzlDE9hRB6DBl0kTvS9XVC6Zjx6o1s+eHXfbg4STiZJqEjv1EDSnwZwLMc896usq43h0EoOPtBe8xHwKihGUk6DA46nasvw0QTvkp+rWI/UTJHGQEjYASMgBFoA4H3NdBpQ5ILkUKDEwbmr4Rb/5Y+hP2mcPZc8Z2cpcdZn7T0lybvKcEp6Xer/JUXywA5cn9g4CkcY+tRyHewl0txGIgYil9mdXA0TeDDwQEM0vkW1GAGNYnn49knrbutYQAWQaZnwuLgJUiCFbqbpIm0rd5bxL5VrCyXETACRsAIGIEUARtQKRoVzxp08Jb+jQZXbCC32xkCp9Sv8sLoYelnb4zXwqm0HPG/KJ3SYNw3VW9bw0DyMMPMHrrRUxVraGr1eE661rA/JxbO2wgYASNgBIzAEgRsQC1BK9Bq4MGyKpbfjR5qcQRbJ2kEgVPqV3m9Vj0qHi0+BkeQ763SjS7Hy9OGNDdKc9IZplyOkr8VDCQHSxw5SXJUHzU0pTK2GtYK9q3iY7mMgBEwAkbACJQQeFAKdNg0Ahpg8eaffSic+Ga3MwROrF/2ww321VXAiSG0xHhi1hSDvznjKZS1FQyYder3n43ooYZmJGmTwa1g3yQ4FsoIGAEjYASMQAkBz0CVUHGYETghAjKgmNFcNKO0RDzxZ6kgB5sM9lot4bE1rTHYGuFx/sZ+HBvHGAEjYASMgBEoIbC6AaXOmJPHeOMdT/ziY7fNDtxKoDjMCBgBI2AE9omA+6h96tWlMgJGwAicEoFVDajQMbGHoFteJD8n1n0kf1Mb108JsPMyAkbACBiBNhBwH9WGHiyFETACRuDSEVh7DxRLkVL3Qp4n6rSYjbooJ5mR23ucLkpr9cJav/VY7YFS+mYZ46yrpZtlZIJWEXAf1apmLNcAAfdRAzh276nte2rpdg9YAwVc24CiSKmxFJfupWENFHtaBBouUTzXzFnxeznTqR3bOgLWb+sa2kS+t9I7M+JzrpZujo/j20Ug7Y/cR7Wrp6uVzH3UVaq+tu+ppbtKEE9Z6FWX8OWCqxHgbR+Dlg8vaR+U5P5DMv/3kmTOsbd/HAHrdxybPcdI7+zPfKf/9eRphLV0e8bqWsomXbuPuhZlX1A53UddkLJWFLW276mlW1E0syogsLUB9U55sifq+0LeTQapYv4gwRhkXYzMTQLZqFC1+hUdR4tzGMqjUJR+b1+jRVtdrD1ioDLRJn0893Kklg7Qt8BpC56rV5AdMAx6dh+1A13upQiqk7NjELcP27S7LdSh0Ca5j2pBGXMyaCDxHpccg0X+uFyvw/Ukxi+9Kz1ve+mYOv6XcJe8LO0AjIdz8oqGN5f/1NDO8TpV/CXKvCY2Kn+VfkVH3eXwk/S/wX/i0xi29/teMVC52As12y4toFu9ruwV+/v+Z4SL+6jKNgys5S6qj7o0ee9bn0vphcFsHyWa1duckiwth+0ZA5XNfdSF2A0PpCzeoNIxvdQf5utwfSb/W11/KG7xQQrx7Qi8xOOSHPLWfi/nRrS/qoy3F1TAVWU+pm6cGatZ/Ya6y4uD/kO1Qcd8QLVmD82Zi3j/7HeOAR/Brvlw8SzdFjhtwfP+NeL8HISL+6g7Ncy2YYm2Vm3vE75bPa4ur+rN4vHLVoWr5DupX7cP/czTnvvo2b4n1KVZui3qyxY8K/8bzZF1S/gECLNOB99rUjgzLLyF/zBKrjDeamFgjTkGmZ8pDQNOjDMasL/lxyBr2oXyfilZqZh2EwgIKwY03wkrGvyLcDX6FQ3/hS9Urr7Ox8Ipjte6vGSY3EMT6S/1vncMVD6W8dHeTf7P5+i2wGkLnpdaD1O5Ay6zfZToJvsn2ivR0Ce5j0oB3uGz9Ly7Pir8D9w/7byPlp7dR11AmxQNKAwljJzHqcxSIsYQHdJjxfVv5FOa9Fn0T+R/qetFEs4A+5nSNz1TI9mZOqfSXtSBFwnOJ30UXrzFZ53uRRhQa+hXPKgfzDpeRJm3qBB7wEBlwEhmiebUiyBe/lTRlXDeAqcteJZkbzFMZXcf5T5qUdVUnbmqPuqa24dYMfaCgcpR1ffU0kV80vsWWG3BM5W5tecPgkAYR3z3iL0/JUOHNzk17o2IoB0sdRrhWcPvlDS8lby9EFlPictBXqonGNX8wS9pJmaRflVGXgY81cWLhThTgZHNdRVuxxhgCH9RocQqui1w2oJnRXlbJnEfpT2YUpD7qIpaqv/Prvsotw/d6qY999FVfY/+ClV0W9SXLXhW/LWbIukMKA0QPx6RigqKG8w+CTiMpOe6uN/oei0e7B06WPakuCoXeLI58jYkgDd84+A17tUi378CDSekQdN/rykolY4GmZENPv8TzVgZFd054gflDOH9LcgY90+wBGTyDXaf8IwPS2UWPR1PNBLA9y+Vsz+RUPHdW71QpE/lx5DC5bqa40PdYraSvH7SBR/0hvtE129pvl3o/X5m9Qt7lSfKhTzIdaMw6uUrXbjJOnJHcvgrHsx4HJX2kNswBJnFe7UlsleAAVg9VDnHXhhFgCfptsBpC54URnwvpv5F8NO76vdY+83/FTf4b6FbhdFXcKcfuFcfFfids39SEd6bbcOCnBfTRx0jr9LM9S277aNU9k36JyqXeF9EG3ElGEz2PegruEm6LbDagidluZT6F4Hv7uqYiqfkKZLBLHs+2OfS04RwllOsdiKZeDGIhufnMS8901AQxp9at86Ptd35EzpOR/sm0ESjq5c3hHMs+SAs94sHfOhkR+kU35/epWdk+yqnVxgnqNDRdvc8/tR+yTErc5RJtHRMA4xJr+vnSBPvga7nHcO5yy3l0+swpKc+wGig6xAHtshEGu4Mgkd1FuMC/Zx+qfPo9eD0yZAembq6FmShjk7yjDJHOYJ/cR0hX13FeqVw8F7l/wgfXVUYiG5UprS8rWEgudEbujzQcyr3FJ3iVsdpC56tYZ/ie9/ngBd63LSPEv+z909Bj7vro4Rt34foudinpvVENEv7lp7/PfmcvY9S2avbnFBfqvqnQJv/h6rb9oir5BtNE/R27z5qCQZT8kSZ07vom8FAspy0jwp1YLa+LMG/lmegawb7tE7MPY8OPAUUA+nBwFn+OLDtDZ25DGrixZcP1/6R0spPQ9kP5PSMPAOQoZfrK5qeYwMzGFQrvB/0pnmkz6JBhkF5s3jeavUDej2T+QAH+aHpZSRe1yjPlP8Wz0GeSZnTfIO8PebEyYEpDzmm6GOqc6rlQ8f0TypHyJc8Bwaq/BgQaXlYKjqoNzmf6IdO16gu4KsLmYt1ReGUF5nAgzqH8cY1aZwrHhlTmRfXEaWfTSMaZOnzieVecie9rloMZmWKeYtnUxiEcna6jDKW7mN0IXxVnLbgSZnkmsK+hPOxYSrbSfoo5XP2/inocq4N4z/ZtwF6pgI020dJtll587pBeXTV9i276aNU5iVtc3X/FOrVRbQRCzGo7p9axCCU9VR9VFV9WYh/Fc8Wsc/bnCl/0YASUN3gME+ocP5oKJV4BrNcg0FunmbOr/TMGnU8x2gVjzIOOoNIH+I6w0XPcaDbvTGSv+9QIn3pHtIVDQLo5XoDQs9TRsXgLYtoSdynHckbDLpOesF98u15rcwj8oB3fJuErinDAEf5RzunyFM0NXzQ04ERpDDyHBgz8pNnL4eeY73ow2Le+T2kndJvrNtFXSk9nfbA0JOfejBqQCkO+XqDGpmg17WojtSkEQ11aNRAzPEo+Umva7S+Kq7HQM9V5RBdcxhIptjmDPSQYzJGp/AtcNqCZzPYByxXa+PE7yR9VFIHptqO2A4NjJVYn8SD/9S9+id4yU22uYrv2y49N99H1cgbMczvSgvmV9NHqazVup/eeQAAEyZJREFU7UPESmkm+6dQp5ppI6LcY/clGIi2qn9qFQPJf7I+KuI9V1+W4L+AZxP1L+C9uH86MKDEiEap2FkonAEUA16U210RqGPv4hM7n8FgOeUnGt460YMUBzwhrltKpWfkwrADDNJwFcuT5QH9LB1p5OA/GDjLHyv8wLBROJgVO9Y0/62fJcOBzHmeoQzggI47o0T3qJ+BkaLw0c5ccWBRy4e8DpbBKQygB3VC/kGjKH+cES3Wi7R8op3Ur+LJ78CQg4dc1G0uz2QHpXR0en19SPj0YYH/aB1Zkka0lHGgpxSDuWelrcJgoUzNYQBGoayT9WaMLqSdrSuii/VmVt8b8WwO+7k6WBMvrE7WRymv2P4N/vupnKLZvH8iP7nJNiyT6aC9V/rq+pjyOsWzZDuQt5RvKENt37KbPkrlrmqbU8yUZrJ/glbuYtqIWgxEt6iet4iBZDpJH7WkvtTiv5DnxdS/tFzx+YFA6Z02cdERcGR5f0wzG7u4AhEV80/F38arT3z8Q9z8OzhCPWMXach/zP0Z5ZRsfC+EI7bfFzFl+UpxNCZT7m9F3kwRJHHw6g+uCOERo9uEjsclfLOkq3pLMucZUJl5i8kBGRHznkYYjuKvOAY10R3NJzIo3SUTdTPFncEN7ve72+TvqB6Sco3xYTkAbumpgwyc2eQZ3TF1ZEmaV8porp5HWQb3hRgskalFDGI9PqjjA1DuOmKCerotcNqCZyhHi9hnEC/zCqtT91FR9+funwBqtA0roFhq75f8bwssNw0qyVvK8Oi+RXXnIvuohe1DCbOpsItoIxZisLSet4jBqfqoqbrRxy3Ev09X8dAi9hVi35E8iJQCiMHoJxrw9cZTiKPDii4dDMawe92VHwYHfJ+WGEkuZm+Ih+7g1LsgN0lpWPnjxMEuYcywMejlJL/4pyK45OAfK20pvgsLFQmsmDkpOTq53M3yzROs6a+QOWZHJ9afehgCU9yeJniDV2pwpgOMJXxi3sfceWt58HHNEUY1+uWNZcnxn/g+1NVS/EGYcEpxy+OPqSM1aRjsHfxH8sxn/EswmJSpYQxivS3Jn8IzRbcaTkmGq/FsGPukuMseQ9tz0j6qof4JsGraME6zor+5mD6qQt60oizpW/bWRy1pH1LMis8X2kYswaDUvvPf6F3DGEz1Pb38epiiW4JVynPqeTWeDWM/Vf5B3AN8oSAYIDx/F64fdCfsa3Ui8S3ct/J/AV3qRDeolGlc5fMz0THTNTB+kEPhbwOP/+r+hcLoGFL3Uh4Gt3Fm4nlBHuSL8Wna9Pk3eYpGXEqkZxpw3Bi/WKHvqP6vgkf/Oe5zMkeZ6HByfDlSHAeGGAWxLjBbM2YkLOEjNstdqBt81LY/Yn2Gy6h+wyAJfaZGYMdO+fAfYNaVur/EpVjl6Y6pIzVpwB09DRz/B13/6mLpS9EdicGcTM1gkBUaubpZ9Cw89x7QbYHTFjxVkFaxzzGu8qvuUh7+i+foo1ronyj6aBtGZOLm2vu5/23C6iSPc/KmQizpW3bRRx3ZPqSYjT0300bM9VFHYlBTz5vBIFMScp2yj8qyH3qPxH/I5NDXKvaHko6EfBDC6ZgoTDrFHZPEATOzOb+qoj/TBT2NeTdgU3jtIDbyHNyVnuV3/1HgS9351gVW7iNd7Enq8tf9bYj7VnfyxTFYZAYiNWaYccLQIg46BsVs4O3LIX/JYajVfBuGN/zdMsaMSeRPvqnDH43ANPyUz2My5zKAPQY0G7SjDl7o+S9dGKroA0ypC1+L7mcueSk7xm50s3xEiMH8XPeu8wx8yIt6GGdBMYZ548wApncKw9BmqWGk6+MmHub0C6/X4s2MK/WJxhejiQ/pDvJXWK3rsEqIj6kjS9L8rbzyTgNd3apc8Jmr37UYLJGpCQwSHfBI/WRwNefG6LbAaQueLWI/h/lY/Nn6KP1/WuifwGWuDYvYjbX3S/63kdcp7mPylvKe7Vto70io+576qNr2oYTZVFgTbURlH1WLwdJ63gQGmZLG+p6MbLQvq8Uq5zfl34Jni9hPYTCMU8WlofF1t6GSzfxzG8sxLIqHTSicuEF6+dmpyWD/bBgHuYoyn1OuY/IGX139yXbBzzraWXxFW6NfDCiMQYy0yQMZFI8sg8NEohwKf6KrdDjG4jpCHrpm6xU0ukZP4lMcZZuti4FuEgPRzMokmuYwQD9yzMSNHgqQ6HCSTjxm64poZnGK+QXZVuGpfJvEPi2rn+fbrBwj6bWmDaPOFdv7pfUxz38Lf5CpKO8W+W3NU+XZrI8S79n2IZYvyFHsn6CRa66NCOWb7KMCzb37p1YxCHJN9j2JjifparBKeI2OZyJNkG2VOthi/UvLWfM8O+isYbIXGimUmZd+cE65gpK7P7SemR0hsDhgVzgzeH1HoGcq2uiAdivclCcNY5XMW8mwBd9QLhoMcI0Xe9EmG9woi+gO9BvjjrkHGYodFDLpOohT2GQdUTx1bFCmuTRRdtFh9A3qb4zjPhWX0tU8i9dkOUJ+LWKATIAxWWeIr6Gbw0o8ZnGa45HH1/AM8jdV//Jy2H+UAXXQhknX1e19Td3ZWi9L5N1alrX5h7I10UdJFvrIgzYglllxLbbPo/1XlLvmrrJVtbuNYlDV9wTZZ/uyGrygkZusL7V8UropnkH+g/qp8EndKX6zMVIqe83z+xDZ3SGgZU50RG+EyYeEyI+ieJvHHiuWDjKwZfkeU5lFJxqU/yhEMkgbpS0yuGfgMTLfM8uTJVfZePtK4zJwwpjTFmed0g/0O5tghCBgzPI+3thQRzh44zfJMVjKKjrqDqdB3ureO4WP1pEoo4i/VDr4dm4qTULD4AqDPV3SGtODG99sG8gY0x5zr5SpGQwoo2TGyOSUycllmbV0NbjV4FTDJ6Wp4SmaprBP5ffzcQhIp4M2TP6L6qOOkfc4pM6TSuU7ex8VMJ7tn0BItM20EZJl1T5K/Eb72bR2tIRB0MlJ+yiVnzakqr6kuE091/JsDfupMpXibEBlqEihrLV/FQevwc9+L/ZS8e2XpUdZZzls771EmbdH5S6HXL9b5qu8aMAxuHtDqDY/peX0yUXplIb6ydrpA6c4jKeT192WMAAUycOA4ZmwmNyXWEt3AHRDAa1h3xA0Fy2K9HrRfVSQ/6L61FNWmFy/W+bdUhshWdxHSdnCwX1URaUXTquOkSqyPCCxAZVBIqXwFoRZqOJANCO398IQOLV+lR/7oNgkXe1oGET8VuniZtjZtCHNjdKc3EiaE64VDCQHb/b4lhhv20ZdLd0og4YiWsG+IUguXhTp1H3UxWtxvACn1q/biM5oaaKfli7cR43/NfoY4dTEGOlBL5EfOgQ0uLrVwwspiOVQdjtD4Az65VRDGsUlDkNoifHEgIplac0ZT6HQZ8dAOmCZAjNPc8ZTFd0SZZ6Z9uzYn7n8u8v+DG3Y7jBsuUBn0K/bCO0dPnc/7T5q0b+yiTGSZ6BGdKbKjIV7q8bsYD/JSBIHXxACp9RvyGvRjNISKMWfpYI/ho53SdKT0Z4bA+XPC5HZjy7X0p0MuBUyOjf2KxTBLAoIBL26jypgs4egU+rXbUQ3C7V4VmNJPRPGk/204t1HLXhxvCb2S3iltDagUjT8bASMgBEwAkbACBgBI2AEjIARmEDgg4m4QZSsY87dZ6kQS1xYXjT7Nlc0dheAgHV7AUqyiEbACBgBI2AEjIARMAJNIFA1AxUG2HzfqNuXIT+nAPGRUR+00IQajxfCuj0eO6c0AkbACBgBI2AEjIARuD4Eag+RYG1o6l7I80SDb2aj7C4bAev2svVn6Y2AETACRsAIGAEjYAROiECtAYVIqbHESXW4NOwuxL+XiECqR+v2EjVomY2AETACRsAIGAEjYAROgkDVHigt1eMjsqnja+i43+9u5V/NUHF8M/umHgWKfhlgOcV06Nr8pnO7jljr9jr07FIaASNgBIyAETACRsAIrINAlQFVyIoDJSYPkSjsrcGQ+plwDdoXHw2+Nr9CmRx0h4B165pgBIyAETACRsAIGAEjYARGEKg6RCJNGwyZhzKCvk7D02fRMPPERys/y8KZuXqj8A/T8LnntfnN5Xet8dbttWre5TYCRsAIGAEjYASMgBGoRaCbgdLAmYMEBsZOyiAaS8GQmTSeQjpO53ua8uBZfN6Kx0P46PnHPH7Cvza/iaz2FWXd7kufLo0RMAJGwAgYASNgBIzAeRGonoHSQPxTifqZDJ9vETn4/8YoWlIEpXsn+l+VbnQG65z8luS9F1rrdi+adDmMgBEwAkbACBgBI2AEtkagag+UBtgsvWNvzAs9x2OvMYCeTQkY0jEThaH1S6DlxLf01LcpFoO4tfkNmF+px7q9UsW72EbACBgBI2AEjIARMAJHITAwoDSY5qCH57q43+h6LcOHpXZvQhgf0O2d4m57T/IQBuUvFfSDrp903SgMA+yVLlz3Qd67x/nftfnN57g/Cut2fzp1iYyAETACRsAIGAEjYAROj0BvQGmAzRI9DCQOfxickid/9aEPCZ//Kl1c3oeh9a3iXocisoyPZYDf6MYR539xF323PFDPvVubX89YD+L9kfJcZMyl6aeexftJUv4p0s3jEgxPplsKBQa6cepicX9dItdsXanhB03qxP8q9JuW2c9GwAgYASNgBIyAETAC2yLQGVAMNJUNxs3BAHtJ9oEPRtiLEeMhLt3jMAlO6uuNJvk/1/Wz0vXLAtfml5ZFvBnY9wab/LPGXJqe55k0Hyn+RnkMjNGcx9Z+yXAO3WI4xT1uUeeDoga5auvKLL8Bc3nE/yr0m5fbfiNgBIyAETACRsAIGIFtEegOkdBgk4Ese5tYrheX5b3T4H/JSXkMWiOfD5U28ulLoPh/8CjuQz0zC/V1amAo7F+F9WnX5kfeuMCX71h1s0/yY8w9lj8ekAEW/5O/N+ZIl7qaNKJhCSMD+apZLtE/FD3LJZe4L8U/zvQdpBPPqJOT6TYKobyZ1eTjyfmHmKMOwLnXd0zHXWn7uhLDp/hFmpCWMjen31RGPxsBI2AEjIARMAJGwAhcJgLRgGKw+ruu3mDQoPfAAJoroga4GEBvlZZjxwdOcRgH5IORwmCe54/TwX8YNGMQdAdOrM1P+TEwZzYD4yidfZo15kibOvGZTSMayvxSefW4pjxO8SwZTqpblfX7WC7lPWVAVdWVWn5Jnlel31hu342AETACRsAIGAEjYAROg8CDkA0D/T81WL2N19Lsg7FAMgyxkmOWB4fxFJd15Uba34q7gWhtfvAMjkMy4mEWMR/kIe/UIRsGwIELss2mAUslZilfLO8BrxMEnFq3s0VaqNtZfhnBtek3K769RsAIGAEjYASMgBEwAlsi8CAwH10CdkTmzMyUHHtivg9GRYzPjRbCGfCnbm1++eEO0bjB2Eldb8ylgeF5SRqMtaIhVuC7RdC5dFtTliW6reEHzbXptxYX0xkBI2AEjIARMAJGwAisgEA0oFjO9kXOL5kpyKMO/MEw4sCEsf0uzHD1y+YCg262KWHW+9fmRx4zM0E1xlwiavdYk4b9T8VT6HJmG/nPpdvR4hyp21F+MeJK9RuL77sRMAJGwAgYASNgBIzACRD4gDw0oP1Vg89nuth8/5suZmMeKrzfzyJ/jWOW6bX4cDgAxhQGEQN4PqSb7gOKhyrks0340xmTtfkxcxTz1uPAIWs6C9UbcwOqoacmTYflMNnpfGfQbW3hanVbyw+6q9PvEnBMawSMgBEwAkbACBgBI3B/BDoDCjYMtHW715Hb4oFx8jgYUOxFYYnWwSl0orsVDbQlI6XfQ7U2P+WHS40k/NGgwnhLHf7UmEvjlqRhlqpUzpTfps+n1G1tQWp1W8svobs6/SZl96MRMAJGwAgYASNgBIzAxgj0BtSa+WhwzCl63Ul6E3w54ptZqc5oC0bXL0qbD4Ax7tbihzEzMJTIr8aYS8uxMA3GUzS4UjYX+Vypi1i2AdYxML2vzM/6TcH1sxEwAkbACBgBI2AEjMDqCMQ9UKsznmOogTPLA9/JePmOS8+fKSxd5jfHYhBfyQ9DhmVeuYvGXBeeG3Pyc5IeSxNTg2AyTZLB1LKyhGw/jwEv8EGvYPezLj5UfJRbwM/6PQphJzICRsAIGAEjYASMgBGoRaD7DlQt8R7oNBhnWSHfnxrMdIUB/qNQRvZ/sUenc4p7ooc3uvpvVBExlYZ4nGgwJH4Wv3stj+yY+WcWAeFt/c6iZAIjYASMgBEwAkbACBiBYxG4RgOKmRBOBJxbYniAqQbnny9NpzR/KM3Bh4UPmDtgFQSEt/W7CpJmYgSMgBEwAkbACBgBI1BC4OoMKEDQIPu1jJpFR4tjPCnpW6Wr3s8U0twoDR8PtjsRAtbviYB2NkbACBgBI2AEjIARuEIEzrYH6sxYs+/qq4UyYAgtMZ7YL8W+LhtPC4Fegdz6XQFEszACRsAIGAEjYASMgBE4ROAqZ6CA4ZgZpUP4xkPEn6VkP8qAGuy1Gk/hmDURsH7XRNO8jIARMAJGwAgYASNgBCIC/x+Dir6STIFV7AAAAABJRU5ErkJggg==\n", "text/latex": [ "$\\displaystyle {{φ}_{(0,0)}}^{3} - \\frac{{{φ}_{(0,0)}}^{2} α \\operatorname{atan}{\\left({{T}_{(0,0)}} γ - T_{eq} γ \\right)}}{\\pi} - \\frac{3 {{φ}_{(0,0)}}^{2}}{2} + \\frac{{{φ}_{(0,0)}} α \\operatorname{atan}{\\left({{T}_{(0,0)}} γ - T_{eq} γ \\right)}}{\\pi} + \\frac{{{φ}_{(0,0)}}}{2} - \\bar{\\epsilon}^{2} δ^{2} \\cos^{2}{\\left(j θ_{0} - j \\operatorname{atan_{2}}{\\left({\\partial_{1} {{φ}_{(0,0)}}},{\\partial_{0} {{φ}_{(0,0)}}} \\right)} \\right)} {\\partial_{0} {\\partial_{0} {{φ}_{(0,0)}}}} - \\bar{\\epsilon}^{2} δ^{2} \\cos^{2}{\\left(j θ_{0} - j \\operatorname{atan_{2}}{\\left({\\partial_{1} {{φ}_{(0,0)}}},{\\partial_{0} {{φ}_{(0,0)}}} \\right)} \\right)} {\\partial_{1} {\\partial_{1} {{φ}_{(0,0)}}}} - 2 \\bar{\\epsilon}^{2} δ \\cos{\\left(j θ_{0} - j \\operatorname{atan_{2}}{\\left({\\partial_{1} {{φ}_{(0,0)}}},{\\partial_{0} {{φ}_{(0,0)}}} \\right)} \\right)} {\\partial_{0} {\\partial_{0} {{φ}_{(0,0)}}}} - 2 \\bar{\\epsilon}^{2} δ \\cos{\\left(j θ_{0} - j \\operatorname{atan_{2}}{\\left({\\partial_{1} {{φ}_{(0,0)}}},{\\partial_{0} {{φ}_{(0,0)}}} \\right)} \\right)} {\\partial_{1} {\\partial_{1} {{φ}_{(0,0)}}}} - \\bar{\\epsilon}^{2} {\\partial_{0} {\\partial_{0} {{φ}_{(0,0)}}}} - \\bar{\\epsilon}^{2} {\\partial_{1} {\\partial_{1} {{φ}_{(0,0)}}}}$" ], "text/plain": [ " 2 2 \n", " 3 φ_C ⋅α⋅atan(T_C⋅γ - T_eq⋅γ) 3⋅φ_C φ_C⋅α⋅atan(T_C⋅γ - T_eq⋅γ) φ_C\n", "φ_C - ─────────────────────────── - ────── + ────────────────────────── + ───\n", " π 2 π 2 \n", "\n", " \n", " 2 2 2 \n", " - \\bar{\\epsilon} ⋅δ ⋅cos (j⋅θ₀ - j⋅atan2(D(φ[0,0]), D(φ[0,0])))⋅D(D(φ[0,0])) \n", " \n", "\n", " \n", " 2 2 2 \n", "- \\bar{\\epsilon} ⋅δ ⋅cos (j⋅θ₀ - j⋅atan2(D(φ[0,0]), D(φ[0,0])))⋅D(D(φ[0,0])) -\n", " \n", "\n", " \n", " 2 \n", " 2⋅\\bar{\\epsilon} ⋅δ⋅cos(j⋅θ₀ - j⋅atan2(D(φ[0,0]), D(φ[0,0])))⋅D(D(φ[0,0])) - \n", " \n", "\n", " \n", " 2 \n", "2⋅\\bar{\\epsilon} ⋅δ⋅cos(j⋅θ₀ - j⋅atan2(D(φ[0,0]), D(φ[0,0])))⋅D(D(φ[0,0])) - \\\n", " \n", "\n", " \n", " 2 2 \n", "bar{\\epsilon} ⋅D(D(φ[0,0])) - \\bar{\\epsilon} ⋅D(D(φ[0,0]))\n", " " ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fe = free_energy_density.subs({\n", " m: m_func(T),\n", " ε: εb * σ(θ),\n", "})\n", "\n", "dF_dφ = ps.fd.functional_derivative(fe, φ)\n", "dF_dφ = ps.fd.expand_diff_full(dF_dφ, functions=[φ])\n", "dF_dφ" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we insert all the numeric parameters and discretize:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0MAAAAWCAYAAAACa0c2AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAUCElEQVR4Ae2di7XctBaGk6wUkEsqADoIUAHQAblUEOggrFTAgg4IFfDoAKgASAdwKwg5HXD/zyPpyLZsbz3smTmR1vLYlvfz35K85dfc//fff+9Ny/37979Q3Qda/tTxl9Pjfb8j0BHoCHQEOgIdgY5AR6Aj0BHoCFw6AprXPJGNX2r5S8tLzW1uYpvvx5MhET/Swd+0/Kr6r2LCvt0R6Ah0BDoCHYGOQEegI9AR6Ah0BK4RAc1zPpHdP2n5WPOcV96H6WSIO0LfieC+J7iEtYx/T3Ywo/OFGR52/uwr9lj7maT0xLqzVEnGn9iq5UfJufEytf9LbL/qmYh+o8XPVvH5K9H8rXUoER110FCgC0E9Vd3+ige5r13NY61/iOl13GqjOQ4ZMvGbdkfBNvZ/kn2/UhGXyPdFjERD26Chg8c/WjytNk9FsoeJvkWe56ldSxfYcUUCe/5wa63ufcaPCv56W6HFj2+9rdq+6iL/fRvEx/e1jNrgNTpX61Mt/96YyT7a5vda3lU79G2zqdoIA+TOxqamyiqFRbZmt2HxPpL6zfG90sSzsNfgcg6Da+218EfxxkXTefocWOTqtPi+JPOuYrLhrym3WZJxdH1NfLHVyr8DnTk3dXaSk3+o8xpPwJ2Kdu75RTUM1n/5/UtYyx6cZOIT28lJmuf7vojrt7ZFD98bLY+2aDmuQvJKYh50x9s6tilPNOhDWLx8M5HDiRK64I+28Xtkq/ahG9mjfWKG7E9imWyrQM9EJxzTNo2Axx+DT9q32JgVB4tMZ+Motq6OycxnExutGDH4xFhPtwf8RGOSF9tQsy19xGkUd+cr2M/6nOqIU2gPNbpb8cqezfae0iU+2mCIp7bBnr71Xor+Gupke5VPufyiL8K+BkvppC+NxpsaeTGv5JrGppjnnNuytzjeztfN8f2c/i3plu2r7a4GlyWde9bX2mvhFw1te9RvtL94nt7T35ayLb4v6btmTGT7ah9Y8dmU2yzxH11fE19stfLvQJeVmzpbZ3OdkBAvERwdkKk+AfdcC0iHZMrZSt0siUzwDwOTaHGeJBu+zcmQaNDLCWw6qGXJEz8nUa8bmbMEEB1a3iRspz50KLa1zGyHlyXB/ws8cb32qeOuVIi99i02ZsXBKJNBZpbwq25IlCY25mA0wxhZKiGWbGtJYTbCPLahZjvW7eWojk6MYaMYOVtJRMMk1vMcvZYNWe19ap/48SOFM31iNCmf8ubuSx7tiT6SXHLlLdFLfpVPVn7RVWG/ZP8l1Ms309h0IbbWxvvQsaYWM2u7s7bjWnta8dfaa+UXXdZ5upV/e8qx+r5kw7VhInurxl7xm3ObJcyOrG8QX9MYadVjpQMjlazc1PFc5WSI5/tI9keJofZBISuZEj0Bg282oYgbno6TpELLJCEk0DGNA3RT3hq/lyeapC+qJ2AiO01ctA0Oswmg6jjZIiNMArQNbqM6L2e6Ft2ij57WyTPHwSgT/0YTM/SpMBCN/NS+FaPZ3RcnE11PIn9M8jx9zRq9WlKTPt9+RhN9Zy8dfLWd1thUwosPWsDNbJdoSXpn7Ut1nCyyZJXYvAdPrU8l/OLJxn4P31vIlC/msamFvloZJfGKdYr/sLEm1ttie63d1eLSwr4cGbX2WvlFZz5P59h/Tlqr70s2XjMmsj177BWPObdZwuzI+gbxNZ3nrXqsdGCkkj1HEA/xGeWYD5B0yUXO8jGH/7D2dup5QxIpCldg9ihfSt8hX9GTLyT+FN5vmZbXVIiGyRmF94fei3iGyujHy6KK92Ju5MfonaOINmsT/LW0jgMd6BP588vEpxeqp7EOJTq2iZFsnH34Q/x0ltc6xntE4Olx2pQ3GNDgR7pT7elTJzq07UjVK/HcRPvXugn2qeJ9Wzqe4rmUuiWbrT7V8l8KDqV2NB2bSo3I4CuO1znGmgy/akmLcalVXMhfa6+VP+c8XejK4WxW35cMu4uYLPlKvSm3WRNw8LHa+Fr5W9Nxt6BJbvowB3AN7FytJlH1yWSKnU/WbX5wQLJ4Z4CE72lKyFKdO7mQKPNyeSrBXGI11Us+VwGaTrKcTPSDGy+P81jUkJhrzUcVSKLe0TItvFBMYTL0t2hvX/YaqsMPdx4o8cSHRvdKsjn2uRYmVujmKv0s+V6zUTyzInp8WY3DlkzsEA0fwWBy+0bbtBts/F3HwscxtG3GSLypwgcm/MSDzlMrL6VjsU76hlgnCIgRcSX+o6K6WYxGBJU7wrqo/+WolQ7aCCU16TwdObVtv91k7fT6cYp+xd3H1bHCiketT7X8TQBaEeLsYwykfCrcQr85VaV/rfg57qyxKa2xrtZqb2286NuSYRrf6zxa57b6uy7l9mgtLreS6rasftXam8OvmOecp+sAqOS24Jfj+5I5l4KJxd8lH3Lq5a8pt8mRWUJr8bc2vlb+1nRLeDg9q7mp4x3l3NPJEEn3LDGDUQr8lfqPtQsdCwkrJzaSniF5VSOIE3JVz4szFkOQYSrigZZkmcSZRPFrLU2Ls4tHgDZ9yFAMNsOX5OBxfvylNYmGT3hJ1JhoTgvYUhZxkhxoOM7kcIid8wM+yoeqD3dLdAzdTA7CZEM0FhsHYeK1xsEkU3Y8lUzaFv4zCWXikIptEUZONo8RTkuRvKmQ0n2H44BRqYxSPulGb1b/K9TlB5ubBP/aBClBbquSb/QH4v1Ubcv3r1XmTDxqfarlX/WlwUEeM2WcZax6oYXHO7cmkub25LD2ZlrGJk/bbH1wvLH73GONOT4ZIJ+9HR8cxyp/ZevsPJ2B9SqpZHPu5MIpFzwfq7+G8/0aYwZ+Vb4v2bAnJimdGf6m2LPrFAdrbpMt28KQ4W9tfK38relGMMhfa24KHxeDGRdDeeC3JIhJBh32ma/zax17om2u1vtPOH+k/Z+1z6QBgdym4gq3aRIhOh7f4pGrpSsoEjkuyNaCfu4c4Mj/ZBf2tiwvJP/blgIlj0nPjZepbTAiSQt3n1TH4MWVhDAhcvHwydxa4ogcYpEaAHlHZprIIPN7yQ8NQTSbNopnKKI1xcEq0/mJLVyBBhva2iy2kpeNkfPxuXinGHB3KFue7GpZfNvldvqhRb5n979KA0NbS8jhJN6kKN4MhuD5TD76vrMpuxCPWp9q+Tf9yiUQftwR8he9YMfGtbEHmuFOqzDMGs/FZhqbBgWNf46Ot/Sddawp9NeK+tnacaFftfaW8q+dp61Yz+hcnx0mQMKD3OV31aUu/s14C/Ar9X2m21XsgsmSsgJ/l0SZ6hUHcmowW81tTMIKiAr8rY2vlb813YCO/DXlphCLlpyQJ6dCHv5AO0+0kEBwEvxARFyZHxXqtMR3EjiR+YkPEyO/PeLba0e6/YTlJ9m+BqzZBMmh4f5gZqgjBC/e/SF5G4p8GjqM6r7QwqQoxjTedhzDlVsCyUR09Kih9m8cUYqPj0KA2YeOZmk1s3FKWBCHkUx8lczPJYd3tPCDiS6JA/bNYqvjuRghf9aeVTeUAnmetcUaXyjmpP1EflW/a4m0v0o0vBfXyCufyH+qtvWNW/wjX41UrE4OLD4djUmO39zBHsYMYecn603bp+S3GJtyfKqlbRKvM481tRik+JvgkhK8U12tvcX86kvJ83QjPzlfhgtqamfkaZ9JJ+fQVqXY9yUDdsZkSe1h9fIvK7c5zLC0otr4Wvlb06W9Ua36weYcQTTcjHlHseJpqecPVfFKFSQPnPz+1Hr0r6yqGxUd58q9P6FxjIlRvD+ir92RvmHCIB3TxP4PycZmlniilq1SOhg4PpKOKjlTxZLLFRq+8Da9A+bxwrfg11S/+F84mYHG69Cx4S6SeIZHWnx9tIbH64mqw+aAq9VG0ZnjYJUpS0he3w0WaUP+fCt+EjAmbbPY5mAkfrCZYae6UDLlBb4GG/jGlYy1GDVQcz4R+KZYYoCfJMTG0Ocoq/E5kZh/wZS2Q4IwlNb41vpUy+/92mM9wWq48DCpa6XWNDa1UlYjp2W8JGt0flHfWBzfa2w+grclLtdgb6m/ivHWebrYfclmDOW8PE0yOafMzp2likp9X9K3JyZLOs9Qn53bnMHGQWVtfK38rek8XmpP5tzU87AWHzkmORgX4e895IeiCh7TIpEgiP6qNYemhQQzXInQ9vB1MxydEjba55E4DOcxjL10MHDgB1dw4sLEz9fzGT4/24xp1raRm0r2GMQoTOjWCkHmEbiR37KTu1jvqz5MhFQ3bRBMctG/VLxdVhtz4rApU/aCAe9njXzDWNVx+5LEIZVEQxKXJEaOgGPgkFPW5OXIWaR1scL/HxeJdMBhRMIU7p4Im9AG3XH66xAbrfEVTEdJl+rOWRhT8HVafGw53qqg54gJZq1Ptfyt8FqTQx9uGZtYl3VsinnOub1XvHYfa3YGbS9c9jK71t4sfo3PlvN0ja+0H8r0HPqP6vz4OhA0+MnyfUnfAZgsqT6s3p2XW+Q2h9ksRbXxtfK3pgOjnNx0wFQx4s4dN3LCjYoHw5HbHxJk38Fua92WCzJCpknc1iNXM1kZFXR0HqGadnhvJ+BWFZJHLbxwzeNaYZFQr5u6kIRmKOOxkwB2xEeiET6dTGC08DW1R55G20zEWMJVbo65eu5ihYkQ9SoMvHH5WjuPYpnuIPbwzojHzWQjPFqscdiUKf3I486Bj6N2RwUsvI34bcYIKQ4nNjkxzEquvJmAugriT4kvKpxq3K/sw//ftHwtrGh7POPKI6lDcce5ssGXCTkOVsgj0bykwt3R1PjA1ZjQBxoZfJTvtT7V8jeCKy3GtS3GnsX2meY011rHJrPAnQmr4nXmsWZPaKpw2dOwBdm19pr5FXP6j+U8vWBqVnXqHMf5o2Ux+76k9GBMlszYvT43t9ndIJuC2vha+VvT4R25pDU39WiQf8B3WxQ4rsQPi2qHq8x+f7rWcW75jv6oSPv8wVjyjy6n/H5f9HRU+Db/NFU0TL5Gf1ipfRJ//m0p1Gsbmck/t4v0Yj80zNqD30vbosPG2Z+CenodW5Wn4yT6JKtBl/bxB7nxH6Qih8R2sIu12w/+IUMFecyCiZNfuJtFAxvFxdFTH2KjbY97+JNP1VltNMUhstPiN5OC4LfHSXXo4sMHMW4mjDyP+H0bCf77Y6yRryXo1nYSc0e72bZi2Vvb0kWSudoOdZzYsQy4s47lap+4B4y1zcn3TUyztC063w42+18sQ3yL7d3JxKeZTNXRZuM25/VzZcaPPdUYSwftaYaB6lb7O8fh0zKz3ds3XYu2yicLf6xT9EXYxzKs29Ll+84qbl6e6Evws4xN1W3C2xivC+0tjreL3eZY4+xK9qHY/tztEn+9Dmf74lil48W4eB2l6xK/au018pvO087+4niLnzEf/pBLgKUKY9no/JnCWDRZ/Vb0xbEW7yGYpPz0dbn+RnzZY690mXIbZ1NxG/A2pta5/oq+OL7ot/DvRGfOTT1OsnU213mIBxmFx+dGdyq0H67eW+XIIO4IcDXDP6q1yCral6LljznjR9joWPGnqUmokIm8mUzH+46O0UApv6mOR9RIfmZfGtMxgGKgYbBANyfu8F8lVnmS/TeyHP2NZGA3fr+LvVoPRdu8J/NYO3zKFp0Uvpw3xRY7kEHnnJaZ3+IfPu040T96J0w0VhtNccCoDJk8mvlMLHzdDlzABP9n/4UkmVaMxD6UV/pF3u+n3fFvjjzR+rb1iPiwP5a2vSc+2hSxY6FtUWiHxI0vNU7vPNJWuVs5bQPwUf6rBex8gX6J1tMMa+ePqf/BIBvpe6v9x8lM9j/xcjeSto3vFO5wTdthNcZgJR20efoJcR/ak+qn2OrQbcnFw3Hm+DTrm5KxyY+eTOyL26fzya8Y582PG5bgJx7L2FTdJrxD8brEXvFvxsvJnfUB1ZvGroi/VRwHt0v8tbQ7h2kOLmf3qyaOVn9FZzpPN4i3H1c4Z8aFfc5/q6WgXeTE2tvmbTgKE69vts7119IHohiO/FW9KbeJ+Jv2DZzP9VcsNfFF5SY/RCpN6eSnOTc9qU//3pegcETBJ2HjCu7wQlE40Dc6Ah0BElOumKduxzZHR7romLwXFgZZ1Q0TMdY6xtW/cFx13G1iwr6a+Ivmoov8OAzjiwaiwLhW2EkOVwhp59NHcQusqmdp5Ve9JcdIuKv+3lW/altFDS7ipa/yGH+4EKY6zh17vmNd6/Imfw0mm8IvkOBt8/fcIRDes7nOg4RRXAHupSPQEZgjwDPg3G04onBy83eQmIgxAfoexc6GYIeOQWe+M4SMCy5HYnzBMBSZVoQd7ce1L9qZv3sZ34kvMqYhU5FfDfUfLequ+ntX/aptHzW40E/DX2uo/3IxafbRpVoDz8Bfg8kZzK1W+bb5Ww1YawEPJwK5ykDS1UtHoCMQIeCSxddR1d6bnOB4xNI/WsYkKJz0dMwf53Gwxxij45uPRkB3qeUMGF8qFNl2lWLn2hfvs3BHkUegSa5eXkpbKvVLPlxluav+3lW/ahtZLS7qpzyC+VwLV7opPD0QnydOtVf0W4vJFbk6mPq2+Xsh8WGeEy4oY9PoMbmh4vTtbR6RmL4bxOFeOgJvJQIasPiAwez9sksAA9tkB49K8Czu1ZZLxvjSQa3BTrw8z8+kmsejk+9Rnsv/Gr/OZXON3rvq7131qybW8HZc5gi+bZi8bf7OI35sjfDmKRrOeU+VM90+XqqdmSUi5lYrL9Fe1IlxZmiv6Ah0BDihcjX/Rn25X8Do7aEj0BHoCHQEOgIdgY5AhIDyJF4n4D1YnoDj6Yf1O0MRb9/sCHQELhwBNxHizhCPyD1TB7/qR+UuHO5uXkegI9AR6Ah0BDoCdwyB/wOB2KMgbjOXxgAAAABJRU5ErkJggg==\n", "text/latex": [ "$\\displaystyle \\left\\{ \\pi : 3.14159265358979, \\ T_{eq} : 1.0, \\ \\bar{\\epsilon} : 0.01, \\ j : 6, \\ α : 0.9, \\ γ : 10, \\ δ : 0.02, \\ θ_{0} : 0.2, \\ κ : 1.8, \\ τ : 0.0003\\right\\}$" ], "text/plain": [ "{π: 3.14159265358979, T_eq: 1.0, \\bar{\\epsilon}: 0.01, j: 6, α: 0.9, γ: 10, δ:\n", " 0.02, θ₀: 0.2, κ: 1.8, τ: 0.0003}" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "discretize = ps.fd.Discretization2ndOrder(dx=0.03, dt=1e-5)\n", "parameters = {\n", " τ: 0.0003,\n", " κ: 1.8,\n", " εb: 0.01,\n", " δ: 0.02,\n", " γ: 10,\n", " j: 6,\n", " α: 0.9,\n", " Teq: 1.0,\n", " θzero: 0.2,\n", " sp.pi: sp.pi.evalf()\n", "}\n", "parameters" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "dφ_dt = - dF_dφ / τ\n", "φ_eqs = ps.simp.sympy_cse_on_assignment_list([ps.Assignment(φ_delta_field.center, \n", " discretize(dφ_dt.subs(parameters)))])\n", "φ_eqs.append(ps.Assignment(φ_tmp, discretize(ps.fd.transient(φ) - φ_delta_field.center)))\n", "\n", "temperature_evolution = -ps.fd.transient(T) + ps.fd.diffusion(T, 1) + κ * φ_delta_field.center\n", "temperature_eqs = [\n", " ps.Assignment(T, discretize(temperature_evolution.subs(parameters)))\n", "]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When creating the kernels we pass as target (which may be 'cpu' or 'gpu') the default target of the target handling. This enables to switch to a GPU simulation just by changing the parameter of the data handling.\n", "\n", "The rest is similar to the previous tutorial." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "φ_kernel = ps.create_kernel(φ_eqs, cpu_openmp=4, target=dh.default_target).compile()\n", "temperature_kernel = ps.create_kernel(temperature_eqs, target=dh.default_target).compile()" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "def timeloop(steps=200):\n", " φ_sync = dh.synchronization_function(['phi'])\n", " temperature_sync = dh.synchronization_function(['T'])\n", " dh.all_to_gpu() # this does nothing when running on CPU\n", " for t in range(steps):\n", " φ_sync()\n", " dh.run_kernel(φ_kernel)\n", " dh.swap('phi', 'phi_temp')\n", " temperature_sync()\n", " dh.run_kernel(temperature_kernel)\n", " dh.all_to_cpu()\n", " return dh.gather_array('phi')\n", " \n", "def init(nucleus_size=np.sqrt(5)):\n", " for b in dh.iterate():\n", " x, y = b.cell_index_arrays\n", " x, y = x-b.shape[0]//2, y-b.shape[0]//2\n", " bArr = (x**2 + y**2) < nucleus_size**2\n", " b['phi'].fill(0)\n", " b['phi'][(x**2 + y**2) < nucleus_size**2] = 1.0\n", " b['T'].fill(0.0)\n", " \n", "def plot():\n", " plt.subplot(1,3,1)\n", " plt.scalar_field(dh.gather_array('phi'))\n", " plt.title(\"φ\")\n", " plt.colorbar()\n", " plt.subplot(1,3,2)\n", " plt.title(\"T\")\n", " plt.scalar_field(dh.gather_array('T'))\n", " plt.colorbar()\n", " plt.subplot(1,3,3)\n", " plt.title(\"∂φ\")\n", " plt.scalar_field(dh.gather_array('phidelta'))\n", " plt.colorbar()" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Name| Inner (min/max)| WithGl (min/max)\n", "----------------------------------------------------\n", " T| ( 0, 0)| ( 0, 0)\n", " phi| ( 0, 1)| ( 0, 1)\n", "phi_temp| ( 0, 0)| ( 0, 0)\n", "phidelta| ( 0, 0)| ( 0, 0)\n", "\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7oAAAF1CAYAAADGJZYlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3df7DldX3n+efLVjD+yAhBsG1w6GQ6ybYzCUEWrc2YSYKEhppJYyruNMlojyHFUAO7cXdSlXatnXXWpYoxcTJxB2HRYW13XFkmxtCb6pEgNRl3dkRpEkTAtLRooKEXBkw0agKh+71/fL/33uP13B997z333Pv9PB9V37rn++uczzk0rzrv8/18Pt9UFZIkSZIkDcULpt0ASZIkSZLWkoWuJEmSJGlQLHQlSZIkSYNioStJkiRJGhQLXUmSJEnSoFjoSpIkSZIGxUJXkiRJkjQoFrqSpKYl+ebIciLJX4ys/+K02ydJkk6eha4kqWlV9bKZBXgU+Hsj2z467fZJ0mokOT/JI0keTfLfT7s90nqx0NWikvxXST6f5CeS/Ock9ybZOe12SZIkaVm+BJwH7AbeleSNSbYm+WSSn03yUJLHkvzclNsprSkLXS0oyUuAfwvcAPwi8FHgVuCjSTLNtkmSJGlpVfXNqvpGVf0R3Xe5S4APAP8Z+GvAc8DPA/uTbJteS6W1ZaGrxVwIBPgg8FLgz4D/FfhbwGum2C5JkiQtQ5Irk3w2yX7gOPBi4O8C1wGnAn9WVZ8FPg+8aXotldaWha4WcxbweFXVzIaq+kvgT4FXTa1VkiRJWlKSvwX8L8DlwP8G/DLdldwXAo/NO/wYfr/TgFjoajGPA2ePdlNO8j3AacDRqbVKkiRJy/Em4P+uqmPAPcC36IahPcd3987bht/vNCAWulrMZ+kC8Rq6LsxbgH8G/L9V9fg0GyZJkqQlfYPu+xvAPwU+XVV/Qlfs/o/AKQBJdgM/DBycRiOlSXjhtBugjauq/qoPvluA8+nGdfxH4K1TbZgkSZKW46PAm5M8ADwM/FK//Vfo5mB5L/A9wCuBn6uqP51KK6UJyMjwS2lBSf4NcKSq3j3ttkiSJGn1kvwy8A+q6ien3RZprXlFV8tSVf9g2m2QJEmSpOVwjK6kiUpyS5Kn+m5T4/YnyfuTHElyf5Lz17uNkrQWkuxKcrjPs31j9v9wks8keTbJry7n3CSnJ7kzycP939PW471I0kI2S9ZZ6EqatA8DuxbZfymwo1+uAm5chzZJ0ppKsgW4gS7TdgJXJNk577CvAf8t8Bsnce4+4K6q2gHc1a9La6KqPmS3ZZ2MzZR1FrqSJqqqPk0XeAvZDXykOncDr0iydX1aJ0lr5kK6uSweqarn6Ga13T16QFU9VVX3AH91EufuBvb3j/fT3Q9VkqZl02Sdha6kadvGd960/mi/TZI2k9Vk2WLnntXfA5X+75mrbKckrcamyboNMRnVGWecUeeee+60myFtavfee+/TVfXKkz3vkp96aT3zteMrf937n30Q+MuRTTdX1c0n8RQZs22Q08GbddLaWEnerUPWrSbLBpeD5p20embd6myIQvfcc8/l0KFD026GtKkl+ZOVnPf0147z2TvOXvHrvmjrl/+yqi5Y8RN0v+adM7J+NvDEKp5vwzLrpLWxkrxbh6xbTZYtdu6TSbZW1bF+WMdTy23zNJl30uqZdatj12VJ03YAeFs/+/IbgK/PdF2RpE3kHmBHku1JTgH20OXbas89AOztH+8Fbl/DNkvSydo0WbchruhKmqbieJ2Y2LMn+Rjwk8AZSY4C/xPwIoCqugk4CFwGHAG+Dbx9Yo2R1LDJZl1VPZ/kWuAOYAtwS1U9mOTqfv9NSV4FHAK+FziR5B3Azqr6xrhz+6e+HrgtyZXAo8BbJvYmJA2AWTfDQldqXAEnJjgUrKquWGJ/AddMrAGSxOSzDqCqDtL9eDe67aaRx/8fXVe9ZZ3bb38GuGhtWyppqMy6ORa6kjjB5H75k6SNwqyT1AKzruMYXUmSJEnSoHhFV2pcURyvTX0XC0laklknqQVm3RwLXUkTH8shSRuBWSepBWZdx0JXalwBxw1ESQNn1klqgVk3x0JXkr/8SWqCWSepBWZdx8moJEmSJEmD4hVdqXEFTlogafDMOkktMOvmWOhK8m5rkppg1klqgVnXsdCVGleUkxZIGjyzTlILzLo5FrpS6wqOm4eShs6sk9QCs26Wk1FJkiRJkgbFK7pS4wrHckgaPrNOUgvMujkWulLzwnEy7UZI0oSZdZJaYNbNsNCVGlfACcdySBo4s05SC8y6OY7RlSRJkiQNild0JdnFRVITzDpJLTDrOha6UuMKA1HS8Jl1klpg1s2x0JXEiTIQJQ2fWSepBWZdx0JXapy//ElqgVknqQVm3Rwno5IkSZIkDYpXdKXGFeG4v3lJGjizTlILzLo5FrqSHMshqQlmnaQWmHUdC12pcY7lkNQCs05SC8y6ORa6UvPC8bKLi6ShM+sktcCsm+GnIEmSJEkaFK/oSo0r4IS/eUkaOLNOUgvMujkWupIcyyGpCWadpBaYdR0LXalxVY7lkDR8Zp2kFph1c/wUJEmSJEmD4hVdSZywi4ukBph1klpg1nUsdKXGdfdbs3OHpGEz6yS1wKyb46cgNa8by7HSRZI2h8lnXZJdSQ4nOZJk35j9SfL+fv/9Sc7vt/9QkvtGlm8keUe/791JHh/Zd9mafiySBsasm+EVXalxTkMvqQWTzrokW4AbgIuBo8A9SQ5U1UMjh10K7OiX1wM3Aq+vqsPAeSPP8zjwiZHzfrOqfmNijZc0GGbdHL/dSpIkrd6FwJGqeqSqngNuBXbPO2Y38JHq3A28IsnWecdcBHy5qv5k8k2WpJO2abLOQlcSxysrXiRps1hl1p2R5NDIctW8p98GPDayfrTfdrLH7AE+Nm/btX33v1uSnLaiNy+pGWZdx67LUuOKOGmBpMFbg6x7uqouWGT/uF/+6mSOSXIK8LPAO0f23wi8pz/uPcD7gF9aToMltcesm2OhK4kTTiolqQETzrqjwDkj62cDT5zkMZcCf1hVT85sGH2c5IPA761VgyUNk1nX8dut1LiZaehXukjSZrAOWXcPsCPJ9v5qxR7gwLxjDgBv62ckfQPw9ao6NrL/CuZ15Zs3ru3NwAMn+dYlNcSsm+MVXUmSpFWqqueTXAvcAWwBbqmqB5Nc3e+/CTgIXAYcAb4NvH3m/CQvoZvF9B/Ne+r3JjmP7vvrV8fsl6R1s5myzkJXalzhpFKShm89sq6qDtJ9wRvddtPI4wKuWeDcbwPfN2b7W9e4mZIGzKybY6EryfvoSmqCWSepBWZdx0JXalwVHHcyKkkDZ9ZJaoFZN8dPQZIkSZI0KF7RlZoXToy93ZkkDYlZJ6kFZt0MC12pcYVdXCQNn1knqQVm3RwLXUneD1dSE8w6SS0w6zoWulLjinDC2wtJGjizTlILzLo5lvuSJEmSpEFZstBN8uIkn0vy+SQPJvln/fbTk9yZ5OH+72kj57wzyZEkh5NcMsk3IGn1jvOCFS9DYdZJw2fWmXVSC8y6znLezbPAT1fVjwLnAbuSvAHYB9xVVTuAu/p1kuwE9gCvBXYBH0iyZRKNl7R6BZyoF6x4GRCzThows26WWScNmFk3Z8l3U51v9qsv6pcCdgP7++37gcv7x7uBW6vq2ar6CnAEuHBNWy1pDYXjq1iGwqyThs6sA7NOGj6zbsayJqPqf7m7F/gbwA1V9dkkZ1XVMYCqOpbkzP7wbcDdI6cf7bfNf86rgKsAXsxLuPgFb1n5u5DEyzntdSs5b+aXP00+67acdhrbf+t9k3wLkhZg1s2ZRNb1z2veSWvolHPOPunvdmbdnGV9ClV1vKrOA84GLkzyNxc5fNxPATXmOW+uqguq6oIXceryWitJEzTprNvyspeuVVMlacUmkXX985p3kjaMk7q9UFX9WZI/oBuj8WSSrf2vfluBp/rDjgLnjJx2NvDEWjRW0mQMravKapl10jCZdd/JrJOGyazrLGfW5VcmeUX/+HuANwF/DBwA9vaH7QVu7x8fAPYkOTXJdmAH8Lm1briktVEVJy3ArJOGzqzrmHXSsJl1c5ZzRXcrsL8fz/EC4Laq+r0knwFuS3Il8CjwFoCqejDJbcBDwPPANVV1fDLNl7QWjk8w2JLsAn4L2AJ8qKqun7f/rwH/BngNXSb9RlX97xNr0MLMOmngJpl1m4hZJw2cWddZstCtqvuBHxuz/RngogXOuQ64btWtk7Sp9V+kbgAupuv+dk+SA1X10Mhh1wAPVdXfS/JK4HCSj1bVc+vZVrNOUgvMOkmtOKkxupKGp4ATkxvLcSFwpKoeAUhyK92tKkYL3QJeniTAy4Cv0V01kKQ1M+Gsk6QNwaybY6ErNS+T7OKyDXhsZP0o8Pp5x/wrujFgTwAvB/5+VZ2YVIMktWqiWSdJG4RZN8NCV2pcd7+1Vf3yd0aSQyPrN1fVzf3j5dyW4hLgPuCngR8A7kzy/1TVN1bTKEkatQZZJ0kbnlk3x0JXEseXd0vthTxdVRcssG85t6V4O3B9VRVwJMlXgB/GWT0lrbFVZp0kbQpmXcdPQdIk3QPsSLI9ySnAHrpuyqMepZ8AJclZwA8Bj6xrKyVJkjQoXtGVGldkYl1cqur5JNcCd9DdXuiW/lYVV/f7bwLeA3w4yRfoujr/WlU9PZEGSWrWJLNOkjYKs26Oha4kTkywc0dVHQQOztt208jjJ4CfmVgDJKk3yayTpI3CrOtY6EqNq4Lj/vInaeDMOkktMOvmWOhKsouLpCaYdZJaYNZ1vK4tSZIkSRoUr+hKjesmLfA3L0nDZtZJaoFZN8dCVxLHsYuLpOEz6yS1wKzrWOhKjSscyyFp+Mw6SS0w6+Z4XVuSJEmSNChe0ZWa51gOSS0w6yS1wKyb4acgiRNkxYskbRaTzroku5IcTnIkyb4x+5Pk/f3++5OcP7Lvq0m+kOS+JIdGtp+e5M4kD/d/T1uTD0PSYJl1HQtdqXEzNxZf6SJJm8Gksy7JFuAG4FJgJ3BFkp3zDrsU2NEvVwE3ztv/U1V1XlVdMLJtH3BXVe0A7urXJWkss26Oha4kTtQLVrxI0mYx4ay7EDhSVY9U1XPArcDuecfsBj5SnbuBVyTZusTz7gb294/3A5cv/x1LapFZ1/FbqiRJ0tLOSHJoZLlq3v5twGMj60f7bcs9poDfT3LvvOc+q6qOAfR/z1ztG5GkRQwm65yMSmpcd2NxuyBLGrY1yLqn53Wzm2/ck9dJHPPjVfVEkjOBO5P8cVV9eiUNldQus26OV3QlORmVpCZMOOuOAueMrJ8NPLHcY6pq5u9TwCfougcCPDnT5a//+9RJvm1JjTHrOha6UuNmbiy+0kWSNoN1yLp7gB1Jtic5BdgDHJh3zAHgbf2MpG8Avl5Vx5K8NMnLAZK8FPgZ4IGRc/b2j/cCt6/qg5A0aGbdHLsuS5IkrVJVPZ/kWuAOYAtwS1U9mOTqfv9NwEHgMuAI8G3g7f3pZwGfSALdd7P/s6o+2e+7HrgtyZXAo8Bb1uktSdJ32UxZZ6ErydmTJTVh0llXVQfpvuCNbrtp5HEB14w57xHgRxd4zmeAi9a2pZKGzKzrWOhKrbMLsqQWmHWSWmDWzbLQlRpX4KRSkgbPrJPUArNujoWuJH/5k9QEs05SC8y6jgPzJEmSJEmD4hVdqXEz09BL0pCZdZJaYNbNsdCVZCBKaoJZJ6kFZl3HQldqXOHsfJKGz6yT1AKzbo6FriRn55PUBLNOUgvMuo6TUUmSJEmSBsUrulLryrEckhpg1klqgVk3y0JXapyz80lqgVknqQVm3RwLXUkGoqQmmHWSWmDWdRyjK0mSJEkaFK/oSo1zGnpJLTDrJLXArJtjoSuJMhAlNcCsk9QCs65joSvJ+61JaoJZJ6kFZl3HQldqXDkNvaQGmHWSWmDWzXEyKkmSJEnSoHhFV5JjOSQ1wayT1AKzrmOhKzXP2fkktcCsk9QCs26Gha4kf/mT1ASzTlILzLqOha7UuMJJCyQNn1knqQVm3Rwno5IkSZIkDYpXdKXWVTcVvSQNmlknqQVm3SwLXUneWFxSE8w6SS0w6zoWulLjCictkDR8Zp2kFph1cxyjK0mSJEkaFK/oSs3zfmuSWmDWSWqBWTfDQleSkxZIaoJZJ6kFZl3HQleSYzkkNcGsk9QCs65joSs1rspAlDR8Zp2kFph1c5yMSpIkSZI0KBa6kjhRWfEiSZvFpLMuya4kh5McSbJvzP4keX+///4k5/fbz0ny75N8McmDSX5l5Jx3J3k8yX39ctmafSCSBsms69h1WZKTFkhqwiSzLskW4AbgYuAocE+SA1X10MhhlwI7+uX1wI393+eBf1JVf5jk5cC9Se4cOfc3q+o3Jtd6SUNi1nUsdCU5lkNSEyacdRcCR6rqEYAktwK7gdEvf7uBj1RVAXcneUWSrVV1DDjWtbH+PMkXgW3zzpWkZTHrOnZdlhpXhKqVL5K0GaxB1p2R5NDIctW8l9gGPDayfrTfdlLHJDkX+DHgsyObr+27/92S5LQVfwiSBs+sm2OhK0mStLSnq+qCkeXmefvH/fI3vwPhosckeRnwceAdVfWNfvONwA8A59FdCXnfilovScszmKxbstBdaNBwktOT3Jnk4f7vaSPnvLMffHw4ySWrbaSkyapVLENh1knDN+GsOwqcM7J+NvDEco9J8iK6L34frarfmW1z1ZNVdbyqTgAfpOs2uGJmnTR8Zl1nOVd0ZwYN/xfAG4BrkuwE9gF3VdUO4K5+nX7fHuC1wC7gA/2gZUkbUX+/Nbsum3XSoE0+6+4BdiTZnuQUunw4MO+YA8Db+hlJ3wB8vaqOJQnwr4EvVtW/GD0hydaR1TcDD6z0I+iZddKQmXWzlix0q+pYVf1h//jPgZlBw7uB/f1h+4HL+8e7gVur6tmq+gpwhDWoyCVNkJd0zTqpBRPMuqp6HrgWuIMuP26rqgeTXJ3k6v6wg8AjdHnxQeAf99t/HHgr8NNjbq3x3iRfSHI/8FPAf7fi949ZJzXBrANOctbleYOGz+pnzqKv0M/sD9sG3D1y2rgByvQDm68CeDEvOdl2S9okkuwCfgvYAnyoqq4fc8xPAv8SeBHd2JC/s66N/O72nMsEsm7Lac4hIw1ZVR2k+4I3uu2mkccFXDPmvP/I+DFtVNVb17iZs9Yy6/rnM++kBmyWrFt2oTt/0HB35Xn8oWO2fdfvA/3A5psBvjenD+i6kLT5TKoLcpZxr7UkrwA+AOyqqkdHvlxNxSSz7tTXnGPWSVM0sOEWq7LWWQfmnbRRmHWdZc26vMCg4Sdn+lL3f5/qty9ngLKkDaRq5csSZu+1VlXPATP3Whv1C8DvVNWjXVvqKabErJOGbYJZt6mYddKwmXWd5cy6vNCg4QPA3v7xXuD2ke17kpyaZDuwA/jc2jVZ0loqVj1pwWL3W1vOvdZ+EDgtyR8kuTfJ2yb6hhdg1knDtgZZNwhmnTRsZt2c5XRdnhk0/IUk9/Xb/gfgeuC2JFcCjwJvAegHI98GPEQ3s981VXV8zVsuaW0UsLpge7qqLlhg33K6vL0QeB1wEfA9wGeS3F1VX1pNo1bArJOGbPVZNxRmnTRkZt2sJQvdxQYN030xHXfOdcB1q2iXpGFY7r3Wnq6qbwHfSvJp4EeBdS10zTpJLTDrJLViWWN0JQ3bBMdyLOdea7cDb0zywiQvAV5PN129JK0px61JaoFZ1zmp2wtJGqgJBVtVPZ9k5l5rW4BbZu611u+/qaq+mOSTwP3ACbpbEK36JuGS9F0G9iVOksYy6wALXUlMdvKBpe611q//OvDrE2uEJE046yRpYzDrZljoSvKXP0ltMOsktcCsAxyjK0mSJEkaGK/oSq0r7OIiafjMOkktMOtmWehKsouLpDaYdZJaYNYBFrqSgIVvqShJQ2LWSWqBWQeO0ZUkSZIkDYxXdCXZxUVSG8w6SS0w6wALXUlgIEpqg1knqQVmHWChK6kAZ+eTNHRmnaQWmHWzLHQlUf7yJ6kBZp2kFph1HSejkiRJkiQNild0JTmWQ1IbzDpJLTDrAAtdSeBYDkltMOsktcCsAyx0JQHxlz9JDTDrJLXArOtY6EqtK+ziImn4zDpJLTDrZjkZlSRJkiRpULyiKzUvjuWQ1ACzTlILzLoZFrqS7OIiqQ1mnaQWmHWAha4kMBAltcGsk9QCsw5wjK4kSZIkaWC8oivJX/4ktcGsk9QCsw6w0JVUOGmBpOEz6yS1wKybZaEryRuLS2qCWSepBWZdxzG6kuZuLr6SRZI2iwlnXZJdSQ4nOZJk35j9SfL+fv/9Sc5f6twkpye5M8nD/d/TVvbmJTXDrAMsdLVMdzxx3+wiSZK+U5ItwA3ApcBO4IokO+cddimwo1+uAm5cxrn7gLuqagdwV78uSVOxmbLOQleSJGn1LgSOVNUjVfUccCuwe94xu4GPVOdu4BVJti5x7m5gf/94P3D5pN+IJC1i02Sdha4kUitfJGmzmHDWbQMeG1k/2m9bzjGLnXtWVR0D6P+eudz3K6lNZl3Hyai0LJe8+rxpN0GT5Ox8klqwuqw7I8mhkfWbq+rmkfVxTz7/a+NCxyznXElaHrMOsNCV5KRSklqw+qx7uqouWGT/UeCckfWzgSeWecwpi5z7ZJKtVXWs7/r31EoaL6kRZt0suy5LkiSt3j3AjiTbk5wC7AEOzDvmAPC2fkbSNwBf77voLXbuAWBv/3gvcPuk34gkLWLTZJ1XdCV5RVdSGyaYdVX1fJJrgTuALcAtVfVgkqv7/TcBB4HLgCPAt4G3L3Zu/9TXA7cluRJ4FHjL5N6FpEEw6wALXUk4qZSkNkw666rqIN0XvNFtN408LuCa5Z7bb38GuGhtWyppyMy6joWuJK/oSmqDWSepBWYd4BhdSZIkSdLAeEVXkr/8SWqDWSepBWYdYKErNe8kbhAuSZuWWSepBWbdHAtdSau9sbgkbQ5mnaQWmHWAha4ksIuLpDaYdZJaYNYBTkYlSZIkSRoYr+hKciyHpCaYdZJaYNZ1LHQl2cVFUhvMOkktMOsAC11Jzs4nqQVmnaQWmHWzHKMrSZIkSRoUr+hKsouLpDaYdZJaYNYBFrqSwECU1AazTlILzDrAQlcSjuWQ1AazTlILzLqOY3QlSZIkSYNioStJkiRJGhS7LktyLIekNph1klpg1gEWupK835qkFph1klpg1s2y0JXkL3+S2mDWSWqBWQdY6EoCA1FSG8w6SS0w6wAno5IkSZIkDYxXdKXGBcdySBo+s05SC8y6ORa6kuziIqkNZp2kFph1gIWuJGfnk9QCs05SC8y6WY7RlTRRSXYlOZzkSJJ9ixz3XyY5nuTn17N9kiRJGh6v6EqaWBeXJFuAG4CLgaPAPUkOVNVDY47758Adk2mJJGF3PkltMOuAZVzRTXJLkqeSPDCy7fQkdyZ5uP972si+d/ZXbg4nuWRSDZe0hmoVy+IuBI5U1SNV9RxwK7B7zHH/DfBx4KlVvpMVM+ukBkwu6zYV804aOLMOWF7X5Q8Du+Zt2wfcVVU7gLv6dZLsBPYAr+3P+UB/pUbSBpZa+QKckeTQyHLVyFNvAx4bWT/ab5t77WQb8Gbgpgm/zaV8GLNOGrRVZt2QfBjzThoss66zZKFbVZ8GvjZv825gf/94P3D5yPZbq+rZqvoKcITuio6kjWx1v/w9XVUXjCw3jzxzFni1Uf8S+LWqOr5m72cFzDqpAV7lAMw7afDMOmDlY3TPqqpjAFV1LMmZ/fZtwN0jx33X1ZsZ/VWfqwBezEtW2AxJG9xR4JyR9bOBJ+YdcwFwaxKAM4DLkjxfVb+7Pk1c1Jpm3ZbTTht3iCRtBOadpEFZ68molnP1ptvYXfW5GeB7c/rAfj+QNpHJ/oJ3D7AjyXbgcbrub7/wHS9ftX3mcZIPA7+3QYrcxawo6059zTlmnTQtA7xasU7MO2kzMetmrbTQfTLJ1v4Xv63MTSCznKs3kjaYSY3JqKrnk1xLN5vyFuCWqnowydX9/mmPy12KWScNyNDGn60x804aCLOus9L76B4A9vaP9wK3j2zfk+TU/grODuBzq2uipImb4FiOqjpYVT9YVT9QVdf1224aV+RW1T+sqt9ei7e0Rsw6aUimNG5tsRmN5x039r7jSX49yR8nuT/JJ5K8ot9+bpK/SHJfv6zmx0PzThoKsw5Y3u2FPgZ8BvihJEeTXAlcD1yc5GG6+2NeD1BVDwK3AQ8BnwSumfYEM5KW5ux8Zp3Ugilm3dgZjb+jbXP3Hb8U2Alc0c94DHAn8Der6keALwHvHDn1y1V1Xr9cvZzGmHfSsJl1nSW7LlfVFQvsumiB468DrlvOi0vSRmHWSZqg3cBP9o/3A38A/Nq8Y2bvOw6QZOa+4w9V1e+PHHc38POraYx5J2lCNlTWrbTrsqQhcRp6SS1YXdYtds/wpXzHjMbAmWOOWfK+471fAv7dyPr2JH+U5D8keeNJtEnSUJl1wNrPuixps7FgldSC1Wfd01V1wUI7k3wKeNWYXe9a5vMvObtxkncBzwMf7TcdA15TVc8keR3wu0leW1XfWOZrShoas26Wha7UuDA+cSRpSCaddVX1pgVfO1loRuNRi85unGQv8HeBi6qq+td8Fni2f3xvki8DPwgcWu37kbQ5mXVz7LosSZI0WQvNaDxq9r7jSU6hu+/4AehmKKUb5/azVfXtmROSvLKf2IUk3083I/IjE3sXkrS4DZV1FrqSHKMrqQ3Ty7qxMxoneXWSgwBV9Twwc9/xLwK39TMeA/wr4OXAnfNurfETwP1JPg/8NnB1VX1t1a2VtLmZdYBdlyUxrNsESdJCppV1VfUMY2Y0rqongMtG1g8CB8cc9zcWeN6PAx9fu5ZKGgKzrmOhK8krs5LaYNZJaoFZB1joSgIDUVIbzDpJLTDrAMfoSpIkSZIGxiu6UuvKMbqSGmDWSWqBWTfLQleSXVwktcGsk9QCsw6w0JWEv/xJaoNZJ6kFZl3HQleSv/xJaoNZJ6kFZh3gZFSSJEmSpIHxiq4kuzUaUbAAABHSSURBVLhIaoJZJ6kFZl3HQldqXWEXF0nDZ9ZJaoFZN8tCV5KBKKkNZp2kFph1gGN0JUmSJEkD4xVdqXHBsRyShs+sk9QCs26Oha4ku7hIaoNZJ6kFZh1goSsJSJmIkobPrJPUArOuY6Ertc7Z+SS1wKyT1AKzbpaTUUmSJEmSBsUrupKctEBSE8w6SS0w6zoWupLs4iKpDWadpBaYdYCFriT85U9SG8w6SS0w6zoWupL85U9SG8w6SS0w6wAno5IkSZIkDYxXdKXWlV1cJDXArJPUArNuloWuJLu4SGqDWSepBWYdYKErNS/4y5+k4TPrJLXArJvjGF1JkiRJ0qB4RVcSlD/9SWqAWSepBWYdYKErCbu4SGqDWSepBWZdx0JXal3hpAWShs+sk9QCs26Wha4kcmLaLZCkyTPrJLXArOs4GZUkSZIkaVC8oivJLi6S2mDWSWqBWQd4RVcS3aQFK10kabOYVtYlOT3JnUke7v+etsBxu5IcTnIkyb6R7e9O8niS+/rlspF97+yPP5zkktW1VNIQmHUdC12pdUU3Df1KF0naDKabdfuAu6pqB3BXv/4dkmwBbgAuBXYCVyTZOXLIb1bVef1ysD9nJ7AHeC2wC/hA/zySWmXWzbLQleQVXUlNmGLW7Qb294/3A5ePOeZC4EhVPVJVzwG39uct9by3VtWzVfUV4Ej/PJIaZtZ1LHQlSZKWdkaSQyPLVSdx7llVdQyg/3vmmGO2AY+NrB/tt824Nsn9SW4Z6Q641DmSdLIGk3VORiXJSQsktWF1Wfd0VV2w0M4knwJeNWbXu5b5/BmzbabFNwLv6dffA7wP+KUlzpHUKrMOsNCVmhfsgixp+CaddVX1pgVfO3kyydaqOpZkK/DUmMOOAueMrJ8NPNE/95Mjz/VB4PeWOkdSm8y6OXZdllq3mgkLnIxK0mYx3aw7AOztH+8Fbh9zzD3AjiTbk5xCN/HKAYD+C+OMNwMPjDzvniSnJtkO7AA+t9rGStrEzLpZXtGVJEmarOuB25JcCTwKvAUgyauBD1XVZVX1fJJrgTuALcAtVfVgf/57k5xH11Xvq8A/AqiqB5PcBjwEPA9cU1XH1/F9SdKoDZV1FrqS7LosqQnTyrqqega4aMz2J4DLRtYPAgfHHPfWRZ77OuC6tWmppCEw6zoWupKcukRSG8w6SS0w6wALXUl4RVdSG8w6SS0w6zoWulLrCjhhIkoaOLNOUgvMulnOuixJkiRJGhSv6EpyLIekNph1klpg1gEWupJwLIekNph1klpg1nUsdCWtxQ3CJWnjM+sktcCsAxyjK4nul7+VLks+d7IryeEkR5LsG7P/F5Pc3y//KcmPTuI9StIks06SNgqzrmOhK2likmwBbgAuBXYCVyTZOe+wrwB/p6p+BHgPcPP6tlKSJElDY9dlqXXFJCctuBA4UlWPACS5FdgNPDT78lX/aeT4u4GzJ9YaSe2abNZJ0sZg1s2y0JUaFyCTG8uxDXhsZP0o8PpFjr8S+HeTaoykdk046yRpQzDr5ljoSoITqzr7jCSHRtZvrqqZ7scZc/zY9E3yU3SF7t9eVWskaSGryzpJ2hzMOsBCV9LqPV1VFyyw7yhwzsj62cAT8w9K8iPAh4BLq+qZtW+iJEmSWmKhK2mSXVzuAXYk2Q48DuwBfuE7Xjt5DfA7wFur6kuTaogk2Z1PUgvMus7EZl1e6pYikjaIWuWy2FNXPQ9cC9wBfBG4raoeTHJ1kqv7w/4p8H3AB5LcN68b9IZn1kmbxASzrgVmnbRJmHWzJnJFd+SWIhfTdV28J8mBqnpo8TMlrb+a6I3Fq+ogcHDetptGHv8y8MsTa8AEmXXSZjLZrBsys07aTMy6GZO6ojt7S5Gqeg6YuaWIpA3IG4uvmFknbSJm3YqZddImYtZ1JlXojrulyLbRA5JcleRQkkN/xbMTaoYkTdRJZd3xb35rXRsnSWtkyawD807SxjKpyaiWvKVIf/uRmwEuuOCCuvPQv51QU6Q2JLl3xSfbxWWlTjrrDv3KP1mPdkmDlnf86spONOtWalm3ijPvpLWVd/zqyr7bmXXA5ArdZd1SRNIGUBDvt7ZSZp20WZh1q2HWSZuFWTdrUl2XZ28pkuQUuluKHJjQa0laraqVL20z66TNxKxbKbNO2kzMOmBCV3Sr6vkkM7cU2QLcUlUPTuK1JGlazDpJLTDrJG1Gk+q6PPaWIpI2qGH9gLeuzDppEzHrVsyskzYRsw6YYKErafPIwLqqSNI4Zp2kFph1HQtdSYMbkyFJY5l1klpg1gEWupIKcHY+SUNn1klqgVk3a1KzLkuSJEmSNBVe0ZUaF8qxHJIGz6yT1AKzbo6FriTHckhqg1knqQVmHWChKwkMREltMOsktcCsAyx0JTlpgaQWmHWSWmDWzXIyKkmSJEnSoFjoSiJVK14kabOYVtYlOT3JnUke7v+etsBxu5IcTnIkyb6R7f9Xkvv65atJ7uu3n5vkL0b23bSqhkoaBLOuY9dlSY7lkNSG6WXdPuCuqrq+/1K3D/i10QOSbAFuAC4GjgL3JDlQVQ9V1d8fOe59wNdHTv1yVZ038XcgafMw6wCv6EqiukBc6SJJm8JUs243sL9/vB+4fMwxFwJHquqRqnoOuLU/b1aSAP818LHVNkjSUJl1Myx0JUmSlnZGkkMjy1Unce5ZVXUMoP975phjtgGPjawf7beNeiPwZFU9PLJte5I/SvIfkrzxJNokSeMMJuvsuiy1rvDKrKThW33WPV1VFyy0M8mngFeN2fWuZT5/xmyb3+Ar+M4rHMeA11TVM0leB/xuktdW1TeW+ZqShsasm2WhK8lp6CW1YYJZV1VvWmhfkieTbK2qY0m2Ak+NOewocM7I+tnAEyPP8ULg54DXjbzms8Cz/eN7k3wZ+EHg0Grei6RNzqwD7LosCWddltSGKWbdAWBv/3gvcPuYY+4BdiTZnuQUYE9/3ow3AX9cVUdn30/yyn5iF5J8P7ADeGS1jZW0uZl1Ha/oSrLrsqQ2TC/rrgduS3Il8CjwFoAkrwY+VFWXVdXzSa4F7gC2ALdU1YMjz7GH756Y5SeA/znJ88Bx4Oqq+tqE34ukjc6sAyx0JUmSJqqqngEuGrP9CeCykfWDwMEFnuMfjtn2ceDja9ZQSVqFjZZ1FrpS6wo44RVdSQNn1klqgVk3y0JXap73w5XUArNOUgvMuhkWupIMREltMOsktcCsAyx0JYGBKKkNZp2kFph1gLcXkiRJkiQNjFd0pdY5aYGkFph1klpg1s2y0JWaV1Anpt0ISZows05SC8y6GRa6khzLIakNZp2kFph1gGN0JUmSJEkD4xVdqXWO5ZDUArNOUgvMulkWupLs4iKpDWadpBaYdYCFriQwECW1wayT1AKzDrDQlUQZiJIaYNZJaoFZN8PJqCRJkiRJg+IVXal1BZzwfmuSBs6sk9QCs26Wha4ku7hIaoNZJ6kFZh1goSsJDERJbTDrJLXArAMsdCVR3m9NUgPMOkktMOtmOBmVJEmSJGlQvKIrta6gykkLJA2cWSepBWbdLAtdSXZxkdQGs05SC8w6wEJXEjhpgaQ2mHWSWmDWAY7RlSRJkiQNjFd0pdZVeWNxScNn1klqgVk3y0JXkl1cJLXBrJPUArMOsNCVBJS//ElqgFknqQVmXcdCV2pe+cufpAaYdZJaYNbNcDIqSZIkSdKgeEVXal3h/dYkDZ9ZJ6kFZt0sC11JUI7lkNQAs05SC8w6wEJXal4B5S9/kgbOrJPUArNujmN0pdZVdb/8rXRZQpJdSQ4nOZJk35j9SfL+fv/9Sc6fyPuU1LYJZ91ikpye5M4kD/d/T1vguFuSPJXkgeWen+SdfX4eTnLJqhoqafMz62ZZ6EqamCRbgBuAS4GdwBVJds477FJgR79cBdy4ro2UpMnbB9xVVTuAu/r1cT4M7Fru+X2e7gFe25/3gT53JWkaNlTWWehKok7UipclXAgcqapHquo54FZg97xjdgMfqc7dwCuSbF37dympdRPMuqXsBvb3j/cDl49tX9Wnga+dxPm7gVur6tmq+gpwhC53JTXMrOs4RlfSJCct2AY8NrJ+FHj9Mo7ZBhybVKMkNWp6E7ScVVXHAKrqWJIz1+j8bcDdI8fN5Keklpl1wAYpdO+9996nk3wLeHrabemdgW0Zx7aMt1Ha8tdXctKf86d3fKp++4xVvO6LkxwaWb+5qm7uH2fM8fN/LlzOMYNg1i3KtoxnW8Y76bybcNaR5FPAq8ac965VvOZSNmx+3nvvvd9Mcnja7ehtpH+7tmU82zKeWde/7JhtS2bdhih0q+qVSQ5V1QXTbguAbRnPtoy3kdqyElU1bozEWjkKnDOyfjbwxAqOGQSzbmG2ZTzbsnYmnHVU1ZsW2pfkySRb+ysUW4GnTvLpFzp/I+fn4Y3y72Uj/du1LePZlrVj1s1xjK6kSboH2JFke5JT6CYSODDvmAPA2/rZl98AfH2m24okDcQBYG//eC9w+xqdfwDYk+TUJNvpJvX73CrbKkkrtaGyzkJX0sRU1fPAtcAdwBeB26rqwSRXJ7m6P+wg8AjdxAIfBP7xVBorSZNzPXBxkoeBi/t1krw6ycGZg5J8DPgM8ENJjia5crHzq+pB4DbgIeCTwDVVdXyd3pMkzbehsm5DdF3u3bz0IevGtoxnW8bbSG3ZcKrqIF0xO7rtppHHBVyz3u2aoo3078W2jGdbxttIbdlUquoZ4KIx258ALhtZv+Jkzu/3XQdctzYtXVMb6d+LbRnPtoy3kdqyqWy0rEv3HVOSJEmSpGGw67IkSZIkaVCmXugm2ZXkcJIjSfZN4fW/muQLSe6bmUo7yelJ7kzycP/3tAm99i1JnkrywMi2BV87yTv7z+lwkkvWoS3vTvJ4/9ncl+SykX2TbMs5Sf59ki8meTDJr/Tb1/2zWaQtU/lstHmZdWbdmLaYdRocs86sG9MWs07TUVVTW4AtwJeB7wdOAT4P7FznNnwVOGPetvcC+/rH+4B/PqHX/gngfOCBpV4b2Nl/PqcC2/vPbcuE2/Ju4FfHHDvptmwFzu8fvxz4Uv+a6/7ZLNKWqXw2LptzMevMugXaYta5DGox68y6Bdpi1rlMZZn2Fd0LgSNV9UhVPQfcCuyecpuga8P+/vF+4PJJvEhVfRr42jJfezdwa1U9W1VfoZuh9sIJt2Uhk27Lsar6w/7xn9PN1ruNKXw2i7RlIRP9bLRpmXVm3bi2mHUaGrPOrBvXFrNOUzHtQncb8NjI+lEW/8c2CQX8fpJ7k1zVbzur+vt49n/PXMf2LPTa0/qsrk1yf98FZqZLybq1Jcm5wI8Bn2XKn828tsCUPxttKhvh34VZtzizbnxbwKzT8m2Efxdm3eLMuvFtAbNucKZd6GbMtvWeBvrHq+p84FLgmiQ/sc6vv1zT+KxuBH4AOA84BrxvPduS5GXAx4F3VNU3Fjt00u0Z05apfjbadDbCvwuzbmFm3cJtMet0MjbCvwuzbmFm3cJtMesGaNqF7lHgnJH1s4En1rMB1d3Xiap6CvgEXXeEJ5NsBej/PrWOTVrotdf9s6qqJ6vqeFWdAD7IXFeNibclyYvoAuijVfU7/eapfDbj2jLNz0ab0tT/XZh1CzPrFm6LWaeTNPV/F2bdwsy6hdti1g3TtAvde4AdSbYnOQXYAxxYrxdP8tIkL595DPwM8EDfhr39YXuB29erTYu89gFgT5JTk2wHdgCfm2RDZsKn92a6z2bibUkS4F8DX6yqfzGya90/m4XaMq3PRpuWWffdzDqzTsNj1n03s86s07SczMxVk1iAy+hmPPsy8K51fu3vp5tJ7fPAgzOvD3wfcBfwcP/39Am9/sfoukf8Fd0vRlcu9trAu/rP6TBw6Tq05f8AvgDcT/c/+tZ1asvfpusWcj9wX79cNo3PZpG2TOWzcdm8i1ln1o1pi1nnMrjFrDPrxrTFrHOZypL+P6AkSZIkSYMw7a7LkiRJkiStKQtdSZIkSdKgWOhKkiRJkgbFQleSJEmSNCgWupIkSZKkQbHQlSRJkiQNioWuJEmSJGlQLHQlSZIkSYPy/wPQr4c+TRaXwwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "timeloop(10)\n", "init()\n", "plot()\n", "print(dh)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "result = None\n", "if 'is_test_run' not in globals():\n", " ani = ps.plot.scalar_field_animation(timeloop, rescale=True, frames=600)\n", " result = ps.jupyter.display_as_html_video(ani)\n", "result" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "assert np.isfinite(dh.max('phi'))\n", "assert np.isfinite(dh.max('T'))" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.7.3" } }, "nbformat": 4, "nbformat_minor": 2 }