{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "9cd6f846",
   "metadata": {},
   "source": [
    "# Capacitance Model"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "3603def1",
   "metadata": {},
   "source": [
    "After solving the Thomas-Fermi equations for n(x), QDFlow builds a capacitance\n",
    "model to determine the dot occupancies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b6c64264",
   "metadata": {},
   "outputs": [],
   "source": [
    "from qdflow.physics import simulation\n",
    "from qdflow import generate\n",
    "import tutorial_helper\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "98ef4861",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a set of default physical and numerical parameters\n",
    "phys = generate.default_physics(n_dots=2)\n",
    "phys.gates[3].peak = 7.5\n",
    "x = phys.x\n",
    "q = phys.q\n",
    "numerics = simulation.NumericsParameters()\n",
    "\n",
    "# Calculate V(x) and solve for n(x)\n",
    "# This may take a few seconds initially due to numba compilation time\n",
    "V = simulation.calc_V(phys.gates, x, 0, 0)\n",
    "K_mat = simulation.calc_K_mat(x, phys.K_0, phys.sigma)\n",
    "delta_x = x[1] - x[0]\n",
    "g0_dx_K_plus_1_inv = np.linalg.inv(phys.g_0 * delta_x * K_mat + np.identity(len(x)))\n",
    "n, phi, converged = simulation.ThomasFermi.calc_n(phys, numerics, V, K_mat, g0_dx_K_plus_1_inv)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "657422d1",
   "metadata": {},
   "source": [
    "The first step is identifying the regions in which charges are confined. This is\n",
    "done by identifying charge \"islands\" for which n(x) is above some relative value.\n",
    "\n",
    "This is the step that is responsible for determining whether the dots should be\n",
    "handled individually, or as a single, combined center dot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "048f9784",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAccAAAEmCAYAAAAX5FuVAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAQkBJREFUeJztnQd0VOW6hr9ACoRA6KFIiUjvUgLoEZWmxoIN5CAiciBwQUVYChwFbBhEVAT7VUC9dBULxSNFpRdp0g8iRWnSEnpIwr7r/Yd/mJnMhEwyZe+Z91lrMzN7/7vMnjDvfN//lQjDMAwhhBBCiJ1CV58SQgghhOJICCGEuIGWIyGEEOICxZEQQghxgeJICCGEuEBxJIQQQlygOBJCCCEuUBwJIYQQFyJdVxDfc/nyZTl06JAUL15cIiIieIsJMRmohXLmzBmpVKmSFCpEm4FQHAMChLFKlSr8eyPE5Pz5559y3XXXBfsyiAmg5RgAYDHq/3glSpQIxCkJIV5w+vRp9QNW/18lhOIYALQrFcJIcSTEvHDag2joXCeEEEJcoDgSQgghLlAcCSGEEBc450gIIX5IDcnKypLs7GzeWxMRFRUlhQsXztNYiiMhhPiQS5cuyeHDh+X8+fO8ryYMuEKqTlxc3DXHUhwJIcSHBT/27t2rrBMUFIiOjmYErIms+WPHjslff/0lNWvWvKYFSXEkxGpcuCDyr3+J7N8v0rixyCuviJQuHeyrIlesRggkciZjY2N5T0xGuXLlZN++fZKZmUlxJCSkMAyRPn1Epk2zvV6xQuTwYZGvvoLPKNhXR67AEnTWz2NltCohVmLyZJGpU0XgEho9WiQyUmTOnKtiSQjxCRRHQqxkNb79tu05XKn//rfIyJG213jEdkKIT6A4EmIVVq0S2bpVpGhRkf79besGDxbB3NYff4isXRvsKyQhzuLFi6Vu3bp5TlHZvn27ig49d+6cWA2KIyFW4eOPbY9du4qULGl7XqyYSOfOtufTpwfv2khY8Nxzz8kLL7yQ51zBevXqSatWreStt94Sq0FxJMQKZGSIfPml7Xnfvs7bunWzPc6cKcKkc+Inli9fLnv27JEHH3zQq/169eolH3zwgSqKYCUojoRYxaUK11RCgkhSkvO2jh1FSpUSOXJEZM2aYF0h8QTmgvHZBWPxYh761ltvlaeeekpZh6VLl5YKFSrIiy++aN8+Y8YM6dChgxQpUuTK2zKkffv20qlTJ/UcnDx5UrlRR+q5cBG1D9b/8ssvlvoboTgSYgV+/NH22L498gSct0VHi9x2m+35smWBvzaSO6iUg4oswVi8rNLz2WefSbFixWTNmjUyduxYefnll2XhwoVq27Jly6R58+ZOaREYv27dOpkwYYJa169fP6lcubKTOKIQQpMmTdT+VoJFAAixkjjCSnTHzTeLfP01fF8iQ4cG9NJI6NCoUSMZNWqUeo4qMu+++64KwoH1t3//flX1xxEI4UcffSSPPfaYHDlyRObPny8bN26USKQYOYD9sL+VoDgSYnaOHxfZsMH2vEMHz+KoiwJcvpzTuiTBA9HEZ88G79xeiqMjFStWlL///ls9v3Dhgt2l6sjDDz8sc+bMkTFjxqi5RYiqK0WLFrVcrVmKIyFm5+efbXNHDRrg28r9mCZNbF+Ep06J7NghUr9+oK+SeAJVWRBVbJGuFY5ERESocnigbNmycgp/Xy5A9NavX68iWHfv3u32uJhzrFGjhlgJ/rwkxOzoIBttHboDX2qtWtmew7VKiI9p2rSpylt0ZciQIapc3oIFC9Tc45IlS3KM2bp1q9rfSlAcCTE769bZHlu2zH2cFk9EthLiYzp16qTSORyZN2+eTJo0SaZOnarmJZ999lnp2bOnk4WJQt8HDx5Uka1WguJIiJlB3uKvv+ZNHPUv899+8/91kbCje/fusm3bNtm1a5d6jfZPvXv3VukeN954o1r30ksvSUJCgopa1UyfPl06duwo1apVEysRYegEFeI3Tp8+LfHx8ZKeni4lSpTgnSZ5B+XiGja0zVmlp9sKjnsCJeQwr4PUDuS4uUQMEv//H7148aLq55iYmOg2eMXqPPvss+peIUI1ry28EKAzbdo0uemmmyTYePP50HIkxMzoeqnIL7tWya7q1W25bZcuifz3vwG5PBJePP/888oC1EE61+LAgQPy73//2xTC6C38aUlIKMw3AqRvIKJ19WqRLVtQ2NLvl0fCi5IlSyqxyys33HCDWqwILUdiLTALsGePzYUYDqxfb3t0qEySKzpPjfOOhISXOL733ntSvXp15S9OSkqStddo0zN79mypU6eOGt+wYUNVwUGTmZkpQ4cOVetRMglVHFDp4dChQzlydDAZjbkI/HLCJPTZYCX1hjMQRcy/4ZdorVpXC3GHKnBdbdtme+6SnO0RiiMh4SeOM2fOlMGDB6vyRhs2bJDGjRur8GJdwcGVlStXSrdu3ZSYoaRR586d1YKcG528iuOMGDFCPX799dcqEuvee+91G6WFGoNz586VpUuXSl/XzgjE/0LRu7dNLJBUjShOdKO4UvcxJNm711YbMybG9oMgL+DHA6DlSEjBMCxEy5YtjQEDBthfZ2dnG5UqVTJSU1Pdju/SpYuRnJzstC4pKclISUnxeI61a9cietfYv3+/er19+3b1et26dfYxCxYsMCIiIoyDBw/m6brT09PVMfBI8sknn8ChahixsYaxe7dh/POfttdNmxrG5cuheVvnzLn6HvPKyZO2fbCcPu3PqwspfPV/9MKFC+o7A4/EfHjz+VjGckRIMEoUOSaSoioDXq/ykPSM9a6Jp7A0PY0HCOVGySS4T/Ux8NyxGj2OiXOjcr07MjIyVLiz40IKAL7qX3/d9vyll2xWFLoAFC0qsnGjrbxaKHLFw6GCbPIKWleVK2d77qGUFyHk2lhGHI8fPy7Z2dkqwdQRvEY1eHdgvTfjkQODOUi4YnWuE8aWL1/eaRwqzqPfmafjpKamqpwpvVSpUsWr90pcQNFtfNFDDHVycZky6KJqez5uXGjeMkScOrpK84ou/ExxJCT0xdHfIDinS5cuqmknKssXhOHDhysLVC9//vmnz64zLJkxw/Z49922PD7NoEG2xwULRDzMO1saiiMJE7755huV8oHi5YOu/L92ty6QWEYcUREeN+no0aNO6/EaHavdgfV5Ga+FEf3GEHTjWCEDY10DfrKyslQEq6fzxsTEqGM4LqQAgThaHBGA42ohoWQa3K4QyFAiI+NqIr83blVAy5EEgSlTptino7wlJSVFHnroIWVIvPLKKx7XBRLLiCO6STdr1kw13tSgSgNet27d2u0+WO84HkD8HMdrYUSrlUWLFkkZuOtcjpGWlqbmOzWoOo9zI5WE+JlNm0T++stmMd55Z87tsCbB99+H1kcBYUREbnw8Osp6ty/FkViIs2fPKgME8SBIpytevLjbdQHHsBAzZswwYmJijClTpqiIo759+xolS5Y0jhw5orb36NHDGDZsmH38ihUrjMjISGPcuHHGjh07jFGjRhlRUVHGli1b1PZLly4Z9957r3HdddcZmzZtMg4fPmxfMjIy7Me54447jKZNmxpr1qwxli9fbtSsWdPo1q1bnq+b0aoFYPx4W+TlnXe63752rW17XJxhOHxmlmf2bNv7Skryft8NG2z7linjjysLSRitejUD4PXXXzdq1KhhREdHG1WqVDFeffVV46efflLRvKdOnbLfs40bN6p1e/futW93XPB9C06ePKm+m/FdXbRoUfV9+t///ldtc7efp3WBjla1lDiCiRMnGlWrVlUfHFI7Vq9ebd/Wtm1bo2fPnk7jZ82aZdSqVUuNr1+/vjFv3jz7Nnyorh+Cuw/jxIkTSgzj4uKMEiVKGL169TLOnDmT52umOBaABx+0fdG/9pr77dnZhlGhgm3Mjz8aIcOrr9re02OPeb8v/jZ1OgdSO0hwxRGpRplng7N4meb03HPPGaVKlVIGyO+//24sW7bM+N///d9riiOMifHjx6vvR21g6O9IGCB169Y1li5dqoyQTp06GTfccIMyTrDfrl271HG++uoru2Hibl2gxdFytVUHDhyoFnf87Cak/+GHH1aLO1BpJy9NSRCZiqryJMDgs1m2LPdGv6gn2qmTyGefifzyi0iHDhISXGkLJLVre78vXNAVK4ocPmyLWM1LXVbiP7LPi8xyCCQLJF3OikQWy9PQM2fOyDvvvCPvvvuu6skIatSoITfffLPb71bXaS9E5iMNzjEWA9NV3333naxYsULatGmj1qH3IyL4EXCD72adDYDvWb2vu3WBxjJzjiQMwRc7gqHQgqlFC8/jdMX/lSslZCiIOALOOxIv2bFjh8rRbteunU+PidQ3x/gMxHXUrl1bbTMzlrMcSRihrUZYPrn1XtPiiKIMmZkiUVFieYu5oOKIQglLl4ZPgXYzUzjWZsEF69x5pCjyiD1QCB4aW4yKUzBjKEPLkZgXtF4C1+oFV6cOeunY6pCGQk1RWMtobIwasmhenB/Q2xHs2+fTSyP5AJ8jXJvBWHDuPIKmxBBI1wh/UO5K1aXDcNVfYRMiyV1cqyjU4kjdunVV6ptjNbETJ06oGtb1TN5SjeJIzIsWuhtvzH0cftXq9JxQcK1qq7FaNVtVoPyAfcH+/b67LhLSoHMRKoQ999xz8vnnn8uePXtk9erV8umnn6pkfMwTvvjii2oecd68efLmm2/miOFACgbEFRXN0NgBgnvfffdJnz59ZPny5bJ582Z59NFHpXLlymq9maE4EnOCX6C6tmhe2jVdmewPKXHMr0sV0HIk+QAdioYMGSIjR45UVl/Xrl1VvmFUVJRMnz5ddu7cKY0aNZLXX39dXn31Vad9EXDTr18/tQ8szbFjx6r1kydPVjnqd999t8obh2sWrQNxTDMTgZDVYF9EqIPC44jkQik5VsvxIgke4oC5RvTOLFw49/FwBaHIfGKi9efZnn3WVi/2ySdtBdbzAyxGCCS+gC5etFnXxO//R1Gfee/evZKYmKgsMWIuvPl8+D+GmNulitJp1xJG0KTJ1R6ImK+zelNnx4jT/ICqOrhvCJpwmCcihOQNiiMxJ5s32x4bN87beJT9091P9L5WBQIPYAXnl8hIkeuusz1nUA4hXkNxJOa2HPMqjgBFyIFLFJ2lwCyHdgtff33BjqXnHRmUQ4jXUByJOdHWX16CcVxdq1YWx5MnMQHmLG75hUE5hOQbiiMxH+fOXbV2vGnXpMVx40axLNpqRPm32LwncLuF6RyE5BuKIzEfv/9+dR7RpYVYnsRx2zaRS5ckbOcbNbQcCck3FEdiPnST31q1vBcD9D9EhKbJ6zZ6xFfzjYCWIyH5huJIQicJHqWytBt2+3aRcBdHHb2LZtFMZybEKyiOJHQsR6DrNVIcbbmOeg5XB/kQQvIExZGYj3AWR1/OOSKgp1Spq9YjIblw6623yqBBg3w2rqCgjmsTHUcQBNiyipi3XVNBxNGKc45ZWVejdH3hVgUoBHDqlE0c69f3zTFJSPL111+bvt5pIKHlSMzFiRMiaWlXexJ6S926VxslWy1i9c8/bQXX0dy5UiXfHFNXyaHlSK5B6dKlpXjx4rxPV6A4EnO6VKtWzV+7JohBXJzNCtMpIVYLxoFL1VeFwvW848GDvjkeCVkc3aXvv/++ajeF4twJCQny0EMPedzviy++kObNmythrVChgvzzn/9UnTw0P//8s0RERKhWVhgXGxurOnigp6MjY8aMUefCcXr37q2KhDuC47Rs2VKKFSsmJUuWlJtuukn2+7H6E8WRmAtYfAUpuo2IVavOO/oyUlVDyzHooPHRuUvngrLkp+nSr7/+Kk899ZS8/PLLSsB++OEHueWWWzyOz8zMlFdeeUX1avzmm29k37598vjjj+cY9/zzz6sekDh+ZGSkPPHEE/Zts2bNUnOMr732mtpesWJFJdAaNEzu3LmztG3bVn777TdZtWqV9O3bV4muv+CcIzFnQEqNGvk/BsRx7VrriaMvg3E0FMegcz7zvMSlxgXl3GeHn5Vi0cW82ufAgQPKOkP/RVhx1apVk6a6brEbHEXu+uuvlwkTJkiLFi1U4+M4eHGuMHr0aCVuYNiwYZKcnKysQ1in48ePV9YiFoBekYsWLbJbj2gphnZiuKYaV74b0G/Sn9ByJOYUiILUFdX/aXbuFAl3y1G7VTnnSPJIhw4dlCBC6Hr06CFTp06V8+fPexy/fv16ueeee6Rq1apKTLUAQmQdQZNkDSxDoN2vO3bskKSkJKfxaIzsOB8Ka7RTp07qXO+8844c9nMrNlqOxFzo9koFEUftkrXqnKM/3KqccwwasVGxyoIL1rm9BQK3YcMGNcf3448/ysiRI5XLc926dWquz5Fz584pwcICES1XrpwSRby+5BIQ5xgJq92hly9fzvN1TZ48Wbl74eadOXOmvPDCC7Jw4UJp1aqV+AOKIwk91yLFMac4otsHfv0XtJg58RoIgbeuzWATGRkp7du3V8uoUaOUKC5ZskQeeOABp3E7d+6UEydOqGCaKlcqMmHO0FvgIl2zZo089thj9nWrV6/OMQ7uXSzDhw9XluW0adMojiQMwC9NbeEUxHLUlhfy+yAKpUuL6UEFG6Sx+HrOEbVmixWzVcnBvc1voBMJG+bOnSt//PGHCsIpVaqUzJ8/X1l4td2Uc4QrNTo6WiZOnCj9+vWTrVu3quAcb3n66aeV2xTRrIhChRW6bds25doFe/fulY8//ljuvfdeqVSpkgoU2r17t5OY+hrOORJz5fnBzVKkiEhCQv6PA+tIz7Xp6FerWMzoQlKihO+OC/cV5x2JF8BK/Prrr+X2229XFt2HH34o06dPl/puikjAjTplyhSZPXu21KtXT1mQ48aN8/p+d+3aVUaMGCHPPfecNGvWTKVo9O/f374d6R+wUh988EGpVauWilQdMGCApKSk+O+zNSzGu+++a1SrVs2IiYkxWrZsaaxZsybX8bNmzTJq166txjdo0MCYN2+e0/avvvrK6NChg1G6dGnEPBsbN27McYy2bduqbY5LSkpKnq85PT1d7YNHkguLFiHw3DDq1Cn4bWrb1nas//s/a9zyr7+2XW+LFr4/tr4XU6f6/tghgq/+j164cMHYvn27eiTmw5vPx1KWIyZhBw8erHzgmDBu3Lixmvh1TDh1ZOXKldKtWzcVHrxx40aVJ4MFpr/jhPLNN98sr7/+eq7n7tOnj4qO0svYsWN9/v7CHl+mMlht3lEnMxfEneyJK5GB4ufoPkJCCUuJ41tvvaVEqlevXsqEh7kPc3vSpEluxyPc94477pBnn31WuQfgC7/xxhvl3XfftY9BqDKisTDxnBs4D6o/6KWEL11fxHeRqhpdes4qblUd9o7KQL6G4khI6IojwoKRT+MoYoUKFVKvUS3BHVjvKnqwND2Nzw1MEJctW1YaNGigIqVyy/vJyMhQSauOCwlQjqOrOFrFcqQ4EmIqLJPKcfz4ccnOzla19xzBa0zUuuPIkSNux2O9N6BWIJJiESWF0kVDhw5V0VKYtHZHamqqvPTSS16dg/jYcrSaW5XiSIipsIw4BhNERmkaNmyoqju0a9dO9uzZYy9l5AgsS8yNamA56hwgco1oVVCtWsFvk/5ckB6BlA7d1zAcxVF3+OCcIyGh51aFS7Nw4cJy9OhRp/V4jTlAd2C9N+Pzii5z9LsHqyQmJkbNSTou5BqgVdOhQ86J6wUBuX16rs3s1mNGBv4wbc8550iIKbCMOCLRFPkvaHuiQWIqXjvW4HME6x3HA5Qb8jQ+r2zatMmpPiDxAXB1QyALF8avGt/cUqvMO+q6p8jvRJ6jr9F/p+npIhcu+P74hIQglnKrwlXZs2dPVUUBfb1QyR2pGIheBaiWULlyZTXnp6suoAgu2qSgAvyMGTNUaSNUWtCcPHlS1QI8dMVq0T3GdFQqXKcoUXTXXXdJmTJl1JzjM888o6pHOBbSJT5yqcIFCIH0BZh3XLbM/OLo6FL1RwseVMmB8KLDAVyrvqzdSkiIYhnLUVdRQPUFpF40adJEWXAoQquDbiByjpXa0VATwgYxRE7kl19+qfqNIeJU891336lafRBP8Mgjj6jXSBPRFitap3Ts2FHq1KkjQ4YMUVUavv/++4C//7AQR1/OzVrFcvTnfCOA4DKdg5DQtRzBwIED1eIOVJF35eGHH1aLJ1DPz11jTg0CaX755Zd8Xi3x2rXoi/lGq+U6+lscAcQRqTIMyiEeuPXWW5XRAY9coAqyz5kzRxVmMSOWE0cSovjDcrRKOod+7/4WR0BxJCbh8OHDqrC5WbGUW5WEMP6wHHU6x7FjtmAUs1uO/kz3oTiSAGMYhmRlZeVYr/s8IqYDkf35xbVfpK+hOJLQtRyLF7/a3cPM1mOg3KqAliPJBYgZpq3i4+NV+hw6ZUDkwBdffKGCIdEMGcKG4iiOda0xrQVX6YIFC1RmAYRv+fLlyl2LYw4aNEgdE1XKAMYiBkTz559/SpcuXVRXkNKlS8t9990n+3RhkCtTYHDBjh49WhVk0S203n//falZs6YUKVJExZ889NBDPvmM6VYl5rIcfW09Yd4ROYQQx2bNxHTgiycQ4qjTY1zyfklgPuJcqk36FXRv8yYA+rPPPlONGtauXasi+1EABT0bUdM6MzNT1aeGKEEUkT0AwUK/R0eGDRumAifRi1G7TXFctKBasWKF2/Pi2BBNpNktW7ZMNVt+9dVXVW1sZAggMBIgNQ9540jJA7jGp556Sgk3AjCRfYD9fUJA+oSEOWxZdQ0yMw2jUCFbW6WDB31783v2tB139GjDlJw4Ybs+LOfP++88339vO0ezZv47h4XxZ8uqs2evfsSBXnDuvILWfHXr1jUuX75sXzd06FC1zh3r1q1T9+zMmTPq9U8//aRef/PNNzmO27Rp0xz7Y+ycOXPU8y+++EK1FnQ8d0ZGhlG0aFHjP//5j3rds2dPIyEhQa13bDlYokQJ4/Tp0+HdsoqEcAEANDmOjCxYk2N36PZXDu4ZU7qTy5UTKVrUf+cpX9726KG9GyGgVatWyt2pgSW3e/duVdcajR/uueceZUnCtYoccp1C5whcr67AzZobmzdvVhXHcNy4uDi1wLV68eJFlWvuWL5TW5GgQ4cOqu41rFR0WEKDiNyaQngD3arEPC5VzIv5qgCARhcx1x0/zEYgXKpA/+iAOOJHuz+KDRCPrs2zZ4N3bl9w8eJF5fbEAgEqV66cEkW8dg2MKYbSjS64W+fI2bNnlYDi2K7gXJ6OAzFFb1/Md/74448qB/7FF1+UdevWqbnLgkBxJMFH11StXNn3xza75RiISFWgv2BQx/XMGRHW+w0Y+B1yDW0wDWvWrHF6vXr1ahXsgs5HJ06ckDFjxtibKGC+z1egzy6a2ZcvX97rWtSYn0RrQiyjRo1SorhkyRJ54IEHCnRNdKsS84ij7h7hD8tx/35b7dZwtRxhQsTF2Z4zKId4/HM8oAJtdu3aJdOnT5eJEyeqMpxwpcKdidd//PGHqiyG4Bxf0b17dxXJighVBNTs3btXWYMItvlLe5bcMHfuXJkwYYKqlrZ//375/PPPVc1tHclaECiOJPjo9AJ/iCOsUcxlZmZeFeFwFEfAeUdyDVCf+sKFC6p29YABA5QwImIVrs0pU6bI7NmzpV69esqCRESqr4iNjZWlS5cqEYbFV7duXRU1C3dubpYkrET01b399tvVPij7CVGvX79+ga+JblUSfLRo+aPLCeYwITx//GFzrZqtr2YgquM4zjviPjAoh1yj/OYHH3yQY3u3bt3U4ojOgQTIZ3R87e64nvYFyJ1EyocnIM6u3HzzzR6PH3BxhLkLsxcmLKKC8IsChboR1YQkTEJM5VbV844QBQTl/OMfYipoORJiSvIsjogieuedd9QkLKoQoEJB0aJFVdIlQm0hjPAbDx06VIXWEmIqcTRjUA5Kax08aHseCItWu1U550iIb8QRliEmY1EN4auvvrJHK2kyMjJk1apVql8iclxQzie3ThiEBGzO0czpHPhRgPzOqCjfNXjODc45EuJbccTkq66H5w7U0IO/GQvq3jnWwyMkV5BacOJEYCxHs4mjdqmi2HqhAMTGURwJ8a045iaMrpQpU0YthHhlNaI6v7/a12jL0Ww/2gIZjONaCIAQ4p9oVRSexYKcEkcaNWqU30OScI9U9VfVFm05QoyQ0gE3ZrgF4wDOOQYMd1GbxFqfi9fiiPp6PXv2lB07dthPhFp8eI5H1OAjxDTBONpigmUKFy4SirVYhkt1HA3dqn4n6soPL0TyI2CRmAtd6q5wHspUei2OTzzxhNSqVUs+/fRTFbXqWKSWENMF4wDM58G1umuXbd7RbOIYaMvx5ElzWdAhBL50kZiu+xwiuZ3fkeYAXs5jx46pzwQl53wujigdhIjVG9AnjxArWI4AgqjF0SwEWhwRC4AfCpgKOX7cP0UXiEpmB46NgIk5KFSokKrCk5cfLF6LY7t27VR7EYoj8anl6O8vajMG5QRaHCGMKECOPEcsFEe/gC/eihUrqiLaaOJLzANSEiGQecFrcfzkk0/UnOPWrVulQYMGdh+75t577/X2kCSc0Qnpvu7jaPZ0DnTGSEuzPQ9kSTu4VnHPadUExMWal7ktYk68Fkck+69YsUIWLFiQYxsDcohpxdFslqNO40DPuUC2j2JQDiF5wuvM4yeffFIeffRROXz4sJrgdFwYqUq8Jlwtx0BHqmoojoT4RxzR8PKZZ55RkaqEFAgEhmj3XqAsRwQAXbwoYTffqNH3mfVVCfGtOKLX1k8//STB4r333pPq1aurQudJSUmydu3aXMej/1idOnXU+IYNG8r8+fOdtqMXWMeOHVVVH7iF0TTTFfQUQ28zjImLi5MHH3xQjvLLpeCcOmUrvu1o0fiLsmWvtmPXwhSO4kjLkRD/zDkix3H48OGyfPlyJTauATno3OwvZs6cqbpUo6ElhHH8+PGqtB26ViMyzJWVK1eq/mOpqaly9913y7Rp06Rz586yYcMGFUwEzp07p3qCdenSRfr06eP2vLCU582bp4Q2Pj5eBg4cqH4kYO6VFAD9AwNl46Kj/XsrEboN1+rWrTbXaq1aElal4zQUR0LyhuEl1atX97gkJiYa/qRly5bGgAED7K+zs7ONSpUqGampqW7Hd+nSxUhOTnZal5SUZKSkpOQYu3fvXpT7MTZu3Oi0Pi0tzYiKijJmz55tX7djxw41dtWqVXm67vT0dDUej8SBJUtQY8kw6tQJzG25+27b+T76KPgfw6232q5l6tTAnve772znbd48sOc1Ofw/SlzJV7PjYJX9Qek6WK0a5Ku0b99eRdC6A+thaToCS/Obb77J83lxTuQq4TwauGmRSIrjt2rVKsc+aOGFRXP69Ok8ny+sCFQwjkb3GTVDxGqwAnI450hInghAnxzfcPz4cRUN6xoIhNdHjhxxuw/WezPe0zGQOIqSUHk9Dty4cL/qxbX/JQmSOJolnQOBSGZwq7I4NiEe8dpyhEBNmTJFFi9e7LYrx5IlSyTcgXXraLHCcqRA5iKO/g7GMZs44n2jcgoqdfi7bJ4rqJAD4NlAIYJA5lgSEsri+PTTTytxTE5OVkEtgSqqW7ZsWVVtwjVKFK91LUNXsN6b8Z6OAZduWlqak/WY23HQ/BkLManluH+/BBVtNUIYA138GxG7WM6ds1mPFEdCfCOOM2bMkFmzZsldd90lgQSuzWbNmimLFRGnAFYrXiN61B2tW7dW2wcNGmRft3DhQrU+r+CciMjFcZDCARAde+DAAa+OQ0wkjsh1hOUUrB8wwUrj0MBSR+wAxJENBAjxjThCpIJVdByuStR1bd68ubRs2VKlciAVo1evXmr7Y489JpUrV1ZzftrKbdu2rbz55pvK0oWw//rrr/Lxxx/bj3ny5EkldIeudIeA8AFYhVgwZ9i7d2917tKlS0uJEiVUlSAIo7tgHGJicURXithYNNuzCVTNmhKW4oj7DXFkri4hvgvIGTJkiLzzzjtB6XTdtWtXGTdunIwcOVKaNGmiEvZ/+OEHe9ANRA5l7TRt2rRRuY0Qw8aNG8uXX36pIlV1jiP47rvvpGnTpko8wSOPPKJeI5dS8/bbb6s8SViOt9xyixJNFA8gFhNHTAGYYd4xWJGqGuY6EnJNIpDPIV5w//33qwo5sKLq16+fowgARSMnCMiBBZqenq4sT6ISbEXQKR3uTVgxWrT8DX4EoUoSvAceij74nQceEJkzR2TiRBEPUwJ+Be/7k09EXn5ZZMSIwJ/fhPD/KCmwWxVBKRBIQgoEcj91Lmgg6/SaISgnWGkcjqX0ABoeE0IKJo7nz5+X2NhYmTx5cl53IeTaLtXixW0WZKAwk1uV4kiI9ecckUqBeTfM37HoNrHcfKNZxBEdQXQnkmDNOdJyJMR34rhz505Veg1pHNWqVVOFv0ePHi1btmzJ6yEICb44BruEnHapItewdOngXAPFkRDfiSNqiSKFYdGiRcpyRO4ghPEf//iHXH/99eo1quOw4TGxhOWocx0DjZ7rhEs1QAU0ckBxJMQ/tVUReYlWUMgbPHbsmHz00UdKFJFvWK5cOZk6dWp+DkvCiWCJI8qnYY4T0bLaigun+UZAcSTE/4XHkcrRoUMHmThxouzfv19VkkHPR0JMKY6OuY7BiFjV4qjdu8EURxRDwEIIKXgqR27RrKg+gyR5QkwrjgDiuGNHcOYdHd2qwQK5tshPRvHzEydsVYMIIf5pWbV792657bbbfHU4EuoEUxyDGZRjBrcqrGe6VgkJjX6OJMQItuUYbHEMplsVUBwJ8Y1bFeXicoNRqsQrwlEcg9nk2BWKIyG+EceMjAzp37+/NGzY0O12BOO89NJLeT0cCWfOnr0aCBJO4ojkf6SPoMlx5coSVCiOhPhGHNEFA93s0TLKHZs3b6Y4Eu+sRgSCxMUF/q455jpeuoQ+bIF1qQajybErFEdCfDPniJZOaWlpubpd0U+REFO7VHXLpiJFbG7Ov/4Kr2AcDcWREN9Yjv/+979z3Q6rkkXJSZ7QtUV1X8FgRGsiIAaNreFavf768Enj0FAcCckVRquS8LMcgzXvSMuRkNASx9WrV3tVDGDbtm0FuSYS6hw7FlzLMdjiGOw0DkDLkZCCi2OPHj1UR47Zs2fLuXPn3I7Zvn27cr3WqFFD1q9fn5fDknAXR9Q5DRbBKCFHy5GQ0JpzhPB98MEH8sILL8g///lPVTu1UqVKUqRIETl16pRqZ3X27Fm5//775ccff/SY7kGI6cQxkJajWeccUYQ9WB1CCLGyOKK4+FNPPaUW1E9dvny5ymu8cOGCNG7cWJ555hlVOu5ahQIIcRJH/QUdDAJdQg4eF9QxdTx3MNH3HqksyDstXjzYV0SItQuPN2/eXC2E5BtYK2axHJHKgQLc/s471JVxUPQ7Pl6CDnJM0brrwgXb50FxJMQJRquS8HSrIlI2JiZwuY5mcqlqGJRDiEcojiSwYH7LDG5VlHALpGvVTJGqGoojIR6hOJLAcuaMbZ4r2JZjoCNWzRSpqqE4EuIRiiMJznwj5ryC3WRXi+Mff/j/XBRHQiwFxZEEFjO4VDW6bNzevf4/l7ZO6VYlJHTFcfHixSrh/1//+pc88cQTTou/ee+996R69eoqxzIpKUnWrl2b63gULqhTp44aj/zL+fPnO203DENGjhwpFStWlKJFi0r79u1l9+7dTmNwvoiICKdlzJgxfnl/IY8ZgnFcxXHPHv+fiwE5hIS2OKJnY8eOHZVAHj9+XBUBcFz8ycyZM2Xw4MEyatQo2bBhg8qxROWev3UhaxdWrlwp3bp1k969e8vGjRulc+fOatm6dat9zNixY2XChAny4Ycfypo1a6RYsWLqmBcvXnQ61ssvvyyHDx+2L08++aRf32vIYiZxrFEjMG5VpIpot2qgipznBc45EuIZw0sqVKhgfP7550YwaNmypTFgwAD76+zsbKNSpUpGamqq2/FdunQxkpOTndYlJSUZKSkp6vnly5fV+3njjTfs29PS0oyYmBhj+vTp9nXVqlUz3n777Xxfd3p6uoFbjcewZ+xYxKsaRo8ewb8Vp07ZrgXLmTP+O8+ePbZzFCmCPzrDNMycabuuW24xwh3+HyWueG05Xrp0Sdq0aSOBBudFzVa4PTWFChVSr1etWuV2H6x3HA9gFerxe/fulSNHjjiNiY+PV+5a12PCjVqmTBlp2rSpvPHGG5KVleXxWjMyMuT06dNOCzHhnGPJkmhE6v95R22ZJiaaq0wbLUdCPOK1OGKecdq0aRJo4MLNzs6WBJc2R3gNgXMH1uc2Xj9e65gomzdjxgz56aefJCUlRV577TV57rnnPF5ramqqElm9oNclMaFbNVDzjlp4IY5mguJIiO/Kx2Eu7uOPP5ZFixZJo0aNVN1VR9566y0JNTDPqcF7jo6OViIJEYxBlRUXhg8f7rQPLEcKpInF8ddf/TvvqI9tpvlGR3FEzVdUCkJhBEJI/sTxt99+kyZNmqjnjoEtAFGc/qJs2bJSuHBhOaob5V4BrytUqOB2H6zPbbx+xDpEqzqO0e/RHXC7wq26b98+qV27do7tEEx3okkc8hzN4FYNVFCOWS3HMmVsj9nZIunpIqVKBfuKCLGuOMK1GAxgrTVr1kxFySLiFFy+fFm9HjhwoNt9WrdurbYPGjTIvm7hwoVqPUhMTFQCiTFaDGHlIWq1f//+Hq9l06ZNar6zfDCb9VoVM1qO/narmtVyxA84FBxH1SL8aKE4EpJ/cQwmcFX27NlTdQVp2bKljB8/XjVf7tWrl9r+2GOPSeXKlZW7Ezz99NPStm1befPNNyU5OVnNG6LlFtzC2tKFcL766qtSs2ZNJZYjRoxQvSq1ACMwB2KJllzFixdXr9Gi69FHH5VS/DKxvjiGs+WoLXgtjjVrBvtqCDEPhsWYOHGiUbVqVSM6Olqldqxevdq+rW3btkbPnj2dxs+aNcuoVauWGl+/fn1j3rx5TtuRzjFixAgjISFBpXC0a9fO2LVrl337+vXrVfpHfHy8UaRIEaNu3brGa6+9Zly8eDHP18ww8SvgnunUiZMnDVOwb5/teqKiDCMry/fHR/qOfs9mTOVp0cJ2bd9+a4Qz/D9KXInAP8EW6FAHrlpEraanp0sJ9PMLVw4eFLnuOpHChW3Fx80QAIL5NvQ1RKI+qtj4ujD45s0icNljfk/Pt5qJO+8U+eEHkUmTRK54YMIR/h8lrpjg24mEZY6jGYQRQKh1AXJ/zDtql6rZ5htdg3IQsUoIsWOSbygSFphtvjEQ846OBQDMCMWRELdQHEn4pnEEImLV7JajY64jIcQOxZEEDrNajlq4/Gk5mlUcteVoxvlQQoIIxZEEjnB0q5o5jQPQrUqIWyiOJDyLjgfCrYpAcLO7VSmOhLiF4kgCh3bdmc1y1MJ18qRIWprvjovi9egLishcsxafpzgS4haKIwkcZnWrxsWJ6FKAvnSt6mMhd9KlQL8pA3KY8kyIHYojCRxmFUfHecfffw+f+UZHyxFFEFBGjhCioDiSwGHWOUdQq5bt8b//DZ9IVRAbK1KkiO050zkIsWbh8ZAHbq3s8xKSoEzbuRMi6ORVuphI1jkxFXUSbde2Z5vvrm3/Ltsxr69svvfrSIVSIocPixz7S6RKCHaaKRyLLgPBvgpiMVhb1Ux1G/EFOisuEJdESPjQ5axIZLFch7C2KnGFblVCCCHEBbpVzeb+wa/cUGT5cpE77rDNv/32m5iOCxdEypcTQY+afXtFyhYwaAhpIbrDB1I6EBFrVnr0EJkzR+SNN0RyafJt6f9XhHgJxdFMYF7kGu4fy3LstEgG5hsrmvM9Fi8mklDN1rZqz0GRClc6deSXvVtt77dyZZGSCWJqcH241uNnzPnZEBIE6FYlgeHoUdtjgomFonZt2+OuXQU/lj6GjoI1MywEQEgOKI4kMFhBHH2ZzqGPQXEkxJJQHElgCDfL0UriqKsD6c+IEEJxJAHCCuJYp47tcft234mjFlwzo+u+HjgQ7CshxDTQciSBwQriWL/+1RJyKBieXy5fFtm92zqWY7Vqtsc//7RdOyGE4kgChBXEsUIFkVKlbAJRkHlHiMz587Zi49ULGPUaCCpVsnUOuXRJ5O+/g301hJgCWo4kMGXxrCCOSKWpV8/2fNu2/B9n69arLlWzduNwBNcIgQR0rRKioDgS/3P2rC3J3uzi6OhaLYg46n0bNBDLoAsWUBwJUVAcif/RVmOxYuauFAO05ViQoBxtOWqhtQIUR0KcoDgS/2MFl6ovLUctjla0HFEhiBBiPXF87733pHr16lKkSBFJSkqStWvX5jp+9uzZUqdOHTW+YcOGMn/+fKfthmHIyJEjpWLFilK0aFFp37697NaRhlc4efKkdO/eXXXUKFmypPTu3VvOwlVIQk8cteWY34hVtObascO64ki3KiHWE8eZM2fK4MGDZdSoUbJhwwZp3LixdOrUSf72EGG3cuVK6datmxKzjRs3SufOndWyVf+yF5GxY8fKhAkT5MMPP5Q1a9ZIsWLF1DEvOnwxQhi3bdsmCxculLlz58rSpUulb9++AXnPIYGVxLFiRVtwCiJWZ8/OX4Nj/O2ggXBiolgGnc5BcSTEhmEhWrZsaQwYMMD+Ojs726hUqZKRmprqdnyXLl2M5ORkp3VJSUlGSkqKen758mWjQoUKxhtvvGHfnpaWZsTExBjTp09Xr7dv344+Dca6devsYxYsWGBEREQYBw8ezNN1p6enq2PgMSwZNQrxqoZx5b6bHvw94Xrr18cfmXf7Tp1q2/fGGw1LsXmz7brLljXCkbD/P0pyYJmuHJcuXZL169fL8OHD7esKFSqk3KCrVq1yuw/Ww9J0BFbhN998o57v3btXjhw5oo6hQVNiuGux7yOPPKIe4Upt3ry5fQzG49ywNO+///4c583IyFCLYyPVPIEgkK5dJeRAl3mrWI6gXz+R116zzTtiDjKv6RiQF+1SbdZMLIV2qx4/LtKokYQUn34q0qJFsK+CWAzLiOPx48clOztbEly+YPF6586dbveB8Lkbj/V6u16X25jyuvbkFSIjI6V06dL2Ma6kpqbKSy+95PV7VO44B5dvyNGkiViCkiVFnnpKZPRoEQ9/W7mSnCzyyitiKeLjRW64wTbXumWLhBTnzgX7CogFsYw4WglYt44WKyzHKrp+ZW7gy2nRIglJSpe2jjiCF18UufNO74NyMGepg3qsBAogrF4tsmmThByNGwf7CogFsYw4li1bVgoXLixHXToH4HUFlP1yA9bnNl4/Yh2iVR3HNLnyRY4xrgE/WVlZKoLV03ljYmLU4jUlSoi0a+f9fsT3REaK3HRTeN1Z9HXk3x8h1opWjY6OlmbNmsnixYvt6y5fvqxet27d2u0+WO84HiDiVI9PTExUAuc4BlYe5hL1GDympaWp+U7NkiVL1LkxN0kIISQEMSzEjBkzVCTplClTVBRp3759jZIlSxpHjhxR23v06GEMGzbMPn7FihVGZGSkMW7cOGPHjh3GqFGjjKioKGPLli32MWPGjFHH+Pbbb43ffvvNuO+++4zExETjwoUL9jF33HGH0bRpU2PNmjXG8uXLjZo1axrdunXL83UzEo4Qc8P/o8QVS4kjmDhxolG1alUjOjpapXasXr3avq1t27ZGz549ncbPmjXLqFWrlhpfv359Y968eU7bkc4xYsQIIyEhQQlvu3btjF27djmNOXHihBLDuLg4o0SJEkavXr2MM2fO5Pma+R+PEHPD/6PElQj8E2zrNdSBqxYpIunp6arKDiHEXPD/KLHsnCMhhBASKCiOhBBCiAsUR0IIIcQFiiMhhBDiAsWREEIIcYHiSAghhLhAcSSEEEJcoDgSQgghLlAcCSGEEBcojoQQQogLFEdCCCHEBYojIYQQ4gLFkRBCCHEh0nUFCT7p6VmSmRmazVKioiJEip2TTCNTzEZURJTER8b77P6b+b36832bCXwG8fH8miPew78ak4EvqJkzj0lIc9sSkbhzYkbuin5Y5n91MSzeq1/ft4no2rUcBZJ4Dd2qJiPUfrm7Jcu8v8kuXMoKm/fq1/dtIsLi/xTxORRHQgghxAWKIyGEEOICxZEQQghxgeJICCGEuEBxJIQQQlygOBJCCCEuUBwJIYQQFyiOhBBCiAsUR0IIIcSq4njy5Enp3r27lChRQkqWLCm9e/eWs2fP5rrPxYsXZcCAAVKmTBmJi4uTBx98UI4ePeo05sCBA5KcnCyxsbFSvnx5efbZZyUr62q1kJ9//lkiIiJyLEeOHPHbeyWEEBJcLCOOEMZt27bJwoULZe7cubJ06VLp27dvrvs888wz8v3338vs2bPll19+kUOHDskDDzxg356dna2E8dKlS7Jy5Ur57LPPZMqUKTJy5Mgcx9q1a5ccPnzYvkBICSGEhCaWKPy4Y8cO+eGHH2TdunXSvHlztW7ixIly1113ybhx46RSpUo59klPT5dPP/1Upk2bJrfffrtaN3nyZKlbt66sXr1aWrVqJT/++KNs375dFi1aJAkJCdKkSRN55ZVXZOjQofLiiy9KdHS0/XgQQ1ishBBCQh9LWI6rVq1SwqSFEbRv314KFSoka9ascbvP+vXrJTMzU43T1KlTR6pWraqOp4/bsGFDJYyaTp06yenTp5WV6giEs2LFitKhQwdZsWKFH94lIYQQs2AJyxHze65uzMjISCldurTHuT+sh+Xnau1BCPU+eHQURr1dbwMQxA8//FAJc0ZGhnzyySdy6623KlG+8cYb3Z4b47BoILaEEEKsQ1Atx2HDhrkNdnFcdu7cGcxLlNq1a0tKSoo0a9ZM2rRpI5MmTVKPb7/9tsd9UlNTJT4+3r5UqVIloNdMCCHEwpbjkCFD5PHHH891zPXXXy8VKlSQv//+22k9IkoRwYpt7sB6BNqkpaU5WY+IVtX74HHt2rVO++loVk/HBS1btpTly5d73D58+HAZPHiwk+VIgSSEEOsQVHEsV66cWq5F69atlchhHhEWHFiyZIlcvnxZkpKS3O6DcVFRUbJ48WKVwqEjTpG6gePp444ePVoJr3bbIhoW6SL16tXzeD2bNm1S7lZPxMTEqIUQQog1scScIyJM77jjDunTp4+a/0OgzcCBA+WRRx6xR6oePHhQ2rVrJ59//rmy7ODORC4kLDjMTULwnnzySSWIiFQFHTt2VCLYo0cPGTt2rJpnfOGFF1RupBa38ePHS2JiotSvX1/lTWLOEcKMSFdCCCGhiSXEEUydOlUJIgQQUaqwBidMmGDfDsGEZXj+/Hn7OswL6rEIkEEk6vvvv2/fXrhwYZUz2b9/fyWaxYoVk549e8rLL79sHwPXLNy/EF8UCmjUqJFK/bjtttsC+O4JIYQEkgjDMIyAnjEMwZwjLFnkXsKCzY309CyZOfOYhDS3LRGJOydm5K7oh2X+VxfD4r369X2biK5dy0l8fKTP/o+S8IDiGAC8/Y8HgczMDM3fLFFRESLFzkmmkSlmIyoiSuIj4312/838Xv35vs0EPoNrCSOgOBLLulXDibz8Z7Y28RI+99/c7zW8/u4ICbEKOYQQQkggoTgSQgghLlAcCSGEEBcojoQQQogLFEdCCCHEBYojIYQQ4gLFkRBCCHGBiU0BQBchYl9HQsyJ/r/JgmFEQ3EMAGfOnFGPbFtFiPn/r6KaFSEsHxcA0Frr0KFDUrx4cdXA2SroPpR//vkn603y3oT03wwsRggjuvygWQEhtBwDAP6zXXfddZb9a8OXnJW+6AIJ703o3BdajMQR/kQihBBCXKA4EkIIIS5QHIlHYmJiZNSoUeqR8N7kBf7NkFCBATmEEEKIC7QcCSGEEBcojoQQQogLFEdCCCHEBYojIYQQ4gLFkci+ffukd+/ekpiYKEWLFpUaNWqoKNVLly453Z3ffvtN/vGPf0iRIkVUFZSxY8fmuHuzZ8+WOnXqqDENGzaU+fPnh9wdfu+996R69erqPSYlJcnatWsllElNTZUWLVqoCk/ly5eXzp07y65du5zGXLx4UQYMGCBlypSRuLg4efDBB+Xo0aNOYw4cOCDJyckSGxurjvPss89KVlZWgN8NIXnEIGHPggULjMcff9z4z3/+Y+zZs8f49ttvjfLlyxtDhgyx35v09HQjISHB6N69u7F161Zj+vTpRtGiRY2PPvrIPmbFihVG4cKFjbFjxxrbt283XnjhBSMqKsrYsmVLyNzjGTNmGNHR0cakSZOMbdu2GX369DFKlixpHD161AhVOnXqZEyePFl97ps2bTLuuusuo2rVqsbZs2ftY/r162dUqVLFWLx4sfHrr78arVq1Mtq0aWPfnpWVZTRo0MBo3769sXHjRmP+/PlG2bJljeHDhwfpXRGSOxRH4hYIXGJiov31+++/b5QqVcrIyMiwrxs6dKhRu3Zt++suXboYycnJTsdJSkoyUlJSQuYut2zZ0hgwYID9dXZ2tlGpUiUjNTXVCBf+/vtvtJkxfvnlF/U6LS1N/QiaPXu2fcyOHTvUmFWrVqnXEMNChQoZR44csY/54IMPjBIlSjj9TRFiFuhWJW5JT0+X0qVL21+vWrVKbrnlFomOjrav69Spk3KvnTp1yj6mffv2TsfBGKwPBeBmXr9+vdN7RN1cvA6V95jXvw2g/z5wTzIzM53uC1zrVatWtd8XPMLNnpCQ4PS3gULl27ZtC/h7IORaUBxJDn7//XeZOHGipKSk2NcdOXLE6YsN6NfYltsYvd3qHD9+XLKzs0P6Pealw8ygQYPkpptukgYNGqh1eO/40VSyZEmP9yUvfz+EmAmKYwgzbNgw1SIrt2Xnzp1O+xw8eFDuuOMOefjhh6VPnz5Bu3ZiThB0s3XrVpkxY0awL4UQv8KWVSHMkCFD5PHHH891zPXXX29/jp6Tt912m7Rp00Y+/vhjp3EVKlTIEX2oX2NbbmP0dqtTtmxZKVy4cEi/x9wYOHCgzJ07V5YuXerUgg3vHS7ntLQ0J+vR8b7g0TWq1/XvhxAzQcsxhClXrpya+8lt0XOIsBhvvfVWadasmUyePDlHw9fWrVurL0XMLWkWLlwotWvXllKlStnHLF682Gk/jMH6UAD3CvfH8T3CzYjXofIe3YHAPQjjnDlzZMmSJSrlxxHck6ioKKf7grlopG7o+4LHLVu2yN9//+30t4Gej/Xq1QvguyEkjwQ7IogEn7/++su44YYbjHbt2qnnhw8fti8aRCQilaNHjx4qpB8pDbGxsTlSOSIjI41x48apaMVRo0aFZCpHTEyMMWXKFJWu0rdvX5XK4RiFGWr079/fiI+PN37++Wenv43z5887pXIgvWPJkiUqlaN169ZqcU3l6Nixo0oH+eGHH4xy5coxlYOYFoojUTls+J3kbnFk8+bNxs0336zEoXLlysaYMWNy3L1Zs2YZtWrVUrmA9evXN+bNmxdyd3jixIlKCPAekdqxevVqI5Tx9LeBvxvNhQsXjP/5n/9R6T740XT//fc7/bgC+/btM+68806VH4scR+TRZmZmBuEdEXJt2LKKEEIIcYFzjoQQQogLFEdCCCHEBYojIYQQ4gLFkRBCCHGB4kgIIYS4QHEkhBBCXKA4EkIIIS5QHAkhhBAXKI6EeODTTz+Vjh07+v3+/PDDD9KkSRNVp5UQYg4ojoS44eLFizJixAgZNWqU3+8PWoShcPfUqVP5WRBiEiiOhLjhyy+/VB0j0NQ3EKC12IQJE/hZEGISKI4kpDl27JjqF/jaa6/Z161cuVK1n3Jtr+UImvnec889OQSsc+fOMm7cOKlYsaKUKVNGNf91bONVvXp1efXVV+Wxxx6TuLg4qVatmnz33XfqOu677z61rlGjRvLrr786HRvnwro9e/b49P0TQvIHxZGEfE/LSZMmyYsvvqjE58yZM9KjRw/Vn7Bdu3Ye91u+fLk0b948x/qffvpJCRgeP/vsM5kyZYpaHHn77beVxblx40ZJTk5W54NYPvroo7JhwwapUaOGem1reGGjatWqkpCQIMuWLfPxHSCE5AeKIwl57rrrLunTp490795d+vXrJ8WKFZPU1FSP49HRPj09XSpVqpRjGxo7v/vuu6pR9N13363Ez9UCxflSUlKkZs2aMnLkSDl9+rS0aNFCHn74YalVq5YMHTpUduzYIUePHnXaD+fbv3+/D985ISS/UBxJWABXaFZWlsyePVsFvsTExHgce+HCBfVYpEiRHNvq168vhQsXtr+Ge9Wxuz2A21QDaxA0bNgwxzrX/YoWLSrnz5/Px7sjhPgaiiMJC+AKPXTokEqX2LdvX65jMZcYEREhp06dyrENUaWOYJxrCobjGGz3tM51v5MnTyo3MCEk+FAcSchz6dIlNd/XtWtXeeWVV+Rf//pXDqvNEQTr1KtXT7Zv3x7Q1BEIeNOmTQN2TkKIZyiOJOR5/vnn1RwiUiUw34d5vyeeeCLXfTp16qSCcgLF6tWrlau3devWATsnIcQzFEcS0vz8888yfvx4+eKLL1TeYqFChdRzRIV+8MEHHvfr3bu3zJ8/X4lqIJg+fboKGIqNjQ3I+QghuRNhOMaTE0LsILr0xhtvlOHDh/v1rhw/flxq166tUk0SExP5CRBiAmg5EuKBN954QyXt+xsECL3//vsURkJMBC1HQgghxAVajoQQQogLFEdCCCHEBYojIYQQ4gLFkRBCCHGB4kgIIYS4QHEkhBBCXKA4EkIIIS5QHAkhhBAXKI6EEEKIOPP/WKkkoTb6cwAAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 300x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# A relative cutoff value is defined in the NumericsParameters dataclass\n",
    "cutoff_value = numerics.island_relative_cutoff * np.max(n)\n",
    "\n",
    "# Calculate the regions correcponding to islands and barriers\n",
    "islands, barriers, all_islands, is_short_circuit = \\\n",
    "        simulation.ThomasFermi.calc_islands_and_barriers(phys, numerics, n)\n",
    "\n",
    "# Plot the results\n",
    "fig, ax = plt.subplots(figsize=(3,3))\n",
    "tutorial_helper.plot_n(fig, ax, x, n)\n",
    "tutorial_helper.plot_islands_and_barriers(fig, ax, x, n, islands, barriers, cutoff_value)\n",
    "ax.legend([\"n(x)\",\"cutoff\",\"islands\",\"barriers\"], bbox_to_anchor=(1.05, 1), loc='upper left');"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "7ba26138",
   "metadata": {},
   "source": [
    "Next, approximate charge occupancies $Z_i$ are calculated by integrating n(x) over\n",
    "each charge island. Note that this will give a real number, not an integer\n",
    "number of charges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e6ed409c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.53816277 1.80652947]\n"
     ]
    }
   ],
   "source": [
    "approx_charges = simulation.ThomasFermi.calc_approximate_charges(phys, n, islands)\n",
    "print(approx_charges)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "e142a09b",
   "metadata": {},
   "source": [
    "Now we build a capacitance model based off the result of n(x).\n",
    "\n",
    "Specifically, we are interested in the energy matrix $E_{ij}$, which gives the\n",
    "energy of an integer charge configuration $\\vec{Q}$:\n",
    "\n",
    "$E=\\sum_{i,j}E_{ij}(Q_i-Z_i)(Q_j-Z_j)$\n",
    "\n",
    "The energy matrix is found by calculating the average Coulomb potential energy\n",
    "between dots i and j. Then an additional term is added along the diagonal to\n",
    "account for the kinetic energy of charges in each dot. The strength of this\n",
    "second term can be adjusted by changing `PhysicsParameters.c_k`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "78e4d188",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Energy matrix:\n",
      " [[0.04923472 0.01370183]\n",
      " [0.01370183 0.04829884]] \n",
      "\n",
      "Energy of charge state [1 1] : 0.05757\n",
      "Energy of charge state [1 2] : 0.01321\n",
      "Energy of charge state [2 1] : 0.03171\n",
      "Energy of charge state [2 2] : 0.01476\n"
     ]
    }
   ],
   "source": [
    "# Calculate the energy matrix\n",
    "energy_matrix = simulation.ThomasFermi.calc_energy_matrix(phys, numerics,\n",
    "                K_mat, n, islands, approx_charges)\n",
    "print(\"Energy matrix:\\n\", energy_matrix, \"\\n\")\n",
    "\n",
    "# Calculate and print out the energies of several charge configurations\n",
    "configurations_to_check = np.array([[1,1], [1,2], [2,1], [2,2]], dtype=np.int_)\n",
    "for config in configurations_to_check:\n",
    "    energy = simulation.ThomasFermi.calc_cap_energy(config, energy_matrix, approx_charges)\n",
    "    print(\"Energy of charge state\", config, \": %0.5f\" % energy)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "d6f67632",
   "metadata": {},
   "source": [
    "QDFlow uses a brute-force search of all integer charge configurations close to\n",
    "`approx_charges` to determine which charge configuration has the lowest energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "6c337aa0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lowest-energy state: [1 2]\n"
     ]
    }
   ],
   "source": [
    "# Calculate the lowest-energy state\n",
    "island_charges = simulation.ThomasFermi.calc_island_charges(numerics,\n",
    "                 energy_matrix, approx_charges)\n",
    "print(\"Lowest-energy state:\", island_charges)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "cdff7ff3",
   "metadata": {},
   "source": [
    "This charge state vector has length equal to the number of charge islands.\n",
    "However, the number of charge islands varies depending on n(x). For example:\n",
    "* A double-dot has 2 islands, as seen above.\n",
    "* There could be 1 island if the voltage of one of the plunger gates is reduced.\n",
    "* There could be no islands if both plunger gates are reduced.\n",
    "* If the barrier is reduced, there is only 1 central, combined charge island.\n",
    "* There might be some other arrangment of islands for extreme / unconventional\n",
    "choices of parameters or layouts with more than 2 dots.\n",
    "\n",
    "When reporting data in a CSD, it is nice to standarize the length of the charge\n",
    "state vector, so that the charge states for all points can be stored together \n",
    "in a single numpy `ndarray`.\n",
    "\n",
    "Our solution is to set the length equal to the number of plunger gates (assuming\n",
    "one plunger gate per dot). Then the charge islands are matched up to the closest\n",
    "corresponding plunger. In the case of combined dots, all charges are placed in\n",
    "the left-most entry, and a seperate variable is used to record that the dots are\n",
    "combined together.\n",
    "\n",
    "For example, a single island with 3 dots will have `island_charges == [3]`.\n",
    "This can be expanded to length 2 in different ways, depending on where the\n",
    "island is located:\n",
    "* If the island is above the left plunger, it will become `[3, 0]`\n",
    "* If the island is above the right plunger, it will become `[0, 3]`\n",
    "* If the island is a single large, combined dot, it will become `[3, 0]`, and\n",
    "`are_dots_combined` will be set to True to distinguish from the first case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "8d09549e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Charge state: [1 2]\n",
      "are_dots_combined: [False]\n"
     ]
    }
   ],
   "source": [
    "# Calculate the charge state\n",
    "charge_centers = simulation.ThomasFermi.calc_charge_centers(phys, n, islands)\n",
    "are_dots_occupied, are_dots_combined, dot_charges = \\\n",
    "        simulation.ThomasFermi.calc_dot_states(phys, islands, island_charges, charge_centers)\n",
    "print(\"Charge state:\", dot_charges)\n",
    "print(\"are_dots_combined:\", are_dots_combined)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "69856d15",
   "metadata": {},
   "source": [
    "Finally, the sensor readout is calculated by calculating the screened potential\n",
    "due to each of the charges at one or more sensor locations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c23d7d33",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Single sensor output: [-2.39504244]\n",
      "Multiple sensor output: [-0.51501873 -0.90571996]\n"
     ]
    }
   ],
   "source": [
    "# Calculate sensor output\n",
    "sensor_output = simulation.ThomasFermi.sensor_from_charge_state(phys, n,\n",
    "        islands, dot_charges, are_dots_combined)\n",
    "print(\"Single sensor output:\", sensor_output)\n",
    "\n",
    "# Calculate output for multiple sensors at a time\n",
    "# Sensors should be provided as a numpy array with shape (num_sensors, 3)\n",
    "# giving the (x, y, z) coordinates of each sensor\n",
    "phys_multiple_sensors = phys.copy()\n",
    "phys_multiple_sensors.sensors = np.array([[-200, 100, 0], [200, 100, 0]])\n",
    "multi_sensor_output = simulation.ThomasFermi.sensor_from_charge_state(\n",
    "        phys_multiple_sensors, n, islands, dot_charges, are_dots_combined)\n",
    "print(\"Multiple sensor output:\", multi_sensor_output)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "qdflow_venv",
   "language": "python",
   "name": "qdflow_venv"
  },
  "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.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}