FigSM5 - purcell_cluster.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "77f428f0-d341-4845-b567-be46fdf1032a",
   "metadata": {},
   "source": [
    "# Code for Fig. SM5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f37c11c6-07e3-4250-8ad6-71411b3148b5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-07T16:05:59.080884Z",
     "iopub.status.busy": "2022-02-07T16:05:59.080340Z",
     "iopub.status.idle": "2022-02-07T16:06:00.146657Z",
     "shell.execute_reply": "2022-02-07T16:06:00.144730Z",
     "shell.execute_reply.started": "2022-02-07T16:05:59.080818Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 1\n",
    "%aimport general_functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "09d458ec-8034-4adf-a118-4f59dd367cae",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-07T16:06:00.161456Z",
     "iopub.status.busy": "2022-02-07T16:06:00.161297Z",
     "iopub.status.idle": "2022-02-07T16:06:00.185665Z",
     "shell.execute_reply": "2022-02-07T16:06:00.185296Z",
     "shell.execute_reply.started": "2022-02-07T16:06:00.161429Z"
    },
    "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",
    "\n",
    "from general_functions import (compute_energies, purcell_factor, complete_dynamics,\n",
    "                               g0_fun, compute_shift)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d29ff92b-ac02-4eb4-b766-8840ec2a2066",
   "metadata": {},
   "source": [
    "## Single impurity (a)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d5f3661-8a29-4657-864b-8de747506c88",
   "metadata": {},
   "source": [
    "### Analytical"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "ed532075-7e85-42d6-ac95-30dfbf0e7de2",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-07T17:50:50.073762Z",
     "iopub.status.busy": "2022-02-07T17:50:50.073195Z",
     "iopub.status.idle": "2022-02-07T17:50:50.090992Z",
     "shell.execute_reply": "2022-02-07T17:50:50.090549Z",
     "shell.execute_reply.started": "2022-02-07T17:50:50.073692Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "pol_lat = np.array([[0, 1, 0]]).T\n",
    "pol_emi = np.array([[1, 1, 0]]).T / np.sqrt(2)\n",
    "\n",
    "d = 1\n",
    "lambda0 = d / 0.3\n",
    "k0 = 2 * np.pi / lambda0\n",
    "aho = 0.1 * d\n",
    "\n",
    "n = 2 ** 10 + 1\n",
    "nmax = 9\n",
    "\n",
    "b1 = 2 * np.pi / d * np.array([1, 0])\n",
    "b2 = 2 * np.pi / d * np.array([0, 1])\n",
    "A = d ** 2\n",
    "\n",
    "G0 = g0_fun(k0, aho)\n",
    "\n",
    "k_x = np.linspace(-np.pi, np.pi, n) / d\n",
    "k_X, k_Y = np.meshgrid(k_x, k_x)\n",
    "\n",
    "k = np.vstack((k_X.flatten(), k_Y.flatten())).T\n",
    "\n",
    "x = np.arange(0, n) - n / 2 + d / 2\n",
    "y = np.arange(0, n) - n / 2 + d / 2\n",
    "x, y = np.meshgrid(x, y)\n",
    "z = np.zeros((n, n))\n",
    "r_lat = [x, y, z]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "ff877fe3-f4e6-4073-bd47-4aeeb373a7c6",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-07T17:50:50.709977Z",
     "iopub.status.busy": "2022-02-07T17:50:50.709441Z",
     "iopub.status.idle": "2022-02-07T17:53:41.284099Z",
     "shell.execute_reply": "2022-02-07T17:53:41.283687Z",
     "shell.execute_reply.started": "2022-02-07T17:50:50.709910Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "J, Gamma = compute_energies(b1, b2, A, nmax, k0, aho, k, G0=G0, polarizations=[pol_lat])\n",
    "omega_k = (J - 0.5j * Gamma).reshape(n, n)\n",
    "\n",
    "omega_q = np.real(omega_k).reshape(n, n)[n // 2, -1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "db6df6ed-e699-4bec-bcb5-399762e06c15",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-07T18:30:56.965064Z",
     "iopub.status.busy": "2022-02-07T18:30:56.964840Z",
     "iopub.status.idle": "2022-02-07T18:30:56.969976Z",
     "shell.execute_reply": "2022-02-07T18:30:56.969464Z",
     "shell.execute_reply.started": "2022-02-07T18:30:56.965040Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "n_pos = 100  # Number of points\n",
    "\n",
    "# The emitter is located in the diagonal (y = -x), and moves from the center\n",
    "# of a plaquette up a distance of 0.5d from the center.\n",
    "\n",
    "x_vector = -np.linspace(0.146, 0.5, n_pos) * d\n",
    "y_vector = np.linspace(0.146, 0.5, n_pos) * d\n",
    "z_vector = np.zeros(n_pos)\n",
    "\n",
    "# Distance to the center of the plaquette. [Note that there is a array dipole at (0, 0)]\n",
    "x_axis_vector_analytical = np.sqrt((x_vector[::-1] - x_vector[0]) ** 2 +\n",
    "                                    (y_vector[::-1] - y_vector[0]) ** 2)\n",
    "\n",
    "r_emi_vector = np.array([x_vector, y_vector, z_vector])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "1c8d16d0-1bff-422a-87dd-c45af367c906",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-07T18:30:58.008749Z",
     "iopub.status.busy": "2022-02-07T18:30:58.008356Z",
     "iopub.status.idle": "2022-02-07T18:31:54.735069Z",
     "shell.execute_reply": "2022-02-07T18:31:54.734555Z",
     "shell.execute_reply.started": "2022-02-07T18:30:58.008731Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0e8dec4d78bf46cb83b97b239283e23a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Computing Purcell:   0%|          | 0/100 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pbar = tqdm(range(n_pos), desc='Computing Purcell')\n",
    "results = Parallel(n_jobs=-1, backend='multiprocessing')(\n",
    "    delayed(purcell_factor)(r_lat, r_emi_vector[:, j], [k_X, k_Y],\n",
    "                            omega_k, omega_q, k0, pol_lat, pol_emi) for j in pbar)\n",
    "\n",
    "results = np.array(results)\n",
    "optical_depth = results[:, 1] / results[:, 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4fc198c-ea97-401b-9832-a48785077bc1",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Numerical"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "f619725d-9816-4a54-b0fd-b8260aa17d2e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-07T17:57:45.681161Z",
     "iopub.status.busy": "2022-02-07T17:57:45.680982Z",
     "iopub.status.idle": "2022-02-07T17:57:45.686057Z",
     "shell.execute_reply": "2022-02-07T17:57:45.685749Z",
     "shell.execute_reply.started": "2022-02-07T17:57:45.681133Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "border = 15\n",
    "max_gamma = 15\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "1f2fef8e-f51d-4509-8610-e60a904739f0",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-07T18:02:57.202257Z",
     "iopub.status.busy": "2022-02-07T18:02:57.201736Z",
     "iopub.status.idle": "2022-02-07T18:29:45.208482Z",
     "shell.execute_reply": "2022-02-07T18:29:45.207977Z",
     "shell.execute_reply.started": "2022-02-07T18:02:57.202191Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b2d06a59d6f844a886c701de16f6334b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Computing dynamics:   0%|          | 0/20 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_pos = 20\n",
    "purcell_numerical = np.zeros(n_pos)\n",
    "\n",
    "x_vector = np.linspace(0, 0.5, n_pos) * d / np.sqrt(2)\n",
    "y_vector = -np.linspace(0, 0.5, n_pos) * d / np.sqrt(2)\n",
    "z_vector = np.zeros(n_pos)\n",
    "\n",
    "x_axis_vector_numerical = np.sqrt(x_vector ** 2 + y_vector ** 2)\n",
    "\n",
    "r_emi_vector = np.array([x_vector, y_vector, z_vector])\n",
    "\n",
    "pbar = tqdm(total=n_pos, desc='Computing dynamics')\n",
    "for i in range(n_pos):\n",
    "    emitter = [r_emi_vector[:, i]]\n",
    "    \n",
    "    results = complete_dynamics(r_pos, k0, pol_lat, dt, tf, psi0, N_x, N_y, emitter, [pol_emi], gamma_emi,\n",
    "                      omega_emi, border, max_gamma=max_gamma)\n",
    "    \n",
    "    purcell_numerical[i] = results['Q_n'][-1]\n",
    "    \n",
    "    pbar.update()    \n",
    "pbar.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "b6f572c0-9fcb-49db-ae55-e2e4945051c1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-07T18:31:54.874378Z",
     "iopub.status.busy": "2022-02-07T18:31:54.874212Z",
     "iopub.status.idle": "2022-02-07T18:31:54.973998Z",
     "shell.execute_reply": "2022-02-07T18:31:54.973553Z",
     "shell.execute_reply.started": "2022-02-07T18:31:54.874361Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEJCAYAAACZjSCSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAApp0lEQVR4nO3deXgUVdr38e+dhLCFXXQUhgFHRzGySWSIjopP2BRHFEUdl3HXkVVQB7kEdUReH0UFQdYRB1EQRAdFRVTQ6MxDWBJABNlcQFEIARXClvW8f3SIERIISbqql9/nuuqq7k5VnTsnnb77VJ06x5xziIhI9IrxOwAREfGXEoGISJRTIhARiXJKBCIiUU6JQEQkysX5HUBpTjjhBNe8eXO/wxARCRsZGRk7nXONK7JvSCaC5s2bk56e7ncYIiJhw8y2VHRfnRoSEYlySgQiIlFOiUBEJMqF5DWC0uTl5bF161YOHjzodyhSxWrUqEHTpk2pVq2a36GIRKWwSQRbt26lTp06NG/eHDPzOxypIs45du3axdatW2nRooXf4YhEpbA5NXTw4EEaNWqkJBBhzIxGjRqppSfio7BJBICSQITS31XEX2GVCERE5EibNm2q1P6eJAIzG2hma8xsrZnd60WZ4eC8884LehnNmzdn586dZf588+bNzJw5s/h5eno6AwYMqFBZnTp1Kr4R8FjlikjFOedYsWIFQ4YMoWXLlvzhD3+o1PGCngjM7GzgTqAD0Aa4zMxOD3a54WDx4sV+h3BEIkhKSmLs2LE+RiQiZfn2228ZOXIkZ511Fu3bt2f06NH89re/rfT/rBctgpbAEufcfudcPvAJcKUH5Vapffv20aNHD9q0acPZZ5/N7NmzycjI4KKLLqJ9+/Z069aNbdu2AYFvxoMGDeLCCy+kZcuWLF++nF69enH66aczbNiw4mMmJCSUWtbbb7/NH//4R9q1a0fnzp3JzMwE4NFHH+W2226jU6dOnHrqqb/6419xxRW0b9+exMREpkyZcsQxhw8fznPPPVf8/KGHHmLs2LE8+OCD/Oc//6Ft27aMHj2a1NRULrvsMgD27t3LrbfeSqtWrWjdujVvvPEGAPfccw9JSUkkJibyyCOPHLXeyipXRMqnsLCQDz74gJ49e9KiRQuGDRvGiSeeyOTJk9m+fTsffPAB/fv3r1QZXnQfXQOMNLNGwAHgUuCIgYTM7C7gLoBmzZod9YD33nsvq1atqtIg27Zty5gxY8r8+YIFCzjllFN49913Adi9ezeXXHIJb731Fo0bN2b27Nk89NBDvPjiiwDEx8fz6aef8txzz9GzZ08yMjJo2LAhv//97xk0aBCNGjUqs6w//elPLFmyBDPjhRde4KmnnuKZZ54BYP369Xz88cdkZ2dzxhlncM8991CtWjVefPFFGjZsyIEDBzj33HO56qqrflXG7bffTq9evRg4cCCFhYXMmjWLZcuW0bp1a55++mneeecdAFJTU4v3GTFiBPXq1ePzzz8H4KeffgJg5MiRNGzYkIKCAlJSUli9ejWtW7cu9Xcpq1wRObrc3FxefvllnnrqKTZu3Ejjxo158MEHufPOO6nqQTmDngicc+vM7EngQ2Av8BmQX8p2U4ApAElJSSE3kXKrVq24//77GTJkCJdddhkNGjRgzZo1dOnSBYCCggJOPvnk4u0vv/zy4v0SExOLf3bqqafy3XffHTURbN26lWuvvZZt27aRm5v7q/71PXr0oHr16lSvXp0TTzyRzMxMmjZtytixY5k7dy4A3333HZs2bfpVGc2bN6dRo0asXLmSzMxM2rVrd9QYABYuXMisWbOKnzdo0ACA1157jSlTppCfn8+2bdv44osvykwEFSlXJJodPHiQqVOn8uSTT/Ldd99xzjnnMGPGDK666iqqV68elDI9uaHMOTcVmApgZv8P2FqZ4x3tm3uw/OEPfyAjI4P58+czdOhQunTpQmJiImlpaaVuf+gPFhMT86s/XkxMDPn5v86DDz30UHFLY9WqVfTv35/Bgwdz+eWXk5qayqOPPnrEcQFiY2PJz88nNTWVhQsXkpaWRq1atejUqVOp/fLvuOMOpk2bxvbt27ntttuO+Ts7547o2vnNN9/w9NNPs3z5cho0aMAtt9xyzHsAjrdckWjknGPWrFkMHTqULVu2cP755zNlyhS6desW9C7WXvUaOrFo3QzoBbzqRblV6YcffqBWrVrceOON3H///SxdupSsrKziRJCXl8fatWsrdOyRI0eyatWq4tNdu3fvpkmTJgC89NJLx9x/9+7dNGjQgFq1arF+/XqWLFlS6nZXXnklCxYsYPny5XTr1g2AOnXqkJ2dXer2Xbt25fnnny9+/tNPP7Fnzx5q165NvXr1yMzM5L333jtmfKWVKyK/WLx4MR07duT666+nQYMGLFy4kP/85z90797dk/tsvBpi4o2iawR5QF/n3E8elVtlPv/8cx544AFiYmKoVq0aEydOJC4ujgEDBrB7927y8/O59957SUxMrHRZjz76KL1796ZJkyZ07NiRb7755qjbd+/enUmTJtG6dWvOOOMMOnbsWOp28fHxXHzxxdSvX5/Y2FgAWrduTVxcHG3atOGWW26hXbt2xdsPGzaMvn37cvbZZxMbG8sjjzxCr169aNeuHYmJiZx66qmcf/75x/x9SitXRGDPnj0MGTKESZMmccopp/Cvf/2Lm266yfv/E+dcyC3t27d3h/viiy+OeE2OT0FBgWvTpo3buHFjyJWrv69Em3nz5rkmTZq4mJgYN3jwYLd3795KHQ9IdxX8zNWdxVHiiy++4LTTTiMlJYXTT/fuNg6/yhUJVfv37+eOO+7g8ssvp2HDhqSlpfHMM89Qu3Zt32IKm9FHpXLOOussvv7666gpVyQUrVmzhmuvvZZ169YxdOhQHn30UeLj4/0OS4lARMQL06ZN45577qFevXq8//77xV3PQ4FODYmIBFFBQQEPPPAAt956K+effz6rVq0KqSQAahGIiARNdnY2N954I/PmzaNv376MGTOGuLjQ+9gNvYhERCLA9u3b6d69O2vWrOH555+nb9++fodUJp0aChPz5s3jf//3fyu0r4aEFvHWt99+y4UXXsimTZt45513QjoJgFoEYSE/P5/LL7+8ePwiEQldX375JSkpKezevZsPP/zQk3lHKiuyWwRpafDEE4F1JW3evJmWLVty5513kpiYSNeuXTlw4MCvJmPZuXNn8aiA06ZN44orruDPf/4zLVq04Pnnn+fZZ5+lXbt2dOzYkR9//BGAr776iu7du9O+fXsuuOAC1q9fD8Att9zC4MGDufjiixkyZAjTpk2jX79+AGRmZnLllVfSpk0b2rRpUzyvwbGGohaR4Fq/fj0XXngh+/bt46OPPgqLJACRnAjS0iAlBYYPD6yrIBls2rSJvn37snbtWurXr188Pn9Z1qxZw8yZM1m2bBkPPfQQtWrVYuXKlSQnJzN9+nQA7rrrLsaNG0dGRgZPP/00ffr0Kd5/48aNLFy4sHgI6kMGDBjARRddxGeffcaKFSuKh7V48cUXycjIID09nbFjx7Jr165K/84iUj5btmyhc+fOFBYW8umnn3LOOef4HVK5Re6podRUyM2FgoLAOjUVkpMrdcgWLVrQtm1bANq3b8/mzZuPuv3FF19MnTp1qFOnDvXq1ePPf/4zEBiaevXq1ezdu5fFixfTu3fv4n1ycnKKH/fu3bvUMUc++uij4kQSGxtLvXr1AI45FLWIBEdmZiZdunRh3759fPLJJ5x11ll+h3RcIjcRdOoE8fGBJBAfH3heSYcPAX3gwAHi4uIoLCwEOGI45sOHny45NHV+fj6FhYXUr1+/zEl2jueW8/IORS0iVWv37t10796d77//ng8//LDMuTlCWeSeGkpOhkWLYMSIwLqSrYGyNG/enIyMDABef/3149q3bt26tGjRgjlz5gCBAQA/++yzY+6XkpLCxIkTgcDNKnv27Cn3UNQiUnVyc3Pp2bMna9eu5d///nfYXBM4XOQmAgh8+A8dGrQkAHD//fczceJEzjvvvAp10ZwxYwZTp06lTZs2JCYm8tZbbx1zn+eee46PP/6YVq1a0b59e9auXUv37t3Jz8+ndevWDB8+vMyhqEWkajjn6NOnD5988gnTpk0L67k2LDB6aWhJSkpyh3riHLJu3TpatmzpU0QSbPr7SrgZM2YMgwYNYtiwYYwYMcLvcDCzDOdcUkX2jewWgYhIECxYsID77ruPXr168Y9//MPvcCrNq6kqB5nZWjNbY2avmlkNL8oVEalqmzZt4tprr6VVq1ZMnz6dmJjw/z4d9N/AzJoAA4Ak59zZQCxwXUWOFYqnsaTy9HeVcHHw4EGuueYa4uLimDdvnq+TyVQlr1JZHFDTzOKAWsAPx3uAGjVqsGvXLn1oRBjnHLt27aJGDTUSJfTdd999rFq1iunTp9OsWTO/w6kyQb+PwDn3vZk9DXwLHAA+cM59cPh2ZnYXcBdQagU3bdqUrVu3kpWVFeSIxWs1atSgadOmfochclSvv/46EyZM4L777qNHjx5+h1Olgt5ryMwaAG8A1wI/A3OA151zr5S1T2m9hkRE/PL111/Trl07WrZsyaeffhoS00seLtR7DXUGvnHOZTnn8oB/A+F514WIRJ2CggJuvPFGzIxZs2aFZBKoLC+GmPgW6GhmtQicGkoB9HVfRMLC6NGjSUtL4+WXXy4eXTjSBL1F4JxbCrwOrAA+LypTYySLSMhbt24dw4YN44orruCGG27wO5yg8WTQOefcI8AjXpQlIlIV8vPzueWWW0hISGDSpEmYmd8hBU3kjj4qIlIJTz/9NMuWLWPWrFmcdNJJfocTVOF/S5yISBXbuHEjjzzyCFdffTXXXHON3+EEnRKBiEgJzjn69u1LzZo1GTduXESfEjpEp4ZEREqYPXs2Cxcu5Pnnn+c3v/mN3+F4ImyGoRYRCbbdu3dz5pln0rRpU5YsWVLqVLGhqjI3lKlFICJSZPjw4WRmZvL222+HVRKoLF0jEBEBVq5cyfjx4+nTpw9JSRX6Yh22lAhEJOo55xg0aBANGzbk8ccf9zscz+nUkIhEvbfeeotPPvmECRMmUL9+fb/D8ZwuFotIVMvNzSUxMZH4+Hg+++wz4uLC8/uxLhaLiFTQhAkT+PLLL3nvvffCNglUlq4RiEjU+vHHH3nsscfo1q0b3bt39zsc3ygRiEjUeuyxx9i9ezdPP/2036H4SolARKLSli1bmDBhArfffjtnn3223+H4SolARKLSP/7xD2JiYnj44Yf9DsV3SgQiEnU2bNjASy+9RJ8+fWjatKnf4fgu6InAzM4ws1Ullj1mdm+wyxURKcvDDz9MrVq1GDp0qN+hhISg95Vyzm0A2gKYWSzwPTA32OWKiJRm5cqVvPbaawwbNozGjRv7HU5I8PrUUArwlXNui8fliogAgYHlGjRowH333ed3KCHD60RwHfBqaT8ws7vMLN3M0rOysjwOS0SiwdKlS3n33Xf5+9//HpVDSZTFsyEmzCwe+AFIdM5lHm1bDTEhIsFw2WWXsWTJEjZv3kxCQoLf4VSpygwx4WWL4BJgxbGSgIhIMKxYsYJ3332XQYMGRVwSqCwvE8FfKOO0kIhIsD3++OPUr1+ffv36+R1KyPEkEZhZLaAL8G8vyhMRKWnNmjXMnTuXAQMGUK9ePb/DCTmeDLXnnNsPNPKiLBGRw40cOZKEhAQGDhzodyghSXcWi0hE27BhA7Nnz6Zfv340bNjQ73BCkhKBiES0J598kho1ajB48GC/QwlZSgQiErG+//57XnnlFe644w7dRXwUSgQiErGee+45CgsL1Ro4BiUCEYlIP//8M5MmTeKaa66hefPmfocT0pQIRCQiTZ48mezsbB544AG/Qwl5SgQiEnFycnIYM2YMXbp0oV27dn6HE/I8uY9ARMRLr7zyCtu3b+fll1/2O5SwoBaBiESUwsJCRo0aRbt27UhJSfE7nLCgFoGIRJQFCxawYcMGZsyYgZn5HU5YUItARCLK6NGjOeWUU+jdu7ffoYQNJQIRiRiff/45CxcupF+/flSrVs3vcMKGEoGIRIwxY8ZQs2ZN7r77br9DCStKBCISEXbs2MGMGTO4+eabNbjccVIiEJGIMGnSJHJycjTUdAUoEYhI2MvJyWHChAlccsklnHnmmX6HE3a8mqGsvpm9bmbrzWydmSV7Ua6IRIdZs2aRmZnJoEGDvCs0LQ2eeCKwDnNe3UfwHLDAOXe1mcUDtTwqV0QinHOOcePG0bJlSzp37uxNoWlpkJICubkQHw+LFkFy+H6/DXqLwMzqAhcCUwGcc7nOuZ+DXa6IRIclS5aQkZFBv379vLuBLDU1kAQKCgLr1FRvyg0SL04NnQpkAf8ys5Vm9oKZ1T58IzO7y8zSzSw9KyvLg7BEJBKMGzeOunXr8te//tW7Qjt1CrQEYmMD606dvCs7CLxIBHHAOcBE51w7YB/w4OEbOeemOOeSnHNJmklIRMpj27ZtzJkzh9tuu42EhATvCk5ODpwOGjEi7E8LgTfXCLYCW51zS4uev04piUBE5HhNnjyZgoIC+vbt633hyclhnwAOCXqLwDm3HfjOzM4oeikF+CLY5YpIZMvNzWXy5MlccsklnHbaaX6HE9a86jXUH5hR1GPoa+BWj8oVkQj1xhtvsH37dvr37+93KGHPk0TgnFsFJHlRlohEh+eff57TTz+drl27Hv/OaWmBnj6dOkXM6Z3K0HwEIhJ2Vq1axeLFixk9ejQxMcd5hjvC7gGoChpiQkTCzvjx46lZsya33HLL8e8cYfcAVAUlAhEJKz///DMzZszghhtuoH79+sd/gAi7B6Aq6NSQiISVadOmceDAgYp3GT10D4CuERQz55zfMRwhKSnJpaen+x2GiISYwsJCzjzzTE444QQWL17sdzhVo4ouXJtZhnOuQp1y1CIQkbCxaNEiNm3axCOPPOJ3KFUjRC5c6xqBiISN8ePH07hxY66++mq/Q6kaIXLhWolARMLCt99+y9tvv83tt99O9erV/Q6naoTIhWudGhKRsDBlyhScc5E1MX2IXLhWIhCRkJebm8sLL7xAjx49aN68ud/hVK0QGLxOp4ZEJOS9+eabZGZmcs899/gdSkRSIhCRkDdhwgRatGhBt27d/A4lIikRiEhI++KLL/jkk0+4++67iY2N9TuciKREICIhbeLEicTHx3Pbbbf5HUrEOmYiMLObzWynmf1oZtPNrI4XgYmI7N27l+nTp9O7d280hW3wlKdFMBzoApwJbAH+X1AjEhEpMnPmTPbs2UOfPn38DiWilaf76B7n3Mqix8PNbOlRty6FmW0GsoECIL+i42GISPRwzjFx4kRat25NsgaGC6ryJIKTzewuYB2wHqhWwbIuds7trOC+IhJllixZwqpVq5g0aRJm5nc4Ea08ieARoDVwA9AKSDCz+cBnwGrn3KtBjE9EotTEiROpU6cON9xwg9+hRLxjJgLn3JSSz82sKYHE0Aq4FChPInDAB2bmgMmHH7PouHcBdwE0a9asHIcUkUi1c+dOZs+ezZ133klCQoLf4US84x5iwjm3FdgKzD+O3c53zv1gZicCH5rZeufcp4cddwowBQLzERxvXCISOf71r3+Rm5urO4k94sl9BM65H4rWO4C5QAcvyhWR8FNYWMjEiRO56KKLSExM9DucqBD0RGBmtQ/de2BmtYGuwJpglysi4en999/nm2++UWvAQ16MPnoSMLfoqn8cMNM5t8CDckUkDE2cOJGTTjqJK6+80u9QwkdaGk3gNxXdPeiJwDn3NdAm2OWISPjbvHkz77zzDg899BDx8fFHblBF8/tGlKLpLn8DTSp6CM1HICIhY9KkScTExJQ++UyIzO8bcg5Nd1kJGnRORELCwYMHeeGFF+jZsydNmzY9coMQmd835BRNd+kC3fQrRIlARELC7Nmz2bVrF3379i19gxCZ3zfkFE13mQk/VPQQ5lzoddlPSkpy6enpfochIh7q0KEDe/fuZe3atWUPKaFrBGUys4yKjuOmawQi4rtly5axfPlyxo0bd/RxhUJgft9IpFNDIuK78ePHk5CQwF//+le/Q4lKSgQi4qtD4wrddNNN1K1b1+9wopISgYj46p///Cc5OTllXySWoFMiEBHf5OXlMWHCBFJSUjSukI90sVhEfPPmm2+ydetWxo8f73coUU0tAhHxzdixY2nRogU9evTwO5SopkQgIr5YsWIF//3vf+nfvz+xsbF+hxPVlAhExBfjxo2jdu3a3HrrrX6HEvWUCETEczt27GDmzJncfPPN1K9f3+9wop4SgYh4bsqUKeTm5tKvXz+/QxGUCETEYzk5OYwfP56uXbvSsmVLv8MRPOw+amaxQDrwvXPuMq/KFZHQ8uqrr7J9+3Zeeuklv0ORIl62CAYC6zwsT0RCjHOOZ599lrPPPpsuXbr4HY4U8SQRmFlToAfwghfliUhoWrRoEZ9//jmDBw8++iij4imvWgRjgL8DhWVtYGZ3mVm6maVnZWV5FJaIeOmZZ57hpJNO4vrrr/c7FCkh6InAzC4DdjjnMo62nXNuinMuyTmX1Lhx42CHJSIeW7t2LQsWLKBfv35Ur17d73CkBC9aBOcDl5vZZmAW8D9m9ooH5YpICBk9ejQ1a9bkb3/7m9+hyGGCngicc0Odc02dc82B64CPnHM3BrtcEQkd27dv55VXXuHmm2/mhBNO8DscOYzuIxCRoBs9ejR5eXncd999focipfB0GGrnXCqQ6mWZIuKvn3/+mYkTJ9K7d29OO+00v8ORUqhFICJBNWHCBLKzs3nwwQf9DkXKoEQgIkFz4MABxowZQ/fu3Wnbtq3f4UgZlAhEJGhefPFFsrKyGDp0qN+hyFEoEYhIUOTl5TFq1CjOO+88LrjgAr/DkaPQnMUiEhSzZs1iy5YtjBs3TsNJhDi1CESkyuXn5zNixAjatGkTmI84LQ2eeCKwlpCjFoGIVLmZM2eyadMm5s6dS8zSpZCSArm5EB8PixZBcrLfIUoJahGISJU61Bpo27YtPXv2hNTUQBIoKAisU1P9DlEOoxaBiFSpGTNm8OWXX/Lmm28Grg106hRoCRxqEXTq5HeIchhzzvkdwxGSkpJcenq632GIyHHKz8/nzDPPpE6dOqxYseKXi8RpaYGWQKdOOi0UJGaW4ZxLqsi+ahGISJWZMWMGX3311S+tgUOSk5UAQphaBCJSJXJycmjZsiX169cnIyNDXUY9phaBiPhu8uTJfPPNNyxYsEBJIMyo15CIVNqePXsYMWIEKSkpdO3a1e9w5DgpEYhIpY0aNYqdO3fy5JNPqjUQhpQIRKRStm3bxrPPPst1111H+/bt/Q5HKsCLyetrmNkyM/vMzNaa2T+CXaaIeOexxx4jLy+Pxx9/3O9QpIK8uFicA/yPc26vmVUD/mtm7znnlnhQtogE0bp16/jnP/9Jnz59+P3vf+93OFJBQU8ELtA/dW/R02pFS+j1WRWR4+KcY8CAAdSpU4fhw4f7HY5UgifXCMws1sxWATuAD51zS0vZ5i4zSzez9KysLC/CEpFKmDt3LgsXLmTEiBE0btzY73CkEjy9oczM6gNzgf7OuTVlbacbykRC2/79+znrrLOoW7cuK1asIC5OtyT5LWxuKHPO/WxmqUB3oMxEICKh7amnnmLLli2kpqYqCUQAL3oNNS5qCWBmNYHOwPpglysiwbF582aefPJJrrvuOi666CK/w5Eq4EUqPxl4ycxiCSSe15xz73hQrohUMecc/fr1IyYmhlGjRvkdjlQRL3oNrQbaBbscEQm+mTNn8u677zJmzBiaNm3qdzhSRXRnsYiUy44dOxg4cCAdO3akX79+focjVUiJQETKZeDAgWRnZzN16lRiY2P9DkeqkBKBiBzTvHnzmDVrFsOHD+ess87yOxypYpqYRkSOaufOnbRu3ZrGjRuzfPly4uPj/Q5JShE29xGISHhxznHbbbexa9cu5s+fryQQoXRqSETKNH78eLLefpuPunSh7YEDfocjQaJEICKlWr16NXMGD+bjmBjOW7AAUlIgLc3vsCQIlAhE5Aj79+/nuuuuo1v16lQ3wwoKIDcXUlP9Dk2CQIlARH7FOccdd9zB+vXr6fz441h8PMTGQnw8dOrkd3gSBLpYLCK/MmrUKF599VVGjhxJh4EDoUOHQEugUydITvY7PAkCdR8VkWLz58/nsssu45prruHVV1/VRPRhpDLdR3VqSEQA2LBhA3/5y19o06YNU6dOVRKIIkoEIsK2bdu49NJLqV69Om+99Ra1a9f2OyTxkK4RiES5n376ia5du5KZmclHH31Es2bN/A5JPKZEIBLF9u3bR48ePdi4cSPz58+nQ4cOfockPlAiEIlSBw8e5KqrrmLp0qXMmTOHlJQUv0MSn3gxVeVvzexjM1tnZmvNbGCwyxSRo8vOzqZHjx68//77TJkyhV69evkdkvjIixZBPnCfc26FmdUBMszsQ+fcFx6ULSKH+fHHH7n00ktJT09n+vTp3HTTTX6HJD7zYqrKbcC2osfZZrYOaAIoEYh4bNu2bXTr1o0NGzbw+uuvc8UVV/gdkoQAT7uPmllzAvMXLy3lZ3eZWbqZpWdlZXkZlkhUWL58Oeeeey5fffUV7777rpKAFPMsEZhZAvAGcK9zbs/hP3fOTXHOJTnnkho3buxVWCJR4ZVXXuGCCy4gLi6O//u//6Nz585+hyQhxJNEYGbVCCSBGc65f3tRpohATk4OgwcP5qabbiI5OZnly5fTtm1bv8OSEONFryEDpgLrnHPPBrs8EQlYvXo15557LqNHj6Z///588MEHqLUtpfGiRXA+cBPwP2a2qmi51INyRaJSfn4+o0aNol/79lz7zTf8d9Qoxo4dS7Vq1fwOTUKUF72G/gto9CoRD6SmpjJgwABqf/45qbGxxB84gD38MJx/voaQljJp0DmRCLBlyxauvfZaLr74Ynbv3s2U668nHjSzmJSLEoFIGPv666+58847Oe2005g3bx6PPvoo69ato1W/fppZTMpNYw2JhBnnHMuWLWP8+PHMnDmTuLg47r77boYMGcJvf/vbwEbJybBokWYWk3JRIhAJEz/99BOvvfYakyZNYtWqVdSuXZv+/fvzwAMPcMoppxy5Q3KyEoCUixKBSAjbtWsX8+bNY86cOSxcuJC8vDzatGnDxIkTuf7666lbt67fIUoECMlEkJmZyZgxYzh8PuVDz51zv3pc2msll9JeKywsPOK1sl4/2mtlHaes7Q+PpeTvcDwOTSNoZsVLTExM8bq0JTY2tnhdcomLiytex8XFUa1ateJ1ySU+Pr54fWipXr068fHx1KhRg+rVq/9qHR8fr+kOj9PPP//MsmXLWLRoEQsXLmTlypU452jevDn33nsvvXv3JikpSfUqVSokJ683M8+CKvlBWvLDtKznJV8v7WcxMTFHHPfw7Q7/+aHn5VVWcjuUcEomqcLCQgoKCn71+PAlmO+BGjVqUKNGDWrWrFm81KpVq3hdu3ZtatWqVbzUrl37iCUhIaF4ffjjcO0bX1BQwObNm1m3bh3r1q1j9erVLFu2jI0bNwJQrVo1kpOT6dy5M5dccgnt27fXh78cVWUmrw/JRNCuXTv38ccfA0d+QB7+QXq018r68D30wSwBhxJEfn5+8ZKXl1fqkpub+6slJyeneJ2Tk8PBgweL1wcOHODgwYPFjw8t+/fvZ//+/cWP9+3bV/zavn37jisxxcfHk5CQQJ06dYqTQ8mlZCIp+fhQ4imZoGrWrEn16tWLl/j4+OJWUlxc3BHvmZL1VvL3zM7OZs+ePezevZsff/yRrKwsduzYQWZmJt9++y1btmzhu+++Iy8vr/hYJ598Mh06dODcc8+lQ4cOJCcnk5CQUGV/Y4l8EZcIkpKSXHp6ut9hiA+cc+Tk5LBv3z727t3Lvn37ih9nZ2eX+vzQ47179xYvJbfNzs5m//79VRZjTEwMhYWFx7VPfHw8J554Is2aNeN3v/sdzZo14/TTT6dly5a0bNmSBg0a/LJxWpp6+8hxq0wiCMlrBBGlqv6pQ+k4QfygMrPi00mNGjWqsuMWFhZy8ODB4uRQsgVSssVyqGVzqKVTUFBAXl4e+fn5vzr9VvJaS1xcXPHprxo1alCnTh3q1q1L3bp1adCgASeeeCJ169YtXys0LQ1SUgI3gcXHB7qAKhlIkIVmIti+PfAP4fcHVWWPU1X/1KF0nKr8oPLwm29MTEzx6aCQHngtNTVQtyXvCFYikCALzUTw/feBDxs/P6iq4jhV9U8dSsepqljCNKEEXadOgfo4VC+6I1g8ELpDTFR0fJTSPqgqoiqOc+ifurK3+YfScaoqlqr6Ox1KKMOHB9ZpaRU/zhNPVHz/qjrOoTuCR4zQaSHxTGi2CKDyH1SV/UZVFcepqtv8Q+k4VRVLVf2dqqKFEkqtSNAdweK50EwETZrAnDn+flBV5XGq4p86lI5TVccIlYQSSqfeRHyg7qMS/iLtor5IBYT0fQRm9iJwGbDDOXd2efZRIhDPhUpPM5EKCvVEcCGwF5iuRCAiEhyVSQRB7zXknPsU+DHY5YiISMWETPdRM7vLzNLNLD0rK8vvcEREokbIJALn3BTnXJJzLimk7/wUEYkwIZMIRETEH0oEIiJRLuiJwMxeBdKAM8xsq5ndHuwyRUSk/ELyhjIzywY2+B1HiDgB2Ol3ECFA9fAL1cUvVBe/OMM5V6ciO4bmEBOwoaL9YSONmaWrLlQPJakufqG6+IWZVfjmK10jEBGJckoEIiJRLlQTwRS/AwghqosA1cMvVBe/UF38osJ1EZIXi0VExDuh2iIQERGPKBGIiEQ53xKBmXU3sw1m9qWZPVjKz83Mxhb9fLWZneNHnF4oR12caWZpZpZjZvf7EaNXylEXNxS9H1ab2WIza+NHnF4oR130LKqHVUUDNv7Jjzi9cKy6KLHduWZWYGZXexmfl8rxvuhkZruL3herzOzhYx7UOef5AsQCXwGnAvHAZ8BZh21zKfAeYEBHYKkfsYZIXZwInAuMBO73O2af6+I8oEHR40ui/H2RwC/X+VoD6/2O26+6KLHdR8B84Gq/4/bxfdEJeOd4jutXi6AD8KVz7mvnXC4wC+h52DY9CUxm45xzS4D6Znay14F64Jh14Zzb4ZxbDuT5EaCHylMXi51zPxU9XQI09ThGr5SnLva6ov98oDYQqT0/yvN5AdAfeAPY4WVwHitvXRwXvxJBE+C7Es+3Fr12vNtEgmj5PcvjeOvidgKtxkhUrrowsyvNbD3wLnCbR7F57Zh1YWZNgCuBSR7G5Yfy/o8km9lnZvaemSUe66B+JQIr5bXDv82UZ5tIEC2/Z3mUuy7M7GICiWBIUCPyT7nqwjk31zl3JnAFMCLYQfmkPHUxBhjinCsIfji+Kk9drAB+55xrA4wD3jzWQf1KBFuB35Z43hT4oQLbRIJo+T3Lo1x1YWatgReAns65XR7F5rXjel+4wJSwvzezE4IdmA/KUxdJwCwz2wxcDUwwsys8ic5bx6wL59we59zeosfzgWrHel/4lQiWA6ebWQsziweuA+Ydts084K9FvYc6Arudc9u8DtQD5amLaHHMujCzZsC/gZuccxt9iNEr5amL08zMih6fQ+DiYSQmxmPWhXOuhXOuuXOuOfA60Mc596bnkQZfed4XvynxvuhA4HP+qO8LX0Yfdc7lm1k/4H0CV8FfdM6tNbO/Ff18EoEr/5cCXwL7gVv9iDXYylMXZvYbIB2oCxSa2b0Eegrs8SvuYCjn++JhoBGBb3wA+S4CR58sZ11cReDLUh5wALi2xMXjiFHOuogK5ayLq4F7zCyfwPviumO9LzTEhIhIlNOdxSIiUU6JQEQkyikRiIhEOSUCEZEop0QgIhLllAhERKKcEoGISJRTIhApJzObbGbnl/L6VDPr4UdMIlVBiUDkGMwstujhHwkMfX24tgTGhRcJS74MMSES6sxsDoHhftsBi8zsDWCjc67AzP4AvAjUK1r/xjm31b9oRSpHLQKR0rUC9jrnLnbOPU5gNrQFZhYHvAIMds61Ak4H1vsYp0ilqUUgchgzqwE0BB4r8XI3AgMf9gLWOeeWFb2+lsDAXiJhSy0CkSMlEpgLOR/AzGoB9Z1zPxCYGzijxLbt0fUBCXNKBCJHagWsLvH8YuDjose7gLMBzKw98BdglZfBiVQ1nRoSOVIrYFmJ55cQmOwE4GVgvpmtAjYAPwPrvAxOpKppPgKRYzCzFcAfnXN5fsciEgxKBCIiUU7XCEREopwSgYhIlFMiEBGJckoEIiJRTolARCTKKRGIiEQ5JQIRkSj3/wFhxP9u7MaRWgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(x_axis_vector_analytical, optical_depth, 'k', label='semi-analytically')\n",
    "ax.plot(x_axis_vector_numerical, purcell_numerical, 'r.', label='numerical')\n",
    "\n",
    "ax.set_xlim(0, 0.5)\n",
    "ax.legend()\n",
    "ax.set_xlabel(r'$r/d$')\n",
    "ax.set_ylabel(r'$P$');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26039d82-5fa3-4d23-861b-e99655f32d76",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Cluster (b-c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "7bb35cdf-27d5-4949-96f7-f7db22183ba2",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-07T23:48:26.929768Z",
     "iopub.status.busy": "2022-02-07T23:48:26.929538Z",
     "iopub.status.idle": "2022-02-07T23:48:26.936582Z",
     "shell.execute_reply": "2022-02-07T23:48:26.936164Z",
     "shell.execute_reply.started": "2022-02-07T23:48:26.929747Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "pol_lat = np.array([[0, 1, 0]]).T\n",
    "pol_emi = [np.array([[1, 1, 0]]).T / np.sqrt(2)] * 2\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 = 15 # Border of adiabatic lossees\n",
    "max_gamma = 15  # Maximum individual free space decay rate in the border\n",
    "r_circles = np.arange(border, border + 5)\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 = [x.flatten() for x in np.meshgrid(x_pos, y_pos)]\n",
    "z_pos = np.zeros(np.shape(x_pos))\n",
    "\n",
    "gamma_emi = 0.002\n",
    "omega_emi = -0.999  # For d / lambda_0 = 0.3\n",
    "\n",
    "n_e = len(pol_emi)\n",
    "signs = [1, -1]  # The emitters are initialized in the anti-symmetric state\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",
    "dt = 2e-3 / gamma_emi\n",
    "tf = dt * 2 ** 12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc31224e-ecd2-4b93-9d70-3788680e82a0",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-07T23:52:25.464326Z",
     "iopub.status.busy": "2022-02-07T23:52:25.464191Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5e165b5cde944860be32a6409854e4ff",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Computing dynamics:   0%|          | 0/100 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_dc = 100\n",
    "dc_vector = np.linspace(0.01, 0.5, n_dc) * d\n",
    "\n",
    "purcell_numerical = np.zeros(n_dc)\n",
    "real_space_emission = np.zeros((n_dc, N_y, N_x))\n",
    "chi_1D = np.zeros(n_dc)\n",
    "\n",
    "pbar = tqdm(total=n_dc, desc='Computing dynamics')\n",
    "for i in range(n_dc - 1, -1, -1):\n",
    "    emitter = [[dc_vector[i] / np.sqrt(2), -dc_vector[i] / np.sqrt(2), 0],\n",
    "               [-dc_vector[i] / np.sqrt(2), dc_vector[i] / np.sqrt(2), 0]]\n",
    "    \n",
    "    shift = compute_shift(emitter[0], emitter[1], pol_emi[0], pol_emi[1], k0, gamma_emi, psi0[-2:])\n",
    "    \n",
    "    results = complete_dynamics([x_pos, y_pos, z_pos], k0, pol_lat, dt, tf, psi0, N_x, N_y, emitter, \n",
    "                                pol_emi, gamma_emi,omega_emi - np.real(shift), border,\n",
    "                                max_gamma=max_gamma, r_circles=r_circles,\n",
    "                                verify=True)\n",
    "    \n",
    "    purcell_numerical[i] = results['Q_n'][-1]\n",
    "    chi_1D[i] = np.average(results['chirality'])\n",
    "    \n",
    "    dt = results['time'][1]\n",
    "    tf = results['time'][-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[i] = population[index_max]\n",
    "    \n",
    "    pbar.update()    \n",
    "pbar.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "401473de-cea8-4591-85e2-3ac3fa5a264e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-08T07:12:44.947556Z",
     "iopub.status.busy": "2022-02-08T07:12:44.946995Z",
     "iopub.status.idle": "2022-02-08T07:12:45.085584Z",
     "shell.execute_reply": "2022-02-08T07:12:45.085069Z",
     "shell.execute_reply.started": "2022-02-08T07:12:44.947488Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEKCAYAAAA8QgPpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAoNElEQVR4nO3dd3xV5eHH8c+TBSEJe4iBMBL2qsiQWZAhMmRFBRU3tpaqrVq1w1GttT9rq1gVUUvBCcioMgRFNioQhjLC3iEskZCEEZI8vz9ACgqS3Nx7nzu+79eLl9zknnO+HC98c87znHOMtRYREZHiinAdQEREgpMKREREPKICERERj6hARETEIyoQERHxiApEREQ8EuU6gDdVrlzZ1q5d23UMEZGgsmLFikPW2irFXS6kCqR27dqkpaW5jiEiElSMMTs9WU6nsERExCMqEBER8YgKREREPKICERERj6hARETEIyoQERHxSEhN4xUR8afjp46z++hu9ufsJ68gj+jIaDrX6gzArqxdRJpIqsZVJToy2nFS31CBiIhcgrWW9QfXs+XwFvo37A/AjZNuZOK6iee9r26Fumy9fysAt//3dubtmEdURBQpFVNoVLkR7Wq043cdfuf3/L6iAhERuYAjJ44wa8ssZm+dzewts8nMySTSRJL9+2xio2O5rv51tKjWgppla3JZ/GWUjipNXEzc2eX/0OkP3NjkRnZm7ST9UDrrD66n0BbyO04XyIDxA0hMSKRXSi+61ulKfEy8qz+qx0woPZGwVatWVleii4incvNyiYqIolRUKZ5f8jyPznmUCqUr0DO5Jz2Te9KuRjsaVG5AhPFs+NhaizGG/MJ8+o/vz/wd8zl26hjREdF0qtWJB696kD71+3j5T3VpxpgV1tpWxV0uYI9AjDF1gT8C5ay1qUVZZu2BtbR6oxXxMfEklEqgfOny3NPyHjrV6sS3x75l7va5VIuvRvX46iSWTaRMdBnf/iFEJCiszFzJ6LTRvL/2fd7s9yZDmg7htha30SmpE20S2xAZEemV7RhjAIiKiGLGTTM4mX+SxbsWM3vrbGZsnsH+3P0A7M3ey5hVY0htnErDyg29sm1f8OsRiDFmDNAXOGCtbXrO13sBI4FI4C1r7d/O+d6kohZIxeSKtv0z7cnJyyHrZBZHThzhhR4vMLjxYObvmE/XcV3Pf39sRT4Y/AE9k3uy9fBWPt36KQ0qN6Bh5YZUj69+9n+2iISegsICxq8dz8ilI1m+dzmxUbHc2PRGfnvVb2lerbmTTIW2kAgTwYS1ExgyeQgATas25cYmNzK06VCSKyb7ZLueHoH4u0A6AznA298XiDEmEtgE9AD2AMuBodba9We+X+QC+alTWMdOHWPr4a3sz91PZnYmu4/uZs/RPdzf9n4aVm7IuNXjuP2j28++v2ypsjSp0oSxA8ZSv1J9jp48SkxkDKWjSnu+A0TEuYLCAiIjIim0hTR9rSkWy72t7mVY82FUiK3gOt5ZGUczmJw+mQ/Xf8jiXYuJMBFkPpRJ1biqZ0+FeUtQFAiAMaY2MP2cAmkHPGWtvebM698DWGufO/PaKwVyKYW2kMzsTDYc2kD6oXTSD6az9uBaptwwhUplKvHMgmd4euHTtKjWgjaJbWhXox2danWiVrlaOlIRCQJHThxh5FcjGff1OL7+5dcklEog42gG1ROqezym4S+7snaxaOcibm5+MwCDJw4mJjKGu6+4m651upY4fzCPgSQCu895vQdoa4ypBDwLXGGM+f33hfJDxph7gHsAkpKSPA4RYSJILJtIYtlEutXt9qPvd6/bneP5x1mWsYz31rzHqLRRxEbFcuSxI8RExrA8YzlV46pSq3wtjzOIiPcdPXmUF798kRe/epGsk1n0b9Cf7LxsEkolkFg20XW8Ikkql3S2PKy1JJVNYuzXYxm/djwNKjXggbYPcGuLW8+bBeYPgXAEcj1wjbX27jOvhwFtrLX3FXfd/pqFVVBYwNoDa9n63VYGNRp0ettvtGJF5gqSKyTTK6UXvev1pkvtLhqoF3FoX84+mo1qxqFjhxjQcABPdH6CK6pf4TqWV5zIP8Gk9ZMYuXQkaXvT+OvVf+X3nX7v0bpC9hRWcbicxrvuwDrmbJvDZ9s+Y96OeRw7dYxBjQYx+YbJAHx3/LuAOr8qEqqstaw5sObsQPif5v6JAQ0H0OryYv/7GBSstXy550vqV6pP5TKVmb5pOot3Leahdg9RJa5oDxkM5gKJ4vQgejcgg9OD6DdZa9cVd92Bch3IifwTLNy5kPiYeNrXbM+eo3uo9VIt2tdszw2NbyC1cSrVE6q7jikSclbsXcEDsx5g+d7lbBixgToV6riO5HdPznuSZxY+Q5noMjzS4REebv/wJc+EeFogfh05MsZ8AHwJNDDG7DHG3GWtzQd+DcwG0oGJnpRHICkdVZqeyT1pX7M9ANER0TzR+QmyTmRx/6z7SfxnIl3HdWXN/jWOk4qEhkPHDnH3x3fT+s3WbD68mVd7v0pSOc/HRIPZn7v+mXW/WkevlF48Of9JGrzSgKnpU32yrZC4Et0Y0w/ol5KSMnzz5s2u4/yk9IPpTFg3gUnrJzHn1jlcFn8Z83fMJzcvl2tSriEqIhDmNYgEj5y8HJJfTubw8cM80PYBHu/8OOVKl3MdKyAs2rmI387+Lb9q/SvuvOLOi74vaE5h+VKgnMIqrkETBjF1w1QSExK542d3cFfLu6hdvrbrWCIBLTM78+yp4NFpo+mY1JEmVZs4ThV4Cm0hcHqm6etpr7P9u+088fMnzpuxFRSnsOTCxqeOZ8oNU2herTnPLnqWuiPrct/MYk9CEwkLeQV5PDX/KWqPrM3c7XMB+EWrX6g8LiLCRJy9TmTDoQ08/8XzNB3V9Oy+K9G6S7wGKbGYyBgGNhrIzJtnsvM3O3m88+O0SWwDnL6528ivRnL05FHHKUXcW5axjJajW/LnBX8mtXEqzao2cx0pqLzU6yXm3zafqIgour3djREzRpCTl+Px+nQKK8BNTZ/KoImDSIhJYHjL4Txw1QNhOzgo4e2vi/7K4/Mep3p8dV7v+zp96/d1HSloHTt1jD9+/kdeXvYyc4bN4eq6V2sMJBQLBCBtbxr/+PIffLjuQwBubn4zr/d5ndjoWMfJRPznteWvkbY3jReveVGD5F6y9fBWkismh/cgejDNwiqJXVm7+OeX/2TdwXV8esunGGM4kHuAqnFVXUcT8bpCW8jIr0ZyWfxlDG021Os3EJT/CetBdGvtNGvtPeXKhfZPJUnlknip10tny2N/zn7qjqxL6sRU1h5Y6zqeiNfsy9nHte9dy4OfPsgnWz4BUHkEoJAokHDz/V+k2OhYHmr3EJ9u/ZTmo5ozdPJQth7e6jidSMl8svkTmo9qzqKdixjVZxTjBoxzHUkuIiROYX0vVMdALuXw8cO88MULjFw6koLCAnb8ZgeXxV/mOpZIsa3Zv4bmrzenWdVmjE8dT+MqjV1HCgthPQbyvXAtkO9lZmcye+tsbv/Z7QBM3zSd7nW76yFYEvDyCvKIiYwB4P017zOw4UBNEvGjsB4DkdOqJ1Q/Wx5bDm/hug+uo/GrjZmSPoVQ+kFBQsv8HfNJeTmFVZmrALip2U0qjyChAglRKRVT+GzYZ8TFxDF44mB6vNODjYc2uo4lcpa1ln9++U+6v92duJg4lUYQCokCMcb0M8a8kZWV5TpKQOlWtxurfrGKV659hbS9abT7dzty83JdxxLh2Klj3DTlJh769CH6N+zPsruX0bByQ9expJhCokDCZRqvJ6IiohjRZgQbf72Rdwa+Q1xMHNbas6cLRFwYtXwUE9ZO4LluzzHp+kkklEpwHUk8EBIFIpdWLb4afer3AeC/G/5LyzdacvfHd3PkxBG3wSSs5BXkAfCbq37DkjuX8FjHx3R9RxBTgYShXim9eKzDY4xdPZbGrzZm+qbpriNJGJiwdgJNXmtCZnYmkRGRtKvZznUkKSEVSBiKjY7lue7PsWz4MqrEVaHfB/149LNHXceSEGWt5fklzzNk8hAui7+M6Mho15HES/T4uzDWsnpLlt29jCfnP0nbxLau40gIyi/M5/5P7mdU2ihubHIjYweM1XVJIUQFEuZKRZXib93/dvb1i1++SE5eDn/o9AciIyIdJpNQ8PSCpxmVNopH2j/Cc92fO/tgIwkNKhA5y1rLNwe+YezqsczbMY93B73L5QmXu44lQeyhdg/RsHJDbmp2k+so4gMh8eOArgPxDmMMY64bw5jrxrA0YyktXm/B59s+dx1LgszB3IOMmDGC46eOU650OZVHCAuJAtF1IN5jjOGOK+4gbXgaVeOqcu1717I7a7frWBIkMo5m0HlsZ8asHsPX+792HUd8TKew5IIaVWnE0ruXMnf7XGqWqwmcHhCNitBHRi4s42gGXcd1ZV/OPmbfMpuralzlOpL4WEgcgYhvxMfEc12D6wCYtWUWLUe3ZNt32xynkkD0w/LoXKuz60jiByoQKZJSkaXYc3QPbd9qy6Kdi1zHkQBz9ORRAGbfMlsXCIYRFYgUSdc6XVl691IqxVai29vd+M+q/7iOJAEgNy8Xay2NqjRi/Yj1Ko8wowKRIqtXqR5f3f0VXWp34c6P72TBjgWuI4lDR08epfPYzvxx7h8BND4WhlQgUizlS5dnxk0zGDdgnM5zh7ET+SfoP74/3+z/ho5JHV3HEUdUIFJs0ZHR3NriVowxpB9M59apt3Ls1DHXscRP8gvzGTJpCPN3zGfcgHH0rtfbdSRxRAUiJZK2N413v3mX3u/1Jvtktus44gcPf/owH238iJd7vayLBMNcSJy0NMb0A/qlpKS4jhJ2hrUYRnRkNMOmDqP7O92ZdfMsKsRWcB1LfKhjUkdio2K5r+19rqOIY8Za6zqD17Rq1cqmpaW5jhGWPt74Mdd/eD2NKjfis2GfUSWuiutI4mXZJ7P15MAQZYxZYa1tVdzldApLvOK6Btcxfeh0EssmEhcT5zqOeNn6g+upM7IOU9KnuI4iASQkTmFJYOiR3IPudbtjjCH7ZDZ5BXlUKlPJdSwpoeyT2QycMJCoiChaX97adRwJIDoCEa8yxmCtpf/4/vR4pwffHf/OdSQpAWstd318F1sOb2FC6oSz90UTARWI+IAxhkc6PMK6g+vo9V4vzc4KYv9a9i8+XP8hf736r/y89s9dx5EAowIRn+iV0osPr/+QFXtXMGjiIPIK8lxHEg9knchiQMMBPNLhEddRJABpFpb41NjVY7njozsY0XoEr/R+xXUc8UChLdSjaEOcp7OwNIguPnX7z27nVMEputXt5jqKFMOjnz1Kz+SedKvbTeUhF6VPhvjc8CuHU7dCXay1LN612HUcuYQp6VN4/ovnWbhzoesoEuBUIOI3ry1/jc7/6czU9Kmuo8hF7MvZxz3T7uHK6lfyp85/ch1HApwKRPzmzivupHVia26ZegsrM1e6jiM/YK1l+LTh5OTl8M7Ad4iOjHYdSQJcSBSIMaafMeaNrKws11HkJ8RGx/LRkI+oFFuJfh/0I+NohutIco4Zm2cwfdN0/tb9bzSq0sh1HAkCmoUlfvfN/m/oMKYDTas25Ys7v8AY4zqScHq21dT0qQxsNFAD52FGs7AkaDSv1pzxg8cTFxOn8ggQuXm5xMXEMbjxYNdRJIjoxwxxok/9PnSp3QWAvdl73YYJc8sylpH0UpJmyEmxqUDEqXe/eZeUl1NYsXeF6yhhqaCwgHtn3EvpqNI0r9bcdRwJMioQcapXSi8ql6nM4ImDOXz8sOs4YWdU2ihWZq7kxWtepGypsq7jSJBRgYhTlctUZtINk9ibvZfh04YTSpM6Al3WiSyemv8U3ep04/rG17uOI0FIBSLOtUlsw7NXP8uU9Cm8tfIt13HCxscbP+bw8cP8vcffNZlBPKJZWBIQHmr/EAt3LaTQFrqOEjaGtRhG2xptqV+pvusoEqRUIBIQIkwEHw/5WD8J+8nB3INUiaui8pAS0SksCRjfl8fEdRN5esHTjtOErvUH11PzxZpMWj/JdRQJcioQCThzt8/lqflP8eXuL11HCUmPzXmMUlGlzl6HI+IpFYgEnOd7PE/NcjW546M7OH7quOs4ISVtbxrTNk3jkfaPULlMZddxJMipQCTglC1Vlrf6vcXGbzfy1PynXMcJKc8tfo5ypcpxX9v7XEeREKACkYDUI7kHw1sO54UvX2DTt5tcxwkJB3IPMH3TdO5rc58uGhSvCIlZWMaYfkC/lJQU11HEi17o+QK96/XWTCEvqRpXlc33bSYhJsF1FAkRup27BIXv7xYrnskvzCcqIiR+XhQf8PR27jqFJQFv0vpJ1HqpFruydrmOErQe+OQB+r7fVxdqilepQCTgtb68Ncfzj/ObWb9xHSUo7cvZx79X/Zvq8dX1oCjxKn2aJODVKl+Lxzs/ztQNU5m5eabrOEHn1WWvkleQxyMdHnEdRUKMCkSCwoPtHqRh5Ybc98l9ujakGPIK8nhz5Zv0qd+HepXquY4jIUYFIkEhJjKGV3u/yo4jO/h8++eu4wSNyesnsz93PyNaj3AdRUKQpmVI0Li6ztVsvX8rtcvXdh0laPRK6cWoPqPomdzTdRQJQToCkaDyfXnsztrtNkiQqBBbgV+2+qUGz8Un9KmSoDN29VjqvlyX9IPprqMEtBe+eIEP1nzgOoaEMBWIBJ0+9fpQJroMj8551HWUgHXkxBGenP8kc7bNcR1FQpgKRIJOlbgq/KHjH5i2aRrzts9zHScgjVs9jmOnjjGijQbPxXdUIBKU7m97P0nlknh0zqOE0u14vMFay+gVo7mqxlW0rN7SdRwJYSoQCUqx0bE83vlx0g+ls/nwZtdxAsryvctJP5TOXVfc5TqKhDhN45WgdVuL2xjQcIAejPQDJ/JP0CmpE9c3vt51FAlxuhuvBD1rLYeOHaJKXBXXUUSCku7GK2FryOQhXPvetRoLAbZ9t40jJ464jiFhQgUiQe/alGtZkbmCaZumuY7i3P2f3E/rN1urTMUvVCAS9G5pfgspFVN4cv6TYf0P5/6c/czaMovURqkYY1zHkTCgApGgFxURxZ86/YnV+1Yza8ss13GceW/NexTYAm772W2uo0iYUIFISBjabCiJCYmMWT3GdRQnrLWMXT2WtoltaVi5oes4EiZCYhqvMaYf0C8lJcV1FHEkJjKGWbfMol7F8HzmxcZvN7LmwBpe6/2a6ygSRjSNV0JOoS0My7vPbvp2E1XjqlK+dHnXUSTIaBqvCLBo5yLq/aseO47scB3F7+pXqq/yEL9SgUhIqVOhDruydvHSVy+5juI3Gw5tYMikIWw5vMV1FAkzKhAJKTXK1mBo06G8tfItvjv+nes4fjFp/SQmrJtAbFSs6ygSZjwqEGPMeGPMO2d+Pe/tUCIl8XD7h8k9lcvoFaNdR/GLyemTaVejHYllE11HkTDj6RHIl9baYdbaYcD/eTOQSEk1r9acLrW7MHrFaAoKC1zH8alt321j9b7VpDZOdR1FwpCn03j7G2MKgdnW2k3eDCTiDU93eZojJ46E/BXZk9dPBmBQo0GOk0g48rRAhgEtgMHGmGRr7d1ezCRSYp1qdXIdwS8SSiWQ2jiV2uVru44iYeiS14EYY24D/sHp013TgRHW2mw/ZCs2XQci5zqQe4CXvnqJX1z5C2qVr+U6jkjA8uV1II8DPYCGwE7gr8XdiIgLJ/JP8H9L/o83VrzhOopP7M3ey6mCU65jSBgrSoEctdaustYesNY+DrTxdSgRb0gql0Tf+n15a9Vb5BXkuY7jdbf99zY6/qej6xgSxopSINWNMfcYYzoZY6oA0b4OJeIt97a6lwO5B5iSPsV1FK/Kycth4c6FdEoKj7EeCUxFKZAngebAX4CNQFNjzExjzHPGmKE+TSdSQj2Te1K3Ql1eT3vddRSvmrd9HnkFeVybcq3rKBLGLjkLy1p73glkY0wNThdKM6A38IFvoomUXISJYETrESzNWMqpglNER4bGAfQnWz4hLjqOjkk6hSXuFHsar7V2D7AHmOn9OCLe92C7B11H8CprLTM3z6R73e6UiirlOo6EsZB4HohIUWz+djMpFVNC4uLCtwe+TalIlYe4pZspSliYtnEa9V+pzxe7v3AdpcSMMXSu1Zm2Ndq6jiJhTgUiYaFrna6UiS7D21+/7TpKib2e9jpL9yx1HUNEBSLhIT4mnsGNBjNh3QRO5J9wHcdj2Sezuf+T+5m0fpLrKCIqEAkft7a4layTWUzbOM11FI/N3T6XU4Wn6F2vt+soIioQCR9da3elRtkavLfmPddRPPbJlk+Ij4mnQ1IH11FENAtLwkdkRCSTb5hMvYr1XEfxiLWWWVtm0b1ud2IiY1zHEVGBSHhpkxi8t3I7eOwgWSez6FG3h+soIoAKRMLQB2s+YFnGMl7s9aLrKMVSNa4qh353iFOFugOvBAaNgUjYST+UzsvLXuZA7gHXUYotMiKS0lGlXccQAVQgEoZSG6dSaAuZmj7VdZQis9bSZWwX/rPqP66jiJylApGw06xqM+pXqs+H6z90HaXINn27iQU7F+j0lQQUFYiEHWMMqY1Smb9jPgdzD7qOUyTzdswDTk9FFgkUKhAJS9c3uZ4OSR2CZhxk3o55JCYkklIxxXUUkbM0C0vC0s8u+xkLbl/gOkaRWGuZv2M+PZN7hsSdhCV06AhEwtrh44fJzct1HeMn5eTl0DO5J/0b9HcdReQ8KhAJW5u/3Uy1F6oxcd1E11F+UkKpBN4Z+A6pjVNdRxE5T8AWiDEmzhgzzhjzpjHmZtd5JPSkVEwhMSGRKRumuI7yk/bn7Mda6zqGyI/4tUCMMWOMMQeMMWt/8PVexpiNxpgtxpjHznx5EDDJWjscuM6fOSU8GGPoW78vn2/7nOOnjruOc0GFtpDGrzXmt7N/6zqKyI/4+whkLNDr3C8YYyKBV4FrgcbAUGNMY6AGsPvM2wr8mFHCSN/6fTmef5z5O+a7jnJBa/av4fDxw1xZ/UrXUUR+xK8FYq1dCBz+wZfbAFustdustXnAeKA/sIfTJQI/kdMYc48xJs0Yk3bwYHDM6ZfA0aV2F8pEl2HG5hmuo1zQol2LAOhcq7PjJCI/FgjTeBP535EGnC6OtsDLwCvGmD7ARZ8AZK19A3gDoFWrVjpRLMVSOqo07w58l2bVmrmOckFLdi8hMSGRpHJJrqOI/EggFMiFJrZba20ucIe/w0j4GdhooOsIF/XF7i/okNRB139IQAqEAtkD1DzndQ1gr6MsEoastUxYN4H4mHj61u/rOs5Z1lr+3uPvVIur5jqKyAUFQoEsB+oZY+oAGcAQ4Ca3kSScGGN4fsnzAVcgxhhuaHKD6xgiF+XvabwfAF8CDYwxe4wxd1lr84FfA7OBdGCitXadP3OJ9KnXhyW7l3D4+A/neLizcOdCVmWuch1D5KL8PQtrqLW2urU22lpbw1r77zNfn2mtrW+tTbbWPlvc9Rpj+hlj3sjKyvJ+aAkLfev3pdAWMnvLbNdRznr404d1/YcEtIC9Er04rLXTrLX3lCtXznUUCVKtE1tTpUwVpm+e7joKAMdOHWPVvlV0qNnBdRSRiwqJAhEpqQgTQe96vdl6eGtA3DZkecZy8gvzaV+zvesoIhcVCIPoIgFhdN/RlIoq5ToGcHr6LkC7mu0cJxG5OB2BiJwRKOUB8MWeL2hUuREVYyu6jiJyUToCETnH7+f8nm1HtjEhdYLTHO8Pep+M7AynGUQuJSSOQDQLS7wlryCPjzZ85PzuvAmlEmhYuaHTDCKXEhIFollY4i3d63bnZMFJluxe4izDrC2zeGLeE5zIP+Esg0hRhESBiHhLp1qdiI6IZs62Oc4yTFw3kdeWv0apyMAZkxG5EBWIyDniY+JpV7Od0wJZmrGUtjXa6gaKEvBUICI/cGvzW+lau6uT60GyT2aTfjCdtolt/b5tkeLSLCyRH7ir5V3Otr0icwUWS+vLWzvLIFJUOgIRuYBTBafY9t02v283MzuThJgEWieqQCTwhUSBaBqveNtNU26i5zs9/b7doc2GcuSxI1QuU9nv2xYprpAoEE3jFW/rlNSJrd9tZft32/2+7QgTEn8tJQzokypyAd3rdgdg7va5ftvm/pz9XDH6CqczwESKQwUicgGNKjeiSpkqLNy10G/bXJaxjNX7VhMbFeu3bYqUhApE5AKMMXSu1ZkFOxb4bZvLMpYRaSK5ovoVftumSEloGq/IRTzS4RGOnTqGtdYvF/Ut37ucplWbUia6jM+3JeINKhCRi2iT2MZv27LWsixjGamNU/22TZGS0ikskZ8wb/s8pqRP8fl2cvJy6Jnckx51e/h8WyLeYgLh8Z0lZYzpB/RLSUkZvnnzZtdxJIT0+6AfGw9tZNN9m1xHEfEZY8wKa22r4i4XEkcgug5EfOXntX7O5sObyczO9Ol2cvNyfbp+EV8IiQIR8ZXOtToDsHCnb6fz9ny3JwMnDPTpNkS8TQUi8hNaVm9JfEy8TwvkVMEpVmaupG75uj7bhogvqEBEfkJURBQdanZg1b5VPtvG2gNrOZF/glaXF/sUtIhTmsYrcgnvDnqXCqUr+Gz9aXvTAHQHXgk6KhCRS/D1nXGX711O+dLlSa6Q7NPtiHibCkSkCB6c/SCJCYk81P4hr697QMMBtKjWQo+wlaCjAhEpglX7VrFg5wKfFEjver29vk4RfwiJQXQ9UEp8rWPNjqzet5rsk9leXe/+nP2sylxFfmG+V9cr4g8hUSC6kFB8rVOtThTaQr7a85VX1/vh+g9p+UZL9uXs8+p6RfwhJApExNeuqnEVESaCxbsWe3W9y/cup1pcNRITEr26XhF/UIGIFEHZUmXpW78vpaNKe3W9yzOW0zqxtQbQJShpEF2kiD4a8pFX15d9MpsNhzZwY5MbvbpeEX/REYhIMVhrvTbgvTJzJRarCwglaKlARIoo+2Q2NV+sycivRnplfS2rt2T2LbPpULODV9Yn4m8qEJEiSiiVQOmo0ize7Z2B9IRSCfRM7km50po9KMFJBSJSDB2TOrJ412K88SC2l5e+zKpM392kUcTXVCAixdApqROHjh1i47cbS7Sew8cP88CsB/h066deSibifyoQkWLomNQRoMTXg3x/B17dwl2CmQpEpBjqV6rPw+0epmnVpiVaz9I9SzEYrrz8Si8lE/G/kLgOxBjTD+iXkpLiOoqEOGMMf+/59xKvZ8nuJTSp2oTypcuXPJSIIyFxBKJ7YYk/FRQWsGLvCo6cOOLR8tZa1hxYo+m7EvRCokBE/Onr/V/T6s1WzNg0w6PljTHseGAHz3V7zsvJRPxLBSJSTC2qtaBC6QrM3T7X43VER0ZTIdZ3j8kV8QcViEgxRUZE0qV2Fz7f/rlH14M8Me8JnlnwjA+SifiXCkTEA93qdGNn1k62H9le7GXHfT2ONQfW+CCViH+pQEQ8cHWdqwH4fNvnxVpuz9E97MrapQF0CQkhMY1XxN8aVm7Ip7d8Svua7Yu13JJdS4D/XZAoEsxUICIeMMbQI7lHsZdbvGsxcdFxtLishQ9SifiXTmGJeCgzO5Mn5z3J5m83F3mZuJg4rmtwHVER+tlNgp8+xSIeyivI4+mFT1OpTCXqVapXpGX+1v1vPk4l4j86AhHxUK3ytUiukMzn24s2kF5QWODjRCL+pQIRKYFudboxf8f8Ij3m9i8L/0K9f9UjryDPD8lEfE8FIlICV9e5mqMnj7Iyc+Ul37to1yLiouOIiYzxQzIR31OBiJRA1zpdiYuOIzM78yffdyD3AAt2LqBnck8/JRPxvZAYRNft3MWVqnFV2ffwPuJj4n/yfe+veZ/8wnxua3Gbn5KJ+F5IHIHodu7iUnxMPNZacvNyL/qesavH0vry1jSp2sSPyUR8KySOQERc6/1+b2KjYply45Qffc9ay1NdnqJUZCkHyUR8RwUi4gXNqzbnH1/+g4yjGSSWTTzve8YYBjQc4CaYiA+FxCksEdfuufIeCmwB/1717/O+nleQx18W/oXdWbsdJRPxHRWIiBckV0zmmuRreHPlm+ddEzJj0wwen/c4aw+sdZhOxDdUICJe8stWv2TP0T3nPep27NdjqR5f3aMbL4oEOhWIiJf0rd+Xtwe8TY/kHszZNod+H/Rj5uaZDGs+TDdPlJCkT7WIl0RFRDGsxTAAjp86TsbRDDomdeTe1vc6TibiGyoQER/o16Af/Rr0cx1DxKd0CktERDyiAhEREY+oQERExCMqEBER8YgKREREPKICERERj6hARETEIyoQERHxiLHWus7gNcaYg8BO1zn8rDJwyHWIAKN9cj7tj/Npf/xYA2ttQnEXCqkr0a21VVxn8DdjTJq1tpXrHIFE++R82h/n0/74MWNMmifL6RSWiIh4RAUiIiIeUYEEvzdcBwhA2ifn0/44n/bHj3m0T0JqEF1ERPxHRyAiIuIRFYiIiHhEBRIkjDG9jDEbjTFbjDGPXeD7DY0xXxpjThpjHnaR0Z+KsD9uNsZ8c+bXF8aYFi5y+lMR9kn/M/tjtTEmzRjT0UVOf7nU/jjnfa2NMQXGmFR/5vO3Inw+uhhjss58PlYbY5645EqttfoV4L+ASGArUBeIAb4GGv/gPVWB1sCzwMOuMwfA/mgPVDjz+2uBpa5zB8A+ied/457NgQ2uc7vcH+e8by4wE0h1ndvx56MLML0469URSHBoA2yx1m6z1uYB44H+577BWnvAWrscOOUioJ8VZX98Ya397szLr4Aafs7ob0XZJzn2zL8UQBwQyjNoLrk/zrgPmAwc8Gc4B4q6P4pFBRIcEoHd57zec+Zr4aq4++Mu4BOfJnKvSPvEGDPQGLMBmAHc6adsLlxyfxhjEoGBwOt+zOVKUf/OtDPGfG2M+cQY0+RSK1WBBAdzga+F8k+Pl1Lk/WGM6crpAnnUp4ncK9I+sdZOtdY2BAYAz/g6lENF2R8vAY9aawt8H8e5ouyPlUAta20L4F/Afy+1UhVIcNgD1DzndQ1gr6MsgaBI+8MY0xx4C+hvrf3WT9lcKdZnxFq7EEg2xlT2dTBHirI/WgHjjTE7gFTgNWPMAL+k879L7g9r7VFrbc6Z388Eoi/1+VCBBIflQD1jTB1jTAwwBPjYcSaXLrk/jDFJwBRgmLV2k4OM/laUfZJijDFnft+S04OpoVqsl9wf1to61tra1trawCTgV9ba//o9qX8U5fNx2Tmfjzac7oef/HyE1N14Q5W1Nt8Y82tgNqdnU4yx1q4zxvzyzPdfN8ZcBqQBZYFCY8xvOD3L4qir3L5SlP0BPAFU4vRPlQD5NoTvwFrEfTIYuNUYcwo4Dtx4zqB6SCni/ggbRdwfqcC9xph8Tn8+hlzq86FbmYiIiEd0CktERDyiAhEREY+oQERExCMqEBER8YgKREREPKICERERj6hARLzMGPNvY0yfi3xvtDGmQ3GWEQlUKhAR7/sZp2+XfSFtOX134OIsIxKQdCW6SAkZY+oDY4ByZ/57mbV2zwXe1wjYZK0tKOoyIoFMRyAiJWCMiQLeBR601jYD6gEbLvL2a4FZxVxGJGDpCESkZAYB6dbaZWder+P0fYQu5BrgjmIuIxKwVCAiJdMcWHHO6yuB+WeOMp7n9DMXdnL6tvLlrbV7z9xm/kfL+CeuiPeoQERK5lugKYAx5kpgKKcfVHQv8JG1dsGZ7/UB5l1iGZGgorvxipTAmQfuzOT0szU2Ap2BJOANYLi1Nv/M+14BJllr519sGWttODzPXkKICkTEB4wx/YHrgMPAc8AcoK1KQkKJCkRERDyiabwiIuIRFYiIiHhEBSIiIh5RgYiIiEdUICIi4hEViIiIeEQFIiIiHlGBiIiIR1QgIiLikf8HKqqh1vXeZ8YAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "P_0 = 1.193\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(dc_vector, purcell_numerical / P_0, 'g--')\n",
    "\n",
    "ax.set_xlim(min(dc_vector), max(dc_vector))\n",
    "\n",
    "ax.set_yscale('log')\n",
    "ax.set_xlabel(r'$d_c/d$')\n",
    "ax.set_ylabel(r'$P_n$');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "790de200-2917-4f40-881b-fe3312944d25",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-08T07:12:49.349556Z",
     "iopub.status.busy": "2022-02-08T07:12:49.348999Z",
     "iopub.status.idle": "2022-02-08T07:12:49.420815Z",
     "shell.execute_reply": "2022-02-08T07:12:49.420438Z",
     "shell.execute_reply.started": "2022-02-08T07:12:49.349489Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEOCAYAAACaQSCZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAe4ElEQVR4nO3deXgW5dn38e+ZhARCEkAS2RcFBKkISBS3VqgbVoXaal+QR61LKVqqrdW6tNbW1tqnj77a4gLWoj5VxLeVKlpArAWpUpWgIIKgYY9sYWkghJCF8/0jS7NB5r7hnjsmv89x5DAz18zk5PImP66Za2bM3REREWlMQrwLEBGRLwYFhoiIBKLAEBGRQBQYIiISiAJDREQCUWCIiEggMQ8MM5tmZtvN7ONDtJuZ/d7Mcs3sIzM7JdY1iYhI5MIYYTwDjDpM+0VAv8qvCcATIdQkIiIRinlguPtCYNdhNhkD/K9XeBdob2ZdYl2XiIhEpilcw+gGbKqxnFe5TkREmpCkeBcAWAPrGnxeiZlNoOK0FW3bth02YMCAWNYlItLsLFmyZIe7Z0Wzb1MIjDygR43l7sDmhjZ09yeBJwGys7M9Jycn9tWJiDQjZrYh2n2bwimpWcDVlbOlTgcK3H1LvIsSEZHaYj7CMLMXgBFAppnlAfcCrQDcfQowG/gakAsUAdfGuiYREYlczAPD3cc10u7A92Jdh4iIHJmmcEpKRES+ABQYIiISiAJDREQCUWCIiEggCgwREQlEgSEiIoEoMEREJBAFhoiIBKLAEBGRQBQYIiISiAJDREQCUWCIiEggCgwREQlEgSEiIoEoMEREJBAFhoiIBKLAEBGRQBQYIiISiAJDREQCUWCIiEggCgwREQlEgSEiIoEoMEREJBAFhoiIBKLAEBGRQBQYIiISiAJDREQCUWCIiEggCgwREQlEgSEiIoEoMEREJBAFhoiIBKLAEBGRQBQYIiISiAJDREQCCSUwzGyUma02s1wzu7OB9nZm9qqZLTOzFWZ2bRh1iYhIcDEPDDNLBB4DLgIGAuPMbGCdzb4HrHT3wcAI4CEzS451bSIiElwYI4zTgFx3X+vuJcAMYEydbRxINzMD0oBdQFkItYmISEBhBEY3YFON5bzKdTU9CpwIbAaWA7e4+8G6BzKzCWaWY2Y5+fn5sapXREQaEEZgWAPrvM7yhcBSoCswBHjUzDLq7eT+pLtnu3t2VlbW0a5TREQOI4zAyAN61FjuTsVIoqZrgZleIRdYBwwIoTYREQkojMBYDPQzs+MqL2SPBWbV2WYjcC6AmXUC+gNrQ6hNREQCSor1D3D3MjObBLwOJALT3H2FmU2sbJ8C/BJ4xsyWU3EK6w533xHr2kREJLiYBwaAu88GZtdZN6XG95uBC8KoRUREoqM7vUVEJBAFhoiIBKLAEBGRQBQYIiISiAJDREQCUWCIiEggCgwREQlEgSEiIoEoMEREJBAFhoiIBKLAEBGRQBQYIiISiAJDREQCUWCIiEggCgwREQlEgSEiIoEoMEREJBAFhoiIBKLAEBGRQBQYIiISiAJDREQCUWCIiEggCgwREQlEgSEiIoEoMEREJBAFhoiIBKLAEBGRQBQYIiISiAJDREQCUWCIiEggCgwREQlEgSEiIoEoMEREJBAFhoiIBBJKYJjZKDNbbWa5ZnbnIbYZYWZLzWyFmb0VRl0iIhJcUqx/gJklAo8B5wN5wGIzm+XuK2ts0x54HBjl7hvN7NhY1yUiIpEJY4RxGpDr7mvdvQSYAYyps82VwEx33wjg7ttDqEtERCIQRmB0AzbVWM6rXFfTCUAHM1tgZkvM7OoQ6hIRkQjE/JQUYA2s8wbqGAacC7QB/mVm77r7p7UOZDYBmADQs2fPGJQqIiKHEsYIIw/oUWO5O7C5gW3muvs+d98BLAQG1z2Quz/p7tnunp2VlRWzgkVEpL4wAmMx0M/MjjOzZGAsMKvONq8AXzazJDNLBYYDn4RQm4iIBBTzU1LuXmZmk4DXgURgmruvMLOJle1T3P0TM5sLfAQcBJ5y949jXZuIiARn7nUvJ3wxZGdne05OTrzLEBH5QjGzJe6eHc2+utNbREQCUWCIiEggCgwREQlEgSEiIoEoMEREJBAFhoiIBKLAEBGRQBQYIiISiAJDREQCUWCIiEggCgwREQkk8MMHzew44OeV+5QAP3f3DTGqS0REmphInlb7E+CH7r7LzDoCvwWuj01ZIiLS1ERySmqvu+8CcPedwJ7YlCQiIk1RJCOMvmb2tcrvDegTg3pERKSJiiQw/gJk1VkWEZEWInBguPuzNZfN7I6jX46IiDRVkcyS+n81F4EhwH8f7YJERKRpiuSU1B53v6FqwcyeiEE9IiLSRDU6S8rMkiu/vb9O00+OfjkiItJUBZlW+66ZDXL3dTVXVk2xFRGRliFIYHwXeN7MbqvbYGZzj35JIiLSFDUaGO6+GBgOnGJm/zCzm8xsqpktD7K/iIg0D0GuYfwK+Bg4GdgK3AO0BUa5+wWxLU9ERJqKICOEbwOnuvtJ7n4lFcGRDjxoZhmxLE5ERJqOIIFxQs0L3O6e7+5jgAXAe7EqTEREmpZG78Nw96JDrJ9qZm8e/ZJERKQpOqKL1u6ee7QKERGRpk2znEREJBAFhoiIBKLAEBGRQBQYIiISiAJDREQCUWCIiEggCgwREQkklMAws1FmttrMcs3szsNsd6qZlZvZ5WHUJSIiwcU8MMwsEXgMuAgYCIwzs4GH2O6/gddjXZOIiEQujBHGaUCuu6919xJgBjCmge2+D7wEbA+hJhERiVAYgdEN2FRjOa9yXTUz6wZcBkw53IHMbIKZ5ZhZTn5+/lEvVEREDi2MwLAG1nmd5UeAO9y9/HAHcvcn3T3b3bOzsrKOVn0iIhJAo0+rPQrygB41lrsDm+tskw3MMDOATOBrZlbm7i+HUJ+IiAQQRmAsBvqZ2XHA58BY4MqaG7j7cVXfm9kzwGsKCxGRpiXmgeHuZWY2iYrZT4nANHdfYWYTK9sPe91CRESahjBGGLj7bGB2nXUNBoW7fzuMmkREJDK601vkCJSUl8S7BJHQKDBEovT+5++T8qsU3ljzRrxLEQmFAkMkSgvWLwBgbu7c+BYiEhIFhkiUTu9+OgBd0rvEuRKRcCgwRKJ0YuaJAKQkpsS5EpFwKDBEolRwoIBxJ43j7J5nx7sUkVAoMESi9OzSZ5nx8QwGdx4c71JEQqHAEIlSwYECHGdjwcZ4lyISCgWGSJT2HNgDwL0L7o1zJSLhUGCIRKngQAEAew/sjXMlIuFQYIhEqaC4MjBKFBjSMigwRKL0q6/+itZJrTXCkBZDgSESpTN7nMnF/S7WCENajFCeVivSHM3+bDYX9rmQawZfE+9SREKhEYZIFNydr8/4Omt2r+HS/pfGuxyRUCgwRKJQXFZM6cFS9h7Yy/x183Gv+5p6keZHgSEShaopta999hpf/d+vsr9sf5wrEok9BYZIFKqm1HbP6A7oXgxpGRQYIlGousu7OjA0U0paAAWGSBT6Z/Zn/jXzOf/48wGNMKRl0LRakShkpGQwovcIyg6WARphSMugwBCJwqodq/hgywd8pddXeGXsKwzMGhjvkkRiTqekRKIw57M5jJ85nrTkNEb3H01mama8SxKJOQWGSBSqptWmJKbw2qevsWbXmjhXJBJ7CgyRKBQUF5CenE5JeQmXvnApr6x+Jd4licScAkMkCgUHCmjXuh1pyWmAZklJy6DAEIlCwYECMlIySExIJLVVqmZJSYugWVIiUXj4woerRxXpyekaYUiLoMAQiULPdj2rv09PSWdPyZ44ViMSDgWGSBSe/vBperfvzcjjRvLcZc/RvnX7eJckEnO6hiEShbvevIsZH88AYHj34fTP7B/nikRiT4EhEoWqWVIAb298m5dXvRzfgkRCoMAQiVBJeQnFZcW0S6kIjMcXP87tb9we56pEYk+BIRKhqndhVI0wNEtKWopQAsPMRpnZajPLNbM7G2gfb2YfVX4tMrPBYdQlEo2qx4JUjTDSU9J1H4a0CDGfJWVmicBjwPlAHrDYzGa5+8oam60DznH33WZ2EfAkMDzWtYlEo3f73mz4wYbqmVHpyekUlRZRfrCcxITE+BYnEkNhjDBOA3Ldfa27lwAzgDE1N3D3Re6+u3LxXaB7CHWJRCUpIYme7XqSkZIBVIwwAApLCuNZlkjMhREY3YBNNZbzKtcdyvXAnJhWJHIEPt7+MfcvvJ8dRTsAGD9oPB9N/Ii2yW3jXJlIbIURGNbAOm9wQ7ORVATGHYdon2BmOWaWk5+ffxRLFAkuZ3MOP53/0+r3endK68SgToNIStB9sNK8hREYeUCPGsvdgc11NzKzk4GngDHuvrOhA7n7k+6e7e7ZWVlZMSlWpDHVs6QqL3rn7clj8nuT+XzP5/EsSyTmwgiMxUA/MzvOzJKBscCsmhuYWU9gJnCVu38aQk0iUauaJVV1DWPt7rXcPPdmVu1YFc+yRGIu5mNody8zs0nA60AiMM3dV5jZxMr2KcDPgI7A42YGUObu2bGuTSQaew7sIbVVKq0SWwEVs6QATa2VZi+Uk67uPhuYXWfdlBrf3wDcEEYtIkeqoLig+nQU/GeWlG7ek+ZOd3qLROixix9jxU0rqpc1wpCWQtM6RCKUnJhMcpvk6mWNMKSlUGCIROjBRQ+SlZrFNUOuAaBNUhvW3LyGY9seG+fKRGJLp6REIvTHD//I3z77W/WymXF8h+NJS06LY1UisafAEIlQ3YveAFNypjDzk5lxqkgkHAoMkQjVfHlSld+/93umL58ep4pEwqHAEIlAaXkpRaVF1TftVclIydAsKWn2FBgiEdhbspcES6h3Sio9RS9RkuZPs6REInBMm2Mou6eMg36w1vr05HS2Fm6NU1Ui4VBgiETIzEi02i9K0ghDWgKdkhKJwIrtK7j+levJ3ZVba/3kiybXuvtbpDlSYIhEIHdXLtOWTqt+xHmVjJQMvUBJmj0FhkgEqh5tXnda7T83/JNbX7+VA2UH4lGWSCgUGCIRqPvypCpLty7l4Xcfrn4Ln0hzpMAQiUBVINQdYVQ/gFD3YkgzpsAQiYDjZKZmkpyYXGt91Y18miklzZkCQyQCP/3KT8m/Pb/eer0TQ1oCBYbIUVB1Smpfyb44VyISOwoMkYBKy0sZ+exIXl71cr2207qdRtk9ZVzY98LwCxMJiQJDJKCV+StZsH4B+0v312tLsAQSExIb2Euk+VBgiAS0ZMsSAIZ1HVavrbismImvTWRu7tywyxIJjQJDJKAlm5eQnpxO32P61mtLSkhi6pKpLP58cRwqEwmHAkMkoCVblnBKl1NIsPp/bZISkmid1FqzpKRZ09NqRQLq3b43J3c6+ZDt6cl6Yq00bwoMkYBmXD7jsO0ZKRnsLt4dUjUikVn/7/W8+PGLR3QMBYZIAAf9YIOnomoa2mWoniUlcbfnwB42FmzkpGNPAiqC4oF/PsC0pdPqvcclUgoMkQBunnMz72x6hw8mfICZNbjNc5c9R0pSSsiVSUtXfrCcV1a/wrw181i0aREfb/+Y1FapFN5dCMDdb97NS5+8xHeHfZe7zr6L7vd0j/pnKTBEAsjZnENGSsYhwwKoDosgoxFpfnbv383Ggo1sKNjAhn9vIG9PHicdexJXDb4qpj/3F2/9gl8u/CUZKRmc0f0MLh94ea1rbVcPvprfnPcberbrecQ/S4Eh0oiyg2Us27aMG7NvbHTbB/75AM8vf57lNy4/bLjIF9/Wwq3MXzefcYPGAXDJC5ewaNOi6vZES+Tc48+tDoyXVr7EoE6D6HdMv6g/G9v3beet9W8xf/18Rvcfzai+o7jhlBs4udPJXDbgsgZvHh3Vd1RUP6shCgyRRqzMX0lxWTHDutS/Ya+uTmmdWJG/gg+3fsgpXU4Jobqmy93ZXbybbYXbKCotqr7h8Q9L/sCOoh18f/j3SUtOi3OVkdt7YC8P/eshHlz0IAf9IBf0uYCOqR25++y72V+2n17tetGzXU+y2mZVvz9l+77tXPHnK3CcbundOKPHGeTtyeP2M2/nGyd+g00Fm5j5yUxuOvUmWiW2qvXzDpQd4EfzfsSC9QtYkV/xGuC2rdoyIHMAo/qOome7nkdl9BCEAkOkEUs2V9zhnd01u9FtR/cfTYIlMPOTmc0+MLYWbmXZ1mVsKNhAQXEBt591OwC3zbuNF1e8yLbCbZQeLAWgW3o38m7NA2DWp7N47dPXeOrDp3hmzDN8udeX4/ZnqKn8YDnb921n+77tbNu3jW2F29i2bxvnHX8eQzoP4aNtH3HlS1eyoWADhSWFXDHwCn597q/pmNoRgItPuLjeMTu06QBAVmoWqyet5h/r/sGb695k8ebF9GzXk6SEil/B05dP58437+SJnCeYdNokVmxfQfvW7XngvAdITkxm/vr59MjowfhB4xl53EiGdRlWL1jCYO4e+g89GrKzsz0nJyfeZUgL8PbGt3lh+QtM/trkQNcmvvrsV9lauJWV31sZQnWxUVpeytbCrWwp3MKWvVvYvHcz1w29jpSkFB7+18P8dtFv2Vq4tXr7pIQkin9STGJCIpPfm8ySLUvonNaZTm070TmtM13Tu3JO73OAipHHO5ve4ZqXr2Hd7nVMOm0SD17wYMUvxnXzeX3N6+zav4td+3exu3g3u/bv4v0b3qdVYivmrZmHu3NBnwsOe1rH3SkqLWJH0Q5SklLonNaZ/aX7mbpkKvn78tlRtIP8onzyi/K5dsi1XDf0OtbsWkPfyfXv4n/kwke45fRb2PDvDfzw9R/SJa0LVw2+itO7n37U+tvdmf3ZbG6ddyuf7vyU9OR0vjnwmzw95unq9qN1itPMlrh74//6aWhfBYbI0fXY+48xac4kVt60khOzTox3OY0qKC7g+eXPc93Q62id1JrJ703mlrm34NT+3bD+lvX0at+LGR/PYG7uXIZ2HsqQzkPoc0wfOrXtFPG/eAtLCvnxGz/miZwn2PiDjfRo14P7F97PfQvv45g2x9ChdQc6tOlAxzYdmf7N6aQlp3HBny7gjbVvMLjTYMaeNJZ9JfsYmDWQcYPG4e4Me3IY+UUVgVBcVgzALcNv4ZFRj7C/dD+pv04l0RLpmNqRrNQsstpmcf3Q6/mvk/+L/aX7eXbZs9XrqwKvsckOR1NJeQlrd6+l7zF9q0cfR5sCQyRG5q2ZR3FZMaP7jw68z+a9m3l88eNMzJ5I94zopzDGSml5Kat3rmbL3i38eeWfeX758xSVFjFr7Cwu7X8p7+W9x5zcOXRN70rX9K50SetCl/QudE7rHJPZX6t3rKZnu560adWG8oPlJFjCIX9Bl5SXMH35dP5n0f+wMn8lhnHNkGuq/yV+5UtX0jqpNR3bdCSrbRaZqZkM7jS4+vrJrv27aN+6fYuexabAEImBZ5Y+w3de/Q7ZXbN557p3YvZLpqi0iPx9+ewv28+AzAEAPP/R86zIX0GrhFaM6D2CM3qcQeuk1rX2c3cKDhSws2gnO/fvpF1KO/pn9qekvIR759/LjqId7Ni/g/x9FaderhtyHXecfQfbCrfR+aHOALROas2VJ13JjafeGOgaTVNR9WfPSMlo0b/8o3EkgRHKRW8zGwX8DkgEnnL339Rpt8r2rwFFwLfd/YMwahOpy925/5/3c8/8ezjv+PN46VsvRfVL6a+f/JW/rvorUy6ZQmqr1FptB8oO8PKql/njh3/k72v/juP06dCH3JtzAXh66dO8teEt3J37Ft5H66TW/PjMH/OLkb8AoM/v+7CxYCNlB8uqj3n90Ot5avRTtEpoxe/e+x3tW7cnMzWTzNRMhnQeQq/2vQDomNqRFy9/kczUTIZ2Hlp9YfaLxMxo37p9vMtocWIeGGaWCDwGnA/kAYvNbJa717wieBHQr/JrOPBE5X8PaXfxbqbkTKm1rmt61+pTBy8sf4GCAwW12nu168VF/S4C4E/L/kRRaVGt9r7H9OXc488F4OkPn66e4VHlxMwTq2d0/GHJH2r+GQEYdOwghncfTkl5CX9a9qd6NZ/S5RSGdhlKUWlRrWe6VO1/atdT+dKxX6KguKD6rW41h+andz+dEzqewM6inczJnVPr2IZxVs+z6N2+N9sKt/HmujdrtQF8pddX6JbRjc/3fM7bG9+ud/wRvUdwbNtj2Viwsfox3VXthjGi9wg6tOnA+n+vZ9nWZbXqN4yRx40kLTmNtbvXsmrHqur9qtpH9B5BSlIKa3evZd3uddXrq/57Vs+zSEpIYvu+7RSVFpGRkkGbpDYkWAIJlhCzWSGl5aUUlxWTnpJO2cEyBjw6gDW713DVyVfx1OinSE5Mjuq4a3ev5bmPnmP59uXM/NZMOqd1JiUphQRL4LZ5t/Ho4kfp1a4Xd519F32O6UPX9K7V+84aN4s2SW0oLClk4YaF/H3t39lX+p/Xv14x8AoSLIHM1Ew6tulIx9SO1Y9dNzP23b3vkKd1khKS+NaXvhXVn0latpifkjKzM4Cfu/uFlct3Abj7AzW2mQoscPcXKpdXAyPcfcuhjtu2V1svuq72L/xzep3Dgm8vAOCEySfw2a7ParVfcsIlvDruVQC6PtSVLYW1Dz/upHFM/+Z0ANIfSKewpLBW+4RTJjD10qkVNf+i/l/GW0+/lYcufIi9B/aS8ZuMeu33nnMvPx/xczbv3Uy3/9utXvuD5z/Ij878Eat3rGbAYwPqtU+9ZCoThk1gyeYlZP+h/ohy+jemM27QOBasX8DIZ0fWa686R/3q6lcZPaP+Ofn518xnRO8RTF8+nfEzx9drz/lODsO6DmNqzlQm/m1ivfbVk1ZzQscTeHDRg9z+xu312jffupku6V24d/693Lfwvnrte+7cQ3pKOrfNu42H/vVQvXa/t+Kz+t1Xv8u0pdNISkgiOTGZVgmt6NCmA599v+L/942v3cgrq18B/hNmndM6kzOh4hTm+JnjmbdmHuUHy9lXuo+S8hIGZg1kxU0Vc9zvX3g/x3c4nrEnjT3ii51zPpvDlTOvpPxgOWUHy3jz6jc5o8cZfJL/CXl78jj3+HN1SkVC1aSvYZjZ5cAod7+hcvkqYLi7T6qxzWvAb9z97crlN4E73D2nzrEmABMqF/sDq2NafNOUCeyIdxFNiPqjNvVHfeqT2vq7e3o0O4ZxDaOhf6LVTakg2+DuTwJPHo2ivqjMLCfafx00R+qP2tQf9alPajOzqGcLhTEWzgN61FjuDmyOYhsREYmjMAJjMdDPzI4zs2RgLDCrzjazgKutwulAweGuX4iISPhifkrK3cvMbBLwOhXTaqe5+wozm1jZPgWYTcWU2lwqptVeG+u6vsBa9Cm5Bqg/alN/1Kc+qS3q/vjC3rgnIiLh0nw+EREJRIEhIiKBKDCaIDMbZWarzSzXzO5soH2Amf3LzA6Y2W3xqDFsAfpkvJl9VPm1yMwGx6POsATojzGVfbHUzHLM7Ox41BmWxvqjxnanmll55f1hzVqAz8gIMyuo/IwsNbOfNXpQd9dXE/qiYmLAGuB4IBlYBgyss82xwKnA/cBt8a65ifTJmUCHyu8vAt6Ld91x7o80/nON8mRgVbzrjmd/1NjuH1RMsrk83nXHu0+AEcBrkRxXI4ym5zQg193XunsJMAMYU3MDd9/u7ouB0oYO0AwF6ZNF7r67cvFdKu7laa6C9EehV/5WANrSwI2wzUij/VHp+8BLwPYwi4uToH0SEQVG09MN2FRjOa9yXUsWaZ9cD8w5TPsXXaD+MLPLzGwV8DfgupBqi4dG+8PMugGXAbWfWNp8Bf07c4aZLTOzOWb2pcYOqsBoegI9JqWFCdwnZjaSisC4I6YVxVfQR+n81d0HAF8HfhnrouIoSH88QsXz6cpjX06TEKRPPgB6uftgYDLwcmMHVWA0PXpMSn2B+sTMTgaeAsa4+86QaouHiD4j7r4Q6GNmmbEuLE6C9Ec2MMPM1gOXA4+b2ddDqS4+Gu0Td9/j7oWV388GWjX2GVFgND1BHqXS0jTaJ2bWE5gJXOXun8ahxjAF6Y++lS8mw8xOoeLCZ3MN0Ub7w92Pc/fe7t4b+Atwk7u/HHql4QnyGelc4zNyGhV5cNjPSChv3JPgPMCjVMysM5ADZAAHzewHVMyA2BOvumMpSJ8APwM6UvEvR4Ayb6ZPKA3YH9+k4vlspcB+4P/UuAjerATsjxYlYJ9cDtxoZmVUfEbGNvYZ0aNBREQkEJ2SEhGRQBQYIiISiAJDREQCUWCIiEggCgwREQlEgSEiIoEoMESOgJn90cwuPkTbVDM7K5J9RJoyBYbIkRlCxaOjGzKciifnRrKPSJOlO71FImBmJwDTgHaV/+3s7nkNbHci8Km7lwfdR6Sp0whDJCAzSwKeA25190FAP2DVITa/CJgb4T4iTZpGGCLBfQP4xN3fr1xeQcUzeBpyIXBthPuINGkKDJHgTgaW1FgeBiyoHEX8lor3DWyg4hHr7d19c+Uj1+vtE065IkeXAkMkuJ3ASQBmNgwYR8WLeW4EXnH3tyrbLgbmN7KPyBeOnlYrElDly2VmU/FuidXAV4CewJPAd9y9rHK7R4G/uPuCQ+3j7i3lfezSjCgwRI6QmY0BRgO7gAeAvwPDFQrS3CgwREQkEE2rFRGRQBQYIiISiAJDREQCUWCIiEggCgwREQlEgSEiIoEoMEREJBAFhoiIBKLAEBGRQP4/vqsbkISnNKQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(dc_vector, chi_1D, 'g--')\n",
    "\n",
    "ax.set_xlim(min(dc_vector), max(dc_vector))\n",
    "\n",
    "ax.set_ylim(0, 1)\n",
    "ax.set_xlabel(r'$d_c/d$')\n",
    "ax.set_ylabel(r'$\\chi_\\mathrm{1D}$');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "1420bd3c-5135-4ead-b0e9-013b29c1e686",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-02-08T07:12:51.872105Z",
     "iopub.status.busy": "2022-02-08T07:12:51.871586Z",
     "iopub.status.idle": "2022-02-08T07:12:51.932625Z",
     "shell.execute_reply": "2022-02-08T07:12:51.932236Z",
     "shell.execute_reply.started": "2022-02-08T07:12:51.872040Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The value of d_c / d = 0.451\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAREAAAEGCAYAAABCR6GtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAVmElEQVR4nO3de6ylVXnH8d9PQBhEiyPIZaDONIwBHA0airaahlaNgyGCtlT8ox0jybQJTbUxrVCa3klNbf2jUWunhUgickmRMt5QIFJCLeBILYIDMsogIwMUlYhBgZGnf+w9ztrPnL323medfd59zvl+ksm8a693v+9zLvPMuu31OiIEAPP1vK4DALC0kUQANCGJAGhCEgHQhCQCoMmBXQewEOxDQzq86zCAZW734xFxZH51WSSRXgLZ3HUQwDL3Vw/O9SrdGQBNSCIAmpBEADQhiQBoQhIB0IQkAqAJSQRAE5IIgCYkEQBNSCIAmpBEADQhiQBoQhIB0IQkAqAJSQRAE5IIgCYkEQBNSCIAmpBEADQhiQBoQhIB0IQkAqAJSQRAE5IIgCYkEQBNSCIAmpBEADQhiQBoQhIB0IQkAqBJZ0nE9vG2v2x7u+17bL+3//pq2zfYvr//94u7ihHAaF22RPZIen9EnCTpdZLOt32ypAsk3RQR6yXd1C8DmFGdJZGI2B0Rd/aPn5S0XdIaSWdJuqx/2mWSzu4kQABjmYkxEdtrJb1a0u2SjoqI3VIv0Uh66ZD3bLa9zfY26alFixXAoM6TiO3DJF0j6X0R8aNx3xcRWyLi1Ig4VTp0egECqOo0idg+SL0EcnlEfLr/8qO2j+nXHyPpsa7iAzBal7MzlnSJpO0R8eGiaqukTf3jTZKuW+zYAIzvwA7v/XpJvyPpG7a/3n/tTyV9UNLVts+T9F1J53QTHoBxdJZEIuJWSR5S/cbFjAXTctA83/fsgkaB6ep8YBXA0kYSAdCEJAKgSZcDq5jYJGMMLeMK5X3ydcatm9Qkv4qMmcwSWiIAmpBEADShOzPTRnUPFqoLULtPrnu2Ujff+09qVXG8Z4L71L6WlnNXNloiAJqQRAA0IYkAaMKYyJKyKpXLH9+oPvyTlbqalnGX+Y61jFKLqTYFndXGVvI9yvIk4zDLHy0RAE1IIgCakEQANGFMZOaUffpav1wa7NPn8ZLcb6/V5V0p87WGxZCvU4thVHy18YnsJ8VxbdzlRZX35fdO8v2rWXnjI7READQhiQBoQndmpuWmei6XTfDcdM/N6tXF8Y7KdSaJaVSTvzYFXbvuJHWTWJ3KZfdmVNdnvt2U5d+9oSUCoAlJBEATkgiAJoyJLCn5x3XQkGOp3v/PTyZ9MpV/UBzn8ZJnK3V5HGFVpS5/LWV97Z6S9MLKdWvL3vPXXV4nj++8MJXz96i0/Mc9amiJAGhCEgHQhO5M53I3ZL4/ktpqVmmwuf6qVHdnKo/bXRgVw7OVutp7R3UPal2fMr5R3ba1xXFetZtjWFWpm2Tae/mhJQKgCUkEQBOSCIAmnY6J2L5U0pmSHouIDf3XVku6Sr0O605Jvx0RP+wqxm61/Hhy//9lxXEe51ibyo8Wx9+rxJT7/7Wl49OaDs7Kc1+S6nK5HAfJU+KParg8LrTyxkFKXbdEPiFpY3rtAkk3RcR6STf1ywBmVKdJJCJu0eDKJkk6S9Jl/ePLJJ29mDEBmEzXLZG5HBURuyWp/3eep5Mk2d5se5vtbdJTixoggH2W7DqRiNgiaYsk2cdGx+FMSe5r1/reeZ3DcRPc5zWpXK4byWsiHiuO89Lw2m5ltXUXuT5fN4+BrKncsxzbyN+Tk1K5/PUftbtbbjCPa/k/SW8WWyKP2j5Gkvp/PzbifAAdmsUkslXSpv7xJknXdRgLgBG6nuK9QtLpko6wvUvSX0j6oKSrbZ8n6buSzukuwllXNvNzFyA3m8up2temulHN/mHyNPIkv061c2vTv9maVC6nmfP78nR1OdyWvwe1aeTl1yVp0WkSiYh3Dal646IGAmDeZrE7A2AJIYkAaLJkp3hXptwXry1Bz8u2T6rU5eXq5fLwPK5QlrdX4snyVGcecyhjyF9nXq4+7H3SYHxHVd4n9T5VsdehqS6P99Sm11f2GAktEQBNSCIAmpBEADRhTGTZyOsc8nhAufA3j0dUlsgfnsYcnjihKOTxkrtSuXxvvmf+6H0Zf44n/5qW2zvW1pTkMZt8z/LcvIakZmV/9D+jJQKgCUkEQBO6MzNnkqZy7YHUuXtTelkqPzBYPHBdcZxO3VBMt95d6x5Ig92Jtakufyq2ttR+fSqXU6q521ZuC5GveXMqrxpyLM3/U7srDy0RAE1IIgCakEQANGFMZNkY9fH5svyNVJd2NiuHZQ5Ppx5WHL/Bg3W3pbGLPeU4za50oVem8ouGHKd7StIhxfHjqW7g+5Cnbdemcn7yXyn/06g90HtloyUCoAlJBEATujNLWtnvGPVg8LJpf0Kqy9PB3993+Ej6BO0RxfEp6W1rU3lbMf16b9p57cD0qdkjhhxL9VnvJ/K5tYdg5S7VUZW62gO0aruerTy0RAA0IYkAaEISAdCEMZElpfYAqFyXpyTLXdHzzmZ50KEYK/hxqiqnV3+a6s5O5VOK40fSGMi96dwTi+M8bZvDu75SN7BbWZ7mzh8NKM8dtfNabfyJnc0AYN5IIgCakEQANGFMZEnJffGyn55/lHk8oBwHyR+fz0/PK9eUpHUiXy+OT62EI2n9+/+3OPVrA3VXPfrOgfJzN79g6HX0yVQuNz57/KlUWY575OX9+dxyK4NJdjZb2WMgGS0RAE1IIgCaDO3O2P6MpBhWHxFvm0pEqMjN6Ny9KeWl7OUUb+1h39LgJ2y/P1h1dLmzWXrbbw0Wj9dDPz9+uc8bqLssPj9Q/pd3/t7Pj2+95M3pnuk+5RTvft+Dstv20lSXv87y3El2hkOp1hL5B0n/qN7eeT+R9K/9Pz/W/r8+C872Rtv32d5h+4Jp3w/A/AxtiUTEf0qS7b+JiF8rqj5j+5ZpBmX7AEkflfRm9T4Z9VXbWyPim9O8L4DJjTMmcqTtX9pbsL1O0pHTC0mSdJqkHRHxnYh4RtKVks6a8j0BzMM4U7x/JOlm29/pl9dK2jy1iHrWSEWHutcaeW15gu3N++L4hSmHsxTlPn65e3me4l2TyuW56UFS5XL116W3pU7uM689+OfHJ8RvDtTt8DUD5efHu/cVNqTr/nsql0vk781jF+W0bZ7izWNB853iRak2sHpMROyOiOttr9e+H929EfH0lOPyHK8NDPJGxBZJWyTJPnboADCA6aq1RC61/WL1HtZxvaRbI2Kxnh+4S9LxRfk4SQ8v0r0BTGDomEhEnCHpdPWSyNsl3Wb707Y32/7FKcf1VUnrba+z/XxJ50raOuV7ApiH6phIRPzU9gmS/jIiftgfVD1D0kdsHx0Rp00jqIjYY/sPJH1R0gGSLo2Ie6Zxr+VrknUij6VyuU4knXtKsS7jkfS2QwaLt161b73HV07/1YG6d8ZbB8oPlQ3Pnem6eYzktrKQluUPjG3ksZ6FWifCVgClcQZWj5Z0h+3/kXSppH+OiI/1WwhTExGfl/T5kScC6NTIKd6I+DNJL5d0iaR3S7rf9sXab9gewEo01qd4IyJsP6JeA3aPenNj19i+ISL+ZJoBolRb5p7VpnjXpro85Vs25dP/FeWm6LmbcWMqF7u2P/fJFwxUXXHiewbPLZe25/CeSOVTiuObU91AF2bUFG/5dY+a4mW392FGJhHbfyhpk3qb1v2bpD+OiGdtP0/S/ZJIIsAKNk5L5AhJ74iIB8sXI+I522dOJywAS8XIJBIRf16p276w4QBYatjZbEmp7fY+6txybCNPX+bxk2KKN/+GlE+n25bq8rmPDzmWpM+mXcZOLHaDzzumHZ7K5S7z+WHfA7vT560AcoA7iuO8u1veLb8cB2GKt8SmRACakEQANCGJAGjCmMiSUtvtPffpa9sGvmbEfYr+/5705LpybCM/AS8/LW9g7CWPw6Sn8N1bjNnsSOtW8pYDpbx14o5yGXxez7Ejlcs1JZOsE1msz6EuDbREADQhiQBoQndmSSuneHMTO0/blrt47Ux1p6dy5dfi8XL399xlun/4+/b7pPDqVL593+Gekwarbj1hsLyh2LPqiXSZgV5Hnnp9ZSqXy+JHPfxrZU/j1tASAdCEJAKgCUkEQBPGRJaN3IevfVw9jTkofwSqrE9TsQP3uSvV5XGZct+APEazs3Ld21Nd+lruLqZmD087mw2EsC5d52upnGMq1b7uSbYCWP7TwbREADQhiQBoQndmSZnkx5V3K/vJkGNp/63EdlbO/YGGSw//HpgCzu/L3a/axsh3pnIR0xN5mrnsouRr5u9JGUOOnQd8j4uWCIAmJBEATUgiAJowJjLTRu3ufuCQY0lKO4cNjHuM+sRqOQ6SxwK+P+S8ucpl/Dm+fG7tVzFPk5bL6/N1yunpfM08JrKzUpeneMuvJcez/Kdxa2iJAGhCEgHQhCQCoAljIktKbYyktvO7NLhOIz/oOi97L6+V13fUngQ3ydL7ST5an39N87YCw+Qd3PJWBeX3YdQ157vsffnrpCVi+xzb99h+zvapqe5C2zts32f7LV3EB2B8XbVE7pb0Dkn/Ur5o+2RJ50p6haRjJd1o++UR8bPFDxHAODpJInufnGc7V50l6cqIeFrSA7Z3SDpN0n8vboSLKTfrJ3lod3lubmLnKcuyWzJqF7Ty07e1adzcfckPfKp1AWpdn1HXLX9tczek/J7k5fL5YVbD7i/VH4ie/9mU5ZW3A9qsDayukfRQUd6l/TvwkiTbm21vs71t/zURABbL1Foitm/U/hv6S9JFEXHdsLfN8VrMdWJEbJG0pXevY+c8B8D0TS2JRMSb5vG2XZKOL8rHSXp4YSICMA2zNsW7VdKnbH9YvYHV9ZLu6DakWTLJA73z1GzZx9+Z6vKDr8Zdxj1qqrM2PjBJXe2hXfnccrl6ji/vMP9gcZy/B7UtDzKWvS8622+3vUvSr0j6nO0vSlJE3CPpaknflHS9pPOZmQFmW1ezM9dKunZI3cWSLl7ciADM16zNzgBYYmZtTAQDcn+/9nH63KfP4wHltfK5eWuAcgyiNg5T++h/rs+x53GESa77bKXuyUpd3p2+XEuT7znfcY78vuW/boSWCIAmJBEATejOLCm1JnZtybk0uDw819W6IbX7jHro9aox63L9qF3Qatct35t3Zastbc/n1rpQK3tKN6MlAqAJSQRAE5IIgCaMicyc+U4JTjK1OOppdDULFV+tftSv5bgxjJoq3lkc54/+18Zl8v2X/zRuDS0RAE1IIgCakEQANGFMZKa19LUnecLcqDUSpUl2Pa8tT6+NV4way6iNn5R1+X2TfJ2TrBOpbU2w/NESAdCEJAKgCd2ZJaU2tVhrfo86t3afSZr8415zlEnia1Fep7ajfMay9xItEQBNSCIAmpBEADRhTGTZmGScoOXcScZWFurXq7b7+yR1k9xjEitvWrdESwRAE5IIgCZ0Z9BgWs34Uded70OxFjIG7EVLBEATkgiAJiQRAE0YE8EUMa6wEtASAdCkkyRi+0O277V9l+1rbR9e1F1oe4ft+2y/pYv4AIyvq5bIDZI2RMSrJH1L0oWSZPtkSedKeoWkjZI+ZvuAjmIEMIZOkkhEfCki9n6e+jZJx/WPz5J0ZUQ8HREPSNoh6bQuYgQwnlkYE3mPpC/0j9dIeqio29V/bT+2N9veZnub9NSUQwQwzNRmZ2zfKOnoOaouiojr+udcpN4OL5fvfdsc58dc14+ILZK29K5z7JznAJi+qSWRiHhTrd72JklnSnpjROxNArskHV+cdpykh6cTIYCF0NXszEZJH5D0togo+yJbJZ1r+2Db6yStl3RHFzECGE9Xi80+IulgSTfYlqTbIuL3I+Ie21dL+qZ63ZzzI+JnHcUIYAydJJGIOKFSd7GkixcxHAANZmF2BsASRhIB0IQkAqAJSQRAE5IIgCYkEQBNSCIAmpBEADQhiQBoQhIB0IQkAqAJSQRAE5IIgCYkEQBNSCIAmpBEADQhiQBoQhIB0IQkAqAJSQRAE5IIgCYkEQBNSCIAmpBEADQhiQBoQhIB0IQkAqAJSQRAE5IIgCaOiK5jaGb7/yQ92HUcfUdIerzrIBJiGg8x1b0sIo7MLy6LJDJLbG+LiFO7jqNETOMhpvmhOwOgCUkEQBOSyMLb0nUAcyCm8RDTPDAmAqAJLREATUgiAJqQRBaI7Q/Zvtf2XbavtX14UXeh7R2277P9lkWM6Rzb99h+zvapqa6rmDb277nD9gWLdd854rjU9mO27y5eW237Btv39/9+8SLGc7ztL9ve3v+ZvbfrmMZFElk4N0jaEBGvkvQtSRdKku2TJZ0r6RWSNkr6mO0DFimmuyW9Q9It5YtdxdS/x0clnSHpZEnv6sfShU+o97WXLpB0U0Ssl3RTv7xY9kh6f0ScJOl1ks7vf2+6jGksJJEFEhFfiog9/eJtko7rH58l6cqIeDoiHpC0Q9JpixTT9oi4b46qrmI6TdKOiPhORDwj6cp+LIsuIm6R9IP08lmSLusfXybp7EWMZ3dE3Nk/flLSdklruoxpXCSR6XiPpC/0j9dIeqio29V/rUtdxTSL34vSURGxW+r9o5b00i6CsL1W0qsl3T4rMdUc2HUAS4ntGyUdPUfVRRFxXf+ci9Rrml6+921znL9g8+rjxDTX26YZU0VX910ybB8m6RpJ74uIH9lzfctmC0lkAhHxplq97U2SzpT0xti3AGeXpOOL046T9PBixTTEVGOawfuO61Hbx0TEbtvHSHpsMW9u+yD1EsjlEfHpWYhpHHRnFojtjZI+IOltEfFUUbVV0rm2D7a9TtJ6SXd0EeMMxPRVSettr7P9fPUGd7cuwn3HtVXSpv7xJknDWnILzr0mxyWStkfEh2chprFFBH8W4I96g5MPSfp6/8/Hi7qLJH1b0n2SzljEmN6u3v/+T0t6VNIXZyCmt6o3e/Vt9bpcXf28rpC0W9Kz/e/ReZJeot4MyP39v1cvYjxvUK9rd1fxO/TWLmMa9w/L3gE0oTsDoAlJBEATkgiAJiQRAE1IIgCakETQCdvX216TXjvd9me7ignzQxLBorO9Sr31Dt/rOha0I4lgamz/cn9/lUNsv6C/T8YGSadLurl/zsb+Piy3qrdtAZYYFpthqmz/raRDJK2StCsi/s72P0n6D0lfUW8l5m+ot+L3KkmHRsSZHYWLeaAlgmn7a0lvlnSqpL/vv/Z6SbdKOlHSAxFxf/T+N/tkNyGiBZ/ixbStlnSYpIMkHWL7KEkPRcQz/Y+50xRe4ujOYKpsb1VvB7N1ko5Rb8eun0XEx20fot6H8X49Ir5t+wpJL6Q7s7TQEsHU2P5dSXsi4lP9/VW/ot7H2V8pSRHxU9ubJX3O9uPqdXE2dBYw5oWWCBaN7YMl/VfM+AOqMRmSCIAmzM4AaEISAdCEJAKgCUkEQBOSCIAmJBEATf4fDNzV/6ojRpUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "dc_target = 0.45\n",
    "index = np.argmin(np.abs(dc_target - dc_vector))\n",
    "print('The value of d_c / d = {:.3f}'.format(dc_vector[index]))\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.imshow(real_space_emission[index], 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
}