{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6aa46c82",
   "metadata": {},
   "source": [
    "# 8. 2D Ising model\n",
    "\n",
    "For a system of interacting spins on a 2D lattice of sites, the Hamiltonian describing the energy of interaction is\n",
    "$$\n",
    "  {\\cal H} = - \\sum_{i,j} J_{ij} s_i s_j -\\mu \\sum_j h_j s_j\n",
    "$$\n",
    "where the index $j$ runs over all spins on the lattice, while the index $i$ runs over those spins with which spin $s_j$ interacts, with the interaction strength $J_{ij}$. In the case of spin-1/2, each $s_j$ can take on values of +1 or -1.  Typical example is a nearest-neighbour interaction, on a square lattice implying interactions with spins directly North, South, East and West of site $j$. The sign of the interaction $J_{ij}$ determines whether the interaction favors (lowers the energy for) spins aligned parallel to their neighbours ($J_{ij}>0$, ferromagnetic), or anti-parallel to their neighbours ($J_{ij}<0$, antiferromagnetic). In the absence of spin-spin interactions, $J_{ij}=0$.  \n",
    "\n",
    "The second sum describes the interaction of spins with an external field: $\\mu$ is the unit of strength (a magnetic moment of a single spin) and $h_j$ represents the strength pf the field at site $j$.  Positive $h_j$ value implies that the alignment along the external field is energetically favourable, $h_j < 0$ favours spins aligning against the field, and $h_j=0$ indicates the absence of the external field.\n",
    "\n",
    "The probability of encountering a particular configuration $\\{s\\}$ of spins on the lattice is given by the Boltzmann factor:\n",
    "$$\n",
    "  P(\\{s\\}) = \\frac{e^{-{\\cal H(\\{s\\})}/k_B T}}{\\sum_{\\{s\\}} e^{-{\\cal H(\\{s\\})}/k_B T}}\n",
    "$$\n",
    "with the sum over all possible configurations of the system in the denominator being the partition function $Z$.  Often the notation $\\beta = 1/k_b T$ (the inverse temperature) is used; with that,\n",
    "$$\n",
    "  P(\\{s\\}) = \\frac{1}{Z} e^{-\\beta{\\cal H(\\{s\\})}}\n",
    "$$\n",
    "In two or more dimensions, the Ising model exhibits interesting phase behaviour, including ferromagnetic phase transitions, and is widely studied.\n",
    "\n",
    "The limitations of finite size of the model lattice of spins are normally overcome by imposing periodic boundary conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "cc4d8af6",
   "metadata": {},
   "outputs": [],
   "source": [
    "%plot --format png -w 900 -h 600"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "76327018",
   "metadata": {},
   "outputs": [],
   "source": [
    "### function to display the spin state of the system, on a rectangular lattice of spins\n",
    "function DisplayIsing (S)\n",
    "  [n,m] = size(S);\n",
    "  N = round((n-1)/2);\n",
    "  [X,Y] = meshgrid(-N:N,-N:N);\n",
    "  quiver(X,Y-0.5*S,0*S,S,0,'color','blue');\n",
    "  axis equal;\n",
    "  xlim([-N-1,N+1]);\n",
    "  ylim([-N-1,N+1]);\n",
    "  axis off;\n",
    "endfunction;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "988b61a4",
   "metadata": {},
   "source": [
    "## Set up a 2D rectangular lattice of spins\n",
    "\n",
    "Initial state will be a random orientation of $2N+1 \\times 2N+1$ spins. Here we assume that all sites can only have the magnetic moment of $+\\mu$ or $-\\mu$. `rand()` produces a random number somewhere in the 0..1 range, and and a shift of -0.5 makes it equally likely for the sign to be up or down."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "487d0326",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4QAAAJYCAIAAAC1p7+MAAAJMmlDQ1BkZWZhdWx0X3JnYi5pY2MAAEiJlZVnUJNZF8fv8zzphUASQodQQ5EqJYCUEFoo0quoQOidUEVsiLgCK4qINEWQRQEXXJUia0UUC4uCAhZ0gywCyrpxFVFBWXDfGZ33HT+8/5l7z2/+c+bec8/5cAEgiINlwct7YlK6wNvJjhkYFMwE3yiMn5bC8fR0A9/VuxEArcR7ut/P+a4IEZFp/OW4uLxy+SmCdACg7GXWzEpPWeGjy0wPj//CZ1dYsFzgMt9Y4eh/eexLzr8s+pLj681dfhUKABwp+hsO/4b/c++KVDiC9NioyGymT3JUelaYIJKZttIJHpfL9BQkR8UmRH5T8P+V/B2lR2anr0RucsomQWx0TDrzfw41MjA0BF9n8cbrS48hRv9/z2dFX73kegDYcwAg+7564ZUAdO4CQPrRV09tua+UfAA67vAzBJn/eqiVDQ0IgALoQAYoAlWgCXSBETADlsAWOAAX4AF8QRDYAPggBiQCAcgCuWAHKABFYB84CKpALWgATaAVnAad4Dy4Aq6D2+AuGAaPgRBMgpdABN6BBQiCsBAZokEykBKkDulARhAbsoYcIDfIGwqCQqFoKAnKgHKhnVARVApVQXVQE/QLdA66At2EBqGH0Dg0A/0NfYQRmATTYQVYA9aH2TAHdoV94fVwNJwK58D58F64Aq6HT8Id8BX4NjwMC+GX8BwCECLCQJQRXYSNcBEPJBiJQgTIVqQQKUfqkVakG+lD7iFCZBb5gMKgaCgmShdliXJG+aH4qFTUVlQxqgp1AtWB6kXdQ42jRKjPaDJaHq2DtkDz0IHoaHQWugBdjm5Et6OvoYfRk+h3GAyGgWFhzDDOmCBMHGYzphhzGNOGuYwZxExg5rBYrAxWB2uF9cCGYdOxBdhK7EnsJewQdhL7HkfEKeGMcI64YFwSLg9XjmvGXcQN4aZwC3hxvDreAu+Bj8BvwpfgG/Dd+Dv4SfwCQYLAIlgRfAlxhB2ECkIr4RphjPCGSCSqEM2JXsRY4nZiBfEU8QZxnPiBRCVpk7ikEFIGaS/pOOky6SHpDZlM1iDbkoPJ6eS95CbyVfJT8nsxmpieGE8sQmybWLVYh9iQ2CsKnqJO4VA2UHIo5ZQzlDuUWXG8uIY4VzxMfKt4tfg58VHxOQmahKGEh0SiRLFEs8RNiWkqlqpBdaBGUPOpx6hXqRM0hKZK49L4tJ20Bto12iQdQ2fRefQ4ehH9Z/oAXSRJlTSW9JfMlqyWvCApZCAMDQaPkcAoYZxmjDA+SilIcaQipfZItUoNSc1Ly0nbSkdKF0q3SQ9Lf5RhyjjIxMvsl+mUeSKLktWW9ZLNkj0ie012Vo4uZynHlyuUOy33SB6W15b3lt8sf0y+X35OQVHBSSFFoVLhqsKsIkPRVjFOsUzxouKMEk3JWilWqUzpktILpiSTw0xgVjB7mSJleWVn5QzlOuUB5QUVloqfSp5Km8oTVYIqWzVKtUy1R1WkpqTmrpar1qL2SB2vzlaPUT+k3qc+r8HSCNDYrdGpMc2SZvFYOawW1pgmWdNGM1WzXvO+FkaLrRWvdVjrrjasbaIdo12tfUcH1jHVidU5rDO4Cr3KfFXSqvpVo7okXY5upm6L7rgeQ89NL0+vU++Vvpp+sP5+/T79zwYmBgkGDQaPDamGLoZ5ht2GfxtpG/GNqo3uryavdly9bXXX6tfGOsaRxkeMH5jQTNxNdpv0mHwyNTMVmLaazpipmYWa1ZiNsulsT3Yx+4Y52tzOfJv5efMPFqYW6RanLf6y1LWMt2y2nF7DWhO5pmHNhJWKVZhVnZXQmmkdan3UWmijbBNmU2/zzFbVNsK20XaKo8WJ45zkvLIzsBPYtdvNcy24W7iX7RF7J/tC+wEHqoOfQ5XDU0cVx2jHFkeRk4nTZqfLzmhnV+f9zqM8BR6f18QTuZi5bHHpdSW5+rhWuT5z03YTuHW7w+4u7gfcx9aqr01a2+kBPHgeBzyeeLI8Uz1/9cJ4eXpVez33NvTO9e7zofls9Gn2eedr51vi+9hP0y/Dr8ef4h/i3+Q/H2AfUBogDNQP3BJ4O0g2KDaoKxgb7B/cGDy3zmHdwXWTISYhBSEj61nrs9ff3CC7IWHDhY2UjWEbz4SiQwNCm0MXwzzC6sPmwnnhNeEiPpd/iP8ywjaiLGIm0iqyNHIqyiqqNGo62ir6QPRMjE1MecxsLDe2KvZ1nHNcbdx8vEf88filhICEtkRcYmjiuSRqUnxSb7JicnbyYIpOSkGKMNUi9WCqSOAqaEyD0tandaXTlz/F/gzNjF0Z45nWmdWZ77P8s85kS2QnZfdv0t60Z9NUjmPOT5tRm/mbe3KVc3fkjm/hbKnbCm0N39qzTXVb/rbJ7U7bT+wg7Ijf8VueQV5p3tudATu78xXyt+dP7HLa1VIgViAoGN1tubv2B9QPsT8M7Fm9p3LP58KIwltFBkXlRYvF/OJbPxr+WPHj0t6ovQMlpiVH9mH2Je0b2W+z/0SpRGlO6cQB9wMdZcyywrK3BzcevFluXF57iHAo45Cwwq2iq1Ktcl/lYlVM1XC1XXVbjXzNnpr5wxGHh47YHmmtVagtqv14NPbogzqnuo56jfryY5hjmceeN/g39P3E/qmpUbaxqPHT8aTjwhPeJ3qbzJqamuWbS1rgloyWmZMhJ+/+bP9zV6tua10bo63oFDiVcerFL6G/jJx2Pd1zhn2m9az62Zp2WnthB9SxqUPUGdMp7ArqGjzncq6n27K7/Ve9X4+fVz5ffUHyQslFwsX8i0uXci7NXU65PHsl+spEz8aex1cDr97v9eoduOZ67cZ1x+tX+zh9l25Y3Th/0+LmuVvsW523TW939Jv0t/9m8lv7gOlAxx2zO113ze92D64ZvDhkM3Tlnv296/d5928Prx0eHPEbeTAaMip8EPFg+mHCw9ePMh8tPN4+hh4rfCL+pPyp/NP637V+bxOaCi+M24/3P/N59niCP/Hyj7Q/Fifzn5Ofl08pTTVNG02fn3Gcufti3YvJlykvF2YL/pT4s+aV5quzf9n+1S8KFE2+Frxe+rv4jcyb42+N3/bMec49fZf4bmG+8L3M+xMf2B/6PgZ8nFrIWsQuVnzS+tT92fXz2FLi0tI/QiyQvpTNDAsAAAAJcEhZcwAACxMAAAsTAQCanBgAAAAddEVYdFNvZnR3YXJlAEdQTCBHaG9zdHNjcmlwdCA5LjI1wZk/DQAAES9JREFUeJzt3cFynbq2QNHN+/9/5jWccqUSOznOBU9pMUYr5dZeBUIT0chxnucLAAAK/1f/AAAAnkuMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKKs7jvoX3OM4Jo82krn2YonBLsQoG/Dk3YiLtSNXbSMuFvOIUeB68/bLt4mmzgUQEqMsbepOeRyv8/zxD+BybyvrPMcusalz8UxilA0M3lGGeY9s9mKJ7cISYyQxClxvWNy8F8DUuQBCYpR1Td0pf55rWNzACp6wxKbOxTOJUfbgybu+qS8PD2GJrc8SYyoxCtxiTNz8UgBT5wKoiFEWNXWn/H2uMXEDK3jOEps6Fw8kRtmGJ+/Kpr48PIoltjJLjMHEKHCXAXHzYQFMnQsgIUZZ0eCd8sO5Bgw7YARmeNQSGzAXvMQoa/rsCevJu53dL9nUW3HqXMCOxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkDnO86x/AwAAD+VkFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGWd1x1L/gNlNHM9dezLWdwaPxTGKUDYx88o4cik1NvRunzgXDiFEAADJilKVNPdg4jtd5/vgHhKbegVPngpHEKBs4T1sL3GjqEps6FwwjRgEAyIhR1vX+LXuYn+dyckPoCUsMWJ8YZQ+iDW41dYlNnQsmEaMAAGTEKIua+qHt97mc3JB4zhIDFidG2YZog1tNXWJT54IxxCgAABkxyqI+/NA24Ovb1LnYztRbcepcMJgYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhS4wHF87e+7mDoX23ErMpgYBe5yHK/zrH8EzGWJMYMYBfjUeX5w8qQAAC4kRoELfBhtwFUsMQYTo8AtHB/CrSwxxhCjAH/yy4mUAgC4lhgFruEzItzKEmMqMQpcz/Eh3MoSYxIxCvAX7ydSCgDgcmIUuIzPiHArS4yRxChwPceHcCtLjEnEKMDfve39CgDgcmIUuJJcg1tZYsxznO5rAAAiTkYBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYpTVHUf9C24zdTRz7cVc2xk8Gs8kRtmAJ+9GXCwAvkSMAjydVwggJEZZmj1yUy7cLlwpICdG2cB52jL3cByv86x/BABbEaMAT+d9DwiJUdblmG1f4mYLlhiwAjHKHsTN+pQNAP9AjALgfQ/IiFEW5Zhtd+JmcZYYsAgxyjbEzcqUDQD/RowC8Hp53wMiYpQVOWabQdwsyxID1iFGWdFn26Ttc1ku2V5cL2AdYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgE+dRxf+zsAXyVG4Vt9GDHKBq4ydYlNnQteYhTgD87zg/3+OF7nWfwagInEKACs7sP3IphBjMK3+n1HccwGF3rOEps6Fw8kRgH+5Je4UQAA1xKjALABX+qZSozCd/t5R3HMBpd7whKbOhfPJEYB/uI9bhQAwOXEKADswZd6RhKjEHjfURyzwR3eVtbgk+ypc/FMYhTg7972fgUAcDkxCgDb8EbEPMfpvgYAIOJkFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGWd1x1L/gNlNHM9dezLWdwaPxTGIUGrYTuJtVBlsQo2zAjkLu7SZ0K5JzEzKPGGVpHrsAv/NsZBIxCoHjeJ3nj38Al3tbWedpicEGxCgbsKPQen95cCvSer8VYRIxyro8dgE+5L2IScQofLefI9uOApezxGAvYpQ92FGo/HJC71ak4mMRU4lRFuWxC/AH3osYQ4zCt/o9su0ocCFLDLYjRtmGHYXv9+EJvVuR7+djEYOJUVbksQvwV96LmEGMsqLPSnRAoX44woC5pnJptjN1iQ1+KoIYBfgyBQBwFTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaPABY7ja3+n5XoB6xCjwF2Ow//hvqjz1J3AKsQocAFxM4CXByAhRgEAyIhR4BaO2RbnMBtYhBgFriFutublAaiIUQAAMmIUuJ5jti04zAZWIEaBy4ibTXl5AEJiFACAjBgFrueYbReuFJATo8CVxM2OXDUgJEYBAMgcpzdiAAAiTkYBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYpTVHUf9C+5xHJNHG8lce5k612v0aDyTGGUDnrwbcbHgVpYY84hRluaxu6l5F+5tonlzAeTEKASO43WeP/4B8FUeHUwiRtnAeXry7uE9soE7WGKMJEZZl8fuvoa9P7zfisPmAliBGIXv9nNkixvgH3h0MIkYZQ+evOtzkg23ssSYSoyyKI/d3Y15f/jlVhwzF8AixCh8q98jW9wA/8CjgzHEKNvw5F2Zk2y4lSXGYGIUuMuA94cPC2DAXADrEKOsaPAZwIdzDRh2wAiwsqmPDniJUdb02RPWk3c7u18ytyLA3cQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCPM5xfO3vG/lwhAFzwWBiFOBxzlOfAasQowAAZMQoAK/X63Ucr/Osf8T/7PdD3xlzwWBiFOCJfKkHFiFGAQDIiFEARn3L/vnQd9JcMJUYBXgoX+qBFYhRgKdzfAiExCgA07y1tciGLYhRgOfSakBOjAI8nSQFQmIUgJlENmzhOC1WAAAiTkYBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYpTVHUf9C+5xHJNHG8lce5k612v0aDyTGAUAICNGAQDIiFEIHMfrPH/8AwCeTIwCAJARowAAZMQofLf3b/Sv1+s8fakH4NHEKAAAGTEKAEBGjMK3+vkb/Rtf6gF4MjEKAEBGjAIAkBGj8K1++Ub/hz8CwBOIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0ZZ0XF87e/kpl4yc+1l99//B1MvGbzEKGs6T0/YCY7jdZ71j+AjU5fYZ3ONHPZliTGFGAX41IdxowAALiRG2YYCWNnUk7ZHscRWZokxmBhlUZ68u1M2i5u6xH6fa+qtOHUuHkiMAvzJL3GjAACuJUbZgwJY39STtoeYusQmzWWJMZUYZV2evPuaVACDTV1iU+f6mSXGJGIU4C/e40YBAFxOjLIBBbCLJ5xIjTR1ic2byxJjJDHK0oZtJM/hwu1i6pWaOte78QPyKGIU4O/e9n4FAHA5McoGFMBGXKwdTb1qI+caORQPd5zuawAAIk5GAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmKU1R1H/QtuM3U0c+3FXNsZPBrPJEbZwMgn78ihYClWGWxBjAL83VvWTI2bqXMBWxCjLM0eCbeyxICcGGUD5zltyzyO13n++AdwubeVNe/RASOJUYC/eH95mBo3U+cCtiBGWdd7AQB3sMSAFYhR9jDp5ObnApg0FyzCEoO9iFGAP/nl+HBq3EydC1ifGGVRPiDCrSwxYBFilG3MOLn5vQBmzAWLsMRgO2IU4FMfHh9OjZupcwGLE6OsyAdEuJUlBqxDjLKiz7bJAdvnhyMMmGuqqbfi1LlelhhsSIwCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxClzgOL72d+BLPlxK1hcziFHgLsfxOs/6RwCwNjEKXOA8HdLAjSwxBhOjALAfXx4YQ4wCt7BTAvBfiFHgGj4jwq0sMaYSowCwGV8emESMAtezUwLwH4lR4DI+I8KtLDFGEqMAsBlfHphEjALXs1MC8B+JUeBKMhRuZYkxz3G6rwEAiDgZBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAIPP/5m0DJPkGpwgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N=3;\n",
    "mu=0.5;\n",
    "S = mu*sign(rand(2*N+1, 2*N+1)-0.5);\n",
    "\n",
    "DisplayIsing(S);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2989692",
   "metadata": {},
   "source": [
    "To visualize the state of the spin lattice, we can use colour to enhance the previous picture.  While we are at it, we can make the function a little more robust by adding error checking:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4c11f260",
   "metadata": {},
   "outputs": [],
   "source": [
    "### function to display the spin state of the system, on a rectangular lattice of spins\n",
    "function DisplayIsing (S,up_colour='blue', dn_colour='red')\n",
    "  [n,m] = size(S);\n",
    "  if (n != m)\n",
    "    error(\"need a square matrix for S\")\n",
    "    endif\n",
    "  N = round((n-1)/2);\n",
    "  if (2*N+1 != n)\n",
    "    error(\"matrix dimensions are not 2N+1 x 2N+1\")\n",
    "    endif\n",
    "  [X,Y] = meshgrid(-N:N,-N:N);\n",
    "  hold on;\n",
    "  quiver(X,Y-0.5*S,0*S,S,0,'color',dn_colour);\n",
    "  ### zero all negative values in S, repaint those that are +ve in up_colour\n",
    "  S = ifelse(sign(S) == 1, S, 0);\n",
    "  quiver(X,Y-0.5*S,0*S,S,0,'color',up_colour);\n",
    "  hold off;\n",
    "  axis equal;\n",
    "  xlim([-N-1,N+1]);\n",
    "  ylim([-N-1,N+1]);\n",
    "  axis off;\n",
    "endfunction;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "22764cad",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4QAAAJYCAIAAAC1p7+MAAAJMmlDQ1BkZWZhdWx0X3JnYi5pY2MAAEiJlZVnUJNZF8fv8zzphUASQodQQ5EqJYCUEFoo0quoQOidUEVsiLgCK4qINEWQRQEXXJUia0UUC4uCAhZ0gywCyrpxFVFBWXDfGZ33HT+8/5l7z2/+c+bec8/5cAEgiINlwct7YlK6wNvJjhkYFMwE3yiMn5bC8fR0A9/VuxEArcR7ut/P+a4IEZFp/OW4uLxy+SmCdACg7GXWzEpPWeGjy0wPj//CZ1dYsFzgMt9Y4eh/eexLzr8s+pLj681dfhUKABwp+hsO/4b/c++KVDiC9NioyGymT3JUelaYIJKZttIJHpfL9BQkR8UmRH5T8P+V/B2lR2anr0RucsomQWx0TDrzfw41MjA0BF9n8cbrS48hRv9/z2dFX73kegDYcwAg+7564ZUAdO4CQPrRV09tua+UfAA67vAzBJn/eqiVDQ0IgALoQAYoAlWgCXSBETADlsAWOAAX4AF8QRDYAPggBiQCAcgCuWAHKABFYB84CKpALWgATaAVnAad4Dy4Aq6D2+AuGAaPgRBMgpdABN6BBQiCsBAZokEykBKkDulARhAbsoYcIDfIGwqCQqFoKAnKgHKhnVARVApVQXVQE/QLdA66At2EBqGH0Dg0A/0NfYQRmATTYQVYA9aH2TAHdoV94fVwNJwK58D58F64Aq6HT8Id8BX4NjwMC+GX8BwCECLCQJQRXYSNcBEPJBiJQgTIVqQQKUfqkVakG+lD7iFCZBb5gMKgaCgmShdliXJG+aH4qFTUVlQxqgp1AtWB6kXdQ42jRKjPaDJaHq2DtkDz0IHoaHQWugBdjm5Et6OvoYfRk+h3GAyGgWFhzDDOmCBMHGYzphhzGNOGuYwZxExg5rBYrAxWB2uF9cCGYdOxBdhK7EnsJewQdhL7HkfEKeGMcI64YFwSLg9XjmvGXcQN4aZwC3hxvDreAu+Bj8BvwpfgG/Dd+Dv4SfwCQYLAIlgRfAlxhB2ECkIr4RphjPCGSCSqEM2JXsRY4nZiBfEU8QZxnPiBRCVpk7ikEFIGaS/pOOky6SHpDZlM1iDbkoPJ6eS95CbyVfJT8nsxmpieGE8sQmybWLVYh9iQ2CsKnqJO4VA2UHIo5ZQzlDuUWXG8uIY4VzxMfKt4tfg58VHxOQmahKGEh0SiRLFEs8RNiWkqlqpBdaBGUPOpx6hXqRM0hKZK49L4tJ20Bto12iQdQ2fRefQ4ehH9Z/oAXSRJlTSW9JfMlqyWvCApZCAMDQaPkcAoYZxmjDA+SilIcaQipfZItUoNSc1Ly0nbSkdKF0q3SQ9Lf5RhyjjIxMvsl+mUeSKLktWW9ZLNkj0ie012Vo4uZynHlyuUOy33SB6W15b3lt8sf0y+X35OQVHBSSFFoVLhqsKsIkPRVjFOsUzxouKMEk3JWilWqUzpktILpiSTw0xgVjB7mSJleWVn5QzlOuUB5QUVloqfSp5Km8oTVYIqWzVKtUy1R1WkpqTmrpar1qL2SB2vzlaPUT+k3qc+r8HSCNDYrdGpMc2SZvFYOawW1pgmWdNGM1WzXvO+FkaLrRWvdVjrrjasbaIdo12tfUcH1jHVidU5rDO4Cr3KfFXSqvpVo7okXY5upm6L7rgeQ89NL0+vU++Vvpp+sP5+/T79zwYmBgkGDQaPDamGLoZ5ht2GfxtpG/GNqo3uryavdly9bXXX6tfGOsaRxkeMH5jQTNxNdpv0mHwyNTMVmLaazpipmYWa1ZiNsulsT3Yx+4Y52tzOfJv5efMPFqYW6RanLf6y1LWMt2y2nF7DWhO5pmHNhJWKVZhVnZXQmmkdan3UWmijbBNmU2/zzFbVNsK20XaKo8WJ45zkvLIzsBPYtdvNcy24W7iX7RF7J/tC+wEHqoOfQ5XDU0cVx2jHFkeRk4nTZqfLzmhnV+f9zqM8BR6f18QTuZi5bHHpdSW5+rhWuT5z03YTuHW7w+4u7gfcx9aqr01a2+kBPHgeBzyeeLI8Uz1/9cJ4eXpVez33NvTO9e7zofls9Gn2eedr51vi+9hP0y/Dr8ef4h/i3+Q/H2AfUBogDNQP3BJ4O0g2KDaoKxgb7B/cGDy3zmHdwXWTISYhBSEj61nrs9ff3CC7IWHDhY2UjWEbz4SiQwNCm0MXwzzC6sPmwnnhNeEiPpd/iP8ywjaiLGIm0iqyNHIqyiqqNGo62ir6QPRMjE1MecxsLDe2KvZ1nHNcbdx8vEf88filhICEtkRcYmjiuSRqUnxSb7JicnbyYIpOSkGKMNUi9WCqSOAqaEyD0tandaXTlz/F/gzNjF0Z45nWmdWZ77P8s85kS2QnZfdv0t60Z9NUjmPOT5tRm/mbe3KVc3fkjm/hbKnbCm0N39qzTXVb/rbJ7U7bT+wg7Ijf8VueQV5p3tudATu78xXyt+dP7HLa1VIgViAoGN1tubv2B9QPsT8M7Fm9p3LP58KIwltFBkXlRYvF/OJbPxr+WPHj0t6ovQMlpiVH9mH2Je0b2W+z/0SpRGlO6cQB9wMdZcyywrK3BzcevFluXF57iHAo45Cwwq2iq1Ktcl/lYlVM1XC1XXVbjXzNnpr5wxGHh47YHmmtVagtqv14NPbogzqnuo56jfryY5hjmceeN/g39P3E/qmpUbaxqPHT8aTjwhPeJ3qbzJqamuWbS1rgloyWmZMhJ+/+bP9zV6tua10bo63oFDiVcerFL6G/jJx2Pd1zhn2m9az62Zp2WnthB9SxqUPUGdMp7ArqGjzncq6n27K7/Ve9X4+fVz5ffUHyQslFwsX8i0uXci7NXU65PHsl+spEz8aex1cDr97v9eoduOZ67cZ1x+tX+zh9l25Y3Th/0+LmuVvsW523TW939Jv0t/9m8lv7gOlAxx2zO113ze92D64ZvDhkM3Tlnv296/d5928Prx0eHPEbeTAaMip8EPFg+mHCw9ePMh8tPN4+hh4rfCL+pPyp/NP637V+bxOaCi+M24/3P/N59niCP/Hyj7Q/Fifzn5Ofl08pTTVNG02fn3Gcufti3YvJlykvF2YL/pT4s+aV5quzf9n+1S8KFE2+Frxe+rv4jcyb42+N3/bMec49fZf4bmG+8L3M+xMf2B/6PgZ8nFrIWsQuVnzS+tT92fXz2FLi0tI/QiyQvpTNDAsAAAAJcEhZcwAACxMAAAsTAQCanBgAAAAddEVYdFNvZnR3YXJlAEdQTCBHaG9zdHNjcmlwdCA5LjI1wZk/DQAAEUNJREFUeJzt3cFynbq2QNHN+/9/5jWccqUSO4nPBU9pMUYr5dZeBUIT0chxnucLAAAK/1f/AAAAnkuMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKKs7jvoX3OM4Jo82krn2YonBLsQoG/Dk3YiLtSNXbSMuFvOIUeB68/bLt4mmzgUQEqMsbepOeRyv8/zxD+BybyvrPMcusalz8UxilA0M3lGGeY9s9mKJ7cISYyQxClxvWNy8F8DUuQBCYpR1Td0pf55rWNzACp6wxKbOxTOJUfbgybu+qS8PD2GJrc8SYyoxCtxiTNz8UgBT5wKoiFEWNXWn/H2uMXEDK3jOEps6Fw8kRtmGJ+/Kpr48PIoltjJLjMHEKHCXAXHzYQFMnQsgIUZZ0eCd8sO5Bgw7YARmeNQSGzAXvMQoa/rsCevJu53dL9nUW3HqXMCOxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkDnO86x/AwAAD+VkFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGWd1x1L/gNmNHmzrY0LmGjjV2rtdr9mw8kRhlAyMfvCOHYlNT78apc8EwYhQAgIwYZWlTDzaO43WeP/4Boal34NS5YCQxygbO09YCN5q6xKbOBcOIUQAAMmKUdb1/yx7m57mc3BB6whID1idG2YNog1tNXWJT54JJxCgAABkxyqKmfmj7fS4nNySes8SAxYlRtiHa4FZTl9jUuWAMMQoAQEaMsqgPP7QN+Po2dS62M/VWnDoXDCZGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAWucBxf+/sups7FdtyKzCVGgdscx+s86x8Bc1lijCBGAT53nh+cPCkAgOuIUeAKH0YbcBVLjLnEKHAPx4dwK0uMKcQowB/9ciKlAAAuJUaBi/iMCLeyxBhKjAI3cHwIt7LEGESMAvzN+4mUAgC4mhgFruMzItzKEmMiMQrcwPEh3MoSYxAxCvAP3vZ+BQBwNTEKXEquwa0sMcY5Trc1AAARJ6MAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCAJARowAAZMQoAAAZMcrqjqP+BfcZOtvQscYONnSswYONHo1HEqNswIN3Iy4WAF8iRgGezisEEBKjLM0euSkXbheuFJATo2zgPG2ZeziO13nWPwKArYhRgKfzvgeExCjrcsy2L3GzBUsMWIEYZQ/iZn3KBoD/QIwC4H0PyIhRFuWYbXfiZnGWGLAIMco2xM3KlA0A/40YBeD18r4HRMQoK3LMNoO4WZYlBqxDjLKiz7ZJ2+eyXLK9uF7AOsQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCfO44vvZ3AL5IjML3+jBilA1cZeoSmzoXiFGAPznPD/b743idZ/FrAAYSowCwvA/fi2AEMQrf6/cdxTEbXOg5S2zqXDyPGAX4o1/iRgEAXEqMAsAOfKlnKDEK3+7nHcUxG1zuCUts6lw8khgF+Jv3uFEAAFcTowCwCV/qmUiMQuF9R3HMBnd4W1mDT7KnzsUjiVGAf/C29ysAgKuJUQDYhzcixjlOtzUAABEnowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxyuqOo/4Ftxk72tDBho41drChY71eo0fjmcQoNGwncDerDLYgRtmAHYXc203oViTnJmQeMcrSPHYBfufZyCRiFALH8TrPH/8ALve2ss7TEoMNiFE2YEeh9f7y4Fak9X4rwiRilHV57AJ8yHsRk4hR+G4/R7YdBS5nicFexCh7sKNQ+eWE3q1IxcciphKjLMpjF+APvBcxhhiFb/V7ZNtR4EKWGGxHjLINOwrf78MTerci38/HIgYTo6zIYxfgr7wXMYMYZUWfleiAQv1whAFzTeXSbGfqEhv8VAQxCvBlCgDgKmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0aBKxzH1/5Oy/UCliFGgdsch//EfVHnqTuBRYhR4AriZgAvD0BBjAIAkBGjwD0csy3OYTawBjEKXETcbM3LAxARowAAZMQocAPHbFtwmA0sQIwC1xE3m/LyAHTEKAAAGTEK3MAx2y5cKaAmRoFLiZsduWpAR4wCAJA5Ti/EAABEnIwCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCirO476F9zjOCaPNtLUucYONnWuyZPxUGKUDXjybsTFgltZYswjRlmax+6m5l24t4nmzQWQE6MQOI7Xef74B8BXeXQwiRhlA+fpybuH98gG7mCJMZIYZV0eu/sa9v7wfisOmwtgBWIUvtvPkS1ugP/Ao4NJxCh78ORdn5NsuJUlxlRilEV57O5uzPvDL7fimLkAFiFG4Vv9HtniBvgPPDoYQ4yyDU/elTnJhltZYgwmRoG7DHh/+LAABswFsA4xyooGnwF8ONeAYQeMACub+uiAlxhlTZ89YT15t7P7JXMrAtxNjAIAkBGjAABkxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKMDzHMfX/r6RD0cYMBfMJUYBnuc89RmwCDEKAEBGjALwer1er+N4nWf9I/5nvx/6zpgL5hKjAI/kSz2wBjEKAEBGjAIw61v2z4e+k+aCocQowFP5Ug8sQIwCPJ7jQ6AjRgEY562tRTbsQIwCPJhWA2piFODxJCnQEaMADCWyYQfHaa0CABBxMgoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCAJARoyzvOOpfcIvjmDrZ2Etmrs1Mnes1ejQeSYwCAJARowAAZMQoBI7jdZ4//gEATyZGAQDIiFEAADJiFL7b+zf61+t1nr7UA/BoYhQAgIwYBQAgI0bhW/38jf6NL/UAPJkYBQAgI0YBAMiIUfhWv3yj/8MfAeAJxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjLOk4vvZ3clMvmbn2svvv/4OplwzEKIs6T0/YCY7jdZ71j+AjU5fYZ3ONHPZliTGEGAX43IdxowAAriNG2YcCWNnUk7ZHscRWZokxlxhlVZ68u1M2i5u6xH6fa+qtOHUunkeMAvzRL3GjAAAuJUbZhAJY39STtoeYusQmzWWJMZQYZWGevPuaVACDTV1iU+f6mSXGIGIU4G/e40YBAFxNjLIDBbCLJ5xIjTR1ic2byxJjIjHK2oZtJM/hwu1i6pWaOte78QPyJGIU4B+87f0KAOBqYpQdKICNuFg7mnrVRs41ciie7Tjd1gAARJyMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCAJARowAAZMQoAAAZMQoAQEaMAgCQEaMAAGTEKAAAGTEKAEBGjAIAkBGjAABkxCgAABkxCgBARowCAJARowAAZMQoqzuO+hfcZuxoQwcbOpa5NjR5Np5IjLKBkQ/ekUPBUqwy2IIYBfi7t6yZGjdT5wK2IEZZmj0SbmWJATkxygbOc9qWeRyv8/zxD+Bybytr3qMDRhKjAH/x/vIwNW6mzgVsQYyyrvcCAO5giQErEKPsYdLJzc8FMGkuWIQlBnsRowB/8svx4dS4mToXsD4xyqJ8QIRbWWLAIsQo25hxcvN7AcyYCxZhicF2xCjApz48PpwaN1PnAhYnRlmRD4hwK0sMWIcYZUWfbZMDts8PRxgw11RTb8Wpc70sMdiQGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUuMJxfO3vwJd8uJSsL0YQo8BtjuN1nvWPAGBpYhS4wnk6pIEbWWLMJUYBYEO+PDCFGAXuYacE4B+IUeAiPiPCrSwxhhKjALAbXx4YRIwCN7BTAvBvxChwHZ8R4VaWGBOJUQDYjS8PDCJGgRvYKQH4N2IUuJQMhVtZYoxznG5rAAAiTkYBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGAQDIiFEAADJiFACAjBgFACAjRgEAyIhRAAAyYhQAgIwYBQAgI0YBAMiIUQAAMmIUAICMGAUAICNGAQDI/D/ENAMkSd+gVgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "DisplayIsing(S);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b50e82b2",
   "metadata": {},
   "source": [
    "## System evolves in time according to the rules encoded in `Evolve()`\n",
    "\n",
    "Just as a demonstration, go through all spins in the system and randomly \"flip\" some of them. Setting the cutoff at the midpoint of the 0..1 interval should flip about a half of them. Going over the entire spin system once corresponds to a single step in the time evolution of the system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "2af059c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "### a single step in the evolution of the spin system\n",
    "function S_out = Evolve (S_in)\n",
    "  [n,m] = size(S_in);\n",
    "  if (n != m)\n",
    "    error(\"need a square matrix for S\")\n",
    "    endif\n",
    "  N = round((n-1)/2);\n",
    "  if (2*N+1 != n)\n",
    "    error(\"matrix dimensions are not 2N+1 x 2N+1\")\n",
    "    endif\n",
    "  S_out = S_in;\n",
    "  for i = [1:n]\n",
    "    for j = [1:m]\n",
    "      ### for demo purposes, one possible \"evolution\" is to randomly flip some spins       \n",
    "      if (rand() > 0.5)\n",
    "        S_out(i,j) = - S_in(i,j);\n",
    "      endif\n",
    "    endfor\n",
    "  endfor\n",
    "endfunction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "48f01104",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4QAAAJYCAIAAAC1p7+MAAAJMmlDQ1BkZWZhdWx0X3JnYi5pY2MAAEiJlZVnUJNZF8fv8zzphUASQodQQ5EqJYCUEFoo0quoQOidUEVsiLgCK4qINEWQRQEXXJUia0UUC4uCAhZ0gywCyrpxFVFBWXDfGZ33HT+8/5l7z2/+c+bec8/5cAEgiINlwct7YlK6wNvJjhkYFMwE3yiMn5bC8fR0A9/VuxEArcR7ut/P+a4IEZFp/OW4uLxy+SmCdACg7GXWzEpPWeGjy0wPj//CZ1dYsFzgMt9Y4eh/eexLzr8s+pLj681dfhUKABwp+hsO/4b/c++KVDiC9NioyGymT3JUelaYIJKZttIJHpfL9BQkR8UmRH5T8P+V/B2lR2anr0RucsomQWx0TDrzfw41MjA0BF9n8cbrS48hRv9/z2dFX73kegDYcwAg+7564ZUAdO4CQPrRV09tua+UfAA67vAzBJn/eqiVDQ0IgALoQAYoAlWgCXSBETADlsAWOAAX4AF8QRDYAPggBiQCAcgCuWAHKABFYB84CKpALWgATaAVnAad4Dy4Aq6D2+AuGAaPgRBMgpdABN6BBQiCsBAZokEykBKkDulARhAbsoYcIDfIGwqCQqFoKAnKgHKhnVARVApVQXVQE/QLdA66At2EBqGH0Dg0A/0NfYQRmATTYQVYA9aH2TAHdoV94fVwNJwK58D58F64Aq6HT8Id8BX4NjwMC+GX8BwCECLCQJQRXYSNcBEPJBiJQgTIVqQQKUfqkVakG+lD7iFCZBb5gMKgaCgmShdliXJG+aH4qFTUVlQxqgp1AtWB6kXdQ42jRKjPaDJaHq2DtkDz0IHoaHQWugBdjm5Et6OvoYfRk+h3GAyGgWFhzDDOmCBMHGYzphhzGNOGuYwZxExg5rBYrAxWB2uF9cCGYdOxBdhK7EnsJewQdhL7HkfEKeGMcI64YFwSLg9XjmvGXcQN4aZwC3hxvDreAu+Bj8BvwpfgG/Dd+Dv4SfwCQYLAIlgRfAlxhB2ECkIr4RphjPCGSCSqEM2JXsRY4nZiBfEU8QZxnPiBRCVpk7ikEFIGaS/pOOky6SHpDZlM1iDbkoPJ6eS95CbyVfJT8nsxmpieGE8sQmybWLVYh9iQ2CsKnqJO4VA2UHIo5ZQzlDuUWXG8uIY4VzxMfKt4tfg58VHxOQmahKGEh0SiRLFEs8RNiWkqlqpBdaBGUPOpx6hXqRM0hKZK49L4tJ20Bto12iQdQ2fRefQ4ehH9Z/oAXSRJlTSW9JfMlqyWvCApZCAMDQaPkcAoYZxmjDA+SilIcaQipfZItUoNSc1Ly0nbSkdKF0q3SQ9Lf5RhyjjIxMvsl+mUeSKLktWW9ZLNkj0ie012Vo4uZynHlyuUOy33SB6W15b3lt8sf0y+X35OQVHBSSFFoVLhqsKsIkPRVjFOsUzxouKMEk3JWilWqUzpktILpiSTw0xgVjB7mSJleWVn5QzlOuUB5QUVloqfSp5Km8oTVYIqWzVKtUy1R1WkpqTmrpar1qL2SB2vzlaPUT+k3qc+r8HSCNDYrdGpMc2SZvFYOawW1pgmWdNGM1WzXvO+FkaLrRWvdVjrrjasbaIdo12tfUcH1jHVidU5rDO4Cr3KfFXSqvpVo7okXY5upm6L7rgeQ89NL0+vU++Vvpp+sP5+/T79zwYmBgkGDQaPDamGLoZ5ht2GfxtpG/GNqo3uryavdly9bXXX6tfGOsaRxkeMH5jQTNxNdpv0mHwyNTMVmLaazpipmYWa1ZiNsulsT3Yx+4Y52tzOfJv5efMPFqYW6RanLf6y1LWMt2y2nF7DWhO5pmHNhJWKVZhVnZXQmmkdan3UWmijbBNmU2/zzFbVNsK20XaKo8WJ45zkvLIzsBPYtdvNcy24W7iX7RF7J/tC+wEHqoOfQ5XDU0cVx2jHFkeRk4nTZqfLzmhnV+f9zqM8BR6f18QTuZi5bHHpdSW5+rhWuT5z03YTuHW7w+4u7gfcx9aqr01a2+kBPHgeBzyeeLI8Uz1/9cJ4eXpVez33NvTO9e7zofls9Gn2eedr51vi+9hP0y/Dr8ef4h/i3+Q/H2AfUBogDNQP3BJ4O0g2KDaoKxgb7B/cGDy3zmHdwXWTISYhBSEj61nrs9ff3CC7IWHDhY2UjWEbz4SiQwNCm0MXwzzC6sPmwnnhNeEiPpd/iP8ywjaiLGIm0iqyNHIqyiqqNGo62ir6QPRMjE1MecxsLDe2KvZ1nHNcbdx8vEf88filhICEtkRcYmjiuSRqUnxSb7JicnbyYIpOSkGKMNUi9WCqSOAqaEyD0tandaXTlz/F/gzNjF0Z45nWmdWZ77P8s85kS2QnZfdv0t60Z9NUjmPOT5tRm/mbe3KVc3fkjm/hbKnbCm0N39qzTXVb/rbJ7U7bT+wg7Ijf8VueQV5p3tudATu78xXyt+dP7HLa1VIgViAoGN1tubv2B9QPsT8M7Fm9p3LP58KIwltFBkXlRYvF/OJbPxr+WPHj0t6ovQMlpiVH9mH2Je0b2W+z/0SpRGlO6cQB9wMdZcyywrK3BzcevFluXF57iHAo45Cwwq2iq1Ktcl/lYlVM1XC1XXVbjXzNnpr5wxGHh47YHmmtVagtqv14NPbogzqnuo56jfryY5hjmceeN/g39P3E/qmpUbaxqPHT8aTjwhPeJ3qbzJqamuWbS1rgloyWmZMhJ+/+bP9zV6tua10bo63oFDiVcerFL6G/jJx2Pd1zhn2m9az62Zp2WnthB9SxqUPUGdMp7ArqGjzncq6n27K7/Ve9X4+fVz5ffUHyQslFwsX8i0uXci7NXU65PHsl+spEz8aex1cDr97v9eoduOZ67cZ1x+tX+zh9l25Y3Th/0+LmuVvsW523TW939Jv0t/9m8lv7gOlAxx2zO113ze92D64ZvDhkM3Tlnv296/d5928Prx0eHPEbeTAaMip8EPFg+mHCw9ePMh8tPN4+hh4rfCL+pPyp/NP637V+bxOaCi+M24/3P/N59niCP/Hyj7Q/Fifzn5Ofl08pTTVNG02fn3Gcufti3YvJlykvF2YL/pT4s+aV5quzf9n+1S8KFE2+Frxe+rv4jcyb42+N3/bMec49fZf4bmG+8L3M+xMf2B/6PgZ8nFrIWsQuVnzS+tT92fXz2FLi0tI/QiyQvpTNDAsAAAAJcEhZcwAACxMAAAsTAQCanBgAAAAddEVYdFNvZnR3YXJlAEdQTCBHaG9zdHNjcmlwdCA5LjI1wZk/DQAAFQdJREFUeJzt3dFVG1e7x+HZZ50CDB2ElEDSAS2EDmyXIEowdGCXYHyVW0gHdglx7nJLOphzoXiWjh1jJA37P3r1PFd8Wfa3mc3o1U8jM2rjOA4AAJDwP+lvAACA4yVGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUYq7v79vrd3f3z/+x3777bfW2mq16vNdAcxlc3ytVqsfjjtYmv9NfwOwCB8+fHj58uXr16/T3wjAdjbH183NzTAMFxcX6W8KtuDKKEfhzZs3rbWTk5N3794Nw7BarVprm9cShmG4v79fX1FYrVYnJyettbOzs0+fPg1fLq+u//u3fx2gj4eHh4uLi/X8ubi4eHh42Bxf669vb2/Xg+urSfXVHIPlcGWUo/Dw8HB3d7darVar1d9//31zc/P+/fthGC4vL09PTy8vL29ubi4uLi4uLu7v729ubl6+fHl5eblardbjfv1/cnt7++rVq+vr66/+uiQF+nj37t0ff/wxzZ/b29vN8XV+fr7++uzs7NtJdX5+PnyZY+HDgP/PlVGOwvX19cXFxdXV1T///PP777//9NNPHz9+/Pjx44sXL+7u7tYz+uTk5Ozs7O7ubhiGd+/eTX9++gdY19fX19fXd3d3X/315IEBx2S1Wq2Hz9u3b4dh+PPPPzfH1/T1ycnJ9ybVeo4FDwG+JUY5It++OfXq1avLy8s+fx1gT9fX17/88sswDFv9A/fNSeU9ehZIjHIU1u9Vra8l/Prrrw8PD1dXV1dXV58/f57ehV/7+eefh2FY/9PS9+/fv3jxYn2xYXJ+fv7IXwd4Puvrnevh8/ifNKk4IGKUo3B7e9ta+/Dhw9u3b9+8eXN2dnZ6enp6evr58+ev/vnUq1ev1r+Xuv7z19fXX11IuLq6euSvAzyf169f//XXX6enp9/7B0IvXry4ubm5v783qTggbRzH9PcAPdzf35+dnZ2dnU3/8+Tk5KurnpNPnz49PDycn59/7y2tx/86wDNZX+b83vB5eHj49OnTNLtMKg6CGAUAIMbb9AAAxIhRAABixCgAADFiFACAGDEKAECMGAUAIEaMAgAQI0YBAIgRowAAxIhRAABixCgAADFiFACAGDEKAECMGAUAIEaMAgAQI0YBAIgRowAAxIhRAABixCgAADFiFACAGDEKAECMGAUAIEaMAgAQI0YBAIgRowAAxIhRAABixCgAADFiFACAGDEKAECMGAUAIEaMAgAQI0YBAIgRowAAxIhRAABixCgAADFiFACAGDEKAECMGK2jtaqLdV2t8Db2XK7zkVFAG/pOlY7L9X6g99zJstO59wl55MRoEYWf+wsfGhDRvzPMMXiEGC3FvNvTegNt455sIwehT5J6IMAPiVEOg4EOACWJ0QpaG8Yx/U08j8KHBkS0oY1D77Eyjl5Rw3eJ0VLMu31M4Wsb92EbORTjMD73O/VeTsNTiFEOg7IBgJLEaAWbr7yLvQrveWiFt7En28jCbb5H/9zv13s4wFOIUQAAYsQoAAAxYhQAgBgxCgBAjBhlS5u/0+7323fWeRv91Fiwr+6v5DPBd7a5dc++jT2nSt8J1nUbGYZBjLI191iaRWob3faQ5dm832fkjvT1VNvG0MCsto0LJkYBgC1NgeglLnsTo+zKAJpFn210PZtD4ELUnjp8plRYl4FZfxuXR4yyPQ06i8g2+tmxSBp0XgX3MzG7Cm7jUolRAGB760D0Epe9iVF2ZQDNots2+nmxeC5E7a/4HvaaY8W3cXna6CkKAIAQV0YBAIgRowAAxIhRAABixCgAADFiFACAGDEKAECMGAUAIEaMAgAQI0YBAIgRowAAxIhRAABixCgAADFiFACAGDEKAECMGAUAIEaMAgAQI0braK3qYl1Xs40HuBgVtKHvw6Hjcp0fDV0HZuFt7HtCHjkxytKpmlnYRjgGtR/p3Y5OiXYmRouoPYC6WW+jzYTa+qdGyami2JiLGC2l5LxjRs4QmPRpqfIPusLbqLa7EaMsWmvDOP77BTubthEozMCcRRvaOJiYXYnRCqTGLKZtHEejHMqKpEa9qWIbmZEYLaXwA7XwofVkG2EyDuNzvw97DFcKCm9jh0NjTYxWsPkoLTb4eh6abYRjsHk977mv7XV+6E1L9Fir8DZ2PDTWxCgAADFiFACAGDEKAECMGAUAIEaMArPa/HV9v7oPwI+IUbYkNWZhG2HDdAMdd9KBIyRG2ZI7Vc6i8DZOh3YMt1gEYG9ilF1JjVnYRgCOmxgFIGx9a3GfCQ7HSYyyPVfyZlF4G9eHVvgAAZiPGGVXUmMWthGA49ZGz4UAAIS4MgoAQIwYBQAgRowCABAjRgEAiBGjAADEiFEAAGLEKAAAMWIUAIAYMQoAQIwYBQAgRowCABAjRgEAiBGjAADEiFEAAGLEKAAAMWIUAIAYMQoAQIwYraO1mou11vXQqq7Ve7nOx8bh6/1IHzoOsY5rdV6u6sBsQ+v8UztyYrQIT/2z6L+N3VZ0hsBEZ8yi8MCkMzFaigfqLPpsox8WZEnSWRiY7E+MsmitDeP47xfszzayTOszcxydouS1oY3DOHi50pEYrWAqNvYR2cY+z77OEJhMqcE+Cg9M+hOjpXigzqLDNqpDiBuH0aWv/RmY7E+MsmjTACo2iTYPp8OhVd1Gyuh5irosenB6D8wvZ4hTpRsxWkHnB2pVthGOjdrYmYHJjMQoAAAxYhQAgBgxCgBAjBgFACBGjLKlzXt4uI8U33KGwBebt47qcBupaQm3rOKwiFG25F6mAE/jVqbwFGKUXbkNMf9perniDIEvfOwTPEKMAkAF04VY4cthEaNsz+UugKfRhfBDYpRdSVK+Z31uOEPgC0kKjxCjAFCE6uUQtdHVCwAAQlwZBQAgRowCABAjRgEAiBGjAADEiFEAAGLEKAAAMWIUAIAYMQoAQIwYBQAgRowCABAjRgEAiBGjAADEiFEAAGLEKAAAMWIUAIAYMQoAQIwYraO1mmt1Xq/roXXex6rbSAlt6HrSdB2YnQ+t53J1n3g6/9SOnBhlR2oDOFD9x5eygUeI0SLWs1Ug7sk2zsI2wkSGwg+JUbYmMoACjDJYCDHKjsbRKAcOT2vDOPZedBxGl0jhe8RoBdNsFYj7sI2zsI0cig6naBvaOHQvXzg0YrSCzVf5HV7xT0v0v7rwrDpvY1W2kYXreYpulqgqhe8RowAAxIhRAABixCgAADFiFACAGDEKAECMGGXZNu+8UulGQZ2Pq+o2UoIbcD6HZ9/VnlOl7wTb3DonZx9ilAMRuVH180ndhLPYNlLC5g3hPffvY3MbS91JKjQwq23jgolRAGBLUyB6icvexCjLVv5jfPrM8fLbCBzDh452GZj1t3F5xCgHot4r78gR1dtGSli/Gepd0bkU3MbE7Cq4jUslRgGA7a0D0Utc9iZGWbzak67b0dXeRkpwIWp/xfew1xwrvo3L00ZPUQAAhLgyCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxWkdrVRfru1zVtfou1/nIKKANXU+anss5tJkW67uNhlhHYhQAiugcvlUp0c7EaBEeOSzH+mx0TrJY/YutZCOWPCgixGgpnv4BnqhPS5UvtvIHSAdiFAAqaEMbh3EQiPtpbRjHf7+gDzFawfTIgbjpbBxHo5wlmoqtp3EYiwViZBupSoyW4ukf4Ik6BOIxFFu9zqY/MVrB5mVRl0jJcjaycJt1+Nyl2HOtzSV6rtVnuZ6mwWWCdSNGAQCIEaMAAMSIUQAAYsQoAAAxYhS+2LwTgbsSAEAXYpQtdS62aYlidWgbYcN0byA3CeI/uFJQnRhlS8dwL9MOnyJwDNsIMAsDszoxyq7qfe5TZN7V20aAZ2JgFiVGWbYpEA2gfdhGlm36RPVit08HnkKMsr3yQdPnAMtvI8BcDMzSxCi7Kjka+h9UtxVL/ryAo2KOFdVGP1oAAEJcGQUAIEaMAgAQI0YBAIgRowAAxIhRAABixCgAADFiFACAGDEKAECMGAUAIEaMAgAQI0YBAIgRowAAxIhRAABixCgAADFiFACAGDEKAECMGAUAIEaM1tFa1cW6LtfzyOruYvdj4/C11vfRN3ScKh3X6rxc4bU6/9SOnBgtwlP/LPpvY7cVnSEw6d8ZygYeIUZLERyzsI1wDPoEogyFHxKjkBFJXp3NMq3PzHF0ipLXhjYO4+CFREditILWhnFMfxOHL7KNfZ59nSEwmVKjp3EYlQ18jxgtxXWFWXTYRnUIcR0CMRK+cHDEKARslmiHKp2WUMAsk1OU5ZheP3gh0Y0YraBz2VRlG+EYbBbGc9dGz7XgcIlRAABixCgAADFiFACAGDEKAECMGGVLmzc96nOHzJ7LsT8/MgiZ7lTV4Z6mm0t0uENWt7WIEKNsKXUvU3fmBFieDvdS9ZEB5YlRdqUO+U/TyxVnCNQlEJmRGGXZfKgUwNNMgVj1nqY+0aoqMcr2Ipe7XGMDOFYatDYxyq7UId+zPjecIVDaOhBlIvsToyyepgF4mtppWPvojlkbPdMDABDiyigAADFiFACAGDEKAECMGAUAIEaMAgAQI0YBAIgRowAAxIhRAABixCgAADFiFACAGDEKAECMGAUAIEaMAgAQI0YBAIgRowAAxIhRAABixGgdrVVdrOtqtvEAF6OC3qdMx/U6H1obOh5a0bWGwRDrSoyydAbCLGwjy+csheMkRotYD3GjfE+2EZhX4amyvlTZ+4IlFYlRDkPJUd6fbWSZnJlwzMQoi9baMI7pb+Lw2UYOwjiWrdLChwb7E6MVTKlh3u3DNgLzKvw6sA1tHMZhGMZh9E49exKjHAaBOAvbyAIVLjbgKcRoBZtzvNhM73lothEiCp+flQ9tGP/za9iBGAUAIEaMAgAQI0YBAIgRowAAxIhRtrT5y9h+MXtnhbex8KFRg1P00GzeOson1JckRiHBPZaAuU2hVuzGn71vZToNZ3cd60WMsiUVxeOmM8QcZ+GcorAMYpRdmeOzsI3QmVfUB2v62Kfn5QzpToxCiAYFZjW9nV3vLvSZIzKlexGjbM/jk8etzxDnCQvnFIVlEKPsyhyfhW2E/jzuDlO/66POkL7aaMcBAAhxZRQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYsQoAAAxYhQAgBgxCgBAjBgFACBGjAIAECNGAQCIEaMAAMSIUQAAYv4PXPJeqO6q1ggAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N=3;\n",
    "rn = rand(2*N+1, 2*N+1);\n",
    "mu=0.5;\n",
    "S = mu*sign(rn-0.5);\n",
    "\n",
    "clf;\n",
    "subplot(1,2,1);\n",
    "DisplayIsing(S);\n",
    "title(\"before\");\n",
    "\n",
    "### apply a single step of evolution of the system, and see the result\n",
    "S_now = Evolve(S);\n",
    "  \n",
    "subplot(1,2,2);\n",
    "hold on;\n",
    "DisplayIsing(S_now);\n",
    "DisplayIsing(0.5*(S_now-S),'green','green');\n",
    "hold off;\n",
    "title(\"after\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "224f6588",
   "metadata": {},
   "source": [
    "## Monte-Carlo and $\\pi$\n",
    "\n",
    "The principle of Monte-Carlo sampling can be well illustrated by the following method of determining $\\pi$, from the ratio of the area of a unit circle, $S_c = \\pi r^2 = \\pi$ and the area of the square that circumscribes it, $S_s = 2\\times 2 = 4$.  Imagine randomly dropping stones anywhere inside the square, but counting how many stones fall inside the circle ($N_{in}$), and how many outside of the circle ($N_{out}$). If the probability of dropping a stone anywhere inside the square is equal, then the fraction of the total number of stones that fall inside the circle will correspond to the fraction that the area of the square that is occupied by the circle, i.e. \n",
    "$$\n",
    "  \\frac{N_{in}}{N_{in}+N_{out}} = \\frac{S_c}{S_s} = \\frac{\\pi}{4}  \n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "db44fca8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4QAAAJYCAIAAAC1p7+MAAAJMmlDQ1BkZWZhdWx0X3JnYi5pY2MAAEiJlZVnUJNZF8fv8zzphUASQodQQ5EqJYCUEFoo0quoQOidUEVsiLgCK4qINEWQRQEXXJUia0UUC4uCAhZ0gywCyrpxFVFBWXDfGZ33HT+8/5l7z2/+c+bec8/5cAEgiINlwct7YlK6wNvJjhkYFMwE3yiMn5bC8fR0A9/VuxEArcR7ut/P+a4IEZFp/OW4uLxy+SmCdACg7GXWzEpPWeGjy0wPj//CZ1dYsFzgMt9Y4eh/eexLzr8s+pLj681dfhUKABwp+hsO/4b/c++KVDiC9NioyGymT3JUelaYIJKZttIJHpfL9BQkR8UmRH5T8P+V/B2lR2anr0RucsomQWx0TDrzfw41MjA0BF9n8cbrS48hRv9/z2dFX73kegDYcwAg+7564ZUAdO4CQPrRV09tua+UfAA67vAzBJn/eqiVDQ0IgALoQAYoAlWgCXSBETADlsAWOAAX4AF8QRDYAPggBiQCAcgCuWAHKABFYB84CKpALWgATaAVnAad4Dy4Aq6D2+AuGAaPgRBMgpdABN6BBQiCsBAZokEykBKkDulARhAbsoYcIDfIGwqCQqFoKAnKgHKhnVARVApVQXVQE/QLdA66At2EBqGH0Dg0A/0NfYQRmATTYQVYA9aH2TAHdoV94fVwNJwK58D58F64Aq6HT8Id8BX4NjwMC+GX8BwCECLCQJQRXYSNcBEPJBiJQgTIVqQQKUfqkVakG+lD7iFCZBb5gMKgaCgmShdliXJG+aH4qFTUVlQxqgp1AtWB6kXdQ42jRKjPaDJaHq2DtkDz0IHoaHQWugBdjm5Et6OvoYfRk+h3GAyGgWFhzDDOmCBMHGYzphhzGNOGuYwZxExg5rBYrAxWB2uF9cCGYdOxBdhK7EnsJewQdhL7HkfEKeGMcI64YFwSLg9XjmvGXcQN4aZwC3hxvDreAu+Bj8BvwpfgG/Dd+Dv4SfwCQYLAIlgRfAlxhB2ECkIr4RphjPCGSCSqEM2JXsRY4nZiBfEU8QZxnPiBRCVpk7ikEFIGaS/pOOky6SHpDZlM1iDbkoPJ6eS95CbyVfJT8nsxmpieGE8sQmybWLVYh9iQ2CsKnqJO4VA2UHIo5ZQzlDuUWXG8uIY4VzxMfKt4tfg58VHxOQmahKGEh0SiRLFEs8RNiWkqlqpBdaBGUPOpx6hXqRM0hKZK49L4tJ20Bto12iQdQ2fRefQ4ehH9Z/oAXSRJlTSW9JfMlqyWvCApZCAMDQaPkcAoYZxmjDA+SilIcaQipfZItUoNSc1Ly0nbSkdKF0q3SQ9Lf5RhyjjIxMvsl+mUeSKLktWW9ZLNkj0ie012Vo4uZynHlyuUOy33SB6W15b3lt8sf0y+X35OQVHBSSFFoVLhqsKsIkPRVjFOsUzxouKMEk3JWilWqUzpktILpiSTw0xgVjB7mSJleWVn5QzlOuUB5QUVloqfSp5Km8oTVYIqWzVKtUy1R1WkpqTmrpar1qL2SB2vzlaPUT+k3qc+r8HSCNDYrdGpMc2SZvFYOawW1pgmWdNGM1WzXvO+FkaLrRWvdVjrrjasbaIdo12tfUcH1jHVidU5rDO4Cr3KfFXSqvpVo7okXY5upm6L7rgeQ89NL0+vU++Vvpp+sP5+/T79zwYmBgkGDQaPDamGLoZ5ht2GfxtpG/GNqo3uryavdly9bXXX6tfGOsaRxkeMH5jQTNxNdpv0mHwyNTMVmLaazpipmYWa1ZiNsulsT3Yx+4Y52tzOfJv5efMPFqYW6RanLf6y1LWMt2y2nF7DWhO5pmHNhJWKVZhVnZXQmmkdan3UWmijbBNmU2/zzFbVNsK20XaKo8WJ45zkvLIzsBPYtdvNcy24W7iX7RF7J/tC+wEHqoOfQ5XDU0cVx2jHFkeRk4nTZqfLzmhnV+f9zqM8BR6f18QTuZi5bHHpdSW5+rhWuT5z03YTuHW7w+4u7gfcx9aqr01a2+kBPHgeBzyeeLI8Uz1/9cJ4eXpVez33NvTO9e7zofls9Gn2eedr51vi+9hP0y/Dr8ef4h/i3+Q/H2AfUBogDNQP3BJ4O0g2KDaoKxgb7B/cGDy3zmHdwXWTISYhBSEj61nrs9ff3CC7IWHDhY2UjWEbz4SiQwNCm0MXwzzC6sPmwnnhNeEiPpd/iP8ywjaiLGIm0iqyNHIqyiqqNGo62ir6QPRMjE1MecxsLDe2KvZ1nHNcbdx8vEf88filhICEtkRcYmjiuSRqUnxSb7JicnbyYIpOSkGKMNUi9WCqSOAqaEyD0tandaXTlz/F/gzNjF0Z45nWmdWZ77P8s85kS2QnZfdv0t60Z9NUjmPOT5tRm/mbe3KVc3fkjm/hbKnbCm0N39qzTXVb/rbJ7U7bT+wg7Ijf8VueQV5p3tudATu78xXyt+dP7HLa1VIgViAoGN1tubv2B9QPsT8M7Fm9p3LP58KIwltFBkXlRYvF/OJbPxr+WPHj0t6ovQMlpiVH9mH2Je0b2W+z/0SpRGlO6cQB9wMdZcyywrK3BzcevFluXF57iHAo45Cwwq2iq1Ktcl/lYlVM1XC1XXVbjXzNnpr5wxGHh47YHmmtVagtqv14NPbogzqnuo56jfryY5hjmceeN/g39P3E/qmpUbaxqPHT8aTjwhPeJ3qbzJqamuWbS1rgloyWmZMhJ+/+bP9zV6tua10bo63oFDiVcerFL6G/jJx2Pd1zhn2m9az62Zp2WnthB9SxqUPUGdMp7ArqGjzncq6n27K7/Ve9X4+fVz5ffUHyQslFwsX8i0uXci7NXU65PHsl+spEz8aex1cDr97v9eoduOZ67cZ1x+tX+zh9l25Y3Th/0+LmuVvsW523TW939Jv0t/9m8lv7gOlAxx2zO113ze92D64ZvDhkM3Tlnv296/d5928Prx0eHPEbeTAaMip8EPFg+mHCw9ePMh8tPN4+hh4rfCL+pPyp/NP637V+bxOaCi+M24/3P/N59niCP/Hyj7Q/Fifzn5Ofl08pTTVNG02fn3Gcufti3YvJlykvF2YL/pT4s+aV5quzf9n+1S8KFE2+Frxe+rv4jcyb42+N3/bMec49fZf4bmG+8L3M+xMf2B/6PgZ8nFrIWsQuVnzS+tT92fXz2FLi0tI/QiyQvpTNDAsAAAAJcEhZcwAACxMAAAsTAQCanBgAAAAddEVYdFNvZnR3YXJlAEdQTCBHaG9zdHNjcmlwdCA5LjI1wZk/DQAAIABJREFUeJzt3b1vG9eax/EzFxfYTpJ704C92HVBVnkpRJf2AmIr3JgppcUVU6qR4ttoi7i5DtW4jLVAXEbMXjVb0ICdbk1h1/GtyMILxAuYruNQ/8BsMdZkNO8vZ+Y5M/P9IAhoipw5c0iRPz1zzhnLtm0FAAAASPiDdAMAAADQXoRRAAAAiCGMAgAAQAxhFAAAAGIIowAAABBDGAUAAIAYwigAAADEEEYBAAAghjAKAAAAMYRRAAAAiCGMAgAAQAxhFAAAAGIIowAAABBDGAUAAIAYwigAAADEEEYBAAAghjAKAAAAMYRRAAAAiCGMAgAAQAxhFAAAAGIIowAAABBDGAUAAIAYwigAAADEEEYBAAAghjAKAAAAMYRRAAAAiCGMAgAAQAxhFAAAAGIIowAAABBDGAUAAIAYwigAAADEEEYBAAAghjAKAAAAMYRRAAAAiCGMAgAAQAxhFAAAAGIIowAAABBDGAUAAIAYwigAAADEEEYBAAAghjAKAAAAMYRRAAAAiCGMAgAAQAxhFAAAAGIIowAAABBDGAUAAIAYwigAAADEEEYBAAAghjAKAAAAMYRRAAAAiCGMAgAAQAxhFAAAAGIIowAAABBDGAUAAIAYwigAAADEEEYBAAAghjAKAAAAMYRRAAAAiCGMAgAAQAxhFAAAAGIIowAAABBDGAUAAIAYwigAAADEEEYBAAAghjAKAAAAMYRRAAAAiCGMAgAAQAxhFAAAAGIIowAAABBDGAUAAIAYwigAAADEEEYBAAAghjAKAAAAMYRRAAAAiCGMAgAAQAxhFAAAAGIIowAAABBDGAUAAIAYwigAAADEEEYBAAAghjAKAAAAMYRRAAAAiCGMAgAAQAxhFAAAAGIIowAAABBDGAUAAIAYwigAAADEEEYBAAAghjAKAAAAMYRRAAAAiCGMAgAAQAxhFAAAAGIIowAAABBDGAUAAIAYwigAAADEEEYBAAAghjAKAAAAMYRRAAAAiCGMAgAAQAxhFAAAAGIIowAAABBDGAUAAIAYwigAAADEEEYBAAAghjAKAAAAMYRRAAAAiCGMAgAAQAxhFAAAAGIIowAAABBDGAUAAIAYwigAAADEEEYBAAAghjAKAAAAMYRRAAAAiCGMAgAAQAxhFAAAAGIIowAAABBDGAUAAIAYwigAAADEEEYBAAAghjAKAAAAMYRRAAAAiCGMAgAAQAxhFJG++OILy7IePHiglHrw4MGLFy+ybuHDhw/37t17/fq1888HDx5cu3bt2rVrzjZD74m6Uyn14sULy7KcZuRrj7t957neDZYqx46iOiHqR2/fvr137969e/ei7knT+TE7BQCgLDYQQSm1t7f3yy+/OLe//vrr9M/95Zdfvv7665s3byqlnj9/btv28+fPlVJ//etf//rXvzp3Bu8JfZi7zV9//fX58+e//vprjvb4jst5rneDpcq6I+fYv/vuu6+//lopNZlM3B9F9c+f/vSnjY0N72+09540nR/T8wAAlIcwCtu27V9//fXu3bvO3yd379799ddfnRh08+ZNNxLdvHnz559/tm3b+acb6ZwQ8/XXX29sbLgbdMLo3t6eG2ucZzk/dR4fvCf0Ye42nR09f/7c2x5fY3ztiTmun3/+2d2gs18nunl/5LTfyYWhvRTsSeeJf/rTn5RSGxsbzhO9O0rju+++29vbc7rR3bsjtH+eP39+8+ZN749896Tp/JieBwCgPJymh1JKPXny5KeffppMJpPJ5Keffvrxxx/v37+vlHLO87q3b9269ejRo2+//dZ55Lfffvvo0SNnCz/++ONoNHI36DzSeaLj9evXTqFUKfXJJ5+8fv06eE/ow4Ktddvzn//5n6GNcdsTc1y3bt1yH/zixYtvv/32iy++eP78+bVr19xT2x8+fHj+/Pknn3zy1VdfBXvpyZMnUf359u3b58+f37p168GDBx8+fAg+4MWLFw88fJtyWv7ixYsvvvjik08++eKLL2K60X18cAuZOj9NzwMAoB1hFEop9eDBg59//vnnn3/+7rvvlFK//PLLp59+qpS6du3arVu33NvXrl1zSm7Ogzc2Npyan1Lq0aNH3iwY6tq1a4k3ou70ctvzX//1X6GNcdsTc1zejTtPfPLkyb179/7yl7/89ttvf//735VSX331lXfYpfMUZ5uTySRmYOWjR4/cTYWmurdv3772cCqgQbdu3fr73//uG2zq658HDx58+umn3tGiwXuCz0p5AwCAshFGoZRSjx49+uyzz5RSTgkwpdFo5NY+E+PLp59+6iQ8pdRPP/3kZFzfPaEPy9EYtz35jitmF8+fP793796TJ0/u37+fOMsnpk9Go9ELD1+OdxLqvXv3fvzxx5s3b04mE/dHwf55/fr1f/zHf1iW9e233yqlnBljvnvSdH7ungcAoAjCKJRSyql3/uUvf3n79m38Iz/99NMPHz785S9/cR4ceg461D/+4z8qpX788UenzvfZZ58F7wl9WJHGpDwuZ6fOee3JZLKxsfFP//RPwYc9ePDgu+++e/To0YcPH27evPn27du3b98+ePAguHEnPjrlWKcWG9yU5eGrYk4mk88+++z169cvXrz4v//7v1u3brk7CvbPo0ePnOlHzghXpyTsuydN52fqeQAAtJEetAojuLU3Z4KOM3lFeWaxOJN7nCnhn3zyifPgTz75xJknriJm5/h+5J39E3VP1J2+rTntOTs78zXG97DQ43KPxftId66SUuq7777zzW1yflOcwQDOY5xRAcFjd+5xWxU1gcmZhuXyTlGybfuXX37xDt/0dXJU/3hnIAXvSdP5UVsGAKA8lm3bWfMrGsmpLIaW8ZRSHz58eP369aeffuqcen7x4sW1a9eiHhzDGUDpfWLwnqg7o9oT35jgcfmOxbtT55HxQw6cwqH7sCdPnty6dcstbb548eJf/uVfnOx469atIie7Y44rsX9Cpen8fFsGACA3wiiQ3+vXryeTiXfEpxtGfWfeAQBAKMIooNPbt2+/++67r776iglAAACkQRgFAACAGGbTAwAAQAxhFAAAAGIIowAAABBDGC1qtVotl0vpVgBAfsvlcrVaSbcCQEsRRouaTCY//PCDdCsAIKfVarW9vb2+vi7dEAAtRRgtZHt7+/j4WLoVAJDf8fHx2toalVEAUgijhZydnf35z3+WbgUA5HRycnLjxo1Op0NlFIAUwqgZLOvjf6G3AZ+6vzG8b3K9D0YW8/l8Pp/v7e1JNwRAq/1RugENd/v27f/93/91/5ntEgN8++pQdicmvKLui1jG1SUsK3Kzpe4XTXFycnL9+vXxeLxcLo+OjkajUafT8T7g9u3bUm0DWuXNmzfSTZBEGC3dlQAajAjexElucHojbz+E5s6y+zQq7ObZby0SZLHX6AqTD7Md9vb2Li4ulFLn5+eDwWBtbS34mDK+Iy0r88X/LGXZyX/6xW026qfB+xPv8f4zx7GkkaeLUjyFLiqji1LeGdNF/NVHGDWG+x4NzanFxVTRDJH9eH1PyHp47tOL9EvUc63CW47da7V11oIPDsr0XMPft3XW6/WcG2tra/1+X7YxCSxl2cl5FEAdEUaLOjw8LH0fxWtR7gjUmn+vJ6ZPLfkyN+/ePzbArQ1cfaTm5gVLqllf6HzvkFqUcpHC06dPq9xdjjqZbdtW0qCb+M1G/TR4f+I9ZdT5Etug5Sl0UcHHhP405Z3Vd1GNEEarFXzzxbwdnR/VYuSollASenrI95DoB+TYsV3+iFLvvrys0B/5esDMqGdUYwBpRIpEdBESEUbbwbbrVRYNZkSjmp47AYc+UbKa6xbd6/PeAAA0DGHUeLpSQvGqW/mFumAsy1q5zNQyK3Aj6xZU6Kn52IcFhR5stmYUfEVy/6Gi5Z1g5FsRAFAZwmitaJzIrFfKJvnafxkprMt7gltJ3K55fZFKaO70l0ttW9X2AAEASIkwWh9uNahWJ9xDRLS/yCHlm7pue56YaQvFZ98nVnyD5dLwnWatEYb+PVP87VS8VFnrtzQAoADCaH044z6VxBKPuiuyweRnJ6aupK3lE1qVzLcFLQ+LeWLRtaK8M+HcvwdMiIAmtAEAIIcwWitFql+56a7I/r7Ob/TZef9T0j2sdjIdka+k+nsfmrDeQsXNIL8CQIMQRhtHeybQUZH9mCZzDYJMc0bbunxkZSGl4BR4nYXSrK9LSUkueLEGIiMAIAXCaPmMnXWUXvEYGnG/90d17iAB3kJpqq6rxTuQIAsA7UMYbRwzlhTNUTiMekrKIBssoFZwcr/40WlpgInDGKTfgQCAuiCMNpFoDpAKRlF7LOncvVFRK3SGvlgLC16Fy4QhsDCSrWxLcXl6oIEIo+VrTYmopNWRSr1AUWVXPwqt3QZbkmbZ/KjHpFw0ynSt+X0BADgIow0iN95OSzU0/dO1HJ4V+GcFeTSRxmZUOZ0rogUMAAUAJCOMigr9tq7bhCf50BMhsRjpY+yBOOx0LfRfxsnsgwIAgDBqmFpdZqmarFP8ckeJjynv1Hb6YQApm5ppdyHTm6r8UyfrAFAqqQDQSoRRUXqna1T4Fd6Akpsdcbt5wiNpkXdL+tRIrAQAJCGMSgteJdzs0/QNiKHtVIOXjKn0ANBKhFHzmBpDlfGjKg1UcXcl7862lfjrGPPnlsFvfgBASQijuCoYFC6LVb9fUx41ZxescGs5xQ8AgFJKqT9INwAm8U6fCjAtiVqX/9WReOPtkmdu/c6yPv73ccemvY8AAMKojFbO5G/lwIjV38/nGthaHXIXCFM+0YQhtjFt+LhcVMXvyYa+lwAA+RBGpZmWTT3N+JhEDWlYmXKPoUz5xODDKqpKptiLQa+u2VP3AAAlIYwaqcpv5bA0bEI9L5HG5okcae5ybL7nxrmsiJcysSnl2zjlCrum/fGGanF5eqCRCKOVSz+JWGjd+xZOmc9xyCkfX7Ani6fPtE+0bZX6Ik+lMH5RMwBASQij1Qp+3SZ+9VZYCqpFQdQrJqsVOfMetc3KRJ1bL+PkfvCiTZJ5FADQPoRRCW7JMypoVvmtHJyrVGfW1dv21RtBlY3d1KW818j9U6Towk+lIrACQOMQRqFUzZNofMuty//H5FFjpW9b7oJuaBxPDPEAAOhCGK1WjpPylZydr2nmCDbbLnyeveKuCK1BFm9Dpk6IKg9LnrIHALQGYVSU9DlHQ87G6h2maV+thpqcqDSOEIgKlPHHntjzRp+yBwA0AmG0vSqIaN54JJhmqt+1SPy1r97WlXTlT9kzyx4AGo0wapjgN2762fRZ5t0bWyxsgMRRqq6sL4FgkVK+wCy00hkAoGyE0TaqPlXE7057Y5qdWVK+fLoe430wf8MAALT7g3QDoFtS9ajKPGFf/peJFT38sRbLMNme/zeMzEpYzluasigANBSVUeOl/w5uxNl5y3PDDru/FtJ0co5pW0WmIuVmxKr4JFEAaC4qo21hRWQIqz4VRxiidlcKQJM4l6eXbgUAnaiMtoJRBdH4WTh1XJo+X1Wy1pEuQ33U8hwlBU4AQABhtPmMSqIuM1uVnpbT4hp7oLzOTHMl1ci9W5b/n+RRAMBVnKY3mGX5v8vj7w/7UWLmyzfHyN24CbU98ZEGUvEq8cDL7hn3wCN34YueJFEAQABh1HgxuTPpnlKrj7kjTpH4axrbc6MZR6QfARQAEIvT9PWR8To0dT8PbriS+lbLNeX1im+StzjK+w0AkANhtHxZLox0RejjY0bdee6PSgblLQAkSO+x1OtS7NqvJpD7HRI3pcm2Iwv8AIDW4zS9SWIGg2baTH2ylLHMiU6GtKToO8q2OV8PAAhFZdQ8UbVPp7xUyQWWguWxlPXCRlZeyxDVn3bYw2T/usi0Jn9kU3OfHwAANB1htHwav30Lbyr3rHn3hnfKTu5NmZxHMrUtNHwnHmOw2GlOiLeLJeBKr89EwAWARiCM1lbYN3EwB4isx56mgFplfCgY9fI9Xcsx1jFkhefRTMNGSZkA0CaEUZO4X9jxp+MjvtT1Jjz76u2UFc30DWjSwNYcc30qUDCCF1zMPzyP5hO/jgSBFQBqjjDaEGVnO10bF7kGZnCn1tWfFtly4p1pQqFReco35CDYtjR/mWg7Xx9TTyWGtpVzeXrbrN8bAPkRRo2R/hLebgH18mEx3/oGflpXOqzQs1OR5+qSo8xZPPfHDx5Nc3GvnC80KRMA2oQwap4038QplhQ1mZaBrQUbUMFUqtxZUFeHiL8xNPzhwRqlANB0hFHDZKwJxX/Tx2Qac6ZvSwmd/x61plXwKbn3lXsLUsMbChbdM+TRqHlLFEoBoNEIo8bI+43rq/BlLfhVn0qNXd3J8PqbIesPRK1LZeALCgCoBcJojcWEp5Tj+QQVX7JUVRKDYk7oJ+b4NA0r44+BItvM+q5I83iRUcIAgLogjNZVyhCQ5qc1DQqhS/Hnk7KjWiiYa3MvFJX8MnE6HgBaiTBaSzGnSvNdNad65A6HeIk6ivMCJbYt6+x+XncAgA9htH6sq0HB2DQTxWmwlnJmjgPXeFpc49qrerdcfB0rS99SrIVwKSYAaAHCaI3Zsf8MMm0GvWydTO/eTetbXfQeDsVRAEAQYbRmgmVRVZ9v9/gZV45Mx5JvTEK9Csnp1WJWe7Y8yiKjANAChNE68SXRfHxbYMBoqaoPiFpKjxUsUJAhjwIAGo0w2iImf6tXWbAMXe4+UXx8j7pgpnvD5M4HaofL0wNNQhg13uVl6LUszOnyhT8ThjwW3HWOQwiNvxprmdVcdNS7u7pg8CgAwEUYNdvlgLkyvrlTTnjSst+6jPvTXsskbzmC76VUedQ7YJTz9QDQUITRSkStUJN65Zqyv4djzpLrzcFGBYpgHoqpZZrQ8lpMUQqqy58iAAARhFGz2bYq+YRmlSfoC172M34Lua8MpGU7FWhepEsujrqz6SmLAkBzEUYrlOsLVXBoncYT9A2IEiYMq62pmLp7qjwKAGi0P0g3oDbm8/lqtcr5ZNsO/0517nd+ZFkNXlKxsQemlLq8XlHWY8z6LDvL5V7zNak8+S5UC71ms1n+DzEAKA2V0VR2dnY6nc5isTg4OOj3++79n3/+ebfbVUp1u93Dw8P8O3BjqGV5Y2sFZdH4NYkMCRDpl14ypMFIL7I4ygl6fVar1e7u7ubm5tHR0cHBwWAwcH+k7UMMAPIijCabTqedTufhw4fL5fLo6MgNo8vlstvtPn36tOgOIgqiZSRRqcRW9h4rWNEzfrhqZTXIlH8nBJtU8UufaXchebS5ZwlETCaTwWCwt7f35ZdfHh0duWFU24cYABRAGE02n897vZ5SqtPpnJ+fu/cvl8v19fWjo6O1tbXRaLS+vi7XxmwSE5vGdFXZEpsi1zpy957vwqSliqp5G305qNAMSom0sL29PefGfD7vdDru/fX9EAPQJIwZTcX9+N7c3HTvXK1W169fHwwGGxsb+/v7Uc+1PNLvUVdoyD12sHaD/BIbbNowyiLq9dKk4f/7xzuc2kGttLDxeHx8fOz8ae1I8yFmXVVVY4GG49fKy7KpNyQZj8c3btwYDodKqdu3b7958yb4mO3t7bOzs+D9UY9PpDGMOniZtXdFsBYbugtegpQspezQlXepjOqzWq3u3bv36tWr4I9CP8Ryf4JVgMuBojFM/kWrBpXRZL1e7927d+pyfJV7/8nJyWw207Ybz2x67edS+cAuQ9nF4yaVctOwlbKCBVEVvRgFUjs6OppOp0op34l4zR9i1XIuTy/dCgAaMGY02WAwODk5GY/H5+fnztCr2Wy2u7v7t7/9bX9/f2tra7FYOHVTDSwrX7E6qvzGd7ir7AsH0NVacM36MgyHw/39/fl8vlgsRqOR8wn25s2bzc1N/R9iAJARp+nTms1mnU7HO/Y//n5Hhtq7WxbVGkZRtvJ6XuQ1LXWniamdt3F5VqvVYrHI+iFm+NlDztSjGQz/RasAYbRcWd9huctCdfwWb0ZBMdjzuY/Le8ZRqlsqyNYxG6/j27jZDP+OJIyiGQz/RasAp+kbohafx96UVruhXlF/J9iBh2XdbOh2pFS5Zmro3lUN3xsAgCKYwGSQZo+Wq3XCsDz/L48duCFI+8GGrjAQ3Gk1XQ0AMAeV0XYx58y4CW3QSNfJ+souVaX3CqspB4NGlV2JngDQZlRGTVFlWVTku9+u4UL6Ltvz/zTS97CB3ZK1Md5rscavRWVdvR18pLc2TEIFgJYgjLaRUdGnLhIvoNok5R1OYvI2LZoDAMpGGC2Ts1pTiit9FSmLZloXvfpv+vas3O72rWlzkrzSF2IzPTL93kOfYmB5GABQGcaMVssbTHUsqtWGkKdF8cGy9Z1eVkE4Ttxy8AFpnlLfPgcApEdltEzp4qbvGzdfKZFUGqN452Sa4l3H/KS9gF1RRdy5iG6Kkw8AAGNRGS1fmZcVMH+eR2OuHpSJeAOq5HsHajx2iqOI51yennXvgbojjFarhGDqW0tI8eUdULxDDE9F8S+9sc3OzK2Act04AGgQwqiwYMRJ+TVrVPQ0Z/lSVVozDDm6kmg/Ol0bTPgzgFQKAPVHGEVR5owTaEkwyf13SNa/GVI+vopuZ1QoADQXYVRS8TO/dY9fRtV36ytT76V815mT/mzda1AAAIxCGDVaTFbLfZkcO9396ZEO6iXxop3xjwEAQC/CaI35okN9A4QJLRfPYekbUM2fDSa8KB/ZtuXURymLAkATsc6oGJNnZ1fGwOvumHN6Og0ty3maf4ksWymLJAoADUVl1Gi+xfBVCVUxvuGD2vl3QjuPGgAgjjAqI+qL33eu1g2goYWrtkWHUheQ8nVy9clM/NUUbEChMRKtXHz09PT03bt3N27cWFtb6/f76+vr0i0CgPw4TW+uUs+cVnS1xhKEtlnL4XjHDNQl19g6hjpoGSxR9tvJ/IuNVWY8Hk+nU6XUu3fv3r9/v7u7K90iACiEMFoPwYvXt1bZMdHAYaw1Um4e5Ur0SimlJpPJ48eP79y5o5Ta29vrdruz2Uy6UQCQH6fpBcScArav3g5+61qeGwZmJr1n0rMuF2Bgh1RMak2A4jXpAk9u3cu+trbm/edyufTd0ypcnh5oAMKo0cr7fC31kzt9UNYVnnxDbCGFzq/AaDTa3d3tdrvL5XJ/f18p1ev1pBsFAPkRRmsjzVrliY+sl3zDBA2pGQsuXKpxjzH5PniA+fab4U8IZ8HRXHtpkuFw2O12p9Npt9u9cePGcDiUbhEAFEIYrVqzv021H1qD+6oMwbXAlI4+THzT5t5F1j827Kb/BqUxn88vLi6cMaNKqdls1u/3ZZsEAEUQRpumXt/TelubsoBXpXq9HDGsq4XqxhxXHS0Wi3fv3jm3z8/PndWdZJsEAEUQRmsj5mwmycBAvhfF9/JVNry1SM3SXVjA8vy/+Ja1b6RtfOfl9/f3V6sVS40CqC+WdqpU8TOMla1qE7+AVH2Xl5JdYDU01ZWkyBpVie2s7zq1zXP9+vXFYiHdCgDIj8oo4sSflrWulvq0VLm0nwjWW3tr53nq4MFqj6HpO5Zho+Px2Js+F4vFaDQSbA8AFEQYLZ+myxWKf/vGxAW3kCbeyPSCl1ottfHBFWRN7ivzW9jOq4A6BoOBO3tJKdXtdjlHD6DWCKMI5x0sGP+YevEdV77Z3GnEjwq1kx5QgfgUXtL0eRSxs7MTev/Tp08rbQcAaEUYLZ9tO1WcsktNZVT40lwpSmktpEmlnPT7zdTCxG4p712hd2knXZuKkmnLrY3CnI4H0EiE0UpoOpNYMBOkfHpUxa5epbLEwqRK0SHVnMHXVSLN11qpF67om9lu3+l5pVjCCUAjEUabI9/1iqKYPmQwlvcsfMxR6O2xqF3E/0i8AUVkuj4TNHLWvXduT6fTwWDQ5pzK5emBuiOMNkprP4xzR58Gn3euYLSGo8jfLa19xxYxHo/Pz8+9a4seHBzINgkAiiCMVkRLoTH9FkJLVimf7lbscrTZkHpY1ulBUc2u7CgqfntkUtJrWnBrNZjvX5rJZPLq1avT01Ol1HA4jJrVBAB1waL3TVb8FHC9VjX3LZ/UzqRSJTu6n2N+BC06nY5zUdBut8ui9wBqjTCKOGaGieDlf2JyjxX2eJ8yDpMLFFXEsj7+1xr379/f3t7udrvPnj0bj8eTyaTT6Ug3CgDyI4w2k+35f+4t5Hh6fethJTU7MSKZ32Pmt7Btvvzyy4ODg/X19cePHyulHj9+TBgFUGuMGa2CyOA22fRgyOBRcbU+/Hq8iJfr+LbHkydPZrPZl19+ubW1dXh4KN0cACiKMIr60TVRvdSwZWyAq0fETOH3OUwtW3L04cOHy+Xy2bNn+/v7nU6n3+8Ph0PpRgFAfpymB/JLHI0qsik0XqfT2dvb++abb1ar1b/9279JNwcACqEyihLVomCVb/mqWhxaqPRL/VdzjKGV2rp3cqnm8/l0On327Fmv1xsMBs7IUQCoL8IoSuG9vJCxkaJIw6zLaT26apkVXA7Kuy/tos7+F1nEoLUricabTqcbGxtPnz5l3hKAZiCMGqcxQ/pacsa5sgsdAQ4mLQFoGMJo6QSrO+Xl2jRbLljqK1JVrSDQm3m5I8NL0aGi1sxP+Vyqp1Bcnh6oOcKooUz7TM2XcoofRXlRI3duM+2lCTJzKTFdTapj4AYAxCCMGqfI5eBL+p62PDfc5fRrfRY+tPEmDJAwoWN9/WBCt7iCb0UAQN0RRgXk+3ZPHzS9kbE8ZUeBgpePqq8aHXjK0RrQZWdn5+DgYDqd+i5G//TpU6EWAYAGhNGaSVMQCi6RoyQiTmh6NqTMFlOAbFJ4ytfboZ1jSMwtr/xfC6PRqNPpDAaDO3fuSLcFALQhjNae+/Vs1IlL8XPNiYJ9ZVrvic97U6lfx6ypt8gBmvMyVa/f7yul1tfXT05OXr58qZTqdruj0Ui6XQCVM9bTAAAgAElEQVRQCFdgKpNlqYjQY+eaDJTjWaid9Pkv8ZEF3y1ptp/jPekd92n+Hy0GGo/H8/l8NBo5MXR7e1u6RQBQCJXR8llW2dfONnDEXtR+zQnTlVUfy9hR+gyXqVQZ9XhzXjUopSaTyatXr5zb/X7//fv3s9nMKZoCmbAeFgxBZRSGKvVa7VbgRklybN+tNUoVDouX7eNfO+8ji38NmrD+QMW63e5yuXT/uVwuuRQTcrNa9wsEE1EZLZNtq9u3yy6LGsj8KSY1nQfj/dKoZjWD8gYix2+2hREzvW63u729ff/+faXU+fm5UuqHH35QSt25c4f6KOJRCoWZCKMoi1ETqoKqaZvJPeBl4LCKunRd9W7cuOEkUaXU5uambGNQO24edW6UURkl8iIrwij00L5mU2WfZFWWG9OLGrtZi4JuytVw0z9Y+9NrbTgcSjcBrZAvU3LeHzkwZhT6VTnxv9ShpbLc47KuHmYL4xeQqKQiX5PkCJd0KapBGK1IvTJTvVrbSJVNsQLQQonBlFPtqBKn6avW1OJW2UdU5Bqq8c9q3mtRCwW7veWv2unp6bt3727cuLG2ttbv99fX16VbBCh1dRwqcTaU0zn0jA+VUYSo0W9JU68FYAf+Axzj8Xg6nSql3r179/79+93dXekWoQZEoiFn+ZESYbQioWsrmnk23LqcCG9CAPL1jwlNkmJ7/o82m0wmjx8/di5Pv7e31+12Z7OZdKNQJ9Qsq2cpy/nP/af3R0KNMghhtDoxbzcD34lGBWV3laisTcpaVixy1Fb5nRZ6IIl7DP1p2U0N3X6OPRryDjTK2tqa95/L5dJ3D1BQ7nn07rPIu16+uOkurRX8UWsxZlRG8N1n+KqcsuiZRFHvHyv2p76HldrPOT5xUza+bUaj0e7urnMdpv39faVUr9eTbhSqk2PQofnR0PwWFhR6dCRRF5XRoubz+Wq1SvNI9zRr6HnnrNW78hg1SDE4vKHihkX1drDyp/ESlxXwNrLI26nKCro5pXpZw+Hwm2++WVtb63a7/X7/6dOnBTc4m81SfoiheYosiZX1ie3MXqE521a2c6d7o+UIo4Xs7Oycnp7u7u6mHLPli1PuPzMl0YoZ9eGh61c2KkIFr7GeiUiIj//jISbBl9raYKty/JET2nij3pAixuNxr9c7PDw8PDwcDoenp6fz+Tzfplar1fb29suXL7e3t51JUWi8qPTjHc5YMB55k5YvgIbmUQNDasEm+UaIRiGGujhNn990Ou10Og8fPlwul0dHR1FXhQ49AcobECUJrmaV+GYr492o6wpJwae3+XdnPp8fHx8vFovFYuHeuVgszs7O8m1wMpkMBoO9vb0vv/zy6OhoMBhoailaJ83CpZkeL8LAWNwShNH85vO5M1Sr0+mcn5+XtyPv93rB390mXUSx+DDHIilNpAODhyz+gla53wa8aQvq9XpPnz4dj8eHh4daNri3t+fcmM/nnU5HyzZRqhwXlE8T+3KEMLclaTZuYO5MlL7ZofOTuCRYJoTRQtyP783NzfhHWpallLLtUn4hy178yOQZJPnaFnyKeKpL5L0mU3xKbtvnn7EvWUkODw/n8/nFxYXzz+l0OhgMos7MpDEej589ezYajUJ/6nx2uUr6EDNcTeNUDr4jdRNVytwZs1l3I8FVjbyzy03oZ7c9WQ/ZPcA0qwp4f7P++Z//uUB7m4AwWshyuXRuxFRGP74NW/kJjtyCRVA7dTFY41std0aPmdqfY2vxTP5jqQzj8fj8/Hy1WrkXXjo4OCiywcPDw9FodO/eveFwGPxpjdJnwcAUhfpWvJQh0oSUGS/qzRN6v3unr1adviDq/c26fft2zkY3BROY8uv1eu/evVNKLZfLbrdb3o7cmR/FZxNnmkSid3p4aOOLzMU2ata/XlEdEjzeBncCokwmk7Ozs9FoNBwOz87OiiwyenR05Mxb4oKitZAjaptQ0/WVWt3bRqX8YGNiVgMNrl3vnRbmvSHe+XVBZTS/wWBwcnLiVCncoVdlSPx91f5mD61gxZe1Uv7UkCKWr0sbk+eiSqdFCpwGfV3Ao9PpvHz5UinV7XYXi0W+0/TD4XB/f38+ny8Wi6jT9DCQ9hqwN3iF7itq+lFiITC0kaaNp0zTmKgKKHFTC6tGp2DMNJvNOp1O1Nj/27dvv3nzpuAugsmpbOnDaMxvcNTDQlfqqfhdqKVLDRxmGtUkA5uqkSF/4VTG+QP4+++/397e3tracgqluacfrVarxWIR9SGm5ROsSuWdpjchc+S4xFGmDok5T62uJuAiF1sK3YjypL3quzrlUgDlhdHa/aJpRxgtl94wWtlLFZM77au3U4ZRLQ3Qu7Xis74MTHjNC6NpWt62MKqUms1m/X5/Pp9Pp9M7d+4Umb0Ur3bfkdqjjGBCCm1J1D9TPivfg/WG0fjdVdzV3nCZ5mHeB2tsau1+0bTjNH0NVP8RGLrHlCd/DTr1Eis+3CQ+xkxRDa7dgSDGarVaLpfj8dj558uXL8sLo6g1EzJ0DJPb5hOaPmvUfvMRRmvJnLQUbIB4kyrTniOFUfb399fX169fvy7dEDRNTLqKKcrWKJMFq8vOjZSHUKMjrR3CaJ2kLDqWveyoujxBX97vpd4tV//5Yc5fC3VHBwYtFotXr15Jt6IVDK8sVqYZnRA6I74Zh9YALO2EnIKjEq2rt+tyvj6oduslFezwNE+Pekzu/WZ6YnDvtX6DFXT//v2TkxPpVrSL+Ozvmi7qZCDvVePpH3MQRmsscYCgxt+z+LxS6od0faNt8f7PeuAV91Lx3RXcQh3fFVocHx/f9pBujkHEU2PjlZ1xQ1eS0rjx0NsQx2n6ShU8e5v+WfySJSr1TLreRQCyPqvsd1f1ooaFmNzm8vz7v//7//zP/7BMPZBe+oVRIYIwCg3s2H8argHjO33rbeV7evrHaOmxHMvvF3l6k2xubko3oRU4zS1LS//7Vqp38LIaiDAqoI6/B7ranGNyVR27S4uCcS1R8csNtHCBTxN0Op179+55L0H89OlTueagdC0cMOqerC9yFN4MWuveaAPCaKWi1pBH9bzVxJQqWKagejkCZfzFDlC24XA4GAykWwFhIkvxV7M77Zct0Lg1lIQwKsO0a7W3FrkqHxPetK194Xq9nnQTYJy6l0KDtAzuNOTqWUhEGK2HBoxrdBjY/kzTwppXzy7pWKrsqCa9HPF2dnYODg6m0+lisfDez2l67cy8QjrSoyBaL4RRGY2MNW3QpNervGPR8iXQmD/ANBqNRp1OZzAY3LlzR7otqEiahBpczp1Qqy5P99MVtcA6o6WL+mLOsay6sb9S9V0KVBxdl1u9LkygRb/fX19fd1a87/f7/X5/bW3t+PhYul3QgEpeqHzdQhyvHSqj9WDgb1XdC1dmtt/MVmWlpfFphvO289v74ODg6Ojo5cuXSqnz8/NvvvlGukWtwOhDoDyEUeB31WdB02ZQBRuTsivKGHbibi3+dWnbLMBOp3NwcLC7u6uUevz4MfOZfIqnRkKnUXg52oAwCg0M/JxImY0KLhevhYG9l1vbcqGIk5OT6XT6/fffr6+vOyXSw8ND6UYhWUnLhfrWdW9SdMsxob5Jh98qjBlFTrbBg/bSfHoFG29ChdKcXi3jmqLae9hOsdPmuX79+tnZWa/X63Q6zKNvhuIRqtlXucx0aOTROqIyWjrbspSd+RejAWMH61gky9TgBrxGQXbE7UxPDGV5/t+YqXsiBoPBbDZzxoxubGyMRiPpFkEPhqWGSp+zGxzHG4/KaCUsmd8QwWna3uRhLHcmuzn1yEwM794K1PEPnuLG4/Hx8fHGxkav13v37p0zeBS6lBQHQy+SXlCzq6EOonlLUBltvnZ+YTeeG/fr8uKaNlWrviaTyYsXL9bX15VSg8Fgf39/Pp8zjal5UsZi72MorKKmCKPls+18l//Wtn99m8q005RHHXqyu+CMovR7L7gXOFKOWKDHtFhbW5NuAnQiQWpET9YUYVRAjrGG+YYnyv5SFt97kbJfmidmysrpn1KNNIHbqFGtWpYsMOqIpGxtbe3u7g4Gg7W1tfl8vlwuKYsCZPpaI4xWitOUUTItZpn+8c1Wo07gna/R4eGhM4Hpt99+6/V6BwcH0i1qjvhAk3uCUXk5ifiFZiCMiuEjxBHaD3ROvVT8ermrw7bzfTIejw8PD/v9vnRDzNWAaen52l/rQ0abMZu+Cr6pG96J2yknvAcvRYNq2Ff/M4Hlmb0U/34INttK/ZbTS1cH2pe/TW3+RVgsFrPZTLoV0ENjaG52Em38ugEtR2W0UlEfFQbWeMw8IS54aUrTWJ4bNRo5GmR48wzkzKDvdrvuPSx9j2ZzSt0xwb3uhXAQRmuDlXGyMiTiE7agV7fbffz4sXQrGqjUFUa1bxZoEsKosCpnx+dIRXyCGijltY58j9R15UwTqs7umNF26vV6JycnzhWYut0uV2BqEncde/KrT8wK/5zBbwDGjFakdnXNmBF+UuMO0zNqiKd4MzR2RRkver4Xy5Cyt4jxeDyfz0ejkRNDt7e3pVuEtHyJiqJpVr7ec/5TZPf6ozIKlK55H5OCYw9M/iuoGpPJ5NWrV87tfr///v372WzG5Hq0RBkXVoU4KqMtonFGc1PV9+MtuFxDGS+TeMm5vi+QRt1ud7lcuv9cLpedTkewPc1AkdJ87gvkS6K8cA1AZRR51OVXP1MBL/Rq7yYMkUwkdZ16w7ulqbrd7vb29v3795VS5+fnSqkffvhBKXXnzh3qo6ZhudDyUB9tEsJodXzXb2zeJOs0R1T8qEX6LV/Oa95L7JI6opYvd++4ceOGk0SVUpubm7KNMVZJ695r3Cy12OLowMYgjKKutNcsTfhUy3dQaa5Tj8YYDoe+2fTr6+vSjYI2BKx4FEQbiTGjwir41NE7+V3L1vJVGXPsN/dMbetqLrSjH1CGHFuu7OvL/LUUGo/Z9GgtJ6nbyiayNwyVUTE5ElKOZ1UpTdt0TZ+qsjdy78XkF6u+qAEzm147lrsHZFEZrVTx1UYNqZmJf76asIxoI6erp2R+CxuM2fQAGobKaCs0YGBllftN3JfvCkCJE7aaFN2adCw1xWz6WuMaS0AQYbQSlqWUUraGjx4+vWqkXgMrDZz7b2CTTMBsegANQxgtn2V5b3/8Ws0eTPk+Lqh5Rco0cuQ531OqT4Qxo0IZMKqUGg6H8/n84uLC+ed0Oh0MBtREc2NkJyCOMIpWEClS1vH7rfiw5oKiLj1Qx84syXg8Pj8/X61W7opOBwcHsk1CEBkXSI8wWj5vEdS2rxRKq8LpzjJk6kzZumyRhQ5Kiqdp3pPhu3Z/g3SMe6kjZzb96empUmo4HO7s7Ei3yFAlrXuvhZmtAqQwm75ytq1su16jCRsgdLlQ5Vk104q4vyDxVTkT5+YHWxic0a9rgn+arvC9Ut7GUB/16nQ67969U0p1u93FYiHdnBZh0XVAO8IoTFF2bjMhx/ANplIUiU14pUx2//797e3tbrf77Nmz8Xg8mUxY2ik3kdIpZVHAhzDaComrV4oX8ET4yoFBRa4U5e1Pb+e3rZNjBN+TMe/D38uitv3xv7Y6PDw8ODhYX19//PixUurx48eEUdMYOzwAMBNjRmVEXUycwZ0igr2tq/+rudxrkR0ZtYArUnLmzvd6vV6vJ90WACiKyihMUf3Fh/JdkbXg7rwlUi1jUn03AACoFyqjUKrFJausq28mzqHJ+lMm5XixvCgAtBBhVEzomfp8gxTzPRENUPx1zz0yhCElqDtGdgKGIIw2hCHJIGUztLfWm8hLSudRw3yLbDA3/gIBjEXGBbIijErSnm/SKJJjRK4M6RW6Xyvitnbajzpff6Y/xvSvtcnfnP7fkdYveg8ADcMEptrLNOmn7Gku4ulAvAHmSDm3yS0nZ51QVf2EMyCrkhaoZ917QC8qo8J0FUfLzgSJS3Km/Gn6x6RXfNxtSpWdHI/fUcr3jB29nawHUrwinrvrIo+UsigANAVhtF3MT70uQ0bBuryFxlJXIU1Tb0nZgMSRu+5qUwVzdvyL5avRBtf/z7ZfYih0YGQnYA7CaJ2k//JOM9QSmRRPbKaxI26bgxWdUDtkXCAHwqg8kWlMXmkyLh+uqqpOqGAvWXdRsEmlrHLANCYAaArCaPPV9LvaO7FGlX8Upo0KqJd8K3mlHEUAAGg2wqgR0k9JSb9B7YhrMJFlURwFgFpjaSfUYI0ek9tmiBxrM1UmdDp/0aYSQFEAIzsBo1AZNYX4yNHyFCmpVtYhjez5MuQbNeF9fJoY2tTfBTQbGRfIhzCKON4cKfIR25ixAc04ELMqrxRHkYKzQL32jFjSZoF2IowaRLw4WmpgCtbGWvspXkY/VzbTv3geTbMwWWvfGwDQQoTRRtGVcnxpQMvKPol7KajUgNuMumZWoV1awbIGqXbB0k4A0BRMYEprPp+vVquy96Kl8qRR8bPzBSfWZJ1cZVTveQUPpC4ZytguRVbz+Xy5XEq3Qh6n1wHTUBlNZWdnp9PpLBaLg4ODfr/v3v/55593u12lVLfbPTw8lGvgFbk/ZSuLwrq+ByouWIbOwjF5VlYmJlR/OUFfktVqtbu72+12l8ul78PKzA+xOiLjArkRRpNNp9NOp/Pw4cPlcnl0dOSGUedj/enTp3p3V2TkqIGXnjcw4OZIXe38hjH6qDk7n8VkMtnc3HSy5t27d93QWdKHGABkQhhNNp/Pe72eUqrT6Zyfn7v3L5fL9fX1o6OjtbW10Wi0vr6ua48lzWSSKn2RGpAGZdHybG1tOTd8Y43K+xADgPQYM5pKp9Nxbmxubrp3rlar69evDwaDjY2N/f39qOdaHqU31LfrYlVJ6+p/BiqyXH/xWVmGXykgq+YdEVydTqfT6cxms93d3dFo5N6f5kPMuqqqJgMNx6+VF5VRv9PT0+l06v5zMBgopdxR/97K6GAwcH7a7/e9T/Gxc51P1FgczTE1Pv43o+4LM9W35eURHzBKWbRs4/H4/fv3jx8/dv+0Vuk+xPJ9ghlL48hOlhpFEd7frNu3bwu2xASEUb/hcDgcDr33TKfT+XyuLsdXufefnJx0u13vfCbDFfnIDL2iYyM/g8uLZeKBL43gy1qLZiPe6empk0R999fuQyyfClIjqRQogjCabDAYnJycjMfj8/Pzvb09pZRztutvf/vb/v7+1tbWYrHw5VctihdHDZzPBAQ19Q8bcziLOu3s7Dj/HI1Gu7u7b9682dzcLPVDDADSsBp2CqY8s9nMGXeV8n7H7du337x5U2S/Bn5PB0/TN6l4puVYQkcyGN5LUc0ru9kh73AWtK9WzIdY8U8wQ7iVS70lzJI2i7ZpzC9ablRG04o6k5X2DFfw+9W8b9yUscOU5paj+NFFjbg1vN+imlc8lBfcyOW2jPt9aYzGn6YHYDjCqOmKnKw3vBSXktRRNKP3DJf83maeKQA0HWG0xkqa1d6YSwpJMaQ3zA/TkUmU2ifKUdLJdM7RAwURRqsS/H5N/Y1rq8v6kOcpaepFWa9+ZOYw0Mqu4aRX/J8K9V0eK9Mbo44HCACoGGG0DgqcqWxGGog5ivTZKGu8LmORV19LKqD3DWBdvV1wOGmqpwfHWAMAmoUwWife7+9mpEwRBi5Q4NBbjTbzGB15XgLO3cM8zgqm0q0Aao8wWge2XX1ZqC7f/NpP4utKhDFPr0vfBmVqeX2HIgAAqkQYlWZZqUo+tq0iBo+iytPuQeYMrjVT6CWd6CtUrNSLMDF7CSiIMCrKSZbeqmd8ynQfmTLCtl69inO6Rqlq3Jp2JFEAgM8fpBuALNwAShJNwQrcSP+UrP3L6xHFt0QDHQUA8CGMinIyZaaIadvKthkwn0Y1k73s7CMpW/LyZe0ZoFScoweMxWl6ab48qlJd9tC9LFOO09Api3/ewFTfz1rTWu4t1tqee/SeoDfT72VRxj0DADwIo0IKX2jbkKXg6zUoM5F7IMHInhjic5zi1/gKyr4ZQt8Gvib5u4VxzwAApRRhtE4K59ePz8549aY0Oyt1LKDhM3LSqPgQKu4oK3A72ID6vnYAgLIRRkWFxsr0lwnNG09TVvvaHCDskhO29r1krZSXmo+tywGjcU2iLAoAUEoRRsWEfhOnCZeVf4WnSS0VtKn65GLH/tMnzWn9shO/9guiFtmvWyJlBj0AIB5htD6CMfTyykyC3/eGxNCos8OGnOI3dgyD9oZ5Z2Wp+HcmY0YBAEopwmh1sp5ST/n4yysz5cuj6Z9icmowYSJXPEMyscP2/L+afQEAEIMwapIChaLyzoc2I09o7JzQZGn4lJ2KFz1IeDc6Ff0c73ZNc/gAAEZh0XtTBdcfTXxGHWqEZbCj11f33i/YOeVFp5hjD6qsB5KjP2kSAHCJymhVcnz7GvyF3bDlRXPIeuzt6StmLAEAMiGMNkrFk5crjh2GXDjekKSVqTca8pYw+M8zAEBuhNGasDynWC8n0X+8fVWOPJq7zNnsaKBx1pFRE5gAADAKYbQqRa7HbVlx/wzIlEfzjSMUD1VGjRMga35kWR97IM37nNlIAAClFBOYKpIUHxOemP3bupGTmYKTdSo4RmODUszUJevyv7J5d3Rld6FveMv6+B8AAB5URg3mfG07i+DkzaOJTzMwbBlS9SzYAO9RiB9L2ULeae7iTaRPAEAswmgldJ2IzHhmM2UeNXP6c3yrjGqwe9mhYJvN7NtQRf4GuHJ2PmX65Ow8AEApRRg1mndt8Kjznu4jo7aRlIfcc6xEg6y86S3VBTCrkqMBuYcOhxyv792Ye4l7AEA7EEbNFvUVnuXUp52unppYGKtsmk4tYkviC1DxMlvVy1BJ1ZtEmfkEAM1CGDVM7kn3eb+hGznVSVD8NUKNlfVtkDlnF1lNAgDQaIRRk5Q61cO2g6XN4P4MmTxUUAVHYUf/s8jeBVeJSr/HnElUL3ItADQFSzsZKfHLOzi/Pua72XmwZdnRmw1Ouwl9TJrLoFe2rlBMA3w3RDS13lzd2IOodaByLS4BADAWlVGTeC+tFPxRzJ2pz4HalmXZdmie8E3EMVZi7dCOrU1WWfot49JWsgvs50yiRbIjk58AoOmojJbJXShU43LfMZtKs33btq/GNV+xM2X503Axq8FXs/cizzW28yXnY7FYKQA0F5VRw+iYPp9my+XN9TY2S2VVpAbZvGlhAknUV/gPLnNGxRQAGoHKaJncL8uyvzWd7Zs9j17XQNIitdsG1H3zHUKRgbyhi/kDAKALldHy5UqKyRtMc2filqoaQ1le6pUdQ5nIzFalZ0QMZRV9AGg0wmhN5PvqjTmheTnnyZ23xHe7Dx1i7ruCJAoADUIYbaXACNSyLxfUtuxg2oDRfBcIbdurBgAQQRg1j9D8jMpO2fsUPM/ezsBU4uAEy1JKWXbltUcu0QQAbUUYlVBZ3IzafsSCpsVP2Xs3anKoLTt2F+zAKpdBvbK7y3eF5BJO3jdnzC8Ic+oBoCkIo20V/RUuUiKtcl+mnUN3iV8yKuRVIPMBAEpGGDWPAd/6Vc5q8l75qeBy8camzCJCe0bj63Lltc793it+kt373JhLkQEAGocwWlvll6xyzGrK1xRduSPl3sXDftSIgvirgGpvtq4/A36nccWlxO0Y8DcbAEALwmhz6UirUrOa4pm2tmjF7dG4u+SMmzhq01fRBAAgI8KoBO3f2WXORNZfPwvbRasWEso0oqCkbin6gvpOoxNDAQB5EUZNlVjX9F25u2S+UaTaa4HudqK2LFWgTX+kFYxn0EVD9GdYJwBAE8Jo/UXFghKKVXonCeVItJaRF5evuD1Fdqcz01MNBQDoQBg1VabKU4WxoIKz9maq+0CC+r1krCoFAO1AGDVYxd/BWQaealn7KbTOGrVBqTzia6Rpc6eSSV1RCQCAdP4g3YBWsizjxtu57cnSMCeoxTzBSnpA1Dn3xCcaJUdTKz4627T3W6LaNRgAkBdhVI5RX7du4SxjBc1OEUnrzvaE5mDvFDnwUjvNUmW+x0r9g8p7gj70Dens3ajfIABAXpymx6Viy5GqOo5KzMXcA7wMZ5ZtK1879Z6lJwUCAPQhjMoxZxRfmpkiwccE7vFF0iLDK1M+xYQRnEVWGChv3dZ6S5y9x8JSANAgnKYvk/N9GfzWjDr5WCX3RKfuL3W7WD6roxyrTeldoCpymITzTtP+Zitps8FdFHkAAKAmqIwiC9+1H6ODrDeP1n1RJJNdGRpRXjjTuMqS9z1DmgQAEEarUItv3MRGeh/gLgIV+yzvz4oPJw09I1/qpTK92zdhPICPtut5lh1hM13dPtPTAQCNwGn6Mgl+jyaegnfPtGY93ZnrtH5Jk+7jt1bN+lAVr0Ll7k7gSlROLTzr0I58V7HPtdYYAKCOqIxWpW2Xkwk7Xt917X33Z97D5f9rNA2oSIW1rMUKLCuk7K0CQzI0Cs6BC30MlVEAaAfCKALiQ0DhcBA8953pWRUI7kswEJU1QqCCCenxbxVv7TP0kcRQAGgHwmi1Kvt+1TLRpKrW5stbdlJZtJrWVzZu1Qga3xLUPgEASinCaHXSfOM26VR+ikOIuZpRyplDubtJ7ylvbyPjt5y4x0ozqO90ufNP5//OwNCy34cNeJ8DAAojjCJALiLkO4Ofm/bxpvma7X1WhvZoD4vBDVZYHQ9B3RQA2oHZ9EWtVqvlcqlnWzGX4S5O6xL3VtR2NO3i4+XgLcv5T8XOWNc7nz3y0LSylH9efLbZ8VHXU0h4UvZDk5rMnmW/1bxkDTafz7V9iLUPb79EdBESURktajKZ/Pbbb4eHh3o2F7PmYsrpIOKlrBL4KqZ24ZqZ3j7KevFS7Q0oKtiTRaY3mfBWRDqr1Wp3d7fb7S6Xy263q+1zDACyoDJayPb29vHxsXQrypGykpqvFJpy44E0c6V8aFlVLvCZlRX4zwN41d0AAArHSURBVFvuLbp17/jORmrwoZlkMplsbm4+fPjw6dOnz549q2y/OUplaZ4S/5ionwbvT7ynglIfXZTInC5KeWf1XVQjVEYLOTs7G4/H5e7DrVEVHMBXxre7b2XyMnaR+iJPKnQF00rGHYZ+oiTsMtOrGXoUjZ9dJN6AFtja2nJurFYr2ZYAaDPL5hO/GCeMRp3eun37drXNAdrozZs30k2osdlsdnx8PBwOh8Oh70d8ggHVaPmHGJXRDE5PT6fTqfvPwWAQ/Oz2acXby9gxgt6TIKa1LT1juxeNMB6P379///jx406nE/xpKz7BAEgjjGYQWjmAuSHJ2IZl0oyjgJFOT0+dJCrdEACtRhgFgJZyFnXa2dlx/vn06VPJ1gBoK8aMAgBSGY/HFxcXw+Gw1+tJt8VcJycne3t70q0w1Gq1cpag6ff7g8FAujkmcrro4uJib2+vPb9oLO0EAEh2enp648aNg4ODxq5nV9hyudze3n758qV0Q8w1mUz6/f7Dhw9PT0+l22KoJ0+eDIfDx48fn5ycSLelOoTRKui8SpMB5vN5gxeCadLRNelYYjTs98tY8/l8a2trfX1duiHm6nQ6Z2dn0q0w2tbWllMQvbi4kG6LoQ4PDzudznQ6vX79unRbqkMYrcJkMvnhhx+kW6HHzs7O6enp7u7ubDbz3v/555/v7Ozs7OyUvvBqmaKOro5ijqUZL5arSb9f4nzJ3vv3zHK5JImq2C6CI6aLOp3ObDbb2dn55ptvhFpnhMR30cXFRavyOhOYSre9vb1YLP785z9LN0SD6XTa6XQePny4XC6Pjo76/b5zv3MtwbrPfog6ujqKOZZmvFiuJv1+mcB7feOdnZ1Op7NYLA4ODvr9frfbnc/n7RnEFiWmi6SbZoqYLppOp7PZrDGfP7nFdNF4PD48PBwOh96lJBuPMFq6Kq7SVBX3q6jT6Zyfn7v3OyWTo6OjtbW10WhU0/JJ1NHVUcyxNOPFcjXp90ucN9kH/5758ssvj4+Pr1+/3u12pVsqJr6LpFtnhPguevbs2Wq1chZwaG0kje+iXq+3v79//fr10KV/m4owimzcX4/NzU33ztVqdf369Tt37iwWi/39/fp+xIQeXU1FHUtjXixo5032wb9nOp3ON998s1wu21wcje8iR8t/p+K7iEVtVVIXDQaDXq+3Wq1a9YtGGNUpxyWaTBY8HKWUO8zF++E7GAycnzpnYaptpk6hR1dTUcfSmBcLZQv+PbO+vl73UrpeTfrztSR0UaJgF3U6nVaVRRVhVK+GXaIpeDjT6XQ+n6vLcYfu/ScnJ91ut+5nqXq9XujR1VHMsTTjxUIFmvS3WUnookR0USK6SBFGkclgMDg5ORmPx+fn53t7e7PZbHd3982bN5ubm/v7+1tbW4vFor5x3Hd00s0pJHgsDXuxULYm/W1WErooEV2UiC5ycAUmZDabzUJPIkTdXy/NOApHzLE06TChkTOUzZnku729vbm56fw9w8VyXHRRIrooEV3kQxgFAITjj5ZEdFEiuigRXUQYBQAAgBiuwAQAAAAxhFEAAACIIYwCAABADGEUAAAAYgijAAAAEEMYBQCgdLPZbGdnx7ntrHOeg/NE76aql7vxQBTCKAAApet2uwcHB87t4+PjfBtxnujdVPVyNx6IQhgFAKB0y+VyOp0qpcbj8WKxcK7Bc3p6evfu3c8//9z553w+Pz093dnZcf45Ho8///xz96fuE72bunv37t27d92nHx0dHR0d3b59e39/f7Va+drgbPDu3bunp6ehT3duOA2bz+fBDXobD+hCGAUAoHQXFxeLxUIpdXh42O12Dw8Pneh5dnb24sWLxWJxenp6cXFxfHw8Go0ODw9ns9n5+fmrV69evXp1fn4+n8/dJzqbch5wdnZ2dnZ2fn4+m80uLi6ePXs2HA7fvHmjlHJ255rNZovF4sWLF2dnZ8fHx8+fPw8+3X3Ku3fvLi4ught021B5/6HJCKMAAAiYTqfdbnexWCwWi263++7dO6VUt9vt9/tKqX6///33389ms9PT09VqdXFx4Xv6y5cvB4PB+vr6+vr6YDB4+fKl8/Rer6eUun79eszjv//++//+7/8OPj0oZoOALoRRAABkLJfLly9fOkHwxo0b3h/N5/Pt7W3nR+vr66FPX1tbc26kSYrv3793r37e6/X+4R/+IdPTgfIQRgEAEHDjxo319fXDw8PDw8ONjQ03Gjqm0+nW1tbh4eFwOIx6ujuxfTabbWxsxO+u1+u55c/t7e0//vGPMU8/Pz/PcURAPn+UbgAAAO3ijBAdDofT6XRnZ6fT6SwWi++//947yvPOnTv7+/vqcujndDrt9/vOE50C59bWljPbaX19fblcHhwc+AaJ+ty/f393d/fo6Gi5XHa73X/913/d3d31Pt3Z13g8fv/+vS8ZhzZeV28Alm3b0m0AAKBFlsvlarVyxmLO5/OLiwtnnGjwYU5wVEotFot+v+99osOpbnrviTebzdbW1tzH+56+Wq2cAaxRAwN8jQe0IIwCAABADGNGAQAAIIYwCgAAADGEUQAAAIghjAIAAEAMYRQAAABiCKMAAAAQQxgFAACAGMIoAAAAxBBGAQAAIIYwCgAAADGEUQAAAIghjAIAAEAMYRQAAABiCKMAAAAQQxgFAACAGMIoAAAAxBBGAQAAIIYwCgAAADGEUQAAAIghjAIAAEAMYRQAAABiCKMAAAAQQxgFAACAGMIoAAAAxBBGAQAAIIYwCgAAADGEUQAAAIghjAIAAEAMYRQAAABiCKMAAAAQQxgFAACAGMIoAAAAxBBGAQAAIIYwCgAAADGEUQAAAIghjAIAAEAMYRQAAABiCKMAAAAQQxgFAACAGMIoAAAAxBBGAQAAIIYwCgAAADGEUQAAAIghjAIAAEAMYRQAAABiCKMAAAAQQxgFAACAGMIoAAAAxBBGAQAAIIYwCgAAADGEUQAAAIghjAIAAEAMYRQAAABiCKMAAAAQQxgFAACAGMIoAAAAxBBGAQAAIIYwCgAAADGEUQAAAIghjAIAAEAMYRQAAABiCKMAAAAQQxgFAACAGMIoAAAAxBBGAQAAIIYwCgAAADGEUQAAAIghjAIAAEAMYRQAAABiCKMAAAAQQxgFAACAGMIoAAAAxBBGAQAAIIYwCgAAADGEUQAAAIghjAIAAEAMYRQAAABiCKMAAAAQQxgFAACAGMIoAAAAxBBGAQAAIIYwCgAAADGEUQAAAIghjAIAAEAMYRQAAABiCKMAAAAQQxgFAACAGMIoAAAAxBBGAQAAIIYwCgAAADGEUQAAAIghjAIAAEAMYRQAAABiCKMAAAAQQxgFAACAGMIoAAAAxBBGAQAAIIYwCgAAADGEUQAAAIghjAIAAEAMYRQAAABiCKMAAAAQQxgFAACAGMIoAAAAxBBGAQAAIIYwCgAAADGEUQAAAIghjAIAAEAMYRQAAABiCKMAAAAQQxgFAACAGMIoAAAAxBBGAQAAIIYwCgAAADGEUQAAAIghjAIAAEAMYRQAAABiCKMAAAAQQxgFAACAGMIoAAAAxBBGAQAAIIYwCgAAADGEUQAAAIghjAIAAEAMYRQAAABiCKMAAAAQQxgFAACAGMIoAAAAxBBGAQAAIIYwCgAAADGEUQAAAIghjAIAAEAMYRQAAABiCKMAAAAQ8/8/b8SkvB0WzAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "clf;\n",
    "subplot(1,2,1);\n",
    "hold on;\n",
    "Pi = 2* acos(0);\n",
    "theta = [0:0.01:2*Pi];\n",
    "### unit circle, area = \\pi*r^2 = \\pi\n",
    "plot(cos(theta),sin(theta),'c-');\n",
    "\n",
    "x_in = [];\n",
    "y_in = [];\n",
    "x_out = [];\n",
    "y_out = [];\n",
    "val = [];\n",
    "\n",
    "L = 1000;\n",
    "for i = [1:L]\n",
    "  x = 2*rand()-1.0;\n",
    "  y = 2*rand()-1.0;\n",
    "  if (x^2+y^2 <= 1.0)\n",
    "    x_in = [x_in x];\n",
    "    y_in = [y_in y];\n",
    "  else\n",
    "    x_out = [x_out x];\n",
    "    y_out = [y_out y];\n",
    "  endif\n",
    "  n_in = size(x_in)(2);\n",
    "  n_out = size(x_out)(2);\n",
    "  val = [val 4*n_in/(n_in+n_out) ];\n",
    "endfor\n",
    "\n",
    "plot(x_in,y_in,'c.');\n",
    "plot(x_out,y_out,'r.');\n",
    "hold off;\n",
    "\n",
    "### circumscribed square, area = 2x2=4\n",
    "xlim([-1 1]);\n",
    "ylim([-1 1]);\n",
    "axis equal;\n",
    "axis('square');\n",
    "title(sprintf(\"after %d iterations, pi = %f\",L,val(end)));\n",
    "\n",
    "subplot(1,2,2);\n",
    "semilogx([1:L],val,'g-');\n",
    "xlabel('iteration count');\n",
    "ylabel('approximate pi value');\n",
    "axis('square');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e3bbc94",
   "metadata": {},
   "source": [
    "## Metropolis algorithm\n",
    "\n",
    "Given a configuration $\\{s\\}$, the evolution of the system is always toward <i>lower</i> total energy.  Therefore, for each spin one needs to:\n",
    "<ul>\n",
    "<li>flip the spin\n",
    "<li>calculate the new total energy of the system\n",
    "<li>if the energy is lowered by the flip, accept the new configuration, otherwise reject it and keep this spin as is\n",
    "</ul>\n",
    "A flip of a spin only affects the spins in the vicinity of this spin, i.e. the interactions are always local (for example over the four nearest-neighbour spins). Periodic boundary conditions take care of the edge spins, for example the East neighbour of the last spin in a row is the first spin in this row, etc.\n",
    "\n",
    "Clearly, such an algorithm should move the system toward an \"all up\"/\"all down\" (depending on the external field) states of the system in the ferromagnetic configuration, and into a \"checkerboard\" patter of up/down spins in the case of the antiferromagnetic interaction. \n",
    "\n",
    "But how does one mimic a non-zero absolute temperature and the related fluctuations for such a system?\n",
    "\n",
    "For any finite temperature, the fluctuations can be modeled by modifying the last step in the above algorithm:\n",
    "<ul>\n",
    "<li>if the energy is lowered by the flip, accept the new configuration, otherwise accept it only with probability $e^{-\\beta{\\cal H}}$\n",
    "</ul>\n",
    "In practice, this is achieved by generating a random number $r\\in[0,1]$, and accepting the new configuration if $r < e^{-\\beta{\\cal H}}$  (even if the energy is increased by the flip).\n",
    "\n",
    "Repeating this evolution step many times mimics the evolution of the spins system toward a thermal equilibrium with an imaginary \"reservoir\" at temperature $T$.  During this evolution,one can monitor the total energy ${\\cal H}$ and the total magnetization $M$ of the system:\n",
    "$$\n",
    "  M = \\mu \\sum_j s_j\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "def27aa9",
   "metadata": {},
   "source": [
    "## Plan of action\n",
    "\n",
    "<ul>\n",
    "<li>Modify `Evolve()` function to implement the Metropolis algorithm.\n",
    "<li>Calculate total energy and magnetization and keep a record of these.\n",
    "<li>Repeat until the values settle to a stable equilibrium value; may need to ignore some early states that may reflect the initial conditions, wait until those are \"forgotten\".\n",
    "<li>Change the temperature and repeat, recording the equilibrium values and the extend of fluctuations for each temperature.\n",
    "<li>Plot magnetization as a function of temperature, find the conditions for ferromagnetic and antiferromagnetic transitions by exploring a range of parameters.\n",
    "<li>As an independent part of the project, modify to consider one of:\n",
    "    <ul>\n",
    "    <li>nearest- and next-nearest-neighbour interactions\n",
    "    <li>triangular lattice\n",
    "    <li>anisotropic spin-spin interactions, e.g. East-West coupling is different in strength from North-South coupling\n",
    "    <li>lattice of spins that can re-orient at their sites, with the Hamiltonian given by\n",
    "        $$\n",
    "          {\\cal H} = - J \\sum_{i,j} \\vec{s}_i \\cdot \\vec{s}_j -\\mu \\sum_j \\vec{h} \\cdot \\vec{s}_j\n",
    "        $$\n",
    "        i.e. dependent on the [cos of the] angle between the two nearest-neigbour spins, and of each spin with the direction of external magnetic field.  Without loss of generality, one can assume that $\\forall j, \\vec{h}_j=\\vec{h}=h\\hat{k}$.\n",
    "    </ul>\n",
    "</ul>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "67d59842",
   "metadata": {},
   "outputs": [],
   "source": [
    "### a single step in the evolution of the spin system; streamline it as much as possible\n",
    "function Sij = Evolve (S,n,i,j)\n",
    "  J = 1; ### >0 ferromagnetic, <0 antiferromagnetic\n",
    "  ### wrap-around at boundary\n",
    "  #printf(\"For (%d,%d) NN are: W=(%d,%d), E=(%d,%d), S=(%d,%d), N=(%d,%d)\\n\",...\n",
    "  #  i,j,mod(i-2,n)+1,j, mod(i,n)+1,j, i,mod(j-2,m)+1, i,mod(j,m)+1);\n",
    "  Snn = [ S(mod(i-2,n)+1,j) S(mod(i,n)+1,j) S(i,mod(j-2,n)+1) S(i,mod(j,n)+1) ];\n",
    "  dE = J * S(i,j) * sum(Snn);\n",
    "  if (dE < 0 )\n",
    "    Sij = - S(i,j);\n",
    "  else\n",
    "    Sij = S(i,j);      \n",
    "  endif \n",
    "endfunction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "7c922561",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4QAAAJYCAIAAAC1p7+MAAAJMmlDQ1BkZWZhdWx0X3JnYi5pY2MAAEiJlZVnUJNZF8fv8zzphUASQodQQ5EqJYCUEFoo0quoQOidUEVsiLgCK4qINEWQRQEXXJUia0UUC4uCAhZ0gywCyrpxFVFBWXDfGZ33HT+8/5l7z2/+c+bec8/5cAEgiINlwct7YlK6wNvJjhkYFMwE3yiMn5bC8fR0A9/VuxEArcR7ut/P+a4IEZFp/OW4uLxy+SmCdACg7GXWzEpPWeGjy0wPj//CZ1dYsFzgMt9Y4eh/eexLzr8s+pLj681dfhUKABwp+hsO/4b/c++KVDiC9NioyGymT3JUelaYIJKZttIJHpfL9BQkR8UmRH5T8P+V/B2lR2anr0RucsomQWx0TDrzfw41MjA0BF9n8cbrS48hRv9/z2dFX73kegDYcwAg+7564ZUAdO4CQPrRV09tua+UfAA67vAzBJn/eqiVDQ0IgALoQAYoAlWgCXSBETADlsAWOAAX4AF8QRDYAPggBiQCAcgCuWAHKABFYB84CKpALWgATaAVnAad4Dy4Aq6D2+AuGAaPgRBMgpdABN6BBQiCsBAZokEykBKkDulARhAbsoYcIDfIGwqCQqFoKAnKgHKhnVARVApVQXVQE/QLdA66At2EBqGH0Dg0A/0NfYQRmATTYQVYA9aH2TAHdoV94fVwNJwK58D58F64Aq6HT8Id8BX4NjwMC+GX8BwCECLCQJQRXYSNcBEPJBiJQgTIVqQQKUfqkVakG+lD7iFCZBb5gMKgaCgmShdliXJG+aH4qFTUVlQxqgp1AtWB6kXdQ42jRKjPaDJaHq2DtkDz0IHoaHQWugBdjm5Et6OvoYfRk+h3GAyGgWFhzDDOmCBMHGYzphhzGNOGuYwZxExg5rBYrAxWB2uF9cCGYdOxBdhK7EnsJewQdhL7HkfEKeGMcI64YFwSLg9XjmvGXcQN4aZwC3hxvDreAu+Bj8BvwpfgG/Dd+Dv4SfwCQYLAIlgRfAlxhB2ECkIr4RphjPCGSCSqEM2JXsRY4nZiBfEU8QZxnPiBRCVpk7ikEFIGaS/pOOky6SHpDZlM1iDbkoPJ6eS95CbyVfJT8nsxmpieGE8sQmybWLVYh9iQ2CsKnqJO4VA2UHIo5ZQzlDuUWXG8uIY4VzxMfKt4tfg58VHxOQmahKGEh0SiRLFEs8RNiWkqlqpBdaBGUPOpx6hXqRM0hKZK49L4tJ20Bto12iQdQ2fRefQ4ehH9Z/oAXSRJlTSW9JfMlqyWvCApZCAMDQaPkcAoYZxmjDA+SilIcaQipfZItUoNSc1Ly0nbSkdKF0q3SQ9Lf5RhyjjIxMvsl+mUeSKLktWW9ZLNkj0ie012Vo4uZynHlyuUOy33SB6W15b3lt8sf0y+X35OQVHBSSFFoVLhqsKsIkPRVjFOsUzxouKMEk3JWilWqUzpktILpiSTw0xgVjB7mSJleWVn5QzlOuUB5QUVloqfSp5Km8oTVYIqWzVKtUy1R1WkpqTmrpar1qL2SB2vzlaPUT+k3qc+r8HSCNDYrdGpMc2SZvFYOawW1pgmWdNGM1WzXvO+FkaLrRWvdVjrrjasbaIdo12tfUcH1jHVidU5rDO4Cr3KfFXSqvpVo7okXY5upm6L7rgeQ89NL0+vU++Vvpp+sP5+/T79zwYmBgkGDQaPDamGLoZ5ht2GfxtpG/GNqo3uryavdly9bXXX6tfGOsaRxkeMH5jQTNxNdpv0mHwyNTMVmLaazpipmYWa1ZiNsulsT3Yx+4Y52tzOfJv5efMPFqYW6RanLf6y1LWMt2y2nF7DWhO5pmHNhJWKVZhVnZXQmmkdan3UWmijbBNmU2/zzFbVNsK20XaKo8WJ45zkvLIzsBPYtdvNcy24W7iX7RF7J/tC+wEHqoOfQ5XDU0cVx2jHFkeRk4nTZqfLzmhnV+f9zqM8BR6f18QTuZi5bHHpdSW5+rhWuT5z03YTuHW7w+4u7gfcx9aqr01a2+kBPHgeBzyeeLI8Uz1/9cJ4eXpVez33NvTO9e7zofls9Gn2eedr51vi+9hP0y/Dr8ef4h/i3+Q/H2AfUBogDNQP3BJ4O0g2KDaoKxgb7B/cGDy3zmHdwXWTISYhBSEj61nrs9ff3CC7IWHDhY2UjWEbz4SiQwNCm0MXwzzC6sPmwnnhNeEiPpd/iP8ywjaiLGIm0iqyNHIqyiqqNGo62ir6QPRMjE1MecxsLDe2KvZ1nHNcbdx8vEf88filhICEtkRcYmjiuSRqUnxSb7JicnbyYIpOSkGKMNUi9WCqSOAqaEyD0tandaXTlz/F/gzNjF0Z45nWmdWZ77P8s85kS2QnZfdv0t60Z9NUjmPOT5tRm/mbe3KVc3fkjm/hbKnbCm0N39qzTXVb/rbJ7U7bT+wg7Ijf8VueQV5p3tudATu78xXyt+dP7HLa1VIgViAoGN1tubv2B9QPsT8M7Fm9p3LP58KIwltFBkXlRYvF/OJbPxr+WPHj0t6ovQMlpiVH9mH2Je0b2W+z/0SpRGlO6cQB9wMdZcyywrK3BzcevFluXF57iHAo45Cwwq2iq1Ktcl/lYlVM1XC1XXVbjXzNnpr5wxGHh47YHmmtVagtqv14NPbogzqnuo56jfryY5hjmceeN/g39P3E/qmpUbaxqPHT8aTjwhPeJ3qbzJqamuWbS1rgloyWmZMhJ+/+bP9zV6tua10bo63oFDiVcerFL6G/jJx2Pd1zhn2m9az62Zp2WnthB9SxqUPUGdMp7ArqGjzncq6n27K7/Ve9X4+fVz5ffUHyQslFwsX8i0uXci7NXU65PHsl+spEz8aex1cDr97v9eoduOZ67cZ1x+tX+zh9l25Y3Th/0+LmuVvsW523TW939Jv0t/9m8lv7gOlAxx2zO113ze92D64ZvDhkM3Tlnv296/d5928Prx0eHPEbeTAaMip8EPFg+mHCw9ePMh8tPN4+hh4rfCL+pPyp/NP637V+bxOaCi+M24/3P/N59niCP/Hyj7Q/Fifzn5Ofl08pTTVNG02fn3Gcufti3YvJlykvF2YL/pT4s+aV5quzf9n+1S8KFE2+Frxe+rv4jcyb42+N3/bMec49fZf4bmG+8L3M+xMf2B/6PgZ8nFrIWsQuVnzS+tT92fXz2FLi0tI/QiyQvpTNDAsAAAAJcEhZcwAACxMAAAsTAQCanBgAAAAddEVYdFNvZnR3YXJlAEdQTCBHaG9zdHNjcmlwdCA5LjI1wZk/DQAAIABJREFUeJzt3b9vG2eeP/DRF9ueZJcHmAKSA+KCqpLdwkppHyBWCwhYK6UExFrcNW6sOMV5cYibTahGZeQD4u4iBVBzBQ3YweGAMwVc1leRhYu9Qky9DvUP6FvMLY/Rj6Eyw+FnSL5eMAKJw5nnGfJ5wrcezjzPwtnZWQIAABH+X3QFAACYX8IoAABhhFEAAMIIowAAhBFGAQAII4wCABBGGAUAIIwwCgBAGGEUAIAwwigAAGGEUQAAwgijAACEEUYBAAgjjAIAEEYYBQAgjDAKAEAYYRQAgDDCKAAAYYRRAADCCKMAAIQRRgEACCOMAgAQRhgFACCMMAoAQBhhFACAMMIoAABhhFEAAMIIowAAhBFGAQAII4wCABBGGAUAIIwwCgBAGGEUAIAwwigAAGGEUQAAwgijAACEEUYBAAgjjAIAEEYYBQAgjDAKAEAYYRQAgDDCKAAAYYRRAADCCKMAAIQRRgEACCOMAgAQRhgFACCMMAoAQBhhFACAMMIoAABhhFEAAMIIowAAhBFGAQAII4wCABBGGAUAIIwwCgBAGGEUAIAwwigAAGGEUQAAwgijAACEEUYBAAgjjAIAEEYYBQAgjDAKAEAYYRQAgDDCKAAAYYRRAADCCKMAAIQRRgEACCOMAgAQRhgFACCMMAoAQBhhFACAMMIoAABhhFEAAMIIowAAhBFGAQAII4wCABBGGAUAIIwwCgBAGGEUAIAwwigAAGGEUQAAwgijAACEEUYBAAgjjAIAEEYYBQAgjDAKAEAYYRQAgDDCKAAAYYRRAADCCKMAAIQRRgEACCOMAgAQRhgFACCMMAoAQBhhFACAMMIoAABhhFEAAMIIowAAhBFGAQAII4wCABBGGAUAIIwwCgBAGGEUAIAwwigAAGGEUQAAwgijAACEEUYBAAgjjAIAEEYYBQAgjDAKAEAYYRQAgDDCKAAAYYRRAADCCKMAAIQRRgEACCOMAgAQRhgFACCMMAoAQBhhFACAMMIoAABhhFEAAMIIowAAhBFGAQAII4wCABBGGAUAIIwwCgBAGGEUAIAwwigAAGGEUQAAwgijAACEEUYBAAgjjAIAEEYYBQAgjDAKAEAYYRQAgDDCKAAAYYRRAADCCKMAAIQRRgEACCOMAgAQRhgFACCMMAoAQBhhFACAMMIoAABhhFEAAMIIowAAhBFGAQAII4wCABBGGAUAIIwwCgBAGGEUAIAwwigAAGGEUQAAwgijAACEEUYBAAgjjAIAEEYYBQAgjDAKAEAYYRQAgDDCKAAAYYRRAADCCKMAAIQRRgEACCOMAgAQRhgFACCMMAoAQBhhFACAMMIoAABhhFEAAMIIowAAhBFGAQAII4wCABBGGAUAIIwwCgBAGGEUAIAwwigAAGGEUQAAwgijAACEEUYBAAgjjPIzr169WlhYePXq1citjx8/fvXq1fWfnyE9VL6tMNLvfve7hYWFx48f525L7969u3fv3ps3b9JfHz9+fPPmzZs3bz5+/LjII6mL3SrfaQ72vWa/A6iKMxjyl7/85eXLl3/5y19Gbk2S5LPPPst+/suXL5MkefnyZXah6aHybYWRkiR58ODBn//85xxt6c9//vNnn3323nvvDVpy2qr/+Mc//vGPf0wfzPfIoIiL3Sr3aab7ZvdKgKoxMsrPvHnz5u///u/fvHmTDq5sb28vLCwsLCzs7+8Pb02Hdr777rt//dd/TR9Jh47SJ9+7d+/du3eXHv/i0waH+vd///eLRxhsTQtNt54bWIKBqxrYq1ev/vEf/zH5a1tKkuRcc0obfDp4efGw9+7dG/ycxsp0nDUZCpq/9JHBAS92q5E1zOhHb968GRww+etw7MLCwvvvv58+cmnXvmb/BSiDMEqWd+/evXz58sMPP/z9738//Pj9+/eTJLl3797f/u3fpo/s7+9///33h4eHh4eH33///XfffXfpAc89bX9/f3Co//iP/7h4hMHWf/u3f/vqq6/SrV999dWXX35Z6okzpS62w0ET+od/+If0h/fff//LL7+8tDl9991329vbwwdMn5weJPXmzZt0oDRJkg8//DANfzkeuVj5QVVH1jDjNN9///3BAV+9evXVV1/97ne/e/ny5c2bN4cj9bmufbFj5nv9AXIQRsny+9///t69e8OfYamPPvooSZKbN2/+zd/8TfrI48eP//SnP/3pT3/6+uuvkyRJvxK9KB12Sp98eHj4+PHjwaH++Z//+eIRBlv/8z//87333kufcOPGjeGBJRi42A4HTei3v/1t+sPNmzdfvnx5aXP68ssvr/N3zmD0tOAP5wyqOrKGGac5fPB0r/39/Xv37n3++ec//fTT4ELSc137Yscc+SIAjIswynh8+eWXv/71r5MkOTeGes729vbLly/v3buXjokOf+Zd8wjpQYZHqmDg+q1oYLg5XRUTh3300Uf//d//nf78/fffv//++/keufY5XVLDHKc5soirOiZA2YRRxiMdyPn888//53/+J+Npjx8//vrrr7/88st379699957w0/OPsJHH3307t27zz//PH2Ca9q41DXbYZHm9Hd/93dJknz33XfpKOOvf/3rfI9kl5Jdw2ueZlpo+p374eHhjRs30gHUizI6JkDpou+goloG978P3wj/2WefpU1l+MEbN24kSZJ+p/ny5cvDw8O0Rd29ezf56129F++mT79zTJ+Zfvk4ONQ//dM/XTzCYOvR0dGHH36YPuHDDz90pzCXurQdDn5I21J6p/m55pQ988O5renB01KKPHLx4Nes4aWnOdh3+JkPHjwY/N/+66+/PlfcoGtf2jEBJmPh7OzsF8ZXSJIkeffu3Zs3bz766KPBN5vp+M1VQy/D0sGhwb6DQ7179+7iEYYLevXq1c2bN69TBHMrox2ea7RFmlN6B9LwvvkeyV3Di6d5sUsOCk2fOfIihHMdE2AyhFEAAMK4ZhQAgDDCKAAAYYRRAADCCKMAAIQRRgEACCOMAsAY9Pv9Xq83/Ein0+n3+1H1gWkhjALAGBweHn777bfpz/1+f319/eDgIP1vbMWg4oRRAChqfX19d3d38OuLFy/u3Lnz9OnTo6OjVqsVWDGoPmEUAIo6Ojr69NNPB7/W6/Xj4+N2u314eFiv1wMrBtX3q+gKAMCsqdVqi4uLrVar2+1ubGxcfMLt27cnXysq6+3bt9FViCSMAsCY7e/vNxqNNIbevXv30jyakT8WFs4v1n3ukcGv53646r8VLOjcgxUpKP01SZLyChouMS3IXya+pgcAIIyRUQAYs08++WRzc/Pk5KTb7a6trUVXByotawgaAMit3W7XarVarXZx0+3btydzmWD2d80KqkJBE2sMlWVkFABKsbq6Gl0FmAKuGQWAmTWx7z8VRG7CKAAAYYTRMAsLpWxdWMjcN/O42YVmbM4utMiZZhdaxmGLFDpqM1cqqTskSZK7xY94LyO6Q3bPzl1o/u6QfdTMrQADwmikfNFlGgPPVXUueC6F8mg5hZKb7jBir4hCC5rCNwcIIIxGync5ysi9Qq5yyS70qq0Fq5qv0IJcRFSSkrpDiDJa5sidQrpDtkq+OUDlCKPzZ4oGoELGNksdrWWyRr8ZV7xfub+jnz7FusMMvRBAGGF0/lRzKGmSSn0FvLzT5Yr3a8S7OEvv8iydCzCdhFGup4Lf/c9SoZTDe1nUqO7gFQaKE0bDFLnAK2PriM+OzM259y3pXMorNPfOhQ7L1cprQrlbfO7B0ZBziekO2UfN3AowIIwCABBGGAUAIIwwGqbINNSTn9Z+xNYChebfOm1nSob8E9cXWTchc7/p6g4lvUSjStUdgDEQRqkuH2ZwfboLMKWE0UghM35WasLQUitTqTNlpHyDo9P4gufrDrPU8ZkxC/4WohhhNFLIkjOVmi6p1DVjKnWmjJTvRvJpfMFL6g4VnOlhGt8dYPKE0VlUxnDEyGNO0cBLwQEooz1TpKQFhGapkegO49Pv93u93vAjvV7v3CPARb+KrgBT4uysip86uceWM85l5Jka7ZkhOd/LanaHfHSH8Tk8PPzpp592dnbSX588eZIkSa/XazQaGxsboVWDSjMyGqas2enPzsqasz17qv3JzwFezgT+I7YWKZSr5Z+4vkh3KFSnanWHQt/R6w7jsL6+vru7O/i13W4nSfL06dO9vb30Z+AqRkYBoKijo6Nmszn49fXr18vLywcHB0mS7O3txdULpoCRUQAYv2+//TZJkpOTk83NzUufsDBkojWjArz7w4yMAsD4ra2tpZeKrq+vX/qEs3m7lIEhw+/+7du3A2tSBUZGw5S1EEtJS60UKbXA1twrMJW0tdDLy9VKfDMnX2p2k846aIGFlEK6Q9Zu8z7x5PLy8uDnfr8fWBOoPiOjkRYW8lzjf525Vmbj7+1qRruZeXmrxgubrZrdgatsbGysr683m81ut7u9vR1dHai0BV8TREk/Wq56+TM+mEfumLF1xAd+7q3lHDbmXMqrElcrqTskSbKQcad5RHdII+WVLShjUwW7Q+Zd/Nlb50S73a7VarVa7eKm27dvv337dvJVKkN2q2akWWoM+RgZjVTGLJm5D1tUCTksZjLHUUVKmyUpozuEDSZmdod8LWh0d/C3UPWsrq5GVwGmgGtG508ZqxaVvft4D+v7TgbmpwmVc6ZT+EIAlSOMzp/cA1BFnlDSgE0Z52Jsaa6EdIeSlNGwR+2ltwDF+Zo+TFkLsZS01EqRUgtszb0CU0lbC728XC1kXaKQ7jCiBWVvrVp3yNpNTgWuy8goAABhhFEAAMIIowAAhBFG4+ReeSh7c0mHjVmAKWbJmZLOlAyFFg4bceicCwjpDv+7sZwyAQaE0blT0lQ2ZXwshXyYFTxTcsv3uoa8G/PTHUaqZq2A6SKMzp2SprKpYKElqWatZkDOqeDHXIvrFao7/FU1awVMF2E0VMgoZTnCln3KsTVzLx+uMyVozvaZ6Q7GPYEJEEaZXVclAnmTOaQ7AFUljIaapdVfQszPmZLPXC0gVEJ3mKnXB6gqYTRO7pWHsjeXdNiYBZhilpwp6UzJUGjhsBGHzrmAkO7wvxvLKRNgwHKgAMC1pJcRn/3814s/T/gvkfIubvYn1WQIowBATmlcC7/XTWqcar6mB4Ax6Pf7vV7v3IO9Xq/f74fUB6aFMBqmvMWQSlqmJXehJZ1pSeeSvdmSMyUpssZS/n2r1x1GnGm+Co3amn8JK93h5w4PD7/99tvhR/r9/vr6+tLSUlSVYCoIo5GmaCGWgoVW6kxLPZc5/AAel0o1klILnaVGojsMrK+v7+7unntwd3d3cXHRyChkE0YjTdHMTgULrdSSM6Uun+MO4twq1UgKFjo/jWR+znSko6OjTz/9dPiRZ8+eLS8v12q1eRgZnae3mvETRiup1PGEXAuxlKWCo5RzNZgz60a/l5XqDpnyn0vm1tzf0ZOt0+l0Op0HDx5EVwSmgDBaSaWOJ1RqIZaQM4WBSnWHTEUrdMUZVe48Z8WzZ89u3brVbDZ7vd6TJ08u3tiUJMnCkMnXsLgF7aeAaX/3x8vUTkyzkq5I8L+GWXFWgRlnJkd3qJIHDx6cnp4mSXJ8fNxoNBYXFy8+56x6f/MwMcPv/u3btwNrUgXCaJiQxZCquC5R7iWjyjmX7M2ukCtJkTWW8u9bve4w4kzzHzfnGkvZ++oOGVZWVtIfFhcXV1dXYysDFbfgLzMAuFS/3+92u8OPjCtZ3r59++3bt2M51CQNr8B07mv62BWYpjrKTGljGCMjowBwuT/84Q9Jkty6dWvwiGFOGDthNMzCQtbXWPm3Fjhu7l2zdyxSaMbmIi9gWac6olSuVFZ3yN4c0R2yh3BmqTtM+2BV6sWLF3M+ZAUT4G56rmXkPQzTODEUlTJL76buMDPu37/f6XSiaxFsYejf4JGRu0yMhj8DjIxyLSPvqZ29iaGYsFl6N3WHWbK1tVWv1we/Pn/+PK4uYQYTUwyuFr30OQMTDoja/rQTRgHgchsbG41GI7oWMOOE0VlUwesa8xUasjj9da5IMAY1LWbp3dQdIgxmaALKI4zOnyn6BrGkCbezD1vNKxIoyRS9m7pDkHa7/fr16yRJbty4cf/+/XlYaB4mzA1MYQrNjJ09FXze+bhzF1rWtPaZm8ubab+0+f250vx0hyINc7q6w2x0hmazubu7e+PGjZWVlZOTk62tregawQwyMgoAlzs8PHz16lU6GtpoNB4+fNjpdHx3D+NlZBQALnfpmvLAeBkZBYDLra2tbW1tNRqNxcXFTqfT6/UMi8LYWZs+TFmrFpWz3FH21pDlc6auUDKUt2pRxr6FCs192Oq1zBFLJWUUWmQ1qekxuIFpeXl5bW1tXDcwTdFy5AuXzTN6lt0LJnjR8Ays9TVFjaEkRkYjVWq+o4KmaMKoaSx0HpTRHUavE5Ov0HIOW5CWOV6bm5uPHj1qtVrdbjd9pNvttlqt+Zz0HkoljM4dn1gQSx+cCtvb27VardFofPzxx9F1gRnnBqZIuSffHCFzsKikT8FSDjtqEHh2CqWc7jDykDkLLXbYWeqDlfyeZjxWV1eXlpZardbqkF6vN1dL1S9c/RZfZ3n6CTSPGW6Bc8XIKNdTzYsDYH7og5PV6XR2d3e73e7ga/okSbrd7tHRUWCtJu/swg/Jdf4eS5JkUknR+MAMcAMTAFyu2Wzu7Oxc88n9fv/09LRWqw0e6XQ6S0tLw48MTMU9KwXvDZrArUUzcPdSMiWNoVRGRgHgcjs7O51O5/T0NP211Wo1Go3V1dVLn3x4ePjTTz+l4bXf729tbdXr9V6vV6/Xr59oYQ4JowBwuWazeXx83O/3BzM6PXr06NJnrq+vd7vdTz/9NP318PDwzp07aQa9e/euMAoZhFEAuNzh4eEPP/xwcHCQJMnGxsbm5uZVzzw6Omo2m4Nf19bW0h/6/X7JdSzLbHwDzlRwNz0AZKnVaicnJ0mS1Ov14fuZsnep1Wrtdntra2t7e7vkCsJ0MzIap8CSKSWtwDRdiyGVt9ZUWQtccbXcyx0l2W0kSZKSlkqq2ApMI1pm3jWWsvedhxWY7t+/v76+/s033zx58iRJksPDw08++eSa+zabzR9//HFvb+/SG5iSJFkYmiFhJm8mzl6lac4Nv/sffPBBYE2qwMjo9Ck4wUtJ88OMWAhneialmaKqzph8L3xId7jOCkxjL7SaZuhULrezs/Po0aOlpaW9vb0kSTKS5TkHBwfZSTRJkrMhY6sxU8K7P8zI6PQ5Oyv0YRYy4fYU9bWCLy+55Wsj2e/X2ai0lHvS+yJtZIq6w0gzdCpXSu+dX1lZWVlZuf5enU6n1+sNrjG1iChkMM9onPQjdOzfQWceNv9XhCUdtkCh+auUfdjSXgcylPQ1ffa+U/Y1fZHukCTJVRXOveN8fE2/ubk5PDHT5ubmuGJl9aeWHMs37GV/TT8blwFUvzGUzdf0zJ8Z+ISEX+jKRq87jHJ6evrkyZPpvSkeqk8YDTVFa1SHfGKNLLSatYIyRHQHbT1JkqdPny4vL29tbfV6vei6wGxyzWicAldZZm3MvePIz7KrNxc5bO5Ci5xpyMtLhhFvV/bW7AaWe8eQw5bUHTKPmnvfIu/LdHnw4EG9Xs+YZJQJcDH/DDMyCgCXazQa6Q+rq6t7e3u/6B4mxu7sr/8d/scMMDIKAOdtbm4+evTo5OTEgCiUTRgFgPO2t7drtVqj0fj444+j6wIzztf0YUZMZpl31uyFhaxdixSafdisXfMeNntzoXnFC7y8BUrlSiPekOx9Mw+buzvkfqsLdYfMQrNaZvZ+eQ87otACL2D1ra6uLi0ttVqt1SG9Xq/T6URXDWaNkdFIIfPw5Ss0ZN2m6+yef2XFcgolt5wtM6TQac9Z45N7Ntbq63Q6u7u73W53eDH6brd7dHQUWCuYScJopClaDGnkukQlHXbk7rm3llQouZWxGFJUu52fRpJ7zoHqW1lZef78ebPZHMx4D5REGJ1C1RylHHnw6gxAlbpwuIFTRrqikRT6Mj33Klaldofpt7Ozc3BwcHJysry8vLi4mH53H10pmDWuGZ1C0zjiV6l8Nj9nSjVd0UiKzAY6usxyDjvzms1mq9VKkuTk5OTHH3/c2tqKrhHMIGGUCrPAElxbKQ137rvD4eHh3t5eekN9OvV9u92OrhTMGmE0TJElgkpaDCn38kJlLXeUubm8M81dqE/u3IpceljOokVlHbakpZKKvERl/d8mu9Apsbi4OPxrr9c79whQnGtGAeBy29vbW1tb9Xq91+s9fPgwSRKLMMHYCaMAcLmNjY16vd5qter1+vLy8sbGRnSNJmG8t62lR/ulI+VTPS8Yv5QwGifztuv8M2gWmHszpND8W6tXqDvpcxtxJ3iRN/PqI5dU6IgPUd1hqvT7/dPT08EiTO12e3V1NbZKkzGuty49zhjT7fS3KS4hjHItcz/BC6UrYyAkqtka1JkZf/jDH5IkuXXr1uCRjDCaJtdarTZ4pNPp1Go1s0FBNmGUayk4Oz2MVEZ6y54SvzyS6Mx48eLF27dvr/nkw8PDn376aTBJ/ubmZq1W63a7jx49mpPxVMjH3fQAcLn79+9fczH69fX13d3dwa+tVqtWqz19+nRvb29/f7+0CsIsMDIaKmT19EoVWnD1lwoWCiOVtCCZ7lCO9G76wa/Pnz+/9GlHR0fNZnPwa6fTSe+7r9Vqx8fHJdcRppswGmqWJnXPN1Vh9tf/BZcADykURsrXhHSHCBsbG41GI9++g4tH79y5M74awQwSRuMUmAo+a2uBueBDCs2/tXqFzvrnconKWaNgxJFLKrSkGeZHlao7lKLIrKK9Xi/94aqR0YWhoH82G68X1zb87n/wwQeBNakCYRQALtfpdE5PT4cfWVxcvE5CXVlZSS827fV6w9/yDxNA59nwu3/79u3AmlSBMAoAlzs+Pt7f30/TZDrSWavVTk9Pj46OsndsNBrPnj1rNpvHx8cPHjyYRF1hai34ywwALrW5ufno0aPBUOjDhw+/+OKL/f39lZWV61xL2m63a7Xa8MyjA7dv377+pFETNvaJcnMc8OIuMzx9b5Ubw2QYGQ1TwXVYcq+nMmOFZmwtdFiuVuIKTLmbUN4qFTlsWd2hnEWhCh12SnS73eEoubi42O12b9y48eOPP15nd9OLwnWYZzRSvmnkrzPXStUKzSek0JFM/l8pxScmGnuhuQ9bwUK5f//+w4cPDw4OWq3WkydP0mzaarXW1taiqwazw9f0YdJPjhzjHCN3zHfY7M3zU2iRKs3GUFCIspaJL9KEMha1L+ew/7tzGd3ByGgB7Xb79evXSZLcuHHj/v376f1Ml37z/otU55vZ4fdx8BdNqV/TD//ddGlBw9VYGPphFtrTZarTGKL4mj5SyGyD+cxPoSPHl2bi47WKyvikKakJ6Q7/t31288HA4uLixx9/nP68u7vbaDRm+8v3iryhZ0nMWr6EEEYBJmoeAtzMSG+H7/f7S0tL6SOPHj2KrRLMHmEUfm7kABTlmJ98Nk1nOqo7TNO55HJ4ePjDDz8cHBwkSbKxsbG5uRldI5hBbmAKU8F1WHKvpzJjhWZsLXRYrlbO2zVia+6lm8o7bFndoUChWd2hyGGnSq1WOzk5SZKkXq93u93o6sCsEUYB4HL3799fX1+v1+svXrxoNpuHh4fFb10CzhFGAeByOzs7jx49Wlpa2tvbS5Jkb29PGIWxc80oAFyu3+/3er1ms5n++vr169m+lR5CGBmNU2S27oytJR02c+uIG34yN5d0okUKLelMyRLRbsvqDlm7jWqZeXcdUdn8Ncradx66w8OHD9vtdnQtYMYZGWU8ypjf+jpLzuQotOBH5KzM5M2UKdgd8s0nNbKvzHx36Ha7P/zwQ3QtSrQw6/PJMxWMjDIeJc3yXUahBas62x+9VFbB7pCv2Y7ca+a7w/379589exZdC5hxwmioSi1OX8aS9uPafZKFTuOZzoYy3q+gxenzfVM/g01nJrrD7u7u7SHR1YEZ5Gv6UCGjClcVWupk75U609krdJ4VbLdVaiQzuP7h9HeHf/mXf/mv//qvwfJL0+Vic8poY2W3veuXe9Uq9swwYXQKVXNh7JkxP2daNVPUbud+XaL5cefOnegqFHJpsLvYPifZYs+VtfDXBxcuZOWznz+HGSaMxqngEkwR6xKVdKJFCrUCU4CIdltWd8jabVTLzLtrSatJZe87D92hVqvdu3evXq8PHnn+/Pn1d+90OrVabUoHVmFihFEAuNzGxkaj0cixY7/f39raqtfr7XZ7e3t7Y2Nj7HWDmSGMAsDlVlZW8u344sWLO3fu7Ozs9Pv9hw8fCqOQQRgFgDGr1+sHBwftdrvb7Q5/yz8xl04daj5RqsnUTmEKzUiTdyWWQiumVG/ZpzJqVOSw7nzKrcgCQiFLJYV0h+zFkCZfoxGFZh119tVqtcXFxVar1Wq1lpeXo6sDlWZkNFLI4iUzv2LKwFVnagWmavLCZiuY7fJ1h6KdZY7H4fb39xuNRvrt/N27dy/9mn5h6PU9m+nWP4PTlhU2/O5/8MEHgTWpAiOjkao0xeFUyne7c5FbqK/zBPLxwmYruBhSSHfwlmY7GxJdFybNuz/MyOgUClmyvejC2BHDMiGM701YsSaUc+gupDswVT755JPNzc2Tk5Nut7u2thZdHag0YXQKhczjXdI6N6UWGkISnbBijSTnuxXSHZgqtVrt+++/b7fbn3zySa1Wi64OVJowGqbQt2N5J78uNEn1VM20HzJLuiSaW6Fve3N3h9yHLVBoka2555+PWTQga+O8WF1dja7C/3HhJpXlmlEAmCnzfOsY08jIKABMqzkZ7JyT05xbwigATLHrD4IGDpeeXT1eezb036v2Zbb5mj5MkenTszYWmPS+epN8z1ShZCjUHfIeuoKNZH4KBRiuuU7XAAAO00lEQVQQRqdPyP/irzOVjUIpYope2PlpmbpDxbk2lNkgjEaqVKwstTIKpZoq1Ujmp1CAYcJopHwTARWcPihmIZa5L5SRKrUg2fy0TN0BCCeMziLDEcyN3BeMwrTzHT0zQxgFppgPY4BpZ2qnMEXW8snaWGAFpuqtODNThZKhUHfIvWP1Gsn8FAowYGQUAIAwwigATBkXjDJLhFEAAMIIo3EiFkUpa5mWAoXmLrWkw47YWORMyRCxglBZXTBrv1GHzd4372ELdYcyduSXW/j5vxk222fHpYTRuZMvMFUzZoUsOcMsKak7FEmG+Y6Z+7BMkbOf/5s6GQvTX+dBZpgwOnfc4goDlVp4gtnT6/V6vV50LaDqTO0UamEhz+fSdZaLHvfH3dnZNb4HHPtn7KgzLeNDPeZMSWKaUEnKaB9nowZHy2iVo18+3eFqT548SZKk1+s1Go2NjY0ih3K7ErPNyCgwzXw5TSW12+0kSZ4+fbq3t5f+DFzFyGio3N8RZn8AV2qR74LHrGDUMA5UkjK6QzWb0PQYORyrO1zl9evXy8vLBwcHSZLs7e0VOZRhUWaekdE4EYuilLVMS4FCc5da0mFHbCxypmSIWEGorC6Ytd+ow2bvm/ewhbpDGTvOh2+//TZJkpOTk83Nzei6QKUZGQWA8VtbW0svFV1fX7/0CQtDw/ZnVyR7w6Kzavjd/+CDDwJrUgVGRgFgzJaXlwc/9/v9S59zNmRS9aIqvPvDjIwCwJhtbGysr683m81ut7u9vZ3vIIZFmRPCaJzsKVHybk0H/q/ctUChGRtHzO6S97AjtpbzAo7YWuSwZCjt3cz6OC+tOyRX98HseJG/YWanFt0hwtHRUbvd/uSTT2q1WnRdoNKEUa6lpLlNi6zHCFEKtsx8w13XWfZprqNfJa2urube1xvK/HDN6AwqaZKlMgotcjszRCnYHfK169GF5josQDhhNNRVYx2lDhjmK7Sgkg6e74Uq+PIazi1JCe/mtRYQylGZaaQ7VM9C5j9/XTA/hNFQV411lDpgmK/Qgio1D//snek8KOndDGkkIXSH6jnL/Def5vnc55kwyjSboisSGKmEF9ZbVZTuAJRPGI0TsahMSesSFSm0pGWfQha48sGcX2nvZjmNr0AbydpY1rJPugNQZcIoAABhhFEAAMIIowAAhBFG42RPSZM9I8rVmxcWMvctMg9L3kILzf1SzpkWKTT/Ybla3jdk9Ob87aB63SG7vefvDlllZp5pzv0AfkYYZRIqNZljSfO0juXg82x+kk3OGT/LqErJZuldA8ojjM6gSs3pmb21pDkiy9jrmru7gTi3IreZV00Z3WHkC1DBl6iCVQIqSBgNlW/cIGQFptkbwyypUCaspBWYytgr1hV1zv0d/bV2B7gGYTRUpUb2KjjiFzJwynSZxtH1EFfUudAEpZYVAMbhV9EVmGMFpoIPmfQ+d6GFUm45Z5p7Lv0KJvbZUOg7+oxGMqLUKesOuZeAyL8uQOa+s3RlBRDIyCgAAGGEUQCYWQuTusRZQeQmjAJAWXq9Xr/fj64FVJowCgCl6Pf76+vrS0tL0RWBShNGAaAUu7u7i4uLOUZGL35BfO6Rwa/nfrj438oWdO7BihSU/lBqQcMluhIgJYwCwPg9e/ZseXm5VqsZGYVsC2em3wCAsep0Os+ePdvb29vc3Hz+/PnFJ9y+fXvilaK63r59G12FSMIoAIzZw4cPb926lSTJixcvVldXt7e3a7VadKWgokx6DwBj9uDBg9PT0yRJjo+PG43G4uJidI2guoRRABizlZWV9IfFxcXV1dXYykDF+ZoeACI1m83T09ONjY1BhC3Vs2fPHjx4UGoR/X5/d3c3SZLV1dVGo1F2Qaenpw8ePJjAq3dwcFCv10stqNPppC9dvV7f2dkpr6BKcTc9AIQ5ODhYXl5+9OhRGkFK1ev11tfXX79+XXZBh4eHq6urT58+PTg4KLWg/f39jY2Nvb29Z8+elVpQkiS9Xm9/fz+9+qI8p6enjUbj+fPn85NEE2EUAAJ1Op21tbXJTP9Uq9WOjo4mUNDa2lo6IFp2dNvZ2anVaq1WK71drFT7+/tra2tll9Ltdtvt9sOHD9vtdtllVYcwCgDl6vf7vV5v+JFOp5NOht/r9caYRDMKGq/sgmq1Wrvd3tzc/OKLL8Ze1sUzOj09HUvqzSio2WxubGwUL2JkQfV6/cGDB1988cUERsqrQxgFgHIdHh5+++23g183NzcPDg62trba7Xa9Xu90OhMoaFxFXKegVqvVarWeP38+lssrh8s6V1Cz2VxaWtrY2DiXjMde0PLycqvV6na7rVar1IJ6vd7KysrS0tJczcDgbnoAKNH6+nq32/3000/TX1utVq1We/r0aa/Xe/LkydOnT3d3d2/dulWv10staIw39Y8s6MWLF/1+f3NzM0mSS+f8z1fWxYJWVlbSKV2LT+OaXVA6LNpsNj/++ONSC1pcXEzPqHh7mCLCKACU6OjoqNlsDn7tdDrpeGGtVjs+Pq7Val988UU6HlZqQYPHC6bD6xS0t7dXsIhLy7pYUKPRWFlZ6ff74331rnrpxnJT0cgzWl1dHUt7mCLCKABM1GAY786dO0mSLC0tlXQD07mCyhNYUK1WK2N1q8AzKq89VJZrRgFgogYXOA6PuilIQaUWVGXCKABMzsrKysnJSZIkvV6v1OsCFaSgaeFregCYnEaj8ezZs2azeXx8XOpKSApS0LSwHCgATFq73S7pYkcFKWjqCKMAAIRxzSgAAGGEUQAAwgijAACEEUYBAAgjjAIAP9Nut9OV5ZMk6XQ6+Q6S7jh8qMnLXXkmSRgFAH6mXq8/evQo/Xl3dzffQdIdhw81ebkrzyQJowDAz/R6vVarlSRJs9nsdrvNZjNJkoODg7t37/7mN79Jf+10OgcHB5ubm+mvzWbzN7/5zWDrYMfhQ929e/fu3buD3Z88efLkyZPbt28/fPiw3+8PVyA92t27dw8ODgaPnNs9/SGtWKfTuXjA4cpTZcIoAPAzp6en3W43SZKdnZ16vb6zs5NGz6Ojo1evXnW73YODg9PT093d3e3t7Z2dnXa7fXx8/MMPP/zwww/Hx8edTmewY3qo9AlHR0dHR0fHx8ftdvv09PTFixcbGxtv375NkiQtLtVut7vd7qtXr46OjnZ3d/v9/qW7D3Y5OTk5PT29eMBBHSJeQn4BYRQAGKHVatXr9W632+126/V6uqJ6vV5fXV1NkmR1dfWbb75pt9sHBwf9fv/09PTc7q9fv240GktLS0tLS41G4/Xr1+nuKysrSZLcunXrqid/8803V+1+0VUHpOKEUQBgtF6v9/r16zQILi8vD2/qdDrr6+vppqWlpUt3X1xcTH8YmRR//PHHwdqYKysr6QGvvztTRxgFAEZYXl5eWlra2dnZ2dm5cePGIBqmWq3W2trazs7OxsbGVbsPbmxvt9s3btzIKGtlZWUw9rm+vt7pdLJ3Pz4+zndSVMSvoisAAFRXeoXoxsZGq9Xa3Nys1Wrdbvebb74Zvsrz448/fvjwYfLXSz9brdbq6mq6YzrGuba2lt7ttLS01Ov1Hj16NLz7Offv39/a2nry5Emv10u/ea/Vaud2T8tqNps//vjjuWR8aeXH+IIwdgtnZ2fRdQAAKqrX6/X7/fRazE6nc3p6ml4nevFpaXZMkqTb7a6urg7vmEpHN4cfydButxcXFzN27/f76QWsV10YcK7yVJYwCgBAGNeMAgAQRhgFACCMMAoAQBhhFACAMMIoAABhhFEAAMIIowAAhBFGAQAII4wCABBGGAUAIIwwCgBAGGEUAIAwwigAAGGEUQAAwgijAACEEUYBAAgjjAIAEEYYBQAgjDAKAEAYYRQAgDDCKAAAYYRRAADCCKMAAIQRRgEACCOMAgAQRhgFACCMMAoAQBhhFACAMMIoAABhhFEAAMIIowAAhBFGAQAII4wCABBGGAUAIIwwCgBAGGEUAIAwwigAAGGEUQAAwgijAACEEUYBAAgjjAIAEEYYBQAgjDAKAEAYYRQAgDDCKAAAYYRRAADCCKMAAIQRRgEACCOMAgAQRhgFACCMMAoAQBhhFACAMMIoAABhhFEAAMIIowAAhBFGAQAII4wCABBGGAUAIIwwCgBAGGEUAIAwwigAAGGEUQAAwgijAACEEUYBAAgjjAIAEEYYBQAgjDAKAEAYYRQAgDDCKAAAYYRRAADCCKMAAIQRRgEACCOMAgAQRhgFACCMMAoAQBhhFACAMMIoAABhhFEAAMIIowAAhBFGAQAII4wCABBGGAUAIIwwCgBAGGEUAIAwwigAAGGEUQAAwgijAACEEUYBAAgjjAIAEEYYBQAgjDAKAEAYYRQAgDDCKAAAYYRRAADCCKMAAIQRRgEACCOMAgAQRhgFACCMMAoAQBhhFACAMMIoAABhhFEAAMIIowAAhBFGAQAII4wCABBGGAUAIIwwCgBAGGEUAIAwwigAAGGEUQAAwgijAACEEUYBAAgjjAIAEEYYBQAgjDAKAEAYYRQAgDDCKAAAYYRRAADCCKMAAIQRRgEACCOMAgAQRhgFACCMMAoAQBhhFACAMMIoAABhhFEAAMIIowAAhBFGAQAII4wCABBGGAUAIIwwCgBAGGEUAIAwwigAAGGEUQAAwgijAACEEUYBAAgjjAIAEEYYBQAgjDAKAEAYYRQAgDDCKAAAYYRRAADCCKMAAIQRRgEACCOMAgAQRhgFACCMMAoAQBhhFACAMMIoAABhhFEAAMIIowAAhBFGAQAII4wCABBGGAUAIIwwCgBAGGEUAIAwwigAAGGEUQAAwgijAACEEUYBAAgjjAIAEEYYBQAgjDAKAEAYYRQAgDDCKAAAYYRRAADCCKMAAIQRRgEACCOMAgAQRhgFACCMMAoAQBhhFACAMMIoAABhhFEAAMIIowAAhBFGAQAII4wCABBGGAUAIIwwCgBAGGEUAIAwwigAAGGEUQAAwgijAACEEUYBAAgjjAIAEEYYBQAgjDAKAEAYYRQAgDDCKAAAYYRRAADCCKMAAIQRRgEACPP/AbXeHSdXw2HaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 8;\n",
    "n = 2*N+1;\n",
    "rn = rand(n);\n",
    "mu = 0.5;\n",
    "S = mu*sign(rn-0.5);\n",
    "\n",
    "clf;\n",
    "subplot(1,3,1);\n",
    "DisplayIsing(S);\n",
    "title(\"initial state\");\n",
    "\n",
    "S_1 = S;\n",
    "\n",
    "M = [ sum(S(:)) ];\n",
    "L = 10000;\n",
    "\n",
    "for k = [1:L]\n",
    "  i = randi(n);\n",
    "  j = randi(n);\n",
    "  S(i,j) = Evolve(S,n,i,j);\n",
    "  M = [M sum(S(:))];\n",
    "endfor\n",
    "\n",
    "subplot(1,3,2);\n",
    "hold on;\n",
    "DisplayIsing(S,'blue');\n",
    "DisplayIsing(0.5*(S-S_1),'cyan','cyan');\n",
    "hold off;\n",
    "title(sprintf(\"after %d iterations\",L));\n",
    "\n",
    "subplot(1,3,3);\n",
    "semilogx([1:L+1],M,'c-');\n",
    "xlabel('iteration count');\n",
    "ylabel('magnetization');\n",
    "axis('square');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb856fd8",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Octave",
   "language": "octave",
   "name": "octave"
  },
  "language_info": {
   "file_extension": ".m",
   "help_links": [
    {
     "text": "GNU Octave",
     "url": "https://www.gnu.org/software/octave/support.html"
    },
    {
     "text": "Octave Kernel",
     "url": "https://github.com/Calysto/octave_kernel"
    },
    {
     "text": "MetaKernel Magics",
     "url": "https://metakernel.readthedocs.io/en/latest/source/README.html"
    }
   ],
   "mimetype": "text/x-octave",
   "name": "octave",
   "version": "4.2.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}