RPM_FW_N5N1_mesolve.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "46c31d44",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import linalg\n",
    "from scipy import integrate\n",
    "from scipy import sparse\n",
    "import matplotlib.pyplot as plt\n",
    "import qutip as qt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ef849016",
   "metadata": {},
   "outputs": [],
   "source": [
    "opstr2fun = {'x': lambda dim: qt.spin_Jx((dim-1)/2),\n",
    "             'y': lambda dim: qt.spin_Jy((dim-1)/2),\n",
    "             'z': lambda dim: qt.spin_Jz((dim-1)/2),\n",
    "             'p': lambda dim: qt.spin_Jp((dim-1)/2),\n",
    "             'm': lambda dim: qt.spin_Jm((dim-1)/2),\n",
    "             'i': qt.identity}\n",
    "def mkSpinOp(dims, specs):\n",
    "    ops = [qt.identity(d) for d in dims]\n",
    "    for ind, opstr in specs:\n",
    "        ops[ind] = ops[ind] * opstr2fun[opstr](dims[ind])\n",
    "    return qt.tensor(ops)\n",
    "def idOp(dims):\n",
    "    return mkSpinOp(dims, [])\n",
    "def zeroOp(dims):\n",
    "    d = np.prod(dims)\n",
    "    return qt.Qobj(sparse.csr_matrix((d, d), dtype=np.float64), \n",
    "                       dims=[list(dims)]*2, type=\"oper\", isherm=True) \n",
    "def mkH1(dims, ind, parvec):\n",
    "    axes = ['x', 'y', 'z']\n",
    "    components = [v * mkSpinOp(dims, [(ind,ax)]) for v, ax in zip(parvec, axes) if v!=0]\n",
    "    if components:\n",
    "        return sum(components)\n",
    "    else:\n",
    "        return zeroOp(dims)\n",
    "def mkH12(dims, ind1, ind2, parmat):\n",
    "    axes = ['x', 'y', 'z']\n",
    "    components = []\n",
    "    for i in range(3):\n",
    "        for j in range(3):\n",
    "            if parmat[i,j] != 0:\n",
    "                components.append(parmat[i,j] * mkSpinOp(dims, [(ind1,axes[i]), (ind2,axes[j])]))\n",
    "    if components:\n",
    "        return sum(components)\n",
    "    else:\n",
    "        return zeroOp(dims)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a1220d7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def point_dipole_dipole_coupling(r):\n",
    "\n",
    "    # H = - \\frac{ \\mu_0 }{ 4 \\pi } \\frac{ \\gamma_j \\gamma_k \\hbar^2}{ r_{jk}^3 }\n",
    "    #       \\left( 3 (\\bold{I}_j \\cdot \\bold{e}_{jk})  (\\bold{I}_k \\cdot \\bold{e}_{jk}) - \\bold{I}_j \\cdot \\bold{I}_k \\right)\n",
    "\n",
    "    dr3 = -4*np.pi*1e-7 * (2.0023193043617 * 9.27400968e-24)**2 / (4*np.pi*1e-30)/6.62606957e-34/1e6 # MHz * A^3\n",
    "\n",
    "    if np.isscalar(r):\n",
    "        # assume r is aligned with z\n",
    "        d = dr3 / r**3\n",
    "        A = np.diag([-d, -d, 2*d])\n",
    "    else:\n",
    "        norm_r = np.linalg.norm(r)\n",
    "        d = dr3 / norm_r**3\n",
    "        e = r / norm_r\n",
    "        A = d * (3 * e[:,np.newaxis] * e[np.newaxis,:] - np.eye(3))\n",
    "\n",
    "    return A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "169ceeb6",
   "metadata": {},
   "outputs": [],
   "source": [
    "g = 2.00231930436256 #electron g-factor\n",
    "beta = 9.274009994e-24 #bohr magneton\n",
    "hbar = 6.62607015e-34/(2*np.pi) #hbar\n",
    "\n",
    "# hfcs: mT -> Mrad/s\n",
    "\n",
    "N5 = g*beta/hbar*1e-9 * np.array([[-0.0994933,   0.00287092, 0],\n",
    "                                  [0.00287092, -0.0874862,  0.], \n",
    "                                  [0.,          0.,         1.75687]])\n",
    "\n",
    "N1 = g*beta/hbar*1e-9 * np.array([[-0.0529598, -0.0586562, 0.0460172], \n",
    "                                  [-0.0586562,  0.564443, -0.564764],\n",
    "                                  [0.0460172, -0.564764,  0.453074]])\n",
    "\n",
    "Is = [1, 1]\n",
    "As = [N5, N1]\n",
    "indE = [0, 1]\n",
    "\n",
    "Is = [1]\n",
    "As = [N5]\n",
    "indE = [0]\n",
    "\n",
    "rFW = np.array([  8.51061  , -14.251621 ,   6.5492562])\n",
    "Dee = point_dipole_dipole_coupling(rFW) * 2*np.pi # Mrad/s\n",
    "Dee = 0 * Dee\n",
    "\n",
    "b0 = 1.4 * 2*np.pi\n",
    "\n",
    "kb = 2. # 1/us\n",
    "kf = 1. # 1/us"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "89bdff6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "dims = [2, 2, *[round(2*I+1) for I in Is]]\n",
    "Ps = 1/4 * mkSpinOp(dims,[]) - mkH12(dims, 0, 1, np.identity(3)) # Ps = 1/4 - S1 . S2 = 1/4 - S1xS2x - S1yS2y ...\n",
    "one = mkSpinOp(dims,[])\n",
    "Pt = one - Ps\n",
    "Hhfc = sum(mkH12(dims, indE[i], i+2, As[i]) for i in range(len(As)))\n",
    "Hdee = mkH12(dims, 0, 1, Dee)\n",
    "K = kb/2 * Ps + kf/2 * one\n",
    "rho0 = Ps/Ps.tr()\n",
    "H0 = Hhfc + Hdee"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "4d15b843",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA+eUlEQVR4nO3deXxU9b3/8feZNXtCWMIWQkRBZHEJFQNi3cCil5YuV1pb0Ype+blVU7sgbbVcbmmrtWgV1KvU2rpwrWsVlbRWFpFWkCgCCrKFJRASIJOFTDIz5/dHSDSQQGYyM2eW1/PxmEeck3O+55NEH/P2ux3DNE1TAAAAMcRmdQEAAADHIqAAAICYQ0ABAAAxh4ACAABiDgEFAADEHAIKAACIOQQUAAAQcwgoAAAg5jisLqArAoGA9u7dq8zMTBmGYXU5AACgC0zTVG1trfr37y+bLbg+kbgIKHv37lV+fr7VZQAAgBDs2rVLAwcODOqauAgomZmZklp+wKysLIurAQAAXeHxeJSfn9/2OR6MuAgorcM6WVlZBBQAAOJMKNMzmCQLAABiDgEFAADEHAIKAACIOQQUAAAQcwgoAAAg5hBQAABAzCGgAACAmENAAQAAMYeAAgAAYk7QAWX58uWaMmWK+vfvL8Mw9PLLL5/0mmXLlqmoqEgpKSk65ZRT9Mgjj4RSKwAASBJBB5T6+nqdeeaZeuihh7p0/vbt23X55ZdrwoQJWrdune666y7ddttteuGFF4IuFgAAJIegn8UzefJkTZ48ucvnP/LIIxo0aJDmz58vSRo+fLjWrFmj++67T9/85jeDvT0AAEgCEZ+D8t5772nSpEntjl122WVas2aNmpubO7zG6/XK4/G0ewEAgOQR8YCyb98+5eXltTuWl5cnn8+nqqqqDq+ZN2+esrOz2175+fmRLhMAAMSQqKziOfYxy6Zpdni81axZs1RTU9P22rVrV8RrBAAAsSPoOSjB6tu3r/bt29fuWGVlpRwOh3r27NnhNW63W263O9KlAQCAGBXxHpTi4mKVlpa2O7Z06VKNGTNGTqcz0rcHAABxKOiAUldXp7KyMpWVlUlqWUZcVlam8vJySS3DM9OnT287f+bMmdq5c6dKSkq0adMmLVq0SE888YTuvPPO8PwEAAAg4QQ9xLNmzRpddNFFbe9LSkokSddcc42efPJJVVRUtIUVSSosLNSSJUt0xx136OGHH1b//v314IMPssQYAAB0yjBbZ6zGMI/Ho+zsbNXU1CgrK8vqcgAAQBd05/ObZ/EAAICYQ0ABAAAxh4ACAABiTlwFFDMQsLoEAAAQBXEVUAIEFAAAkkJcBRTTJKAAAJAM4iygxPyKaAAAEAYEFAAAEHPiKqAEAn6rSwAAAFEQVwFF9KAAAJAU4iqgMMQDAEByiLOAwioeAACSQVwFlAO1XjU2Mw8FAIBEF1cB5cpH3tVF971jdRkAACDC4iqgfM/+dw3yrLO6DAAAEGEOqwsIxq3OV5TlflXSnVaXAgAAIiiuelAAAEByIKAAAICYQ0ABAAAxh4ACAABiDgEFAADEHAIKAACIOXEZUHgmDwAAiS0uA8rza3fr4z01VpcBAAAiJK42amv1479+KMnQjl9fYXUpAAAgAuKyB8UQQzwAACSyuAwoNgIKAAAJLS4DCj0oAAAktrgMKPSgAACQ2OIyoPQxDsmtJqvLAAAAERKXq3hWum/XtkBfSV+3uhQAABABcdmDIkmn2PZZXQIAAIiQuA0oEjvKAgCQqOI8oFhdAQAAiIS4DigBEgoAAAkpzgOK1RUAAIBIiPOAQkIBACARxXVAAQAAiSmuA0ogELC6BAAAEAHxHVAY4gEAICHFdUAxA36rSwAAABEQ1wHlv1/bqD+v3ml1GQAAIMziOqC8/MEu/fzlj60uAwAAhFlcBxRDzEEBACARxXVAsYlVPAAAJKK4DiiG1QUAAICIiOuAQg8KAACJKc4DCnNQAABIRHEdUJgkCwBAYkqIgGKyoywAAAklrgNK6xDPY8u3aemGfRZXAwAAwsVhdQHd0RpQ5r3xiSRpx6+vsLIcAAAQJnHdg2KwigcAgIQU5wEFAAAkorgOKOyDAgBAYorrgPJb52O6wf5a23tW8wAAkBhCCigLFixQYWGhUlJSVFRUpBUrVpzw/Kefflpnnnmm0tLS1K9fP33/+99XdXV1SAV/0QX29ZrtfKbtvS9AQAEAIBEEHVAWL16s22+/XbNnz9a6des0YcIETZ48WeXl5R2ev3LlSk2fPl0zZszQhg0b9Pzzz+v999/X9ddf3+3ijxWgBwUAgIQQdEC5//77NWPGDF1//fUaPny45s+fr/z8fC1cuLDD81evXq3BgwfrtttuU2Fhoc4//3zdeOONWrNmTbeLPxb5BACAxBBUQGlqatLatWs1adKkdscnTZqkVatWdXjNuHHjtHv3bi1ZskSmaWr//v3661//qiuuCP+eJQQUAAASQ1ABpaqqSn6/X3l5ee2O5+Xlad++jndyHTdunJ5++mlNmzZNLpdLffv2VU5Ojv7whz90eh+v1yuPx9Pu1RUmz+YBACAhhDRJ1jDa70BimuZxx1pt3LhRt912m37xi19o7dq1evPNN7V9+3bNnDmz0/bnzZun7Ozstld+fn6X6mKOLAAAiSGogNKrVy/Z7fbjeksqKyuP61VpNW/ePI0fP14/+tGPNHr0aF122WVasGCBFi1apIqKig6vmTVrlmpqatpeu3bt6lJ9LDMGACAxBBVQXC6XioqKVFpa2u54aWmpxo0b1+E1DQ0Nstna38Zut0vqPFC43W5lZWW1e51YSzv0oAAAkBiCHuIpKSnR448/rkWLFmnTpk264447VF5e3jZkM2vWLE2fPr3t/ClTpujFF1/UwoULtW3bNr377ru67bbbdO6556p///5h+iGOJhMCCgAACSHopxlPmzZN1dXVmjNnjioqKjRy5EgtWbJEBQUFkqSKiop2e6Jce+21qq2t1UMPPaQf/vCHysnJ0cUXX6zf/OY3YfshbAooIBv7oAAAkCAMMw4mbng8HmVnZ6vmp5nKch8/GXdY45PyyqVRA7J10el9VDJxqAVVAgCAL2r7/K6p6cJ0jfbi+lk8rYyjYzvr99TowX9ssbgaAADQXQkRUGxMPgEAIKEkSEAJtHsfB6NWAADgBBIyoLDcGACA+JYgAaV9IqEHBQCA+JaQAYUeFAAA4ltCBJSr7aWaYPuo7T37oQAAEN+C3qgtFt3hfEGSNLjxGYsrAQAA4ZAQPSjHogcFAID4lpABZcE/t+qfn1ZaXQYAAAhRQgzxHOuhf34mSdrx6yssrgQAAIQiwXpQGNoBACARJFRAccrf7j37oQAAEJ8SKqAYx23YZlEhAACgWxI6oLCaBwCA+JRQAYVn8gAAkBgSKqAYx7ynBwUAgPiUUAHl2B4U8gkAAPEpwQJK+0TiJ6EAABCXEiqgHLsPCkM8AADEp4QKKMf2oJiBTk4EAAAxLaEDCj0oAADEp4QKKMfug/Kzlz/Ww0efywMAAOJHQgWUY3tQXl9foXvf+tSiagAAQKgSKqDwsEAAABJDQgWUY3tQAABAfEqKgBJgz3sAAOJKQgWUYyfJtmI1DwAA8SWhAorN6DiIsKMsAADxJaECiqGOd2YjnwAAEF8SLKB0zM8cFAAA4kpCBZRjn2bcamd1g6rqvFGuBgAAhMphdQHh1NkqnssfXCFJ2vHrK6JZDgAACFFC9aBcYVutUcY2q8sAAADdlFAB5YfOv+pv7p91+n2T2bIAAMSFhAooJ+NjsiwAAHEhqQJKs7/jSbQAACC2JGRA6Ww1T7OfHhQAAOJBQgaUzra851mCAADEh4QMKJ0+NJBJsgAAxIUEDSidbHkf5ToAAEBokiqg0IMCAEB8SNCA0nEQWb75gDbsrYlyNQAAIFgJtdV9q84CSsn/fSiJLe8BAIh1CdmDYnQyxAMAAOJDQgaUznpQWjX5CDAAAMSypAwo7CgLAEBsS8qAwmoeAABiW0IGlJPNQeGZgQAAxLaEDCgn60Ex6UEBACCmJWVAoQcFAIDYlqAB5WRDPCQUAABiWUIGFMNgkiwAAPEsIQPKyYZ45r62SYvfL49SNQAAIFgJGVAutJVpiLGn0++/+uFe/eSF9fIzGQUAgJiUkM/imeP8kyRpcOMzJzyv2R+Q3WaPRkkAACAICdmD0lXMRQEAIDaFFFAWLFigwsJCpaSkqKioSCtWrDjh+V6vV7Nnz1ZBQYHcbreGDBmiRYsWhVRwcFhuDABAPAp6iGfx4sW6/fbbtWDBAo0fP16PPvqoJk+erI0bN2rQoEEdXnPllVdq//79euKJJ3TqqaeqsrJSPp+v28WfjE2mAjI6/T5zUAAAiE1BB5T7779fM2bM0PXXXy9Jmj9/vt566y0tXLhQ8+bNO+78N998U8uWLdO2bduUm5srSRo8eHD3qu4iuwIKnKCTiB1lAQCITUEN8TQ1NWnt2rWaNGlSu+OTJk3SqlWrOrzm1Vdf1ZgxY/Tb3/5WAwYM0NChQ3XnnXfqyJEjnd7H6/XK4/G0e4Xi5Bu2hdQsAACIsKB6UKqqquT3+5WXl9fueF5envbt29fhNdu2bdPKlSuVkpKil156SVVVVbrpppt08ODBTuehzJs3T7/85S+DKa1DJwsoa3YcVP+cVI0ckN3tewEAgPAJaZKsYbSf12Ga5nHHWgUCARmGoaefflrnnnuuLr/8ct1///168sknO+1FmTVrlmpqatpeu3btCqVM2U8SUP7rz2v1H39YyVwUAABiTFA9KL169ZLdbj+ut6SysvK4XpVW/fr104ABA5Sd/XkvxfDhw2Wapnbv3q3TTjvtuGvcbrfcbncwpXXoZDvKtmI/FAAAYktQPSgul0tFRUUqLS1td7y0tFTjxo3r8Jrx48dr7969qqurazu2efNm2Ww2DRw4MISSu+5kQzyt2A8FAIDYEvQQT0lJiR5//HEtWrRImzZt0h133KHy8nLNnDlTUsvwzPTp09vOv+qqq9SzZ099//vf18aNG7V8+XL96Ec/0nXXXafU1NTw/SQdONkQT6tmPwEFAIBYEvQy42nTpqm6ulpz5sxRRUWFRo4cqSVLlqigoECSVFFRofLyzx/El5GRodLSUt16660aM2aMevbsqSuvvFJz584N30/Ria4O8TAHBQCA2GKYcbAZiMfjUXZ2tmp+mqksd+cbrx1rbOND2q/ck57377suUZ+slO6UCAAAjtH2+V1To6ysrKCuTehn8XR1iMdHDwoAADEloQOKzeha8CCeAAAQWxI7oHR1FQ89KAAAxJSEDihdHeL59/aDWr+7JsLVAACArgp6FU886WoPyg+f/1CStH3e5Z3uiAsAAKInoXtQurrMuBWTZQEAiA0JHVC6OsTTqskX3PkAACAyCChfQEABACA2JHRACXYBcbOfgAIAQCxI6IAS7BwUZqAAABAbEjqgGEFGjn9sqtTanQcjVA0AAOiqBF9mHFxAueul9ZKkHb++IhLlAACALkroHhSbAjKCnCgrSV6fPwLVAACArkrogPKi+x79yfmboK9r8BJQAACwUkIHFEm6wL4+6Gv8JtNlAQCwUsIHlBbBBY4AAQUAAEslRUAJdrJsgO1QAACwVFIEFIeCm1NCDwoAANZKioDS1acatyKgAABgraQIKEH3oDDEAwCApZIioAT70EB6UAAAsBYBpQM/f+VjPbpsa4SqAQAAJ5MkASW4IZ4VW6o0741PdKSJDdsAALBCUgQURwjb3UtSYzMBBQAAKyRFQLEboQUU5qIAAGCN5AgoQQ7xtGLLewAArJEUASXYZcatlm7Yr6Ub9oW5GgAAcDIOqwuIhmBX8bT62csfS5LW3zNJmSnOcJYEAABOICl6UIJ9Fs+x6r1MlgUAIJqSJKB0b2tYH1vLAgAQVUkRUIxu9qCQTwAAiK6kCCjdHeJhNQ8AANGVJAGle10g/gABBQCAaEqKgHK38yldZ3/D6jIAAEAXJUVAOcf2mX7h/HM3WqAHBQCAaEqKgNLKCHGohxEeAACiK6kCiku+kK777uP/0vf/+G+ZTJYFACAqkiqguNUc0nUHar3656cHVOsNLeAAAIDgJFVA6faGbX56UAAAiIakCiihPpOnldfHlvcAAERDUgWU7vag7D3cqEP1TWGqBgAAdCYpnmbcqrs9KN9cuEqStPVXl8tuM8JREgAA6EBS9aB0N6C0qmtksiwAAJGUVAHFZoQnoNQ3EVAAAIikpAoo4epB4dk8AABEFgElBAE2bAMAIKKSKqAYYXqmDh0oAABEVlIFFIZ4AACIDwSUEDDEAwBAZCVVQOnuRm2t6EEBACCykiqghKsH5fk1u/XaR3vD0hYAADgeASUEi97drlueWaeaI6E9HRkAAJxYUgWUcA3xtGpgwzYAACIiqQKKPUzLjFsxFwUAgMhIqoASrq3uW/2+dIvmvrYxrG0CAIAke5pxuId4XvhgtyTp6uICFfRMD2vbAAAks5B6UBYsWKDCwkKlpKSoqKhIK1as6NJ17777rhwOh84666xQbttt6WqUXf6wt1vvDX+bAAAks6ADyuLFi3X77bdr9uzZWrdunSZMmKDJkyervLz8hNfV1NRo+vTpuuSSS0Iutrsec/1eL7ruDnu7zEUBACC8gg4o999/v2bMmKHrr79ew4cP1/z585Wfn6+FCxee8Lobb7xRV111lYqLi0MuNhzOtG2TwjxZ1hcI79ARAADJLqiA0tTUpLVr12rSpEntjk+aNEmrVq3q9Lo//vGP2rp1q+6+u2u9F16vVx6Pp90rnFLlDWt79KAAABBeQQWUqqoq+f1+5eXltTuel5enffv2dXjNli1b9NOf/lRPP/20HI6uzcmdN2+esrOz2175+fnBlHlSjjBPljWMsDYHAEDSC2mSrHHMJ7JpmscdkyS/36+rrrpKv/zlLzV06NAutz9r1izV1NS0vXbt2hVKmZ1yiA3WAACIZUEtM+7Vq5fsdvtxvSWVlZXH9apIUm1trdasWaN169bplltukSQFAgGZpimHw6GlS5fq4osvPu46t9stt9sdTGlBCXcPyuMrtuu9rdW65eLTwtouAADJKqgeFJfLpaKiIpWWlrY7XlpaqnHjxh13flZWltavX6+ysrK218yZMzVs2DCVlZVp7Nix3as+ROFeavzGx/t039LNqvQ0hrVdAACSVdAbtZWUlOjqq6/WmDFjVFxcrMcee0zl5eWaOXOmpJbhmT179uipp56SzWbTyJEj213fp08fpaSkHHc8mhxGINwLeSRJXh+reQAACIegA8q0adNUXV2tOXPmqKKiQiNHjtSSJUtUUFAgSaqoqDjpnihWi8RmbZLU5CegAAAQDoZpmjG/Rtbj8Sg7O1s1P81Ulrv7S2Yu8d6rreaAMFTW3nfHDtKpfTL0/fGFYW8bAIB40/b5XVOjrKysoK5NqmfxtLKHeZJsq6f/1dJzNOXM/uqVEblJvgAAJLqkeppxK0eEhnhaVdc1RbR9AAASHQElAtj6HgCA7iGgREDsz+oBACC2JWVAidQcFAAAEB7JGVCMyAYUelAAAOiepAwotgj3oPhNU3GwehsAgJiVlAEl0kM8Vz76nr6+YJV8bNwGAEBIkjKgnGfbqFOMvRFrv8kXUNmuw9peVR+xewAAkMiSMqDc7HhVb7vvVEQeyPMFAUZ5AAAISVIGlFZuNUe0/V0HG/TJPk9E7wEAQCJK6oDilC+i7V//1Bp9Zf4KVdY2RvQ+AAAkGgJKFGytZC4KAADBSPKAEtkdZVux9T0AAMFJ6oDiMqLTg9LMcmMAAIKS1AElWkM8dKAAABCcpA4okX5oYKsd1fXauJfVPAAAdJXD6gKs5IpSD8rc1zdJkj74+UTlpruick8AAOJZUveguCK8D8qx9hw6EtX7AQAQr5I6oNijNMQDAACCE1cBZa96h7U9e4S3uj9WfZNP9d7oDCsBABDP4iqg7D/n9rC2Zzei24Py7cdWa8Tdb+lIEz03AACcSFwFlHD3d0RrFc+xNlawogcAgBOJq4AS7ohikzUblHib6UEBAOBE4iqgmM1NYW3PYVFACUR36gsAAHEnrgKKUnuEtTlW8QAAEJviKqCcdeE39I++12lZ5hVhac9uUQ/K/67Ypvve+lSmSVcKAAAdiauAYrfbdcnM3yv19Ilhac+qIZ5lmw/ooX9+xmRZAAA6EVcBpY0ZnqGZ+5wL9VPHs2FpKxR1jeyJAgBAR+IzoATCE1Bchl8zHX+L2lONj8VkWQAAOhaXAcWw28PaXorCuzqoq556b4fufuXjrs9F8dZKnorIFgUAQAyIy6cZp2WGdzVPippUq7SwttkVb3y8T5J02ci+Gjek1/EnfPBnqf6AlHuKVLVZevdBqalW+uFmKTMvytUCABA9cRlQhhdP0Y5VZ8jl86i/b3e323MbTeHfpjYIh+o7earyq7d0fHzPGun08KxkAgAgFsXlEI/N6dLgn76ngxf9NiztpVo0xNOq2X/MaiLTPPE8m4PbpT0fRLYoAAAsFJc9KK3MQHgmt1o1B6XVn97boQ17a3TX5cNlGIb0l29Ih3Z2fsHS2S1fbyuTcgujUiMAANEUlz0orQL+8AQUq1bxtFpXflj/u2K7ynYdbjmw9W3p4NaTX7j7fSlMvwMAAGJJXAcUBTqZuxEklxEbH/ID/nGL9NiFXb/gxRukh4qkI4ciVhMAAFaI64BipOWGpR1HTDyTx1Sfna9Je9cFd9mhHdKWv0ekIgAArBLXAWXUuZdqaZ/va1XKl7vVjsPiIR5JSpM39ItNa7bsBwAgUuI6oBg2mybdNF/p517VrXZcMRBQujcPhi1pAQCJJa4DSquAurezbCwM8XTrycqHdko73g1fMQAAWCyulxm3Cvi7FzCsXsXzY8dzOs+2MfQG3vlVy9cb/ikNOCc8RQEAYKGECCimr3v7mDgNa3tQbnK8Gp6GylcTUAAACSEhhnjOuGCq9ihPVWZ2SNff53xUP3f8OcxVWcBv7YZzAACES0IElLT0LOX9bKP2XvJgyG3McLwht8U7ynZb9RbmogAAEkJCBBRJcjgcLc+w6YYc1YWpmq4pMPbpVKP7Dztss+4v0pOXS/s3hK9NAAAskBBzUFqZ3RziSIniU41tCmip68dyR2IX271lUt6I8LcLAECUJEwPiiQ5B57VrevdCs/W+V27V1Nkwokk9kUBAMS7hAooZwwdptcvekMvDv99SNdH86nG3dr35GS6OdQFAIDVEiqgSNIVXx6nwf3zQro2YQIKPSgAgDiXcAFFkoxAaEMnKUb0AkpEd6+lBwUAEOcSMqA47aH9WNF6Jo9bTUqNZBhaeb/03Hclv/XPGAIAIBQJtYqnVf7gU0O6zhWVSbKmXnH9XKfbdkXuFod2tLw2viyN+lbk7gMAQIQkZA9K9sDhOvTludo54uagrotGD4pT/siGky+q2x+d+wAAEGYJGVAkqcdFt6rp9KlBXeOM2LLfzzmi+WBC5qIAAOJUwgYUSQr4vEGdH419UBwRXb1zDDOK9wIAIIxCCigLFixQYWGhUlJSVFRUpBUrVnR67osvvqiJEyeqd+/eysrKUnFxsd56662QCw6GX0ZQ50djiMceydU7xyKgAADiVNABZfHixbr99ts1e/ZsrVu3ThMmTNDkyZNVXl7e4fnLly/XxIkTtWTJEq1du1YXXXSRpkyZonXr1nW7+JMZNHysnvFdpGX+0V06/2v2d3W+bX3E6vmabaW+74hOOJMk/f1u6U9TJG90nzEEAEB3GaYZ3ESFsWPH6pxzztHChQvbjg0fPlxTp07VvHnzutTGiBEjNG3aNP3iF7/o0vkej0fZ2dmqqalRVlZWMOWqtrFZ28ve0eg3u76aZWjjn9QkZ1D3ORmHfPosZXpY2+yyKQ9IRddac28AQNLqzud3UD0oTU1NWrt2rSZNmtTu+KRJk7Rq1aoutREIBFRbW6vc3NxOz/F6vfJ4PO1eocpMccpuDy5s9FTo9+tMjurD3maXNTdad28AAEIQVECpqqqS3+9XXl77reTz8vK0b9++LrXxu9/9TvX19bryyis7PWfevHnKzs5ue+Xn5wdT5nH8/uDmYmQaDd26X0eis8dKJ5iLAgCIMyFNkjWM9pNPTdM87lhHnn32Wd1zzz1avHix+vTp0+l5s2bNUk1NTdtr167u7Rtyyqix2mnmab+Z06XzsyLQ22E3ojg59lgEFABAnAlqJ9levXrJbrcf11tSWVl5XK/KsRYvXqwZM2bo+eef16WXXnrCc91ut9xudzClnVBGerpss9Zr5ycfKO/lr5z8fONI2J+354zm6p1jmRbeGwCAEATVg+JyuVRUVKTS0tJ2x0tLSzVu3LhOr3v22Wd17bXX6plnntEVV1wRWqXdlJbilt3oWk9CmoLbP6UrIvpwwJOhBwUAEGeCfhZPSUmJrr76ao0ZM0bFxcV67LHHVF5erpkzZ0pqGZ7Zs2ePnnrqKUkt4WT69Ol64IEHdN5557X1vqSmpio7OzuMP8rJ+dP7dum8dCO8k0qvti/VhAguXz6pTX+T6iqliXOkICcMAwBghaADyrRp01RdXa05c+aooqJCI0eO1JIlS1RQUCBJqqioaLcnyqOPPiqfz6ebb75ZN9/8+bNxrrnmGj355JPd/wmCcPqQU7To9MeV49mkb+z9XafnpYa5B+W/nU+Gtb2g7Vnb8hp8vnS6NT1YAAAEI+h9UKzQnXXUHflk3bs6/ZXLO/3+b5q/rYX+r3b7PlLLzrFbU64OS1vddvl90rk3WF0FACBJRG0flERhO8kwR2oYh3ic0Xw44MkEmCwLAIgPSRlQ8vufeMXRbY6X9ajzfoVjKY+l+58ci9U8AIA4kZQBJbVXgerG/0R7T/tup+dcZl+jXmHYUdYdUz0oMVQLAAAnEPQk2USRMfEu7d38qbTl6U7P6WtUq8oMfaVRieP/VGh0bYfdqGj0SIfLpZxBVlcCAMAJJWUPSiub48RzUXobNSG3naV63eZ4WVPsq0NuI+xW3CfNHyVVb7W6EgAATiipA0pu3kA95Zuo1YHhHX9ftSG3naKmkK+NuM1vWl0BAAAnlNwBJd2lIdc+oprzf9Hh93sYoQcUtxHDAYXVPACAGJfUAUWSxp/aS4WDCzv8Xm53Akosrd45TsxvfQMASHJJH1AkqV/+KR0ez+3GKp5YDiiNzazmAQDENgKKpMxUtypu+EifTVzU7nhPw6NQehtKHP+nnzv/Eqbqwm/N9mqrSwAA4IQIKEf1G1Agd78R7Y5NtH+gp52/kqGuPw04Qw26zfGyzrNtCneJYbPqswP62sPvqqou/E9tBgAgHAgoX5Dbb7A+CeS3OzbevkFDjd1dbiNd4X0SciRcaV+mqyp+q//793arSwEAoEMElC9IT01Rxm0r9fGU19sdH26Ud3JFB22E8Tk+kTLYtl/THO+o6f0/655XN6ixmVU9AIDYQkA5xsBeOeo7cHC7Y8Nsu7p8fariZ9jEVrtbT67aoRc/2GN1KQAAtENA6UCv3v1U4+7X9v60Lgzx2OXXq67Zet09O5KlhZX96Nya2S+v1zcWvKuahthdeQQASC4ElI7Y7Mq+dYUOXnyvJOlS+zrlG/tPeEmBsV+jbfE1p8N2dIWSaUoflB/Wonfjq34AQOIioHQmo7eah17R9vYl191KO8EE2BzVRaOqsDKOWUL95KoduvmZD3SkiTkpAABrEVBOIDM3Tw/5viZJ6mV49A37ig7PMxTo1q6zVrEdE1BqjjTr9Y8qtOCdzyyqCACAFgSUE0hzOXTe9fP18ei7JEnX2t/qYE8UUy+67tHjrt9Fv8Bu+i/H63rL9ePjdsz9w9uf6bon32dOCgDAMgSUkxgzOFeFl96gWjNVp9r26nr7knbf763DOtsWvz0Ow2y79W3728cdf/uTSj22Yqv2HD5iQVUAgGRHQOmC9Kxcec7/uSRpluNZ3WB/TVmq11zHE/q9c4HF1XWf2+j42TwP/3Orxv/6bX2463B0CwIAJD0CShcNuPQmNZ59nWyGqdnOZ/RRyg36nuMfOt++werSus2uE0+KnfPaRj3773IFAjwFGQAQHQ6rC4gbhiH3lN/pud2Z+nrVo3Kbsb9jbFc5TvKsobU7D2ntzkOyGdK0Lw2KUlUAgGRGD0oQDJtN3755jvQf91tdSlidrAel1U9eWK/Lfr9clZ7ECWcAgNhEQAmB+6wrtSXvcqvLCBuH/HKpayt2Pt1fq/9ZskmH6psiXBUAIJkZpmnG/MQCj8ej7Oxs1dTUKCsry+py2tR+9p4y//IVq8sIiwbTrUu992qvenXpfJshLb6xWF8anBvhygAA8ao7n9/0oHSDe+BZ+ihQqGbTbnUp3ZZmeDXRvrbL5wdM6T8feU9fe2glQz4AgLAjoHSDKyVVKy/+q34/+lWrSwmLnkZN0Nd8uLtG9/xtg/6+cT+rfAAAYUNA6aabLjxVP5w6Tsv9o7TH7Gl1Od3iUsf7oZzMkvX7dP1Ta/TSuj1qbOY5PgCA7iOghIHdbtOZs/4p48ZlVpfSLaEGlFY/fP5DnTVnqT7afTg8BQEAkhYBJUyy05zqm9dfm+2nqko5VpcTErv8cnYzpDQ2B/SNBav0qyWbdLghwit9qrdK7z4oVX7S8vXAZulfj0mNnpNfCwCIaaziCbOA368jdTVK/32h1aWEpN50a2IQq3lO5D+LBqqgZ5qun3CKUpxhnEj80fPSsl9Ltfukprrjv597ijSgSLrkbiknP3z3BQAEpTuf3+wkG2Y2u13p2bkqG/oD7a8P6LI9f7C6pKCkG1593b5SD/undrut59fuliQ1NPk1aURfnZWf070GN7wkvT1Xqj7JwxkPbmt5ffqmdMqXpa/Mk3LYARcA4gk9KBEUCJiyzcmxuoygLfB9Vb/1fTvs7d73n2dqbGGu8nPTgruw0SPtWSP9+euh3bjfmdJ1SyVnSmjXAwBCwj4oMcpmM/TgaYt0f8+7rS4lKLaTPJsnVHc+/6Em/Pafeu2jvQoqF794Q+jhRJIqPpT+p6+0MTGWgwNAMiCgRNht3/2mSm4t0Uv+8fpX4HSry+kSe4QCSqtbnlmn0fcs1b+3HzzxsuSKj1qGdDa/GYa7mtL/XS39Y45Uuz8M7QEAIokhniiprvPK29Ss/g8OsLqUk3rNf54e812hj8whUbnfszecp/NOyZVhGO2/8b+XtAzthFv+WGnG0vC3CwBohyGeONAzw63+uRkqP7NEu4Zda3U5J/Qf9tV61f1znWNsjsr9vvO/q3Xa7De0dudB1TQ0t0xwfW9BZMKJJO36l7T6EampPjLtAwC6jVU8UTbo6y3zUV7/vVeHPHX6nvk3iyvq3CX2D/SBb2hU7uULmPrmwvdktxla3+NHSqvfHdkbvvkTaf/H0tceiux9AAAhIaBY5PLbF8ofMKX/7mF1KZ1yKLrb1tsU0ECzMvLhpNW6P0v9z5LO+h4rfAAgxjDEYxHDMOSw27TxO//We5e9bnU5Heru1vfB+pnjL1rmLonqPfX6D6Xnr43uPQEAJ0VAsdgZw4apuPh8PXXmX/SHgb+zupx2hhm7NNG2RlJ05lFf5wjHap0QbH5D8kV4W34AQFBYxRNjDi/6luyeXco8/InVpbT5f00/0BuBsRFrv8j4VNPs7+hKh4UPWzTs0refloZNtq4GAEgwrOJJIDnX/VWZP1itR/r/j2Ybt1ldjiRpvO3jiLb/gvuX1oYTSTL90rPfsbYGAEAbJsnGIsPQjTfcLF/A1EeLDii7aq0KvNFZ8tuRFKM5Iu3aFFAfHYpI26ExpT8USZffJw25yOpiACCp0YMSowzDkNNu0+gbHlHBrPf1QMpNeso30ZJabAooVY1hb3eu4wmtTrk17O12S/Vn0p+nWl0FACQ95qDEicZmvxqa/HI+N02u6k1yN1RE9f5Npl1XNM3TFnNg2NrckXJV2NoKuwvvkr50vZTe0+pKACBuMQclCaQ47cpNdynzuhfl/uHHeuJLr+mWnAVRu7/L8GuqfWVY2rrQVqYHnX8IS1sR886vpKWzra4CAJIWc1DijWFIdodmXDFBknRwQz817FyrHe+9pDG2zRGbLyJJPeUJSzv3Oh9Vb6MmLG1F1IfPSoPPl0Z8Q3KlWV0NACQVelDiXO6IizXw8h8pdcbf9PpX3lVtTssTk/ebOWG/17cd72iBc76yFNozbPqrSt+wLY+PcNLqlZtV+7dZ2u9pVByMhnau4aAUz/UDSDrMQUk0Pq/krdOL63brvY3bde/ea8J+i/ubv6UH/d8I+rpHnffrMnuEHgAYYVO9c3TaWefLJ7u+d16BDtY36dzBucpOc1pdWkvw2POBlD1AWvcXqfcw6dVbpYLx0s53pf5nS1v/KZ16iWRzSMO/Ku37SBp6mTSoWLK7JRv/rwIg/Lrz+U1ASXBN1eWqr/hEpf8qU9rulfoPs/v7jfzdf7Z+1HyjDqlrf4s8HdQ42wb93rWw2/e20r3NV+ph/9S296f0SpfXF9C3igbqswN1unbcYOWkOlXQM11OuyHDMCJXjGlK3lpp2W8kh1taEeIuxI6Uls3pHCnS8CnS6VeEt04ASY2Agi7xB0w1b1shb/k6vbPPpcNb3tM15ishtRUwDX2veZZWBUae9NxnnXNVbN8Y0n1izR1N/09/DxSpVp3PSTmld7p2Hzyi6cUFstsNTTqjrzbsrdHFp/fRgJzU7gWXxhqpZre06CtSwCc1N4TeVkd6D5dGTJUKL5B6FEpZ/cLbPoCkQkBBSAIBUzs2va/9G1bI4w0oY8db6mtWaUhge5eurzKzNLf5e3otcJ58Hcy3/q797/q6faXG2KzbZC4S/uK7RD/zzQj6ulSnXYYhXTairzbu9ejKL+XL7bDpy0N767PKOo0amK1eGe7jLwz4pc/+Idns0tPfksxAGH6KLnBlShf/TMo/VxpwTnTuCSChEFAQNof3blXVrk2qbwqosuwtpR76RMX+tfLKpTTD2+E16wODtdQ/RtXKVrWZqTNt23Su7ZOECyZf9IJ/ghb4vqqt5oCwtZmb7pLPH9A5BT20r6ZRXx3ZUwUVb2p0Vr3y11n4IElXhpQ3Qjr9P6TxsfH4BQDxgYCCiDnoqdP2HTvUK92hPRtWyldTobztLyrDaWiAd6vV5VlqbeA0zW3+nj40hygQxgVxaWrUjY7X1EO1mu4oDVu7YZHeu2XOyrn/JWX0lTJ6W10RgBgW9YCyYMEC3XvvvaqoqNCIESM0f/58TZgwodPzly1bppKSEm3YsEH9+/fXj3/8Y82cObPL9yOgxCDTVKPngD7cvF0HVz+tMzLqVVtdoTTfIdmOVGuwsd/qCqPmSd8k/Y/ve2ru5rZC2arTTx3PKteojY/VTj1Pk4qukQaMkQqKra6mPW+dJFPav7Fl5dLhnZIzVar4UPI3S2m5UlN9y0RjX6OU2qPlidYpWS1Dar2GtjxAMveUlu/3GdHSnj0GVm0BcSSqAWXx4sW6+uqrtWDBAo0fP16PPvqoHn/8cW3cuFGDBg067vzt27dr5MiRuuGGG3TjjTfq3Xff1U033aRnn31W3/zmN7t0TwJK/PD5AzrY0KRco1779uxUVZ1X6//9toYNO0P6cLGa0vspf8/r8rj6akTTR6pTirKMI1aX3W0+06b5vm/qlcA47TL7SOraRNjBRoVyVau5zj9qgHFA2UaYJ71GSU3KAFX2OV+fDfyGUnP7qc7VW6lOu3pnupXqtOtIs18pTrsyUxxy2m0yJDkdNrkdNjltNjX5AzIMye2wd+2GgUDL0uj9G1qWSX/295ZNDNc/3zJH5/Au6cjBlonE4ZKaK+UMktyZ0tCvSP6mlodKNtZI+edJhk1yuMJ3PyABRDWgjB07Vuecc44WLvx8yejw4cM1depUzZs377jzf/KTn+jVV1/Vpk2b2o7NnDlTH374od57770u3ZOAklj8AVM2v1fVlXvkc2bo/TX/0hmnFKjyveeUPqRYnrXPydWzUIXbntGGnIuUf3CVGty9NKp5vdWld8mmwCC9HThLHwaGaI/ZSz7ZZcpQipqUZTTIIb+KbRs00bZWp9j2WV1uRJQFhmhtYKhWB4ar3OyjzeZAmTJ0bHBrXdBkmpLLblOPdKccNptM01S626HBmX4ZkkY0fSSf7Brn+5cMw6ZTGz5UtnevnGZT1H+2zpiuDKnHYAVSc2UvvECm6ZdRMF5mc4OMwi/LbDwsI7NvlIo5ugzd7pI8e1p6iA580jJEV7lJyshr2QunR+HRrwUtx/ucIe3/WOp5asveOv3PaumF6je65fsF41vOzz9XqvykZfL0/o+lIRe3hMWC8VLVlpbraiukzH5SU52Ukt1SUySX3iMmRS2gNDU1KS0tTc8//7y+/vWvtx3/wQ9+oLKyMi1bdvweGxdccIHOPvtsPfDAA23HXnrpJV155ZVqaGiQ03l8l6nX65XX+/mETI/Ho/z8fAJKEjJNU81+U6bp1+7tm5U/cJB2b1qtzJ79tXXZ0+ozdpq2rP2H+hSO0pF3H1HKsEuU+umLqs8bowHbntfGvCkaUPmO9mWO1vk1f9N615k6u3mdttoGa0hgh+rMFGUY4X9SM45XaeaowsyVJB02M3REbjXLLkOm3PLJKZ/canlUQ45Rpz7GIfU0aq0sOaz8pqEqo4dMU6qy9VQgYGqrbbAMSW6HJH+TbM4UNfoNuRwONTc3Kc0pBXzNSrWb8vualWoPSAGf3HZTR/wO+VyZSmk+pCybVxnNB+W0G8puPiC3Yie4SVKz4ZJk6ogtUwHDpoPOfkoJ1Mnj7C23v161jp5yBLzy2jMkBeSzpbTM6zIM+eWQadhkypAhqfUD69ioYxhHN0s2JMM0ZZdPdrPlFTBscgUa5TOcSvHXq9mWojR/jXyGWyn+WjU4cpTdvF/1jh7Kajqgemeuspsr1eDIUpqvRg2OHGU0V6vO0UMZvkOqd+Qo3XdIdY5cpfsOq96RrTSfR7WOXKX5PS3H/YflcfRSSqBONc4+SvHXqcmWKkOSz+Zq+Wo45Dv6u3GYzQrIdvRnNGQaNhky5TccbT93wGg/jGxIMr+Y8iUZ+vyr0XaeefT7rb898+ixz89vPdZ23Dh67AsRoeXcz8857ljreW33MlXX2Kxz71oa0ud3UIPmVVVV8vv9ysvLa3c8Ly9P+/Z1/H+C+/bt6/B8n8+nqqoq9et3/D4L8+bN0y9/+ctgSkOCMgxDLochyaYhQ0dIkk4pmihJ6j34Vy3vTz+z5eRxE49edZMkyTTnqd/R/3iHmaYO1nh0dnaWGmurVeBO14aP/qVTRo3T+o1rNejUkfpk5Ys65ZyLVP7BUg0cfZG2v/OU+n1pqra+97L6jbxINf9+VqmjrlD9h68oc+RkNa5bLPPUifJtfF2u/HPk3vmOjN5DZe7boB4ZKTJq9yrT5pWvuVkuu1QXcEmuDDV7G+Rwp8nmrZHTblOqz6MDroHq1bRHn7pHaYR3ncqdp2hw81Z9Zhus081t+sw2WKcGdqjMGK5RgU+029ZfBeYe7TZ7a6BxQPvNHsozDqnc7KNBRmXb7++gmaFco67tvcdMU1YXhpEaTHenq7ZC1cc4rD7G4bC2GU/shqk8HZQMqa/Z8nW0+VnLN1sfodX6K28dmfKf5P2xI1hRWoEerNaeLqf/oCQpx1clServ7dqWBlHn3Xb0a+v7XZKk3Ka9x3xteap8z6Y9kqReTbslSXneHZKkvjrajhJjH6hQeLyhr8MJaVbfsRtNmaZ5ws2nOjq/o+OtZs2apZKSkrb3rT0oQDC++O+XYRjqmZMtSUrJ6iVJGvGlCyVJo85pmeA59vKWxwL0vvx6SVLe91pCcsGws1saGdNyviZ8teXrhd9q+Trlv47e5e4T1pNx9Oux/700+/wadHTuxdmS6hubNMTlUOORBg1NTdfhQ1Ua0qOnajyHNSIjW1UHq9U3J1tVB3arV88B2rt3m3IHnKq9+8rVs89A7dy1WT36F6pm92Zl9Rui7bs2Kid/pA7t+VTZfYdoS/lG5Qw4TXV7PlFG/6Gq37Ve6fmjdWjXJmUNOE0Hd29Rdt8h2r13o9L7n6HavZuV3neIGvZsUGrf0+Tdv0WuPqfJu+8TOXsPUdPBcjl65Kv54C7ZcgYocGCrbD2HyHdop2wZfZTadFABGWquq5bdZqjJdCiluUa2QJNsgSY57HY1yaFmOWXaXQoEAjLsDtnsDjU4c9XbVqvDypDRcFC19hylBWp1xEhXWopTtT6HegQOqcGWLrevXofNNGWoTkdMl1wOh7x+KcveLI/frjSHTQ0+UzmOZlU3+JSdla3DdQ3qnWbTgZoGZWRmyVN/RH0yHDpwuFbZPfvo0MFDysnO0KFDh5WdmSlPbY1ycvvo4OFDysruodqqvcrOy5enqkKp/U5X3d5P1LtHjuqrdsmV1UeNngOy5wyUr2qb0rJ7yV+zV67UDHkbapWWliG/Z4+cDqcamk05nS41eY/IbTfU5PPJ7XTqiN+Qy+VWg09KS3HL0ySluN2qazaV7fCpue6QHOnZamz0yp7WQ40NNXJm95P/8F75e54qW/Vnqut9tlKqPlZt7ihlH/pYrr5DZdu/Ub7ew+U++Km8uacr5fBm1fU6S5kHP5Yn+3Sl7l+r+t5nKqP6Y9XnnqGMQ5/oUM+zlV71oepzRyjz0AbV5wxTzsGPVNXjTGUfXK+67NOUefgT1WQNU2bNpzrUY7SyDm9SbUah0uu2qzGtn1Ib98vnzFLqkX1qdOcqpbFKDan9lOqtks+RKkdzvWTYZPM3yiZTpmnKZvplM088l+i4gQDDkN9wKmA45Dccspk+NdtS5DCb5bWlyRU4Iq8jUzazWQE5lNp8WHWuXkprqlJDSp5SG/fL4+yj9OaDanDlKr25WkecuUpranmf1nxQDc6W4/WuXkpvPqQjzhylHf2a2nxIR5w9lNZULa8jUyk+jxodWXL66uSzueX0H5HfsMtpNsswW+Zf+W0uGYGATEMyTKkljRoyAj6Zhk02MyDjC7+Hls4is/UXcLTL44t9Jp/3bRiG2g2xtvS62I7+c8v5Mj/vjTHb+l/a9Ysc/f4X+2aOPf7FY5JkU32jV9I9J/z7dSYmh3iOxRwUAADiT3c+v4PavMHlcqmoqEilpe33ZigtLdW4ceM6vKa4uPi485cuXaoxY8Z0KZwAAIDkE/TuUiUlJXr88ce1aNEibdq0SXfccYfKy8vb9jWZNWuWpk+f3nb+zJkztXPnTpWUlGjTpk1atGiRnnjiCd15553h+ykAAEBCCXoOyrRp01RdXa05c+aooqJCI0eO1JIlS1RQUCBJqqioUHl5edv5hYWFWrJkie644w49/PDD6t+/vx588MEu74ECAACSD1vdAwCAiIjaHBQAAIBoIKAAAICYQ0ABAAAxh4ACAABiDgEFAADEHAIKAACIOQQUAAAQcwgoAAAg5hBQAABAzAl6q3srtG526/F4LK4EAAB0Vevndiib1sdFQKmurpYk5efnW1wJAAAIVnV1tbKzs4O6Ji4CSm5uriSpvLw86B8Q4eXxeJSfn69du3bxXCSL8beIHfwtYgt/j9hRU1OjQYMGtX2OByMuAorN1jJVJjs7m3/ZYkRWVhZ/ixjB3yJ28LeILfw9Ykfr53hQ10SgDgAAgG4hoAAAgJgTFwHF7Xbr7rvvltvttrqUpMffInbwt4gd/C1iC3+P2NGdv4VhhrL2BwAAIILiogcFAAAkFwIKAACIOQQUAAAQcwgoAAAg5sR8QFmwYIEKCwuVkpKioqIirVixwuqSktLy5cs1ZcoU9e/fX4Zh6OWXX7a6pKQ1b948felLX1JmZqb69OmjqVOn6tNPP7W6rKS0cOFCjR49um1DsOLiYr3xxhtWlwW1/HdiGIZuv/12q0tJOvfcc48Mw2j36tu3b9DtxHRAWbx4sW6//XbNnj1b69at04QJEzR58mSVl5dbXVrSqa+v15lnnqmHHnrI6lKS3rJly3TzzTdr9erVKi0tlc/n06RJk1RfX291aUln4MCB+vWvf601a9ZozZo1uvjii/W1r31NGzZssLq0pPb+++/rscce0+jRo60uJWmNGDFCFRUVba/169cH3UZMLzMeO3aszjnnHC1cuLDt2PDhwzV16lTNmzfPwsqSm2EYeumllzR16lSrS4GkAwcOqE+fPlq2bJkuuOACq8tJerm5ubr33ns1Y8YMq0tJSnV1dTrnnHO0YMECzZ07V2eddZbmz59vdVlJ5Z577tHLL7+ssrKybrUTsz0oTU1NWrt2rSZNmtTu+KRJk7Rq1SqLqgJiT01NjSSF9DAuhI/f79dzzz2n+vp6FRcXW11O0rr55pt1xRVX6NJLL7W6lKS2ZcsW9e/fX4WFhfr2t7+tbdu2Bd1GzD4ssKqqSn6/X3l5ee2O5+Xlad++fRZVBcQW0zRVUlKi888/XyNHjrS6nKS0fv16FRcXq7GxURkZGXrppZd0xhlnWF1WUnruuef0wQcf6P3337e6lKQ2duxYPfXUUxo6dKj279+vuXPnaty4cdqwYYN69uzZ5XZiNqC0Mgyj3XvTNI87BiSrW265RR999JFWrlxpdSlJa9iwYSorK9Phw4f1wgsv6JprrtGyZcsIKVG2a9cu/eAHP9DSpUuVkpJidTlJbfLkyW3/PGrUKBUXF2vIkCH605/+pJKSki63E7MBpVevXrLb7cf1llRWVh7XqwIko1tvvVWvvvqqli9froEDB1pdTtJyuVw69dRTJUljxozR+++/rwceeECPPvqoxZUll7Vr16qyslJFRUVtx/x+v5YvX66HHnpIXq9XdrvdwgqTV3p6ukaNGqUtW7YEdV3MzkFxuVwqKipSaWlpu+OlpaUaN26cRVUB1jNNU7fccotefPFFvf322yosLLS6JHyBaZryer1Wl5F0LrnkEq1fv15lZWVtrzFjxui73/2uysrKCCcW8nq92rRpk/r16xfUdTHbgyJJJSUluvrqqzVmzBgVFxfrscceU3l5uWbOnGl1aUmnrq5On332Wdv77du3q6ysTLm5uRo0aJCFlSWfm2++Wc8884xeeeUVZWZmtvUyZmdnKzU11eLqkstdd92lyZMnKz8/X7W1tXruuef0zjvv6M0337S6tKSTmZl53Dys9PR09ezZk/lZUXbnnXdqypQpGjRokCorKzV37lx5PB5dc801QbUT0wFl2rRpqq6u1pw5c1RRUaGRI0dqyZIlKigosLq0pLNmzRpddNFFbe9bxxGvueYaPfnkkxZVlZxal91feOGF7Y7/8Y9/1LXXXhv9gpLY/v37dfXVV6uiokLZ2dkaPXq03nzzTU2cONHq0gDL7N69W9/5zndUVVWl3r1767zzztPq1auD/uyO6X1QAABAcorZOSgAACB5EVAAAEDMIaAAAICYQ0ABAAAxh4ACAABiDgEFAADEHAIKAACIOQQUAAAQcwgoAAAg5hBQAABAzCGgAACAmENAAQAAMef/AxQYao19zR0tAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[0.5523135378447771, 0.44558911812086077]"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "oris = [np.array([0, 0, 1]), np.array([1, 0, 0]) ]\n",
    "tmax = 20. / kf # us\n",
    "tlist = np.linspace(0, tmax, int(np.ceil(1000*tmax)))\n",
    "\n",
    "opt = qt.Options()\n",
    "opt.atol = 1e-10\n",
    "opt.rtol = 1e-8\n",
    "\n",
    "yields = []\n",
    "for ori in oris:\n",
    "    b0vec = ori * b0\n",
    "    Hzee = mkH1(dims, 0, b0vec) + mkH1(dims, 1, b0vec)\n",
    "    Heff = H0 + Hzee - 1j * K\n",
    "    Leff = -1j*qt.spre(Heff) + 1j*qt.spost(Heff.conj().trans())\n",
    "    sol = qt.mesolve(Leff, rho0, tlist, e_ops=[Ps], options=opt) \n",
    "    ps = sol.expect[0]\n",
    "    yr = kb * integrate.simps(ps, tlist)\n",
    "    yields.append(yr)\n",
    "    plt.plot(tlist, ps)\n",
    "\n",
    "ax = plt.gca()\n",
    "ax.set_xlim([0,5])\n",
    "plt.show()\n",
    "yields #change by adding electron dipole interaction and "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6787a415",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n"
   ]
  }
 ],
 "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}