FigSM7 - decays.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0b3366ef-660a-40a7-a966-2e77584e4253",
   "metadata": {},
   "source": [
    "# Code for Fig. SM 7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5745e2d6-816a-4237-9036-dd9cd9370e83",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-03T13:14:02.181477Z",
     "iopub.status.busy": "2022-02-03T13:14:02.181477Z",
     "iopub.status.idle": "2022-02-03T13:14:03.326549Z",
     "shell.execute_reply": "2022-02-03T13:14:03.326549Z",
     "shell.execute_reply.started": "2022-02-03T13:14:02.181477Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 1\n",
    "%aimport general_functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "36460638-b3ad-443a-8da5-d0f584069031",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-03T13:14:03.328549Z",
     "iopub.status.busy": "2022-02-03T13:14:03.327549Z",
     "iopub.status.idle": "2022-02-03T13:14:03.405591Z",
     "shell.execute_reply": "2022-02-03T13:14:03.405591Z",
     "shell.execute_reply.started": "2022-02-03T13:14:03.328549Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from tqdm.auto import tqdm\n",
    "from joblib import Parallel, delayed\n",
    "from scipy.interpolate import CubicSpline\n",
    "import copy\n",
    "\n",
    "from general_functions import (complete_dynamics, compute_shift)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b1b8bd4-5627-47ed-9607-5177c5ffdfad",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Single layer"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe74ef19-36b2-4f4b-9c3d-c6e5223e3e83",
   "metadata": {},
   "source": [
    "### Individual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "27d14f2d-57f1-4c71-9e92-d9eee80464b4",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-03T13:08:50.337726Z",
     "iopub.status.busy": "2022-02-03T13:08:50.337726Z",
     "iopub.status.idle": "2022-02-03T13:08:50.352715Z",
     "shell.execute_reply": "2022-02-03T13:08:50.352715Z",
     "shell.execute_reply.started": "2022-02-03T13:08:50.337726Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "pol_lat = np.array([[0, 1, 0]]).T\n",
    "pol_emi = [np.array([[1, 1, 0]] / np.sqrt(3)).T]\n",
    "\n",
    "d = 1 # Wave-lenght of transition\n",
    "lambda0 = d / 0.3 # Wave-lenght of transition\n",
    "k0 = 2 * np.pi / lambda0\n",
    "\n",
    "border = 10\n",
    "max_gamma = 10\n",
    "\n",
    "N_x = 60\n",
    "N_y = N_x\n",
    "x_pos = np.arange(0, N_x) - N_x / 2 + d / 2\n",
    "y_pos = np.arange(0, N_y) - N_y / 2 + d / 2\n",
    "x_pos, y_pos = np.meshgrid(x_pos, y_pos)\n",
    "x_pos = x_pos.flatten()\n",
    "y_pos = y_pos.flatten()\n",
    "z_pos = np.zeros(len(x_pos))\n",
    "r_pos = [x_pos, y_pos, z_pos]\n",
    "\n",
    "gamma_emi = 0.002\n",
    "omega_emi = -0.999  # For d / lambda_0 = 0.3\n",
    "\n",
    "dt = 2e-3 / gamma_emi\n",
    "tf = dt * 2 ** 10\n",
    "\n",
    "psi0 = np.zeros(len(x_pos) + 1)\n",
    "psi0[-1] = 1\n",
    "\n",
    "emitter = [[0, 0, 0]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "116589d9-160d-4fa7-8cfa-d41707301ec3",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-03T13:08:51.040576Z",
     "iopub.status.busy": "2022-02-03T13:08:51.039574Z",
     "iopub.status.idle": "2022-02-03T13:09:21.199640Z",
     "shell.execute_reply": "2022-02-03T13:09:21.199640Z",
     "shell.execute_reply.started": "2022-02-03T13:08:51.040576Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "results = complete_dynamics([x_pos, y_pos, z_pos], k0, pol_lat, dt, tf, psi0, \n",
    "                            emitter= emitter, pol_emi=pol_emi, gamma_emi=gamma_emi,\n",
    "                            omega_emi=omega_emi, border=border, max_gamma=max_gamma)\n",
    "\n",
    "population_ind = np.abs(results['psi'][:, -1]) ** 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24333b05-7856-4f42-bc28-029c8d11a556",
   "metadata": {},
   "source": [
    "### Collective"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "4ce1c301-1200-4b8d-8ba1-0a8a1f29377f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-03T13:09:21.200642Z",
     "iopub.status.busy": "2022-02-03T13:09:21.200642Z",
     "iopub.status.idle": "2022-02-03T13:09:21.215640Z",
     "shell.execute_reply": "2022-02-03T13:09:21.215640Z",
     "shell.execute_reply.started": "2022-02-03T13:09:21.200642Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "de = 14.14 * d\n",
    "emitter = [[-de / np.sqrt(2) / 2, de / np.sqrt(2) / 2, 0],\n",
    "           [de / np.sqrt(2) / 2, -de / np.sqrt(2) / 2, 0]]\n",
    "\n",
    "n_e = 2\n",
    "signs = [1, 1]\n",
    "psi0 = np.zeros(N_x * N_y + n_e, dtype=complex)\n",
    "for i in range(n_e):\n",
    "    psi0[-i - 1] = signs[i] / np.sqrt(n_e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "4c58bad9-3df3-48a3-9d90-c658349e34e2",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-03T13:09:21.216642Z",
     "iopub.status.busy": "2022-02-03T13:09:21.216642Z",
     "iopub.status.idle": "2022-02-03T13:09:52.400729Z",
     "shell.execute_reply": "2022-02-03T13:09:52.400729Z",
     "shell.execute_reply.started": "2022-02-03T13:09:21.216642Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "results = complete_dynamics([x_pos, y_pos, z_pos], k0, pol_lat, dt, tf, psi0, \n",
    "                            emitter=emitter, pol_emi=pol_emi * 2, gamma_emi=gamma_emi,\n",
    "                            omega_emi=omega_emi, border=border, max_gamma=max_gamma)\n",
    "\n",
    "population_coll = np.sum(np.abs(results['psi'][:, -n_e:]) ** 2, axis=-1)\n",
    "\n",
    "population = np.abs(results['psi'][:, :-n_e].reshape((-1, N_y, N_x))) ** 2\n",
    "index_max = np.argmax(np.sum(population, axis=(1, 2)))\n",
    "real_space_emission = population[index_max]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "e9428a69-4385-492d-ac1e-8e0d8dacafd9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-03T13:09:52.402733Z",
     "iopub.status.busy": "2022-02-03T13:09:52.401731Z",
     "iopub.status.idle": "2022-02-03T13:09:52.508730Z",
     "shell.execute_reply": "2022-02-03T13:09:52.508730Z",
     "shell.execute_reply.started": "2022-02-03T13:09:52.402733Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'population')"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEOCAYAAABmVAtTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAzgUlEQVR4nO3dd3RVRdfA4d9OQ3oNLUSKFGlqICAiWFAQFAwifFLEgpEX6aIiiMoLFlCkChaQ4gtKF0GkiYABRXroBiFIEaQkQKSXzPfHJCGEkNzAbUn2s9ZZyal3e9eRnTNzZo8YY1BKKaUAfDwdgFJKKe+hSUEppVQSTQpKKaWSaFJQSimVRJOCUkqpJJoUlFJKJXFZUhCRCSJyVES23WC/iMgoEdktIltEpIarYlFKKeUYVz4pTAIap7G/CVAhYekIfO7CWJRSSjnAZUnBGBMBxKZxSBjwP2P9DhQQkRKuikcppVT6PNmnEAQcSLZ+MGGbUkopD/Hz4GdLKttSrbkhIh2xTUwUhpqBBQqQ+447XBmbUkplKRs2bDhujAlM7zhPJoWDQHCy9VLAodQONMaMBcYCVPb1Nb/8+y+FIyLwzZXL9VEqpVQWICL7HDnOk81H84DnEt5CqgOcMsYcTu+k+KJFKXrlCptff931ESqlVDbjyldSpwKrgUoiclBEXhKRTiLSKeGQBUA0sBsYB3R25Lp5S5ZkR0AAeSdNwsTHuyR2pZTKrlzWfGSMaZPOfgN0yfCFRfinbVuiJ03i8JIlPNA4rbdelVJKZUSmHNFc9/PPebtYMT4cMcLToSilVJaSKZPCbbfdRo8ePYhdvJhty5Z5OhyllMoyMmVSAOjyxBOsBXb06OHpUJRSKsvItEkh3113satsWepu28Yf21Itr6SUUiqDMm1SACg2YAClgMVdMt5frZRS6nqZOinkb9uWE3nzUi0igt27d3s6HKWUyvQydVLA1xe/rl2pAYx8911PR6OUUple5k4KQN633uKDV17hi5kz2bt3r6fDUUqpTC3TJwXy5OHVfv3wEeHjDz/0dDRKKZWpZf6kAATlyMG+3LmRCRM4cOBA+icopZRKVZZIChQuTMGgIF6Jj+ejwYM9HY1SSmVaWSMpiJCjZ0+qAzvHjePw4XSLrSqllEpF1kgKAG3bciVfPjpevsyQIUM8HY1SSmVKWScp5MqF70sv8TTw/eefc/ToUU9HpJRSmU7WSQoAXbtydNgw/r5wgaFDh3o6GqWUynSyVlIoV46SPXvydOvWjBkzhuPHj3s6IqWUylSyVlIAOHeO4UWKUO/MGUbofAtKKZUhWS8pBARQbO5chgQGMmrUKE6cOOHpiJRSKtPIeknB1xc6daL6sWME/fsvw4YN83RESimVaWS9pADw0ksQEMCwO+5gxIgRHDt2zNMRKaVUppA1k0LRotCqFY3++Qc5c4bBOspZKaUckjWTAkDXrvg+/DAvt2jBmDFjOHjwoKcjUkopr5d1k0KdOvDDD3QdMoT4+Hjef/99T0eklFJeL+smhQRlgbdatWL8+PFER0d7OhyllPJqWTspXLoEtWrx5unT+Pn58d///tfTESmllFfL2knB3x+ee46cCxbQ98UXmTJlCtu3b/d0VEop5bWydlIAeOUVuHyZ1/PkIU+ePLyrczkrpdQNZf2kUKECPPEEub7+mt7du/Pdd9+xZs0aT0ellFJeKesnBYBXX4WTJ+lVrx5Fixald+/eGGM8HZVSSnmd7JEUGjSAAwfI1bgx/fv3JyIigh9//NHTUSmllNeRzPYXc2hoqFm/fv1Nn38pLo6qoaH4+/uzefNm/Pz8nBidUkp5JxHZYIwJTe+47PGkkKhVK/yffZZBgwaxY8cOvv76a09HpJRSXiV7JYXKlWH+fFpUr06dOnV49913OXv2rKejUkopr5G9ksIrr4CfHzJ6NB9//DGHDh3SiXiUUiqZ7JUUSpSANm1gwgTqV6/Ok08+yeDBgzly5IinI1NKKa/g0qQgIo1FJEpEdotIn1T25xeRH0Rks4hsF5EXXRkPAD17wpkzMHEiH3/8MefOndMBbUoplcBlSUFEfIExQBOgCtBGRKqkOKwLsMMYczfwEDBURAJcFRMAISEwbRqEh1OpUiW6du3KV199xebNm136sUoplRm48kmhNrDbGBNtjLkITAPCUhxjgLwiIkAeIBa47MKYrGeegbx5AXj33XcpUKAAr776qg5oU0ple65MCkHAgWTrBxO2JTcaqAwcArYCPYwx8SkvJCIdRWS9iKx32tSa8+ZBy5YULFCAgQMHsnz5cubNm+ecayulVCblyqQgqWxL+af4Y0AkUBK4BxgtIvmuO8mYscaYUGNMaGBgoHOiO3kSZs+GRYv4z3/+Q+XKlXnttde4cOGCc66vlFKZkCuTwkEgONl6KewTQXIvAt8ZazewF7jThTFd1bo1BAXBxx/j5+fHsGHD2LNnD6NHj3bLxyullDdyZVJYB1QQkbIJncetgZTtM/uBRwBEpBhQCXDP9GgBAbZQ3ooVsG4djRs3pkmTJgwcOBCnNVEppVQm47KkYIy5DHQFFgM7gRnGmO0i0klEOiUc9h5QV0S2Aj8Dbxpjjrsqpuu8/DLkzw9DhgAwdOhQzpw5o6+oKqWyrWxXEO86n38OBQva5iSge/fujBkzhsjISKpXr+68z1FKKQ9ytCCeJoUUYmNjqVChAtWqVWPFihXYt2WVUipz0yqpGREXB4MGwcGDFCpUiEGDBhEREcE333zj6ciUUsqtNCkAxMbCO+/AJ58AEB4eTu3atXn99dc5efKkZ2NTSik30qQAUKYMtG8PY8fCkSP4+Pjw2WefcfToUfr37+/p6JRSym00KSTq2xfOn4fhwwGoWbMmnTp1YvTo0URGRno2NqWUchNNCokqVrQ1kcaMsc1JwAcffEDhwoXp3Lkz8fHXVd9QSqksR5NCcm+9BXXrwokTABQsWJCPP/6Y1atXM2nSJM/GppRSbqCvpKYjPj6eBx54gKioKKKioihUqJDbPlsppZxFX0m9Ffv3w5IlAEmdzidOnKBfv34eDkwppVxLk0JqunSBdu3sDG3AXXfdRbdu3fjyyy9Zu3ath4NTSinX0aSQmr594fhx2+mcYMCAAZQsWZLw8HAuXbrkweCUUsp1NCmkpm5daNwYPv4Y/v0XgHz58vHZZ5+xdetWhiQU0FNKqaxGk8KNDBwIMTEwcmTSpieffJJWrVoxcOBAoqKiPBicUkq5hiaFG6lVC1q0SBqzkGjUqFHkzJmTjh076tgFpVSWo0khLTNnwrBh12wqXrw4Q4cOJSIignHjxnkoMKWUcg1NCmnxSfh61q2zTUkJXnzxRRo0aEDv3r35+++/PRScUko5nyaF9OzfD3XqJFVQBRARvvzySy5evEiXLl3IbAMAlVLqRjQppOf2221NpFGj4J9/kjaXL1+eAQMGMHfuXL777jsPBqiUUs6jScERAwfCxYv2ZzK9evUiJCSELl26cPy4+6aWVkopV9Gk4Ijy5aFjRzvfwq5dSZv9/PyYOHEisbGxdO3a1YMBKqWUc2hScNS770KxYrB9+zWb7777bvr378/06dOZPn26h4JTSinn0CqpGXHxIgQEXLf58uXL3H///ezevZvt27dTvHhxDwSnlFI3plVSXSEgAIyBn3+2PxP4+fnx9ddfc/bsWV5++WV9G0kplWlpUsioqVPh0Udh0aJrNt95550MGjSI+fPn8/XXX3soOKWUujXafJRRFy9C5cqQOzds2gS+vkm74uPjadCgAZs2bWLr1q3cfvvtnotTKaWS0eYjVwkIgA8+gK1b4dtvr9nl4+PDxIkTuXLlCi+99JLWRlJKZTqaFG7G//0f1KgBb78N585ds6ts2bIMGzaMpUuX8vnnn3soQKWUujmaFG6Gj48te+HjA9HR1+1++eWXady4MW+88QY7duzwQIBKKXVzNCncrIcftgPZqla9bpeIMHHiRPLkyUObNm04f/68BwJUSqmM06RwK/z9bfPRDz9ct6t48eJMmjSJLVu20KdPHw8Ep5RSGadJ4VYNHQphYbBhw3W7Hn/8cbp3787IkSNZsGCBB4JTSqmM0VdSb1VcHFSoAJUqwS+/gMg1u8+fP8+9997L4cOH2bJli452Vkp5hL6S6i758sH778PKlTB79nW7b7vtNqZOncrp06d5/vnn9TVVpZRX06TgDB06wF13weuvw9mz1+2uUqUKw4cPZ8mSJYwYMcL98SmllINcmhREpLGIRInIbhFJtbdVRB4SkUgR2S4iv7gyHpfx9YVPP4WSJeEG8yp07NiR5s2b06dPHzak0v+glFLewOE+BRHxBYoBfonbjDH70zl+F9AQOAisA9oYY3YkO6YA8BvQ2BizX0SKGmOOphWH1/UpJGfMdX0KycXExBASEoKfnx8bN26kQIEC7otNKZWtObVPQUS6AUeAn4AfE5b56ZxWG9htjIk2xlwEpgFhKY5pC3yXmFzSSwheTwSOHoXhw6+popqocOHCzJgxg4MHD/LCCy9oNVWllNdxtPmoB1DJGFPVGFM9YbkrnXOCgAPJ1g8mbEuuIlBQRFaIyAYRec7BeLzX1KnQqxfMm5fq7jp16jBkyBDmzp3L0KFD3RycUkqlzdGkcAA4lcFrp9aOkvJPYz+gJvAE8BjwjohUvO5CIh1FZL2IrD927FgGw3Czzp2hWjXo0SPVTmeA7t278/TTT9OnTx9WrVrl5gCVUurGHE0K0cAKEekrIr0Sl3TOOQgEJ1svBRxK5ZhFxpgzxpjjQARwd8oLGWPGGmNCjTGhgYGBDobsIf7+MGYM7NsHgwaleoiIMH78eMqWLcszzzzD0aOZu9VMKZV1OJoU9mP7EwKAvMmWtKwDKohIWREJAFoDKdtU5gL1RcRPRHIB9wI7HQ3eaz3wADz7LHz8MURFpXpI/vz5mTVrFrGxsbRr144rV664OUillLqeX/qHgDFmAICI5LWr5rQD51wWka7AYsAXmGCM2S4inRL2f2GM2Skii4AtQDzwlTFm203+t3iXTz6xTw3589/wkLvvvpvRo0cTHh7OwIEDGTBggBsDVEqp6zn0SqqIVAMmA4USNh0HnjPGbHdhbKny6ldSb4Ixhg4dOjBp0iTmzp3Lk08+6emQlFJZkLPLXIwFehljShtjSgOvAeNuJcBsY/duaNYMDh9OdbeI8NlnnxEaGsqzzz7Lzp2Zv/VMKZV5OZoUchtjlieuGGNWALldElFWEx8PP/1k30a6gZw5c/Ldd9+RM2dOwsLCOHnypPviU0qpZBx++0hE3hGRMgnL28BeVwaWZVSsCO+8AzNnpjrvQqLg4GBmzZrF3r17teNZKeUxjiaFDkAg8B0wJ+H3F10VVJbzxht27ELnzvDvvzc8rH79+owaNYoFCxbQv39/NwaolFKWQ0nBGHPCGNPdGFPDGBNijOlhjDnh6uCyjIAAGDcO/v7bvpWUhk6dOvHyyy/zwQcfMGvWLDcFqJRSVppvH4nICGNMTxH5getHI2OMcfurMpn67aM5c+CxxyBXrjQPu3DhAg8//DCbN29m9erV3HVXehVFlFIqbY6+fZReUqhpjNkgIg+mtt8Y4/ZS15k6KSQ6e9YWzMt94776w4cPU6tWLXx8fFi7dq3O2KaUuiVOeSXVGJNY+P8eY8wvyRfgHifEmf2cOwc1ath+hjSUKFGCH374gZiYGMLCwjh7gzpKSinlTI52ND+fyrYXnBhH9pEzJzzxBHz+OSxdmuahISEhTJ06lXXr1ulUnkopt0gzKYhIm4T+hLIiMi/ZshyIcU+IWdD779tXVV96CeLi0jz0ySef5JNPPmHWrFm8/fbbbgpQKZVdpVf76DfgMFAESF78/19svSJ1M3LmhEmToF49eO01+2ZSGl599VV27drFoEGDqFixIi+88IJbwlRKZT9pJgVjzD5gH3Cfe8LJRu67D15/HVatgvPn4bbbbnioiPDpp58SHR1Nx44dKVOmDA899JD7YlVKZRuOTsdZR0TWichpEbkoIldEJO12D5W+996DX35JMyEk8vf3Z8aMGVSoUIHmzZuzdetWNwSolMpuHO1oHg20Af4EcgLhwKeuCirbCAgAPz84dsz2M6RTsbZAgQIsWLCA3Llz06RJE/bv3++mQJVS2YWjSQFjzG7A1xhzxRgzEXjYdWFlM99/b+sjjR6d7qGlS5dm0aJFnD59msaNGxMbG+v6+JRS2YajSeFswuxpkSLysYi8ilZJdZ7wcPua6htvwPb0p6ioXr06c+fOZc+ePTRr1oxz5865IUilVHbgaFJoj509rStwBjv38tOuCirbEYHx4yFfPmjb1nY8p+PBBx/km2++YfXq1bRu3ZrLly+7IVClVFbnaEG8fcaYc8aYOGPMAGNMr4TmJOUsxYrBxImwZQsMHOjQKS1btmTUqFHMmzePLl264MgsekoplZY0X0kVka2kUggvkTFGK7U50xNPwIQJ0LSpw6d07dqVQ4cOMWjQIAoXLsyHH37owgCVUlldeoPXHP/XSTnHiwnTVFy+DDEx9gkiHR988AExMTEMGjSIvHnz0rdvXxcHqZTKqhwZvKY8oUULOHQIfv0VcuRI89DEeZ7PnDnDW2+9RZ48eejWrZubAlVKZSWODl77V0TiEpbzOnjNDV56CTZsSLeaaiJfX18mTZpE8+bN6d69OxMmTHBxgEqprMjRjua8xph8Cctt2DeP0n+pXt28sDB49VX49FM7v7MD/Pz8mDZtGo0aNeLll19mxowZLg5SKZXVODx4LTljzPdAA+eGoq4zeDDUqWP7GRwYvwCQI0cO5syZw/3330+7du2YP3++i4NUSmUljjYftUi2tBSRwaTxVpJykoAAmDUL7r4bMjCXQq5cuZg/fz4hISG0bNmSRYsWuTBIpVRW4uiTQrNky2PY0tlhrgpKJRMUZCupVq9u1x0ci5AvXz4WLVpElSpVCAsLY8GCBS4MUimVVaT3SioAxpgXXR2ISoOIfUW1SxcoVcrWSXJAoUKFWLp0KY0aNeKpp55i9uzZNM3AGAilVPbjaPNRORH5QUSOichREZkrIuVcHZxKxtfXzu/cvz9k4K/+QoUK8dNPP3H33XfTokUL5s2b58IglVKZnaPNR98CM4ASQElgJjDVVUGpVIjAF1/Y/oV27WDXLodPLViwIEuWLEnqY/j+++9dF6dSKlNzNCmIMWayMeZywjIF7Wh2v1y5YM4c8Pe3pTBOnHD41AIFCrBkyRJq1qxJq1atmD17tgsDVUplVo4mheUi0kdEyohIaRHpDfwoIoVEpJArA1QplCljE0NsLGzblqFT8+fPz+LFi6lVqxbPPPMMU6ZMcU2MSqlMSxyprCkie9PYbYwxbutfCA0NNevXr3fXx3mv06chT56bOvXff/+lefPmLFu2jE8//ZSuXbs6OTillLcRkQ3GmND0jnP07aOytx6ScqrEhDBmjH1NNQP/sOfNm5cff/yR1q1b061bN06ePEm/fv0QERcFq5TKLBx9+8hfRLqLyKyEpauI+Ls6OJWO+HhYuhR69IAff8zQqbfddhuzZs2iffv2vPPOO7z++us6H4NSyrEnBeBzwB/4LGG9fcK2cFcEpRzk4wOTJ8NDD8H//R+sWAG1ajl8up+fH5MmTaJAgQIMGzaMkydPMnbsWHx9fV0WslLKuzna0VzLGPO8MWZZwvIikO6/PiLSWESiRGS3iPRJ47haCZVXWzoauEqQJ499SihWzE7Ss2dPhk738fFh5MiRvPvuu0yYMIGWLVvqnM9KZWOOJoUrInJH4krCwLUraZ0gIr7AGKAJUAVoIyJVbnDcR8BiR4NWKRQrBgsX2r6F5cszfLqIMGDAAEaNGsXcuXN55JFHOH78uAsCVUp5O0ebj97AvpYanbBeBkiv9EVtYLcxJhpARKZh6yXtSHFcN2A2Djx5qDRUqmQHtBUseNOX6NatG0FBQbRr1466deuycOFC7rjjjvRPVEplGY4+KfwKfAnEJyxfAqvTOScIOJBs/WDCtiQiEgQ8BXyR1oVEpKOIrBeR9ceOHXMw5GwoMSEsXw7t28OlSxm+RIsWLfj555+JiYnhvvvuY926dU4OUinlzRxNCv8DygLvJSxlgcnpnJPa+40pX28ZAbxpjEmzKcoYM9YYE2qMCQ0MDHQs4uzsjz9gyhTo0CFDJbcT1a1bl99++408efLw0EMP6ZwMSmUjjjYfVTLG3J1sfbmIbE7nnINAcLL1UsChFMeEAtMS3o8vAjwuIpcTJvFRN+uVV2wJjH79IH9+O3tbBscgVKpUidWrV9O0aVPCwsIYOXIkXbp00bEMSmVxjj4pbBKROokrInIvtkkpLeuACiJSVkQCgNbANSU6jTFljTFljDFlgFlAZ00ITtK3r53fecwYh0ttp1SsWDFWrFhB06ZN6datG507d+bSTTRJKaUyD0efFO4FnhOR/QnrtwM7RWQrtszFXSlPMMZcFpGu2LeKfIEJxpjtItIpYX+a/QjqFonARx/ByZPw11+2Gckn47Ov5s6dmzlz5tCvXz8GDx5MVFQUM2fOpHDhwk4PWSnleY7WPiqd1n5jzD6nRZQOrX2UQVeu2ATh42PnY8iZ86YvNXnyZMLDwwkODuaHH36gcuXKTgxUKeVKjtY+cuhPR2PMvrSWWw9XuYyvr00IR47YuRg++yz9c26gffv2/PLLL5w+fZo6deqwcOFCJwaqlPIGGW9PUJlTwYJw5512Ss9bSAx16tRh3bp1lCtXjqZNm/LRRx9pzSSlshBNCtlFQADMmgXNmt1yYggODmbVqlW0atWKPn368PTTTxMXF+fEYJVSnqJJITtJmRim3vyMqrlz52bq1KkMHz6cefPmUatWLbZv3+7EYJVSnqBJIbtJTAyvvw4NG97SpUSEnj17snz5cuLi4qhduzbTpk1zUqBKKU/QpJAdBQTAkCFQpAhcuAC3+A95/fr12bhxIzVq1KBNmzb07NmTixcvOilYpZQ7aVLI7saNgzZtoE8fW2X1JpUoUYJly5bRs2dPRo4cSf369YmOjk7/RKWUV9GkkN117gydOtmBbp0731StpET+/v4MHz6cWbNmERUVRUhICNOnT3disEopV9OkkN35+Ng3kd58E7744qarqyb39NNPExkZSZUqVWjdujUdO3bk7NmzTgpYKeVKmhSUHfE8eDB8+CEsWAD7bn08YpkyZYiIiKBPnz589dVX1KpVi23btjkhWKWUK2lSUFf17WvLbpcvb/sXTp26pcv5+/szaNAgFi9ezPHjx6lVqxajRo0i/haaqJRSrqVJQV2rWDH7c9gwWxbjjz9u+ZINGzZk8+bNNGjQgB49evDYY49x8ODBW76uUsr5NCmo1D34oC2gV7curFp1y5crXrw48+fP58svv2T16tVUr16dqbcweE4p5RqaFFTqQkNh9WoIDIRHH4VvvrnlS4oIHTt2JDIyksqVK9O2bVtat25NbGysEwJWSjmDJgV1Y+XKwW+/wb33wvPPw+7dTrls+fLliYiI4IMPPmD27NlUq1aNuXPnOuXaSqlbo0lBpa1wYfjpJ1i82HZAg52j4Rb5+fnx1ltvsXbtWgIDA2nevDmtW7fm6NGjt3xtpdTN06Sg0hcQAI88Yn9fsMA2Lf31l1MuHRISwvr163nvvfeYM2cOVapU4ZtvvtFy3Ep5iCYFlTH+/rB3L9SqBStXOumS/rz99tts2rSJChUq8Oyzz9KsWTN9Q0kpD9CkoDKmYUNYswYKFYKHH7avrjrpr/oqVaqwatUqhg8fzvLly6lSpQqjRo3i8uXLTrm+Uip9mhRUxlWqBGvXwpNPwmuvgROn5fT19aVnz55s3bqV++67jx49elCrVi1+//13p32GUurGNCmom5M/P8yeDfPmQZMmdtu5c067fLly5Vi0aBEzZ87k2LFj3HfffXTs2JGYmBinfYZS6nqaFNTNE7GzuInArl32FdYpU5x4eaFly5bs3LmT1157jQkTJlCpUiXGjx+vpTKUchFNCso58uWzzUrt20OHDnD6tNMunTdvXj755BM2bdpE5cqVCQ8Pp169eqxbt85pn6GUsjQpKOcoXhyWLoW334avv4aQEHDyP9rVq1cnIiKCiRMnsmfPHmrXrs3zzz/P33//7dTPUSo706SgnMfPD957D5Yvh/PnYdIkp3+EiPDCCy/w559/0qdPH6ZNm0bFihUZOHCgztmglBNIZhskFBoaatavX+/pMFR6TpyAHDkgVy7Ytg0KFIBSpZz+MXv37uXNN99k5syZlCpVisGDB9OmTRt8fPTvHaWSE5ENxpjQ9I7T/3OUaxQsaBOCMfDcc1C9Okye7LQxDYnKli3LjBkziIiIoFixYjz77LPUrl2bpUuXOvVzlMouNCko1xKB6dOhShWbHJo3h3/+cfrH1K9fn7Vr1/L1119z/PhxGjZsyKOPPqqd0UplkCYF5XoVKkBEBAwdCkuW2AThgqk5fXx8eO6554iKimLEiBFs2bKF2rVr07JlS/5wwmRBSmUHmhSUe/j6Qq9eEBkJzzwDd95pt7ughEWOHDno0aMHe/bs4b///S+LFy+matWqhIeH85eTCvkplVVpUlDuVakSfP65fVMpJsauf/qpU8pxp5Q3b1769+9PdHQ03bt3Z/LkyVSoUIHw8HCio6Od/nlKZQWaFJTnnD9v52jo3h3q1IFNm1zyMYGBgQwfPpw9e/bQqVMnpkyZQsWKFenQoQN79uxxyWcqlVlpUlCeExQEixbB1Klw4ICdp6FXL5c0KQGUKlWKTz/9lOjoaLp06cLUqVOpVKlS0rgHpZQmBeVpItC6NfzxB3TsCPv22aYlFypZsiQjR44kOjqabt26MX36dO68807atGnDxo0bXfrZSnk7lyYFEWksIlEisltE+qSyv52IbElYfhORu10Zj/JiBQrYvoYZM+z6rl3QoAFs2OCyjyxRogTDhw9n79699OrVix9//JGaNWvSsGFDlixZorO/qWzJZUlBRHyBMUAToArQRkSqpDhsL/CgMeYu4D1grKviUZmEr6/9+ddfsH27neEtPBxcOHdz8eLFGTJkCAcOHOCjjz5i+/btPPbYY4SEhPDtt9/qJD8qW3Hlk0JtYLcxJtoYcxGYBoQlP8AY85sx5kTC6u+A8+sgqMypUSP7tNCrly2wV6ECDB/u0o/Mnz8/vXv3Zu/evYwfP56LFy/Srl07ypcvz4gRIzh16pRLP18pb+DKpBAEHEi2fjBh2428BDhvCi+V+eXPD598Yge6PfAA7N9/dZ8LXmFNlCNHDjp06MC2bduYN28ewcHBvPrqqwQFBdGlSxd27Njhss9WytNcmRQklW2pNtKKyMPYpPDmDfZ3FJH1IrL+2LFjTgxRZQqVKsEPP9gEAbBsGVSrBnPmOL2WUnI+Pj40a9aMlStXsn79elq1asX48eOpWrUqjz76KHPnzuWKC5OTUp7gyqRwEAhOtl4KOJTyIBG5C/gKCDPGpDrXojFmrDEm1BgTGhgY6JJgVSaQ2N8gCX9vtGgBdevCihUu/+iaNWsyceJEDhw4wIcffkhUVBTNmzenfPnyDBkyRKcJVVmGK5PCOqCCiJQVkQCgNTAv+QEicjvwHdDeGLPLhbGorOThh2HrVhg3zjYpPfywLZ3hBoGBgfTt25e9e/cya9YsSpcuTe/evQkKCqJt27YsX75c31pSmZpL51MQkceBEYAvMMEY84GIdAIwxnwhIl8BTwP7Ek65nF69b51PQV3j3DmbHHLnhpdesgPffv3V9kFIai2Yzrd161bGjh3LlClTOHnyJOXLlyc8PJwXXniBYsWKuSUGpdLj6HwKOsmOylq++QaefRbq1bNTgzZq5LbkcO7cOWbNmsW4ceNYuXIlfn5+PPnkk4SHh9OwYUP8XDwoT6m06CQ7Knt6+mkYPRr27oXGje1c0d9+69K3lRLlzJmT9u3bExERwc6dO+nRowcRERE8/vjjlCpVitdee43IyEhtXlJeTZ8UVNZ04YJ9ahgyxJbN2LLFPjHEx4Mbp+q8cOECCxYsYPLkycyfP59Lly5RvXp12rdvT7t27ShZsqTbYlHZmzYfKQU2CRw+bIvvxcXBPffYWkuvvALBweme7kwxMTHMmDGD//3vf/z+++/4+PjwyCOP0LZtW5o3b06BAgXcGo/KXrT5SCmwTwVBCWMmT52ySeGjj6BsWWjVClaudOlYh+QKFy7MK6+8wurVq4mKiqJfv37s3r2bF198kaJFi9KsWTOmTJlCXFycW+JRKjX6pKCyn7/+gs8+g6++ghMn7IjpqlU9EooxhnXr1jFjxgxmzJjBgQMHyJEjB02aNOGZZ56hadOm5MmTxyOxqaxFm4+USs/Zs7B4MTz1lF3v1s3+DA+Hu91fsDc+Pp41a9Ywffp0Zs6cyaFDh7jtttto1KgRzZs3p2nTpujgTXWzNCkolRHGwMsvw5QptpM6NNQmhzZtIF8+t4cTHx/Pr7/+ysyZM/n+++85cOAAPj4+3H///TRv3pywsDDuuOMOt8elMi9NCkrdjNhY+9bSuHF21PSbb8LgwbbDWsRtYx6SM8YQGRnJ999/z/fff8+WLVsAqFatGmFhYYSFhVGzZk183PhWlcp8NCkodSuMgfXroXhx+5bSwoXQtSu0bQvt2sGdd3ostL179zJ37lzmzp1LREQE8fHxFC1alMcee4zGjRvTqFEjihQp4rH4lHfSpKCUM0VEwHvv2Qqt8fFQo4ZNDp07w223eSysmJgYFi5cyMKFC1m8eDExMTGICLVr16ZJkyY0btyY0NBQfBOLCapsS5OCUq5w+DBMm2abmP75xxbk8/GB1avtREAe/Av9ypUrbNiwgYULF7Jo0SLWrFmDMYbChQvTqFEjHn30UR555BFKly7tsRiV52hSUMrVTpyAggXtk0OpUnDkCDz4oC218dRT4OHRyjExMfz0009JTxFHjhwBoFy5cjRo0CBp0aJ92YMmBaXcxRiIjITZs+3yxx92+4cfQt++Hg0tkTGGHTt2sGzZMpYtW8aKFSs4efIkAFWrVk1KEPXq1dP+iCxKk4JSnrJzp00ODRrYSYDWr4fnn4emTeGJJ+w2D1dMvXLlCps2bWLZsmX8/PPPrFy5knPnzgFQuXJl6tWrR/369alXrx5lypRBPPDWlXIuTQpKeYvVq+Gdd2xn9aVLtsnpscdg2DAoUcLT0QG2cN+6detYtWoVK1eu5Ndff+XUqVMABAUFJSWIevXqUa1aNe24zoQ0KSjlbeLi4Kef4McfYfly2LEDcua0YyL27YNHH4X77oMcOTwdKVeuXGH79u2sXLkyaTl0yM6mmzt3bkJDQ7n33nuTlqDE+lLKa2lSUMqbGXN1INx//gPjx9s5H3LlsrPGhYVBp06ejTEZYwx//fUXq1atYu3ataxZs4bIyEguXboE2KeJ5EkiNDSU3LlzezhqlZwmBaUyk7g4WLHCPkksXWpfb52XMKV5jx5QvrxNFtWqgZc03Zw/f57IyEjWrFmTtERHRwPg4+NDlSpVCAkJISQkhBo1anDPPfeQP39+D0edfWlSUCozO3/eDoo7e9ZWcP3rL7s9f3471egrr9hOay9z7Ngx1q1bx5o1a9iwYQMbN27k8OHDSfvLlStHjRo1kpJFSEgIxYsX92DE2YejSUEnjVXKGyWOks6Vy04tun+/nfshIsIuif/Q/vmnLeR3771Qp4796cHxEYGBgTz++OM8/vjjSduOHDnCpk2b2LhxI5s2bWLTpk3MmjUraX/x4sW55557qFatWtJSuXJlcuXK5Yn/hGxPnxSUyowS+yTWrrUlvzdtsm82ga3VNGsW1K4N//5rX3/NmdOz8aZw6tQpIiMjk5LF1q1b2blzJxcuXABARLjjjjuuSRTVqlWjYsWK+Pv7ezj6zEmbj5TKTs6ftwPo1qyB33+H4cNtMb+hQ22l1ypVICTE1myqUcO+5eThsRIpXb58mT179rBt27Zrlj///JMrV64A4O/vT8WKFalcuTKVKlVKWipWrKjTmaZDk4JSCtatg7lz7ZPEhg22FIevL5w+bZuovvkG/v7bTlNavbpNJF42UO3ChQtERUVdkyj++OMPoqOjk5IFQNGiRa9LFJUqVaJcuXL6dIEmBaVUag4fhl27bI0msJMITZt2dX+hQnYk9syZdj0qCooVAy/8K/zixYtER0cTFRXFrl27iIqKSlqOHTuWdJyfnx/lypXjjjvuoFy5ctct2WW6U00KSinHxMTYCYW2bbM/c+WyzU9gm5127rQF/6pVs/NIPPggNG/u0ZDTc+LEiWsSxa5du4iOjmbPnj1JI7UTFS1aNNVkUa5cOUqWLJllRm9rUlBK3brFi21fxbZtdomKgv/7P5g0yXZ2BwfbJqdKla4utWpBuXKejvyGYmNjiY6OTnXZt28f8fHxScf6+fkRFBTE7bffTnBwMLfffnvSkrieP3/+TFEbSpOCUsr54uPt2Ik8eWzndu/eNlH88Yd9bRagXz94/31bWrxZM5sgki9Vq9r6T17o0qVL7N+/PylJ7N+/n/3793PgwIGkn5cvX77mnLx5816TMIKDgwkODqZEiRKULFmSEiVKUKhQIY8nDk0KSin3OnvWjpsoUABKl7b1nF54AaKj4cAB+2QB8OWX0LGjTSavvQZly8Ltt9slONg2U+XL58n/khu6cuUKR44cuSZJpEwcR48eve68gICAa5LEjX4WLlzYZclDk4JSyntcuGCfJKKjbT9FcLAdY/Gf/9htcXFXj12wAJo0gZ9/hv79ryaLxJ8PPmhHdnupc+fOcfjwYQ4dOpT0M7XfE+ezSC4gIIDixYtTrFgxihYtmrQEBgZes564LSAgwOG4dESzUsp75Mhh6zlVqHB1W+3a9lVZgFOn7NPEgQN2O9gnC39/O/Zi9my4eNFu37LFvj775ZcwYIAdwV2ihF1KloSePe3TyrFj9pxixdw6JiNnzpxJHdVpSZ48kieNw4cPc/ToUQ4dOkRkZCRHjx5NKjyYUoECBa5LFkWKFKFw4cIUKVLkmsVR+qSglPJ+8fFw9KhNGtWr2zEWS5fC1Klw6JB91fbwYZsIYmNtUujbFwYPtuMuCheGwEC7LF1qk83ixbB799XtRYte/elFjDHExcVx9OhRh5bY2NhrOsuT0eYjpVQ2c+mSfSoQgY0bbRNVYrI4etQ+kfz0kz32uedg8uRrz8+fHxKbdd54w55fqNDVpXRp6NrV7t+82ZY7L1jQ7suXzysG/sXHx3Py5EmOHz9OTEwMx48f5/jx43To0EGTglJK3dDly3aMRmLCOHbM9n0895zd/9//2smQTpywTx+xsbZTfPt2u79+fVi16ur1fH1tefNly+x69+52BHm+fDbZ5Mtn+1NatrT7V6+2CSz5/pw5XZZYtKNZKaWc7dIl2/QEtmzIwYM2WSQmjsBAO/8FQIsWNoHExdnl7Flb7nz+fLs/KMg2fSX3zDNXR5jXr28TRL58kDu3fQ34kUegbVu7f9QoO9AwT56r+8uWhTJlbH9MbKzdFhAAItrRrJRSTpe8hlLNmna5ke++u3b90qWrneVgS4mcOGGbtOLi7M/kHfGBgXb/4cNw5oytV1WkiE0KFy9eTT7J9e4NH31km8ASO5d9fW1ycJBLk4KINAZGAr7AV8aYwSn2S8L+x4GzwAvGmI2ujEkppTzC3//apFK3btrHp0wqKa8VE2MTRWLCOHPGPn2AfToYOfLqvvPnYdgwh8J0WfORiPgCu4CGwEFgHdDGGLMj2TGPA92wSeFeYKQx5t60rqvNR0oplXGONh/5uDCG2sBuY0y0MeYiMA0IS3FMGPA/Y/0OFBCREi6MSSmlVBpc2XwUBBxItn4Q+zSQ3jFBwOHkB4lIR6BjwuoFEdnm3FCdrghw3NNBpENjdA6N8dZ5e3yQNWIs7chFXJkUUnuvKmVblSPHYIwZC4wFEJH1jjwCeZLG6Bwao3N4e4zeHh9krxhd2Xx0EAhOtl4KOHQTxyillHITVyaFdUAFESkrIgFAa2BeimPmAc+JVQc4ZYw5nPJCSiml3MNlzUfGmMsi0hVYjH0ldYIxZruIdErY/wWwAPvm0W7sK6kvOnDpsS4K2Zk0RufQGJ3D22P09vggG8WY6UY0K6WUch1XNh8ppZTKZDQpKKWUSuJVSUFEGotIlIjsFpE+qewXERmVsH+LiNRw9Fw3xtguIbYtIvKbiNydbN9fIrJVRCJFxCXDsh2I7yEROZUQQ6SIvOvouW6M8Y1k8W0TkSsiUihhn8u/w4TPmSAiR280JsbT96ID8Xn0PnQwRm+4F9OL0RvuxWARWS4iO0Vku4hcV/TIqfejMcYrFmxn9B6gHBAAbAaqpDjmcWAhdnxDHWCNo+e6Mca6QMGE35skxpiw/hdQxMPf4UPA/Js5110xpji+GbDMXd9hss95AKgBbLvBfk/fi+nF57H7MAMxevRedCRGL7kXSwA1En7Piy0f5LJ/G73pSeFWymI4cq5bYjTG/GaMOZGw+jt27IW73Mr34DXfYQptgKkuiCNNxpgIIDaNQzx6L6YXn4fvw8QY0vsOb8Rd92JGY/TUvXjYJBQKNcb8C+zEVn5Izmn3ozclhRuVvHDkGEfOdVeMyb2Ezd6JDLBERDaILd3hqfjuE5HNIrJQRKpm8Fx3xYiI5AIaA7OTbXb1d+goT9+LGeHu+zAjPHkvOsxb7kURKQOEAGtS7HLa/ehN8yncSlkMh8plOIHDnyMiD2P/Z6yXbPP9xphDIlIU+ElE/kj4S8Wd8W0EShtjToutUvs9UMHBc50hI5/TDPjVGJP8LzlXf4eO8vS96BAP3YeO8vS9mBEevxdFJA82KfU0xsSl3J3KKTd1P3rTk8KtlMVwV7kMhz5HRO4CvgLCjDExiduNMYcSfh4F5mAf7dwanzEmzhhzOuH3BYC/iBRx5Fx3xZhMa1I8rrvhO3SUp+/FdHnwPnSIF9yLGeHRe1FE/LEJ4RtjTGoTLTjvfnR1J0kGOlP8gGigLFc7RKqmOOYJru1MWevouW6M8XbsCO26KbbnBvIm+/03oLEH4ivO1UGLtYH9Cd+n13yHCcflx7b15nbnd5gihjLcuJPUo/eiA/F57D7MQIwevRcdidEb7sWE7+R/wIg0jnHa/eg1zUfmFspi3OhcD8X4LlAY+EzsBNyXja1cWAyYk7DND/jWGLPIA/G1BF4RkcvAOaC1sXePN32HAE8BS4wxZ5Kd7vLvMJGITMW+HVNERA4C/QH/ZDF69F50ID6P3YcZiNGj96KDMYKH70XgfqA9sFVEIhO2vYVN/E6/H7XMhVJKqSTe1KeglFLKwzQpKKWUSqJJQSmlVBJNCkoppZJoUlBKKZVEk4JSSqkkmhSUugkiUkpEnkn4/Uqy8sqRiWMulMqMvGbwmlKZzCNAFWA6cM4Yc49nw1HKOXTwmlIZJCL1gLnASeBfoIIxJrdHg1LKSbT5SKkMMsasAtZhC83dg+creCrlNJoUlLo5lYAoTwehlLNpn4JSGSQihYFTxphL6RzXHygEnDTG9HdLcErdIn1SUCrjypJOTXoRCcJW2zyJLWWsVKagTwpKZdwf2FLL24AbTcP4HtADCOTaSU6U8mr69pFSLiAirwH5sHMaRBljPvVwSEo5RJOCUkqpJNqnoJRSKokmBaWUUkk0KSillEqiSUEppVQSTQpKKaWSaFJQSimVRJOCUkqpJJoUlFJKJdGkoJRSKsn/AyzHGEzPTJ76AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(results['time'] * gamma_emi, population_ind, 'k', label='individual')\n",
    "ax.plot(results['time'] * gamma_emi, population_coll, 'r--', label='collective')\n",
    "ax.set_ylim(0, 1)\n",
    "ax.set_xlim(0, tf * gamma_emi)\n",
    "ax.set_xlabel(r'$t\\Gamma_a$')\n",
    "ax.set_ylabel('population');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "945456ac-7b8f-48ff-b550-811c316c9332",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-03T13:09:52.509731Z",
     "iopub.status.busy": "2022-02-03T13:09:52.509731Z",
     "iopub.status.idle": "2022-02-03T13:09:52.632731Z",
     "shell.execute_reply": "2022-02-03T13:09:52.632731Z",
     "shell.execute_reply.started": "2022-02-03T13:09:52.509731Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAREAAAEGCAYAAABCR6GtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAYIElEQVR4nO3dfbTVVZkH8O9XMd/QDFDeR2ykpYA2TQ7RaGWJiUWipqm9DEwaq8aZsKkmSNdIM7pgjYy9LFOHzMEVJLIqEa1UsEixQMlRvHYhrghxE3DISBnzBXnmj3M8Z++He37n3LvvOb9zLt/PWiz2Pvuc89scLg97P3uf/aOZQUSkpw7IuwMi0toUREQkiYKIiCRREBGRJAoiIpKkX94d6A3kYQYclXc3RPq4bTvN7Gj/aJ8IIoUAMj3vToj0cV/b0tWjms6ISBIFERFJoiAiIkkUREQkiYKIiCRREBGRJAoiIpJEQUREkiiIiEgSBRERSaIgIiJJFEREJImCiIgkURARkSQKIiKSREFERJL0kUOJmt1BvfQ+r/XS+4j0Ho1ERCSJgoiIJFEQEZEkyonkIvzYs/Il3cmBKF8i+dBIRESSKIiISBJNZ+rCT1H8x3xohbJ/7Z9dm6/3dLrj+6epkPScRiIikkRBRESSKIiISBLlRBrC5yCODMoDMp77gmvz9RczruHzJ7219V4kppGIiCRREBGRJLkFEZIjSf6cZDvJp0jOKD4+gORykhuLv78lrz6KSHV55kT2APiimT1G8ggAvya5HMA0AA+Y2VySMwHMBPCVHPvZA1n7QoA4DzIqo+1517bR1fcEZb/Xozt/tf61B2W0icRyG4mY2TYze6xYfhFAO4DhAKYAuK34tNsAnJtLB0WkJk2REyE5CsA7AKwBMNjMtgGFQAPgmAqvmU5yLcm1wEsN66uIxHJf4iXZH8APAVxhZi+QrOl1ZjYfwPzCewyz+vWwN/iPOZyynBg39Q+nEiPitt1+WrQmKPupj39u6EVX1/Kv9FyuIxGSB6EQQBaZ2Y+KD+8gObTYPhTAc3n1T0Sqy3N1hgC+C6DdzK4PmpYBmFosTwVwV6P7JiK1y3M6cyqATwF4kuTjxce+CmAugCUkLwXwOwAX5tM9EalFbkHEzFYBqJQAOaORfekdB1UoA5lLvP3dcycH5f7uZfeMjuvbw63ta+K2aPkXiJdqs/IlQLxlvrdOXpO+qilWZ0SkdSmIiEgSBRERSZL7PpH9Q8ZRAKNc04Sg7LaJ4ARX/8bJ5XKn3yeyLqM/R7i6z5/4YwREKtNIRESSKIiISBJNZ+rCf6wZH7NvCmYWb568PWoaePAfovqmQ8aWK1ec7t6n071xuPE3a/kXiJeAs06Y16nxopGIiCRSEBGRJAoiIpJEOZGGyFhC3Y6Kjjp4V1QfhJ1R/fXL15fKW/a49d8rPuHe7ZagXO3YgPCoAn/CvH+t7O80EhGRJAoiIpJEQUREkignUhfV9mEEOZGdrunxcnHLHXGe47CL4rNkz8J9pfKaGbuitie2T4jqmDsxqCxxF63W31DWSfDaN7I/0khERJIoiIhIEk1nGsJPF3YETTvipqWDy+VBcVPHxOOj+siBW0vlv8Uvo7aX5sTLths3v71cWXyq68/Drh5ukfc3HD+mwvMATV/2TxqJiEgSBRERSaIgIiJJlBOpi4wlXQDxVvLNcdPuICeywr3rniOj+oqrysu2bxu8IWr7eyyI6rfcflmpvGnnafEbr/i96187KguXcf12+VqXhqs9V1qJRiIikkRBRESSaDrTa8Lhuf9Y/dA9vKG2/1bsr8vFtlGubWBU23vK4aXyQ1PfG7V9HIui+hfw9VJ51tI5UdvuIRfFl9n9naDil3HDKYw/8Nnz3wCWvkgjERFJoiAiIkkUREQkSa45EZK3onAL6+fMbFzxsQEA7kDhtk6bAXzMzP6YVx9rFy5hVrsZVPjcP7i2MH8SL+miLc6JYF65+MSe+Fu7gy+Nt9P/A24slb95+Iyo7bObb457MOgzQe1brn9hDsdviR/u6mFOREu6fVXeI5EFACa5x2YCeMDMRgN4oFgXkSaVaxAxswex7/LEFAC3Fcu3ATi3kX0Ske7JeyTSlcFmtg0Air8f09WTSE4nuZbkWuClrp4iIg3QsvtEzGw+gPkAQA6znLuD7DvDZT33RdcWvtZvP3fHBrQFX+n/AaOm+/9qSlQf9c7NpXK4ZwQAFg2MT4b/2D13lyuTP+L6EJ6K5rfLH+/qYU5He0b6qmYciewgORQAir/73U4i0kSaMYgsAzC1WJ4K4K4c+yIiVeS9xHs7gNMBDCLZCeBqAHMBLCF5KYDfAbgwvx72lqwlX98W5pn9t2R9/cly8d6T4qaX4+nN/BvKy7qjxm6O2i6LbmwFfP/D55bKH//B0vh9LwinPvHW+n1z5F2ms4o0vekrcg0iZnZJhaYzGtoREemxZpzOiEgLURARkSQtu8Tb3Krd1CnruVk6XD3MQbj3WfnOuD63XLz6G1+Lmk4auC6qX/KLci771Y/GxwRM+9Id5cq8d7n++FPjw9yPPzYg6/Q3nYLWSjQSEZEkCiIikkRBRESSKCeSi6z5f8jnAo7MaPdb0N3+k4Xjy69yp8Z/ckG832P5+84slacujG/+Pey6Z0vlD659KL7GSr8vJHyt3+My2NXD/utG4a1EIxERSaIgIiJJNJ1piO4s+YbTEP/X409IC+v+28D+ucF0YvHJUcufXh4S1c+5s7zE2zbtL6O2hz+1qlx53F1i9pi4vjRcAl7jnuxPQQtPSXPfVt5H+PlpapM3jUREJImCiIgkURARkSTKieSi1htf+2MCPL9sGvJftW+v3LY0zpFsH/nWUvmY3++O2vYOeb1cWRm/zYl3PhZf8dqzy5Wr/I+a3yIfnoqWdcfAarLyJVoqrgeNREQkiYKIiCTRdKbpZB347Kc3WcP8rKmO55aDO8vLr3s73LLtZUH5iuujpvbt/xzVx815tFRuW39m1IaFm10fwh2rfvk362S47kxRNH2pB41ERCSJgoiIJFEQEZEkyok0tWpz+DAfUO2vMsyfZG2fB6Jv2J7m7gu2OjxFPs6BYG6cI2nrV27/wPfuidp+1v8zUR03Lwgq1U6N999YDoWfg8+lZC35Kl/SUxqJiEgSBRERSaIgIiJJlBNpaa9VKAP75gP6VSh3VQ+/iu9uKj6hfPr7Adv/L2ra63Mk18wuFX+2e3bUdOJNbot827RyZdW3XH/8nyU8FS3r2IDu5Dm6cxN2CWkkIiJJFEREJEnF6QzJuwFYpXYzO6cuPZIeqjYcD2XdRBwAngvKx7q2BaXS3hOmRS0DdsZLr8/3m12ufGN21NaOuD7uoWCL/Ls/H19y9U2uDwMqlIHsrwJkbZH3y9xSq6yRyDwA/wngGRQ+4e8Uf+0G0FbvjpGcRHIDyQ6SM+t9PRHpmYojETP7BQCQ/Hcze2/QdDfJB+vZKZIHAvg2gDMBdAJ4lOQyM/tNPa8rIt1XS07kaJKlU2pIHgfg6Pp1CQAwHkCHmW0ys1cBLAYwpc7XFJEeqGWJ9wsAVpLcVKyPAjC9bj0qGA5ga1DvBBDdPZrk9HI/3lzn7rSCakuQWX/VPh8Q5g58viS48dWuOFfx/Gmfi+rjfhXkOQ6ZHb+Ny5G0vVyuD1jlciuD4vfFrgVBJSsn4o9DyDopTku8PZWVWB1qZtvM7F6SowGcUGxab2av1Llf7OKxKMlrZvMBzAcAcljFBLCI1FfWf0+3knwLCidp3gtglZlVO/Szt3QCGBnURwB4tsJzRSRHFXMiZnY2gNNRCCLnAVhN8kckp5P8izr361EAo0keR/JNAC4GsKzO1xSRHsjMiZjZyySPBzDbzP5YTKqeDeAGkkPMbHzW63vKzPaQ/EcA9wE4EMCtZvZUPa7Vd1Sb04eDyO7sifA5kfA6LuewOt6u3vae8n6PEx9y29zdPpEwR/L8nrjtgPVue/2QaUHtDte/k4Ky3xLv/yw6SrE31JJYHQLgEZL/A+BWADeZ2Y3FEULdmNlPAPykntcQkXRVl3jN7CoAbwPwXQDTAGwkeS0KeQoR2c/V9C1eMzOS2wFsR2FcPADAD0kuN7N/qWcHpVbdOQWtO21eOEXwp7I7qxaUiu2fmxY1feAmf9LZ7HLlmviEtH2+Hbw6KE/4mLvok0E56wQ0IP5z+8/A/9MIp4Oa6oSqBhGSnwcwFcBOALcA+LKZvUbyAAAbASiIiOzHahmJDAJwvpltCR80s70kJ9enWyLSKqoGETP714y29kptIrJ/0Mlm+43uLGdm3YUv9JyrZ5zKfvN3opaf7Y5Pex/3vWCL/J7sU+SjU+ZXuc3Np4U/0j5n47fBhzc290cI9HSL/P6XL9GhRCKSREFERJIoiIhIEuVE9ktZxwT69qw2v33ebysPcxIdcdPC5VG17YQzS+UT57gt8kNcjuSKci7jgGsOjJr2jhhTrnT6/ZDrUDufP8k6SjHMn1T7bPsejUREJImCiIgk0XRmv1RtiN3TIXjWN379cuvDcfWq8pSgfe3ZcdvsuHp1cJLdP7m3Hbfn6VJ5+3lvjRuXnub6EE5v/JYn/2fJmrJUel615/aNqY5GIiKSREFERJIoiIhIEuVEpJtq3RIPxHkFnxMZ7OprysWl7q57u8ZE1VOtnNsYuHBV1LYsuLPImYvjZeQ/TRsSv+/i8BQ0v2zbnS3y4WeSdYSAf27fWA7WSEREkiiIiEgSTWckQdZuVt/uD032N50KboqFO+OmlfG3hT/45YdK5QXXXRS1Tf3FklJ54fs+EbWdf0P8vq/1C665MLo3GoBfu3rW4dZh36sdgh3+k2vUHVjqSyMREUmiICIiSRRERCSJciLSi7KWKH2uwC+THpPRtiauzju+VJw2Ib551Zs+Wr5N9EQ8ELV9beDVUf2rM79ernS6S678a/dAmO/x+Z1qp8qHsj6j1lzy1UhERJIoiIhIEgUREUminIjk5IWMNr9F3n8tf1G5eEF8B7yP37O0VF7y4Y9Ebee5/Sebx44qlefPmxFf4ip3ivy94RZ53/dwz0u1097Cf3JZJ6S1jlxGIiQvJPkUyb0kT3Fts0h2kNxA8qw8+icitctrJNIG4HwA/xU+SHIMgIsBjAUwDMAKkm8zs9cb30URqUUuQeSNO+eR9E1TACw2s1cAPEOyA8B4AL9qbA+l8cIpgp/O+PrmoHx33DT586XiJ3YuippuHvjZqP4h/KT8ju8cFbXdf8GUqI7O4Ge1zZ+QtjEo+xPSjnD1cPm62hb5rAOzm0ezJVaHA9ga1DtR4dbzJKeTXEtyLfBSQzonIvuq20iE5AoAQ7poutLM7qr0si4es66eaGbzAcwvXGtYl88RkfqrWxAxs4k9eFkngJFBfQSAZ3unRyJSD822xLsMwPdJXo9CYnU0gEfy7ZLUR9YWb7+N3OcVwq/e+y3y5RuHvzYqvmn4jO3fjOpzDp9VKp+BFVHbjkvjk9ee6DehXJnnLtkWLvH62bf/c4bLuD7P4XMk1XImzSGvJd7zSHYCeDeAH5O8DwDM7CkASwD8BsC9AC7XyoxIc8trdeZO7HPyTKntWgDXNrZHItJTzbY6IyItptlyIiLI3hIPZG8zD/Ipu+NjAnafGx+l+PXlXyiVL8MtUdt78GBUf3JSedv73rWHu2sOLBfbBro2d4xBtB2hWg5E+0REZD+gICIiSTSdkZx0Z3juh/m1ng62OW5aEd/oatMl5e3r/337tKjNn4o2cXB5Cfjns98f9+aaYMnZXRK7R7kHsm6C1Zrf6tVIRESSKIiISBIFERFJopyINKFqp56Hd8TzW+KzTo1/OK4uLm9R3zjq7VHTYXPi/MS7gqXa4wd2RG3tE4KT4Xe6Sy52Ny7fE9b9qfGtSSMREUmiICIiSRRERCSJciLShHwOpDt3hguf63+8/euCr//PvSxqeWLIhKi+a8ZRpfJgl8s49qL1pfKWzhMqXgIAsP3QCn0FWvWfo0YiIpJEQUREkrTm+EmkxC/jhlvF/Y/3ka4evjY+GR5f+mRU3dKvPE058PL4nKxdrxxVuXv+lOFoOtM3/vlpJCIiSRRERCSJgoiIJOkbkzLp46ot+YbC7eoDXJv/cQ/r7sbbe1bG9bmnl4qbXh4bt43KuETmt/l9Y2t89d/TSEREkiiIiEgSTWekD6t286fwx99PkTbG1c5gajTv5LhtclDe7d5ms79meLJZ8x6+3B0aiYhIEgUREUmiICIiSZQTkRZU67d4q+VEwlPR/Hv6+rpyMdq6DmDx6MqX2O3fJ1xKrnaKfWvkTDQSEZEkuQQRkteRXE9yHck7SR4VtM0i2UFyA8mz8uifiNQur5HIcgDjzOxkAL8FMAsASI4BcDGAsQAmAbiR5IE59VFEapBLTsTM7g+qqwFcUCxPAbDYzF4B8AzJDgDjAfyqwV2UPsHnFFwuIzNn4p8bHhuwLm7anbXV/vmMemtuc/eaISfyaQA/LZaHA9gatHUWH9sHyekk15JcG99pXUQaqW4jEZIrsO+RLABwpZndVXzOlSiE4zdOhGEXz7eu3t/M5gOYX3ifYV0+R0Tqr25BxMwmZrWTnIrChuEzzOyNINAJYGTwtBEAnq1PD6VvyloW7c70Iesbtn6K0h6U/elp3Vnibc3pTV6rM5MAfAXAOWYWzkWWAbiY5MEkjwMwGsAjefRRRGqT12azGwAcDGA5SQBYbWafNbOnSC4B8BsUwvLlZvZ6xvuISM7yWp05PqPtWgDXNrA7IpJA295lPxbmK/xRAD5fEf5TyTphvtpW9j9ntLWmZljiFZEWpiAiIkkUREQkiXIish/JykF0Jz/h93OEr612/ED43Kz3aR0aiYhIEgUREUmi6YxITWqdarTm1vUUGomISBIFERFJoiAiIkmUExFJ1ppLs71FIxERSaIgIiJJFEREJImCiIgkURARkSQKIiKSREFERJIoiIhIEgUREUmiICIiSRRERCSJgoiIJFEQEZEkCiIikkRBRESSKIiISBIFERFJQjPLuw/JSP4vgC1596NoEICdeXfCUZ9qoz5lO9bMjvYP9okg0kxIrjWzU/LuR0h9qo361DOazohIEgUREUmiINL75ufdgS6oT7VRn3pAORERSaKRiIgkURARkSQKIr2E5HUk15NcR/JOkkcFbbNIdpDcQPKsBvbpQpJPkdxL8hTXllefJhWv2UFyZqOu20U/biX5HMm24LEBJJeT3Fj8/S0N7M9Ikj8n2V78O5uRd59qpSDSe5YDGGdmJwP4LYBZAEByDICLAYwFMAnAjSQPbFCf2gCcD+DB8MG8+lS8xrcBnA1gDIBLin3JwwIU/uyhmQAeMLPRAB4o1htlD4AvmtmJACYAuLz42eTZp5ooiPQSM7vfzPYUq6sBjCiWpwBYbGavmNkzADoAjG9Qn9rNbEMXTXn1aTyADjPbZGavAlhc7EvDmdmDAJ53D08BcFuxfBuAcxvYn21m9lix/CKAdgDD8+xTrRRE6uPTAH5aLA8HsDVo6yw+lqe8+tSMn0VosJltAwr/qAEck0cnSI4C8A4Aa5qlT1n65d2BVkJyBYAhXTRdaWZ3FZ9zJQpD00VvvKyL5/faunotferqZfXsU4a8rtsySPYH8EMAV5jZC2RXH1lzURDpBjObmNVOciqAyQDOsPIGnE4AI4OnjQDwbKP6VEFd+9SE163VDpJDzWwbyaEAnmvkxUkehEIAWWRmP2qGPtVC05leQnISgK8AOMfMXgqalgG4mOTBJI8DMBrAI3n0sQn69CiA0SSPI/kmFJK7yxpw3VotAzC1WJ4KoNJIrtexMOT4LoB2M7u+GfpUMzPTr174hUJyciuAx4u/bg7argTwNIANAM5uYJ/OQ+F//1cA7ABwXxP06UMorF49jcKUK6+/r9sBbAPwWvEzuhTAQBRWQDYWfx/QwP6chsLUbl3wM/ShPPtU6y9texeRJJrOiEgSBRERSaIgIiJJFEREJImCiIgkURCRXJC8l+Rw99jpJO/Jq0/SMwoi0nAkD0Vhv8Pv8+6LpFMQkboh+TfF81UOIXl48ZyMcQBOB7Cy+JxJxXNYVqFwbIG0GG02k7oieQ2AQwAcCqDTzOaQ/BaApQB+icJOzA+gsOP3DgCHmdnknLorPaCRiNTbvwE4E8ApAP6j+NipAFYBOAHAM2a20Qr/my3Mp4uSQt/ilXobAKA/gIMAHEJyMICtZvZq8WvuGgq3OE1npK5ILkPhBLPjAAxF4cSu183sZpKHoPBlvPeb2dMkbwdwhKYzrUUjEakbkn8HYI+Zfb94vuovUfg6+0kAYGYvk5wO4Mckd6IwxRmXW4elRzQSkYYheTCAh63Jb1At3aMgIiJJtDojIkkUREQkiYKIiCRREBGRJAoiIpJEQUREkvw/ONWxIkKaG3EAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.imshow(real_space_emission, origin='lower', aspect='equal', cmap='jet',\n",
    "          extent=[min(x_pos), max(x_pos), min(y_pos), max(y_pos)], vmin=0)\n",
    "\n",
    "ax.set_xlabel('x/d')\n",
    "ax.set_ylabel('y/d');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1aa52b77-00a0-4391-a704-8be4e6e1a190",
   "metadata": {},
   "source": [
    "## Cluster"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be8d7003-a358-4de8-99b4-1f66dfee11bb",
   "metadata": {},
   "source": [
    "### Individual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "75b23624-ea0d-4561-90de-c2903e50b7eb",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-03T13:09:52.633732Z",
     "iopub.status.busy": "2022-02-03T13:09:52.633732Z",
     "iopub.status.idle": "2022-02-03T13:09:52.647730Z",
     "shell.execute_reply": "2022-02-03T13:09:52.647730Z",
     "shell.execute_reply.started": "2022-02-03T13:09:52.633732Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "dc = 0.45 * d\n",
    "emitter = [[-dc / np.sqrt(2), dc / np.sqrt(2), 0],\n",
    "           [dc / np.sqrt(2), -dc / np.sqrt(2), 0]]\n",
    "\n",
    "n_e = 2\n",
    "signs = [1, 1]\n",
    "psi0 = np.zeros(N_x * N_y + n_e, dtype=complex)\n",
    "for i in range(n_e):\n",
    "    psi0[-i - 1] = signs[i] / np.sqrt(n_e)\n",
    "\n",
    "shift = compute_shift(emitter[0], emitter[1], pol_emi[0], pol_emi[0], k0, gamma_emi, psi0[-2:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "fcb3528c-ed60-457a-a418-e9bd0034c5b0",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-03T13:09:52.648732Z",
     "iopub.status.busy": "2022-02-03T13:09:52.648732Z",
     "iopub.status.idle": "2022-02-03T13:10:24.106674Z",
     "shell.execute_reply": "2022-02-03T13:10:24.106674Z",
     "shell.execute_reply.started": "2022-02-03T13:09:52.648732Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "results = complete_dynamics([x_pos, y_pos, z_pos], k0, pol_lat, dt, tf, psi0, \n",
    "                            emitter=emitter, pol_emi=pol_emi * n_e, gamma_emi=gamma_emi,\n",
    "                            omega_emi=omega_emi - np.real(shift), border=border,\n",
    "                            max_gamma=max_gamma)\n",
    "\n",
    "population_ind = np.sum(np.abs(results['psi'][:, -n_e:]) ** 2, axis=-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3508be8-1f79-48e3-aa1b-8a7a8aa1986d",
   "metadata": {},
   "source": [
    "### Collective"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "2abeef6f-edcb-4cbf-ac00-2ed6f00c25ad",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-03T13:10:24.107674Z",
     "iopub.status.busy": "2022-02-03T13:10:24.107674Z",
     "iopub.status.idle": "2022-02-03T13:10:24.121676Z",
     "shell.execute_reply": "2022-02-03T13:10:24.121676Z",
     "shell.execute_reply.started": "2022-02-03T13:10:24.107674Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "dc = 0.45 * d\n",
    "n_plaquette = 5\n",
    "emitter = [[n_plaquette + -dc / np.sqrt(2), n_plaquette + dc / np.sqrt(2), 0],\n",
    "           [n_plaquette + dc / np.sqrt(2),  n_plaquette -dc / np.sqrt(2), 0],\n",
    "           [-n_plaquette + -dc / np.sqrt(2), -n_plaquette + dc / np.sqrt(2), 0],\n",
    "           [-n_plaquette + dc / np.sqrt(2),  -n_plaquette -dc / np.sqrt(2), 0]]\n",
    "\n",
    "n_e = 4\n",
    "signs = [1, 1, 1, 1]\n",
    "psi0 = np.zeros(N_x * N_y + n_e, dtype=complex)\n",
    "for i in range(n_e):\n",
    "    psi0[-i - 1] = signs[i] / np.sqrt(n_e)\n",
    "\n",
    "shift = compute_shift(emitter[0], emitter[1], pol_emi[0], pol_emi[0], k0, gamma_emi, psi0[-2:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "a616bb5b-203e-472c-8681-3f513d3efc8d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-03T13:10:24.123674Z",
     "iopub.status.busy": "2022-02-03T13:10:24.122675Z",
     "iopub.status.idle": "2022-02-03T13:10:56.558540Z",
     "shell.execute_reply": "2022-02-03T13:10:56.558540Z",
     "shell.execute_reply.started": "2022-02-03T13:10:24.123674Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "results = complete_dynamics([x_pos, y_pos, z_pos], k0, pol_lat, dt, tf, psi0, \n",
    "                            emitter=emitter, pol_emi=pol_emi * n_e, gamma_emi=gamma_emi,\n",
    "                            omega_emi=omega_emi - np.real(shift), border=border,\n",
    "                            max_gamma=max_gamma)\n",
    "\n",
    "population_coll = np.sum(np.abs(results['psi'][:, -n_e:]) ** 2, axis=-1)\n",
    "\n",
    "population = np.abs(results['psi'][:, :-n_e].reshape((-1, N_y, N_x))) ** 2\n",
    "index_max = np.argmax(np.sum(population, axis=(1, 2)))\n",
    "real_space_emission = population[index_max]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "1baa4e74-f84f-4517-bde8-7c89282aaf95",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-03T13:10:56.560543Z",
     "iopub.status.busy": "2022-02-03T13:10:56.560543Z",
     "iopub.status.idle": "2022-02-03T13:10:56.666540Z",
     "shell.execute_reply": "2022-02-03T13:10:56.666540Z",
     "shell.execute_reply.started": "2022-02-03T13:10:56.560543Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'population')"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEOCAYAAABmVAtTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAnR0lEQVR4nO3deXxV9Z3/8dcnCSENBAKyE0gCggoKmLC44VIHIShQq+MyUxFa6zgVx6oztrVarR2ptj8dW3e0jkqlCi7VsVa0IiIiCqggYd+JgOw7AZJ8f3+cS7wkITkh995z7837+XicR3LW++E+jnl7zvd8v8ecc4iIiACkBF2AiIjED4WCiIhUUiiIiEglhYKIiFRSKIiISCWFgoiIVIpaKJjZs2a22cwWHmO9mdkfzWyFmS0ws4Jo1SIiIv5E80rhOWBYLeuLgB6h6XrgiSjWIiIiPkQtFJxzM4DttWwyCnjBeWYD2WbWMVr1iIhI3YJsU+gMrA+bLwktExGRgKQF+NlWw7Iax9wws+vxbjFxAhTmnn46lqI2chERv+bNm7fVOde2ru2CDIUSoEvYfA6woaYNnXMTgAkA/c3c3ydOpG3v3tGvUEQkSZjZWj/bBfm/228Co0NPIZ0B7HLObfSz4561vv5tIiJST1G7UjCzvwDnA23MrAS4G2gC4Jx7EngbGA6sAPYDY/0ee9/69XVvJCIi9Ra1UHDOXV3HegfceDzHLt3o64JCRETqKeFaaxcAa3r0CLoMEZGklHChcBjYXVoadBkiIkkp4UKhI9Bi/vygyxARSUoJFwqdgLbFxUGXISKSlBIuFBzAvn1BlyEikpQSLhTKgZS9e4MuQ0QkKSVeKJjRdOfOoMsQEUlKiRcKKSlk7toVdBkiIkkp4UJhS8uWFGVlBV2GiEhSSrhQyMjM5OuNG1m3bl3QpYiIJJ2EC4WsjAz+F1j74otBlyIiknQSLhSapKczBuCLLwKuREQk+SReKGRkcAio0EipIiIRl3ChgBlbUlNJ3bw56EpERJJO4oUCsPs736GJ+iqIiERcQobCzlat2HvwYNBliIgknYQMhVevuIIRFRV47+kREZFISchQyMnJ4cCBA+zYsSPoUkREkkpChkLhli28C2xYtCjoUkREkkpChkKHpk0ZAmxbuDDoUkREkkpChkL2yScDsGfp0oArERFJLgkZCq1OPx2AwytWBFyJiEhySchQSMvLoxxI1aB4IiIRlZChQJMmzM/KYvP+/UFXIiKSVBIzFIDxF13Eg2lpQZchIpJUEjYUcnJyWL9+vTqwiYhEUMKGQtHXX/Ppvn3s1hhIIiIRk7Ch0KZ5c3oDG7/6KuhSRESSRsKGQuaJJwKwXR3YREQiJmFDIbtXL0Ad2EREIilhQ+GEU08FoHT16oArERFJHgkbCuldu/JRejpr9+4NuhQRkaSRsKFA06bcMXAgr5aXB12JiEjSSNxQALp3787KlSuDLkNEJGkkdCjcvGABf/76a0pLS4MuRUQkKUQ1FMxsmJktNbMVZvbzGta3NLP/M7P5ZlZsZmPrc/ysrCx6AGvWrIlUySIijVrUQsHMUoHHgCKgF3C1mfWqstmNwCLnXF/gfOBBM0v3+xkZ+fm0B1YvXx6ZokVEGrloXikMBFY451Y55w4BLwGjqmzjgCwzM6A5sB0o8/sBLXr1Ig3Y+OWXkalYRKSRi2YodAbWh82XhJaFexQ4BdgAfAXc7JyrqHogM7vezOaa2dwtW7ZULs8KvYFtl97VLCISEdEMBathWdUhTYcCXwKdgH7Ao2bWotpOzk1wzvV3zvVv27bttx9wyin8NTubVZs3R6xoEZHGLJqhUAJ0CZvPwbsiCDcWeM15VgCrgZN9f0KPHjxz9tnM3L69obWKiAjRDYU5QA8zyw81Hl8FvFllm3XAhQBm1h44CVhVnw/Jy83l61X12kVERI4haq8uc86Vmdk4YCqQCjzrnCs2sxtC658EfgM8Z2Zf4d1u+plzbmt9PufeKVMYuHs3O3fuJDs7O7L/CBGRRiaq77N0zr0NvF1l2ZNhv28ALmrQZ7RqRactW1i1ahUFBQUNOZSISKOX0D2aAVK7dqUTaLgLEZEISPhQyDzxRDoBy9WBTUSkwaJ6+ygW0nNzSQfWLV4cdCkiIgkv4a8UOPdcnu/aldUrVgRdiYhIwkv8UDjrLD688EK+Wrs26EpERBJe4oeCc/Rq3579Gzdy4MCBoKsREUloiR8KBw/yn/ffz43AWl0tiIg0SOKHQkYGh7Oz6QKsXr066GpERBJa4ocCQE4OXYBVGu5CRKRBkiIU0rp1I89MHdhERBooKULBunQhJyVFHdhERBooKUKBK6/kpb59WbZ0adCViIgktOQIhcGDWTt0KKtWr6aszPfbPEVEpIrkCIVDhxiUkUHLsjLWrFkTdDUiIgkrOUJh3TpG3X03FwPLli0LuhoRkYSVHKGQkwN47/5UKIiIHL/kCIWMDFybNnRPT2epGptFRI5bcoQCYDk59MjMVCiIiDRA0oQCOTl0NVMoiIg0QPKEwu23M/PSS9mwYQO7d+8OuhoRkYSUPKEweDAZF18MqLFZROR4JU8o7NhB4caNtAHdQhIROU7JEwrLlpE7bhxnql1BROS4JU8o5OYCUHDCCSxZsiTgYkREElPyhEK7dtC0Kae1aKFQEBE5TskTCikp0LUr3Zs0YdmyZRoYT0TkOCRPKADk5tLp8GEOHjyoV3OKiByH5AqFhx5iw333AbB48eKAixERSTzJFQqnnUZ+UREAixYtCrgYEZHEkxZ0ARFVUkLLN9+kT4cOulIQETkOyXWlsHIl3Hgjwzp2pLi4OOhqREQSTnKFQqivwumtW7N48WIqKioCLkhEJLEkVyh07gwpKZyUkcH+/fv1ak4RkXpKrlBo0gQ6dyanvByAhQsXBlyQiEhiiWoomNkwM1tqZivM7OfH2OZ8M/vSzIrN7MMGf2huLq1CQ2erXUFEpH58P31kZqlA+/B9nHPr6tj+MWAIUALMMbM3nXOLwrbJBh4Hhjnn1plZu3r/C6qaNIm0Fi3octppCgURkXryFQpmdhNwN/ANcKT11gF9atltILDCObcqdIyXgFFAeAeCfwFeOxIuzrnN9aq+Jl26ANC7d2/dPhIRqSe/t49uBk5yzvV2zp0WmmoLBIDOwPqw+ZLQsnA9gVZmNt3M5pnZaJ/1HNvixfCLXzAoP58lS5ZQHmpfEBGRuvkNhfXArnoe22pY5qrMpwGFwMXAUOAuM+tZ7UBm15vZXDObu2XLlto/dd06uP9+BmZnc/DgQVauXFnPskVEGi+/bQqrgOlm9jfg4JGFzrmHatmnBOgSNp8DbKhhm63OuX3APjObAfQFjnqfpnNuAjABoH///lWD5WihvgonZ2QA3hNIPXtWyxkREamB3yuFdcB7QDqQFTbVZg7Qw8zyzSwduAp4s8o2bwCDzSzNzDKBQUDDxqfo2hWAzqHbRmpsFhHxz9eVgnPu1wBmluXNur0+9ikzs3HAVCAVeNY5V2xmN4TWP+mcW2xm7wAL8Bqwn3HONax1ODMT2ral6caN5OfnKxREROrB79NHpwITgdah+a3AaOdcrX9xnXNvA29XWfZklfnfA7+vR811y82FzZv1BJKISD35vX00AbjVOZfrnMsFbgOejl5ZDTR9Orz+Or1792bZsmUcOnQo6IpERBKC31Bo5pz74MiMc2460CwqFUVCs2ZgRt++fTl8+LDerSAi4pPfUFhlZneZWV5ouhOI3/ddfv45/OAHDMrJAWDevHkBFyQikhj8hsIPgbbAa8Drod/HRquoBtu9G158kbwDB8jKylIoiIj45Pfpox3Af0S5lsjp1g2AlDVrKCgo4PPPPw+4IBGRxFBrKJjZw865n5rZ/1G9NzLOuZFRq6whOneG9HRYuZLCwkIef/xxysrKSEtLrrePiohEWl1/JSeGfv6/aBcSUampkJ/vhcL3v09paSmLFi2iT5+6hmsSEWncam1TcM4duRnfzzn3YfgE9It6dQ3Rpw+YUVBQAKixWUTED78NzdfWsGxMBOuIvMmTYcoUevbsSfPmzRUKIiI+1NWmcDXeOw/yzSx83KIsYFs0C4uUlJQUTj/9dDU2i4j4UFebwixgI9AGeDBs+R688Yri14IFcNNN8PDDFBYW8tRTT6mxWUSkDrX+hXTOrQXWAmfGppwIatIEZsyA4mIKCws5cOAAS5Ys4dRTTw26MhGRuOWrTcHMzjCzOWa218wOmVm5me2OdnEN0q0bpKTA8uUUFhYCamwWEamL34bmR4GrgeXAd4DrgEeiVVRENG0KeXmwbBk9e/akWbNmCgURkTr4vsHunFthZqnOuXLgf81sVhTrioyePWHpUlJTU9XYLCLig98rhf2ht6d9aWa/M7NbiOdRUo8YPBhOPBGAgoICvvjiC8pDb2QTEZHq/IbCNXhvTxsH7MN79/Jl0SoqYu64w+uvABQWFrJ//36WLl0acFEiIvHL74B4a0O/HgB+Hb1yoie8sblXr14BVyMiEp9qvVIws6/MbMGxplgVedy2bfOeQpowgZNPPpnMzEw1NouI1KKuK4VLYlJFtLRqBd98A4sXk5qaSr9+/dTYLCJSCz+d1xJXSkrlE0jg3UJ69tlnKS8vJzU1NeDiRETij9/Oa3vMbHdoKk2IzmtHnHQSLFsGwKBBg9i3bx8LFy4MuCgRkfjkKxScc1nOuRahKQPvyaNHo1tahPTsCatXw8GDnHmmN1rHJ598EnBRIiLxye8jqUdxzv0V+G5kS4mS886D666DffvIz8+nXbt2CgURkWPw9UiqmX0/bDYF6E8Nr+eMSxde6E2AAWeeeaZCQUTkGPxeKYwIm4biDZ09KlpFRVx5OezcCXihsHz5crZu3RpsTSIicchv57Wx0S4kqk4/3WtbeOWVynaF2bNnc8klif3ErYhIpPl9+qibmf2fmW0xs81m9oaZdYt2cRHTvTsUFwPQv39/0tLSmDUr/sfzExGJNb+3jyYBk4GOQCdgCvCXaBUVcb16wfLlcOgQmZmZFBQU8NFHHwVdlYhI3PEbCuacm+icKwtNfyZRGpoBevf22hVC/RXOO+88PvvsM/bv3x9wYSIi8cVvKHxgZj83szwzyzWz24G/mVlrM2sdzQIjondv72foFtL555/PoUOHmD17doBFiYjEH78v2bky9PPfqiz/Id4VQ3y3L5x0Etx7L4Tez3zOOeeQkpLC9OnT+e53E6O7hYhILPh9+ig/2oVEVUYG3HVX5WyLFi0oKCjgww8/DLAoEZH44/fpoyZm9h9m9kpoGmdmTaJdXETt3AlhndbOP/98Zs+ezYEDB4KrSUQkzvhtU3gCKAQeD02FoWWJ45FH4OyzYe9ewGtsPnToEJ9++mnAhYmIxA+/oTDAOXetc25aaBoLDKhrJzMbZmZLzWyFmf28lu0GhEZevdxv4fXWty84B6ERUsPbFURExOM3FMrNrPuRmVDHtfLadjCzVOAxoAjoBVxtZtXegxna7gFgqt+ij0ufPt7P+fMByM7Opl+/fgoFEZEwfkPhv/AeS51uZtOBacBtdewzEFjhnFvlnDsEvETN4yXdBLwKbPZZy/HJzYWWLStDAb5tVygtLY3qR4uIJAq/ofAx8BRQEZqeAuoaarQzsD5sviS0rJKZdQYuBZ6s7UBmdr2ZzTWzuVu2bPFZcrWDeFcLVULh4MGDalcQEQnxGwovAPnAb0JTPjCxjn2shmVVe0E/DPzMOVfrrSjn3ATnXH/nXP+2bdv6q7gm99/vNTiHDB48GDPjgw8+OP5jiogkEb+d105yzvUNm//AzOYfc2tPCdAlbD4H2FBlm/7AS2YG0AYYbmZloZf4RN5ZZx01m52dzYABA5g6dSr33HNPVD5SRCSR+L1S+MLMzjgyY2aD8G4p1WYO0MPM8s0sHbgKeDN8A+dcvnMuzzmXB7wC/CRqgQBQWgqTJh11C6moqIhPP/2Ubdu2Re1jRUQShd9QGATMMrM1ZrYGrz3hPDP7yswW1LSDc64MGIf3VNFiYLJzrtjMbjCzGyJQe/2ZwZgx8JdvB3gtKirCOce7774bSEkiIvHEnKt7sFMzy61tvXNubcQqqkP//v3d3LlzG3IAyM6Gf/wDgPLycjp06MCwYcOYOLGuZhIRkcRkZvOcc/3r2s7v2Ecx+6Mfdf37w0svQUUFpKSQmprK0KFDmTp1KhUVFaSk+L14EhFJPo3vL+CAAbBrF6xYUbmoqKiILVu2MG/evAALExEJXuMMBYAvvqhcNHToUMyMt99+O6CiRETiQ+MLhd69Yd06uOKKykVt2rRh4MCB/P3vfw+wMBGR4DW+UEhNhS5dvCeRwhQVFfHZZ5+xdevWgAoTEQle4wsFgA8/hNGjoaysctHw4cNxzjF1anTH5RMRiWeNMxRKSmDiRFi0qHJRYWEhbdu21S0kEWnUGmconBHqnD1rVuWilJQUioqKePvttzl8+HBAhYmIBKtxhkK3btChA8ycedTiyy67jB07dmiAPBFptBpnKJjBOedUC4WLLrqI5s2b88orrwRUmIhIsBpnKACce6433MWePZWLMjIyGDFiBK+//jplYY3QIiKNReMNhXHj4MsvISvrqMWXXXYZW7duZcaMGcHUJSISoMYbClbTO4C8/gqZmZm8+uqrMS5IRCR4jTcUAH7xC+82UpjMzEyGDx/Oa6+9Rnl5rS+EExFJOo07FDIyvMbmXbuOWnz55ZezadMmZoU9sioi0hg07lA45xxwDj755KjFw4cPJyMjgylTpgRUmIhIMBp3KAwa5I2FVOXR1KysLIqKipg8ebKeQhKRRqVxh0Lz5nD66fDRR9VWXXPNNXzzzTe89957ARQmIhKMxh0KAGPHQlFRtcUXX3wxrVu35oUXXgigKBGRYPh6HWdS+8lPalycnp7O1VdfzZ/+9Cd27dpFy5YtY1yYiEjs6UoBYN8+KC6utnj06NGUlpaqwVlEGg2FAsAPfgAjR3pPIoUZMGAAJ598Ms8//3xAhYmIxJZCAWDIEFi1CpYvP2qxmTF69GhmzpzJypUrAypORCR2FAoAw4Z5P995p9qqa665hpSUFJ555pkYFyUiEnsKBfDer9CzJ9Tw1rWcnBxGjhzJM888Q2lpaQDFiYjEjkLhiGHDYPp0OHCg2qpx48axdetWXn755djXJSISQ+aqNK7Gu/79+7u5c+dG/sArVsDu3dCvH6QcnZXOOXr37k1mZiZz5szBjjHCqohIvDKzec65/nVtpyuFI048EQoKqgUCeA3O48aNY968eXz66acBFCciEhsKhXCzZsEvf1njqtGjR9OiRQseffTRGBclIhI7CoVwc+fC+PHVHk0FaN68OWPGjGHy5Mls2rQpgOJERKJPoRBu1Cjv52uv1bj6xhtv5PDhwzz99NMxLEpEJHYUCuFyc2HAAHjllRpX9+zZk6FDh/LEE09w8ODBGBcnIhJ9CoWqLr/cu420dm2Nq2+99VY2btzIxIkTY1yYiEj0KRSquuwy6NwZjjGsxZAhQygsLOT+++/XC3hEJOlENRTMbJiZLTWzFWb28xrW/6uZLQhNs8ysbzTr8aV7d1i/Hr773RpXmxl33HEHK1eu5JVj3GYSEUlUUeu8ZmapwDJgCFACzAGuds4tCtvmLGCxc26HmRUB9zjnBtV23Kh1XquqvBwOHoTMzGqrKioqOPXUU0lNTWX+/Pmk1NC3QUQknsRD57WBwArn3Crn3CHgJWBU+AbOuVnOuR2h2dlAThTr8W/fPsjLgwcfrHF1SkoKv/rVr1i4cCEvvvhibGsTEYmiaIZCZ2B92HxJaNmx/AioPiJdEJo183o4//nP1d6xcMQVV1xBYWEhd955pwbKE5GkEc1QqGmAoBr/wprZBXih8LNjrL/ezOaa2dwtW7ZEsMRaXHMNLFsGs2fXuDolJYXf//73rFu3Tr2cRSRpRDMUSoAuYfM5wIaqG5lZH+AZYJRzbltNB3LOTXDO9XfO9W/btm1Uiq3mn//Zu2L405+OuckFF1xAUVER9913H9u3b49NXSIiURTNUJgD9DCzfDNLB64C3gzfwMy6Aq8B1zjnlkWxlvrLyoIrr4SXXoI9e4652QMPPMCuXbsYP358DIsTEYmOqIWCc64MGAdMBRYDk51zxWZ2g5ndENrsV8AJwONm9qWZxeCxonq47TZvyIsankA64rTTTuPaa6/lkUceYc2aNbGrTUQkCvQ+hQgoKSmhR48ejBgxgsmTJwddjohINfHwSGpy2LED/vM/4ZNPjrlJTk4Od9xxB1OmTOHdd9+NYXEiIpGlK4W67NsHOTkwdKjXvnAMpaWl9OnTB4AFCxaQkZERqwpFROqkK4VIadYMfvhDePVVKCk55mYZGRk89thjLF++nHvvvTeGBYqIRI5CwY+bbvI6sT38cK2bDRkyhLFjx/LAAw8wZ86c2NQmIhJBCgU/8vK8x1OfesprY6jFQw89RIcOHRg7dqx6OotIwlEo+HX77TBiBBw4UOtm2dnZPP300xQXF3PbbbfFqDgRkchQKPjVty9MmgSdOtW56fDhw7n11lt5/PHHmTJlSgyKExGJDIVCfS1YcMzXdYb77W9/y6BBg/jRj37EymO8sEdEJN4oFOrr17+GsWNh8+ZaN0tPT+fll18mNTWVK664Qu90FpGEoFCor/HjvXYFH4+d5ubm8vzzz/P5559zyy23xKA4EZGGUSjU10knwY9/7D2JtKzuMfxGjhzJ7bffzhNPPMFjjz0WgwJFRI6fQuF43HMPfOc74PP//sePH88ll1zCzTffzNSpU6Nbm4hIAygUjkf79nDffdCnj/cu5zqkpqYyadIkevfuzWWXXcasWbNiUKSISP1p7KMY2rRpE4MHD2bLli1MmzaNgoKCoEsSkUZCYx/Fynvvwc9qfItoNR06dOD999+nZcuWXHTRRRQXF0e5OBGR+lEoNNTMmfC738Fbb/navGvXrrz//vukp6fzT//0TyxZsiTKBYqI+KdQaKhf/hJOOw2uu67WUVTDnXjiifzjH/+goqKCwYMH8/nnn0e5SBERfxQKDZWe7r1nYd8+uPTSOsdGOqJXr17MnDmTzMxMLrjgAr2cR0TigkIhEnr1ghdfhHnz4Nlnfe/Wo0cPZs6cSV5eHkVFRTz00EMkWsO/iCQXhUKkjBwJH30EP/lJvXbr0qULH3/8MZdeeim33XYb1157rYbcFpHAKBQi6eyzwQxWrYIPPvC9W/PmzZk8eTL33nsvEydO5Nxzz2X16tVRLFREpGYKhWj4t3+DSy6B6dN975KSksJdd93FX//6V5YuXUrfvn15/vnndTtJRGJKoRANf/4z5OfD8OEwbVq9dh01ahQLFiygoKCAMWPGMHLkSNatWxelQkVEjqZQiIb27b0w6N4dLr4Y3nyzXrvn5uby/vvv8+CDDzJt2jR69erFww8/TFlZWZQKFhHxKBSipV07r12hTx/44x+hnreBUlNTufXWWykuLua8887jlltuoV+/frzzzjtRKlhERKEQXW3aeFcMr7ziNUBv3w71fNlOXl4eb731Fq+99hqlpaUUFRUxZMgQZsyYofYGEYk4hUK0NWsG2dlQUeF1bhs4EBYurNchzIxLL72URYsW8T//8z/Mnz+f8847j3POOYe//e1vCgcRiRiFQqykpMDtt8OmTdC/v/fmNp+9n49IT0/npz/9KWvWrOGRRx6hpKSESy65hH79+vHcc8+xf//+KBUvIo2FQiGWLr4YvvoKRo2Cu++G3r3rfdUAkJmZybhx41ixYgXPP/88ZWVljB07lk6dOnHTTTexYMGCKBQvIo2BQiHW2rWDl1/22hp69IC8PG/5zp31PlSTJk0YPXo0CxcuZPr06QwfPpwJEybQt29fTj31VMaPH8/KlSsjWr6IJDe9ZCcelJXBKadAbi7ccIN3JdGkyXEdauvWrUyePJlJkybx8ccfA3DKKacwcuRIRowYwaBBg0hLS4tk9SKSAPy+ZEehEA8OHoQ//AEefRTWr/f6OYwZA//+715QHKc1a9bwxhtv8OabbzJjxgzKysrIyspi8ODBXHDBBVxwwQX069eP1NTUyP1bRCQuKRQSUXk5TJ0KEyZ4L+159VXvquHrr73bS716eY+2HoedO3fy7rvvMm3aNKZPn87SpUsBaNmyJQMGDDhq6ty5M3acnyMi8UmhkOi+/tq7YkhLg3vugV//2hs647zzYPBgbzrxxOMOiQ0bNjB9+nRmzJjBZ599xldffVXZY7pt27b06tXrqOmUU06hQ4cOCguRBKVQSCYbNnhDZbzzjvf6z23bIDMTdu+G1FSvc9yBA167RLdu0KpVvcOitLSU+fPnM2fOHObPn8+iRYtYtGgRO8MawLOyssjLy6s2de3alY4dO9KuXTuaHGdbiIhEV1yEgpkNA/4ApALPOOfur7LeQuuHA/uBMc65Wt9N2ShDIVxFBSxd6k3f+563bPBgLyyOaNkShgyBKVO8+ddf9/Zr186b2rf3tqkjOJxzbNq0icWLF1NcXMzKlStZvXo1a9asYfXq1ezZs6faPieccAIdOnSonNq2bUvr1q1p1apVtZ+tWrWiRYsWpKen6wpEJMoCDwUzSwWWAUOAEmAOcLVzblHYNsOBm/BCYRDwB+fcoNqO2+hDoSZlZV5ILF8Oq1d773No08brCwFeY3XVkVZHjoQ33vB+HzoUDh2CrKxvp3PPhX/5F2/9c895ne/S072notLTcd26sbNTJ9asWsXuadPYuns33+zYwTc7drBp2zZWbt/Oqq1b2bp5MxX79lEB1aYjZ15aWhrNmzc/asrKyqJ58+Y0a9aMpk2b0rRpUzIyMip/9ztlZGSQnp5OamoqaWlplT+P9Xv4TwWVJBO/oRDNZxMHAiucc6tCBb0EjAIWhW0zCnjBeck028yyzayjc25jFOtKPmlpXke43r1rXv/pp/DNN7B587c/u3b9dn3z5t6ydetgzx5vatLEC4WKCvjhD6sN6Gc330yrhx+m1cknez20q7rzTvjNb7we3B07Vlv9xVVXMfOss2iydi0/fugh3O7duD17KgPjj9268XLz5nTeto0n16/H4YWIcw4H3AC8ApwBHBmD1oX9vBaYCnwXmFjD+iuBWcAI4NEa1o8CFqakcCXw3xUVGOBCIeGA72dmsjYtjX89dIhbDh48aj3A5W3asDUtjdH79zNm796j/u3OjKu7duVAWhpjtm/ne7t2gdlRn391jx44M67bvJkLd+06av/SlBSu69kTgBs2buSs3bsrjwuwKzWVW0Lrx5WU0Dfs850ZW5o04Vc9egDw07VrOWnfvqOOX5KRwf2h9bevXEnXsJ73Bqxs1oyHu3cH4M6lS2lfZTyv4qwsnszPB+A3ixeTffjwUevnZWfzXOj8+11xMd8pLz9q/SetWzMpJweAPy5YQNVo/qBNG17r1In08nIeLC6mqqnt2/O3Dh1ocfgw/71oUbX1b3TsyPvt2tH24EHuWrKk2vrJOTnMPOEEcg4c4PZly6qtn9i1K3NataL73r3cXEM/oKfz8liYnU3v3bv5cQ0vy3qsWzeWZ2Vx+o4dXFvDsPgP9ejBusxMzti2jatKSqqtH3/SSWzOyOD8LVv43oYN1dbffcop7EpPZ+imTRR980219XWJ5pXC5cAw59x1oflrgEHOuXFh27wF3O+cmxmafx/4mXNubpVjXQ9cH5o9Fah/N+DYagNsDbqIOqjGyFCNDRfv9UFy1JjrnGtb10GieaVQ07V31QTysw3OuQnABAAzm+vnEihIqjEyVGNkxHuN8V4fNK4aoznMRQnQJWw+B6h6reNnGxERiZFohsIcoIeZ5ZtZOnAV397+PeJNYLR5zgB2qT1BRCQ4Ubt95JwrM7NxeO19qcCzzrliM7shtP5J4G28J49W4D2SOtbHoSdEqeRIUo2RoRojI95rjPf6oBHVmHCd10REJHo0dLaIiFRSKIiISKW4CgUzG2ZmS81shZn9vIb1ZmZ/DK1fYGYFfveNYY3/GqptgZnNMrO+YevWmNlXZvalmUWlW7aP+s43s12hGr40s1/53TeGNf5XWH0LzazczFqH1kX9Owx9zrNmttnMauwTE/S56KO+QM9DnzXGw7lYV43xcC52MbMPzGyxmRWb2c01bBO589E5FxcTXmP0SqAbkA7MB3pV2WY48He8/g1nAJ/63TeGNZ4FtAr9XnSkxtD8GqBNwN/h+cBbx7NvrGqssv0IYFqsvsOwzzkXKAAWHmN90OdiXfUFdh7Wo8ZAz0U/NcbJudgRKAj9noU3fFDU/jbG05VC5bAYzrlDwJFhMcJVDovhnJsNZJtZR5/7xqRG59ws59yO0OxsvL4XsdKQ7yFuvsMqrgb+EoU6auWcmwFsr2WTQM/FuuoL+Dw8UkNd3+GxxOpcrG+NQZ2LG11ooFDn3B5gMdC5ymYROx/jKRQ6A+vD5kuo/g8/1jZ+9o1VjeF+hJfeRzjgXTObZ97QHUHVd6aZzTezv5vZkQGT4u47NLNMYBjwatjiaH+HfgV9LtZHrM/D+gjyXPQtXs5FM8sDTgc+rbIqYudjPL2styHDYvgaLiMCfH+OmV2A9x/jOWGLz3bObTCzdsB7ZrYk9H8qsazvc7wxUPaaN0rtX4EePveNhPp8zgjgY+dc+P/JRfs79Cvoc9GXgM5Dv4I+F+sj8HPRzJrjhdJPnXO7q66uYZfjOh/j6UqhIcNixGq4DF+fY2Z9gGeAUc65bUeWO+c2hH5uBl7Hu7SLaX3Oud3Oub2h398GmphZGz/7xqrGMFdR5XI9Bt+hX0Gfi3UK8Dz0JQ7OxfoI9Fw0syZ4gfCic+61GjaJ3PkY7UaSejSmpAGrgHy+bRDpXWWbizm6MeUzv/vGsMaueD20z6qyvBmQFfb7LLxRZGNdXwe+7bQ4EFgX+j7j5jsMbdcS715vs1h+h1VqyOPYjaSBnos+6gvsPKxHjYGei35qjIdzMfSdvAA8XMs2ETsf4+b2kWvAsBjH2jegGn8FnAA8bt749mXOG7mwPfB6aFkaMMk5904A9V0O/LuZlQEHgKucd/bE03cIcCnwrnMufLD/qH+HR5jZX/CejmljZiXA3UCTsBoDPRd91BfYeViPGgM9F33WCAGfi8DZwDXAV2b2ZWjZHXjBH/HzUcNciIhIpXhqUxARkYApFEREpJJCQUREKikURESkkkJBREQqKRRERKSSQkHkOJhZjpldGfq9PGx45S+P9LkQSURx03lNJMFcCPQCXgYOOOf6BVuOSGSo85pIPZnZOcAbwE5gD9DDOdcs0KJEIkS3j0TqyTk3E5iDN9BcP4IfwVMkYhQKIsfnJGBp0EWIRJraFETqycxOAHY55w7Xsd3dQGtgp3Pu7pgUJ9JAulIQqb986hiT3sw64422uRNvKGORhKArBZH6W4I31PJC4FivYfwNcDPQlqNfciIS1/T0kUgUmNltQAu8dxosdc49EnBJIr4oFEREpJLaFEREpJJCQUREKikURESkkkJBREQqKRRERKSSQkFERCopFEREpJJCQUREKikURESk0v8HZVaEtwj8pXMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(results['time'] * gamma_emi, population_ind, 'k', label='individual')\n",
    "ax.plot(results['time'] * gamma_emi, population_coll, 'r--', label='collective')\n",
    "ax.set_ylim(0, 1)\n",
    "ax.set_xlim(0, tf * gamma_emi)\n",
    "ax.set_xlabel(r'$t\\Gamma_a$')\n",
    "ax.set_ylabel('population');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "289db1f4-286c-4a84-834a-f44a4e1dceb7",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-03T13:10:56.667541Z",
     "iopub.status.busy": "2022-02-03T13:10:56.667541Z",
     "iopub.status.idle": "2022-02-03T13:10:56.746541Z",
     "shell.execute_reply": "2022-02-03T13:10:56.746541Z",
     "shell.execute_reply.started": "2022-02-03T13:10:56.667541Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAREAAAEGCAYAAABCR6GtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAWXElEQVR4nO3dfawdxXkG8OcJOJjUpMYB/N3YDY74uEGKQg0VtKKFCBMhDKhEjvrhKFQWFVWTKm0xdVVwW1rStPmjSoh7CyiWCjhIQHC+IDaKS53EgKEIbIzjCzbxBduUEssgMGB4+8dZ4zmvz5k9e+bs7jnnPj/J8s6Zvbvj6+vXM+/MztLMICLSrQ/U3QARGWwKIiKSREFERJIoiIhIEgUREUlybN0N6AXyQwZMrbsZIkNuzytmdrL/dCiCSCOALKu7ESJDbuULrT7VcEZEkiiIiEgSBRERSaIgIiJJFEREJImCiIgkURARkSQKIiKSREFERJIoiIhIEgUREUmiICIiSRRERCSJgoiIJFEQEZEkCiIikmRINiUSOWxScPxOba2YSNQTEZEkCiIikkRBRESSKCciA26SK8d+pMvKkfg2dGo4cjbqiYhIEgUREUmi4YxIR2JDlvCfkT+vyJBlMIc36omISBIFERFJoiAiIkmUE5ESVbEEPeW6RaZmY3mP4zv8OgA41OZ4cKknIiJJFEREJEltQYTkXJI/JrmN5FaSX8w+n0ZyHckd2e8n1tVGEclXZ07kEIAvm9kTJE8A8DjJdQA+D+AhM7uZ5HIAywFcV2M7pWtlLUGP5TKKXDfWPp+vCO/pcyDTItd505UPdHh/YFDWjdTWEzGzPWb2RHb8GoBtAGYDWAxgdXbaagCX19JAEelIX+RESM4D8EkAjwCYbmZ7gEagAXBKm69ZRnIzyc3AG5W1VUSa1T7FS3IKgHsAfMnMDpDs6OvMbBTAaOMas6y8Fk4E3S7VLjJFmnfNIteKTakWGQKEQxb/T8GXTwiOp0fa4/8c+1z5zTbHg6vWngjJSWgEkDvM7N7s430kZ2b1MwG8XFf7RCRfnbMzBHAbgG1m9rWgai2ApdnxUgD3V902EelcncOZ8wD8IYCnST6ZffbXAG4GcDfJqwH8AsBV9TRPRDpRWxAxs40A2iVALqyyLRNPkRxIbOcwX+fH+EXG/OF185aDd3puLHeS97UfjpRPdXXh9yHvexBO8frv32Aug++L2RkRGVwKIiKSREFERJLUvk5E+lE4Vvc/ImFuwI/3U5Zph/mAvOuE943lEXxdXo4kNDtS9utEwu/JDlfnv3/d7gzfv9QTEZEkCiIikkTDmQmj213G/FRn6ECkrqgi08ydTh37a/of99hydT+c+VRw7J+yCFcqDMe0bRHqiYhIEgUREUmiICIiSZQTmTC6XdrudZsHKXIPz0/NvtrhdX0+wpfDfM/prm5Oc3FqWHBPa+yPfU/8930wdisrQj0REUmiICIiSRRERCSJciKSw4/3izzeH1uHUca6kNj9W90z3PJwYXPVVJf3eD04nuzvE7bvxUgdUGx5/2BQT0REkiiIiEgSDWcmpCJPksaWbfvrnNDyrNbXiQ1R/Lmxbn+szg9n/I5kQXky21YBAHYFx6/8n6s80Oa4VXk4hjAh9UREJImCiIgkURARkSTKiUgLsWnIMA8Sy4HkXcfnK4pMfYZfW2S393Oai3OCHcr8v4RxVz4YFnyeY1twnLfbe2xXtsHMl6gnIiJJFEREJImGMxNSr7rNr+Vcd1pwHNtVDCj2yuXYi7jDey5wdW6D5ZHILTa48sFwA+bYKl7/Au+8J4kHn3oiIpJEQUREkiiIiEiSWnMiJG8HcCmAl81sJPtsGoBvA5iHxmLjz5rZL+tqo3jdPn2bt+w9vG5s+tef66eZg2ncGWc1V/ml7HuD4/2u7qArN+VBnvaVgSLL+wdzSteruyfyLQCL3GfLATxkZgsAPJSVRaRP1RpEzOxhNG+YCQCLAazOjlcDuLzKNolIMXX3RFqZbmZ7ACD7/ZRWJ5FcRnIzyc3AG5U2UESOGNh1ImY2CmAUAMhZ/pVkA6zbN9VV1YbYS7FjeY4ib4bLOzdsg39T3SeOHJ7vqua58qbg+En/53zClX2HOeTXhsQMRx4k1I89kX0kZwJA9nuRVUgiUrF+DCJrASzNjpcCuL/GtohIjrqneO8CcAGAk0iOA7gBwM0A7iZ5NYBfALiqvhZWocguY0U2O065T+zrOn3CF2j+8Sqy+bJfVu5fKh4MWSZ/urnqoshlN7rys2HBD1f8EGVXcOyHV+GfrZcvOR8MtQYRM/tcm6oLK22IiHStH4czIjJAFEREJMnATvEOryLLv3txj14q8vLq2NRxbEoXAII8yB+4qguC402ubr0r7388KOxylX6bg0mRuthU9vBN6XrqiYhIEgUREUmiICIiSZQTqV3eWovYuWGuIO+l13WPzYu8lPujrjyvubgqOPYv114THPucyH7/5rqwTX6diC/HtmR8p83xxKCeiIgkURARkSQazvQ1PwSIPUFbRJEl8L1aWp93nXBp+2XNVTOal73fcM2Rl2+vvNk9wB1O4x58xN1jhyuHQ5Yi0+dFhmbDTz0REUmiICIiSRRERCSJciJ9J7aE2gvH5kV2SC9rGtL/OMV2SHNvo8PvB5dxj/7vfbypuBK3vH+85bqPNdWNLL8zKO1y9/C5jPB7NPEe4e8V9UREJImCiIgk0XCmr+U9FdvtbmVl8cOFsOx3J/uT5uI1wfEqOO7x2znXvX848lvuOrgnOPa7k/mnb0MT7+nbXlFPRESSKIiISBIFERFJopzI0Mhbit2rl2KF1/Ev044tHf9sc3Geq171TFC4z1WuaC5OCY437nTn+guH/CuMwmXvyoF0Sz0REUmiICIiSRRERCSJciIDLTaOj72NLmX8H14ntowcAP44OHTL3L/jr3t3cHxDc5V/MffGcBn8T1xlbLe3vF3kQ8qRdEo9ERFJoiAiIknaDmdIfheAtas3s8va1Ukdijzxm3LdsOyHA4uaiyPBEObWcXfurc3FqTceOT7VnbpxnftgW6R94bStb5//cQ+/Jxq+dCvWE/kXAP8KYCca3+3/yH69DmBL2Q0juYjkdpJjJJeXfT8R6U7bnoiZ/RcAkPx7M/vtoOq7JB8us1EkjwHwDTTelzgO4DGSa83smfhXikjVOsmJnEzy1w8XSM4HcHJ5TQIALAQwZmbPm9nbaLxNZHHJ9xSRLnQyxfvnADaQfD4rzwOwrLQWNcwGsDsojwM4JzyB5LIj7fjVkpszCPLG9Me3OW4l9gh/eB+XvDh7fnN5LCy4HEg4/QsAc4LjzX4p+5grh1O1sZdrF8kD5eWUlDNpJ5ZYnWlme8zsAZILAJyWVT1rZm+V3C62+KwpyWtmowBGAYCc1TYBLCLlivVEbid5IoANAB4AsNHMirycI8U4gLlBeQ6Alyq6t4gU0DYnYmaXALgAjSByBYBNJO8luYzkr5XcrscALCA5n+QHASwBsLbke4pIF6I5ETM7SPJUADea2S+zpOolAL5OcoaZLSyjUWZ2iOSfAngQwDEAbjezrWXca3jF1nf4DqX/MZgWHPu8wnlHDs//VHPVRj+qXHnk8KQbm6sud6feGm5leKerjP2YpuQqwtxQbHd8iekksToDwKMk/wfA7QC+aWa3ZD2E0pjZDwD8oMx7iEi63CleM/sbAB8HcBuAzwPYQfImNOfTRWSC6ugpXjMzknsB7EWj3zcNwD0k15nZX5XZQOlWrDs+zZVPceV57c89LZjG3eiv+4+uHOxmNsVVHbUMPnyKN2+5emzX9iLCIUze8KVXO8MNn9wgQvLPACwF8Aoak/1/aWbvkPwAGq9ZVxARmcA66YmcBOBKM3sh/NDM3iN5aTnNEpFBkRtEzOxvI3Xb2tWJyMSgnc2GVmwZ96uuHBnj+2ncXWHhK+5ktyt79K12d7tyHS/Ujr083f/TqGqd5eDRpkQikkRBRESSKIiISBLlRIaWz3OEj/T7N9f5x/2PrCPc8t8fa6oZ4V8EpYuav8z/NK0K8xx3uEqfAwnb4Nve7daOvRS2ya9jKbLeZPioJyIiSRRERCSJhjMDLZyW9H+VsSlL3x2f11ycfGQTuZGvPNdUdUOwX9TKGe6p3b1+iBIOYfYhrldDlrwdykLvtDluVe70RVf+/sM/vFFPRESSKIiISBIFERFJopzIQIlNLXp+Gvf0NscATvpIc/nc4HhGc9XKVUEe5Bq/F8AuVy6ylL1XuYMqchA+fxPbIW34qSciIkkUREQkiYYzfSe2g5bvKoe7jvmhjh/OBPVT3fBliTs1HM5scHVN76B62lX6p4NjP16+vXWsSi2yW1nYvtg0sr/O8E/5qiciIkkUREQkiYKIiCRRTqTvhGNm/3Rt7NzZru6c5uLU6UeOT2uuwl5XXhMcr/f3XBcc+zxGbOl47GVaXt6LpIosba9Cp/mS4aSeiIgkURARkSQKIiKSRDmR2sUe2fe5gOmRsn+rnSuPBMfnNlcdtVr9e+Ej/ve6yrHg2OcufA4ntuw9L58Sij2G769TxZvqilx3+NaFeLX0REheRXIryfdInu3qric5RnI7yYvraJ+IdK6unsgWAFcC+PfwQ5JnoLF+8kwAswCsJ/lxM3u3+iaKSCdqCSKH35xH0lctBrDGzN4CsJPkGICFAH5WbQvLFpsGDOv8EMULX8TtXjI1z536enC8ydWNuXLTcvYXXV1sOjPWdc97GVTsSdgiS+LrHj7Uff/q9VtidTaA3UF5HEcvgAAAkFxGcjPJzcAblTRORI5WWk+E5HoctRsFAGCFmd3f7stafGYtPoOZjQIYbdxrVstzRKR8pQURM7so/6yjjAOYG5TnAHipNy0SkTL02xTvWgB3kvwaGonVBQAerbdJvRDLgfjpy9iu7Ke6cjClOtlV+aXtoQf8BxtceUdwXGQKNbY83ed3YnkOv6VAkUfvpWp1TfFeQXIcwG8C+D7JBwHAzLai8br4Z9D4Ub9WMzMi/a2u2Zn7ANzXpu4mADdV2yIR6Va/zc6IyIDpt5zIBHFsm2OgeSm7z4m4vMKU+cGxO/VZV25aeuHfRrfDlV9Ge7Fl+TE+z+EVWR/zWofXyaN8Si+oJyIiSRRERCSJhjOliD2ZCzQPU/yQJSy7l0z5p2TDad057tQtrnzwmaDg17nHplvzpnFDsWXw/uv8E7+xH0U/fIntmBZeZ+K9SKoO6omISBIFERFJoiAiIkmUE6lFOI73uYEPt6+b4t5ctz849mmOg/6ZxG2Rk/NyOKEwH+HzOUWu473Z5jhPkXtoSrcM6omISBIFERFJoiAiIkmUE6lELFfg8wqz29f5v61wqft+v5R9lyuHay38+glf7vTHwre9SC6jyG7vsfsW2WFeyqCeiIgkURARkSQaztRiUptjz3XV9/tp23BL2nFX94Qrh8MF/9fe7dAithy9TLEXnYdDs7z2VPGiq+GnnoiIJFEQEZEkCiIikkQ5kdr5v4LY4/N++nJXcOzfVBeb6vR5GD/F2+mj9nlTur3KOcRyIEWmdIu8sU86pZ6IiCRREBGRJBrO1CLsRsdWj+6KfB3Q/DRu3kbIsV2+YkOq2Mu18nYO63ZT57zNlsMhTGwoJlVQT0REkiiIiEgSBRERSaKcSCX8OD0cx/spyjAH4b/OT6keaHOcp8jTtl6RHdSLnBv+uf3XdbssX1O6VVBPRESS1BJESH6V5LMknyJ5H8mpQd31JMdIbid5cR3tE5HO1dUTWQdgxMzOAvBzANcDAMkzACwBcCaARQBuIXlMTW0UkQ7UkhMxsx8FxU0Afi87XgxgjZm9BWAnyTEACwH8rOImliwcm/vxvt+hLOTXgoRfm7c7WbdrNoq8Ha+I2FqQXuUulAOpQj/kRL4A4IfZ8WwAu4O6cTTvF/g+kstIbia5GXij5CaKSDul9URIrgcwo0XVCjO7PztnBRr/hd5x+MtanO934ml8aDYKYLRxnVktzxGR8pUWRMzsolg9yaUALgVwoZkdDgLjAOYGp80B8FI5LSxTXjc6NsUbGz4UmeosMr1axdCiqutK1eqanVkE4DoAl5lZOBZZC2AJyeNIzgewAMCjdbRRRDpT12KzrwM4DsA6kgCwycyuMbOtJO8G8Awa/5Vea2bv1tRGEelAXbMzp0bqbgJwU4XNEZEEWvZei9juZbHdymN5jl7mGDp9oXcvKUcyqPphildEBpiCiIgkURARkSTKiVSiyHi/qryHSG+oJyIiSRRERCSJhjO16LdhSaw9/dZW6TfqiYhIEgUREUmiICIiSZQTkRaUB5HOqSciIkkUREQkiYKIiCRREBGRJAoiIpJEQUREkiiIiEgSBRERSaIgIiJJFEREJImCiIgkURARkSQKIiKSREFERJIoiIhIEgUREUmiICIiSWhmdbchGcn/BfBC3e3InATglbob4ahNnVGb4j5qZif7D4ciiPQTkpvN7Oy62xFSmzqjNnVHwxkRSaIgIiJJFER6b7TuBrSgNnVGbeqCciIikkQ9ERFJoiAiIkkURHqE5FdJPkvyKZL3kZwa1F1PcozkdpIXV9imq0huJfkeybNdXV1tWpTdc4zk8qru26Idt5N8meSW4LNpJNeR3JH9fmKF7ZlL8sckt2V/Z1+su02dUhDpnXUARszsLAA/B3A9AJA8A8ASAGcCWATgFpLHVNSmLQCuBPBw+GFdbcru8Q0AlwA4A8DnsrbU4Vto/NlDywE8ZGYLADyUlatyCMCXzex0AOcCuDb73tTZpo4oiPSImf3IzA5lxU0A5mTHiwGsMbO3zGwngDEACytq0zYz296iqq42LQQwZmbPm9nbANZkbamcmT0M4FX38WIAq7Pj1QAur7A9e8zsiez4NQDbAMyus02dUhApxxcA/DA7ng1gd1A3nn1Wp7ra1I/fi9B0M9sDNP5RAziljkaQnAfgkwAe6Zc2xRxbdwMGCcn1AGa0qFphZvdn56xAo2t6x+Eva3F+z+bVO2lTqy8rs00Rdd13YJCcAuAeAF8yswNkq29Zf1EQKcDMLorVk1wK4FIAF9qRBTjjAOYGp80B8FJVbWqj1Db14X07tY/kTDPbQ3ImgJervDnJSWgEkDvM7N5+aFMnNJzpEZKLAFwH4DIzeyOoWgtgCcnjSM4HsADAo3W0sQ/a9BiABSTnk/wgGsndtRXct1NrASzNjpcCaNeT6zk2uhy3AdhmZl/rhzZ1zMz0qwe/0EhO7gbwZPZrVVC3AsBzALYDuKTCNl2Bxv/+bwHYB+DBPmjTZ9CYvXoOjSFXXX9fdwHYA+Cd7Ht0NYCPoDEDsiP7fVqF7TkfjaHdU8HP0GfqbFOnv7TsXUSSaDgjIkkUREQkiYKIiCRREBGRJAoiIpJEQURqQfIBkrPdZxeQ/F5dbZLuKIhI5Ugej8Z6hxfrboukUxCR0pD8jWx/lckkfyXbJ2MEwAUANmTnLMr2YdmIxrYFMmC02ExKRfIfAEwGcDyAcTP7J5L/BuA7AH6KxkrM30Vjxe+3AXzIzC6tqbnSBfVEpGx/B+DTAM4G8M/ZZ+cB2AjgNAA7zWyHNf43+896migp9BSvlG0agCkAJgGYTHI6gN1m9nb2mLu6wgNOwxkpFcm1aOxgNh/ATDR27HrXzFaRnIzGw3i/Y2bPkbwLwAkazgwW9USkNCT/CMAhM7sz21/1p2g8zv4JADCzgySXAfg+yVfQGOKM1NZg6Yp6IlIZkscB+In1+QuqpRgFERFJotkZEUmiICIiSRRERCSJgoiIJFEQEZEkCiIikuT/AYZG9t5nDxkTAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.imshow(real_space_emission, origin='lower', aspect='equal', cmap='jet',\n",
    "          extent=[min(x_pos), max(x_pos), min(y_pos), max(y_pos)], vmin=0)\n",
    "\n",
    "ax.set_xlabel('x/d')\n",
    "ax.set_ylabel('y/d');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67046e9c-b4b8-473e-9d5f-13e7dcd37252",
   "metadata": {},
   "source": [
    "## Bilayer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "11473fab-8ed9-4085-9d86-5c8ab7ccd5a8",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-03T13:14:07.685260Z",
     "iopub.status.busy": "2022-02-03T13:14:07.684260Z",
     "iopub.status.idle": "2022-02-03T13:14:07.704259Z",
     "shell.execute_reply": "2022-02-03T13:14:07.704259Z",
     "shell.execute_reply.started": "2022-02-03T13:14:07.685260Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "data = np.load('divergence_energy_bilayer_shift.npy', allow_pickle=True).item()\n",
    "omega_X_int = CubicSpline(data['z_vec'], data['omega_X'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c2cf0e1-c359-44a1-a1a5-4739ad7d8ac4",
   "metadata": {},
   "source": [
    "### Individual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "3d08de55-58d6-4cc3-9282-50419de6bfcb",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-03T13:14:08.321297Z",
     "iopub.status.busy": "2022-02-03T13:14:08.320298Z",
     "iopub.status.idle": "2022-02-03T13:14:08.340300Z",
     "shell.execute_reply": "2022-02-03T13:14:08.340300Z",
     "shell.execute_reply.started": "2022-02-03T13:14:08.321297Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "pol_lat = np.array([[1, 1, 0]] / np.sqrt(2)).T\n",
    "pol_emi = [np.array([[0, 1, 0]]).T]\n",
    "\n",
    "border = 10 # Border of adiabatic lossees\n",
    "gamma_max = 10  # Maximum individual free space decay rate in the border\n",
    "\n",
    "d = 1 # Wave-lenght of transition\n",
    "lambda0 = d / 0.3 # Wave-lenght of transition\n",
    "k0 = 2 * np.pi / lambda0\n",
    "\n",
    "border = 10\n",
    "max_gamma = 10\n",
    "\n",
    "db = 0.1 * d\n",
    "shift = np.array([0.5, 0.5, db]) # Displacement vector\n",
    "\n",
    "N_x = 50\n",
    "N_y = N_x\n",
    "x_pos_A = (np.arange(0, N_x) - N_x / 2 + d / 2)\n",
    "y_pos_A = (np.arange(0, N_y) - N_y / 2 + d / 2)\n",
    "x_pos_A, y_pos_A = [x.flatten() for x in np.meshgrid(x_pos_A, y_pos_A)]\n",
    "z_pos_A = np.zeros(np.shape(x_pos_A))\n",
    "n_A = len(x_pos_A)\n",
    "\n",
    "x_pos_B = x_pos_A + shift[0]\n",
    "y_pos_B = y_pos_A + shift[1]\n",
    "z_pos_B = z_pos_A + shift[2]\n",
    "n_B = len(x_pos_B)\n",
    "\n",
    "x_pos = np.append(x_pos_A, x_pos_B)\n",
    "y_pos = np.append(y_pos_A, y_pos_B)\n",
    "z_pos = np.append(z_pos_A, z_pos_B)\n",
    "r_pos = [x_pos, y_pos, z_pos]\n",
    "\n",
    "emitter = [[0.5, 0, db / 2]]\n",
    "n_e = 1\n",
    "\n",
    "psi0 = np.zeros(len(x_pos) + 1)\n",
    "psi0[-1] = 1\n",
    "\n",
    "gamma_emi = 0.002\n",
    "omega_emi = omega_X_int(db)\n",
    "\n",
    "dt = 2e-3 / gamma_emi\n",
    "tf = dt * 2 ** 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ef60f4f7-7cf1-465a-8c8a-3ecd69612f05",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-03T13:14:08.981904Z",
     "iopub.status.busy": "2022-02-03T13:14:08.981904Z",
     "iopub.status.idle": "2022-02-03T13:15:16.320669Z",
     "shell.execute_reply": "2022-02-03T13:15:16.320669Z",
     "shell.execute_reply.started": "2022-02-03T13:14:08.981904Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "results = complete_dynamics([x_pos, y_pos, z_pos], k0, pol_lat, dt, tf, psi0, \n",
    "                            emitter=emitter, pol_emi=pol_emi * n_e, gamma_emi=gamma_emi,\n",
    "                            omega_emi=omega_emi, border=border, max_gamma=max_gamma)\n",
    "\n",
    "population_ind = np.sum(np.abs(results['psi'][:, -n_e:]) ** 2, axis=-1)\n",
    "\n",
    "population = np.abs(results['psi'][:, :n_A].reshape((-1, N_y, N_x))) ** 2\n",
    "index_max = np.argmax(np.sum(population, axis=(1, 2)))\n",
    "real_space_emission = population[index_max]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c85dd850-7480-4fee-a61d-0e59fa8ccb85",
   "metadata": {},
   "source": [
    "### Collective"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "956da10e-0cf9-49f5-94a1-41fc8c3689b5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-03T13:16:32.085011Z",
     "iopub.status.busy": "2022-02-03T13:16:32.085011Z",
     "iopub.status.idle": "2022-02-03T13:16:32.100005Z",
     "shell.execute_reply": "2022-02-03T13:16:32.100005Z",
     "shell.execute_reply.started": "2022-02-03T13:16:32.085011Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "n_plaquette = 4\n",
    "emitter = [[n_plaquette, n_plaquette + 0.5, shift[2] / 2],\n",
    "           [-n_plaquette, -n_plaquette + 0.5, shift[2] / 2]]\n",
    "n_e = 2\n",
    "\n",
    "signs = [1, 1]\n",
    "psi0 = np.zeros(n_A + n_B + n_e, dtype=complex)\n",
    "for i in range(n_e):\n",
    "    psi0[-i - 1] = signs[i] / np.sqrt(n_e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "a995c6ab-6a64-488b-9f97-a33277e2fa5c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-03T13:16:33.738350Z",
     "iopub.status.busy": "2022-02-03T13:16:33.738350Z",
     "iopub.status.idle": "2022-02-03T13:16:33.761348Z",
     "shell.execute_reply": "2022-02-03T13:16:33.761348Z",
     "shell.execute_reply.started": "2022-02-03T13:16:33.738350Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Distance between emitters: 15.264 d\n"
     ]
    }
   ],
   "source": [
    "angle_target = 57.6 / 180 * np.pi\n",
    "counter = 0\n",
    "while True:\n",
    "    previous = copy.deepcopy(emitter)\n",
    "    theta = np.arctan((emitter[0][1] - emitter[1][1]) / (emitter[0][0] - emitter[1][0]))\n",
    "    if theta < angle_target:\n",
    "        desp = -1\n",
    "    else:\n",
    "        desp = 1\n",
    "    if counter % 2 == 0:\n",
    "        emitter[1][1] += desp\n",
    "    else:\n",
    "        emitter[0][1] -= desp\n",
    "    theta_new = np.arctan((emitter[0][1] - emitter[1][1]) / (emitter[0][0] - emitter[1][0]))\n",
    "    if np.abs(theta_new - angle_target) > np.abs(theta - angle_target):\n",
    "        emitter = previous[:]\n",
    "        break\n",
    "    else:\n",
    "        counter += 1\n",
    "        pass\n",
    "\n",
    "d_e = np.sqrt((emitter[0][0] - emitter[1][0]) ** 2 + (emitter[0][1] - emitter[1][1]) ** 2)\n",
    "print('Distance between emitters: {:.3f} d'.format(d_e))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "7a20f013-8019-46e7-93a8-7d731285f4e1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-03T13:16:36.153663Z",
     "iopub.status.busy": "2022-02-03T13:16:36.152664Z",
     "iopub.status.idle": "2022-02-03T13:17:44.197360Z",
     "shell.execute_reply": "2022-02-03T13:17:44.197360Z",
     "shell.execute_reply.started": "2022-02-03T13:16:36.153663Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "results = complete_dynamics([x_pos, y_pos, z_pos], k0, pol_lat, dt, tf, psi0, \n",
    "                            emitter=emitter, pol_emi=pol_emi * n_e, gamma_emi=gamma_emi,\n",
    "                            omega_emi=omega_emi, border=border, max_gamma=max_gamma)\n",
    "\n",
    "population_coll = np.sum(np.abs(results['psi'][:, -n_e:]) ** 2, axis=-1)\n",
    "\n",
    "population = np.abs(results['psi'][:, :n_A].reshape((-1, N_y, N_x))) ** 2\n",
    "index_max = np.argmax(np.sum(population, axis=(1, 2)))\n",
    "real_space_emission = population[index_max]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "94332e3c-0204-4131-8cb0-f0264eb82849",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-03T13:17:44.198361Z",
     "iopub.status.busy": "2022-02-03T13:17:44.198361Z",
     "iopub.status.idle": "2022-02-03T13:17:44.513362Z",
     "shell.execute_reply": "2022-02-03T13:17:44.513362Z",
     "shell.execute_reply.started": "2022-02-03T13:17:44.198361Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'population')"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEOCAYAAABmVAtTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAlsUlEQVR4nO3deXxV9Z3/8dcnIWFXliSyBAhgEoUgq5CKC2pl3Gha6yiKC+g8LK1OrXajzK8uj3amYxendtRB6lir1aJTtIoCLnVBtFAWkU2BsAQCIgQIsgkJ+f7+OJd4yXqS3Jtzl/fz8biPm3PO99z74c7pvD3ne77fY845REREAFKCLkBERGKHQkFERKopFEREpJpCQUREqikURESkmkJBRESqRS0UzOwJM9tlZqvr2W5m9jszKzazlWY2Ilq1iIiIP9E8U3gSuLSB7ZcBuaHXbcD/RLEWERHxIWqh4JxbAOxtoEkR8JTzLAK6mFnPaNUjIiKNC7JPoTewLWy5NLROREQC0ibA77Y61tU554aZ3YZ3iYnuMDJ78GDatGsXzdpERBLKsmXLypxzmY21CzIUSoE+YcvZwI66GjrnZgIzAUaZuSfvuouCW2+NfoUiIgnCzEr8tAvy8tHLwE2hu5AKgf3OuU/97Hhw/froViYikqSidqZgZn8GxgEZZlYK3AukATjnZgBzgcuBYuAwMMXvZ1fs3BnpckVEhCiGgnPuuka2O+D25nz28V27mlWTiIg0LO5GNK8yY/6ZZwZdhohIQoq7UHBt2rCrvDzoMkREElLchUImMHLZsqDLEBFJSHEXCqdWVXHe5s1BlyEikpDiLhRcSgqdjh4NugwRkYQUd6FQlZrKKRUVQZchIpKQ4i4UaNOGrs5RVVkZdCUiIgknLkMBYP/WrQEXIiKSeOIuFCpOPZU0oOz48aBLERFJOHEXCqlpaTigrKws6FJERBJO3IVCGvAYULloUdCliIgknLgLhTYpKdwG2Nq1QZciIpJw4i4UUtu2BaBSM6WKiERc/IVCejoVAHv2BF2KiEjCibtQANiXkkKKQkFEJOLiMhQ+a9eOo4cOBV2GiEjCCfIZzc32kwsvZPv27XwYdCEiIgkmLs8UevXqxY4dO4IuQ0Qk4cTlmcJXy8u5ctcuKo4dIy09PehyREQSRlyeKWQDXwM+W7cu6FJERBJKXIZCu5wcAPasWRNsISIiCSYuQ6Fzbi4An69fH3AlIiKJJS5DoeuZZwJwRI/lFBGJqLgMhW6DBrEJ2LdvX9CliIgklLgMhZRu3bggO5v5XbsGXYqISEKJy1AAjVUQEYmGuBynAPDTsjIOFhcHXYaISEKJ21DoW1WF7d8fdBkiIgklbi8fVWZkkHn8OF988UXQpYiIJIy4DQXr1YtM4NPS0qBLERFJGHEbCm379iUV2K3HcoqIREzc9im0HzKEBcCBbduCLkVEJGHE7ZnCKd/8JhcAGyoqgi5FRCRhxG0odOvWjfT0dI1VEBGJoLgNBTtyhFXOccY77wRdiohIwohqKJjZpWa2zsyKzWxaHdtPNbM5ZvaRma0xsym+P7x9e/pXVtL+s88iWrOISDKLWiiYWSrwCHAZMAi4zswG1Wh2O7DWOTcUGAf8xsz8PUrNjM/btSO9vDxiNYuIJLtonimMBoqdc5ucc8eAWUBRjTYO6GxmBnQC9gKVfr/gcKdOdDx0KFL1iogkvWiGQm8g/H7R0tC6cA8DZwI7gFXAnc65qpofZGa3mdlSM1u6e/fu6vUVXbvS7fhxDhw4EPHiRUSSUTRDwepY52os/xOwAugFDAMeNrNTau3k3Ezn3Cjn3KjMzMzq9eVDh7II+PTTTyNVs4hIUotmKJQCfcKWs/HOCMJNAV5wnmJgM3CG3y8onzqVO0G3pYqIREg0Q2EJkGtm/UOdxxOBl2u02QpcDGBmpwH5wCa/X9CrVy9AoSAiEilRCwXnXCVwB/Aa8DHwvHNujZlNNbOpoWY/A84xs1XA34AfO+fK/H5H33/8g3Lg0KpVkS1eRCRJRXXuI+fcXGBujXUzwv7eAYxv7ue379aNDsChLVua+xEiIhImbkc0A9hppwFwTNNni4hERFyHAllZADiNahYRiYj4DoXQ7ampe/YEXIiISGKI71Do0IHFBQX84+BBnKs5BEJERJoqvkMBeH/KFP7v2DH2798fdCkiInEv7kOhV8+epKOxCiIikRD3oXDxY4+xGIWCiEgkxH0opGdlkYVCQUQkEqI6eK01tO/Xjw7A9m3bGm0rIiINi/tQSO/tzcZdrlHNIiItFveXj04MYDtcUhJwISIi8S/+Q2HYMP7crx8lGsAmItJi8R8Kgwbx2rhxLNdUFyIiLRb/oeAc+ZmZHNqxg4qKiqCrERGJa3Hf0UxlJT/59a/5Ati+fTs5OTlBVyQiErfi/0whLY1jnTuTBWzdujXoakRE4lr8hwLgMjPJRKEgItJSCREKbXr2JAso0W2pIiItkhChkNqjBz1TU3WmICLSQgkRCkyezKzsbIWCiEgLJUYoXHklH40YoctHIiItlBihcPAghR078llJiZ7AJiLSAokRCvPn86M//Ylehw+zb9++oKsREYlbiREKmZkAugNJRKSFEiMUQjOlagCbiEjLJFQoZKIzBRGRloj/uY8AunbFpabS24wtetiOiEizJcaZQkoK9sQTLOnbl02bNgVdjYhI3EqMUAC46SaODR7M5s2bg65ERCRuJU4orF/PRe3bs2nTJo1VEBFppsQJhfvuY/Jbb3Hw4EHKysqCrkZEJC4lTihkZdHx0CEA9SuIiDRTQoVC2pEjtEWhICLSXAkVCuCNVVBns4hI80Q1FMzsUjNbZ2bFZjatnjbjzGyFma0xs3eb/WWhUBjUvbvOFEREmsl3KJhZqpn1MrO+J16NtQceAS4DBgHXmdmgGm26AI8CX3PODQb+uan/gGqFhfDqqzBwoEJBRKSZfIWCmf0r8BnwBvBq6PVKI7uNBoqdc5ucc8eAWUBRjTbXAy8457YCOOd2NaH2k2VlweWXk5mbq1AQEWkmv9Nc3AnkO+f2NOGzewPbwpZLgTE12uQBaWb2DtAZeMg591QTvuNLVVXw6qt8pWNH/rxtGxUVFaSlpTXro0REkpXfy0fbgP1N/GyrY13NUWVtgJHAFcA/AT81s7xaH2R2m5ktNbOlu3fvrufbDK65hgu2bqWqqkqzpYqINIPfM4VNwDtm9ipw9MRK59yDDexTCvQJW84GdtTRpsw5dwg4ZGYLgKHA+vBGzrmZwEyAUaNG1T1c2Qyyssg4UfCmTQwcOLCxf5eIiITxe6awFa8/IR3vMs+JV0OWALlm1t/M0oGJwMs12rwEnGdmbcysA97lpY/9Fl9LVhanHvUyS/0KIiJN5+tMwTl3P4CZdfYW3UEf+1Sa2R3Aa0Aq8IRzbo2ZTQ1tn+Gc+9jM5gMrgSrgcefc6mb+WyAri3Y7d5Kenq5QEBFpBl+hYGYFwNNAt9ByGXCTc25NQ/s55+YCc2usm1Fj+VfAr5pQc/169MA++oicnBwNYBMRaQa/l49mAnc75/o55/oB3wd+H72ymmnaNJg3j4EDB7Jhw4agqxERiTt+O5o7OufePrHgnHvHzDpGqabmy80FIC8vjwULFuCcw6yum6BERKQufs8UNpnZT80sJ/T6f0DsXZ/Zvh1mzGBYjx4cOnSI7du3B12RiEhc8RsKt+DNNfcC8GLo7ynRKqrZNm6Eb3+bs1K8f9b69esb2UFERML5vftoH/DdKNfScj16AJDTrh0A69at46KLLgqyIhGRuNJgKJjZb51z3zOzOdQejYxz7mtRq6w5evYEoOvRo3To0EFnCiIiTdTYmcLTofdfR7uQiOjUCTp0wD77jLy8PNatWxd0RSIicaXBPgXn3LLQn8Occ++Gv4BhUa+uqcy8S0g7d5Kfn69QEBFpIr8dzTfXsW5yBOuInDffhIcfJi8vjy1btnD06NHG9xEREaDxPoXr8J550N/Mwuct6gw0ZRrt1tO/PwD5+flUVVWxceNGBg0a1MhOIiICjfcpfAB8CmQAvwlbfwBvvqLY8+67sGAB+ZddBnh3ICkURET8aTAUnHMlQAnwldYpJwIWLIB77iFv6lRAYxVERJrC7+M4C81siZkdNLNjZnbczD6PdnHNEhqrcMoXX9CjRw91NouINIHfjuaHgeuADUB74F+A/45WUS0SCgXdgSQi0nR+QwHnXDGQ6pw77pz7A3Bh9MpqgROhEBqroMtHIiL++Z0l9XDo6WkrzOyXeJ3PsTdLKpwUCmeccQZlZWWUlZWRkZHR8H4iIuL7TOFGvKen3QEcwnv28jejVVSL9O4N+/bBLbcwePBgANasafBZQCIiEuIrFJxzJc65I865z51z9zvn7g5dToo9KSnQpQuYUVBQAMDq1c1/wqeISDJpbPDaKuqYCO8E59xZEa8oEmbMgPJyev34x3Tp0kWhICLiU2N9Cle2ShWR9uabsGYNNm0aBQUFCgUREZ/8DF6LP336wPz54BwFBQXMmjVLj+YUEfHB7+C1A2b2eej1RUwPXgPIzoZDh6C8nIKCAsrLy9mxY0fQVYmIxDy/T17rHL5sZl8HRkejoIjo08d737aturN5zZo19O7dO8CiRERin+/Ba+Gcc38FYvc5l336eA/c2bu3+rZU9SuIiDTO15mCmV0VtpgCjKKBu5ICV1gIn38OZmQAPXr0UCiIiPjgd0TzhLC/K4EtQFHEq4mUGh3KugNJRMQfv30KU6JdSMTdfTdkZUHottSZM2dSVVVFSkqzrpiJiCQFv3cfDTCzOWa228x2mdlLZjYg2sW1yKJF8PrrgHemcPjwYbZs2RJsTSIiMc7vfzY/CzwP9AR6Af8H/DlaRUVEnz6wbRtA9R1IK1fG5sPiRERihd9QMOfc0865ytDrT8RyRzN4oVBaWj2AzcxYsWJF0FWJiMQ0v6HwtplNM7McM+tnZj8CXjWzbmbWLZoFNlufPvDFF7B7Nx07duSMM87gww8/DLoqEZGY5vfuo2tD79+qsf4WvDOG2OtfyM2FM8+E8nLIymL48OEsWLAg6KpERGKa37uP+ke7kIi7/HLvFTJixAieffZZdu/eTWZmZoCFiYjELr93H6WZ2XfN7C+h1x1mlhbt4iJpxIgRALqEJCLSAL99Cv8DjAQeDb1GhtbFtquugnvvBWD48OEALF++PMiKRERimt8+hbOdc0PDlt8ys48a28nMLgUewnuU5+POuf+sp93ZwCLgWufcX3zW1LitW73OZqBLly70799foSAi0gC/ZwrHzWzgiYXQwLXjDe1gZqnAI8BlwCDgOjMbVE+7B4DX/BbtW04ObN5cvThixAiFgohIA/yGwg/xbkt9x8zeAd4Cvt/IPqOBYufcJufcMWAWdc+X9K/AbGCXz1r8y8mBLVvAeUMqRowYwcaNG9m/f3/Ev0pEJBH4DYX3gceAqtDrMeDvjezTG9gWtlwaWlfNzHoD3wBmNPRBZnabmS01s6W7d+/2WTLQv793+eizz4Av+xXU2SwiUje/ofAU0B/4WejVH3i6kX3qevZlzVHQvwV+7Jxr8FKUc26mc26Uc25Uk24nLSiA8ePhyBEARo0aBcCSJUv8f4aISBLx29GcX6Oj+W0fHc2lQJ+w5Wyg5jMxRwGzQs9OzgAuN7PK0EN8Wu6CC7xXSGZmJgMGDGDRokUR+XgRkUTj90zhQzMrPLFgZmPwLik1ZAmQa2b9zSwdmAi8HN7AOdffOZfjnMsB/gJ8J2KBUI/CwkIWL14cza8QEYlbfkNhDPCBmW0xsy14/QkXmNkqM6tz6lHnXCVwB95dRR8Dzzvn1pjZVDObGoHa/SkshDvuqF4cM2YM27dvp7S0tNVKEBGJF34vH13anA93zs0F5tZYV2ensnNucnO+w0cR8Mkn1YuFhd4Jz+LFi8nOzo7KV4qIxCtfZwrOuZKGXtEuskXy8mDDhurFoUOHkp6ern4FEZE6JP6zKfPyvJHNoTuQ2rZty/Dhw9WvICJSh+QIBYDi4upVhYWFLF26lIqKioCKEhGJTYkfCsOGwa23Qtu21avGjBnDkSNHWLVqVXB1iYjEoMQPhfx8ePzxL88YgLFjxwKwcOHCoKoSEYlJiR8KAFVV3hPYQvr27Uu/fv147733gqtJRCQGJUcoXHYZXHHFSavOP/98FixYgHM1Z94QEUleyREKOTmwfv1Jq84//3x27drF+hrrRUSSWXKEQl4elJXB3r3Vq8477zwAFixYEFRVIiIxJzlCITfXew8bxJaXl0dWVpZCQUQkTHKEQn6+9x423YWZVfcriIiIJzlCYeBAuOceb8xCmPPPP5+tW7eyZcuWQMoSEYk1yREKbdrA/ffD0KEnrR43bhwAb731VgBFiYjEnuQIBYDPP4ca8x0VFBTQo0cPXn/99YCKEhGJLckTCo884j1bYf/+6lVmxvjx43nzzTepqqoKsDgRkdiQPKFQUOC9r1170urx48ezZ88ePvzwwwCKEhGJLckXCqtXn7T6q1/9KoAuIYmIkEyh0K8fdOhQKxROO+00hg4dqlAQESGZQiElBQYPhjVram0aP34877//PgcPHgygMBGR2JE8oQDwy1/CL35Ra/X48eOpqKjQQDYRSXrJFQrjxsHZZ9dafe6559K+fXvmzp3b+jWJiMSQ5AqFI0fg+edr3YHUrl07LrnkEl5++WVNpS0iSS25QqGyEiZOhL/8pdamoqIitm3bxooVK1q/LhGRGJFcodC5szeN9vLltTZdeeWVmBkvvfRSAIWJiMSG5AoFgJEjYdmyWquzsrI455xzFAoiktSSLxRGjIDSUti1q9amoqIiVqxYwdatWwMoTEQkeMkXCiNHeu91TGtRVFQEoLMFEUlayRcKhYWweTOMH19rU15eHkOGDOG5554LoDARkeAlXyi0awc5OWBW5+aJEyfy/vvv6xKSiCSl5AsFgLffhm99C+oYk3DttdcC8Pzzz7d2VSIigUvOUCguhpkzvfcaBg4cyNlnn82sWbMCKExEJFjJGQpjx3rv779f5+aJEyeybNkyNmzY0IpFiYgELzlD4YwzoEuXekPhmmuuwcx45plnWrcuEZGAJWcopKTAOefABx/UuTk7O5uLL76YJ598Uo/pFJGkEtVQMLNLzWydmRWb2bQ6tk8ys5Wh1wdmNjSa9Zzk3HMhNRWOHq1z86233kpJSQlvvfVWq5UkIhK0qIWCmaUCjwCXAYOA68xsUI1mm4ELnHNnAT8DZkarnlqmTYOVK6Ft2zo3f/3rX6dr16488cQTrVaSiEjQonmmMBoods5tcs4dA2YBReENnHMfOOf2hRYXAdlRrOdkJ8Yp1DNVdrt27Zg0aRIvvPAC+/btq7ONiEiiiWYo9Aa2hS2XhtbV51ZgXhTrqe2++7y+hXrccsstHD16lKeeeqr1ahIRCVA0Q6GuIcN1/me5mV2IFwo/rmf7bWa21MyW7t69O3IVduwIixbBp5/WuXn48OEUFhby8MMPq8NZRJJCNEOhFOgTtpwN7KjZyMzOAh4Hipxze+r6IOfcTOfcKOfcqMzMzMhVePHF3vvf/lZvkzvvvJPi4mI9qlNEkkI0Q2EJkGtm/c0sHZgIvBzewMz6Ai8ANzrn1kexlroNGwbdusGbb9bb5Jvf/Ca9e/fmd7/7XevVJSISkKiFgnOuErgDeA34GHjeObfGzKaa2dRQs3uA7sCjZrbCzJZGq546paR4ZwtvvFFvh3NaWhrf+c53eOONN1hb49nOIiKJxuLtQfWjRo1yS5dGMDteeQUWL4bp06F9+zqblJWV0adPH2688UZmzmy9u2ZFRCLFzJY550Y11i45RzSHu/JK+NnP6g0EgIyMDKZMmcKTTz7Jtm3b6m0nIhLvFAoAX3wB777bYJNp07wB2Q888EBrVCQiEgiFAsAjj8C4cVBSUm+Tvn37MnnyZH7/+9+zffv21qtNRKQVKRQAJkzw3ht5NvP06dOpqqril7/8ZSsUJSLS+hQKAHl5UFAAjTybOScnh5tuuonHHnuMkgbOKkRE4pVC4YTrr/em0t68ucFm9913H2bG9OnTW6kwEZHWo1A44brrvPdGLiH16dOHH/zgBzz77LMsXry4FQoTEWk9GqcQ7qOP4KyzvpxBtR4HDhwgNzeXgQMHsnDhQqyR9iIiQdM4heYYOrTRQADo3LkzP//5z/nggw/0yE4RSSgKhZqmT4fvf7/RZlOmTKGwsJC77rqLsrKyVihMRCT6FAo17d4NM2ZAeXmDzVJTU3n88cfZv38/d999d+vUJiISZQqFmr79bTh8GJ5+utGmgwcP5ic/+QlPP/008+fPb4XiRESiSx3NdRkzxjtTWLsWUlMbbHr06FFGjBjB3r17+eijj8jKyopubSIizaCO5pb4wQ9g/Xp48cVGm7Zt25ZZs2axb98+Jk+erCe0iUhcUyjU5aqr4O67vVHOPgwZMoQHH3yQefPm8dBDD0W5OBGR6NHlowhxznHVVVcxZ84cXn/9dS666KKgSxIRqabLR5Hw8cdw113g45KQmfHHP/6R/Px8rr76aoqLi1uhQBGRyFIoNOTDD+G3v4Vnn/XV/JRTTmHOnDmkpKQwYcIE9u7dG936REQiTKHQkIkTYcQI+Ld/825T9WHAgAHMnj2bTZs2cdlll3HgwIEoFykiEjkKhYakpMB//Rds3Qr33+97twsuuIDnn3+eZcuWUVRUxJEjR6JYpIhI5CgUGnP++XDrrfCb33gT5vlUVFTEH//4R9555x2uuOIKnTGISFxoE3QBceFXv4KsLBg4sEm7TZo0CYCbb76Ziy66iHnz5pGRkRGNCkVEIkJnCn507Qr/8R/QqRMcPdqkXSdNmsSLL77I6tWr+cpXvsLatWujVKSISMspFJqipASGDGn0sZ01TZgwgTfffJMDBw5QWFjIyy+/HKUCRURaRqHQFL16QWYmTJ7sPbqzCcaOHcvSpUvJz8+nqKiIH/7whxxt4lmHiEi0KRSaIi0N/vpXyM6GCRPgk0+atHt2djYLFizgW9/6Fr/+9a8ZPXo0q1atik6tIiLNoFBoqsxMmD8f2rSBCy+EdeuatHv79u2ZMWMGc+bMYefOnYwcOZJp06Zx6NChKBUsIuKfQqE5Bg6Et9+GkSOhR49mfcSVV17J6tWrmTRpEg888ABnnnkmzzzzDMePH49wsSIi/ikUmmvQIHjlFTj1VDh0CGbPbvJHZGZm8oc//IGFCxfSvXt3brjhBs466yxmz56tKbhFJBAKhUh46CG4+mqYMqXRx3jWZezYsSxbtoznnnuOqqoqrr76agYPHsyjjz7KwYMHI1+viEg9FAqR8KMfwfTp3iM8Bw/2OqObOCV5SkoK11xzDatXr+aZZ56hc+fO3H777WRnZ/Pd736XJUuWEG/TnItI/FEoREKbNvDv/w6LFkG3bvCNb8C99zbro1JTU7n++utZvHgxf//737n88suZOXMmo0eP5owzzuD+++/nww8/VECISFToITuRVlEB//u/cNFFkJcHK1d6z2W46irvltZmKC8vZ/bs2fzpT3/i3XffxTlHr169uOKKK7jkkks477zz6NHMDm8RSQ5+H7KjUIi2733P63PIzPT6Ha69FsaO9c4ummHXrl3MmzePV155hddee616or3c3FzOO+88Ro8ezbBhwxgyZAgdOnSI4D9EROKZQiFWVFXB3Llef8OcOXDkiDdVxsqV3vbiYujTB9q2bfJHV1RUsHz5ct577z3ee+89Fi5cWP1gn5SUFPLy8hg6dCh5eXmcfvrp5Obmcvrpp5ORkYGZRfJfKSIxLiZCwcwuBR4CUoHHnXP/WWO7hbZfDhwGJjvnljf0mXEXCuEOHoR587wH9tx8s9cZ3bMn7N3rBcXQoXD66d503eee2+SPd85RUlLCihUrql8rV66kpKTkpFtcTz31VLKzs+nduze9e/emV69e1e8ZGRl069at+pXWzEteIhJbAg8FM0sF1gOXAKXAEuA659zasDaXA/+KFwpjgIecc2Ma+ty4DoWaqqrghRdg2TJYuhRWr4adO73nQj/4oHdWkZHhDZA78crK8jqyx4/3wuXFF72xEqecAu3be2cc2dleh3dlJRw+zLGUFDaXllK8cSMbNmxg48aNlJaWsn37drZv387OnTvrHRfRuXNnunfvTrdu3ejcuTOdOnU66dWxY8eT/m7Xrh3p6em0bdu2+r2+v9PT00lNTa31SklJ0ZmMSIT5DYVoPk9hNFDsnNsUKmgWUASEzx1dBDzlvGRaZGZdzKync+7TKNYVO1JSvH6Gq6/+ct3Bg3DsmPd3RQVMneoFxc6d3lxLCxZ4Hdjjx8OOHXDDDbU/9+GH4fbbYc0aGDaMdCAfyE9Lg9RUryP8+uu9u6UmTMBlZlJlRqVzHDdj8eTJrM3OpuPy5Yx/6SUqy8s5vmcPx6uqOF5Vxf/LyGDe8eOMLi/np6HHlLqw12S8/yN/DZgeKqki9DoAXA9sAa4BvhvaJ/wzrgb2pKQw2YzJVVVeQJhhZhhwY9eufJGWxs1HjlB04jGpoRAx4IZevXApKdxSXs7Fhw4RHi/HUlK4rV8/zIzbdu/mnIMHq/cFOJCayt0DBmBm3L5jB8PC9zdjd1oa9+bmYmbcuWUL+TWmJylt145fnH46AD/auJF+NZ66t7FDBx4cMACAezZs4LQTkyKGaljbqROP5uQA8PNPPqFLZeVJ+y8/5RSe6NsXgF+tXUv7GmH+965deTY7G4DfrVpFzWh9p3t3ZvfqRXpVFb9Zs4aa5mdl8eppp3FKZSX/8fHHtbb/tWdP/paZSebRo9xTxxQvz/fuzXvdu5N95Ag/3rCh1van+/ZlSdeuDDh0iO8VF9fa/nhODitPPZVBn3/O1M2bT95oxiMDBrC+c2eG79vHzVu31tr/v3JzKenQgcI9e5hYWlpr+y/y89nVvj0X7N7NN7Zvr7X93kGD2J+ezvidO7l8585a238yZAhHUlP52o4dXLRrV63tdw8dSpUZ/7xtG2P37Dlp27GUFH501lkA3FBSwqh9+07afqBNG35aUADArZs3M2T//pO2727bln8/80wAvlNcTF6NMUzb27fnV/n5Xh3r19PX5yOEa4rmmcLVwKXOuX8JLd8IjHHO3RHW5hXgP51zC0PLfwN+7JxbWuOzbgNuCy0WAKujUnTkZABlQRfRCNUYGaqx5WK9PkiMGvs55zIb+5BoninUdf5fM4H8tME5NxOYCWBmS/2cAgVJNUaGaoyMWK8x1uuD5KoxmoPXSoE+YcvZwI5mtBERkVYSzVBYAuSaWX8zSwcmAjUfOfYycJN5CoH9SdOfICISg6J2+cg5V2lmdwCv4d2S+oRzbo2ZTQ1tnwHMxbvzqBjvltQpPj56ZpRKjiTVGBmqMTJivcZYrw+SqMa4G7wmIiLRownxRESkmkJBRESqxVQomNmlZrbOzIrNbFod283MfhfavtLMRvjdtxVrnBSqbaWZfWBmQ8O2bTGzVWa2wsyiMizbR33jzGx/qIYVZnaP331bscYfhtW32syOm1m30Lao/4ah73nCzHaZWZ1jYoI+Fn3UF+hx6LPGWDgWG6sxFo7FPmb2tpl9bGZrzOzOOtpE7nh0zsXEC68zeiMwAEgHPgIG1WhzOTAPb3xDIbDY776tWOM5QNfQ35edqDG0vAXICPg3HAe80px9W6vGGu0nAG+11m8Y9j3nAyOA1fVsD/pYbKy+wI7DJtQY6LHop8YYORZ7AiNCf3fGmz4oav+/MZbOFKqnxXDOHQNOTIsRrnpaDOfcIqCLmfX0uW+r1Oic+8A5d2L8+iK8sRetpSW/Q8z8hjVcB/w5CnU0yDm3ANjbQJNAj8XG6gv4ODxRQ2O/YX1a61hsao1BHYufutBEoc65A8DHQO8azSJ2PMZSKPQGtoUtl1L7H15fGz/7tlaN4W7FS+8THPC6mS0zb+qOoOr7ipl9ZGbzzGxwE/dtrRoxsw7ApcDssNXR/g39CvpYbIrWPg6bIshj0bdYORbNLAcYDiyusSlix2M0p7loqpZMi+FruowI8P09ZnYh3v8Yw+fAHuuc22FmWcAbZvZJ6L9UWrO+5XhzoBw0b5bavwK5PveNhKZ8zwTgfedc+H/JRfs39CvoY9GXgI5Dv4I+Fpsi8GPRzDrhhdL3nHOf19xcxy7NOh5j6UyhJdNitNZ0Gb6+x8zOAh4Hipxz1VMlOud2hN53AS/indq1an3Ouc+dcwdDf88F0swsw8++rVVjmInUOF1vhd/Qr6CPxUYFeBz6EgPHYlMEeiyaWRpeIDzjnHuhjiaROx6j3UnShM6UNsAmoD9fdogMrtHmCk7uTPmH331bsca+eCO0z6mxviPQOezvD/BmkW3t+nrw5aDF0cDW0O8ZM79hqN2peNd6O7bmb1ijhhzq7yQN9Fj0UV9gx2ETagz0WPRTYywci6Hf5Cngtw20idjxGDOXj1wLpsWob9+AarwH6A48at4c+ZXOm7nwNODF0Lo2wLPOufkB1Hc18G0zqwSOABOdd/TE0m8I8A3gdedc+AMLov4bnmBmf8a7OybDzEqBe4G0sBoDPRZ91BfYcdiEGgM9Fn3WCAEfi8BY4EZglZmtCK2bjhf8ET8eNc2FiIhUi6U+BRERCZhCQUREqikURESkmkJBRESqKRRERKSaQkFERKopFESawcyyzeza0N/Hw6ZXXnFizIVIPIqZwWsiceZiYBDwHHDEOTcs2HJEIkOD10SayMzOBV4CyoEDQK5zrmOgRYlEiC4fiTSRc24hsARvorlhBD+Dp0jEKBREmicfWBd0ESKRpj4FkSYys+7AfudcRSPt7gW6AeXOuXtbpTiRFtKZgkjT9aeROenNrDfebJvleFMZi8QFnSmINN0neFMtrwbqewzjz4A7gUxOfsiJSEzT3UciUWBm3wdOwXumwTrn3H8HXJKILwoFERGppj4FERGpplAQEZFqCgUREammUBARkWoKBRERqaZQEBGRagoFERGpplAQEZFqCgUREan2/wHhT+UaxovtfwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(results['time'] * gamma_emi, population_ind, 'k', label='individual')\n",
    "ax.plot(results['time'] * gamma_emi, population_coll, 'r--', label='collective')\n",
    "ax.set_ylim(0, 1)\n",
    "ax.set_xlim(0, tf * gamma_emi)\n",
    "ax.set_xlabel(r'$t\\Gamma_a$')\n",
    "ax.set_ylabel('population');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "9512cda3-c7d8-40d9-8720-32773a5badfd",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-03T13:18:14.914833Z",
     "iopub.status.busy": "2022-02-03T13:18:14.914833Z",
     "iopub.status.idle": "2022-02-03T13:18:15.008342Z",
     "shell.execute_reply": "2022-02-03T13:18:15.008342Z",
     "shell.execute_reply.started": "2022-02-03T13:18:14.914833Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAREAAAEGCAYAAABCR6GtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAUv0lEQVR4nO3df5CdVX3H8c9Hy8+iAoIkJGhSCSPxx4w2RRy1kxatwWEMMGKhThtHbIYpVLFWIcVRaEu1qIxatDQFBmZAI1WRoCUYGCOiAqbUiUDABBNhIQQptWKpQOTbP+6N7P3es7s3e/be5+7m/ZrJZM95nnOfbzabb8793vOcxxEhAJis5zQdAIDpjSQCoApJBEAVkgiAKiQRAFV+q+kApoK9b0j7Nx0GMMNtezQiDs69MyKJtBLI8qaDAGa4835a6uXtDIAqJBEAVUgiAKqQRABUIYkAqEISAVCFJAKgCkkEQBWSCIAqJBEAVUgiAKqQRABUIYkAqEISAVCFJAKgygzZTwQYtD0mMebpKY9iGDATAVClsSRi+zDb37K90fZdtt/X7j/Q9lrbm9q/H9BUjAAm1uRMZIekD0TEkZKOlnS67YWSzpZ0U0QskHRTuw1gSDWWRCJiW0Tc0f76cUkbJc2RtFTSFe3TrpB0fCMBAujJUBRWbc+T9GpJt0k6JCK2Sa1EY/tFY4xZrt/szvyCgcSJ3UWpaJr/qTy/h9f5RWpTWO0L2/tJ+oqkMyMif9fHFBErI2JRRCyS9u1fgADG1WgSsb2HWgnkqoj4art7u+3Z7eOzJT3SVHwAJtbkpzOWdKmkjRFx4ahDqyUta3+9TNK1g44NQO+arIm8XtKfSvqR7R+2+/5G0sclXW37VEn3SzqpmfCw+8g1kOcVzpmX2oek9tbCmO2TjGd6aSyJRMQtkjzG4WMGGQuAyWu8sApgeiOJAKhCEgFQZSgWmwH9M5mFY0d1D9nvdzvbv8wnbJhELDNj8RkzEQBVSCIAqpBEAFShJoIZJtcdSj/iuQbyms7mG1L9Q5LentpnRup4sHCdvCBtZi4+YyYCoApJBEAVkgiAKtREMMOVNg86srP5std2NA/+zv1dI372jRennitSO9VVJJVvypt5mIkAqEISAVCFJAKgCkkEQBUKq5hh8o906WEBizubX+9cOPZOXdU14tOnrUg981K7tBvaxtTORd7SDXg7ejhnuDATAVCFJAKgCkkEQBVqIphh8g14r+o+5ZLO5j+99M872h/8nwu6x4ykTYdmLe5sP3x9IZbDUzvfpLdPYcz0e2oeMxEAVUgiAKqQRABUoSaCaay0CXNaF3L0gq4z/uHU93e0n9KeHe1f7X9g98senfq25hNKN+CtTu1cA8lrQsbqG27MRABUIYkAqEISAVCFJAKgCoVVTGOlxVqdu5TN/f6mrjNWbPl0R9ur087tBxVedlFq35oXga0rDHostXsprA7/4rKMmQiAKiQRAFVIIgCqUBPBNFZYFPb2zsVlD2xx1ynXzv+jzo4r0wkfL1zqPbm2si61S0/Ay4vh/i+1p9/CshJmIgCqNJpEbF9m+xHbd47qO9D2Wtub2r8f0GSMAMbX9EzkcklLUt/Zkm6KiAWSbmq3AQypRmsiEXGz7Xmpe6me3Un3CrXefJ41uKgwvHKN4U1dZ3zo387r7Dix+1WOP+uGzo5z0wnHbSlc+7rUzpsHldas5Hin3xqQXjQ9Eyk5JCK2SVL799J23bK93PZ62+ulJwYaIIBnDWMS6UlErIyIRRGxSNq36XCA3dYwJpHttmdLUvv3RxqOB8A4hjGJrJa0rP31MknXNhgLgAk0Wli1/UW1iqgH2R6R9FG1lvpcbftUSfdLOqm5CDFc0lPmLprbdcaJPrejfUW8o/tljk3tNf+VOq4uXDsXTvPT7CbzT6mXxWal3dsmMtgCbtOfzpwyxqFjBhoIgEkbxrczAKYRkgiAKtyAhyFReu+ffzw7F5c985Hum+t8Qmf7aK8qvO7dqf3d1J5XGJMXl2WTuZmu9M8v9/VSE8k39g0WMxEAVUgiAKqQRABUoSaCIfah1O58739e3gdZ0nnfSZsua23hdXOd4ZAJjkvd/1RyDaRUM8m1iry2JF9X6l7jkV/j8cKYZjETAVCFJAKgCkkEQBWSCIAqFFbRkFxUfE33KWekxWQPd+4bc96sXESVdNGXJn5dbUztvGt8LwvHetnZrJdC6kSvmwurw7c7GjMRAFVIIgCqkEQAVKEmgimQF2eV6gNzUjvVC9a8tnvIkrRQ7NNv7myfWVpIdnhqby+ck2sgk7mBbYI/T/F1cz2jFFteTNbLU/OarZMwEwFQhSQCoApJBEAVkgiAKhRWMQm5kJofUviqwpgjO5unpULkksJdsLNSIfWT+YQHC9eZzMKxfE6pMDxvgnNKxc1cOM3twm3IExZSWWwGYIYhiQCoQhIBUIWaCJJc73he4Zx8I9m81M4LvqSuxViP5uOXdw85/r2d7YtXpxNeWbjO1vGvW5RrOqUxuWaTr1Oqb+S+XO8oLXIbvprHRJiJAKhCEgFQhSQCoAo1kd1KXttQ2tU81wfmFs5Jaz661mbc0T3k8Bd2tr+8Ll0m1T8k6eLckWsxpZpC/jOV1onkm+dyHaK0/mRrauc1H8N3Y9ygMBMBUIUkAqAKSQRAFZIIgCpjFlZtXyepsJ12S0S8rS8RYQrlQupEi8Qk6SWpXVp4lYuKeSHWUd1DuhaXpeLrGxZ3j1mVdy7L8Zd+fPNCt1LBM8efC6mlR2LmHcd2j6JpL8b7dGbnPZMnSpol6cp2+xR1l6oB7KbGTCIR8W1Jsv13EfH7ow5dZ/vmfgdme4mkz0h6rqRLIuLj/b4mgF3XS03kYNu/s7Nhe76kg/sXkmT7uZI+J+lYSQslnWJ7YT+vCWByells9n5J62z/pN2eJ2l53yJqOUrS5oj4iSTZXiVpqaS7+3zdaaSXG+VyPSMvssqLxKTuBVylG8vytVNN4eQF3UNWXdjZ3v+v0vFS+W1eauc6RP7zSN31jJ8Wzsk1kFwjmczu77uv8QqrsyNiW0Sssb1A0svah+6JiCf7HNccSQ+Mao9I6nimgO3l+k0ye0GfwwEwlvFmIpfZPkDSOklrJN0SEb3sNTcVXOjr+K8qIlZKWilJ9qFjfooEoL/GK6wea3tvSYslnSDpk7bvVyuhrImI+/sY14ikw0a150p6qI/XAzBJ49ZEIuJXtg+XdG5E/He7qHqspItsz4qIwoKAKfEDSQva13tQ0smS/qRP1xpCpRvj8l9VL2s+cp1kglqGpO66Q6lukiekaUPlWwpD9JrO5kHp8M9vn3hMV33mR4UxI6n9SOGc/GekBlKjl8LqLEm32/5PSZdJ+ueI+LztPfsVVETssH2GpBvU+oj3soi4q1/XAzB5E37EGxEflnSEpEslvUvSJtvnq3yP+JSJiH+PiCMi4qURcX4/rwVg8nq6dyYiQtLD7V871JrjfsX2BX2MDcA0MOHbGdvvlbRMrbsfLpH0wYh42vZzJG2S9KH+hghgmPVSEzlI0okR0bFqJyKesX1cf8KaSUpF0nxjXG6Xipm5L7d72aE8L6oq/fU/kdqlp8Glguxp6fDFn+0e8va0c9mX8wmlHeLzDmlbUzv/eSRulBu8CZNIRHxknGMbpzYcANMN+4kAqEISAVCF3d6r5ZpH/paWagq5npEXjpXqA/k6eYFUaeFYPqcUSzYvtQsLuua+q7PdVd94V+7oXlzWVa/ZXIhlQ2rnhWPUO4YBMxEAVUgiAKqQRABUoSayS0prPvJNbnm9RmmzoFwDya+b1zpI3TWEfE7przL3Teavu7Dxz1+n9plpQ+U16YY8SVpyferIf55eajrUQIYRMxEAVUgiAKqQRABUIYkAqEJhdVwTLSSTJr55rnRjXN4ZLBcVSwXEHEt+3V52Q8s7gx1ZGJNuejvjj7tPyTuXnZEKqUtuK7zu1tQu3TyH6YiZCIAqJBEAVUgiAKpQE9klpbpDlhdRlXYSz315TOlGuXztvHlQSa6t5DrK1u4hB3XWQD56UfcjgM7T/6aeHO+NPcSCmYKZCIAqJBEAVUgiAKqQRABUobBaLd9Nmxd0leRCZC6alp6bnouvWS9j8o5p3YvN4o2dhdT4dvernvfY1tRz3QSxYSZjJgKgCkkEQBWSCIAq1ER2SWnhWJYXVZUWqJXqF7t67YlurpOk93Q2l7ywo3n59d03192W1pYdfVF0v+wZ+Ql3vXxfMFMxEwFQhSQCoApJBEAVaiLj6tdNYxO9bul4rjvkm+ne0T3k6501kPwUumUnXt015Kw4t7PDI4VYSrvRY3fFTARAFZIIgCqNJBHbJ9m+y/YzthelYytsb7Z9r+23NBEfgN41VRO5U9KJkv5ldKfthZJOlvRySYdKutH2ERHx68GHCKAXjSSRiNgoSXbXrllLJa2KiCclbbG9WdJRkr4/2AgHrZfHRebHWS7ubF4yv3vIuZ3Nr/0gTew+1T3kgpM+mnouL8TCLmV41rDVROZIemBUe0TFh8FKtpfbXm97fW9bBQLoh77NRGzfKGlW4dA5EXHtWMMKfYV111JErJS0snWtQ4vnAOi/viWRiHjTJIaNSDpsVHuupIemJiIA/TBsi81WS/qC7QvVKqwukHR7syHVKu3cnr/tuSaSNw+SWnXoUc5IN/adXRjy4c7m0i3f7GgfNv/H3WO+vCl19LLJEnZnTX3Ee4LtEUmvk/QN2zdIUkTcJelqSXdLWiPpdD6ZAYZbU5/OXCPpmjGOnS/p/MFGBGCyhu3TGQDTzLDVRGaAvAlRLzWRfM7i7iGz0uuuT8cf7R7yqff9RUf7Yzqzoz3yugWF2K5MbTYcwviYiQCoQhIBUIUkAqAKSQRAFQqrU66Xb2kuVubbg+4ojDm2s3nrho7m3j+f2zViTz3V0f7ApZ9Pr5EXlknSI6nNzXYYHzMRAFVIIgCqkEQAVKEmMuXy0+1+UTgnL0jbnNqFG6AfXtfZnru4o3naCz7WNeQv7/vXzo735DM25A5RA8GuYiYCoApJBEAVkgiAKtREplyuKfTyNLsDU/vBwpi0OdDFnc2r9M7uIcfl3SbXpXZeEyJ113SA8TETAVCFJAKgCkkEQBWSCIAqFFaHQl58VroBLz2ZLq1P+9kbX9w95J7bUsfG1C4thAN2DTMRAFVIIgCqkEQAVKEm0ohcA9me2qVnmKeFY+npdvrlfxTG5NpKroGUFpZxAx52DTMRAFVIIgCqkEQAVKEm0oiJ6g7PL/SlNR+/zHWUrYUxj+/idYFdx0wEQBWSCIAqJBEAVUgiAKpQWB0K+6T2Y4Vzch8LxzAcmIkAqEISAVClkSRi+xO277G9wfY1tvcfdWyF7c2277X9libiA9C7pmoiayWtiIgdtv9R0gpJZ9leKOlkSS+XdKikG20fERG/bijOAcn1jFJNJKPegeHQyEwkIr4ZETv/5dwqaW7766WSVkXEkxGxRa39u45qIkYAvRmGmsi7JV3f/nqOpAdGHRtR+b542V5ue73t9dITfQ4RwFj69nbG9o2SZhUOnRMR17bPOUetufxVO4cVzo/S60fESkkrW69zaPEcAP3XtyQSEYVH2z/L9jJJx0k6JiJ2JoERSYeNOm2upIf6EyGAqdDUpzNLJJ0l6W0RMfq9yGpJJ9vey/Z8SQsk3d5EjIP19CR+AcOhqU9nLpK0l6S1tiXp1og4LSLusn21pLvVeptz+sz/ZAaY3hpJIhFx+DjHzpd0/gDDAVBhGD6dATCNkUQAVCGJAKhCEgFQhSQCoApJBEAVkgiAKiQRAFVIIgCqkEQAVCGJAKhCEgFQhSQCoApJBEAVP7up2PRl+2eSfjqFL3mQpEen8PX6jXj7i3hbXhIRB+fOGZFEpprt9RGxqOk4ekW8/UW84+PtDIAqJBEAVUgiZSubDmAXEW9/Ee84qIkAqMJMBEAVkgiAKiSRUWx/wvY9tjfYvsb2/qOOrbC92fa9tt/SYJi/Yfsk23fZfsb2onRsGONd0o5ns+2zm46nxPZlth+xfeeovgNtr7W9qf37AU3GuJPtw2x/y/bG9s/B+9r9A42XJNJpraRXRMSrJP1Y0gpJsr1Q0smSXi5piaTP235uY1E+605JJ0q6eXTnMMbbvv7nJB0raaGkU9pxDpvL1fqejXa2pJsiYoGkm9rtYbBD0gci4khJR0s6vf09HWi8JJFRIuKbEbGj3bxVrWcBS9JSSasi4smI2CJps6SjmohxtIjYGBH3Fg4NY7xHSdocET+JiKckrVIrzqESETdLeix1L5V0RfvrKyQdP8iYxhIR2yLijvbXj0vaKGmOBhwvSWRs75Z0ffvrOZIeGHVspN03rIYx3mGMqVeHRMQ2qfUPV9KLGo6ni+15kl4t6TYNON6mnsXbGNs3SppVOHRORFzbPucctaaKV+0cVjh/IJ+N9xJvaVihr+nP8ocxphnB9n6SviLpzIj4Rfv51gOz2yWRiHjTeMdtL5N0nKRj4tlFNCOSDht12lxJD/Unwk4TxTuGxuIdxzDG1KvttmdHxDbbsyU90nRAO9neQ60EclVEfLXdPdB4eTsziu0lks6S9LaIeGLUodWSTra9l+35khZIur2JGHs0jPH+QNIC2/Nt76lW4Xd1wzH1arWkZe2vl0kaawY4UG5NOS6VtDEiLhx1aLDxRgS/2r/UKkA+IOmH7V8Xjzp2jqT7JN0r6dimY23HdIJa/8M/KWm7pBuGPN63qvWp131qvR1rPKZCjF+UtE3S0+3v7amSXqjWpxyb2r8f2HSc7VjfoNZbwg2jfmbfOuh4WfYOoApvZwBUIYkAqEISAVCFJAKgCkkEQBWSCBphe43tOalvse2vNxUTJockgoGzvY9aaxcebDoW1COJoG9s/157b5a9bf92e8+LV0haLGld+5wl7T1cblFrWwNMMyw2Q1/Z/ntJe0vaR9JIRHzM9mclfU3S99RaVfmHaq0W/pKkfSPiuIbCxSQwE0G//a2kN0taJOmCdt/rJd0i6WWStkTEpmj9b3ZlMyGixm53Fy8G7kBJ+0naQ9Letg+R9EBEPNW+ZZ2p8DTH2xn0le3Vau1iNl/SbLV23/p1RFxse2+1bsj7g4i4z/YXJT2PtzPTCzMR9I3tP5O0IyK+0N5j9Xtq3Zr+SkmKiF/ZXi7pG7YfVestzisaCxiTwkwEA2N7L0nfjWn0cGxMjCQCoAqfzgCoQhIBUIUkAqAKSQRAFZIIgCokEQBV/h/PQYA9NCgO6AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.imshow(real_space_emission, origin='lower', aspect='equal', cmap='jet',\n",
    "          extent=[min(x_pos), max(x_pos), min(y_pos), max(y_pos)], vmin=0)\n",
    "\n",
    "ax.set_xlabel('x/d')\n",
    "ax.set_ylabel('y/d');"
   ]
  }
 ],
 "metadata": {
  "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}