{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "e87edff3-c6c0-434f-876b-28f12d376a14",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.integrate\n",
    "from classy import Class\n",
    "from classy import CosmoComputationError\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from tqdm import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "445fc9a0-6bb9-45aa-8ea2-726bada97199",
   "metadata": {},
   "outputs": [],
   "source": [
    "def class_compute(class_params):\n",
    "    '''\n",
    "    A function to handle CLASS computation and deal with potential errors.\n",
    "\n",
    "    Args:\n",
    "        :class_params (dict): Dictionary of CLASS parameters\n",
    "\n",
    "    Returns:\n",
    "        :cosmo (CLASS): Instance of the CLASS code\n",
    "        :isNormal (bool): Whether error occurred in the computation\n",
    "    '''\n",
    "    cosmo = Class()\n",
    "    cosmo.set(class_params)\n",
    "    try:\n",
    "        cosmo.compute()\n",
    "        isNormal=True\n",
    "    except CosmoComputationError as e:\n",
    "        if \"DeltaNeff < deltaN[0]\" in str(e):\n",
    "            # set YHe to 0.25. Value is from https://arxiv.org/abs/1503.08146 and Plank 2018(Section 7.6.1) https://arxiv.org/abs/1807.06209\n",
    "            warnings.warn(f\"Adjust YHe to 0.25 due to CLASS CosmoComputationError for cosmology {class_params}.\")\n",
    "            class_params['YHe'] = 0.25\n",
    "            cosmo.set(class_params)\n",
    "            cosmo.compute()\n",
    "            isNormal=False\n",
    "        else:\n",
    "            raise e\n",
    "    return cosmo, isNormal\n",
    "\n",
    "\n",
    "def get_class_linear(k, As, Om, Ob, h, ns):\n",
    "    \"\"\"\n",
    "    Compute linear P(k) using class for the cosmology of interest\n",
    "    \n",
    "    Args:\n",
    "        :k (np.ndarray): k values to evaluate P(k) at [h / Mpc]\n",
    "        :As (float): 10^9 times the amplitude of the primordial P(k)\n",
    "        :Om (float): The z=0 total matter density parameter, Omega_m\n",
    "        :Ob (float): The z=0 baryonic density parameter, Omega_b\n",
    "        :h (float): Hubble constant, H0, divided by 100 km/s/Mpc\n",
    "        :ns (float): Spectral tilt of primordial power spectrum\n",
    "        \n",
    "    Returns:\n",
    "        :plin_class (np.ndarray): Linear power spectrum at corresponding k values [(Mpc/h)^3]\n",
    "        :isNormal (bool): Whether error occurred in the computation\n",
    "    \"\"\"\n",
    "    \n",
    "    Oc =  Om - Ob\n",
    "    redshift = 0.0\n",
    "    \n",
    "    class_params = {\n",
    "        'h': h,\n",
    "        'omega_b': Ob * h**2,\n",
    "        'omega_cdm': Oc * h**2,\n",
    "        'A_s': As*1.e-9,\n",
    "        'n_s': ns,\n",
    "        'output': 'mPk',\n",
    "        'P_k_max_1/Mpc': k.max() * h,\n",
    "        'w0_fld': -1,\n",
    "        'wa_fld': 0,\n",
    "        'Omega_Lambda': 0,  # Set to 0 because we're using w0_fld and wa_fld instead\n",
    "        'tau_reio':0.0561,\n",
    "        'z_max_pk': 3.0,  # Max redshift for P(k) output\n",
    "    }\n",
    "    \n",
    "    cosmo, isNormal = class_compute(class_params)\n",
    "    plin_class = np.array([cosmo.pk_lin(kk*h, redshift) for kk in k]) * h ** 3\n",
    "    \n",
    "    # Memory cleanup for class\n",
    "    cosmo.struct_cleanup()\n",
    "    cosmo.empty()\n",
    "\n",
    "    return plin_class, isNormal\n",
    "\n",
    "\n",
    "def get_sig(R, As, Om, Ob, h, ns):\n",
    "    \"\"\"\n",
    "    Compute the bulk flow variance when the field is smoothed with a top\n",
    "    hat filter of side length R.\n",
    "    \n",
    "    sigma^2 = H_0^2 f^2 / (2 \\pi^2) \\int_0^\\infty dk W^2(k,R) P(k)\n",
    "    \n",
    "    where W(k,R) is the Fourier transform of the top hat filter\n",
    "    \n",
    "    Args:\n",
    "        :k (np.ndarray): k values to evaluate P(k) at [h / Mpc]\n",
    "        :As (float): 10^9 times the amplitude of the primordial P(k)\n",
    "        :Om (float): The z=0 total matter density parameter, Omega_m\n",
    "        :Ob (float): The z=0 baryonic density parameter, Omega_b\n",
    "        :h (float): Hubble constant, H0, divided by 100 km/s/Mpc\n",
    "        :ns (float): Spectral tilt of primordial power spectrum\n",
    "        \n",
    "    Returns:\n",
    "        :sigma (float): Bulk flow RMS [km/s]\n",
    "        \n",
    "    \"\"\"\n",
    "    \n",
    "    k = np.logspace(-4, 2, 300) # (h / Mpc)\n",
    "    plin_class, isNormal = get_class_linear(k, As, Om, Ob, h, ns) # (Mpc/h)^3\n",
    "    x = k * R\n",
    "    integrand = plin_class * (3 * (np.sin(x) - x * np.cos(x)) / x ** 3) ** 2 # (Mpc/h)^3\n",
    "    \n",
    "    # Check units\n",
    "    sigma = scipy.integrate.simpson(integrand, x=k) # (Mpc/h)^2\n",
    "    H0 = 100 # h km/s/Mpc\n",
    "    f = Om ** 0.55\n",
    "    sigma = np.sqrt(sigma / 2) * H0 * f / np.pi # (Mpc/h) * h km/s/Mpc = km/s\n",
    "\n",
    "    return sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "0e09f7ec-92e7-474c-b9a1-46094ab7385f",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 20/20 [00:11<00:00,  1.73it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "58.40857421174422\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAg6UlEQVR4nO3deXhddb3v8fd3Z9iZ57FpmgRawLZAoekAFB9ksIAei16UekUrcsRzDh5FPc+5crk+F6+Xe5zxeh2OICgqUFBQUAYZRChTQ1pK6UhT0rRpk2ZopmYefvePvRpCSdukGdbO3p/X86xnrf3ba+18f3naz175rcmcc4iISGQJ+F2AiIhMPoW7iEgEUriLiEQghbuISARSuIuIRKBYvwsAyMnJcaWlpX6XISIyo2zYsKHJOZc72nthEe6lpaVUVlb6XYaIyIxiZjXHek/DMiIiEUjhLiISgU4Y7maWYGYVZvaGmW01s2967Vlm9rSZ7fLmmSO2udnMqsxsp5mtnMoOiIjIe41lz70XuNg5dzawCLjczJYDXweedc7NA571XmNm84HVwALgcuBnZhYzBbWLiMgxnDDcXchh72WcNzlgFXCP134PcJW3vApY65zrdc5VA1XA0sksWkREjm9MY+5mFmNmm4AG4Gnn3Hog3zlXB+DN87zVi4B9Izav9dqO/swbzKzSzCobGxsn0AURETnamMLdOTfonFsEzAaWmtnC46xuo33EKJ95h3Ou3DlXnps76mmaIiJyksZ1toxzrhX4O6Gx9INmVgjgzRu81WqB4hGbzQYOTLTQ0exv7eZ7f93B/tbuqfh4EZEZayxny+SaWYa3nAhcCuwAHgXWeKutAR7xlh8FVptZ0MzKgHlAxSTXDcDhngF++txuXq5qmoqPFxGZscZyhWohcI93xksAeNA59xczewV40MyuB/YCHwdwzm01sweBbcAAcKNzbnAqip+Xl0JGUhwV1Yf4eHnxiTcQEYkSJwx359xm4JxR2puBS46xzW3AbROu7gQCAWNJaRYVew5N9Y8SEZlRZvwVqsvKsqhp7uJge4/fpYiIhI0ZH+5LSrMAqKjW3ruIyBEzPtwXzEojKT5G4S4iMsKMD/fYmACLSzIV7iIiI8z4cIfQuPvOgx20dPb5XYqISFiIiHBfWpYNQGVNi8+ViIiEh4gI97NmpxMfE6CiutnvUkREwkJEhHtCXAyLijM07i4i4omIcAdYWpbFlgPtHO4d8LsUERHfRVS4Dw45NmrcXUQkcsL93JJMAgav6VYEIiKRE+4pwVgWFqWzXuPuIiKRE+4AS0uz2LSvlZ7+KbkJpYjIjBFZ4V6WRd/AEJtr2/wuRUTEVxEV7kduIqZxdxGJdhEV7pnJ8ZyWn6JxdxGJehEV7hAamtmw5xADg0N+lyIi4psIDPdsOvsG2VbX7ncpIiK+ibxw18M7REQiL9wL0hOYk5WkcBeRqBZx4Q6hcffX9hxiaMj5XYqIiC8iNtxbuvqpajzsdykiIr6IyHBfVhYad9cpkSISrSIy3OdkJZGfFuQ1hbuIRKmIDHczY0lpFhXVh3BO4+4iEn0iMtwhNDRT397DvkPdfpciIjLtIjbcjzw0e72eqyoiUShiw31eXgoZSXE6311EotIJw93Mis3sOTPbbmZbzezLXvutZrbfzDZ505UjtrnZzKrMbKeZrZzKDhxLIGCUl2TpDpEiEpVix7DOAPA159xGM0sFNpjZ0957tzvnvj9yZTObD6wGFgCzgGfM7DTn3LQ/QWNZWRbPbD/IwfYe8tMSpvvHi4j45oR77s65OufcRm+5A9gOFB1nk1XAWudcr3OuGqgClk5GseO1tEz3mRGR6DSuMXczKwXOAdZ7TV80s81mdreZZXptRcC+EZvVMsqXgZndYGaVZlbZ2Ng4/srHYMGsNJLiYxTuIhJ1xhzuZpYCPATc5JxrB34OnAosAuqAHxxZdZTN33OyuXPuDudcuXOuPDc3d7x1j0lsTIDFJZkadxeRqDOmcDezOELBfq9z7mEA59xB59ygc24IuJN3hl5qgeIRm88GDkxeyeOztDSLHfUdtHb1+VWCiMi0G8vZMgbcBWx3zv1wRHvhiNU+Cmzxlh8FVptZ0MzKgHlAxeSVPD5Hxt1f29PiVwkiItNuLGfLXAB8GnjTzDZ5bf8d+KSZLSI05LIH+AKAc26rmT0IbCN0ps2Nfpwpc8TZxRnExwSoqG7msvn5fpUhIjKtThjuzrkXGX0c/fHjbHMbcNsE6po0CXExLCrOoEJ77iISRSL2CtWRlpRlsmV/G529A36XIiIyLaIi3JeWZTM45Ni4V3vvIhIdoiLcF5dkEjBdzCQi0SMqwj0lGMvConQ9mUlEokZUhDvAktIsNu1rpXfAtxN3RESmTdSE+9KyLPoGhthc2+Z3KSIiUy5qwn1JqW4iJiLRI2rCPSs5ntPyUzTuLiJRIWrCHUJDMxtrWhgYHPK7FBGRKRVV4b6kNIvDvQNsr+vwuxQRkSkVVeF+5CZiemi2iES6qAr3wvRE5mQl6aCqiES8qAp3CO29v7bnEEND73l+iIhIxIi+cC/NoqWrn92Nh/0uRURkykRfuA+Pu2toRkQiV9SFe0l2EnmpQY27i0hEi7pwNzOWlmVRUX0I5zTuLiKRKerCHWBZWRb17T3UtnT7XYqIyJSIynBfonF3EYlwURnup+Wlkp4YR4UuZhKRCBWV4R4IGEtKs3RQVUQiVlSGO4TG3fc0d9HQ3uN3KSIiky5qw/3IuHvFHu29i0jkidpwXzArjaT4GA3NiEhEitpwj4sJsLgkU+EuIhEpasMdQveZ2VHfQWtXn9+liIhMqugOd2/cfd2uJp8rERGZXFEd7uWlWRRnJfKbV/b4XYqIyKQ6YbibWbGZPWdm281sq5l92WvPMrOnzWyXN88csc3NZlZlZjvNbOVUdmAiYgLGZ88v47U9LWyubfW7HBGRSTOWPfcB4GvOufcBy4EbzWw+8HXgWefcPOBZ7zXee6uBBcDlwM/MLGYqip8MnyifTUowll+9tMfvUkREJs0Jw905V+ec2+gtdwDbgSJgFXCPt9o9wFXe8ipgrXOu1zlXDVQBSye57kmTmhDHx8tn8+c3DnBQFzSJSIQY15i7mZUC5wDrgXznXB2EvgCAPG+1ImDfiM1qvbajP+sGM6s0s8rGxsaTKH3yfPb8Ugad47ev1Phah4jIZBlzuJtZCvAQcJNzrv14q47S9p4bpzvn7nDOlTvnynNzc8daxpQoyU7m0vflc+/6Gnr6B32tRURkMowp3M0sjlCw3+uce9hrPmhmhd77hUCD114LFI/YfDZwYHLKnTrXryijpaufP72+3+9SREQmbCxnyxhwF7DdOffDEW89CqzxltcAj4xoX21mQTMrA+YBFZNX8tRYVpbF/MI07n6pWk9oEpEZbyx77hcAnwYuNrNN3nQl8G3gMjPbBVzmvcY5txV4ENgGPAnc6JwL+7EOM+NzK8p46+BhXqrSfd5FZGazcNhLLS8vd5WVlX6XQe/AIBd8+2+cWZTOr64L2xN8REQAMLMNzrny0d6L6itUjxaMjeHa5SU8t7OR3Y2H/S5HROSkKdyP8qllJcTHBPi1LmoSkRlM4X6U3NQgqxbN4g8bamnr6ve7HBGRk6JwH8V1F5TR3T/I2tf2+l2KiMhJUbiPYv6sNM47JZt7Xt7DwOCQ3+WIiIybwv0YPreijANtPTy5td7vUkRExk3hfgwXn5FHSXYSd79Y7XcpIiLjpnA/hpiAcd35pWzc28rre1v8LkdEZFwU7sdxdXkxqbrXu4jMQAr340gJxnLNkmIef7OOurZuv8sRERkzhfsJrDm/lCHd611EZhiF+wkUZyXxwfkF3Fexl+6+sL//mYgIoHAfk+svLKO1q5+HX6/1uxQRkTFRuI9BeUkmZxalc/eL1QwN+X8XTRGRE1G4j0HoXu+l7G7sZF1Vk9/liIickMJ9jD505ixyU4O6qElEZgSF+xjFxwb4zPISnn+rkaqGDr/LERE5LoX7OPzXZXOIjw1wty5qEpEwp3Afh+yUIB87p4iHN9bS0tnndzkiIsekcB+n6y4oo6d/iPt1r3cRCWMK93E6vSCVFXNz+M3LNfTrXu8iEqYU7ifhcytKqW/v4Yktute7iIQnhftJuOi0PE7JSeauF6txThc1iUj4UbifhEDAuO6CUt7Y18rGva1+lyMi8h4K95P0sXNnk5YQy90v6aImEQk/CveTlByM5ZNL5/DklnpqW7r8LkdE5F0U7hOw5vxS4mKMW/64RWPvIhJWFO4TMCsjkVuufB/Pv9XIb1/VwzxEJHwo3Cfo2uUlXHR6Lrc9tl33nBGRsHHCcDezu82swcy2jGi71cz2m9kmb7pyxHs3m1mVme00s5VTVXi4MDO+e/VZJAdj+fLaTfQN6MImEfHfWPbcfw1cPkr77c65Rd70OICZzQdWAwu8bX5mZjGTVWy4yktN4D8+diZbD7Tzo2fe8rscEZETh7tz7gXg0Bg/bxWw1jnX65yrBqqApROob8ZYuaCA1UuK+fnzu6moHuuvS0RkakxkzP2LZrbZG7bJ9NqKgH0j1qn12t7DzG4ws0ozq2xsbJxAGeHjGx+ez5ysJL7ywCbae/r9LkdEotjJhvvPgVOBRUAd8AOv3UZZd9RzBJ1zdzjnyp1z5bm5uSdZRnhJDsZy+zWLqG/v4dZHt/pdjohEsZMKd+fcQefcoHNuCLiTd4ZeaoHiEavOBg5MrMSZ5dw5mXzxA3N5eON+/rI5qrouImHkpMLdzApHvPwocORMmkeB1WYWNLMyYB5QMbESZ54vXjyXRcUZ3PLHLdS39fhdjohEobGcCnk/8ApwupnVmtn1wHfN7E0z2wx8APgKgHNuK/AgsA14ErjROTc4ZdWHqbiYALdfs4i+gSH+7fdvMDSkq1dFZHpZOFw2X15e7iorK/0uY9LdX7GXmx9+k298eD7XryjzuxwRiTBmtsE5Vz7ae7pCdQqtXlLMpe/L5ztP7mBHfbvf5YhIFFG4TyEz4zv/5UzSEuK4ae0megeiboRKRHyicJ9i2SlBvnv1meyo7+AHT+nqVRGZHgr3aXDxGflcu3wOd657m5ermvwuR0SigMJ9mtxy5XzKcpL52u/foK1LV6+KyNRSuE+TxPgYfnTNIho7evnGI1tOvIGIyAQo3KfRWbMzuOnSeTz6xgEe2bTf73JEJIIp3KfZP180l/KSTP7Hn7awv7Xb73JEJEIp3KdZTMC4/ZpFDA05vvrAJgZ19aqITAGFuw+Ks5K49SMLWF99iF+ue9vvckQkAincfXL14tlcsbCA7z+1kw01eriHiEwuhbtPzIz/89EzKcpI5DN3VbD+7Wa/SxKRCKJw91FmcjwPfOE8CtITWPOrCtbtiownUomI/xTuPstPS+CBL5xHaXYy199TybPbD/pdkohEAIV7GMhJCbL2huWcnp/KP/1uA09uqfO7JBGZ4RTuYSIjKZ57P7+MM4vSufG+13WRk4hMiMI9jKQlxPGb65exuCSTmx7YxIOV+/wuSURmKIV7mEkJxnLPdUtZMTeHf//DZn77ao3fJYnIDKRwD0OJ8THc+ZlyLjkjj2/8aYsudBKRcVO4h6mEuBh+fu1irlhYwP9+bDs/fa7K75JEZAZRuIex+NgA/++T57Bq0Sy+99ed/PCpnYTDA81FJPzF+l2AHF9sTIAffmIRwdgAP/5bFT0DQ9x8xRmYmd+liUgYU7jPADEB49sfO4v42AB3vPA2vf2D/M9/WEAgoIAXkdEp3GeIQMD41qqFJMTG8MsXq+kbHOK2q85UwIvIqBTuM4iZccuH3kdCXAw/ea6K3v4hvnv1WcTG6NCJiLybwn2GMTP+beXpBGMD/ODpt+gdGOL2axYRH6uAF5F3KNxnqH+9ZB4JcTHc9vh29h7q4vZrzmZuXqrfZYlImNDu3gz2+fefwi8+vZj9rd186Mcv8uuXqhnSY/tEhDGEu5ndbWYNZrZlRFuWmT1tZru8eeaI9242syoz22lmK6eqcAlZuaCAJ2+6kAvm5nDrn7ex5lcV1Lf1+F2WiPhsLHvuvwYuP6rt68Czzrl5wLPea8xsPrAaWOBt8zMzi5m0amVUeakJ3LWmnNs+upDKPS2s/NEL/GXzAb/LEhEfnTDcnXMvAEc/5HMVcI+3fA9w1Yj2tc65XudcNVAFLJ2cUuV4zIxPLSvhsS+toDQnmS/e9zpfeWATbd39fpcmIj442TH3fOdcHYA3z/Pai4CR96mt9drew8xuMLNKM6tsbNTj5SbLKbkpPPRP5/GVS0/j0TcOcMWPXuDl3U1+lyUi02yyD6iOdkXNqEf4nHN3OOfKnXPlubm5k1xGdIuNCfDlS+fx0D+fTzAuhk/9cj23PbaNnv5Bv0sTkWlysuF+0MwKAbx5g9deCxSPWG82oMFfnywqzuCxL63g2mUl3LmumlU/eYltB9r9LktEpsHJhvujwBpveQ3wyIj21WYWNLMyYB5QMbESZSKS4mP51lUL+dV1SzjU1cdVP32JXzy/m0GdMikS0cZyKuT9wCvA6WZWa2bXA98GLjOzXcBl3mucc1uBB4FtwJPAjc45jQWEgQ+cnsdfb3o/F5+Rx388sYNP3vkq+w51+V2WiEwRC4f7g5eXl7vKykq/y4gKzjke2rifWx/dCsA3P7KAj51bpFsIi8xAZrbBOVc+2nu6QjXKmBlXL57NE1++kPmFaXzt929wzR2vUrnn6LNdRWQmU7hHqeKsJO6/YTnfumoh1U2dXP2fr/C5X7/G1gNtfpcmIpNAwzJCV98A97xcw38+v5u27n7+4exZfPWy0yjLSfa7NBE5juMNyyjcZVhbdz93vvA2d3kPA/lE+Wz+9eJ5zMpI9Ls0ERmFwl3GpbGjl58+V8V96/eCwaeXl/AvF51KdkrQ79JEZASFu5yU2pYu/u8zu3hoYy2JcTFcf+EpfP7CMlIT4vwuTURQuMsEVTV08MOn3+LxN+vJSIrjXy46lc+cV0pCnG74KeInhbtMijdr2/jeUzt54a1G8tOCfOmSeXyivJg4PcNVxBcKd5lUr77dzPf+upMNNS2UZCfxjyvKuOqcIg3XiEwzhbtMOuccz+1s4Pand/Hm/jaS42NYdU4R1y4rYf6sNL/LE4kKCneZMs45Nu1r5Xev7uUvmw/QOzDEuXMyuHZ5CVeeWahxeZEppHCXadHa1ccfNtRy7/q9VDd1kpkUx8fLi/nUsjmUZOuCKJHJpnCXaTU05Hh5dzO/e7WGp7cfZHDIceG8HK5dXsIlZ+QRqwOwIpNC4S6+qW/rYe1re1lbsY/69h4K0xNYvWQOq5cWk5+W4Hd5IjOawl18NzA4xDPbG7h3fQ3rdjURGzAum5/PtctLOO+UbAIB3XJYZLyOF+6x012MRKfYmACXLyzg8oUFVDd1ct/6Gn6/oZYnttRTkJbA5QsLuPLMQspLMhX0IpNAe+7im57+QZ7cUs9jb9bx/FuN9A0MkZca5PKFBVyxsJClZVnEKOhFjknDMhL2DvcO8LcdDTy+uY7ndjbQOzBETko8KxeE9uiXlWXpQKzIURTuMqN09g7w952NPL6ljr9tb6C7f5Cs5HhWLsjnioWFnHdqtm55IILCXWaw7r5Bnn+rgcffrOfZ7Qfp7BskPTGOD87P58qzCrng1BziYxX0Ep0U7hIRevoHeeGtRp7YUs8z2w7S0TtAakIs75+Xy4p5OayYm0NxVpLfZYpMG50tIxEhIS6GDy4o4IMLCugdGOSlqiaeeLOedbuaeOzNOgBKs5O8oM/lvFOzSU/UzcwkOmnPXWY85xxVDYdZt6uJF6uaePXtZrr6BgkYnF2cwYVzc1gxL5dz5mRorF4iioZlJKr0DQzx+t4WXqxqYt2uJjbXtjLkICUYy/JTsljhhf2pucmY6VRLmbkU7hLV2rr6eeXtpuE9+5rmLgAK0xNYMTeHpWVZLC7JpCxHYS8zi8JdZIR9h7pYt6uJdbsaeXl3M23d/QBkJsVx7pxMzi3JZHFJJmfPziAxXrcslvClcBc5hqEhx+7Gw2yoaWFDTQsb97awu7ETgJiAMb8wjcUlmZwzJ4PFJZkUZSRq717ChsJdZBxaOvt4fV8LG2ta2VDTwqZ9rXT3DwKQnxbk3DmZXuBnsrAojWCs9u7FH1N2KqSZ7QE6gEFgwDlXbmZZwANAKbAH+IRzrmUiP0dkOmUmx3PxGflcfEY+ELqj5Y76DjbubWFjTQsb9rbwxJZ6AOJjA5xRkMqCWWnMn5XOgllpvK8gTcM54rsJ7bl74V7unGsa0fZd4JBz7ttm9nUg0zn33473Odpzl5mmoaOHjTWtbNzbwpb9bWw90D48dh8wOCU3JRT4hWks8EI/Mzne56ol0kzZsMwxwn0ncJFzrs7MCoG/O+dOP97nKNxlpnPOsb+1m60H2tl2oN2bt3GgrWd4nVnpCcN796E9/TSN4cuETOUVqg54yswc8Avn3B1AvnOuDsAL+LwJ/gyRsGdmzM5MYnZmEisXFAy3H+rs88I+tHe/9UAbz+44yJF9qoykOE7LT2VeXkpo8pZzU4MKfZmQiYb7Bc65A16AP21mO8a6oZndANwAMGfOnAmWIRKespLjQ7dDmJcz3NbVN8CO+o7hvftdBw/z5zcO0N4zMLxOWkLscNDPHRH6hekJCn0Zk0k7W8bMbgUOA59HwzIi4+Kco/FwL1UHD7Or4TC7GjrYdfAwVQ2Hae7sG14vOT6GuXkpzM1LZV5+aG//lNwUZmcm6tYKUWhKhmXMLBkIOOc6vOUPAv8LeBRYA3zbmz9ysj9DJFqYGXmpCeSlJnD+3Jx3vdd8uJeqhlDoV3nBv25XIw9trB1eJyZgFGUkUpqTTGl2EqXZyZTmJFGSnUxxZpJuixyFJjIskw/80fsTMRa4zzn3pJm9BjxoZtcDe4GPT7xMkeiVnRIkOyXIslOy39Xe1tVPVWMH1U1d7GnqZE9zJzXNXbxe00JH7ztDPAGDoszEUOBnJ1MyHP7JFGcl6jz9CKWLmEQijHOOQ5197GnuZE9TFzXNnVQ3e/OmTjp63h38BWkJ3sHgRIoyE5mdmTj8ujA9UXv9YUz3cxeJImY2vLe/uCTrXe8552jp6vf28jupbuqitqWL2pZu1lcfom5TN0Nu5GcdCf9Q4BdlHBX+GQna8w9TCneRKGJmZCXHk5Ucz7lzMt/zfv/gEPVtPezzAr+2pZv9Ld3UtnRRUX2Iurb3hn9OSpBZ6QkUpCdQmJ7ozRMoSEtgVkYieWlBfQH4QOEuIsPiYgIUZyUd83GFR8K/1gv82pZu6tt6ONDWzduNnbxc1fyu8f4jclLiKUhPoCAtkcL0BAozjnwBJJKfFiQvLYGUoOJoMum3KSJj9u7wzx51nY6efg6291DX5k2tPdS3d1PX1kNtSxev7Tk0fKuGkZLiY8hLDZKXmkBuWnB4OS81SF7aO8sZSXE6138MFO4iMqlSE+JITYhjbl7qMdfp6hug3gv/g+09NHT00tDeS0NHaHnbgXb+3t5DZ9/ge7aNjwmQOxz4QXJSjkzxZHvL2Snx5KQESUuIjdovAoW7iEy7pPhYTskNXYB1PJ29A17we18AHaEvgMb20HJ1Uyev7WmhpauP0U78i4sxspOD5KTGk50cCv1cL/xD7UGyk+PJTI4nKyk+ou7mqXAXkbCVHIylLBhLWU7ycdcbGBziUFcfzYdDU9PhXm/qo9lbbu7so6rhMI2He+kbGBr1cxLiAmQnB8lMjiMzKXTg+cj8yDTc7q0TrlcGK9xFZMaLjQkMX+F7Is45DvcODH8JNHf20dLZx6Eub97ZT0tXH82dfdQ0d9HS2TfqQeIjUhNiyUiKIyMxPjRPiicjMW705aQ4MhLjSE+MI3aKvxQU7iISVcxs+LhA6Qn+Ijiib2CI1q7QF8ChztDU0hn6Amjt6qetO/SF0NrVT21LN61dfbR197/rtNGjHflSuHxBAbd8aP4k9e4dCncRkROIjw2Ql5ZAXtqJ/zI4YmjI0dEzQGt3Hy1d/cOB39LZR2t3//CXQkF64pTUrHAXEZkCgYCRnhRHelIcJaOfNTq1P3/6f6SIiEw1hbuISARSuIuIRCCFu4hIBFK4i4hEIIW7iEgEUriLiEQghbuISAQKi2eomlkjUON3HT7KAZr8LsJH6r/6r/6fnBLnXO5ob4RFuEc7M6s81kNuo4H6r/6r/5Pffw3LiIhEIIW7iEgEUriHhzv8LsBn6n90U/+ngMbcRUQikPbcRUQikMJdRCQCKdynmJkVm9lzZrbdzLaa2Ze99iwze9rMdnnzzBHb3GxmVWa208xW+lf95DGzGDN73cz+4r2Omv6bWYaZ/cHMdnj/Ds6Lsv5/xfu3v8XM7jezhEjuv5ndbWYNZrZlRNu4+2tmi83sTe+9H5uZjasQ55ymKZyAQuBcbzkVeAuYD3wX+LrX/nXgO97yfOANIAiUAbuBGL/7MQm/h68C9wF/8V5HTf+Be4B/9JbjgYxo6T9QBFQDid7rB4HPRnL/gfcD5wJbRrSNu79ABXAeYMATwBXjqUN77lPMOVfnnNvoLXcA2wn9g19F6D893vwqb3kVsNY51+ucqwaqgKXTWvQkM7PZwIeAX45ojor+m1kaof/sdwE45/qcc61ESf89sUCimcUCScABIrj/zrkXgENHNY+rv2ZWCKQ5515xoaT/zYhtxkThPo3MrBQ4B1gP5Dvn6iD0BQDkeasVAftGbFbrtc1kPwL+HRga0RYt/T8FaAR+5Q1L/dLMkomS/jvn9gPfB/YCdUCbc+4poqT/I4y3v0Xe8tHtY6ZwnyZmlgI8BNzknGs/3qqjtM3Y81XN7MNAg3Nuw1g3GaVtxvaf0F7rucDPnXPnAJ2E/iw/lojqvze2vIrQkMMsINnMrj3eJqO0zdj+j8Gx+jvh34PCfRqYWRyhYL/XOfew13zQ+9MLb97gtdcCxSM2n03oz9iZ6gLgI2a2B1gLXGxmvyN6+l8L1Drn1nuv/0Ao7KOl/5cC1c65RudcP/AwcD7R0/8jxtvfWm/56PYxU7hPMe8I913AdufcD0e89SiwxlteAzwyon21mQXNrAyYR+jAyozknLvZOTfbOVcKrAb+5py7lujpfz2wz8xO95ouAbYRJf0nNByz3MySvP8LlxA67hQt/T9iXP31hm46zGy593v7zIhtxsbvI8uRPgErCP05tRnY5E1XAtnAs8Aub541YptbCB0138k4j5CH8wRcxDtny0RN/4FFQKX3b+BPQGaU9f+bwA5gC/BbQmeGRGz/gfsJHV/oJ7QHfv3J9Bco935nu4Gf4N1RYKyTbj8gIhKBNCwjIhKBFO4iIhFI4S4iEoEU7iIiEUjhLiISgRTuIiIRSOEuIhKB/j/gQm2IsntfPgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "As = 2.105  # 10^9 A_s\n",
    "h = 0.6766\n",
    "Om = 0.3111\n",
    "Ob = 0.02242 / h ** 2\n",
    "ns = 0.9665\n",
    "\n",
    "As = 1.7988  # 10^9 A_s\n",
    "h = 0.68\n",
    "Om = 0.335\n",
    "Ob = 0.049\n",
    "ns = 0.97\n",
    "\n",
    "all_R = np.linspace(50, 1000, 20)\n",
    "all_sigma = np.empty(len(all_R))\n",
    "\n",
    "for i in tqdm(range(len(all_sigma))):\n",
    "    all_sigma[i] = get_sig(all_R[i], As, Om, Ob, h, ns)\n",
    "plt.plot(all_R, all_sigma)\n",
    "\n",
    "print(get_sig(500, As, Om, Ob, h, ns))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9581f953-0bb4-4c71-a1ba-32ade3a95e27",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "symreg",
   "language": "python",
   "name": "symreg"
  },
  "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}