2D_SCFT.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code for Chapter 6: SCFT Numerics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import some libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pyfftw\n",
    "from scipy import integrate\n",
    "from numba import jit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define some helper functions.\n",
    "\n",
    "`get_fselector` defines which part of the copolymer sequence belongs to monomer type A (denoted by 1's) and which to type B (denoted by 0's). `sf` selects according to a head-to-tail perspective, while `sft` selects according to a tail-to-head perspective. These correspond to what the propagator functions q and qt 'see' as they 'walk' in their respective directions along the contour of the copolymer sequence. \n",
    "\n",
    "`get_k` gets the corresponding wavenumber/frequency grid.\n",
    "\n",
    "`solve_phi` solves for the respective monomer densities according to eqn 5.47a and eqn 5.47b."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_fselector(fa,fb,Ns):\n",
    "    ka = int(fa*Ns)\n",
    "    kb = int(fb*Ns)\n",
    "    sf = np.concatenate((np.zeros(ka+1),np.ones(kb+1)),axis=0).astype(int)\n",
    "    sft = np.concatenate((np.ones(kb+1),np.zeros(ka+1)),axis=0).astype(int)\n",
    "    return sf,sft\n",
    "\n",
    "def get_k(Nx):\n",
    "    freqs = np.fft.fftfreq(Nx,(1.0/float(Nx)))\n",
    "    kx, ky = np.meshgrid(freqs,freqs)\n",
    "    k = kx + ky\n",
    "    k2 = kx**2+ky**2\n",
    "    return k2\n",
    "\n",
    "def solve_phi(q,qt,Nx,Na,ds,V):\n",
    "    ax1 = np.linspace(0.,V,num=Nx)\n",
    "    phi = np.zeros((Nx,Nx,2))\n",
    "    for j in range(Nx):\n",
    "        for k in range(Nx):\n",
    "            # here is where the flipping of the backward propagator occurs.\n",
    "            tmp1 = q[j,k,:]*qt[j,k,::-1]\n",
    "            \n",
    "            # eqn 5.47a and eqn 5.47b\n",
    "            phi[j,k,0] = integrate.trapz(tmp1[:Na+1],dx=ds)\n",
    "            phi[j,k,1] = integrate.trapz(tmp1[Na:],dx=ds)\n",
    "    Q = 1./(V*V) * integrate.trapz(integrate.trapz(q[:,:,-1],x=ax1),x=ax1)\n",
    "\n",
    "    return phi/Q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The main SCFT function. The algorithm of section 6.1.1 is referred to as simply 'the algorithm' here.\n",
    "\n",
    "The modified diffusion equation is solved according to section 6.1.2: solving the MDE. Note that here, both the forward and backward propagator functions are updated and stored in the same direction, so a flipping of one of them will be necessary to account for the opposite directions (i.e. the head-tail perspectives)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def scft(Nx,fa,V,gam,XN,it,wisdom=None):\n",
    "    # define number of monomers for each type.\n",
    "    N = 20\n",
    "    fb = 1.-fa\n",
    "    Na = int(fa*N)\n",
    "    Nb = int(fb*N)\n",
    "    print('fa = %.3f' %fa)\n",
    "    print('fb = %.3f' %fb)\n",
    "    print('')\n",
    "\n",
    "    # define copolymer sequence contour step-sizes (ds), and contour resolution (Ns).\n",
    "    ds = 1./(N)\n",
    "    Ns = int(N)\n",
    "\n",
    "    # define time step-size (gamma_1 and gamma_2 of eqn 6.1a)\n",
    "    eps = gam\n",
    "\n",
    "    # define the Flory-Huggins parameter such that X_ab == constant always.\n",
    "    Xab = XN/N\n",
    "    \n",
    "    # random seed for reproducibility.\n",
    "#     np.random.seed(555)\n",
    "\n",
    "    # initialize empty arrays and initial conditions.\n",
    "    w = np.zeros((Nx,Nx,2))\n",
    "    C = 0.1 # scaling factor\n",
    "    w[:,:,:] = C*np.random.uniform(-1.,1.,(Nx,Nx,2))\n",
    "    phi = np.zeros((Nx,Nx,2))\n",
    "    q = np.zeros((Nx,Nx,Ns+1))\n",
    "    qt = np.zeros((Nx,Nx,Ns+1))\n",
    "    \n",
    "    # import pyfftw wisdom if it is defined.\n",
    "    if wisdom != None:\n",
    "        pyfftw.import_wisdom(wisdom)\n",
    "    \n",
    "    # initialise pyfftw arrays.\n",
    "    qpp = pyfftw.empty_aligned((Nx,Nx), dtype='complex128', n=16)\n",
    "    qppt = pyfftw.empty_aligned((Nx,Nx), dtype='complex128', n=16)\n",
    "    qp = pyfftw.empty_aligned((Nx,Nx), dtype='complex128', n=16)\n",
    "    qpt = pyfftw.empty_aligned((Nx,Nx), dtype='complex128', n=16)\n",
    "    v = pyfftw.empty_aligned((Nx,Nx), dtype='complex128', n=16)\n",
    "    vt = pyfftw.empty_aligned((Nx,Nx), dtype='complex128', n=16)\n",
    "    \n",
    "    # define pyfftw FFT plans.\n",
    "    fft_object = pyfftw.FFTW(qpp, qp, threads=1, axes=(0,1))\n",
    "    ifft_object = pyfftw.FFTW(qp, v, direction='FFTW_BACKWARD', threads=1, axes=(0,1))\n",
    "    \n",
    "    fft_object = pyfftw.FFTW(qppt, qpt, threads=1, axes=(0,1))\n",
    "    ifft_object = pyfftw.FFTW(qpt, vt, direction='FFTW_BACKWARD', threads=1, axes=(0,1))\n",
    "\n",
    "    # get the monomer configuration in the polymer.\n",
    "    f, ft = get_fselector(fa,fb,Ns)\n",
    "\n",
    "    # get wavenumber and spatial grid.\n",
    "    freqs = np.fft.fftfreq(Nx,(1.0/float(Nx)))\n",
    "    kx, ky = np.meshgrid(freqs,freqs)\n",
    "    k = kx + ky\n",
    "    k2 = kx**2+ky**2\n",
    "    \n",
    "    # get spatial grid.\n",
    "    dx = np.linspace(0,V,num=Nx)\n",
    "\n",
    "    # calculate xi and phi at first iteration. Step 1 of algorithm.\n",
    "    pressure = (w[:,:,0] + w[:,:,1])/2\n",
    "    pressure = pressure-pressure.mean()\n",
    "    xi = Xab*N-pressure\n",
    "\n",
    "    # calculate the initial phi values.\n",
    "    phi[:,:,0] = (w[:,:,0] - xi)/(Xab*N)\n",
    "    phi[:,:,1] = (w[:,:,1] - xi)/(Xab*N)\n",
    "\n",
    "    # for each subsequent iteration...\n",
    "    for i in range(it):\n",
    "        # store the old phi values for equation 6.1c, step 4 of the algorithm.\n",
    "        phi_old = np.copy(phi)\n",
    "\n",
    "        # define initial conditions (eqn 5.30)\n",
    "        q[:,:,0] = 1\n",
    "        qt[:,:,0] = 1\n",
    "\n",
    "        # solve the MDE.\n",
    "        for s in range(1,Ns+1):\n",
    "            # solving MDE: substep 1 (eqn 6.5a)\n",
    "            qpp[:,:] = np.exp(-0.5*ds*w[:,:,f[s]])*q[:,:,s-1]\n",
    "            qppt[:,:] = np.exp(-0.5*ds*w[:,:,ft[s]])*qt[:,:,s-1]\n",
    "            \n",
    "            # solving MDE: substep 2 (eqn 6.5b)\n",
    "            qp[:,:] = np.exp(-4*np.pi**2*ds*k2/(V**2))*fft_object(qpp)\n",
    "            qpt[:,:] = np.exp(-4*np.pi**2*ds*k2/(V**2))*fft_object(qppt)\n",
    "            \n",
    "            # solving MDE: substep 3 (eqn 6.5c)\n",
    "            q[:,:,s] = np.exp(-0.5*ds*w[:,:,f[s]])*ifft_object(qp).real\n",
    "            qt[:,:,s] = np.exp(-0.5*ds*w[:,:,ft[s]])*ifft_object(qpt).real\n",
    "        \n",
    "        # get the phi values using helper function.\n",
    "        phi = solve_phi(q,qt,Nx,int(fa*Ns),ds,V)\n",
    "\n",
    "        # get the pressure term (2nd term of xi, equation 5.47e).\n",
    "        pressure = (w[:,:,0] + w[:,:,1])/2\n",
    "        pressure = pressure-pressure.mean()\n",
    "\n",
    "        # equation 6.1c.\n",
    "        delta_phi = phi_old.sum(axis=2)-phi.sum(axis=2)\n",
    "        delta_phi = delta_phi-delta_phi.mean()\n",
    "\n",
    "        # update w(i) to w(i+1), equation 6.1a. This completes substep 4 of the algorithm.\n",
    "        w[:,:,0] = w[:,:,0]+gam*(Xab*N*(phi[:,:,1]-0.5)+pressure-w[:,:,0])-eps*(delta_phi)\n",
    "        w[:,:,1] = w[:,:,1]+gam*(Xab*N*(phi[:,:,0]-0.5)+pressure-w[:,:,1])-eps*(delta_phi)\n",
    "\n",
    "        # print every 100th iteration.\n",
    "        if i%100==0:\n",
    "            print('it = %i' %i)\n",
    "        \n",
    "    wisdom = pyfftw.export_wisdom()\n",
    "\n",
    "    return phi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define parameters used and run the SCFT simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fa = 0.700\n",
      "fb = 0.300\n",
      "\n",
      "it = 0\n",
      "it = 100\n",
      "it = 200\n",
      "it = 300\n",
      "it = 400\n",
      "it = 500\n",
      "it = 600\n",
      "it = 700\n",
      "it = 800\n",
      "it = 900\n"
     ]
    }
   ],
   "source": [
    "# Nx = spatial resolution; fa = fraction of copolymer sequence belonging to monomer type A.\n",
    "# gam = 'temporal' step-size (see eqn 6.1a)\n",
    "# it = number of iterations; V = length of domain\n",
    "# XN = Flory-Huggins parameter multiplied by the polymerisation index\n",
    "Nx = 128\n",
    "fa = 1.0 - 0.3\n",
    "gam = 0.2\n",
    "it = 1000\n",
    "V = 8.0\n",
    "XN = 30.0\n",
    "\n",
    "# run simulation.\n",
    "phi = scft(Nx,fa,V,gam,XN,it)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot resulting morphology."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATYAAAD8CAYAAAD9uIjPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOy9e5Al2Vkf+PtOZt5n3arq6uru6emZkUYgOzCwEqAQ7BKBMRJGZh2SIxaDIPAKB8QsGwi87HrXaNkVDnlxyLsOGzaCYJkF2WAwAmRjz+7OWgiwgtjAYA2gQEgCNBqNpnu6p1/V9byvfHz7x/edk+fmzXsrb9Wt7urq/EXcyLz5PJm/POd7nO98h5gZNWrUqHGWYB50AWrUqFFj2agbtho1apw51A1bjRo1zhzqhq1GjRpnDnXDVqNGjTOHumGrUaPGmUOlho2IfpiIPk1Ef0JEv0xErZMuWI0aNR4eENGHiOgWEf3JjP1ERP87Eb1IRH9MRF/t7XsPEX1Of+9ZRnkObdiI6AqAHwLwFmb+CgABgHcv4+Y1atQ4M/jnAN4xZ/9fA/BG/T0D4KcBgIg2APwYgK8F8FYAP0ZE545bmKqmaAigTUQhgA6A68e9cY0aNc4OmPl3AGzNOeRdAH6BBb8HYJ2ILgP4FgAfY+YtZr4H4GOY30BWQlihwK8S0T8G8AqAAYDfYObfKB5HRM9AWmJ0OviaL/2SEBmADAQASO2SpS2NOcSI5fbDNJJlEiIbBwAAM5brmhgwiYyOoExHSWTefztywg2g4Mn1ktXJgpf8IQBk1/Ml2/VAllmglw4J+gjghtwoaiRohzEAoGkSAECDZBkhRUCZPJ/e2YBh9Pp+kWiygIV/i6H4ChjsttkRKPM4G3OIsXI2sJzFEXgs+wPLWQKQ5Sy1/NxHzuw2szzOQqQITxFnV68muLuVHefS+Ja/0uW7W2mlY//gj0efBjD0Nj3LzM8ucLsrAK56/6/ptlnbj4VDGzZVC98F4GkA2wB+jYi+m5l/0T9OH/JZAHjTfxLx889vYshAn+VL2ssaAIC7WRcA8Fq8ji+MLgAA/nTvkixvX8Lgiz0AwMor8vl0XsvQuicvP9qTj870dTmKgVg+PEqUoCQFMvkA3ZIZs4aOEVFeIQL96gPj1jmSV8TNENyUmpCsyHK8ppV8I8DBFbnG4PVSu69c2cKXb9wAAHxp5xYA4MlIBNpj4Q7WzUCeTytQhxgtLUdE2lCAYLQaBbrPeEq2qVBlMq9qZJD3kbpGLEPMsm3olrM5ezXewMvDTQA5Z5+9+RjiL8p+x9mtDM172ohbzgby34xiYCzbHGdplnOV5hXtyJyFumxE4LZy1i1wdj7AweOTnD3+xBa+YgZnF8JdnDd9eb5TwNnbvvX2odc5DHe3UvzHjz5V6djg8ueGzPyWY9yu7MF5zvZj4dCGDcDbAXyBmW8DABH9awD/GYBfnHUCESEiQgpGC/Khxir5ekYa/WHQx6VoFwCw32nKcr2Jl0ZSpP64I9dK84+C9RVEJpfCpNoCYq0QcQJKCw1bxrM/J+NVEiPX4jAAXOWQ8mStEElHKke8KvsG52TZv0QYXJHn27goz/T61S083twBAFwI9wAA64FUjBbFiFT6R8phAKkUE0WDKa0c8v9Ywjq/jq2YTktjRFqZWspZl6Tinw/20W9IY2c5O9ho4vOxvoekDQCgzACq2dl3G/qcGcuZXB9p6nGmDVXGoGLD5ngq4SwwgBVCZZz1lLONac7OXRR+Xte7N5OzLo3vG2cGNNG4Te2nok54NDDyhvM+4BqAJ73/T0BcWtcAfGNh+8ePe7MqDdsrAL6OiDoQU/RtAF447KQAhAbl0qdJ0vC0INK6Zwa4oI3HsCnFGPdCxPphf1FbsQO0we7jlWtzqJUlMggG1nRVzS0OgMRK/0IDVwYikfaA17AZVznSplw/6YaIV7S85+S4wUUpR/+JBKuPSUV4el0k/FPtLVyKpJLYytExI13GaOn7iGxjTeQamWJlOQ78SuIqGllN1iBVwWPvGYDRIjk+1TLaco85cBV+pNrrqBdilMq7upoqZ9wC0yRnTctdZBCEso9GspHiNNfU0lzLxqyGzdfYrAAq4Szt5JwNCpwNnkjRuyzP8oZzdwEAT3bunQrOfJRxhiU1RgxGzNVM0SXgOQDvJaIPQzoKdpj5BhF9FMA/9DoM/iqA9x33ZlV8bL9PRB8B8IcAEgB/BDU5a9So8XBjWRobEf0yRPPaJKJrkJ7OCACY+f8A8DyAbwXwIoA+gL+t+7aI6B8A+IRe6gPMPK8TohKqaGxg5h/TglYCQdRlwwz1z6KlEieDaFaZGSKDdUqLVM1aBilPSr5XsIGDUMLmsoYcnzbkmGaTEB7ItnComtsohYknNTVKS6S/fTZjnGOZVZPIQoO0pVK/LcvxinGa2lBcgxhelmdZfWwPX7JxBwDwdFek/+XGTm7OmNyckXeROqnf8CS+lfq+j2YZJqg9p8y8CXznO8Qktf6chnJm3QmrZogY8j7iSN8PGyRqmltcM+dwEMzgbNcgVI0qGMp1zSiFsVq2LinL5mpsbDlTbTtrBMisdm056xmM1lVDuyDL4ePqEnlsD19a4OyJxj1shPsApjmLkJ0azpYFBjuuj30t5u88ZD8D+IEZ+z4E4ENLKYiiUsNWo0aNs4mTbDgfJE6sYQtAyEDInCYgL9BK/9TESCE9g2lopZNxoQWB+nlCynC1sQ4A2G9Lh0LSlmLHWwEau3J8QzW3YBgiGGu3fGylP+duCSuhyEr8XPpnkWoXLULSUm1yRfaNzhGG5+Xc5JJI8c0L1q92F6/rqG+tKctL0TbWgwMAnp/GhXvkmqwLHaD5Un1ZnQVyLc9vw7YEVqMW36hAfW3KWUbkNJlUNaU0Msg6qrUoZ4HJcC20nEmHgtWi4ns5Z1Ffte1BAGM588NEZmkTASGzmpqqUWnLTHE2XicMN+UasXK2cUE6d55e33Ka2pXmPQDS87luCr415axF2anibBndBwzpLDqLqDW2GjUeYdQa2wIgjeUxxIjYbQSQv8gWUmQkPaRQKYkQMKo5RBorFFKKbiTS82pLOk7urqwAAMZrDTTvqfTfFY0gOmCEQ9lm7OUTdp1KXiEBiMKSWqmvqkrSBuKeSv01Ke/4fILmedEwr6yLpvZUT7SzJ1rbE1IfEB+NDW3pelIfkF61qNCbdlj80zwElB+X8mxncJnfxt7Lam7zOINJnOabaq+hXMM+l2h2TZOgG4mGdLUlmtv2isS6je800LynISZ7ljODcKg9t7EuUy8gW2GDbdl4mprPmdWu11VL20gRKmdPnrOcCU9PtbdwuSE9oJeibQDCWR7eYQN0TzFnxwQDiM/o1AAnqrEFoDwIXF+gVedBLBUFyM1E03clsqZoRClWAmnY1hrSUFxvrwEAXlvtYW9XTJ3hjlw53DeuQ0FPg4kBrXMgG/Su31UWAhqHilSH9sddRroqrWK0Jhe5tNrHY12pHE90pCI8pjFPF6NdnA/U6azmZ88MvQZtMkygQbnT2QZ2zop/mmfO+JWjuG1eZcmvbaZ6xRbmLIA3iiJv4LqhvLfVSDi70V0FANxc6+FgR160cZzlnUCBN+LEcubqdAlniXKWdBnJmgrDVXUVrB3g8RURNI93hKvHm8LdZrjnxarN5qz5EHB2VDC4NkVr1KhxxsBAejbbtRPsPCCaDCakSSc1mHP1SbUAw4wAVkuwWk6Cjg4cXQvFrNiIRMJe6qzg1qoMwdreFM1tt9/CaCiaQKYhIIgNKLUmlxZDi8MBA025V6jLTmeE9bZoGptt0cQutPZxoSHrm5FI+g2npfWxqmandTp3KSmV+vJseWBn0ZSpijLJX7a/TAvwzZu5AaAVOEMmnQUAEITqZjCxew+9QN9jU97VpXYPt9fElbDVF84OBk0MBsIZD23QrgGUs1zLlhWOGIj0nm15x+3OCGuOM/k+hDPhajOS+9twjo1g/6HmbHkjD84mao2tRo1HFuQSHZw1nFDngUq1knCCCS3AhV7YsXeJ06SMahANStHSToauStY1dfBebOzhckt8KNuxSP+9uIX9RMYxDhId9pOEiLNC0KT14ZkMLc3o0AlFM1yJRliPRDs8F/X1ngN3X+dgtlqJGaCrZXRDx2g6oNOOxzSev2aiTCcwHjQgM9N3M3O41QKcwSS5L45yLdtyZjW2c6FoUZvRPnY0BGR7RcJ3duIW9uNJzsZpgDidDPy1nIUmQzvSb6KEs3XlbCM8cJz1AtlntbSeGbixsBMhHQ8RZ8eFdB7UDduRUOx1m6wsak7oHgMGYGO9tGEz2VTD1jNSCdaDPjbVAbzXEC9yP22ir57lUabxbhwgUwIz1/uXO7ptihq77AQjZ/7m9xy6dVuejvuf5h0E1pSm8sohz3v0DBCHmTPLQBXO/AYuUM6sGyEwA/eOitytB30cRJr0wHKWNdBPpzlLuNCwwTZsOWctY7kYO85sI9YxI3dfO4LgrHJ2FEgcW92w1ahR44whqzW2xeF3TZdqASrJMif9ATjTyJ7HiFQtjzQTgZW+QzPEkDXhoV1mkRt7GnMu/dOCyRDo9Q1xngCyxJRySxO79Yaea83OhmfC2NCIyDNdihkgjhImcBxUdUrb8lThzNe2jTvXvr8MLbYmuaY+Yqv5RjlXgWhpQ44QK2djjzOLopYdIJvgCpjkx+esgclynCXOjotaY6tRo8aZA4OmBP5ZwYk1bNbJaaVbuRZgfTeyiFhyCALwwj4YsacJAHA5pFqcYIyRblMtzYROClmNbZ5UCsCedpj73Rrunvm9nd/Pi0aXZ819M/YzichMZX5w9zxCBoiT9tOUhYDM5yzXtoucBWAvXEczhChnMccY6/Vjk2vWLsOL7ktBbtxwEQFl+ffhadv3m7PDcFp9az5qU/QIKOt1m6gs1sTQ4zPKnFljP9wYPNHIyTY7KDvFWK+R6QeeYuzIcjn7D2s0Cvc0lGeRDcirrF6lACYrRpkJMyveadEKcj9RlTN424zrpcvN03QBzmKOS7nKZjRsxmvYfM7yxi5z+06Ss5MyQ+8XGOQEyllDbYrWqPGIQgJ0T6+QPQ5OvGErc047FEemzzBP7eDjtCClU+TS2Q4NybxzfelfTGBpYTUy/7oGLvekK20Av5vfbpsd4+Q7m/1tRVSV+kUn8iJmTpUxiGVlmsvZhDY1zzy1yzmcUTqlVaecj1ktg69JSxkfbc6OirrzoEaNGmcKzLP9mA87qky/9xcB/Iq36Q0A3s/MP7HIjcrSr1T1u6VOKltfiJX+nI91s1oaM7wJ3Nxa5mlm/vV9+N6GoqSX/eWhAP6zlI0fPAn/zP2Q6PM4m9C25/hKMzcqIde67Zl+ZolilokiX5PlmkaZdibbD+fMHUuz98m208/ZosiWpLER0TsA/CSEip9l5g8W9v9TAH9F/3YAXGTmdd2XAviU7nuFmd953PJUmczlzwC8WQsQAHgVwK8f98Y1atR4sJDOg+Mbbdou/BSAb4ZMp/cJInqOmT/j7sX8w97xPwjgq7xLDJj5zccuiIdFn+ptAD7PzF886g3nagFAqd8t74ebDD8wvkbmaQbzcky51G9zyjgh1T0pXuaTcfvOcC/aYdr2lJbqc1bg04BLtbgI05il35RrbDVni2KJnQdvBfAiM78EADrF3rsAfGbG8d+JBSaHOgoWbdjeDeCXl3HjWYN5pz6uEpMnP7Yk1gqTDd5RcNhg56rmZr7vbFSOMs785y7jopSzksbuuKg5OxpmdaotiCsArnr/r0HmDp0CEb0OwNMAftvb3CKiFyADxT/IzP/muAWq3LARUQPAOzFjMlMiegbAMwDw1JW6T6JGjdOOBUcebGrjY/EsM9v5hctax1nS6t0APsI8MVPzU8x8nYjeAOC3iehTzPz5qgUrwyIt0F8D8IfMfLNspz7kswDwNW9qVhLBRal4qAY3NXEBEMDOb+mbSdOYlU75sGDZeckEDzt3GVL/JCfbOM58lxYzXQoWJZyVaXFFzEt/XXO2PMwKgi7BHWZ+y4x91wA86f1/AsD1Gce+G4X5RZn5ui5fIqKPQ/xvx2rYFjGwvxNLMkNr1Kjx4CGD4E2l3yH4BIA3EtHTatm9G8BzxYM0wuIcgP/gbTtHRE1d3wTw9Zjtm6uMShobEXUgPR7/VdULZ4UQjSrwj63kf8MMv04JAhxt6EjV4U/HlaQPYhq0Wfc8ac5KtbgCjsrXofcuPe5o3B2HszKNtEp5l/mdMGgii8qRr8OcENF7AXwUEu7xIWb+NBF9AMALzGwbue8E8GGdFd7iywD8DBFlEEXrg35v6lFRqWFj5j6A80e5wVErT9n+yhXnhHG/GrFlzUY0D+XC4uxxlt/7ZBuxRTlLmZ1Qrnrust4fM5YWoMvMzwN4vrDt/YX/f7/kvN8F8JVLKYSH2stfo8YjC1pagO5pw4k0bBIfY+eaLJcIk07nai+3ynHLUNWX7qCtUKYyaZ3Omcx2Ec1groQvMQsfRs6q3qsq7idnWckph3G2jDfGWJ7GdtpQa2w1ajzCqBNNLggnuWY4iicDOxfXBGbhtARWzpOph0n60v3zrleiJfjR9ynSiX0TwaxlErvmrGTfcjgr48qiKmcZyyzuxwWD6kSTi4A1620Amqo4vqN0nuP6tHzsi6KsclQxWTJk0wPBvWOm9h1WEM7T+UyXxx9+ZNP/zG7sas4E94szYJq3Sc6W854ZeZbps4az+VQ1atSogHrC5IXAEMmVgac1Al8bmBozOG3qPCxaQBWpX2a6WKleJukzTGsJ1kBZJKjAT7wISLp1QNN6s71Xnt7HmkE1Z8fjbNKYPBpnFkGBs2V1Hiww8uChQq2x1ajxCKPW2BYAM2PI2QyNwNPgitKCpn04MmvS0V7+UWcJWjQh4CzJX+aTATTZopPqsyW9bMvXZZ+m+vGcvmWxSH7GjLI02va/vy5LctvmcTbP73Y/OTtq8saHnbN5nRBVwUy1xlajRo2zBek8qGepqowMwJDt/JK6zU604cJAvIlbLNiU+nAW9d3cr/kcfalf5puZ2manouN8j5X4MXJJH3sSPoadIMUmZcyXVtqWmRP5tHQZAp6U/nYyFcPs5tu0bywCO1+c5Syyj0nypABK/W7H8bcdlbN5M6aX4UFwNnb/82dclDOrNUecISBe0oiBR3jOg6MgA6HPASKvOzyfGzQ/yndY223zTJ3T4pwumjF+ZUid6Z1NOZnzSsKI3fGyjEEY2kmf9R3EHGCMwK3LMp8E2n6UZVlQ/Umg7eTF/oTQANCiGA2223SJzM0Fapf5HKGEuS4FTyjVnJVzVgyvSNk4XiwiSkonhAaABlJEJbN6HQXSeVD72GrUqHHGUI88WAApCNtZEy1K0FJJY2dWTz0toFHB1AmInHR90FpAVakPyPPFqg3Zo8bWhOHcdPEl/lCl+ZAjtzzImrpf92WR/g8mNIEifInvNDQjOkeLdGnG+bpbTnMW1ZxV5uwga+b8Oa7yYzJPs7Mw3sz19n/OkSyt5tY1I7QoXkpK73rkQY0aNc4k6pngF0DCAW6nq+jQCF0zAiCaAAB0dClS0nYkyMIwu9nEfS3gsCSSRVhH8jJDB6o6nX1nc+w5mWWbLIccOKlvJf5B1kSfm24dAPbSFvayFgCgn8q2ftYAAAzSCInVHDJZZky5k9moT8YkaBp5571gCADomDEAYCUYOn56ZgAg1wgAoKvLpmpwJ8nZUXFaONtNlaesOcXZvi5HWYhRJucm2XRvZKicNU2CdiAlsFx1AsvTEL1giHgJVZcZiLO6YauMMYd4NT6H9aCPIUuFWTVSqVLb0CF1Ee6Zc07LPwCTlaXgnPbNmwfplC5zOltTxnc2D3nShBly4D5+24htpx3spF0AwE7aBgDcizvYinXbWI+PtfGLGxglaopq9yUzgbRhCwMpRytM0Imkcqw2hIO1SJbrUR/nwr5s0+V60HeNXGxsQyicdZEAJ8TZUYXRojgpzrYSWd5LOtiJhT/L2V6sDV0cYRhrw5b6pugkZ40wQSeSkvSUq3MN4WQ96uNc1Mc4W0LDhjqOrUaNGmcQ9ciDBTDKQnxucAmb0T42wz0AwDDcBwCsqwRcNUNvQJzGbRN7GQ58LaAQO+VpAfOc0lW1gEXNmXlOZ9+UsVK/r9L1QJ3Ke1kbuyr9t5IVAMCdpIdb4x4A4NZIlrcHK9g66Mi5fe1EGGjAzNDAjOW51LoHmACjZVNmuZEBLTVLO6IFrHREC9js9HGhLbxcau4CAM5HB9iMhLM4kH1j5SwzQ6S2A6Lm7FDO9g/k+GQg16JhABrJ9SnVpafccijl4AaDlbNIOet2RGve6Apng6xsiunFcJbDPSrpoUS0TkQfIaI/JaLPEtF/etIFq1GjxklDTNEqv0OvRPQOIvozInqRiH6kZP/3ENFtIvqk/r7P2/ceIvqc/t6zjCerqrH9JIB/x8zfptNrdeYd3E8b+NS9x3G5s4srbZFatgs8DuWWMQKkEL9Opk5qmEQ0AGBSCyhxTgOYGFt6mBYwSwOYJfnLwgTKpD5Q7nQeMjmpv6sOZuujuZus4E6yCgC4Gcvy1cE6Xj1YAwDc2hWNoL/dhtmRa0R7Uv6uvDIEI0D9yjBpXlZWpz0HGuDcMEha6tfpSMfDvVWh795aDzfWRdO41JPlE91t7DfV2d0QzjZC60kzSDHQm47tG3ScuRELzHnnAiYxK6fbIrgfnFlNbTuTd+VzdmMsPN0YrpVyFmxPctbUVxYMgWAs9yRf4dV3lYX6LUdA2pK3aTnb7Ynfbnt9BTfXexjEx9fYgPIxq4uCiAIAPwWZye4agE8Q0XMls039CjO/t3DuBoAfA/AWiBL5B3ruveOU6dAvjIhWAXwDgJ8DAGYeM/P2cW5ao0aNBw/pFQ0q/Q7BWwG8yMwvMfMYwIcBvKtiMb4FwMeYeUsbs48BeMeRH0pRRWN7A4DbAP4ZEb0JwB8A+DvMfOAfRETPAHgGAILz63jp6gXcO9/G/ppI/0FXJMyoKcs0ovzuRtQQGUfqiTKIFuAyGczQAqqg6Ls5alYIYDofVwZM9ab1s9BpatupSP2tVKT6zXgNV4cbAIBX+ucAAFe317F9V/aHd+Qdde4RGuL6QnSg72OkyzE7Tc1KfzDcO7Lxn1lISDWqNmnJMu7KOxivNnFwXu71+U3VCDba2FnVXryuaAvDpva+MiENCxLejBFAnHzGCwWxAbwZTWrZfhjIolr2w8pZONTnjBkmKWps7GnZsi0LCVk0i7MGhhsRstES5gNdLEB3k4he8P4/y8zP6voVAFe9fdcAfG3JNf4LIvoGAH8O4IeZ+eqMc69ULdQsVGnYQgBfDeAHmfn3iegnAfwIgP/ZP0gf8lkAaD3+JLdfamL7IMTgksZdbQjxyWpOiAsO9EuhJs5EZSlEu09UlgVDQKpUjlnOZ7ueORNGlmPmiXgnQEwZa3raynF9LBXii8MNvLS3CQC4enddrnGzg/ZNObe5Jddq7maI+tqQDbRhHuvzJpkbtEj+NEe2ktgGLjTIwtwsBYCkrdHvOwaNXVkf7UpZ7xyE6F/UyPlUOVvRzoNW3tAEoTcG2G12vRgwTggVOhYW4GyRhuxUcWaF0HCaM0oKrhZPGDnuQnJmaaqcpY4zQmPXQMPajo0FTNE7zPyWGfvKLsKF//8XgF9m5hERfT+AnwfwTRXPXRhVlJ9rAK4x8+/r/49AGroaNWo8xLC9olV+h+AagCe9/08AuD5xL+a7zGyb4/8TwNdUPfcoOFRjY+bXiOgqEf1FZv4zAG8DMHcK+mAMdK8xgkGAwVACF1+K1Rla0sPiMlGEnCc8NNZMSb0xdPb4XAsI1GV92Dymi6IYqe4HdLqlN47QSv1+locH3E0mpf4XBiLxX9zdxNXbsg2vigm48hqhdVeu19yVe4cHKcKh2CxkpX6s/5MMyDJbuGlo5wEb49azhpQx7OfLqG90XU3FQYjBUDoSPp/IcakXne7GNZLlJ5se6+hxVgzePclw0JPg7Itqfv75zsWT58y+HGPAqrEFkZSRHU8hwgHDJCXnHwFLCtD9BIA3EtHTAF4F8G4A3+UfQESXmfmG/n0ngM/q+kcB/EMi0peLvwrgfcctUNVe0R8E8EvaI/oSgL993BvXqFHjwYKZShWNxa/DCRG9F9JIBQA+xMyfJqIPAHiBmZ8D8ENE9E6Iv2ILwPfouVtE9A8gjSMAfICZt45bpkoNGzN/EtIdWwkmYXRupwjGBiaRF9fXYUKvQCSgIXaal8s1RWmef4o1QJEyL221wPpLAiIngas6peehLF30xGQenuMZ8J3PxvPTiE9xO+3gZiKhANdGIoxe2jsPALh6awPmqvhyOteljO07GZo78szhgeZPGyagkUp7K/VT9TqnnvRXUMZgq9bq+yAiINBA3qGU0TR1OQoRjOQTCNQZHYwDmNhyJg70L9g+AWI3ntFyFlGChq5bLXseZ5kfkVoxaHceTpqzF3cvABDOSDnrLoEzKpn+3XFnDMhyFqjGFlnOUgSj0HVCHBfLCtBl5ucBPF/Y9n5v/X2YoYkx84cAfGgpBVGcyMgDShmNnRiUhnnuQf2YD4w2cME5NANpxNqBdBg0TeyladGKw9lUr9ukSWrNoNwkPa456ps0fvxT6jmegXyAdAzjVQ4xvW8nq3htJJXk5X1pzK0pQ9da6LwqD9G5reMDd1JEe/LswVCel4YJKFabI7XmjFYSZsBWDq8iu8/UNvSGAKPvw1YWdQtQM5VOCAAm1ej4lEHavU/K2b6RBu6VMENLObMD61sUu0HzuVDKOYtKODsJc3QZnN1NV6Y5uyMdBXSt5Rq0zq1jcJbO7hBx3AUm5y9ULpQzjkMxaedcpyrO8siDeqxojRqPMOqGbRFkjGB/XAhDkIV1iu43OnilKdJwJZLOkpXAT3OUJ0B0CQ9pMglgxuxM0OOYpIeZM/Ze7r66HKtG088i9LPJ+Kc7yQquDkRDe+WeLNPXcqdz+45cpbmtmuleDNOXZzYjNZjGcS7tfakPAGkKrjBTkZiiKvWd5qbXyjKQm1Ipsw/unay7lLODRhdfjEQbsZz1gqHL/pFrbht6rqkAACAASURBVOkUZ75JWoxtqzqCxMdxORty4DizIR634lXH2dVt+TbTG8LnymuE9u0jcpbm3E1xlnGuzioo8TizmpnPWZaVmrKLok40WaNGjTOJ5UwKc/pwMj62jGGG4jeL7HtT4ZzqhrRhsNeRrvVrLc1c0Oi7vGA9LxeYTVVtszH4fhvndzuJB8FkKvOin2bkcnVFLvPDjmps14fruLonUv/gtmzr3pRStu6wczo7H83BGKRS3y4RJ06T4kT9Ns6vluXamy+9C9KfiZw/xjqk4XxouQZhfE3CBorqtWyAb9oMsN1WztrybOvRACuawLKrwdUtE596zvpZcypbx2ujVcfZ/i3xu3Vu3QfOLIwNrp7PGdJswq96VDADSZ1ocgEwg0axzD7lZnoV0lpRXknintz+Vk9jh9pr2IhkpNZ6IA1cN8szupaZNxO9bcDCoxGqRqxnmDZn7LRqQ46wp72+d2J5lhuDVdy5J/FgjVvynHnMU4poTztE+nnFmKgcAJAkYGvOZHbpVZKyypEW/htyMWTs8uPYEQueSZN5jY02bJGNf3PDewySFYn5ek0HzW+293G+oSmplLNeNsDQuRK0gSv0kgKew//UcLaWc3ZbOGvf8Tg7OCHOLCx3Pmd2DkTb0GUZEPJSGjYpWq2x1ahR4wyh9rEtCmaRYIbcnKGsKnUYaSqXlkF8V9b3zolJcHNtBY93xATYjsR865mhS1WdaTiB9W+nYERHzK1f5nyeh5TZm08yl/qApIq2Juhtm3jwYAXpljinO1tSxuaOapz7qQsPMGMNExjFudSPRQvgNMulfjop/Zl5Ko6tFGxcunDYsZcq7dkE0zqRIdfJEKhpFOl4xeYOOc7210XbubG6iost0dg2NZloP+ijq6Nn7Hybdm7TFOze/aLzIpw0Z6/t93LO7pZw1n8QnBVcC4n26SxJY+O6YatRo8ZZQ915sAiYgSQBGZKxigDMSDW2gXhZov0AjR15qY0t2ba10cWtFfF3XGqINN0I9vMklSxSMSMr/XNH8WEO6YWlfXGMIfJ5Ja0WMnRjDFtuApbbQyn/vZ0uGltyXGNbI/JtdPogyX0zY12mKaDOZrZd/EnipL4LE8gmta6JbT5saEeaurQ4NrzC/5RZFQ6XYMIL6DXqYwt0bGmjadDclm2je8rZ+S7urOiEM015BxtZEz1zujk7yBpTnG3vdnLOdk4pZxnL0y1BYZMY77phq1GjxpkCTSQ4OEs4IY0N4CwDJSnIpiEIrOammSUGARr7VqNRv81OE3c3RPrvdjTZYZTPrG17tGwPV0TFfA7SS7Zoxo+yCT+KYwwzwM2+HXthHoD4a+7FGpg7kPInuw20rUZ6oL24mlPNjJJ8HKEXzDnVm+YH4RalfpbND9C1w3nIG1JlwxD0kAnNTXvfJjjT4TzByGotKaICZwc7Tdw9p1PPteUdXHwIONvL2o6zraEs450mepaz/QJn4/R0cBbYe9U+tnk4IY2NhSQiQAfB2yhs+3GYUeqZpWry7AXY7msm1558bP1mE0OdJDg2+vH7nQg2Ov6QgfFHhb1WyrkDeozJStLPGti2c0kOpEEOdwJEe3JuaJNE+oOjXXiA52guOJ0nnM1F80Ye+tDySyRFucPae3352ESfMy0jjbWBGwaILGeaz9/shrinnNm8/P2sgVgnFbYNSuam0vKLf/KcWR+S5cyWy+ds23K2GyDan8HZODldnC3DFEVtitaoUeOsYXnhcKcOJ2aKIs3AlIF0dmunvluNbZy6/P062ADhvnHzZ+7Y2bOzhpP6VtqmsDMk5dL5qGEfZUhLxGGGXLoVy9NPm9gdi/QfaPkb+4RwoFrIKDdnAEk46MYAWlMj9YI33TIrl/r+MUAexuHDzheQcT4k0zmkvevbbVqeCc6cQ1zLn2QINH+/TVke7RMOBpOc7act56SPzSRnKefvd5mclSFD3mlQNEX7aRP7sZS7r5xF+4Swf7o4c7umOFtSgG7dK1qjRo2zBK47D44AzuTnghR1ME2aawNGUyeHTnMjHOhM506apo3cEW39Nao5ZZ6w8UMIFqFq2pU9eb3U/SeksJlEJkMH+lkDB4n4AdO+Jm4c5M9ln9NqPkjSKT+Me1/wJH2ZneAPz5kHf79+vHkyDd1njLsXueDdbLpsnn/UcTa0nBn0C5wNs8hxZt/ZSXPm516rytmeapjpgWxr9z3O4tPBmRv7W+RsSSZkbYouBBYyvUHW5KvqENXefjx2EtlgCEDj3fqxNBSjLMRYTT47q5X7WEG5WTOnNMtIPunfN3aOaFkO0gj7YykvDTVeb5jPTkSJbczz1DPO9HNpiLJJc0aXeQ/b7MrBJQ5p8gfD23PKGji/00AKBC5meS3hzMQ5Z6yZeS1ngzRynOXug9PHWV8nHS7lLD4Zzsq4sliYsyXgrPaKVvpyiOhlIvqUTk3/wuFn1KhR47SDWRq2Kr/DQETvIKI/I6IXiehHSvb/t0T0GSL6YyL6LSJ6nbcv1bblk0T03DKebRGN7a8w851KR7JINcp4SqJZByyn7Nad9B8zjErPYSJFG2WRk7K+1Ac0RommJaA1VYKJXBKHYzK99OQyZXIz+qTOIW3LGGIUy7oZSdnMmN1zUVpYZl53lCfB55ozBcyT/P5+Mod8lHbcqK9ZF1KOT3Bms4Ak1szGQpzNwoPgbDAWjc3nzEalPAycLQPLCPcgogDATwH4Zsh0ep8goueY2Z/N7o8AvIWZ+0T0XwP4XwF8h+4bMPObj10QD2fTc1ijRo1KYK72OwRvBfAiM7/EzGMAHwbwrsn78L9nZo1/wO9B5g89MVTV2BjAb5CkHPgZb2p7ByJ6BsAzANCirp6VO1fzK1lndS79rVQ0CWBi9YnonJajLJySum4uRJqW0kdBMaX0zOOsIxyTIQTjLESc6ixCWn6T5IVyYQK+dHdR6YWlB56YsOVoT8gZ5xpAwW9jc+QVTvALMLEkrzxWYwtihhnr+NE414ZmcZYyIVMt+zRylmtqp5MzNlhagAaDkFXvFd0suKGe9dqBKwCuevuuAfjaOdf6XgD/r/e/pddOAHyQmf9N1ULNQtWG7euZ+ToRXQTwMSL6U2b+Hf8AfchnAWAt2CyxD8s+FLsuC0rlB+SZPRM2Ux9nVRw1mj0rqTR+vI+trM4xnQVIUx2srOWnLHe+O2v5MLOlSkqbJYKZyyuJq0zTnOWxXLKg1OWIdMOXMqaas6rgLI9fewBYwKi9w8yzpuAsI6v00kT03ZCpPP+yt/kpbV/eAOC3iehTzPz56kWbRqU3yszXdXkLwK9DVM8aNWo8zFhe58E1AE96/58AcL14EBG9HcCPAngnsybsw0T78hKAjwP4quM9WIWGjYi6RNSz65Ap6P/k0CsfJsk8CUisEtL7+S80Y5Obn5iMTzoNyLwPQKS+/hhlfRuCeU6MLJv5/jjjaSe0Nfn932HnTD3E7HtOlFkxwVkmvywz8gOdbc5KL3hEzorr8845CXDF33x8AsAbiehpImoAeDeAid5NIvoqAD8DadRuedvPEVFT1zcBfD0Av9PhSKhiil4C8Os6HCcE8C+Z+d8d98Y1atR48FhGHBszJ0T0XgAfhUxr8SFm/jQRfQDAC8z8HID/DcAKgF/TtuQVZn4ngC8D8DNElEEUrQ8WelOPhEMbNlUP33TcG51lmJkivgJOeLxkjXIci7Nlo8THdmjIxxLAALJsOfdh5ucBPF/Y9n5v/e0zzvtdAF+5lEJ4qMeK1qjxqIKB4pwKZwUn17CZGe674hAe5O+WDZzXzxgN3gXD0KT/IVjWQLklwVDmJt+wvedsDvlm5mlqXoroqdPs3JO+/+WQXrVK0n8WXxMXqjmbfcHlcna/UI8VPSrIffWF7eRmrrIVg4N8rHyolSQ0qasUZk7kU9lnctSkhYZoinEDduUI3Lym9j8jDOWD1uB7ZAHchMMcFBrzWeWaUzlOAjSrHEXOvHIXOYOZ5swQV+KsDPeTs0DTM50Jzo6KumGrUaPG2UK1caAPI06mYSOVLL66XZB8TORSsmRhvswaIkUbqgFFlCHSAXxW6lozJ6DydDeLZoWwWkLV3nWrhUQa2dkwCaJA1vuRjr0MyUl9pwWodJ/IaU/+clK7IaI84aCzl3TMraFK4QATZqjTxOZ8zAtylkYEjiY5a5rE42g2Z8cxxpbBWSOU72pQxllRUztNnC1Ta6s1tho1apwpMNwkPmcNJ9uwmVzCO6lvfTShAVtNLdJlE+CmSLd2JHM3toOxk7LOP3JCYsZqeinSKa0iIJ7SQqwm2Q5itCNZv9eSsqUN457LamzwNDjnKzH5kqigLRyCUqd0Yd/hFyFdeOWZxVlAOWd22QCylryPlmpATRMvzNlRc68dlbOmSdCynDVLOLNchbnGdF85m9W5UKWTZyHUDdsC0MpBBLJEOKezDuYNDTL9aFL9mJIWQG352LqR5MhvmsR9jM4k9SpLlSQ3y0hYGHiO6EhztLZM7Mpoy8stnRi4HSBt6Ecc2QHMaooa494HJdNmnrNlmF2221nmjVzvkI+zYM7kFdTkFcVrzKi4zedMy22fLWkDpM/c0XfQDmK0SNYfJGf+fS1ntsFtBzFWLGftU8aZ36gVOQOW27jVpmiNGjXOHOqGbQEQgCDIf0AuZXQSXg4N0oZm8Gip9O8yGm3RglYiGSPbMWM0aFLa2qhxX24FqlIHC6rWBmYqh34AcpkhfBPKrhfL0zFj9BpDeTzVOJNOA2lTHez6nEHkvYPYarLeXBBkOxfUme05rGdqAcDUOEM5YY7U98I43DbfIW4KZbOcRQHSpmps+mxJhxEpZ6v6DnzOXNjHA+DM17JteayW3TFjdPUbc5y1HxBnZWbnYZwtw4KsA3Rr1KhxFlEH6C4KK1kC65/JpT4AZFGAtKUaW1el/0qGjY5I0bVoAADoBCMnZQOXPtoGWS5n5vB5sNc3zIgKDugWSbl6wRBrkWgrK11Z7q20EHdU+utzZg1NbBgFQCzrdr5VCgyYC74YNhOzEwElWgAwqQlYeD6cMqnv/ruODePK4TQ0x52WPwqQqUYTe5ytd+SZe6Fq2R5nUz6208SZ1TD1mzvoPSScBQGW5vSve0UXAcnLD8Mpc8Z+KGk7cB/RuCeHZGsxNtqSPXg1zM0a+zFGlMe2AWLWuJ4wr7Iso+K463FeIY2bnNmaNWMt4wg9bdjWO9Igb6/2EPckp36yLaWMmtpgDANQpK/eTiSdZbn4tBlr4blAZlQWAFOT6wIznM2FykKByflxPHmcaRm5oTNNNQ3ijpwbW85WE2zoM69HfX0f05yZEmfO/eTMNqw+Z1Z4ntPy766uzOcsnPyWHxhnZkmmKOakaHrIUZuiNWo8qqiWa+2hxAl1HhAoDIEwABekftaUZdI2iFdUY1uXt9tZHeJ86wAAcC5Uzc0MnPRvFE3Rklv7YQJVQwbccZSbEymspgG3bOj+hpodtlwdM3LlteW/vjbCeF2edXxPnjPalxIHwxAU2xzood4680yW/MmsYJ6lBfjHIMvKQwHKpD4gmpnei7zOHec2cJqa/E+6AcYFzlqrI/fMa6FoPl0zcqaf5axRomVbnCRnARW07BLONtv7AIAba6vVOLOTfi+DszJ4mhu5zjfrMrCcLUNlo7rzoEaNGmcQtca2AAhOW7MaW9ZSTa2jswStGIzOqQN6Q6Told4BNpsiPddUmrZM7AXo5lIfEL9M1ZABM2N/NoNZez0/hCD3501qAV0zduW92JLyb6we4LWNFgBgdE+evbGvIQTDAKQzOgU6SzxnWT6RiN6Hk7w8rvQ20t3vzrIaQeDpsL52UIyYd1qA5zdqiG+JoxCs65aztKOzuncNxuvaeXBOOLvYO8BmQzW2INfYch9bgTNUD/NYBmeRDdB1nQjTnJ1vSvkfCGcTD1zCWVFjs5wtS9G6v/MH3TfUGluNGo8q6jg22NmeXwDwKjP/9UMOFk2tEYHb2svUkWW8ovOFrhsMz6tvbUMk52PdXZyPRHr2jIZ70AitgsYWWeGFyZ414ORCCUQ7VOlf8LF1zQjrgTzDhpb/cncXWxsyv+poR5bRvmoSo9DNy+lmVmd2egh5SzdI2T6X9e9MBHiWPHPZuEPftwb1q3maGgDhzGpqbdXULGdr5DhrbUgv8GPdXZxviJbaCw7nLCCa4uykUJWzC1r+Cc62ZdnYewCckReGU9DcrO96WRk+ltUrSkTvAPCTEKX8Z5n5g4X9TQC/AOBrANwF8B3M/LLuex9krtEUwA8x80ePW55FNLa/A+CzAFYPPZIIaDaQtaPcjOkJMcN1qVzDC0B2UeKHHlvbk2VrFxuhNAyrgVQccURr5Thip8Esk8bus6bNhENa++NdxDwzIj3HOsJbsCMPRugZKe9mJM/yeLuDO+srAICrFxryzHuyDIcGRs0aO0EvODd1XWUx5GKm8tRB/sxGeuS8wZde+igqNGx+5w488zPtWiEk+4ZrytkmIVHOnljf1efcwWaoDZu+A3EflHf0+MzY9x0QnQrOHmt3pzgb7MsyGN1nzmxDZhuwUHmynC1LOCyhYVOl56cAfDNkKr5PENFzhUlZvhfAPWb+UiJ6N4B/BOA7iOgvQWa1+nIAjwP4TSL6C8x8rMydlbqgiOgJAP85gJ89zs1q1KhxJvFWAC8y80vMPAbwYQDvKhzzLgA/r+sfAfA2klb7XQA+zMwjZv4CgBexhHmLq2psPwHgfwDQm3UAET0D4BkAaDXWRFvrRhhrwONwXSTQ4KJImuFjCS5siqR8auUeAOBitIcNlf4dEs2gRQladtyhZ84A4ixedJzhovAd0tbCiErMm66R8m4EUv4LjQ6eWNkGAOxtNgEAO4M1ueY4gEls0OY0BcZJei/Thj1e01n7waFcMi4mHwOaj12cGq8bhS4IN1OXQdqOck2twNngcoIN5exJy1ljD2vBpPugRTEiTJugwOnkzJqklxq72F5pAwD2LlTgzHvvS+fMBU5Pcuc6d5aU4WMBU3STiF7w/j/LzM/q+hUAV7191wB8beF8d4xO17cD4Lxu/73CuVcql2oGDm3YiOivA7jFzH9ARN846zh9yGcBYLX7+BntRK5R4wyBsciQqjvM/JYZ+8ouUmwDZh1T5dyFUUVj+3oA7ySibwXQArBKRL/IzN8984zAIFlpIO6FGJ5TqX/JSn3RvlYu7eOpVZH6V1qi2WxGe07qW2naojSXtnp5fxiV76eRfb6vrRppxkl4L5NHIfW0IcqHVOkxvt/GlnsYyN5L0Q76HfHP9Ndl+efqo+mPV0CpHWOoFyP/waykJ8AF8ur4RH8IViHUAMzTvpfA5L6ewpAgboQuYDrt5H61kfpBc01NnrNz6QCvWxPOnlDOLka7TuPxOWs5v6g+kl0+IM7sHeZxNizh7M/Gsq8fd8s5s2ZEFc6Y846E0oea1rJdkk8b6B4FyJphnrz0uFiOCnINwJPe/ycAXJ9xzDUiCgGsAdiqeO7CqDJh8vsAvA8AVGP7u3MbNQBZQBidizBcN65y9B/XD+pxMVue3tjC0927AIBLkTiizwf7WFWHbteLgyozZ/zlsuGnxckrXebio4rmTQspYu0FtJVlHAYYslSKkTXtUll+ITXoc0cuQvYDJpeVNtJB58HAwIzkuvlIBTVr0iyPhfLNmmLD5iVItIPZ3SiQRpD3fNrOnbUAgwuTnEWWs/NbeKorDdvlxg4AYD04mOKsRWmp28BfLhuHcab5I11DeBhnA+1AGZ6X9/NSSsvhrDC2dEIY+UsvvRfgJY9ohkjawdIatiX1in4CwBuJ6GkAr0I6A76rcMxzAN4D4D8A+DYAv83MTETPAfiXRPRPIJ0HbwTwH49boDqOrUaNRxlLaNjUZ/ZeAB+F9Pd+iJk/TUQfAPACMz8H4OcA/AsiehGiqb1bz/00Ef0qgM8ASAD8wHF7RIEFGzZm/jiAjx92XBYB/YsGg03CUE3P9mVxqj99fgsA8IaVO7jSFOl/IRSNbT04QEfNmaYb25enko5U+hrPhDmOOVPEPPMGbCbCCKQ8ghTsOjhsCEHKBrF20cd2rGUv7+N/Wa91EEi8VBYZpJr5pKnqRXRgEA40bc5YTRjVAijNJsIO3LI4s1SQzwVqpb4d+5m2zESMGgAMNsm5C5qXxcR8elM06zes3MFTTeHPcWb6U5xF4FITVP4/JJy1VEPSEBJmwst62YNANLelcOajhDObDSfPimMwXgnc3BPHxpK84cz8PIDnC9ve760PAfzNGef+OIAfX05JBLXGVqPGIwriOm3RQsgawP5TwPjSGOsXRFN73bpoZ69fEen/VHMLF0IJHbAhEl0ao+sSAubO54gmNbWTDhcAcm3CT0Ht7ku5DwdQ6a+Bn6kes2qGyOwMSnbYSju/vp01/eVQlgetrhtHG99VLWDHIDqQa4QDDXRVLcDE+ThFpCVfpzcjlk1znWma66Qt++KucZracFOuMb4Yo3dJ+Hj9OeVMfaFPtrZwKbK+NZt9Zeg4a3ijDCxnkb6rB81Z5riTxVE4s3PHvhxtAKjOmYktZwyT2tAPvaivZetnlQX5bFmZcufy4HUJo3VC2qj8WuajTjRZo0aNs4ZaY1sEzQz8hj6unNvDUz2R+k92bG+ahAk8Fu5MSH0A6JjY+T38ntBij9pRhuIEM+ZpTEsmQjlsyI4Vt5M9bvYLSexFkEJ62+xbtrnBAjCamjq7FcjyWmsdt1ck/jleleDQ8ZZBY1c0gWhfzg2H+uzjfGiP8TQ2O7+lLWoWkptSLpXEFXkevFVgtKHaxKYOb9vYxZM94eh1HfGn2R7QS9EOzqt2vW5sttycM80ajsYp5Mx+TzaFd1XOXFna8zkbr3mc7ShnB8r3UHO6jdkbbzr9XGwT83qcJaox+pyNz6fImtPnHwl1w1YdzSjBX7h8C1c623i8mVcKAG5kwbrpTzRoANClZKJyANJhYGi6chQxz/k8q4LYfWUVpfzY6TgpAIi8MEM/nACahnrqOsgm5rcEgLVoiFdbUqlu9qSy7J5vY7Sj40v31LxRMycY5Zc3tg/JL4etJA0g1UqQrEjZkp46s9fGWF/VBAQ9cQtc6ezgcS9GDYBrzDa8cByfs+gsc2aH0lI6l7NbqzLGdOd8p4QzbeSHBH1tIDuclPMEG6z3Sn3OusrZqs4dsTbG+dU+bje8/EhHRe1jq1GjxplE3bBVRzuI8eVrN3Ap2s01NDU7bTBkz3M6+6aML/UBkbCzHM/Lmi28/NqTYQSzAkABADQdTiCiUMWyagFujk1PY7OhEr1giA1N2HixLRrbndUVbJ2X0ILdviZAHGgWjmEIxFJGOzM5mJwI5lC/2IhhdN7MVlvKcV5nlTrf7uOCJsbc1NQ9Fxu72NCxn5azdf3fpfGEpgaIideqOcOdjmhsN3s93Dsv9qPlbDgQDS4b5nOT5pzBdRrknGUIlLO2crba9jhr7uOLYbzYy5mBMpP4LODkvrIaNWrUeEA4GR8bxfiS1q2J4TZWynVJJJCMJ5x2OheDcAMcfWzhPD9N2XGznNKAaAFlE4hA9xZDQHRAjK6rf8So9KfMm5lc3kfXjLCmGtJmJNrTVquLvRWR+ttj0QL2E3G+HMQNjHQimETnqMyYXECqDSdphgm6kdxjRef9XG/odHlh303Acs7mwTMDlzCy54ZKzeesGIT7KHN2ubmD7US07CJngyRynMVZHvgbKFe2c6kdxmirRrais9WvR/kUh2thHx8LlqOx1aboAmhQgivRFro09mYF0rGfsGljsrlOZ783raxyyP/FetMOw2GVpbTXDdDKYk0cjTfizJloxp2Xz7HZ0I+5xXZ8Zewaf2sCboZ72MukYes3pXL0MzFrRlmIkabPsdHxGSi/l5ataRI0jVTSjlbIXpALm64zq/KUQ3lDJuflo0BqzoD5nB2ETfS1y3Kvqe4D7RXopw3E2vWZZPkoFMuVFUpNk7gJpztuHlTL3QAdM0KD6s6Deag7D2rUeJRRN2zVEVCG8+YALUqctG8Uct/7IwqCEqfzUdNG30/Mi3EDGWR2PKCNybNmDmUIrONaNY2IU6cpDdUEHHLksk0MM1mOVeIPuYFM4wTSEldp4BI9pvkcn86UyjWOFuXrskxcWu9iTGGE3LQs6yiYx9lpQ0B0bM5anCzEWcyh09iqcuZzZZcNShHZTo7jom7YatSocZZAOLu9oiejsYGxZkZoUObkUtE3Y5BrahH5juXZTucHjbJMEhZ+RgnoUVbqRy44VHcxO39KDKsVZYg1W0tLJ6ccY+QkfGx0hnL9n8G4MY4pl0h/bzIVK/3djE0uVCN1M7Xns0ll3lyccq0yzqp2FFg8KC27LARE/mfH5iymbEJ7A8o5szzFHLp1i5SN48oiQuq22RCThvtOUkSUTZ1zJNQ+tsUQgNE1mVYEQRWzEzidlaOIWZUFwIRz2lUXN2lwPqQq9dYBSfVjK0ymH/OYDTJthFKIyZObn+QqSVbSsBmvYXOxWFS4J2VenBZ0W80ZcPKcWWRsHFcWVTgLlmVD1g1bjRo1zhzqhq06iAgtjW8qM10sysYRLkPq267/RUMIqo4/tCgzTafCCoASJ3V+rpW8Kdilrx5zLp3t+HZ3NT0mBeWpdea8noB4SivL9+VR98Ecc9MdX2J2Tjyzh5qzHJnHWV7uae4Czy4scuamVoBo1Mty0NSmaI0aNc4eHtWGjYhaAH4HQFOP/wgz/9jcc5AHbhbHDJb5Zsr25duO7qNZVJofFX4gaL4tH6dYdFJLaGiuCbhr2BmVdFvm7U/ZagkWuee37CnLJLr1nU1qW3bftJY2z5/mP+f0vWvOgBmcKbISVcl/k8WJ4n2NGvDmID0O+P70ihLRBoBfAfB6AC8D+HZmvlc45s0AfhrAKuSV/Tgz/4ru++cA/jKAHT38e5j5k/PuWUWjHQH4JmZ+E4A3A3gH1mAl2gAAIABJREFUEX1dtUeqUaPGqQZX/B0PPwLgt5j5jQB+S/8X0QfwXzLzlwN4B4CfIKJ1b/9/z8xv1t/cRg2oNv0eA9jXv5H+5j4qIQ/hKEr0w3wz+b7T0ZtWFXN9NxYlPpzck5N5UtmT7M53Y3vWpl99dY2toG2VaGdl5a45U5wQZ/O0i1mcLetN3ycf27sAfKOu/zxkQqi/5x/AzH/urV8nolsALgDYPsoNK/nYiCgA8AcAvhTATzHz75cc8wyAZwDgySsB/PGCPs5SxSiD/wxlpo5D0QbwQjbKKk6+7+hfYlkqoTI+qpib+b6aMznOnn9/OKNlvffqRdskohe8/88y87MVz73EzDcAgJlvENHFeQcT0VsBNAB83tv840T0fqjGx8yjedeo1LDpPH9vVtXw14noK5j5TwrHPAvgWQD46jc1z6hLskaNM4TFzMw7zPyWWTuJ6DcBPFay60cXKRIRXQbwLwC8h9k5XN8H4DVIY/csRNv7wLzrLDqv6DYRfRxiA//JvGMXGTmwDKlfNhpgGThO2eaNVJhr8igCBM4BnZ+XIys1QmdcvwRlGnXVc08jZ8so00lzZlHG3XE4OwoIyzNFmfntM+9DdJOILqu2dhnArRnHrQL4fwD8T8z8e961b+jqiIj+GYC/e1h5Dn2TRHTBOvGIqA3g7QD+9LDzatSocfph5xY97HdMPAfgPbr+HgD/dqocRA0Avw7gF5j51wr7LuuSAPwNHKJUAdU0tssAfl79bAbArzLz/z3vBMLh0mcR6boMyT5LuzmsnPPuXfUZisdV0gZQ7q/xEUwFBSyORcbi3m/OFkWVez5oziyOw51Zpmfz/tD0QQC/SkTfC+AV6IzwRPQWAN/PzN8H4NsBfAOA80T0PXqeDev4JSK6AGlaPgng+w+7YZVe0T8G8FWLP4tgUQoWqRDzTLEiiuaBm72o4jXKP+Kyj/3w5y07pmrFuR84arWpyt0ivB0HjxJnR8Z9aNiY+S6At5VsfwHA9+n6LwL4xRnnf9Oi96xHHtSo8aiizu6xOJZptsyT8EVNrOp5Wclpc6VtiaP4MI1gkXdw2LHLNO2WGaIxq1wnwVkZas6Oibphq1GjxllDnWjyBDBPopVJ7qKkLz2mLMp7jobgR9+n3qi+qcDIkpxn/lcxTxNYZhjCg8b94gyYzVvN2fJQm6JLRlkFmVUxitv9ilD8+Cf2VSkI+2aIX5biUBaeX3EKos+vNNPR7Kf/gy/DIpwV9y+bs/LYsJqzhbCccaCnErUpWqPGo4y6YVsODpP6VtpPbNNzMpcGZlrCT6eIqSj9PeQpfIDYzU6kIyiYnUZg3ODmVI+nabNnjskjE/k+PBpAVc78fUXO3LHgpXA2lTQTNWeLYpkjD04bao2tRo1HGFQWHnAGcN8atqLUnyXxiz6ZjPNJNMokvd0We5fPZ2+yQbizJa2feSHwJs7ItTfrTPbmmHQJGG2ZyZt0I09QWAwCLvPhLKoFHGXG9KMmbzxpzvy050XOZPvsLBlBQdU4jDM7D6rjhwraG3AmOFsItY+tRo0aZxG1KXpEVJH6vsQv+tEy5BNlFLWzFIRY/ST53I35fJv+tHTFac8sAnA+VR3nGpudZ9Potok5Up3mps9E+TyU+W2ymT6co2gBR5H6xXOragHL4MweF7vzoP/zSWhifQ8Z59PS+dwVOfOnnPOnF7TLeZxNT51HE9qbfdKHlbMjo27YFodfQaZCNrwwDt+EsQ5gv2LYhmysH52tEDEbNzntkCM9JoD9nP3JhfP72lmW7AebTz7rTyhsJ6ptkVTNBmf5pMKwS3bLtGDyzDJPAQCUnVrn9P3kbKyiIeYAMU9OLpyymeDNh89ZpOZkg9KJiaCB+Zz5M0wZ11g/nJwdB7XGVqNGjbOHumFbDBkmTRYL35SxUj9Wdds3O63EH7Jx0n5Y0M6GWZSvT2xr6DVyjaA4W7pvyhQlfcvETlNrGZnNu0XeNj2+pcdHyNCyos8+3yzzFBAtoMTEmYejzrtZFYdpanZ7Fc5ikMdVNc6sxpZzFs50H0SUItAytozlJHb8dc3IbStyZjW4FlL3zHZu1YeNs2OD6yFVNWrUOGOo49iOiHlO55h96Z/7ZoacdwIAIvEPVAPrZ00AwIEu97KWW99PW7ItbWGQigweZXYZIiuEftjQgdCkaBqR5u1ApHsvGKKjmtpKMAQArJoBOqoJWI2gS1abS5A6n5zny1GpP7djwdMCljlO0cc8B/QiHQWzOMu161xLs5wNlYM+C0+7acvxuJcJZ/20iX4mx1vuEg4QZ5PJGJ2WTYyGcmZ56gRjx5XdNo+zmBI0zwhnx8KccdQPM06kYWN4vWYllQMQ562tHEPPhOlnoW6T5W7Wwl7WBpBXhO20AwC4l3RxL5b1rXFXjo9b2BtLxRkmco1RHCLNjCsbkMdBhUGGZiiVpBvJR99rDHGuMQAArEd9AMC5sI+NUGYhXA9k2zCQY7o0Rs8M9Tm1siDN71ZaWTwTB5hwTs+rLCdl3hzWoAGTnNneziHnZmdfG7EDbjjOdlXg7KTCz1bSxb1EONuJ5ZidcQt7sRw3iFUYJSHSTJ6fVSgZow1QkKGlnHWUs9XGEGuR5UyWG+EB1pSrImc9M0QXY33Oh4ezZWcmrjW2GjVqnC08ygG6RPQkgF+ATK2VQeYT/MkqF/fDA3ypL//ZaWpDlYAjDrDLk2bKdtrFVrICALijy9vjHgDgtUEPdway7V5fpH//oImsL49FI53gdkyg1Hbf60KFJ4eMrKWSuKVhAp0Yq13RwDbaIukvtPdxqbkLALjY2AMAbIayXA/6ztFtTZ3UxABEq8jKtIAKzul54QQpZ5U0gEVNmrKQjsM4O9BOAKulbacd3E2Vs1i4spzdGq3gtnK2dSCaW7/fRDLQT9FxZqD+/hz6KrKIwQ0pCLXloGYnRq8jnG12DmTZ2selpnB0sSHcbWSybxj0EQeybrVt4fB0c7Zs3I/OAyLaAPArAF4P4GUA387M90qOSwF8Sv++wszv1O1PA/gwgA0AfwjgbzHzeN49q+jGCYD/jpm/DMDXAfgBIvpLVR6oRo0apxuUVfsdEz8C4LeY+Y3QCY9nHDdg5jfr753e9n8E4J/q+fcAfO9hN6wymcsNAHYW5z0i+iyAKwA+M+88K/nLpD4gfjUr9fvqT9vLWthVTW1LJf7NeA03xmsAgOsDWb66L8s7OysYb4uGF+7INZp7BHWjINC5ok3MjhwqaGxZYJA15E/Slmsk3Qbu9kSbuLsu5XhtrYc7q+In2u7Ivr2WlPVyY9uFKaTBfv4S9B5d1QKMpwUUndNl84UeZ2KQqppamS+0GNJRxtlIn3eXm87naTm7nfTw2ki5Gq4DAK5bzva6GGzLewuUs2jPoGk5E+UJJs4HaDvOrNIdGqQNfc6WhomsNHGnwNmN1VXcWRVNzXL2eGtbniWKJoKBAWDVDB8azpZiQTLuV+fBuwB8o67/PICPQyY9PhQ65d43Afgu7/y/D+Cn5523EAtE9HrIjFW/X7LvGSJ6gYheuHM3Le6uUaPGKcQC84pu2vqtv2cWuM0lO+mxLi/OOK6l1/49Ivobuu08gG1mto6JaxDFai4qdx4Q0QqAfwXgv2Hm3eJ+Zn4WMv08vupNDbaSvyj18940M6GpAcB21sGtZBUAcGMskv7a8Bxe2T8n69si9ft3RPpGd0Os3BPxGe1psGWfEaojyCS6jDnv7NJysGZx4BDIQpXcOsYmaRHiFR2WtSpl2z/fwOcuaDk3tJd2VbTFUSfEqCF+ptT5WgwA8c/l4sNyMx1WkFEeTuDPS7lob9si/plZabozF3A7yVkMmtDUAPGn3VbObsbCz7XROby8vyHrO7Jt745ou+GdCN1t5Uy/ognONHbEpAyy8rGgKWUhQTvPC5zpUK017W0938TepqzfK3DWbzcQNzUYOAz1fTycnB0L1RW2O8z8llk7ieg3IX74In50gdI8xczXiegNAH6biD4FYKqtQYVSV2rYiCiCNGq/xMz/umopM86TCvrhAYB1Oos9sZ1JQ3UrWcWrY2nEvtA/DwB4efc8rt+RRg6vyUfZvS0fRXOL0djTEIC+fjDDDMFY1kkbNso4V7kLlYQNgUNtjLSBy5oGya42bDuybOwaDPekYbu5r7FZj2nk/LnIVSofpvARBsZ+8GmeWseZNRIko4XSMi4WTnAYZnUWyD4dSeCN/Sxy1s/CiQYNAG4nq7g2lkbslYEsX9o7j1eVs/SmNCidW2Wc6bhMjzOTZLZg07nCrDAilHIWd3Q8quVsx2C4K5y9diCfev+ixtedi5z7IGvmjczDxNlxscwAXWZ++8z7EN0kosvMfENndb814xrXdfkSEX0cYh3+KwDrRBSq1vYEgOuHledQU1Rt3J8D8Flm/ieHHV+jRo2HBCzCo8rvmHgOwHt0/T0A/m3xACI6R0RNXd8E8PUAPsPMDODfA/i2eecXUUVj+3oAfwvAp4jok7rtf2Tm5+edlGrGi+mxnyIlDzhyUv+uhnHcGK87Te1z2xcAAK+9to7GdZGy7Zsi+VpbIu0aexnCA00WOJSlSTJQLOtkpT/z7JzTBuBAtS2rBUQGYV+2hQeyjPoGYV/NDg1N2B+LCfZiEiDRAGCs5pe24xldJhG2DmZ226wmGXimzrJHFs4aXTArtGMeZ3uFzp0b8TpeHkxydv3WOsLrotmtWM7uysUaexkiy9lAOYtTx5XjLMumObMvhghs9E8JZ5HjLEDUl/sP+hqSMhKCPpcaJCXjh4ucRaxlPGWc8bIC0O5PHNsHAfwqEX0vgFcA/E0AIKK3APh+Zv4+AF8G4GeIKIO8zg8ys+2g/HsAPkxE/wuAP4IoWnNRpVf0/wNmBObUqFHjocb9GHnAzHcBvK1k+wsAvk/XfxfAV844/yUAb13knic3pIpl+E1R6tthN3tZ22ls1un8ymADn9/ZBAC8dkN8bc1rETo35BpOU9vVrBr7CcxIu+XHqqXFKZCo1znLJWxRnWaTt9VkpX8gS4oCmJGU14zkFQXjEEGs2+wysc+0gpf1WqGxWUMyl23CLq2PpmEyl9TS5v+KwZNOaSlkaQCo/F+OrClOuuJzNvbG6wIyRtdq15azLw7OO86u3xS/WnStic51KV/7boGzgwRmUMKZ5Sr1NDaF5c5xRgSyzvpQykahmeZsFCIYFzhL7bP08JIdqqUcRJROceZzd9o4Ozakoi7nWqcMJ9awxWCkkJ40YDp9zW7Wwh3tTbs2kkbspb3zuHFLKkfjVTmucx3o3JmsHOG+uLXNIIEZa69VLEtKSipJibPVfWJEnjdYG7jYgCIdvRBbcykD6QBtGxFPbpC2wYGRXr8vBJoWJ4jdoHo7ENulQOIYke0tRl5xJp3SyzFvDotZs3DJIpFzNnLCSDsMso4zQa8OpaPg5f0NvHpbOAtfleM6r1LO2Y6mC9rXxmwYwwzncGa58jt8FKWcxdYkDUDaePlmbRXOXg51rGgQu4QIxbRVDZOeGs4yXpoh+ugOqapRo8bZRT0IfhEwI1WTJk8/lGfrAICtZAU3Y9HYXjkQje3Vu2sI1Oncec2aMika2yJFwwNZBn2RojROQCMNSkjV/EyzfN3GrDFPq9y+KWrNGu1EIGNybU/NTUoZoV6P0slrsQmQqUm03xCN5uVGjF4omlpPw+knkiKyNXXylDluDgWXrLI8VQ6wnLTUE/MVeB0GszjbTru4pR0m1/qipV29uw5zXfZ3buTmZ9Nyptp10JcljeJJTQ0Q10EZZzNARC70w2rZlKTOLLXcBUnmjV4ocEYBMuV7L5JxrK8cgzOLbEZ6I9m3HM6WiXr6vRo1apwtMGpTdBEwcl9N0bdmM0BsJSt4dSBS/1WNTk9utdG7Oe10ntLUhir9x3HeUZDIMZxNa2wAZjtJDblRCFbSU2Ccz4d0ycxu0g+Lhs0QEgBZZIM95ZXe6fRwrTWZ080mQuyaUe7DUS0g9gJAIzdrOR/ZZzMv1fdE7jXPtyblIC+8Q8JsLGf3ki5u6BjQq7uyHN/qoHtrkrPmTppragcFzuLEaWwTnFl+2PO1zdCyGcg1No3gL+Usm54Oxr5bnzM73vR2t4drTeGoyFnPDOZyZt+3PwJhUVT1hy5La5MA3bPZsp1c54GaNNacsZlubU/orbiHG30xa7a3xInbvhWguaU9h17PZzDQSuE3aIBUEFs5XG9aOl1J5qnbKXKzVCsJs3GNoq0khFy4GVepdGhQYNCKtDFq6TCdXgPXe/J851uSHudclCc97Bmb9z/Pxe/PoARIL5r9oAPYGbdy82bR3rZixPrEyBC2y3wWKZepWJNF3hr3cP1AGrT/v71rjZHkusrfuVXV3dPz3PV4109hI0VOogAJskJCJB62IxkREUAhsiUig0D+wyO8ScgfhEAKEsLkBwqsnJeEZRMtAawk4MjBwB8UxSQRdjARlmMlix/r8e7svPtRdfhxz711q7q6unqmZvox95NGNV1dXXVvfVX3nHvueVy9oqfczcuB9VFrXksXd+zUs07OzOjrmBEqcybb0Axsiixn/QXDWRMvL+tp6XUt3U+XM5OFt5yzNBqhDs7yMN4GZVP18U5Yz2mmDX4q6uFxiuE1tjGRwEh/idET6X9N0kJvdJZweVtL/WBDT3maV4CmxBGGe+JmcdCzCwQZqQ8A/T5YpH/q2pGkmoBtzAixJMZeIkdrMEZsJYsHSN0NzNSVAr0NQoVoR2IhJcC7u6Gwc0b89EQLONdKE1QuK4l7NYkpuY+EzLRDms1spzZG2xo11amS0MaNMkinOhpay5ZYWMOZiQvtLuGVbd0Xel1z1roCNLeynAUHfZC44VAv69oxjDOjoVlNpIyzGOmiQS2cSYzphsL2mp49vLqk+3lj6xoAYD1sTQ9nqMk05m1sHh4e84da4kCnEscysCUgdFmhC2UN0FvWAC3Sf38JO5uS+UFSDzW33NhPkfid2JH2uYWCfj+V+rJgoF07Chxzh2kASqW/NZKVBy0g3E8jFCjnaqCCAKGsJES7JrME0LkiRumzUshkUYqYNNq2MIzRjtro2RRPkbOIkNpuhqNM4rtG58FFhDSDh4ky6CLV2IyN7Upft//y/jK2r2nOWldMPxnRrtiUJKKADvpAgXYNnDBnCVmXkjLOGqK5NTaDAc5eW9SzinON7RPlrGjRQH+u2Szmp6IeHh5zBa4l7fdU4hhtbIQeBziQepGm7uemlMvb2GtDXdOSr5FJOGhsayL9uz3HkTO3muY441rbTBwPSPrSFaQ4HpDmSFJ56moB3Bc7jbHXKBPzGID39W9NhonGtkJDbDe7ksftithvrrUWsBdpbair9KpbjxVix83DbIetniVI836NQtFqaPqdOR9JOwLHNcdwprW01/cWgU3NZ0ObnhDtJQj3He0a2q52nJyRe+wozsy+Is4kttTYBqeJs4Hf5jirL6TKa2yVod099MJBaojWL8TVrmQy3W0h2hJ3iR0xsu4nUPblMKmH0qB2Ng90YvzUkuzLAegHPE/WCDtC6igu51epT1TqL8XpNczLFJvY0j6oZwKvTRLFBNG2+PBtSd8PpNbm4oKd5pn70+MAsbgR5A3S+v9sQsOj+EsZxI67R2wHttAWOTbC6GpXC6OtAs7C/QRkEnv20oHtODlzsw0dibNOjrP9BNGO5mx/W9+DaeasFsznuOanoh4epxk0ymNgRnFMGhshBqHLgZX+xslzu6u3vd0Glna0BAsPxPjcSZNEZlIP2TjQ7JbjJJXSToRAobQfloqZVOpq4GpuIuFZrkVA6vluHECdFDs2C4ik4gkPAkQmMeW23m7t675v9ZtWgzXuMDHITgdjZ06VVHQZqIJ8xbAYrqYmiwfOVHRPSkFtCWfd3QYWd0Vj2085U3nXjiLOzD2ug7MkVdlq5Ww/sIs/4ZbUdjhBzobVoMhu9XW5jhRIqS/x3MFrbB4epxQE9g664yIGoYfAhufsS16sra62U9Be4NT/lOXxbpxNDQ0AcZLaadxcXeaYfDYIV/JXKZjhHiOagM4V6Nhu5PxktI+izBJxmrdNb9nWyAwlPfXBvpb4e/0GOokJqTL2mhAxssWtY2Yb2zguytw8imA1NwSOjU1CqhzOQiniFHQNZ0k1zhx7WmkGj6qcmSpPjrZdyFmSu34NnJlnum7OXJRxFjP5xYMRGDmwEdEnAbwHwGVmfkuVkzKAJLd4sC/TmoOu/hzsK1tyzVaVijlNEuhmUk2yg5Z9IVz/p0wDcitsJYZoUgUPYcLpFMfuS2yGXbjXN1ub5si8JIl9+cMDeUkOpFZCr2mnebbQMgiJZHRNnCZlYxDHT2RYVNEocbaxuaZ0uIizvY7eqn2FwJTJk75RP7HzsEqcme8zDeLxOTPHOwNcIWc2qN74u1XjLOgM58zcn0lxVitOYGAjorMA/hbAbQBeBPB+Zr6aO+bHATzk7HojgPuY+R+I6NMAfhSArMXjF5j5GyhBlXv+aQD3VjjOw8NjlmBsbFX+joYPAfgyM78BwJflc7YpzE8x81uZ+a3Qld/3AHzJOeR3zfejBjWgWjGXf5cK8GMhBiF2MkWYqWhH3CJUF5AcfrYwLsXpNCVjYLZSd/AOF05BzXcVwkU44VRry2sBbjuCIL0GZbUQcgziNrNEn9OCzWa20tXn7cRhqqk5mpKZDqZbRlTSdjdrRBUkhZpA7polnAUOZ6ZvFCfpFLAKZyNSSU0LZ0FHzjXFnNWBE1oVfS+AH5P/PwPgX6ErTw3D+wD8EzPvHfaCdVcN8/DwmBmwFjRV/o6G88z8MgDI9tyI4+8D8Ghu358Q0X8R0UOm/mgZals8IKIHATwIAOdvChFD6QWERF+iK0U0+lJwI+iRrfpt7WqJY/dwbSEDzpslUqZihe3sT0SKF9nb7EFpDqzUhpO2y1ZSMvamxNHYjGbak8wfcYCO3BejBQCp1D8puBnOY5FxhZz1xW2lS6mmNoqzPBzOBhYNpomzeHY4OzIY4wxa60T0tPP5AjNfMB+I6EkANxT87iPjNEkqxX8fgCec3R8G8AqABoAL0NreH5Wdp7aBTTp5AQDu+P4WA9ogbdRmU1A4kW3YT+PUbEGJUdIhP00ZMsCVTmeKpi5FxxhfqaDgwS2ZXpHzctv+mRdHQrJ6cWCN9e6AMgnkr1vImZSsC2PYF38kZ8NMBIfhrArq4syYRZLp5qy28a26TNlg5juHfcnM9wz7joheJaIbmfllGbgul1zn/QD+nplNfgaj5QFAh4g+BeB3RjXWT0U9PE4xiLnS3xHxOIAH5P8HAPxjybH3IzcNlcEQpAN+fxrAs6MuOHJgI6JHAfwHgDuI6JKUqS8Hp+4DMSvErLQmwMoKeGIA+b9MoB2PjPEcC5wcbspTlobZPWdec3H6Rbn+MmuP9TqNwOMgZi4IjFejOUswGc7G5K4WzpJsf6eRs1pwMja2jwJ4NxH9L4B3y2cQ0Z1E9LA5SBYpbwXwb7nfP0JEzwB4BsA6gD8edcEqq6L3V2y8h4fHLMH15TvWy/DrAO4u2P80gF92Pr8I4OaC4+4a95onFlKl5jHxU5mdzmM64TnL4rRGHnh4eMwx/MA2BijV0AJTDZvM8rwssSukebWM2ULBicOs2ZZxSElNRXF/RecyZd1sLjDYftl+mkU7ldhwm0nAZp1wHmpXozachUqyZDicWa48ZyeKIs6ODEa9NtEpwrFqbIoSRLJ+HioZ4GTLIcAmpYxZnidyanw6D2f+xbH1JdPc9y7Ipqo5JGmkxntJiWxf7HIMEYy7k+kfhzLYq8QOGoGtQTmZByx/XUXpCzzAWQAkQQFnRZgBzixmkLN6RAgfakFtFuCnoh4epxWME1k8mASOZWAjMAIkCMB2itOQXPNRpKVer8EQR+5Uc1PkpHUu0Qhsmue0XgEb6UlqbClU6r1uDyKnNkJuS4Pt5pCsdiPZbsChbmMjiK0ma7ZK7tdJIiDAXNLVQiLhKs9ZfxhnJuuJ4WVGOOMg215Ws8VZLfA2Ng8Pj7mDH9jGQwBGRH20SEv9hUBHSDQj/bnTYCQNkY6h0XIUOF+7MwjSNOGUrf7NGTuJ/O/GAjpSfZjtJiP5iwzMytFGclpCRhuwUl/CkAKFJJKMDxKySy3dj4Wwh6bRhuT+BJRK/3RbjnErHqmiRYP8NR27aFNsSilnqSZTyplS6b08JGeVc+iVceZewxxexlk0W5wdHbU4304ljmkqqo3Q+iXJDmztht5ebSfotyR3fFMesEgBoWQ/DdKXheR/Ntn8TNLAGGlGVJM/XyGNGXSmN6VTl/zLoZwpDBW8REoe30C2SgGh/B+al4QQy8At5R4QNvW9WIo6aEmUdYOMQZrtyrHbmkDMxMEhzcUB0cDCl7vgF5hrislACyPdtoVA51synG22E8Qt3U8jlDhQNjbTcqbUIGfuekEFzkZONatw5g6wR+SsHXZPlDPYuFQz7YXdBlTT4gGjPKHEDMNPRT08TjO8xjYeAjAixFbKtSXb4kpTJ5V/qR2j385qbEmowCI9yUjTQAF9I4Fln8m8EAS2IlEGdorjyNG8cbpwCqN/R0TZKajdV7BoAK2psJX6sm0o9HVJSvTbur3tBX0PFsOuvR+pQbo/YIg+SmUqM+Vxi/TGKLhX5lpybZezJcm2uNTQW2730W9LKiPDWUOBI72vkDPXpABgaDr3KpwZuNzVyVmk0BdNLc/ZctQ5Uc7KoN096hiQTiakahLwGpuHx2kFA+z92KqDwIgoQYNia68x0n+loTW25mIXvWUpaLJg7BoKwb5I/wNj/whAsTTTViFybG2c1bzI0eAyXyUFGpp7LuSN35TdFwRWYySjGYSp3YZFI0maotG0CL1F/dv+irbTrC7ovi+HB2grfT9ajhZgDfkZ+3pWAziKRmBsPom1AXGqqRm3nALO1pq6nFhzsYvekuapp2t5AAAIuUlEQVSsbzhrKgRGYwsczvpGu84moaTAMX2bfW4leHdBYVjkwXFxtkDo5zhbbul7UMRZAJ4IZ7VGu/rIg+ogOKuiMq1ZlrpmZxs6jfnq4j4ur7QBAL0l82AphHv6fyUPG/UCIJDQnlA3l4uyoLorYc4LY3cNeRoy4TeOH5Z9Sc1LFKQGcZh2yEvCUQi2A5qUZmsr9Jblp8v6Hqy19ACxEh7Y+xIhNUQHOUO02+T8atq4q2tFCIigTHFfO8ClnLVlYFuNpN2LB9go4+zA4UzujRE07JgR7AqpucdEac4vd0FhCIZxlk5Fx+MsbqacdVfkp8LZ2QX9vBZxFlE8Ec7qizyAt7F5eHjMGYbZO+cAx6axReIPZaY1Ro1fk4q76+1dvLamVZrumraydzcVwl0j/XXTqBmn1YRymhj3AQSyT2wFDDhT1lRClkq4fAFk1+gcpFMXOFMtAEAkWkAzQCyuK/22PkdvmdBd0+1YXdZ9vq65qz+He860RrQAShDlfKEC0KFdBlyYqVBiZ3upb5S4bSGCcfdIOVtWWsteE41tvb2LDcPZqmg5m47G1hHOemkxYmNGMNyx88il2pwa4CzT6wI+89ENFASH5qy3OMjZypLu87RxprjmqajX2Dw8POYLXOxVMAc4xlhRRsuV/oGWgGciLQGvb+3gpRWtyWye1W7ejU2FaE/Lo0Ckv+ontiKSXZoOjZbGqQaQiIYApEv6ZdXHDVytzjU6G9tMYFwZwlTam21Db5NmaN0gukv6d50zAJ/VRubrF3Wf1xs7AIDVYB+LIv0bxl6DxBqgxzU2K0dDSPLZOqAGalm6Tp+ppiExkUgsZ6aNZ0Npf2sHq6uas62zehHB5SwUG5vqhU6t0Wx7CL2UM+OqEceDnLl8Ga3LdmoIZ0ZjE3va2JytpZydW9JcTRtngOasPgddr7F5eHjMG06zuwcR3QvgY9DC4mFm/mjp8QAiAiIelP6rormda27jhuVtAMDWul5p61xrIdoZlP4Qjc2s4LnSiiRki/t2h9XsyCVtmERV5ITdOE6crqYGAFHoSH0dMJk0jV0tRHdZH985o89xsJ5gdU1rN+fbWwCA9UhL/xW1b++LiaVtUDKwsqaIrMQ2GsGolTVl3QOGS2LXhcBEPEWO20feLroaSD+a23hVNJmtdW0X7WwuINqVNhq7aC+xBTArcQbA1POznJVpQFU5C4NSzszKbmdNOLt+kLMzYhOeOs6O4EJiwDhC/rsxQEQ/B+APAbwJwNul1kHRcYXjDBHdDuAxAGcBfA3AB5i5W3bNkQMbEQUA/hK6uswlAF8loseZ+b/LfqegiTcPwSJJ5IESQ3S0jZvb1wAAV87oge3VcxHCXd2koGvcBdi6dLjnBuRlMdMZJy0Om2ywbjWiYdHJNPiSkFIDRmftHqCnX7ygX5J4UW97K6Hzcogx+foD3LiiX47zTT2AmwGirTr2fhj/sey0sFqsoSr5vuhlKTJIm4HH3NPI5cws+Ei716Md3LSY5WzjfIRwTzjriOE/jmzd0XzptswAZzlT4LiAszLkIhrG5mw5xMGZLGeqgLOzoR7Ip42zeqaiJ5Zo8lkAPwvgr4cdMGKc+VMADzHzY0T0VwB+CcDHyy5YZYHl7QCeZ+YXZJR8DMB7q/TGw8NjusFxXOnvSNdgfo6ZvzXisMJxRmqJ3gXgohz3GejaoqWoMhW9GcB3nc+XAPxQ/iAiehDAg/Kx8z23vjKyqOkMYh3ARtWDX5DtE8fTljoxVr9mDGP17duynQHO7jjqCbZx9Ykn+eJ6xcNbROROIS8w84WjtsHBsHHmOgCbzNZwcQkFJfryqDKwFWm9A/ME6eQFACCip5n5zgrnnin4fs0e5rVvuUHmUGDme+toCwAQ0ZMAbij46iPMXFb53Z6iYB+X7C9FlYHtEnR1ZoNbALxU4XceHh6nBMx8zxFPMWyc2QCwRkShaG2Vxp8qNravAngDEd1ORA0A9wF4fOxme3h4eAxH4TjD2rHxKQDvk+MeADBSAxw5sMko+avQZofnAHyWmb854md1zr2nCb5fs4d57dvM9IuIfoaILgF4J4AvENETsv8mIvoiMHKc+X0Av0VEz0Pb3D4x8po8p7FiHh4epxe1xtN6eHh4TAP8wObh4TF3qHVgI6J7iehbRPQ8EX2oznNPCkR0KxE9RUTPEdE3ieiDk25T3SCigIi+TkSfn3Rb6gIRrRHRRSL6H+HunZNuU10got+UZ/FZInqUiFqTbtO0obaBzQmJ+AkAbwZwPxG9ua7zTxB9AL/NzG8C8A4AvzIn/XLxQWiD7TzhYwD+mZnfCOAHMCf9I6KbAfw6gDuZ+S3QUV33TbZV04c6Nba5DL1i5peZ+Wvy/zb0CzLS83lWQES3APhJAA9Pui11gYhWAPwIZPWMmbvMvDnZVtWKEMACEYUA2vB+pQOoc2ArComYmwEAAIjoNgBvA/CVybakVvwFgN8DRtR8my18L4DXAHxKptgPE9HipBtVB5j5/wD8GYDvAHgZwDVm/tJkWzV9qHNgO1Tow6yAiJYA/B2A32DmrUm3pw4Q0XsAXGbm/5x0W2pGCOAHAXycmd8GYBfAvNh8z0DPhG4HcBOARSL6+cm2avpQ58A2t6FXRBRBD2qPMPPnJt2eGvEuAD9FRC9Cmw7uIqK/mWyTasElAJeY2WjWF6EHunnAPQC+zcyvMXMPwOcA/PCE2zR1qHNgm8vQK0mb8gkAzzHzn0+6PXWCmT/MzLcw823QfP0LM8+89GfmVwB8l4hMBoy7AZTmD5whfAfAO4ioLc/m3ZiThZE6UVtqcGbuE5EJiQgAfLJC6NUs4F0APgDgGSL6huz7A2b+4gTb5DEavwbgERGyLwD4xQm3pxYw81eI6CJ0Jtk+gK9jhsKrTgo+pMrDw2Pu4CMPPDw85g5+YPPw8Jg7+IHNw8Nj7uAHNg8Pj7mDH9g8PDzmDn5g8/DwmDv4gc3Dw2Pu8P+L6SmbKffNJwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# calculate density difference\n",
    "U = phi[:,:,0] - phi[:,:,1]\n",
    "\n",
    "# plot.\n",
    "fig = plt.figure()\n",
    "plt.imshow(U ,origin='left',extent=[0,V,0,V], vmax=1.0, vmin=-1.0)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# output arrays of the monomer densities.\n",
    "np.savetxt(\"phi0_fA=0.8_it=500.txt\", phi[:,:,0])\n",
    "np.savetxt(\"phi1_fA=0.8_it=500.txt\", phi[:,:,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}