{ "cells": [ { "cell_type": "markdown", "id": "cac2fea9", "metadata": {}, "source": [ "# Comparing the splitstep and plane-wave methods" ] }, { "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:03<00:00, 992.25it/s] \n", "100%|██████████| 19999/19999 [00:19<00:00, 1038.26it/s]\n", "100%|██████████| 10000/10000 [00:09<00:00, 1033.70it/s]\n", "100%|██████████| 19999/19999 [00:18<00:00, 1084.60it/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": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGfCAYAAACneiONAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAVbNJREFUeJzt3XtcVHX+P/DXDHdU8IICGiYVheaFREXM1tpYsWxXNjNzK82fX93cbHWpXHVJSy1cyzLTImt1c8tw2cpKDTPMzCQM8H6/BoLDRYSBUa5zfn+Yo8PMwAxzObfX8/HgoZz5zJn3mRnmvOZzPudzNIIgCCAiIiKSOa3YBRARERG5AkMNERERKQJDDRERESkCQw0REREpAkMNERERKQJDDRERESkCQw0REREpAkMNERERKQJDDRERESkCQw0REREpgndb7rRq1Sq89tpr0Ol0GDBgAN5++20MGTLEZvuMjAy8+OKLOHfuHKKiovDPf/4TDz74oOn2zz77DGlpacjLy0NFRQX27t2LmJgY0+0VFRVYsGABvvnmGxQUFKBr165ISkrCokWLEBwcbFfNRqMRxcXF6NChAzQaTVs2m4iIiDxMEARUV1eje/fu0Gpb6YsRHJSeni74+voKa9asEQ4fPixMnTpV6Nixo1BSUmK1/Y8//ih4eXkJS5cuFY4cOSKkpKQIPj4+wsGDB01t1q1bJ7z88svC+++/LwAQ9u7da7aOgwcPCg8//LDw5ZdfCqdOnRKysrKEqKgoYezYsXbXXVhYKADgD3/4wx/+8Ic/MvwpLCxsdV+vEQTHLmgZFxeHwYMHY+XKlQCu9oBERETg2WefxZw5cyzajx8/HgaDAZs2bTItGzp0KGJiYpCWlmbW9ty5c4iMjLToqbEmIyMDTzzxBAwGA7y9W+9wqqqqQseOHVFYWIigoCA7tpSIiIjEptfrERERgcrKylaPzjh0+Km+vh55eXmYO3euaZlWq0VCQgKys7Ot3ic7OxvJyclmyxITE7Fx40ZHHtpCVVUVgoKCbAaauro61NXVmX6vrq4GAAQFBTHUEBERyYw9Q0ccGihcXl6OpqYmhIaGmi0PDQ2FTqezeh+dTudQe3vrWLRoEaZNm2azTWpqKoKDg00/ERERbX48IiIikj7Znf2k1+sxevRo9OnTBy+99JLNdnPnzkVVVZXpp7Cw0HNFEhERkcc5dPgpJCQEXl5eKCkpMVteUlKCsLAwq/cJCwtzqH1LqqurMWrUKHTo0AGff/45fHx8bLb18/ODn5+fw49BRERE8uRQT42vry9iY2ORlZVlWmY0GpGVlYX4+Hir94mPjzdrDwDbtm2z2d4WvV6PkSNHwtfXF19++SX8/f0duj8REREpm8Pz1CQnJ2PSpEkYNGgQhgwZguXLl8NgMGDy5MkAgIkTJ6JHjx5ITU0FAMycORMjRozAsmXLMHr0aKSnpyM3NxerV682rbOiogIFBQUoLi4GABw/fhzA1V6esLAwU6C5fPkyPvroI+j1euj1egBA165d4eXl5dyzQERERLLncKgZP348ysrKMH/+fOh0OsTExCAzM9M0GLigoMBscpxhw4Zh/fr1SElJwbx58xAVFYWNGzeib9++pjZffvmlKRQBwGOPPQYAWLBgAV566SXk5+cjJycHAHDbbbeZ1XP27Fn06tXL0c0gIiIihXF4nhq50uv1CA4ONp0KTkRERNLnyP5bdmc/EREREVnDUENERESKwFBDREREisBQQ0RERIrAUENERESKwFBDRNQGP5+rwEc//QKVnEBKJAsOz1NDRETAuLRsAEBkSDvcfVuIyNUQEcCeGiJygdqGJtX2WJwtN4hdAhH9iqGGiJxSWHEZ0S9m4tlP9opdChGpHEMNETnlPz/9AgDYdOCCyJUQkdox1BAROUGdB92IpImhhmSptqEJj7y7G73mbMbzGftVO56DiIiuY6ghWfpf3nnk/nLJ9H8O1pSG3HMVYpdARCrGUEOyVNdoNPu90cieGil45NfTnFWFvYREksFQQ7Kk1YhdAanZm9tOmP6/LvsX1DU2iVgNEV3DUEOyxExDYqkw1OOtrJOm30+W1mDVd6dFrIiIrmGoIVnSaBhrSBz1zQ59AkDOmYsiVEJEzTHUkCwx00iH2s48s/beU9lTQCRZDDUkS+ypISKi5hhqSJYYaUhKBE7BRyQJDDWkCOtzCsQugVTCWqDmjAJE0sBQQ7JTebkeKRsPmS379+5zqLrSIFJFpHZqG1dEJFUMNSQ7Czcdsbr83z+e82whBECF45usDRT2fBVEZAVDDcnO6TLrl0RoNFqeakvkCeyoIZIGhhoiIifx8BORNDDUkOyo7GCH5Klth67hO5BIshhqSDFUtm8lkahtCBGRnDDUEBE5iXmaSBoYaoiIHGCto4a9hETSwFBDRE4pqLgsdglERAAYakiGOKZBOj7ZU4Cth0vELkN0vEwCkTQw1BBRm63IOil2CR6nuskGiWSEoYaIiIgUgaGGiMhJHChMJA0MNSQ7tjr/Oa6BPMHaZIMMNUTSwFBDROQkZhoiaWCoIcXg9PXkbk1GAe/uOG2xXG2XiiCSKoYakh1bZ5/U1DV6uBJSm0/zzuODXWctlh/TVWPj3iIRKiKiGzHUkGL8e/c51DcaxS5DNQorLuNCVa3YZXjU6fIam7fN2rDPc4UQkVUMNaQoFYZ6sUtQjQfe+kHsEjzOW8tDnERSxlBDRG2ixsN9Xlp+ZBJJGf9CSXZa+q7M07rJndhPQyRtDDVE5FIvZOzn2UBEJAqGGlIU7kvFl5F3HtmnL4pdhlvwsk9E0sZQQ7KT+8slsUugVlyubxK7BLfgXEhE0sZQQ0RERIrAUEOKwqNPRETq1aZQs2rVKvTq1Qv+/v6Ii4vDnj17WmyfkZGB6Oho+Pv7o1+/ftiyZYvZ7Z999hlGjhyJLl26QKPRYN++fRbrqK2txTPPPIMuXbqgffv2GDt2LEpKStpSPhG5mVLDJcfUEEmbw6Fmw4YNSE5OxoIFC5Cfn48BAwYgMTERpaWlVtvv3r0bEyZMwJQpU7B3714kJSUhKSkJhw4dMrUxGAwYPnw4/vnPf9p83L/97W/46quvkJGRge+//x7FxcV4+OGHHS2fZE5f29Di7bnnKjxUibqV19SJXQIRkQWN4OC5l3FxcRg8eDBWrlwJADAajYiIiMCzzz6LOXPmWLQfP348DAYDNm3aZFo2dOhQxMTEIC0tzaztuXPnEBkZib179yImJsa0vKqqCl27dsX69evxyCOPAACOHTuG3r17Izs7G0OHDm21br1ej+DgYFRVVSEoKMiRTSYJefJfOfjhZLnN2328NDj5yoMerEides3Z3OLt708chN/1CfVQNZ6zIusk3th2wubt55aM9mA1ROrgyP7boZ6a+vp65OXlISEh4foKtFokJCQgOzvb6n2ys7PN2gNAYmKizfbW5OXloaGhwWw90dHR6Nmzp8311NXVQa/Xm/2Q/LUUaACenULuxXcXkbQ5FGrKy8vR1NSE0FDzb2ChoaHQ6XRW76PT6Rxqb2sdvr6+6Nixo93rSU1NRXBwsOknIiLC7scjIucodfI9jqkhkjbFnv00d+5cVFVVmX4KCwvFLok8gJdJICJSL29HGoeEhMDLy8virKOSkhKEhYVZvU9YWJhD7W2to76+HpWVlWa9NS2tx8/PD35+fnY/BimDQjsIiIjIDg711Pj6+iI2NhZZWVmmZUajEVlZWYiPj7d6n/j4eLP2ALBt2zab7a2JjY2Fj4+P2XqOHz+OgoICh9ZDyqe0TLPjeCniU7Pw46mWxxJJjUahx2mUul1ESuFQTw0AJCcnY9KkSRg0aBCGDBmC5cuXw2AwYPLkyQCAiRMnokePHkhNTQUAzJw5EyNGjMCyZcswevRopKenIzc3F6tXrzats6KiAgUFBSguLgZwNbAAV3towsLCEBwcjClTpiA5ORmdO3dGUFAQnn32WcTHx9t15hOph9LGcjy19mcAwOMf5MjqzBqlvQ5EJA8Oh5rx48ejrKwM8+fPh06nQ0xMDDIzM02DgQsKCqDVXu8AGjZsGNavX4+UlBTMmzcPUVFR2LhxI/r27Wtq8+WXX5pCEQA89thjAIAFCxbgpZdeAgC8+eab0Gq1GDt2LOrq6pCYmIh33nmnTRtNysVdKRGRejk8T41ccZ4aZWhtfhRAWXOF3Li9Utqu1l6H1U/GYuSd9o+bk4tV353Ca1uP27xdSq8RkVK4bZ4akq5SfS32FVaKXQaRonFIDZG0OXz4iaRpyKtXB1F/NWM4+t0ULHI1pHaq6P4lIslhT43C/MxrHxERkUox1CiASoZFEYnOaGz5b40X+iQSF0ONzAmCgIlr9ohdBpHi7T5djte/sX0xSwCoutLyVeSJyL0YamSu8nJDqxd5JCLn/W3DvlbbcBwxkbgYahSGB6KU65eLBrFLoFaszykQuwQiVWOoIZKJEa/tELsEVdPY0Q/zwa6zHqiEiGxhqCHFaWgyil2CSxw8XyV2CXQDzlFDJH0MNaQ4n+cXiV2CS/x+5S6xSyAikhWGGlKcS5frxS5B9ZQ4ywA7aoikj6GGSEbKqjkPilg0PP5EJHkMNTLHz1lLCuwkMBn8yrdil0BEJFkMNTKnxG5+Z/E5ISJSJ4YaUhxB0X01wNcHL4hdgh2U/RoQkTQx1CgMrwOl/J6a6R/nQ1dVK3YZpBL62gZU1/LyDyQPDDUKs3jzUVTy7B/FqzDwNfY0NY5fq280ov9L36DfS9+gqZWLeRJJAUONAi3efFTsEkSlxp2P9CjvRVDj++rGq44b6htFrITIPgw1CnSmrEbsEkRlz3T25G7K+1av9vdVXYMyZuomZWOokTk1fnskIs8b/Mq3POxJksdQI3PWBsUq7zuyYxj0yB3U9L66WFOHRZuO4ERJtdnyrYd1IlVEZB9vsQsgcjU17HvEPG29+Y5OLdTwvrrm758ewLdHS/EvXnWcZIY9NQqk9FOaSTyCIGDkmzvtaOeBYshtDhbxCvEkTww1CqT2/YmaDhN4mprP6uW1n4ikj6GGFEcNZ6moYRulhs84kfQx1Mjc2YsGy4Uq7/vnF2oi5zA0k1wx1MjYqdIaPPzObrHLIFIH7ufV/n2JZIChRsZ2ny4XuwQi1bA30zQ0cZI6IrEw1CiQ2r9MqWFAp9KvRC5nGbnnxS7Babb+hFTwp0Uyx1BDilNTy2vUuAv3aa278XpJcsXDTCRXDDUKpPYPpDe/PSF2CYql5reWGnoAieSOoYZk4+OcX8QugeykxPBjb6Rh9CESD0ONAil1vMU/Pj8kdgmSIfVTbv/ycT6qLjeIXQa1ka1OKbX3ApP0MdQoED94SApe/+a42CW4FI8+EUkfQw2RDMmhN06nrxW7BJeyt3eM4YdIPAw1ROQWRjVfKErGvj54AReqrAdSBjaSOoYaBSqtroPAY1CKtvWQjq+xh6llhz7943yxSyBqM4YaBSqrrsMrm4+KXQa50Yrtp/Dt0VKxyyAreOo3kXgYahTqg11nxS6B3Cy/4JLHH5O9Q+rGl5+kjqGGiMgO7IEBcn+p4FgpkjSGGiIiOzDSAJ/lFyEjr1DsMohsYqghkqC3vj0pdgl0g0NFVThyQS92GZLwxb5isUsgsomhhkiC7Ll+Fcc3eM5Db+8SuwQisgNDDRG5hVqHoGTk8vAMkVgYakiRKgz1YpegemrtSTp38bLYJbiVWl9XkgeGGlKkBV8eFrsEReL+jIikjKGGFOlMWY3YJbid1A/vZB0rxWkVvA5EJB0MNaRIaphKQw6HAR5b/ZPYJRCRirQp1KxatQq9evWCv78/4uLisGfPnhbbZ2RkIDo6Gv7+/ujXrx+2bNlidrsgCJg/fz7Cw8MREBCAhIQEnDxpfkrriRMnMGbMGISEhCAoKAjDhw/Hd99915bySQU4QZg0lFXXiV0CEamIw6Fmw4YNSE5OxoIFC5Cfn48BAwYgMTERpaXWr0Oze/duTJgwAVOmTMHevXuRlJSEpKQkHDp0yNRm6dKlWLFiBdLS0pCTk4N27dohMTERtbXXrxT70EMPobGxEdu3b0deXh4GDBiAhx56CDqdrg2bTUpnlEM3BhERuZTDoeaNN97A1KlTMXnyZPTp0wdpaWkIDAzEmjVrrLZ/6623MGrUKLzwwgvo3bs3Fi1ahIEDB2LlypUArvbSLF++HCkpKRgzZgz69++PdevWobi4GBs3bgQAlJeX4+TJk5gzZw769++PqKgoLFmyBJcvXzYLRzeqq6uDXq83+yEiIiLlcijU1NfXIy8vDwkJCddXoNUiISEB2dnZVu+TnZ1t1h4AEhMTTe3Pnj0LnU5n1iY4OBhxcXGmNl26dMEdd9yBdevWwWAwoLGxEe+99x66deuG2NhYq4+bmpqK4OBg009ERIQjm6oIa1R8UUupD6IlIiLXcyjUlJeXo6mpCaGhoWbLQ0NDbR4G0ul0Lba/9m9LbTQaDb799lvs3bsXHTp0gL+/P9544w1kZmaiU6dOVh937ty5qKqqMv0UFqpvQqyFm46IXQK5kcATrEkEfN+RlHmLXYA9BEHAM888g27duuGHH35AQEAAPvjgA/z+97/Hzz//jPDwcIv7+Pn5wc/PT4RqyR0u1zeKXQJBHmdcEZF6OdRTExISAi8vL5SUlJgtLykpQVhYmNX7hIWFtdj+2r8ttdm+fTs2bdqE9PR03H333Rg4cCDeeecdBAQE4MMPP3RkE0imLl1uELsEydHwutFERGYcCjW+vr6IjY1FVlaWaZnRaERWVhbi4+Ot3ic+Pt6sPQBs27bN1D4yMhJhYWFmbfR6PXJyckxtLl++Ou24VmterlarhdFodGQTSKaamthFQERELXP48FNycjImTZqEQYMGYciQIVi+fDkMBgMmT54MAJg4cSJ69OiB1NRUAMDMmTMxYsQILFu2DKNHj0Z6ejpyc3OxevVqAFfHy8yaNQuLFy9GVFQUIiMj8eKLL6J79+5ISkoCcDUYderUCZMmTcL8+fMREBCA999/H2fPnsXo0aNd9FTIT5WKei+aHDzuwV4MIiL1cTjUjB8/HmVlZZg/fz50Oh1iYmKQmZlpGuhbUFBg1qMybNgwrF+/HikpKZg3bx6ioqKwceNG9O3b19Rm9uzZMBgMmDZtGiorKzF8+HBkZmbC398fwNXDXpmZmfjHP/6B3/72t2hoaMCdd96JL774AgMGDHD2OZClrYd1WLbthNhleIzgYKhRw2BGNWyjXAmCAA1PwSPyuDYNFJ4xYwZmzJhh9bYdO3ZYLBs3bhzGjRtnc30ajQYLFy7EwoULbbYZNGgQtm7d6nCtSvXqlqNil+BR3H2TnAgCpxUgEgOv/SRTXir7xHT0rBsefiIxlSr48hA8A46kjKFGprRa7rSJpGruZwfELoFIlRhqZEp9mcbBgcKqe348g+N47PPLxctil0CkSgw1MqVV2V7b0S5vdpGTmPj2IxIHQ41MqS7UiF0AERFJHkMNyYLDA4XVlflIYhydgkAqSvS1YpdA5BSGGpmS50dm23EsB5H7xb2a1XojIgljqCFFOqarFrsERTpSrBe7BFlgBCcSB0MNyYJMe/PbxO5DFx5+Ts6U1eCP7+z27IPKlJLfrwreNFIAhhqSBTXNoPzOjtNil2DVgfNVYpdARNQihhqShR9Olotdgses3H5K7BLISRwDRiQOhhoiueIZXpKl5MNPRFLGUCNT3J8RSRdDDZE4GGpk6sgFnoWiVFcamuxq19jk2T0nD6kQkdQx1BDJ1L92nUX26Ytil0EqU9doFLsEIpsYaohk7IX/7Re7BLJCrjMK22N/YSX+l3de7DKIrGKokaFP9hSIXQKR3a7U23c4TUmUG2muej6DYZqkiaFGhuZ+dlDsEkgiPNkhoGnj8PRH0tQ3YZ+CO2qIJI2hhojc6jAvrUBEHsJQo3BGI78ykmvw7CcikjqGGoVTY9c/ERGpE0ONwuUXVIpdApHqsFeLSBwMNURERKQIDDVERC7Gs5+IxMFQQyQhlwz1DrVX8iRvRESOYqghWfs2eYTN2/ILLnmwEtc4VVYjdgnkAhpecZZIFAw1JGu3dWtv87aH35HfmV9aB/eG7Kdxv92nyx2+T4m+zg2VEFFrGGqIJMRLy6/4UvOn93PadL+yagYbIk9jqCGSEC8et1CMXy4axC6BSHUYaogkRCvhv0iOSXZMI2fzJvI4CX+EEqlPWy8aSdLTxFBD5HEMNUQyxt4TIqLrGGpItvx9rr59p997q8iVqAOH+ziGgZPI8xhqSPJsTTCXl/I7AMDsxDsQF9nZkyVJBq8xRER0HUMNSd7WwyVWl7fz8wYAaDQadGnv68mSVIk9D0QkdQw1JHmFFZdbbaOUHe7J0mqxSyAiki2GGpI8ox2JRSmhJuXzQ2KXQC6y6UCx2CUQqQ5DDUmePWfGKmVsiT0BjuQh/edCsUsgUh2GGpI8NfXUGOqbHGqvlO2Wq6dH8Mw7IilhqCHJs3X2k1kbD9RB1NycB6LFLoGIbsBQQ5Jn1+EnlaYalW42EZFVDDUkeZxunoiI7MFQQ5J3/tIVO1rJP/jU1DWKXQIRkawx1JDkfZp/vtU2Sjj8tPr702KXQEQkaww1pAgKyDSoutIgdgktUkJwJNvyfqkQuwQipzHUyExtg2On/KqFPWdIKZFKN5vcYOy72WKXQOQ0hhqZiX4xU+wSJIn7dvK0V/7YV+wSiKiZNoWaVatWoVevXvD390dcXBz27NnTYvuMjAxER0fD398f/fr1w5YtW8xuFwQB8+fPR3h4OAICApCQkICTJ09arGfz5s2Ii4tDQEAAOnXqhKSkpLaUTzIzsk9oq22U0GOh0WjELqFFEi/Po+4I7YDH424WuwwiasbhULNhwwYkJydjwYIFyM/Px4ABA5CYmIjS0lKr7Xfv3o0JEyZgypQp2Lt3L5KSkpCUlIRDh65f42bp0qVYsWIF0tLSkJOTg3bt2iExMRG1tbWmNp9++imefPJJTJ48Gfv378ePP/6IP/3pT23YZJIbH6/W36YKyDQkI/6+Xqb/56UkiFgJEd3I4VDzxhtvYOrUqZg8eTL69OmDtLQ0BAYGYs2aNVbbv/XWWxg1ahReeOEF9O7dG4sWLcLAgQOxcuVKAFd7aZYvX46UlBSMGTMG/fv3x7p161BcXIyNGzcCABobGzFz5ky89tprePrpp3H77bejT58+ePTRR9u+5SQbRy/oW22j1jE17D0Rx41Pe5f2fqLVIabzly6LXQKRBYdCTX19PfLy8pCQcP2biVarRUJCArKzrQ8yy87ONmsPAImJiab2Z8+ehU6nM2sTHByMuLg4U5v8/HwUFRVBq9XirrvuQnh4OB544AGz3p7m6urqoNfrzX5Ifkr0tThTbhC7DIIyDvGR6wz/53dil0BkwaFQU15ejqamJoSGmo9xCA0NhU6ns3ofnU7XYvtr/7bU5syZMwCAl156CSkpKdi0aRM6deqEe++9FxUV1k9DTE1NRXBwsOknIiLCkU0liThdVmNXO3sOURERkbLJYk9gNBoBAP/4xz8wduxYxMbGYu3atdBoNMjIyLB6n7lz56Kqqsr0U1hY6MmSycPmPdjb6nKlH5ZS+OZJFg/7EUmTQ6EmJCQEXl5eKCkpMVteUlKCsLAwq/cJCwtrsf21f1tqEx4eDgDo06eP6XY/Pz/ccsstKCgosPq4fn5+CAoKMvsh5bqtW3ury3nZKCIi9XAo1Pj6+iI2NhZZWVmmZUajEVlZWYiPj7d6n/j4eLP2ALBt2zZT+8jISISFhZm10ev1yMnJMbWJjY2Fn58fjh8/bmrT0NCAc+fO4eabeVqlGn098x672im9p6a8pg4/nbkodhmqw44aImly+PBTcnIy3n//fXz44Yc4evQopk+fDoPBgMmTJwMAJk6ciLlz55raz5w5E5mZmVi2bBmOHTuGl156Cbm5uZgxYwaAq3NzzJo1C4sXL8aXX36JgwcPYuLEiejevbtpHpqgoCA8/fTTWLBgAb755hscP34c06dPBwCMGzfO2edA8YwK7K7oHW7Z87ZwzJ0Wy+S06W09pPHY6p9cWwgRkUx5O3qH8ePHo6ysDPPnz4dOp0NMTAwyMzNNA30LCgqg1V7PSsOGDcP69euRkpKCefPmISoqChs3bkTfvtdn45w9ezYMBgOmTZuGyspKDB8+HJmZmfD39ze1ee211+Dt7Y0nn3wSV65cQVxcHLZv345OnTo5s/2qkPDm98ic+Rv4estiCJUZjQPfie/sHmyxTJDRDDZK7lTa8HMBxg/uKXYZLiP1iRKJ1MrhUAMAM2bMMPW0NLdjxw6LZePGjWuxR0Wj0WDhwoVYuHChzTY+Pj54/fXX8frrrztcr9qdKTNgX2ElhkR2FrsUt7K2n1FyUJCTv396UFGhhq6qa2yCn7dX6w2JPER+X92pTZQ+tgSwHmBUsNkeUVR5Bc9l7Be7DJKYukaj2CUQmWGoIQWxTDBGGaWaf+8+J3YJNr2y+YjYJUgKDz4RSRNDjUrIZ9duzpGhC1Z7alxXiqrVNfAbOVmS0XcGUgmGGpWQ64ePI3XLdBNJhjhOmEiaGGpIMayduq6GsURERHQVQw1JmkOHn+xcRo7j82hucVI/s98//r84kSohohsx1JBiWOuUOaGr9nwhpGiLkvrijrAOZsvuvi1EpGqI6EYMNSohp0no2sraNj6Sli1CJcqj1iEkTVYOaar1uSCSA4YatVB+plHHNpLHCIKAxOU7LZeLUItk8ckgiWGoIcXg56v7qPG5bTIKOFVa49Q6nL0/ETmGoUZGqq40iF2CpPFEJ5KaN789IXYJ7sVjcSQxDDUy8vrW42KX4HGOfGaqYdwQkaTwT44khqFGRgoqLrf5vmr47LEyppPI9RTYJZhfcEnsEohcgqFGJRT4OWyBE+2RK6nl3XSypBoPv7Nb7DKIXIKhhhRDLTshko/NBy6IXUKrDhZViV0Ckcsw1JCkfX1IZ39jGaeaK/VNYpdAzail44/XsSIlYaghSfv37nN2tzXKeC+0ePMRp+6feci9PQI8tHcdn4nrHPn7JPIEhhoZqalrbPN91XBmkLXZX+XCoR4pK57+KN9FldA1avibAQCNE+dlK/6UdZIdhhoZaXRip62GL9py7qkhEgsPP5GSMNSQYvj5eIldgmJpVLjnY0Ymkh+GGpKd/0wZYnX5b6K6ergS9eCYmuv4VBBJF0ONnPDTFPdEheAeG+HFS6tB53a+Hq6ISN7U2AtHysVQIydOfPgoJQ619gHMj2ciIvViqJGRo8V6sUuQPH7pJFdhxyiR/DDUyEh9k1HsEkTXWmaZnRjtkTrU5tJlXiH+GlvjiyI6B3i4Etfg9wBSEoYaUpRR/cLELqFNpLxjWfXdKewrrBS7DI+zNU9NROdAq8sz/jzMneW4DXs3SUkYalRCKZ9brX0AK2U7peS1rcfFLkEy/vybW/Db6G5WbwsL9vdwNa7hzOR7RFLDUKMSHB4gbXx9pMfaUaapv7mFZwsRSRhDDclKa7sT7nCIiNSLoYZkpbUeDblGGrnWrWTW3mtKfJ34PYCUhKGGFEWuH9A8/ERicfZPhrNNk5Qw1KiETPf1FpSyHTcSBAEVhnqxy6BmrO2seXjT0rYjJWKXQGTCUKMSavkuJcczOcprGGjkQn7vrtY5+9lQWl3nkjqIXIGhhiTLaHT8mzK/SJOrqOWLAJGSMNSQZH2+t8gl6zmuq3bJetyFQUw++FpZYvgjKWGoIcnKK7hksawt+5SDRVXOF0Oqo5bxr8xppCQMNSqx7JvjaFDBtaOsfZPm2RnkKnIcs9Ua/nWQkjDUqMSB81X4cPc5sctwGrv/yWO4t7cPvzSQhDDUqMip0hqxS3BIWz4rrX2TlvpHLnOactzcxfrFLpXMUN8kdglEJgw1MnHFBR8c9Yo4/MSzn8gzrF2l29aVu9VsydfHxC6ByIShRibu/ud2p9ehhl5iq5lGBdstB0szufOTmkuGevzl43yxyyByGYYameCMs1e1pSfGqIY0JwPv7DiN0upascuwm7W3jdLeSv9k0CSFYaghRbE2Od+czw6KUIn9XDX1/tELepesx53qG5VwCFQ5SvTyCZlE9mCoIQmz/FocHuzf4j3UPKTmgbd+ELuEVsmpp6PJSrGtlS+n7SNSIoYakpXnRt4hdgnkBLkcCrxc34hBi791+H5NVi7tQUSew1CjInKfhC6kvS+CA3xabMOzn6RNLm/BzEO6Nt3PWqj5Yp9rLvdBRK1jqCHJar4D/Nvvbm/1Pq4an+JJrqxY6sFV2tVdZ+tt1Nrz66W1vOPM9H0uqMg9XPX3IvX3HalHm0LNqlWr0KtXL/j7+yMuLg579uxpsX1GRgaio6Ph7++Pfv36YcuWLWa3C4KA+fPnIzw8HAEBAUhISMDJkyetrquurg4xMTHQaDTYt29fW8pXLTnu8G/UKdBX7BLISXI5/GSLj3fLH5mrHh/ooUqk5ZjELxpL6uFwqNmwYQOSk5OxYMEC5OfnY8CAAUhMTERpaanV9rt378aECRMwZcoU7N27F0lJSUhKSsKhQ4dMbZYuXYoVK1YgLS0NOTk5aNeuHRITE1Fbazkyf/bs2ejevbujZRPk922qebn2RrK/3h/l8lrcSeZZ0yFGmYw5sTYz9cIxdyLIv+XDnzERHfGcHT2KSsOz2kgqHA41b7zxBqZOnYrJkyejT58+SEtLQ2BgINasWWO1/VtvvYVRo0bhhRdeQO/evbFo0SIMHDgQK1euBHB1R7t8+XKkpKRgzJgx6N+/P9atW4fi4mJs3LjRbF1ff/01vvnmG7z++uuObynJzoVmp5vau/NXUUawILPcKlnW3msT43t5vA654NuOpMKhUFNfX4+8vDwkJCRcX4FWi4SEBGRnZ1u9T3Z2tll7AEhMTDS1P3v2LHQ6nVmb4OBgxMXFma2zpKQEU6dOxX/+8x8EBrZ+fZW6ujro9XqzH5KXnSfKmi1Rc1xRBjX1ShGR5zkUasrLy9HU1ITQ0FCz5aGhodDprJ8toNPpWmx/7d+W2giCgKeeegpPP/00Bg0aZFetqampCA4ONv1ERETYdT+SLu4QSQ74PiUSjyzOfnr77bdRXV2NuXPn2n2fuXPnoqqqyvRTWFjoxgpJSuTWFW5t/EZbyW3bpUrug+rtpY6tJDVxKNSEhITAy8sLJSUlZstLSkoQFhZm9T5hYWEttr/2b0tttm/fjuzsbPj5+cHb2xu33XYbAGDQoEGYNGmS1cf18/NDUFCQ2Q/Jm1I/gFO+ONR6I5HUNjh/dfgbccwPEbmTQ6HG19cXsbGxyMrKMi0zGo3IyspCfHy81fvEx8ebtQeAbdu2mdpHRkYiLCzMrI1er0dOTo6pzYoVK7B//37s27cP+/btM50SvmHDBrzyyiuObAKRpBw4X4mv9heLXYZN49Ksj5UjIpIib0fvkJycjEmTJmHQoEEYMmQIli9fDoPBgMmTJwMAJk6ciB49eiA1NRUAMHPmTIwYMQLLli3D6NGjkZ6ejtzcXKxevRrA1W7eWbNmYfHixYiKikJkZCRefPFFdO/eHUlJSQCAnj17mtXQvn17AMCtt96Km266qc0bT/KideKQgCAIkjykkHvukkvXd/W0fddt58GiKpetC5DPeBNnypTi+4xILRwONePHj0dZWRnmz58PnU6HmJgYZGZmmgb6FhQUQKu93gE0bNgwrF+/HikpKZg3bx6ioqKwceNG9O3b19Rm9uzZMBgMmDZtGiorKzF8+HBkZmbC37/lixeSuti7r+jVxfLsOEGQzw6VxOfMe0VO80G56m9CTttMyuZwqAGAGTNmYMaMGVZv27Fjh8WycePGYdy4cTbXp9FosHDhQixcuNCux+/Vqxf/iNpALc/YmJgeSP7vfrNlatl2qZPLn60zg7flso1ESiSLs5+IAPu/VVq7/o5aQrA6ttL92KtHJE8MNSqi5s9p7uylQS5hwZky+V4jEg9DjQyUVlteA6st5P5hy0MCJAfW3mtq6SkkEhtDjQwUVlxxyXoammR+0TlnBm/KPtLZR+r7TqnXd41TA4WtvNfkst1EcsdQoyJbDuqw+cAFsctoM6cOCXCnQg5x7XEyvv2IPIOhRhZc95H4zPp8l62LyFFVVxrELsEuzp3SbW2ZVGONa8LbBz+cdcl6iJzFUCMDZdX1YpcgCc5MaibZfYqLSf0w2yNp2Tim04tdRqtcPVBY2q+K8zYflG8PMCkLQ40MzP3sgNglSIJMTpyhVvz35/Nil9Aqp2YFtpKg1RKqicTGUCMDly7Lo8teyqTeg0HS4vqeGr7/iDyBoYZkw7lrP7mwEAlTy3a6m+vH1LR9fURkP4Yakg1ndjQ85k+OkMskgURkjqGGZMOZ/cy/eHYGOcCZiR4bjPKYD6q+0YifzlwUuwwil2KoIUn678+FFsucOvtJomMaXN0jcP6SayZqdCdZ9II4UeNjg3taLJPi4aeXvjqMmrpGscsgcimGGpKk2Z9anvHlzM7QKMGdijskvPE9iiqlH2ykzpncFRnSDpv/OtxsmRRD9fqcArFLIHI5hhqSDWd2NJUqOoMs/5dLYpcge9au9O6IIH8fs9+l2FNDpEQMNSQbWgd2NCNu72r2e3lNnWxms3UW95/Oa36mXaCvl0P3b96ryNeEyDMYakg2orq1t7vtojF9LZYdKqpyZTmkYKXVtWa/757zW4fu33z8l3Qvk0CkLAw1JAtfz7wHHQN9nVqHHManqoEcXofnM8zHdDn73mOkIfIMhhqShfBgf7FLkA32CjivycmR5c2DW32jPE7zdoZRLaPxSdIYakgWHD2dW4pnm5A8VBicv4Bs87froMXfKn6n/9WBYrFLIGKoIXlw9GQUq50VEjvuIQgCDnKcj+Qs++a40+uwNnnfkQvSvzq5Mw6c53uZxMdQQ7Lg6HWfrGcaaaWajLzz+Cy/SOwyqJnjumqn12Ht7dqo8J4aHvUkKWCoIVlwdOI9a+NKpDaT7QYrsyaTcjk7ToeIWsdQQ5JjLZA4c4VuqXLXTo7fmJ3jireatVUYFf7CcBwbSQFDDUmOtX29wz01rinFrfYVVopdQou+3O+egZ8KzKeW1LCNzSh9IDTJA0ONxKnxg8LaN1pHx8NY+1Kswv2MU/76yV63rHf7sVJJn+Kswj85l+DzRlLAUCNx9U3S/fB3F2uHZRw9+8nbyh2cucq3nEj9MMDpMgPe/PaE2GXY5JJ5fqT9EriF1N93pA4MNRKnxsGF1vYpjo6publLoIuqkR85DN34X955sUuwyRV/cjJ4CVzuo5941W8SH0ONxKW74QyZ3afLXb5OV7J6+MnBThZrvTIq6aghJ3FGZiL5YqiRuB9PuT6A/On9HJev05WshxrnEwkzDdnDJT01zEVEomCokTglnsrcGhUecSMJccWp1xxfQiQOhhqJU2GmYfc/icol44T5FiYSBUONxG07UiJ2CR7nrp4aqQVEd9XDHapzXPH08SUgEgdDDUmO+874kliqIUly17uEYZPI/RhqSHLcd/hJWnsVRixpckUPmqPzKhGRazDUkOS4q6NGat+U1TIZoBqFBflbLOPLTeR+DDUkOe668J/EMg1JlEsuaKnR4J6oEOdXREQOYaghyXFXqFHjdbSkSsqdFo5eZ8wWqfUMEqkBQw1Jjrt2BmrZx8hhO+VQo7PCgy0PQRGRezHUkOQcLq5yy3qldmVodw2Ifj5jP365aHDLutWg+dl3PTu37Tpicx/s7Ypy3OI/2efELoHILRhqSHKe/ijfJetZ/WSs2e8T1+xBWXWdS9YtdYs3HxW7hBZJ+fBTcICP2e//mx7fpvV0bufrinLc4sUvDotdApFbMNSQ5EWHdWjT/UbeGWaxbOX2k86WIwu1DU1ilyBbB85Xmv4f0t4X3TrwMBKRXDDUkORt/us9LltXncQOQbkLB6m2TUOTEYb664HQlc+jGl6TU6U1YpdAKsdQQ5LnpdCZzHgylvRsOXjBbeue8+kB1DUquwftUJF7xsMR2YuhhlRFKt+Wvz9RJnYJZEV1baPZ7658u5wpN2Dd7l9cuEYiao6hhiSl8nK9W9fvrjlwHDVpzR63rl/qs9eWVtehwuDe19oVXH2GWlHlFZeuT2qqaxvELoFUjqGGJOXzvUVuXb9W6nt7FXnigxyxS/A4qYRqd+FZVSS2NoWaVatWoVevXvD390dcXBz27Gn5W2dGRgaio6Ph7++Pfv36YcuWLWa3C4KA+fPnIzw8HAEBAUhISMDJk9fPUjl37hymTJmCyMhIBAQE4NZbb8WCBQtQXy/9b3rOaGxSx6DWG3l7uTdnM9NIx5ELerFLsPCfbPPDQ8qOIETK4/AeZMOGDUhOTsaCBQuQn5+PAQMGIDExEaWlpVbb7969GxMmTMCUKVOwd+9eJCUlISkpCYcOHTK1Wbp0KVasWIG0tDTk5OSgXbt2SExMRG1tLQDg2LFjMBqNeO+993D48GG8+eabSEtLw7x589q42fLwxb5isUvwOHdnDoYaasnxkmqz313dsaLwjhoi0Tkcat544w1MnToVkydPRp8+fZCWlobAwECsWbPGavu33noLo0aNwgsvvIDevXtj0aJFGDhwIFauXAngai/N8uXLkZKSgjFjxqB///5Yt24diouLsXHjRgDAqFGjsHbtWowcORK33HIL/vCHP+D555/HZ5991vYtl4HnMvaLXYICMdW0RhAETHTzmB8iIndwKNTU19cjLy8PCQkJ11eg1SIhIQHZ2dlW75OdnW3WHgASExNN7c+ePQudTmfWJjg4GHFxcTbXCQBVVVXo3Lmzzdvr6uqg1+vNfui6s+XSnEbf3T0pCj073KVOl9VgJ8/OAuC+S1kQkXs4FGrKy8vR1NSE0NBQs+WhoaHQ6XRW76PT6Vpsf+1fR9Z56tQpvP322/jzn/9ss9bU1FQEBwebfiIiIlreOJW57/UdYpfgEc1nI+bhp9Zx/pzr+FQQyYvszn4qKirCqFGjMG7cOEydOtVmu7lz56Kqqsr0U1hY6MEqqa24Q3WNH06W40xZ22Z3ZeeE+wiMSURu5VCoCQkJgZeXF0pKSsyWl5SUICzM8jo7ABAWFtZi+2v/2rPO4uJi3HfffRg2bBhWr17dYq1+fn4ICgoy+5ETtV67x+jiVNN8B62mU7offc/24VuyEzMIkaw4FGp8fX0RGxuLrKws0zKj0YisrCzEx1u/km18fLxZewDYtm2bqX1kZCTCwsLM2uj1euTk5Jits6ioCPfeey9iY2Oxdu1aaLWy62RyyJkyaY55cTdXz+PR/JuxeiINUF6j7CkPPMHVmUajgnegWr+QkTR4O3qH5ORkTJo0CYMGDcKQIUOwfPlyGAwGTJ48GQAwceJE9OjRA6mpqQCAmTNnYsSIEVi2bBlGjx6N9PR05ObmmnpaNBoNZs2ahcWLFyMqKgqRkZF48cUX0b17dyQlJQG4HmhuvvlmvP766ygruz6I0VYPkdypqEPBjKsPPzVfn0YCTyw/9NVLDYef9Fca4O/jJXYZpFIOh5rx48ejrKwM8+fPh06nQ0xMDDIzM00DfQsKCsx6UYYNG4b169cjJSUF8+bNQ1RUFDZu3Ii+ffua2syePRsGgwHTpk1DZWUlhg8fjszMTPj7+wO42rNz6tQpnDp1CjfddJNZPTw7QVlcf/hJeu+PpFU/il1Ci9Sw47WXs++fpWP7Y/anB1xUjTw0cmAcicjhUAMAM2bMwIwZM6zetmPHDotl48aNw7hx42yuT6PRYOHChVi4cKHV25966ik89dRTbSlVtiTQoSCKoIA2vSVtav7xKoXn9ZiuuvVGJAnO7p69vSTwhvMwKfyNkXope2CKjKnh2Ls1gb4uDjXN9krenKimVWp971kjwY4+p3EOIlIyhhqJUuu3HZcPFG62vg93/4KqK7yScEt4+Om6NU8Ndur+zQ/FSCEkuXu26EsG/n2ReBhqJEqlmcYi1PS/Kdip9Y2803wgeX2TEfM+O+jUOkkdnhrWC/G3dnFqHa4eIyYHS7ceE7sEUjGGGpIUY7MLk//v6WFOrS/5d7dbLPue3e9kBz8f5z8eb+na3ux3NfTAlurrxC6BVIyhhiRFX2vede3r7dxb1NqppVI8I4rE19YZmFsyJNL8+nQf/VSAEn2tyx9HSvjXRWJiqJEoNXyja65EX4uXvzri9sdR4REByfrXrrNil2Ay/4vDHnmcv23Y55HHEUt1LcfUkHgYaiRLfanmxY2HPPI4ajoBavepcofv48mOrEWb3B9i7VXfZGy9kQscvaD3yOOI5fylK2KXQCrGUEOScaHKM93yarr+058+yBG7BPlijx6R7DDUSFTVFfVdt8fVp3PbJGKmkdLhFjK352yFRx6nvtEzPUJEasRQI1Fj31XfFZZLqz1z1oSYPTVSOtxiS01do9glKJqnDnMRqRFDDUlGmcdCjUceRpYMdY0Yl6a+QG2Nu/oNG5qUf1yrkcGNRMJQo2IT1+xR5enNUrhSt1S9/8MZsUsgBcjIOy92CaRSDDUS9GMbzlhpi50nyqCvVf6hhl5dAs1+Z6SxrUYF7wdyvwuVPAOKxMFQI0G7PBRqAHV0Ez82pKfZ7+ypsY1z+Fynxl5MIrljqFG5JhXsxZpvI8fUkCe9+/hAsUsw8dS4NVXOHkqSwFBDitf8G3dpdR0qDOo7Zd4eTc0vvqViruqouTsqxDUrcpKhrhGDX/lW7DKI3IqhhhTPWmfU45yUzkLVlQZ8mP2L2GVIhpeXa3obpNJnUXjpsuceq8Jzj0V0I4YaCfLkoXwpHXz6/YDubllvg5VxQ2JMVb/rpOfGSt3ow93n7Gr3xb4i9xYiMzd1Cmy9kR2kMoZL48F49flevpdIHAw1JBk+Lvpm3NyD/cLdsl5HFFVewRP/Eqd3aMGXnrlQo5L06xGMxwZHuGRd1sZwXRLh8KdEshWRWzHUSJDgwf4TqZzgUXW5AZ/lu+fbXe/wILes1xFFMrjIn1j7vL98nCfSI1+XeUhn9vuipL7w8XLNx6O1HhIxDn8y05AaMNRI0NEL1R57LE8GqJY8/ZHljm3oLZ1FqMQ9ZHF6sEhf5bcc1Il+Fl7z958rz5Czdk2zIwq/UjeRWBhqJGjniTKPPZZU9rXZZy5aLEt7IlaEStSrUsQzwjx2MVM7ufL6YO38vF22Lmfw8BOpAUONyklrV3Jd7/AgdAz0FbsMl5Hq83yjZdtOiPbYEss0ol701H08u028MCqJgaFG5RoapTkvCSfI86wfTnqud9AaqRwGvUaRmcbDFn0l/SvSk/Iw1Kjcva/vQKm+VuwyLHCn4lkrsk6K+vhS66kh523ILRS7BFIhhhqJOe/BCbKuWb+nwOOP2RpPdP976sKh5TV1eGz1Tx55LLkaveIH1DY0iV0GEckcQ43EzFi/V+wSJMETHTWeOq32ne9Oe+RxWrLwqyOSPgPrdJkBXx+6IHYZJhJ+qtpEEARsOSid55fIXRhqJGZfYaXYJUiDi3tqfhvdzaXrc0R9k/g9EGt+PIsTJTVil9EiJV8w/ud/JIj6+FlHS/GGiAPBiTyFoYY8On26vVxd0eKkvi5eo/2k8q2fh3ess3aWjtbFn4xdO/hZLPPk3Dx7Cy957LGIxMRQQ6KrutxgsczVQ2qCA3xcu0IHNDZJI9VUXLY+D82Jkmr8fE78nV6jSF01o5bvFOVxx7+XLcrjEikZQw2J7g+rdlksc3VPjVgToH17pEQyZ4FMXvszDhVVWSx/PmO/CNVYmvPZQRwutqzP3c5buYRFO1/3v19yfxE/SBIpDUMNiX769C8XLc/4cseVjR8e2MNiWb2b5+mZLoHrGt3o45xfzH6vutKAA+c9HyRsWZp5XOwSMP3eWxHR2TVX6Fa7BiUPlCJJYqghrPruFKprLQ8BialzOzfMJmzlKNC9r33n+sf51eX6RjRI5NCTLWIdepGyF0beIXYJimGtZ5DInRhqCHWNRryy+ajYZZhZOOZOl6/TWrworqp126nOk9bscct6nfHJnkKU3DDZ4oUqaU286OkIaG2wrpbTWbvMH9/ZLXYJpDIMNRLyxb4i0R7753MVoj12c/94sDfCgwNcvl5b4cVdZ6FIYfCtNX/95OpcSAcldNjpmp0nyvC9By/oyrlbiJSFoUZCZqbvE+2xpXSQxF1jfGxll0Y3hJolXx9z+Tpd5XCxHgDw+5WWA7SlwJM9XJ48rbpPeJDFMk+MOblkqMcqCUwASeQJDDUkKk8ec7e1Axuz8keX7lyMRgFp30t3J1JT14h12efELkMSZv/vgMce670nYy2W3fvaDrfP9JzyxSG3rp9IShhq6CqRumoeetvK6dxu6qoZcXtXq8uPl1Rj10nXXAdKEAQkvfOjS9blTvO/OCx2CS366cxFjzxOvQfPzrF2RlVR5RWbPYiuIsXDjETuwlBDAIAz5QaPT35m6xuqu4Zpjo29yeZt+89XOv2NWRAEpH1/RlKnSMvVY6t/wiXD1ckCy2vq8K9dZ/HujtP40/s/4bZ5W/D//v0z8pyc56VEIlend/chKKPIU1pfrrecsZnIXcSZkYwkae5nB/HauAEee7w6G3PEuGtMjZdWg/kP9cHCTUcsblv+7Un07R6MhD6hbVp3k1HAq1uO4l+7zjpbJv3qrkXbbN62/Vgpth8rRerD/fDY4Ig29e5Zex/88S7LuYzczd3jesS+TMdb357E3Ad7i1sEqQZ7aiRizqeeO7ZvS0beeY8+3pV669ciGn5biNse866eHW3e9n/rclFa7fi396rLDXjo7V0MNCKY+9lB/Oa171BYYTmBY2vOlhksli3zYKg31VFuWYerCIKAi4Y6t63fHu/tPCPq45O6MNRIRPrP0phK31POlhusfhNPnzYUUaEd3Pa4rXXFD3klC6V2HpYQBAE/nCzDgIXf4OgFvSvKozYorLiCe5Z+h21HSuy+T/qeAhyx8pqJMUfNQ2/vQlGl5aUaXOG1rcdR2yD+rL7u2j6i5hhqJMDdZz9I0fSPrF8+4HY3BhrA9mndNxryalarZ2Ud11Xj9yt34cl/SW+CPbWaui4XJ0qqW20nCALmfHbQAxVZ+vuoaKvL8910Hah3dkjjLLyaWo6rIc9gqJGA01a6wcXyvge6ivN+uYRjOus7Hy83f1P2tnP9D729C+PSduPjnF9Q13j9MJkgCFj13SkkLt+JQ0XsnZGakW/ubPVLwqwN+zxTjBXT773V6nKlf63Zc9YzZ7MRMdRIQMpGcb41WvPKFvdeLiHvlwqMfdf21Ok+Xu4NNQNu6mh325/PXcI/Pj+EO1Iyseyb4yisuIyRb+7Ea1vFv+gi2Tb8n9/ZPONm04FifLGv2MMVte6vn+xt07iglnhyZubWvPjFYcldX46UiaFGZI1NRvx0RjqXKABsD+B11pf7izH23ewW2/h6ufctqdVqcDb1QYfv9/b2U7hn6Xc4WVrjhqpsz6EjJbd1a+/wfaxNOOduRZVX0Gf+VrNxHKfLajD9ozzMWL/X4/XYa0mma2ehltq1x9b+eE7sEkgFGGpElnlY12qbRUl9sWHaUCTFdHfP1aubGbn8e5eur7S6FqlfHzVdc6gl3m4ONYD7JvdzxtqnBtvVrnuwv5srse22ro6HmuiwDghp7/73rDV3L9mO3i9moteczbh/2ff4+lDLf2vp04Z6pK4OftZn0th84ALOlLkmNHtqDp4nhvbEP8f2Q1iQP3p1sZxc8EZvbDvh8bmwSH0YakR0sqS6xW+OZ1MfxLklo/Hk0JsRd0sXLH/sLvz499+2ut6eVmYudURhxRWnx9YIgoCjF/QY/142hryShfe+b319m54d7tRjypm1s24+/8swi2WzR0Xjvjs806szrtlkhYl92zaHj5gh8kqD/b2OQ2/p4sZKrls3ZYjN23735k6nJ+NrMgqIezXLofvcE+X4NAqPx/XE4qR+GD+4J36adz92vHAffph9X4v3+X8f5jr8OESOaFOoWbVqFXr16gV/f3/ExcVhz56WuzkzMjIQHR0Nf39/9OvXD1u2bDG7XRAEzJ8/H+Hh4QgICEBCQgJOnjxp1qaiogKPP/44goKC0LFjR0yZMgU1Ne45FOAJRZVX8Ls3d9q8/cTiB6zuDAJ8vfDD7Ptwf3Q3m/ddPTEWva1cPO+mTvZf+fqVLUcxavlOGB2YGEwQBBRWXMZ7359G5NwteOCtH5Bz1v5Da317BNvd1lmPDrI9u7BU3NWzk9XlYVZ6a76eeQ/CXdSL86e4ntg/fySWjO2Pj6bEYWSfUKQ9EYukmNYnpru52bf1djZ6JdTsrp6d8IcB3a3e1mQUMP697DYdAhYEAafLanDrvC2tN77BVzOG4z9T4hx+vDvCLM9UjOgciHNLRtvsjdp5ogxf7ZfemCZSDoc/cTZs2IDk5GSkpaUhLi4Oy5cvR2JiIo4fP45u3Sx3tLt378aECROQmpqKhx56COvXr0dSUhLy8/PRt29fAMDSpUuxYsUKfPjhh4iMjMSLL76IxMREHDlyBP7+Vz+oH3/8cVy4cAHbtm1DQ0MDJk+ejGnTpmH9+vVOPgWeVVPXiA93n2txsOmZVx9scb6MiM6B+NdTg9Frzmart0eHBeHrmffgyX/l4Idfr2m0/v/iMOy2EJv3seaYrhq33PABGRnSDjV1jWgyCtBfaXDp1a373+S5QAMASx8ZgP/menayQVeIu6UzfnN7V3yyx3xeo6AAH3w6fRgWfHkYRZeuIPl3t6NrBz+MWeXYdahG9w/Hq3/sZ/p9eFQIhjvwLf7Z30bh+Yz9AIDXHumPkPZ+mPdgNP62Yb+pTQc/b1TXSesU339NGuTRx7u1hUN5+QWV6D0/EwDwSOxN6BMehO4dA9Ap0AeRXdvBW6tFcIAPKi/X49xFA/YVVmHNrrMOzwVz3x1dMXFYL/T79W/v8bie+DinwK773n1bF0wY0tPm7QdfTsTYd3dbvZTFs5/sxbOf7MVXM4ajb48gSR4OJvnSCA5OkhIXF4fBgwdj5cqVAACj0YiIiAg8++yzmDNnjkX78ePHw2AwYNOmTaZlQ4cORUxMDNLS0iAIArp3747nnnsOzz//PACgqqoKoaGh+Pe//43HHnsMR48eRZ8+ffDzzz9j0KCrHz6ZmZl48MEHcf78eXTvbvmtp66uDnV112fS1Ov1iIiIQFVVFYKCLHsx2upUaTU++qkARkFAo1FAU9Ov/xqNv/579XejUcAPp8pRb+PSANccejkR7e38dmsroJxbMhoAUFhxGX/5OB//d08kxvz6LXvLwQv4y8f5DmyhZ5x85QH4eGA8zY0+2VOAuSLMV/LE0J745nAJSquvvj9H9gnF6omDzF7P50fejhm/jTJbtnvOb9G949Xetq/2F+PZX8copYzujf+75xarjzXyze9xosT+Hs1r7x1bNh0obvGQ6VuPxWBm+j6LdZ0pq8HUdbno0SkQayYNQml1HYYt2W53Xe709oS78HsbPSfuUlPXiAfe2onCCs9PSndkYSICfS0/YxqajJj72UH8r9nM4qkP97P4O8lNSUBIe79WHyt9T0GrcwL17RGEHh0DEODjBR8vLby9NNBqNNBocPVfSHMcHFl3a7f2eHLozS5dp16vR3BwsF37b4dCTX19PQIDA/G///0PSUlJpuWTJk1CZWUlvvjiC4v79OzZE8nJyZg1a5Zp2YIFC7Bx40bs378fZ86cwa233oq9e/ciJibG1GbEiBGIiYnBW2+9hTVr1uC5557DpUvXU39jYyP8/f2RkZGBP/7xjxaP+9JLL+Hll1+2WO7qUPP9iTKXnGWQOeseRIc5Vtep0mqcKq1BOz9vpO8pxOaDF3BHaAds/dtvWrzf5fpGfH1QhysNTUjZeMiZsl3iWi+SGI7p9Bi1/AePPd7/uzsS83/fB299exJvfnsCMREd8flfhkGj0eCSoR5Ltx7DI7ERiL356qGna6EmMqQdvnv+XtN6ahuaMGr5Tgzs2QlvjI+x+XiCICByrn2HI9KeiMWovmF2tTUaBeT+cgl/27APE4ZE4JM9hXg+8XaM7tcdT3yQg8GRnfBCovWJ5q5pbDJiyKtZqPj1wpVieHvCXXiof7goO83ahiYMWvwtajzUa/XE0J74+6hodPD3sdnm2ut6saYO0z/OR+d2vvj5Hwm4a+E3AIDNf70HdY1NuK2b/ZNkGuoa8fuVu3BGQvNxkfv85vauWPf/bI8bawu3hZri4mL06NEDu3fvRnx8vGn57Nmz8f333yMnJ8fiPr6+vvjwww8xYcIE07J33nkHL7/8MkpKSrB7927cfffdKC4uRnh4uKnNo48+Co1Ggw0bNuDVV1/Fhx9+iOPHzQ/ZdOvWDS+//DKmT59u8bie6qk5V25ARl4hvLRaeGs18NJqLP710mqhr23AsQt6lNXUIapbB3Tt4IdRfcPQo2MA/H28nK5DEASU19Sja4fWvz01v9+RC3p88MNZtPfzxoWqWrTz80K/HsGI6BwIDYBGowCtRgOjIOBKfRPqGo2ob7z6b4WhHlqtBjlnLkLA1cntQoP8UVZdh4YmI4ICfGCoa0RwgA8G3twJJVW1uLVbe/h5a3Fr1/YY2LOTKFPTN9fYZETuL5dwsqQaBRWXcd8d3QANEOjrjS7tfBEU4IP8gkvo1yMYQf4+uFLfhOq6BtzUKRCGuka08/PGcV01tBrAx0uLD3adQa8u7RAeHIDR/a++r0v0tejWwQ8ajQaNTUbsOVeBmIiOVr81X3O5vhHZpy8i7pYuFj14giDYtTOurm3Al/uLEdmlHXp2CcRNnQJRYajH29tP4qZOgejWwQ+/ub0rggNs7+zcqepyA+qbjGjn5wV/by9cbmjCcZ0eGo0GhrpGnCipwSVDPToG+iDI3wcHiipRU9sIHy8tqmsb0bNLIMKD/ZF77hJ0+lpEh3XAlYYm7C2ohI+XBl3a+WF0/3B0DPTBmTIDdp4owy1d2+GpYZHo0911nwVttb+wErtOlWPXyXJcNNShvtEIHy8tNBogNMgfPTsHwsdLC61Gg87tfODn7YUu7X0RHOCDjoE+OFlSg0uXGyBAQJd2vvD38YK/jxdqahuh1QJ3hAbh1m7t4Oft2OfMoaIqRHQORHCAj6mH2dfbud7US4Z6/HCqHKdKqnHRUA99bSP8vLW4qVMAAn290GgU0NgkwCgIEISrkxIKv/7fHgIEXO3bIbH0CmmHR2JdO2bRkVCj2FF8fn5+8PNzbAffFr1C2rX6jdQTNBqNw4Hm2v3u7B6MN1v4tq8G3l5aDL2lS4tnwNx3x/UxY77eWgQHXg0B1wbD3jhwcnFSPzQXGnR9IK+3lxbDbm29dyrQ1xv397Z+1pG9vQsd/H3weJx5d3Dndr5Y8Ps77bq/u117Hq9p7+eN2Js7m36/J8r8bK9HB0dYXc/kuyPtery/3h/lYIXuNSCiIwZEdMQz993Wpvvf+Fy50o0D950NM9d0audrc5A0kSs49E4NCQmBl5cXSkrMLxxXUlKCsDDr3dZhYWEttr/2b2ttSktLzW5vbGxERUWFzcclIiIidXEo1Pj6+iI2NhZZWdfnQDAajcjKyjI7HHWj+Ph4s/YAsG3bNlP7yMhIhIWFmbXR6/XIyckxtYmPj0dlZSXy8q5fBHH79u0wGo2Ii3P8VEQiIiJSHocPPyUnJ2PSpEkYNGgQhgwZguXLl8NgMGDy5MkAgIkTJ6JHjx5ITU0FAMycORMjRozAsmXLMHr0aKSnpyM3NxerV68GcLULfdasWVi8eDGioqJMp3R3797dNBi5d+/eGDVqFKZOnYq0tDQ0NDRgxowZeOyxx6ye+URERETq43CoGT9+PMrKyjB//nzodDrExMQgMzMToaFXj/sXFBRAq73eATRs2DCsX78eKSkpmDdvHqKiorBx40bTHDXA1YHGBoMB06ZNQ2VlJYYPH47MzEzTHDUA8PHHH2PGjBm4//77odVqMXbsWKxYscKZbSciIiIFcXieGrlyZPQ0ERERSYMj+29e+4mIiIgUgaGGiIiIFIGhhoiIiBSBoYaIiIgUgaGGiIiIFIGhhoiIiBSBoYaIiIgUgaGGiIiIFEGxV+lu7tocg3q9XuRKiIiIyF7X9tv2zBWsmlBTXV0NAIiIiBC5EiIiInJUdXU1goODW2yjmsskGI1GFBcXo0OHDtBoNC221ev1iIiIQGFhIS+pIGF8neSDr5U88HWSDzW9VoIgoLq6Gt27dze7tqQ1qump0Wq1uOmmmxy6T1BQkOLfLErA10k++FrJA18n+VDLa9VaD801HChMREREisBQQ0RERIrAUGOFn58fFixYAD8/P7FLoRbwdZIPvlbywNdJPvhaWaeagcJERESkbOypISIiIkVgqCEiIiJFYKghIiIiRWCoISIiIkVgqCEiIiJFUHWoOXfuHKZMmYLIyEgEBATg1ltvxYIFC1BfX2/W7sCBA7jnnnvg7++PiIgILF261GJdGRkZiI6Ohr+/P/r164ctW7Z4ajNUa9WqVejVqxf8/f0RFxeHPXv2iF2SqqSmpmLw4MHo0KEDunXrhqSkJBw/ftysTW1tLZ555hl06dIF7du3x9ixY1FSUmLWpqCgAKNHj0ZgYCC6deuGF154AY2NjZ7cFFVZsmQJNBoNZs2aZVrG10kaioqK8MQTT6BLly4ICAhAv379kJuba7pdEATMnz8f4eHhCAgIQEJCAk6ePGm2joqKCjz++OMICgpCx44dMWXKFNTU1Hh6U8QjqNjXX38tPPXUU8LWrVuF06dPC1988YXQrVs34bnnnjO1qaqqEkJDQ4XHH39cOHTokPDJJ58IAQEBwnvvvWdq8+OPPwpeXl7C0qVLhSNHjggpKSmCj4+PcPDgQTE2SxXS09MFX19fYc2aNcLhw4eFqVOnCh07dhRKSkrELk01EhMThbVr1wqHDh0S9u3bJzz44INCz549hZqaGlObp59+WoiIiBCysrKE3NxcYejQocKwYcNMtzc2Ngp9+/YVEhIShL179wpbtmwRQkJChLlz54qxSYq3Z88eoVevXkL//v2FmTNnmpbzdRJfRUWFcPPNNwtPPfWUkJOTI5w5c0bYunWrcOrUKVObJUuWCMHBwcLGjRuF/fv3C3/4wx+EyMhI4cqVK6Y2o0aNEgYMGCD89NNPwg8//CDcdtttwoQJE8TYJFGoOtRYs3TpUiEyMtL0+zvvvCN06tRJqKurMy37+9//Ltxxxx2m3x999FFh9OjRZuuJi4sT/vznP7u/YJUaMmSI8Mwzz5h+b2pqErp37y6kpqaKWJW6lZaWCgCE77//XhAEQaisrBR8fHyEjIwMU5ujR48KAITs7GxBEARhy5YtglarFXQ6nanNu+++KwQFBZn9zZHzqqurhaioKGHbtm3CiBEjTKGGr5M0/P3vfxeGDx9u83aj0SiEhYUJr732mmlZZWWl4OfnJ3zyySeCIAjCkSNHBADCzz//bGrz9ddfCxqNRigqKnJf8RKi6sNP1lRVVaFz586m37Ozs/Gb3/wGvr6+pmWJiYk4fvw4Ll26ZGqTkJBgtp7ExERkZ2d7pmiVqa+vR15entlzrtVqkZCQwOdcRFVVVQBg+vvJy8tDQ0OD2esUHR2Nnj17ml6n7Oxs9OvXD6GhoaY2iYmJ0Ov1OHz4sAerV75nnnkGo0ePtvis4uskDV9++SUGDRqEcePGoVu3brjrrrvw/vvvm24/e/YsdDqd2esUHByMuLg4s9epY8eOGDRokKlNQkICtFotcnJyPLcxImKoucGpU6fw9ttv489//rNpmU6nM/tDBmD6XafTtdjm2u3kWuXl5WhqauJzLiFGoxGzZs3C3Xffjb59+wK4+nfh6+uLjh07mrW98XWy5++LnJeeno78/HykpqZa3MbXSRrOnDmDd999F1FRUdi6dSumT5+Ov/71r/jwww8BXH+eW/rc0+l06Natm9nt3t7e6Ny5s2peJ0WGmjlz5kCj0bT4c+zYMbP7FBUVYdSoURg3bhymTp0qUuVE8vTMM8/g0KFDSE9PF7sUaqawsBAzZ87Exx9/DH9/f7HLIRuMRiMGDhyIV199FXfddRemTZuGqVOnIi0tTezSZMVb7ALc4bnnnsNTTz3VYptbbrnF9P/i4mLcd999GDZsGFavXm3WLiwszOIsgGu/h4WFtdjm2u3kWiEhIfDy8uJzLhEzZszApk2bsHPnTtx0002m5WFhYaivr0dlZaVZL8CNr1NYWJjFWWvN/77IOXl5eSgtLcXAgQNNy5qamrBz506sXLkSW7du5eskAeHh4ejTp4/Zst69e+PTTz8FcP15LikpQXh4uKlNSUkJYmJiTG1KS0vN1tHY2IiKigrVvE6K7Knp2rUroqOjW/y5NkamqKgI9957L2JjY7F27VpoteZPSXx8PHbu3ImGhgbTsm3btuGOO+5Ap06dTG2ysrLM7rdt2zbEx8e7eUvVydfXF7GxsWbPudFoRFZWFp9zDxIEATNmzMDnn3+O7du3IzIy0uz22NhY+Pj4mL1Ox48fR0FBgel1io+Px8GDB80+iLdt24agoCCLD3hqm/vvvx8HDx7Evn37TD+DBg3C448/bvo/Xyfx3X333RZTIpw4cQI333wzACAyMhJhYWFmr5Ner0dOTo7Z61RZWYm8vDxTm+3bt8NoNCIuLs4DWyEBYo9UFtP58+eF2267Tbj//vuF8+fPCxcuXDD9XFNZWSmEhoYKTz75pHDo0CEhPT1dCAwMtDil29vbW3j99deFo0ePCgsWLOAp3W6Wnp4u+Pn5Cf/+97+FI0eOCNOmTRM6duxodnYGudf06dOF4OBgYceOHWZ/O5cvXza1efrpp4WePXsK27dvF3Jzc4X4+HghPj7edPu1U4VHjhwp7Nu3T8jMzBS6du3KU4Xd7MaznwSBr5MU7NmzR/D29hZeeeUV4eTJk8LHH38sBAYGCh999JGpzZIlS4SOHTsKX3zxhXDgwAFhzJgxVk/pvuuuu4ScnBxh165dQlRUFE/pVou1a9cKAKz+3Gj//v3C8OHDBT8/P6FHjx7CkiVLLNb13//+V7j99tsFX19f4c477xQ2b97sqc1Qrbffflvo2bOn4OvrKwwZMkT46aefxC5JVWz97axdu9bU5sqVK8Jf/vIXoVOnTkJgYKDwxz/+0exLgyAIwrlz54QHHnhACAgIEEJCQoTnnntOaGho8PDWqEvzUMPXSRq++uoroW/fvoKfn58QHR0trF692ux2o9EovPjii0JoaKjg5+cn3H///cLx48fN2ly8eFGYMGGC0L59eyEoKEiYPHmyUF1d7cnNEJVGEARBnD4iIiIiItdR5JgaIiIiUh+GGiIiIlIEhhoiIiJSBIYaIiIiUgSGGiIiIlIEhhoiIiJSBIYaIiIiUgSGGiIiIlIEhhoiIiJSBIYaIiIiUgSGGiIiIlKE/w8KnzyapVYsIgAAAABJRU5ErkJggg==", "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/0k/8rf6137x12n9xmg_yll_84gw0000gn/T/ipykernel_6071/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": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGdCAYAAAD60sxaAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAXIxJREFUeJzt3Xl8VNX5P/DP7FknIUASAmETZZFNUCFWxYUSbGq10tbtp6ioXy3aCnWjX4qK7RdqrWgVl1YktpVStO4giyBYJKAiUWSTTQKECRBIJuus5/fHnXvn3pk7k4UkM5N83q/XvCZz75nJuRkyeXjOc84xCCEEiIiIiLoQY6w7QERERNTRGAARERFRl8MAiIiIiLocBkBERETU5TAAIiIioi6HARARERF1OQyAiIiIqMthAERERERdjjnWHWgvfr8f5eXlSE9Ph8FgiHV3iIiIqBmEEKipqUFeXh6MxvbL03TaAKi8vBz5+fmx7gYRERG1wuHDh9GnT592e/1OGwClp6cDkH6Adrs9xr0hIiKi5nA6ncjPz1f+jreXThsAycNedrudARAREVGCae/yFRZBExERUZfDAIiIiIi6HAZARERE1OUwACIiIqIuhwEQERERdTkMgIiIiKjLYQBEREREXQ4DICIiIupyGAARERFRl8MAiIiIiLocBkBERETU5TAAIiIioi6HARARJZ5tbwBlm2PdCyJKYC0KgB5//HEYDAbNbciQIcr5xsZGTJ8+Hd27d0daWhqmTJmCiooKzWuUlZWhqKgIKSkpyM7OxkMPPQSv16tps379eowZMwY2mw2DBg1CcXFx66+QiDqX/euA934JvFYI+Dyx7g0RJagWZ4DOPfdcHDt2TLlt3LhROTdjxgx88MEHePPNN7FhwwaUl5fjuuuuU877fD4UFRXB7XZj06ZNeP3111FcXIw5c+YobQ4ePIiioiJcfvnlKC0txQMPPIA777wTq1atOsNLJaJOYffy4Neehtj1g4gSmkEIIZrb+PHHH8e7776L0tLSsHPV1dXo2bMnlixZgp/97GcAgN27d2Po0KEoKSnB+PHj8dFHH+HHP/4xysvLkZOTAwB4+eWX8cgjj+DEiROwWq145JFHsHz5cnz77bfKa99www2oqqrCypUrm31hTqcTGRkZqK6uht1ub/bziCiOCQEsGA44jwA3LgUGXxXrHhFRG+uov98tzgDt3bsXeXl5GDhwIG6++WaUlZUBALZu3QqPx4OJEycqbYcMGYK+ffuipKQEAFBSUoIRI0YowQ8AFBYWwul0YseOHUob9WvIbeTXiMTlcsHpdGpuRNTJHCuVgh9LKjDwslj3hogSWIsCoHHjxqG4uBgrV67ESy+9hIMHD+KSSy5BTU0NHA4HrFYrMjMzNc/JycmBw+EAADgcDk3wI5+Xz0Vr43Q60dAQOd09b948ZGRkKLf8/PyWXBoRJYJ9a6X7QVcA9ZVSIXTD6dj2iYgSkrklja+6KphuHjlyJMaNG4d+/fph2bJlSE5ObvPOtcSsWbMwc+ZM5bHT6WQQRNTZVEkZZ2SfC/zrBsCxHbjpTeCcSbHtFxElnDOaBp+ZmYlzzjkH+/btQ25uLtxuN6qqqjRtKioqkJubCwDIzc0NmxUmP26qjd1ujxpk2Ww22O12zY2IOpmMfCDvPKDH2YC9t3Sspjy2fSKihHRGAVBtbS3279+PXr16YezYsbBYLFi7dq1yfs+ePSgrK0NBQQEAoKCgANu3b8fx48eVNmvWrIHdbsewYcOUNurXkNvIr0FEXdiEh4C71wMjfgak95KOORkAEVHLtSgAevDBB7FhwwZ8//332LRpE37605/CZDLhxhtvREZGBqZNm4aZM2fik08+wdatW3H77bejoKAA48ePBwBMmjQJw4YNwy233IKvv/4aq1atwuzZszF9+nTYbDYAwD333IMDBw7g4Ycfxu7du/Hiiy9i2bJlmDFjRttfPRElLjkDxACIiFqhRTVAR44cwY033ojKykr07NkTF198MTZv3oyePXsCABYsWACj0YgpU6bA5XKhsLAQL774ovJ8k8mEDz/8EPfeey8KCgqQmpqKqVOnYu7cuUqbAQMGYPny5ZgxYwaee+459OnTB6+++ioKCwvb6JKJKCHJK3YYDNK9nRkgImq9Fq0DlEi4DhBRJ3PqAPBiAdB9EHDvZ9KK0P/4KdBzKDCd22IQdRZxuw4QEVFMOMsBb2Nw9ef0POmeRdBE1AotGgIjIooZ5zHp3h4IfDL6ABfPlB4LERwaIyJqBgZARJQYnEelezkAsqUBEx+LXX+IKKFxCIyIEkNNSAaIiOgMMAAiosQgZ4DSVQFQjUPaDuP0odj0iYgSFgMgIkoMNdJ+gcr0dwD4+AngtULg27di0yciSlisASKixJA9VLrP6BM8ltpduueGqETUQgyAiCgx/OT58GNJmdI9AyAiaiEOgRFR4krOlO4bqmLZCyJKQAyAiChxJXeT7hkAEVELMQAiovh3+hDwhzzg+bHa4/IQWGNVR/eIiBIcAyAiin8NpwFPHeCq1R7nEBgRtRKLoIko/skZHnnIS5bZD5jwKJCW3eFdIqLExgCIiOKfnOGRMz6y1B7A5bM6ujdE1AlwCIyI4p+cAZJrfoiIzhADICKKf/I6P6FDYABwch9QtgVw1XRsn4gooTEAIqL4F2kIDAD+eR3w2iTg+K6O7BERJTgGQEQU/+x5QJ8LgKyB4eeUmWBcDZqImo9F0EQU/8b9j3TTw8UQiagVmAEiosTGxRCJqBUYABFRYuMQGBG1AgMgIop/C8cDz5wLOL4NP6fsCF/VkT0iogTHGiAiin/Oo4DLCZht4eeUGiBmgIio+RgAEVF883ml4AfQXweobwEw4RGg1+gO7RYRJTYGQEQU3xqrg18nZYSf7ztOuhERtQBrgIgovsnZH0sKYLLEti9E1GkwA0RE8c1dK93b0vXPe13AqQOAzw30GtVx/SKihMYMEBHFN3mPL2ua/vnT3wMvjgde/0mHdYmIEh8zQEQU30w2oPf5QEYf/fNyYOSuBYQADIaO6xsRJSwGQEQU3/qMBe5aG/m8LRAA+b3ScJglqWP6RUQJjUNgRJTY1ENjcr0QEVETGAARUWIzmgBLqvS1PGOMiKgJDICIKL5tfFbaBmP9/Mht5GEwFzNARNQ8DICIKL7VnQCcRwB3XeQ28hR5DoERUTOxCJqI4ps8DT7SOkAAMPZ2qZ09r2P6REQJjwEQEcW35gRAF93XMX0hok6DQ2BEFN/kYa1ICyESEbUCM0BEFN/kwmZblACo/hRQe1zaLT49p2P6RUQJjRkgIopvzRkCW/ck8OI44MvXOqZPRJTwmAEioviW1V/a6DSle+Q2cnAkB0tERE1gAERE8e36fzbdxipPg2cARETNwyEwIkp8XAiRiFqIARARJT4OgRFRCzEAIqL41XAaWDAcWDge8Psit5OnyHMlaCJqJtYAEVH8anQC1YcBc5K06WkkzAARUQsxACKi+NXcRRCzBgDjfwlk5Ld/n4ioU2AARETxS1kEMcoaQACQNRCYPK/9+0NEnQZrgIgofgWmtYtoq0ATEbUCAyAiilvfHjgKADjgNERvKATgLAdOfBe9WJqIKIABEBHFrfXbDwAADtU08VHl9wHPDAUWXgA0VndAz4go0TEAIqK45TUlYb+/F46IntEbmsyAySZ9zanwRNQMLIImori1Nf0KPOsYCQC4tanG1hSgwQW469u9X0SU+JgBIqK4ZTG14CPKkirde+rapzNE1KkwACKiuGUxNVH8rGZNke6ZASKiZmAARERxa8rp17DS+gh+ZtrQdGNLIADyMAAioqadUQA0f/58GAwGPPDAA8qxxsZGTJ8+Hd27d0daWhqmTJmCiooKzfPKyspQVFSElJQUZGdn46GHHoLX69W0Wb9+PcaMGQObzYZBgwahuLj4TLpKRAmop9eBIcbDSEczghprYAjMzSEwImpaqwOgL774Aq+88gpGjhypOT5jxgx88MEHePPNN7FhwwaUl5fjuuuuU877fD4UFRXB7XZj06ZNeP3111FcXIw5c+YobQ4ePIiioiJcfvnlKC0txQMPPIA777wTq1atam13iSgBJYkGAEA9kpps+4FnLNbafwrRbUB7d4uIOoFWBUC1tbW4+eab8be//Q3dunVTjldXV2PRokV45plncMUVV2Ds2LFYvHgxNm3ahM2bNwMAVq9ejZ07d+Kf//wnRo8ejauuugpPPvkkFi5cCLfbDQB4+eWXMWDAAPz5z3/G0KFDcd999+FnP/sZFixY0AaXTESJwiYaAQD1wgYhRMR2jR4f7j8wHtOO/xyHkwZ3VPeIKIG1KgCaPn06ioqKMHHiRM3xrVu3wuPxaI4PGTIEffv2RUlJCQCgpKQEI0aMQE5OjtKmsLAQTqcTO3bsUNqEvnZhYaHyGkTUNdj8cgbIBpfXH7FdZZ1b+dpqZmkjETWtxesALV26FF999RW++OKLsHMOhwNWqxWZmZma4zk5OXA4HEobdfAjn5fPRWvjdDrR0NCA5OTksO/tcrngcrmUx06ns6WXRkRxRskAIQmNHh+SLCbddqdq3UiCC3bUQzScAjLyOrKbRJSAWvRfpcOHD+PXv/413njjDSQlNT0m35HmzZuHjIwM5Zafnx/rLhHRGbL4AhkgYUOjJ1oGyIUHzP/B50nTkbbl2Q7qHRElshYFQFu3bsXx48cxZswYmM1mmM1mbNiwAX/5y19gNpuRk5MDt9uNqqoqzfMqKiqQm5sLAMjNzQ2bFSY/bqqN3W7Xzf4AwKxZs1BdXa3cDh8+3JJLI6I4VGvpDofohloko9ETeZPTU3VuNIjAVhicBk9EzdCiAOjKK6/E9u3bUVpaqtzOP/983HzzzcrXFosFa9euVZ6zZ88elJWVoaCgAABQUFCA7du34/jx40qbNWvWwG63Y9iwYUob9WvIbeTX0GOz2WC32zU3Ikpsr5+7CONdC7Ff9EajN3oAVA8pADIwACKiZmhRDVB6ejqGDx+uOZaamoru3bsrx6dNm4aZM2ciKysLdrsd999/PwoKCjB+/HgAwKRJkzBs2DDccssteOqpp+BwODB79mxMnz4dNpv0AXbPPffghRdewMMPP4w77rgD69atw7Jly7B8+fK2uGYiShA+f3DmV4M7cgBUWedGA5gBIqLma/PNUBcsWACj0YgpU6bA5XKhsLAQL774onLeZDLhww8/xL333ouCggKkpqZi6tSpmDt3rtJmwIABWL58OWbMmIHnnnsOffr0wauvvorCwsK27i4RxTG/aup7tBqgU7VueAJDYEYGQETUDAYRbXGNBOZ0OpGRkYHq6moOhxElotoTcLwwGYfqzbjePQeLb78Alw/O1m165+tfwrznA7xsfRb1uRcg5Z6PO7izRNRWOurvd5tngIiI2oTLidzGfUgxSBMfXFGLoF1IBzNARNR8XDGMiOKTuxZAcBuMqENgdW4cET2wzDsB1f0nd0j3iCixMQAiovjkljI59YHanmjT4Cvr3NgveuNh7//g+Hm/7pDuEVFiYwBERPEpsKu7nAFqiBAAub1+1DR6lce+zlnWSERtjAEQEcUnjxwAyRkg/SGw0/XyPmACKWiEsf44wCCIiJrAAIiI4pOcARJyDZB+BqiyVgqA7KjHzqQ7MPJfFwA+t25bIiIZAyAiik8GE6pNWahEOgBEXAn6VGAneDlTBEAJnoiIIuE0eCKKT6Oux293nI3l3xwDALgiDIHJmSEvzHALE6wGH+Bp6LBuElFiYgaIiOKWvxlbYahXi+Z2GETUXAyAiChuabbCiDAEpm4jzxjjEBgRNYUBEBHFp0+fxm+O/BrXGDcCiFwErUoSKWsGMQNERE1hAERE8enEHpzj+hY9DdUAIk+D1x0CczMAIqLoWARNRPEpkMVpQPSVoNUZoLX+85A18DzkpelvmkpEJGMARETxSVkHKBAAeSNkgFQR0ALvzzHswvOR1yun/ftHRAmNQ2BEFJ8CU9mVlaCbMQsMAHx+rgJNRE1jAERE8SmwFUajHABFnAWmfiRg8NRzHSAiahIDICKKT3IGqInd4NUZoN+bX0Phe+cBnz3X/v0jooTGAIiI4pM5GY2woQFWAFFmgfn1ZoFxHSAiio4BEBHFp3s34sbsd/CtGAigmesAcSVoImomBkBEFLfUwY3b1/Q6QPLO8VwHiIiawgCIiOKWenhLCECI8BleQrMVhjwEVtvufSOixMYAiIjij6sWWPwjPFY9B2Z4lcM68Y9m2js3QyWi5uJCiEQUf9x1wKHPMAYGeGFSDvuFgBEGTVPdvcA4BEZETWAGiIjiT2ANIBdsgCrg0VvjUF0DdET0RFnOlcCAS9q7h0SU4JgBIqL4E1gDqNFg0xwOXfUZ0A6LfSPOwobRP8EtBf3bs3dE1AkwA0RE8ScwhCWvAi3TC4B8Ice4EwYRNQcDICKKP8pO8Emaw00NgQGAz+cHvK526xoRdQ4MgIgo/njkDJBVc7ipIbBcVOK2tWOBeX3atXtElPgYABFR/PF7AXNSWAZI6KyFqF4rqBFWGOEHfG7A5w1vTEQUwACIiOLP0KuB2RX4pflxzeHQep/QYw3qmiEP9wMjosgYABFR3PKHrfkTHgCp64JcsMAvf6xxLSAiioIBEBHFLb8/dIZX9K0wAAO8pmTpS64GTURRMAAiovhTugRYcj2u9n+sOay3FUZoUOQxyhuicgiMiCJjAERE8adiB/DdSvQXRzWHddcBCimM9jADRETNwJWgiSj+BFaCVvb2CvDpLAQUukP8oYwLkNlvFGBLb7/+EVHCYwBERPEnkL2pF9p1gJozBLZ64KMYVTik3bpGRJ0Dh8CIKP4EAqA60YytMEKGwLgVBhE1BwMgIoo/gSnsDSJ0JejwpqFBkTJzTC9dREQUwACIiOJPoAaothU1QJMP/B8wtztQsrD9+kdECY8BEBHFH58bgF4NUPSFEIHA4ol+L2eBEVFUDICIKP7cuQb+2Sex3j8aAGAMLAitNwQWuj2Gm+sAEVEzMAAiorjkM5iUbS3MJuk+2krQhkCQ5DYEAiBmgIgoCgZARBSX1PU+5kAKSHcvsMAsMItR+jhzKRkgBkBEFBnXASKi+LNsKswwohsm4TTssJiMAHxKsKMmB0VmkwFuH+BSMkAcAiOiyBgAEVF8EQLY+S7MAEz4IQDAYoqcAZJrgEyBLJESADEDRERRcAiMiOJLYAo8ANRDmgZvijIEJh+yBOqETplzgAGXArkj2rmjRJTImAEioviiCoAaIU2DNxvlIujw5soQWCBI2pN6PvCLae3cSSJKdMwAEVF8CdTuCHOSMgtMHgKLtg6QJcpMMSKiUAyAiCi+BGp3hDkZgDS93RjI7uitBC1vfRFtmIyIKBQDICKKL4H1e4QlBQBgNBhgMsjBTXhz9SwwAOjR8D0wvx/wzLnt31ciSlisASKi+BKoAfIHMkAmgwFGQ7QhMOmYvA6QF2agsUraDoOIKAJmgIgovvS7CJh9AhXXfwQAMBqDqzzrZ4CkezkD1BCYOQZ3HXeEJ6KIGAARUXwxGACzFT5LKgBpCEzOAIXu+wUEa4DkWWCNBnkHeQF4Xe3fXyJKSAyAiCguyQXPJoMBgdEt/a0wlBogqVGjnAECuB8YEUXEGiAiii/7Pga+/jfSup0HoB+MxmARdLRp8HIGyAsTYLICPncgAMrqoI4TUSJhBoiI4svxXcD2ZUhyfAFAmt5ukGeBRdkLTF4HyOcXQGAGGbfDIKJIWhQAvfTSSxg5ciTsdjvsdjsKCgrw0UcfKecbGxsxffp0dO/eHWlpaZgyZQoqKio0r1FWVoaioiKkpKQgOzsbDz30ELxe7WyN9evXY8yYMbDZbBg0aBCKi4tbf4VElFgCQYvPJM0CMxqkG9DUEJhqHaD8cUD/SwCjqQM6TESJqEUBUJ8+fTB//nxs3boVX375Ja644gpcc8012LFjBwBgxowZ+OCDD/Dmm29iw4YNKC8vx3XXXac83+fzoaioCG63G5s2bcLrr7+O4uJizJkzR2lz8OBBFBUV4fLLL0dpaSkeeOAB3HnnnVi1alUbXTIRxTVPaAAULILWDYACWaHgdhkCuHkZcNuHQPezOqDDRJSIWlQDdPXVV2se/+EPf8BLL72EzZs3o0+fPli0aBGWLFmCK664AgCwePFiDB06FJs3b8b48eOxevVq7Ny5Ex9//DFycnIwevRoPPnkk3jkkUfw+OOPw2q14uWXX8aAAQPw5z//GQAwdOhQbNy4EQsWLEBhYWEbXTYRxa1AAKSsA2Q0KCtBN2cvML1hMiKiUK2uAfL5fFi6dCnq6upQUFCArVu3wuPxYOLEiUqbIUOGoG/fvigpKQEAlJSUYMSIEcjJyVHaFBYWwul0KlmkkpISzWvIbeTXiMTlcsHpdGpuRJSAAgGQ16zOAEmnou0GLw+B6U2VJyIK1eIAaPv27UhLS4PNZsM999yDd955B8OGDYPD4YDVakVmZqamfU5ODhwOBwDA4XBogh/5vHwuWhun04mGhgZEMm/ePGRkZCi3/Pz8ll4aEcUDuQbImARAWgjRGGUrDF9IEbQQAnj7buCpgcD2tzqgw0SUiFocAA0ePBilpaXYsmUL7r33XkydOhU7d+5sj761yKxZs1BdXa3cDh8+HOsuEVFrBLbC8LVwKwyzesNUdx1QXwm4mAkmIn0tXgfIarVi0KBBAICxY8fiiy++wHPPPYfrr78ebrcbVVVVmixQRUUFcnNzAQC5ubn4/PPPNa8nzxJTtwmdOVZRUQG73Y7k5OSI/bLZbLDZbBHPE1GC+Hkx4K7DiaMNAL6G0WhQtsLQ3Q0+ZAjMLwBYAp8VnAZPRBGc8TpAfr8fLpcLY8eOhcViwdq1a5Vze/bsQVlZGQoKCgAABQUF2L59O44fP660WbNmDex2O4YNG6a0Ub+G3EZ+DSLq5CxJQGp3eEyBITCDAaZoRdDKVhiqWWDyOkCeyMPmRNS1tSgDNGvWLFx11VXo27cvampqsGTJEqxfvx6rVq1CRkYGpk2bhpkzZyIrKwt2ux33338/CgoKMH78eADApEmTMGzYMNxyyy146qmn4HA4MHv2bEyfPl3J3txzzz144YUX8PDDD+OOO+7AunXrsGzZMixfvrztr56I4pY8m8vU1DR4vXWArNI+YvDUtX9HiSghtSgAOn78OG699VYcO3YMGRkZGDlyJFatWoUf/vCHAIAFCxbAaDRiypQpcLlcKCwsxIsvvqg832Qy4cMPP8S9996LgoICpKamYurUqZg7d67SZsCAAVi+fDlmzJiB5557Dn369MGrr77KKfBEXcXq2YC7Dpa8WwEARmNwFli0rTCCK0EjOATGDBARRdCiAGjRokVRzyclJWHhwoVYuHBhxDb9+vXDihUror7OZZddhm3btrWka0TUWXyzDKitgHHyNQAAkxHKVhg+nTV+5KBIHiYT6iEwNzNARKSPe4ERUXwJZG08xuatAyQXRlvUs8Ay+gC9RgMZXA6DiPRxN3giih9CKFkbrzG8CDrqNHiTqgh61A3SjYgoAmaAiCh++DyA8AEA3MbgVhiGKAshhq4EzYWgiag5GAARUfxQzdrymKSZodJu8M2YBWbkVhhE1HwMgIgofsiztoxmeGGRvlTVAOkthOgTOusAHf4ceHYkUPzj9u8zESUk1gARUfyQV262pCiZHZPRAJMh8vCWvF6QxaTaDV74gapDgIH/xyMifQyAiCh+ZA0AHtwHeBvhPxgMgAxRhsCEXhE0V4ImoiYwACKi+GE0AWk9AQA+/xEA0hpAwWnw4U+Rj5nU0+CVAIh7gRGRPuaHiSguKUNgTRRB+0KKoP0CgJUBEBFFxwwQEcWP8lJg2z+AnkPgh7T9jcloQKC+Wdn4VE1/CCywFYbfC3jdgNna7l0nosTCDBARxY+T3wFfvArs/lDJ7Bg0m6GGP0XZC8yoyhJZUoMNmAUiIh3MABFR/HDXSveWVCXb0/zd4OXNUAVgsgA9zgHMNikLREQUggEQEcUPeRq8NUVT3BxtN3h5bSDNStAGA3DfF+3dWyJKYBwCI6L44QmuAyQHNkbVNHi9VZ6FMgSmygARETWBARARxY/ARqiwpipDW9qtMMKfol4wUf2YiCgaBkBEFD88OitBGwwwybPAotQAKStBy22W3Qo8Nxr4/rN27TIRJSYGQEQUP1Q1QL7AFhdGY7AIOtpWGMFp8IETznLg9EGg4XQ7dpiIEhWLoIkofhT+HrhkJpCUCf/nUuBiNCC4FYbOGFjYbvDKvHhuh0FEkTEAIqL4kdxNugHw+U8B0M4C0yuCDk6DNyjHhBAwKAFQXTt2mIgSFYfAiCguBYugDUqBs+4QWOCY2Rj8OPP5RXA7DDcXQiSicMwAEVH8+OwvQH0lMObW4EKIUXaDVw+JWVQZIL8AN0QloqgYABFR/Nj2T+DkHmDQlfCJbABSBii4G3xIAKR6LBdBK8cZABFRFBwCI6L4oUyDT1WGtoyqrTDkmWEydU20XAQtHRdAeg7QbQBgS2/PHhNRgmIGiIjih7IQYopqCAwRt8LQZIBUAZDPL4BLfiPdiIh0MANERPEjwlYYxgirPEceAmvnfhJRwmMARETxwe8DvI3S19ZUZcq7UbMbfMhTVI81RdCMgIioCQyAiCg+qIuVLSnKlHdTM4ugTaE1QPs+Bl6+GHhvert1mYgSF2uAiCg+KOv1GABLsnYILMJK0OrHJoMBBoO0VpBPCOn1HNsBS2qHdJ+IEgsDICKKDyndgfu2Spkgg0E1BBZ5N3j1Y3mozCeElD3iStBEFAUDICKKDyYz0GOQ8lCodoNvzhCY0WiAyWCAD0K7EjT3AiMiHawBIqK4pDcLLHQrDL8qSwQABnWgZEmWHnArDCLSwQwQEcWHk3uB0iVAt/7A2KnKoofRt8KQ7uUhMrkQ2u9HsPaHK0ETkQ5mgIgoPpzYA2x8Btj2DwDa7I6yG3xoEbRqqrz6XpMBYgBERDqYASKi+KBaBBEI2Q0+YhG0PEyGQFvp3icEYE0FUntKr+f3AUZT+/afiBIKAyAiig/KNhjS0JVPtRu8nNkJ3QpDfqhkgIyqdilZwEP72rvXRJSgOARGRPEhSgbIEGEWmFIoLdcARcgUERGFYgBERPFBnq0VmL6uFDirMkC+JmeBybvGMwIiougYABFRfJAXLLSmAYCyEKLJYFBqfMJ3g5fu5aEveT9UJVO09GbglQnASQ6FEZEWa4CIKD7INUDyEJhSAxQyu0tFRJoFFsgewfENUFUGNJxuz54TUQJiAERE8eHimcDom6SZWwhmgAzq3eD92qf4QobAwgIlrgVERBEwACKi+GDvJd0C5OEtkzoAamIhRHmoTA6MlO0w3NwPjIi0WANERHHJr5kGHzgWYS+w0FlgSq2QsiEqM0BEpMUMEBHFh6/+AdSdAIZdA3Q/S5nJZTAEi5xDJ3cF1wGS7+VZYIEGgYJqZoCIKBQDICKKD1++BpR/BWQPBbqfpWR31Ashhq0DpKoTAtSBEofAiCg6DoERUXyQh6kCK0H71dPglSEw7VPUQRIQzATJw2dI7gakdOc2GEQUhhkgIooP8kKIFu1WGMaoW2FEmgUWaFD0Z+lGRBSCGSAiig/KQojSsJVPqe8JboURvht8sI363hcSKBERhWIARETxwa3dC0wow1vBIa7QITB1obTUVr9WiIgoFAMgIoo9vx/wNkhfh+wGbzREHgJrsgZoz0qg+MfAx0+0Z++JKAGxBoiIYk+9Tk8gA+RTrQMUaTd4EToEFpopqj8JfP9fwJzUPv0mooTFAIiIYs+cBNy1ThoGsyQD0C5yaIywy3twCCykBkhuZ+VWGESkjwEQEcWeyQz0Hqs5pC5wDg6BIaSNdhZY+ErQgQCI6wARUQjWABFRXFJvhWEKfFJFGgKTa4CU2WIiJAPEAIiIQjADRESxV30E2P4WYM8DRv4CgHand0Po+j4B/pCVoMNmi1m5FxgR6WMGiIhi7+Re4OPHgI3PKoeU4a1oW2H4tUNgSjulBkjeC6y2nTpORImqRQHQvHnzcMEFFyA9PR3Z2dm49tprsWfPHk2bxsZGTJ8+Hd27d0daWhqmTJmCiooKTZuysjIUFRUhJSUF2dnZeOihh+D1ejVt1q9fjzFjxsBms2HQoEEoLi5u3RUSUfxTtsFIUQ75AxuaarbCaGohxLC9wFIBkxUw2dqn30SUsFoUAG3YsAHTp0/H5s2bsWbNGng8HkyaNAl1dcHx9RkzZuCDDz7Am2++iQ0bNqC8vBzXXXedct7n86GoqAhutxubNm3C66+/juLiYsyZM0dpc/DgQRQVFeHyyy9HaWkpHnjgAdx5551YtWpVG1wyEcWdkEUQAe00+LAtLgKEar8wIJgJUmaB2fOA350AHtrbTh0nokTVohqglStXah4XFxcjOzsbW7duxaWXXorq6mosWrQIS5YswRVXXAEAWLx4MYYOHYrNmzdj/PjxWL16NXbu3ImPP/4YOTk5GD16NJ588kk88sgjePzxx2G1WvHyyy9jwIAB+POfpT18hg4dio0bN2LBggUoLCxso0snorihbIORqhwK7vQOGKA/BCbHOYaQWWBcCZqImnJGNUDV1dUAgKysLADA1q1b4fF4MHHiRKXNkCFD0LdvX5SUlAAASkpKMGLECOTk5ChtCgsL4XQ6sWPHDqWN+jXkNvJrEFEnI8/SUmWAhGqVZ6MyC0z7NJ9qrSAgcrE0EVGoVs8C8/v9eOCBB/CDH/wAw4cPBwA4HA5YrVZkZmZq2ubk5MDhcCht1MGPfF4+F62N0+lEQ0MDkpOTw/rjcrngcrmUx06ns7WXRkQdzR1eA6QMgRkMkOOZiLvBBwIkebq8ZsHE/9wF1DqAnzwPdOvf1j0nogTV6gzQ9OnT8e2332Lp0qVt2Z9WmzdvHjIyMpRbfn5+rLtERM0lD4FZVENgqlWejaHr+wT4QzJAupuhHtoEHPwUqD/VHj0nogTVqgDovvvuw4cffohPPvkEffr0UY7n5ubC7XajqqpK076iogK5ublKm9BZYfLjptrY7Xbd7A8AzJo1C9XV1crt8OHDrbk0IoqFsbcDt74HnH+7cki9yGHY9PYAeaaYMdJWGEAwq8TFEIlIpUUBkBAC9913H9555x2sW7cOAwYM0JwfO3YsLBYL1q5dqxzbs2cPysrKUFBQAAAoKCjA9u3bcfz4caXNmjVrYLfbMWzYMKWN+jXkNvJr6LHZbLDb7ZobESWIbv2AgZcBPQcrh3yqGV6RtsLwhW6FYdQLgLgfGBGFa1EN0PTp07FkyRK89957SE9PV2p2MjIykJycjIyMDEybNg0zZ85EVlYW7HY77r//fhQUFGD8+PEAgEmTJmHYsGG45ZZb8NRTT8HhcGD27NmYPn06bDZprY577rkHL7zwAh5++GHccccdWLduHZYtW4bly5e38eUTUbxSFjk0AvDrz+5SF0oDEWaBKfuBcTFEIgpqUQD00ksvAQAuu+wyzfHFixfjtttuAwAsWLAARqMRU6ZMgcvlQmFhIV588UWlrclkwocffoh7770XBQUFSE1NxdSpUzF37lylzYABA7B8+XLMmDEDzz33HPr06YNXX32VU+CJOqtdHwC1FcDAy4HuZwHQ1vcIg/R1eA2QdG8IWQjR51c1UvYDYwaIiIJaFACFzsDQk5SUhIULF2LhwoUR2/Tr1w8rVqyI+jqXXXYZtm3b1pLuEVGi2vIK8P1/gSmLVAGQdMpkDM4Ci7QXWOhu8JoMEIfAiEgHN0Mlothz6yyE6A9mgEyBACf0P2F+v3YWmFG3BigFMJgAn7s9ek5ECYoBEBHFnke7FYZ6tpfRAAh5L7CwDFCgjVwDpLcO0NV/AX7yQnC5aCIiMAAionigLIQoZYDUQ1gmo0EJdMK3wghZB0hvCMxoao8eE1GCO6OtMIiI2oQ8QysQAKmLnY3G4EKIQmiHwYK7wQfbAiEZICIiHQyAiCj2lAAoDUBwgUNAyurI09wB7TBYaA2QnAHSzBbbvw5YejOw4al26DgRJSoOgRFRbHndwQJlmxQAaTJABoMyzR2QhrdMIbvDh22FoY6SnMeA3R8C3sZ2uwQiSjwMgIgotowm4JZ3pZlgNmkFd79mCAwwqjJC/mYNgalePxBUwcWFEIkoiAEQEcWW0QScdbnmkDqDo94KA9Buh9GsImgrV4ImonCsASKiuOPza4fA1AGQ+pxfvV0GIq0DlC7dMwAiIhVmgIgotpzHgO8+AtJ7AYOvAqAtdDYaDUqAI53TGwLTZoC86hfgEBgR6WAARESxdXwn8OEMIGe4KgDSbnJqNESYBRZWBB04rskABQIgZoCISIVDYEQUWyFT4IHgEJac0dHWAImwr8OKoNU1QLbAEJjPA/i8bdp1IkpczAARUWzJQ1O28ABIjntUywBp6nvkQMcQWgStzgAldwP+1wGYk7gdBhEpGAARUWzpZIDkBI48BGaIOASmbWfSywAZDIAluY07TUSJjkNgRBRbbp0MkNAOgQHB4Ea7FUbIEJhBZxYYEZEOBkBEFFvyEJg8XR3BAMaoGvuSv4y2FYbZpLMOEACsni1th3FiT1v2nIgSGAMgIootnQxQaHEzEBwG05sGbwgplg7LAO1bK22H4Tzalj0nogTGGiAiiq2xtwP9LwG6D1IO+UKmwQPBYEizEKLSDpr2mq0wANVU+Lo27DgRJTIGQEQUWznDpJuKL2RoCwjWA6lHt+Svo26FAXAxRCIKwyEwIoo7/kAGRx0AGXWCm+B0+cAQmN5WGAAXQySiMMwAEVFs7V4OeBqkYbD0HADhK0EDwSV8/DqzwEyhK0GHZYACBdaumrbuPRElKGaAiCi21j4J/GeatCVGgFwDpN4DTM7u6G+FEWgTqQiaGSAiCsEAiIhiSy5MtgWnwfv94esARRsCM4YuhBgWAKVK96wBIqIADoERUWy5A8NSqpWgQ3d5V3+tDoDkXd/NxiaKoC99ULpZUtq060SUuBgAEVFsRdkLTHchRNUUdyVTZGyqCDq1LXtMRJ0Ah8CIKHa8LsDvkb7WZICaNwTmDQmA5PY+7oRBRE1gAEREsaOuybFG3g0eUG+FET4LzBxSA+QPzQA5tgPv3AusndtWPSeiBMcAiIhiR67/MScDpuCIvP40+PBZYF6fdqgs4hBY3Qng6yXAnpVt2n0iSlysASKi2EnpDvz8dcDv1RzWC4CU7E60DFCkImh5o1U31wEiIgkDICKKHVs6cO61YYd9uitBS/dCpwbIqKwELT+fCyESUXQcAiOiuBO6wKH0dfgQmBzomE2hRdAhAVCSXbpvdGo3EyOiLosBEBHFTuV+4Nv/AEe3ag6HTm8HggXR6uxO6KapEYugbYEASPgAT32bdZ+IEhcDICKKnQOfAG/dAfz3Gc1hZSsMQ/QaICUDFBj7Uoqgw2qAUgGDSfq60dl2/SeihMUAiIhiRw5G5AxNQGhmR/210BkCkzdBVYqgVYslApDSR6wDIiIVFkETUey4AgFQkjYAkoMc/WnwqgyQMlvMqGkfVgQNAPd9IWWCuB0GEYEBEBHFkpyNiZQB0tsKI0oGyBipCBoA0rLboMNE1FlwCIyIYqdRPwOkZHZUs8CC2Z3g+Ja8EGJoBiisCJqIKAQDICKKHZd+DZCIUgTtU2+GGrJnmJwJ0s0AbX0dePeXwP5P2qLnRJTgGAARUexEygDJCyGqhsDMehmg0N3gDVFqgL7fCJS+AVR82yZdJ6LExhogIoqdCQ8D1UeAvDGawz6d3eDlIMerCm78oQshRhsCUy+GSERdHgMgIoqdsy7XPexXiqCDx+S1ftTZnbCtMKIVQcvDbC4GQETEITAiikP+KDVAcuEzoF4IMTQDpPOizAARkQoDICKKDSGAb98G9n4M+DyaUz6drTDMOmv8hLYzRVoJGlAthMgAiIg4BEZEseKuBd66Xfr6t8cAk0U5FTq7C9CvAQouhNiMImhbhnTPAIiIwAwQEcWKPBRlNAOWZM0pOX4xqAIgudBZPQssUgYI0CmE5hAYEakwA0REsaFeA0gV6ADhKzxLX0sPvNGGwFSv4xMCRqhet/8lwMxdQFJGm10CESUuBkBEFBsR1gACgtmbZtcAybPAVAGTzy9gMale1Joi3YiIwCEwIooVJQOUHnbKF20WWDOKoAHtpqlERKEYABFRbDRWS/e28CEpOcZRB0DNmQWmbh9WCO11AytnAe9OB7yuM+4+ESU2BkBEFBuulg2B6a0D5A0URIeuAyS9RsiLGs3A5peA0n8Ggy8i6rJYA0REsdH/EuAnzwPpvcJO6Q2Bhe4FJoQIZooiFEFrGI3ScJvLKdUfpWW32aUQUeJhAEREsdHjbOmmI7gSdPBY6Cww9RCXHBypN0/VWwvIZ7XD5HIyA0REHAIjotiod3vxXulRVNd7ws7pzgIzaWuA1BkevaEyvSLo7+utAIDTp46fafeJKMExACKimFi3+n38Z9nr+Puaz8PO+QL1O0a9GiCdDJAmAIqyGvRxr7TgYlVlxRn2nogSHQMgIoqJYd+9iL9b/4hMx8awc3pbYYTOAosUAMlrAYUGQI0eH077U6Vzdafa4AqIKJExACKimLC6pTqcSl9q2Dn9GiA5AySlhzQBkCE8AxQ6BFZV70GVkL6XqD99pt0nogTX4gDo008/xdVXX428vDwYDAa8++67mvNCCMyZMwe9evVCcnIyJk6ciL1792ranDp1CjfffDPsdjsyMzMxbdo01NbWatp88803uOSSS5CUlIT8/Hw89dRTLb86IopbSb4aAMAJb/jqzHJwY4yyEnTkDJD+EFhVgxvPeH+BcY0vYMeA29rgCogokbU4AKqrq8OoUaOwcOFC3fNPPfUU/vKXv+Dll1/Gli1bkJqaisLCQjQ2Niptbr75ZuzYsQNr1qzBhx9+iE8//RR33323ct7pdGLSpEno168ftm7dij/96U94/PHH8de//rUVl0hE8SjZK60D5PAkh53T3w0+MAvMpw2AjAbtpqmRiqBP13lwEhmoQBZqfJwAS9TVtfhT4KqrrsJVV12le04IgWeffRazZ8/GNddcAwD4+9//jpycHLz77ru44YYbsGvXLqxcuRJffPEFzj//fADA888/jx/96Ed4+umnkZeXhzfeeANutxuvvfYarFYrzj33XJSWluKZZ57RBEpElKB8XqSKOgBAhScp7LRfpwjaEmEWmNmo/X9csAha+5rVDW7l63qX9ww6T0SdQZvWAB08eBAOhwMTJ05UjmVkZGDcuHEoKSkBAJSUlCAzM1MJfgBg4sSJMBqN2LJli9Lm0ksvhdVqVdoUFhZiz549OH1af+ze5XLB6XRqbkQUp1Tr8JQ3hgdA0fYC8wQCIDkTFBL/RBwCO13vwVmGo5ht/gfO3f+3M7wAIkp0bRoAORwOAEBOTo7meE5OjnLO4XAgO1u7AqvZbEZWVpamjd5rqL9HqHnz5iEjI0O55efnn/kFEVH7aJD+I+MUyah2C4iQ4argOkDBY6ErQfubyACFDYHVu5FtqMKd5o8w+PhHbXQhRJSoOs0ssFmzZqG6ulq5HT58ONZdIqJIUrvjd7gH8703wecXaPD4NKf1M0DaGiCvP3ymmNROPwNUXe9BlUgDACR5mSEm6uratBIwNzcXAFBRUYFevYL7+1RUVGD06NFKm+PHtauwer1enDp1Snl+bm4uKiq0C5XJj+U2oWw2G2w2W5tcBxG1L78tE/90XQo5SVPT6EWKNfhx1Jzd4OV7s0n7/zhlHSCdDJAcACV7nYAQgCEkeiKiLqNNM0ADBgxAbm4u1q5dqxxzOp3YsmULCgoKAAAFBQWoqqrC1q1blTbr1q2D3+/HuHHjlDaffvopPJ7gEvlr1qzB4MGD0a1bt7bsMhHFQL3HB3V8UtOoLUqOuhu8P3QWmDaIUYbAQqfB13tQBWkdIDO8gLvuTC+DiBJYiwOg2tpalJaWorS0FIBU+FxaWoqysjIYDAY88MAD+P3vf4/3338f27dvx6233oq8vDxce+21AIChQ4di8uTJuOuuu/D555/js88+w3333YcbbrgBeXl5AICbbroJVqsV06ZNw44dO/Dvf/8bzz33HGbOnNlmF05EsdNwbA8uNX6Nfgappq+mUbsfmO46QKGzwOQMUMgYWMR1gOo9aIANLhHINDVwMUSirqzFQ2BffvklLr/8cuWxHJRMnToVxcXFePjhh1FXV4e7774bVVVVuPjii7Fy5UokJQVnerzxxhu47777cOWVV8JoNGLKlCn4y1/+opzPyMjA6tWrMX36dIwdOxY9evTAnDlzOAWeqJMwfvsW/m59Bm94r8T/eqehNmRaektWgjYZ9TNAoUNgVQ1uAAZUIw3ZqJICoExOliDqqlocAF122WVhMzbUDAYD5s6di7lz50Zsk5WVhSVLlkT9PiNHjsR///vflnaPiBKAr17ai0sekgobAmvGXmDeSAGQvBBiyDpApwO7zleJVGQbqpgBIuriuBwqEXU4UV8FAEpRcnOGwJRZYIFzwWnwIUNgOhkgIQSq6qWFEKd5HkSWPR3v9buoTa6FiBITAyAi6niNUvalCnIApM0ABWa6R88A+cKDJECdAQoGQPVuHzyB9odFDmo8FsBkaZNLIaLE1GnWASKixGFqlBdC1B8Ck4fZ1WscKjVAviYyQDpF0Kfr3Zo2ddwKg6jLYwaIiDqcxSUFQJUiHQDCiqD1prhHqgEKnwYfeA3VEFh1gzTEZjMbMca/HRMNX8FbWgPz6Ovb5HqIKPEwA0REHS7JLRVBVyIDQOQaIP11gAJbYSgLITY9BNYYWGm6R5oNIwwHMM38Efx7P26biyGihMQMEBF1LCGwMn8GduzdD6elB+DWGwKT7jUZIFPzMkB6RdCNHiloSrWZUF2XCQDw155oowsiokTEDBARdSyDAZvSCvFX39XIzMgEoDMEFm0vsCYWQtTbC8zllTJASRYT6syB1eTrGAARdWUMgIiow8kBT15mMgDAGToLTGcILNJeYBFngakyQK5ABshmNqLekiU9r/5kG1wJESUqBkBE1LGqj6L/6c9wluEoemVIK8TXhtQAKQsh6s0CU4bApKAm4jpAqoUQXV45ADLBZZMyQObGSiDKoq5E1LkxACKijnVgPR46ORtzzP9AboaUAYq0FYYhyiywYJDU/CJom9kIt607AMDo9wCNVW1ySUSUeBgAEVHHCtTenIQd3VOtAIJFyjI5e2My6MwCC5yU1wMKDYD0iqDlDFCSxQRrUjKcIiXQFw6DEXVVDICIqGPJAZDIQGaKtBqzXKQs8+vWAEkfV2EZoNB1gAKfanpF0DazESlWM37qfgJvX/EJkDWwTS6JiBIPAyAi6liBrEulsCMzJZgBUm+y7NebBWYKrQFqYghMZxq8zWJEqtWE/aI3ThkyAaOpzS6LiBIL1wEioo4VyABVigxkJgf343J5/UiySAFJcBp88GlhNUARFkIMFkHrZYBMSl1RnUubdSKiroUZICLqWHIABLsyBAYE63SASENgwQyQECLyVhh66wCpMkDJFhMmGL/GRfueBna+32aXRUSJhQEQEXUoERgCOynsSE+yKFkelyeYkVEyQDo1QADgF1EWQjTorAOkmgafbDFhjPE7XOBYChzc0FaXRUQJhgEQEXUo1+WP4Q+em3BYZCPZYoLNLA17qWeCBZb40a0BAqQ1gCIthBjcDT54TD0NPtlqQqWwSye4GjRRl8UaICLqUHXn/BR/80m7wNvMRiRZjGjw+DQzwfRmeKkzPT6/ULJELckAyTVGlULahBXcD4yoy2IAREQdqsEj78tlhNFoCAQlHk0GKJjdCT5PXQ/k9Qv4Iq0DpKwXpD8N3mAAjotM6UTNsTa5JiJKPBwCI6KOU30Exr2rMdBQrmRjbGbpY6hRkwGS7tXBjTob5PMFM0Dh0+ADbfSmwZulIuhjkPYDg7Oc22EQdVEMgIio4+xfh7wVUzHb/E8kBwIgORByqWuAdNYBMhoNSsG0R1UDFLYQojwEFmE3+GSLCRUiEAD5XED9qba6OiJKIAyAiKjjVB8FADhElhIA2SxyEbRqFliEKe7q1aCDCyFqP8aUImjdWWBGJFlN8MCM04ZM6aTz6BlfFhElHgZARNRxAsHGMZGlBD66Q2BNrPLs9QlVG+230M0AKesAmZTAa0byH4CH9gO5I878uogo4bAImog6jrMcAOBAFpItUuQSfQhM+3T1atCRMkB6CyE2qoqg5QBorz8PSO1x5tdERAmJGSAi6jiBAKhcdEeyNVADpJMB8unUAAHa/cB8ETJAukNgnuA0ePn7Nni4FQZRV8YAiIg6jpwBEllIMofWAIUvhBg6BKbOAPkiZYB0i6C1s8AA4Gz3buCjR4Etr5z5dRFRwmEAREQdw1UDuKoBBAKgkAyQ7kKIkWqA/P6ICyHqZ4CCQ2DykFsvXzmw5SVg1wdtcHFElGhYA0REHcNoBq57FZu+3om6Hclh0+DlDJAQwTV+DGE1QMFZYJEWQjQpu8EHj6lXgpaHwI6J7tLJQFaKiLoWZoCIqGNYkoGRP8eW3BsBSCtBA8FZYHKWxu3zK2sTysGRzKTaEb6phRDlITC/X8DtU02DD3w/LoZI1LUxACKiDiUXO4cthBjI0jS4g0NhySEBkG4NkCH6EJhblQqyWUwwm4ywmoyoEN2kg94GLoZI1AVxCIyIOkbZZqChCua6JADBwEfOBMkLIcqzsywmAywm/SnuXp+6CDr6StDqBRblbFOSxQinzwpvai7MdQ7g9EEgtXsbXSgRJQJmgIioY2x+EfjX9Rh8cg0AqPYC064EXe/WZojUTLqzwPQLpeUMkJxZMhmDAZVcB+Sy95eeVLn/DC+OiBINAyAi6hiVBwAA5abeANRDYPIsMO0QWIo1PEFtNqlmgUUIgIwG7UKILtVGqDL5e9en95MOVB1q9WURUWLiEBgRtT8hgFNSAHTE0AuAKgNkiZABsuplgPT2AgtZKygQJMlT6dWrQMvk733w3F+j5zX/B6Rw+Iuoq2EGiIjaX40D8NQBBhOOIBsAkGwNmQUmZ4A8kYfAzOpZYP7g0JZapAyQekaZ/HW1pYe0HUbofHsi6vQYABFR+zsVqLHJ7Is6rxRshK8DFCiCdnsBACm6GSBVDVBg5nroLDBLIAPk8ck1QOEZIPl7czsMoq6LARARtT+5yLj7WUrQYYuwEGK0ITC9DJA85CVLs1kAALUuKZAKboMRfD35tRtdXmD1bGDJDUBd5ZlcIRElGNYAEVH7kzNAWWehoUI7xGUL2Qoj2iwws0muAQoWQYdumJqWJH2s1TZKAVCjEnDpZIC8fmDHu0D1YaByH6fCE3UhDICIqP2ddyuQMwLIGojGb6T9wCJthSEHLHpDYGaddYBC9wJLswUCoJAMUJI5vAaoweMDsgYGAqC9QN9xZ3qlRJQgOARGRO2vxyBg5M+BPmOVACd0IcSwDJDONPjmrAOUnhQaAOlkgAIF2A1uH5AzXDp47JszuUIiSjAMgIioQ4XO8pJrc1yhNUBNzgLTD4BSQzJAjVHWAWr0+IBeo6SDx74+k8siogTDAIiI2pfjW2DT88CRrRBCBDNA1uC2FEBwvZ5oQ2DaWWD6AZA8BOb2+uHy+pRNVm2qgEozCyxvdKCf2wE/Z4URdRUMgIiofe1dLc20KnkBbp8f/pCd3uXaHE+grqc+MA2+qVlgXl/0AAgA6lw+1Sww1UKIVlUGqPsgwJIirVPELTGIugwGQETUvuShpV6j0OgO7syuDIGpanNcXp8yBKafAQrOAvNHyACZjAblubWNXtUQmF4GyA8YTUDuCMCWATiPtP46iSihcBYYEbUvVQBU4/IAAKwmo7IxqXp2VqPHr+wF1lQNkLIVhs4qzmk2M+rdPtS4PKhukL5nRrJFOa8EQIHvhZv+DSRlckVooi6EGSAiaj81DuD0QQAGoNconKx1AwB6pFmVJkajAdZAMNTo8QWLpPUyQIFFD33qafAmnQBItRZQZZ0LANA9Nfg9k9VDYACQ3I3BD1EXwwCIiNrPgfXSfa9RQEoWTtRIwUjPdJummXo/sPpou8HrzAILXQgRANJVM8EqA0FXd1XQJdcfyfVGCiFYCE3URTAAIqL2s3+ddH/WFQAQOQBSDUlFGwLTWwfIbAz/GEtTrQV0sjaQAUoLfs9uKVIwdKrOHXzSJ/8HPDMU2PVBCy6QiBIVAyAiaj9HvpDumwiA5OGpk7Uu1HuangXm8fuV2V0Ws34NEADUNHpRGQhy1ENgvTKSAADHqhshAsXUaHQCNceCQRsRdWoMgIio/fxyC3DbCiD/QgDAidpGAEDPNG0A1CtTDkga0BCYKRZtFlidy6tkb3LtSWHt5A1Raxq9OF0n1x0Fv2e2Xfra5fXjdL1UJC0Hadi/ThoKI6JOjQEQEbUfsxXo/wPALAUcJ2ukYCQ0A9QrIxkAUF7ViIZAXU60vcCOnG5Q2qhnd8nk7TDKqxqU2WLdUoPtbGaTUoh9rFp6LfT/AWBNk/YFO7SpFRdLRImEARARtT2fV7qFOFGrPwSWlxHMANV7mq4B+v5knfS8zGQYIkyDB4DvK6V26UlmzTpAQDDoclRLWSlYU4Hh10lff/X3Ji6QiBIdAyAianu73gcWnAuULNQcjlQDlBsIgL4/Wa+MPunVAOUEhru+r6wHAPTOTNb99vJ+YIcC7XqEDLmpv2e5HAABwJip0v3Od4GG07qvTUSdAwMgImpbPi/EhqeAWocmiBBCBAOgNG3dTl4gkDlwslY5ppcBGpybpvu8UPIssLJTUgCkLoCWyYXQDnkIDAB6j4Wr+1DA2wjvZ8/rXx8RdQoMgIioTZWtWQjDiV04LdJQ7P+RcrzOHVzksEe6NiCRgxF5oUSryQizKfzjaVB2uuZx78zwAmgguA6QTL0GUPB7SsHTMVUGqLrRi0dOX4M3vFfif767AFX17rDnEVHnENcB0MKFC9G/f38kJSVh3Lhx+Pzzz2PdJSKKwnngC/TY/AcAwJ+9P8cfNxxHeZWUYZGzP6lWU9gih3IwItMb/gKk7SzkYAmIkgEKCYCyUsOHwJSp8FXBAGjBmu/wbv1I/K93GtaW+fG793bovj4RJb64DYD+/e9/Y+bMmXjsscfw1VdfYdSoUSgsLMTx48dj3TUi0nOoBMYlv0AKXPjKPBrf9Z6CBo8P8z7aDSBy/Q8gBTyZKcFZWnozwGRn5wSzQE0Ngcl66GSA5Bogh1MKgPYdr8U/Nh8CADw8eTAMBuDDr4+gfPVzQC0/d4g6m7gNgJ555hncdddduP322zFs2DC8/PLLSElJwWuvvRbrrhFRKFctfH//KdK8p7DD3w+N1y3GY9eOhNEAfPB1ObYcqMSR01I9jl4ABGizQJEyQAAwOCdYBxSpCDo0AxStBuhYdQOEEHh61R74/AITh+bgl5cNwjWj8jDP/CryNs2BWDgOWD8fOLGHawQRdRJxuRu82+3G1q1bMWvWLOWY0WjExIkTUVJSovscl8sFl8ulPHY6nQCAL/56H1KTwz/8tva4BieSBwIAetftxMhTKyP2p7R7EY6lDAYA5NbvwXmVyyO23Z41CUdShwMAejYcwAUn34nYdkfmFTiUfh4AIKvxMMYf/3fEtrszL8UBu7SYXKbrGH5Q8c+IbfdmXITvMn4AAEjznMSEY4sjtj2QfgF2dbsMAJDsrcaV5a9EbHsobTS2Z00CAFh9dZh09IWIbY+knovS7j8GAJj8bvzo8DMR25anDMHWntcqj68+NC9i2+PJA7El+3rlcVHZn2ASHt22lbZ8fJZ7i/J48uFnYfPX6batsvbChl53KI9/eOQFpHirdNvWWHpgXe97lMdXHH0Zds9J3bb1ZjtW9/mV8njCsUXo5irXbes2pWBF/m+Uxxc7/o6ejd/rtvUZLHi/X/D3o6BiCXIb9uq2BYB3+s1RNvu88Pib6FMXeWjn/X6z4DVKQcrYE++iX22p5rzF34gUXzVsvnq8MuQ1CIMJVfVuXOwZD4PwYOM5j2LBMOl368YL++KNLWWYuexr1DRK79OQXLvu983LSMKuY9LvbXMyQAZDcFZYqPSk0Bqg8KArx54Ek9GARo8fP3u5BFsPnYbRIGV/AOA3kwZj+s6rMcq/H0MbDgPr5wHr58FjsKHK1gvVlmz866w/wWeUMldjTr6PXvV7AAACqqn5gZ/7qt73w2eUPotGnFqlvAfBtsHnrMu7C25TKgBg2Ol16Fv7TcSfx4Zet6PBnAEAGFz1Xwyo2Rqx7cbc/4daSw8AwKDqzRjk3ByxbUn29ai29QIADKj5EoOrNkZs+3nPKTiVlA8A6Fv7NYad/iS8UeDyvurxE+WzN69uJ0acWh3xdb/ufhUcgc/enPrvMLpyRcS232b9EEdTzwUA9Gg4iLEn34vYdlfmZShLHw0A6NZ4BBeeeCti2+8yLsZB+/kAALvbgYKKpRHb7rePw76MAgBAqqcSFzv+EbHt9+ljsCfzUgBAkteJCcci/+f+cNpI7OwmLdRp8dXjivK/RmxbnjIU27sXAgCMfg9+eHRhxLYVyYNQ2uPHyuPCw89GbHsyqb/mc3rikYURP3tP23rj8+yfK48vP/oKrP4G3bZOazZKcm5SHl9a/hr8tfqfp20tLgOgkydPwufzIScnR3M8JycHu3fv1n3OvHnz8MQTT4Qdv6DyXdht4euEvHK0Pz7xSwmwKcYvcZc18i/A34/mYblf+gD9kXEb7onS9q3yLLzpk/6Hepnxa0y3/idi2w/K0/EPXzcAwDjDLvzK9nbEtmuP2VDsywYAjDTsxwO2yIHVJocBxd7eAIBBhiN4MErbrxxuFHv7AwB64wQeSXo3YttdjloUe88BAGTBid8mRf6AOVRxCsV7pEAwCS78Lun9iG0dxx0o3jtaefx4UuS9mD7xjULxvnHK44dtK5BicOm23eIfguIDlyiP77OtRA+DU7ft1/6BKD54hfJ4mnUN8o0ndNvu9fdG8aHJyuObrGtxjvGobtsjogeKy36iPL7W+glGGw/otq0U6Sg+PEV5XGj5FKNMO3XbNggrrjt6o/L4UstGjDKV6rYFgGvKb4X8F+h8yyaMNG2J2Pb6YzegAVJgca55M0aaP43Y9uPNX+KwkH5Pl+N2FJzTC3+9Yaxy/jeTBmPVjgocDdQBnd+vmxJghDorOw1rd0tDTdefnx/xew7Pk/7Y53dLgdWsn8TOsSchM8WCqnoPbGYjRvTOCGuTZDHhkcmDMf+j3dh6SJqtdmtBf5wTCLDys1Lw2LQpuGlxHi5xf4afmTZgnHE3bHChZ+P36Nn4PV7bfBgikEg/3/IJLojyc73tcBHqAz/XP5k/wbgoP9f7Dl+Ok5D6/IR5A8ab10Rs++CRi3BESJ8Nj5o3Yrw58u/PnCNj8Z2Qfra/Nn2G/2eJ/Pk0/8hwlArpd+su02ZMtSyL2Pa5o2ejJLCJ7P8zfY47orT965G+WKf57I3c9p9HcvBh4LP3KuM23BOl7X+OZuJNnxQ0Sp+9kduuOJqMv/syAcifvZHbrj9qxGKfFDSOMBzAjChtt5R7sdgrBY1nGY7iwShtvymvw2JvXwDyZ2/ktvvKK7HYKwWN3eDEb6O0/Y/vEizeLf2O2eDG76K0XeG7EIsDn9MA8FiUtut9o7D4u9HK4wdtbyE1ymfv4r0XKI9/aXsHPSN89n7jH4DF+y5SHt9hfQ8Zno4ZcjYIEX/53PLycvTu3RubNm1CQUGBcvzhhx/Ghg0bsGVL+IeMXgYoPz8fa567B6nJ4f/729nzRzid3A8AkF23B2dX6vyPJWBP94k4mToIANCjfj8Gn/w4Ytu9WRNwPG0IAKBbwyEMPbkqYtsDmRfBkS7947M3lmP4iQ8jtv0+YxzK7aMAAGmu4xhxPHLwccQ+BoczpD9CyZ7TGO2I/CFXnj4ChzKlgMLmrcF5jsi/ABWpQ3Cwm5RZMvsacf6xNyK2PZFyNvZnSf+7Mfq9uKA88sJyp5L7Y2/3YPAx7shrUP9vWK0qqQ/29Pih8viCo/+AAfq7d9dYc7Cr51XK4/OO/RsWf6Nu2zpLd+zIDv5PaJTjP7D5anXbNpgzsD3nWuXx8OPvI9lTpdvWbUrF17nBoGbYiRVIdVfqtvUabdjW6xfK48En1yDdXaHb1m8w4atewQDo7MpPkBEhswQAX/a6SclEDDz9X3RrKIvYdlvuL+APZDX6V5Wge/3BkH5a0WjOQIMlA8fSRsBjSkay1YRheXZcMqhH2OytEzUufLLnOKrrPbhpXF9ljZ5QzkYPPtl9HD8Y1EN33R61d7YdwYAeaRidnxmxzak6N8qrGpDfLQUZKeGrRct2ljux+UAlhvfOwAX9u4UtrOiobsSG747j6OkGwOeB3eVAhqscqZ6T2NmzSGl3TuXH6FGvH9wCwJbetykZoLMrP0F23R4A0sevASKsrceUAgAYeOq/yKvdHvF1v+x1MxotUrDU/3QJ+tRsi9j2q9zrUW/tDgDIr/4C/aq/iNj265wpqLFJwW1v5zYMqNLPvgPA9uxrUJ0k/acrt+ZbDDqtDe7UV7ez51U4ldwfAJBduxvnnIr+2Xsi9WwA0mfvkJORs0V7sy5DRdpQAEC3hjKceyJypv5At4tRnj4CAJDRWB718/T7zPE4Ypcy9WmuCoyuiPwf1cP2sTiUKWXqUzynMOZY5Kz+0fRRONhN+sNv8zpxQXnkz1NH2jDsy5oAALD4GjDuaHHEtsdTz8F33a8EIH32XnTkbxHbViYPxK6ehcrji8teitj2dFK+5jPyosOvwhghA+S09cI3qs/IcUeKYYmQAaq19kRp7s+UxxeU/xPempP44a9fQXV1Nex2/YxxW4jLAMjtdiMlJQVvvfUWrr32WuX41KlTUVVVhffei/yPVeZ0OpGRkdHuP0AiIiJqOx319zsui6CtVivGjh2LtWvXKsf8fj/Wrl2ryQgRERERtUZc1gABwMyZMzF16lScf/75uPDCC/Hss8+irq4Ot99+e6y7RkRERAkubgOg66+/HidOnMCcOXPgcDgwevRorFy5MqwwmoiIiKil4rIGqC2wBoiIiCjxdOkaICIiIqL2xACIiIiIuhwGQERERNTlMAAiIiKiLocBEBEREXU5DICIiIioy2EARERERF0OAyAiIiLqchgAERERUZcTt1thnCl5gWun0xnjnhAREVFzyX+323ujik4bAFVWVgIA8vPzY9wTIiIiaqnKykpkZGS02+t32gAoKysLAFBWVtauP8B443Q6kZ+fj8OHD3epPdB43bzuroDXzevuCqqrq9G3b1/l73h76bQBkNEolTdlZGR0qX84MrvdzuvuQnjdXQuvu2vpqtct/x1vt9dv11cnIiIiikMMgIiIiKjL6bQBkM1mw2OPPQabzRbrrnQoXjevuyvgdfO6uwJed/tet0G09zwzIiIiojjTaTNARERERJEwACIiIqIuhwEQERERdTkMgIiIiKjLSdgA6A9/+AMuuugipKSkIDMzs1nPEUJgzpw56NWrF5KTkzFx4kTs3btX0+bUqVO4+eabYbfbkZmZiWnTpqG2trYdrqB1Wtq/77//HgaDQff25ptvKu30zi9durQjLqlZWvO+XHbZZWHXdM8992jalJWVoaioCCkpKcjOzsZDDz0Er9fbnpfSIi297lOnTuH+++/H4MGDkZycjL59++JXv/oVqqurNe3i7f1euHAh+vfvj6SkJIwbNw6ff/551PZvvvkmhgwZgqSkJIwYMQIrVqzQnG/O73o8aMl1/+1vf8Mll1yCbt26oVu3bpg4cWJY+9tuuy3sfZ08eXJ7X0aLteS6i4uLw64pKSlJ06Yzvt96n18GgwFFRUVKm0R4vz/99FNcffXVyMvLg8FgwLvvvtvkc9avX48xY8bAZrNh0KBBKC4uDmvT0s8MXSJBzZkzRzzzzDNi5syZIiMjo1nPmT9/vsjIyBDvvvuu+Prrr8VPfvITMWDAANHQ0KC0mTx5shg1apTYvHmz+O9//ysGDRokbrzxxna6ipZraf+8Xq84duyY5vbEE0+ItLQ0UVNTo7QDIBYvXqxpp/65xFpr3pcJEyaIu+66S3NN1dXVynmv1yuGDx8uJk6cKLZt2yZWrFghevToIWbNmtXel9NsLb3u7du3i+uuu068//77Yt++fWLt2rXi7LPPFlOmTNG0i6f3e+nSpcJqtYrXXntN7NixQ9x1110iMzNTVFRU6Lb/7LPPhMlkEk899ZTYuXOnmD17trBYLGL79u1Km+b8rsdaS6/7pptuEgsXLhTbtm0Tu3btErfddpvIyMgQR44cUdpMnTpVTJ48WfO+njp1qqMuqVlaet2LFy8Wdrtdc00Oh0PTpjO+35WVlZpr/vbbb4XJZBKLFy9W2iTC+71ixQrxv//7v+Ltt98WAMQ777wTtf2BAwdESkqKmDlzpti5c6d4/vnnhclkEitXrlTatPRnGUnCBkCyxYsXNysA8vv9Ijc3V/zpT39SjlVVVQmbzSb+9a9/CSGE2LlzpwAgvvjiC6XNRx99JAwGgzh69Gib972l2qp/o0ePFnfccYfmWHP+YcZKa697woQJ4te//nXE8ytWrBBGo1HzYfrSSy8Ju90uXC5Xm/T9TLTV+71s2TJhtVqFx+NRjsXT+33hhReK6dOnK499Pp/Iy8sT8+bN023/i1/8QhQVFWmOjRs3TvzP//yPEKJ5v+vxoKXXHcrr9Yr09HTx+uuvK8emTp0qrrnmmrbuaptq6XU39RnfVd7vBQsWiPT0dFFbW6scS4T3W605nzsPP/ywOPfcczXHrr/+elFYWKg8PtOfpSxhh8Ba6uDBg3A4HJg4caJyLCMjA+PGjUNJSQkAoKSkBJmZmTj//POVNhMnToTRaMSWLVs6vM+h2qJ/W7duRWlpKaZNmxZ2bvr06ejRowcuvPBCvPbaaxBxskTUmVz3G2+8gR49emD48OGYNWsW6uvrNa87YsQI5OTkKMcKCwvhdDqxY8eOtr+QFmqrf4/V1dWw2+0wm7Vb/8XD++12u7F161bN76XRaMTEiROV38tQJSUlmvaA9L7J7Zvzux5rrbnuUPX19fB4PGEbRq5fvx7Z2dkYPHgw7r33XlRWVrZp389Ea6+7trYW/fr1Q35+Pq655hrN72dXeb8XLVqEG264AampqZrj8fx+t0ZTv99t8bOUddrNUEM5HA4A0Pyxkx/L5xwOB7KzszXnzWYzsrKylDax1Bb9W7RoEYYOHYqLLrpIc3zu3Lm44oorkJKSgtWrV+OXv/wlamtr8atf/arN+t9arb3um266Cf369UNeXh6++eYbPPLII9izZw/efvtt5XX1/j3I52KtLd7vkydP4sknn8Tdd9+tOR4v7/fJkyfh8/l034fdu3frPifS+6b+PZaPRWoTa6257lCPPPII8vLyNH8IJk+ejOuuuw4DBgzA/v378dvf/hZXXXUVSkpKYDKZ2vQaWqM11z148GC89tprGDlyJKqrq/H000/joosuwo4dO9CnT58u8X5//vnn+Pbbb7Fo0SLN8Xh/v1sj0u+30+lEQ0MDTp8+fca/O7K4CoAeffRR/PGPf4zaZteuXRgyZEgH9ahjNPe6z1RDQwOWLFmC3/3ud2Hn1MfOO+881NXV4U9/+lO7/kFs7+tW/9EfMWIEevXqhSuvvBL79+/HWWed1erXPVMd9X47nU4UFRVh2LBhePzxxzXnYvF+U9uZP38+li5divXr12sKgm+44Qbl6xEjRmDkyJE466yzsH79elx55ZWx6OoZKygoQEFBgfL4oosuwtChQ/HKK6/gySefjGHPOs6iRYswYsQIXHjhhZrjnfH97khxFQD95je/wW233Ra1zcCBA1v12rm5uQCAiooK9OrVSzleUVGB0aNHK22OHz+ueZ7X68WpU6eU57eH5l73mfbvrbfeQn19PW699dYm244bNw5PPvkkXC5Xu+3H0lHXLRs3bhwAYN++fTjrrLOQm5sbNnOgoqICABL+/a6pqcHkyZORnp6Od955BxaLJWr7jni/9fTo0QMmk0n5ucsqKioiXmNubm7U9s35XY+11ly37Omnn8b8+fPx8ccfY+TIkVHbDhw4ED169MC+ffvi4g/imVy3zGKx4LzzzsO+ffsAdP73u66uDkuXLsXcuXOb/D7x9n63RqTfb7vdjuTkZJhMpjP+N6RoUcVQHGppEfTTTz+tHKuurtYtgv7yyy+VNqtWrYq7IujW9m/ChAlhs4Ei+f3vfy+6devW6r62pbZ6XzZu3CgAiK+//loIESyCVs8ceOWVV4TdbheNjY1tdwGt1Nrrrq6uFuPHjxcTJkwQdXV1zfpesXy/L7zwQnHfffcpj30+n+jdu3fUIugf//jHmmMFBQVhRdDRftfjQUuvWwgh/vjHPwq73S5KSkqa9T0OHz4sDAaDeO+99864v22lNdet5vV6xeDBg8WMGTOEEJ37/RZC+htns9nEyZMnm/we8fh+q6GZRdDDhw/XHLvxxhvDiqDP5N+Q0p8WtY4jhw4dEtu2bVOmdG/btk1s27ZNM7V78ODB4u2331Yez58/X2RmZor33ntPfPPNN+Kaa67RnQZ/3nnniS1btoiNGzeKs88+O+6mwUfr35EjR8TgwYPFli1bNM/bu3evMBgM4qOPPgp7zffff1/87W9/E9u3bxd79+4VL774okhJSRFz5sxp9+tprpZe9759+8TcuXPFl19+KQ4ePCjee+89MXDgQHHppZcqz5GnwU+aNEmUlpaKlStXip49e8bdNPiWXHd1dbUYN26cGDFihNi3b59meqzX6xVCxN/7vXTpUmGz2URxcbHYuXOnuPvuu0VmZqYyO++WW24Rjz76qNL+s88+E2azWTz99NNi165d4rHHHtOdBt/U73qstfS658+fL6xWq3jrrbc076v8mVdTUyMefPBBUVJSIg4ePCg+/vhjMWbMGHH22WfHRUAva+l1P/HEE2LVqlVi//79YuvWreKGG24QSUlJYseOHUqbzvh+yy6++GJx/fXXhx1PlPe7pqZG+fsMQDzzzDNi27Zt4tChQ0IIIR599FFxyy23KO3lafAPPfSQ2LVrl1i4cKHuNPhoP8vmStgAaOrUqQJA2O2TTz5R2iCw1onM7/eL3/3udyInJ0fYbDZx5ZVXij179mhet7KyUtx4440iLS1N2O12cfvtt2uCqlhrqn8HDx4M+zkIIcSsWbNEfn6+8Pl8Ya/50UcfidGjR4u0tDSRmpoqRo0aJV5++WXdtrHS0usuKysTl156qcjKyhI2m00MGjRIPPTQQ5p1gIQQ4vvvvxdXXXWVSE5OFj169BC/+c1vNNPFY62l1/3JJ5/o/l4AEAcPHhRCxOf7/fzzz4u+ffsKq9UqLrzwQrF582bl3IQJE8TUqVM17ZctWybOOeccYbVaxbnnniuWL1+uOd+c3/V40JLr7tevn+77+thjjwkhhKivrxeTJk0SPXv2FBaLRfTr10/cddddLf6j0BFact0PPPCA0jYnJ0f86Ec/El999ZXm9Trj+y2EELt37xYAxOrVq8NeK1He70ifSfK1Tp06VUyYMCHsOaNHjxZWq1UMHDhQ83dcFu1n2VwGIeJkrjMRERFRB+ky6wARERERyRgAERERUZfDAIiIiIi6HAZARERE1OUwACIiIqIuhwEQERERdTkMgIiIiKjLYQBEREREXQ4DICIiIupyGAARERFRl8MAiIiIiLocBkBERETU5fx/JrYQTExWI1wAAAAASUVORK5CYII=", "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 our Bragg computation function as we did before, but we've modified it to allow for passing in a list of detuning values." ] }, { "cell_type": "code", "execution_count": 6, "id": "f343ce53-5fcd-4059-8fe2-240581d1c741", "metadata": {}, "outputs": [], "source": [ "from numba import jit, float64\n", "from mwave.integrate import bloch_rhs, omega_fnc_gaussian, phase_fnc_constant\n", "from scipy.integrate import solve_ivp\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(kvec, phi0, tfinal, delta, omega, sigma, omega_mod=None, mod_phase=0.0, phase=0.0, method='DOP853', atol=1e-10, rtol=1e-10, dense=False, max_step=0.1):\n", " \n", " # Compute t0 from tfinal\n", " t0 = tfinal/2\n", "\n", " # Convert passed arguments to floats\n", " tfinal = np.float64(tfinal)\n", " delta = np.float64(delta)\n", " omega = np.float64(omega)\n", " sigma = np.float64(sigma)\n", "\n", " # Determine if we should do multifrequency or single frequency\n", " if omega_mod is None:\n", " # Integrate and return\n", " return solve_ivp(lambda *x: bloch_rhs(x[0], x[1], kvec, delta, omega_fnc_gaussian, np.array([omega, sigma, tfinal/2]), phase_fnc_constant, np.array([phase])), [0, tfinal], phi0, method=method, atol=atol, rtol=rtol, dense_output=dense, max_step=max_step)\n", " else:\n", " # Integrate and return\n", " return solve_ivp(lambda *x: bloch_rhs(x[0], x[1], kvec, delta, multi_omega_fnc, np.array([omega, sigma, tfinal/2, omega_mod, mod_phase]), phase_fnc_constant, np.array([phase])), [0, tfinal], phi0, method=method, atol=atol, rtol=rtol, dense_output=dense, max_step=max_step)\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", " phases = deltas*t\n", "\n", " if idx == 3:\n", " phases += dphase\n", "\n", " soly = np.full((len(deltashift), len(ifr.kvec)), np.nan, dtype=np.complex128)\n", " for i in range(len(deltashift)):\n", " if idx == 1 or idx == 2:\n", " sol = gbragg(ifr.kvec, make_phi(ifr.kvec, k_init/2), 6*sigma, deltas[i], omega, sigma, phase=float(phases[i]))\n", " soly[i,:] = sol.y[:,-1]\n", " elif idx == 3 or idx == 4:\n", " sol = gbragg(ifr.kvec, make_phi(ifr.kvec, k_init/2), 6*sigma, deltas[i], omega, sigma, phase=float(phases[i]), omega_mod=omega_m, mod_phase=omega_m*t)\n", " soly[i,:] = sol.y[:,-1]\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/0k/8rf6137x12n9xmg_yll_84gw0000gn/T/ipykernel_6071/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/0k/8rf6137x12n9xmg_yll_84gw0000gn/T/ipykernel_6071/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": "lab", "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.11.0" } }, "nbformat": 4, "nbformat_minor": 5 }