Gaussian.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "yDOMJ1RpdP8F"
   },
   "source": [
    "# Optimization of a State-to-State Transfer for a Quantum Charger-Battery Driven by Laser Field Using Krotov's Method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "tp3TsOLZlINM"
   },
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pip! install drill"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "id": "YwFSPufgdP8T"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "D:\\ProgramData\\Anaconda3\\lib\\site-packages\\qutip\\__init__.py:93: UserWarning: Old version of Cython detected: needed 0.29.20, got 0.29.15.\n",
      "  .format(_cy_require, _Cython.__version__)\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "import qutip\n",
    "import numpy as np\n",
    "import scipy\n",
    "import matplotlib\n",
    "import matplotlib.pylab as plt\n",
    "from src import krotov\n",
    "\n",
    "\n",
    "import math\n",
    "from src.krotov import Integrals as inte\n",
    "import scipy.integrate as integrate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "MO1oqm_IdP8V"
   },
   "source": [
    "$\\newcommand{tr}[0]{\\operatorname{tr}}\n",
    "\\newcommand{diag}[0]{\\operatorname{diag}}\n",
    "\\newcommand{abs}[0]{\\operatorname{abs}}\n",
    "\\newcommand{pop}[0]{\\operatorname{pop}}\n",
    "\\newcommand{aux}[0]{\\text{aux}}\n",
    "\\newcommand{opt}[0]{\\text{opt}}\n",
    "\\newcommand{tgt}[0]{\\text{tgt}}\n",
    "\\newcommand{init}[0]{\\text{init}}\n",
    "\\newcommand{lab}[0]{\\text{lab}}\n",
    "\\newcommand{rwa}[0]{\\text{rwa}}\n",
    "\\newcommand{bra}[1]{\\langle#1\\vert}\n",
    "\\newcommand{ket}[1]{\\vert#1\\rangle}\n",
    "\\newcommand{Bra}[1]{\\left\\langle#1\\right\\vert}\n",
    "\\newcommand{Ket}[1]{\\left\\vert#1\\right\\rangle}\n",
    "\\newcommand{Braket}[2]{\\left\\langle #1\\vphantom{#2}\\mid{#2}\\vphantom{#1}\\right\\rangle}\n",
    "\\newcommand{op}[1]{\\hat{#1}}\n",
    "\\newcommand{Op}[1]{\\hat{#1}}\n",
    "\\newcommand{dd}[0]{\\,\\text{d}}\n",
    "\\newcommand{Liouville}[0]{\\mathcal{L}}\n",
    "\\newcommand{DynMap}[0]{\\mathcal{E}}\n",
    "\\newcommand{identity}[0]{\\mathbf{1}}\n",
    "\\newcommand{Norm}[1]{\\lVert#1\\rVert}\n",
    "\\newcommand{Abs}[1]{\\left\\vert#1\\right\\vert}\n",
    "\\newcommand{avg}[1]{\\langle#1\\rangle}\n",
    "\\newcommand{Avg}[1]{\\left\\langle#1\\right\\rangle}\n",
    "\\newcommand{AbsSq}[1]{\\left\\vert#1\\right\\vert^2}\n",
    "\\newcommand{Re}[0]{\\operatorname{Re}}\n",
    "\\newcommand{Im}[0]{\\operatorname{Im}}$\n",
    "\n",
    "This first example illustrates the basic use of the `krotov` package by solving\n",
    "a simple canonical optimization problem: the transfer of population in a two\n",
    "level system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "jUiMagksdP8W"
   },
   "source": [
    "## Hamiltonian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import dill \n",
    "dill.settings['recurse'] = True\n",
    "\n",
    "H0=dill.load(open(\"H_0\", \"rb\"))\n",
    "He=dill.load(open(\"H_e\",\"rb\"))\n",
    "Hde=dill.load(open(\"H_de\",\"rb\"))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "g=0.2, T=100,target= np.sqrt(2), np.sqrt(2),0,5 * np.sqrt(2), nt=100 it is good"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter.      J_T    ∫gₐ(t)dt          J       ΔJ_T         ΔJ  secs\n",
      "0     1.00e+00    0.00e+00   1.00e+00        n/a        n/a     3\n",
      "1     1.00e+00    6.56e+00   7.56e+00   3.35e-04   6.56e+00     7 **\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:matplotlib.legend:No handles with labels found to put in legend.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQMAAAEvCAYAAABBv6SGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3hUZfr/8feTDgRCh0AIiZQgoJQEBFEXZV0Vd1eQJrqKCwoCimBBYlt1Y0HEBVRQ1hVQqqwiKCqyVkR/kqC0AIYOgQBJIAVC6ty/PzLJF0iAJFPOTHK/rmuuzMwzc859GOYz5znlOUZEUEopH6sLUEp5Bg0DpRSgYaCUstMwUEoBGgZKKTsNA6UUAH5WzLRx48YSERFhxayVqtE2btyYJiJNymuzJAwiIiJISEiwYtZK1WjGmAMXanNKN8EYc7Mx5ndjzG5jzBRnTFMp5V4Oh4Exxhd4C7gF6AgMN8Z0dHS6Sin3ckY3oSewW0T2AhhjlgK3AdsdmWhaWhq///47oaGh5zwfHBxM06ZNAdi3bx/nH05dt25dmjRpgoiwb9++MtMNCQmhUaNG2Gw29u/fX6a9QYMGNGjQgMLCQg4ePFimvWHDhtSvX5/8/HySk5PLtDdu3Jh69eqRl5fH4cOHy7Q3bdqU4OBgzpw5Q0pKSpn2Zs2aUadOHU6fPs2xY8fKtIeGhlKrVi2ys7NJTU0t096yZUsCAwPJysoiLS2tTHurVq3w9/cnIyODEydOlGkPDw/Hz8+PEydOkJGRUaY9IiICHx8f0tPTyczMLNN+2WWXAZCamkp2dvY5bcYYIiMjATh27BinT58+p93X15fWrVsDcPToUXJycs5p9/f3p1WrVgAcOXKE3Nzcc9oDAgIICwsDIDk5mfz8/HPag4KCaNGiBQAHDx6ksLDwnPbatWvTvHlzAA4cOEBRUdE57Z76f+/QoUNcccUVNGzYsMx7KkVEHLoBg4F3z3p8N/Dmxd4THR0tl3LrrbeKj4+PAOfc7rzzztLX1KlTp0z7mDFjRETEZrOVaQPk0UcfFRGRzMzMctufe+45ERFJTk4ut3369OkiIrJjx45y2+fOnSsiIvHx8eW2L1myREREvvnmm3LbP/30UxERWblyZbnt3333nYiILFy4sNz2jRs3iojInDlzym1PSkoSEZFp06aV256SkiIiIs8++2y57dnZ2SIiMmnSpHLbS9x3331l2urWrVvaPmzYsDLtLVq0OOfzP789KiqqtP26664r0372/6tu3bqVae/bt29pe9u2bcu0/+Uvfyltb9asmVf93+vRo0c536KygAS5wPfSGWsGppznypz9ZIwZDYyG4l+fS5k2bRobNmwgMzOT8ePH07VrV4DSXxaAuXPnlkn39u3bl95fsGBBmel27FjcgwkKCiq3vUuXLkBxSpfXHhMTAxT/QpfX3qtXL6D4F7S89quuugqAyy+/vNz2kuXs1q1bue1RUVEA9O7du9z2kl/W66+/vtz2Zs2aAXDLLbeU/sqdLSQkBICBAwfSpk2bMu2BgYEADB8+vLTW8owcOZJrr732nOf8/f1L748bN47+/fuf016rVq3S+5MmTWLo0KHntNerV6/0/pNPPsmoUaPOaW/UqFHp/RdeeKHMmk/JsgO8+uqrZdZcStYqAGbNmlVmzcNT/+899NBDdO7cuczrK+1CKVHRG9AbWHPW41gg9mLvqciagYhISkqKdO/eXYKDgyUjI6NC71GqJjl27JgA8vrrr1fo9VxkzcAZexPigXbGmEhjTABwB7DKCdOlefPmzJ07l1OnTjFv3jxnTFKpaiUxMRHAKWsGDoeBiBQCDwJrgB3AhyKS6Oh0S0RHR9OnTx/eeOONMht0lKrptm3bBnhIGACIyOci0l5E2ojIi86Y5tkefvhh9u7dy+rVq509aaW82rZt22jYsGHpXhBHeMW5CQMHDiQsLIyZM2daXYpSHiUxMZFOnTphTHnb8SvHK8LAz8+PBx54gG+++abc/bNK1UQiwrZt25yzJwEvCQOAO++8E4ClS5daXIlSnuHw4cNkZmbWvDCIjIykd+/eLFmyxOpSlPIIztx4CF4UBlB8oMuWLVtKd6coVZOVhEGnTp2cMj2vCoOhQ4fi4+OjawdKUbzxMDQ09JwjLx1hyXgGVdWsWTP69evH4sWLeeGFF/DxqViWZWRksH37dnbt2kVeXh7GGEJDQ7nsssto3749fn5e9c+gFAA7duwoPcTZGbzuWzBq1CjuuOMOPvnkE26//fYLvi4nJ4dFixaxePFivv/++zJnmJWoU6cOvXv3pn///gwePLj0rDilPF1qaipt27Z12vS8qpsAMHjwYNq1a0dcXFy5X/CioiLmzp1Lu3btGD16NEeOHOHpp5/m008/JSkpieTkZA4ePMjPP//MBx98wL333ktKSgqPPPII4eHh3HDDDSxfvpyCggILlk6piktPT3daFwFw/ESlqtwqeqLShcybN08AWb169TnPb9++Xa666ioBpHfv3vLtt9+KzWar0DSTkpIkLi5OWrduLYC0bt1a5syZI7m5uQ7VqpQrFBQUnHPac0VxkROVvDIM8vPzpXXr1hIdHS2HDh2SnJwcee655yQwMFAaNmwoixYtqnAInK+wsFBWrlxZGiqtW7eWhQsXSlFRkUM1K+VMqampAsisWbMq9b5qFwYiIosWLRJjjPj6+krTpk0FkGHDhsnRo0cdnrZI8QAVa9aska5duwogPXv2lPj4eKdMWylH7dy5UwBZtGhRpd53sTDwum0GJe6880727NnD448/TkxMDN988w1Lly49ZwALRxhj+NOf/sTGjRtZsGABBw4coGfPnowbN67c4cCUcqeSgVscHursbBdKCVfenLFm4G4ZGRny8MMPi4+PjzRv3lyWL19e5a6IUo767LPPBJBffvmlUu+jOq4ZuFtISAgzZsxgw4YNtGjRgiFDhjBgwAAOHTpkdWmqBkpPTwdw6t4EDYNKio6O5pdffuG1115j7dq1dOzYkRkzZpQZD08pV3JFN0HDoAr8/Px49NFHSUxM5Nprr2XSpEn06NGDH3/80erSVA1x4sQJfHx8SgewdQYNAwdERkayevVqli9fTnp6Otdeey133HEHe/bssbo0Vc2lp6fToEGDCh+SXxEaBg4yxjB48GB27tzJs88+y6pVq+jQoQNjx47VgViUy5w4ccK5exLQMHCa2rVr8/zzz7Nnzx7uv/9+/vOf/9C2bVvuvvtuNmzYYHV5qprRMPACoaGhzJ49m7179zJhwgQ++eQTrrrqKnr06MHs2bPLvaSZUpXl9PMS0DBwmbCwMF5//XUOHz7Mm2++SV5eHuPHj6d58+b85S9/4YMPPtCDl1SV6ZqBF6pXrx7jx49ny5Yt/Pbbb0yYMIHNmzdzzz330KRJE2666Sbefvttjhw5YnWpyotoGHi5rl278tprr7F//35+/vlnHnnkEfbu3cvYsWNp2bIlvXv3Ztq0aezdu9fqUpUHKywsJDMzU7sJ1YGPjw+9evVi6tSpJCUlkZiYSFxcHAUFBUyePJk2bdoQExPD9OnTdY1BlXHy5EnAyecloGFgOWMMHTt25KmnniIhIYF9+/Yxbdo0jDE89thjtGrViltvvZWVK1fqUY4KcNFJSmgYeJyIiAgee+wx4uPj2blzJ7GxsWzatIkBAwbQtm1bZs6cyenTp60uU1nIFeclgIaBR4uKiiIuLo4DBw7w8ccfExYWxsSJE2nTpg2zZ8/WodlqKI9cMzDGDDHGJBpjbMaYGGcVpc7l5+fHwIED+fHHH1m3bh3t27dn/PjxREdHk5CQYHV5ys08MgyAbcDtwA9OqEVVwDXXXMP333/PihUrSE9P56qrruL555/HZrNZXZpyE4/sJojIDhH53VnFqIoxxjBgwAASExO56667eO655xg0aBDZ2dlWl6bcoOSMxXr16jl1um7bZmCMGW2MSTDGJKSmprprttVa/fr1WbBgATNmzGDVqlXccMMNZGZmWl2WcrETJ044/YxFqEAYGGP+Z4zZVs7ttsrMSETmikiMiMQ0adKk6hWrcxhjePjhh1mxYgWbNm3illtu0TWEas4V5yVABa6oJCJ/dPpcldP99a9/ZdmyZQwdOpRBgwbxxRdf4Ovra3VZygVccSgy6K7FauX222/nnXfeYe3atTz77LNWl6NcxCPDwBgz0BiTDPQGVhtj1jinLFVVo0aN4v777+ell15i5cqVVpejXKBkm4GzObo3YYWIhIlIoIg0E5GbnFWYqrpZs2YRHR3NqFGj0I211U9GRobnhYHyTEFBQbz//vtkZWUxceJEq8tRTmSz2cjMzKR+/fpOn7aGQTVVcvLT4sWLWb16tdXlKCc5deoUNptNw0BVTmxsLJ06dWL8+PHk5uZaXY5ygpLRsTQMVKUEBAQwc+ZMDhw4wJtvvml1OcoJNAxUlfXr14/+/fvz4osvlh7TrryXhoFyyKuvvkpWVhZxcXFWl6IcpGGgHNKpUydGjhzJ7NmzOXz4sNXlKAdoGCiHPfXUU9hsNqZOnWp1KcoBGgbKYREREYwYMYK5c+eSkpJidTmqikrCwJkXXC2hYVCDxMbGUlhYyLRp06wuRVVRRkYGwcHB+Pld8hzDStMwqEHatGnD3/72N95++23S0tKsLkdVQUZGhku6CKBhUOM8/vjjnDlzhjlz5lhdiqoCDQPlNJ06daJ///688cYbelSiF9IwUE71+OOPk5qayvvvv291KaqSNAyUU/3hD38gOjqa6dOn66jKXkbDQDmVMYZHH32UpKQk1qzR8Wi8iYaBcrpBgwbRvHlz3njjDatLURXkyrEMQMOgxgoICOCBBx7giy++YNeuXVaXoyrAlWMZgIZBjTZmzBj8/f1denpzYWEh69atY9myZSxevJgvvviCpKQkvU5kFbjyUGSowFDpqvpq3rw5Q4cOZd68ecTFxVG3bl2nTfv06dPExsayaNGi0msDni0wMJAuXboQExNDTEwM0dHRdOjQgYCAAKfVUN1oGCiXGj9+PIsWLWLJkiWMHj3aKdNMSkpi0KBBJCYmMnz4cAYOHEjHjh3x9fUlLS2N3bt3s3XrVjZu3MgHH3zA7NmzAfD396dDhw506NCBqKgoIiMjiYyMJDw8nFatWtX4oNAwUC7Vq1cvunTpwpw5c7j//vsxxjg0vYMHD9K7d2+MMXz55Zf86U9/Oqc9KiqKPn36lD622WwkJSXx22+/sXnzZhITE/ntt9/46KOPztntaYwhPDycqKgoevXqxXXXXUffvn1r1IViNAyUSxljeOCBBxg7diy//PILvXr1qvK0CgoKGD58OAUFBSQkJNC+fftLvsfHx6d0bWD48OHnTOvQoUPs27ePQ4cOsX//fnbv3s22bduIi4vDZrMRHh7OAw88wKRJkwgKCqpy3d7C1WGAiLj9Fh0dLcpzZGVlSXBwsNxzzz0OTSc2NlYAWbJkiZMqK19WVpYsX75c+vXrJ4B07txZtmzZ4tJ5eoKZM2cKIGlpaVWeBpAgF/heahgoEREZO3asBAYGVvk/2pYtW8QYI6NGjXJyZRf3+eefS7NmzSQwMFDWrFnj1nm72/PPPy+AFBQUVHkaFwsD3bWoABg7dix5eXnMnz+/Su9/8sknCQkJ4dVXX3VuYZdwyy23sGXLFjp06MDtt99OQkKCW+fvTq4cywD0OANld8UVV3DNNdcwZ86cSp+v8OOPP/LZZ5/xxBNPuOSCoJfStGlTvvjiC5o0aUL//v05cOCA22twB1ceigyOX3h1mjFmpzFmizFmhTHGdZUqlxs3bhx79uxh7dq1FX6PiDBlyhRCQ0OZMGGCC6u7uNDQUNasWcOZM2cYM2ZMcR+4mvHoMADWAp1F5EogCYh1vCRlldtvv52mTZtWauCTtWvXsn79ep555hlq167twuourX379rz88susWbOGDz74wNJaXMGjw0BEvhKRQvvD/weEOV6SskpgYCD33Xcfn376Kfv27avQe+Li4ggLC2PkyJEurq5ixo0bR58+fZg4cSLHjh2zuhyn8ugwOM9I4AsnTk9ZYNy4cfj6+vLaa69d8rU//PAD69atY/LkyQQGBrqhukvz8fHh3XffJTs7mxdeeMHqcpzK8jAwxvzPGLOtnNttZ73mKaAQWHSR6Yw2xiQYYxJSU1OdU71yupYtW3LPPffw3nvvXfKXNS4ujqZNm3Lfffe5qbqK6dChA6NGjWLu3Lns3bvX6nKcJj093bUbaC+0z7GiN2AE8DNQu6Lv0eMMPNvOnTvFGCOxsbEXfM13330ngEydOtWNlVVccnKyBAUFyd133211KU6Rm5srgPzzn/90aDq46jgDY8zNwBPAX0Ukx8FcUh4iKiqKQYMG8dZbb5V7xmFhYSETJkygdevWPPTQQxZUeGktW7bkoYceYuHChSQmJlpdjsNKLprbuHFjl83D0W0GbwJ1gbXGmE3GmLedUJPyAM888ww5OTmMGzeuzG66f//732zZsoXXXnuNWrVqWVThpT3xxBPUqVOHl156yepSHFbStW7SpInrZnKhVQZX3rSb4B3i4uIEkIULF5Y+t2vXLmnUqJH07dtXbDabhdVVzKOPPio+Pj6yZ88eq0txyNq1awWQ77//3qHpoIcjq6qYMmUKffr0Ydy4cbz66qt8+OGH9OjRA4C33nrL4dOd3eGRRx7Bz8/P6y8pV3IFLE/uJqhqzNfXl4ULF9K5c2eeeOIJhg0bRlhYGPHx8XTs2NHq8iqkRYsW3HvvvcybN4+jR49aXU6VuaOboGGgLioiIoL169eze/duFi9ezE8//URkZKTVZVXK5MmTKSgoYMaMGVaXUmWpqakYY1y6a1HDQFVImzZtGD58uFPHSXSXNm3aMHjwYN5++22ysrKsLqdK0tLSaNiwoUtHdtIwUDXC5MmTyczMZO7cuVaXUiWpqamu3ZOAhoGqIaKjo+nXrx//+te/yMvLs7qcStMwUMqJJk+ezJEjR1i06IJHzXustLQ0l+5JAA0DVYPceOONdOvWjalTp1JUVGR1OZWiawZKOZExhilTppCUlMSKFSusLqfCbDYb6enpGgZKOdOgQYNo27YtL7/8steMhnTy5EmKioo0DJRyJl9fXyZPnsyvv/5aqeHdrOSOow9Bw0DVQPfccw8tWrTwmhOY3HKSEhoGqgYKDAzkscce4/vvv2f9+vVWl3NJGgZKudDo0aNp3LgxL774otWlXJJ2E5RyoTp16jBp0iS++OILNm7caHU5F6VrBkq52Pjx4wkJCfH4bQepqakEBwe7/OKyGgaqxgoJCeHBBx/k448/Zvv27VaXc0HuOPoQNAxUDTdx4kRq167Nyy+/bHUpF+SOow9Bw0DVcI0bN2bMmDEsWbLEY4dV1zBQyk0ee+wxfH19mTp1qtWllCs1NVW7CUq5Q4sWLRg5ciTz58/n8OHDVpdzjvz8fI4ePUqLFi1cPi8NA6UoPr25qKioQpeVc6ctW7ZQUFBA9+7dXT4vDQOlgMjISP72t7/xzjvvcPz4cavLKbVhwwYAevbs6fJ5aRgoZRcbG0tubi6vv/661aWUio+Pp0mTJoSHh7t8XhoGStlFRUUxbNgw3njjDY9ZO4iPj6dHjx5uuUaFhoFSZ3nuuefIzc3llVdesboUsrOz2b59u1u6CKBhoNQ5oqKiGDFiBLNnz3b6noXKDqby66+/IiKlV7FyNQ0Dpc7z7LPPYrPZeP7556s8DZvNxvr165kyZQq9e/emefPm+Pj44O/vT/369enevTt33XUXc+bMYdeuXeVOo2TjobvCwM+RNxtj/gncBtiA48C9InLEGYUpZZWIiAgefPBBZsyYwb333svVV19d4ffabDb++9//EhcXx9atW/Hz86NXr178+c9/pkWLFhQVFZGVlcXu3bv57rvvWLx4MQCdO3fmzjvvZMSIEaXHFMTHxxMREeGWow8Bx67CDNQ76/4E4O2KvE+vwqw8XXZ2toSHh0vHjh0lNze3Qu9JTEyU3r17CyBRUVHy3nvvycmTJy/4epvNJklJSTJr1iy5+uqrBRA/Pz8ZPny4xMbGSuPGjWXIkCHOWiQRufhVmJ12mXUgFphTkddqGChvsHr1agHk6aefvujrCgoK5MUXXxR/f39p1KiRzJ8/XwoLCys9v127dsnEiROlXr164ufnJw0aNJClS5dWtfxyuTQMgBeBQ8A2oElF3qNhoLzFiBEjBJBZs2aV275z587SX/Vhw4bJ8ePHHZ6nzWZzeBoX4lAYAP+zf9HPv9123utigecvMp3RQAKQEB4e7rKFVcqZ8vLyZMCAAaVrCEeOHBGbzSZ79uyRiRMnip+fn4SEhMjixYutLrVC3NVNaA1sq8hrdc1AeZO8vDwZOnSoAOLj4yMNGzYsvT969Gg5evSo1SVW2MXCwNG9Ce1EpGS/yF+BnY5MTylPFBAQwLJly3j++edZvHgxKSkpxMTEcP3119O+fXury3MaUxwWVXyzMR8BURTvWjwAPCAilzxSIyYmRhISEqo8X6VU1RhjNopITHltDq0ZiMggR96vlPIcegSiUgrQMFBK2WkYKKUABzcgVnmmxqRSvMHxUhoDaS4ux110WTxPdVkOqPiytBaRck92sCQMKsoYk3ChLZ/eRpfF81SX5QDnLIt2E5RSgIaBUsrO08NgrtUFOJEui+epLssBTlgWj95moJRyH09fM1BKuYnHhoEx5mZjzO/GmN3GmClW11NVxpj9xpitxphNxhivOiHDGPOeMea4MWbbWc81NMasNcbssv9tYGWNFXWBZXnOGHPY/tlsMsb0t7LGijDGtDLGfGuM2WGMSTTGPGx/3uHPxSPDwBjjC7wF3AJ0BIYbYzpaW5VDrheRrl64G2s+cPN5z00BvhaRdsDX9sfeYD5llwXgX/bPpquIfO7mmqqiEHhURC4HegHj7d8Nhz8XjwwDoCewW0T2ikg+sJTigVeVG4nID8CJ856+DVhgv78AGODWoqroAsvidUQkRUR+td/PBnYALXHC5+KpYdCS4qHUSiTbn/NGAnxljNlojBltdTFO0ExEUqD4PybQ1OJ6HPWgMWaLvRvhFV2eEsaYCKAb8AtO+Fw8NQzKu5aUt+726CMi3Snu8ow3xlxndUGq1BygDdAVSAGmW1tOxRljgoGPgIkikuWMaXpqGCQDrc56HAZ45fUYxH4dCRE5DqyguAvkzY4ZY0IB7H8946KEVSAix0SkSERswL/xks/GGONPcRAsEpGP7U87/Ll4ahjEA+2MMZHGmADgDmCVxTVVmjGmjjGmbsl94E8UDybrzVYBI+z3RwArLazFISVfHruBeMFnY4qvwPofYIeInH25aIc/F4896Mi+m2cG4Au8JyIvWlxSpRljLqN4bQCKR5Va7E3LYYxZAvSl+Iy4Y8A/gE+AD4Fw4CAwREQ8fsPcBZalL8VdBAH2A2NK+t2eyhhzDbAO2ErxcIMAT1K83cChz8Vjw0Ap5V6e2k1QSrmZhoFSCtAwUErZaRgopQANA6WUnYaBUgrQMFBK2WkYKKUAB6+1WFWNGzeWiIgIK2atVI22cePGtAtdN8GSMIiIiECvwqyU+xljLnjxIu0mKKUADQOllJ2GgVIKsGibQXkKCgpITk4mNze3TFtQUBBhYWH4+/tbUJlSnuG734/z3e+pbEnO4MTpfAqK/u+M466t6vPWXd0dmr7HhEFycjJ169YlIiKC4vEbiokI6enpJCcnExkZaWGFSlkjIyefZ1cmsmrzEWr5+3JFyxC6tKqPr4/Bx/5duaxJHYfn4zFhkJubWyYIAIwxNGrUiNTUVIsqU8o6aafyGPDWeo5m5vLIje0Z27cN/r6u6d1XKgyMMfWBd4HOFI8OMxL4HVgGRFA8WsxQETlZlWLOD4JLPa9UdVZQZGPcol9Jzc5j2ZjeRLd27eDNlY2YmcCXItIB6ELxmO3eelENpTzaPz/bzoZ9J3h18JUuDwKoRBgYY+oB11E8GCMiki8iGXjpRTWU8mTf7DzG+z8f4L5rIrmtq3suGVKZNYPLgFRgnjHmN2PMu/YRf512UY0Ljceo4zSqmiQzp4ApH20lqlldHr85ym3zrUwY+AHdgTki0g04TSW6BMaY0caYBGNMQnkbA4OCgkhPTy/zxS/ZmxAUFFSJUpXyXi98tp300/m8NqQLgX6+bptvZTYgJgPJIvKL/fF/KQ6DY8aYUBFJudjFG0RkLjAXICYmpsxPfVhYGMnJyeXuNSg5zkCp6u6bncf46NdkHry+LVeEhbh13hUOAxE5aow5ZIyJEpHfgX7AdvttBPAKDlxUw9/fX48jUDVaZk4BsR8Xdw8e6tfW7fOv7HEGDwGL7Fc52gv8neKuxofGmFHYL97g3BKVqhle+Gw7aafyefeeHm7tHpSoVBiIyCYgppymfs4pR6ma6esd1nUPSuiJSkpZ7OTpfKZ8vJUOza3pHpTwmMORlaqp/rEqkZOn85n/d2u6ByV0zUApC32+NYVVm48woV87OrWwpntQQsNAKYscz87lqRVb6RIWwti+bawuR8NAKSuICLEfbSUnv4jpQ7u67EzEyrC+AqVqoGXxh/h653Em39yBtk2DrS4H0DBQyu0OpJ/mhc+2c3WbRvz96giryymlYaCUGxXZhEc/3Iyvj2HakC74+HjOWB26a1EpN3rnhz0kHDjJv4Z1oWX9WlaXcw5dM1DKTbYdzuRfa5O49YpQBrhpjILK0DBQyg1yC4qYtGwTDWoH8OLAzh45lJ92E5Ryg6lf7mTX8VMsGNmT+rUDrC6nXLpmoJSLrduVyrz1+7n36gj+0L7ca556BA0DpVwoIyefx5Zvpm3TYKbc0sHqci5Kw0ApF3pmZSLpp/KZMawrQf7WnYRUERoGSrnIyk2H+XTzESbd2J7OLa09CakiNAyUcoGUzDM888k2ols3YMx1l1ldToVoGCjlZDab8PjyLRTahNeHdsHPA05CqgjvqFIpL7LwlwP8uDuNp269nNaNHL8gqrtoGCjlRHtTT/HS5zv4Q/sm3Nkz3OpyKkXDQCknKbIJjy3fTKCfL68OvtIjjzK8GD0CUSkneXfdXn49mMHMO1krH3cAAA79SURBVLrSrJ73XQFM1wyUcoJdx7KZvjaJmzo1469dWlhdTpVoGCjloMIiG4/9dwt1AnyJG3CF13UPSmg3QSkHvbd+H5sPFXcPmtQNtLqcKtM1A6UcsDf1FNO/SuLGjt7bPSihYaBUFdlswpSPthLo58OLAzxzjILKqHQYGGN8jTG/GWM+sz+ONMb8YozZZYxZZr8oq1LV3pL4g2zYf4Knb+1IUy/ce3C+qqwZPAzsOOvxVOBfItIOOAmMckZhSnmyo5m5vPL5Tq5u04ghMWFWl+MUlQoDY0wYcCvwrv2xAW4A/mt/yQJggDMLVMoT/WPVNvKLbLx8u/fuPThfZdcMZgCTAZv9cSMgQ0QK7Y+TAc8b6VEpJ/oq8ShrEo8x8Y/tvercg0upcBgYY/4MHBeRjWc/Xc5L5QLvH22MSTDGJKSmplayTKU8w6m8Qv6xKpEOzety37WRVpfjVJVZM+gD/NUYsx9YSnH3YAZQ3xhTcrxCGHCkvDeLyFwRiRGRmCZNPHccOKUuZvpXv3M0K5eXbr/CI66P6EwVXhoRiRWRMBGJAO4AvhGRu4BvgcH2l40AVjq9SqU8wNbkTBb8tJ+7rgqne3gDq8txOmdE2xPAI8aY3RRvQ/iPE6aplEcpLLIRu2ILjYIDmXyzZw9sWlVVOhxZRL4DvrPf3wv0dF5JSnme938+wLbDWbwxvBv1gvytLsclqlenRykXOJJxhulf/c4f2jfhz1eGWl2Oy2gYKHURIsKzK7dhE4irBoccX4yGgVIX8eW2o/xvx3Em3diOVg1rW12OS2kYKHUBmTkFPLsqkY6h9RjZp3odU1AeHc9AqQv45+rtnDidz7x7e3jNcOeOqP5LqFQVfPf7cf67MZkH/nCZV1wNyRk0DJQ6T2ZOAbEfb6Vt02AeuqGd1eW4jYaBUud5dtU2UrPzeH1oF4+/WKozaRgodZZVm4+wctMRJvRrx5Vh9a0ux600DJSyO3Qih6dWbKVrq/qM69vG6nLcTsNAKSC/0MaDi38F4I3h3WrE3oPz6a5FpYCXPt/B5uRM3v5bdLU/uOhCal78KXWejzYmM/+n/YzsE8nNnZtbXY5lNAxUjfbbwZPErtjK1W0aEdu/ep6aXFEaBqrGOnQih9EfbKRZvUDeurN7tRu5qLJ0m4GqkU6ezmfEvA3kFRSx6L6raFBHL/ehYaBqnNN5hYxaEE/yyTMsHHUV7ZvVtbokj1Cz14tUjZOTX8jI+fFsTs5k1h1d6RnZ0OqSPIaGgaoxcvILuW9BAvH7T/D60C7c3Ln6jlpUFdpNUDVCRk4+f58fz+ZDGbw2pAu3ddVr/ZxPw0BVe8kncxg5P579aTnMviu6Rh9LcDEaBqpaS9h/gjEfbCS/yMb8kT24uk1jq0vyWBoGqloSEeb/tJ+XPt9BWIPavDsihjZNgq0uy6NpGKhq58TpfKZ8tIWvth/jj5c35bUhXahfW48juBQNA1WtfLntKE9/spXMMwU8fevljLomsloPb+5MGgaqWjiccYbnVyXy1fZjdG5Zj4X3XUWH5vWsLsurVDgMjDGtgPeB5oANmCsiM40xDYFlQASwHxgqIiedX6pSZeUWFPHuur289e0eAJ64uQP3XRtZ488zqIrKrBkUAo+KyK/GmLrARmPMWuBe4GsRecUYMwWYQvHFWJVyGZtN+HTLEV798ncOZ5zhpk7NeObPHQlrUDPHInCGCoeBiKQAKfb72caYHUBL4Dagr/1lCyi+IKuGgXKZn3an8cqXO9mSnEnH0HpMG3Kl7jJ0giptMzDGRADdgF+AZvagQERSjDFNnVadUmfZdjiTqV/uZN2uNFqEBDF9SBcGdmuJj49uIHSGSoeBMSYY+AiYKCJZFd1Sa4wZDYwGCA8Pr+xsVQ22J/UUr69NYvWWFOrX9uep/pdzd+/WNWoYc3eoVBgYY/wpDoJFIvKx/eljxphQ+1pBKHC8vPeKyFxgLkBMTIw4ULOqIZJP5jDzf7v46Ndkgvx9mXBDW+677jLqBflbXVq1VJm9CQb4D7BDRF4/q2kVMAJ4xf53pVMrVDXO8axc3vx2N0s2HMQYw71XRzLu+jY0Dg60urRqrTJrBn2Au4GtxphN9ueepDgEPjTGjAIOAkOcW6KqKdJP5fHOD3tZ8NN+Cm3C0JgwHrqhHS3q17K6tBqhMnsTfgQutIGgn3PKUTVRRk4+/163l3nr95NbUMSAri2Z0K8dEY3rWF1ajaJHICrLZOYU8O6PxSFwOr+QW68IZeIf29G2qQ5DZgUNA+V2GTn5vPfjPuat3092XiH9r2jOw/3aE9VcQ8BKGgbKbdJP5fGfH/fx/s8HOJVXyC2dmzOhXzsuD9VzCDyBhoFyueNZucz9YS+LfjlIbmER/a8I5aEb2uqJRB5Gw0C5zOGMM7zz/R6Wxh+iyCbc1qUF465vo9sEPJSGgXK65JM5vPXtHv678RAAg7qHMbZvG1o30r0DnkzDQDnN8axc3vhmN0vjD2IwDOvRirF929JSjxPwChoGymE5+YXM/WEvb3+/h8IiYWiPVjx4fVs9WMjLaBgoh6zdfoxnV24jJTOXW68MZfJNUdod8FIaBqpK0k/l8czKbXy+9Sgdmtdl1vBu9IjQS5V5Mw0DVWk/JKXy6PLNZOYU8PhNUYy+7jIdZqwa0DBQFVZkE2b+L4lZ3+ymfbNg3h/ZUw8YqkY0DFSFZOTkM2HpJn5ISmVoTBgv3NZZBxepZjQM1CUlHcvm/vcTSMnI5eXbr2B4Tx2pqjrSMFAX9VXiUR75cDO1AnxZMroX0a0bWF2SchENA1Uum02Y+fUuZn69iy5hIbxzdwzNQ4KsLku5kIaBKiPtVB6Tlm1i3a40BkeHETdAtw/UBBoG6hzfJ6Xy+PLNZJwp4MWBnbmzZ7heq7CG0DBQQPHeghdX72D5xmTaNg1mge42rHE0DGq4wiIbizcc5PW1SWTnFjKubxsm9Gun3YIaSMOghsovtPHxr8nM/m4PB0/kcHWbRjzz5466NlCDaRjUMIdO5LAs/hBL4w+SdiqfK8NC+MdfYrihQ1PdNlDDaRjUAIdO5PC/HcdYvSWFhAMnMQZuiGrKiKsjuLZdYw0BBWgYVDsiwuGMM/x6MIMN+9L5aU86e1NPA9C+WTCP3xTFX7u0oFVDvXS5OpeGgRfLLShiX9ppko5ls/NoNtuPZJF4JIu0U3kA1AnwpWdkQ+7sGc4fL2+mFyVRF6Vh4OFy8gtJPnmGQydyOJCew4H00+xLz2F/2mkOncxB7Jew9fMxtG0aTN+oJlwZFkL38AZ0aF4XPz21WFWQU8LAGHMzMBPwBd4VkVecMd3qrqDIxvHsPI5m5nIsK5cjGWdIySz+ezjjDIdPniH9dP457wkO9COicW2uDAthYLeWtG0aTLtmwVzWOJgAP/3iq6pzOAyMMb7AW8CNQDIQb4xZJSLbHZ22N8ovtJGRk0/66XxOnC7+m5adR9qpPFKz8zhuv6Vm55J+Or/0l71ELX9fWjaoRWhIEJ1ahBDWoBZhDWrRqmFtWjesTcM6AbrBT7mEM9YMegK7RWQvgDFmKXAb4FFhUGQTCopsFBTZKCyy37cJhfbn8gpt5NtvJfdzC4vILbCRW1DEmfwicvKLyMkv5HR+IadyCzmVV0hWbiFZZwrIOlNAxpkCcvKLyp2/n4+hcXAgjesG0CIkiK6tQmhaN4hm9YIIDQmieUgQLUJqUa+Wn37ZlSWcEQYtgUNnPU4GrnJ0ossTDpF4JKv0sU2EIlvxLb/IRl7Jl9R+yyso/vKWfKELiux/7SFw/i9wVQX4+RAc6EedQF/qBvoTHORHWIPa1GvhR4PaAdSv5U+DOgE0qB1Ao+AAGtUJoHFwICG1/PHx0S+58lzOCIPy/oeX+eoZY0YDowHCwy89OEb8/hN8ue1o6WMfH4OPMfj5GPx9fagV4EuQvw+1/H0JDvSjUR1fAv19CPT1IcCv+ObnU/zX39fg5+ODn68hwLf4r79v8fP+vj74+foQ4OtDoN//vTfIz5cAv+LpB/kXz6+Wv69ukFPVlhEHfzKNMb2B50TkJvvjWAAReflC74mJiZGEhASH5quUqjxjzEYRiSmvzRk/c/FAO2NMpDEmALgDWOWE6Sql3MjhboKIFBpjHgTWULxr8T0RSXS4MqWUWznlOAMR+Rz43BnTUkpZQ7eGKaUAJ2xArNJMjUkFDlTgpY2BNBeX4y66LJ6nuiwHVHxZWotIk/IaLAmDijLGJFxoy6e30WXxPNVlOcA5y6LdBKUUoGGglLLz9DCYa3UBTqTL4nmqy3KAE5bFo7cZKKXcx9PXDJRSbuKxYWCMudkY87sxZrcxZorV9VSVMWa/MWarMWaTMcarTsgwxrxnjDlujNl21nMNjTFrjTG77H+94kqsF1iW54wxh+2fzSZjTH8ra6wIY0wrY8y3xpgdxphEY8zD9ucd/lw8MgzOGjDlFqAjMNwY09HaqhxyvYh09cLdWPOBm897bgrwtYi0A762P/YG8ym7LAD/sn82Xe1H0nq6QuBREbkc6AWMt383HP5cPDIMOGvAFBHJB0oGTFFuJCI/ACfOe/o2YIH9/gJggFuLqqILLIvXEZEUEfnVfj8b2EHxmCIOfy6eGgblDZjS0qJaHCXAV8aYjfYxHbxdMxFJgeL/mEBTi+tx1IPGmC32boRXdHlKGGMigG7ALzjhc/HUMKjQgCleoo+IdKe4yzPeGHOd1QWpUnOANkBXIAWYbm05FWeMCQY+AiaKSNalXl8RnhoGyUCrsx6HAUcsqsUhInLE/vc4sILiLpA3O2aMCQWw/z1ucT1VJiLHRKRIRGzAv/GSz8YY409xECwSkY/tTzv8uXhqGFSLAVOMMXWMMXVL7gN/ArZd/F0ebxUwwn5/BLDSwlocUvLlsRuIF3w2pni03P8AO0Tk9bOaHP5cPPagI/tunhn834ApL1pcUqUZYy6jeG0AiseOWOxNy2GMWQL0pfiMuGPAP4BPgA+BcOAgMEREPH7D3AWWpS/FXQQB9gNjSvrdnsoYcw2wDtgK2OxPP0nxdgOHPhePDQOllHt5ajdBKeVmGgZKKUDDQCllp2GglAI0DJRSdhoGSilAw0ApZadhoJQC4P8DoqrXVz3+BTsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "omega=1\n",
    "ampl0=0.00001\n",
    "#T=8\n",
    "mu=0.5\n",
    "g=0.2*omega\n",
    "#ampl0=1\n",
    "T=20\n",
    "nt=100\n",
    "tlist = np.linspace(0,T, nt)\n",
    "def hamiltonian(omega, ampl0, mu,g):\n",
    "      \"\"\"Two-level-system Hamiltonian\n",
    "\n",
    "      Args:\n",
    "          omega (float): energy separation of the qubit levels\n",
    "          ampl0 (float): constant amplitude of the driving field\n",
    "      \"\"\"\n",
    "      HA = -1j*qutip.Qobj(H0(g))\n",
    "      H1 =-1j*qutip.Qobj(Hde(mu,g))\n",
    "      H2= -1j*qutip.Qobj(He(mu))\n",
    "      def guess_control(t, args):\n",
    "          return ampl0*krotov.shapes.flattop(\n",
    "              t, t_start=0, t_stop=T, t_rise=0.005, func=\"blackman\"\n",
    "          )\n",
    "      def guess_control2(t, args):\n",
    "          return ampl0*krotov.shapes.flattop(\n",
    "              t, t_start=0, t_stop=T, t_rise=0.005, func=\"blackman\"\n",
    "          )\n",
    "      \n",
    "      return [HA, [H1, guess_control], [H2, guess_control]]\n",
    "def S(t):\n",
    "    \"\"\"Shape function for the field update\"\"\"\n",
    "    return krotov.shapes.flattop(\n",
    "        t, t_start=0, t_stop=T, t_rise=0.05 * T, t_fall=0.05 * T, func='sinsq'\n",
    "    )\n",
    "def plot_iterations(opt_result):\n",
    "    \"\"\"Plot the control fields in population dynamics over all iterations.\n",
    "\n",
    "    This depends on ``store_all_pulses=True`` in the call to\n",
    "    `optimize_pulses`.\n",
    "    \"\"\"\n",
    "\n",
    "    fig, [ax_ctr,ax] = plt.subplots(nrows=2, figsize=(4, 5))\n",
    "    n_iters = len(opt_result.iters)\n",
    "    EEnergy=np.zeros(nt)\n",
    "    for (iteration, pulses) in zip(opt_result.iters, opt_result.all_pulses):\n",
    "        controls = [\n",
    "            krotov.conversions.pulse_onto_tlist(pulse)\n",
    "            for pulse in pulses\n",
    "        ]\n",
    "        objectives = opt_result.objectives_with_controls(controls)\n",
    "        dynamics = objectives[0].mesolve(\n",
    "            opt_result.tlist, e_ops=[]\n",
    "        )\n",
    "        if iteration == 0:\n",
    "            ls = '--'  # dashed\n",
    "            alpha = 1  # full opacity\n",
    "            ctr_label = 'guess'\n",
    "            pop_labels = ['0 (guess)', '1 (guess)']\n",
    "        elif iteration == opt_result.iters[-1]:\n",
    "            ls = '-'  # solid\n",
    "            alpha = 1  # full opacity\n",
    "            ctr_label = 'optimized'\n",
    "            pop_labels = ['0 (optimized)', '1 (optimized)']\n",
    "        else:\n",
    "            ls = '-'  # solid\n",
    "            alpha = 0.5 * float(iteration) / float(n_iters)  # max 50%\n",
    "            ctr_label = None\n",
    "            pop_labels = [None, None]\n",
    "        ax_ctr.plot(\n",
    "            dynamics.times,\n",
    "            controls[0],\n",
    "            label=ctr_label,\n",
    "            color='black',\n",
    "            ls=ls,\n",
    "            alpha=alpha,\n",
    "        )\n",
    "    EField=np.transpose(np.array(opt_result.optimized_controls))\n",
    "    EEnergy[0]=(np.square(EField[0]))*(T/nt)\n",
    "    a=0\n",
    "    for i in range (1,nt):\n",
    "      a+=np.square(EField[i-1])\n",
    "      EEnergy[i]=(np.square(EField[i])+a)*(T/nt)\n",
    "      \n",
    "    \n",
    "    ax.plot(tlist,np.transpose(EEnergy))\n",
    "    plt.legend()\n",
    "    plt.show(fig)\n",
    "    \n",
    "    \n",
    "    \n",
    "\n",
    "\n",
    "\n",
    "H = hamiltonian(omega,ampl0,mu,g)\n",
    "pulse_options = {\n",
    "    H[1][1]: dict(lambda_a=0.1, update_shape=S), \n",
    "\n",
    "    \n",
    "}\n",
    "target=np.array([1,0, 2*np.sqrt(2),0,2*np.sqrt(2),1/2,0,0,0,1/2,0,0,1/2,0,1/2])\n",
    "#target=np.array([1,np.sqrt(2), 2*np.sqrt(2),np.sqrt(2),2*np.sqrt(2),1/2,0,0,0,1/2,0,0,1/2,0,1/2])\n",
    "r=np.pi/4\n",
    "r2=np.pi/8\n",
    "theta=np.pi/4\n",
    "objectives = [\n",
    "    krotov.Objective(\n",
    "        #initial_state=qutip.Qobj(np.array([1,0,0,0,0,1/2*np.cosh(2*r2),0,-1/2*np.sinh(2*r2)*np.cos(theta),-1/2*np.sinh(2*r2)*np.sin(theta),1/2*np.cosh(2*r2),-1/2*np.sinh(2*r2)*np.sin(theta),1/2*np.sinh(2*r2)*np.cos(theta),1/2*np.cosh(2*r2),0,1/2*np.cosh(2*r2)])), target=qutip.Qobj(qutip.Qobj(np.array([1,0,0,0,0,1/2*np.cosh(2*r),0,-1/2*np.sinh(2*r)*np.cos(theta),-1/2*np.sinh(2*r)*np.sin(theta),1/2*np.cosh(2*r),-1/2*np.sinh(2*r)*np.sin(theta),1/2*np.sinh(2*r)*np.cos(theta),1/2*np.cosh(2*r),0,1/2*np.cosh(2*r)]))), H=H\n",
    "        initial_state=qutip.Qobj(np.array([1,0,0,0,0,1/2,0,0,0,1/2,0,0,1/2,0,1/2])), target=qutip.Qobj(target), H=H\n",
    "      )\n",
    "  ]\n",
    "\n",
    "opt_result = krotov.optimize_pulses(\n",
    "  objectives,\n",
    "  pulse_options=pulse_options,\n",
    "  tlist=tlist,\n",
    "  fieldcoupling=mu,\n",
    "  propagator=krotov.propagators.expm,\n",
    "  chi_constructor=krotov.functionals.chis_ss,\n",
    "  info_hook=krotov.info_hooks.print_table(J_T=krotov.functionals.J_T_ss),\n",
    "  check_convergence=krotov.convergence.Or(\n",
    "      krotov.convergence.value_below('9e-3', name='J_T'),\n",
    "      krotov.convergence.check_monotonic_error,\n",
    "      #krotov.convergence.check_monotonic_fidelity,\n",
    "\n",
    "  ),\n",
    "    store_all_pulses=True,\n",
    "    overlap=inte.coherent_overlap,norm='gaussiano'\n",
    ")\n",
    "plot_iterations(opt_result)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "target=np.array([1,0.6*np.sqrt(2), 2*np.sqrt(2),-0.6*np.sqrt(2),np.sqrt(2),1/2,0,0,0,1/2,0,0,1/2,0,1/2]) for T=30 good\n",
    "target=np.array([1,0.6*np.sqrt(2), 2*np.sqrt(2),-0.6*np.sqrt(2),np.sqrt(2),1/2,0,0,0,1/2,0,0,1/2,0,1/2]) for T=30 good\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot in phase space the last state and the target for the battery"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAI/CAYAAABEVcwAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dfbRlZ10n+O8vdasqJCBIUiLmxYoSGQtbG6yOOCDSRDBBh9IxWV1xjUQWmlEJ4OisnuCswe5ov+DqRcaXtHQkaaMCCSsgFhqNQBx1bIy5gQgUMU2ZhqQ6UYokEyAJqVTqmT/OiVwvt1I3dc/e+7x8PmvdVeec/Zyzf3ffc6u+9fyevU+11gIAwGQdN3QBAADzSMgCAOiAkAUA0AEhCwCgA0IWAEAHhCwAgA4sDV3AaieffHLbvn370GUAABvwqU8ln//80FX04ZbPtda2rbVl6kLW9u3bs7y8PHQZAMAx+vCHk+/5nqGr6Et95khbtAsBgIlpLfmpn0oeemjoSoYnZAEAE/P+949ahQhZAMCEHD6cvP71yYMPDl3JdBCyAICJ+MQnknvvHbqK6SFkAQAT8cEPJo89NnQV00PIAgAm4n3vS770paGrmB5CFgCwYY8+mtx889BVTJdjDllVdVpV/UlV3VZVe6vqjWuMqar6laraV1Ufq6oXbKxcAGAa/dVfJVu2DF3FdNnIxUgPJfnZ1tpHquppSW6pqg+01j65Ysy5Sc4cf31Hkl8f/wkAzJHbb7cea7Vjnslqrd3TWvvI+PYXktyW5JRVw3Yl+a028pdJnlFVzz7magGAqfTgg8mhQ0NXMV0msiarqrYneX6Sm1ZtOiXJXSvu789XBjEAYMY99JCZrNU2HLKq6qlJ3pPkp1trqz8KstZ4SlvjNS6qquWqWj5w4MBGSwIAerZ5c3Kc0+n+kQ0djqranFHAekdr7b1rDNmf5LQV909NcvfqQa21K1prO1trO7dtW/ODrAGAKXbiicmmTUNXMV02cnZhJbkyyW2ttbceYdieJK8en2X4wiQPtNbuOdZ9AgDT6aSTRrNZfNlGzi58UZIfSfLxqrp1/NjPJTk9SVprb0tyfZJXJtmX5KEkr9nA/gCAKfXiFyePPDJ0FdPlmENWa+3/zdprrlaOaUled6z7AABmw9d+bfKsZyV33jl0JdPDEjUAYCLOPTepJ5x+WSxCFgAwEd/3fcnTnjZ0FdNDyAIAJuK7v9u6rJWELABgIr7qq5LXvjbZunXoSqaDkAUATMyllyZLG7l2wRwRsgCAiTnppORNb0pOOGHoSoYnZAEAE/UzP5Mcf/zQVQxPyAIAJuopT0kuu8xslpAFAEzcq1+dXHjhYgctIQsA6MSv/mrykpcsbutQyAIAOrFpU/Le9ybf9E2L+eHRQhYA0JmnPCX50IeSr/maUehaJEIWANCpk09OPvzh5Bu/cbFah0IWANC5005LPvrR5BWvSE48cehqNm7z5uSZz3ziMUIWANCLE05I3ve+5N/9u9HtqqErOjYnnpi86EXJbbc98TghCwDoTVXy+tcnf/VXo/bhLM1qbd48WmP2b/9tcuONo3VmT0TIAgB697znjWaCfvVXR2FlmsPW5s2jtWSvfnXyqU8lb3jD+mbhhCwAYBBLS8lrXpPs3z+6Qvy2bdMVtpaWRuHqgguS229P3v725JRT1v98IQsAGNTmzcmP//gobP3SL43ORnzqU0chZwhPe1qyZUty/vnJJz+ZXH11cvrpT/51hCwAYCps2ZL81E8lf//3yQc/mLzxjaOZoxNOSLZu7W6/VclXfdVo1ur7vi+54orkwIHkne9MzjhjA6/bWptclROwc+fOtry8PHQZAMCUuO225D3vSX7nd5K//dtRGDruuOThh5NHH31yr1U1akk+/vwTTxwFqx/+4eTss598mKuqW1prO9fcJmQBALPi0UeTO+4YrZG6/fbRtbf27k0+/elRaHrsseTw4VGYam3UclxaSp71rOS5z02e//xkx47R7ec+N3nGMzZWzxOFrIG6nQAAT97mzV8OSEfyeNBaWhr2WlxCFgAwVzZtmo7PSbTwHQCgA0IWAEAHhCwAgA4IWQAAHRCyAAA6IGQBAHRAyAIA6ICQBQDQASELAKADQhYAQAeELACADghZAAAdELIAADogZAEAdEDIAgDogJAFANABIQsAoANCFgBAB4QsAIAOCFkAAB045pBVVVdV1Wer6hNH2P7Sqnqgqm4df7352MsEAJgtSxt47m8m+bUkv/UEY/68tfb9G9gHAMBMOuaZrNbanyW5b4K1AADMja7XZH1nVf11Vf1hVT2v430BAEyNjbQLj+YjSb6+tfbFqnplkvclOXOtgVV1UZKLkuT000/vsCQAgH50NpPVWvt8a+2L49vXJ9lcVScfYewVrbWdrbWd27Zt66okAIDedBayquprq6rGt88a7+vervYHADBNjrldWFXvSvLSJCdX1f4kP59kc5K01t6W5LwkP1lVh5I8nGR3a61tuGIAgBlwzCGrtXbBUbb/WkaXeAAAWDiu+A4A0AEhCwCgA0IWAEAHhCwAgA4IWQAAHRCyAAA6IGQBAHRAyAIA6ICQBQDQASELAKADQhYAQAeELACADghZAAAdELIAADogZAEAdEDIAgDogJAFANABIQsAoANCFgBAB4QsAIAOCFkAAB0QsgAAOiBkAQB0QMgCAOiAkAUA0AEhCwCgA0IWAEAHhCwAgA4IWQAAHRCyAAA6IGQBAHRAyAIA6ICQBQDQASELAKADQhYAQAeELACADghZAAAdELIAADogZAEAdEDIAgDogJAFANABIQsAoANCFgBAB4QsAIAOCFkAAB3YUMiqqquq6rNV9YkjbK+q+pWq2ldVH6uqF2xkfwAAs2KjM1m/meScJ9h+bpIzx18XJfn1De4PAGAmbChktdb+LMl9TzBkV5LfaiN/meQZVfXsjewTAGAWdL0m65Qkd624v3/8GADAXOs6ZNUaj7WvGFR1UVUtV9XygQMHOi4JAKB7XYes/UlOW3H/1CR3rx7UWruitbaztbZz27ZtHZcEANC9rkPWniSvHp9l+MIkD7TW7ul4nwAAg1vayJOr6l1JXprk5Kran+Tnk2xOktba25Jcn+SVSfYleSjJazayPwCAWbGhkNVau+Ao21uS121kHwAAs8gV3wEAOiBkAQB0QMgCAOiAkAUA0AEhCwCgA0IWAEAHhCwAgA4IWQAAHRCyAAA6IGQBAHRAyAIA6ICQBQDQASELAKADQhYAQAeELACADghZAAAdELIAADogZAEAdEDIAgDogJAFANABIQsAoANCFgBAB4QsAIAOCFkAAB0QsgAAOiBkAQB0QMgCAOiAkAUA0AEhCwCgA0IWAEAHhCwAgA4IWQAAHRCyAAA6IGQBAHRAyAIA6ICQBQDQASELAKADQhYAQAeELACADghZAAAdELIAADogZAEAdEDIAgDogJAFANABIQsAoAMbCllVdU5V3V5V+6rqkjW2/2hVHaiqW8dfP7aR/QEAzIqlY31iVW1KcnmSlyfZn+TmqtrTWvvkqqHXttYu3kCNAAAzZyMzWWcl2ddau6O1djDJNUl2TaYsAIDZtpGQdUqSu1bc3z9+bLUfqqqPVdV1VXXaBvYHADAzNhKyao3H2qr770+yvbX2rUk+mOTqNV+o6qKqWq6q5QMHDmygJACA6bCRkLU/ycqZqVOT3L1yQGvt3tbaI+O7v5Hk29d6odbaFa21na21ndu2bdtASQAA02EjIevmJGdW1RlVtSXJ7iR7Vg6oqmevuPuqJLdtYH8AADPjmM8ubK0dqqqLk9yQZFOSq1pre6vq0iTLrbU9Sd5QVa9KcijJfUl+dAI1AwBMvWpt9TKqYe3cubMtLy8PXQYAwFFV1S2ttZ1rbXPFdwCADghZAAAdELIAADogZAEAdEDIAgDogJAFANABIQsAoANCFgBAB4QsAIAOCFkAAB0QsgAAOiBkAQB0QMgCAOiAkAUA0AEhCwCgA0IWAEAHhCwAgA4IWQAAHRCyAAA6IGQBAHRAyAIA6ICQBQDQASELAKADQhYAQAeELACADghZAAAdELIAADogZAEAdGBp6AI61Vpy4EBy++2jr3vvTR55JHn44eTgwWTr1uQpT0m2bEme9azkuc8dfT3zmUNXDgDMuPkJWYcPJx/+cPKnf5rcckvyyU8mn/nM6PHjj08OHUq+9KXksce+8rlLS6Mxxx03GrNlS7J9e/It35K84AXJ2Wcnz39+UtX7twUAzKZqrQ1dwz+yc+fOtry8vL7BBw8mN96YvPOdye/93mjm6uGHR4FqUrZsGX0df3xy3nnJ7t3Ji1+cbNo0uX0AADOpqm5pre1ca9vszWS1lvzhHyZve1vywQ8mmzcnX/jC6PEuHDw4+vriF5Pf+I3kHe8YPf7935/85E8m3/Vd3ewXAJhps7Pw/fFwtWNH8i/+RfL+949mrT7/+e4C1mqPPTYKdF/4QnLNNcm55yYvfGHyF3/Rz/4BgJkx/SGrteQDH0i+9VuT889P/uZvRrNKQ2stefDB5Kabkle8YtRCvOmmoasCAKbEdIes5eXRwvMf/MHkE58YhZpp9NBDyX/5L8nLXpa89KXJf/2vQ1cEAAxsOkPWY48lv/ALyUtektx66/SGq5VaG4WtP//z0ZmIV1zRXxsTAJg603d24bd9W1veujXZu3cUWmbViSeOQuK11yZPe9rQ1QAAHXiiswunL2Rt3tyWq5JHHx26lI3bujV59rNHZ0F+4zcOXQ0AMGFPFLKmr1146NB8BKxkdHX5O+8ctQ9vuWXoagCAHk1fyJo3hw+PLvnw8pcnn/700NUAAD0RsvrywAPJd393cv/9Q1cCAPRAyOrL4cPJ3/3d6JpajzwydDUAQMeErD4dPDg6a/KCC1zeAQDmnJDVt4cfTv74j0cfCwQAzK0NhayqOqeqbq+qfVV1yRrbt1bVtePtN1XV9o3sb248+GDy+tePzqQEAObSMYesqtqU5PIk5ybZkeSCqtqxathrk9zfWntOksuSvOVY9zd37r03ufLKoasAADqykZmss5Lsa63d0Vo7mOSaJLtWjdmV5Orx7euSnF1VtYF9zo8HH0ze9KbZ+MggAOBJ20jIOiXJXSvu7x8/tuaY1tqhJA8kOWkD+5wvjzySvPWtQ1cBAHRgIyFrrRmp1afMrWdMquqiqlququUDGyho5jz0UHLVVUNXAQB0YCMha3+S01bcPzXJ3UcaU1VLSZ6e5L7VL9Rau6K1trO1tnPbBgqaSXffnXz2s0NXAQBM2EZC1s1JzqyqM6pqS5LdSfasGrMnyYXj2+clubFN2ydSD23r1uTGG4euAgCYsGMOWeM1VhcnuSHJbUne3VrbW1WXVtWrxsOuTHJSVe1L8jNJvuIyDwvvC19Ifv/3h64CAJiwmraJpZ1VbXnoIvq2bZuWIQDMoKq6pbW2c61trvg+De691+cZAsCcEbKmwebNrpcFAHNGyJoGmzYJWQAwZ4SsadDaaDYLAJgbQtY0OHQoeepTh64CAJggIWsatJaccMLQVQAAEyRkTYNv+7bkOD8KAJgn/mUf2pYtyQ/8wNBVAAATJmQNbevW5OUvH7oKAGDChKyhPfpo8u3fPnQVAMCECVlDe/GLk6WloasAACZMyBrSCSckl146dBUAQAeErKFs2pR813cl3/mdQ1cCAHRAyBrK5s3JL//y0FUAAB0RsoawZUuye3fy3OcOXQkA0JFqrQ1dwz9SVQeSfGboOp7AyUk+N3QRU8hxWZvjcmSOzdocl7U5Lkfm2Kytr+Py9a21bWttmLqQNe2qarm1tnPoOqaN47I2x+XIHJu1OS5rc1yOzLFZ2zQcF+1CAIAOCFkAAB0Qsp68K4YuYEo5LmtzXI7MsVmb47I2x+XIHJu1DX5crMkCAOiAmSwAgA4IWUdRVedX1d6qOlxVRzxLoao+XVUfr6pbq2q5zxqH8CSOyzlVdXtV7auqS/qscQhV9cyq+kBVfWr851cfYdxj4/fKrVW1p+86+3S090BVba2qa8fbb6qq7f1X2b91HJcfraoDK94nPzZEnX2rqquq6rNV9YkjbK+q+pXxcftYVb2g7xqHsI7j8tKqemDF++XNfdc4hKo6rar+pKpuG/+b9MY1xgz2nhGyju4TSf7nJH+2jrH/vLX2T4c+ZbQnRz0uVbUpyeVJzk2yI8kFVbWjn/IGc0mSD7XWzkzyofH9tTw8fq/809baq/orr1/rfA+8Nsn9rbXnJLksyVv6rbJ/T+J349oV75O391rkcH4zyTlPsP3cJGeOvy5K8us91DQNfjNPfFyS5M9XvF8W5YNxDyX52dbaNyd5YZLXrfG7NNh7Rsg6itbaba2124euY9qs87iclWRfa+2O1trBJNck2dV9dYPaleTq8e2rk/zAgLVMg/W8B1Yes+uSnF1V1WONQ1jE3411aa39WZL7nmDIriS/1Ub+MskzqurZ/VQ3nHUcl4XUWruntfaR8e0vJLktySmrhg32nhGyJqcl+eOquqWqLhq6mClxSpK7Vtzfn69888+bZ7XW7klGv/xJvuYI446vquWq+suqmucgtp73wD+Maa0dSvJAkpN6qW446/3d+KFxe+O6qjqtn9Km3iL+vbJe31lVf11Vf1hVzxu6mL6Nlxo8P8lNqzYN9p5Z6mMn066qPpjka9fY9H+21n5vnS/zotba3VX1NUk+UFV/M/6fx8yawHFZazZi5k9nfaLj8iRe5vTx++UbktxYVR9vrf3tZCqcKut5D8zl++Qo1vM9vz/Ju1prj1TVT2Q02/eyziubfov4flmPj2T08S5frKpXJnlfRu2xhVBVT03yniQ/3Vr7/OrNazyll/eMkJWktfY9E3iNu8d/fraqfjejdsBMh6wJHJf9SVb+7/vUJHdv8DUH90THpar+vqqe3Vq7Zzwd/dkjvMbj75c7qur/yeh/X/MYstbzHnh8zP6qWkry9Mx/W+Sox6W1du+Ku7+RBVirtk5z+ffKRq0MFq2166vqP1bVya21uf9Mw6ranFHAekdr7b1rDBnsPaNdOAFVdWJVPe3x20lekdHC8EV3c5Izq+qMqtqSZHeSuT6TLqPv78Lx7QuTfMWMX1V9dVVtHd8+OcmLknyytwr7tZ73wMpjdl6SG9v8X8DvqMdl1ZqRV2W01oTRcXr1+IyxFyZ54PEW/SKrqq99fC1jVZ2V0b/v9z7xs2bf+Hu+MsltrbW3HmHYYO8ZM1lHUVU/mORXk2xL8gdVdWtr7Xur6uuSvL219sokz0ryu+P391KSd7bW/miwonuwnuPSWjtUVRcnuSHJpiRXtdb2Dlh2H/59kndX1WuT3Jnk/CSp0WUufqK19mNJvjnJf6qqwxn9RfjvW2tzGbKO9B6oqkuTLLfW9mT0F+RvV9W+jGawdg9XcT/WeVzeUFWvyujsqfuS/OhgBfeoqt6V5KVJTq6q/Ul+PsnmJGmtvS3J9UlemWRfkoeSvGaYSvu1juNyXpKfrKpDSR5OsnsB/rOSjP6T+iNJPl5Vt44f+7kkpyfDv2dc8R0AoAPahQAAHRCyAAA6IGQBAHRAyAIA6ICQBQDQASELAKADQhYAQAeELACADghZAAAdELIAADogZAEAdEDIAgDogJAFANABIQsAoANCFgBAB4QsAIAOCFkAAB0QsgAAOiBkAQB0QMgCAOiAkAUA0AEhCwCgA0IWAEAHhCwAgA4IWQAAHRCyAAA6IGQBAHRAyAIA6ICQBQDQASELAKADQhYAQAeELACADghZAAAdELIAADqwNHQBq5188slt+/bt6x5/553JgQPd1bPSyScnX//1/ewLABbePfckd9/d/36PPz553vPWNfSWW275XGtt21rbpi5kbd++PcvLy+se/+M/nrz97R0WtMKuXf3tCwAW3i/+YvLmNyet9bvfM85I1plFquozR9o28+3CQ4fmc18AwGwTsqZ0XwCw8PqewZqwmQ9Zjz46n/sCADJM0JrQPmc+ZJnJAgCmkZA1pfsCAGabkDWl+wKAhWdN1rCELABgoqzJGhGyAIBpJGRN6b4AgNk28yHrscfmc18AsPCsyRqWmSwAYKKsyRoRsgCAaTTzIUu7EADmlHbhsIQsAGAazXzI6rOFJ2QBAOs18yHr8OH+9iVkAQDrNfMhS7sQAObUUGuynF044uxCAJhTFr4Pq892YZ/7AgBm28yHLO1CAGAarStkVdU5VXV7Ve2rqkvW2P6SqvpIVR2qqvNWbbuwqj41/rpwUoU/TsgCgDk172uyqmpTksuTnJtkR5ILqmrHqmF3JvnRJO9c9dxnJvn5JN+R5KwkP19VX73xsr+sz+CjXQgArNd6ZrLOSrKvtXZHa+1gkmuS7Fo5oLX26dbax5KsjiHfm+QDrbX7Wmv3J/lAknMmUPc/cAkHAGAarSdknZLkrhX3948fW4+NPHddzGQBANNoPSGr1nhsvc3KdT23qi6qquWqWj5w4MA6X3rETBYAzKl5X5OV0ezTaSvun5rk7nW+/rqe21q7orW2s7W2c9u2bet86RGXcACAObUA18m6OcmZVXVGVW1JsjvJnnW+/g1JXlFVXz1e8P6K8WMTo10IAEyjo4as1tqhJBdnFI5uS/Lu1treqrq0ql6VJFX1z6pqf5Lzk/ynqto7fu59SX4ho6B2c5JLx49NjJksAJhTM94uXFrfvtr1Sa5f9dibV9y+OaNW4FrPvSrJVRuo8Si1dfXKX0nIAgDWyxXfnwQhCwBYr5kPWWayAIBpNPMhq8/gM+MnOQAA69HjJRymVmv9Bp/WzGYBwNyrtS7z+eTNdMh67LGJHYd1qXJBUgBgfWY+ZB3X43dw3HFCFgCwPjMdsg4d6n8m69Ch/vYHAMyumQ5ZfbcLzWQBQI/6/Ee+g/3OfMjqs11oTRYAsF4zHbK0CwGAaTXTIcvZhQAwx4ZqF06IkPUkCFkA0CNrsoajXQgATKuZDllmsgCAaTXzIatPQhYAsF4zHbL6bhc+vk8AoAfWZA1HuxAAmFYzH7IWYZ8AwOyZ6ZDVd+vO2YUA0CPXyRpO3+3Cx/cJAPTAmqzhaBcCANNqpkPWEK077UIAYD1mOmSZyQKAOWZN1nCELABgWglZM7BPAGD2zHTIsiYLAJg4Zxf2P6vUmpksAOiNNVnDeeyxUfDpe58AAEezrpBVVedU1e1Vta+qLllj+9aquna8/aaq2j5+fHNVXV1VH6+q26rqTZMsXrsQAJhWRw1ZVbUpyeVJzk2yI8kFVbVj1bDXJrm/tfacJJclecv48fOTbG2t/ZMk357kf308gE1C3zNZ2oUAsAB6XJN1VpJ9rbU7WmsHk1yTZNeqMbuSXD2+fV2Ss6uqkrQkJ1bVUpKnJDmY5PMTqTxCFgDMtQVYk3VKkrtW3N8/fmzNMa21Q0keSHJSRoHrwST3JLkzyX9ord23wZr/gXYhADCt1hOy1oqRq+ePjjTmrCSPJfm6JGck+dmq+oav2EHVRVW1XFXLBw4cWEdJI2ayAIBptZ6QtT/JaSvun5rk7iONGbcGn57kviQ/nOSPWmuPttY+m+QvkuxcvYPW2hWttZ2ttZ3btm1bd/FCFgDMuSFahj2uybo5yZlVdUZVbUmyO8meVWP2JLlwfPu8JDe21lpGLcKX1ciJSV6Y5G8mUnlGrbu+Q5Z2IQD0ZN7XZI3XWF2c5IYktyV5d2ttb1VdWlWvGg+7MslJVbUvyc8kefwyD5cneWqST2QU1v5za+1jkyreTBYAMK2W1jOotXZ9kutXPfbmFbe/lNHlGlY/74trPT4pfYesw4eFLABgfWb6iu/ahQDAxPnswtGs0uHD/e1PuxAAejTva7KmmXYhADCtZjpkHTrU/0yWdiEAsB4zHbKcXQgAc26GW4YzHbKGmFV69NH+9wkAC2mogGXhe3LwYP/7FLIAgPWY6ZDlA6IBgGklZD1JQ8yeAQCzZ6ZD1hCL0Ps8mxEAFpo1WcMRsgCAaTXTIWuIwCNkAUCPXMJhGEIWAMwx7cLhCFkAMMeOO26YoHXcZOKRkPUkueI7APTETNZwzGQBwBwzkzWcIQJPn5+VCAALTcgajnYhAMwx7cLhmMkCgDlmJms4Q8wqmckCgJ4MNZMlZFn4DgBzzUzWcIZo3QlZANATIWs4ZrIAYI5Z+D4cC98BYI5NaEZpqP0KWU+She8A0JOqYWazzGRpFwLAXDOTNRztQgCYY0MtfN+0aSIvI2TNwD4BYCFZ+D4cIQsA5tgiXMKhqs6pqtural9VXbLG9q1Vde14+01VtX3Ftm+tqg9X1d6q+nhVHT+RyiNkAcBcG2omq692YVVtSnJ5knOT7EhyQVXtWDXstUnub609J8llSd4yfu5Skt9J8hOttecleWmSRydSeYQsAJhrQy1877FdeFaSfa21O1prB5Nck2TXqjG7klw9vn1dkrOrqpK8IsnHWmt/nSSttXtbaxO7CMIQi9AtfAeAnixAu/CUJHetuL9//NiaY1prh5I8kOSkJN+UpFXVDVX1kar6lxsv+cvMZAHAHJvxD4heWseYtb7D1fM5RxqzlOTFSf5ZkoeSfKiqbmmtfegfPbnqoiQXJcnpp5++jpJGhCwAmGMLcJ2s/UlOW3H/1CR3H2nMeB3W05PcN378T1trn2utPZTk+iQvWL2D1toVrbWdrbWd27ZtW3fxQhYAzLGhrvjeY8i6OcmZVXVGVW1JsjvJnlVj9iS5cHz7vCQ3ttZakhuSfGtVnTAOX9+d5JMTqTzWZAHAXJvxmayjtgtba4eq6uKMAtOmJFe11vZW1aVJlltre5JcmeS3q2pfRjNYu8fPvb+q3ppRUGtJrm+t/cFEKo+ZLACYa/MespKktXZ9Rq2+lY+9ecXtLyU5/wjP/Z2MLuMwcUIWAMyxGV/47orvM7BPAFhIMz6TNdMhy5osAJhjQy189wHRZrIAYK4twBXfp5aZLACYY0OFLDNZZrIAYK5Z+D6cIQKPmSwA6ImF78MZIvCYyQKAngw1k6VdaE0WAMw1C9+Ho10IAHPMwvfhaBcCwBzTLhyOswsBYI5Z+D4ca7IAYI65hMNwrMkCgDlmJms41mQBwBwTsrsA3YAAAAxTSURBVIajXQgAc8zC9+FY+A4Ac8xM1nDMZAHAHLPwfThCFgDMsSFmsqqErETIAoC5dtxx/f/DW+VjdRJrsgBgrmkXDsdMFgDMsaHahWayhCwAmGtDzGRZkzXiiu8AMMdcwmGxCFkA0BML34dj4TsAzDHtwuFYkwUAc2yodqGZLCELAOZa1TDtQjNZQhYAzLVFWPheVedU1e1Vta+qLllj+9aquna8/aaq2r5q++lV9cWq+t8nUvWYkAUAc2zer5NVVZuSXJ7k3CQ7klxQVTtWDXttkvtba89JclmSt6zaflmSP9x4uf+YkAUAc2wB2oVnJdnXWrujtXYwyTVJdq0asyvJ1ePb1yU5u2oUA6vqB5LckWTvRCpewXWyAGCOLcDC91OS3LXi/v7xY2uOaa0dSvJAkpOq6sQk/0eSf73xUqeDkAUAPRliJivpdSZrrTi3+js+0ph/neSy1toXn3AHVRdV1XJVLR84cGAdJY2YyQKAOTbUmqwJ7XdpHWP2Jzltxf1Tk9x9hDH7q2opydOT3JfkO5KcV1W/lOQZSQ5X1Zdaa7+28smttSuSXJEkO3funOoYI2QBQE9m/Irv6wlZNyc5s6rOSPLfk+xO8sOrxuxJcmGSDyc5L8mNrbWW5Lu+XHP9qyRfXB2wNsLCdwCYY0Nc8T3pbyartXaoqi5OckOSTUmuaq3trapLkyy31vYkuTLJb1fVvoxmsHZPpLqj0C4EgDk2xEzW4/udgPXMZKW1dn2S61c99uYVt7+U5PyjvMa/Oob6po6QBQA9GeqzC32sjnYhAMy1GZ/JErIAgOm0ANfJmlpmsgBgji3AFd+nlpAFAHNsET4gmi8TsgCgJ0Nd8V27cLjAI2gBQA+Guhjpos9kDRl0hCwA6IGF78M4fHi4y2cMcRFUAFg4C/AB0VNpqKAjZAFAT1wnaxitDTeTpV0IAD2w8H0YQ7ULH983ANCxIdZkWfhuJgsA5p524TCsyQKAOTdUy0q7ULsQAObeEP/YL/pMlnYhANAJM1mLuW8AWChmsvpnJgsAFoCQ1T8L3wFgAQw1ozIBMx2yLHwHgDlnJqt/Q7XstAsBoEdmsvpnJgsAFoCZrP6ZyQKABSBk9c9MFgDQCe3CYfbr7EIA6NEQHxK96DNZyXAzWQBATyx8BwDogDVZAAAd6DtkTfAjZWY2ZA15hp+zCwGgJ0PMZG3aNJGXWVfIqqpzqur2qtpXVZessX1rVV073n5TVW0fP/7yqrqlqj4+/vNlE6kaAFgMEwo8T8rS0kRe5qghq6o2Jbk8yblJdiS5oKp2rBr22iT3t9aek+SyJG8ZP/65JP9Ta+2fJLkwyW9PpGoAYDH0fXZha73OZJ2VZF9r7Y7W2sEk1yTZtWrMriRXj29fl+TsqqrW2kdba3ePH9+b5Piq2jqJwrULAWAB9D2T1XPIOiXJXSvu7x8/tuaY1tqhJA8kOWnVmB9K8tHW2iPHVioAsHBmuF24nldZa8XZ6rmcJxxTVc/LqIX4ijV3UHVRkouS5PTTT19HSQDAQpjzduH+JKetuH9qkruPNKaqlpI8Pcl94/unJvndJK9urf3tWjtorV3RWtvZWtu5bdu2J/cdAADza87bhTcnObOqzqiqLUl2J9mzasyejBa2J8l5SW5srbWqekaSP0jyptbaX0yk4jFrsgBgAcxzyBqvsbo4yQ1Jbkvy7tba3qq6tKpeNR52ZZKTqmpfkp9J8vhlHi5O8pwk/1dV3Tr++pqJVD4QH+UDAD0aImT1uCYrrbXrk1y/6rE3r7j9pSTnr/G8X0zyixusEQBYVPM8kzWttAsBYAEIWcPQugOAOdd3yDp8uL8rvgMADMZMFgBAByY0q7Ruhw8LWdZkAcAC6DtkmckCABZC3+3CqoldZV7IAgCmV98zWRM8q25mQ5Z2IQAsgL5D1gQ/K3FmQxYAsAD6bhcKWQDAQjCTBQDQASGrf9ZkAcACELIAADqweXO/+xOyAICF0PdM1gQX2gtZAMD00i7snzVZALAA+m4XmskazgQvBAsAHI12IQBAB7QL+6ddCAALYMuWfvdnJmtE6w4A5px2IQBAByx8BwDogJDVP2uyAGABbNrU7/ogIQsAWAibNk30jL917W9ChCwAYHotLZnJ6pt2IQAsADNZAAAd6HtN1gQvGSFkAQDTq+92oZAFACyEeT+7sKrOqarbq2pfVV2yxvatVXXtePtNVbV9xbY3jR+/vaq+d1KFW5MFAAtgntuFVbUpyeVJzk2yI8kFVbVj1bDXJrm/tfacJJclecv4uTuS7E7yvCTnJPmP49cDADi6OW8XnpVkX2vtjtbawSTXJNm1asyuJFePb1+X5OyqqvHj17TWHmmt/bck+8avBwBwdPM8k5XklCR3rbi/f/zYmmNaa4eSPJDkpHU+d+ZoFwJATya4Rmpdeg5Za8XH1THjSGPW89xU1UVVtVxVywcOHFhHScPqM1ADwEKb83bh/iSnrbh/apK7jzSmqpaSPD3Jfet8blprV7TWdrbWdm7btm391QMA823O24U3Jzmzqs6oqi0ZLWTfs2rMniQXjm+fl+TG1lobP757fPbhGUnOTPJXkykdAJh7fbcLN2+e2EsdNa611g5V1cVJbkiyKclVrbW9VXVpkuXW2p4kVyb57aral9EM1u7xc/dW1buTfDLJoSSva609NrHqAYD5NsGZpb73t65Xaq1dn+T6VY+9ecXtLyU5/wjP/TdJ/s0GalzTkOuirMkCgJ7MebsQAGAYM9wuFLIAgOk1w+3CmQ1Z2oUAsADMZAEAdEDIAgDoQN/tQiELAFgIc/6xOlPJmiwAWACbNvX3ocFVEw11MxuyAIAF0Ge7sMpM1uP6CrYAwED6bBced5yZrES7EAAWgnYhAEAHlpb6DVnahcPRogSAHmkX9k+7EAAWgHYhAEAHtAsBADpgJgsAoANCVv+syQKABdBnu/Dx/U3IzIYsAGABbNqUHD7cz77MZAEAC0O7sH/ahQCwALQLAQA60Ge78PH9TYiQBQBML+1CAIAOVPW7Tke70JosAFgYx/UYV8xkAQALo8/ZDSFrpM+TDQCAgUww+ByVdqF2IQAsDDNZAAAdsCYLAKADfbYLJ7ivalO2sKmqDiT5TA+7OjnJ53rYD/3xM51Pfq7zx890/izyz/TrW2vb1towdSGrL1W13FrbOXQdTI6f6Xzyc50/fqbzx890bdqFAAAdELIAADqwyCHriqELYOL8TOeTn+v88TOdP36ma1jYNVkAAF1a5JksAIDOLGTIqqpzqur2qtpXVZcMXQ8bU1WnVdWfVNVtVbW3qt44dE1MRlVtqqqPVtXvD10LG1dVz6iq66rqb8a/r985dE1sXFX9b+O/ez9RVe+qquOHrmlaLFzIqqpNSS5Pcm6SHUkuqKodw1bFBh1K8rOttW9O8sIkr/MznRtvTHLb0EUwMb+c5I9aa/9Dkm+Ln+3Mq6pTkrwhyc7W2rck2ZRk97BVTY+FC1lJzkqyr7V2R2vtYJJrkuwauCY2oLV2T2vtI+PbX8joL+5Thq2KjaqqU5N8X5K3D10LG1dVX5XkJUmuTJLW2sHW2v83bFVMyFKSp1TVUpITktw9cD1TYxFD1ilJ7lpxf3/8gzw3qmp7kucnuWnYSpiA/zvJv0xyeOhCmIhvSHIgyX8et4DfXlUnDl0UG9Na++9J/kOSO5Pck+SB1tofD1vV9FjEkLXWR3k7xXIOVNVTk7wnyU+31j4/dD0cu6r6/iSfba3dMnQtTMxSkhck+fXW2vOTPJjEmtgZV1VfnVE36IwkX5fkxKr6X4atanosYsjan+S0FfdPjanNmVdVmzMKWO9orb136HrYsBcleVVVfTqjlv7Lqup3hi2JDdqfZH9r7fFZ5usyCl3Mtu9J8t9aawdaa48meW+S/3HgmqbGIoasm5OcWVVnVNWWjBbo7Rm4JjagqiqjdR63tdbeOnQ9bFxr7U2ttVNba9sz+h29sbXmf8czrLX2d0nuqqrnjh86O8knByyJybgzyQur6oTx38VnxwkN/2Bp6AL61lo7VFUXJ7kho7Mgrmqt7R24LDbmRUl+JMnHq+rW8WM/11q7fsCagK/0+iTvGP8H944krxm4HjaotXZTVV2X5CMZnen90bj6+z9wxXcAgA4sYrsQAKBzQhYAQAeELACADghZAAAdELIAADogZAEAdEDIAgDogJAFANCB/x87vAU8ObfC8wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x720 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.]]\n"
     ]
    }
   ],
   "source": [
    "psi=opt_result.states[0]\n",
    "circle1=plt.Circle((1/np.sqrt(2)*np.real(psi[nt-1][2]), 1/np.sqrt(2)*np.real(psi[nt-1][4])), 0.25, color='r')\n",
    "circle2=plt.Circle((1/np.sqrt(2)*target[2], 1/np.sqrt(2)*target[4]), 0.25, color='b')\n",
    "\n",
    "fig, ax = plt.subplots(2,figsize=(10,10)) \n",
    "\n",
    "x=[1/np.sqrt(2)*np.real(psi[nt-1][2]),1/np.sqrt(2)*target[2]]\n",
    "y=[1/np.sqrt(2)*np.real(psi[nt-1][4]),1/np.sqrt(2)*target[4]]\n",
    "\n",
    "overlap1=np.exp(-np.abs(1/np.sqrt(2)*target[2]+1j*1/np.sqrt(2)*target[4]-(1/np.sqrt(2)*np.real(psi[nt-1][2])+1j*1/np.sqrt(2)*np.real(psi[nt-1][4])))**2)\n",
    "\n",
    "ax[0].scatter(x,y)\n",
    "ax[0].add_patch(circle1)\n",
    "ax[0].add_patch(circle2)\n",
    "#ax[0].set_ylim([-0.4, 0.3])\n",
    "#ax[0].set_xlim([0.6, 2.3])\n",
    "\n",
    "\n",
    "x=[1/np.sqrt(2)*np.real(psi[nt-1][1]),1/np.sqrt(2)*target[1]]\n",
    "y=[1/np.sqrt(2)*np.real(psi[nt-1][3]),1/np.sqrt(2)*target[3]]\n",
    "\n",
    "\n",
    "overlap2=np.exp(-np.abs(1/np.sqrt(2)*target[1]+1j*1/np.sqrt(2)*target[3]-(1/np.sqrt(2)*np.real(psi[nt-1][1])+1j*1/np.sqrt(2)*np.real(psi[nt-1][3])))**2)\n",
    "\n",
    "circle3=plt.Circle((1/np.sqrt(2)*np.real(psi[nt-1][1]), 1/np.sqrt(2)*np.real(psi[nt-1][3])), 0.25, color='r')\n",
    "circle4=plt.Circle((1/np.sqrt(2)*target[1], 1/np.sqrt(2)*target[3]), 0.25, color='b')\n",
    "\n",
    "ax[1].scatter(x,y)\n",
    "ax[1].add_patch(circle3)\n",
    "ax[1].add_patch(circle4)\n",
    "#ax[1].set_ylim([-1, 1])\n",
    "#ax[1].set_xlim([-1, 1])\n",
    "\n",
    "plt.show()\n",
    "print(1-overlap1*overlap2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Energy vs time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1f3fa4c04c8>]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD4CAYAAADFAawfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAeMklEQVR4nO3deXxU5aHG8d+bfSGBhIQ9EEAWQfaIAkrdiyu1Sou2VqUtVVur1t5q6221tZterbVqVdyoSl1Ra62tUhEUlX2VPYGQBAgJxOz7zHv/mAnGmJBJyMyZzDzfz2c+mZlzJvNwZvJw5p2zGGstIiISvCKcDiAiIsemohYRCXIqahGRIKeiFhEJcipqEZEgF+WPX5qWlmYzMzP98atFRELSunXrDltr01ub5peizszMZO3atf741SIiIckYs6+taRr6EBEJcipqEZEgp6IWEQlyKmoRkSCnohYRCXIqahGRIKeiFhEJcipqEZHjZK1lxe7DPLY8xy+/3y87vIiIhIs1uSXc985OVu0tYVBKPNdMzyQuOrJLn0NFLSLSCbUNLm5fvJk3Nh4gPSmWX18ylrlTM4iN6tqSBhW1iEiHFVfUMf+5tWzML+XHZ4/g+q8MJz6m6wu6iYpaRKQDdh+q4NqFazhcWcej35rMrJP6+/05VdQiIj4qra7nO0+vptFteWn+NCZk9ArI86qoRUR8YK3ltsWbOVxZx2vXz2DcoJ4Be25tnici4oNFq/J4Z+shbps1OqAlDSpqEZF27Sys4O63tvGVkenMmzE04M+vohYROQaX2/KTlzeSFBfNfXMmEBFhAp5BRS0icgxvbNjP1gPl/PKiE0lPinUkg4paRKQNtQ0u/rRkF+MG9uTi8QMcy6GiFhFpw7Of5LK/tIafnz/akSGPJipqEZFWlFbX8/DSbM4Ylc70E9IczaKiFhFpxV+X5VBR18hts0Y7HcW3ojbG3GKM2WqM+dQY84IxJs7fwUREnFJYVsvCj3O5dNJATuyf7HSc9ovaGDMQ+DGQZa09CYgE5vo7mIiIU/66LBu323LLOSOdjgL4PvQRBcQbY6KABOCA/yKJiDhnf2kNL67OZ05WBhmpCU7HAXwoamvtfuA+IA84CJRZa99tOZ8xZr4xZq0xZm1xcXHXJxURCYCHl+4G4MazTnA4yed8GfpIAWYDQ4EBQKIx5tst57PWLrDWZllrs9LT07s+qYiIn+UdqeaVtQVcMTWDAb3inY5zlC9DH+cAe621xdbaBuA1YLp/Y4mIBN5flu4mMsJww5nBszYNvhV1HnCqMSbBGGOAs4Ht/o0lIhJY2w6U89r6Aq46dQh9k4NrwzZfxqhXAa8C64Et3scs8HMuEZGAsdZy15tb6ZUQw41njXA6zpf4dOIAa+2dwJ1+ziIi4og3Nx1gdW4Jf/z6OHomRDsd50u0Z6KIhLWqukZ+//Z2xg/qyTeyMpyO0yqdiktEwtpDS7M5VF7Ho9+e4uiBl45Fa9QiEra2HSjnqRV7uGzyICYPTnE6TptU1CISlmobXPzk5Y30jI/hjgtPdDrOMWnoQ0TC0p+W7GJHYQXPXHMyqYkxTsc5Jq1Ri0jYWbnnCE98uIcrTxnMmaP7OB2nXSpqEQkrZdUN3PryJoakJnDHBcE95NFEQx8iEjZcbstNL22gqKKWl38wjcTY7lGB3SOliEgXeGDJLpbtLOZ3l57EpCDeyqMlDX2ISFj4z6cHefj9bOaenMGVUwc7HadDVNQiEvK2Hyzn1pc3MTGjF7+ePRbP8eW6DxW1iIS0/aU1XPPMapLionns21OIjYp0OlKHqahFJGSVVTdwzdOrqa5zsXDeyfTrGVyHL/WVvkwUkZBU1+hi/nNryT1Sxd/mTWV0P+fPJt5ZKmoRCTmNLjc3/n0Dq/aW8ODciUwfnuZ0pOOioQ8RCSlut+W2xVt4d9sh7rx4DLMnDnQ60nFTUYtIyLDWcve/trF4fQE3nzOCa2cMdTpSl1BRi0jI+Mt72TzzUS7XzsjkprOD75RanaWiFpGQsPCjvTzw311cNnkQv7xwTLfbVvpYVNQi0u29vqGAu/65jXPH9OWey8YF7ZlaOktFLSLd2tIdh/jpK5uZNqw3D10xiajI0Ku10PsXiUjY2Jhfyg8XbWBM/2SeuDqLuOjut9ehL1TUItIt7TtSxXcXriEtKYanrzmZHt3kkKWdoaIWkW6npKqea55Zg8taFl47lfSkWKcj+ZWKWkS6lQaXm+ueX8eB0hqeujqL4ek9nI7kd6H7WUFEQtJv/rmN1d5dw6cMSXU6TkBojVpEuo0XVufx3Mp9/GDmsJDYNdxXKmoR6RbW7SvhV//4lJkj0/nZrNFOxwkoFbWIBL3S6npu/PsG+veM56G5k4gMsR1a2qMxahEJatZafvbqZoor61h8/XR6JkQ7HSngtEYtIkHt+VV5vLvtED/76mjGD+rldBxHqKhFJGjtKCzn7re2MXNkOt89LTQOWdoZKmoRCUoNLje3vLSJ5Lho7p8zIeQOtNQRGqMWkaD0+PIcth8s5/GrpoT8noft0Rq1iASd3Ycq+Mt72Vw4rj9fHdvP6TiOU1GLSFBxuS23Ld5MQmwkd10y1uk4QcGnojbG9DLGvGqM2WGM2W6MmebvYCISnp79JJf1eaXcefGYsB/yaOLrGPWDwH+stZcbY2KABD9mEpEwVVRRy/3v7mLmyHS+Fka7iLen3aI2xiQDM4FrAKy19UC9f2OJSDi69z87qWt0cdfFoXXOw+Ply9DHMKAYeMYYs8EY86QxJrHlTMaY+caYtcaYtcXFxV0eVERC28b8Ul5dV8C804YyLAwOXdoRvhR1FDAZeNRaOwmoAm5vOZO1doG1Nstam5Went7FMUUklLndljvf3Ep6Uiw3njXC6ThBx5eiLgAKrLWrvLdfxVPcIiJdYvH6Ajbll/Lz80eH9Cm1OqvdorbWFgL5xphR3rvOBrb5NZWIhI2aehf/985OJg3upS8Q2+Drf103Aou8W3zsAa71XyQRCSfPfLyXooo6HvnW5LDeTfxYfCpqa+1GIMvPWUQkzJRW1/PoshzOHt2HkzPD47RanaE9E0XEMY8uz6GyrpH/mTWq/ZnDmIpaRBxRWFbLwo9yuXTiQEb3S3Y6TlBTUYuIIx58bxdua7nl3JFORwl6KmoRCbj8kmpeWVvAlVMHk5GqI1K0R0UtIgH3yPvZREQYbjjzBKejdAsqahEJqPySal5dV8AVJ2fQNznO6TjdgopaRALq0eU5RBjDdWcMdzpKt6GiFpGA2V9awytr8/nmyRn07xnvdJxuQ0UtIgHz6LJsAK1Nd5CKWkQC4lB5LS+vKWBOVgYDe2ltuiNU1CISEE+t2Euj2811M7U23VEqahHxu7LqBhat3MfFEwYwuLe2m+4oFbWI+N2zn+RSVe/iuq9obbozVNQi4lc19S6e+TiXs0b34cT+OqZHZ6ioRcSvXlqTR0lVPTdoS49OU1GLiN80uNw88eFeTs5MIUvHm+40FbWI+M1bmw+wv7RGY9PHSUUtIn5hreXx5XsY0acHZ47q43Scbk1FLSJ+8eHuw+worOD7M4fpXIjHSUUtIn7x+Ac59E2OZfbEAU5H6fZU1CLS5T7dX8ZH2Ue4dsZQYqMinY7T7amoRaTLLfhgDz1io7jylMFORwkJKmoR6VL5JdX8a8tBrpiaQXJctNNxQoKKWkS61FMr9mKAeacNdTpKyFBRi0iXKa2u56U1+VwycYBODNCFVNQi0mWeX7mPmgYX3z99mNNRQoqKWkS6RG2Di4Uf72PmyHQdfKmLqahFpEu8vmE/hyvr+MFMrU13NRW1iBw3t9vyxId7GDsgmenDezsdJ+SoqEXkuP13+yH2FFcxf+YwjNHu4l1NRS0ix8Vay2PLc8hIjefCcf2djhOSVNQiclzW5H7G+rxSvn/6MKIiVSn+oKUqIsflseU5pCbGMGdKhtNRQpaKWkQ6bWdhBUt3FHH1tEziY3TwJX9RUYtIpz3+QQ7x0ZF8Z9oQp6OENBW1iHTKgdIa3tx4gLlTM0hJjHE6TkjzuaiNMZHGmA3GmLf8GUhEuocFH+wB4HvaXdzvOrJGfROw3V9BRKT7OFxZx4tr8rh00kAG9tLBl/zNp6I2xgwCLgSe9G8cEekOnl6xl7pGN9edobOLB4Kva9R/Bn4GuP2YRUS6gbKaBp77ZB8XnNSf4ek9nI4TFtotamPMRUCRtXZdO/PNN8asNcasLS4u7rKAIhJcnl+5j4q6Rm44U2vTgeLLGvUM4BJjTC7wInCWMeb5ljNZaxdYa7OstVnp6eldHFNEgkFNvYunVuzlzFHpjB3Q0+k4YaPdorbW/txaO8hamwnMBZZaa7/t92QiEnQWrdpHSVU9PzzzBKejhBVtRy0iPqltcPH4B3uYPrw3WZmpTscJKx0qamvtMmvtRf4KIyLBa9GqPIor6rjp7BFORwk7WqMWkXbVNrh4bHkO04b15pRhOjFAoKmoRaRdf29amz5Ha9NOUFGLyDE1rU2fOiyVU7U27QgVtYgc06JVeRRV1HHT2SOdjhK2VNQi0qaqukb++n4204f3ZppOWusYFbWItOmZj/ZypKqe//nqKKejhDUVtYi0qrS6nsc/2MO5Y/oyaXCK03HCmopaRFr12PI9VNY1cut5Gpt2mopaRL6kqLyWhR/vZfaEAYzul+x0nLCnohaRL/nze7tpdFluPkdr08FARS0iX7DrUAUvrs7j26cOITMt0ek4gopaRFr4/dvbSYyN4sc6pkfQUFGLyFEf7i5m2c5ibjzrBFJ1ZvGgoaIWEQBcbsvv/rWdjNR4rp6e6XQcaUZFLSIAvLounx2FFdw2azSxUZFOx5FmVNQiQll1A/f8ZydThqRw4bj+TseRFlTUIsJ97+6ktLqeu2efhDHG6TjSgopaJMxtKSjj+VX7+M60TMYM0M4twUhFLRLG3G7LL//xKb0TY7nlXO3cEqxU1CJh7OW1+WzML+UXF4ymZ3y003GkDSpqkTBVWFbL797eztShqVw6aaDTceQYVNQiYchayx2vb6HB5ebey8brC8Qgp6IWCUP/2HiA93YU8dPzRul4Ht2AilokzBRX1HHXP7cyeXAvrp0x1Ok44gMVtUgYsdbyi9e3UF3v4t7LJxAZoSGP7kBFLRJGnl+Vx5Jth/jZV0dxQp8eTscRH6moRcLEjsJy7n5rG2eMSmeehjy6FRW1SBioqXfx4xc2kBwXzX1zJhChIY9uJcrpACLif795ayu7DlXy7LyppPWIdTqOdJDWqEVC3N9X5fHC6nyuP2M4M0emOx1HOkFFLRLC1u0r4c43P2XmyHR+et4op+NIJ6moRUJUYVkt1z2/noG94nlo7iRtiteNaYxaJARV1zfyg+fWUl3XyKLvnULPBB1wqTtTUYuEmAaXmxsWrWfL/jIWXJXFyL5JTkeS46SiFgkh1lpuX7yFZTuL+cPXx3HOmL5OR5IuoDFqkRBy7zs7Wby+gFvOGckVUwc7HUe6iIpaJEQ8sGQXjy7L4cpTBvPjs09wOo50oXaL2hiTYYx53xiz3Riz1RhzUyCCiYjvHliyiwff282cKYP4rU5QG3J8GaNuBG611q43xiQB64wxS6y12/ycTUTaYa3lz//dzYPv7ebyKYO457Lx2j08BLVb1Nbag8BB7/UKY8x2YCCgohZxkNtt+c1b21j4ca5KOsR1aKsPY0wmMAlY1cq0+cB8gMGD9SWGiD/VNbr4ycub+Nfmg3z3tKHcccGJKukQ5vOXicaYHsBi4GZrbXnL6dbaBdbaLGttVnq6jicg4i9l1Q1c+8wa/rX5ID8/fzT/e6FKOtT5tEZtjInGU9KLrLWv+TeSiLRl96EKvv/sWvaX1nD/nAlcNmWQ05EkANotauP5+vgpYLu19k/+jyQirVmy7RA3v7iB+JgoXpx/KlOGpDodSQLElzXqGcBVwBZjzEbvfb+w1r7tv1gi0qTB5ea+d3fy+PI9jB/Uk8evmkL/nvFOx5IA8mWrjxWABsBEHJBfUs2NL2xgY34pV54ymF9dNIa46EinY0mA6VgfIkHIWstr6/dz1z+3AvDIlZO5cHx/h1OJU1TUIkGmsKyWX7y+haU7ijg5M4U/fWMiGakJTscSB6moRYKEy215YXUe9/xnBw0uN7+6aAzXTM/UpneiohYJBlsKyvjfN7awqaCMacN684evjyMzLdHpWBIkVNQiDioqr+X+d3fx8rp80nrE8uDciVwyYYAOqiRfoKIWcUB1fSNPfriXx5bn0OByM2/GUG46ZwTJcTpllnyZilokgOoaXbywKo+H38/hcGUd55/Uj9tmjdYwhxyTilokAOob3SxeX8DDS7PZX1rDKUNTefyqydq7UHyiohbxo9oGF6+szefRZTkcKKtlQkYv/njZOE47IU3j0OIzFbWIH1TUNrBoVR5PrdhLcUUdU4ak8IfLxjNzhApaOk5FLdKFispreebjXJ5fuY+K2kZOOyGNB785kWnDe6ugpdNU1CJdYPehCp74cA9vbDhAg9vNrLH9uP6M4Ywf1MvpaBICVNQinWStZUX2YZ78cC/LdxUTFx3B3KkZfPe0oQzpra04pOuoqEU6qK7RxZsbD/DUir3sKKwgrUcst547km+dOoTUxBin40kIUlGL+OizqnqeX7mPv32yj8OVdYzqm8S9l49n9sQBxEbp0KPiPypqkXbkHq7iqRV7eWVdPrUNbmaOTOf7pw/VJnYSMCpqkTZsLijlseU5/PvTQqIjIpg9cQDfO30Yo/olOR1NwoyKWqQZay2f5BzhkWXZfJR9hKS4KK7/ynCumZ5Jn+Q4p+NJmFJRi+Ap6GU7i3lo6W7W55XSJymWX1wwmiumDiZJB0oSh6moJaxZa3l/ZxF//u9uNheUMbBXPL/92klcPmWQzk0oQUNFLWGpaQ36gf/uYnNBGRmp8dx72XgunTyQ6MgIp+OJfIGKWsKKtZaPso9w/5KdbMgrZVCKClqCn4pawsYnOUd4YMkuVueWMKBnHL+/dByXTxlETJQKWoKbilpCmrWWT/Yc4S/v7WblnhL6Jsfy60vGMndqhnZSkW5DRS0hqelLwoeXZrM+r5T0pFjuvHgMV0wdrC8JpdtRUUtIqW1w8caG/Tz90V52HapkYK947v7aSczRVhzSjamoJSTsO1LFi2vyeXlNPkeq6jmxfzL3z5nAJRMH6EtC6fZU1NJtVdY18u7WQhavL+Cj7CNEGDhrdF/mnZbJtGE6UL+EDhW1dCsVtQ0s31XM21sO8t72Iuoa3QxKiefWc0cyJyuDfj21m7eEHhW1BDVrLdlFlXy4+zDv7yxi5Z4jNLgsvRNj+ObJGVwyYQCTB6cQEaG1ZwldKmoJKtZacoorWZP7GWv2lvBxzhEKy2sBGJaWyLUzhnLumL5MHpxCpMpZwoSKWhxVVFHL1gPlbMovZWN+KZvyS/msugGAtB4xnDK0N6eNSOO0E9LISE1wOK2IM1TUEhC1DS72Hq5iZ2EFOwor2FFYztYD5RRX1AFgDIzsk3R0bXnq0FSGpiXqC0ERVNTShay1fFbdwN7DleQUV7GnuIqc4kqyiyrZd6QKt/XMFx1pGJ7eg9NHpDF2QE/GDkjmpIE96RGrt6NIa/SXIR3idluKKurIK6n2XI5Usa+kmtwj1eQerqKspuHovNGRhiG9ExndL4mLx/dneJ8ejO6XzLD0RG3bLNIBKmr5goraBg6V11FYVsuBshoOltZyoLSG/c0u9Y3uo/NHGBjQK54hvRO4eEJ/MnsnMjQtkeHpPRiUEk+UClnkuPlU1MaYWcCDQCTwpLX2j35NJV3CWktNg4vymkY+q67ns+p6SqsbKKmqp6SqniOVdRyurKe4oo7iyjqKymupqnd96fek9YhlYEo8Y/onc96YvgxKTWBwagIZKfEMSknQ0edE/KzdojbGRAKPAOcCBcAaY8yb1tpt/g4XDKy1uNwWtwW3tVgLLms9192e+zwX77ze62735/O63Pbzad7HuNxNty2Nbs/tRrel0eWmwdV02019o5t6l5uGRjd1jZ7btY0uahvc1Da4qGlwUVPvorreRXV9I5V1LirrGqisbaSitpHGpoHhViTHRdG7RyzpSbGMGZDMGaPS6ZccR1/vZWCvePr2jNVR5kQc5ssa9VQg21q7B8AY8yIwG+jyor74oRXUNnx5ja5l1VhrvzjNen7apsLEU5LWfn6f29ovzOP2lqRtut6siJuXbzCKMBAfHUmc95IQE0lCbBQJ0ZEM7BVDUlwSibGRJMdFkxwfTXJcNL0SPJeUhBhSEmJITYzRmrBIN+FLUQ8E8pvdLgBOaTmTMWY+MB9g8ODBnQozPD2Repe71WkG0/KOL1w1xjNHhPn8ujHGexsijPHcbyDy6P3m6LQIAxERxjMfEBlhjv6eyIjP52+67nnM5/dHNLs/svnvjfDeF2G+8JjICPP5xRiiIiOIjDBERRiiIg1RERFERRpiIiOIiYogOjKC2CjPReO+IuHFl6JubUPWL61rWmsXAAsAsrKyOrUu+ue5kzrzMBGRkObLqlkBkNHs9iDggH/iiIhIS74U9RpghDFmqDEmBpgLvOnfWCIi0qTdoQ9rbaMx5kfAO3g2z3vaWrvV78lERATwcTtqa+3bwNt+ziIiIq3Q5gMiIkFORS0iEuRU1CIiQU5FLSIS5Ezz3bG77JcaUwzs6+TD04DDXRinqyhXxyhXxyhXx4RiriHW2vTWJvilqI+HMWattTbL6RwtKVfHKFfHKFfHhFsuDX2IiAQ5FbWISJALxqJe4HSANihXxyhXxyhXx4RVrqAboxYRkS8KxjVqERFpRkUtIhLkHCtqY8wsY8xOY0y2Meb2VqbHGmNe8k5fZYzJDECmDGPM+8aY7caYrcaYm1qZ5wxjTJkxZqP38it/5/I+b64xZov3Ode2Mt0YY/7iXV6bjTGTA5BpVLPlsNEYU26MubnFPAFZXsaYp40xRcaYT5vdl2qMWWKM2e39mdLGY6/2zrPbGHN1AHL9nzFmh/d1et0Y06uNxx7zNfdDrruMMfubvVYXtPHYY/7t+iHXS80y5RpjNrbxWH8ur1a7IWDvMWttwC94DpeaAwwDYoBNwJgW89wAPOa9Phd4KQC5+gOTvdeTgF2t5DoDeMuBZZYLpB1j+gXAv/GckedUYJUDr2khno32A768gJnAZODTZvfdC9zuvX47cE8rj0sF9nh/pnivp/g513lAlPf6Pa3l8uU190Ouu4Cf+vA6H/Nvt6tztZh+P/ArB5ZXq90QqPeYU2vUR0+Ya62tB5pOmNvcbOBv3uuvAmcbY1o7LViXsdYetNau916vALbjOWdkdzAbeNZ6rAR6GWP6B/D5zwZyrLWd3SP1uFhrPwBKWtzd/D30N+BrrTz0q8ASa22JtfYzYAkwy5+5rLXvWmsbvTdX4jlrUkC1sbx84cvfrl9yef/+vwG80FXP56tjdENA3mNOFXVrJ8xtWYhH5/G+qcuA3gFJB3iHWiYBq1qZPM0Ys8kY829jzNgARbLAu8aYdcZzIuGWfFmm/jSXtv+AnFheAH2ttQfB84cG9GllHqeX2zw8n4Ra095r7g8/8g7JPN3Gx3gnl9fpwCFr7e42pgdkebXohoC8x5wqal9OmOvTSXX9wRjTA1gM3GytLW8xeT2ej/cTgIeANwKRCZhhrZ0MnA/80Bgzs8V0J5dXDHAJ8Eork51aXr5ycrndATQCi9qYpb3XvKs9CgwHJgIH8QwztOTY8gKu4Nhr035fXu10Q5sPa+W+Di0zp4ralxPmHp3HGBMF9KRzH9U6xBgTjeeFWGStfa3ldGttubW20nv9bSDaGJPm71zW2gPen0XA63g+gjbn5EmIzwfWW2sPtZzg1PLyOtQ0/OP9WdTKPI4sN+8XShcB37LegcyWfHjNu5S19pC11mWtdQNPtPF8Ti2vKODrwEttzePv5dVGNwTkPeZUUftywtw3gaZvRy8Hlrb1hu4q3jGwp4Dt1to/tTFPv6axcmPMVDzL8IifcyUaY5KaruP5MurTFrO9CXzHeJwKlDV9JAuANtd0nFhezTR/D10N/KOVed4BzjPGpHg/6p/nvc9vjDGzgNuAS6y11W3M48tr3tW5mn+ncWkbz+fUya7PAXZYawtam+jv5XWMbgjMe8wf35D6+C3qBXi+Oc0B7vDe9xs8b16AODwfpbOB1cCwAGQ6Dc9Hks3ARu/lAuA64DrvPD8CtuL5tnslMD0AuYZ5n2+T97mbllfzXAZ4xLs8twBZAXodE/AUb89m9wV8eeH5j+Ig0IBnDea7eL7TeA/Y7f2Z6p03C3iy2WPned9n2cC1AciVjWfMsuk91rR10wDg7WO95n7O9Zz3vbMZTwH1b5nLe/tLf7v+zOW9f2HTe6rZvIFcXm11Q0DeY9qFXEQkyGnPRBGRIKeiFhEJcipqEZEgp6IWEQlyKmoRkSCnohYRCXIqahGRIPf/6hlhVmk02PUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "psi=opt_result.states[0]\n",
    "Energy=np.zeros(nt)\n",
    "for i in range (0,np.shape(psi)[0]):\n",
    "    Energy[i]=1/2*(np.real(psi[i][9])+(np.real(psi[i][2]))**2+np.real(psi[i][14])+(np.real(psi[i][4]))**2)-1/2\n",
    "    #Energy[i]=1/2*(np.real(psi[i][5])+(np.real(psi[i][1]))**2+np.real(psi[i][12])+(np.real(psi[i][3]))**2)-1/2\n",
    "plt.plot(opt_result.tlist,Energy)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "Quantum object: dims = [[15], [1]], shape = (15, 1), type = ket\\begin{equation*}\\left(\\begin{array}{*{11}c}1.0\\\\0.072\\\\0.045\\\\0.025\\\\0.025\\\\\\vdots\\\\0.0\\\\0.0\\\\0.500\\\\0.0\\\\0.500\\\\\\end{array}\\right)\\end{equation*}"
      ],
      "text/plain": [
       "Quantum object: dims = [[15], [1]], shape = (15, 1), type = ket\n",
       "Qobj data =\n",
       "[[1.        ]\n",
       " [0.07199861]\n",
       " [0.04486285]\n",
       " [0.02484165]\n",
       " [0.02543346]\n",
       " [0.5       ]\n",
       " [0.        ]\n",
       " [0.        ]\n",
       " [0.        ]\n",
       " [0.5       ]\n",
       " [0.        ]\n",
       " [0.        ]\n",
       " [0.5       ]\n",
       " [0.        ]\n",
       " [0.5       ]]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inte.opt_result.states[0][nt-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function to get the states at each t for optimized field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import simpson\n",
    "\n",
    "\n",
    "def field_discrete(lista,t,dt):\n",
    "    return lista[int(t/dt)]\n",
    "\n",
    "def states(psi0,g,mu):\n",
    "    HA = qutip.Qobj(H0(g))\n",
    "    H1 =qutip.Qobj(Hde(mu,g))\n",
    "    H2=qutip.Qobj(He(mu))\n",
    "    states=np.zeros((len(opt_result.tlist),15),dtype=complex)\n",
    "    integral=np.zeros(nt)\n",
    "\n",
    "    for time_index in range(len(opt_result.tlist)-1):\n",
    "        if time_index==0:\n",
    "            integral[time_index]=0\n",
    "        else:    \n",
    "            integral[time_index]=simpson(np.array(opt_result.optimized_controls[0])[:time_index], np.array(opt_result.tlist)[:time_index])\n",
    "    \n",
    "    dt=opt_result.tlist[1]-opt_result.tlist[0]\n",
    "    for i in range (0,len(opt_result.tlist)):\n",
    "        #A=HA+H1*field_discrete(opt_result.optimized_controls[0],opt_result.tlist[i],dt)+H2*field_discrete(integral,opt_result.tlist[i],dt)\n",
    "        A=HA+H1*opt_result.optimized_controls[0][i]+H2*integral[i]\n",
    "        '''if i==len(opt_result.tlist)-1:\n",
    "            A=HA+H2*opt_result.optimized_controls[0][i]+H1*(opt_result.optimized_controls[0][i]-opt_result.optimized_controls[0][i-1])/dt\n",
    "        else:\n",
    "            \n",
    "            A=HA+H2*opt_result.optimized_controls[0][i]+H1*(opt_result.optimized_controls[0][i+1]-opt_result.optimized_controls[0][i])/dt'''\n",
    "\n",
    "        if i==0:\n",
    "            states[i,:]=psi0.trans()\n",
    "        else:\n",
    "            states[i][:]=((-A * dt).expm())(qutip.Qobj(states[i-1][:])).trans()\n",
    "\n",
    "    return states"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9706997877916896"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1-np.abs(inte.coherent_overlap(opt_result.states[0],qutip.Qobj(target)))**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Krotov Optimization Result\n",
       "--------------------------\n",
       "- Started at 2022-01-24 17:18:42\n",
       "- Number of objectives: 1\n",
       "- Number of iterations: 3\n",
       "- Reason for termination: Reached convergence: Loss of monotonic convergence; error decrease < 0\n",
       "- Ended at 2022-01-24 17:18:51 (0:00:09)"
      ]
     },
     "execution_count": 118,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "opt_result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function to get the real field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def field():\n",
    "    integral=np.zeros(nt)\n",
    "    for time_index in range(len(opt_result.tlist)):\n",
    "        if time_index==0:\n",
    "            integral[time_index]=0\n",
    "        else:    \n",
    "            integral[time_index]=simpson(np.array(opt_result.optimized_controls[0])[:time_index], np.array(opt_result.tlist)[:time_index])\n",
    "        \n",
    "    return integral"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Energy vs time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Electric field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-30.65462645179146\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXwU9f3H8dcnm4tcBHJxhJBw3yQQ5EatqEhVxKtKBYpWindrL4+2av3Z/lqt/YlSFRTrgQW0avFWtBVEEBIIEO4QAwlXDiD3ne/vj9lIigkkZDez2Xyej8c8sjszu/PJZPPe2e9857tijEEppZR38rG7AKWUUu6jIa+UUl5MQ14ppbyYhrxSSnkxDXmllPJivnYX0FBkZKSJj4+3uwyllGpXUlNT840xUY0t86iQj4+PJyUlxe4ylFKqXRGRA00t0+YapZTyYhrySinlxTTklVLKi2nIK6WUF9OQV0opL+b2kBeRaSKyR0QyROQ+d29PKaXUKW4NeRFxAIuAy4AhwI0iMsSd21RKKXWKu/vJnwdkGGMyAURkOTAD2OnKjeQePkjmB09S6hdFiX8kJf7R5If0B0cADh8I9HPQyd9BkL+DiOAAeoQH0q1zJ0ICPOoyAaVUR2QMbH4FQmJg4DSXP727U64nkN3gfg4wtuEKIjIfmA8QFxd3ThspPJJBcvbL+Erdt/PKTAAb6waxtm4Yq+tGc8B0+87jIkP86R8dyoCYEAZ3D2NYz84M7BaKn0NPVSil2kDxUVh1N+z7GIZe7ZaQF3d+aYiIXAdcaoz5sfP+bOA8Y8xdja2fnJxszvmK17paKM2H4iNw8iAmay3s/zdSsA+A6ujhFPW5nMNRk8j06c3hwioy80rYm1tCxrFiSqtqAfD39WFQt1CG9ujM0B5hDO0RxsCYYIJMBVQWQ3U5VJdCTSXUVkNdtbXthnx8T02+/uAIAN8A8OsEvoHgF2TdFzm331Up1f7tXAWr7rKy5OJHYMyt4HNuB5gikmqMSW5smbuP5HOAXg3uxwKH3bIlHweExlhTj0RkyJXW/JMHYecq/Ha8TcSGPxIBDA/sDL3GQdcE6BpM3fBgCkvKOF5wjLKT+dSUFOC77TihacWESwkBlIG4+M1QHBAQAv4hEBAKAWHWz8DO1tQpHALDoVMXawrqCp26nrrtG+DaepRSbefYTnhzHnQbAVcvgch+btuUu4/kfYG9wEXAIWATMMsYs6Ox9Vt1JN8chYcg60s4sA4OboCSo1BZAsZ5JB4QZoVrpy6YoAjKfDuTXxvM0Up/ssv8OVDiQ04JlJoAKvGnGgfhIUF07xJM97BOdA8PpHvnQHqG+RMZ7MBhaqyj/ZpK51RufRKoKoXqMutnZQlUFlmfEiqLoKIQyk9CxUmoq2n6d/ELdgZ/wzeBLqemwHAIDLPeMALqf4Zak2+gfopQyi7GwEvTIW8X3JkKwRGtfkrbjuSNMTUicifwMeAAljYV8G2ic08Y+QNrqmcM1FRYTSsOv29nCxDsnHpz6kRCZU0tWfll7MstJjOvlP15JXydV0LW3jJKKmuAOqACf4cPcRFBxEeE0ycqmITIYPpEBtMnKoTIEH/kbCFrjPUmUH4Cyo87f56AsuPW/bIG88uOw9HtznVOnnrTaoo4rE8Q/sHgH2T99Au2mpP8g8C3k3W7fmp43zfwv5fVP84vyHqsXyfruRvsS6VUA2mvw8Gv4IqFLgn4s3HrkXxLuf1I3o2MMeSXWO38WQWlfJNf9u3trIIyqmpOnRQOC/Slb3QIfaOsaUBMCANiQukZ3gkfn1YeYRtjfSqoOAkVzk8GFYXWvKpia15VKVSVOKeyU58sqstPnXOornDeLrPOO7SUj5/15hEQeqpJqv6TRWBn65NGUIT1CSQoEkKird4FwVHg0F5PykuVHYdnkqFrX7j543Nugz+dnW3yHYaIEBUaQFRoAGP7/Pe7c12d4XBh+bdH/vvzSsjMK2XN3jzeTM35dr0gfwf9Y0IZFBPKwG6hDOkRxpAeYYQFtuCoWMQZpmGu+tWsE8v1bwA15c43gPo3gtPeHKrKnD+dTVFVpdabS2Wx9QI//o2zSepE4584xAdCe0B4LwiPg4h+ENEXIgdYk56LUO3ZZ49Yn7Yv/6vLAv5s9EjeZkUV1ew7VsLeY8XsOVrM3mPF7D5azPHSqm/XiesaxPCenRke25kRPTszLLZzy4LfExljnYMoO271iio5Zk3FR6AwB05mw4ksKDr1Jog4rKDvNgx6jIKeo6H7CKuJSClPdzQdnpsE426DaX906VOf6UheQ94DGWPIK65kx5Eidh4uIv1QIemHC8k+Xg5YB+v9o0NI6tWFUb3DSY7vSp/I4LO387dH1eVwPBPy9sCxHdZ0dBsUHbKW+/haPRR6T4C48dbPoK721qxUY5ZdB9lfwz1brc4RLqQh7yVOlFax7VAhW7NPsvngCbYcPElhudVe3jXYn+TeXRjXJ4JxfSIY1C209e37nqz4KBzaDDmbrJ5Sh1KhthIQ6DYcEqZAnwut0PcPsrta1dFlrYO/T4epD8Okn7n86TXkvVRdnSEzv4SUrBNsyjrBpqzjHDxeBkB4kB/j+0QwoV8kE/tGkOCtR/r1aiqt0M/6Er75ArI3WqHv8LeO8PtdBP0uhujBbdN9tK4OCrOhIMN6QyrNtZqlqkpPndswdVazFVjnIhx+zovo/KzbDv8GP+tvO+/X9wbzqX+M479/iuPU7frnrX8e34BTvav8gtqsbbjDMgZevNhqhrxrs1sOOjTkO5BDJ8vZsL+A9ZkFfJWRz+HCCgBiu3Ti/AFRTBkQxcR+kd4/bk91ORz4CvZ/bk25zuGSwmKh/1ToNxUSznfNCeq6OijYZ32qyNlkvdnk77OCvCG/ICtc/QKtrqg+9X8DsU5C1zmvq6irhdoq602qtsa6fS49nJpFTvVwCo6CzrHQpTd0iYeogRA1WD8Jtdau92DFD60uk6PnumUTGvIdlDGGrIIyvszIZ83ePL7KyKe0qhY/hzAmvisXDIzie4Ni6Bvl5Uf5YB1FZayGfZ9C5hdWjx8fX+iZDPGTIH4ixI6xunqeTVWpdV1CTgocXG+9mZQft5YFdrZOCEcPgcj+1oni0O5WF1H/4HOvv67OCvraamfo1zR4U6ix3hjq55na/75fV+N8s6g8dWHetz2giqGswPqUUZpnnfAuOgTU54JYvZt6JEGvsdY+ihmm3Vybq7Yanp0IGLhtvdv2m4a8AqCqpo7UAyf4z95cvtiTx+6jxQDERwTxvUExTB0SzXnxXfH19gHaaqut5pyMT+GbNXA47VR3zs69IGqQ1X3TP8g6+jbmVO+fE1mQv9dqagHokgC9J0LcOCsEI/q1/+aPmiprOJC8XVaPkGPp1htayVFruX8oxI21fu/4ydAjUS9+a8rnj8GaP8ONy2HgZW7bjIa8atShk+V8vusYq3flsn5/AVW1dXTu5Mf3BkVzyZAYzh8YRZB/Bzhiqyy2ej0cToO83ZC72zqarW87R6ymjNAYq7mn+wjryLZHEoR+d3RTr2SMdY4he6P1yeXAOmtfgdUEFTfO+YloMnQfqaEP1hvji5fAiOth5nNu3ZSGvDqr0soa1u7L49OduXy++xgnyqoJ8PVhcv8opg3rxsWDY+gc1AH/cevqAGOdxFT/rSTPCvustfDNWsjfY833Cz51pN97IvQc1fEuYqsqhecmW01rt62zmvHcSENetUhNbR2bsk7w8Y6jfLzjKEcKK/D1Ecb3jWD68O5cOrQbXYP97S5TeZqSXGfof2kd7def7HYEWG35vSdYU6/zWnd+oj14/+ew6QWY+67VndfNNOTVOTPGsDWnkI/Sj/JR+hGyCspw+AgTnIE/bWg3umjgq8aUHbfC/uB6K/yPbLXOZfj4QvdE6D3eGvK713nWiWlvkL8PPv0d7PkAxt0B0/7QJpvVkFcuYYxh55Ei3t92hA+2nwr8if0iuXxEd6YN69b+h1tQ7lNRZLXpH/zKCv9DqVZzBkB4b+sitm7DIWaodb9zrHVl6Jl6ftVUWYPx1fcQKsu3bteP0lpR6BzCu8g5zlKFdZ6lpupUt9Vvr1UQ5/UFvlYvGB+/U91dfQOsk/D1U/13QdR3LzXGOlmdtswatXXyvTDhrjY7N6Ehr1zOGMOOw0W8v/0I7207TPbxcvwdPlwwMIoZiT25aHA0gX7ajq3OoLrCOrrP/hoOpVg9eY5ncqr7JlZgBoQ4w7WT8xqCSiukK53B3RT/EGu00/oRUP2DrefzDXBeGOa8aEx8rG0ac6r7aX1X1doqZ5dT54is1WWnRm6tKj7VywqsN4dRc+CCByAkyl17rVEa8sqtjDGkZZ/k3a1W4OcWVxIS4MulQ7txVVIPJvSNxOHNQywo16kqtcYpKsyxpqJD/z0Uto+vM6T9Tn3JT2B4gwu6Iq0hrDt1cf/JXmOsNwA49SnApusHNORVm6mtM2zILOBfaYf4cPtRiitriAkLYEZiT64e1ZNB3Vw4BLJSCtCQVzapqK7l8925vLU5h//syaOmzjC0RxjXjIplRmIPIkI6WLc6pdxEQ17ZrqCkklVbD/PPzTmkHyrCzyFMHRzD9WN6MaV/lDbnKNUKGvLKo+w+WsQbKTm8veUQx0ur6NE5kOuSe/GDMb3oEa5fAKJUS2nIK49UVVPH6l3HWL4pm7X78hDggoHRzB7Xm/MHRHn3ePhKuZCGvPJ42cfLWLEpmxUp2eQVVxLXNYjZ43pz/ZhedO6kfe+VOhMNedVuVNXU8fGOo7y6/gAbs44T5O/gmlGx/GhiPH2jQuwuTymPpCGv2qUdhwt5aV0Wq9IOU1VbxyVDYlhwQV9Gxbn2+zGVau9sCXkReRi4FchzznrAGPPBmR6jIa8ak1dcyavrs3h5/QEKy6s5L6Erd3+vPxP7RXj/l50o1Qx2hnyJMeaJ5j5GQ16dSWllDcs3ZbNkTSZHiyoY3bsLP53an0n9IjXsVYd2ppBv519hozqS4ABfbpmUwBe/uoBHZwzl8MlyZr+4kRsWb2DLwRN2l6eUR3L3kfyPgCIgBfi5MeY7/4kiMh+YDxAXFzf6wIEDbqlHeZ/KmlqWb8xm4Wf7KCitYvrwbvx62iB6R3j5WOVKncZtzTUishpo7PvPHgQ2APlYQ8o9CnQ3xtx8pufT5hp1Lkoqa3hhbSaL12RSU2u4ZXICd17Yj+CADvDVhUrhAb1rRCQeeM8YM+xM62nIq9Y4VlTBnz7azVubDxEdGsBvLh/CFSO6a3u98nq2tMmLSPcGd2cC6e7allIAMWGBPHl9Im/dPoGYsEDu/scWbv77JnJOnGHMcaW8nDtPvP5ZRLaLyDbgQuBnbtyWUt8aFdeFd+6YyG8vH8LX3xznkr+u4eWvsqir85xrQpRqK3oxlPJqOSfKePDtdL7Ym8ekfpH8+doROgia8jrahVJ1WLFdgvj7vDH8YeZwNh88waX/t4Z/pR2yuyyl2oyGvPJ6IsKssXF8eM9kBsaEcs/yNO775zbKq2rtLk0pt9OQVx1G74hgls8fxx0X9mVFSjZXLVpHRm6x3WUp5VYa8qpD8XX48MtLB/HyvPPIL6lkxjPr+Cj9qN1lKeU2GvKqQ5oyIIr3755Mv5hQFryWyl8+2aO9b5RX0pBXHVa3zoGsmD+O65NjefrzDG59JYXSyhq7y1LKpTTkVYcW6OfgT9eM4PczhvLvPblc99x6jhSW212WUi6jIa86PBFhzvh4XvzRGA4UlHLVonWkHyq0uyylXEJDXimnCwdG8+ZtE/AR4QfPr2ftvryzP0gpD6chr1QDg7uH8c4dE+nVNYh5L23inS164ZRq3zTklTpNTFggKxeMJzm+Cz9dkcaSNZl2l6TUOdOQV6oRYYF+vHzzeXx/eHce+2AXf/xwF540zpNSzaXfqqBUEwJ8HSy8MYkuwX48/0UmJ0ureWzmMHwdemyk2g8NeaXOwOEjPDpjGF2DA1j42T5OlFWx8MYkAv0cdpemVLPoIYlSZyEi3HvxAB66Ygif7DzG3KUbKaqotrsspZpFQ16pZpo3MYGnbkgk9cAJbnh+A3nFlXaXpNRZacgr1QIzEnvywtxkvskv5drnviIrv9TukpQ6Iw15pVrogoHRvH7rWIrKq7nm2a9Iyz5pd0lKNUlDXqlzkBTXhX/eNoGgAAc3Lt7AZ7uO2V2SUo3SkFfqHPWJCuGt2ybSPyaEW19J4e/rvtG+9MrjaMgr1QpRoQH849ZxXDQ4hoff3clv3kmnurbO7rKU+paGvFKtFBzgy/M3jWbB+X1Z9vVB5i7dyInSKrvLUgrQkFfKJXx8hPsuG8QT140kJesEVzzzJTsO63DFyn4a8kq50LWjY1nxk3HU1BquefYrHcVS2a5VIS8i14nIDhGpE5Hk05bdLyIZIrJHRC5tXZlKtR9JcV14965JjIgN56cr0vjtO+lUVNfaXZbqoFp7JJ8OXA2saThTRIYANwBDgWnA30REB/tQHUZUaADLfjyWWycn8OqGA1zzrF44pezRqpA3xuwyxuxpZNEMYLkxptIY8w2QAZzXmm0p1d74OXx48PtDeGFOMjknyrn86S9ZtfWw3WWpDsZdbfI9gewG93Oc875DROaLSIqIpOTl6detKe8zdUgM7989if4xIdz9jy386s2tlFXV2F2W6iDOGvIislpE0huZZpzpYY3Ma/QqEWPMYmNMsjEmOSoqqrl1K9WuxHYJYuVPxnPHhX15IzWHK57+kp2Hi+wuS3UAZx1P3hgz9RyeNwfo1eB+LKCfU1WH5ufw4ZeXDmJC30h+tiKNqxat41fTBnLzxAR8fBo7LlKq9dzVXLMKuEFEAkQkAegPbHTTtpRqVyb2i+Sjn07h/IFR/M/7u5j70kaOFVXYXZbyUq3tQjlTRHKA8cD7IvIxgDFmB7AS2Al8BNxhjNE+ZEo5dQ32Z/Hs0fxh5nA2ZR1n2v+t4dOdOsiZcj3xpAGVkpOTTUpKit1lKNWmMnJLuGf5FnYcLuKHY+P4zfeH0Mlfexyr5hORVGNMcmPL9IpXpWzWLzqEt26fwE+m9GHZ1we5/Om1pB/SIRGUa2jIK+UBAnwd3D99MMt+PJaSyhpm/m0dz3+xn7o6z/mkrdonDXmlPMjEfpF8dM8ULhoUwx8/3M1NL37N0UI9KavOnYa8Uh6mS7A/z940ij9dM5wtB09y2VNr+GTHUbvLUu2UhrxSHkhE+MGYON67exI9wjsx/9VUfvPOdsqrtJOaahkNeaU8WN8o66Ts/Cl9eG3DQa585kt2HdErZVXzacgr5eECfB08MH0wr95yHifLq5mxaJ1+n6xqNg15pdqJyf2j+OieyUzqF8nD7+7k1ldSOK5fM6jOQkNeqXYkIiSAF+cm89AVQ1izN5/LnlrDV/vz7S5LeTANeaXaGRFh3sQE3rp9AsH+vvzwha954uM91NTW2V2a8kAa8kq1U8N6dubduyZx3ehYnvl3Btc/v57s42V2l6U8jIa8Uu1YcIAvf752JE/dkMjeYyVMX7iWD7cfsbss5UE05JXyAjMSe/LB3ZPpExnMbcs285t3tuuXhytAQ14prxEXEcQbC071qb9q0ToycovtLkvZTENeKS/i7+vDA9MH89K8MeQWV3LF0+tYmZKtfeo7MA15pbzQhQOj+fCeyST2CudXb27jpyvSKK6otrssZQMNeaW8VExYIK/9eCw/v3gA7249zBVPf6nj1HdAGvJKeTGHj3DXRf1ZPn88FdV1XP23r3j5qyxtvulANOSV6gDOS+jKB/dMZlL/SB5atYMFr6VSWKbNNx2BhrxSHUTXYH9enJvMg9MH89muXL7/9Fq2HDxhd1nKzTTklepARIRbp/ThzdsmAHDdc+tZsiZTm2+8mIa8Uh1QYq9w3r97MhcNjuaxD3Zx6yvafOOtNOSV6qA6d/LjuZtG87vLh/DF3lymL1xLWvZJu8tSLtaqkBeR60Rkh4jUiUhyg/nxIlIuImnO6bnWl6qUcjUR4eZJCbyxoL75RnvfeJvWHsmnA1cDaxpZtt8Yk+icFrRyO0opN7KabyYxpX8UD63awZ3/2KIXT3mJVoW8MWaXMWaPq4pRStknPMifJXOS+dW0gXy4/Yhz7JsSu8tSreTONvkEEdkiIl+IyOSmVhKR+SKSIiIpeXl5bixHKXU2Pj7C7Rf047Ufj+VkWTVXLVrHR+lH7S5LtcJZQ15EVotIeiPTjDM87AgQZ4xJAu4FXheRsMZWNMYsNsYkG2OSo6Kizu23UEq51IS+kbx71yT6Roew4LVUHv94N3V12k7fHvmebQVjzNSWPqkxphKodN5OFZH9wAAgpcUVKqVs0SO8Eyt/Mo6H/rWDRf/ez56jxfz1B4mEBvrZXZpqAbc014hIlIg4nLf7AP2BTHdsSynlPgG+Dv549XB+P2Mo/96Tx8y/fUVWfqndZakWaG0XypkikgOMB94XkY+di6YA20RkK/AmsMAYc7x1pSql7CAizBkfz2u3jKWgpJKZf1vHpiz9d24vxJP6wyYnJ5uUFG3RUcpTHSgoZd5Lm8g5Uc7j141gRmJPu0tSgIikGmOSG1umV7wqpZqtd0Qwb90+gcS4cO5ZnsYzn+/TC6c8nIa8UqpFwoP8efWW87gqsQdPfLKXh1btoFZ73niss/auUUqp0wX4Onjy+kSiwwJZvCaT/JJKnrw+kUA/h92lqdNoyCulzomPj/DA9MFEhQTw2Ae7OFm2iSVzkgkO0FjxJNpco5RqlVun9OHJ60eyIbOAOUs3UqRj3ngUDXmlVKtdPSqWRbNGsS3nJLOWbOB4aZXdJSknDXmllEtcNrw7i2cns/dYCbOWbKCgpNLukhQa8kopF7pwUDQv/WgM3+SX8sMXvtYjeg+gIa+UcqmJ/SJ5ca4V9LOWbOCEBr2tNOSVUi43qX8kL8xNJtN5RF9Yridj7aIhr5Ryi8n9o1g8ezT7couZ99JGSitr7C6pQ9KQV0q5zQUDo3n6xiTSsk8y/9UUKqpr7S6pw9GQV0q51bRh3Xn82pGsyyjgrn9soaa2zu6SOhQNeaWU210zOpbfzxjKpzuP8cDb23VQszak1x8rpdrEnPHx5JdUsfCzfUSFBvDLSwfZXVKHoCGvlGozP5van7ziShb9ez+RIQHMm5hgd0leT0NeKdVmRIRHZwyloKSS37+3k25hgVw2vLvdZXk1bZNXSrUpX4cPC29MIrFXOD9dkcaWgyfsLsmracgrpdpcoJ+DJXOSiQkL5Mcvp3CwoMzukryWhrxSyhaRIQG8NG8MNXWGH/19I4VlelWsO2jIK6Vs0zcqhMWzR5N9vIw7Xt+sfejdQENeKWWrsX0ieOyq4XyZkc+j7+20uxyvo71rlFK2u35ML/YeK+aFL7+hX0wos8f1trskr9GqI3kReVxEdovINhF5W0TCGyy7X0QyRGSPiFza+lKVUt7s/umDuXBgFA+v2sH6/QV2l+M1Wttc8ykwzBgzAtgL3A8gIkOAG4ChwDTgbyKiX+OulGqSw0dYeGMS8RFB3PH6ZnJOaI8bV2hVyBtjPjHG1I8fugGIdd6eASw3xlQaY74BMoDzWrMtpZT3Cw30Y/GcZKpr6vjJq6mUV+mola3lyhOvNwMfOm/3BLIbLMtxzvsOEZkvIikikpKXl+fCcpRS7VHfqBCeujGRnUeKuP+tbTqYWSudNeRFZLWIpDcyzWiwzoNADbCsflYjT9XoX8oYs9gYk2yMSY6KijqX30Ep5WW+NyiGn188gHfSDrN0XZbd5bRrZ+1dY4yZeqblIjIXuBy4yJx6y80BejVYLRY4fK5FKqU6ntsv6MfWnEL+8MEuhvUIY2yfCLtLapda27tmGvBr4EpjTMOzJKuAG0QkQEQSgP7AxtZsSynVsfj4CH+5fiS9uwZxx+tbOFZUYXdJ7VJr2+SfAUKBT0UkTUSeAzDG7ABWAjuBj4A7jDF6BkUp1SJhgX48N3s0ZVU13PZaKlU1ekVsS7W2d00/Y0wvY0yic1rQYNljxpi+xpiBxpgPz/Q8SinVlAExofz52hFsPniSP3ywy+5y2h0d1kAp5fEuH9GDWyYl8PevsvhX2iG7y2lXNOSVUu3CfZcNYkx8F+7753b2Hiu2u5x2Q0NeKdUu+Dl8eGbWKIIDfFnwWiollTVnf5DSkFdKtR8xYYE8MyuJAwVl/PpNvVCqOTTklVLtyrg+Efzq0oG8v/0IL375jd3leDwNeaVUuzN/Sh8uHRrDHz/czcZvjttdjkfTkFdKtTsiwuPXjSSuaxB3vr6Z3GK9UKopGvJKqXYpLNCPZ28aRVFFNXe+voVq/erARmnIK6XarUHdwvjj1cPZ+M1x/vzRbrvL8Uga8kqpdm1mUixzxvdmydpv+GD7EbvL8Tga8kqpdu833x9CUlw4v3xjKxm5eqFUQxrySql2z9/Xh7/9cBSBfg5+8moqxRXVdpfkMTTklVJeoXvnTjw9K4msgjJ+8cZWvVDKSUNeKeU1JvSN5P7LBvHxjmM8+8V+u8vxCBrySimvcsukBK4Y2YMnPt7Dmr36vdEa8kopryIi/Oma4QyICeWuf2zhYEHZ2R/kxTTklVJeJ8jfl+dnj8YYw/xXUyir6rgjVmrIK6W8Uu+IYBbemMSeY8X8+p/bO+yJWA15pZTXumBgNL+8dCDvbj3MkrWZdpdjCw15pZRXu+38vkwf3o3//XA3X3TAE7Ea8kopryYiPH7tSOtE7OubycovtbukNqUhr5TyesEBviyZk4yPj3DrKykd6qsDNeSVUh1Cr65BLJo1isz8Uu5dkUZdXcc4EduqkBeRx0Vkt4hsE5G3RSTcOT9eRMpFJM05PeeacpVS6txN7BfJg9MH88nOYzz12T67y2kTrT2S/xQYZowZAewF7m+wbL8xJtE5LWjldpRSyiXmTYznutGxPPXZPj5K9/6hiVsV8saYT4wx9Y1bG4DY1peklFLuIyL8z8xhJMWFc+/Krew+WmR3SW7lyjb5m4EPG9xPEJEtIvKFiExu6kEiMl9EUkQkJS+v43VvUkq1vQBfB8/fNJrQQF9ufSWFE6VVdpfkNmcNeRFZLSLpjUwzGqzzIFADLHPOOgLEGX12NO8AAA5aSURBVGOSgHuB10UkrLHnN8YsNsYkG2OSo6KiWv8bKaVUM0SHBfL87GSOFVVy+7LNXvsdsWcNeWPMVGPMsEamfwGIyFzgcuCHxnndsDGm0hhT4LydCuwHBrjv11BKqZZL7BXO/149nPWZBTz63k67y3EL39Y8WESmAb8GzjfGlDWYHwUcN8bUikgfoD/QMa8pVkp5tKtHxbL7aDGL12QyqFsYs8bG2V2SS7Uq5IFngADgUxEB2ODsSTMF+L2I1AC1wAJjzPFWbksppdzi19MGsedoMb/7Vzr9okM4L6Gr3SW5jHjSyGzJyckmJSXF7jKUUh1QYXk1Mxeto7C8mn/dOZHYLkF2l9RsIpJqjElubJle8aqUUkDnTn4smZtMVW0d819J9Zox6DXklVLKqW9UCAtvTGLX0SJ++cY2rxiDXkNeKaUauHBgNPdNG8T724/wzOcZdpfTaq098ep21dXV5OTkUFFRYXcpbhMYGEhsbCx+fn52l6KUAuZP6cPuo8X85dO9DOwWyiVDu9ld0jnz+JDPyckhNDSU+Ph4nD14vIoxhoKCAnJyckhISLC7HKUU1tAHf7x6OJl5JfxsRRpv3T6Rgd1C7S7rnHh8c01FRQURERFeGfBgvZgiIiK8+pOKUu1RoJ+D52cnExzgy49f2dRuhz7w+JAHvDbg63n776dUe9WtcyDPzx7NscJK7nh9MzXtcOiDdhHySilll6S4Lvzh6uF8tb+Axz7YZXc5LaYh3wwOh4PExESGDRvGFVdcwcmTJ8/5ueLj48nPz3dhdUopd7t2dCw3T0zgpXVZrEzJtrucFtGQb4ZOnTqRlpZGeno6Xbt2ZdGiRXaXpJRqYw9MH8TEfhH85u10Nh88YXc5zebxvWsaeuTdHew87NoB/of0COOhK4Y2e/3x48ezbdu2b+8//vjjrFy5ksrKSmbOnMkjjzwCwFVXXUV2djYVFRXcc889zJ8/36V1K6Xalq/Dh2duHMWVi75kwaupvHvXJGLCAu0u66z0SL4Famtr+eyzz7jyyisB+OSTT9i3bx8bN24kLS2N1NRU1qxZA8DSpUtJTU0lJSWFhQsXUlBQYGfpSikX6BLsz5I5yZRU1rDgtVQqa2rtLums2tWRfEuOuF2pvLycxMREsrKyGD16NBdffDFghfwnn3xCUlISACUlJezbt48pU6awcOFC3n77bQCys7PZt28fERERttSvlHKdQd3C+Mt1I7lt2WZ++046f7pmhEf3kNMj+Waob5M/cOAAVVVV37bJG2O4//77SUtLIy0tjYyMDG655Rb+85//sHr1atavX8/WrVtJSkrSfvBKeZHLhnfnru/1Y2VKDq+sP2B3OWekId8CnTt3ZuHChTzxxBNUV1dz6aWXsnTpUkpKSgA4dOgQubm5FBYW0qVLF4KCgti9ezcbNmywuXKllKv9bOoALhoUzaPv7eTrTM9tjtWQb6GkpCRGjhzJ8uXLueSSS5g1axbjx49n+PDhXHvttRQXFzNt2jRqamoYMWIEv/3tbxk3bpzdZSulXMzHR/jrDYnERQRx+7LNHD5ZbndJjfL4Lw3ZtWsXgwcPtqmittNRfk+lvE1GbglXLVpHQmQwbywYT6Cfo81r0C8NUUopN+kXHcJff5DI9kOFPPD2do8bg15DXimlWuniITHcc1F/3tp8iJe/yrK7nP+iIa+UUi5wz0X9mTo4hkff38UGDzoRqyGvlFIu4OMjPPmDkfSOCOIODzoRqyGvlFIuEhbox+LZyVTV1PHjl1M84svANeSVUsqF+kWHsHBWEruPFvHzlVupq7P3RKyGfDPUDzU8dOhQRo4cyZNPPkldnfXlASkpKdx9990AVFZWMnXqVBITE1mxYgVr165l6NChJCYmUl7uGR/dlFLud+HAaB6YPpgP04/yf6v32lpLq8euEZFHgRlAHZAL/MgYc1iswRyeAqYDZc75m1u7PTvUD2sAkJuby6xZsygsLOSRRx4hOTmZ5GSre+qWLVuorq7+dt0FCxbwi1/8gnnz5tlWu1LKHrdMSmDvsWIWfp5Br65BXJfcy5Y6XDFA2ePGmN8CiMjdwO+ABcBlQH/nNBZ41vnz3H14Hxzd3qqn+I5uw+Gy/2326tHR0SxevJgxY8bw8MMP88UXX/DEE0+wdOlSbrrpJvLy8khMTOS2225j5cqVfPzxx6xevZply5a5tm6llEcTEf7nquEcPlnBr/65DV+HMDMpts3raHXIG2MaDvAeDNQ3QM0AXjHWlQEbRCRcRLobY460dpt269OnD3V1deTm5n47Lzo6mhdeeIEnnniC9957D4D169dz+eWXc+2119pVqlLKRv6+PiyZk8wtL2/i5yu34iPCjMSebVqDS4YaFpHHgDlAIXChc3ZPoOH3ZOU45x057bHzgfkAcXFxZ95QC4643c3TrmpTSnmmTv4OXpibzM1/38TPVqRxorSKOePj8fFpm+GJm3XiVURWi0h6I9MMAGPMg8aYXsAy4M76hzXyVN9JRmPMYmNMsjEmOSoq6lx/jzaVmZmJw+EgOjra7lKUUu1AkL8vS380hvMHRPHwuzu5YckGsvJL22TbzTqSN8ZMbebzvQ68DzyEdeTe8ExDLHC4RdV5oLy8PBYsWMCdd97p0V8UoJTyLPVB/0ZqDo++t5NpT63h+uRefG9QNOP6RLhtYDNX9K7pb4zZ57x7JbDbeXsVcKeILMc64VrYXtvj678Zqrq6Gl9fX2bPns29995rd1lKqXZGRLg+uRdT+kfx2Ae7WJmSzSvrD9DJz8EPx8bxm8uHuH6brW1bFpF/AgOxulAeABYYYw45u1A+A0zD6kI5zxiT0vQz6VDDHeH3VEqdUlFdy/rMAv69O5f4iGBunpRwTs9zpqGGXdG75pom5hvgjtY+v1JKeatAPwcXDozmwoHuO7+nV7wqpZQXaxch7+3dFb3991NK2cfjQz4wMJCCggKvDUJjDAUFBQQGBtpdilLKC7nkYih3io2NJScnh7y8PLtLcZvAwEBiY9v+cmellPfz+JD38/MjIeHczjgrpVRH5/HNNUoppc6dhrxSSnkxDXmllPJirb7i1ZVEJA/rqtlzFQnku6gcV9K6Wkbrahmtq2W8sa7exphGR3j0qJBvLRFJaerSXjtpXS2jdbWM1tUyHa0uba5RSikvpiGvlFJezNtCfrHdBTRB62oZratltK6W6VB1eVWbvFJKqf/mbUfySimlGtCQV0opL9buQl5EponIHhHJEJH7GlkeICIrnMu/FpH4Nqipl4j8W0R2icgOEbmnkXUuEJFCEUlzTr9zd10Ntp0lItud2/3Ot3OJZaFzn20TkVFurmdgg/2QJiJFIvLT09Zps/0lIktFJFdE0hvM6yoin4rIPufPLk08dq5znX0iMrcN6npcRHY7/05vi0h4E48949/cDXU9LCKHGvy9pjfx2DP+/7qhrhUNasoSkbQmHuvO/dVoPrTZa8wY024mwAHsB/oA/sBWYMhp69wOPOe8fQOwog3q6g6Mct4OBfY2UtcFwHs27bcsIPIMy6cDHwICjAO+buO/6VGsizls2V/AFGAUkN5g3p+B+5y37wP+1MjjugKZzp9dnLe7uLmuSwBf5+0/NVZXc/7mbqjrYeAXzfhbn/H/19V1nbb8L8DvbNhfjeZDW73G2tuR/HlAhjEm0xhTBSwHZpy2zgzgZeftN4GLnN836zbGmCPGmM3O28XALqCnO7fpYjOAV4xlAxAuIt3baNsXAfuNMa250rlVjDFrgOOnzW74OnoZuKqRh14KfGqMOW6MOQF8ivWdxm6ryxjziTGmxnl3A9DmY1Q3sb+aozn/v26py5kB1wP/cNX2musM+dAmr7H2FvI9gewG93P4bph+u47zn6EQiGiT6gBn81AS8HUji8eLyFYR+VBEhrZVTYABPhGRVBGZ38jy5uxXd7mBpv/x7NpfADHGmCNg/ZMCjX0Jp537DeBmrE9gjTnb39wd7nQ2Iy1tounBzv01GThmjNnXxPI22V+n5UObvMbaW8g3dkR+eh/Q5qzjFiISAvwT+Kkxpui0xZuxmiRGAk8D77RFTU4TjTGjgMuAO0RkymnLbdlnIuIPXAm80chiO/dXc9n5WnsQqAGWNbHK2f7mrvYs0BdIBI5gNY2czrb9BdzImY/i3b6/zpIPTT6skXkt2mftLeRzgF4N7scCh5taR0R8gc6c20fLFhERP6w/4DJjzFunLzfGFBljSpy3PwD8RCTS3XU5t3fY+TMXeBvrY3NDzdmv7nAZsNkYc+z0BXbuL6dj9U1Wzp+5jaxjy35znny7HPihcTbcnq4Zf3OXMsYcM8bUGmPqgCVNbM+u/eULXA2saGodd++vJvKhTV5j7S3kNwH9RSTBeRR4A7DqtHVWAfVnoK8FPm/qH8FVnO19LwK7jDFPNrFOt/pzAyJyHta+L3BnXc5tBYtIaP1trBN36aettgqYI5ZxQGH9x0g3a/Loyq791UDD19Fc4F+NrPMxcImIdHE2T1zinOc2IjIN+DVwpTGmrIl1mvM3d3VdDc/hzGxie835/3WHqcBuY0xOYwvdvb/OkA9t8xpzx9lkd05YPUH2Yp2lf9A57/dYL3qAQKyP/xnARqBPG9Q0Cesj1DYgzTlNBxYAC5zr3AnswOpRsAGY0Eb7q49zm1ud26/fZw1rE2CRc59uB5LboK4grNDu3GCeLfsL643mCFCNdeR0C9Z5nM+Afc6fXZ3rJgMvNHjszc7XWgYwrw3qysBqo61/ndX3JOsBfHCmv7mb63rV+drZhhVe3U+vy3n/O/+/7qzLOf/v9a+rBuu25f5qKh/a5DWmwxoopZQXa2/NNUoppVpAQ14ppbyYhrxSSnkxDXmllPJiGvJKKeXFNOSVUsqLacgrpZQX+3+8PCH3yoPHwQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Efield=field()\n",
    "plt.plot(opt_result.tlist,Efield,label=\"Real\")\n",
    "plt.plot(opt_result.tlist,opt_result.optimized_controls[0],label=\"Diff\")\n",
    "plt.legend()\n",
    "print(Efield[nt-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.+0.j]]\n",
      "[[-0.09224555+0.j]]\n",
      "[[0.02034789+0.j]]\n",
      "[[-0.18902043+0.j]]\n",
      "[[0.17992509+0.j]]\n",
      "[[-1.77635684e-15+0.j]]\n",
      "[[0.+0.j]]\n",
      "[[0.+0.j]]\n",
      "[[0.+0.j]]\n",
      "[[6.10622664e-15+0.j]]\n",
      "[[0.+0.j]]\n",
      "[[0.+0.j]]\n",
      "[[-1.88737914e-15+0.j]]\n",
      "[[0.+0.j]]\n",
      "[[5.55111512e-15+0.j]]\n"
     ]
    }
   ],
   "source": [
    "for i in range(0,15):\n",
    "    print(opt_result.states[0][i]-psi[nt-1][i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8\n",
      "7\n",
      "6\n",
      "5\n",
      "4\n",
      "3\n",
      "2\n",
      "1\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "tlist=np.linspace(0,1,10)\n",
    "for time_index in range(len(tlist) - 2, -1, -1):\n",
    "    print(time_index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "colab": {
   "collapsed_sections": [],
   "name": "Non-dissipative, Time=pi/g, g=0.2*omega.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}