SimulateBleaching.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 222,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline\n",
    "from _do_calculation import solve_mbes_general, solve_mbes, solve_mbes_without_cavity, get_pulses_from_cavity_field, solve_mbes_wait, solve_mbes_only_decay, solve_obes\n",
    "from _mbes import mbes, mbes_soc, mbes_soc_no_cavity, mbes_soc_det, obes_soc, mbes_masing\n",
    "import matplotlib\n",
    "# matplotlib.use('NbAgg')\n",
    "import numpy as np\n",
    "from scipy.stats import logistic\n",
    "from scipy.optimize import basinhopping\n",
    "from scipy.special import erf\n",
    "from scipy import interpolate\n",
    "from IPython.display import clear_output, Javascript, Markdown\n",
    "import matplotlib.colors as colors\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.io\n",
    "import gc\n",
    "import datetime\n",
    "import time\n",
    "import random\n",
    "import matplotlib\n",
    "from _create_pdf import create_spins, create_spins_inhomo_g\n",
    "from _helper_functions import get_thermal_photons, tsallis, cauchy, find_nearest, odeintz, find_nearest_sorted_left\n",
    "from _steady_state import cav_steady, a_steady, a_steady_maser, a_steady_maser_temp, sz_steady_maser_from_a2, sz_steady_maser_from_a2_temp,  a2_steady_maser\n",
    "from importlib import reload\n",
    "import _spin_params as spin_pars\n",
    "import _sys_helpers\n",
    "_sys_helpers.copy_connection_file()\n",
    "# matplotlib.pyplot.switch_backend('nbagg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 243,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.143362827573153\n"
     ]
    }
   ],
   "source": [
    "gcoll = spin_pars.gcoll\n",
    "q = spin_pars.q\n",
    "quali = spin_pars.quali\n",
    "fr = spin_pars.fr\n",
    "kappa = np.pi*fr/quali  # HWHM\n",
    "delta = spin_pars.delta  # detuning cavity to spin central frequency ws-wc\n",
    "deltac = spin_pars.deltac  # detuning drive to cavity wc-wp\n",
    "gperp = spin_pars.gperp  # HWHM\n",
    "gpar = spin_pars.gpar  # 0.001 #HWHM\n",
    "gammaq = spin_pars.gammaq\n",
    "nbins = spin_pars.nbins  # 701#5000#20\n",
    "nbins = 22001\n",
    "g_number = spin_pars.g_number  # how many bins for inhomogeneous coupling\n",
    "#g0 = gcoll/sqrt(nspins)\n",
    "print(kappa)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 244,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pulse_rwa(t, args):\n",
    "    wgrund = args[0]\n",
    "    fcomps = args[1]\n",
    "    pulse = 0\n",
    "    for i in range(len(fcomps)):\n",
    "        pulse = pulse+(fcomps[i])*np.sin((i+1)*wgrund*t)\n",
    "    return pulse\n",
    "\n",
    "\n",
    "def pulse_rwa_gauss(t, args):\n",
    "    wgrund = args[0]\n",
    "    fcomps = args[1]\n",
    "    nfcomps = len(fcomps)\n",
    "    tend = np.pi/wgrund\n",
    "\n",
    "    width = 2*tend/(nfcomps+5)\n",
    "    t0 = 3*width/2\n",
    "    pulse = 0\n",
    "    for i in range(len(fcomps)):\n",
    "        pulse += fcomps[i]*np.exp(-4*(t-t0)**2*np.log(2)/width**2)\n",
    "        t0 += 0.5*width\n",
    "    pulse *= logistic.cdf(t, loc=0.01, scale=.002) * \\\n",
    "        (1-logistic.cdf(t, loc=tend-0.01, scale=.002))\n",
    "    pulse *= 0.46971947398573455\n",
    "    return pulse\n",
    "\n",
    "\n",
    "def pulse_rwa_gauss_const(t, args):\n",
    "    wgrund = args[0]\n",
    "    fcomps = args[1]\n",
    "    maxampl = args[2]\n",
    "    tend = np.pi/wgrund\n",
    "\n",
    "    spline = 0\n",
    "    for i in range(len(fcomps)):\n",
    "        spline += (fcomps[i])*np.sin((i+1)*wgrund*t)\n",
    "    pulsei = 0\n",
    "    pulseq = 0\n",
    "    pulse = 0\n",
    "    pulse = maxampl*logistic.cdf(t, loc=0.02, scale=0.004) * \\\n",
    "        (1-logistic.cdf(t, loc=tend-0.02, scale=0.004))\n",
    "\n",
    "    pulsei = pulse*np.cos(spline)\n",
    "    pulseq = pulse*np.sin(spline)\n",
    "\n",
    "    return pulsei, pulseq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 245,
   "metadata": {},
   "outputs": [],
   "source": [
    "def g_inhomo_function(x, sigma):\n",
    "    return np.exp(-x**2/sigma**2)\n",
    "\n",
    "\n",
    "x_start = 5\n",
    "x_stop = 20  # µm\n",
    "x = np.linspace(x_start, x_stop, g_number)\n",
    "sigma = 5\n",
    "g_temp = g_inhomo_function(x, sigma)\n",
    "g_list = g_temp/np.sum(g_temp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 246,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEMCAYAAADEXsFmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmUXGd55/HvU71Vt7q71JvUrV2WJWzLNrYsY+MFzBYnLBPCkIAHQgxkHGYmYIYZZyYzZJIMCUkgyQQ8EHASB46TwzIMcWIgNofNaGywkRdky8SWZW2tbkkt9b5UL1XP/FG3WqV2t7u61bduLb/POTpdde+tqufKVv/qve/73tfcHRERkVjUBYiISHFQIIiICKBAEBGRgAJBREQABYKIiAQUCCIiAigQREQkoEAQERFAgSAiIoHqqAtYivb2dt+yZUvUZYiIlIzHHnvstLt35HNsSQXCli1b2Lt3b9RliIiUDDM7ku+xumQkIiKAAkFERAIKBBERARQIIiISUCCIiAgQYiCY2QNm1mdmH51nX9zM/t7M9gQ/42HVISIi+QmzhfB+4I4F9t0K/Iu73wg8GzwXEZEIhRYI7t79ErtvAr4RPL4PeFVYdYiEaTqV5qs/OcbXHutmJpWOuhyR8xLVxLRWYCB4PAi0LXSgmd0G3AawadOm8CsTyZO78+GvPMk39/UC8OihM3zi7S+PuCqR5YuqU7kfWB08TgTP5+Xud7n7bnff3dGR1+xrkYL43r+c4pv7evnIG3bwgVdv46t7u/nRwTNRlyWybFEFwoPAG4PHbwyei5SUz//wBdavruff3bSND79+O+2NtfzVnheiLktk2cIcZfRXZDqVbzWze83sCjPLdjJ/AbjMzPYAlwXPRUrGsf5xHj3Uz7uv3UxNVYx4TRXvuHojP3j2FKdHJ6MuT2RZQutDcPd/O8/mJ4N9E8AtYX22SNge2H8CgDde1jm77U2XreMz3z/IA/tP8K5rNkdVmsiyaWKayDJ8e/9JLu5qZnPbqtltF3c1sbmtge/97FSElYksnwJBZIkmplI8cWyAV+1oP2e7mXHdtnYePdSvIahSkhQIIkv0xLEBplPONVtbX7Tvum1tjEzO8HTPcASViZwfBYLIEv3k0ABmcNXmFwfCtRdkptQ88oKGn0rpUSCILNFPDvdzUWczifqaF+3raKpjY2s9+7qHIqhM5PwoEESWwN3Z1z3IFRtXL3jM5etXs+/4YAGrElkZCgSRJTg+OMFwcoad65oXPObyDQmO9U/QPzZVwMpEzp8CQWQJngk6iy95iUC4bEMCgKeO67KRlBYFgsgS7O8ZJmZwcefCgXDp+kwgPK1AkBKjQBBZgmd6h9navor62qoFj2mO19CViPP8qdECViZy/hQIIkvwTM8wl6xLLHrc9rVNPHdypAAViawcBYJInsYmZzg+OMFFnU2LHrtjTSPPnxollfYCVCayMhQIInk6dHoMgAvaVy1yJOxY28TkTJpj/eNhlyWyYhQIInk62JfpE9i2pnHRY7evzRyjy0ZSShQIInk62DdGzGBzW8Oix25fm7msdEAdy1JCFAgieXqhb5QNLQ3UVS88wiirsa6azuY4L/SNFaAykZWhQBDJ08G+MS7oWLz/IGtLewOHzygQpHQoEETykE47h06Psq1j8f6DrC1tqziiQJASokAQyUPP0ATJ6fSSWgib21ZxenSKkeR0iJWJrBwFgkgesn0BF7Tn30LY2p7pfD5yRkNPpTQoEETycDSYT5DPCKOs7HrL6keQUqFAEMnDsYFxaqqMtc3xvF+TDQ+1EKRUKBBE8tDdP8GGlgaqYpb3axpqq1nTVDc7w1mk2CkQRPJwtH+cDS31S37dlnaNNJLSoUAQycOxgXE2tubff5C1ubVBl4ykZCgQRBYxkpxmcHyajS1LD4T1LfWcGplkciYVQmUiK0uBILKIY/0TAGxsXfolow1BiPQOJle0JpEwKBBEFnFsIHPJZ1kthNWZEOkemFjRmkTCoEAQWUR2TYPl9CFkO6K7B9SPIMVPgSCyiO6BCRrrqmlpqFnya7sScapixvFBtRCk+CkQRBaRHXJqlv8chKzqqhidzXFdMpKSoEAQWcSx/uUNOc1a31LPcQWClIBQA8HMbjWzh83sITPbNWffBWb2QzP7gZl938w2hFmLyHK4O90DE8ualJa1YXW9+hCkJIQWCGbWAnwIuAl4N/DpOYf8e+Bv3P0m4IvAB8OqRWS5hiammZhOzY4WWo4NLfWcGE4ynUqvYGUiKy/MFsI1wB53n3L3Q0CjmdXl7N8PrA4etwKnQqxFZFl6gvkDXYnzCYQG0g4nhjQXQYpbdYjv3QoM5DwfCrb1Bs+/AzxgZu8H6oBXzPcmZnYbcBvApk2bQitWZD69Q5lr/12r87/L6VzrW87ORTifvgiRsIXZQujnbAsAIBFsy/oT4KPufjnwe8DH53sTd7/L3Xe7++6Ojo6wahWZV0/wrX7debUQNBdBSkOYgfAIcIOZ1ZjZJmDU3Sdz9htwOnh8ikzrQaSonBiaoCpmdDTVLX7wAjoTmdZFry4ZSZEL7ZKRuw+Y2WeBBwEHbjezK4A3uPsngT8APm9mM0AN8Bth1SKyXL2DSdY21S1pHYS56qqraG+snb38JFKswuxDwN3vBu6es/nJYN9+4IYwP1/kfPUMTdB1HiOMsroS9WohSNHTxDSRl9A7lKQrsfwO5azORFyjjKToKRBEFuDu9A4lWbciLYQ4PbqfkRQ5BYLIAs6MTTE1k16xFsJwcoaxyZkVqEwkHAoEkQX0rsCktKzssNUTw7psJMVLgSCygOyooHXnMSkta3boqVZOkyKmQBBZQHZU0Eq0ELpm5yKoH0GKlwJBZAE9QxPUVBltq2rP+73WNmcCQSONpJgpEEQW0DuYpDMRJ3Yek9Ky4jVVtK2qnb0VhkgxUiCILKB3aGJFLhdlZeYi6JKRFC8FgsgCegaTrFuBIadZXYm4ZitLUVMgiMwjlXZODidX5LYVWV2Jeg07laKmQBCZx+nRSWbSvqIthM5EnMHxaSamUiv2niIrSYEgMo+VHHKapaGnUuwUCCLz6A3uO9S5on0IwWxl9SNIkVIgiMzjbAthZTuVAQ09laKlQBCZx8nhJLVVMVpXYFJaVra1oaGnUqwUCCLzODGcZG2iDrPzn5SWFa+poqWhRi0EKVoKBJF5nBhK0tm8cpeLsroS9ZxUIEiRUiCIzOPEcJLOFRxhlNWpyWlSxBQIInO4e9BCqFvx9+5MxDU5TYqWAkFkjqGJaSZn0rN3KF1JXc1x+semSE5rcpoUHwWCyBxhTErLyo40OqlWghQhBYLIHNlLOp2Jlb9klA0Z9SNIMVIgiMyRHQUUxiWjs3MRFAhSfBQIInNkv72vaQoxEHTJSIqQAkFkjpPDSdob66itXvl/Ho111TTFq9VCkKKkQBCZIzMHYeX7D7IyC+Xo9hVSfBQIInOENUs5qzNRrxaCFCUFgsgcmRZCeIHQ1azZylKcFAgiOZLTKQbHp0NtIaxNxOkbnWQ6lQ7tM0SWQ4EgkiM7YSyMIadZXYk47nBqZDK0zxBZjlADwcxuNbOHzewhM9s1z/7/YmbfNbMfmNlrw6xFJB9hzlLO0roIUqyqw3pjM2sBPgRcC6wH7gFuyNn/C0DC3V8XVg0iS3UyxFnKWWfXVlY/ghSXMFsI1wB73H3K3Q8BjWaW+6/sV4B40EK4x8wSIdYikpcTIc5Szupq1trKUpzCDIRWYCDn+VCwLWsdkA5aCI8Avz3fm5jZbWa218z29vX1hVasCGS+ta+qraIpXhPaZzTXV1NfU6VAkKITZiD0A6tznieCbbn77w8e3w9cPt+buPtd7r7b3Xd3dHSEUqhI1smQh5wCmFlmcppuXyFFJsxAeAS4wcxqzGwTMOruucMqfgDsDh7vBp4PsRaRvIQ9ByGrMxFXC0GKTmiB4O4DwGeBB4EvAR82syvM7I7gkC8Al5jZ94H3AR8PqxaRfJ0YSobaf5ClQJBiFNooIwB3vxu4e87mJ4N9k8B7wvx8kaVIpZ1TI5Ozo4DC1JWIc3I4SSrtVMUs9M8TyUdeLQQzuz2fbSKl7MzoJKm0hzpLOaszUc9M2jkzqslpUjzyvWT0a/Nsu3UF6xCJ3IkCzFLO6mrWXAQpPi95ycjMbgH+DbDVzP4pZ1cTcCbMwkQKrRCzlLM6cyanvXxj6B8nkpfF+hAeBnqBduDPcraPAPvCKkokCrP3MQpxlnKWbl8hxeglA8HdjwBHgFcWphyR6JwYSlIdM9pXhR8IrQ211FbFNBdBikq+ncpvM7MDZjZkZsNmNmJmw2EXJ1JI2SGnsQKM+onFjLWJOk6qD0GKSL7DTj8BvMXdfxZmMSJROjGcZG1z+K2DrK7menUqS1HJd5TRSYWBlLtCzVLO6kzEZ0c2iRSDxUYZvS14uNfMvgLcC8wOnHb3r4dYm0jBuDsnhpK8ekfh7pfVlYhz//4k7o6ZJqdJ9Ba7ZPSWnMfjwM/lPHdAgSBlYWRyhvGpVEFmKWd1JuJMzaQZGJ+mdVVtwT5XZCGLjTJ6b6EKEYlSIdZBmOvsQjkTCgQpCnl1KpvZp+fZPATsdfd/XNmSRAqvZzAzH2D96vAnpWV1Js4ulLNzndaHkujl26kcB64ADgR/Liez2M37zewvQqpNpGBmZykXMBC0lKYUm3yHnV4IvNbdZwDM7C+BbwNvAJ4KqTaRgukZnCBmsLapcMNO2xvrqIqZboMtRSPfFsJ6YFXO81XAOndPkTPqSKRU9QxmJqVVV4W5ZtS5qmLGmqY6tRCkaCxlYtqTZvYDwIBXAR83s1XAd0KqTaRgeocmCjrCKCszF0H3M5LikFcguPvfmNm3gFeQCYT/5u49we47Fn6lSGnoGZxg5/rCd+x2JeI8e2Kk4J8rMp+XbB+b2UXBz11AF3AMOAp0BttESp670zuULOgIo6zO4PYV7l7wzxaZa7EWwkeA2zj31tdZDrx2xSsSKbD+sSkmZ9KRXDLqSsQZn0oxMjlDc7ym4J8vkmuxiWm3BT9fU5hyRAqvZ7BwC+PMdXZdhKQCQSKX7+2vG8zso2Z2V/B8u5m9OdzSRAqjZ6jwk9KyNBdBikm+Y+z+FpgCrguedwN/EEpFIgXWG8xS7lodzSgj0MppUhzyDYRt7v4JYBrA3SfIjDYSKXk9Q0lqq2O0RXA/oTVNcczUQpDikG8gTJlZPZmOZMxsG5qQJmWiZ3CCdYl4JLegrq2O0d5Yp9nKUhTynZj2u8D9wEYz+3vgeuDWsIoSKaTeoWQkHcpZXYm4WghSFPINhPcA3wS+BrwA3O7up0OrSqSAegYneOW2tsg+v7M5zpEz45F9vkjWUjqV48C/Aj4NfN7Mbg+tKpECmUmlOTkczaS0rK5EfHakk0iU8r11xffM7EHgauA1wAeAncCnQqxNJHSnRiZJezRzELI2tDQwkpxhaGKaRL3mIkh08l0g57tk7nD6I2APcLW7nwqzMJFC6IlwyGnWhpZMGHUPjJOo10I5Ep18LxntIzMP4VIyi+NcGow6Eilp3QOZQNjY0hBZDRuCzz7Wr8tGEq18Lxn9RwAzawTeS6ZPoRMo3GoiIiE41p/pzM1+S4/CxtazLQSRKOV7yeg3gRuBq4AjwN1kLh2JlLRjA+N0NNURr6mKrIZEfQ2NddWzrRWRqOQ77LQe+HPgsewymvkws1vJ3C3VgQ+6++PzHPP7wLvc/cJ831dkpRzrn2BjhK0DADNjQ0u9WggSuXwvGX1yqW9sZi3Ah4BrySzBeQ9ww5xj1gI7lvreIiule3CcXZtaoi6DDS0NCgSJXJgLyF4D7HH3KXc/BDSa2dw+h98B/ijEGkQWNJNK0zOYjLT/ICvTQpjQQjkSqTADoRUYyHk+FGwDMrfQBhrdfd9LvYmZ3WZme81sb19fXziVSkXqHUqSSnukI4yyNrTUMzo5w+D4dNSlSAULMxD6gdU5zxPBtqzfAz622Ju4+13uvtvdd3d0dKxshVLRZoecthZDIGRqUMeyRCnMQHgEuMHMasxsEzDq7rl3SL0A+IyZ3Q90mdmnQ6xF5EWOBdfsi6GFoKGnUgzyHWW0ZO4+YGafBR4kM8rodjO7AniDu3/S3V+ZPdbMnnf3D4VVi8h8uvvHiVm0s5SzZienKRAkQqEFAoC7301mzkKuJ+c5TkNOpeC6ByboStRTUxVmQzk/ifoamuKaiyDRiv5fgkhEjg2MF8UIo6yNLQ0KBImUAkEq1rH+idlLNcVgQ0v97K00RKKgQJCKlJxOcXIkOduZWww2tzVwtH+cdFpzESQaCgSpSEf7x3GHre2roi5l1pb2VUzOpOkd1nKaEg0FglSkQ6fHANjSVjyBsDWo5XBQm0ihKRCkImV/6W4pshYCnA0rkUJTIEhFOnxmjNZVtUW1ZGVnc5y66phaCBIZBYJUpEOnx9jSVjwjjABiMWNL2yoOn1EgSDQUCFKRDp8eL6rLRVlb2ht0yUgio0CQijMxleLEcHK2E7eYbGlfxbH+CVIaeioRUCBIxclekinGFsLWtlVMpdL0DGrGshSeAkEqTrbTtpjmIGRppJFESYEgFedQMbcQgprUsSxRUCBIxTnUN0Z7Yx2NdaHe7HdZ1jTV0VBbpRaCREKBIBXn+b5RLlxTfK0DADNjW0cjz58ajboUqUAKBKko7s7zJ0fZsbYp6lIWtH1NIwdOKhCk8BQIUlFODCcZmZxh+5rGqEtZ0Pa1TZwYTjI0MR11KVJhFAhSUZ4LvnlvL+IWwo61mbDSZSMpNAWCVJQDJ0cAiruFsCYTVtlaRQpFgSAV5flTo7StqqWtsS7qUha0oaWe+pqq2daMSKEoEKSiPHdyhAuLuHUAmZvcXbimkQOn1EKQwlIgSMVwdw6cKu4RRlkaaSRRUCBIxTg5PMlIcobta4u7hQAaaSTRUCBIxXimdwiAizqbI65kcdmRRs+pY1kKSIEgFWP/8WEALu4q/ktGO9clAHimZzjiSqSSKBCkYuzvGWZLWwNN8eJZNnMha5vraG+s5enjQ1GXIhVEgSAV45ne4dlv3sXOzNi5LsHTaiFIASkQpCIMJ6c52j/OJeuKv/8ga+e6Zg6cHGFyJhV1KVIhFAhSEbLX4kspEC5dn2Am7Tx3QsNPpTAUCFIR9geBsLOUAiG4vPV0j/oRpDAUCFIR9h8foqOpjjVN8ahLydvG1nqa4tXsVyBIgYQaCGZ2q5k9bGYPmdmuOft+y8weCfbdaWYWZi1S2R4/OsCVG1dHXcaSmBmXrkuwr1uBIIURWiCYWQvwIeAm4N3Ap+cc8g/ufo27Xw+sBV4bVi1S2c6MTnL4zDi7NrdEXcqS7dq8mmd6hpmYUseyhC/MFsI1wB53n3L3Q0Cjmc3eYtLdD+QcOwXMhFiLVLAnjg4CsGtT6QXCVZtbmEk7P+0ejLoUqQBhBkIrMJDzfCjYdg4zuwnoAn4435uY2W1mttfM9vb19YVRp5S5x48OUB0zLltfGnMQcl25MRNijx0ZWORIkfMXZiD0A7kXbRPBtllmdjnwR8A73N3nexN3v8vdd7v77o6OjtCKlfL1+NEBLu5qpr62KupSlqxlVS3bOlbxuAJBCiDMQHgEuMHMasxsEzDq7pPZnWZ2IXA38E53Px1iHVLBZlJp9nUPsWtTaXUo59q9uZXHjg6QTs/7nUlkxYQWCO4+AHwWeBD4EvBhM7vCzO4IDvkLMi2IL5rZD8zsTWHVIpVr3/EhxqdS7N7yoquVJeOqzS0Mjk/zwmlNUJNwVYf55u5+N5lWQK4ng31vDvOzRQB+dPAMANdta4u4kuW7emsmzH70Qj8Xrin+O7VK6dLENClrDx88zUWdTUW9hvJitrQ1sH51PQ8d0JVVCZcCQcpWcjrF3sMDXLetPepSzouZccOF7Tx88DQp9SNIiBQIUraeODrI5Ey6pC8XZV2/vZ3h5Az7NB9BQqRAkLL1wwN9VMeMV1xQuh3KWdcHofbQ87psJOFRIEjZ+u7PTvKKra00l8AKaYtpa6xj57pmfvicAkHCo0CQsnT0zDjPnRzldRevjbqUFfO6i9ey90g/Z0YnFz9YZBkUCFKWvvOzkwC8/uI1EVeycm7euZa0nz03kZWmQJCy9J2fnWT7mkY2t62KupQVc0lXMxta6nlgvwJBwqFAkLJzaiTJj184w89f2hl1KSvKzLh5Zyf/78BpRid1c2BZeQoEKTv3/bSXtMMvXrEu6lJW3C9c2slUKs0DT5+IuhQpQwoEKTv/+ORxdq5rLsvbPFy1uYXNbQ187bHuqEuRMqRAkLJysG+Ufd1DvPWK9VGXEgoz4+27NvCjF85wrH886nKkzCgQpKz83Y+PUFNl/OKV5Xe5KOttV23ADP7v42olyMpSIEjZGJ+a4WuPdfPzl3axpikedTmhWb+6nhu3d/ClR48yNZOOuhwpIwoEKRv3PtHDSHKG97xyc9SlhO5912/h5PAk33yqJ+pSpIwoEKQspNLOX+15gZ3rmtm9uSXqckL36h0dbF/TyF/vOcQCq8+KLJkCQcrCfT/t4dDpMT742gsxs6jLCZ2Z8es3bmV/zzA/eLYv6nKkTCgQpOTNpNLc+b0DXNTZxM9dUl6T0V7K23ZtYHNbA5944FmttywrQoEgJe9Ljx7lYN8YH379DmKx8m8dZNVUxfjIG3bws95h7tunvgQ5fwoEKWkDY1P86bef47ptbdy8s3zubJqvt1y+jkvXN/Pxb/2M4eR01OVIiVMgSEn72DefYXRyht99y86K6DuYKxYz/vCtl9E3Mskn73826nKkxCkQpGT981O9fP3x4/yHm7bxss7yu01Fvl6+cTW3XreVe358hD0H1MEsy6dAkJJ05MwYv/0PT3H5hgQffN32qMuJ3B03v4yXrW3iw19+khNDyajLkRKlQJCSM5yc5v1f3AvAnbdcSU2V/jeur63is+/eRXI6xW/cs5cx3R5blkH/kqSkjE7O8P4v/ITDp8f4y3ddVVYL4JyvbR2NfOqdV/J0zzAf+LvHmJxJRV2SlBgFgpSMwfEpbr37UR4/Osin3nklr9zWFnVJRef1l6zlj992GXsOnOa9f/sTRjTySJZAgSAl4flTI7z1Mw/x0+5B7rzlSt50eVfUJRWtX969kT//lZfz6KF+fvlzP+LQ6bGoS5ISoUCQopZOO198+DBvufMhRidTfPm2a3njZQqDxbxt1wbuvvVqeoeSvPnTe/g/e4/pnkeyKAWCFK2fHO7n7Z97mN/9p/28Ymsr933weq7a3Bp1WSXjVTs6+Ofbb+SSdc3c8bV9vOPzP+bp40NRlyVFzErpW8Pu3bt97969UZchIUqnnQef6+MLDx/mwef66Giq47dufhlvv2pDRU48WwnptPPVvcf4xAPP0j82xWte1sFtr9rGtRe06u+0ApjZY+6+O69jFQgSNXfnqeND3P/0Cb6xr5ej/eN0NNXx3uu38N7rtlJfWxV1iWVhaGKae350mLsfOkz/2BSbWhv4pSvXc/POTi7ualI4lCkFghS15HSK50+N8tiRAX5yuJ9HD/VzamSSqpjxygvaeMfVG7l5Zye11bqiGYaJqRTfeqqXrz/RzcMHz+AOa5rquGF7O7s2tfDyDat5WWeT/v7LRNEEgpndCtwGOPBBd388Z18c+BtgE3AUeL+7v+QUSwVCaXB3Bsan6Rmc4MRQkt6hCXqGkhw8NcqBU6McOTNG9m7NXYk4V29p5cbt7bzhkrWsbqiNtvgKc2o4yYPP9fHgc308fPAM/WNTANRWxdjS3sCWtlVs7VjF1rZVdCbirGmKs6a5jtaG2oq6s2wpK4pAMLMW4LvAtcB64B53vyFn/weADnf/mJn9D+CUu3/upd5TgZA/dyeVdmbSmZ8pd1KpzPO0B9tTzkw6Pft8JnV23+R0muRMisnpFMnpNBPTKZLB4+R0iuRMiompFEMT02f/jJ99PDPn/vzVMWNTWwM71jSxY20j29c2ceWm1WxoaYjob0jmcne6Byb4afcgT3UPcbBvjEOnRznaP8506sX/Pdsb62iur6Y5XkNzfQ1N8czjpng19TVVxGuqiNfEqKuuoi74Ga+JEa+porY6Rk0sRlXMqK6yzM9Y9meM6qpzn2f3m6FLW0u0lECoDrGOa4A97j4FHDKzRjOrc/fJYP9NwJ8Ej+8D7gBeMhCW68137mFiKjNrc/Z/az/72N1zHmd3+9nHOf8WsgGae3z21Wdfm/sazznuxe9x9jNy3vec98nv88ipO5V2wl4vpbY6Rn1NFYn6GlY31JCor2H96noS9ZnH7Y11rFsdpzNRz7pEnLbGOqr0jbKomRkbWxvY2NrAmy9fN7t9JpWmdyjJyeEkp0YmORX8PD06yUhyhuHkNH0jkxzsG808n+cLwcrXCgbELAgJMhtiwWOzYB8E2+1FrwHLHB+8JrZI2OTuOucxtvBx52y3ebfP3TDfa1obavnqB165YG0rJcxAaAUGcp4PBdt659k/CMw77dTMbiNz2YlNmzYtq5Dta5qYmknP/k1n/8LNLOfxudtnj5t9jc1z3Nl9s49n/2vmbDvnM2zOa1/8P+Hc4170eTl1M08N1TEjZsE3rCqjys79BlYVi2WOOWfb2WNiMaOuOvNNLp7zrS77jS9eXaXLBRWkuio2GxT5mk6lmZwJWpNBy3JyJvg5nWJyJn1OC3YmnWYmlduqTefsy/ycTqVnv1i5++yXo7Sf/aLknrsP0n72i116nteAk06f3ZbrnC+CnPNkvofBa3zefb6M1+Q+aYqH+av6rDA/pR9YnfM8EWybb//cfbPc/S7gLshcMlpOIf/rHVcs52Uiskw1VTFqqmI01hXmF5msjDCHETwC3GBmNWa2CRjNuVwE8CDwxuDxG4PnIiISkdACwd0HgM+S+UX/JeDDZnaFmd0RHPIF4DIz2wNcFjwXEZGIhNqec/e7gbvnbH4y2DcB3BLm54uISP4080RERAAFgoiIBBQIIiICKBBERCSgQBAREaDE7nZqZn3AkajrWIJ24HTURUSgEs9b51w5Su28N7t7Rz4HllQglBoz25vvTaXKSSWet865cpTzeeuSkYiIAAp0eTD8AAAFFElEQVQEEREJKBDCdVfUBUSkEs9b51w5yva81YcgIiKAWggiIhJQIITMzL5oZt/Jeb7LzB4ys4eDNafLhpn9hZn9OPjzX3O2l+05Q2bt8ODcHjKzXVHXEwYzuzI4vx+a2ffM7ALLuNPM9pjZN8ysNeo6w2JmO8xs2sxuKOfzViCEyMwu49xFggDuBN5NZgnRDwVrT5eLz7j7tcB1wC+a2bZge9mec3AuHyJzbu8GPh1pQeHpBX7e3V8F/Cnw+8DNQIO73wh8FfitCOsL2+9wds2Wsj1vBUK4/gfw8ewTM6sDVrn7oWCt6T3A1VEVt9Lc/UDwMw2kgFS5nzM5a4e7+yGgMTjnsuLuJ9x9JHg6BcyQCcFvBNvuA14VQWmhM7NXACeA7mDTTZTpeSsQQmJmNwHPASdzNreRWT86a8G1pEuZmf0qcNDdD1P+57zQ2uFlycxWAX8IfJIXr4teruf9UeCPc56X7XlrwdPzYGZVwEPz7PomcD3wTs69ZNRPZv3orAXXki5WL3XO7v4xM3s98GvAW4LtJX/Oi1hs7fCyYWY1wFeAP3L3Z8xs7rroAwu+uESZ2ZuAve5+xsyym8v2vBUI58HdU8C1c7ebWRPwr4EvA/XATjP77+7+h2Y2Hqwx3QvcQOZabMlY6JwBzOwa4GPALwQr4uHuyVI/50U8AvxB8MuyixevHV4WzCwG/B1wr7vfG2x+EPgl4F7Kd130K4CbzOw6Mkv9XgT8JZnzLbvz1jyEkJnZFuCv3f31wfPdwKcAC7bPXWK0ZJnZ08HD7I2//pO7P1bO5wxgZu8Dfh1w4HZ33xtxSSvOzN5OZt3z7Lk9BdxOZsDA5cAw8B53PxNJgQVgZl8A/hp4mDI9bwWCiIgA6lQWEZGAAkFERAAFgoiIBBQIIiICKBBERCSgQBAREUCBICIiAQWCVCwz+yUzczO7aM72PzOzZ8zszjnb683sweD2HQSvvSdnf7WZ9ZnZN3K2jc55j1vN7H8vUE9tcHtp3UFAIqFAkEp2C5mZt+/MbjCzC4Dr3f0Sd//gnOPfB3w9uH0HwBhwqZnVB8/fABxfbjHB3WC/C7xjue8hcj4UCFKRzKwReDXwfjLBgJm9jMx9aTab2RPBnT1zvQv4xznb/hl4U/D4FuBLeX7+B8zsyeDPITP7frDr3uBzRApOgSCV6q3Ad9x9HzBmZrvc/Vngi8DvuPuV7j6WPdjMaoELglt65/oy8E4zi5O5t80jc/bX5/zifxL4nwDu/jl3v4LM2hDdwJ8Hxz9Nea0XISVEgSCV6hYyq10R/LwleHwZ8NN5jm/n3HUdAAgCZUvw+m/N87oJd78i+4fMokm5PgV8z93vC94vBUwFd8wVKSgFglQcM2sDXgHcH2z6CvAOy9zwfiewf56XTQDxBd7yn8gsK5nX5aKcOm4FNvPi24HXAcmlvJfIStBoBqlEbwe+lV23wN0PmdkJMms1TLv7+NwXuPuAmVWZWdzd5/6yvhsYcvengpXyFmVmVwH/GbgxWHI0u70N6HP36WWdmch5UCBIJboFuNzMDudsayNzff/peV+R8W0yofGd3I3u3k3m0s9S/CaZpRe/H6zEtdfdfx14DfNfehIJndZDEMmTmV0JfMTdfzXEz/g68NtBB7dIQakPQSRP7v4EmW/0VWG8fzCS6V6FgURFLQQREQHUQhARkYACQUREAAWCiIgEFAgiIgIoEEREJKBAEBERQIEgIiKB/w+AXIVKlCd4KAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cooperativity is: 8.276906\n"
     ]
    }
   ],
   "source": [
    "plt.close('all')\n",
    "matplotlib.interactive(False)\n",
    "spins, gs, pdf = create_spins_inhomo_g(\n",
    "    gcoll, g_list, q, delta, gammaq, nbins, holes=False, inhomo=True, hyperfine=True)\n",
    "plt.figure()\n",
    "plt.plot(spins[:len(pdf)], pdf/max(pdf))\n",
    "plt.xlabel('$\\Delta f$ (MHz)')\n",
    "plt.ylabel('weight')\n",
    "plt.show()\n",
    "# %qtconsole\n",
    "coop = gcoll**2/kappa*np.sum(pdf/(gperp+1j*2*np.pi*spins))\n",
    "print(\"Cooperativity is: %f\" % coop.real)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 247,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEMCAYAAAAxoErWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd8VfX9+PHX+2aThDDC3hvZIwwFBERwoCLugaNaqQvtt7baoa2tWqvV/lpxtNa6FScyBAERUGSaMGRvQggrjJAQyLzv3x/nRgJJSEjuzckN7+fjkcfJ/Zxzz3mfR2vefLaoKsYYY0xJPG4HYIwxpvqyJGGMMaZUliSMMcaUypKEMcaYUlmSMMYYUypLEsYYY0plScIYY0ypLEkYY4wplSUJY4wxpQp1O4DKio+P19atW7sdhjHGBJWkpKSDqtqgrOuCPkm0bt2axMREt8MwxpigIiLJ5bnOmpuMMcaUypKEMcaYUlmSMMYYUypLEsYYY0plScIYY0ypLEkYY4wplSUJY4wxpbIkYUwgJC+BjTPdjsKYSgv6yXTGVDuq8N7VkJ8Nj2yC2MZuR2RMhVlNwhh/O5riJAiA5MXuxmJMJVmSMMbf0jaf/P3QNvfiMMYPLEkY428HfUnCEwZHdrgbizGVZH0Sxvjb0d0QFg2Nu8ORnW5HY0ylWE3CGH87tg9iG0HdVpCe4nY0xlSKJQlj/C1zH8Q0hphGcGy/M9rJmCBlScIYf8vc5wx7jWkEBTmQne52RMZUmCUJY/zt2H4nSRTOjzh2wN14jKkESxLG+FNOJuQec2oRMQ2dsmP73Y3JmEqwJGGMP2X6EkJsE6dfomiZMUHItSGwIvIkcCNQ+F/QCFUtEJEFQASQA6xR1QnuRGhMBRzb5xxjrSZhaga350k8o6rvl1B+varurvJojKmsTF+SiGkMkXEQGnkycRgThNxubnpURL4XkYeKlCnwkYjME5GL3ArMmAoprDXENgIRpzZhzU0miLmZJCYCPYGRwFUicqGv/HpVHQzcAbwmIrGnf1FExotIoogkpqWlVV3ExpTl2AEICYfIOs7nmEaQZaObTPAKaJIQkRARWVrCzxOqekgdJ4DJQF8AVT3oO6YAq4H2p99XVV9X1QRVTWjQoEEgX8GYs3PsAEQ3dGoR4JtQZ0nCBK+A9kmoagEwsKRzIlJHVdNFRIBhwNu+32NVNcNXg+gOJAcyRmP8KusAxBT5h0tMQ9i1xL14jKkkNzuu/ykinQABFqjqTBEJA+aLyAkgDHhSVQ+7GKMxZ+fYAWf4a6GYRnD8EBTkQUiYe3EZU0GuJQlVvbOEsjx8zU7GBKWsNGjS4+TnwmGwWWlQu6k7MRlTCW6PbjKm5vB6nZpETKOTZYW/21wJE6QsSRjjLyeOgBY4HdeFfkoSNgrPBCdLEsb4S+FQ19M7rsFqEiZoWZIwxl8KE0HRmkS0JQkT3CxJGOMvhU1KRfskwiKd5TlsroQJUpYkjPGXkpqb4OQOdcYEIUsSxvjL6UtyFIpuaDUJE7QsSRjjL6cvyVEopqHVJEzQsiRhjL+cviRHIVu/yQQxSxLG+EthTeJ0MQ0hNxNys6o+JmMqyZKEMf5y7Aw1icLzxgQZSxLG+ENBvtPcVHRxv0KWJEwQsyRhjD9kHQD1lryIn826NkHMkoQx/pCxxznWblb8nC3yZ4KYJQlj/CEj1TmW1NwUHQ/iseYmE5QsSRjjDxl7nWNJNQlPCEQ3gMy9VRuTMX5gScIYf8hIhZAIqFWv5PO1m56sbRgTRCxJGOMPmXuhdpPis60LxTWHo5YkTPBxLUmISIiIvCAic0VkgYh08ZX3EZFFIrJYRO50Kz5jzkrGnpKbmgrVbu7UJFSrLiZj/MC1Pa6B8cBmVf31aeUTgXFAKrBURKaq6pEqj86Ys5GRCs0SSj8f1xxyj0F2OkTVrbq4jKkkN5ubrgdaich8EXlZRMJFJAKIVtUdqpoLLAT6uRijMWVTdTquS5ojUSjOV8uwJicTZNxMEs2Avao6HMgG7gLqA+lFrkn3lZ1CRMaLSKKIJKal2d7BxmXHD0NBzpmbm+JaOMeju6smJmP8JKDNTSISAiwq4dQM4DAwy/d5FnAN8DYQV+S6ON91p1DV14HXARISEqyR17jraIpzPFNNojCBZFiSMMEloElCVQuAgSWdE5FaQAKwtfCoqtkiclxEWgJ7gcHAnwMZozGVlp7sHOu2Kv2amEbgCbOahAk6bnZcPw+8JSL34tQWbvOVPwxMAgR41TqtTbV3xJck6pwhSXg8zhBZ65MwQca1JOH74391CeWJwKCqj8iYCkpPhsg4iKpz5uviWlhNwgQdm0xnTGUdST5zLaJQnZYnm6aMCRKWJIyprPTkM/dHFKrbxplPkXci8DEZ4yeWJIypDFVI31W+mkS9ts7xiNUmTPCwJGFMZRw7APnZULd12dfWa+McD28PaEjG+JMlCWMqI70cI5sK/VST2BG4eIzxM0sSxlTGkZ3OsU7Lsq+NqgsRcVaTMEHFkoQxlXFoq7PrXGFT0pmIONcdtpqECR6WJIypjINbnFpEaET5rq/XxmoSJqhYkjCmMg5tgfodyn99vbbOWk8FeYGLyRg/siRhTEV5vXBoG8SfRZKo2wa8+ScXBTSmmrMkYUxFZe6BvONQv135vxPf0TmmbQ5MTMb4mSUJYyrq4BbneDbNTQ06Oce0jf6Px5gAsCRhTEUd2uocz6a5KaoOxDaBtE2BickYP7MkYUxFHdwC4THOH/2z0aATpG0ITEzG+JklCWMq6sB6aNDZmf9wNhp0dvokvN7AxGWMH1mSMKYiVGH/OmjU9ey/26Az5GXZCCcTFCxJGFMRmfvgxOGKJwmwfgkTFCxJGFMRB9Y5xwolicIRTtYvYao/SxLGVMR+X5Jo2OXsv1urHsQ2PXkPY6ox1/a4FpEQ4Dmgly+O+1V1vYgsACKAHGCNqk5wK0ZjSrV/vfOHvla9in2/SU/Yu9q/MRkTAK4lCWA8sFlVf13CuetV1XaMN9VXRTutCzXpCVtmQ24WhEf7Ly5j/MzN5qbrgVYiMl9EXhaRcF+5Ah+JyDwRucjF+IwpWX6OM2O6sklCvdbkZKo9N5NEM2Cvqg4HsoG7fOXXq+pg4A7gNRGJPf2LIjJeRBJFJDEtLa3qIjYGYP9a8OZB094Vv0eTHs7RmpxMNRfQ5iZfv8OiEk7NAA4Ds3yfZwHXAKjqQd8xRURWA+2BlUW/rKqvA68DJCQkaECCN6Y0qSucY7O+Fb9H7WZQqz7sXeWfmIwJkIAmCVUtAAaWdE5EagEJwNbCo4gIEKuqGb4aRHcgOZAxGnPWUldAdAOIa17xe4hY57UJCm42Nz0P3OQbzdQf+A9O0povIt8Dc4EnVfWweyEaU4I9K6Bpn7NfjuN0TXrCgQ2Ql+2fuIwJANdGN6nqEeDqEk5Vog5vTIDlZDozpbuOrfy9mvV1NiDa9yO06F/5+xkTADaZzpizsXc1oE5NorKa+xJDyrLK38uYALEkYczZ2J3oHJv5IUnENoK6rS1JmGrNkoQxZ2PXUqjXDqLj/XO/FgMgZbmzqqwx1ZAlCWPKy+uFXUug1fn+u2eL/nBsP6TbID5TPVmSMKa80jZCdjq0GuS/e7YY4BxTlvvvnsb4kSUJY8or2TcvtKUfaxINu0B4rFNDMaYasiRhTHntWuLsZ123tf/u6QmBVhfAju/8d09j/MjNVWCNqZbmbzrAB0uTiY+J4KERHWhaJ8rpWE5e4vxBr+wkutO1udBZEfZoKsQ1Iysnn5fnb2X9ngxGdmnErQNaIv5+pjHlVGaSEJHzgXHAEKAJcAJYi7P+0vuqejSgERpThb78cQ8PfriSxrUjST+Ry/xNB/ji/kE0zd8NmXv82x9RqO1Q57jjO7K73sBt/1vGypR0WtarxeNT1pJy5Di/u+w8/z/XmHI4Y3OTiHwF/ByYDVyKkyS6AI8DkcBUEbkq0EEaUxUOZGTz28/X0LdVXRb8Zhhf3D+IrJwCfvnxKnTrN85F7Uf4/8ENuzqL/e34lr/P3sSKXem8cksfFvx6GLcMaMl/vt3O0u2H/P9cY8qhrD6J21T1blWdpqp7VDVfVY+p6gpVfVFVhwGLqyBOYwLuhTmbyMkv4IXrexIZFsJ5TWrz+OjzWL7jMHtXzoR6bf3bH1HI44HWQ8jbOp+3Fm3nlgEtubx7E0SEJ0Z3oUW9KJ76cj1qcymMC86YJAqX7S5KROKlSANpSdcYE2z2Hc1m8opUbh3QijbxJ3eKuyGhBX2bR1Nn/zIK2gwPXABthxKWtY/ukWk8dknnn4qjwkOYcFEH1u3JYN7GA4F7vjGlKKu5aaCILBCRySLSW0TW4vRH7BeRS6smRGMC750lOylQ5a5BbU4p93iEx3tmUotsvtceAXv+9trOOk7/13oXcbXCTjk3tnczmtWJ4s1FOwL2fGNKU1Zz08vAX4FJwDzg56raGLgQeDbAsRlTJXLyC5i0fBejujSiZf1axc73yllJPiE8vS6evAJvQGJ4eWU+W7U5gwoSi50LC/FwU78WLNp6iJ0HswLyfGNKU1aSCFXVOar6KbBPVZcCqOrGwIdmTNWYvzGN9ON53Ny/ZYnnZctsMhv0YctRYeqqPX5/fsrh40xdvYcDTYYRtnsxZGcUu+aGfi0I8QiTlu/y+/ONOZOykkTRfzadOO2c9aKZGmHa6lTqR4czuH0Ji/Yd2Qn711Kn9xjOa1KbVxdspcDr3//rv7FwOx6BjkOuc/aX2Dav2DWNakdyUeeGfL4i1e/PN+ZMykoSPUUkQ0QygR6+3ws/d6+C+IwJqMzsPOZuOMAVPZoQGlLCfw4bZwIgnUfzwPB2bE/LYva6fX57/qFjOXz0QwpjezcjvvMQiKoLm2eVeO3Y3s04eCyHJdtsOKypOmWNbgpR1dqqGquqob7fCz+Hnem7xgSDWWv3kZvv5apezUq+YOMMZ32lem25rFsT2sZH88r8rX4bjvrWop3kFngZf2E7CAmF9iNhyxzwFhS79qLODYmNCGXqqlS/PNuY8ihrdFO9M/1U5sEi8rxv5NQCEdkrIhN85X1EZJGILBaROyvzDGPKMm31HlrUi6JPyzrFTx4/DLsWQ+fRAIR4hHuHtWPdngwWbEqr9LMzs/N4d8lOLu3amPYNY5zCTpfB8UOQXHz6UWRYCJd0a8ystfvIziueRIwJhLKam5KARN8xDdgMbPH9nlSZB6vqo6o6zDchLw2Y7Ds1EWcZkGHAQyJStzLPMaY0BzKzWbT1IGN6Nit5baTNs0C90Onyn4oKh6O+7IfaxKTlu8jIzue+Ye1OFna8BMJqwbovSvzOmF5NyczJZ77NmTBVpKzmpjaq2hZnWY4rVTVeVesDV3Dyj3qliEgf4ICqpopIBBCtqjtUNRdYCPTzx3OMOd2Xq/fiVbi6d9OSL1g7GeJaQNPePxWFhXgYf2FbkpKPsGzH4Qo/Oye/gDcW7mBw+3h6NC9SiwmPdhLF+qlQkF/sexe0iyc+JoIp1uRkqkh5lwrvp6ozCz+o6lfAUD/FMA74wPd7fSC9yLl0X9kpRGS8iCSKSGJaWuWr/ebcNHX1Hro0qU37hrHFT2YddEYZdbu22KqvN/ZrQXxMBK/M31rhZ7+/dBcHMnO4v2gtolDXsXD8ICR/X+xUiEe4smcT5m9K4+iJvAo/35jyKm+SOCgij4tIaxFpJSJ/AMocYiEiISKytISfJwrPA2OAz31fOQzEFblFnK/sFKr6uqomqGpCgwYNyvkKxpy042AWq1PSGdOrlFrEui9AC6DHDcVORYaF8PMhbVi45SCrU9JL+PKZHT2Rx8R5WxjSIZ4LShp222EUhMc4NZkSjOnVjNx8r19HWRlTmvImiZuBBsAXvp8GvrIzUtUCVR1Yws9TvktGAImqmuG7Phs4LiItRSQMGAzYvo7G76at2oMIXFVakljzmTOqqVHXEk/fOqAlcVFhPD9741n3Tby2YBvpx/N47NLOJV8QFuV0YG+YBgXFaws9m8fRqn4tpgVgYp8xpytXklDVw6r6sKr2VtU+qvpLVa14g+xJ44D3Tyt7GGcZkG+BV1X1iB+eY8xPVJWpq1Pp37oeTeKiil9wJBlSlkL360q9R2xkGI+M6siirYfOqn9gw94M3li4nWv7NKdbs7jSL+x+PZw4AptnFzslIlzVsymLtx3kQGZ2uZ9tTEWUNQT2dREpcdKciESLyF0icmtFH66qt6vq9NPKElV1kKpeoKpvVvTexpRmTepRtqdlcXXvUuZGrPnUOXYrPUkA3DqgFb1a1OGpLzew72jZf6xz8gt49LMfiYsK4/HRZWwi1G4ExDSGlaf/G8oxpldTvAozftxb5nONqYyyahKvAk+IyAYR+VREXhWRN0VkIc4+ErHAZwGP0hg/mrJyD+EhHi7v1qT4Sa8XVrwLrYdA3VZnvE+IR3jh+h5k5xVw/wdJ5OSXPndBVXly2nrWpB7lr9d0p250+JmDDAmFXjc7E+syi/c9tG8YS5cmtQOylpQxRZU1BHaVqt6AMwz1FZwhqdNwVoPtqar/UtWcKojTGL/IL/Ay/cc9DO/coNiS3ADs+BbSk6HvneW6X/uGsTx/XQ9W7ErnvvdXlJgoVJUX5mxi0vJd3DesHZd0bVy+YHuNczrPV39U4ukxvZqyKiWd5EO2MqwJnPL2SRxT1QWqOklVp6jqpkAHZkwgLN52iLTMHK4ubRmOpLchqh50vqLc97yiR1Oevrob8zYe4NrXFrNi18lutJTDx/nFe0m8Mn8bN/dvwW9GdSp/sPHtoeX5TpNTCZ3jV/Z0Ot2tA9sEUqjbARhTlaasSiU2MpThnRsWP3kszVmrqf94CIs8q/uOG9iKRrUjeezzH7nm1cU0jYskLNRD8qHjhId4eHz0edw1qA0eTwkzu8+k9ziY+oCzTEfrQaecaloniv5t6jFlVSoPXtS+5FnjxlSSJQlzzjiRW8DstfsY3aMJkWEhxS9Y/SF486DvHRW6/8gujbig3XAmr0wlaedh8rzKDQktuLZPcxrHnV3S+UnXa2DO47D8P8WSBDhNTn/4Yi3r92bQtekZRksZU0HlThIi4lHVwGzLZUwV+HrDfrJyC0oe1VSQD8vfgFaDoMFZNAmdJjoilNsGtuK2gWfu9C638FrQ5w5YPBHSU6BOi1NOX96tCX+auo5pq/ZYkjABUdYQ2LdFZIaItAL+XUUxGRMQnyXtpklcJAPaFFvpBTZOh6O74PwHqj6wsvT7OaDwwxvFTtWNDmdYpwZ8sTKV/ABtrWrObWV1XBcA1wFPUMIaSsYEi5TDx1m4JY3rE5xtQItZ8grUawsdL6364MpSp4XTkZ70NuQeL3b6xn4tOZCZw3w/LF9uzOnKShIpqnoC+A3gp/qzMVXv06TdANyQ0Lz4yZTlsPsHGHAfeEroq6gOBtwL2enwY/HhsMM7NaBhbAQf2f7XJgDKmifxpO94RFUTqiQiY/yswKt8mpjCkA4NaF63VvELlrwMkXHQ65aqD668Wl3gLFm+6F/FlhAPDfFwXd/mzN90oFwzv405G2fsuBaR28t5n1Wq+qMf4jHG777bnMbeo9n88YouxU+mbYb102DwLyEipuqDKy8RGPIIfDwO1k0utjrtjf1a8OqCbXyamMKEER1cCtLURGU1N7Up508Jez8aUz1MWr6L+tHhjDivUfGTC19wVl0dWA07rE/XaTQ0OA8W/sNZPqSIVvWjuaBdfT5OTKHA65/9t42BMmoSqvrn0s6JSJiq2q4nplpLOXycuRv2M/7CdoSHnvZvokPbnMX8Bt4PMUGwL4nHA0N+BZPvgU0z4bxTZ4XfOqAVD3y4gnkbDzCySwkJ0ZgKKO9+EgCI4yIReQPYHaCYjPGbd5fsRES4/fwSxl0sfBFCwmHQw1UeV4V1vQbqtoFvnytWm7ikayOaxkXy5vc7XArO1ETlShIiMkBE/gUk4yzwtxAoZccUY6qHrJx8Pvohhcu6NaZpndP2jTi0zVk4L+EuiClhiY7qKiQUhj4G+36E9VNOORUa4uH2C1qzZPshNuzNcClAU9OUNZnuGRHZAvwVWAP0BtJU9R3bDMhUd5+v2E1mdj53DW5T/OQ3f4HQSBj0y6oPrLJ63ODsmjfv6WI7193UrwVRYSG8tchqE8Y/yqpJjAf2A68B76vqIcB6xUy15/Uqby3aSa8WdejTsu6pJ3cnOv8Kv2ACxAZh270nBEb8EQ5vK7YpUZ1a4VzbtxlTVu3h4DFbxd9UXllJojHwDHAVsFVE3gOiRMQWBjTV2ux1+9hxMIu7T69FqMKcJyC6IVzwoDvB+UPHS6HFAKdv4rRZ2D8b1Ia8Aq/VJoxflDWZrkBVv1LV24H2wFScHelSReTDyjxYRJ4XkQW+n70iMsFXvkBElviOEyvzDHNuUlUmzttKm/hoLu9+2u5zm76CXYth2G8hItadAP1BBC7+M2TudSbYFdGuQQyXd2vCO4uTOXrcBiCayin36CZVzVbVz1T1WqADUHyH9rOgqo+q6jBVHQakAZOLnL7ed25CZZ5hzk3fbDjA+r0Z3D+s3anrNOXnwNdPQP0O0Ke880SrsVbnQ7drYdE/4UjyKaceGN6eYzn5vLXYahOmcsrquG5Z0g/O5Ln5RcpqVzQAEekDHFDVVF+RAh+JyDwRuaii9zXnJlVl4vytNK8bVXxJ8MUT4dBWuPRvEFLC1qXBaORfQDzOnhNFdGlam4vPa8Rbi3aSmW21CVNxZfUtvIPzR/tMW14p8DbwbgVjGAd8UOTz9ap6UERaAHNFJEFVM4t+QUTG43Sq07Jlywo+1tRE3205yOqUdP46tjthIUX+DXQkGb57Ac67Ejpc7F6A/hbX3JlgN+9p2L4A2g776dSEi9oz5pVFvLskmQeGt3crQhPkREvYO9dvNxcJARaVcGqGqj7lO78Z6K2qxQZ2i8gnwLOqurK0ZyQkJGhiYqLfYjbBy+tVrnz5e46eyOObR4YSEVpkRdePboVt8+DBH5w/rDVJXja80h9CI+De752jz91v/8DynYf57jfDqRsd7mKQproRkaTyLNx6VjOuz5av43tgCT9P+S4ZASQWJgjfjO7avt9jge44E/iMKdO01XtYtyeDX4/qdGqC2DwbNn4JQx+teQkCnP24R/8DDm52aktFPHppZ7Jy8nl5/laXgjPBLqBJohzGAUUHeofi9HV8D8wFnlTVw65EZoJKTn4BL8zZRNemtbmqZ9OTJ7KPwpf/Bw06B8cifhXV4WLocSN8/w/Yv+6n4k6NY7m2T3PeW5JMyuHiGxYZUxZXk4Sq3q6q04t8zlPVvqo6WFUHqOrHbsZngsf7S3ex+8gJfntZZzxFRzTN/oMzTHTMqxBaw5tbLnnW2Rdj2gTwFvxU/KtRHRGBF+dscjE4E6zcrkkYU2mHjuXw0jdbGNIhniEdiqzmumUurHzPWcCveV/3Aqwq0fXhsuchNcnZjtWnSVwUdw1uw5RVe0hKttV0zNmxJGGC3nOzNpKVk8+friyyqVD2UZj+kNPMNPS37gVX1bpd6+yHPe8p2Lfmp+IHhrenUe0I/jh1re03Yc6KJQkT1JKSD/NJ4m7uHtKG9g19M6hVYcavTzYzhUW6G2RVEoErX4KouvD5PZB3AoCYiFCeuKIL6/Zk8MEyGwtiys+ShAla+QVenpiyjsa1I3nooiJbdq6eBGs+gWG/OzeamU4XXR+ufhXSNsDcJ38qHt29CYPbx/P32ZtIy7TF/0z5WJIwQevNRTtYvzeDJ67oQnSEb17owa1OLaLVYGdP6HNV+4thwL2w7N9O3wwgIvx5TFey8wp4esZ6lwM0wcKShAlKWw8c44U5mxnZpRGXd2/sFObnwGc/c0YxXfO6s6T2ueziJ519JybfA+kpgLP434PDOzB11R5mrd3nangmOFiSMEGnwKv85rPV1AoP4Zmx3RDxDXmd9Ttnx7Yxr0JcszPf5FwQFgU3vOdsTPTpHU4SBe4f3o6uTWvz+JQ1HM7KdTlIU91ZkjBB583vd7ByVzp/vqorDWN9ndIr3oXE/8EFD0Hny90NsDqJb+/0T6QmwezfAxAW4uHFG3py9EQef5y61uUATXVnScIElXV7jvL32ZsY1aXRyZnVuxNhxiPQdrjTxGJO1eUqZxe+H96A1c781M6Na/PwiA58+eNepq5KLeMG5lxmScIEjaycfCZ8uJK60WH87doeTjNT5n74eBzENoHr3rR+iNKMeNLpzJ/+kJNUgXuHtqNvq7r8fvIadhzMcjc+U21ZkjBB449T17HzUBb/uqk39aLDnTkAH98KJ9Lhpg+gVj23Q6y+QkLhhnchtjFMugnSdxEa4uGlm3sTGuLhwQ9XkJNfUPZ9zDnHkoQJCp8n7ebzFbt58KIODGxb31mbaPI9zr+Kr3kdGnd3O8TqL7o+3PIp5OfChzdCdgbN6kTxwvU9Wbcng2dnbnQ7QlMNWZIw1d6a3Uf5/RdrGNCmHg9d5Ns8Z84TsGE6XPKM0+ZuyqdBR7jhHUjbBJ/dBQX5jOzSiLsHt+HtxTuZtnqP2xGaasaShKnWDh7L4RfvJVI/OpxXbu1DaIgHlv0Hlr7iTBYbeL/bIQafdsNh9Iuw9Wunj0KVxy7tTL/WdXn0s9Ws2X3U7QhNNWJJwlRbeQVe7v9gBYeycvnPbQnEx0TA2s/hq8ecRewu+auzVpE5ewk/c5YtWfUBzHmc8BDhtXF9qVcrnPHvJXIgM9vtCE01YUnCVEuqyh+nrmP5jsM8f10PujePg01fweTx0PJ8uOa/NpKpsoY+Bv3Hw5KXYdE/iY+J4L93JHDkeC73vW8d2cZhScJUS68u2Mak5bu4b1g7xvRqBtsXwCd3QOMecMvHEF7L7RCDnwhc+hx0u85ZCDDxTbo2jeOF63uSlHyERz5ZjdeWFT/nhbodgDGn+zxpN3+fvYmrezXlN6M6QcpymHQL1G8H4z6HyNpuh1hzeDxw9WuQk+ls8yoeruh7J6lHTvDsVxtpEBvBH6/ocnLpE3POcS1JiEgC8DKQA2QBN6pqpoj0ASYCAryuqm+7FaOpegu3pPHY5z8yqH19nr+uJ56UpfDB9c6s8gwrAAAYLElEQVT4/tum2FyIQAgNd+ZQfHIbTH8YgPEX3sH+jBzeXLSDxrUj+cXQdi4HadziZnPTb4HHVHUosBwY5yuf6Pt9GPCQiNR1JzxT1RJ3HuYX7yXRvmEMr43rS/iuhfD+NRDbCO780jmawAiLdBYD7DAKpj+MrHiHx0efxxU9mvDsVxv5JDHF7QiNS9xsbloH1PH9XhdYIyIRQLSq7gAQkYVAP2COOyGaqrIqJZ073/qBxrUjeffu/tRO+daZTV23Ddw+1RJEVShMFL4ahacgjxdvuIujJ/J47PMfCQ/xcHVvW133XONmTeJz4CURWYuTCKYC9YH0Itek+8pMDbY29Si3/28Z9aLD+fCegTRM/QY+uhniO8CdMyxBVKWwSLjxfeh0Ocz8NRHf/53Xx/VlYJv6/OqTVUy3yXbnnIDWJEQkBFhUwqkZwGXANaqaJCK/A/4Pp6kprsh1ccDhEu47HhgP0LJlS3+HbarQ+j0Z3Pa/ZcRGhvHhPQNovPVj+PKX0KQX3DbZ2avZVK3QCKdGMW0CLHiWqOOH+N8dz3Dn20n88uNVhHqEy7o3cTtKU0UCmiRUtQAYWNI5EbkcSPN9PAC0V9VsETkuIi2BvcBg4M8l3Pd14HWAhIQEG6MXpJKSj/Czt5YTHRHKhz/vT/PVL8GCZ52tN69/ByJi3A7x3BUSCmNecQYKLHmZWscP8eZtL3PHO6t4cNJKXsz3WtPTOcLNPonfAp+ISDbg5WTH9cPAJJzRTa+q6hGX4jMB9P2Wg4x/L5GGsRG8f1dfmi/6nbNxUK9b4cp/QUiY2yEajwdGPQ3R8TD3SWKOHeCdm9/mnk+28suPV5GRncft57d2O0oTYKIa3P8QT0hI0MTERLfDMGdhzrp9PPjhSto2iOa9cefRYPYDsGU2DPk1XPS4LbVRHa3+CKY+CHVbkXPDRzw4+yhfr9/PIyM78uBF7W0eRRASkSRVTSjrOptxbarUh8t2cd8HK+jStDaf3NCYBh9fCVvnwuh/wIgnLEFUVz1vgjumwfHDRLw9kteGZHNNn2a8+PVm/vLlegpsZnaNZUnCVAmvV3n2qw38/os1DOkQz6SRudR+bxRk7oXbvoB+d7sdoilLqwvg53OhVn1C37uaF9qv5a5BbXhr0U7ufT+J47n5bkdoAsCShAm47LwCJkxayX++3c6tA1ryvy6rifroOohuCOPnQ9uhbodoyqt+OydRtLoAz7QH+KO8wV9Gd+CbDfu58T9LOZBhq8fWNJYkTEAdyMzm1jeWMWPNXp64pBVPy6uEfPVraDfC+WNTr63bIZqzFVUXxk2GCx6CxP9x+6YHePe65mxLO8bVryxiw94MtyM0fmRJwgTMil1HuHLi96zbc5R3rqrD3RvGI6smwYWPws2TbKG+YBYSCqOeguvfhv3rGDzvWr68ykOBKte8utgm3dUgliRMQHy4bBc3/Wcp4aEe5o46zNAFNzj9D+M+g4v+YHtB1BRdx8I98yAyjrYzbuKb/kl0axLDhEkreWbGevILvG5HaCrJkoTxq5z8An43+Ud+/8UaBreOZm6n6TT/5n5o2AXuXehMlDM1S8POcM986DKGmO//yse1nuPBhFr8d+EOxv1vGQeP5bgdoakESxLGb3YczOLa1xYzaXkKf+pXwP9yHyVi5Zsw8AH42UyIa+52iCZQImvDdW/CVS/jSU3k19vu4sMLj7ByVzpXvPQ9y7YfcjtCU0GWJIxfTF6xmyteWsjuQ1nM7r+Kn62/Gzlx2Nkk6NK/2gzqc4EI9LkNxn8Lcc24YPkDLO45i7iwfG7+71L+8fVma34KQpYkTKUcy8nnVx+v4lefrObCRrksa/4SnX583tmX4L4l1rx0LmrQEX7+DQy4j/rr3uariN/xcKejvPTNFm7+71JS00+4HaE5C7Ysh6mwlbuO8KtPVpN86Bj/7rqBkSkvId4CuOw56D3OZk8b2DbfWU02I5XN7e7ixs1DKfBE8Ow1PRjdw1aSdZMty2ECJie/gOdmbeTa1xZTL3cPK1q/yqitTyONuzud031uswRhHO2Gw32Lofc4Om59g2XxT3FxnT088OEKHvxwBYezct2O0JTBahLmrKzZfZRHPl3F1v0Z/L82P3DVwf8i4oGRf4a+dzkrhxpTki1fw7QJ6LEDrGx2K3fuGE54VCzPjO3OJV0bux3dOae8NQk3lwo3QSQnv4BX5m/jlflb6VdrHyubv0fc3hVOn8MV/4Q6LdwO0VR3HUbC/UuQOU/QZ+W7JNadxzPczS/ey2Vs72b86cou1KkV7naU5jRWkzBlWrr9EH/4Yg370g7yctM5DEv/HImIhUv/Bj1utKYlc/aSF8P0X8LBTWyOv5g794wlJ6oRj19xHlf3amZLj1eB8tYkLEmYUh3OyuXZmRv4NCmF22NX8IfQ94k4sR/63A4jnoRo237cVEJ+Liz+F3z7dwo84bwZfgt/OzSYge0b8tSYbrRtYDsTBpIlCVNhqsrnK1J5ZsZ6GmQn85/4T2iTsRwa93D2fWjRz+0QTU1yaBvMeAS2zyc9ph2/OXYz3+Z344Fh7bl3WFsiQm0Jl0CwJGEqZG3qUf4yfT2bdu7i6XozuSJ7BhJWy9kQKOEuW3PJBIYqbJoJs38PR3byY8wgHjx0HSH12/L46PO4qHNDa4LyM0sS5qykZebw4pxNTE7cwfjIeUwInUx4/jGk9zi46AmIaeh2iOZckJ8DS16B717AW5DLx6FX8nTGaPp2bMkTo8+jQ6NYtyOsMap9khCRBOBlIAfIAm5U1UwRWQBE+MrXqOqEM93HkkTl5OZ7eXvxDl76ZgtD8pfxTMwn1MvZDW2Hw6inoXE3t0M056KMvfDNX2D1h5wIq8v/yxvLu7nDuWlgO355cQcbBeUHwZAkPgMmquq3IvIksF9VX/MliXGqurs897EkUTGqyux1+/jbVxupd3gVz9aeTKecHyG+E1zyjDO01ar3xm2pSfD1n2DnQg6FN+UvWdfwXfgQHrq4E7cMaGn9FZUQDPMk1gF1fL/XBdb4flfgIxHJBZ5W1XluBFeTLd52kOdmbSJv9yr+Fj2ZgRGJENoALv4H9LnD2VDGmOqgWV+4Yzps/Yb6c//Ev3JfZkfoLP444zr+t7A/j1zSiTE9m+Hx2D9oAsXNmkQPYDqQCWQAF6pqvojEq+pBEWkBzAUSVDXztO+OB8YDtGzZsm9ycnIVRx+c1qYe5blZG9mzdTW/j/qCEd7FaGQcMuhh6P8LiLAhh6Ya83phzafo/KeR9F2sCu3B01lXc6xRPx67tDPDOjWwzu2zUC2am0QkBFhUwqkZwGXABFVNEpHfAfmq+vfTvv8J8KyqriztGdbcVLYdB7N4cc4mVq1ZzaMRX3AFC5GwKOT8++H8ByGqTtk3Maa6yM+BxDfRhS8iWWkkenrw3ImrkVYX8OilnUhoXc/tCINCtWhuUtUCYGBJ50TkciDN9/EA0F6cfwbEqmqGiMQC3QGrJlTQtrRjvDJvKz+uXs6DodN5KWIREhKK9L8fBv8fRMe7HaIxZy80Agbeh/S5AxLfpO+if/Kp9y8s29+d5/4zlvB2g3h4REf6t7Fk4Q9uNjcNBZ4DsgEvMA4naSwFTgBhwD9U9eMz3cdqEsVtPZDJxHlb2fLjEh4MncplnmUQGon0vRMGPQS1m7odojH+k3sckt5Cv/8nknWA5dKNf+RcjbQazMMjOzKwra0MUJJq0dxUFSxJnLRpXyYT520hde1CJoRN5SJJwhseg6f/eBh4P8Q0cDtEYwLntGSxVjowMWc0R1uO4qGRnTi/bX3rsyjCksQ5JHHnYV7/diu5m+Zwb9hXDJS1eCPr4jn/fuh/D0TVdTtEY6pO3glY9QHeRRPxpO8kmaa8mjea7U2v4O6hnRjZpTEhNhrKkkRN5/UqX2/Yz1vfbqRV6peMD5tJO1LxxjRxkkPCXTZayZzbvAWwfire7/+JZ99qDkpd/pt7KYvqXMnNF3bj2j7NiQw7d+dZWJKoobLzCpi8IpVPv1vJkPRp3Bn2NfU4irdRdzwXTICuYyHUZqMa8xNV2L7AaYbasYDjEsXHeRcyLeIKRgw6n3EDW52TM7gtSdQwaZk5TFq+i0WLv+OqnBlcF7qQCHLxth/pJIc2F9oMaWPKsmcVuvRVdO1k8OYzv6AXk+QymvS+nDsGtaZ9w3NnbShLEjXEyl1HeG/RVvLXTeNWzxwGeDbi9YQjvW5Gzn8AGnRyO0Rjgk/mfkh6i/xlbxB6Io2t2oy380exr/XV3DToPIZ3bljj+y0sSQSxnPwCvly9ly+/T6Rn2lRuCZlPQzlCXu2WhA34OfQaZxv+GOMP+Tmw7gvyFr9K2P7VZBDNZ/lDmB8zmqGDBnN9QgviosLcjjIgLEkEoT3pJ/hwaTKbl3/FmLyZXBqSiAeloN3FhA64x1l0z/ZzMMb/VCFlOd5l/4b10/FoHsu8nfmMi4nqOZYbBnagW7M4t6P0K0sSQSKvwMs3Gw7w1ZKVNE2ewg2eBbTx7CMvvA6hfW9D+t0F9dq6HaYx546sg7DqA3KWvUlExk6OaAyfFwwhKX4MQ84fxFW9mhITEfyLYFqSqOZ2Hszikx92sO+H6VyWN4fhIasIxUt20wFE9r/TGaUUFuV2mMacu7xe2PkducveJHTzDDyaz1LveXzBCCJ6jOG6gR3p3iwuaCfoWZKohrLzCpi9bh8LFi+h454pXBuykIaSTk5kA8L63IKnz+0Q397tMI0xpzt2AF35PrnL3yIicxfHNIovCwaQWPcyep5/CVf2ahZ0w2gtSVQTqsqKXenMXL4e77qpXOr91hmhRAi57S4mst+d0GEkhNTMzjFjahSvF3YtJjfpfWT9FMIKTrDD24gpOpSD7cYyvH9fhnZqQFiIx+1Iy2RJwmUph48zNXEH+5KmMej4N4zwrCRc8jleux2RCePw9L4FYhu7HaYxpqJyjqEbppG17D1i9i7Gi7CkoAuzQi8iqucYrurXka5Na1fb5ihLEi44eiKPmT+msmHZ13Q+MJPRIUuJk+NkR9TH0/06wvvcDE162aQ3Y2qaI8kUrJpETuL71MpK4YSGM9fbh6TaI2je70qu7NOGRrUj3Y7yFJYkqkh2XgHzN+xnReIi6u38kivle5rLQfI8keR2HE10wi3QZphtCWrMuUAVdi0lZ9XH6NopROYd4ajWYpa3P1sbXkLbfpdxSfdm1It2v//CkkQA5eQX8N2mNJJ+WETcjhmM0sW08+zFi4djzYYQ2/9WpPNoW2DPmHNZQR5s/5bMxEmEb/2KiIIs0jSOmd6B7Gp2OV0SRjCqW2NiI93pj7Qk4Wd5BV4WbUnjh+WLid0+nYu9S2jv2YMXDxmNBhDb9zpCuoyxPRuMMcXlnUA3zybzh4+olfwNoZrLbo3na+1HWvNL6DJgJCPOa0JUeNVNlrUk4Qe5+V6WbjtIUtISojZPZ4R3MR08qXjxcLRhf2onXE9Il6sgpmFAnm+MqYGyj6IbvuTois+J2f0doZrHAa3DN9qP/c0voW3CKIZ1aUrtANcwLElUUFZOPt9u2sfmxHnE7ZrDUO8PtPXsw4uH9Ab9iE24nrCuYywxGGMqLzsD76ZZHEn6jNiUBYRrDoc1hrnaj5TGF9Oiz2WM6Nac+jERfn90tU8SItIReN33cQXwiKqqiPQBJgICvK6qb5/pPv5IEoezcpm/NpnUFbNosvcbhksS8ZJBvoSS3uh84nqNIazrVRDbqFLPMcaYUuVm4d3yNelJnxO9cy4R3uNkaC3me3uxo/5Q4nuP5qKe7Wlaxz8rMQRDkpgC/E1Vl4rIa8BkVf1aRBYB44BUYCkwQlWPlHafiiaJPeknmL9qE+mrptPu8HdcKKupJTlke6I51nI4dfuMJaTjKIisXcE3NMaYCsrLRrfPJz3pcyK2f02t/HRyNYSl3i5sqD2YiK6jOb9PLzo2iqnwPIzyJgk3x2V2BAr/ui8HhovId0C0qu4AEJGFQD9gjj8frF4vO14azY0FKwkVL5kRDTjR/gai+o4lsvUQIm1nN2OMm8IikU6XUbfTZc42rCnLyVo1le4bZ3Jh1muw/DXWLm3N4QdmUL9R84CG4maSWANcKiIzgEuBI0B9IL3INem+slOIyHhgPEDLli3P+sHi8dC6Yw8yoy+gbp+xxDbpDZ7qP43eGHMO8oRAq/Op2+p8GPM3OLiFzNVTid+2lPoNmwX88QFtbhKREGBRCadmAG/h9D3EAluBPcDzwBJV7e37/j+Br1R1dmnPcHsynTHGBKNq0dykqgXAwDNcMlacBrV3cfokskXkuIi0BPYCg4E/BzJGY4wxpXOtuUlEbgHuARR4T1XX+k49DEzCGd306pk6rY0xxgSWa0lCVT8EPiyhPBEYVPURGWOMOZ311hpjjCmVJQljjDGlsiRhjDGmVJYkjDHGlMqShDHGmFIF/SqwIpIGJLsdx1mKBw66HUQVs3c+N9g7B49WqlrmBjhBnySCkYgklmemY01i73xusHeueay5yRhjTKksSRhjjCmVJQl3vF72JTWOvfO5wd65hrE+CWOMMaWymoQxxphSWZJwiYi8IyJzi3zuIyKLRGSxiNzpYmh+JyL/FJGlvp/fFimvye98p++9Fvn2ba9xRKS37/2+E5F5ItJWHBNFZKGIfCki9dyOMxBEpKOI5InI4Jr+zpYkXCAi3YE6pxVPxNnbexjwkIjUreq4AugVVR0IXACMEZF2vvIa+c6+93gI573GAS+5GlDg7AUuVdULgRdw9n65BKilqkOAT4BHXYwvkJ4AvvX9XqPf2ZKEO/4I/LXwg4hE4NvbW1VzgcK9vWsEVd3iO3qBAqCghr/zAGChqub69muP8b1vjaKq+1Q10/cxF8jHSYxf+sqmAxe6EFpAiUh/YB+w21c0jBr8zpYkqpiIDAM2A/uLFJdrb+9gJyK3AdtUdSc1+53r4ezZXuior6xGEpFo4Bng75z67unUzPd+HPhbkc81+p1d23SoJitjb+9BwE2c2tx0GIgr8jnOVxY0zvTOqvqUiFwM3AFc6SsP+nc+g8Oc+r9vTXq3U4hIGPAx8KyqrheRou8ex6nJMuiJyGggUVUPOTsvA6f+713j3tmSRACUtre3iMQC1wIfAVFAVxH5g6o+E+x7e59pP3MRGQA8BVymqid819fk/cyXAU/7/oA2AY6pao7LMfmdiHiA94EpqjrFV/wtMBaYAlzOyXb7mqIXMExELgC6A52B13DetUa+s82TcImItAbeUNWLfZ8TgH/h7O39hqq+6V50/iUihfuXFy6C9oiqJtXwd74L+DnOHu4P+7blrVFE5DrgbaDw3dbg7FE/EegBZAC3q+ohVwIMMBF5G3gDWEwNfmdLEsYYY0plHdfGGGNKZUnCGGNMqSxJGGOMKZUlCWOMMaWyJGGMMaZUliSMMcaUypKEMcaYUlmSMKYIERkrIioinU8rf1FE1ovIxNPKo0TkW9+yJPi++16R86EikiYiXxYpO3baPe4UkZdLiSfctxS3rY5gXGFJwphT3Ywzg/imwgIRaQsMUtUuqjrhtOvvAib7liUByAK6iUiU7/NIILWiwfhWyP0GuLGi9zCmMixJGOMjIjHAUOBunGSBiHTCWYunlYis9K14WtStwNTTyr4CRvt+vxmYVM7n3ysiq3w/O0Rkvu/UFN9zjKlyliSMOelqYK6q/ghkiUgfVd0EvAM8oaq9VTWr8GIRCQfa+pY+L+oj4CYRicRZz2fZaeejiiSDVcBfAFT136raC2dfjd3AP3zXr6Xm7LVhgowlCWNOuhlnZzF8x5t9v3cHVpdwfTyn7okBgC/JtPZ9f2YJ3zuhqr0Kf3A2oSrqX8A8VZ3uu18BkOtbRdiYKmVJwhhAROoD/YFZvqKPgRvF2TSgK7CuhK+dACJLueU0nC09y9XUVCSOO4FWFF82PQLIPpt7GeMPNmLCGMd1wMzCfR9UdYeI7MPZ5yJPVY+f/gVVPSIiISISqaqn/wF/Eziqqmt8uxGWSUT6Ar8Ghvi2ei0srw+kqWpehd7MmEqwJGGM42agh4jsLFJWH6e/YG2J33DMwUkkc4sWqupunGajs/EgztaX8327niWq6s+B4ZTcbGVMwNl+EsZUgoj0Bn6lqrcF8BmTgd/5OtGNqVLWJ2FMJajqSpx/+YcE4v6+EVRTLEEYt1hNwhhjTKmsJmGMMaZUliSMMcaUypKEMcaYUlmSMMYYUypLEsYYY0plScIYY0ypLEkYY4wp1f8Hl9RyDAg2djMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.close('all')\n",
    "matplotlib.interactive(False)\n",
    "\n",
    "initsteady = (0.001, 0.001)\n",
    "fstart = 50\n",
    "deltaarr = np.linspace(-fstart*2*np.pi, fstart*2*np.pi, 2000)\n",
    "etasteady = 0.01\n",
    "losasteady = np.zeros(np.shape(deltaarr))\n",
    "losareal = np.zeros(np.shape(deltaarr))\n",
    "losaimag = np.zeros(np.shape(deltaarr))\n",
    "\n",
    "newinit = initsteady\n",
    "for i, delt in enumerate(deltaarr):\n",
    "\n",
    "    ar, ai = a_steady(etasteady, delt, gs, spins, gpar, gperp, kappa, newinit)\n",
    "    newinit = (ar, ai)\n",
    "    losasteady[i] = ai**2+ar**2\n",
    "    losareal[i] = ar\n",
    "    losaimag[i] = ai\n",
    "\n",
    "fig = plt.figure()\n",
    "plt.xlabel('$\\Delta f$ (MHz)')\n",
    "plt.ylabel('|A|² (dB)')\n",
    "plt.plot(deltaarr/2/np.pi, 10*np.log10(losasteady))\n",
    "\n",
    "ar, ai = cav_steady((fr*2*np.pi+deltaarr), fr*2*np.pi,  etasteady, kappa)\n",
    "losasteady = ai**2+ar**2\n",
    "\n",
    "\n",
    "plt.plot(deltaarr/2/np.pi, 10*np.log10(losasteady), '')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Power Scan Steady State Resonant"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 312,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEJCAYAAABhbdtlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VPW9//HXdyYrWYGEJSSEVVZZg6CAVbELWhVspXah2rq0tXax3mvVLuLPe7W3t9beK2rrdmu9VuRasbaKCijKorIbQGQPW1gSsi+TTGa+vz/OhCQkA0GSnAx5Px+P8cx8v2dmPgeTeeec75zvMdZaREREWuJxuwAREem8FBIiIhKWQkJERMJSSIiISFgKCRERCUshISIiYSkkREQkLIWEiIiEpZAQEZGwotwu4GylpaXZAQMGuF2GiEhEWb9+faG1Nv1060V8SAwYMIB169a5XYaISEQxxuxrzXo63CQiImEpJEREJCyFhIiIhKWQEBGRsBQSIiISlkJCRETCUkiIiEhYEX+exNm49957m7VNmzaNK664gpqaGu6///5m/TNmzGDGjBmUlZXxm9/8pln/zJkzmT59OgUFBTzyyCPN+mfNmsUFF1zAwYMHefzxx5v1z5kzh3HjxrFnzx6efvrpZv1z585lxIgRbNu2jeeff75Z/80338ygQYPYtGkTCxcubNZ/2223kZmZyZo1a3j11Veb9d9xxx2kp6ezYsUKFi9e3Kz/7rvvJjk5mWXLlrFs2bJm/ffddx+xsbG88cYbrFy5sln/gw8+CMCiRYtYu3Ztk76YmBjmzZsHwIIFC8jNzW3Sn5SUxD333APAc889x/bt25v09+zZkzvvvBOAp556ir179zbpz8jI4Pbbbwdg/vz55OfnN+kfOHAgt9xyCwAPP/wwx48fb9I/bNgwbrjhBgAeeughysvLm/SPGTOG66+/HoB58+ZRW1vbpH/SpEnMnj0b0M9eV/rZs9bSp28G3731+9QFLM88+QSHDx/GWosFsNA3qz9XXvctAkFLQmwUg9MTm9Xvli4dEiLSdVhr8fmDVNTU4fMHCATtiRu+IC98tI8af5AteUUUFFdhgxC0lqC1UBrg+8+vx1cXoDw3H39RKUHLiQ/6wJE6Fj64FH/AkpS3j+iqIiwWa533rjkU4JH8twHol3+I2NrSJrV9WOjl8cNOsE0fmsbzN03uyH+aUzK2fisiVE5OjtUZ1yJdi88foKC8hoKKGkqqaimu9FNcVRu6+SmudO6XVPkp99VR7vNTWesEQ2vFRHmIi/IQF+0N3TzERjnLuGgvsVEeYqI8RHk8RHkN0fVLr4doryHK6yHa4ywb9zdp9xi8jW/GkJYUy7is1Hb813MYY9Zba3NOt572JESkU6mtC3KopJr9RVXsL6riaKmPo2U+jpbXOPfLfZRU+Vt8rtdjSI2PpntCDN27RZPVoxvJcdEkxUWREOslMTaaxFgviXFRJMREkRgXRWKsc+sWE3UiAGK8Hjwe08Fb3jkpJESkwwWDlvzSanYerWDH0XL2FFSyr6iSA0XVHC6tpvEf/B4D6Umx9EmOo3/Pbkwa2J3eSXH0TokjLTGG7t1Ct4QYkmKj8BigtgKqi6G6BGrKIVALwQAE/RDwO49rfVBZDXU+8PtC6/ghWBdatw6sBU8UeLyhZZSzTm0V+KugthL81U6btYB1ljYYWgac17L1rxeEYLBRe7DRzTbcxzZta/K40Xp37YHouHb9f6WQEJF25fMH+ORwGbkHStiaX8aOo+XsOlZBZW0dcdSSRBUDEusYkhTgsvQ6+vX30y+6gl6eUlKDxcTXHMdUFTofspUWKhp9aJ74QA99CAdqwVfi3P8sPNHgjW4IBnA+1IN1oZvfWSemG0QnhJbxTpsxgAHjce4bj/M6UbHO0njAeJ3XNZ6Gpalf1t9wlo1f68T9xjfTUGM7UkiISNvx+8g/coitu/I4uH8vFcf2Qlk+fSlkKMe50FtJd6+PxOgq4rxVeGzow7wOKA7d6nmiICEdEntBtzTwxjR8+NZ/cNb/de8Jffh6YyAuFeK7N9xik5x2TxR4Q+t7YyAqzvmAj4qFqHgnHIwOMZ1MISEiZyYYhKLdkL8J8jdSe2A9dUV5RPlKiLE+MoCMxqt7PdTGp+NJzSQ6KQsTlwKxyc6Hd1yycz8upeEWm+yEQ3x38OhULrcpJEQkvIAfCnfA4Vw4kguHP8YezsXUOueI1BDDlmA2e+15VHmTSe7Ri169M8jOzKRvv2w8qVl4kvoQ5412eUPks1JIiEgDa51Q2PEW7HwbDqyBQA0Awag4DscN4aPAVD7w92eLHUR8xihmjMrgc+elM7Jvsr4RdA5SSIh0dX4f5K2AHW86wVCy32nvNYrq8d9ldVUmCw50Z1lBMqbKy/Shacwc3Yd/Hd6LXknt+80acZ9CQqQrqix09ha2vwG73wV/JUR3g0GXYKfewdaEyTy72c8/PzxMbV2Q8f1Tue/qflw5pi9pibFuVy8dSCEh0lUE/LB9MWx4DnYtAywkZcDY62HYTILZ03h7RwmPL99N7sE8EmOjuH5SFnOnZDO0d5Lb1YtLFBIi57qiPbDhL7DxBag85gTD9Dth5NXQZwwBC//MzWf+/DXsPFbBgJ7deGDWaGaP70dirD4iujr9BIicS+oHnvd/6NwOfOiEhPHA0C/CxBthyOXgjcJay3s7CnjojU/ZfrSc83on8l/Xj+PLYzLwagBaQhQSIucCa+HjF2HJfc7eAjgnoPWfAjk3wehrIbnh7IVth8t48I1trNhZSHbPbjz2jQnMHN1H306SZhQSIpGuLB/+8RPnm0lZU+Dy+5xlz8HNziCurg3wh6U7eHrlXhJjo/jVl0cyd0o2MVE6aU1appAQiVTWwqa/wpv3OHMWfek/4IJbw56lvGpXIXe/ksuBomq+lpPFPVcMJ7VbTAcXLZFGISESiSoL4dUfOHsP2VPh6kedPYcW1NQF+N1b23lqxV4GpSXw4i1TuHBwzw4uWCKVQkIk0hzaAC/NhcqC0+497Cmo4EcvbmRrfhlzp2TziytHEBfd/jOHyrlDISESSTY8D6/fCYm94aa3IGN82FWXbTvKTxZsItpreOrbOXx+ZO8OLFTOFQoJkUhQVwOL74L1f4ZBl8BXnoWElg8ZWWt54r3d/Odb2xmdkcKf5k4kIzW+I6uVc4hrIWGMmQ08CAy01sY1al8OxAI1wGZr7Y/cqVCkk6guhhfmwME1MO0OuOxXYS824w8E+fnLubyy8RBXjc3gt18ZQ3yMDi/JZ+fmnsT7wHhgSwt911lrD3ZwPSKdj68Unr8WjmyG656DUbPCr+oPcPtfN7B02zF+9vnz+NFlQzC6iI6cJddCwlp7HGjph9gCC4wxtcC/WWvf6ejaRDoFXxn871ecgPja8zBsZthVy31+bn5uHWvyinhg1mjmTsnuwELlXNYZxySus9YWGmOygKXGmBxrbbnbRYl0qJoKeOE655tMc547ZUBU1tTx7WfXsPlgKX/42jiuGdevAwuVc127hoQxxgusaqHrdWvtAy09x1pbGFoeMMZ8DAwBNp70urcCtwL079+/TWsWcV1tJfx1DhxcC199BkZcFXZVnz/Azc+tI/dgKY99YwJfGt2nAwuVrqBdQ8JaGwCmtHZ94xx7SrLWlhljkoDzgX0tvO6TwJMAOTk5to3KFekcXv4u7P8Arn0KRs0Ou5o/EOS2FzbwwZ7jPPK1sQoIaReuTdhijJlujFkKZBhjlhpjrsUJrXeNMSuBpcA8a22RWzWKdLjyI84V4qb/C5z/1bCrWWuZ99pW3vn0WGha78wOLFK6EjcHrlcAl7fQNbGjaxHpNPa+7yyHX3HK1f68Oo8XPtrP9z43SIPU0q409aNIZ7LnPYhLhT5jwq6yfPsxHvjnJ3xhZG9+/sXhHVicdEUKCZHOwlrYsxwGXhz2ZLn8kmp++tImzuudxCNfG6frP0i7U0iIdBZFe6DsIAz6XIvd/kCQH724EX9dkCe+NZEEXVpUOoB+ykQ6iz3LneXAS1rs/t3b21m/r5j//vp4BqYldFhZ0rVpT0Kks9j7HiT3a/G6EB/sPs6f3tvDNyb35+qxGS08WaR9KCREOoNg0Plm06BLml1ytKq2jp//LZfsnt345ZUjXClPui4dbhLpDI7kOrO9Dmw+HvHbN7ezv6iKBbdOoVuMfmWlY2lPQqQz2Puesxx4cZPmtXlF/Hl1HjdcmM2UQbrkqHQ8hYRIZ7BnOaQPh+S+J5rqAkF+9eoW+qXGc9eXdD6EuEMhIeK2uhrY90GzQ00vrtnPp0fK+eWVI/R1V3GNQkLEbQfXQl11k/MjiitreXjJDi4c1FMT94mrFBIibtvzHhgPZE890fT7JTso99Vx39UjdXU5cZVCQsRte5ZDxgSITwVgb2Elf12zn29O7s/wPsnu1iZdnkJCxE2+Mji0vsmhpkeW7CDG6+FHlw11sTARh0JCxE37VoMNnBi0/vRIGf/IzefGqQNIT4p1uTgRhYSIu/a+B1FxkDUZcPYiEmOi+N7Fg1wuTMShkBBxU95KJyCi49h+pJy3th7lu9MGktotxu3KRACFhIi7SvZB+jAA/vT+buKjvdx40QB3axJpRCEh4paaCvCVQnIGh0qqeW1TPtdfkEX3BO1FSOehkBBxS/lhZ5ncj6dX7AHg5ukai5DORSEh4payQwBUxfVi4doDXD02g36p8S4XJdKUQkLELWX5ALy530tlbYAbNBYhnZBCQsQtoT2Jpz72MTYzhbFZqS4XJNKcQkLELaWH8Mf2YFtBLd+aku12NSItUkiIuKUsnyP0ILVbNFfputXSSSkkRFxSV3KQndXJfGVCJnHRXrfLEWmRQkLEJXUlh8gPduerEzPdLkUkLIWEiBv81cT5Swgk9WNEX00HLp2XQkLEBXv27ARg4EBNBy6dm0JCxAWrN+YCMHb0SJcrETk1hYRIBwsGLbt27gAgpdcAd4sROQ2FhEgH27C/mHjfUedBcl93ixE5DYWESAd7ffNhMr1F2LhUiElwuxyRU1JIiHSgYNCyePMRRiZWYJL7uV2OyGkpJEQ60MYDxRwp8zEgugSSdZa1dH4KCZEO9HruEWK8HlL8BQoJiQgKCZEOYq3l7U+OcOmQZDxVBZCiM62l81NIiHSQnccqOFhczZUDjdOgPQmJAK6FhDHmD8aYD0O3uxu1TzDGrDLGrDbG3OhWfSJtbdm2YwBMTa91GhQSEgHc3JN4zFo7BbgIuMYYMzjU/ijwLeAS4MfGmO4u1SfSppZtO8rofsn0DBY6Dfp2k0QA10LCWrsztAwCASBgjIkFEqy1e621tcAKYJJbNYq0laLKWjbsL+ay4b1PXJFOexISCVwfkzDGzAV2W2vzgJ5ASaPuklDbyc+51RizzhizrqCgoGMKFTkL7+04RtDCjOG9oPQQxCZDbJLbZYmcVlR7vrgxxgusaqHrdWvtA8aYy4EbgKtC7UVASqP1UkJtTVhrnwSeBMjJybFtWrRIO1i27RhpibGc3y8FVh3SXoREjHYNCWttAJjSUp8xZjLwADDTWlsdWt9njKkyxvQHDgPTgPvbs0aR9hYMWlbuKuSy4b3weAyU5SskJGK4ebjpGSAJeNUYs9wYMzHU/hPgReA94HFrbbFbBYq0ha35ZZRU+Zk+NM1pUEhIBGnXPYlTsdaODtO+DpjaweWItJsVu5xxs6lD0iDgh4qjkKwT6SQyuD5wLXKuW7WrkOF9kuiVFAflRwCrPQmJGAoJkXbk8wdYm1fs7EWAc6gJdI6ERAyFhEg7WrO3iNq6INNOjEfoHAmJLAoJkXa0alch0V7D5IE9nAaFhEQYhYRIO/pwz3HGZ3WnW0zoOyJl+RCdAHEpp36iSCehkBBpJ1W1dWzJL2PSwEbTj5WFTqQzxr3CRM6AQkKknWzaX0IgaMkZ0KOhUedISIRRSIi0k3X7ijEGJvRvvCeRr282SURRSIi0k7V5RQzrnURKfLTTEKhzzpPQnoREEIWESDuoCwTZsK+YnAGN9iIqj4ENQIr2JCRyKCRE2sGnR8qprA0w6eTxCNDhJokors3dJHJO+u1vYdIk1sUMAGBidnd4911Yuxa+PMxZR4ebJIJoT0KkLU2aBHPmUPrG2/RNiaPfxg9hzhynvbT+RDrtSUjk0J6ESFu69FLsSy/x7auu5bwvzsH8fhEsXAiXXgpvLYWoOIjXZdslcmhPQqSNHRw3hb+Mm8nMRU/BD37gBAQ0nCOhE+kkgigkRNpY3sv/5Fsb36Dgp3fBE084YxIAxXshJcvd4kTOkEJCpC29+y4T/vX7/OtXf0GPh3/jHGqaMweWvg1Ht0LfsW5XKHJGNCYh0pbWruX+b/6aupyL8HqMc6hp4UJY9ipE1UK/CW5XKHJGtCch0oZKb7+DhclDmZTdaHD60kvhmjHO/Yzx7hQm8hmdNiSMMR5jTGzofq/2L0kkcq3fXwTQdFI/gPyNEN8DUrNdqErks2vNnsQi4GVjzHjg39q5HpGIti6vmCiPYVxWatOO/E3OXoS+2SQRpjUhcQS4BrgJ0J9BIqewLq+YUf1SiI/xNjTWVsGxT3SoSSJSa0Jii7U2CPwLUN7O9YhErJq6AJsOljQdjwA4usWZ2E8hIRHotN9ustY+Glr6gK+2e0UiEWrLoVJq64LNxyMObXCW+maTRKDThoQx5tetfK3l1tr3z7IekYi1Nq8YoOn04OAMWif2hqS+LlQlcnZac57Evla+VsnZFCIS6dblFTMwLYG0xNimHfkbIWOCBq0lIrUmJMqtta+0eyUiESwYtKzfV8TlI3o37agph8IdMPor7hQmcpZaM3D9y3avQiTC7SmsoLjK3/QiQwCHPwasBq0lYumMa5E2cMrxCFBISMRqzeGm4caY3BbaDWCttWPauCaRiLMur5ieCTEMTEto2nFogzPza2K6O4WJnKXWhMRe4Kr2LkQkkq3bV8TE7O6Ykwen8zdCxjh3ihJpA60JiVprbWu/4STS5Rwr87HveBXfmnzShARVRc41JCZ8253CRNpAa0LikDGmfyvWK7HWlp1tQSKRZt2+MOMRhzc5S41HSARrTUgkAH/GGYMIx4bW+cvZlyQSWdbmFREb5WFURkrTjhOD1jrcJJGrNdNyXNoRhYhEqjV7i5jQvzsxUSd9WfDQBugxCOK7t/xEkQigr8CKnIUyn59th8u4YGCP5p3104OLRDCFhMhZWL+vmKCFySeHRMUxKDvoTMchEsEUEiJnYc3eIqI8hvH9dRKdnJtaM3DdLowxfwCmhB6+aq39Tah9ORAL1ACbrbU/cqdCkdNbs7eIMZknXWQIQiFhoK/ONZXI5uaexGPW2inARcA1xpjBjfqus9ZeooCQzqy6NkDuwRIuGNizeeee5dB7FMQmdXhdIm3JtZCw1u4MLYNAIHQD5+u0C4wx7xhjLnOrPpHT2XigGH/ANh+PKM6D/R/AqNmu1CXSllw73FTPGDMX2G2tzQs1XWetLTTGZAFLjTE51tryk55zK3ArQP/+rTnPT6TtrdlbhDEw8eST6HIXOssxczq+KJE21q4hYYzxAqta6HrdWvuAMeZy4AYazQ1lrS0MLQ8YYz4GhgAbGz/ZWvsk8CRATk6ObafyRU5pzd4iRvZNJjkuuqHRWvh4AQyYDqn6A0YiX7uGhLU2QMPgdBPGmMnAA8BMa211qM0ASdbaMmNMEnA+rb8ynkiHqakLsGF/MddPOikIDq6Dot0w7Q53ChNpY24ebnomtHw1NHPmnUAu8K4xphqIBuZZa4tcqk8krA37SvD5g0wdkta0I3cBRMXByGvcKUykjbkWEtba0WG6JnZoISKfwYqdBXg9himDGg1a19XClr/B8CshLtm94kTakE6mE/kMVu4qZEL/VJIaj0fsWgLVxTDmevcKE2ljCgmRM1RcWcvmQ6VMG3LS1eY+fhES0mGwvrkt5w6FhMgZWrW7EGth2tBG4xHVxbDjLRj9VfC6/s1ykTajkBA5Qyt3FpIUF8XYzEbXj9i6CAK1MFaHmuTcopAQOQPWWlbsLOSiwT2J8jb69fn4JUgfDn3HulecSDtQSIicge1HyzlUUs3nzuvV0Fi0Bw58CGO+BuZUF3AUiTwKCZEzsGTrUQAuH9EoJD5eABhNwyHnJIWEyBlYsu0o47JS6ZUc5zT4q2HtMzD0C5CS6W5xIu1AISHSSkdKfeQeLOXzI3s3NOa+BFWFcNHt7hUm0o4UEiKttGSbc6jpC/UhEQzCB485g9UDprtYmUj7UUiItNKST44yoGc3hvRKdBp2vg2FO+DCH2nAWs5ZCgmRViipquWD3YV8YVQfTH0gfDAfkvvBqFnuFifSjhQSIq3w+ubD+AOWq8dmOA35GyFvBUz5AXijT/1kkQimkBBphb9vzGdIr0RGZYRmd109H2KSYMK33S1MpJ0pJERO42BxFWvyipg1LsM51FRywJmGY+INEJdy+hcQiWAKCZHT+PumfACuGdfPafjoj85y8vddqkik4ygkRE7BWsurGw8xMbs7WT26ga8U1j8Ho2ZDapbb5Ym0O4WEyCmszStm57EKrpsYOpv64wVQW66T56TLUEiInMLzH+4jKS6q4VDTjjchbRhkjHe3MJEOopAQATi+G/5yDVQcO9F0rNzHm1sOc93ELOJjvFBbBXmrYMjlLhYq0rEUEiIAq/4Ae5bD5v870fTSmgP4A5ZvTunvNOSthEANDJnhTo0iLlBIiFQVQW4oHLb9A4CaugAvfLSfaUPSGJwemoZj11KIiofsqS4VKtLxFBIiG5+HumoYeQ3s/xDKj/LKhkMcKfNx68WDGtbbtRQGTofoOPdqFelgCgnp2oIBWPM0ZE+DS+4BLIFPXuOJ5bsZm5nC9KFpznpFe6Bot8YjpMtRSEjXtn0xlO6Hybc616juOZTja19mf1EVt182tGEyv13LnKVCQroYhYR0bWv+BMmZMOxKMIbA8C/Ts3ANk3pZZgxvdInSXUuh+0DoOdi9WkVcoJCQruvYNtj7Pky6CbxRACwOXICXIPOG7cfjCe1F1NU462kvQroghYScO+pq4eWbYMdbrVv/oz+BNxYm3ABAabWfX33k5Zi3NyNLljest281+KsUEtIlKSTk3LHxL7DlZXjlVig/cup1q4ud61Offx0k9ARg/js7KfHV4Rl5NWbPu+Arc9bdtRS8Mc43m0S6GIWEnBv81fD+76DXKKjzwT/vAGvDr7/xBWfvYPKtAGzNL+XZVXlcNzGTtEnXQaC2YY9k1zLIvghiEjpgQ0Q6F4WEnBvWPgPlh+GK/4TLfgnb34DchS2vGwzA2qeg/4XQdyx1gSB3vZxL924x3HvFCMicBIl9YNvfofQgFGzToSbpshQSEvlqymHl72HQpTBgKky5DbImw+K7mh92qi6GF78OxXnOesCTK/awNb+MB64ZRWq3GPB4YMSXYefSE2dgM+TzHbtNIp2EQkIi30d/hKrjcNmvnMceL1zzuHPY6R8/bTjsdGQLPHkJ7H4HrnwYRlzFpgMlPLJkB1ec34eZ5/dteM0RVztnYb/3W+crsunDOnyzRDqDKLcLEDkr1cWw6lEYdgVkTmxoTxsCM34Nb93rDFAbD7z2Y4hPhe+8AVkXUFrl54cvbKB3chwPzR7T9HWzp0J8D6gughFXQf1JdSJdjPYkpONsfhn+dnPbvubq+VBTCpf+onnf5O9D1hQnHF65BfpNgO+9D1kXEAhafrZwE0fLfDz69fGkdItu+lxvFAy/0rmv8QjpwhQScvZ2LYVnZzrfMDqVDx5zpuKuq2mb960ogA+fgFHXQp/Rzfs9Xpj1OHTrAVN+CN/+OyQ6Z1H/++vbWPbpMe67aiTj+3dv+fUn3QyDL3NuIl2UDjfJ2TNe2L/aCYsRV7W8TukhyN/g3C87BD0GtbxeOL4yKNjujDPU1TjXddi6yBk3uPTe8M/rORh+tq3J4aLnVufx7Kq9fGfqAOZeOCD8czPGwdxFZ1anyDlGISFnb8B06NbT+dAOFxLb32i4X3rwzEKitgr+OA1K9jXvy/kupA099fMbBcRLa/dz32tb+fzI3vzyypGtr0Gki3ItJIwxPwNmAV5gN3CTtdZvjJkAPAoY4Elr7Z/dqlFayRvlhEPu/zmHnKLjm6+z/Q2ITYaaMickzsTqR52AuOq/nHCJinPOgI6Oh7TzWv0yC9cd4O5XNvO589KZ/43xeD0ajBY5HTfHJOZbay+21tZf5usLoeWjwLeAS4AfG2PCHDCWTmXUbPBXws4lzft8pbB3BYz9uvP4TEKi9JBzadGRs2DijTDwYsi6wDkUlD6sVd86stYy/52d3PVyLtOGpPGnuROJjfK2vgaRLsy1kLDW1gIYZ8J+D7DLGBMLJFhr94b6VwCT3KpRzkD2NOiW5hxyOtnOJRD0w+ivQEL6mYXEsv/nnCH9+fs/U1n+QJBfvLqF3729g9nj+/HMDZOIi1ZAiLSWq99uMsb8AtgB9AAOAD2BkkarlITaTn7ercaYdcaYdQUFBR1Sq5yGNwpGXg073nTGEBr79HVI6AWZOZCS2fqQOLgechfAhT+E7gPOuKT8kmq+9qcP+OtH+/nBJYP5/ZyxxETpC30iZ6Jdf2OMMV5jzIct3H4FYK39d+A8YC9wI1AEpDR6iZRQWxPW2iettTnW2pz09PT23AQ5E6NmO5Pm7Xy7oa2uxtmTGPYl5yuprQ0Ja+HNuyGxN0z/2RmXsuSTo1z53yvYfqSc+d8Yz8+/NLzhKnMi0mrtOnBtrQ0AU1rqM8bEWWt91lprjCkFqqy1PmNMlTGmP3AYmAZ8tuMM0vGypzqHk7YuglGznLa8FVBbDsO/7DxOyYJd7zghcKoP7S1/g4Nr4JrHIDap1SUcr6jh/n98wmsf5zO8TxKPf3MCg9ITz2KjRLo2N78C+7AxZhSh8QjgvlD7T4AXcb7d9Li1ttil+uRMebww8hpnGu7aSmdq7U9fh+gEGPg5Z52UTGeAu7rYOcmtJbVVsOQ+6DMGxn6jVW/tDwT560f7+cPSHVTU1PHTy4dy2yVDdHhJ5Cy5FhLW2h+GaV8HTG2pTyLAqNmw9mnnWgzGDx98AAAPD0lEQVQjZ8H2xTBkBkTHOf0pmc6y9GD4kPjgMSg7CF95ypmR9RSCQcuSbUf5jzc/ZU9BJRcO6sl9V49keJ/kNtwoka5LJ9NJ2+p/oTOOsHURpGY713ionwMJnBlVwQmJvmNafo1tf3dO0Mu+KOzb+ANB/pmbzx+X72H70XIGpyfwzA05XDa8l8YeRNqQQkLaVv0hpw1/gaS+zpQdQ7/Q0N94T6IlwSAU7oJJN7XYfaikmv9bd4CFaw+QX+rjvN6JPPK1sXx5TAbRXh1aEmlrCglpeyNnwZonncNO2Rc1PayUkO6cLV0WJiRK9zvzMTU6k7q02s+ybUd57eN83tvhfOV52pA07r9mNDOG98KjM6dF2o1CQtpe/ynO5T8rjjQ91ATOGENyv/B7EgU7ADgWl82yNft5c8sRVu8uxB+w9E2J4/ZLhzAnJ4usHt3aeSNEBBQS0h7qDzmt+ZNzMaCTnXSuhLWWo2U1rMkrwvvhcq4EPv/8UUqppH+Pbnx36kC+OLoP4zJTtdcg0sEUEtI+Lr0Hhs2E7tlNmoNBS3V8X7x57/PHpTvYfLCU3EOlFJQ715h4OG4bZd5UfnrVBVw0OI3zeidqIFrERQoJaXOBoKWgNp5DUeM49HE+eYWV7C6oYNexCvYUVPJ96+d271EeW7aNgb1SuXhoOmMyUxiXlcqYt/8b4xnFd6YOdHszRASFhLSStZaKmjqKKms5XllLUUUtxytrTtwvrKghv8RHfmk1R0p91AVtk+dndo9ncHoikwf2ZHLNGLxbFvHxHWPo1mtg4zeBwu3OwLeIdAoKiXOAtZa6oCUQdJa1dUFq6gL4/M6yxh/E5w9QU9ewbHzf5w9Q7qujosZPha+Oipq60OPQzec8rg0EW3z/+GgvPRNjyEiNZ9KAHvRNiSMjNZ5+qfH0TY0ju0cC8TGNZl7ddRi2QLfqw0CjkKgsdM7ETh/Wvv9gItJqXTYk/mvpTj5YOB8shP6DBZKyR5M6bDIBfy2H3nm+yXMskDBgLAmDxhPwVVK4amFDX+gP5/hBE4nNHEVdVSnl614FW//qjpjBk/H2HkpdeSG1uW826rPOawy6kLrUbIKlR/HufAeLEwI29Dol/aZQ1a0P0eX59Mj/yOlr9Pr5fS7CF5dGYsVBehVuaLbdh/peTE1sKknleaQfzz3R7jEGr8dQM2ImCcndiTu+g7gDm0jwGDI8Tl+Ux3D5179Hv/Qe5OV+xJZ1K4n2ePAEDJQCpXDfjfcRGxvLG2+8wZLXVjZ7/wcffBBSslh0NJO1Dz8Lif9s+LcJVDIvHkg7jwULFpCbm9vkuUlJSdxzzz0APPfcc2zfvr1Jf8+ePbnzzjsBeOqpp9i7d2+T/oyMDG6//XYA5s+fT35+fpP+gQMHcssttwDw8MMPc/z48Sb9w4YN44YbbgDgoYceory8vEn/mDFjuP766wGYN28etbW1TfonTZrE7NmzAbj33uaXXJ02bRpXXHEFNTU13H9/8ynLZsyYwYwZMygrK+M3v/lNs/6ZM2cyffp0CgoKeOSRR5r1z5o1iwsuuICDBw/y+OOPN+ufM2cO48aNY8+ePTz99NPN+ufOncuIESPYtm0bzz//fLP+m2++mUGDBrFp0yYWLlzYrP+2224jMzOTNWvW8Oqrrzbrv+OOO0hPT2fFihUsXry4Wf/dd99NcnIyy5YtY9myZc3677uv4Wdv5cowP3vAokWLWLt2bZO+mJgY5s2bB+D6z15n02VDYuWuAsoqG36J64dGjxVUUOMtxATrSPb5Q33mxAqlJT7qDpfh8VcTWxs4aY46g6/KD2U+PDW1mEDTQy4AtXVBogJBbND50D/xdGPwGOgWF0V8ajzWE09ZjBdjDIbQXHjGMHpIGsl9M6k+FuRgZdyJPoPBGLhi2kDSM/pTvD+aT1bvwBN6XRNa/uqbE8nK6seOzcm88/aBE/2e0Ibc8b0LQ7+osHjx7mb1Xz02w/lFPdSNnZ/1wj0p/ZxloOakf5xKiCe0J6Ep4EU6A2Nt8w+ySJKTk2PXrVvndhlypv5jAIy6Fr78+4a2xT+HDc/DvYdadcU5EfnsjDHrrbU5p1tP8xiIO1q6rkTBdkgbqoAQ6UQUEuKOlKzmIVG4U4PWIp2MQkLccfKeRE2FM59TozmbRMR9CglxR0om1JSCr9R5XOjM2aQ9CZHORSEh7jgxZfghZ1kfEmkKCZHORCEh7kjJcpb1h5wKtoMnCnpoOg6RzkQhIe6o35Oov65E4Q7oMQi80e7VJCLNKCTEHYm9nT2HxnsSGrQW6XQUEuIOjxeSMpyQCPiheK8GrUU6IYWEuKf+a7BFeyBYp0FrkU5IISHuScmE0gPOoSaAdB1uEulsFBLinpRMKMuHgk+dxz2HuluPiDSjkBD3pGQ6h5n2vg/JmRCb6HZFInIShYS4p/5cif0f6lCTSCelkBD31J8rEfRr0Fqkk1JIiHvqQwK0JyHSSSkkxD1xyRCb4tzXnoRIp6SQEHfVX8pUZ1uLdEoKCXFXSibEd4eENLcrEZEWRLldgHRxF/4QSg7okqUinZRCQtw16BK3KxCRU9DhJhERCUshISIiYSkkREQkLIWEiIiE5drAtTHmZ8AswAvsBm6y1vqNMcuBWKAG2Gyt/ZFbNYqIdHVu7knMt9ZebK2dGnr8hUZ911lrL1FAiIi4y7WQsNbWAhhjTKiOXfVdwAJjzDvGmMvcqk9ERFwekzDG/ALYAfQADoSar7PWTgNuAJ4wxiS18LxbjTHrjDHrCgoKOq5gEZEuxlhr2+/FjfECq1roet1a+0BoHQPMB7Zaax8/6fkLgYestRtP8R7lwPa2q9oVaUCh20WcJW1D53EubIe2of1lW2vTT7dSuw5cW2sDwJSW+owxcdZan7XWGmNKgapQYCRZa8tCexDnA/tO8zbbrbU5bVt5xzLGrNM2uO9c2AY4N7ZD29B5uDktx8PGmFE0jEfcF6rnXWNMNRANzLPWFrlYo4hIl+ZaSFhrfxima2KHFiIiImGdCyfTPel2AW1A29A5nAvbAOfGdmgbOol2HbgWEZHIdi7sSYiISDuJ6JAwxniNMb8zxiw1xiw3xowMtU8wxqwyxqw2xtzocpmnZIyZZ4zZFqp/eehrwxG1DfWMMc8ZY5Y2ehwx22CM+Zkx5v1QvX8xxkSH2iNpG/5gjPkwdLu7UXskbcPs0O+D76T2SNqGG0N1rjLGTHC7nrNmrY3YG/AD4NYW2lcBA4EYYAPQ3e1aT7EN84BvRfI2hOo9H/g7sDQStwGIaXT/L8CVEbgNQ0NLD/ABMDgCt6EnEAfsOqk9IrYB6B6qLyZU70q3azrbW0TvSQDXAdnGmHeNMfONMTHGmFggwVq71zpTf6wAJrlb5mndZYxZaYz5MUCEbsOvgQfrH0TaNtgWpomJwG3YGVoGgQAQiMBtOG6tPXkvIpK2YTKwwlpba63dCySG6o9YkR4S/YDD1tpLAR/wXZy/REoarVMSauusHgXGAp8HrjbGXEyEbYMx5hKc6VWONmqOqG2AFqeJibhtADDGzAV2W2vziNBtOEkkbUMPoLjR49JQW8Tq9Ne4PtXUHkAR8Gbo8ZvAtcCfgZRG66WE1nNNa6YnAaqNMa/gnCeyhgjaBmAqcD2Q2qi9iAjaBmvtA9bafzfGPIgzTcyNwLNE2DYYYy7HmffsqlB7xP1/aKG9023DKRTR9PegM9faKp0+JOypp/boBuTgnLGdg3Mc02eMqTLG9AcOA9OA+zuq3pacZhtSrbUlocMclwB/jqRtCE2f8hVgARAPjDLG/CL0gRsR2wAtTxMTSf8fAIwxk4EHgJnW2urQ+hG1DWHW73TbcAofAf8W+uJDX6DCWlvjck1nJaLPkzDGdAf+Bye5i4C51tpKY0wO8F+AAZ621j7rYpmnZIz5MzAMp9bl1tq7Q+0Rsw31jDEDcGq9PPQ4YrbBGPMY0HiamO9Z5yJYkbQNW0J36yeVu9Nauz7CtmE6zhQ9FwGrgcetta9E2DZ8F7gZ57IHP7HWrnO5pLMS0SEhIiLtK9IHrkVEpB0pJEREJCyFhIiIhKWQEBGRsBQSIiISlkJCRETCUkjIOc0YEzDGbDLGbDHG/F/oBEw36vipMebbLbQPqD+/wRhziTGmNFRvbmh2415n+D4LjDFD26puEYWEnOuqrbXjrLWjgVrg++39hvXTvTd6HIUzr9hfW/H0FaF6xwBrgXCX+Q3nCeCuM3yOSFgKCelKVgBD4MT1I7aEbj8Ntd3VaCbeR4wx74TuzzDG/G/o/heMMR8YYzaE9kwSQ+15xphfG2NW4sxO3NhlwAZrbV1o3YnGmI+NMR8QJgRC07QkEZoszjjXHXnOGPN26L2uNcb81hiz2RjzZmgaiPptvDwUTCJnTSEhXULoQ3MmsNkYMxH4Ds60zlOAW4wx44H3gemhp+TgTPMcjTNX0ApjTBrwS+Bya+0EYB3ws0Zv47PWTrPWLjjp7acC6xs9/h/gx9baC1sodboxZhOwH7gcZ5LBeoOBK4FrgP8F3rXWng9Uh9rrpwnfhTOzsMhZU0jIuS4+9KG7DueD9xmcD/1F1tpKa20F8ApOOKwHJoYmLazBuXBPTqhvBU6gjARWhV7zBiC70Xu9FKaGvkABgDEmBUi11r4X6nv+pHXrDzdl4YTJbxv1LbbW+oHNgJeGGZA3AwMarXcMyDjVP4pIa2mXVM511dbacY0bQodymglN6JeHs5exGsgFLsX5C35baLnEWvv1MO9VGa4GnKutgTNBXWsnTHsN+FujxzWhOoPGGL9tmHgtSNPf5bjQe4qcNe1JSFf0PjDLGNPNGJMAzMbZU6jv+5fQcgXOQPem0Afyh8BUY0z9uEY3Y8x5rXi/bYTGQqy1JUCpMWZaqO+bp3jeNGD3GW2Z4zxg62d4nkgz2pOQLsdauyE0RfuaUNPT1tqNofsrgF8AH4SmnfeF2rDWFhhjbgRebHRJyl/iXM3uVBbT9LDSd4BnjTFVwFsnrVs/JmFwrmp285lsmzGmN87e0+EzeZ5IOJoqXKQDGGMWAXfVX4e6Hd/nDqDMWvtMe76PdB063CTSMe7GGcBubyXAcx3wPtJFaE9CRETC0p6EiIiEpZAQEZGwFBIiIhKWQkJERMJSSIiISFj/H+/pjoa7sNphAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2013075576076097.8\n",
      "8.266555645753396e-14\n"
     ]
    }
   ],
   "source": [
    "plt.close('all')\n",
    "matplotlib.interactive(False)\n",
    "\n",
    "powers = np.arange(-30, 30, 0.1)\n",
    "detuning = 0\n",
    "drs = 10**(powers/20.)\n",
    "\n",
    "initsteady = (0.001, 0.001)\n",
    "losasteady_p = np.zeros(np.shape(powers))\n",
    "losareal_p = np.zeros(np.shape(powers))\n",
    "losaimag_p = np.zeros(np.shape(powers))\n",
    "\n",
    "newinit = (losareal[int(len(losareal)/2)], losaimag[int(len(losareal)/2)])\n",
    "for i, dr in enumerate(drs):\n",
    "\n",
    "    ar, ai = a_steady(dr, detuning, gs, spins, gpar, gperp, kappa, newinit)\n",
    "    losasteady_p[i] = ai**2+ar**2\n",
    "    losareal_p[i] = ar\n",
    "    losaimag_p[i] = ai\n",
    "    newinit = (ar, ai)\n",
    "\n",
    "fig = plt.figure()\n",
    "plt.xlabel('Power (dBm)')\n",
    "plt.ylabel('|T|²')\n",
    "\n",
    "data = scipy.io.loadmat(\n",
    "    \"/mnt/samba2/experiments/Quantum InterConnect/Cryogenics/Dilution Fridge/Measurements/20220223_DCR_n-diamond_33/vnaPowerScans/powerScan_S12_onres_7_mfile-v6\", simplify_cells=True)\n",
    "\n",
    "datax = data[\"measPars\"][\"vna\"][\"measPower\"]\n",
    "datay = np.array([trace['realpart']**2+trace['imagpart']\n",
    "                 ** 2 for trace in data['traces']])\n",
    "\n",
    "\n",
    "pows = powers-27\n",
    "\n",
    "\n",
    "plt.plot(pows, 10*np.log10(losasteady_p/drs**2))\n",
    "\n",
    "plt.plot(datax, 10*np.log10(datay)+64)\n",
    "\n",
    "\n",
    "val1 = 10*np.log10(1/kappa**2)\n",
    "val2 = 10*np.log10(1/(kappa+kappa*coop.real)**2)\n",
    "mean = (val1+val2)/2\n",
    "\n",
    "plt.plot([pows[0], pows[-1]], [val1, val1], '--', color=(0.3, 0.3, 0.3))\n",
    "plt.plot([pows[0], pows[-1]], [val2, val2], '--', color=(0.3, 0.3, 0.3))\n",
    "\n",
    "\n",
    "idx = find_nearest(10*np.log10(losasteady_p/drs**2), mean)\n",
    "\n",
    "plt.plot(pows[idx], 10*np.log10(losasteady_p[idx]/drs[idx]**2), 'rx')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "eta2 = (drs[idx]*1e3)**2/nbins*spin_pars.nreal\n",
    "print(eta2)\n",
    "k1 = 25e3\n",
    "print(eta2*scipy.constants.Planck*fr*1e6/(2*k1))\n",
    "\n",
    "# scipy.io.savemat(\"/home/andreas/wenz_daten.mat\", {'power_sim': pows, 'y_sim': 10*np.log10(\n",
    "#    losasteady_p/drs**2), 'power_data': datax, 'y_data': 10*np.log10(datay)+64})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 257,
   "metadata": {},
   "outputs": [],
   "source": [
    "%qtconsole"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.4"
  },
  "toc-autonumbering": false,
  "toc-showcode": true,
  "toc-showmarkdowntxt": true
 },
 "nbformat": 4,
 "nbformat_minor": 4
}