Figure1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6d29583a-541e-4716-829c-a8856f24459f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import mrcfile \n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import fsc_utils as fsc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3031cda7-efbb-4c0f-ab4a-654f8d9a0c20",
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap_a = [plt.get_cmap('tab20c').colors[idx] for idx in [0, 4, 8, 12]]\n",
    "cmap_b = [plt.get_cmap('tab20c').colors[idx] for idx in [2, 6, 10, 12]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "8d51510d-04cf-4fb0-b71b-1fe165ef8439",
   "metadata": {},
   "outputs": [],
   "source": [
    "# You will need the 3-D volume and mask from: https://www.ebi.ac.uk/emdb/EMD-11657\n",
    "\n",
    "map_file = ''\n",
    "mask_file = ''\n",
    "    \n",
    "with mrcfile.open(map_file) as mrc:\n",
    "    x = mrc.data\n",
    "    voxel = mrc.voxel_size.x\n",
    "    mrc.close()\n",
    "    \n",
    "with mrcfile.open(mask_file) as mrc:\n",
    "    mask = mrc.data\n",
    "    voxel = mrc.voxel_size.x\n",
    "    mrc.close()    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c3ab250d-0aeb-41ea-90e8-839f8545d8d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "\n",
    "snr = 30\n",
    "B_signal = 100\n",
    "B_noise = 0\n",
    "\n",
    "x = x * mask\n",
    "\n",
    "x = fsc.fourier_downsample(x, 1.5)\n",
    "r = x.shape[0]//2\n",
    "\n",
    "y1 = fsc.generate_noisy_data(x, voxel, snr, B_signal, B_noise, return_noise=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "58cac62f-e993-414e-b96a-b282ee630ab8",
   "metadata": {},
   "outputs": [],
   "source": [
    "s1 = y1[:, :, ::2]\n",
    "s2 = y1[:, :, 1::2]\n",
    "S2 = fsc.phase_shift_3d(np.fft.fftshift(np.fft.fftn(s2)), 0.5, 0, 0)\n",
    "s2_shift = np.fft.ifftn(np.fft.ifftshift(S2)).real\n",
    "\n",
    "s3 = y1[:, ::2, :]\n",
    "s4 = y1[:, 1::2, :]\n",
    "S4 = fsc.phase_shift_3d(np.fft.fftshift(np.fft.fftn(s4)), 0, 0.5, 0)\n",
    "s4_shift = np.fft.ifftn(np.fft.ifftshift(S4)).real\n",
    "\n",
    "s5 = y1[::2, :, :]\n",
    "s6 = y1[1::2, :, :]\n",
    "S6 = fsc.phase_shift_3d(np.fft.fftshift(np.fft.fftn(s6)), 0, 0, 0.5)\n",
    "s6_shift = np.fft.ifftn(np.fft.ifftshift(S6)).real"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "226de8dd-d676-4e4f-ab80-3c9bb9ebd862",
   "metadata": {},
   "outputs": [],
   "source": [
    "c1 = fsc.two_volume_fsc(s1, s2_shift, r)\n",
    "c2 = fsc.two_volume_fsc(s3, s4_shift, r)\n",
    "c3 = fsc.two_volume_fsc(s5, s6_shift, r)\n",
    "\n",
    "c_avg = np.mean([c1, c2, c3], axis=0)\n",
    "\n",
    "freq = fsc.get_radial_spatial_frequencies(s1, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "917f3f37-2aab-4d72-a432-faf1254bb16b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAasAAAEcCAYAAAB53pugAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABFfUlEQVR4nO3dfVzV9f3/8ccLCETB6xQFCk20UKOQiSRTs7bUctZmViurdWEXq7bK1tzcVr/WWle7qG3pStvWhZldmmbWzCuKWMRXFHRiqSnmVV4hhijw/v1xDmdEfD4gnMPn8zm87rfbuck553M+n9d5fj7w8nMtxhiUUkopN4twugCllFKqKdqslFJKuZ42K6WUUq6nzUoppZTrabNSSinletqslFJKuZ42KxV2RGSWiPyqDabzDxH5bWs/KyJjRKQsuNW1nohcKyK5rfj8EhG5Jpg1qfYryukCVPshIluB3kBNvZcHGmO+COZ0jDE3B3N8qmkich8wwBhzVd1rxpjxzlWkwo2uWam2NtEYE1fvEdRGJSKRrfhs2P7nrbHvFs7fV4UfbVbKcSISIyJ/EpEv/I8/iUiM/71vbIoSESMiA/w//0NEnhKRt0XkCHBuw81zInKRiKwRkYMi8qGInFnvva0icq+IrAWONPwDLj5/FJE9InJIRNaKyJB6g3QTkcUiclhE8kXktHqfPV1E3hOR/SKyUUSmtDCfwfXGs1tEftGM3MaISJn/u+0CnhWR+0TkFRF5XkTKgWtFpIuIzBGRnSKyQ0R+a9XwReTPIrJdRMpF5BMR+bb/9XHAL4DLRKRCRIr8r68QkRv8P0eIyEwR+dyf5b9EpIv/vRT/PL1GRLaJyJci8suWZKXClzYr5Qa/BEYAZwHpwHBg5gl8/ofAg0A80LCxZQBzgZuAHsBsYGHdH3W/K4ALga7GmOoG4/4uMAoYCHQFLgP2Nfjs/UA34FN/HYhIJ+A94EWgl3+4v4nI4BP4XohIPPBv4B2gLzAAWOZ/u6ncEoDuwKnANP9rk4BX/N/lBeCfQLV/vGf7v+8NFuV87J9Wd//3WiAiHYwx7wC/A+b715bTG/nstf7HuUB/IA74S4NhcoBBwHnAr0XkDIs6VDukzUq1tTf8azgHReQN/2tXAv/PGLPHGLMX3x//qScwzjeNMR8YY2qNMUcbvHcjMNsYk2+MqTHG/BOowvdHvs4TxpjtxpjKRsZ9HF8TPB0QY8wGY8zOeu+/Zoz5j7/JvYDvjznARcBWY8yzxphqY0wh8Cow+QS+V914dhljHjfGHDXGHDbG5Pvfayq3WuA3xpiqet8tzxjzhjGmFugMjAd+aow5YozZA/wRuLyxQowxzxtj9vm/z+NADL7m0hxXAn8wxmw2xlQAM4DLG6zJ3m+MqTTGFAFF+BqwUoAeYKHa3sXGmH83eK0v8Hm955/7X2uu7TbvnQpcIyK313stusH4LT9vjHlfRP4C/BU4RUReB6YbY8r9g+yqN/hX+NYY6qabJSIH670fBTxn90UakQx8ZvFeU7ntbaR51/+upwInATtFpO61CCzyEJG78a119QUMvmbXs+mvYFlrFL4DbupYZamUrlkpV/gC3x/OOqf4XwM4AnSse0NEEhr5vN2tA7YDDxpjutZ7dDTGzGvm5zHGPGGMGQYMxrc58B674etNd2WD6cYZY25pxmcbjuc0i/fscoPGv1f917bjW8vsWa/GzsaYb2yq9O+fuheYAnQzxnQFDgF1Xa6p2zc0Vms1sLuJzykFaLNS7jAPmCkiJ4tIT+DXwPP+94qAwSJyloh0AO47wXE/DdwsIln+gyU6iciF/n1BTRKRb/k/exK+xnmUrx96b2URMFBEporISf7Ht1qwH2YRkCAiP/UfUBEvIln+9+xya5J/c+a7wOMi0tl/EMRpIjK6kcHj8TWXvUCUiPwa35pVnd1AiohY/U2ZB9wpIv1EJI7/7eNquI9QqUZps1Ju8FugAFgLrAMK/a9hjCkF/h++gww20eAAiqYYYwrw7bf6C3AA30EQ157AKDrja3gH8G262gc81ozpHsZ3sMLl+NYqdgEP49vP02z+8XwHmOgfxyZ8BymATW4n4Gp8m0XX4/uOrwB9GhluKbAEKMWXw1G+vrlwgf/ffSJS2Mjn5+LbBLoK2OL//O2NDKdUo0RvvqiUUsrtdM1KKaWU62mzUkop5XrarJRSSrmeNiullFKu59hJwT179jQpKSmtHs+xY8eIjo5ufUFhSvOxptnY03ysaTb2WpPPJ5988qUx5uSGrzvWrFJSUigoKGj1eFasWMGYMWNaX1CY0nysaTb2NB9rmo291uQjIp839rrnNwOmp+vlw+xoPtY0G3uajzXNxl4o8vF8szp8+LDTJbia5mNNs7Gn+VjTbOyFIh/PN6vNmzc7XYKraT7WNBt7mo81zcZeKPLxfLNSSikV/jzfrIJxRGE403ysaTb2NB9rmo29UOTTZLMSkbn+21AXW7wvIvKEiHzqv+V3RtCrtNG9e/e2nJznaD7WNBt7mo81zcZeKPJpzprVP4BxNu+PB1L9j2nAU60vq/kKCxu7wLOqo/lY02zsaT7WNBt7ocinyfOsjDGrRCTFZpBJwL+M7/LtH4lIVxHp0+DW30FXWVnJBRdcwKFDh+jSpQsA9e52GngeERFBVFQUMTExREdH06FDB2JjY4mNjSUq6n9fPzo6mk6dOtGxY0cSEhJISUmhX79+9OnTh4gIz28tVUopTwvGScGJfP2+NmX+177RrERkGr61L/r27cuKFSsA6N+/P/Hx8RQVFQHQo0cPBg8ezKpVq3xFRkWRk5NDYWEh5eW+u4mffvrprF69Ogjl24uLi+PUU0+lX79+nHbaaZx77rkcOXKE3r17ExkZSXZ2Nlu2bGHXLt8dudPS0qipqWHjxo0AJCYmkpSURH5+fmB8mZmZ5OXlUVVVBUBOTg6lpaXs2bMHgCFDhlBVVcWmTZsASE5Opnfv3oGTqDt37kxGRga5ublUV/vuXTdq1ChKSkrYt28f4DvP4fDhw1RWVrJixQpSUlLo3r174H883bp1Iz09nZUrV2KMQUQYPXo0RUVFHDhwAICMjAz279/P1q1bgZbNp8zMTHbv3s327b5FJDU1lZiYGIqLfVuVe/XqxcCBA8nN9d2mKiYmhuzsbAoKCqioqAAgKyuLsrIyduzYAcCgQYOIjIxk/fr1ACQkJNCvXz/y8vIAiI2NJSsri/z8fCorKwEanU8xMTGBZdDp+VR39JSb5lNsbGwgHyfnk5t+n+rmU1RUFOXl5a6YT275fao/nyIiIqisrGzRfLLSrPtZ+desFhljhjTy3mLgIWNMrv/5MuBnxphP7MaZmZlpWnMFi5qaGj744IPA84bfwxiDMYba2lqqq6s5duwYVVVVVFVVUVlZSWVlJTU1NYFhq6qqOHLkCBUVFXzxxRds3bqVLVu2sHfv3kanP3jwYObMmUNWVlaj7yullDpxIvKJMSaz4evBWLMqA5LrPU/Cd2fUkIqMjGTUqFGsXLmS0aMbuwt3cOzdu5d169axdu1aPv30UzZv3syaNWsoKSnhnHPO4ac//SkPPPAAHTt2DFkNrRHqfLxMs7Gn+VjTbOyFIp9gNKuFwG0i8hKQBRwK9f6q+kJ9p+OTTz6ZsWPHMnbs2MBrlZWV3HfffTz22GP84Q9/YNmyZSxZsoQ+fRq7G7iz9E7Q1jQbe5qPNc3GXijyac6h6/OAPGCQiJSJyPUicrOI3Owf5G1gM/Ap8DRwa9CrtK+vLScH+LbhPvzww+Tn55OamkpRURHnnHNOYHutmziRj1doNvY0H2uajb1Q5NOsfVah0Np9Vm6xd+9eLrroIv7zn//Qo0cPJkyYQElJSaBxde7cma5duzJlyhRmzJhBTEyMwxUrpZR7We2z8vwx2XVH0jjl5JNP5v3332fChAns27eP5557jsLCQo4cOcKRI0fYuXMnGzZs4P777yczM5NPPrE97iTonM7HzTQbe5qPNc3GXijy8Xyzqjss1EmdOnXijTfe4B//+AezZ88mNzeX/fv3U15ezvbt21m6dCkDBgyguLiYrKws/vznP7dZbW7Ix600G3uajzXNxl4o8nHs5ovh5qSTTuKaa675xuvx8fEkJSVRVFTEzJkz+dOf/sRPf/pTYmNjmTZtmgOVKqWU93h+n1V5eTmdO3cOQkVt46mnnuLWW29FRJg3bx6XXXZZSKfntXzakmZjT/OxptnYa00+YbvPav/+/U6XcEJuueUWHnzwQYwxTJ06lYULF4Z0el7Lpy1pNvY0H2uajb1Q5OP5ZlV36RIvmTFjBnfffTfHjx/n4osv5re//S21tbUhmZYX82krmo09zceaZmMvFPl4vll5kYjw6KOP8sADDwDwq1/9iu9///uB63+diKqqKj744APd4auUCmueb1b9+/d3uoQWERFmzpzJokWL6NKlC2+++Sbjxo0LXCiyKbt37+b+++/nlFNOIScnh1NOOYXp06cHLk5Zx6v5tAXNxp7mY02zsReKfDzfrOLj450uoVUmTJhAQUEBycnJ5OXlMXXqVNtNgjU1NTzyyCOceuqp3HfffezZs4ekpCQqKip4/PHH6devH4888kjgcidezyeUNBt7mo81zcZeKPLxfLMKh5PzBgwYwJIlS+jSpQuvvvoq99xzT6PDbdmyhXPPPZd7772XqqoqJk6cyPLly9m2bRsFBQVMmTKF6upq7r33Xn784x9TU1MTFvmEimZjT/OxptnY05OCw9jgwYN57bXXOOmkk/jDH/7AD3/4Q5YsWcKxY8fIy8vjpptuYujQoaxevZqEhATefvttFi5cyJgxYxARhg0bxvz585k/fz4xMTE89dRTfP/732fbtm2BW6EopZRXef48q3Xr1jF06NAgVOQOzz//PNdcc01gU2BMTMzXbkh26aWX8tRTT9GjRw/LcaxevZpJkyYFDrro2LEjZ555JkOHDmXIkCGceeaZ5OTkfO1Oye1RuC07wab5WNNs7LUmH6vzrDzfrGpra8PutvObNm1i3rx5zJ8/n/Xr15OQkMDUqVO55pprGDx4cLPGsWHDBmbOnMnHH38cuKtofWeffTZ///vfycz8xjLRboTjshNMmo81zcZea/IJ22a1YsUKxowZ0/qCXGrXrl307NmzxWtBK1asYOjQoRQVFVFcXExJSQlLlixh+/btREREcNttt/Hggw8SFxcX5MrdL9yXndbSfKxpNvZak0/YXsEi3CUkJLR6c12PHj0YO3Ysd9xxB7Nnz2bDhg1Mnz4dEeGJJ57grLPOIj8/P0gVK6VU8Hm+WbX3/S5NaSyfTp068eijj1JQUMCZZ57JZ599xsiRI3nggQfa1cEYuuzY03ysaTb2QpGP5zcDqtapqqril7/8JY8//jgAN954I7Nnz9Y7oSqlHBG2mwELCwudLsHVmsonJiaGxx57jHfffZcOHTrw9NNP8+STT7ZRdc7SZcee5mNNs7EXinw836xacj299qS5+XznO9/hH//4BwB33nkn7777bgircgddduxpPtY0G3uhyMfzzUoFz2WXXcbMmTOpra1lypQplJaWOl2SUkoBYbDPqqKiol0edt1cJ5pPbW0tkydP5vXXX2fEiBHk5uYSGRkZwgqdo8uOPc3HmmZjrzX5hO0+q927dztdgqudaD4RERHMnTuXvn378tFHH/HEE0+EqDLn6bJjT/OxptnYC0U+nm9WjV2dQf1PS/Lp2rUrs2fPBuCXv/wln376abDLcgVdduxpPtY0G3uhyMfzzUqFxkUXXcSVV15JZWUlN954Y8juZKyUUs3h+WaVmprqdAmu1pp8/vznP9OrVy9WrFjBU089FcSq3EGXHXuajzXNxl4o8vF8s4qJiXG6BFdrTT49evTgb3/7GwD33HMPmzZtClZZrqDLjj3Nx5pmYy8U+Xi+WRUXFztdgqu1Np8f/OAHXHXVVVRWVjJ16lSqq6uDVJnzdNmxp/lY02zshSIfzzcrFXpPPvkkSUlJ5Ofn8/DDDztdjlKqHfJ8s+rVq5fTJbhaMPLp2rUrzz77LAD33Xdf2FxqRpcde5qPNc3GXijy8fxJwdXV1XoFZBvBzOeOO+7gySefZMCAARQWFhIfHx+U8TpFlx17mo81zcZea/Jp1UnBIjJORDaKyKci8vNG3u8iIm+JSJGIlIjIj1pUZQvk5ua21aQ8KZj5PPLII5x55pl8+umn3HTTTTj1H51g0WXHnuZjTbOxF4p8mmxWIhIJ/BUYD6QBV4hIWoPBfgysN8akA2OAx0UkOsi1Kod16NCBl19+mU6dOjFv3jzmzp3rdElKqXaiOWtWw4FPjTGbjTHHgJeASQ2GMUC8+G6CFAfsB9rksDE9hNResPMZNGhQ4Jyr22+/nfXr1wd1/G1Jlx17mo81zcZeKPJpcp+ViEwGxhljbvA/nwpkGWNuqzdMPLAQOB2IBy4zxiy2G6/efNHbrr32Wv75z3/y7W9/m5UrV+rNGpVSQWG1z6o5e8Aa+yvUsMNdAKwBxgKnAe+JyGpjzNduaiIi04BpAH379mXFihUA9O/fn/j4eIqKigDfyaiDBw9m1apVviKjosjJyaGwsDBwn5TMzEx2797Nxo0b6dixI6mpqcTExASO7+/VqxcDBw4MbDuNiYkhOzubgoICKioqAMjKyqKsrIwdO3YAvrWGyMjIwNpCQkIC/fr1Iy8vD4DY2FiysrLIz8+nsrISgOzsbLZs2cKuXbsASEtLo6amho0bNwKQmJgYOOwbIC4ujszMTPLy8qiqqgIgJyeH0tJS9uzZA8CQIUOoqqoKnISbnJxM7969qWvunTt3JiMjg9zc3MB5T6NGjaKkpIR9+/YBkJ6ezuHDhykuLqZjx46kpKTQvXv3wJF83bp1Iz09nZUrV2KMQUQYPXo0RUVFHDhwAICMjAz279/P1q1bvzGffvCDH7Bw4UJWr17NzJkz+c53vmM7n+quFeam+bRx40ZqampcMZ82b94MEPT5BCf2+1R/Pm3dupXjx487Pp/c9PtUN59EhLPPPtsV88ktv0/155MxhhEjRrRoPlkyxtg+gGxgab3nM4AZDYZZDHy73vP3geF24x02bJgJhuXLlwdlPOEqlPnMmTPHACYhIcEcOnQoZNMJFV127Gk+1jQbe63JBygwjfSM5uyz+hhIFZF+/oMmLse3ya++bcB5ACLSGxgEbG7GuJWHXXvttWRlZbFr1y7uv/9+p8tRSoWxZp1nJSITgD8BkcBcY8yDInIzgDFmloj0Bf4B9MG32fD3xpjn7cYZrH1WlZWVxMbGtno84SrU+RQWFpKZmUlERARr164lLa3hgaLupcuOPc3HmmZjrzX5tOo8K2PM28aYgcaY04wxD/pfm2WMmeX/+QtjzHeNMUONMUOaalTBVFZW1laT8qRQ55ORkcG0adOoqanhkUceCem0gk2XHXuajzXNxl4o8vH85ZbqdhKqxrVFPvfccw8iwrx58wI7S71Alx17mo81zcZeKPLxfLNSzjvttNO48MILOXbsGE8//bTT5SilwpDnm9WgQYOcLsHV2iqfO+64A4CnnnoqcLiz2+myY0/zsabZ2AtFPp5vVpGRkU6X4Gptlc/555/P6aefzo4dO3j99dfbZJqtpcuOPc3HmmZjLxT5eL5ZeflyP22hrfIREW6//XbAd/8rL9Blx57mY02zsReKfDzfrJR7XH311XTu3Jnc3NywueeVUsodPN+sEhISnC7B1doyn7i4OK6//noAfvvb37bZdFtKlx17mo81zcZeKPLx/M0Xq6qq9ArINto6n507d3LaaadRWVlJQUEBw4YNa7NpnyhdduxpPtY0G3utyadVJwW7Wd3FFlXj2jqfPn36cNttvgvy/+pXv2rTaZ8oXXbsaT7WNBt7ocjH881Kuc/PfvYz4uLiWLJkCR988IHT5SilwoDnm5Ven8ueE/n07NmTO++8E3D32pUuO/Y0H2uajb1Q5OP5fVbKnQ4ePEi/fv04ePAgy5cvZ8yYMU6XpJTygLDdZ1V3cy/VOKfy6dq1Kz/5yU8AmDVrliM1NEWXHXuajzXNxl4o8vF8s6q7c6VqnJP5XH/99YgIr7/+euCOq26iy449zceaZmMvFPl4vlkp90pOTuaCCy7g2LFjPP98m901RikVhjy/z0rPd7DndD6vvvoqkydPZsiQIaxduxYRcayWhpzOxu00H2uajT09z6oRW7ZscboEV3M6n4kTJ3LyySdTXFzMxx9/7GgtDTmdjdtpPtY0G3uhyMfzzWrXrl1Ol+BqTucTHR3N1KlTAZgzZ46jtTTkdDZup/lY02zshSIfzzcr5X511wucN28eR44ccbgapZQXeb5ZpaWlOV2Cq7khn7S0NLKzszl8+DCvvfaa0+UEuCEbN9N8rGk29kKRj+ebVU1NjdMluJpb8rnyyisBXNWs3JKNW2k+1jQbe6HIx/PNauPGjU6X4Gpuyefiiy8G4J133nHNpkC3ZONWmo81zcZeKPLxfLNS3pCYmMiIESM4evQo77zzjtPlKKU8xvPNKjEx0ekSXM1N+Xz/+98H3LMp0E3ZuJHmY02zsReKfDzfrJKSkpwuwdXclM8ll1wCwKJFizh27JjD1bgrGzfSfKxpNvZCkY/nm5VeUNKem/IZMGAAQ4cOpby8nPfff9/pclyVjRtpPtY0G3t6IVvleW7bFKiU8gbPN6u4uDinS3A1t+VT16zeeOMNxw//dVs2bqP5WNNs7IUiH89fyFZ5izGG1NRUPvvsM1asWMHo0aOdLkkp5SJheyHbvLw8p0twNbflIyJMnjwZgBdeeMHRWtyWjdtoPtY0G3uhyMfzzaqqqsrpElzNjfnUXdj25ZdfdvQmdm7Mxk00H2uajb1Q5NOsZiUi40Rko4h8KiI/txhmjIisEZESEVkZ3DJVOBk8eDDDhg3j0KFDvPXWW06Xo5TygCb3WYlIJFAKfAcoAz4GrjDGrK83TFfgQ2CcMWabiPQyxuyxG2+w9llVV1cTFRXV6vGEK7fm88QTT/CTn/yECy+8kEWLFjlSg1uzcQvNx5pmY681+bRmn9Vw4FNjzGZjzDHgJWBSg2F+CLxmjNkG0FSjCqbS0tK2mpQnuTWfK664gqioKN555x12797tSA1uzcYtNB9rmo29UOTTnNaXCGyv97wMyGowzEDgJBFZAcQDfzbG/KvhiERkGjANoG/fvqxYsQKA/v37Ex8fT1FREQA9evRg8ODBrFq1yldkVBQ5OTkUFhZSXl4OQGZmJrt372bz5s3s2bOH1NRUYmJiKC4uBqBXr14MHDiQ3NxcAGJiYsjOzqagoICKigoAsrKyKCsrY8eOHQAMGjSIyMhI1q/3rTQmJCTQr1+/wM7C2NhYsrKyyM/PD+xryc7OZsuWLYGbjaWlpVFTUxO4kGNiYiJJSUmBk+Ti4uLIzMwkLy8vsF03JyeH0tJS9uzx9fghQ4ZQVVXFpk2bAEhOTqZ3797UrYl27tyZjIwMcnNzqa6uBmDUqFGUlJSwb98+ANLT0zl8+HAgn5SUFLp3705hYSEA3bp1Iz09nZUrV2KMQUQYPXo0RUVFHDhwAICMjAz279/P1q1bWz2ftm/3LUJ186mkpIThw4fz4Ycf8txzz5GZmdnm82nHjh2BzN0wnwBXzacvvvgi8F319+nr8+nYsWMkJSW5Yj658e9eVVUV/fr1a9F8smSMsX0AlwLP1Hs+FXiywTB/AT4COgE9gU3AQLvxDhs2zATD8uXLgzKecOXmfF555RUDmLPOOsuR6bs5GzfQfKxpNvZakw9QYBrpGc3ZDFgGJNd7ngR80cgw7xhjjhhjvgRWAenNGHerDRkypC0m41luzueiiy6iW7durFmzJvC/y7bk5mzcQPOxptnYC0U+zWlWHwOpItJPRKKBy4GFDYZ5E/i2iESJSEd8mwk3BLfUxukhpPbcnE9MTAxXXXUVAA899FCbT9/N2biB5mNNs7HnyKHrxphq4DZgKb4G9LIxpkREbhaRm/3DbADeAdYC/8G32bA46NU2om47tGqc2/O55557iI6O5uWXX2bdunVtOm23Z+M0zceaZmMvFPk06zwrY8zbxpiBxpjTjDEP+l+bZYyZVW+YR40xacaYIcaYPwW9UhWWkpOTmTZtGsYY7r//fqfLUUq5lOevYJGcnNz0QO2YF/KZMWMGHTp04NVXX2XNmjVtNl0vZOMkzceaZmMvFPl4vln17t3b6RJczQv59O3bl1tuuQWA3/zmN202XS9k4yTNx5pmYy8U+Xi+WemV2+15JZ97772Xjh07snDhwsC5K6HmlWycovlY02zshSIfzzcrFR569+7NtGnTAJg1a1YTQyul2hvPN6vOnTs7XYKreSmfumY1b968wNn2oeSlbJyg+VjTbOyFIh+9+aJylZEjR/Lhhx8yd+5cfvSjHzldjlKqjYXtzRfrroGlGue1fG644QYAnnnmmZBPy2vZtDXNx5pmYy8U+Xi+WdVdeFI1zmv5XHrppcTFxfHhhx+yYUNoL4LitWzamuZjTbOxF4p8PN+sVHiJi4vjiiuuAGDOnDkOV6OUcgvP77Oqra0lIkJ7rhUv5pOfn8+IESM4+eSTKSsrIzo6OiTT8WI2bUnzsabZ2GtNPmG7z6qkpMTpElzNi/kMHz6cIUOGsHfv3pDe9t6L2bQlzceaZmMvFPl4vlnV3RxNNc6L+YgI1113HQDPPvtsyKbjxWzakuZjTbOxF4p8PN+sVHi66qqriIqKYsmSJezcudPpcpRSDvN8s0pPb5N7PHqWV/M5+eSTueiii6itreW5554LyTS8mk1b0XysaTb2QpGP55vV4cOHnS7B1bycT91Jwc8++yyhOBDIy9m0Bc3HmmZjLxT5eL5Zbd682ekSXM3L+YwfP55evXrx3//+l/z8/KCP38vZtAXNx5pmYy8U+Xi+WanwddJJJ3H11VcDoT3QQinlfp5vVikpKU6X4Gpez6duU+BLL73EV199FdRxez2bUNN8rGk29kKRj+ebVffu3Z0uwdW8nk9aWhrDhw+nvLychQsXBnXcXs8m1DQfa5qNvVDk4/lm1VY36vOqcMjnyiuvBGD+/PlBHW84ZBNKmo81zcZeKPLxfLNS4W/y5MmICG+//TaHDh1yuhyllAM836y6devmdAmuFg759O3bl1GjRnHs2DHefPPNoI03HLIJJc3HmmZjLxT5eL5Z6cl59sIln8suuwwI7qbAcMkmVDQfa5qNPT0puBErV650ugRXC5d8fvCDHxAZGcm7777L/v37gzLOcMkmVDQfa5qNvVDk4/lm5dQtTrwiXPLp1asXY8eOpbq6mtdffz0o4wyXbEJF87Gm2dgLRT6eb1Yi4nQJrhZO+dRtCnzppZeCMr5wyiYUNB9rmo29UOTj+ZsvqvZj//799O7dm9raWnbu3EmvXr2cLkkpFWRhe/PFoqIip0twtXDKp3v37owbN47a2lqefvrpVo8vnLIJBc3HmmZjLxT5eL5ZHThwwOkSXC3c8rnjjjsAePLJJzl69GirxhVu2QSb5mNNs7EXinw836xU+3L++edz1llnsXv3bp5//nmny1FKtZFmNSsRGSciG0XkUxH5uc1w3xKRGhGZHLwS7WVkZLTVpDwp3PIREaZPnw7AY489Rm1tbYvHFW7ZBJvmY02zsReKfJpsViISCfwVGA+kAVeISJrFcA8DS4NdpJ1gnXMTrsIxnylTppCcnMzGjRtZtGhRi8cTjtkEk+ZjTbOxF4p8mrNmNRz41Biz2RhzDHgJmNTIcLcDrwJ7glhfk7Zu3dqWk/OccMznpJNO4s477wTg0UcfbfF4wjGbYNJ8rGk29kKRT1QzhkkEttd7XgZk1R9ARBKBS4CxwLesRiQi04Bp4Lve24oVKwDo378/8fHxgSNIevToweDBg1m1apWvyKgocnJyKCwspLy8HIDMzEx2795NRUUFK1asIDU1lZiYGIqLiwHfSaQDBw4kNzcXgJiYGLKzsykoKKCiogKArKwsysrK2LFjBwCDBg0iMjKS9evXA5CQkEC/fv3Iy8sDIDY2lqysLPLz86msrAQgOzubLVu2sGvXLsB3S4uamho2btzoCy8xkaSkpMCdbuPi4sjMzCQvL4+qqioAcnJyKC0tZc8eX58fMmQIVVVVbNq0CYDk5GR69+5N3aH+nTt3JiMjg9zcXKqrqwEYNWoUJSUl7Nu3D/Bd7uTw4cOBfFJSUujevXvgasjdunUjPT2dlStXYoxBRBg9ejRFRUWBnaMZGRns378/sOC1Zj5t3+5bhII1nyZNmsSvfvUrcnNzeemll7jkkktOeD4dP348sAw6PZ/q7qzqpvlUXV0dyEd/n74+n44dO0Z5ebkr5pMb/+5VVVVRWVnZovlkyRhj+wAuBZ6p93wq8GSDYRYAI/w//wOY3NR4hw0bZoLh888/D8p4wlU453PvvfcawFx66aUt+nw4ZxMMmo81zcZea/IBCkwjPaM5mwHLgOR6z5OALxoMkwm8JCJbgcnA30Tk4maMu9Xi4+PbYjKeFc753H777URFRfHqq6+2aLNDOGcTDJqPNc3GXijyaU6z+hhIFZF+IhINXA587Zatxph+xpgUY0wK8ApwqzHmjWAX2xg9Oc9eOOeTmJjI5ZdfTm1tLX/+859P+PPhnE0waD7WNBt7jpwUbIypBm7Dd5TfBuBlY0yJiNwsIjcHvSKlTsBdd90FwDPPPKM3ZlQqjDXrPCtjzNvGmIHGmNOMMQ/6X5tljJnVyLDXGmNeCXahVnr06NFWk/KkcM/n7LPP5txzz6WiooJnnnnmhD4b7tm0luZjTbOxF4p8PH8h29raWiIi9EIcVtpDPosWLWLixIkkJyezefNmoqKac5Br+8imNTQfa5qNvdbkE7YXsq07zFM1rj3kM2HCBAYNGsT27duZPXt2sz/XHrJpDc3HmmZjLxT5eL5ZKRUREcFDDz0EwMyZMwPnbSilwofnm1VzN/m0V+0ln4svvphx48Zx8OBB7r333mZ9pr1k01KajzXNxl4o8vH8Piul6mzatIkhQ4Zw7NgxcnNzGTlypNMlKaVOUNjus6q73IlqXHvKJzU1lZ/97GcA3HrrrYFL51hpT9m0hOZjTbOxF4p8PN+s6q6ZpRrX3vKZMWMGKSkprF27lldffdV22PaWzYnSfKxpNvZCkY/nm5VS9XXs2DGwdvXXv/7V4WqUUsHi+X1WFRUVxMXFBaGi8NQe86moqCAxMZHy8nLWrFlDenq65XDtLZsToflY02zstSafsN1ntXv3bqdLcLX2mE9cXBzXXnstAH/5y18sh2uP2ZwIzceaZmMvFPl4vlnV3ddFNa695vPjH/8YgBdeeMHyrqXtNZvm0nysaTb2QpGPniygwtLAgQO54IILWLp0KXPnzmXs2LE8+uijbNiwgWHDhjFy5Eg6duzodJlKqWbyfLNKTU11ugRXa8/53HbbbSxdupSZM2d+7Q6kRUVFzJ07l6ioKMrLy5k2bZqDVbpXe152mqLZ2AtFPp7fDBgTE+N0Ca7WnvMZP348/fv3p6qqiri4OO6++26WL1/O448/zkUXXUR1dTU33XQTd911FzU1NU6X6zrtedlpimZjLxT5eL5ZFRcXO12Cq7XnfCIjI1m4cCFPP/0027Zt47HHHmPMmDHcddddvPXWW9xzzz1ERUXxxz/+kYsvvvhra1+qfS87TdFs7IUiH883K6XsDB48mBtuuIFu3bp9470JEybw3nvv0b17dxYtWsQtt9yCU6dyKKXseb5Z9erVy+kSXE3zsdarVy/GjBnDe++9R2xsLM8++yxPPPGE02W5hi471jQbe6HIx/MnBVdXV+sVkG1oPtbqZzN//nwuv/xyIiIiWLJkCd/97ncdrs55uuxY02zstSafsD0pODc31+kSXE3zsVY/m8suu4yZM2dSW1vLpZdeyuLFix2szB102bGm2dgLRT6eb1ZKBcv999/PZZddRnl5ORdddBG/+MUvmrxyu1KqbXi+WekhpPY0H2sNs4mIiODFF1/kd7/7XeDuwyNHjmTOnDns27fPoSqdo8uONc3GXijy8fw+K6VCYeXKlVx++eXs2rUL8B0Gn52dzRlnnMGgQYNITEwkJiaGmJgYEhMTOfPMMxERh6tWyvus9ll5vlkVFBSQmfmN76X8NB9rTWVz8OBBFixYwCuvvMKyZctsTxzu06cPEyZM4Pzzzyc9PZ3U1FQqKyt5//33effdd+nSpQt33303PXr0CMVXCQlddqxpNvZak49Vs/L84SwVFRVOl+Bqmo+1prLp2rUrN954IzfeeCP79u3jk08+YePGjWzcuJHdu3dz/PhxqqqqWLduHTt27GDOnDnMmTMH8G0Gqa2t5fjx44HxzZ49m4ceeojrr7+eyMjIkH63YNBlx5pmYy8U+Xi+WSnVFnr06MF3v/vdRg9pN8awdu1aFi9ezEcffcTatWv5/PPPiYiI4JxzzuG73/0uq1at4v333+emm27ihRde4J133iE2NtaBb6KUN3l+M2BlZaX+0tvQfKyFMptDhw4hInTu3BnwNbQFCxbwk5/8hF27dnH99dfzzDPPhGTawaLLjjXNxl5r8gnb86zKysqcLsHVNB9rocymS5cugUYFICJMmTKFJUuW0KFDB+bMmcOzzz4bsukHgy471jQbe6HIx/PNaseOHU6X4GqajzUnsjnrrLP461//CsCtt97KmjVr2ryG5tJlx5pmYy8U+Xi+WSnlNddddx3XXXcdR48e5eKLL6a0tNTpkpRyPc83q0GDBjldgqtpPtaczOYvf/kLw4cP5/PPP+ecc87hww8/dKwWK7rsWNNs7IUiH883Ky8cAuwkzceak9nExsaybNkyJkyYwL59+zjvvPN47bXXHKunMbrsWNNs7IUin2Y1KxEZJyIbReRTEfl5I+9fKSJr/Y8PRSQ96JVaWL9+fVtNypM0H2tOZxMXF8ebb77JjTfeyNGjR5kyZQoLFy50tKb6nM7HzTQbe6HIp8lmJSKRwF+B8UAacIWIpDUYbAsw2hhzJvAA8PdgF6pUOIqKimL27Nn8/Oc/p6amhilTprBy5Uqny1LKdZqzZjUc+NQYs9kYcwx4CZhUfwBjzIfGmAP+px8BScEt01pCQkJbTcqTNB9rbslGRPjd737HTTfdRFVVFRMnTqSwsNDpslyTjxtpNvZCkU9zrmCRCGyv97wMyLIZ/npgSWNviMg0YBpA3759WbFiBQD9+/cnPj6eoqIiwHe1gMGDB7Nq1SpfkVFR5OTkUFhYSHl5OQCZmZns3r2bnTt3smvXLlJTU4mJiaG4uBjw3aly4MCBgfuqxMTEkJ2dTUFBQeBSIFlZWZSVlQUOsxw0aBCRkZGBVdiEhAT69etHXl4e4NvPkJWVRX5+PpWVlQBkZ2ezZcuWwAVP09LSqKmpYePGjb7wEhNJSkoiPz8f8G36yczMJC8vj6qqKgBycnIoLS1lz549AAwZMoSqqio2bdoEQHJyMr1796buJOrOnTuTkZFBbm5u4BYWo0aNoqSkJHB18PT0dA4fPhzIJyUlhe7duwf+CHbr1o309HRWrlyJMQYRYfTo0RQVFXHggO//HRkZGezfv5+tW7e2ej5t3+5bhNw0nzp16hRYBp2eT5s3b+bSSy9lx44dLFq0iPPOO48nn3ySoUOHOjafunTpEshHf5/+N58AkpKSKC8v198ni/nUu3dvKisrWzSfLBljbB/ApcAz9Z5PBZ60GPZcYAPQo6nxDhs2zATD8uXLgzKecKX5WHNjNlVVVeaCCy4wgElJSTE7duxwrBY35uMWmo291uQDFJhGekZzNgOWAcn1nicBXzQcSETOBJ4BJhlj2t/Nf5QKgujoaF555RWGDx/O1q1bGTduHAcPHnS6LKUc15xm9TGQKiL9RCQauBz42iFLInIK8Bow1RjTpmc46vW57Gk+1tyaTVxcHIsXL+b0009n3bp1TJw4MbD5pS25NR830GzshSKfZl3IVkQmAH8CIoG5xpgHReRmAGPMLBF5BvgB8Ln/I9WmkQsR1qc3X1TK3rZt2xg5ciRlZWVMnDiR1157jagovVGCCm+tupCtMeZtY8xAY8xpxpgH/a/NMsbM8v98gzGmmzHmLP+jze5KVrcDTzVO87Hm9mxOOeUUli5dSrdu3Xjrrbe46aabaM5/LoPF7fk4SbOxF4p8PH8FCyc2j3iJ5mPNC9mkpaWxePFiYmNjmTt3LjNmzGizaXshH6doNvZCkY/nm5VS4S47O5sFCxYQGRnJww8/zEMPPeR0SUq1Oc/ffLGqqoqYmJggVBSeNB9rXsvmxRdf5KqrrsIYwxNPPMHtt98e0ul5LZ+2pNnYa00+YXvzxS1btjhdgqtpPta8ls0Pf/hD/v5335XM7rjjjpDfvNFr+bQlzcZeKPLxfLOqO4NaNU7zsebFbG644Qb++Mc/AnDLLbfw+eefN/GJlvNiPm1Fs7EXinw836yUam9++tOfcsUVV1BVVdWmB1wo5STPN6u0tIYXgFf1aT7WvJzNQw89RExMDPPmzeOjjz4KyTS8nE+oaTb2QpGP55tVTU2N0yW4muZjzcvZnHrqqdx1110A3HXXXSE5/8rL+YSaZmMvFPl4vlnVXeVXNU7zseb1bH7+85/Tq1cv8vLyWLBgQdDH7/V8QkmzsReKfDzfrJRqrzp37swDDzwAwN13382XX37pcEVKhY7nm1ViYqLTJbia5mMtHLK57rrrGDFiBGVlZfzwhz8M6uaXcMgnVDQbe6HIx/PNKimpzW5K7Emaj7VwyCYqKooFCxZw8skn89577/Gb3/wmaOMOh3xCRbOxF4p8PN+s9IKS9jQfa+GSTVJSEvPnzyciIoIHH3yQhQsXNv2hZgiXfEJBs7GnF7JVSjXq3HPP5fe//z3gu9KF/jFV4cbzzSouLs7pElxN87EWbtlMnz6dqVOncuTIEcaPH09xcXGrxhdu+QSTZmMvFPl4/kK2Sqn/OX78OJMnT2bhwoX06dOH3Nxc+vfv73RZSjVb2F7INi8vz+kSXE3zsRaO2Zx00knMnz+fc889l507dzJ+/HgOHz7conGFYz7BotnYC0U+nm9WVVVVTpfgapqPtXDNpkOHDrz55psMHTqU0tJSbrvtthaNJ1zzCQbNxl4o8vF8s1JKfVN8fDzz58+nY8eO/Otf/+L55593uiSlWsXz+6yqq6uJiooKQkXhSfOx1h6ymTNnDjfccANxcXEUFhaSmpra7M+2h3xaSrOx15p8wnafVWlpqdMluJrmY609ZHPddddx+eWXU1FRQWZmJldffTVvvfUWx44da/Kz7SGfltJs7IUiH883qz179jhdgqtpPtbaQzYiwqxZsxg1ahTl5eU899xzfO9732PkyJFNHnjRHvJpKc3GXijy8XyzUkrZ69KlCytXrmTjxo08+OCDJCUlUVBQwJQpU6iurna6PKWaxfPNasiQIU6X4Gqaj7X2ls3AgQP5xS9+wfLly+nZsyfvvPMOt9xyi+W9sNpbPidCs7EXinw8v4dQDyG1p/lYa6/ZDBgwgIULFzJ27FieeeYZ9u7dy4gRI0hLS+Occ86hZ8+eQPvNpzk0G3t66HojNm3a5HQJrqb5WGvP2WRnZ/Piiy8SERHBm2++yYwZM5g0aRJ9+vRh0qRJvPLKK5SUlDhdpmu152WnOUKRj+fXrJRSLXPJJZdQXFzM6tWr2bBhA0VFRaxatYqFCxeycOFCevbsya9//WtuuOEGYmNjnS5XtXOeb1bJyclOl+Bqmo81zQbOOOMMzjjjjMDz3bt3M2/ePObMmUNxcTF33HEHDz30ENdccw3nnXceI0eO1MaFLjtNCUU+nj8puKKiQq+AbEPzsabZWKutrWX+/Pk8+uij/N///V/g9ejoaCZNmsT06dMZPnw4xhhKSkpYunQpSUlJjBs3ji5dujhYedvQZcdea/KxOinY82tWBQUFjBkzxukyXEvzsabZWIuIiKBPnz588sknLFu2jKVLl7Js2TLWrFnDggULWLBgAdnZ2ezbt+9rJ4BGRUUxatQo0tLS6N69O927d6d///4MGTKEU089lYgIz+8mB3TZaUoo8mlWsxKRccCfgUjgGWPM7xu8L/73JwBfAdcaYwqDWqlSqs2JCOeffz7nn38+ADt27ODJJ59k1qxZgStr9+zZkwsvvJCtW7eyevVq3n//fd5///1vjKtTp06cc845XHDBBYwbN460tDR8fzqUalqTzUpEIoG/At8ByoCPRWShMWZ9vcHGA6n+RxbwlP/fkOvcuXNbTMazNB9rmo29xvJJTEzk97//Pb/85S9ZuHAhCQkJjB49OnAduP3797Ns2TK++OILDhw4wN69eyktLaW4uJhdu3bx3nvv8d577zF9+nS6d+/O8OHDGTZsGNHR0Rw+fJgjR44QHR1Np06diIuLIz09nZycnG/UcvToURYvXsySJUs4cuQItbW11NTUUFlZyVdffUVtbS0XXngh119/PT169Ah8bufOnfz73//mvffeY+3atWRlZXHxxRczduxYqqqq2LZtG4cOHWLYsGF06NDhG9+/traWDz/8kHfffZcBAwaQlJT0tZo+//xzUlNTv7YGWV5eTkFBQeAE7OjoaE477TQSExO/saZZWlrKCy+8wPr16zn99NNJT0+nf//+VFVV8dVXX1FVVUVERAQRERF0796djIyMr43j2LFjbN68mb1797J3716MMQwZMoQBAwYQGRlpO79ramo4evQoHTt2REQ4fvw4xcXFFBQUEBsby+TJkxvN5MCBAxQVFbF161auvfZaIDS/W03usxKRbOA+Y8wF/uczAIwxD9UbZjawwhgzz/98IzDGGLPTarx680Wl2pfdu3ezbNky3nnnHd577z127drVrM9FRERw1llnkZiYSHx8PDU1NSxZsoTy8vImP9uhQwcuueQSDh06RFFRETt27Gh0uMjISGpqagLP4+Pj+d73vsfEiROJiIjg0KFDlJaWMn/+fLZt2wb47h125ZVXcumll7Jo0SLmzZvHwYMHOeWUU7jqqqvIyspiwYIFvPrqq1RWVn5jmrGxsaSkpNCtWze6dOnCnj17+OSTT5qVSZ2kpCQuu+wyzjjjDJYuXcrSpUsbzSU2NpbBgwczdOhQhg4dSlJSUqCxb926lfz8fAoKCqioqCAyMpIuXbpw5MiRr50v1bt3b+68806+/e1v89FHH5Gbm0tBQQHbt28HfGvhhw8fplOnTif0HRqy2mfVnGY1GRhnjLnB/3wqkGWMua3eMIuA3xtjcv3PlwH3GmMKGoxrGjANoG/fvsNeeOEFAPr37098fDxFRUUA9OjRg8GDB7Nq1SrAtx08JyeHwsLCwIzIzMxk9+7dbNiwgbi4OFJTU4mJiQncyrtXr14MHDiQ3NxcAGJiYsjOzg7MEICsrCzKysoCC/CgQYOIjIxk/XrfSmNCQgL9+vULbO6IjY0lKyuL/Pz8wMKXnZ3Nli1bAr94aWlp1NTUsHHjRsD3P9GkpCTy8/MB3+2eMzMzycvLCywIOTk5lJaWBq6nNWTIEKqqqgLnKiQnJ9O7d2/qmnvnzp3JyMggNzc38L+1UaNGUVJSwr59+wBIT0/n8OHDrF27lri4OFJSUujevTuFhb6ts926dSM9PZ2VK1dijEFEGD16NEVFRRw4cACAjIwM9u/fz9atW1s9n+oWaDfNp+Li4sD/Sp2eT5s3bwZw1XwqLS0NbKYL9nwyxlBeXs7x48dZvHgxNTU1xMbGcsYZZ7Br1y6+/PJLysvL2bJlC2vWrGn0slCpqamMHz+erKwsPvvsM2pqaoiJiSEzM5P169fz4osv8p///Odrn+nQoQPDhw9n/PjxxMTEUFRURF5eHqWlpcTExNCrVy+ioqLYsmXLN6ZX5+STTyYlJYVPPvmE2trar70XFxcXyKO+008/nYSEBKqrqzlw4EBgzbOh+Ph4Ro4cyZlnnskXX3zB/v37KS0t5aSTTqJDhw707NmTo0ePUllZybZt2xq9Bl9KSgqdOnWia9euGGMoKysLNNimREdHf+0ix0lJSZxxxhl8+eWXXzvQpr6YmBjOOOMM0tLSmDx5Mt26deP48ePk5OS06PcpKSmpxc3qUuCCBs1quDHm9nrDLAYeatCsfmaMsfxvQrDWrFasWKE7Om1oPtY0G3tuyefIkSMUFhayf/9+Dh8+zNGjR8nJyeH0009v8rOlpaW8/fbbJCYmkp6ezmmnndbo5rCvvvqK2NjYQHP+7LPPePnll1m9ejUdO3aka9eugX1zI0eOZNWqVSQnJ/PYY4/xwQcfcP755/OjH/2IwYMHs3r1av71r39RXFzM+PHjufrqq+nfv/83pnnw4MHAZsdDhw4RGRnJmDFjmn1qQG1tLR999BEvvfQSZWVljB07lgsvvJB+/fp9Y9gDBw5QXFzMunXrWLduHXv27KFjx4507NiRXr168a1vfYvhw4eTkJDA8ePHOXToENHR0YHNecYY3n33XR5//HF27NhBdnY2OTk5jBgxgtTU1G9k2pplx2rNCmOM7QPIBpbWez4DmNFgmNnAFfWebwT62I132LBhJhiWL18elPGEK83HmmZjT/OxptnYa00+QIFppGc0Z80qCigFzgN2AB8DPzTGlNQb5kLgNnxHA2YBTxhjhtuNN1hrVrW1tWFzOGwoaD7WNBt7mo81zcZea/Jp8c0XjTHV+BrRUmAD8LIxpkREbhaRm/2DvQ1sBj4FngZubVGVLaDXL7On+VjTbOxpPtY0G3uhyKdZ51kZY97G15Dqvzar3s8G+HFwS2ueuh3VqnGajzXNxp7mY02zsReKfHQ9VimllOt5vlmlp6c7XYKraT7WNBt7mo81zcZeKPLxfLM6fPiw0yW4muZjTbOxp/lY02zshSIfzzerupMpVeM0H2uajT3Nx5pmYy8U+Xi+WSmllAp/jt3PSkT2Ap8HYVQ9gS+DMJ5wpflY02zsaT7WNBt7rcnnVGPMyQ1fdKxZBYuIFDR2Apny0XysaTb2NB9rmo29UOSjmwGVUkq5njYrpZRSrhcOzervThfgcpqPNc3GnuZjTbOxF/R8PL/PSimlVPgLhzUrpZRSYU6blVJKKdfzTLMSkXEislFEPhWRnzfyvojIE/7314pIhhN1OqEZ2ZwuInkiUiUi052o0UnNyOdK/zKzVkQ+FJF2c+G3ZmQzyZ/LGhEpEJEcJ+p0SlP51BvuWyJSIyKT27I+pzVj+RkjIof8y88aEfl1iyfW2B0Z3fYAIoHPgP5ANFAEpDUYZgKwBBBgBJDvdN0uyqYX8C3gQWC60zW7MJ9zgG7+n8frsvO1YeL4377tM4H/Ol23m/KpN9z7+G6jNNnput2UDzAGWBSM6XllzWo48KkxZrMx5hjwEjCpwTCTgH8Zn4+AriLSp60LdUCT2Rhj9hhjPgaOO1Ggw5qTz4fGmAP+px8BSW1co1Oak02F8f/VAToB7emIrOb83QG4HXgV2NOWxblAc/MJCq80q0Rge73nZf7XTnSYcNRev3dznWg+1+NbQ28PmpWNiFwiIv8FFgPXtVFtbtBkPiKSCFwCzKL9ae7vVraIFInIEhEZ3NKJeaVZSSOvNfwfXnOGCUft9Xs3V7PzEZFz8TWre0NakXs0KxtjzOvGmNOBi4EHQl2UizQnnz8B9xpjakJfjus0J59CfNf6SweeBN5o6cS80qzKgOR6z5OAL1owTDhqr9+7uZqVj4icCTwDTDLGtJd7lp/QsmOMWQWcJiI9Q12YSzQnn0zgJRHZCkwG/iYiF7dJdc5rMh9jTLkxpsL/89vASS1dfrzSrD4GUkWkn4hEA5cDCxsMsxC42n9U4AjgkDFmZ1sX6oDmZNOeNZmPiJwCvAZMNcaUOlCjU5qTzQAREf/PGfh2pLeXZt5kPsaYfsaYFGNMCvAKcKsx5o02r9QZzVl+EuotP8Px9ZwWLT9RrSy2TRhjqkXkNmApviNQ5hpjSkTkZv/7s/AdiTMB+BT4CviRU/W2peZkIyIJQAHQGagVkZ/iO2qn3Km620ozl51fAz3w/a8YoNq0gytqNzObH+D7T+BxoBK4rN4BF2Gtmfm0W83MZzJwi4hU41t+Lm/p8qOXW1JKKeV6XtkMqJRSqh3TZqWUUsr1tFkppZRyPW1WSimlXE+blVJKKdfTZqWUUsr1tFkppZRyPW1WSilLIvKkiBSKyLeaMWx/EZkjIq+0RW2qfdFmpZRqlIh0wncvtJuAi5oa3n+riOtDXphql7RZqbAnIl1F5NYGr33YjM9VWLx+h4hsEJEXglWjG4hIrIisFJFIAGPMEaAPsAJ4ot5wQ0VkUYNHL4txRovIKhHxxKXdlHvpAqTag67ArcDf6l4wxpzTivHdCow3xmype8F/sU4xxtS2YrxOuw54re52FyLSA+gIHAYCt8AwxqyjGWta/mGPicgy4DIgrJq7alu6ZqVcTUQ6ichi/83bikXkMhFJEZH/isg/RWStiLwiIh39w78hIp+ISImITPOP5vf4bm2xRkQe9Q9XUW8ajX3Gqp5Z+G7jvVBE7vSvYf0N3317kv3DXCUi//FPb3bdmoqI/FJENorIv0VknohM93+X4nrjny4i9zUxnhT/dJ/21/yuiMT637van0mRiDwnIg+IyE/qjf9BEbnD4utdCbxZ7/lM4DGgBEizy8U/7h7+fM4WkRn13nrDP26lWq6p+97rQx9OPvBd9fvpes+7ACn4bvI20v/aXGC6/+fu/n9jgWJ8V1NPAYobjLei3s/f+EzDYRp8divQ0z/eWmBEvffOAN4CTvI//xtwNTAMWIdvTaUzvrsDTG9Ym/+1+6zG4/85BagGzvI/fxm4ChgMbAR61n0v/7CF/ucRwGd136/Bd4oGdtV7ngLk47vB3l+AG1sxDyOBvU4vS/rw9kM3Ayq3Wwc8JiIPA4uMMatFpBuw3RjzgX+Y54E78K0F3CEil/hfTwZSgV1NTKOxzzT3njufG2M+qvf8PHyN6WP/7UZigT34GsfrxpivAESkqXuOWY2nzhZjzBr/z5/gay7dgFeMMV8CGGP2A/tFZJ+InA30Bv7PNH5zyZ7AwXrPfwv8P2OMEZEN+BphixhjakTkmIjEG2MOt3Q8qn3TZqVczRhTKiLD8N2r7CEReRf4F9+8fbYRkTHA+UC2MeYrEVkBdLAbf0s+08CRhqME/mmMmfG1F333EGvsfjzVfH1zfN20Gx1PPVX1fq7B18zEYhrPANcCCfjWQhtTWTdtETkL+D6QIyJ/9b++zuJzzRUDHG3lOFQ7pvuslKuJSF/gK2PM8/jWnDL8b50iItn+n68AcvFtIjzgbzqnAyP87x8G4i0mYfWZlloGTK47Ok5EuovIqcAq4BL/EXfxwET/8LuBXv79PTH878AFq/E0Ne0p/gMjEJHu/tdfB8YB38J3o7xvMMYcACJFpAPwMDDR/O8OuOm0Ys3KX89eY8zxlo5DKV2zUm43FHhURGqB48At/tc3ANeIyGxgE/AUvjWMm0VkLb59Nx8BGGP2icgH/gMZlhhj7qk3/nca+0xLGWPWi8hM4F0RifDX/GNjzEciMh9YA3wOrPYPf1xE/h++/UNbgP/ajcf/Watpl4jIg8BKEakB/g+41viOyFsOHDT+I/0svAuMAjoZY5bVG+9u/4Eu3f2bFk/Uufju5K1Ui+mdgpXniEgKvv1XQ5yupaX8R/xVGGMea4NpReA7WvFSY8wmm+HOBu4yxkwN8vRfA2YYYzYGc7yqfdHNgEqFMRFJw3fk4TK7RgVgjPk/YHndIfJBmn408IY2KtVaumallFLK9XTNSimllOtps1JKKeV62qyUUkq5njYrpZRSrqfNSimllOtps1JKKeV62qyUUkq53v8HObm7USmrBRAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 504x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(7,4))\n",
    "# plt.plot(freq, c1, color=cmap_a[0], linewidth=2, label='split row')\n",
    "# plt.plot(freq, c3, color=cmap_a[2], linewidth=2, label='split layer')\n",
    "# plt.plot(freq, c2, color=cmap_a[3], linewidth=2, label='split column')\n",
    "plt.plot(freq, c_avg, color='k', linestyle='-', linewidth=2, label='average')\n",
    "plt.rc('axes', labelsize=12)\n",
    "plt.rc('xtick', labelsize=10)\n",
    "plt.rc('ytick', labelsize=10)\n",
    "plt.rc('legend', fontsize=11)\n",
    "plt.grid(which='major', linestyle='--')\n",
    "plt.xlabel('spatial frequency 'r'(${\\AA}^{-1}$)')\n",
    "plt.title('Fourier shell correlation')\n",
    "# plt.legend()\n",
    "# plt.savefig('')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eba5dc90-d121-4fc8-b96e-af5f14a49df9",
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}