{ "cells": [ { "cell_type": "markdown", "id": "cac2fea9", "metadata": {}, "source": [ "# Comparing the splitstep and plane-wave methods" ] }, { "cell_type": "markdown", "id": "1201dfaa-5042-486d-859f-d9f0af332183", "metadata": {}, "source": [ "**What you'll learn:**\n", "- How to run a split-step (real-space lattice) simulation of a full interferometer sequence\n", "- How to compare split-step results against the plane-wave (Bloch) method\n", "- How momentum width affects the agreement between the two approaches" ] }, { "cell_type": "markdown", "id": "feaf49d6-98c3-4825-a86c-e3d2e0892a64", "metadata": {}, "source": [ "I've seen that as I make the wavefunction narrower in momentum space the splitstep simulation moves towards the plane wave simulation.\n", "\n", "I will investigate this effect in more detail by simulating an ensemble of atoms with the plane wave simulation that have a momentum spread equivalent to the wavefunction width in the splitstep simulation.\n", "\n", "Lets start with the splitstep simulation:" ] }, { "cell_type": "code", "execution_count": 1, "id": "d657d8fd-fcba-4c54-a792-d105228751bb", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 3107/3107 [00:02<00:00, 1043.72it/s]\n", "100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 19999/19999 [00:17<00:00, 1170.78it/s]\n", "100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 10000/10000 [00:08<00:00, 1157.81it/s]\n", "100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 19999/19999 [00:17<00:00, 1137.69it/s]\n" ] } ], "source": [ "import numpy as np\n", "from mwave.integrate import make_kvec, make_phi\n", "from matplotlib import pyplot as plt\n", "from mwave.integrate.xintegrate import splitstep, x_to_kvec, kvec_to_x, xspace_to_pspace, make_continuous_phi, pspace_to_xspace\n", "from numba import jit\n", "\n", "# Define fixed parameters\n", "n_init = 0\n", "n_bragg = 4\n", "omega = 19.4\n", "sigma = 0.259\n", "tau_factor = 3\n", "T = 10\n", "Tp = 5\n", "target_phase = 0*np.pi/2\n", "mod_phase = 0\n", "\n", "# Compute dependent parameters\n", "delta = 4*n_bragg\n", "t_bragg = 2*tau_factor*sigma\n", "t_spacing = T - t_bragg\n", "tp_spacing = Tp - t_bragg\n", "omega_m = 8*n_bragg-target_phase/(2*T*n_bragg)\n", "\n", "tseg0 = 0\n", "tseg1 = t_bragg\n", "tseg2 = 2*t_bragg + t_spacing\n", "tseg3 = 3*t_bragg + t_spacing + tp_spacing\n", "tseg4 = 4*t_bragg + 2*t_spacing + tp_spacing\n", "\n", "# Define kvec\n", "kvec, n0_idx, nf_idx = make_kvec(n_init, n_init + n_bragg)\n", "\n", "# Make splitstep parameters\n", "klaser = np.sqrt(2)\n", "xsig = 15\n", "xextentneg = 1.2*(T+2*t_bragg)*-n_bragg*2*klaser\n", "xextentpos = 1.2*((T+2*t_bragg+T)*n_bragg*2*klaser + (2*t_bragg+T)*2*n_bragg*2*klaser)\n", "xvec = np.arange(xextentneg,xextentpos,0.1)\n", "psi0 = np.exp(-np.square(xvec)/(2*xsig**2))/np.sqrt(xsig*np.sqrt(np.pi))*np.exp(1j*n_init*2*klaser*xvec)\n", "ckvec = np.fft.fftfreq(len(xspace_to_pspace(psi0)), np.diff(xvec)[0])*2*np.pi/klaser\n", "\n", "# Compute with splitstep\n", "@jit(nopython=True)\n", "def OmegaProfile(t):\n", " if t < tseg1: # Bragg 1\n", " return omega*np.exp(-np.square(t-(tseg1-t_bragg/2))/(2*np.square(sigma)))\n", " elif t < tseg2 - t_bragg:\n", " return 0.0\n", " elif t < tseg2: # Bragg 2\n", " return omega*np.exp(-np.square(t-(tseg2-t_bragg/2))/(2*np.square(sigma)))\n", " elif t < tseg3 - t_bragg:\n", " return 0.0\n", " elif t < tseg3: # Bragg 3\n", " return 2*np.cos(omega_m*(t-tseg3-t_bragg) + mod_phase)*omega*np.exp(-np.square(t-(tseg3-t_bragg/2))/(2*np.square(sigma)))\n", " elif t < tseg4 - t_bragg:\n", " return 0.0\n", " elif t < tseg4: # Bragg 4\n", " return 2*np.cos(omega_m*(t-tseg3-t_bragg) + mod_phase)*omega*np.exp(-np.square(t-(tseg4-t_bragg/2))/(2*np.square(sigma)))\n", " else:\n", " return 0.0\n", "\n", "# Set x grid, potential, and phi0\n", "@jit(nopython=True)\n", "def Vfnc(t):\n", " return OmegaProfile(t)*2*(np.square(np.cos(klaser*xvec-delta*t/2)) - 0.5) # minus 0.5 to remove AC stark shift\n", "\n", "psi1 = splitstep(xvec, tseg1, 0.0005, Vfnc, np.copy(psi0), klaser, tstart=tseg0, store_hist=False)\n", "psi2 = splitstep(xvec, tseg2, 0.0005, Vfnc, psi1, klaser, tstart=tseg1, store_hist=False)\n", "psi3 = splitstep(xvec, tseg3, 0.0005, Vfnc, psi2, klaser, tstart=tseg2, store_hist=False)\n", "psi4 = splitstep(xvec, tseg4, 0.0005, Vfnc, psi3, klaser, tstart=tseg3, store_hist=False)" ] }, { "cell_type": "code", "execution_count": 2, "id": "09017626-27e4-47e8-9fdd-7b31573bebe3", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGfCAYAAACneiONAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAATvhJREFUeJzt3Qd4FGX+B/BfeqihhF4DRIo0aRFEsXCA4h0oInIqyJ8D5UThUDhAigcoiqJ0ERWVE4SLIiAgLQiCiYEktFBDCQmENEIqpO//eV/cJZvMJrvZMu/M+/08zxJ29t3Zd2dmZ37zVjeDwWAgAAAAAI1zVzsDAAAAAI6AoAYAAAB0AUENAAAA6AKCGgAAANAFBDUAAACgCwhqAAAAQBcQ1AAAAIAuIKgBAAAAXUBQAwAAALqAoAYAAAB0wbMyb1q5ciV99NFHlJiYSF26dKHly5dTr169LKYPDg6m2bNnU2xsLAUGBtKHH35ITz31lOn1zZs30+rVqykyMpLS0tLo2LFj1LVrV9PrbNncuXNpz549FBcXR/Xq1aOhQ4fS/Pnzyc/Pz6o8FxcXU0JCAtWoUYPc3Nwq87UBAADAxdhsTllZWdS4cWNyd6+gLMZgo40bNxq8vb0Na9euNZw+fdowbtw4Q61atQxJSUmK6X///XeDh4eHYdGiRYYzZ84YZs2aZfDy8jKcOnXKlGbdunWG//znP4YvvviCzUNlOHbsmNk6WNpnn33WsG3bNsPFixcNISEhhsDAQMOwYcOsznd8fDxfNx7YBjgGcAzgGMAxgGOANLcN2HW8Im7sH1sipqCgIOrZsyetWLHCVALSrFkzeuONN2j69Oll0o8YMYJycnJo+/btpmUPPvggL4lhpTMlsZKcgICAMiU1lkp/XnrpJb5uT8+KC5wyMjKoVq1aFB8fTzVr1rThGwMAAIBaMjMzeZyRnp5eYe2MTdVP+fn5vIpoxowZpmWsKKh///4UFham+B62fMqUKWbLBg4cSFu2bCF7sCCFBSeWApq8vDz+MGJFVwx7D4IaAAAAbbGm6YhNDYVTU1OpqKiIGjRoYLacPWfta5Sw5baktzYfrD3N+PHjLaZZuHAhj+iMDxblAQAAgH65a7EYavDgwdShQwd69913LaZjpUmsNMf4YNVOAAAAoF82VT/5+/uTh4cHJSUlmS1nzxs2bKj4HrbclvTlYVVIgwYN4j2YfvrpJ/Ly8rKY1sfHhz8AAABADjaV1Hh7e1P37t0pJCTEtIw1FGbPe/furfgetrxkembv3r0W05dXQjNgwACeh23btpGvr69N7wcAAAB9s3mcGtbod/To0dSjRw8+Ns2SJUt4D6QxY8bw10eNGkVNmjThbVqYSZMmUb9+/Wjx4sW82mjjxo0UERFBa9asMRuHho0/w8aRYc6fP8//stIc9jAGNLdv36bvvvuOP2cPho1Zw0qPAAAAQG42BzWsi3ZKSgrNmTOHN/ZlXa937dplagzMgpOSg+P06dOHNmzYQLNmzaKZM2fywfdYz6eOHTua0rCSF2NQxLzwwgv8Lxtwj7WbiYqKovDwcL6sTZs2Zvm5cuUKtWzZsjLfHQAAAHTE5nFqtIqV7LBeUMau4AAAAKCv67fmej8BAAAAKEFQAwAAALqAoAYAAAB0AUENAAAA6AKCGgAAANAFBDUAAJVwNDaNvvvjKknSgRRAn+PUAAAA0fDVYXwzBPhXo4fa+GOTAAgAJTUAYLfcgiJpSyyupOaonQUA+BOCGgCwS3zabWo3exe98f0xbEkAUBWCGgCwy3//uMr/bj95A1sSAFSFoAYAwA5yVroBiAlBDWi2Dcdzn4VSy+k76O3gE9K25wAAgHsQ1IAm/RB5jSKu3jL9H401xRARm6Z2FgBAYghqQJPyCovNnhcWo6RGBM/92c1ZKiglBBAGghrQJHc3tXMAMvt07wXT/9eFXaW8wiJV8wMAdyGoAU1CTANqScvJp6UhMabnMcnZtPLXS9ghAAJAUAOa5OaGsAbUkV+q6pMJv3xTlbwAgDkENaBJiGnEIVvPM6VjT7JNACAsBDWgSSipAQCA0hDUgCah8glEYsAQfABCQFADurAhPE7tLIDEATVGFAAQA4Ia0Jz02/k0a0u02bJvQmMp406BankCucnWrghAVAhqQHPmbT+juPyb32NdnheQsH2TUkNhNfIBAGUgqAHNuZSSo7i8sLhsV1sAV0BBDYAYENQAANgJ1U8AYkBQA5ojWWWH8GS7oLvhCAQQFoIa0A3Jrq2gEtmaEAFoCYIaAAA7IZ4GEAOCGgAAGygV1KCUEEAMCGoAwC5xabexBQFACAhqQHPQpkEc3x+Jo92nk0h2mCYBQAwIagCg0paFxEi39aQbbBBAQxDUAAAAgC4gqAEAsBMaCgOIAUENaI6lwn+0awC1BhtEUAMgBgQ1AAB2wjg1AGJAUAO6geHrwdmKig302YFLJPtUEQCiQlADuul9kp1X6PK8gFx+jLxGXx6+Umb5ucQs2nLsuip5AoB7ENSAbnwTGkv5hcVqZ0Ma8Wm36UZGLsnkUmq2xdcmbzru0rwAQFkIakBX0nLy1c6CNJ5ceohk4+mOMWoARIagBgAqRcbqPg93nDIBRIZfKGhOeffK6NYNah17AKA+BDUA4FBTg0+gNxAAqAJBDegKetaqLzjyGoVdukl6hGmfAMSGoAY0J+LqLbWzABW4nV+ky22EsZAAxIagBgAAAHQBQQ3oCsZ1BQCQV6WCmpUrV1LLli3J19eXgoKC6MiRI+WmDw4Opnbt2vH0nTp1op07d5q9vnnzZhowYADVrVuXjxZ7/HjZQaxyc3Pp9ddf52mqV69Ow4YNo6SkpMpkHwCcTK/BJdrUAOgsqNm0aRNNmTKF5s6dS1FRUdSlSxcaOHAgJScnK6YPDQ2lkSNH0tixY+nYsWM0dOhQ/oiOjjalycnJob59+9KHH35o8XP/9a9/0c8//8wDpIMHD1JCQgI9++yztmYfNC4zt6Dc1yNi01yWF5mlZuepnQUAgDLcDDbOxMZKZnr27EkrVqzgz4uLi6lZs2b0xhtv0PTp08ukHzFiBA9atm/fblr24IMPUteuXWn16tVmaWNjYykgIIAHP+x1o4yMDKpXrx5t2LCBnnvuOb7s3Llz1L59ewoLC+Prq0hmZib5+fnxddWsWdOWrwwCefmrcDoUk2rxdS8PN4p57ymX5klGLafvKPf1L0b1oL90aEB6sywkhj7Ze8Hi67EfDHZpfgBkkGnD9dumkpr8/HyKjIyk/v3731uBuzt/zoILJWx5yfQMK9mxlF4J+8yCggKz9bDqrObNm1tcT15eHt8QJR+gfeUFNAx6p4AzYfA9ALHZFNSkpqZSUVERNWhgfgfGnicmJiq+hy23Jb2ldXh7e1OtWrWsXs/ChQt5ZGd8sNIkAHANGwuANQNtagDEptveTzNmzOBFVcZHfHy82lkCF8A0CQAA8vK0JbG/vz95eHiU6XXEnjds2FDxPWy5LektrYNVfaWnp5uV1pS3Hh8fH/4Auei0gAAAABxdUsOqgLp3704hISGmZayhMHveu3dvxfew5SXTM3v37rWYXgn7TC8vL7P1nD9/nuLi4mxaD+if3mKaA+eTqffCEPr9YvltiUTDhmbQI71+LwApS2oY1p179OjR1KNHD+rVqxctWbKE924aM2YMf33UqFHUpEkT3qaFmTRpEvXr148WL15MgwcPpo0bN1JERAStWbPGtM60tDQeoLBu2saAhWGlMOzB2sSwLuHss+vUqcNbP7PeViygsabnE8hDb205Xvn6KP/74pfhmupZo7f9AAA6DWpYF+2UlBSaM2cOb6TLul7v2rXL1BiYBSesR5RRnz59eFfsWbNm0cyZMykwMJC2bNlCHTt2NKXZtm2bKShiXnjhBf6XjYXz7rvv8v9/+umnfL1s0D3Ws4n1oFq1apV93x50B5dSAAB52TxOjVZhnBo5xkdhtFSiYcv3Fel7VbQf1rzcnQbcb327Oa1Y+etF+mj33ZJkJSLtIwC9cNo4NSCu5MxcOh6frnY2AHQNTWoAdFb9BGLq9f7dRtQ/T+xLnZr6qZ0dkJwUxb8AIByU1OjMUcx9BAAAkkJQowOSNIsCUF1xcfm/NUz0CaAuBDU6CGhGrT2idjYAdC/0Uip9vMfyZJZMxp3yZ5EHAOdCUKNx6bcLKpzkEQDs969NxytMg6H5ANSFoEZnUBGlX1dv5qidBajAhvA4bCMAFSGoAdCIfh8dUDsLUnOzohzmy8NXXJIXAFCGoAZ0p6ComPTg1LUMtbMAJWCMGgDxIagB3fkp6jrpwV9XHFY7CwAAmoKgBnTn1u18tbMgPT2OMoBGwADiQ1ADoCEpWXlqZ0Fabqh/AhAeghqNw3m2LB0WEpj0fG+f2lkAABAWghqN02Mxv72wTQAA5ISgBnTHoOuyGqJfTt0g8el7HwCAmBDU6AzmgdJ/Sc2E9VGUmJGrdjZAEpm5BZSVi+kfQBsQ1OjMgh1nKR29f3QvLQc9vFxNxvZr+YXF1PndPdTp3T1UVMFkngAiQFCj08BGZjJefMSjv50g43FVctbxnPxCVfMCYA0ENTp0OSWbZGbNcPbgbPq7q5f9uMor0MdI3aBvCGo0Tsa7RwBQZzgBVHuC6BDU6LBRrP7ukW2DQA9wXNnnZnYezd9+hi4kZZkt3306EQcXCM1T7QwAOJoMhVdqdlsvfaGThQzHldG/fzxJ+84m01eYdRw0BiU1OqT3Ls2g7pABAz79zYp0LskOOMmp65ghHrQJQY0OyX49QfWT88jcqxdzPwGID0EN6I4MvVRk+I6iwRYHEB+CGo27cjOn7ELJy/5RUgNg528IIRxoFIIaDbuYnE3PrgpVOxsAckBRjez3S6ABCGo0LPRSqtpZAJCGtTFNQREGqQNQC4IaHZL9ZkqGBp16n4lcy4IjrpHWWfoJSfDTAo1DUAO6k52LOWqcBdc02+ZL0ipUM4FWIajRIdlPSJ/uu6B2FnRL5kNLhhJAAK1DUAOasT78qtpZAImDH2tDGoQ+AOpBUKNDem1v8c5P0WpnQRiid7n95/ooyrhdoHY2oJIsFUrJXgoM4kNQo0M48YAIPt5znvQEtU8A4kNQA6BBWiiNS8zMJRlLxxD8AKgHQQ0AOEWxzBNFadgvp27QjQzlgBQBG4gOQY0OJWfl8dmUQb92RydiH7uYLBf0Ceuj1M4CQKUhqNGhlKw8em/HWbWzAU60bP9F2nc2GdtYQOj6DaAeBDU69eXhK2pnAZwsKu6Wy7cxSgDlhgJgEB2CGgAAK6AEhijiahraSoHQENQAAFhBkiY15docdZ2CI+PVzgaARQhqAAS0dF+M2lmAEqKvZ9CZG5nYJkS09XgCtgMIC0ENgEbnr0L7Btd5evlhF34aAFQWghoAcApZukCXFhyB6hkAtSCoAV1Ky8lXOwvSk7UkKfbmbdIzWfcraAOCGtCludtOq50FXcL1DABEhqAGdOlySjbpnejVOyHnkumSBPsBAMSBoAZ0SYZph7RQDfDCmj/UzgIASKRSQc3KlSupZcuW5OvrS0FBQXTkyJFy0wcHB1O7du14+k6dOtHOnTvLjFI6Z84catSoEVWpUoX69+9PMTHmXVovXLhAQ4YMIX9/f6pZsyb17duXfv3118pkHySAyRTFmbIDAEDYoGbTpk00ZcoUmjt3LkVFRVGXLl1o4MCBlJysPA9NaGgojRw5ksaOHUvHjh2joUOH8kd0dLQpzaJFi2jZsmW0evVqCg8Pp2rVqvF15ubemyn26aefpsLCQtq/fz9FRkbyz2XLEhMTK/vdQceKtVCMAQAA6gY1n3zyCY0bN47GjBlDHTp04IFI1apVae3atYrply5dSoMGDaKpU6dS+/btaf78+dStWzdasWKFqZRmyZIlNGvWLF4S07lzZ1q3bh0lJCTQli1beJrU1FRecjN9+nT+emBgIH3wwQd0+/Zts+CopLy8PMrMzDR7AAAAgH7ZFNTk5+fzUhJWPWRagbs7fx4WFqb4Hra8ZHqGlcIY01+5coWXtpRM4+fnx6u1jGnq1q1Lbdu25cFOTk4OL7H5/PPPqX79+tS9e3fFz124cCFfj/HRrFkzks1aiSe1FL0RLQAAqBzUsBKToqIiatCggdly9txSNRBbXl5649/y0rCJ5Pbt28err2rUqMHb5rASo127dlHt2rUVP3fGjBmUkZFhesTHyzcg1rztZ9TOAjiRAR2sQQU47kBknqQBrIrq9ddf5yUzhw4d4o2Jv/zyS/rrX/9KR48e5Q2MS/Px8eEP0Ifb+YVqZwE00uMKAORlU0kN63nk4eFBSUlJZsvZ84YNGyq+hy0vL73xb3lpWOPg7du308aNG+mhhx7ibXJWrVrFg5tvv/3Wlq8AGnXrdoHaWRCOG+aNBgCofFDj7e3N27CEhISYlhUXF/PnvXv3VnwPW14yPbN3715T+oCAAB68lEzDGvWyXlDGNKxBMM+su3l22XP2+aB/RUUoIgAAAAdXP7Hu3KNHj6YePXpQr169eM8l1niX9YZiRo0aRU2aNOENdZlJkyZRv379aPHixTR48GBe2hIREUFr1qwxtZeZPHkyLViwgPdqYkHO7NmzqXHjxrzrN8OCG9Z2hn0uG8+GldB88cUXvJExW6esMiQqvSiysd4DpRgAAPKxOagZMWIEpaSk8OCCNeTt2rUrb7BrbOgbFxdnVqLSp08f2rBhA++yPXPmTB64sK7aHTt2NKWZNm0aD4zGjx9P6enpfGA9tk7WINhY7cWev/POO/T4449TQUEB3X///bR161Y+Xo2Mdp9OpMV7L5AsWLsqm9JL0IhWhu+o5eOV3bABgAYaCk+cOJE/lBw4cKDMsuHDh/OHJezHP2/ePP6whJUM7d69uzLZ1aX3d54lmeDyDVrCYnDENACuh7mfNMpDsjOmrb1uUP0EakrW8fQQ6AEHIkNQo1Hu7nIFNQBaMmPzSbWzACAlBDUaJV9MY2NDYem2j2ugHY91rt6822MTAFwLQY1GuUt21ba1yBtF5KAmtAEDUAeCGo2SLqhROwMAACA8BDWgz4bCcsV8oPEhCESRlJmrdhYA7IKgRqO0ecqsPLTlAHC+oPfNR38H0BoENaBL5xKz1M6CLp1JyFQ7C5og200HgCgQ1IAmaLQ037lVFy7eJpdTsumZVaGu/VCN0vPxquOvBjqAoAY0QaYRlFcduEQiOnktQ+0sAACUC0ENaMKhmFSSxYr9F9XOAtgJbcAA1IGgBkCr0MNLWHqufgIQGYIajcL1DEBcCGoA1IGgRqPO3EAvFL26U1BkVbrCItcWB6BKBQBEh6AGQKO+OnyFwi7dVDsbIJm8wmK1swBgEYIaAA2b+sMJtbMAOhpR2Bon4tPph8hramcDQBGCGg36/kic2lkAsNqdfOuq0/REvyHNXW8HI5gGMSGo0aAZm0+pnQUQhCsLBNwq2Tz9udXyDdin44IaAKEhqAEApzqNqRUAwEUQ1OhccTFuGcEx0PsJAESHoEbnZCz6BwAAOSGo0bmouHS1swAgHZRqAagDQQ0AAADoAoIaAAAHQ+8nAHUgqAEQyK2cfJvS63mQNwAAWyGoAU3bN6Wfxdei4m6R1lxMyVY7C+AAbphxFkAVCGpA09rUr27xtWdXaa/nl7uNV0OU0zhf6KVUm9+TlJnnlLwAQPkQ1AAIxMMdt/ii+fsX4ZV6X0oWAhsAV0NQAyAQD9Rb6MbVmzlqZwFAOghqAATiLvAvEm2SbVOI0bwBXE7gUyiAfCo7aSSIpwhBDYDLIagB0DCUngAA3IOgBjTL1+vu4Tvh0dZqZ0UKaO5jGwScAK6HoAaEZ2mAuchZf+F/pw1sS0EBdUhGmGMIAOAeBDUgvN2nkxSXV/Px5H/d3NyobnVvF+dKPih5AADRIagB4cWn3ZbmghuTnKV2FgAANAtBDQiv2IqIRS9BzayfotXOAjjI9pMJ2JYALoagBoRnTc9YvbQtsSaAA23YeDRe7SwASAdBDQhPppKanPwim9Lr5Xtr1Wv90PMOQCQIakCzvZ/M0rgkJwDmpj/ZDpsEQCAIakAf1U+SRjWSfm0AAEUIakB4GG4eAACsgaAGhHft1h0pyiyy8wrVzgIAgKYhqAHh/Rh1TYrqpzUHL6mdBQAATUNQA7qgg5iGMu4UkMj0EDiCZZFX07B5QPMQ1GhMboFtXX5lYU0PKT2S9GuDEwz7LAzbFTQPQY3GtJu9S+0sCAnXdnC1957piI0OoIegZuXKldSyZUvy9fWloKAgOnLkSLnpg4ODqV27djx9p06daOfOnWXusufMmUONGjWiKlWqUP/+/SkmJqbMenbs2ME/j6WpXbs2DR06tDLZB40Z0KGBFCUWbGJOkQmePZdq26AGvRjUQu1sAIC9Qc2mTZtoypQpNHfuXIqKiqIuXbrQwIEDKTk5WTF9aGgojRw5ksaOHUvHjh3jgQh7REffm+Nm0aJFtGzZMlq9ejWFh4dTtWrV+Dpzc3NNaX788Ud6+eWXacyYMXTixAn6/fff6e9//7ut2QcN8vKo+DDVQUwDGuLr7WH6f+Ss/qrmBQDsCGo++eQTGjduHA8uOnTowAORqlWr0tq1axXTL126lAYNGkRTp06l9u3b0/z586lbt260YsUKUynNkiVLaNasWTRkyBDq3LkzrVu3jhISEmjLli08TWFhIU2aNIk++ugjeu211+i+++7jn/3888/bmn3QoLM3MitMI2ubGpSeqLTdS/y/bnUfktG1W7fVzgKAfUFNfn4+RUZG8uoh0wrc3fnzsDDlRmZsecn0DCuFMaa/cuUKJSYmmqXx8/Pj1UzGNKxE6Pr16/yzHnjgAV5N9eSTT5qV9pSWl5dHmZmZZg/QnqTMXLqcmqN2NkAnVXzgOH0//BWbE7Qd1KSmplJRURE1aGDexoE9Z4GJEra8vPTGv+WluXz5Mv/77rvv8hKd7du38zY1jz76KKWlKXdDXLhwIQ+OjI9mzZrZ8lVBEJdSsh1WRQUAAPqmiStBcXEx//vOO+/QsGHDqHv37vT111/zhpWsEbKSGTNmUEZGhukRHx/v4lyDK818qr2U1VI6/3rCQrUfgA6CGn9/f/Lw8KCkpCSz5ex5w4YNFd/DlpeX3vi3vDSsuolh7WiMfHx8qFWrVhQXF6f4uez1mjVrmj1Av9rUr17pyTABAEDCoMbb25uXkoSEhJiVorDnvXv3VnwPW14yPbN3715T+oCAAB68lEzD2r+wXlDGNOwzWZBy/vx5U5qCggKKjY2lFi3QrVJGv0x62Kp0ei+pSc3Ooz8u31Q7G9JB73YAnVQ/se7cX3zxBX377bd09uxZmjBhAuXk5PDeUMyoUaN41Y8R67W0a9cuWrx4MZ07d463i4mIiKCJEyfy11kV0uTJk2nBggW0bds2OnXqFF9H48aNTePQsFIW1uuJdSPfs2cPD27Y5zLDhw931LbQrWIdFle0b1S25G3ekPvLLNPSV69slcYLa/5wdFYAADTJ09Y3jBgxglJSUvhgeawhb9euXXnQYmzoy6qDWC8loz59+tCGDRt4A9+ZM2dSYGAg76rdseO90TinTZvGA6Px48dTeno69e3bl6+TDdZnxLpze3p68rFq7ty5w3tH7d+/nzcYhvL1//Qg7Zr0CHl7aqIJlRk3G+6J72/sV2aZQUMj2Oi5UGnT0Tga0bM56YXoAyUCyMrmoIZhpSzGkpbSDhw4UGYZK00pr0SFnSDmzZvHH5Z4eXnRxx9/zB9gm8spOXQ8Pp16BdTR9aZTus7oOVDQkn//eEpXQQ3clVdYRD6e9wYiBFCb9m7doVL03raEUfqKEnxtl7iefofeCj6hdjZAMHmFd3umAogCQQ3oSNkIplhDUc03obEkqvd2nFE7C0JB5ROAmBDUSEI7l3ZztjRdUCypcWhu5JVXgDtyKEtD9wwgCQQ1ktDqyceWfGv0K4IGoZ0wgJgQ1ICuu67L0JYIAADuQlAD+ql+snIZ2A7b0dyCoZ3Mnq//RxAOKwABIKgB3VAqlLmQmKVGVkDH5g/tSG0b1jBb9lAbf9XyAwD3IKiRhJYGoXPkd3xudZgqedEbWXv7FClUacq6LQC0AEGNLPQf08jxHcFlWHusgUt+K7sc+wAbA4SFoAZ0AxcbbFtHl9JcTM62ax32vh8AbIOgRkMy7hSonQWhoaMTiObTfRdI11AXB4JBUKMhH+8+T7Kx5ZwpQ7shAKHgJweCQVCjIXFptyv9XhnOPQptOgEcT4dFglFxt9TOAoBDIKiRhA7Pw2VgoD1w6PEkyeaMScqiZ1eFqp0NAIdAUAO6IctFCLRjx8kbJLpT1zPUzgKAwyCoAaH9Ep0oRVRzJ79I7SyAhKWbDOaxAj1BUANC+yY01uq0xRq+Ci3Yccau9++Kdm6JAKr2SmwLp25p/f4+AVwBQY2GZOcVVvq9MvQMUhr9VZclUgpe+y7KYXkBeX4zjJsd/bJ132UdNAdBjYYU2nHR1nAhhhQlNQBqQfUT6AmCGtANHy8PtbOgW24SXvkQIwNoD4Ia0Jz/ju2luPyRwHouz4ss0Kam5LZQcUcAQLkQ1GgJzqb0cKA/PWwhePFwd6M61bxdvlsAtEzGUjjQLwQ1WmLHyccgyQkYp2cAAHkhqNGQswmZamdBeLjpBEdBwSiA9iCo0ZD8omKSXUUlMdMGtnNRTuRy6zZmiK+ofVGzOlVIi1C6CXqCoAZ0ZVCnhqRFIl9YVv56kY7Hp5NsLI1T06xOVcXlwa/2IS1C6SboCYIaSYh80XTkCVgv31MkH+0+r3YWhPHqI63o8Xb1FV9r6OdLsg2+ByAaBDWS0EtDYb3C/hGPUi3TuEdaobcQgMAQ1ICmVHRPie6pAADyQlADuirR0GpBulbzLduxpsf9hDY1oCcIakBXtHqCRvUTqMXenwxGmwaRIKiRhEav9br9HqUvCmk5+WpnA6y4WKN6s6y9Z5Jw7IAwENRIQpaSAC325EjNRkCjFdo7upx/bkjOynNQTgDsh6AGhFVcbPudslarn0A8stwIAOgJghoQ1k/HrjtkPecTs0hkCMS0A/uqLAR/IBIENSCsyLhbZZZVpiDm1PUMh+QH5CLL3E8o3AQ9QVAjicV7zlOBBHNHKd1Jo3cGOOz40mEIIEnsBpJAUCOJk9cy6NvQWNI6FP+Dy+Bqb+V2woYCcSCokcjF5GzS+7lS6U5a9FOu/u795dWirvJkl3qWk1+kdhYATBDUaMQdB5w48nVR/YTeT6DeLN2WZu6W2Qe/nFM7CwAmCGo04qEP99u9DhlKiRVDHgm+txYs2oWLn2hu5eTTP9dHqZ0NAIdBUKMRGHG28m1qimWI5jRg1YFLlJyVS1qhdNjo7VD6EIEm6AyCGtAVpcH5pm8+RSJz1ND7Z29kkujyC/VQBaofSZnaCTIBrIGgBgRW9ra4kZ9vue+QudHtk0sPkei0VNJRpJBZg46+H4AeIagBTXlrQFu1swB20EpV4O38QuqxYJ/N7ytSmNoDAFwHQY1EtD4InX91b/Kr4lVuGoxjIzatHIK7ohMr9T6loGbrccdM9wEAFUNQA5q5AP7rL/e5rH2KK7lJFLiKnbt7LB1GFW1fD/eyb5y08TiJylG/F9GPO5BHpYKalStXUsuWLcnX15eCgoLoyJEj5aYPDg6mdu3a8fSdOnWinTt3lvlBzJkzhxo1akRVqlSh/v37U0xMjOK68vLyqGvXrvzHePy4uCcLEWnxgl9S7areamcBJKl+ssTLs/xT5soXu5GMzgk+aSzIw+agZtOmTTRlyhSaO3cuRUVFUZcuXWjgwIGUnJysmD40NJRGjhxJY8eOpWPHjtHQoUP5Izo62pRm0aJFtGzZMlq9ejWFh4dTtWrV+Dpzc8u2zJ82bRo1btzY1myDBu+mSmfX2pDszScCSUs0HmvapFgjbU6URqaeN+R+qulbfvVn12a16C0rShT1Br3aQLNBzSeffELjxo2jMWPGUIcOHXggUrVqVVq7dq1i+qVLl9KgQYNo6tSp1L59e5o/fz5169aNVqxYYbrQLlmyhGbNmkVDhgyhzp0707p16yghIYG2bNlitq5ffvmF9uzZQx9//HFlvy9oyI1S3U2tvfhLFCOUobG4VVhKx9qo3i3VyIom4LADTQY1+fn5FBkZyauHTCtwd+fPw8LCFN/DlpdMz7BSGGP6K1euUGJiolkaPz8/Xq1Vcp1JSUk8mPrvf//Lg6iKsGqqzMxMswdoy28XUkotkTlc0QeZSqUAQPCgJjU1lYqKiqhBgwZmy9lzFpgoYcvLS2/8W14aVprzyiuv0GuvvUY9evSwKq8LFy7kwZHx0axZMxu+KYgIF0TQAhynAOrRRO+n5cuXU1ZWFs2YMcPq97C0GRkZpkd8fLxT8wji0FpRuFL7DVm+u6i03qjeWnJ8S5CJTUGNv78/eXh48Kqgktjzhg0bKr6HLS8vvfFveWn279/Pq6J8fHzI09OT2rRpw5ezUpvRo0crfi5LW7NmTbMHaJteT8Cztt5rNC+a3AL7Z4cvCW1+AECYoMbb25u6d+9OISEhpmXFxcX8ee/evRXfw5aXTM/s3bvXlD4gIIAHLyXTsPYvrBeUMQ3rGXXixAnehZs9jF3CWU+s9957z5avACCUk9fS6ecTCSSq4auV28oBAIjI09Y3sO7crHSElZL06tWL91zKycnhvaGYUaNGUZMmTXibFmbSpEnUr18/Wrx4MQ0ePJg2btxIERERtGbNGlMx7+TJk2nBggUUGBjIg5zZs2fzbtus6zfTvHlzszxUr16d/23dujU1bdrU/q0AmuBuR5UAa5clYpVCROwtJ3Tbd9z3PHU9gxxJwF2gyJ5sinicAcjC5qBmxIgRlJKSwgfLYw152UB4u3btMjX0jYuL4z2ijPr06UMbNmzgXbZnzpzJAxfWVbtjx45mY8+wwGj8+PGUnp5Offv25etkg/UBGFl7rWhZt2zvOHatx7UGrGXPsaKl8aAc9ZvQ0ncGfbM5qGEmTpzIH0oOHDhQZtnw4cP5o7w7m3nz5vGHNdhoxvgR2U6W086Qrk1oyv9OSPndRaeVa589jbe18h0B9EgTvZ8AbLmrVJp/R5YgWI5v6Xwo1QPQJgQ1EpG5ph8XezFoJViwJ5s41gDUg6BGA5Kzys6BJePJFlUCoAVKhYKylBQCqA1BjQbEp91xyHoKiopJ0+xpvKn5kM46ol87Rc+fQxoKKxxrWvneAFqHoEYiO08l0o6TN0ir7KoSwEUFXHa0KRx/2PoALoGgRhMcd0p8fUOUw9YFYKuMOwUSdOlWWmbQdfD25aErDlkPgL0Q1GhASla+2lkQgj2Dmgl7TXEw0avZnlsdRucSM0m2hsJi7xX77Til3RJg0BcENRowY/NJtbMgBI10nIEK/O/oNeG3kV2jAitE0LIE1QBqQ1CjAbdua6PIXmSil2CA3ktqcPwBuAKCGpBk7ieSgizfU3ttauzKDgBYCUENSHGhQZ0/uOpYAwD1IKgBzbDnOvMVemeATcda5Y+2gmJtjAeVX1hMf1y+qXY2ABwKQQ0I6X9H4x3b+0nQNg2OLhG4dssxAzWS7KUgduTxhZ7NNVH99O7Ppyk7r1DtbAA4FIIaENK0H0869GJYLOBFxRn6f3KQrqeLH9iIzp64K8C/Gu14s6/wQfWG8Di1swDgcAhqQIoLTbpEPciirt5SOwuapzTTuy1q+noJX1IDoEcIakAz3G240PS7r57Z89TsPM2MZmsvXD8d39OuqreHTe8vXaqIfQLgGghqQDMC61e3Ou38IR3LLIu+nuHgHIFeJWflmj0Pnf64Te8v3f5L3GkSAPQFQQ1owi+THqZaVb3tWocW2qfKQAv74e1g8zZd9h57CGkAXANBDWhCIz9ftbOgGSgVsF+RnS3L3RS6T+tdsSyt8UFoCGpAE2ztzi1ibxPQhrQc+yeQLX249liwT/cX/Z9PJqidBQAENaANtnZGUWzC4CZeicoptPMRzuI9550yeN+ZG+LPTm6Pk9fQZg3Uh5Ia0OW8T8oxjVhRTXDkNdocdV3tbEAp5xOz7N4mSodroc5LatAWGkSAoAY0wdaB95TalYg2ku0mhVGTQb/sbacDABVDUAPCUQpI7JmhW7aLHO6Y7eOIQ01pFcU63zFoxwYiQFADwlG61ttcUkPiOx6fTiLbdsI5DT91GJ+WJcN3LEXvDaFBGxDUCE7GE4XSHa2t7WGUboolvM7Y5c3vjzllvfvPJQvdxVnCn5xDYLuBCBDUCC6/SNyTvyurZWzt/eSp8AZ7ZvnWEtGrAS6l5NCn+y6Qrsf5EXsXSHncgRwQ1AhOxsaFStcUW9vUtKhblWSlhaYbP0ReI1E54iengV3gcN/9gVm/QX0IagS30Qk9ZEIvpZLmqp9sLGRRKpWRpKAG7IQRmQG0C0GN4H6/6PgA5O9fhJP2ghr7IxLENGDd8SdHaRmAHiGoEZweuzJXRMIaNxCII7peo30JgDoQ1AhOwpgGxf+g8vEnxjoAwHYIagS390wSycZZJTWiBYjOyg8uqHZuP0HWAQC2Q1ADEvX4EiyqASE56yhBsAngfAhqQKLeJ2LdPyPEEpMjStBsHVcJABwDQQ0Ix1kFNaLdKcsyGKCMGtb0LbMMuxvA+RDUgHCcNfGfYDEN6HlCSzc3ejjQ3xHZAQAbIKgBaYIaGefREpXIZVS2zjOmlZJBABkgqAHhOOtiIMs1RgvfUwt5tFcjv7JVUADgXAhqQDinEzKcsl7RZoZ2VoPot4NP0NWbOU5Zt4y975rXqdw8YjOeak+i+m9YrNpZAHAKBDUgnNe+i3LIeta83N3s+ai1RyglK49ksGDHWRKZyNVPflW8zJ7/MKF3pdZTp5o3iWr21tNqZwHAKRDUgPDaNaxRqfcNuL9hmWUr9seQDHILitTOgmadvJZu+r9/dW+qXwPVSABagaAGhLfjzYcdtq48waqgnAWNVCunoKiYcvKLnLIdZdgnF5Oz1c4CSA5BDQjPQ6cjmaEzlnh2nrrhtHVP//Ek5RXquwQt+rpz2sMBWAtBDUhFlLvlgxdS1M4CKMjKLTR77sjD5XJqDq0LvYrtDuBECGpAKOm38zU5Bo6tRq894tT1iz56bXJWHqXlOHdfi9hD7Xr6HdKzrNwCtbMAkkNQA0L56dh1p67fXfSrvURe+jKcZCNKUO0s6FUFmgxqVq5cSS1btiRfX18KCgqiI0fKv+sMDg6mdu3a8fSdOnWinTt3lrkbmjNnDjVq1IiqVKlC/fv3p5iYe71UYmNjaezYsRQQEMBfb926Nc2dO5fy88W/07NHYZEcjVpL8vRwbpyNmEYcZ25kkmj+G2ZePaTvEARAf2y+gmzatImmTJnCg4qoqCjq0qULDRw4kJKTkxXTh4aG0siRI3lQcuzYMRo6dCh/REdHm9IsWrSIli1bRqtXr6bw8HCqVq0aX2dubi5//dy5c1RcXEyff/45nT59mj799FOedubMmaRnW48nkGycXY6CoAbKcz4py+y5owtWdF5QA6C9oOaTTz6hcePG0ZgxY6hDhw48uKhatSqtXbtWMf3SpUtp0KBBNHXqVGrfvj3Nnz+funXrRitWrDCV0ixZsoRmzZpFQ4YMoc6dO9O6desoISGBtmzZwtOw93/99dc0YMAAatWqFf3tb3+jt99+mzZv3kx69lbwCbWzoEOofqoI+02ygQoBAHQd1LDqnsjISF49ZFqBuzt/HhYWpvgetrxkeoaVwhjTX7lyhRITE83S+Pn58WotS+tkMjIyqE6dOhZfz8vLo8zMTLMH3HMlVcxh9J1dkqLT3uEOdSklm35D7yynTmUBAAIENampqVRUVEQNGjQwW86es8BECVteXnrjX1vWefHiRVq+fDm9+uqrFvO6cOFCHhwZH82aNbPyW8rhsY8PkIyjEaP6qWIYP+cehDQA2qK53k/Xr1/n1VHDhw/n1WCWzJgxg5fmGB/x8fEuzSdUDi6ojnEoJpUup1RudFcUTjiPAWESgDhBjb+/P3l4eFBSUpLZcva8YcOy8+wwbHl56Y1/rVkna2fz2GOPUZ8+fWjNmjXl5tXHx4dq1qxp9tASWefuKXZwVFP6Ai1Tl+7nP7dcfQtWQlENgH6DGm9vb+revTuFhISYlrFeSex5797KM9my5SXTM3v37jWlZ920WfBSMg1r/8J6QZVcJyuhefTRR/nns0bDrC2Pnl1OEbPNi9bG8Sh9ZyxPSEOUmq3vIQ+0GNO4SXAEynpDBmLwtPUNrDv36NGjqUePHtSrVy/ecyknJ4f3hmJGjRpFTZo04W1amEmTJlG/fv1o8eLFNHjwYNq4cSNFRESYSlrc3Nxo8uTJtGDBAgoMDORBzuzZs6lx48a863fJgKZFixb08ccfU0rKvSHmLZUQaZ1EBQpOrX4qvT52vKkNJ315yVD9lHmngHy9PNTOBkjK5qBmxIgRPKhgg+Wxhrxdu3alXbt2mRr6xsXFmZWisKqiDRs28C7bbFwZFriwrtodO3Y0pZk2bRoPjMaPH0/p6enUt29fvk42WJ+xZIc1DmaPpk2bmuUHvRP0xfHVT+JdRIau/J1EJsOF11XHz6JhnWnajydJJoVoGAdaCmqYiRMn8oeSAwfK9qphjXrZwxJ29zxv3jz+UPLKK6/wh0wEKFBQRc0qlTokLTIIuF3PJZoP8Abisje88/QQ4IBzMRF+YyAvfTdM0TAZ6t6VVPV2cFBT6qrkiYFqKiTrsadEwII+u2EMItAzBDWCkvVux+ENhUut79vQq5RxBzMJl7vNUP1ksvaVng6tihEhSHL2aNG3cvD7AvUgqBGUpDFNmaCmc1M/u9Y34H7zhuT5RcU0c/Mpu9YJcnilT0vq3bquUG3EtGDR7nNqZwEkhqAGhFJcamLyH17rY9f6pvzlvjLLDmIKALCCj5f9p8dW9apLVwKbnJmndhZAYghqQCiZueZF196e9h2iSl1LRewRBeqr7AjM5ekVYD4/3Xd/xFFSZi7pGX5doCYENYKS4Y6uNHay/8/PZ5z+ORLWCAjrq8NXSBRztp52yef8a9Nx0rOsUjcmAK6EoEZY8kU1s7dEu+RzZOoAFXox1eb3uLIga/525wex1mLtrVzh7I1M0rNrt+6onQWQGIIaEMaNDNcUy8s0/9PfvwxXOwvahRI9AM1BUCOojDvyzdvj6O7cFqkY04hU3QLmjlxJc8kmyS90TYkQgIwQ1Ahq2GfyzbCcnJWn+5IakapbLMnOK1Q7C7rmqmouABkhqAFhpLgsqHHJx2hSTl4hDV8tX0CtxFnlhgVF+q/XKkTgBipBUCMxNrKojN2bRZipW1RfHLqsdhZAB4Ijr6mdBZAUghoB/V6JHiuVnQMmM1f/VQ0t61Y1e46QxrJsCY4HcL4b6egBBepAUCOgwy4KamQpJn6hV3Oz5yipsQxj+NwjYykmgNYhqJFckQRXsdLfEW1qwJU+e7GbdO3WpBw9FISAoAaku+NmvazScuTrMm+NotKTb0nMUQU1DwX6kyiNwHu+t0/tbAA4FYIa0D2lwqgXMShdGRl3CujbsKsu2Sda4OHhmNIGUcos4m/ddt1npbnuswBKQlAjIFdW5YtU+fTXLo2dst4ChXZDagxVfzjGdW2lSvo2NNaqdFuPX3d6XrSkaW3zBuaVJUobLjcXhlc/HcOxBOpAUAPC8HLQnXFpT3VqRGq7nn6HXvpKnSkL5m5zzUSNetKpiR+90LOZQ9al1IbrlgrVn4LEVgBOhaBGQAYXlp+I0sEj43YBbY5yzt1d+0Y1SW3XNTDJn1rXvH+ujyS17YpONHs+f2hH8vJwd1oJiRrVn4hpQAYIagR09kaWyz7LlQFUeV77ruyF7cFWdUgvNNE9WKVb+Z2nElXvhVf6+HNkDzmlOc3O6HymbgC1IKgREBsUz1VEudaGXb5ZZtnql7qrkhdZpavYI8xlk5mqMD9YNR9PEgGqn0AGCGokJ9alxLzKqFZVb9ILUbdzSYv3XlDtswWLaVSd9NR5XPudMDEqqAFBjeQKCsUclwQD5LnWoRjXlQ6KXA1qpMuYxsXm/yz+jPSgPwhqJPfoxwcoOTOXRIOLimstC4khNYlWUgP22xQRj80ILoegRjDXXDhAltGGI3EkGlcU/7tq4tDU7Dx6Yc0fLvksrRq87BDlFhSpnQ0A0DgENYKZuOGY2lkQgitK/13VrXbVr5dIbfN+PiN0D6xLKTn0S/QNEoXAm6pS2L7feUqc7QvgLAhqBHM8Pl3tLIjBwSU1j7erT2rJL1K/BGLt71foQlI2iUzPE8Yffae/qp8fcjaZPlGxITiAqyCoAZcOn24tR+dowdCOJPtdP6p3rO+l4+7gM2O9Gj5llrlybJ5j8bdc9lkAakJQA0KMJlyao5vU+FXxIrUUFokR1aTdVh6H5kJSFh2NVf+iV6hSUc2gJb+p8rkjPg9T5XMB9AxBDajubysPO72kRq0B0PadSRKmF8iYr49S9PWMMsvfDj5BIpi++RSdTiibP2e7pjCFRTVv5x8vEVfVDyQB9AZBDajeffrqzdsumdn42W5NyizLd/I4PRMEmNeopPXhV82eZ9wpoJPXXB9IWLJo13m1s0ATHm1Nzeo4ZoZu2RXouaEUCAlBDdDKXy9SVm7ZKiA11anmhNGEFWqBHv3oV3KW2/mFVCBI1ZNoVS8imzqgrdpZ0A2lkkEAZ0JQA5RXWEzv7Tgr1JaYN+R+h69TKbxIyMh1Wlfn0WuPkGi+PxJPSSUGW7yRIdbAi64OAZUa67pjOGuHeWZVqONWBmAFBDUC2Xr8umqffTQ2jUTxzlPtqZFfFYev11Lw4qxeKCI0vlXy5vd3x0I6JVC1U8nJXA+6cEJXjN0CoC8IagQyaeNx1T7bIEEbH0uxS6ETgpoPfjlHojqdkMn//nVF2QbaInBlCZcru1V3aFRTlTYnt3LyaaUAA0ACuAKCGpCmzt3SBWzIit8denEpLjbQ6oOXhB6XZV1YrNrZEMK0H0667LM+f7l7mWWPfnTA6SM9z9oa7dT1A4gEQQ2oWlTz9HKF7txOKqrpd189xeXnk7LocIxj5oFiF6ihq34n0c3ZeppE9sflmy75nHwX9s5R6lF1Pf2OxRJERxGxmhHAWRDUAHc5Ncflg59ZukN1Vg/zYd2bWnztxLV0u++Y2ftXH7wsVBdprWITgLJqE+OEoF8dvkKfHbhEf//iD2ozcyf93zdHKdLOcV5KNphWk7OroIpVHtKa9QIEcBV1RiQDIc3YfIo+Gt7Fpb2uXNmmxsPdjeY83YHmbT9T5rUl+2KoY2M/6t+hQaWrtt7feZZffMExHpi/1+Jr+88l88fCZzvRCz2bVap0T+k4eOaBsmMZab1dj9rTdCzdF0MznmqvbiZAGiipEcT0H11Xt29JcOQ1l37enXzliR77tvF32mc+0LyWxdf+sS6CkrNyKzXNA6tGQ0CjTiD+yEe/Unxa2QEcK3IlJafMssUuDOpN+Ugtmw9HYaWHN3PySE2f/3ZZ1c8HuSCoEcTGo2IMpe/KE7nSnfjG8Q9SYIMaqhXF93ovhJKtrJZgF4xDMSnUZd4eOnvjbo8icL34tDv08KJfae+ZJKvfs/FIHJ1R2GdqjFHDAmLWtsYZPtp9nnIL1B/V11nfD6A0BDUCcHbvBxFN+E55+oD7nBjQMNaU9Pd6P6TCXlnnE7N4l+iXvxJvgD1ZjVsXwSfntOb3xuaZUsO/B7VTXB7lpHmgVh0Qoxdedi7a1YBrIKgRwCWFYnC1fOGComLWwPNcYpbFdi/O5Gnl+tnd8/DVoXyupLzCIrMLIptWYuCS3yj6OkpnRDPg098qvEmYvOm4qvNKKdH7bc2RK67pzQaAoEYAs7aoc9eo5L2dzp0uIfJqGg37zPLQ6V4ezg1qujS13KZGaUTgd36KprazdtHiPed5uw120WRF+iCuvh/+arHHzfaTCbT1eAKJOMpzZdoFlceVIzNXZPbW08LNLwf6hKBGZawb9R+XxZmioLwGvPbadiKBhn0WVm4abw/nHpKszcSVhU/Z/L7l+y/ydhsxydkuHUNHJG3qV3fIgHOuaL/RYc5us3Ycl1KyeZXnxA13p4gQ0Qe7zul67rGvf8eAj+B8CGpUtut0YoVp5g/tSJvGP0hDuzZ2zuzVpQxYctCh62M9ihb+ctY051B5PJ0c1DhzcD97fP1KT6vSNfbzJbW0qWd7UNOuYQ3yr+78Y1bJQx/sp/azd1HL6TvoicUH6Zfo8n9rrJG6K9TwUR5JY8fJG3Q5xTFBs6vG4Hnpweb04bBO1LCmL7WsW3ZwwZI+2XvB5WNhgXwQ1KgoJimr3DtHVqIQ+8FgevnBFhTUqi4teeEB+v3fj1e43uYKI5fa2pvE3rY1rF0D6xE04vMw3qPo84MVr2/7G31JVkq9bn76Z58yy6YNakePtXVNqc7wUoMVDuzYQHNB5J0C60sdH2xVl1xh3dheFl/7y6e/2T0YHxv3Juj9EJve83Cg7cMovBjUnBYM7UQjejanP2Y+QQemPkaHpj1W7nv+79sImz8HwOlBzcqVK6lly5bk6+tLQUFBdORI+cWcwcHB1K5dO56+U6dOtHPnzjIXwDlz5lCjRo2oSpUq1L9/f4qJiTFLk5aWRi+++CLVrFmTatWqRWPHjqXsbOdUBbgCKxpnJzBLLix4UvFiUMXbg584nmhX3+J714zqTu0VJs9rWruKTW1rBi35jc9jZC22H1m7gM8PXqKAGTvpyaWHKPyK9VVrHZv4kas838Py6MKieKB5bcXlDRVKa36Z9DA1clApzt+DmtOJOQPog2Gd6buxQTSgQwNa/VJ3Gtq14oHpWpS6W69moVRCZmy//q1LY4sBCbsRqEwVMPv9sWq21jPNz68V+XliX/rv2CCbP69twxqKU0GwGzFLpVFsFvafT4jXpgn0w+YzzqZNm2jKlCm0evVqHtAsWbKEBg4cSOfPn6f69cteaENDQ2nkyJG0cOFCevrpp2nDhg00dOhQioqKoo4dO/I0ixYtomXLltG3335LAQEBNHv2bL7OM2fO8ECIYQHNjRs3aO/evVRQUEBjxoyh8ePH8/VpCZtM8NvQ2HIbm15+/6lyx8tgJ46vXunJi9WVtGtYk1/kXv4qnA79OafRhn8EUZ82/hbfo4T1UGpV4gQZ4F+N55+deDPvFDh0duvOTV0X0DCLnutC/4tw7WCDjhDUqg49cl89+v6I+bhGNat40Y8T+tDcbafp+q07NOUv91G9Gj40ZKVt81AN7tyI3n+mk+l530B//rDWG48H0tvBJ/j/P3quM/lX96GZT7Wjf226u4xhF7ysPLG6+H41uodLP691OVV5UXHp1H7OLv7/57o35bN7N65VhWpX9aKAetXI092d/Kp4UfrtfIq9mUPH4zNo7eErNo8Fw0r8RvVpSZ3+/O2xkpf14XFWvfehNnVpZK/mFl8/9Z+BvEOA0lQWb3x/jD9YMNWxSU0hq4NBu9wMNg6SwgKZnj170ooVK/jz4uJiatasGb3xxhs0ffr0MulHjBhBOTk5tH37dtOyBx98kLp27coDI/bxjRs3prfeeovefvtt/npGRgY1aNCAvvnmG3rhhRfo7Nmz1KFDBzp69Cj16HH35LNr1y566qmn6Nq1a/z9peXl5fGHUWZmJs8nWzcr7XGUi8lZ9N0fcXxQN3aRLyr6829x8Z9/7z5nJR6HLqZSvoWpAYyi/zOQqlt5d2spQGF3SgwrNfnn+ij6x8MBNOTPu+ydp27wZaKJee9J8nJBe5qSvj8Sx0ekdTXWDmHP6SRKzrp7fLKSkDWjepjtz7cH3EcTHw80WxY6/XF+cWPY3S67MDCzBrenfzzcSvGzBnx6kC4kWV+iaTx2LGG9h8qrMl36QleatPF4mXWxtiJsHJkmtavS2tE9+Hfv88F+EsHykQ/QXy2UnDgLuzl4culvvKrX1c7MG0hVvcueY1i1F/s9/FBqZHE2FUXp30nErP48YLVmkMOKxgRigU2TWlWoipcHPwd4eriRu5sbny6F/xW0HRwoa12/Om8y4Ujs+u3n52fV9dumoCY/P5+qVq1KP/zwAy9tMRo9ejSlp6fT1q1by7ynefPmvGRn8uTJpmVz586lLVu20IkTJ+jy5cvUunVrOnbsGA90jPr168efL126lNauXcuDnlu37kX9hYWFvBSHVW0988wzZT733Xffpf/85z9lljs6qGHdJh3Ry2DX5Id5CYutAdXF5GxexL/xSDztOHWD2jaoQbv/9Ui572PdXX85lcjbG8zaEk1qM5YiqeFcYiYNWnLIZZ/3fw8F0Jy/duDz4Xy67wJ1bVaLt51hJ202geOi3efoue7NqHuLu1VPxqCGlZL9+vajpvXkFhTx6sFuzWvTJyPu/W5KYz9vVhVoDVbFNKhjQ6vSsiA94uot+tem4zSyVzNecvT2wPtocKfG9NKX4dQzoDZNHag80JwRazTKBjpM+3PiSrUCmqc7N1Llosn2YY8F+3iA46pgmg3+V8PXq8L9ejM7jyasj+IdE46+058emLeHv77jzYf5uE1t6ls/SGZOXiEfqPKyQONxgfOwkuR1/2e53ZhQQU1CQgI1adKEVyn17t3btHzatGl08OBBCg8PL/Meb29vXq3EqqCMVq1axQOOpKQkvq6HHnqIr5u1qTF6/vnn+YmGVXe9//77fB2siqskVt3F1jNhwgTVSmpiU3MoODKePNzd+cBubPC40n/Za5m5BXTuRialZOdRYP0avGqAXUDYHYqvl4fd+WC7MTU7n6/X1vex4eK/PHSFlxDdyMilaj4e1KmJH6/mYqd6VtLE7phYaRSr62cTUeYX3v3LLkisqiz88k0+gBj7zg1q+lJKVh6/82PVIuykxorLu7WoTUkZuTyS9/F050Xw7KKsxtD0ShdYdjJnjbfj0m7TY23r8+nC2R1t3Wre/HtExd3i26WmrxffDll5BdS0dlX+/VhgyUYZZl+F3W1+efgytaxbjRr5VeFVOsYeKfVr+PDjmn3ekdg0HtQo3TWXDEDDLt3kDcVLl+CxfWfNxZiND8K60wfUrUbN61bleWb7bfn+GP5/lid2ImL7SA1s7qz8omJ+3Pl6etDtgiI6n5jJvxvbtqykiQV8tap68W1/8no6H6GWbees3EL+nVh7oojYW5SYmct7XLGA/VhcOh/3qG41H74P2PvZhZW162hVrxq90ieAOjR23Lmgsk7Ep9Phi6l0OCaVz9PESnPZd2O7lv2WWMN/9pz9ButU8yIfTw+qW92b7y/2nWLY9rldQAYy8GOVnU/Yg20jd3eitg1qUuv61fj7bMFG1WbnAPY5xhJmb0/7SlPZfmQl1heTsuhmTj5l5hbycwFr71fV24OfawqLDPxcw65Mhj+Pc2uvUmwb3C3bAbW09K/Gq03VCmp024rPx8eHP1yxAyu6I3UFdgGwNaAxvu/+xn70aTl3+zJgXclZ75fyesDwQOdP7OTuV9XLrDFsyYaTrFdIaewCVfLz+rSuuHSKBTxPtFfudWRt6QK7M38xyLw4mN2Bz/3r/SQC43Y0YsFb9xZ1TM8fDjTv7fV8z2aK6xnzUIBVn/fmE4Ekki7NavHH64+1qdT7S24rRyrZcN/eYMaodjVvi42kARzBpiPV39+fPDw8eAlLSex5w4bKxdZseXnpjX8rSpOcnGz2Oqt+Yj2iLH0uAAAAyMWmoIZVJXXv3p1CQu6NgcAaCrPnJaujSmLLS6ZnWA8mY3rW24kFJiXTsKImVpVlTMP+sjY7kZH3JkHcv38//2zWcBkAAADA5uon1uiXNQxmvZB69erFu3Sz3k2sizUzatQo3u6GdeFmJk2axBv9Ll68mAYPHkwbN26kiIgIWrNmjakInTUiXrBgAQUGBpq6dLMeTcbGyO3bt6dBgwbRuHHjeI8p1qV74sSJvGeUUs8nAAAAkI/NQQ3rop2SksIHy0tMTOQ9lFj3atYFm4mLiyN31jrtT3369OFjycyaNYtmzpzJAxfW88k4Ro2xoTELjNi4M6xEpm/fvnydxjFqmPXr1/NA5oknnuDrHzZsGB/bBgAAAKBS49RolS2tpwEAAEB712/M/QQAAAC6gKAGAAAAdAFBDQAAAOgCghoAAADQBQQ1AAAAoAsIagAAAEAXENQAAACALiCoAQAAAF3Q7SzdpRnHGGSD+AAAAIA2GK/b1owVLE1Qk5WVxf82a9ZM7awAAABAJa7jbGTh8kgzTQKb0TshIYFq1KjBJ9GsKCpkwU98fDymVBAY9pN2YF9pA/aTdsi0rwwGAw9o2ATWJeeWlLqkhm2Ipk2b2vQedqDo/WDRA+wn7cC+0gbsJ+2QZV/5VVBCY4SGwgAAAKALCGoAAABAFxDUKPDx8aG5c+fyvyAu7CftwL7SBuwn7cC+kryhMAAAAOgbSmoAAABAFxDUAAAAgC4gqAEAAABdQFADAAAAuoCgBgAAAHRB6qAmNjaWxo4dSwEBAVSlShVq3bo178qdn59vlu7kyZP08MMPk6+vLx+WetGiRWXWFRwcTO3ateNpOnXqRDt37nThN5HTypUrqWXLlnybBwUF0ZEjR9TOklQWLlxIPXv25FOP1K9fn4YOHUrnz583S5Obm0uvv/461a1bl6pXr07Dhg2jpKQkszRxcXE0ePBgqlq1Kl/P1KlTqbCw0MXfRh4ffPABnypm8uTJpmXYT2K4fv06vfTSS/z3wq5J7FoSERFhep11Vp4zZw41atSIv96/f3+KiYkxW0daWhq9+OKLfJThWrVq8WtcdnY2ScMgsV9++cXwyiuvGHbv3m24dOmSYevWrYb69esb3nrrLVOajIwMQ4MGDQwvvviiITo62vD9998bqlSpYvj8889NaX7//XeDh4eHYdGiRYYzZ84YZs2aZfDy8jKcOnVKpW+mfxs3bjR4e3sb1q5dazh9+rRh3Lhxhlq1ahmSkpLUzpo0Bg4caPj666/57+L48eOGp556ytC8eXNDdna2Kc1rr71maNasmSEkJMQQERFhePDBBw19+vQxvV5YWGjo2LGjoX///oZjx44Zdu7cafD39zfMmDFDpW+lb0eOHDG0bNnS0LlzZ8OkSZNMy7Gf1JeWlmZo0aIFvyaFh4cbLl++zK9NFy9eNKX54IMPDH5+foYtW7YYTpw4Yfjb3/5mCAgIMNy5c8eUZtCgQYYuXboY/vjjD8OhQ4cMbdq0MYwcOdIgC6mDGiUsMGEHidGqVasMtWvXNuTl5ZmW/fvf/za0bdvW9Pz55583DB482Gw9QUFBhldffdVFuZZPr169DK+//rrpeVFRkaFx48aGhQsXqpovmSUnJ7MxrwwHDx7kz9PT03lwHxwcbEpz9uxZniYsLIw/Z0GMu7u7ITEx0ZTms88+M9SsWdPsNwf2y8rKMgQGBhr27t1r6NevnymowX4SA7uu9O3b1+LrxcXFhoYNGxo++ugj0zK273x8fPjNNnPmzBn++zp69KjZzbubm5vh+vXrBhlIXf2kJCMjg+rUqWN6HhYWRo888gh5e3ublg0cOJAXs9+6dcuUhhUDlsTSsOXgeKx6MDIy0mybswlL2XNsc3V/O4zx98P2UUFBgdl+YlW0zZs3N+0n9pcVsTdo0MDst8NmID59+rTLv4OesWpAVs1X+lyF/SSGbdu2UY8ePWj48OG8GvaBBx6gL774wvT6lStXKDEx0Wz/sUkeWdV7yd9TrVq1+HqMWHp2fgwPDycZIKgp4eLFi7R8+XJ69dVXTcvYQVTyhMsYn7PXyktjfB0cKzU1lYqKirDNBVJcXMzbaDz00EPUsWNHvowd/+xmgJ1kLf02rPl9gf02btxIUVFRvB1UadhPYrh8+TJ99tlnFBgYSLt376YJEybQm2++Sd9++63Z76G8a01iYiIPiEry9PTkNxqy/J50GdRMnz6dN4Qr73Hu3LkyDbQGDRrEo+Rx48aplncArZYCREdH84sniCU+Pp4mTZpE69ev543qQdwbg27dutH777/PS2nGjx/Pr0WrV69WO2ua4kk69NZbb9Err7xSbppWrVqZ/p+QkECPPfYY9enTh9asWWOWrmHDhmV6axifs9fKS2N8HRzL39+fPDw8sM0FMXHiRNq+fTv99ttv1LRpU9NydvyzqsL09HSz0pqSvw32t3SvtdK/L7APq15KTk7mF0wjVtLJ9teKFSt4qQD2k/pYj6YOHTqYLWvfvj39+OOPZr8H9vtgaY3Y865du5rSJCcnm62D9SRkPaJk+T3psqSmXr16vO6+vIexjQwroXn00Uepe/fu9PXXX/O6x5J69+7Nf/ysbYDR3r17qW3btlS7dm1TmpCQELP3sTRsOTge23dsf5Xc5uwuhz3HNncd1tGABTQ//fQT7d+/nw+NUBLbR15eXmb7ibVFY124jfuJ/T116pTZiZj9dlh31NIneKicJ554gm/j48ePmx6szQXr9mv8P/aT+ljVbekhES5cuEAtWrTg/2e/LxaYlPw9sbZnrK1Myd9Teno6D2SN2G+TnR9Z2xspGCR27do13t3tiSee4P+/ceOG6VGydTnr0v3yyy/zrqusK3HVqlXLdOn29PQ0fPzxx7x3x9y5c9Gl28nYfmCt/r/55hve4n/8+PG8S3fJXjTgXBMmTODdSw8cOGD227l9+7ZZV2HWzXv//v28S3fv3r35o3SX7gEDBvBu4bt27TLUq1cPXbqdrGTvJ+wncbrbs+vIe++9Z4iJiTGsX7+eX2u+++47sy7d7DzHhh85efKkYciQIYpduh944AHeLfzw4cO8xxu6dEuCjbHB4jqlR0lsPADW1Y5dRJs0acIPrNL+97//Ge677z4+dsr9999v2LFjhwu/iZyWL1/OL5hsm7Mu3mxcBnAdS78d9rsyYifbf/7zn3xYBHaCfuaZZ8xuGpjY2FjDk08+ycd/YmPUsHGiCgoKsCtdGNRgP4nh559/5kE+u9a0a9fOsGbNmjLdumfPns1vtFkadkN+/vx5szQ3b97kQUz16tX50Ahjxozh3fll4cb+Ubu0CAAAAMBeumxTAwAAAPJBUAMAAAC6gKAGAAAAdAFBDQAAAOgCghoAAADQBQQ1AAAAoAsIagAAAEAXENQAAACALiCoAQAAAF1AUAMAAAC6gKAGAAAASA/+HwqfPJrAF0dfAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(xvec, np.abs(psi4)**2)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 3, "id": "2b01c0dc-3c4b-4d09-97fc-4d41c120a22c", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/var/folders/d7/6t4hylss1p9_9vftvfm4fbm80000gn/T/ipykernel_41528/630846714.py:11: DeprecationWarning: `trapz` is deprecated. Use `trapezoid` instead, or one of the numerical integration functions in `scipy.integrate`.\n", " cpops[i] = np.trapz(pops[i0:i1], shifted_ckvec[i0:i1])\n" ] }, { "data": { "text/plain": [ "(np.float64(0.3613023794837458),\n", " np.float64(0.12050147253410211),\n", " np.float64(0.3599777025771895),\n", " np.float64(0.12058717738992883))" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "phi = xspace_to_pspace(psi4)\n", "\n", "shifted_ckvec = np.fft.fftshift(ckvec)\n", "cpops = np.zeros_like(kvec, dtype=float)\n", "pops = np.fft.fftshift(np.abs(phi)**2)\n", "\n", "for i in range(len(kvec)):\n", " i0 = np.argmin(np.abs(shifted_ckvec - (kvec[i] - 1)))\n", " i1 = np.argmin(np.abs(shifted_ckvec - (kvec[i] + 1)))\n", "\n", " cpops[i] = np.trapz(pops[i0:i1], shifted_ckvec[i0:i1])\n", "\n", "cpops /= np.sum(cpops)\n", "\n", "state2idx = lambda nstate: np.argmin(np.abs(kvec - 2*nstate))\n", "arridx = np.array([state2idx(2*n_bragg), state2idx(n_bragg), state2idx(0), state2idx(-n_bragg)])\n", "\n", "pA, pB, pC, pD = cpops[arridx]\n", "pA, pB, pC, pD" ] }, { "cell_type": "markdown", "id": "98bd0174", "metadata": {}, "source": [ "I will copy these values to a cell below for comparison with the planewave method!" ] }, { "cell_type": "code", "execution_count": 4, "id": "b1d8b214-82b5-48ca-a1d6-8dd611f67dbe", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGdCAYAAAD60sxaAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAV1VJREFUeJzt3Ql4VNX5P/DvzGSyb4R9X0RZZBNUwLqhFKTUaqGtVv+KirZa9FeliqU/ihbbH1atSxXFuoCtWou27ohsQkUWFUXZZZMthD37OpP7f865c+7cmbl3kmCSuZP5fp4nzHYy3JtJJm/e877nuDRN00BERESUQNyxPgAiIiKi5sYAiIiIiBIOAyAiIiJKOAyAiIiIKOEwACIiIqKEwwCIiIiIEg4DICIiIko4DICIiIgo4SShhaqtrUV+fj6ysrLgcrlifThERERUD2J95pKSEnTq1Alud9PlaVpsACSCn65du8b6MIiIiOgU7N+/H126dEFTabEBkMj8qC9gdnZ2rA+HiIiI6qG4uFgmMNTv8abSYgMgNe0lgh8GQERERPHF1cTlKyyCJiIiooTDAIiIiIgSDgMgIiIiSjgMgIiIiCjhMAAiIiKihMMAiIiIiBIOAyAiIiJKOAyAiIiIKOEwACIiIqKEwwCIiIiIEg4DICIiIko4DICIiIgo4TAAIqL48+UrwL61sT4KIkqUAOj++++Xu7OaP/r27Ws8XllZiSlTpqB169bIzMzExIkTcfjw4ZDn2LdvH8aPH4/09HS0a9cO99xzD3w+X8iYFStWYOjQoUhJSUHv3r0xf/7873qeRNRS7FoOvP0r4MWxgL8m1kdDRImSATrzzDNx6NAh42PVqlXGY3fddRfeffddvP7661i5ciXy8/MxYcIE43G/3y+Dn+rqaqxevRovvfSSDG5mzpxpjNmzZ48cM2rUKGzYsAF33nknbr75Znz44YeNcb5EFO+2vR+8XlMRyyMhojjm0jRNa0gG6K233pKBSbiioiK0bdsWr776Kn7yk5/I+7Zt24Z+/fphzZo1GDFiBD744AP88Ic/lIFR+/bt5Zi5c+fi3nvvxdGjR5GcnCyvv//++9i0aZPx3FdffTUKCwuxaNGiep9YcXExcnJy5HFlZ2fX+/OIyMHE29VjA4DiA8DPXwP6jIv1ERFRI2uu398NzgDt2LEDnTp1Qq9evXDttdfKKS1h/fr1qKmpwejRo42xYnqsW7duMgASxOXAgQON4EcYO3asPNnNmzcbY8zPocao57BTVVUln8f8QUQtzKENevDjzQB6XRzroyGiONagAGj48OFyykpkYp555hk5XXXBBRegpKQEBQUFMoOTm5sb8jki2BGPCeLSHPyox9Vj0caIgKaiwj7dPXv2bBkxqo+uXbs25NSIKB7sXKZf9r4EKD+uF0JXnIz1URFRHEpqyOBx44Lp5kGDBsmAqHv37liwYAHS0tIQS9OnT8fUqVON2yJgYhBE1MIU6hlntDsT+OfVQMFG4JrXgTPGxPrIiCiR2uBFtueMM87Azp070aFDB1ncLGp1zEQXmHhMEJfhXWHqdl1jxDxgtCBLdIyJMeYPImphcroCnc4C2pwOZHfW7yvJj/VREVGiBUClpaXYtWsXOnbsiGHDhsHr9WLZskCKGsD27dtljdDIkSPlbXG5ceNGHDlyxBizZMkSGaz079/fGGN+DjVGPQcRJbCL7gF+sQIY+BMgq6N+XzEDICJq4gDo7rvvlu3t3377rWxj//GPfwyPx4Of//znsu5m8uTJchrqo48+kkXRN954owxcRAeYMGbMGBnoXHfddfjqq69ka/uMGTPk2kEigyPceuut2L17N6ZNmya7yJ5++mk5xSZa7ImIDCoDxACIiJq6BujAgQMy2Dl+/LhseT///POxdu1aeV147LHH4Ha75QKIoitLdG+JAEYRwdJ7772H2267TQZGGRkZmDRpEmbNmmWM6dmzp2yDFwHPE088gS5duuD555+Xz0VECUyt2OFy6ZfZzAARUTOtAxRPuA4QUQtzYjfw9EigdW/gtk/0FaH/8WOgbT9gCrfFIGopip26DhARUUyIqS5fZXD156xO+iWLoImoqafAiIhipviQfpkdCHxyugDnT9Vvi0S2mhojIqoHBkBEFB+KD4YGQCmZwOj7YnpIRBS/OAVGRPGhJCwDRET0HTAAIqL4ygCp2h+hpEDfDuPk3pgdFhHFJwZARBQfRLBjbn8Xlv4BeHEssOmNmB0WEcUn1gARUXxo1y9Y/KxktNYvuSEqETUQAyAiig8/ejLyvtRc/ZIBEBE1EKfAiCh+pakAKHQTZiKiujAAIqL4ldZKv2QAREQNxACIiJxPdHn9qRPw5DDrKbBKZoCIqGEYABGR84kan5oyoKo09H5OgRHRKWIRNBE5n8rwqCkvJbc7cNFvgcx2MTksIopfDICIyPlUjY/K+CgZbYBR02NySEQU3zgFRkTxkwFSNT9ERN8RAyAicj61zk/4FJhwbCewbx1QVdLsh0VE8YsBEBHF7xSY8PIE4MUxwJGtzX5YRBS/GAARkfOJHeC7nAPk9Yp8zOgEC2SJiIjqgUXQROR8w3+pf1jhYohEdAqYASKi+MbFEInoFDAAIqL4xikwIjoFDICIyPnmjAAePRMo2BRlR3huh0FE9ccaICJyvuKDQFUxkJQSpQaIRdBEVH8MgIjI2fw+PfixWweo20jgonuBjkOa/dCIKH4xACIiZ6ssCl5PzYl8vNtw/YOIqAFYA0REzqayP950wOON9dEQUQvBDBAROVt1qX6ZkmX9uK8KOLEb8FcDHQc366ERUfxiBoiInE3t8ZWcaf34yW+Bp0cAL/2oWQ+LiOIbM0BE5GyeFKDz2UBOF+vHVWAkMkWaBrhczXp4RBSfGAARkbN1GQbcssz+8ZRAAFTr06fDvKnNdmhEFL84BUZE8c08NabqhYiI6sAAiIjim9sDeDNCO8aIiOrAAIiInG3V4/o2GCserHsarIoZICKqHwZARORsZUeB4gNAdZn9GNUizykwIqonFkETUXy0wdutAyQMu1Efl92p2Q6LiOIbAyAiiv8A6Lzbm+1wiKhl4BQYETmbmtayWwiRiOgUMANERM6mCptVobOV8hNA6RF9t/is9s12aEQUv5gBIqL4nwJb/gDw9HDg8xeb7bCIKL4xA0REzpbXQ9/oNL21/RgVHKlgiYioDgyAiMjZrnq57jHJqg2eARAR1Q+nwIgo/nEhRCJqIAZARBT/OAVGRA3EAIiInKviJPDYAGDOCKDWbz9OtchzJWgiqifWABGRc1UWA0X7gaRUfdNTO8wAEVEDMQAiovhfBDGvJzDiV0BO12Y5LCKKfwyAiCgOFkGMsgaQkNcLuGx2sxwSEbUMrAEiIucKtLVr0VaBJiI6BQyAiMixNu0+KC93F7uiD9Q0oDgfOPpN9GJpIqIABkBE5FgrNu6Wl3tL6nirEkHPo/2AOecAlUXNc3BEFNcYABGRY/k8qdhV2xEHtLbRB3qSAE+Kfp2t8ERUDyyCJiLHWp91CR4vGCSvX1/X4OR0oKIKqC5vjkMjojjHDBAROZbX04C3KG+GfllT1mTHQ0QtBwMgInIsr6eO4ufwDJDADBAR1QMDICJyrIknX8Si5HvxE8/Kugd7AwFQDafAiKiJA6AHH3wQLpcLd955p3FfZWUlpkyZgtatWyMzMxMTJ07E4cOHQz5v3759GD9+PNLT09GuXTvcc8898Pl8IWNWrFiBoUOHIiUlBb1798b8+fO/y6ESURxq6ytAX/d+ZKEeQU1yYAqsmlNgRNSEAdBnn32GZ599FoMG6QWKyl133YV3330Xr7/+OlauXIn8/HxMmDDBeNzv98vgp7q6GqtXr8ZLL70kg5uZM2caY/bs2SPHjBo1Chs2bJAB1s0334wPP/zwVA+XiOJQqlYhL8uRWufYd2uGYVn2j6G16tkMR0ZECRkAlZaW4tprr8Vzzz2HVq1aGfcXFRXhhRdewKOPPopLLrkEw4YNw7x582Sgs3btWjlm8eLF2LJlC15++WUMGTIE48aNwwMPPIA5c+bIoEiYO3cuevbsib/85S/o168fbr/9dvzkJz/BY4891ljnTURxIEWrlJflWgo0sdihjcoaP+7YPQKTj/wU+1P7NOMRElFCBUBiiktkaEaPHh1y//r161FTUxNyf9++fdGtWzesWbNG3haXAwcORPv27Y0xY8eORXFxMTZv3myMCX9uMUY9BxElhpRalQFKQZWv1nbc8TL9jychOYmljUTUBOsAvfbaa/jiiy/kFFi4goICJCcnIzc3N+R+EeyIx9QYc/CjHlePRRsjgqSKigqkpaVF/N9VVVXyQxFjiaiFZICQKrM8qV6P5bgTpdVIRRWyUQ6t4gSQ06mZj5SI4k2D/lTav38/fv3rX+OVV15Bamrdc/LNafbs2cjJyTE+unbtGutDIqLvyOsPZIC0FFTWRMsAVeHOpH/j09QpyFz3OL/uRNS4AZCY4jpy5IjszkpKSpIfotD5r3/9q7wusjSijqewsDDk80QXWIcOHeR1cRneFaZu1zUmOzvbMvsjTJ8+XdYgqQ8RrBFRfCv1tkaB1gqlSJMZIDsnyqpRoQW2wmAbPBE1dgB06aWXYuPGjbIzS32cffbZsiBaXfd6vVi2bJnxOdu3b5dt7yNHjpS3xaV4DhFIKUuWLJHBTf/+/Y0x5udQY9RzWBHt8uI5zB9EFN9eOvMFjKiag11aZ1T6ogdAok5IcDEAIqLGrgHKysrCgAEDQu7LyMiQa/6o+ydPnoypU6ciLy9PBiF33HGHDFxGjBghHx8zZowMdK677jo89NBDst5nxowZsrBaBDHCrbfeiqeeegrTpk3DTTfdhOXLl2PBggV4//33G3K4RBTn/LXBzq+Kan/UIuiKQADEDBARxWQzVNGq7na75QKIoihZdG89/fTTxuMejwfvvfcebrvtNhkYiQBq0qRJmDVrljFGtMCLYEesKfTEE0+gS5cueP755+VzEVHiqDW1vkerARJF0DWBKTA3M0BEVA8uLdriGnFMdIGJYmhRD8TpMKI4VHoUBU9dhr3lSbiqeibm3XgORvVpZzn05pc+R9L2dzE3+XGUdzgH6bcubfbDJaL4+v3d6BkgIqJGUVWMDpU7ke7SGx+qohZBVyErMAXGDBAR1QdXDCMiZ6ouDdkGI+oUWFk1DmhtsMB3EYp6XNZsh0hE8YsBEBE5U3W5sQaQEK0NXhRBi06xab5f4shZv262QySi+MUAiIicKbCru8oAVdgEQNW+WpRU+ozb/pZZ1khEjYwBEBE5U40KgFKiToGdLFf7gGlIRyXc5UcABkFEVAcGQETk7AyQpmqArDNAx0v1AEjsA7Yl9SYM+uc5gD+4OSoRkRUGQETkTC4Pijx5OI4sedNuJWhRAG3OFJmDJyIiO2yDJyJnGnwVfrf5dLz/9SF5s8pmCkxlhnxIQrXmQbLLD9Tom6gSEdlhBoiIHKu2HlthmFeL5nYYRFRfDICIKD62wrCZAjOPUR1jnAIjorowACIiZ/rvI/jNgV/jCveqqEXQpiSRsWYQN0QlorowACIiZzq6HWdUbUJbV1HUNnjLKbDAIopERHZYBE1EzhTY1V0FNfXJAC2rPQt5vc5Cp0zrTVOJiBQGQETk8HWAAgGQzyYDZIqAHvP9FP3PPRudOrZvpoMkonjFKTAicqZAK7uxEnQ9usAEvzklRERkgwEQETl6K4xKFQDZdoGZb2lwiakzrgNERHVgAEREzs4A1bEbvDkD9MekFzH27bOAT55opoMkonjFAIiInCkpTWZ/KpAcvQus1qoLjFthEFF0DICIyJluW4Wft3sTm7Re9V8HSAVAgQ4yIiI7DICIyLHMwU21v+51gNTO8VwHiIjqwgCIiBzLPL0l4hwtrONLv9+8FYaaAittngMkorjFAIiInKeqFJj3A9xXNBNJ8Bl3W8Q/IW3v3AyViOqLCyESkfOIIua9n2AoXPDBEzLd5Yar7r3AuBUGEdWBGSAicuwaQFVySisY8FitcWiuATqgtcW+9pcCPS9onuMkorjFDBAROXYNoEpXIKNjs+qzYL7ra+00rBzyI1w3skfTHyMRxTVmgIjIeQJTWGoV6GgBkD/sPu6EQUT1wQCIiBy8E3ygrb2eU2CCX7TL+6qa9viIKO4xACIixwZAlYFVoOs7BdYBx3HDsmHA7C5Nf4xEFNcYABGR89T6gKTUiAyQVht9rSARMLlRC/irAX+wfZ6IKBwDICJynn6XAzMO41dJ90et9wm/z1gHyNRJRkRkhQEQETlWbcSaP5EBkLkuqApe1Kq3Na4FRERRMAAiIscyT2/J23VshSHWDPJ50vSr3BCViKJgAEREzrPhVeDVq3B57dKQu622wggPimrcakNUToERkT0GQETkPIc3A98sQg/tYN3rAIUVRtcwA0RE9cCVoInIsStBG3t7WWx8qoTvEL835xzkdh8MpGQ18UESUTxjAEREzhOo3ynXkhs8Bba4128xeGzfpj0+Iop7nAIjIscGQGVhGaD6TIFxKwwiqg8GQETkPIEW9oqwDFB9tsIwOses0kVERAEMgIjIsTVApadQA3TZ7v8DZrUG1sxp4oMkonjGAIiInEdsZWFZAxR9IURj8USxlQbXASKiKBgAEZHz3LwEtTOOYUXtEHnT7bKfAgvfHqOa6wARUT0wACIiR/K7PMa2Fkked50rQbsCQVK1K7AQIjNARBQFAyAiciRzvU9SIAVkuRdYoAvM69bfzqqMDJBeSE1EZIXrABGR8yyYhCS40QpjcBLZ8MoMkN8IdsxUUJTkcaHaD1QZGSBuhUFE9hgAEZGziIBmy1vyzcmD78u7vB77DJCqAfIEskRGAMQMEBFFwSkwInJkC7xQjpSQ4Ma6Bki/1LNEwImk9kDPC4EOA5vneIkoLjEDRESODYAqobfBJwXqe6IthKjqhLZnnA38bHLzHCsRxS1mgIjIWQK1O1pSqtEFpqbAoq0DpDJAVlkiIqJwDICIyFkCtTtaUprR3u4OZHesVoJWW19EmyYjIgrHAIiInCWwfo/mTZeXbpcLnsAiP1GnwAJZojYV3wIPdgcePbP5jpmI4g5rgIjIkTVAtYEMkAh+RBBkPwWmhawD5BNva5WF+nYYREQ2mAEiImfpfh4w4ygOX/WBvCniGrXKs3UGCCEZoIpA5xiqy7gjPBHZYgBERM4iop2kZPi9GfKm25QBCt/3y1wDpLrAKl1qB3kN8FU122ETUXxhAEREjqQKnuUUWOCdynIrDKMGSB9UqTJAAvcDIyIbrAEiImfZuRT46l/IbHWWmA+THWCqCDpaG7zKAPngATzJgL86EADlNe/xE1FcYAaIiJzlyFZg4wKkFnxmtLe7VBdYlL3A1DpAMnMU6CDjdhhE1CgB0DPPPINBgwYhOztbfowcORIffKAXKgqVlZWYMmUKWrdujczMTEycOBGHDx8OeY59+/Zh/PjxSE9PR7t27XDPPffA5wvt1lixYgWGDh2KlJQU9O7dG/Pnz2/IYRJRC1gHyO/Ru8BEYieQ3KljCsy0DlDX4UCPCwC3p/mOm4habgDUpUsXPPjgg1i/fj0+//xzXHLJJbjiiiuwefNm+fhdd92Fd999F6+//jpWrlyJ/Px8TJgwwfh8v98vg5/q6mqsXr0aL730kgxuZs6caYzZs2ePHDNq1Chs2LABd955J26++WZ8+OGHjXneRORUNeEBULAI2jIACmSFgttlaMC1C4Ab3gNan9Z8x01ELbcG6PLLLw+5/ac//UlmhdauXSuDoxdeeAGvvvqqDIyEefPmoV+/fvLxESNGYPHixdiyZQuWLl2K9u3bY8iQIXjggQdw77334v7770dycjLmzp2Lnj174i9/+Yt8DvH5q1atwmOPPYaxY8c25rkTkYMDIGMdILcogq7HQohu+2kyIqJGqwES2ZzXXnsNZWVlcipMZIVqamowevRoY0zfvn3RrVs3rFmzRt4WlwMHDpTBjyKCmuLiYiOLJMaYn0ONUc9hp6qqSj6P+YOI4jcA8iWZM0Coczd4NQVm1SpPRPSdA6CNGzfK+h5Rn3PrrbfizTffRP/+/VFQUCAzOLm5uSHjRbAjHhPEpTn4UY+rx6KNEQFNRUVwl+hws2fPRk5OjvHRtWvXhp4aETmpBsidKi/FzFZwCixyuD+sCFp2iv3nF8BDvYCNbzTfcRNRyw6A+vTpI2tz1q1bh9tuuw2TJk2S01qxNn36dBQVFRkf+/fvj/UhEdF32ArD38CtMNQUmOwCE6tAlx8HqpgJJqJGWgdIZHlEZ5YwbNgwfPbZZ3jiiSdw1VVXyeLmwsLCkCyQ6ALr0KGDvC4uP/3005DnU11i5jHhnWPitug6S0vT3xCtiIyU+CCiOPfT+TKAOXpQBEJfyfoftRWG5W7wYVNg8rY3LSSbRETU6OsA1dbWyvobEQx5vV4sW7bMeGz79u2y7V3UCAniUkyhHTlyxBizZMkSGdyIaTQ1xvwcaox6DiJq4bypQEZr1HhSg7vBRyuCNrbCMHWBqXWAAtkkIqLvlAES00zjxo2Thc0lJSWy40us2SNa1EXdzeTJkzF16lTk5eXJoOaOO+6QgYvoABPGjBkjA53rrrsODz30kKz3mTFjhlw7SGVvRF3RU089hWnTpuGmm27C8uXLsWDBArz//vsNOVQiinOqm8s8BVbvdYCS9X3EUFPWbMdLRC04ABKZm+uvvx6HDh2SAY9YFFEEP9///vfl46JV3e12ywUQRVZIdG89/fTTxud7PB689957snZIBEYZGRmyhmjWrFnGGNECL4IdsaaQmFoT7fXPP/88W+CJEsXiGXIKzNvpenlTTIGpLrBoW2EEV4I2TYExA0REjREAiXV+oklNTcWcOXPkh53u3btj4cKFUZ/n4osvxpdfftmQQyOiluLrBUDpYbgvu0LeFHGN2gpDBjdhVFCkpsnkbWMrDGaAiMga9wIjImcJZG1q3PVbB0gVRnvNXWA5XYCOQ4AcLodBRNa4GzwROYcIcAJZG587sgg6aht8YApM3h58tf5BRGSDGSAicg5/DaD55dXqQAYoZDd4iy6w8JWguRA0EdUHAyAicg5T11aNJ8W0G3w9usDUFBgjICKqBwZAROQcqmvLnQQfvBE1QFYLIaqAJ2QdoP2fAo8PAub/sNkOnYjiC2uAiMg51MrN3nQjsyOmwMRaQIJVcketF+RV6wCJ21otULgXcPFvPCKyxgCIiJwjrydw907AV4naPZpFDZDVbvAWRdBcCZqI6sAAiIicw+0BMtvKq/7aA/LSFdIGH/kp6j7VKSanyYwAiHuBEZE15oeJyJGMKbA6iqCDNUCmTrFkBkBEFB0zQETkHPkbgC//AbTti1qMNTI7gfpmY+PTuqfAAlth1PoAXzWQlNxsp0BE8YEZICJyjmPfAJ89D2x7z8js6FNgUXaDV3uBGRkgEQAFNkMVOA1GRBaYASIi56gu1S+9GUa2p/67wavNUDXA4wXanAEkpehZICKiMAyAiMh5bfDJog0ewSmwKLvBq7WBQlaCFgHT7Z8111ETURziFBgROYearvKmG4GN29QGb7XKs7rLGygUsloskYgoHAMgInKOwEaoSM4wprZCt8KI/BTzgonm20RE0TAAIiJHZoCCbfBiN3jUWQNkrAStxiy4HnhiCPDtJ81z7EQUVxgAEZEja4D8gS0uxBSYux5bYQTb4AMPFOcDJ/cAFSeb4cCJKN6wCJqInGPsH4ELpgKpuaj9VA9cxMyWsRVGbT12gzf64tNDN1glIjJhAEREzpHWSv+QgcyJiC4wf9Q2+MCgQLeYywiAAnVFREQmnAIjIkcKFkGLGqAoU2CB+5LUctEqC6S2w1DTakREJswAEZFzfPJXoPw4MPT64EKIUXaDN0+JqSJofZx5CowBEBFFYgBERM7x5cvAse1A70vh19oZGaDgbvBhAZDptiqCNu5nAEREUXAKjIgc2AYv1gGCKQBSBc6hw8010aoIWr9fA7LaA616AilZzXDgRBRvmAEiIgcuhJhumgLTO8GstsIIyQCZAiBZA3TBb/QPIiILzAARkeO3whAfDZsCa57DJaL4xQCIiJyh1g/4KvXryRlGy7t5Ciw8sDHfDimCZgRERHVgAEREzmDu1vKmGy3vnnoWQatWeeP+nUuBuecDb09p6iMnojjEGiAicgZjvR4X4E0LnQKzWQnafFsESmKYiH1k9kg8X8FGWVBNRBSOARAROUN6a+D29XomyOUyTYHZ7wZvvq2mysTnyU/lStBEFAUDICJyBk8S0Ka3cVN1fNV3CkxkisRYP7TQlaC5FxgRWWANEBE5klUXWPhWGMHtMvTbgURRYCHENP0Gt8IgIgvMABGRMxzbAWx4FWjVAxg2yVj0MPpWGPqlmiJThdDyflX7w60wiMgCM0BE5AxHtwOrHgW+/EdEdsfYDT68CNrUKm++DMkAMQAiIgvMABGR4xZBjNgN3rYIWk2T6beNQEncn5wBZLTVn0+sMeT2NNOJEFE8YABERA7bBiMjJNsjprVUZid8Kwx108gAGbVCGpCeB9yzs9kOn4jiC6fAiMjxGaCQ4marQmlVA2STKSIiCscAiIicQXVrBdrXjQJnUwbIX2cXmNo1nhEQEUXHAIiInKFGTYFlygu/eR0gt91u8AiZ+lL7oRqZoteuBZ69CDjGqTAiCsUaICJyVg2QmgIzaoDCurtMNLsusED2CAVfA4X7gIqTzXMORBQ3GAARkTOcPxUYco3euWXKAIlprYjAJsC8XYZ+GRYocS0gIrLBAIiInCG7o/4RNr2lb4VRv4UQ1VSZCoyM7TBUdomIKIA1QETkSMEpsLr3AgvvAjNqhYwNUdVO80REOmaAiMgZvvgHUHYU6H8F0Po0o5NLxDSqyDm8uSu4DpC6VF1gCCmoZgaIiMIxACIiZ/j8RSD/C6BdPxkAqeyOeSHEiHWATHVCQjBQ4hQYEUXHKTAicgY1TRVYCdoIgGQNEKJuhaE2QTXGqYFprYD01twGg4giMANERM5aCDHQuWWs8hx1Kwy7LrDAgPF/0T+IiMIwA0REDlsIMT1k1WfzVhiRu8EHx5gvjS4wIiIbDICIyGEZoPSQ7I5YCFFNcYVPgZkLpfWx1rVCREThGAARUeyJBX18FZa7wYusjt0UWJ01QNsXAfN/CCz9Q/OcBxHFDdYAEVHsmdfpCWSAVAAkghu73eCDbfDhXWCBAeXHgG8/BpJSm/oMiCjOMAAiotgTAcoty/VpMG9axCKHwfV9NJspsLAaIDUukE3iQohEFI4BEBHFnicJ6DzMtsA5OAUW+mnBIAk2K0EHAiBuhUFEYVgDRESO3wpDFEJHmwJTNUBGt5ixECIDICKyxgwQEcVe0QFg4xtAdidg0M8idnpXU1x2CyGqxyO6xdRmqNwLjIjCMANERLF3bAew9D5g1eOR01vRtsIwOsX028Y4owZI7QVW2vTnQEQtNwCaPXs2zjnnHGRlZaFdu3a48sorsX379pAxlZWVmDJlClq3bo3MzExMnDgRhw8fDhmzb98+jB8/Hunp6fJ57rnnHvh8vpAxK1aswNChQ5GSkoLevXtj/vz53+U8iSgutsFID+mMj9gKo66FECP2AssAPMmAJ6Wpz4CIWnIAtHLlShncrF27FkuWLEFNTQ3GjBmDsrLACq4A7rrrLrz77rt4/fXX5fj8/HxMmDDBeNzv98vgp7q6GqtXr8ZLL70kg5uZM2caY/bs2SPHjBo1Chs2bMCdd96Jm2++GR9++GFjnTcROXgRxPA2+IgtLgKMxRKNLrDQz5VTar8/Ctyzo8lPgYhacA3QokWLQm6LwEVkcNavX48LL7wQRUVFeOGFF/Dqq6/ikksukWPmzZuHfv36yaBpxIgRWLx4MbZs2YKlS5eiffv2GDJkCB544AHce++9uP/++5GcnIy5c+eiZ8+e+Mtf9D18xOevWrUKjz32GMaOHduY509EjtoGI8Nip3fABespMBXnGCtB20yVERE1ag2QCHiEvLw8eSkCIZEVGj16tDGmb9++6NatG9asWSNvi8uBAwfK4EcRQU1xcTE2b95sjDE/hxqjnoOIWhjVpm7KAAW3wnDBbXSBhX5asFBadYFZZ4qIiBqtC6y2tlZOTX3ve9/DgAED5H0FBQUyg5ObmxsyVgQ74jE1xhz8qMfVY9HGiCCpoqICaWn6QmlmVVVV8kMRY4kozqbAki2mwFwuqHjGdjf4QICk2uVDFkz89y1AaQHwoyeBVj2a8CSIKCEyQKIWaNOmTXjttdfgBKJAOycnx/jo2rVrrA+JiBo6BaYWLgxb5dmo7bHZC0xlgCw3Q927GtjzX6D8BF8PIvpuAdDtt9+O9957Dx999BG6dOli3N+hQwdZ3FxYWBgyXnSBicfUmPCuMHW7rjHZ2dmW2R9h+vTpckpOfezfv/9UTo2IYmHYjcD1bwNn32i5yGFEe3tYp5jRBWa1ZYbKKnE1aCI61QBIpJtF8PPmm29i+fLlslDZbNiwYfB6vVi2bJlxn2iTF23vI0eOlLfF5caNG3HkyBFjjOgoE8FN//79jTHm51Bj1HNYEe3y4jnMH0QUJ1p1B3pdDLTtY9ylsj2eKFthmBdLlGPdVgEQ9wMjou9YAySmvUSH19tvvy3XAlI1O2LKSWRmxOXkyZMxdepUWRgtgpA77rhDBi6iA0wQbfMi0Lnuuuvw0EMPyeeYMWOGfG4RxAi33nornnrqKUybNg033XSTDLYWLFiA999/vyGHS0RxzFjkUPyZVmvd3WUulLbtAjP2A+NiiER0igHQM888Iy8vvvjikPtFq/sNN9wgr4tWdbfbLRdAFEXJonvr6aefNsZ6PB45fXbbbbfJwCgjIwOTJk3CrFmzjDEisySCHbGm0BNPPCGn2Z5//nm2wBO1VFvfBUoPA71GAa1Pi6jv0VyaTQ0QQneDNzJApkHGfmCBQmsiooYGQOEdGFZSU1MxZ84c+WGne/fuWLhwYdTnEUHWl19+yReJKBGsexb49mNg4gumAAhGdke989TWczf4kAwQp8CIyAI3QyWi2Ku2WAjR2OfLBY/L+o8wVRQdvhVGRBG0ywP4q5v2HIgorjAAIiLn7AUWWAjR3O0lYhrNZZcBCoxRNUBW6wBd/lfgR08Fl4smImIARETOWggxI2IKS0yBqXimtq51gKymwNyepj12Ikq8rTCIiBqF6tAKBEDmYmeR3VE1PuJu8zRYcDf44Fj5+dwLg4jqwACIiBwUAGWGLHCosjqqzV0+ptnXAKkMUEi32K7lwGvXAisfatJTIKL4whogIootX3WwQDklMzID5HIZbe5qessTtjt8xFYY5iip+BCw7T3AV9kMJ0NE8YIBEBHFlqjRue4tvRMsJTuihkcshOg2ZYRq6zUFZnr+QFCFKi6ESERBDICIKPYB0GmjQu4yZ3DMW2EI5tmtehVBGwshMgAioiDWABGR45iLmN1hAZD5MaMGKPBOZr0OUJZ+yQCIiEyYASKi2BI1Ot98AGR1BPqMk3eZ4xfZBeauawosNAPkMz8Bp8CIyAIDICKKrSNbgPfuAtoPMAVAoZucmjNAIV1gEUXQgftDMkCBGiBmgIjIhFNgROSoFnjzFJbK6ITWAGkR1yOKoM01QCmBKTB/DeD3NdlpEFF8YQaIiGJLdWelRAZAKu4xLQMUUt+jAh1XeBG0OQOU1gr43wIgKZXbYRCRgQEQETkuA6QSOGoKLHQdoOCnmneMN1+GZIDE53rTmu74iSgucQqMiJwRAJkzQKoGyBT4qOAmdCuMsCkwtRI0t8IgojowACIiZ0yBqXZ1UwCjanrkdVfdW2EkeSzWARIWz9C3wzi6vYlOgojiDQMgInJcBii8uNk8DWbVBq8es80A7Vymb4dRfLCJToKI4g1rgIgotobdCPS4AGjdO3IKzCIDFLIQojEOoTVA5q0wQlrhy5rmHIgo7jAAIqLYat9f/zAxpsDMNUCB6+bZLXU96lYYAhdDJKIwnAIjIsepDWRwzAGQ2yK4CbbLu+y3whC4GCIRhWEGiIhia9v7QE2FPg2W1d5yJWhBxUKhNUCh3WLGStARGaBAgXVVSROeCBHFE2aAiCi2lj0A/HuyviVGWA2QeQ8wld2x3goD0YugmQEiojAMgIgotlRhssrSmNrbPfWcAlPBUbAIOjwAyghtuSeihMcpMCKKreqSiJWgw3d5twuA1K7vSe46iqAvvFv/8KY31VkQUZxhAEREjt0LzHIhRFOLu5EpUrvG15UBIiIK4BQYEcWOrwqorbHIANVvCkxlgIy9wFQNUFj8Q0QUjgEQEcWOuSYn2X43+NCtMCK7wIwpMFUoHZ4BKtgIvHkbsGxWE5wEEcUjBkBEFPv6n6Q0wBOckbdug4/sAvMFUj3uuqbAyo4CX70KbF/UVGdCRHGGNUBEFDvprYGfvgTU+kLutgqAjOxOtAyQXRG02mhVBVxElPAYABFR7IjW9zOvjLhb7eUV2gUWulGquQZIjVPrBkVkgLgQIhGF4RQYETlO+AKH+vXIKTAV6CR5wougwwKg1Gz9srI4dDMxIkpYDICIKHaO7wI2/Rs4uD7k7vD2dkElg8zZnfBNU22LoFMCAZDmB2rKm+JMiCjOMAAiotjZ/RHwxk3Ax4+G3G1sheGKXgNkZIACc19GEXREDVAG4PIEs0BElPAYABFR7KhgRGVobDI75uvm2EaNU5ugGkXQpsUSJXE/64CIyIRF0EQUO1XFoTU6ASrIsW6DN2WAjG4xd/S9wITbP9MzQdwOg4gYABFRTFWVRM8AWW2FESUDZOwGb1XonNmucY+diOIap8CIKPZTYGEZICOzY+oCC2Z3aiMWQgzPAEUUQRMRhWEARESxnwILywBpUYqg1RpBVnuGqUyQZQZo/UvAW78Cdn3UyCdBRPGIARAROS8DpBZCNE2BJVllgMJ3g1dTYFYZoG9XARteAQ5vauyzIKI4xCJoIoqdi6YBRQeATkNtpsAiM0Aq6DFPdRkLIUabAjMvhkhECY8BEBHFzmmjLO9WAYza2sK81o85uxOxFUa0Img1zaam3YgooXEKjIgcvBWGRQYoUPgcuhBieAbI4kmZASIiEwZARBQbIsjZ9B9gx1LAXxPyULC93aoGSLMdZxRKW2aAAjvCMwNERKwBIqKYqS4F3rhRv/67Q4DHa9vdZVcDFFwIsR5F0Ck5+iUDICJiBoiIYkYVI7uTAG9ayEMqflGrP5sLnc1dYHYZIPkc4UEQp8CIyIRF0EQU+zWATIGO1QrP+nV3ZAYoPAAyPY/IDrlhet4eFwBTtwKpgUwQESU0BkBE5Kg1gMzZm3rXAKkuMFPAJB7zBjaAl5LT9Q8iIk6BEVHsM0CB4mQTf7QusHoUQYdvmkpEFI5dYEQUG5VFocXJJirGMQdA9ekCM4+PKIT2VQOLpgNvTQF8VY15JkQUhxgAEVFsM0D1nAKzWgfIFyiIDl8HSH+OsCcVxdZrnwE2vBwMvogoYbEGiIhiQxQl/+hJIKtjvabAwvcCExumGpkimyLoEKJASEy3icBL1B9ltmuCkyKieMEAiIhio83p+kfUlaDtu8DMU1wqODJvnmq1FpA/ORseGQAxA0SU6DgFRkQxUV7tw9sbDqKoPHQVaNsuMGMdIC0iw2M1VWZVBP1tebK8PHniSCOeCRHFIwZARBQTyxe/g38veAl/X/JpxGP+QP2OOaMT3gVmzvCEBEBRVoM+4tMXXCw8frjxToSI4hIDICKKif7fPI2/J/8ZuQWrIh6z2gojvAvMLgBSawGFB0CVNX6crM3QHys70bgnQ0RxhwEQEcVEcrVeh3PcrwclddcAqQxQbWQAZF4vyGU9BVZYXoNCTf+/tPKTjXkqRJQIAdB///tfXH755ejUqZPcp+ett94KeVx0ZsycORMdO3ZEWloaRo8ejR07doSMOXHiBK699lpkZ2cjNzcXkydPRmlpaciYr7/+GhdccAFSU1PRtWtXPPTQQ6d6jkTkQKn+Enl51Be5OrMKbsxTYPXPAFlPgRVWVONR388wvPIpbO55QyOfDRG1+ACorKwMgwcPxpw5cywfF4HKX//6V8ydOxfr1q1DRkYGxo4di8rKSmOMCH42b96MJUuW4L333pNB1S9+8Qvj8eLiYowZMwbdu3fH+vXr8fDDD+P+++/H3/72t1M9TyJymDSfvg5QQU3oRqj2u8G7Q9YBMoIkV+imqXZF0CfLanAMOTiMPJT42QBLlOga/C4wbtw4+WFFZH8ef/xxzJgxA1dccYW87+9//zvat28vM0VXX301tm7dikWLFuGzzz7D2WefLcc8+eST+MEPfoBHHnlEZpZeeeUVVFdX48UXX0RycjLOPPNMbNiwAY8++mhIoEREccrvQ4ZWJq8erkmNeFgtYmjOAHltusCSzBuAhRRBhz5nUUW1cb28ytdIJ0JE8apRa4D27NmDgoICOe2l5OTkYPjw4VizZo28LS7FtJcKfgQx3u12y4yRGnPhhRfK4EcRWaTt27fj5EnrufuqqiqZOTJ/EJFDmdbhya9MbdBeYDWBAEhlgsLiH9spsJPlNTjNdRAzkv6BM3c913jnQkRxqVEDIBH8CCLjYyZuq8fEZbt2oSuwJiUlIS8vL2SM1XOY/49ws2fPlsGW+hB1Q0TkUBX6HzLFWhqKqjWZPbZeBwi2K0HX1pEBipgCK69GO1chbk76AH2OfND450REcaXFdIFNnz4dRUVFxsf+/ftjfUhEZCejNX6PW/Gg7xqZqamo8dcjAxRaA6TWAzJ3iunjrDNAYsHFQi1TXk8N1B8RUeJq1ErADh06yMvDhw/LLjBF3B4yZIgx5siR0FVYfT6f7AxTny8uxeeYqdtqTLiUlBT5QUTOV5uSi5erLoRK0pRU+pCenNSg3eDVZZI5TWReB8giA6QCIFmALR43PT8RJZZGzQD17NlTBijLli0z7hO1OKK2Z+TIkfK2uCwsLJTdXcry5ctRW1sra4XUGNEZVlMTXCJfdIz16dMHrVq1asxDJqIYKK/xG8GPCoDqvRt8WABkDpJCpsDC2+BFBgj6OkBJ8AHVehE2ESWmBgdAYr0e0ZElPlThs7i+b98+2Yp655134o9//CPeeecdbNy4Eddff73s7Lryyivl+H79+uGyyy7DLbfcgk8//RSffPIJbr/9dtkhJsYJ11xzjSyAFusDiXb5f/3rX3jiiScwderUxj5/IoqBikPbcaH7K3R36TV9JZU1da8D5LHJAIXNgdmuA1RegwqkoEpLCqlDIqLE1OApsM8//xyjRo0ybqugZNKkSZg/fz6mTZsm1woS7eoi03P++efLtnexoKEi2txF0HPppZfK7q+JEyfKtYMUUcS8ePFiTJkyBcOGDUObNm3k4opsgSdqGdyb3sDfkx/FK75L8b++ySgNa0tvyErQ5ixRSBt8+ErQsg3ehSJkoh0K9QAol80SRImqwQHQxRdfHNGxYSayQLNmzZIfdkTH16uvvhr1/xk0aBA+/vjjhh4eEcUBf7m+F5eakoqYAqvHXmA+uwBILYQYtg6QaIOX/6eWIbvBmAEiSmxcDpWImp1WXigvVVFyfabAjC6wwGPBNviwKTCLDJD4o62wXF8IcXLN3cjLzsLb3c9rgjMjonjBAIiIml+lXn9TiEzLDFCg0z16BshYCNEuAxQMgMqr/agJjN+vtUdJjRfweBv/vIgobrSYdYCIKH54AgFQcWB39vAASE2zm9c4NGqA/HVkgCyKoEULvFkZt8IgSnjMABFRs/NW6QHQcS1LXoYXQVu1uNvVAEW2wSNiCqyoQp9iS0lyY2jtRox2fQHfhhIkDbmq8U+OiOICM0BE1OxSq/Ui6OPIiVoDZL0OUGArDGMhxLqnwCoDK023yUzBQNduTE76ALU7ljbBmRFRvGAGiIial6ZhUde7sHnHLhR72wDVVlNgiMwAeeqXAbIqgq6s0YOmjBQPispy5fXa0qNNcXZEFCeYASKi5uVyYXXmWPzNfzlyc3Ktp8Ci7QVWx0KIVnuBVfn0DFCq14OypMBq8mUMgIgSGQMgImp2KuDplJsmL4vDu8AspsDs9gKz7QIzZYCqAhkgUQNU7s3TP6/8WKOfFxHFDwZARNS8ig6ix8lPcJrrIDrm6CvEl4bVABkLIVp1gRlTYLXR1wEyLYRY5VMBkAdVKXoGKKnyeHCujYgSDgMgImpeu1fgnmMzMDPpH+iQo2eA7LbCECvL22WAgkFS/YugRQaoOqW1vO6urQEq9QUZiSjxMAAiouYVqL05hmy0zkgOKVJWVPbGvBBicB2g2pD1gMIDIKsiaJUBEjVAyalpKNbSA8fCaTCiRMUAiIhiEwBpOchN94YUKSu1ljVAbusMUPg6QIF3NasiaJEBSk9Owo+r/4D/XPIRkNer8c+PiOICAyAial6BrMtxLRu56cEMkHmT5eBu8KYMkCe8BqiOKTCLNvgUrxsZyR7s0jrjhCsXcHua6CSJyOm4DhARxSQDdFxkgNK8IdNUYooqtA0e9jVANgshBougrTJAHqOuqKwqNOtERImFGSAiik0ABJEBCg2Aok+BBTNAIltkuxWG1TpApgxQmteDi9xf4bydjwBb3mmSUyQi52MARETNSgtMgR3TspGV6jWyPFWBTq2QDJBFDZAgYhvbhRBdFusAmdrgRQA01P0Nzil4DdizsilOkYjiAAMgImpWVaPuw59qrsF+rZ0MRkRQEt4JFljix7IGSK0BZLcQYnA3eFi2wacle2T9kcTVoIkSFmuAiKhZlZ3xYzznzzICklSvGxU1/pBOMKsOL3OmRwQ/KkvUkAyQqjES9UcS9wMjSlgMgIioWYlgRxCBj8jW6EFJTUgGKJjdCX6euR5I1P/47dYBMtYLsm6DF/HREU3fgwwlh5rgDIkoHnAKjIiaT9EBuHcsRi9XvpGNEUGJUBmSAUJEcGPOBongR2WAItvgA2Os2uDFFJjXg0PQ9wNDcT63wyBKUAyAiKj57FqOTgsnYUbSyzIQEVQgpDq17NYBEpkdFevUmGqAIhZCVFNgNrvBi//3sBYIgPxVQPmJJjlVInI2BkBE1HyKDsqLAi3PCIBSApeqUDlkCiwsuDGvBh1cCNFtXQRt2QXmRmqyBzVIwkmxEKJQrB8TESUWBkBE1HwCwcYhLc8IfCynwOpY5VnU9wTHoO4MkLEOkJ4BEu5K+xNwzy6gw8DGPksiigMsgiai5iNqbkQGCCID5K7HFFjop5tXg7bLAFkthKiCK1UDJOyo7QRktGnsMySiOMEMEBE1ewCUr7WW6/EIqRYZoOBWGGEZINN+YEYNUNi7mOUUWCC4kjVAgf9XdaMRUWJiAEREzZ8B0vKQmhReAxS5EGL4FJg5A+S3ywBZFkGHdoEJp1dvAz74LbDu2cY+SyKKAwyAiKh5VJUAVUXBACgsA2S5EKJdDZDoArNZCNE6AxScAlNTbh39+cC6Z4Ct7zb6qRKR87EGiIiahzsJmPA8Vn+1BWWb0yLa4FUGSGx0qoKXsBmwkC4wu4UQVQbIvBWGeSVoNQV2SGsdkpUiosTCDBARNQ9vGjDop1jX4efGStDmLjCVpan210Ilb1RwFJkBqnshRDUFJi7Fc6r/S2WcuBgiUWJjAEREzUoVO0cshBjI0lRUB6fC1JioNUCu6FNgKvhR9UZJHjeSPW4c1lrpd/oquBgiUQLiFBgRNY99a4GKQiSVpYYEPioTpBZCVN1ZXo8L3rAWL/M6QMEi6OgrQZsXWFTZJvF/FvuT4cvogKSyAuDkHiAjMCVGRAmBGSAiah5rnwb+eRX6HFsibwb3AgtdCbo8kAEKz/6Er/FjGwCFZYBUZkncrwIqVQdUld1D/6Tjuxr7bInI4RgAEVHzOL5bXuR7OodNgbktp8DSkyMT1EnGOkCmvcDCu8CMIuhAAGTaCFVR/3d5Vnf9jsK9jXeeRBQXOAVGRE1PZGNO6AHQAVfH0AyQ1yYDFMjSmKk1f0JXgnZZBkmqld68CrSi/u89Z/4aba/4PyCd019EiYYZICJqeiUFQE0Z4PLgANrJu9KSw7rAVAaoxn4KTBVB6ytBB6e26pMBMneUqetF3jb6dhjh/fZE1OIxACKipnciUGOT2w1lPj3YiFwHKFAEXe2Tl+nJddQABVrlw7vARPG0UBMYoBZYtJoC43YYRImLARARNT1VZNz6NCPoSLFZCDHaFJhVBkhNeSmZKV55WVrlC9sGI/h86rkrxZjFM4BXrwbKjjfmGRORw7EGiIiaLwOUdxoqDodOcQWnwOruAhNr+Agi+FFTXOEbpmam6m9rpZW+kMxSSqDY2vzcFSI42vwWULQfOL6TrfBECYQBEBE1vbOuB9oPBPJ6ofJrfT8wu60wVMCSHi0DZFoHKHwvsMyUJMsMkNp81fx/ymxUXq9AALQD6Da8kU+ciJyKU2BE1PTa9JbbYKDLMCPACV8IMSIDZNEGX591gLJSwwMgiwxQoABbtty3H6DfeejrxjxjInI4BkBE1KzCu7xUbY7q1oo6BRZSA2QdAGWEZYBUZsmqCFoGYx0H63ce+qoxT5OIHI4BEBE1rYJNwOongQPr5U7vRgYokIUxtsIIZGqiTYGFdoFZB0BqCqzaVyuzP1VhRdcRXWCdhgSOcyNQG9w2g4haNgZARNS0dizWO63WPCU3Jg0kboJTYIEMkGhbF4FNeaANvq4uMFEHFC0AEsqq/KYuMNNCiKoLTARArXsD3nR9nSJuiUGUMBgAEVHTUlNLHQejsjq4M7sxBWaqzREZGzUFZp0BCnaBqZWerfYCU58rOsGCU2BWGaBawO0BOgwEUnKA4gONddZE5HDsAiOiZguASqpq5NVkj9vYmNTcnSWClYp61gAZW2FYrOIsskAikBL/X1GF/n/mpOnrA5mfW/1fuOZfQGouV4QmSiDMABFR026BcXIPAJcMgI6VVsu722QmB9+E3C4ZEKkpKaNI2ioDFFj00G9ugw9bCDF8LaDjZVXyeuuM5MiFEAP/F9JaMfghSjAMgIio6exeoV+KTqv0PBwt0YORtlkpIcPM+4GVR9sN3qILLHwhRCHL1Al2PBB0tTYFXar+SNUbGcS0GguhiRICAyAiajq7luuXp10iL2wDINOUVLQpMKt1gJICdUGWGaAqH46VBjJAmcH/s1W6HgydKNODI+mj/wMe7Qdsffc7nDARxQsGQETUdA58Vq8ASE1PiWClvKbuLrCaWtHirhc3e5Osa4CEEjkFVh0xBdYxJ1VeHiqqlK35UmUxUHIoGLQRUYvGAIiIms6v1gE3LAS6nitvHi2tlJdtTdkYoWOuCkgqUBHoFIvWBVZW5TOyNx2y9c+12hBVBEAnA+PamP7Pdtn6dRFEnSyvCQnSZACkgiIiarEYABFR00lKBnp8D0jSA45jJdWWGaCOOWnyMr+wEhWBupxoe4EdOFlhjDF3d4Vvh5FfWGF0i7XKCI4TLfGqEFsEXZI4zuRMfV+wvau/+7kTkaMxACKixuf36R9hjgbqccIDoE7GlFQFysO2yrCqAfr2WJn+eblpcNm0wctxx8uMgMi8DpA56Coo0rNSSM4ABkzQr3/x9waeMBHFGwZARNT4tr4DPHYmsGZOyN12NUAdAgHQt8fKjdknqxqg9oHprm+Pl8vLzrl6EBNO7Qe2NzDOPP0V/n/mqwBIGDpJv9zyFlBxsn7nSkRxiQEQETUuvw/ayoeA0oKQIEIUGxsBUGZo3Y7I5Ai7j5Ua91llgPp0yLT8PLsusH0nyiMKoMMLoQvUFJjQeRiqWvcDfJXwffJkfc6WiOIUAyAialT7lsyB6+hWnNQyMb/2B8b9ZaLFPTC91SYr2TIYUQslioURkwKLI5r1bpcVcrtzoHjabh0gxbwGUPD/TDM6wZSiSh/uPXkFXvFdil9+cw4Ky01t8kTUojg6AJozZw569OiB1NRUDB8+HJ9++mmsD4mIoije/RnarP2TvP4X30/x55VHZCGyoLI/GcmeiEUOVTCiWE1/CaLgWQVLUTNAYQFQXkbkFJjRCl8YDIAeW/IN3iofhP/1TcayfbX4/dub+XoTtVCODYD+9a9/YerUqbjvvvvwxRdfYPDgwRg7diyOHDkS60MjIit718D96s+Qjip8kTQE33SeKDM+sz/YFrX+RwU8uenBLi2rDjDl9PZZ9Z4CU8xbb4TXABUU6wHQziOl+MfavfL6tMv6QNRWv/fVAeQvfgIo5fsOUUvj2ADo0UcfxS233IIbb7wR/fv3x9y5c5Geno4XX3wx1odGROGqSuH/+4+R6TuBzbXdUTlhHu67chBE09a7X+Vj3e7jOHCy3DYACs8C2WWAhD7tg3VAneuZAYpWAyQ6z0R90iMfbperS4/u1x6/urg3rhjcCbOTnken1TOhzRkOrHgQOLqdawQRtRCO3A2+uroa69evx/Tp04373G43Ro8ejTVr1lh+TlVVlfxQiouL5eVnf7sdGWmRb37r21yBo2m95PXOZVsw6MQi2+PZ0Ho8DqX3kdc7lG/HWcfftx27MW8MDmQMkNfbVuzGOcfetB27OfcS7M06S17Pq9yPEUf+ZTt2W+6F2J2tLyaXW3UI3zv8su3YHTnn4Zuc78nrmTXHcNGhebZjd2edg62tLpbX03xFuDT/WduxezOHyPMTkv1lGHPwKduxBzLOxIbWP5TXPbXV+MH+R23H5qf3xfq2Vxq3L98723bskbReWNfuKuP2+H0Pw6MFFrILczylKz7pcJ1x+7L9jyOlVm+LDleY3BErO95k3P7+gaeQ7iu0HFvibYPlnW81bl9ycC6ya45Zji1PysbiLv9j3L7o0AtoVZVvObbak46FXX9j3D6/4O9oW/mt5Vi/y4t3ugd/PkYefhUdKnbAzpvdZxqbfZ575HV0KbOf2hHP63PrQcqwo2+he+mGkMe9tZVI9xchxV+OZ/u+CM3lkbUy59eMgEurwaozfovH+us/Wz8/txteWbcPUxd8hZJK/XXq2yHb8v8VrfBbDxXXOwMkTkd1hdmtA6SYt8FQxOeKtnqxA/1P5q7B+r0nZcAmsj/Cb8b0wZQtl2Nw7S70q9gPrJgtP2pcKShM6Ygibzv887SH4Xfrmauhx95Bx/Lt8romNn9VAl/3DzvfAb9bfy8aeOJD4zUIjg1+zvJOt6DakyGv9z+5HN1Kv7b9eqzseCMqknLk9T6FH6NnyXrbsas6/D+UetvI672L1qJ38VrbsWvaXYWilI7yes+Sz9GncJXt2E/bTsSJ1K7yerfSr9D/5EeRgwKn90WbHxnvvZ3KtmDgicW2z/tV63EoCLz3ti//BkOOL7Qduynv+ziYcaa83qZiD4Yde9t27Nbci7Eva4i83qryAM49+obt2G9yzsee7LPl9ezqAow8/Jrt2F3Zw7EzZ6S8nlFzHOcX/MN27LdZQ7E990J5PdVXjIsO2f9xvz9zELa00hfq9PrLcUn+32zH5qf3w8bWY+V1d20Nvn8wtBPT7HBab2xoo79PC2P3P2479lhqj5D36dEH5ti+955M6YxP2/3UuD3q4LNIrjU1G5gUJ7fDmvbXGLcvzH8RtaXW76cJEQAdO3YMfr8f7du3D7lf3N62TU+nh5s9ezb+8Ic/RNx/zvG3kJ0SuU7Iswd74KNaPQE20f05bkm2/wH4+8FOeL9WfwP9gftL3Bpl7Bv5eXjdr/+FerH7K0xJ/rft2Hfzs/APfyt5fbhrK/4n5T+2Y5cdSsF8fzt5fZBrF+5MsQ+sVhe4MN/XWV7v7TqAu6OM/aKgGvN9PeT1zjiKe1Pfsh27taAU831nyOt5KMbvUu3fYPYePoH52/VAMBVV+H3qO7ZjC44UYP4O/c1IuD/Vfi+mj/yDMX/ncOP2tJSFSHcFA1+zdbV9MX/3Bcbt21MWoY1L/wUb7qvaXpi/J7ASMIDJyUvQ1X3UcuyO2s6Yv/cy4/Y1yctwhvug5dgDWhvM3/cj4/aVyR9hiHu35djjWhbm759o3B7r/S8Ge7ZYjq3QkjHh4M+N2xd6V2GwJzRQMbsi/3rjN9DZ3tUY5FlnO/aqQ1ejAnpgcWbSWgxK+q/t2KVrP8d+Tf85fR83YuQZHfG3q4cZj4sg4sPNh3EwUAd0dvdWRoAR7rR2mVi2TZ9quups/ZeplQGd9F/2XVulIzmwiapVcCOm1ArLa+RGqwM7659jJjZEvfeyPnjwg20y+BGuH9kDZwQCrK556bhv8kRcM68TLqj+BD/xrMRw9zakoEoGpuLjxbX7oQUS6Wd7P8I5Ub6uN+wfj/LA1/XhpI8wPMrX9fb9o3AM+jH/IWklRiQtsR1794HzcEDT3xt+m7QKI5Lsf35mHhiGbzT9a/trzyf4f17796cHDwzABk3/2brFsxaTvAtsxz5x8HSsCWwi+/88n+KmKGP/dqAbloe899qPfflAe7wXeO8dJ9977cf++2AuXvdnmN577ccuPJiGv/tzTe+99mNXHHRjnl8PGge6duOuKGPX5fswz6cHjae5DuLuKGO/zi/DPF8303uv/did+ccxz6cHja3ke2+Ur4P/Aszbpv+MpaAav48ydqH/XMwLvE8L90UZu8I/GPO+Cb5P353yBjKivPfO23GOcftXKW+irc1779e1PTFv53nG7ZuS30ZOTfNMObs0YyMc58jPz0fnzp2xevVqjBypR9PCtGnTsHLlSqxbt65eGaCuXbtiyRO3IiMt8q+/LW1/gJNp3eX1dmXbcfpxi79YAra3Ho1jGb3l9Tblu9Dn2FLbsTvyLsKRzL7yequKveh37EPbsbtzz0NBlv7Nl12ZjwFH37Md+23OcORnD5bXM6uOYOAR++DjQPZQ7M/Rfwml1ZzEkAL7N7n8rIHYm6sHFCm+EpxVYP8DcDijL/a00jNLSf5KnH3oFduxR9NPx648/a8bd60P5+TbLyx3Iq0HdrQOBh/DD4i/hCKDVqEwtQu2t/m+cfucg/+AC/obb7iS5PbY2naccfusQ/+S2QsrZd7W2Nwu+JfQ4IJ/I8UfbMk2E39tb2wf/EtowJF3kFZjnS0Sf8V/1SEY1PQ/uhAZ1cctx4qsy5cdf2bc7nNsCbKqD1uOrXV58EXHYAAkvn9zbDJLwucdrzEyEb1OfoxWFftsx37Z4WeoDWQ1ehSuQevyPWHHmYzKpBxUeHNwKHMgajxpcsqqf6dsXNC7TUT3lqj9+Wj7ERSV1+Ca4d2MNXrCFVfW4KNtR/C93m0s1+0xe/PLA+jZJhNDuuq/xKyIrTJEAbYIlHJM9UXhtuQXY+3u4xjQOQfn9GgVsbCiWChx5TdHcFCsPu2vQXZVgfxaZ9Qcw5a2441xZxxfijbl1sGtsK7zDUYGSLxe4n1H5H8EV+DSPLbGky6v9zrxMTqVbrR93s87XotKrx4s9Ti5Bl1KvrQd+0WHq1Ce3Fpe71r0GboXBfZqs/BV+4koSdGD287FX6JnoXX2XdjY7goUpep/dHUo2YTeJ0ODO/PZbWk7Tv7MC+1Kt+GME9Hfe49mnG689/Y9Zp8t2pF3MQ5n9pPXxff3mUftM/W7W50v3/uEnMr8qO+n3+aOwIFsPVOfWXUYQw7b/6G6P3sY9ubqmfr0mhMYesg+q38wazD2tNJ/8af4inFOvv37aUFmf+zMu0he9/orMPzgfNuxRzLOwDetLzXee8878Jzt2ONpvbC1rZ4tEs7f94zt2JOpXUPeI8/b/zzcNhmg4pSO+Nr0Hjn8wHx4bTJApcltsaHDT4zb5+S/DF/JMXz/18+iqKgI2dnWGeMWGwCJKTBR7/PGG2/gyiuDX8RJkyahsLAQb79t/81qDoBycnKa/AtIREREjae5fn87sgg6OTkZw4YNw7Jly4z7amtr5W1zRoiIiIioxdQACaIFXmR8zj77bJx77rl4/PHHUVZWJrvCiIiIiFpkAHTVVVfh6NGjmDlzJgoKCjBkyBAsWrQoojCaiIiIqEXUADUG1gARERHFn+JErgEiIiIiakoMgIiIiCjhMAAiIiKihMMAiIiIiBIOAyAiIiJKOAyAiIiIKOEwACIiIqKEwwCIiIiIEg4DICIiIko4jt0K47tSC1yLFSWJiIgoPhQHfm839UYVLTYAOn78uLzs2rVrrA+FiIiITuH3uNgSo6m02AAoLy9PXu7bt69Jv4BOjJxF0Ld///4m3UPFaXjefL0TAb/P+X2eCIqKitCtWzfj93hTabEBkNutlzeJ4CeRAgFFnDPPO3Hw9U4sfL0TS6K+3u7A7/Eme/4mfXYiIiIiB2IARERERAmnxQZAKSkpuO++++RlIuF58/VOBPw+5/d5IuD3eUqTfn1dWlP3mRERERE5TIvNABERERHZYQBERERECYcBEBERESUcBkBERESUcOI2APrTn/6E8847D+np6cjNza3X54h675kzZ6Jjx45IS0vD6NGjsWPHjpAxJ06cwLXXXisXnRLPO3nyZJSWlsIpGnp83377LVwul+XH66+/boyzevy1116DU5zK63LxxRdHnNOtt94aMkasFD5+/Hj5fdSuXTvcc8898Pl8iNfzFuPvuOMO9OnTR36Pi9VU/+d//keurGrmtNd7zpw56NGjB1JTUzF8+HB8+umnUceL792+ffvK8QMHDsTChQsb/LPuBA057+eeew4XXHABWrVqJT/EOYWPv+GGGyJe18suuwzxfN7z58+POCfxeS399bZ6/xIf4v0qnl7v//73v7j88svRqVMneXxvvfVWnZ+zYsUKDB06VHbB9e7dW34PfNf3DEtanJo5c6b26KOPalOnTtVycnLq9TkPPvigHPvWW29pX331lfajH/1I69mzp1ZRUWGMueyyy7TBgwdra9eu1T7++GOtd+/e2s9//nPNKRp6fD6fTzt06FDIxx/+8ActMzNTKykpMcaJb4V58+aFjDN/XWLtVF6Xiy66SLvllltCzqmoqCjkazNgwABt9OjR2pdffqktXLhQa9OmjTZ9+nQtXs9748aN2oQJE7R33nlH27lzp7Zs2TLt9NNP1yZOnBgyzkmv92uvvaYlJydrL774orZ582b5muXm5mqHDx+2HP/JJ59oHo9He+ihh7QtW7ZoM2bM0Lxerzz3hvysx1pDz/uaa67R5syZI79Xt27dqt1www3yHA8cOGCMmTRpkvyeMb+uJ06c0Jykoectvk+zs7NDzqmgoCBkTEt8vY8fPx5yzps2bZLf9+LrEU+v98KFC7X//d//1f7zn//I950333wz6vjdu3dr6enp8ne7+Pl+8skn5XkvWrTolL+WduI2AFLEN0N9AqDa2lqtQ4cO2sMPP2zcV1hYqKWkpGj//Oc/5W3xxRYv0GeffWaM+eCDDzSXy6UdPHhQi7XGOr4hQ4ZoN910U8h99fnGjLfzFgHQr3/966g/mG63O+TN9JlnnpFvtlVVVVpLeb0XLFgg3yxqamoc+Xqfe+652pQpU4zbfr9f69SpkzZ79mzL8T/72c+08ePHh9w3fPhw7Ze//GW9f9bj8bzDiQA+KytLe+mll0J+IV5xxRWakzX0vOt6j0+U1/uxxx6Tr3dpaWlcvd5m9XnfmTZtmnbmmWeG3HfVVVdpY8eObbSvpRK3U2ANtWfPHhQUFMjUqCL2CROpszVr1sjb4lJMM5x99tnGGDFe7Eeybt06xFpjHN/69euxYcMGOZUSbsqUKWjTpg3OPfdcvPjiizKt7ATf5bxfeeUVeU4DBgzA9OnTUV5eHvK8Yvqkffv2xn1jx46VG05u3rwZsdZY349i+ktMoSUlJTnu9a6urpbfk+afS3F+4rb6uQwn7jePV6+bGl+fn/VYO5XzDie+l2tqaiI2jBTTB2I6V0yD3nbbbXJHbac41fMW077du3eXGz1fccUVIT+fifJ6v/DCC7j66quRkZERN6/3qajr57sxvpYtfjPUcOIHRDD/slO31WPiUnwjmYlfGuINRo2JpcY4PvFD1K9fP1k/ZTZr1ixccsklshZm8eLF+NWvfiXfdET9SLye9zXXXCPfNMXc89dff417770X27dvx3/+8x/jea2+H9RjLeH1PnbsGB544AH84he/cOTrLY7P7/dbvg7btm2z/By71838c6zusxsTa6dy3uHE97P43jb/IhD1HxMmTEDPnj2xa9cu/O53v8O4cePkLwaPx4N4PG/xi10E6IMGDZLB/COPPCLfv0QQ1KVLl4R4vUV9y6ZNm+T7t5nTX+9TYffzLf4wraiowMmTJ7/zz44jA6Df/va3+POf/xx1zNatW2XxY0tS3/P+rsQ3z6uvvorf//73EY+Z7zvrrLNQVlaGhx9+uEl/ITb1eZt/6YtMjyiQvPTSS+UbxWmnnYaW/nqLNwxRMNm/f3/cf//9MX+9qfE8+OCDsmhd/PVvLggWGQLz97wIGsT3uhgnvvfj0ciRI+WHIoIf8Ufcs88+K4P7RCACH/F6imytWUt8vZuTowKg3/zmN7KqPZpevXqd0nN36NBBXh4+fFj+IlTE7SFDhhhjjhw5EvJ5oiNIdNaoz4/leX/X43vjjTdk2vz666+vc6xIH4s3l6qqqibbT625ztt8TsLOnTvlm4T43PDOAfH9IMT7611SUiL/OszKysKbb74Jr9cb89fbipiCE3+pqq+7Im7bnaO4P9r4+vysx9qpnLciMiAiAFq6dKn8hVfX95H4v8T3vBN+IX6X81bE97II2sU5JcLrLf44EcGuyNrWxWmv96mw+/kW0/iiw098Hb/r95BBS7Ai6EceecS4T3QEWRVBf/7558aYDz/80HFF0Kd6fKIoOLwbyM4f//hHrVWrVpoTNNbrsmrVKvk8okvEXARt7hx49tlnZRF0ZWWlFq/nLb6vR4wYIV/vsrIyx7/eoqDx9ttvDylo7Ny5c9Qi6B/+8Ich940cOTKiCDraz7oTNPS8hT//+c/y+3PNmjX1+j/2798vv1/efvttLZ7PO7z4u0+fPtpdd93V4l9v9TtOnMuxY8fi8vU2q28RtOjONROdr+FF0N/le8g4Hi1O7d27V7aDqpZucV18mFu7xQ+JaL0zt0qKVjnxzfH111/L6nmrNvizzjpLW7dunfyFKVqIndYGH+34REusOG/xuNmOHTvkD4boIgonWqafe+452UYsxj399NOyDVEsNRCv5y1awGfNmiWDhz179sjXvFevXtqFF14Y0QY/ZswYbcOGDbLNsm3bto5rg2/IeYs3ftERNXDgQPk1MLfHivN14ustWlrFG/z8+fNl0PeLX/xC/pyq7rzrrrtO++1vfxvSBp+UlCR/4Yl28Pvuu8+yDb6un/VYa+h5i3MS3XxvvPFGyOuq3vPE5d133y2DI/E9v3TpUm3o0KHye8YJAf2pnrd4jxeB/65du7T169drV199tZaamirbn1vy662cf/75sgsqXLy83iUlJcbvZxEAieVrxHXxO1wQ5yzOPbwN/p577pE/32LpB6s2+GhfyxYfAIn2P/HFDP/46KOPItY6UcRfCr///e+19u3byy/epZdeqm3fvj1i7QXxC0YEVeIvrRtvvDEkqIq1uo5P/CCEfx0E8Uu9a9euMlIOJ4Ii0RovnjMjI0OuOzN37lzLsfFy3vv27ZPBTl5ennytxfo54gfKvA6Q8O2332rjxo3T0tLS5BpAv/nNb0LaxePtvMWl1c+F+BBjnfp6i7U+unXrJn/Bi7/uxLpHishkiZ/38Nb+M844Q44XLbPvv/9+yOP1+Vl3goacd/fu3S1fVxEACuXl5TKYF0G8CAjFeLE+SkN/KTjtvO+8805jrHg9f/CDH2hffPFFi3+9hW3btsnXePHixRHPFS+v90c270nqXMWlOPfwzxHvUeLrJP5wNf8er8/Xsr5c4p9TmqgjIiIiilMJsw4QERERkcIAiIiIiBIOAyAiIiJKOAyAiIiIKOEwACIiIqKEwwCIiIiIEg4DICIiIko4DICIiIgo4TAAIiIiooTDAIiIiIgSDgMgIiIiSjgMgIiIiAiJ5v8DJrYQTHAP5q4AAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(shifted_ckvec, np.fft.fftshift(np.abs(phi)**2))\n", "psig = 1/xsig/klaser\n", "phi_analytic = np.exp(-np.square(shifted_ckvec)/(2*psig**2))\n", "plt.plot(shifted_ckvec, np.max(np.abs(phi)**2)*np.square(phi_analytic), linestyle='--') # Factor of 4 \n", "plt.xlim([-1,1])\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "d491bfae-7599-476b-b19d-d3f1d3f199ae", "metadata": {}, "source": [ "### Simulating a distribution of momenta to recover the splitstep result" ] }, { "cell_type": "markdown", "id": "ccab9960-4676-4d3d-9fb3-0f2deb9fc7ac", "metadata": {}, "source": [ "Now I want to perform the planewave simulation at several different detuning values to see if I can get agreement with the splitstep method." ] }, { "cell_type": "code", "execution_count": 5, "id": "67fe9c2c-10d4-4841-b949-c6fd95373a8f", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from mwave.integrate import make_kvec, make_phi\n", "from mwave.numeric import NumericBraggInterferometer\n", "from scipy.integrate import solve_ivp\n", "import numpy as np\n", "from numba import jit, float64\n", "\n", "# Set interferometer parameters\n", "nbragg = 4\n", "T = 10\n", "Tp = 5\n", " \n", "# Initialize kvec and tree\n", "ifr = NumericBraggInterferometer(-2*nbragg, 4*nbragg, distance=4) # Should pass in wavefunction initialization functions\n", "\n", "ifr.split(nbragg) # Should pass in the relevant function here as well! Function will inspect input function arguments to ensure correct number and naming\n", "ifr.propagate(T)\n", "ifr.split(nbragg)\n", "ifr.propagate(Tp)\n", "ifr.split([3*nbragg, -nbragg])\n", "ifr.propagate(T)\n", "ifr.split([3*nbragg, -nbragg])" ] }, { "cell_type": "markdown", "id": "8b32e016-ab9c-48d4-81c3-31a7e7512152", "metadata": {}, "source": [ "We define a batched Bragg propagator built on top of :py:func:`mwave.integrate.propagate`. Instead of looping over detunings and calling ``solve_ivp`` once per atom, we hand the whole batch to ``propagate`` in one shot.\n", "\n", "The Bragg pulse for atom :math:`i` needs the laser phase :math:`\\delta_i t` (the laser drive phase already accumulated by the absolute pulse start time :math:`t`). Per-atom phase offsets are not directly supported by ``propagate``'s shared ``phase`` callable, so we absorb them via a change of frame: define :math:`\\psi_i[j] = e^{i\\delta_i t \\cdot j}\\,\\phi_i[j]`, where :math:`j` is the momentum-grid index. The transformed wavefunction satisfies the standard equation with ``phase=0``. After integration we recover :math:`\\phi_i[j] = e^{-i\\delta_i t \\cdot j}\\,\\psi_i[j]`. This is exact (no global-phase ambiguity), and reduces the wall-clock cost of each beamsplitter from ``len(deltas)`` scipy calls to a single batched ``propagate`` call." ] }, { "cell_type": "code", "execution_count": 6, "id": "f343ce53-5fcd-4059-8fe2-240581d1c741", "metadata": {}, "outputs": [], "source": [ "from numba import jit, float64\n", "import mwave.integrate as mi\n", "from mwave.integrate import omega_fnc_gaussian, phase_fnc_constant, make_phi\n", "\n", "@jit(float64(float64, float64[:]))\n", "def multi_omega_fnc(t, args):\n", " omega, sigma, t0, mod_freq, mod_phase = args\n", " return 2*np.cos(mod_freq*t + mod_phase)*omega*np.exp(-np.square(t-t0)/(2*(sigma**2)))\n", "\n", "def gbragg_batched(kvec, phi0_single, tfinal, deltas, omega, sigma, t_off,\n", " dphase=0.0, omega_mod=None, mod_phase=0.0,\n", " atol=1e-10, rtol=1e-10, tol=1e-10):\n", " \"\"\"Batched Bragg propagator across many detunings.\n", "\n", " The Hamiltonian seen by atom ``i`` has detuning ``deltas[i]`` and a\n", " constant laser phase ``deltas[i] * t_off + dphase``. We absorb the per-atom\n", " constant phase via a momentum-grid phase ramp before/after integration so\n", " that the underlying call to :py:func:`mwave.integrate.propagate` only needs\n", " a shared ``dphase``.\n", " \"\"\"\n", " deltas = np.asarray(deltas, dtype=np.float64)\n", " natoms = len(deltas)\n", "\n", " # Frame phase: V[i, j] = exp(i * delta_i * t_off * j) where j is grid index\n", " j_idx = np.arange(len(kvec))\n", " frame_phase = np.exp(1j * np.outer(deltas * t_off, j_idx))\n", "\n", " # Build batched initial state\n", " phi0_b = np.broadcast_to(phi0_single, (natoms, len(kvec))).astype(np.complex128).copy()\n", " phi0_b *= frame_phase\n", "\n", " if omega_mod is None:\n", " result = mi.propagate(\n", " kvec, phi0_b, np.float64(tfinal), deltas,\n", " omega_fnc_gaussian, np.array([omega, sigma, tfinal/2]),\n", " phase_fnc_constant, np.array([dphase]),\n", " backend='numba', tol=tol,\n", " )\n", " else:\n", " result = mi.propagate(\n", " kvec, phi0_b, np.float64(tfinal), deltas,\n", " multi_omega_fnc, np.array([omega, sigma, tfinal/2, omega_mod, mod_phase]),\n", " phase_fnc_constant, np.array([dphase]),\n", " backend='numba', tol=tol,\n", " )\n", "\n", " # Undo the frame transformation\n", " return result.phi_final * np.conj(frame_phase)\n", "\n", "\n", "def deltalookup(v, n_bragg):\n", " return 4*n_bragg + 4*(v/0.0035) # The modification to delta is 4 times the velocity defined in units of recoil velocities\n", "\n", "def cpops(target_phase=np.pi/2, dphase=0.0, deltavals=np.array([0.0])):\n", " bs_lookup_dict = {}\n", " \n", " omega_m = 8*nbragg-target_phase/(2*T*nbragg)\n", " \n", " def calc_bs(deltashift, k_init, k_final, klattice, t, x, idx):\n", " \n", " if idx in bs_lookup_dict:\n", " if k_init in bs_lookup_dict[idx]:\n", " kf_idx = np.argmin(np.abs(ifr.kvec - k_final))\n", " return bs_lookup_dict[idx][k_init][:,kf_idx]\n", " \n", " sigma = 0.259\n", " omega = 19.4\n", "\n", " deltas = 4*nbragg + deltashift\n", "\n", " # All beamsplitters share `t` as the absolute pulse-start time. The\n", " # per-atom constant phase is `deltas[i] * t` (with an extra `dphase`\n", " # for the third pulse), which gbragg_batched handles via the\n", " # momentum-grid frame transformation.\n", " extra_dphase = dphase if idx == 3 else 0.0\n", "\n", " if idx == 1 or idx == 2:\n", " soly = gbragg_batched(\n", " ifr.kvec, make_phi(ifr.kvec, k_init/2), 6*sigma,\n", " deltas, omega, sigma, t_off=t, dphase=extra_dphase,\n", " )\n", " elif idx == 3 or idx == 4:\n", " soly = gbragg_batched(\n", " ifr.kvec, make_phi(ifr.kvec, k_init/2), 6*sigma,\n", " deltas, omega, sigma, t_off=t, dphase=extra_dphase,\n", " omega_mod=omega_m, mod_phase=omega_m*t,\n", " )\n", "\n", " if idx not in bs_lookup_dict:\n", " bs_lookup_dict[idx] = {}\n", " if k_init not in bs_lookup_dict[idx]:\n", " bs_lookup_dict[idx][k_init] = soly\n", " else:\n", " raise RuntimeError('Array should not have been created but it was!')\n", " \n", " kf_idx = np.argmin(np.abs(ifr.kvec - k_final))\n", " return soly[:,kf_idx]\n", " \n", " propfnc = lambda deltashift, t, k: np.exp(-1j*t*k**2)\n", " fnclst = [lambda deltashift, k_init, k_final, klattice, t, x: calc_bs(deltashift, k_init, k_final, klattice, t, x, 1), propfnc, lambda deltashift, k_init, k_final, klattice, t, x: calc_bs(deltashift, k_init, k_final, klattice, t, x, 2), propfnc, lambda deltashift, k_init, k_final, klattice, t, x: calc_bs(deltashift, k_init, k_final, klattice, t, x, 3), propfnc, lambda deltashift, k_init, k_final, klattice, t, x: calc_bs(deltashift, k_init, k_final, klattice, t, x, 4)]\n", " fnclst = fnclst\n", "\n", " # # Do some filtering of the wavefunctions to only keep the ones of interest, i.e. those at the right momentum state\n", " # Instead of calling in parallel have one call for each momentum state\n", " ifr.set_operation_funcs(fnclst)\n", " popfnc = ifr.get_population_func([4*nbragg, 2*nbragg, 0*nbragg, -2*nbragg], lambda x: np.zeros_like(deltavals), lambda x: np.ones_like(deltavals,dtype=np.complex128), lambda x: np.zeros_like(deltavals,dtype=np.complex128))\n", " \n", " return popfnc(4*nbragg, [deltavals]), popfnc(2*nbragg, [deltavals]), popfnc(0*nbragg, [deltavals]), popfnc(-2*nbragg, [deltavals])" ] }, { "cell_type": "markdown", "id": "e5b199e5-e5d1-4f63-b99c-6f9a08ac21f2", "metadata": {}, "source": [ "Take the momentum wavefunction width we used in the splitstep simulation. Lets make sure to compute the planewave simulation with a range of detunings that can capture this momentum width." ] }, { "cell_type": "code", "execution_count": 7, "id": "16bd3432-dd6a-4972-9388-8f525c0002ad", "metadata": {}, "outputs": [], "source": [ "psig = 0.04714045207910317 # Comment this line out if you have computed psig directly from xsig in the previous section\n", "pw_kvals = np.linspace(-5*psig, 5*psig, 50)\n", "pw_pops = cpops(target_phase=0*np.pi/2, dphase=0.0, deltavals=pw_kvals*4)" ] }, { "cell_type": "markdown", "id": "f37b6227-7f71-4ef3-b76f-a78abd502b35", "metadata": {}, "source": [ "Now that we have computed our planewave population at each detuning value, we can compute final populations weighted by the wavefunction value at each detuning." ] }, { "cell_type": "code", "execution_count": 8, "id": "737a6576-d381-45e5-99c4-381d7c3d1a45", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/var/folders/d7/6t4hylss1p9_9vftvfm4fbm80000gn/T/ipykernel_41528/2428742450.py:5: DeprecationWarning: `trapz` is deprecated. Use `trapezoid` instead, or one of the numerical integration functions in `scipy.integrate`.\n", " wf_phi /= np.sqrt(np.trapz(wf_phi*np.conjugate(wf_phi), wf_kvals))\n", "/var/folders/d7/6t4hylss1p9_9vftvfm4fbm80000gn/T/ipykernel_41528/2428742450.py:16: DeprecationWarning: `trapz` is deprecated. Use `trapezoid` instead, or one of the numerical integration functions in `scipy.integrate`.\n", " pops_weighted_planewave = np.array([np.trapz(wf_pop*pw_pops_interp[:,i], wf_kvals) for i in range(4)])\n" ] }, { "data": { "text/plain": [ "array([0.36490007, 0.12204806, 0.36423925, 0.12212211])" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "wf_kvals = np.linspace(-5*psig,5*psig)\n", "psig2 = psig/np.sqrt(2)\n", "wf_phi = np.exp(-np.square(wf_kvals)/(2*(psig2)**2))\n", "# Normalize\n", "wf_phi /= np.sqrt(np.trapz(wf_phi*np.conjugate(wf_phi), wf_kvals))\n", "\n", "# Compute population\n", "wf_pop = np.abs(wf_phi)**2\n", "\n", "# Interpolate pops\n", "from scipy.interpolate import interp1d\n", "pw_pops_fnc = interp1d(pw_kvals, pw_pops, kind='cubic')\n", "pw_pops_interp = pw_pops_fnc(wf_kvals).T\n", "\n", "# Compute weighted populations\n", "pops_weighted_planewave = np.array([np.trapz(wf_pop*pw_pops_interp[:,i], wf_kvals) for i in range(4)])\n", "pops_weighted_planewave" ] }, { "cell_type": "markdown", "id": "06b176e3-67d6-473f-b9d2-0922aaff8152", "metadata": {}, "source": [ "Lets compare this to the splitstep method, where we get the following population outputs" ] }, { "cell_type": "code", "execution_count": 9, "id": "67d657a5-b634-4d63-a697-7917337a9a63", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0.36130238, 0.12050147, 0.3599777 , 0.12058718])" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pops_splitstep = np.array((0.3613023794840803, 0.1205014725343373, 0.35997770257669187, 0.12058717738980357))\n", "pops_splitstep" ] }, { "cell_type": "markdown", "id": "224757fb-ab60-4271-9c81-554dd049b20d", "metadata": {}, "source": [ "This don't seem to agree at first glance beyond 0.002. However if we normalize the populations to 1 then we see" ] }, { "cell_type": "code", "execution_count": 10, "id": "44ccfb22-63a9-403a-9fa9-167c1eec390b", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0.37490651 0.12539491 0.37422758 0.125471 ]\n", "[0.3754303 0.12521341 0.37405382 0.12530247]\n" ] } ], "source": [ "print(pops_weighted_planewave/np.sum(pops_weighted_planewave))\n", "print(pops_splitstep/np.sum(pops_splitstep))" ] }, { "cell_type": "markdown", "id": "e446aeae-98e9-464b-8484-8134f36f48ed", "metadata": {}, "source": [ "Now we are getting agreement to roughly 0.0005 or better.\n", "\n", "Note that this agreement is noticably better if we do not simulate the final interfering Bragg pulse. In this case we get agreement to 0.0001. I think this could be because I am not treating the free evolution phase of the shifted momentum states properly?\n", "\n", "I am including the population results for the case where I compare the first three Bragg pulses below:" ] }, { "cell_type": "code", "execution_count": 11, "id": "79a74768-0a66-444b-a3f5-7f08d7d64f17", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 1.3680e-05, -1.0859e-04, 1.7970e-05, 7.6940e-05])" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "p3_pw = np.array([0.24660472, 0.25327332, 0.25324507, 0.24687689])\n", "p3_ss = np.array([0.24659104, 0.25338191, 0.2532271, 0.24679995])\n", "p3_pw - p3_ss" ] }, { "cell_type": "markdown", "id": "8a78c92f-e701-4d46-aaae-5c802b29b0b3", "metadata": {}, "source": [ "I think this is good enought to proceed with for now." ] } ], "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.12" } }, "nbformat": 4, "nbformat_minor": 5 }