{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "df27beb5-e42d-4fa8-a47b-cbd4b7fd61a4",
   "metadata": {},
   "source": [
    "### Simulating 127-qubit kicked-Ising model\n",
    "\n",
    "In this tutorial we will use the ``pauli-prop`` package to classically simulate the time dynamics of a 127-qubit kicked Ising model on a heavy-hex lattice and approximate the single-site magnetization, $\\langle Z_{62} \\rangle$. The Hamiltonian considered is:\n",
    "\n",
    "$H = -J\\sum\\limits_{\\langle i,j \\rangle} Z_iZ_j + h\\sum\\limits_iX_i$\n",
    "\n",
    "where $J>0$ describes the coupling of nearest-neighbor spins, $i<j$, and $h$ is the global transverse field. A first-order Trotter decomposition of the time-evolved operator will be implemented as a quantum circuit, $U$, over $20$ Trotter steps. The coupling constant, $J$, will be fixed at $J=-\\frac{\\pi}{2}$ such that $U$ is Clifford any time $h\\mod\\frac{\\pi}{2}=0$.\n",
    "\n",
    "Workflow:\n",
    "\n",
    "- Create quantum circuits implementing the Trotterized Ising model\n",
    "    - Vary $h$ between Clifford points, $0.0$ and $\\frac{\\pi}{2}$.\n",
    "- Choose a site in the middle of the lattice (qubit $62$) and estimate the magnetization after $20$ Trotter steps for each circuit\n",
    "- Observe the magnetization roll off from $1.0$ to $0.0$ as $h$ moves from one Clifford point to another."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab545dd3-7771-4226-921a-e797ad4749c1",
   "metadata": {},
   "source": [
    "#### Create circuits\n",
    "\n",
    "First, we use ``FakeSherbrooke`` from ``qiskit-ibm-runtime`` to get the edge indices for the coupling map of the 127-qubit Eagle QPU. Once we have the edge indices, we create the Trotter circuits for the varying values of $h$ we want to simulate. The circuits will implement a first-order Trotterization of the Hamiltonian using 20 Trotter steps and 5420 gates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "92612eff-e570-45a0-be1c-16d22dfa9690",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from qiskit import QuantumCircuit\n",
    "from qiskit_ibm_runtime.fake_provider import FakeSherbrooke\n",
    "\n",
    "# Get the edges in the device coupling map\n",
    "backend = FakeSherbrooke()\n",
    "edges = backend.coupling_map.get_edges()\n",
    "\n",
    "# Num Trotter steps\n",
    "num_steps = 20\n",
    "\n",
    "# Create Trotter circuits with varying global fields\n",
    "hs = [i * np.pi / 32 for i in range(17)]\n",
    "circuits = []\n",
    "for h in hs:\n",
    "    circuit = QuantumCircuit(backend.num_qubits)\n",
    "    for _ in range(num_steps):\n",
    "        circuit.rx(h, [i for i in range(backend.num_qubits)])\n",
    "        for edge in edges:\n",
    "            circuit.rzz(-np.pi / 2, edge[0], edge[1])\n",
    "    circuits.append(circuit)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2dfebe6e-3956-40eb-b19c-97ea1c20516c",
   "metadata": {},
   "source": [
    "#### Simulate the time dynamics with Pauli propagation\n",
    "\n",
    "Next, we specify the observable, $O$, we want to measure. In this case, we will estimate $\\langle Z_{62} \\rangle$.\n",
    "\n",
    "To study the magnetization between Clifford points, we will approximately propagate $O$ to the beginning of the circuit, $U$, resulting in a new observable, $O^\\prime \\approx U^{\\dagger}OU$.\n",
    "\n",
    "During propagation, the number of terms in the evolved observable will grow exponentially in the number of gates, so it is usually necessary to limit how large it may grow. Here, we limit the number of terms in $O^\\prime$ to $10^5$ and truncate terms with coefficients smaller than $10^{-5}$. Truncating terms results in some error.\n",
    "Once we have $O^\\prime$, we can easily estimate the expectation value, $\\langle0|U^{\\dagger}OU|0\\rangle = \\langle0|O^\\prime|0\\rangle$, by summing the coefficients of the diagonal terms in $O^\\prime$, i.e. terms consisting only of ``I`` or ``Z`` on all qubits.\n",
    "\n",
    "We set ``frame='h'`` to specify we want to evolve the observable in the Heisenberg framework ($U^{\\dagger}OU$). Use ``frame=s`` for Schrödinger evolution ($UOU^{\\dagger}$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f0cfea45-2e5a-4e9e-9dd4-4b5502128e3a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finshed 17 simulations in 116s\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "\n",
    "from pauli_prop import propagate_through_circuit\n",
    "from qiskit.quantum_info import Pauli, SparsePauliOp\n",
    "\n",
    "# Z_62 observable\n",
    "id_pauli = Pauli(\"I\" * backend.num_qubits)\n",
    "observable = SparsePauliOp(id_pauli.dot(Pauli(\"Z\"), [62]))\n",
    "\n",
    "# Estimate the single-site magnetization for varying magnetic fields\n",
    "times = []\n",
    "approx_evs = []\n",
    "for circuit in circuits:\n",
    "    st = time.time()\n",
    "    propagated_obs = propagate_through_circuit(observable, circuit, 100_000, 1e-5, frame=\"h\")[0]\n",
    "    times.append(time.time() - st)\n",
    "    approx_evs.append(float(propagated_obs.coeffs[~propagated_obs.paulis.x.any(axis=1)].sum()))\n",
    "print(f\"Finshed {len(circuits)} simulations in {sum(times):.0f}s\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d214da4c-d5cf-43f2-97a7-380858d4a953",
   "metadata": {},
   "source": [
    "#### Visualize\n",
    "\n",
    "Finally, we plot the magnetization as a function of $h$. When $h=0$ the circuit does not affect the initial state in the computational basis, and we know $\\langle0|O|0\\rangle = 1.0$. As $h$ moves further from $0$, the magnetization decays monotonically until $h$ reaches the next Clifford angle, $\\frac{\\pi}{2}$, and the magnetization converges to $0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "01a3e294-f5c3-485c-b027-a14d9a6384c9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x123f68b00>]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAG0CAYAAAA7Go31AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAARFdJREFUeJzt3Qd4VFXex/H/TDqQBEJIJXSQEor0AHYQFVlAWbGirOjKgquw7AIqIjbs67uKsLAquIKgKIjCRhFFBQJIld6RGkICKSSkz/ucA4kJJJAyyb1z5/t5nmHm3rmTey5JZn451eZwOBwCAABgQXajCwAAAFBVCDoAAMCyCDoAAMCyCDoAAMCyCDoAAMCyCDoAAMCyCDoAAMCyPMXN5efny/Hjx8Xf319sNpvRxQEAAGWgpgFMS0uTiIgIsdtLr7dx+6CjQk5UVJTRxQAAABVw5MgRqV+/fqnPu33QUTU5Bf9RAQEBRhcHAACUQWpqqq6oKPgcL43bB52C5ioVcgg6AAC4lit1O6EzMgAAsCyCDgAAsCyCDgAAsCyCDgAAsCyCDgAAsCyCDgAAsCyCDgAAsCyCDgAAsCyCDgAAsCy3nxm5KuTlO2TdwdOSkJYpIf6+0rVxkHjYWTAUAAC3rtH56aefpH///nolUjWl86JFi674mhUrVkjHjh3Fx8dHmjVrJrNmzRIjxW47Ib1e/V7umblGnpi3Wd+rbbW/KoNV3P4k+XLzMX2vtgEAgMlqdNLT06V9+/bypz/9Se64444rHn/w4EHp16+fPPbYYzJnzhxZvny5DB8+XMLDw6Vv375S3VSYGfHxRrk4ZsSnZOr90+7vKLdEhzv9nJO/2iEnUjIL94UH+sqk/q2dfi4AAFyNzeFwmPLPf1Wjs3DhQhk4cGCpx4wbN06WLFki27ZtK9x39913S3JyssTGxpZ59dPAwEBJSUmp1KKeqhZF1dwUDRxFqYarsEBfWTnuRqc1Y5UWrAq+elUEKwAAzKCsn9+mqtEpr7i4OOndu3exfaom58knnyz1NVlZWfpW9D/KGVSfnNJCjqLCiHq++5TvpG5NH6nh7SE1fTylpren1PDxkFo+nlLD21Nq+Xjo+5o+RZ4vOFZvn3/s7WHXNTklpVTHhbCjnu/TOoz+QQAAt+XSQSc+Pl5CQ0OL7VPbKrycO3dO/Pz8LnnNlClTZPLkyU4vi+p4XBan0rL1raoVBCsVwGKa1q3y8wEAYEYuHXQqYsKECTJmzJjCbRWKoqKiKv111eiqsnh+QBtpHFxT0rPyJCM7V9KzciU9O08ysnLl7IV9Z7NyJSM778JzuZKRlafv0y/cl6excfwXv8q1zetJdGSAtIkIlBah/uLtaao+6AAAVBmXDjphYWFy8uTJYvvUtmqrK6k2R1Gjs9TN2dQQctUJWHU8dlymj8593RpWqilJdanKzMmXH/ckyGMfb7zi8b8lZch/k34r3FZNXleF+evgEx0ZKNERgXrb18ujXOVgCD0AwBW4dNCJiYmRpUuXFtu3bNkyvb+6qQ95NdJJdQ5WH/dFw07Bx796vrJhQHXS9vP20H1vrhSsgv195OnbWsmOE6my7ViKvqVm5srWYyn6JnKksOzNQ2pJWxV89C1AWoUH6L5CJWGkFwDAVZhq1NXZs2dl3759+vHVV18tb731ltxwww0SFBQkDRo00M1Ox44dk48++qhweHl0dLSMHDlSD0n//vvv5a9//aseiVXW4eXOGnVlRAgoGHUlpQSri0ddqW/10TPndMjRwef4+QB0Ov3SPkMqjzWtV0sHnzYRAToEtY4IkFX7EhnpBQAwXFk/v00VdNTkfyrYXOzBBx/UEwE+9NBDcujQIX1c0deMHj1aduzYIfXr15eJEyfq48rK2UGnupt1Khus1LdfvbZo8FG3hLTfR6YVpa6jtAkJq2IIPQAAlgk6RqiKoFPdqiJYJaRmyvbjqYW1P+rxseRzZXrtJ490Z6QXAKBKucU8OjhPhRpnB4uQAF99u6FlSOG+uWt/k6cW/j45Y2WH2gMAUNUYZ4wyaxxcy6lD7QEAqGoEHZR7CP2VGsV+3ntKsnLzqqlUAACUjqCDcg+hVy4Xdt5bsV/6v7NSfj2aXG1lAwCgJAQdlIsayaWGkKvRVUWpmp7p93eU9+7rKHVresuek2dl0Hur5dXYXZKZQ+0OAMAYjLqywKgrs430UvPyPLd4uyzeclxvNwupJa8PbidXN6hjcKkBAFbB8PIyIuhUndht8fLMom2SeDZLT0D4yDVNZHSfFuVebgIAgIp+ftN0hSpzS3SYLBt9rQy6OlLUHIP//umA3Pavn2XDb6eNLhoAwE0QdFCl6tT0ln8O6SAzh3aWEH8fOXAqXQZPj5MXv94h57LpuwMAqFoEHVSLPq1DZdno6+TOjvVFNZb+Z+VBXbvzyyFqdwAAVYegg2oTWMNL3ryrvXzwUGcJDfCRg4npcte/42TyV9slIzvX6OIBACyIoINqd2PLUPl29HVyV+fztTsfrjokt/7fz7LmQJLRRQMAWAxBB4YI9POS1wa3l1nDuug5eH5LypC7Z6yRSV9uk/QsancAAM5B0IGhrr8qRL4Zfa3c0zVKb8+O+01u+b+fZPX+xEvm7YnbnyRfbj6m79U2AABXwjw6zKNjGj/tOSUTvtgqx5LP6e37uzeQ8be2kpV7T8nkr3bIiZTfV0VXtUBqOQo1UzMAwP2kMmFg2RB0zCUtM0de+d8umbP2sN4OquEtpzOyLzmuYK0ttRwFYQcA3E8qEwbCFfn7eslLg9rK3OHdJLK2b4khRylI56qmh2YsAEBpCDowpR7NguXFgW0ve4yKN6o5S625BQBASQg6MK3UzJwyHacWFgUAoCQEHZiWWhXdmccBANwPQQem1bVxkB5dVdDxuCTqeXUcAAAlIejAtDzsNj2EXCkt7Dzcq7E+DgCAkhB0YGpq6LgaQh4WWLx5ysvjfLiZ/uN+2Zdw1qDSAQDMjnl0mEfHJagh5Gp0lep4rPrkXBXmLw+8v1a2H0/VC4R++ucYaVi3ptHFBABUEyYMLCOCjus6nZ4td8+Ikz0nz0pkbT/59LEYfQ8AsL5UJgyE1QXV9JaPh3eTxsE19bIR981cIwmpDDUHAPyOoAOXppqx5gzvJvXr+MmhpAy57z9rJelsltHFAgCYBEEHLi+itp/MHd5dwgJ8ZW/CWXng/XWSklG2yQYBANZG0IElNKhbQ+Y80k2Ca3nLjhOpMvTDdXqBUACAeyPowDKa1qul++zUruElW44ky8Oz1ktGdq7RxQIAGIigA0tpGRYg//1TN/H39ZR1h07Lox9tkMycPKOLBQAwCEEHltO2fqDMGtZVanh7yMp9ifKXORslOzff6GIBAAxA0IEldWpYR95/sIv4eNrl+10J8uT8TZKbR9gBAHdD0IFlxTStKzOGdhZvD7ss3Rovf1/wq+Tnu/X8mADgdgg6sLTrWtSTd++9Wi/8uXDTMXl60VZx88nAAcCtEHRgeTe3CZO3h3QQtcj5J+uOyOSvdhB2AMBNEHTgFvq3j5DXBrfXj2etPiSvxu4m7ACAGyDowG0M7lRfXhwYrR9P/3G/vPP9PqOLBACoYgQduJX7uzeUZ/q10o/fWrZHZvy03+giAQCqEEEHbmf4NU1k7M0t9OOXl+6S/8YdMrpIAIAqQtCBWxp1Y3MZeUNT/Xjil9vl01+OGF0kAEAVIOjAbY29+Sr5U8/G+vG4L36VLzcfM7pIAAAnI+jAbdlsNpl4eyu5r1sDUQOwxny6RWK3xRtdLACAExF0IO4edl4YEC13dIyUvHyHPP7JRvlhd4J+HLc/SdfyqHu1DQBwPZ5GFwAwmt1uk9fubCdZufmy5NcT8sjs9RLg5yWn07MLjwkP9JVJ/VvLLdHhhpYVAFA+1OgAKvF72PXsye0iAyU331Es5CjxKZky4uONErvthGFlBACUH0EHuMBus0lCWmaJzxU0XKnlI2jGAgDXQdABLlh38LTEp2aV+ryKNydSMvVxAADXQNABLiitNqeixwEAjEfQAS4I8fd16nEAAOMRdIALujYO0qOrbJc5xtNuk3r+PtVYKgBAZRB0gAs87DY9hFwpLeyoEVmD3lul59oBAJgfQQcoQs2TM+3+jhIWWLx5StX0vHpHW+nYoLakZebKn2b9Iu+t2CcONaUyAMC0bA43f6dOTU2VwMBASUlJkYCAAKOLA5NQQ8jV6CrV8Vj1yVHNWqrGJys3T55bvEM+WXdYH9evbbi8/sd2UsObuTcBwIyf3wQdgg4qYM7a3+S5xdslJ88hLcP8ZcYDnaVB3RpGFwsA3EZqGT+/aboCKuC+bg3lk0e6S3AtH9kVnyb9310pP+89ZXSxAAAXIegAFdS5UZB8/XgvaR9VW1LO5ciDH6yTGT/tp98OAJiI6YLO1KlTpVGjRuLr6yvdunWTdevWXfb4t99+W6666irx8/OTqKgoGT16tGRmMqEbqofqtDz/0e7yx071Ra0M8fLSXfLEvM1yLjvP6KIBAMwWdObPny9jxoyRSZMmycaNG6V9+/bSt29fSUgoeSjv3LlzZfz48fr4nTt3yvvvv6+/xlNPPVXtZYf78vXykNcGt5PnB7TR8+ws3nJc7py2Wo6eyTC6aADg9kzVGVnV4HTp0kXeffddvZ2fn69raR5//HEdaC42atQoHXCWL19euO9vf/ubrF27VlauXFniObKysvStaGcmdQ46I8MZ1hxIkpFzNkpSerYE1fSWd++9Wno0DTa6WABgOS7XGTk7O1s2bNggvXv3Ltxnt9v1dlxcXImv6dGjh35NQfPWgQMHZOnSpXLbbbeVep4pU6bo/5iCmwo5gLN0b1JXFj/eS9pGBsrp9Gx54P118sHKg/TbAQCDmCboJCYmSl5enoSGhhbbr7bj4+NLfM29994rzz//vPTq1Uu8vLykadOmcv3111+26WrChAk6/RXcjhw54vRrgXuLrO0nnz0WI4OujtTz8Tz/9Q7522dbJDOHfjsA4LZBpyJWrFghL7/8srz33nu6T88XX3whS5YskRdeeKHU1/j4+OgqrqI3oCr67bx1V3uZeHtrPdHgFxuPyV3/jpPjyeeMLhoAuBXTBJ3g4GDx8PCQkydPFtuvtsPCwkp8zcSJE+WBBx6Q4cOHS9u2bWXQoEE6+KjmKdW/BzCSzWaTh3s1lo/+1FXq1PCSX4+myB/eXSlrDyQZXTQAcBumCTre3t7SqVOnYh2LVVhR2zExMSW+JiMjQ/fjKUqFJYU+ETCLns2CZfGoXtIqPEASz2bLff9ZK/+NO8TPKAC4U9BR1NDymTNnyuzZs/VoqhEjRkh6eroMGzZMPz906FDdx6ZA//79Zdq0aTJv3jw5ePCgLFu2TNfyqP0FgQcwg6igGvLFiB7Sv32EXgF94pfbZfznW/XaWYrqyxO3P0m+3HxM36ttAEDlmWolwiFDhsipU6fk2Wef1R2QO3ToILGxsYUdlA8fPlysBueZZ57RzQPq/tixY1KvXj0dcl566SUDrwIomZ+3h/zr7g7SNjJAXvnfLpm//ojsPpkmQ7pEyb+W75UTKZnFVkuf1L+1Xk0dAGCReXSMwKKeMMJPe07J459s0ktHlMR24X7a/R0JOwBghXl0AHdybYt6svAvPfRMyiUp+Otj8lc7aMYCgEog6AAGOZmapfvrlEY9o5qz1h08Xa3lAgArIegABklIy3TqcQCASxF0AIOE+Ps69TgAwKUIOoBBujYO0qOrSu6lc75DsnpeHQcAqBiCDmAQtTSEGkKulBZ21PPqOABAxRB0AAOpoeNqCHlYYPHmqVo+ngwtBwCrTRgIuCMVZvq0DtOjq77+9bjMWXtYIgJ9CTkA4ATU6AAmoJqnYprWlb/3vUrPrbMn4awcSkw3ulgA4PIIOoCJ1K7hLd2b1NWPv9keb3RxAMDlEXQAk+nb5vzabgQdAKg8gg5gMqq/jrLxcLIkpDJZIABUBkEHMBk1AqtDVG39+NsdJ40uDgC4NIIOYEJ925yv1aH5CgAqh6ADmLifTtz+JEk5l2N0cQDAZRF0ABNqUq+WNA+ppVc3/2FXgtHFAQCXRdABTIrmKwCoPIIOYPKgs2L3KcnMyTO6OADgkgg6gElFRwZIZG0/OZeTJz/vTTS6OADgkgg6gEnZbDbp05rJAwGgMgg6gAs0Xy3feVJy8/KNLg4AuByCDmBiXRrVkTo1vORMRo6sO3Ta6OIAgMsh6AAm5ulhl96tzjdffbudWZIBoLwIOoCLNF99uz1eHA6H0cUBAJdC0AFMrlfzYKnh7SHHUzJl67EUo4sDAC6FoAOYnK+Xh1x/VT39mNFXAFA+BB3ApWZJpp8OAJQHQQdwATe0DBEvD5vsSzgr+0+dNbo4AOAyCDqACwjw9ZKYpsH6Mc1XAFB2BB3ARfRtUzBLMs1XAFBWBB3ARajlIGw2kS1HkiU+JdPo4gCASyDoAC4ixN9XOjaoox9/u4PmKwAoC4IO4JLNVwQdACgLgg7ggsPM1xw4LckZ2UYXBwBMj6ADuJCGdWtKyzB/yct3yPKdCUYXBwBMj6ADuJibCycPpPkKAK6EoAO4aD+dn/aeknPZeUYXBwBMjaADuJjW4QFSv46fZObky497ThldHAAwNYIO4GJsNlthp+Rvab4CgMsi6AAuqCDofLfzpOTk5RtdHAAwLYIO4II6NawjdWt6S2pmrqw9cNro4gCAaRF0ABfkYbfpJSEURl8BQOkIOoCLKuynsyNe8vMdRhcHAEyJoAO4qB7N6kotH085mZolW44mG10cADAlgg7gonw8PeT6q+rpx99sP2l0cQDAlAg6gAsrOszc4aD5CgAuRtABXJiq0fH2sMuBxHTZl3DW6OIAgOkQdAAX5u/rJT2b1dWPGX0FAJci6AAWab6inw4AXIqgA7i43q1DxW4T2XosRY6eyTC6OABgKgQdwMUF1/KRzg2D9ONvqdUBgGIIOoAF3NyGWZIBoCQEHcBC/XR+OXRaks5mGV0cADANgg5gAVFBNaR1eIColSCW70wwujgAYBoEHcByo69ovgKAAgQdwCL6Rp/vp/PzvkQ5m5VrdHEAwBRMF3SmTp0qjRo1El9fX+nWrZusW7fusscnJyfLyJEjJTw8XHx8fKRFixaydOnSaisvYBZXhfpLw7o1JDs3X37cfcro4gCAKZgq6MyfP1/GjBkjkyZNko0bN0r79u2lb9++kpBQcp+D7Oxs6dOnjxw6dEgWLFggu3fvlpkzZ0pkZGS1lx0wms1mo/kKAC5ic5hoJUBVg9OlSxd599139XZ+fr5ERUXJ448/LuPHj7/k+OnTp8vrr78uu3btEi8vrzKdIysrS98KpKam6nOkpKRIQECAE68GqH4bfjstd06LE38fT9kwsY94e5rqbxkAcBr1+R0YGHjFz2/TvAuq2pkNGzZI7969C/fZ7Xa9HRcXV+JrFi9eLDExMbrpKjQ0VKKjo+Xll1+WvLy8Us8zZcoU/R9TcFMhB7CKq6PqSD1/H0nLypXV+xONLg4AGM40QScxMVEHFBVYilLb8fElV8MfOHBAN1mp16l+ORMnTpQ333xTXnzxxVLPM2HCBJ3+Cm5Hjhxx+rUARrHbbdKndcHkgcySDACmCToVoZq2QkJCZMaMGdKpUycZMmSIPP3007pJqzSqw7Kq4ip6A6ykoJ/Osh0nJU9NrAMAbsw0QSc4OFg8PDzk5Mnif4Wq7bCw82/cF1MjrdQoK/W6Aq1atdI1QKopDHBHMU3qir+vpySezZJNh88YXRwAMJRpgo63t7eulVm+fHmxGhu1rfrhlKRnz56yb98+fVyBPXv26ACkvh7gjlQH5BtbhujHjL4C4O5ME3QUNbRcDQ+fPXu27Ny5U0aMGCHp6ekybNgw/fzQoUN1H5sC6vnTp0/LE088oQPOkiVLdGdk1TkZcGe/DzM/KSYaWAkA1c5TTET1sTl16pQ8++yzuvmpQ4cOEhsbW9hB+fDhw3okVgE1Yuqbb76R0aNHS7t27fT8OSr0jBs3zsCrAIx3XYt6umbn8OkM2RWfJq3C6YsGwD2Zah4dM4/DB1zN8Nm/yHc7E+TJ3s3lyd4tjC4OALj3PDoAnOvmIs1XAOCuCDqARfVuFSp2m8jOE6ly5HSG0cUBAEMQdACLCqrpLV0bB+nHjL4C4K4IOoCFscgnAHdH0AHcoJ/O+t/OyKm03xezBQB3QdABLCyytp+0jQwUNbbyu510Sgbgfgg6gMX1bVOwyCfNVwDcD0EHcJN+Oqv3JUlaZo7RxQGAakXQASyuWUgtaRJcU7Lz8uWH3aeMLg4AVCuCDmBxNputyOSBNF8BcC8EHcCN+ums2JUgmTl5RhcHAKoNQQdwA+3r15bQAB9Jz86T1fsTjS4OALhO0Dl9+rTk5+c7pzQAqoTdbpObW19ovtrGMHMA7qNCQWfHjh3yyiuvSI8ePaRevXoSEhIiQ4cOlc8//1zS09OdX0oATht9pebTyct3GF0cADBX0Nm9e7f87W9/k+bNm0v37t3ll19+kccee0xOnjwpS5culYYNG8rzzz8vwcHBcuutt8q0adOqtuQAyqVbkyAJ9POSpPRsWX/otNHFAQBzBZ3Vq1fr2pp//etfkpiYqGtvVC2OCjZdu3aVF154QbZs2SI7d+6UW265RT8PwDy8POxyU8sQ/fib7TRfAXAPNodDTQ7vvlJTUyUwMFBSUlIkICDA6OIAVSp2W7w89vEGvTTEynE36KHnAGDlz+9Kd0Y+d+6cHDt27JL927dvr+yXBuBk17WoJ75edjmWfE62H081ujgAUOUqFXQWLFig++z069dP2rVrJ2vXri187oEHHnBG+QA4kZ+3h1zbvJ5+/C2TBwJwA5UKOi+++KJs2LBBNm/eLB9++KE8/PDDMnfuXP2cm7eIAaYffUU/HQDuwLMyL87JyZHQ0PMzrnbq1El++uknGTRokOzbt4+2f8CkbmoVInabyO6TafL+zwekdUSgdG0cJB5qJwBYTKVqdNT8Ob/++mvhdlBQkCxbtkyPvCq6H4B5rDmQJJ7287/6LyzZKffMXCO9Xv1eYredMLpoAGCuUVdHjx4VT09PCQs7XxVe1KpVq6Rnz55idoy6gjtRYWbExxvl4l/6grqcafd3lFuiww0oGQCYcNRV/fr1Sww5iiuEHMCdqNmQJ3+145KQoxTsU88zazIAK7FXdmh50XWuli9fLm+++aaeKRmAuaw7eFpOpGSW+ryKN+p5dRwAWEWlgo5aCkJVHSlTpkyRl156SY+2mj59uowdO9ZZZQTgBAlpmU49DgAsP+oqNzdXateurR+rJR/UMhHe3t4yZswY6dChg7PKCMAJQvx9nXocAFi+RkcNLVfhRomIiNBrYClpaWnOKR0Ap1FDyMMDfQs7Hl9M7VfPq+MAwCoqVaPzwQcfyEMPPaRrcWrWrCnt27eXzp07S0JCgrz++uvOKyWASlPz5Ezq31qPulKhpqQux+p55tMB4LbDy++55x55+umnJTo6utj+HTt2yN69e3VTVmRkpHTp0kU8PDzEFTC8HO44xFyNrrq4Y/Jfb2omY/pcZVi5AKAqPr/LFXTsdrvUq1dPj666OOwoagRWRkaG1KpVS1wFQQfuSA0hV6OrVMfjr7Ycl+92JsjADhHy9t1XG100ADB2Hh3VyfjGG2+Ubdu2XfLcqVOnCjsnAzAv1TwV07SuDOgQKSNvaKb3xW6Pl7NZuUYXDQCcqlxBR61fNWvWLB101G3r1q2XHFN0Xh0A5tchqrY0Ca4pmTn5snQry0AAcOOgo1q5VN8btUL5TTfdVGLYYTFPwLWo39k7O9XXj7/YeNTo4gCA8cPLVV+dOXPmSO/evXXYYQFPwLUNvDpS1N8oaw6cliOnM4wuDgAY13R1cdjp06ePrt0h7ACuK7K2n8Q0qasfL9p0zOjiAIBxTVfFXmy3y8cff1wYdjZv3uy8kgGoVnd2vNB8tenYJb/rAOAWQWfJkiV6KFdJYefmm2+WO++809nlA1BNbokOkxreHnIwMV02Hk42ujgAUP1B59ZbbxUfH59Lv4jdLv/9739lwIAB/CUIuKiaPp467Cif0ykZgLsFnfj4eMnKyir9C12o2VmzZo3ePnDggHNKCKDam6++3nJcMnPyjC4OAFRf0FmwYIEEBQXJoEGD5MMPP9STA15s3bp1smjRImnTpo1e9wqAa1EdkiMCfSU1M1eW70wwujgAUH1BZ9SoUbJlyxa55ppr9KSB9evXl169esnLL78sjzzyiISHh8vAgQP1gp6vvPJKiUEIgLnZ7TY91Fyh+QqAFZRrrauikpKS5Ouvv5alS5dKo0aNdP+cmJgYl5swkLWugOL2JZyV3m/9qJeJWDPhJqnnf2m/PACw5KKeVkTQAS41YOoq2XIkWZ7p10qGX9PE6OIAQPUt6gnA+gZ3PN989cVGJg8E4NoIOgAucXu7CPHysMmOE6my80Sq0cUBgAoj6AC4RJ2a3nJTy1D9mIU+Abgygg6AEhWsaL5w03HJzcs3ujgAUCEEHQAluq5FPQmq6S2JZ7Pk532JRhcHACqEoAOgRN6edvlD+wj9+PMNNF8BcE0EHQBXXBLi2x0nJeVcjtHFAYByI+gAKFV0ZIC0CK0l2bn5snTrCaOLAwDlRtABUCo10/kdF2p1aL4C4IoIOgAua9DVkWK3iaz/7Yz8lpRudHEAoFwIOgAuKzTAV3o2C9aPP2emZAAuhqAD4IoGX5hTR00emJ/v1svjAXAxpgs6U6dO1auh+/r6Srdu3WTdunVlet28efN0f4KBAwdWeRkBd3Nz6zCp5eMpR8+ck18OnTa6OADgmkFn/vz5MmbMGJk0aZJs3LhR2rdvL3379pWEhITLvu7QoUMyduxYueaaa6qtrIA78fP2kNvahunHLPQJwJWYKui89dZb8sgjj8iwYcOkdevWMn36dKlRo4Z88MEHpb4mLy9P7rvvPpk8ebI0adKkWssLuOOcOku2npBz2XlGFwcAXCvoZGdny4YNG6R3796F++x2u96Oi4sr9XXPP/+8hISEyMMPP1ym82RlZUlqamqxG4Ar69IoSOrX8ZOzWbny7Y54o4sDAK4VdBITE3XtTGjo+RWTC6jt+PiS31RXrlwp77//vsycObPM55kyZYoEBgYW3qKioipddsAd2O1F5tSh+QqAizBN0CmvtLQ0eeCBB3TICQ4+P/S1LCZMmCApKSmFtyNHjlRpOQEruePqSH2/cu8pOZmaaXRxAOCKPMUkVFjx8PCQkydPFtuvtsPCzneCLGr//v26E3L//v0L9+Xn5+t7T09P2b17tzRt2vSS1/n4+OgbgPJrFFxTOjesoycPXLTpmPz5ukt/xwDATExTo+Pt7S2dOnWS5cuXFwsuajsmJuaS41u2bClbt26VzZs3F97+8Ic/yA033KAf0yQFVI3fm6+OisPBnDoAzM00NTqKGlr+4IMPSufOnaVr167y9ttvS3p6uh6FpQwdOlQiIyN1Pxs1z050dHSx19euXVvfX7wfgPP0axcuz321XfacPCvbj6dKdGSg0UUCANcIOkOGDJFTp07Js88+qzsgd+jQQWJjYws7KB8+fFiPxAJgnEA/L+nTOlSW/HpC1+oQdACYmc3h5nXPani5Gn2lOiYHBAQYXRzAJfywK0GGzfpF6tb0ljVP3SReHvwBAsCcn9+8OwEot2uaB0twLR9JSs+WFbtPGV0cACgVQQdAuXl62GVgh4jChT4BwKwIOgAq5M4LK5ov35kgyRnZRhcHAEpE0AFQIa3CA/QtOy9fvvr1hNHFAYASEXQAVNidHc/PlPz5BpqvAJgTQQdAhQ3oECkedptsPpIs+0+dNbo4AHAJgg6ACqvn7yPXtainH9MpGYAZEXQAVModF5qvFm48Jvn5bj0tFwATIugAqJTerULF39dTjqdkypoDSUYXBwCKIegAqBRfLw+5vd35OXU+33jM6OIAQDEEHQCVNrjT+ear/207IelZuUYXBwAKEXQAVFrHBnWkUd0akpGdJ99sjze6OABQiKADoNJsNpvc0fH8TMlqRXMAMAuCDgCnGHT1+ear1fuT5HjyOaOLAwAaQQeAU0QF1ZBujYPE4RBZuIlOyQDMgaADwOkLfarmK4dKPABgMIIOAKe5NTpMfL3scuBUumw5mmJ0cQCAoAPAefx9veSWNmH6MQt9AjADgg4ApyoYffXVr8clKzfP6OIAcHMEHQBO1bNZsIQG+EhyRo78sCvB6OIAcHMEHQBO5WG3ycALQ81ZEgKA0Qg6AJzuzgvNV6pGJ+lsltHFAeDGCDoAnK5FqL+0jQyU3HyHfLXluNHFAeDGCDoAqsSdHWm+AmA8gg6AKtG/fYR42m2y9ViK7DmZZnRxALgpgg6AKlG3lo/c0DJEP2ahTwBGIegAqPJOyYs2HZO8fJaEAFD9CDoAqswNLetJ7RpecjI1S1btSzS6OADcEEEHQJXx8fSQP7SP0I9pvgJgBIIOgGpZEmLp1hMy/5fDErc/iWYsANXGs/pOBcAdnUg+p2dLzslzyLjPt+p94YG+Mql/a7klOtzo4gGwOGp0AFSZ2G0n5C9zNl5SgxOfkikjPt6onweAqkTQAVAlVLiZ/NUOKamRqmCfep5mLABViaADoEqsO3haTqRklvq8ijfqeXUcAFQVgg6AKpGQlunU4wCgIgg6AKpEiL+vU48DgIog6ACoEl0bB+nRVbbLHOPtYZc2EQHVWCoA7oagA6BKqCHlagi5UlrYyc7Ll2GzfpHUzJxqLRsA90HQAVBl1Dw50+7vKGGBxZunVE3PuFtaSoCvp2z47YzcN3OtnEnPNqycAKzL5nA43HpsZ2pqqgQGBkpKSooEBFCFDlQFNYRcja5SHY9VnxzVrKVqfLYfT5EH3l8np9OzpWWYv3w8vJsE1/IxurgALPT5TdAh6ACG2nsyTe79z1o5lZYlTevVlDnDu19SAwQAFf38pukKgKGah/rLp3+OkYhAX9l/Kl3u+necHD2TYXSxAFgEQQeA4RoH15T5f46RBkE15PDpDLlrepwcSkw3ulgALICgA8AUooJq6JqdJvVqyvGUTF2zo5q1AKAyCDoATEP1zZn/aIzumJyQliVDZqzRHZYBoKIIOgBMpZ6/j3zySHdpGxmoR2PdM2ONbD6SbHSxALgogg4A06lT01vmPNJNOjWsI6mZuXL/f9bKL4dY/BNA+RF0AJhSgK+XfPSnrtK9SZCczcqVoe+vk1X7Eo0uFgAXQ9ABYFo1fTxl1rCucm2LenIuJ08vF/HDrgSjiwXAhRB0AJiar5eHzBzaSfq0DpXs3Hx59L/rJXbbCaOLBcBFEHQAmJ6Pp4e8d19H6dcuXHLyHDJy7ib5cvMxo4sFwAUQdAC4BC8Pu/zr7qvlzo719dpZT87fLJ/+csToYgEwOYIOAJehFgJ9fXA7ua9bA1Gr9P3j81/lo7hDRhcLgIkRdAC4FLvdJi8OjJY/9Wyst5/9crvM+Gm/0cUCYFIEHQAux2azycTbW8nIG5rq7ZeX7pL/+26vOFQ1DwAUQdAB4LJh5+99W8rYm1vo7X9+t0de+2Y3YQeAuYPO1KlTpVGjRuLr6yvdunWTdevWlXrszJkz5ZprrpE6deroW+/evS97PADrGXVjc3mmXyv9eNqK/TL5qx067KgOy3H7k/ToLHWvtgG4H08xkfnz58uYMWNk+vTpOuS8/fbb0rdvX9m9e7eEhIRccvyKFSvknnvukR49euhg9Oqrr8rNN98s27dvl8jISEOuAUD1G35NE/Hx8pCJi7bJrNWHZF9CmuxLSJf41MzCY8IDfWVS/9ZyS3S4oWUFUL1sDhPV86pw06VLF3n33Xf1dn5+vkRFRcnjjz8u48ePv+Lr8/LydM2Oev3QoUPLdM7U1FQJDAyUlJQUCQgIqPQ1ADDOZ+uPyD8W/ColvanZLtxPu78jYQewgLJ+fpum6So7O1s2bNigm58K2O12vR0XF1emr5GRkSE5OTkSFBRU6jFZWVn6P6foDYA13NGxvgT6eZX4XEH4UU1bNGMB7sM0QScxMVHXyISGhhbbr7bj4+PL9DXGjRsnERERxcLSxaZMmaITYMFN1RgBsIZ1B09L8rmcUp9X8eZESqY+DoB7ME3QqaxXXnlF5s2bJwsXLtT9dUozYcIEXc1VcDtyhJlVAatISMt06nEAXJ9pOiMHBweLh4eHnDx5sth+tR0WFnbZ177xxhs66Hz33XfSrl27yx7r4+OjbwCsJ8Tf16nHAXB9pqnR8fb2lk6dOsny5csL96nOyGo7Jiam1Ne99tpr8sILL0hsbKx07ty5mkoLwIy6Ng7So6sKOh5fTO1Xz6vjALgH0wQdRQ0tV3PjzJ49W3bu3CkjRoyQ9PR0GTZsmH5ejaRSTU8F1HDyiRMnygcffKDn3lF9edTt7NmzBl4FACPXwlJDyBVbKX10nr29tT4OgHswVdAZMmSIboZ69tlnpUOHDrJ582ZdU1PQQfnw4cNy4sSJwuOnTZumR2sNHjxYwsPDC2/qawBwT2rouBpCHhZYcvPUqbNZ1V4mAMYx1Tw6RmAeHcCa1BByNbpKdTxWfXK2HkvWa2J5e9jl8xE9pG39QKOLCKAaPr9N0xkZAJxJNU/FNK1buN29SZD8cuiMLNtxUkbO3Shf/7WXBPiWPOcOAOswVdMVAFTlIqBvDG4v9ev4yeHTGfKPz35lAVDADRB0ALiNwBpe8u69HcXLwyax2+P1ulgArI2gA8CtdIiqLRNuPb/a+ctLd8qWI8lGFwlAFSLoAHA7w3o2klvahElOnkP310nJKH3ZCACujaADwC3767w6uJ1EBfnJ0TPnZOyCLfTXASyKoAPALalVzt+7t5Mebq5GYn2wiv46gBURdAC4LTWXztP9zvfXmbJ0p2w6fMboIgFwMoIOALc2NKah9GsbLrn5Dhk1d5MkZ2QbXSQATkTQASDu3l9nyp1tpWHdGnIs+ZyM/Yz+OoCVEHQAuD01Q/LUezvq/jrf7UyQ//x80OgiAXASgg4AiEh0ZKBMvLDy+auxu2TDb/TXAayAoAMAF9zfrYHc3u58f53H526UM+n01wFcHUEHAIr217mjrTQOrinHUzLlb59tkfx8+usAroygAwBF+Bf01/G0y/e7EmTGzweMLhKASiDoAMBFWkcEyHP92+jHr3+zW345dNroIgGoIIIOAJTgnq5RMqBDhOTp/jqbJOlsltFFAlABBB0AKKW/zsuD2kqTejUlPjVTxnxKfx3AFRF0AKAUNX085b37OoqPp11+3HNKpv243+giASgngg4AXEbLsAB5fsD5/jpvfrtb1h5IMrpIAMqBoAMAV3BX5yi54+pIUS1Xj3+ySRLprwO4DIIOAJShv86Lg6KlWUgtSUjLktHzN9NfB3ARBB0AKIMa3uf76/h62eXnvYky9Yd9RhcJQBkQdACgjFqE+ssLA6L1439+t0dW7080ukgAroCgAwDl8MfOUTK4U33dX+eJeZvlVBr9dQAzI+gAQDmpWp0WobV0yHly/iY9qSAAcyLoAEA5+Xl76P46fl4esmpfkrzz/V4dduL2J8mXm4/pe8IPYA42h8Ph1r+NqampEhgYKCkpKRIQEGB0cQC4kC82HtUzJitBNb3ldHp24XPhgb4yqX9ruSU63MASAtZV1s9vanQAoILu6Fhfejatqx8XDTlKfEqmjPh4o8RuO2FQ6QAoBB0AqCDVPLX/1NkSnyuoKp/81Q6asQADEXQAoILWHTwt8amlj7pS8eZESqY+DoAxCDoAUEEJaZlOPQ6A8xF0AKCCQvx9nXocAOcj6ABABXVtHKRHV9kuc0zdmt76OADGIOgAQAV52G16CLlSWthJPpctn64/Uq3lAvA7gg4AVIKaJ2fa/R0lLLB481RYgK90bFBb8vJFJnyxVZ5ZtFWyc/MNKyfgrpgwkAkDATiBGkKuRlepjseqT45qrrLbRN5bsV/e+Ha3qHdatU/NqBxcy8fo4gJu8/lN0CHoAKhiy3ee1AuAns3KlYhAX5kxtLNERwYaXSzApTEzMgCYxE2tQmXRyB7SOLimHE/JlMHTV8viLceNLhbgFgg6AFANmoX4y6KRPeW6FvUkMydf/vrJJnk1dhezJgNVjKADANUk0M9LPnioizx2XVO9PW3Ffnl49i+Sci7H6KIBlkXQAYBqHpI+/taW8n93dxAfT7us2H1KBk1dJfsSSl4zC0DlEHQAwAADOkTK5yN66M7JBxLTddhRnZYBOBdBBwAMokZeLX68l3RtFCRpWbky/KP1MvWHfeLmg2EBpyLoAICB1Jw6Hw/vJvd3b6Dn2nn9m90y6pNNkpGda3TRAEsg6ACAwbw97fLiwLby8qC24uVhkyW/npA7p8XJkdMZRhcNcHkEHQAwiXu7NZC5j3SX4FresvNEqgyYukri9icZXSzApRF0AMBEujQKksWjekl0ZICcTs+W+99fKx/FHaLfDlBBBB0AMJmI2n6y4LEeMqBDhJ5Q8Nkvt8v4z7dKVm6e0UUDXA5BBwBMyNfLQ94e0kGeuq2lXhx0/vojcs+MNZKQmqmfVwFINWt9ufmYvmeGZaBkLOrJop4ATO7HPafk8bkbJTUzV8ICfOWhHo1kdtwhOZFyPvQo4YG+Mql/a7klOtzQsgLVhdXLy4igA8AVHExMl0c+Wl/qDMq2C/fT7u9I2IFbSGX1cgCwDrXy+YLHYvSyESUp+It18lc7aMYCiiDoAICL2HkiTbJy80t9XsUb1Zy17uDpai0XYGaeRhcAAFA2CWm/98m5nLGfbZGYpnWldXiAtI4IkFbhAXrldMAdEXQAwEWE+PuW6bhjyedkwYajxfZF1vbToacg/Kj7+nX8xGYr6N1zeao5TNUUqbClytG1cZBeib2qVPf5jDgn53PToDN16lR5/fXXJT4+Xtq3by/vvPOOdO3atdTjP/vsM5k4caIcOnRImjdvLq+++qrcdttt1VpmAKgO6oNCja6KT8ks7JNTlPoIqefvI5P/0EZ2xafJjhOpsuN4qg4+BbdlO35fIT3A11PX9hQNQM1D/PWSFEXFbjuh+/5U1yiv6j6fEefkfNXHVKOu5s+fL0OHDpXp06dLt27d5O2339ZBZvfu3RISEnLJ8atXr5Zrr71WpkyZIrfffrvMnTtXB52NGzdKdHR0mc7JqCsArkR9gIz4eKN+7CjjqKuUjJzzoedC8FH3+xLSJCfv0rd/tdZWsxD/wuBzNjNX3v5uzyXBqqpGeRVcX3Wdz4hzcr6OTjmfSw4vV+GmS5cu8u677+rt/Px8iYqKkscff1zGjx9/yfFDhgyR9PR0+frrrwv3de/eXTp06KDDUlkQdAC4Gmf8tZydm6+Hqv8eflL0vZqrpzz8fT1lxPVNxdvDLnabTTw9bLp5wsN2/l5t6/12u3jYRTzsdvG028SunlPHFbmpT8ZH/7teEs9ml3q+EH8fmf9ozCVNIGVsgSuxeeWP/46TU2lZlz3nZ49des6Knm/wdPc9n01EwgJ9ZeW4Gyt9PpcLOtnZ2VKjRg1ZsGCBDBw4sHD/gw8+KMnJyfLll19e8poGDRrImDFj5MknnyzcN2nSJFm0aJFs2bKlxPNkZWXpW9H/KBWmCDoAXElV9H9QHweqeaug1mfl3kRZ/9sZp5UZKPDJI911h/nqCDqm6aOTmJgoeXl5EhoaWmy/2t61a1eJr1H9eEo6Xu0vjWrmmjx5spNKDQDGUKGmsh8UF1Mdk+vXqaFvN7cJ03P3lCXodG0UJJF1/CQ33yF5+fk6hKnb+e3iN7Uv3+GQ3LwL90WeS8vMkTMZOVc8n4+Hqimyl9hPqcDl/oR3FHmlOm9JTXgXK6h9qqyC/wN3P19CGUcQOoNpgk51mTBhgq4FurhGBwBQsVFeo/u0cEroUmt23TNzzRWPm/Wnbk4LeWU9538fds45OV/5frYsNWFgcHCweHh4yMmTv48IUNR2WFhYia9R+8tzvOLj46OruIreAAClj/Iq7e98tV89r45zxfMZcU7OJ07/HrpM0PH29pZOnTrJ8uXLC/epzshqOyYmpsTXqP1Fj1eWLVtW6vEAgLJTTRmqg7Ny8QdXwbZ63llzo1T3+Yw4J+cTp38PXSboKKpJaebMmTJ79mzZuXOnjBgxQo+qGjZsmH5eDT1XTU8FnnjiCYmNjZU333xT9+N57rnnZP369TJq1CgDrwIArEON4lLDgdVImaLUdlUM9a7u8xlxTs4XLtXJNKOuCqih5QUTBqph4v/617/0sHPl+uuvl0aNGsmsWbMKj1fz7DzzzDOFEwa+9tpr5ZowkOHlAHBl7jCrrtWvMc9i53O54eVGIegAAGDdz29TNV0BAAA4E0EHAABYFkEHAABYFkEHAABYFkEHAABYFkEHAABYFkEHAABYFkEHAABYFkEHAABYlqe4uYKJodUMiwAAwDUUfG5faYEHtw86aWlp+j4qKsroogAAgAp8jqulIErj9mtd5efny/Hjx8Xf319sNucuNqbC05EjRyy9hpY7XCfXaB3ucJ3ucI3ucp3ucI2VuU4VX1TIiYiIELu99J44bl+jo/5z6tevX2VfX33TrPwD6k7XyTVahztcpztco7tcpztcY0Wv83I1OQXojAwAACyLoAMAACyLoFNFfHx8ZNKkSfreytzhOrlG63CH63SHa3SX63SHa6yO63T7zsgAAMC6qNEBAACWRdABAACWRdABAACWRdABAACWRdCphKlTp0qjRo3E19dXunXrJuvWrbvs8Z999pm0bNlSH9+2bVtZunSpWO06Z86cKddcc43UqVNH33r37n3F/xdX/F4WmDdvnp5Re+DAgWK1a0xOTpaRI0dKeHi4Hg3RokULl/iZLe91vv3223LVVVeJn5+fnp119OjRkpmZKWb1008/Sf/+/fVssOpnb9GiRVd8zYoVK6Rjx476+9isWTOZNWuWmFl5r/GLL76QPn36SL169fSEczExMfLNN9+I2VXke1lg1apV4unpKR06dBCrXWNWVpY8/fTT0rBhQ/0zq36fP/jggwqXgaBTQfPnz5cxY8boIXEbN26U9u3bS9++fSUhIaHE41evXi333HOPPPzww7Jp0yb9wahu27ZtEytdp3pDVdf5ww8/SFxcnP7guPnmm+XYsWNilWsscOjQIRk7dqwOdmZX3mvMzs7WHxzqGhcsWCC7d+/WITYyMlKsdJ1z586V8ePH6+N37twp77//vv4aTz31lJhVenq6vi4V6Mri4MGD0q9fP7nhhhtk8+bN8uSTT8rw4cNNHQTKe43qw1T9vKogvmHDBn2t6sNVvdeaWXmvs+gfIUOHDpWbbrpJzC69Atd41113yfLly/Xvo3rv+eSTT/QfIxWmhpej/Lp27eoYOXJk4XZeXp4jIiLCMWXKlBKPv+uuuxz9+vUrtq9bt26OP//5zw4rXefFcnNzHf7+/o7Zs2c7rHSN6rp69Ojh+M9//uN48MEHHQMGDHCYWXmvcdq0aY4mTZo4srOzHa6kvNepjr3xxhuL7RszZoyjZ8+eDleg3sIXLlx42WP+8Y9/ONq0aVNs35AhQxx9+/Z1WOUaS9K6dWvH5MmTHa6iPNepvn/PPPOMY9KkSY727ds7rHSN//vf/xyBgYGOpKQkp52XGp0KUH/tqr8aVLNM0TWz1LaqxSiJ2l/0eEX9pVna8a56nRfLyMiQnJwcCQoKEitd4/PPPy8hISG6hs7sKnKNixcv1tX/qukqNDRUoqOj5eWXX5a8vDyx0nX26NFDv6ageevAgQO6VuC2224Tq3DF9x5nLNasFns06/tOZXz44Yf651TVQlrR4sWLpXPnzvLaa6/pGmTVZK5qzs+dO1fhr+n2i3pWRGJion7DVx8ARantXbt2lfia+Pj4Eo9X+610nRcbN26cbpu9+I3Wla9x5cqVukpVNQO4gopco3oj/f777+W+++7TH/z79u2Tv/zlLzq0mvUNtiLXee+99+rX9erVS6+EnJubK4899pipm67Kq7T3HrVitPrwUH2TrOaNN96Qs2fP6iYQK9m7d69uav355591/xwrOnDggH6PVX3sFi5cqH8/1XtPUlKSDnkVQY0Oqswrr7yiO+uqH1b1Q2sF6q/EBx54QPdXCQ4OFiv/RaxqrGbMmCGdOnWSIUOG6M6B06dPFytRfcpUTdV7772n+/SoTq1LliyRF154weiioYJUv6vJkyfLp59+qn+GrUKFeBXM1bWpWg4rv/fYbDaZM2eOdO3aVdeuvvXWWzJ79uwK1+pYMxJWMfUB5+HhISdPniy2X22HhYWV+Bq1vzzHu+p1Fv2LSgWd7777Ttq1aydWucb9+/frDrqqo2PRX0xF/YWlOs41bdpUXP37qEZaeXl56dcVaNWqla4dUE1E3t7eYjYVuc6JEyfq4Ko65ypqNKTqPPnoo4/qYKeavlxdae89anSS1Wpz1B9W6nupRriatRa5Mn9krV+/XnewHjVqVOF7j6qJVO893377rdx4443i6sLDw3WTVWBgYLH3HnWdR48elebNm5f7a7r+b7EB1Ju8+itX9QovoH7g1Lbq11AStb/o8cqyZctKPd5Vr1NRbavqL+LY2Fjd1mpm5b1GNT3A1q1bdbNVwe0Pf/hD4YgWNcrMCt/Hnj176uaqghCn7NmzR78JmTHkVPQ6VR+yi8NMQbizyjKArvjeUxFqZM6wYcP0vRplZjUqmF783qOaWdVoJPVYTaVgBT179pTjx4/rpsei7z3q97R+/foV+6JO69bsZubNm+fw8fFxzJo1y7Fjxw7Ho48+6qhdu7YjPj5eP//AAw84xo8fX3j8qlWrHJ6eno433njDsXPnTt1b3svLy7F161aHla7zlVdecXh7ezsWLFjgOHHiROEtLS3NYZVrvJgrjLoq7zUePnxYj5YbNWqUY/fu3Y6vv/7aERIS4njxxRcdVrpO9XuorvOTTz5xHDhwwPHtt986mjZtqkdJmpX6Xdq0aZO+qbfwt956Sz/+7bff9PPq+tR1FlDXVaNGDcff//53/d4zdepUh4eHhyM2NtZhlWucM2eOfn9V11b0fSc5OdlhZuW9zou5wqirtHJeozq+fv36jsGDBzu2b9/u+PHHHx3Nmzd3DB8+vMJlIOhUwjvvvONo0KCB/mBXw1rXrFlT+Nx1112nPwCL+vTTTx0tWrTQx6vhnkuWLHFY7TobNmyof5gvvqlfSCt9L10t6FTkGlevXq2nQFDBQQ01f+mll/SweitdZ05OjuO5557T4cbX19cRFRXl+Mtf/uI4c+aMw6x++OGHEn/HCq5L3avrvPg1HTp00P8n6nv54YcfOsysvNeoHl/ueCt9L10t6PxQgWtUgbx3794OPz8/HXrUlA8ZGRkVLoNN/eO8SicAAADzoI8OAACwLIIOAACwLIIOAACwLIIOAACwLIIOAACwLIIOAACwLIIOAACwLIIOAACwLIIOAACwLIIOAACwLIIOALd2/fXXy5NPPlnpr6NW03nrrbekcePGUqNGDRk4cKCkpKQ4pYwAKo6gA8BprrvuOrHZbPrm7e0trVq1krlz54o7+Pvf/y7Tpk2T2bNny88//ywbNmyQ5557zuhiAW6PoAPAKVSNxqZNm+SNN96QEydOyO7du+WWW26RoUOHysGDB8XK1q5dq2tz5s+fL9dee6106tRJHnnkEVm6dKnRRQPcHkEHgFPs3btX0tLSdLgJCwvTTTgPP/yw5OXl6dBTXrGxsdKrVy+pXbu21K1bV26//XbZv39/sSanv/71r/KPf/xDgoKC9DkvrkFR5bnvvvukZs2aEh4eLv/85z+v2FSVn58vU6ZM0eX38/OT9u3by4IFCy5bVhXubrrpJunYsWPhvtDQUElMTCz3dQNwLoIOAKdQTTV16tSR1q1b6+2jR4/K008/LT4+PtKuXbtyf7309HQZM2aMrF+/XpYvXy52u10GDRqkg0gB1UykQoyqUXnttdfk+eefl2XLlhU+r16/atUqWbx4sd6vmpQ2btx42fOqkPPRRx/J9OnTZfv27TJ69Gi5//775ccffyzx+KysLFmyZIkuW1GZmZkSGBhY7usG4GQOAHCCsWPHOux2u6NmzZoOX19fh3p78fPzc3z44YdO+fqnTp3SX3Pr1q16+7rrrnP06tWr2DFdunRxjBs3Tj9OTU11eHl5OT777LPC55OTkx01atRwPPHEE4X71Ncp2M7MzNTPr169utjXffjhhx333HNPieVSx6pyqWtW115w8/b2dvTt21cfs2jRomLnBFB9PJ0dnAC4J1VTMnLkSN2clJycLGPHjpWePXvKQw89VOGmsGeffVbX1qgmoIKanMOHD0t0dLR+fHFNkWqeSkhI0I8PHDggOTk50rVr18LnVQ3LVVddVeo59+3bJxkZGdKnT59i+7Ozs+Xqq68u8TV79uzRtUqbN28utr9fv376+pVff/1VOnToUM7/AQDOQNAB4LSgozrgNmvWTG+/9957OoiofY0aNdL71Ae/6rTbrVs33X9HBRbVNFSS/v37S8OGDWXmzJkSERGhg446XoWOAl5eXsVeo0Z7FW3aKq+zZ8/qe9UUFRkZWew51QRXktTUVAkODi68buW3337TQe3OO+8sDDrnzp3T1686aqumtIKwBqBq0UcHQKWp2hNVi1P0w1v11WnatGmx4eUTJ06UV155RYcd1eemtJCTlJSkOzA/88wzupOvGqZ+5syZcpWpSZMmOgj98ssvhfvUvDaqBqY0qswq0KhaIxVcit6ioqJKfI0KOerrqlFnBV566SW57bbbCvsrqaCjwp7qL6RqvL788styXQuAiqNGB4BTOiKrUNGiRYti+1VIWbhwoTz11FN6W43IUh2UVY2JGlVVGtWpWY20mjFjhm6OUsFj/Pjx5SqTv7+/PPjgg3p+GzUqKyQkRCZNmqQDlqr5Ke01qslNBTBVM6RGfakQowJKQECA/noXu/HGG3XHYxXg7r77bpkzZ4589dVXsm7dOv28qslRI88effRRva2a09RIMgDVgxodAE5ptmrevLmeJLCo3r176xCkRmApqnbl9OnTuq/Mxc1ORakwMm/ePP3aguat119/vdzlUjVHMTExemi6KotqOlK1Q76+vqW+5oUXXtA1T2r0lTpWhTMVzNRw85KoYeSzZs3SkwW2adNG1qxZIytXriysAdq2bZt07ty58PitW7fq4wBUD5vqkVxN5wLgxo4dO6abcxYtWqT7rqgh3NXdT0UNWVd9b958803dR6g6vP/++xIfH69rshTVqfm7777TNVYAqh41OgCqnGq++eMf/yjvvPOOrhmZMGGCrjmpamqm5k8++URPNKhqndTkgcqAAQOkuqj+OW3bttWPc3NzdV8mQg5QfajRAWBZKugMHz5cd2xWzWpqaQbVnFUQPABYH0EHAABYFk1XAADAsgg6AADAsgg6AADAsgg6AADAsgg6AADAsgg6AADAsgg6AADAsgg6AADAsgg6AADAsgg6AABArOr/ATN/+dVE0aYlAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.xlabel(r\"$R_x$ angle $\\theta_h$\")\n",
    "plt.ylabel(r\"$\\langle Z_{62} \\rangle$\")\n",
    "plt.plot(hs, approx_evs, marker=\"o\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}