{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3343f22a",
   "metadata": {},
   "source": [
    "# Another Power Spectrum Example"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d207e54d53df97c",
   "metadata": {},
   "source": [
    "This notebook shows an example of how to generate a Gaussian random vector field with a power spectrum model, and compute its 1D and 3D power spectra. It also shows how to compare the 1D power spectrum computed from the 3D power spectrum with the one computed directly from the power spectrum model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "284bcc5d3c786a31",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2026-08-11T14:28:39.467661Z",
     "start_time": "2026-08-11T14:28:39.108204Z"
    },
    "execution": {
     "iopub.execute_input": "2026-08-11T18:23:31.164668Z",
     "iopub.status.busy": "2026-08-11T18:23:31.164492Z",
     "iopub.status.idle": "2026-08-11T18:23:31.604519Z",
     "shell.execute_reply": "2026-08-11T18:23:31.604068Z"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.integrate import quad_vec\n",
    "from kspace import FourierAnalysis, GaussianRandomField, PowerLawBetaModel"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93358dcbd3df9b20",
   "metadata": {},
   "source": [
    "We set up our power spectrum model using the `PowerLawBetaModel` class, which has the following functional form:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "eb5dab45e0343cb7",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2026-08-11T14:28:39.480426Z",
     "start_time": "2026-08-11T14:28:39.472892Z"
    },
    "execution": {
     "iopub.execute_input": "2026-08-11T18:23:31.605691Z",
     "iopub.status.busy": "2026-08-11T18:23:31.605595Z",
     "iopub.status.idle": "2026-08-11T18:23:31.607103Z",
     "shell.execute_reply": "2026-08-11T18:23:31.606800Z"
    }
   },
   "outputs": [],
   "source": [
    "l_min = 10.0  # minimum or \"dissipation\" scale\n",
    "l_max = 200.0  # maximum or \"injection\" scale\n",
    "alpha = -11.0 / 3.0\n",
    "f_rms = 10.0  # normalization of the field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "80371b74a68cfa29",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2026-08-11T14:28:39.488552Z",
     "start_time": "2026-08-11T14:28:39.481989Z"
    },
    "execution": {
     "iopub.execute_input": "2026-08-11T18:23:31.607941Z",
     "iopub.status.busy": "2026-08-11T18:23:31.607847Z",
     "iopub.status.idle": "2026-08-11T18:23:31.609360Z",
     "shell.execute_reply": "2026-08-11T18:23:31.609031Z"
    }
   },
   "outputs": [],
   "source": [
    "power_spec = PowerLawBetaModel(l_min, l_max, alpha)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e960f7d1613bea15",
   "metadata": {},
   "source": [
    "and we normalize it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "365c28b6",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2026-08-11T14:28:39.496036Z",
     "start_time": "2026-08-11T14:28:39.490013Z"
    },
    "execution": {
     "iopub.execute_input": "2026-08-11T18:23:31.610283Z",
     "iopub.status.busy": "2026-08-11T18:23:31.610211Z",
     "iopub.status.idle": "2026-08-11T18:23:31.611730Z",
     "shell.execute_reply": "2026-08-11T18:23:31.611448Z"
    }
   },
   "outputs": [],
   "source": [
    "# Renormalize the power spectra to have the desired RMS value\n",
    "power_spec.renormalize(f_rms)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb964a4e8700632b",
   "metadata": {},
   "source": [
    "Next, we set up the Gaussian random field generator, and generate a scalar field realization:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0b20fe06",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2026-08-11T14:28:39.504634Z",
     "start_time": "2026-08-11T14:28:39.498374Z"
    },
    "execution": {
     "iopub.execute_input": "2026-08-11T18:23:31.612637Z",
     "iopub.status.busy": "2026-08-11T18:23:31.612582Z",
     "iopub.status.idle": "2026-08-11T18:23:31.614070Z",
     "shell.execute_reply": "2026-08-11T18:23:31.613741Z"
    }
   },
   "outputs": [],
   "source": [
    "# Parameters for the box and grid\n",
    "le = np.array([0.0, 0.0, 0.0])\n",
    "re = np.array([750.0, 750.0, 750.0])\n",
    "ddims = [256] * 3\n",
    "width = re - le"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "412cb60050fb751a",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2026-08-11T14:28:40.503439Z",
     "start_time": "2026-08-11T14:28:39.505445Z"
    },
    "execution": {
     "iopub.execute_input": "2026-08-11T18:23:31.614872Z",
     "iopub.status.busy": "2026-08-11T18:23:31.614827Z",
     "iopub.status.idle": "2026-08-11T18:23:32.534942Z",
     "shell.execute_reply": "2026-08-11T18:23:32.534480Z"
    }
   },
   "outputs": [],
   "source": [
    "# This makes a Gaussian random field with the specified power spectrum\n",
    "g = GaussianRandomField(le, re, ddims, power_spec, seed=20)\n",
    "s = g.generate_scalar_field_realization()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16e368eb6cd76b58",
   "metadata": {},
   "source": [
    "Now we can show how to take the power spectrum of the field in 3D, derive its 1D power spectrum along the x-direction, and compare it to theactual power spectra of both fields and compare them to their input power spectra. First, we create an instance of the `FourierAnalysis` class, which will help us with these tasks.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "6e25c55263590dad",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2026-08-11T14:28:40.522021Z",
     "start_time": "2026-08-11T14:28:40.515910Z"
    },
    "execution": {
     "iopub.execute_input": "2026-08-11T18:23:32.536138Z",
     "iopub.status.busy": "2026-08-11T18:23:32.536073Z",
     "iopub.status.idle": "2026-08-11T18:23:32.537698Z",
     "shell.execute_reply": "2026-08-11T18:23:32.537380Z"
    }
   },
   "outputs": [],
   "source": [
    "# Give the FourierAnalysis class the same width and dims as\n",
    "# the GaussianRandomField created above\n",
    "fa = FourierAnalysis(width, ddims)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4da24271569275a",
   "metadata": {},
   "source": [
    "We can take the scalar field `s` and compute its 3D power spectrum $P_{\\rm 3D}({\\bf k})$ using the `make_powerspec` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f95d08bfdc03d7fb",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2026-08-11T14:28:40.759596Z",
     "start_time": "2026-08-11T14:28:40.525095Z"
    },
    "execution": {
     "iopub.execute_input": "2026-08-11T18:23:32.538599Z",
     "iopub.status.busy": "2026-08-11T18:23:32.538532Z",
     "iopub.status.idle": "2026-08-11T18:23:32.739195Z",
     "shell.execute_reply": "2026-08-11T18:23:32.738819Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(256, 256, 256)\n"
     ]
    }
   ],
   "source": [
    "P_3D = fa.make_powerspec(s)\n",
    "print(P_3D.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c05f10a55bd9841c",
   "metadata": {},
   "source": [
    "To get the 1-dimensional power spectrum for just one direction, integrated over the other two directions, we can use the `integrate_kspace` method of the `FourierAnalysis` class, and choose the axes we want to integrate over. If we choose `axis=(1,2)`, this will give us the 1D power spectrum along the x-direction. This is equivalent to performing this integral:\n",
    "\n",
    "$$P_{\\rm 1D}(k_x) = \\frac{1}{(2\\pi)^2}\\displaystyle\\int{P_{\\rm 3D}({\\bf k})}dk_ydk_z$$\n",
    "\n",
    "we also use the `average_symmetric_k` method of the `FFTArray` class to average the $-k$ and $+k$ components of the power spectrum, since it is symmetric:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d331f2aa784d5d94",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2026-08-11T14:28:40.775032Z",
     "start_time": "2026-08-11T14:28:40.761295Z"
    },
    "execution": {
     "iopub.execute_input": "2026-08-11T18:23:32.740280Z",
     "iopub.status.busy": "2026-08-11T18:23:32.740214Z",
     "iopub.status.idle": "2026-08-11T18:23:32.744284Z",
     "shell.execute_reply": "2026-08-11T18:23:32.743940Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(128,)\n"
     ]
    }
   ],
   "source": [
    "P_1Ds = fa.integrate_kspace(P_3D, axis=(1, 2)).average_symmetric_k()\n",
    "print(P_1Ds.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae4371684f7a44df",
   "metadata": {},
   "source": [
    "It will be illustrative to show how the fully numerical integration compares to a simple semi-analytic integration, which we can do using the `quad_vec` integration method from SciPy. Since the scalar field is isotropic, $P_{\\rm 3D}(k) = P_{\\rm 3D}({\\bf k})$, and $k = \\sqrt{k_x^2+k_\\perp^2}$, where $k_\\perp^2 = k_y^2+k_z^2$. Then,\n",
    "\n",
    "$$P_{\\rm 1D}(k_x) = \\frac{1}{(2\\pi)^2}\\displaystyle\\int{P_{\\rm 3D}(k)}dk_ydk_z = \\frac{1}{(2\\pi)^2}\\displaystyle\\int{P_{\\rm 3D}(k)}2\\pi k_\\perp dk_\\perp$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c2a6ec685fee6ea",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2026-08-11T14:28:41.191427Z",
     "start_time": "2026-08-11T14:28:40.780480Z"
    },
    "execution": {
     "iopub.execute_input": "2026-08-11T18:23:32.745231Z",
     "iopub.status.busy": "2026-08-11T18:23:32.745163Z",
     "iopub.status.idle": "2026-08-11T18:23:32.855802Z",
     "shell.execute_reply": "2026-08-11T18:23:32.855325Z"
    }
   },
   "outputs": [],
   "source": [
    "kmin = 0.0\n",
    "kmax = 1000.0  # a really big number to approximate infinity\n",
    "\n",
    "def P1_int(k_perp, k_x):\n",
    "    k = np.sqrt(k_x ** 2 + k_perp ** 2)\n",
    "    return 2.0*np.pi*k_perp*power_spec(k)\n",
    "\n",
    "k_x = fa.kx[:128, 0, 0] # Remember that we only chose the positive k_x\n",
    "\n",
    "P_1D = quad_vec(P1_int, kmin, kmax, args=(k_x,))[0]/(2.0*np.pi)**2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "415e96dc27d496ea",
   "metadata": {},
   "source": [
    "Now we can plot them against each other and see the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "184c12ac6eaeecca",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2026-08-11T14:28:41.405242Z",
     "start_time": "2026-08-11T14:28:41.200019Z"
    },
    "execution": {
     "iopub.execute_input": "2026-08-11T18:23:32.857032Z",
     "iopub.status.busy": "2026-08-11T18:23:32.856965Z",
     "iopub.status.idle": "2026-08-11T18:23:33.008289Z",
     "shell.execute_reply": "2026-08-11T18:23:33.007862Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$P_{\\\\rm 1D}(k_x)$')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj0AAAG4CAYAAACuBFb3AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAASrlJREFUeJzt3Qd4lFX69/E7hYReQyeACIJBIAoBURCQKOguioq9IOuLrkZFsKzs+hd0VVRWQAFFQYprwwYuKqhA6CAd0aBIEZBeBBJK+nvdJ5nJzGQSJmXq8/1c11zJPPNk8mQgzI/73OecsNzc3FwBAAAIceH+vgAAAABfIPQAAABLIPQAAABLIPQAAABLIPQAAABLIPQAAABLIPQAAABLiPT3BQSSnJwc2bdvn1SrVk3CwsL8fTkAAMADuuRgamqqNGrUSMLDi67nEHocaOCJjY315PUFAAABZs+ePdKkSZMiHyf0ONAKj+1Fq169uvf/dAAAQJmdPHnSFC1s7+NFIfQ4sA1paeAh9AAAEFzO1ZpCI7OITJw4UeLi4iQhIcFXfy4AAMDHwthw1Lk8VqNGDTlx4gSVHgAAQuz9m0oPAACwBHp6AAA+lZ2dLZmZmbzq8FiFChUkIiJCyorQk9/Tozf9RQQAeG8tlQMHDsjx48d5iVFiNWvWlAYNGpRpHT16ehzQ0wMA3rN//34TeOrVqyeVK1dmEVh4HJZPnz4thw4dMsGnYcOGpX7/ptIDAPA6raTbAk+dOnV4xVEilSpVMh81+OjfodIOddHIDADwOlsPj1Z4gNKw/d0pSz8YoYd1egDAZ9jXEP78u0PoEZGkpCRJSUmRNWvWlPkFBQAAgYnQAwAALIHQAwBAgA7nzJ49WwLBvffeK/379w/6n4nQAwBAMQ4fPiwPPvigNG3aVKKjo81aMX369JHly5d7fYr/Nddcc87zzpw5I7Vr15aYmBhJT0+XQDBy5EiJj48v9c/kLUxZd7M4YccXvpeI6PKdYRAeFibVK1aQWlWipHaVClKrsn6MMh9rVbYdz7uvH/VWsULZV58EgEDz9Oc/yq8HU/16Da3rV5OXb2rv0bk33XSTZGRkyIwZM6RFixZy8OBBWbBggRw9etSr16jhyhOff/65tG3b1qxno1WUW2+9VQJVAw9/Jm8h9OQ3MuvNtrhRemaOhIfnlPuLfTojWw6cPOvx+ZUqROQFIZeQlHcsSmq7BKaalStIdCRBCUBg08CzYXdwrMqsawstXbpUFi1aJD169DDHmjVrJp07d3Y654knnpAvv/zSVFo6deokY8eOlQ4dOtirHhpGHn30UfP5sWPH5J577pHx48fLa6+9JmPGjJGcnBwZMmSI/Otf/3IaCpo1a9Y5h5Xeffddueuuu0zo0c9dQ48+z+TJk+Xrr7+Wb7/9Vho3bmy+73XXXWce1//w33///bJw4UKzYrZWtB566CFzPe689957MnToUNm3b5+pfNnodVarVk169+4tzz33nP17q2nTppkhMtef6Y8//pAnn3zSXJe+dhdeeKEpQnTp0kW8gdATwM5kZsve42fMzVNVoyNNSDKByBaMnKpIeQGqlkOVKTKCUU4AcPtvatWq5qah5dJLL3V6k7e5+eabzeJ5c+fONf9xfvvtt80b/9atW82wk9q+fbt5fN68eebzAQMGyI4dO+SCCy6QxYsXy4oVK+Rvf/ubJCYmlugNX59r5cqV8sUXX5jQo2Fk165dJpg50hDy6quvyujRo03YuvPOO815en0auJo0aSKffvqpWThSr0VDkK58fMstt7j9eTXA/e9//zOf2xYN1FD13Xffmdfpp59+Mj/r/PnzzeP6urhKS0szQVJDmD6XVoHWr19vrsdbCD1e1insF1mb28bj42WVlp5lbnuOeR6UqleMdK4euQlJtsf18xqVKkhEuIfrJexaKdKsq+fHASCAREZGyvTp02Xw4MEyadIkueSSS8wb9W233Sbt27eXZcuWyerVq82bvi0Q/ec//zEh6bPPPjPhQekb+dSpU00lJC4uTnr16iW//vqrfPPNNxIeHi6tW7eWV155RZKTk0sUevQ5tUemVq1a5r72GmlVRStKjrTKcvvtt5vPX3rpJXnjjTfMdfft29ds5mmrzKjzzjvPBKlPPvnEbejRgHfHHXeY72MLPe+//76pEPXs2dNUczQo6mtX3HDWhx9+aPqldLkYWzhs2bKleBOhx4sei/xMHov8Ql7OvE0mZeeVEdXfI/4nT1f4WMZl3SjjsgaIv508m2Vuvx897dH5Wq2sWSl/WK1QRakgJMX9OlEabnxdzvR4Vir2HFawsNSysSLzR4r0eFqk13Dv/nAAUEba0/OXv/zFDHOtWrXKVGy0ajJlyhQ5deqUqVi4bq2hzcVahbFp3ry5CTw29evXN1spaOBxPKbhyR0NNvr9lVZxfv75ZzMspX1Gr7/+uv08HebSobZnn33W6bk1oNlUqVLF7E/l+L10SEkD1O7du821aw+Tu0ZkGw2BCQkJsnfvXlOp0WBoG77y1MaNG+Xiiy+2Bx5fIPS4MaBjE4muXLVML2yzU5vk79u/MJ9rwGnfpIZ8WfUWueLgB3Jn2sfmuAaiZdkXeaXi4025uSJ/ns40tx1yyu05Wsn6LDrvF7HS4ufllQW/yewqN8tjlb6WW4+/m3fS4pcls9kVUqHF5b68fAB+pk3EwXYNFStWlKuuusrc/u///k/+3//7fzJixAjT+6LDQNrz40o3x7TRaoojDQfujhU1tKMBS8OI43NpH4yGDtceHg1D2mit11rc98/J/14ff/yxCUra59O1a1cTznQY7Icffijy9dCwoj1L2t9z9dVXmxCmw1ul2U/Llwg9bmZvjbyubbG7tHqmnciyw3kVDRG59uDbcu3JT0TO/Gk/I6f3SJnSMUmOncqQP09nyLFTmfLnqQw5djoj76P9eEZ+yMiQ46dLv+eIL2mQ0wqXBj71j8iP5P6zc6RWepr9HH383XdPyPl1l0hcw+rSpmE1adMg72PdqtEsVw+EKE9nTQUyHaLSISwd7tLmXx3K0WqOt2g1xZU2Leswm2Pzs3rxxRfNY46hpzjLly+Xyy67zAQ4G8cqVVE0+I0bN84EL+1Fio2NtT8WFRVlf08tilafNMxpY7evqj2EHjezt8pNt6F5H/ODj2PgkcSREt5tqOj/A2pWjvL4KbOyc+TEmUx7SHIKRo6BSUNS/rHU9CzxB9uQni341ApzDjx5j+fKLwdSzU02FHxtnSpRBSGoQTW5sGF1aVmvKtP4AfiUTkvXvhVtMtY3aa2CrF271gxvXX/99ebNXqsjOhtJj2ljss5q0qrHDTfcYGZyeYP2wsyZM8c0AF900UVOj+nMMP3enoaJVq1amYqNVo60n+e///2v6bPRz4ujfT1aIdKZYfr1jjQA7ty50wxhaZO0vm6uTeDaY6T9RfrajRo1ylTMNmzYII0aNTKvqTcQerxNg8/y150DT6VaBYGohHSmVZ2q0ebmqYysHDmuwSg/HGm1qKiQZAtROr2+PGiweSDyK6fA82duVaceJ3eOnsqQ5duOmpuNNk+fF1PFHoL0Y5uG1aVRjYpUhQB4hTbkamOxTkHX6ofu8K0VDe1p+ec//2n+7dFmZK22DBo0yIQRbd694oorTI+Ot2jI0N4cnSXmSo/p0JE2F+ssq3N54IEHTNjQYTL9eTSMaNVHe5eKo0UC7XfSgOc6rV6P64wybdjWKf22KeuOtBqks70ef/xxufbaayUrK8tU0HTkxVvCcnWOGwxbpefEiRPlMLwlzk27rhJHljr4+MLZzGyHClJmEUNueY/pRw0pGq5c2Zq2Xbk2d5dFtYqRcmED5+ExHa+vEk2mBwLF2bNnzf/8tXqg/TEIDb179zYLI+psMH/+HfL0/Zt3BW9yDTxa4bFVfGzHAzT46GrQDWtUMjdPaHbWdYUcQ1LMxjelbUpB4EkNqybVcvNWYbUFofIIPqlns2T178fMzUYnEDSrXdkegvTjhQ2rSWytyhLu6XR7AIBbf/75p2ne1tubb74pwYLQ4y26Do1j4LFVdhyDkH6MvTQk1qvRkmjlqEhza1Ir/+dPGVNwQuJIqdZtqOQuHSthC0bag0/FFpfJ1yeay44jpyQ7p/yKjlq/1Cn4epv38wH78cpREdK6QUEIsn2sVtF5ZgMAQIqdvaXBR9cW0jWGggWhx1s0yOg6NItfdh7Kcmxu1sdDIPCU5OcP6z5UJKzg53+s10B5LH84bduhtLym5v0nzcct+0+aYbPypL1Kuvy94xL42isUH1tTureKkSsuqCsdmtT0fPFFALCg33//XYIRPT3e7umx+orEZfz5D6emyy8HTsov+1NlS/5HDUcZ2TleWyn7t4rtpFtLDUB5IcjTIT4ARaOnB2VFT4+X1ukpV0W9sVsh8JTDz1+3WrTUrVZXureqaz+WmZ0jO4+cMpUgx8rQ/hNny2el7MwbZdzmAfL15v3meKt6Vc331xB0aYs6TJsHgCBFpcfblR74jE7L37Jf1/zJqwjpR93N+WxmjstK0c/b79uCj+ssswHpz7qtBEVFhkuX82rLFSYE1ZUL6ldlujzgASo9CIRKD6GnFC8agoc2R+86espeEdpyIFUu2TNdHsz8r9O6Qe4XTjy3+tWj86tAdaV7yxiz/1ghVh/iBAg9KAcMbwHnoA3JLepWNbdr2zXMP9pJ0hedJ9GL8io+pQ086uDJdPls3R/mptPk2zeuYQKQ3rQ5usKSVwo3sys2XQUAnyvYghWwkOiej+etm+QgM6qmZHUdUurNEHWa/KY/Tsj4hdvk5kkr5b7nx+cFHsmfraZBRzkuW6CPa8UHAMqZbgWhe2OVl549e8pjj+l82+BF6IE1afBw3BpEdyHOOC7P1PxWvh16hawa3lteHdBe/tq+odSsXLo1fJaktzSVI7v5I+X0v2MLr9/EEBfgmaL+g+Dl/zjo9gm6FtnLL+f/JyafbjiqxwOV7p91//33+/syAgqhB9bjbqVsm/yKTIMaFeWWTrEy4Y5LZN0zV8nspMtl2FUXSKdmtUq0ho8OlTkGn8rZJ+2fL232sOxswz9IgEeSR4lM61tQMXX8fdbj+rgXaeOsLsSnC/IFuoyMvPXN6tatK5UrV/b35QQUQg+sxd1K2f/4Pe+jjT7u8D9H2+KFj/ZuJZ89eJlsePYqmXRXR7mjS1NpUquSR8FHm6Ud6f27f71Mev1nkVw/YZlMXbZTDqV6NuUesBz9ffTzULHupq4biepu4O6MHDlS4uPjnY7p0JIOMTlWjHRjTt1ZXDcjrVmzpjz//PNmo80nn3zS7IiuO5Lr5pyO9uzZI7fccos5X8/R3d0dFwe0Pe+LL75odii3rZDsOrx1/Phxs7mofm8Ncbo7+1dffWXfTV43Gm3cuLEJSu3atZOPPvpIQg0rMsNaymGl7OoVK0jfixqYm+45pmsGLdl6WJb8dkRWbj9q9iBzpNPhHZulld7X4xqItA9Iby98nSKXt4yR6+MbS5+29dkaA3D8vdXfV8ctfJa/7jxE7eWh4oiICBNW7rjjDrNzuYaT0li4cKH52iVLlsjy5cvlvvvukxUrVphd2X/44QeZOXOmCSZXXXWVOU93de/Tp4907dpVli5dKpGRkfLCCy9I37595ccffzQ7lasFCxaYWcfff/+92++bk5Mj11xzjaSmpprd188//3xJSUkxP5dtZlTHjh3lH//4h3ke3Tn97rvvNud17txZQgWhx9uLEyLw9Bou0qJn4X8gNfiUcC80Hc+3zQ679/LzJD0rW9b9/qcs/u2wLN16RK449L7T+j+O0+NdN13VrceW/nbE3P41K1wS4+pL//jG0uOCumZ9IMDSHP9jolwDjw82b77hhhtMNWfEiBHy7rvvluo5tFKjO5KHh4ebisyrr74qp0+fln/+85/m8eHDh5veoWXLlsltt91mQpAGlilTptj7h7QSpFUf3ezz6quvNseqVKlizrGFIFfz58+X1atXy5YtW+SCCy4wx1q0aGF/XCs8TzzxhP3+I488It9++6188sknIRV6+JdURJKSkkzi1aYvWISXVsqOjoyQy1rGyPBrLpRv+kc6BZ4xuXfIxenvOPX46OO6YKKr9Kwc+frH/TL4vbXS+aX58s9Zm2X1zmOSo8nIT82cgN9psHGZdWnu+yDw2Ghfz4wZM0x4KI22bduawGOjQ006lGSjlZc6derIoUOHzP1NmzbJtm3bpFq1alK1alVz0+CklZnt27fbv06fo6jAozZu3GgqR7bA40r/0//vf//bPI8+v34fDT27d++WUELoAbw9lKYSR8pDz0yQCXdcLNsuGCyjs283h8dl3eh25WdHx09nyoc/7JZb3l4p0158wDRtHp77sl+aOYFAm3Vp7rs2N3uRDkPpcJNWZBxpkNHhbkc6NOWqQgXn2aBavXF3TKs7Ki0tzQw7aWhxvG3dutUMtdlopac4lSoV3384evRoef31183wVnJysvke+nPamqJDBcNbgI+G0nTR9L+2b2Rux0+3l+8W95EVuxuI7Dzm0VNpRei+7Jnm87o/jJJ3N+2TzK5D5M7ML6TashfyTtJeJXdDd0Aozrq0BSDbcR9VfHT4SYe5bA3DtplSBw4cMMHHNgylwaGsLrnkEjPEVa9evTLtFNC+fXv5448/TFhyV+3R/iJtkL7rrrvMfQ1dem5cXJyEEio9gLe5CSA1K0fJ1dfcIJ880FWWP32l/KNvG2nToPhFEbUi5Dg0dt/ZGXLrwisKAo8ukNibdX8Qgkox69KbdAjozjvvNL05jgv3HT582PTo6LCT9onOnTu3zN9Lv09MTIwJJNrIrHtPaS+PNlNriPFUjx49TJXqpptuMs3O+jx6ffPmzTOPt2rVyhzXpmodutNm6oMHD0qoIfQAfta4ZiV5sOf5Mu+xK2TeY93l7z3Ol0Y1nDfTK2rdH9ctNHqvvlimLd8pJ84ULqsDoTJU7DTr0hZ8zjHrsrzpVHPbEJS68MIL5c033zRhp0OHDqZp2LExuLR0+rjO9GratKnceOON5vvojC/t6Slp5efzzz+XhIQEMzVdKzhPPfWUfQLPM888Y6pKOqSlAU6n5+s0+FDDhqMO2HAUgUIbltf8fkxmb9wn32zeXyjEbIi+3ynw6KwwbZK2qVQhQvpf3EjuvrS5xDVi81yEyC7rbN5raWfLYZd1enqAABQeHiZdWtQxt5HXxcniXw/Llxv3yfwtB2VQ7uxi1/1RulbQR6v3mJuuIn1312ZmXSGdXQYELS/NuoR1EHqAAKdB5eq2DcwtfdFrEr3Is3V/bNbu+tPcYqpGya0JsXJHl2ZmSA0ArIaeHiBY7Fop0Yuet9892e1fMu3yhTI+PG+2RXHr/qgjaRkyMXm7dH9loVn/R1eRNuv+AIBFEHqAIG3mrJ74lAy7urU88M/xktJ2mMfr/mjO+XPLYrln6mrpPWaxTFm6Q06czu8ZYoFDACGMRmYHNDIjKBTTzPlLdFt5f9UumbV+r5zKcL+tymORn8ljkV+Y2V62obDKURHyRuxiSdz7Zl6w0vWFAC80oeommOdaKA9w58yZM2aj1bI0MlPpAUKombNNg+ryQv92suqfveX569tKq3rOu7vr0JcGHttQmDY/q3uyZ+UFHrX4Zfl9/Xwv/xCwGtuqw7rPFFAatr87ritYlwSVHpcNR3UFynMlRSBY6Oqwq3YcM9Wfb38+IFk5uSboFLUJqrJVgHq2risPXHG+XNqitn2FWaAs9u/fL8ePHzerC+v6M/y9gqf/jmng0f3IdKPVhg0bFjrH00oPoacULxoQjA6ePCsfrd5tbjec+tQp+Ng4DnnZdGhSQx7ocb70adtAIsLDWCsFZXrz0q0aNPgAJaWBRxdNdBeWCT2lQOiBFWRm58j3KQel+xcJUi03tcgFDl01r1NZ3mg4T9pvm+S8Kq7jvkj0A8EDWlV3txknUBQd0tId6IvC4oQA3P/jEREu1574WMQh8Lhb4NBVzLH10v7UpLw780eaBRAr9XrCeSNINjyFB/TNq7g3MMBbaGQGrMbdbtX5HJubz7XhaaXF/5bT/44tvBEkq+MCCFCEHsBKPNitWoNP78rbPdrwtHL2Sfvnqd2ecR7yAoAAQ+gBrMTD3aonDk+SF/pfJM3qVHYbfLT/x5He75h8kTz75U+y/8QZ7/8cAFAKzN5yQCMzLMPD3aqzc3LNVPdJi7fLj3+cMMdcp7y7zvyKigg3e3w92PN8aaR7fLEzNgAvY/aWF180wIpTjVfuOCp75rwktx5/95xr/CgNP5OafCtXHpzGbC8AXsWKzADKja6LcVnkb06B55Ws28wUd8ceH8cNT9vnpOQFHjV/pJz4/tW8z11ne7HfFwAfoacHQKn6gQYMeU1uuLixvJNT0NzsuOGp62yvGstfZLYXAL+ip8cBw1uAB1x6dHYcTpMJC7fJH5sWyOqcwju8F9UDdKr7M1Kl95O85ADKjOEtAN7h0gDdom5VGXNrvLw89AG58ZLGojtVeDLbK2FxO3ntu1/l5FlW5gXgGwxvASgXJvzcEi8LHu8pN13SJG+frvxKj2Ozs9L7urP7+IXb5IpXk+XtxdvlbGY2fxIAvIrhLQcMbwHl5/cjp+TnT56Tvxx626PZXvWrR8ujvVvJLZ1izVYZAOAphrcA+FXzUz86BZ5zzfY6eDJdZs3+TK4as1i+3LhXcnJyC56MGV4AygH/nQLgs9lef2nX0GkrC8fZXo9FfiafRT8vfY9/LEM+3ijXvrFUkn85JLlLx4pM6yuSPIo/KQBlElLDW8ePH5fExETJysoytyFDhsjgwYM9/nqGtwDvz/ba/McJGf3dr3L6t6X2wKPVHg08rsNehWZ+DZrHhqYACrHkiszZ2dmSnp4ulStXllOnTslFF10ka9eulTp16nj09YQewHdW7Tgqr877RdbvPm7uuwYc1/6fOXXvl4vveE6a1Cq8HxgAazvpYegJqeGtiIgIE3iUhh/NcyGU6YCQcmmLOvL5g5fJuwM7SZsG1Qrt4O7a8PzInp5y5WuLZdQ3WyTtt6Xun5TeHwDFCKjQs2TJEunXr580atTILHs/e/bsQudMnDhRmjdvLhUrVpQuXbrI6tWrCw1xdejQQZo0aSJPPvmkxMTE+PAnAFAS+nve+8L68s2j3eX12+Llmxq3uV3TxzbDKyMrRyqteFWqfvBXWfP+s+a+nW5vQe8PgGAJPTokpYFFg407M2fOlGHDhsmIESNk/fr15tw+ffrIoUOH7OfUrFlTNm3aJDt37pQPP/xQDh48WOT302qQlsQcbwB8Lzw8TK6PbyzJl25wu6aPDn3Zen8ei/zCfJ6w7XWZ+soQ+frH/XnNzuznBeAcAranR/8HOGvWLOnfv7/9mFZ2EhISZMKECeZ+Tk6OxMbGyiOPPCJPP50/S8TBQw89JFdeeaUMGDDA7fcYOXKkPPfcc4WOs8s64AeOG5GKyJnIGlIp68Q5m5tde390pph0G+q76wbgdyHX05ORkSHr1q0zs7NswsPDzf2VK1ea+1rVSU1NNZ/rD67DZa1bty7yOYcPH27Os9327Nnjg58EgNteHIfAo8Gl0jO7JbXbM4XW9Cmu9+ereg/IgXYP8gIDCO7Qc+TIETM7q379+k7H9f6BAwfM57t27ZLu3bubYS/9qBWgdu3aFfmc0dHRJhE63gD4f00fW6WmWuKTefdFZHaNu+1T3Ivaz+vh3T2k138WyfgFv0n6juXuvxfNzoBlRUoI6dy5s2zcuNHflwGgNHoNF2nRs/A6PBqAYi+V/s26St1tR+Slb7ZI94PvF9n7MynzOslOfkmil34hKW2HyYUDnjXD5U5DaBqw9PsBsJSgCT06C0unpLs2Juv9Bg0alOm5tXFab1pJAuBHroHH5fjlLWNkzsVrJXyB+54eHQJrEnZY7opcYO7H/TxG3t/1p7S/baS0/32ac7Ozu4AFIKQFzfBWVFSUdOzYURYsyPvHzNbIrPe7di3bP1xJSUmSkpIia9asKYcrBeA1u1ZK+IKC3p/lzR+Wy3Lederx0cDzflbvgvtp0yR2cttCPUMEHsB6Air0pKWlmeEp2xCVTjvXz3fv3m3u63T1yZMny4wZM2TLli3y4IMPmmnugwYN8vOVA/BH78/l974oyU/0lEMdHnTaz+uZrPuKbHZe12qIZF/2GH9ggAUF1JT1RYsWSa9evQodHzhwoEyfPt18rtPVR48ebZqX4+Pj5Y033jBT2csD21AAwbmfl9q457h88vkn8uGBxvZjG6Lvdwo8OhSmO723bVRdnr/+IunYrJZPLxuAd1hy763y6OnZunUr6/QAQUr/Ofvfpn3yytxf5Lq0T5w3K3VZ70envze/JFGevqaNxFSNLjZQAQhshB4vvmgAAlvm4jFSIfm5IhcwXJbdVrpF/GwC0AcVbpTHr7pA7rq0mUSufJ3ZXUAQCrnFCQHAI7tWOgWe2XXuN0Najj0+GniUVoLuzPxCRs5Jkff+M5StLIAQFzRT1gGgRM3OOi09caT07zZUYn47Is/+r4rIsbygY6v0KL3/QORXUutMQSUorfszUpUhLiDk0NNDTw8Qmlx6c3RH9qnLd8qyBXNkWUarQnt42WhFKCUyTnpdfZ3crUNeEQ4Fcfp9gIBET48XXzQAwWv/iTPywtdbzO7s7mZ3zci+2uzkruFnUd07zSyvzufVZjVnIIARerz4ogEIfr9/+aI03/BqsefYZnqNj10s/Q6/XfDAoHnM8AKC8P2bnh4A1rNsrFPgOZ5bVWq67OXl1O9zuOCx3I6DJMxdvw9DX0DAY/YWAGvRcOKyJcXZYdvly5j73Z7uOtU9bN00OfrtK84n6Uam0/qKJI/y2mUDKDtCT34jc1xcnCQkJPB3CrDYVha6i3uDGhXl+odHy66Ln7JvZaH9PY5O5layz/iqs/IlWTb9X5KelV3Q66N0xpiGKgABidlbDujpASykiOGorJ3LZf2yedJ5+xuFHkvO7iC9IjbZ75+UalJdUgtOyA9RAHyLxQkBoDhFrMMTuXe1U+BxrPho4NHgY+MYeM70+D8CDxDgGN4CgGL6fdbftl7ejLzbKfik5lZyes00GHVf3kGWLvif2f+rEIa8gIBA6AGAYvp9el9YXwY+OU4WNH7I3sxcLexMoWbncekjpPvSu+WT1x+X3UdPFzxIkzMQMOjpcUBPD4Di+n2OfvyQ1PnlgyI3MrX5T87tUuXKJ+X+8C8lYmHBPmCs7wN4Bz09JcDsLQBOiliHxzHwvJZzR6GNTG2eCP9Ibku+wjnwaOWI/bwAv2J4S0SSkpIkJSVF1qxZ498/DQBBM/R169DXpHebembFZlvw0aEvG8cK0NJmD0vGpUN8f80AnDC85YDhLQAlGfrSpuW5Px2QEf/7WZqlbZK1uW3c7uc1OGOYpNZLkFcHtJcOsTXdPheA0mN4CwC8wSGkhIWFybXtGsr8YT2kdeerzc7trj0+ev+z6Oel15EP5IY3l8vLc3+Rs5n5ixqyijPgU+y9BQBlVKNSBXmx7nyRCh8X2eSs+3hJpsikxddJ7Q0T5f6M9wpWcW7Rk4oP4AP09ABAOa/vs6RpknTKnFyoyVmDjw5/2QOPiGT0etY58LCmD+A1hB4AKOcm5yv+9pJ8mXS5LK53V6Hg41j90dWdE3+4RFZsP5J3gCEvwKtoZM6fsq637Oxs2bp1q5w4cUKqV6/u3VceQOhxaUzOzM6RSYu2y/LkOfJWxGtu1/TRUKQzwCadt0T67p9U8MCgeQx5AeXcyEzoKcWLBgAlcXjuy1L3h1FFPl5okUM2LgVKhNlbABAIlo11CjyOG5jaOAae97N6y4sn+uTN8HJErw9QZvT0AIAPNzA99vCv8n7VQW5PP5NbQe6KXCARK16XfuOXyeY/TuQ9QK8PUC4Y3nLA8BaAcpc8Km9ausOQVXZOruyf+BdpcnR5kV+mvT5Tcq+X6Rcsl26/Tyh4gF4foBB6ekqB0APAK1xXXtbKjUMFqKiNS8/Z68OKzoBBTw8ABArXdXgcAk9O75HyWe+lMjr79kJf5hh4Fud0kMk515sqkcGQF1Bi9PQAgB/X9AnvPlQGX9FCrn94tEyrfK85fCY3qtCX9QjfJEe/fUVun7xKjn/3akFw0qEzmpwBj9DT44DhLQA+42ZoStf1SZl8n3Q48EWRX1bskBfDXbCokx4uOUOlJ39xwri4OElISPDlnxEAK3Ozu3qFla87BZ5zTW//MuZ+OdHx4bw7DHcB50SlxwGVHgB+o1Ua3XU9X2avEfJKal+JXPW6PB1ZsJGp6wyv/1W9RT5uu0qarn+14AFmeMFiTrIis/deNADw1fT2VTuOSu77N0nXnA1uv4QZXoAwvAUAQafX8LwqjcO09Ev3vecUeFyHvByHu9ZW6CjbWw8ueJAhL8AJPT0AECTT23+96HHpGT6t0M7tNp0y18ns8U/Khz/sltylDmsBMcMLMAg9ABAk09tbD3hW5j3WXZrVqVzklzwe/qFc801XCVvgvP2Fu8ZpwGoi/X0BAIBzDHm16GkPLQ03T5LbT04tdjVnx/t7LvmHxLKKM2BQ6QGAQGer0rgMdx3s/LTcXOPDIoe7TudGSY+VHeSNBb/lreRMjw8sjtADAEE63FX/2uEy5+FucrrzI5Kc3aHQ6ZXDMuTdyFdkzPdbZebYYfT4wPIY3gKAIB7uqhQVIc/X+V4kYpNThUcDjzk9YpNsCb9XKqXm3S/U48MqzrAQKj0AEEIzvD6ufp/EpU93qvxUyg9ASofCXjjeRzKychjuguUQegAghIa8bnnsNXmqb2v5f9lPu924VE1ZtlM+HDOU4S5YDttQOGBFZgBBy2WY6o85L0mTda+UbhVnIMiw4WgJsOEogJAa8lo21inwFLeKsw53/ds23GULT0CIotLjgEoPgFDbuDS390j5oMKNsv/rUfJkxEdOp2rDs/b/qEua1pTprZZL9eUv5g2XacM0ECSo9ACAFbn0+IR1Hyp3XdpMBl7WvNCpOsNrWoW8itCN+17LCzyu21ZQ+UEIYco6AIT4tHYd7qr3w6gip7T/Fn6XVAjLH94SkSVNk+SyJl0kUhcz1JlhVH4QIpi9BQAWWsVZh7vm/HWdLMkpmNLuGHh0qvs9Wy+XT8Y9zuwuhBxCDwBYbLjr1oSmEvP3r+SsFJ7SrpWfDdH3yx2pU4tezBAIUoQeALDCcNegeU7T0uN2vCsVxWGV5iJmd43Ovl0+jrop7w57dyHIEXoAwIJT2h2HvM5WqFHkl+lGpU9/sVnmTvoHw10IeoQeALASlx4fHbqq+K/dcrJxT7enP13hY/kx+j655sAkp68xIYqhLgQZQg8AWLjHxwx5LRsr1fcusp/iun1F9bAz9s/Hh98pqxrdw1AXghKhBwCs3OPjUvnJ6ThI3ui61KzU7Eqnur92+i+yber97oe6qPwgwBF6AMDKPT4ulZ/wfuPkH33byA0XNy70Jbq2z2/Rd8ldkfPtxzJ7jch7DpqcEQQIPQBgda6zu5aNldY/veZU4SlqTZ/rNibIn9+9SpMzggKhBwBQ5GKGZ3s+K4+2mFuoz8e2ps/MYzdLrRX521compwRwAg9AIACLsNdFXs+Lu+0WCqV8retcOXY5LzhgsfsjdFm09Pkgq0vgEDALusO2GUdAKSg4mPr1XGo/KRJJakqBUHHdcf2yS2WylX73ip4QIfNHNcIAryAXdYBAKVnW4fHZU2fIw9vlzWRl7htcv41+h7nwMNQFwIMw1sAAI/X9Gn+yzuSkLXefkpmboT98+iwLPvn70TdIzvb3M9QFwJKSIWePXv2SM+ePSUuLk7at28vn376qb8vCQBCdk2fjfVvlFbp/5X03EinL9H736U2ly8nPFl4Vhdr+cCPQqqnZ//+/XLw4EGJj4+XAwcOSMeOHWXr1q1SpUoVj76enh4AOAdtTtYAk1/5Sfn0OYn7ecy5XzY9X2kI0uqRhimgnHj6/u0cz4Ncw4YNzU01aNBAYmJi5NixYx6HHgDAOWhYadHT3uTsGHi0wuM4xGXzflZvqfvzAemzP3//Lg1N+hyKJmdYdXhryZIl0q9fP2nUqJGEhYXJ7NmzC50zceJEad68uVSsWFG6dOkiq1evdvtc69atk+zsbImNjfXBlQOAtZucU7s9IzfX/dLtej53RS4oCDy2qs+eVUxrh7VDz6lTp6RDhw4m2Lgzc+ZMGTZsmIwYMULWr19vzu3Tp48cOnTI6Tyt7txzzz3yzjvvFPv90tPTTUnM8QYAKHmTc7XEJ+XzdquLXM/HZmqleyX1bBYrOMMvAranRys9s2bNkv79+9uPaWUnISFBJkyYYO7n5OSYSs4jjzwiTz/9tD3IXHXVVTJ48GC5++67i/0eI0eOlOeee67Q8XONCQIAil/PJzW3klRzWLjQ5mRuJacFDe07vdueByiFkFunJyMjwwxZJSYm2o+Fh4eb+ytX5u3wq/nt3nvvlSuvvPKcgUcNHz7cvEC2m87+AgCUfT2ffe0ecnu6Y+DZGf8UKzjDp4Im9Bw5csT06NSvX9/puN7XmVpq+fLlZghMe4F0BpfeNm/eXORzRkdHm0ToeAMAlHE9HxGnDUvd0RWcr1p9sfz48UiGuuAzITV7q1u3bmbIq6S0h0hvGqoAAGWY1aV03618Jy7/lyxZvU76Zc4rtILzzxXuluhfstyv4MxQF6xc6dHp5xEREWYdHkd6X6enl0VSUpKkpKTImjVryniVAGBhGlRcqj41KlUoFHhsHKe3f1XvAcm4dIjzCs4sZAirhp6oqCiz2OCCBQvsx7Sqo/e7dqX5DQACbhXn2Euden3m1H1ABqQ/W2gFZ/XT3hPy+RuPOw91sVM7Qnl4Ky0tTbZt22a/v3PnTtm4caPUrl1bmjZtaqarDxw4UDp16iSdO3eWcePGmWnugwYN8ut1AwBc2IantOqTv4LzXy9/TGpM+5dE7y68gOHTFT4WcbdqiG0hQ4a7EGqhZ+3atdKrVy/7fQ05SoPO9OnT5dZbb5XDhw/Ls88+a5qXtVF53rx5hZqbS4qeHgDw/grOYcvGyhW7J55zBWcn9PnACuv0+AN7bwGAl2h/jkOD8474J+W6DQmyRu4qckHDQ12GS71rni5YA4g9u2CVdXoAAEHMpcG5Rf9nZEGXDUUGHp3SnvjDxbJnzktMaUdoDm8BAKyzWWn91aOKPFWntK/JvUOi1zGlHeWHSk9+T09cXJzZ4gIA4OMVnIvg2O/zY5uhrN6MMqOnxwE9PQDgI7oOj87McvBR9b/JrqOn82ZyOcjIjZAL0v8rM9uuki7b3yh44Nr/iHQe7KsrRgCjpwcAEPhr+Tj0+dz46H+kdYNqhU6NCsuW36Lvdg48LRNFvnkiLzwB3q70ZGZmmmnjp0+flrp165q1dIIdlR4ACJyd2rXCo4HH1bbql0rLk6sKDmh4Yh0fSzvpjdlbqamp8tZbb0mPHj3MkzZv3lwuvPBCE3qaNWsmgwcPZisHAECZ+3yWNnvYDGll5kY4nar/TXcKPI7r+ADn4HHoGTNmjAk506ZNk8TERLOTua6WvHXrVlm5cqWMGDFCsrKy5Oqrr5a+ffvKb7/9JsGCRmYACKwp7d0HvSift/tBKrhUesLCxGm/rqyuLvt1AeUxvHX77bfLM888I23bti32vPT0dBOMdK+sv/3tbxJMGN4CgMAc6tJ3KsfAk5UbLi3T35dxTRZJ/yPvFDxAc7MlnfRweIvZW6V40QAAvlu9+VC9blLv0LJCp2nwiQzLcW5u3jaflZst6CQrMgMAgn6oq2WiU+DRoGPjGHg2RnfKCzxKp8KvnuzDC4alFid8/fXXy+NpAAAomNKuQ1W2ICMiO+OfknY5HzkFH9vQV3z62oIDTGeHN0PP5s2b5YEHHpDs7LyGs5SUFNMDBABAqenCgw7Nzef1/5d813m985CWS3MzFR94fe+tKVOmyNixY82sLe2J+f333+Xpp/P/ogbJ7C292UIbACAw9+uKXfdKkc3NRVZ8Th3Jex5YXrk0Mq9Zs0ZGjRolu3fvlj///FMWLlxo1u0JNjQyA0BwNDenxfaUqnsWFXm6VnycApAuYKhYxDAk+bSReejQofL3v/9d1q5dKx9//LH0799fli9fXh5PDQBAoeZmx8Dj7r/uToFHFzDcs4q1fFC6Ss+hQ4ekXr16RT6+d+9eueWWW4Iu+FDpAYAAp7OydMjKw4rPx9Xvk5s6NpEKyc8VHGQtn5Dj1UrPgAEDiux/0VWZGzduLAsWLCjNUwMA4FlzswcVn+tO/Nc58DCzy9JKFXpq1qwpjz76aKHjR48eNVtUqIoVK5b96gAA8GA6u1Z8HJuabSqHZdg/zzm/t/NaPuzXZTmlCj3vvfeefP/99zJ16lT7sS1btkjnzp2lSpUq5Xl9AACUqOJzJjeq0Ok5uSLh2xc49/nAckpd6fn888/lySeflNWrV8u3334rXbt2NQ3Mc+bMKf+rBADAg4rPoS7D5Z3wmwudGu5QBcrpnR942KTUcjxep+fGG2+U+Ph4+61du3YyYcIEufbaa+Xs2bMyfvx4GTRokAQj1ukBgCCu+Og6PDpclThSdIrNkNwPijxdq0DLthyUq/a9lXdAv65KTN7zIOR5PHtLqzobN26UTZs2yZEjR6RWrVrSoUMHc/+mm26SpKQkiYuLkwoVKkiwYvYWAAQpW3+Ow1o+S3PipXv4xuK/jk1KQ4JXd1nXKekagBxvO3bskMjISGnTpo0JQsGI0AMAQS55VF71xhZmirGr1mXS7M8VBQeYyh60PH3/LtU2FDolXW9/+ctf7MfS0tLslSAAAPzW56PDVQ5r+YzOvl2Swj93msmlnAIPW1ZYgseNzLrFRHGqVq0q3bp1M8NctmoQAAD+3qj0+vhGhQKPo0P1uztPZdcFEGHt0JOQkGB2Utd9toqiZaXJkyfLRRddZGZ3AQDgt4pP/n5bF2x+zWnquiNt8Kh3cGnBARYvDGkeD2+lpKTIiy++KFdddZVZeLBjx47SqFEj87luMqqP//zzz3LJJZfIq6++amZ1AQDgV/ML1uPZF3O5NDrivD2S44KGqU16SjXHio9td3eEjBI3Mp85c0a+/vprWbZsmezatcvcj4mJkYsvvlj69OljqjzBikZmALBGY7O+8zkGHtf7ZvHC2EsJPUHCq7O3QhWhBwBCf5PSbdUvlZYnVxV5elq3Z6Rqxci8KpH2BlHxsfaGo6FGFyfUNYa0bwkAENpbVhQXeDJzI2Tm2j0Fw2JaJWLl5pBR5tAzZcoUM7Sle27pTXt9pk+fLsFEZ5xpT1JxTdoAgNDasuLHip0KnVYhLFvuOzuj8NezQWlIKNU6PY6B56233pKxY8eaBmYdKduwYYNZvVk/D9ZtKQAAIb5lRctEaX+OxQudsEFpSChTT48OB82aNUuaNGnidHzPnj1mr65gq5zQ0wMA1uvxmRw9UAac/VxqhaUVH3h0yCv+TpH+b/roQhFQPT2nT58uFHhUbGyseQwAgEBfvPCWTk2KDDx/5laVzXtPFPT4bPxAZPZDPrxYlKcyhR5do6co0dHRZXlqAAB8snhhjeUvFnmahqF2W8Y6H9TgY9vgFNbp6dF9turVq1fouI6YaYkJAIBgWbzQY/T3WLPSk5WVJYcOHSp0O3z4sGRkFL3PCQAAfqerLduGuRy8nHmbmbruzqnuz+R9wjT2oOS1dXp++OEHbz01AADlO8zl0ONzaYs6Zuq6qzO5UTJrw17nNXzYnDSoeG1F5qZNm55zZ/ZAw+wtALAw7dPZs8rzIS/b1hYamDQ8IeDfv8vU03PLLbe4Pa456tixY2V5agAAfM/DwLO79mXS1HFz0ioxebPCENDKFHrmz58v//3vf6Vq1aqFQs+SJUvKem0AAPi2x0fX4dHZWfneirxbBmbOlMphzn2qTY+tcK746Lo/uvAhFZ/QDT09e/aUatWqyRVXXFHosfbt25flqQEA8D3bwoMafBJHyoDUdKn8Q9ETc9Jie0pVx4oPm5OGVk/Prl275Mcff5T69etL586dJVQ2HNVbdna2bN269ZxjggAA6/X35OSKhIcVnKLvnmEO9w36ewK6p6dEoeejjz6Se++9VzIzMyUsLMxsNDp37lypW7euhAIamQEA9tCj09Lz7al9mcQ6DmkVRzc2pb8n+LeheO655+SOO+6QX375Rb777jtz7OmnC69xAABAyKzh0zLR88Bj6+9JHuXVy0PplKjSExUVZYZ/mjdvbu5r+OnYsaOcOnVKQgGVHgBAcZuTbojqJBdnrHX7Ih2u313qHlxacEDX/9HwhOCs9OgKzJUrV7bfb9OmjeTk5MiBAwfKdrUAAAT65qQtE4sMPMop8Kgdi7x8cfD6iswzZsyQFStWSFpa3o60kZGR7KgOAAhdOg1d+3Rss7Q8pbO52Jg0eENP9+7d5YUXXpBu3bpJzZo1pVWrVnL27Fl59913JTk5WVJTU713pQAABELFx1NsTBrcoWfx4sVmvOzXX3+V999/X2644Qbp0aOHvPXWW9K7d2+pVauWXHjhhd67WgAA/EXX4HHwa9UuRZ66vPnDzhuTUvEJrb23du7cKWvXrpUNGzbISy+9JMGIRmYAQLE0wOiwlW3fraJOy+4gvSI2OR9kDZ/gWqcn1BF6AAAlndFVIszoCp7ZWwAAWF5p+ntsmNHlV4QeAADK2N+zpuUQOZlb6dxfx4wuvyL0AABQlhWbE0dKp2a1pHrYmXN/HTO6/IrQAwBAadfv0R4dbZBdULAxaVFSuz3jPKMLPkfoAQCgLBx2Yi/OttXzCs5lmMsvCD0AAJTHMNc5FNrCgqZmnyP0AABQHsNcJZ3RpdUenf4Onwm50KOrROvK0AMGDPD3pQAArFTxcZnRtTL84uK/Rhc41PV+Zj/k3WtD6IaeIUOGyHvvvefvywAAWHxGV/OOfYo8dXOlhIIVnTd+QPDxkZALPT179pRq1ar5+zIAABaf0dVwzctFntbuzBrnAxp8GOqyVuhZsmSJ9OvXTxo1aiRhYWEye/bsQudMnDhRmjdvLhUrVpQuXbrI6tWr/XKtAACUdUZXoaEuprJbJ/ScOnVKOnToYIKNOzNnzpRhw4bJiBEjZP369ebcPn36yKFDh0r1/dLT081+HY43AAC8MaMrp4idLnPO710w1MVUduuEnmuuuUZeeOEF04zszpgxY2Tw4MEyaNAgiYuLk0mTJknlypVl6tSppfp+o0aNMhuU2W6xsbFl/AkAACi8TcXRht0lPKzwK6NbfodvX1BwoNOgvNCE0A89xcnIyJB169ZJYmKi/Vh4eLi5v3LlylI95/Dhw82OrLbbnj17yvGKAQCW5VjtaZkodfYvdXtaWJjLENfaaQxxeVGkBIkjR45Idna21K9f3+m43v/ll1/s9zUEbdq0yQyVNWnSRD799FPp2tV9ao6OjjY3AAC80tRcJSavV+ccsiRcIh2HuLRSRMXHuqHHU/Pn5/+lAQDA3zoPFtm3IW92Vr4dOQ2lRfh+p9MiJafwas2EHusOb8XExEhERIQcPHjQ6bjeb9CgQZmeWxuntUcoISGhjFcJAICL/m+KxN+Z93nLxEKBxy1Wa7Z26ImKipKOHTvKggUFDV85OTnmflHDV55KSkqSlJQUWbPGZd0EAADKK/hc+5+CWVrnwmrNoT+8lZaWJtu2bbPf37lzp2zcuFFq164tTZs2NdPVBw4cKJ06dZLOnTvLuHHjTO+OzuYCACDgh7pOHTFVnDPNekmlXcluT9tSpYtc6Lhasy00oczCcnN1wlxgWLRokfTq1avQcQ0606dPN59PmDBBRo8eLQcOHJD4+Hh54403zCKF5UHX6dGp6zqTq3r16uXynAAAOJkzVGRdCZda0VWe6fEp8/t3QIUef9GeHr3p7LCtW7cSegAA3rFrpci0viX7Gl2756/j+BMpBqGnFKj0AAC8Trea0EZlB1p+cFqzx7G3R4e6dM0fprGX+f07aBqZAQAICS6rNR+p2c5t4MnKDXfenkIrROzNVSaEHgAA/Lhac8zxzW5PiwxzWbtHsTdX6MzeCoSeHgAAAmm1ZifszVUmNDI7oKcHAOBTsx9yWq35TG6UVArLcH+uY3+PhibY0dMDAEAQrdace35ikYHnrEQ59/foLDCUGMNbAAD4O/hEREtYMWv3VBSXMMTeXKVCIzMAAP6kVZuSLlZItadUCD1sOAoACJTZXCWh1R6UCI3MDmhkBgD4teKjQWbxy3IqJl6qHNl47q/RTUx1Ty+LO+nh4oT09AAAECgVH72lHZIqngx32XZi101Mmc3lEYa3AAAIsv6e3ZXimM1VCoQeAACCrL+n6ZkU5wP093iE0AMAQADvzeURZnN5hNDD7C0AQCBhNpfXMHvLAbO3AAAB09uju6rnO1CtrTRI/fncXzdoXl5ospiTHs7eotIDAECA78TuUeBR9PYUi9ADAEAg0mnoug6Pbc8tT9DbUyxCDwAAgUoXHizpas1a7WFDUrcIPQAAhNJsLq32aD9Q8ihvXVHQIvQAABCKs7kY6iqE0MOUdQBAKK7d02mQJWdyFYcp6w6Ysg4ACFg6XKXVG0/ovlzaAK0VIgvsy3WSKesAAFiv2pMZVoF9uYrA8BYAACHU21MhN9P5AGv32BF6AAAIFuzLVSaEHgAAggX7cpUJoQcAgGCijcm6x1b+UNf+am3P/TVMXzci8z4AAICgqvjoLe2QNFw31bOv2bHI8lPYqfQAABCMdKsJTwOPWvyy5benIPSwOCEAwCr9PTsWiZWxOKEDFicEAARdtUf32SqJQfNCbpiLxQkBAAh1VHtKhOEtAACstgv7rpViRYQeAACCGdUejxF6AAAIdlR7PELoAQAg2FHt8QihBwCAUEC155wIPQAAhAKqPedE6AEAIFRQ7SkWoQcAgFBBtadYhB4AAEIJ1Z4iEXoAAAglVHuKROhhw1EAQKih2uMWoUdEkpKSJCUlRdasWeP+VQIAIJhQ7XGL0AMAQCii2lMIoQcAgFBEtacQQg8AAKGKao8TQg8AAKGKao8TQg8AAKGMao8doQcAgFBW0mpPp0F5XxOCCD0AAIQ6D6s9Z+pdLLJ2mkjyKAlFhB4AAEKdh9WeSoc25H2y+GWRXSsl1BB6AACwghaeVXvsdiySUEPoAQDACpqVsLcnBKs9hB4AAKyihbWrPYQeAACsWO1p3Mly1R5CDwAAVtJruEjHv4nsXWu5ag+hBwAAK9m1UmTdVM/PD6FqD6EHAAAraVbChuYQqvYQegAAsJoWJWxoDpFqT8iFnq+++kpat24trVq1kilTpvj7cgAACDzNrFntCanQk5WVJcOGDZOFCxfKhg0bZPTo0XL06FF/XxYAAIGnhfWqPSEVelavXi1t27aVxo0bS9WqVeWaa66R7777zt+XBQBA4GlmvWpPQIWeJUuWSL9+/aRRo0YSFhYms2fPLnTOxIkTpXnz5lKxYkXp0qWLCTo2+/btM4HHRj/fu3evz64fAICg0sJa1Z6ACj2nTp2SDh06mGDjzsyZM83w1YgRI2T9+vXm3D59+sihQ4dK9f3S09Pl5MmTTjcAACyjmbWqPQEVenQ46oUXXpAbbrjB7eNjxoyRwYMHy6BBgyQuLk4mTZoklStXlqlT89Yb0AqRY2VHP9djRRk1apTUqFHDfouNjfXCTwUAQABrYZ1qT0CFnuJkZGTIunXrJDEx0X4sPDzc3F+5Mu/F79y5s/z0008m7KSlpcncuXNNJagow4cPlxMnTthve/bs8cnPAgBAwGhmnWpP0ISeI0eOSHZ2ttSvX9/puN4/cOCA+TwyMlJee+016dWrl8THx8vjjz8uderUKfI5o6OjpXr16k43AAAsp4U1qj1BE3o8dd1118nWrVtl27Ztcv/993v0NdpDpMNlCQkJXr8+AAACTjNrVHuCJvTExMRIRESEHDx40Om43m/QoEGZnjspKUlSUlJkzZo1ZbxKAACCVIvQr/YETeiJioqSjh07yoIFC+zHcnJyzP2uXbv69doAALBctafToLyvCSIBFXq0+Xjjxo3mpnbu3Gk+3717t7mv09UnT54sM2bMkC1btsiDDz5oprnrbC4AAOCjak/jTiJrp4kkjwqqlzwsNzc3VwLEokWLTBOyq4EDB8r06dPN5xMmTDDbS2jzsjYrv/HGG2aRwrLQnh69aaO09gPpTC6amgEAlpQ8Km/oylOD5vm94qPr7OnSM+d6/w6o0ONvnr5oAACErF0rRab1DcnQE1DDWwAAIMh6e3YEzywuQg8AACj9TK4gmsVF6GGdHgAAiqz2/FmpqYTKLC56ehzQ0wMAgIM5Q0XW5e1v6U5O7fMl/Nj2vIDUa7j4Cz09AACg9HatLDbwKBN4gmiIi+EtAABQdkHQ0EzoAQAAZZ/FFQTVHkIPjcwAAJTPflwBXu2hkdkBjcwAALhfoTmrUUeJ3LdOAnGxQhqZAQBA2fUaLtLxb54FngCv9jC8BQAAyjSLK1h6ewg9AACg/BqaA7jaQ+gBAADl29AcoNUeQg+ztwAAKN9qz3k9AnJrCmZvOWD2FgAARdDKzbS+4jEfzuJi9hYAAPBPtadTYG5CyvAWAAAo396etdPy1vcJMIQeAABQ/tWeAGxmJvQAAADvzOQKsKnrhB4AAFCyas95PYOy2kPoAQAAntMQs3NRUFZ7CD2s0wMAQKn6enQT0mCq9rBOjwPW6QEAwENzhnq+J5eGJN241EtYpwcAAHjHruDchJThLQAAYIltKQg9AADAu1PXdy6m0gMAAIJUEG5LQaUHAABYYlsKQg8AALDEthSEHgAAYIltKQg9LE4IAIAltqUg9IhIUlKSpKSkyJo1a/z2BwEAQFDaFTzbUhB6AACAJbaliPTLdwUAAKGj13CRtEMS6ckqzX5cqJBKDwAA8N22FH5cqJDQAwAAfLsthZ/6egg9AADAt1PX/dTXQ+gBAACW2JaC0AMAACyxLQWhBwAAWGKhQkIPAACwxEKFhB4AAOCfWVw+rvYQegAAgH9mcfl4oUJCDxuOAgDgn2qPjxcqDMvNzc312XcLcCdPnpQaNWrIiRMnpHr16v6+HAAAgtecoedepTlxpEi3oT57/6bSAwAA/LMtReyl4kuEHgAA4J8hrj2rxJcIPQAAwD8NzfNHMnsLAACE3kKFpyJc+m30cWZvAQCAUFqoMOvKEVLl//bkNS/b6OOs0wMAAEKmrydxpEReMSzvc52tZQs++rgPKz1MWXfAlHUAAMqZVnLcBZuijpcCU9YBAID/NSsi2PiwwmPD7C0AAGAJhB4AAGAJhB4AAGAJhB4AAGAJhB4AAGAJhB4AAGAJhB4AAGAJhB4AAGAJhB4AAGAJkf6+gECSm5trX84aAAAEB9v7tu19vCiEHgepqanmY2xsrDf/bAAAgJfex2vUqFHk42w46iAnJ0f27dsn1apVk7CwMKcXKiEhQdasWePxC1+S8z05V1OshrE9e/ZI9erVxapK+ucQatfkje9VHs9Z2ucozdeV5+8Wv1el/3MIpevi90qC/j1LKzwaeBo1aiTh4UV37lDpcaAvVJMmTdy+UBERESX6gyvJ+SU5V8+zcugp6Z9DqF2TN75XeTxnaZ+jNF/njd8tfq8C7/fKl79b/F5JSLxnFVfhsaGR2UNJSUleO7+kz21lgfha+fKavPG9yuM5S/scpfk6fres8Xvly+vi90os83vF8FaQ0FKhptgTJ04E5P/IgGDE7xVgrd8tKj1BIjo6WkaMGGE+AuD3Cghk0QH6nkWlBwAAWAKVHgAAYAmEHgAAYAmEHgAAYAmEHgAAYAmEHgAAYAmEnhCjS3737NlT4uLipH379vLpp5/6+5KAkHHDDTdIrVq1ZMCAAf6+FCBoffXVV9K6dWtp1aqVTJkyxaffmynrIWb//v1y8OBBiY+PlwMHDkjHjh1l69atUqVKFX9fGhD0Fi1aZPb3mTFjhnz22Wf+vhwg6GRlZZn/lCcnJ5vFC/U9asWKFVKnTh2ffH8qPSGmYcOGJvCoBg0aSExMjBw7dszflwWEBK2i6obEAEpn9erV0rZtW2ncuLFUrVpVrrnmGvnuu+/EVwg9PrZkyRLp16+f2QlWd3KfPXt2oXMmTpwozZs3l4oVK0qXLl3MX5LSWLdunWRnZ5udboFQ58vfLcCqlpTx92zfvn0m8Njo53v37vXZ9RN6fOzUqVPSoUMH85fCnZkzZ8qwYcPM8t3r16835/bp00cOHTpkP0crORdddFGhm/5lstHqzj333CPvvPOOT34uwCq/W4CVnSqH3zO/yoXf6Ms/a9Ysp2OdO3fOTUpKst/Pzs7ObdSoUe6oUaM8ft6zZ8/mdu/ePfe9994r1+sFrP67pZKTk3NvuummcrtWwEq/Z8uXL8/t37+//fEhQ4bkfvDBBz67Zio9ASQjI8MMSSUmJtqPhYeHm/srV6706Dn07+G9994rV155pdx9991evFrAWr9bAMr+e9a5c2f56aefzJBWWlqazJ0711SCfIXQE0COHDlienDq16/vdFzv60wsTyxfvtyUF3WcVUv1etu8ebOXrhiwzu+W0n+8b775Zvnmm2+kSZMmBCaghL9nkZGR8tprr0mvXr3M+9Pjjz/us5lb5vv77DvBJ7p16yY5OTm82oAXzJ8/n9cVKKPrrrvO3PyBSk8A0enlERERZp0dR3pfp58D4HcLCFQxQfAeRugJIFFRUWahpgULFtiPadVG73ft2tWv1wYEM363AH7PFMNbPqaNW9u2bbPf37lzp2zcuFFq164tTZs2NVP9Bg4cKJ06dTINX+PGjTNTBAcNGsTvLMDvFuBXacH+HuazeWKwT3fVl931NnDgQPsrNH78+NymTZvmRkVFmel/q1at4tUDzoHfLcD7koP8PYy9twAAgCXQ0wMAACyB0AMAACyB0AMAACyB0AMAACyB0AMAACyB0AMAACyB0AMAACyB0AMAACyB0AMAACyB0AMAACyB0AMgJD3xxBPSv39/f18GgABC6AEQknTn5/j4eH9fBoAAQugBEJI2bdpE6AHghNADIOT88ccfcuTIEXvoOX78uPTr10+6desmBw4c8PflAfATQg+AkBzaqlmzpjRv3lw2b94sCQkJ0rhxY0lOTpYGDRr4+/IA+AmhB0BIhp4OHTrIhx9+KD169JCnnnpKJk2aJBUqVPD3pQHwo7Dc3Nxcf14AAJS3AQMGyMKFC83nX3/9tXTt2pUXGQCVHgChWem58cYb5ezZs6afx9Xll18uP/zwg/n8vvvuk7Fjx/rhKgH4GpUeACElNTVVatSoIevWrZMNGzbI0KFDZcWKFdK2bVv7OfPmzZO3335bunfvLlu2bJHJkyf79ZoB+AahB0BIWbZsmfTq1UvS0tIkOjpahg0bJrNnz5bVq1dLTEyM/byOHTuaZmcNQPT6ANZAIzOAkBvaatOmjQk8avTo0dK6dWsz3JWRkWGOrVmzRo4dO2YqQgQewDqo9ACwlL1798q1115rqj833XSTvPfee3LRRRf5+7IA+ACVHgCWcebMGbn55ptl/Pjxct5558nw4cPl3//+t78vC4CPUOkBAACWQKUHAABYAqEHAABYAqEHAABYAqEHAABYAqEHAABYAqEHAABYAqEHAABYAqEHAABYAqEHAABYAqEHAABYAqEHAACIFfx/SyEVXfNM62oAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.loglog(k_x, P_1D, label=\"Semi-Analytic\", lw=4)\n",
    "ax.loglog(k_x, P_1Ds, 'x', label=\"Numerical\", mew=2)\n",
    "ax.legend()\n",
    "ax.set_xlabel(\"$k_x$\")\n",
    "ax.set_ylabel(r\"$P_{\\rm 1D}(k_x)$\")"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "-all"
  },
  "kernelspec": {
   "display_name": "py313",
   "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.13.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}