Successful compilation + HyREC FULL mode

* Successful compilation of the python wrapper based on settings in Makefile and setup.py. Do not mess with these anymore! * Turned on HyREC FULL mode so that we're no longer solving using the RECFAST mode, which was dumb.
hongwanliu · Sep 25, 2019 · 035830c · 035830c
1 parent dc9cb02
commit 035830c
Show file tree

Hide file tree

Showing 7 changed files with 376 additions and 73 deletions.
diff --git a/Makefile b/Makefile
@@ -16,9 +16,13 @@ vpath .base build
 ###### LINES TO ADAPT TO YOUR PLATEFORM ################
 ########################################################
 
+# Successfully compiled with gcc 9.2.0, in the python3 environment,
+# with all the openmp lgomp flags turned on. Note that the GCCPATH
+# seems to have to be specified manually in python/setup.py.
+
 # your C compiler:
 # CC       = gcc
-CC        = /usr/local/Cellar/gcc/9.1.0/bin/gcc-9
+CC        = /usr/local/Cellar/gcc/9.2.0/bin/gcc-9
 #CC       = icc
 #CC       = pgcc
 

diff --git a/hyrec/history.h b/hyrec/history.h
@@ -17,7 +17,7 @@
 #define FULL      3    /* All radiative transfer effects included. Additional switches in header file hydrogen.h */
 
 /** here is the switch **/
-#define MODEL RECFAST     /* default setting: FULL */
+#define MODEL FULL     /* default setting: FULL */
 
 /***** Switches for derivative d(xe)/dt *****/
 

diff --git a/hyrec/hydrogen.h b/hyrec/hydrogen.h
@@ -39,12 +39,8 @@ double rec_HRecFast_dxedlna(double xe, double nH, double H, double TM, double TR
 
 /************* EFFECTIVE MULTI LEVEL ATOM *******************/
 
-// #define ALPHA_FILE  "Alpha_inf.dat"     /* Contains the effective recombination coefficients to 2s and 2p */
-// #define RR_FILE     "R_inf.dat"         /* Contains the effective transfer rate R_{2p,2s} */
-
-// HL: Changed to variables in the hope that they can be accessed. 
-char * ALPHA_FILE = "Alpha_inf.dat"
-char * RR_FILE    = "R_inf.dat"
+#define ALPHA_FILE  "Alpha_inf.dat"     /* Contains the effective recombination coefficients to 2s and 2p */
+#define RR_FILE     "R_inf.dat"         /* Contains the effective transfer rate R_{2p,2s} */
 
 
 /* Boundaries and number of elements of temperature tables */

diff --git a/notebooks/.ipynb_checkpoints/recfast_fudge-checkpoint.ipynb b/notebooks/.ipynb_checkpoints/recfast_fudge-checkpoint.ipynb
@@ -1,6 +1,276 @@
 {
- "cells": [],
- "metadata": {},
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%load_ext autoreload"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# import classy module\n",
+    "%autoreload 2\n",
+    "from classy import Class\n",
+    "\n",
+    "%matplotlib inline\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "did it set correctly???????\n",
+      "b'./hyrec/fake_Alpha.dat'\n"
+     ]
+    },
+    {
+     "ename": "CosmoComputationError",
+     "evalue": "\n\nError in Class: thermodynamics_init(L:388) :error in thermodynamics_recombination(ppr,pba,pth,preco,pvecback);\n=>thermodynamics_recombination(L:2601) :error in thermodynamics_recombination_with_hyrec(ppr,pba,pth,preco,pvecback);\n=>thermodynamics_recombination_with_hyrec(L:2736) :could not open fA with name ./hyrec/fake_Alpha.dat and mode \"r\"",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mCosmoComputationError\u001b[0m                     Traceback (most recent call last)",
+      "\u001b[0;32m<ipython-input-2-bd4272b0fed1>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     10\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0;31m# run class\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 12\u001b[0;31m \u001b[0mLambdaCDM_new\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcompute\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
+      "\u001b[0;32mclassy.pyx\u001b[0m in \u001b[0;36mclassy.Class.compute\u001b[0;34m()\u001b[0m\n",
+      "\u001b[0;31mCosmoComputationError\u001b[0m: \n\nError in Class: thermodynamics_init(L:388) :error in thermodynamics_recombination(ppr,pba,pth,preco,pvecback);\n=>thermodynamics_recombination(L:2601) :error in thermodynamics_recombination_with_hyrec(ppr,pba,pth,preco,pvecback);\n=>thermodynamics_recombination_with_hyrec(L:2736) :could not open fA with name ./hyrec/fake_Alpha.dat and mode \"r\""
+     ]
+    }
+   ],
+   "source": [
+    "from classy import Class\n",
+    "\n",
+    "LambdaCDM_new = Class()\n",
+    "# pass input parameters\n",
+    "LambdaCDM_new.set({\n",
+    "    'recombination': 'hyrec', \n",
+    "    'Alpha_inf hyrec file': './hyrec/fake_Alpha.dat'\n",
+    "})\n",
+    "LambdaCDM_new.set({'output':'tCl,pCl,lCl,mPk','lensing':'yes','P_k_max_1/Mpc':3.0})\n",
+    "\n",
+    "# run class\n",
+    "LambdaCDM_new.compute()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "did it set correctly???????\n",
+      "b'../hyrec/Alpha_inf_orig.dat'\n",
+      "did it set correctly???????\n",
+      "b'../hyrec/Alpha_inf.dat'\n"
+     ]
+    }
+   ],
+   "source": [
+    "# create instance of the class \"Class\"\n",
+    "from classy import Class\n",
+    "\n",
+    "LambdaCDM = Class()\n",
+    "hyrec_mod = Class()\n",
+    "# pass input parameters\n",
+    "\n",
+    "hyrec_mod.set({\n",
+    "    'recombination': 'hyrec',\n",
+    "    'Alpha_inf hyrec file': '../hyrec/Alpha_inf.dat',\n",
+    "    'output':'tCl,pCl,lCl,mPk',\n",
+    "    'lensing':'yes','P_k_max_1/Mpc':3.0\n",
+    "})\n",
+    "LambdaCDM.set({\n",
+    "    'recombination': 'hyrec',\n",
+    "    'Alpha_inf hyrec file': './hyrec/Alpha_BB_n_250_fine.dat',\n",
+    "    'output':'tCl,pCl,lCl,mPk',\n",
+    "    'lensing':'yes','P_k_max_1/Mpc':3.0\n",
+    "})\n",
+    "\n",
+    "# run class\n",
+    "LambdaCDM.compute()\n",
+    "hyrec_mod.compute()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/Users/hongwan/anaconda/envs/python3/lib/python3.6/site-packages/ipykernel_launcher.py:13: RuntimeWarning: invalid value encountered in true_divide\n",
+      "  del sys.path[0]\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "Text(0, 0.5, '$(C_\\\\ell - C_\\\\ell^\\\\mathrm{ref})/C_\\\\ell^\\\\mathrm{ref} [\\\\%]$')"
+      ]
+     },
+     "execution_count": 10,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEMCAYAAAAS+xsDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFNZJREFUeJzt3X+MZWWd5/H3x25AmTH8bKQF2m6kJy4yCWKl0YizZkagMdHGVTPN7GDvLG7HDST+mslCDIFRJ+rssuwadTY9oxGJiviDoWfHWeSXa8YsSLWiggjdMs7SQgCBYYZFQOC7f9xTWpRV3aeq79O36vb7ldzce57zVN3vU+emP33Oee45qSokSWrpeaMuQJI0/gwbSVJzho0kqTnDRpLUnGEjSWrOsJEkNWfYSJKaM2wkSc0ZNpKk5paPuoDF4vDDD6/Vq1ePugxJWjK2bdv2s6pa0aevYdNZvXo1k5OToy5DkpaMJP/Yt6+H0SRJzRk2kqTmDBtJUnOGjSSpOcNGktScYSNJas6wkSQ1Z9hIkpozbCRJzRk2kqTmDBtJUnOGjSSpOcNGktScYSNJas6wkSQ1Z9hIkpozbCRJzRk2kqTmDBtJUnOGjSSpOcNGktScYSNJas6wkSQ1Z9hIkppbtGGTZH2SO5PsSHL+LOsPSPLFbv3NSVbPWL8qyWNJ/nhv1SxJmt2iDJsky4BPAGcAxwNnJTl+RrdzgEeq6jjgUuCjM9ZfCvxd61olSbu3KMMGWAfsqKq7q+op4Apgw4w+G4DLutdfBn4vSQCSnAncDdy+l+qVJO3CYg2bo4B7pi3v7Npm7VNVTwOPAocl+Q3gPwF/uhfqlCT1sFjDJrO0Vc8+fwpcWlWP7fZNks1JJpNMPvjggwsoU5LUx/JRFzCHncAx05aPBu6do8/OJMuBg4CHgZOBtyb5c+Bg4NkkT1TVx2e+SVVtAbYATExMzAwzSdKQLNawuQVYm2QN8FNgI/AHM/psBTYB/wd4K3BDVRXw2qkOSS4GHpstaCRJe8+iDJuqejrJecA1wDLg01V1e5IPAJNVtRX4FHB5kh0M9mg2jq5iSdKuZLAzoImJiZqcnBx1GZK0ZCTZVlUTffou1gkCkqQxYthIkpozbCRJzRk2kqTmDBtJUnOGjSSpOcNGktScYSNJas6wkSQ1Z9hIkpozbCRJzRk2kqTmDBtJUnOGjSSpOcNGktScYSNJas6wkSQ1Z9hIkpozbCRJzRk2kqTmDBtJUnOGjSSpOcNGktScYSNJas6wkSQ1Z9hIkpozbCRJzRk2kqTmDBtJUnOGjSSpOcNGktTcog2bJOuT3JlkR5LzZ1l/QJIvdutvTrK6az81ybYkP+ief3dv1y5Jeq7lfTolObRHt2er6p/2sJ6p91sGfAI4FdgJ3JJka1X9cFq3c4BHquq4JBuBjwK/D/wMeGNV3ZvkBOAa4Khh1CVJWpheYQPc2z2yiz7LgFV7XNHAOmBHVd0NkOQKYAMwPWw2ABd3r78MfDxJquq70/rcDjw/yQFV9eSQapMkzVPfsLmjql6xqw5Jvrur9fN0FHDPtOWdwMlz9amqp5M8ChzGYM9myluA7xo0kjRafcPm1UPq09dse1A1nz5JXs7g0Nppc75JshnYDLBq1bB2yiRJM/WaIFBVTwyjzzzsBI6Ztnw0g8N4s/ZJshw4CHi4Wz4auAp4e1X9eBc1b6mqiaqaWLFixRDLlyRNt6DZaEleleSGJN9KcuawiwJuAdYmWZNkf2AjsHVGn63Apu71W4EbqqqSHAz8LXBBVX2rQW2SpHnqFTZJjpzR9F7gTcB64IPDLqqqngbOYzCT7A7gyqq6PckHkryp6/Yp4LAkO7p6pqZHnwccB1yY5NbuccSwa5Qk9ZeqmadCZumU/DWwDfjPVfVEki3AJPAs8EdV9Zq2ZbY3MTFRk5OToy5DkpaMJNuqaqJP377nbM4EbgX+Z5KzgXczCJoDgRaH0SRJY6T3OZuq+hvgdOBg4KvAnVX1sap6sFVxkqTx0PeczZuS/D1wA3AbgxP2b07yhSQvbVmgJGnp6/s9mw8x+B7NC4CvVdU64L1J1gJ/xiB8JEmaVd+weZRBoLwAeGCqsaq2Y9BIknaj7zmbNzOYDPA08AftypEkjaO+ezZfr6qTdtUhyXd210eStG/qGzb/Ksn3d7E+DC4XI0nSr+kbNi/r0eeZPSlEkjS+dhs2Sd4C/LSqbtoL9UiSxlCfPZu3A8uSfG2qoao+2a4kSdK46RM27wH+I4O7XkqSNG+7DZvu1sx/shdqkSSNqb6Xq3lXkr/sXl/YtiRJ0rjp+6XOlwL3dK9f2KgWSdKY6hs2BbwgyQnAixvWI0kaQ7sNmyQBHmLwxc2zgQtaFyVJGi+7DZsa3MpzLfA94JvAbwMkOSXJuUmOneqbZE2rQiVJS1ffKwhcB+wPHM7gkBrACmAdsC7JQ8DngfcBZw27SEnS0tYrbKrqslnarkqyFXglcALwWuCu4ZYnSRoHvcKmm+78eFVdMr29qp4Bvt09JEmaVd/DaGcDJ85sTPIOYEVVfXioVUmSxkrfqc8/r6rHZ2m/HPjDIdYjSRpDvcMmycqZjVX1JIO7d0qSNKe+YXMJcHWSl0xvTHIE8OzQq5IkjZW+s9G+lORAYFuSm4BbGQTV24CL25UnSRoHffdspqY/HwtcCewHPAGcVVWfa1SbJGlM9J36/Grgpqr6Z+CzbUuSJI2bvns2mxgcQrsiyb9LcmTLoiRJ46XvOZt3AiR5GXAG8JkkBwE3Av8L+Fb3BU9Jkn5N73M2AFX1o6q6FHgLcDrw9wwmCdw87MKSrE9yZ5IdSc6fZf0BSb7Yrb85yepp6y7o2u9Mcvqwa5MkzU/fczbPAzYC/xaYAJ4CDgAeBL7GkC++mWQZ8AngVGAncEuSrVX1w2ndzgEeqarjkmwEPgr8fpLju1pfzuDeO9cl+S33vCRpdPru2dzI4G6dFwArq+qYqjqCwcU3bwI+nGSYVxJYB+yoqrur6ingCmDDjD4bgKkLhH4Z+L3u3jsbgCuq6smq+gdgR/f7JEkj0vfaaK+vql/MbKyqh4GvAF9Jst8Q6zqKX92GGgZ7NyfP1aeqnk7yKHBY137TjJ89aoi1/dJTTz/Lu674botfLUl7xZEHPZ+L3vjy5u+z27BJ8hbgpzz3H/BfM1sY7YHM9hY9+/T52cEvSDYDmwFWrVo1n/q6X1r8+MHH5v1zkrRY/OKZWf95HLo+ezZvB5Yl+dpUQ1V9sl1JwGBv5Jhpy0cD987RZ2eS5cBBwMM9fxaAqtoCbAGYmJiY91/8gOXL+Pp7/vV8f0yS9jl9ztm8B7gDuH3ao7VbgLVJ1iTZn8EJ/60z+mxl8P0fgLcCN3S3sN4KbOxmq61hcEtr77cjSSO02z2bqrob+JO9UMv093w6yXnANcAy4NNVdXuSDwCTVbUV+BRweZIdDPZoNnY/e3uSK4EfMrgi9bnORJOk0cpgZ2A3nZJ3ASdU1X9IcmFVfbB9aXvXxMRETU5OjroMSVoykmyrqok+fftOfX4pv5od9sIFVSVJ2mf1DZsCXpDkBAZflJQkqbfdhk33RcmHGEwpPpvBFzslSeptt2HTzfBaC3wP+Cbw2wBJTklybpJjp/p2s78kSXqOvlcQuA7YHzicX31BcgWDy8CsS/IQ8HngfQz5OmmSpKWv7y0GLpul7aokW4FXAicwuE7aXcMtT5I0Dvpe9flC4PGqumR6e/f9lW/jlyYlSbvQ9zDa2cCJMxuTvANYUVUfHmpVkqSx0nfq88+r6vFZ2i8HhnlrAUnSGOodNklWzmysqicZXBJGkqQ59Q2bS4Crk7xkemOSI4Bnh16VJGms9J2N9qUkBwLbktwE3MogqN4GXNyuPEnSOOi7ZzM1/XkNcCWwH/AEcFZVfa5RbZKkMdF3NhoAVfUvwGcb1SJJGlO992wkSVqoBYVNkjcOuxBJ0vha6J7Nnw21CknSWFto2GSoVUiSxtpCw2b395KWJKnjBAFJUnOGjSSpuYWGzf1DrUKSNNYWFDZVdeqwC5EkjS8Po0mSmjNsJEnNGTaSpOYWermaU5Kcm+TYaW1rhleWJGmczOuqz9OsANYB65I8BHweeB9w1rAKkySNjwWFTVVdlWQr8ErgBOC1wF3DLEySND56hU2SC4HHq+qSqbaqegb4dveQJGlOffdszgZOnNmY5B3Aiqr68FCrkiSNlb4TBH5eVY/P0n458IdDrIckhya5Nsn27vmQOfpt6vpsT7Kpazswyd8m+VGS25N8ZJi1SZIWpnfYJFk5s7GqngSeHm5JnA9cX1Vrgeu75edIcihwEXAyg4kKF00Lpf9SVS8DXgG8JskZQ65PkjRPfcPmEuDqJC+Z3pjkCODZIde0Abise30ZcOYsfU4Hrq2qh6vqEeBaYH1VPV5VNwJU1VPAd4Cjh1yfJGmeep2zqaovJTkQ2JbkJuBWBkH1NuDiIdf0oqq6r3vf+7pAm+ko4J5pyzu7tl9KcjDwRuC/D7k+SdI89Z2Nlqq6LMlXgTcDLwf+H3BWVU1O69PrpmpJrgOOnGXV+/uVPeudQn/53kmWA18APlZVd++ijs3AZoBVq1b1fGtJ0nz1nY12Y5KvAFdX1WenGpPsn+R3gU3AjcBn+vyyqnr9XOuS3J9kZbdXsxJ4YJZuO4HXTVs+GvjGtOUtwPaq+m+7qWNL15eJiQnvPipJjfQ9Z7MeeAb4QpJ7k/wwyd3AdgZXDbi0qj4zpJq2MggvuuerZ+lzDXBakkO6iQGndW0k+RBwEPDuIdUjSdpDfc/ZPAF8Evhkkv2AwxlMh/6nBjV9BLgyyTnA/2VwXogkE8A7q+odVfVwkg8Ct3Q/84Gu7WgGh+J+BHwnCcDHq+qvGtQpSeopPU+zjL2JiYmanJwcdRmStGQk2VZVE336eosBSVJzho0kqTnDRpLUnGEjSWrOsJEkNWfYSJKaM2wkSc0ZNpKk5gwbSVJzho0kqTnDRpLUnGEjSWrOsJEkNWfYSJKaM2wkSc0ZNpKk5gwbSVJzho0kqTnDRpLUnGEjSWrOsJEkNWfYSJKaM2wkSc0ZNpKk5gwbSVJzho0kqTnDRpLUnGEjSWrOsJEkNWfYSJKaW3Rhk+TQJNcm2d49HzJHv01dn+1JNs2yfmuS29pXLEnanUUXNsD5wPVVtRa4vlt+jiSHAhcBJwPrgIumh1KSfwM8tnfKlSTtzmIMmw3AZd3ry4AzZ+lzOnBtVT1cVY8A1wLrAZL8JvBe4EN7oVZJUg+LMWxeVFX3AXTPR8zS5yjgnmnLO7s2gA8ClwCPtyxSktTf8lG8aZLrgCNnWfX+vr9ilrZKciJwXFW9J8nqHnVsBjYDrFq1qudbS5LmayRhU1Wvn2tdkvuTrKyq+5KsBB6YpdtO4HXTlo8GvgG8Gnhlkp8wGNsRSb5RVa9jFlW1BdgCMDExUfMfiSSpj8V4GG0rMDW7bBNw9Sx9rgFOS3JINzHgNOCaqvqLqnpxVa0GTgHumitoJEl7z2IMm48ApybZDpzaLZNkIslfAVTVwwzOzdzSPT7QtUmSFqFUefQIBofRJicnR12GJC0ZSbZV1USfvotxz0aSNGYMG0lSc4aNJKk5w0aS1JxhI0lqzrCRJDVn2EiSmjNsJEnNGTaSpOYMG0lSc4aNJKk5w0aS1JxhI0lqzrCRJDVn2EiSmjNsJEnNGTaSpOYMG0lSc4aNJKk5w0aS1JxhI0lqzrCRJDVn2EiSmjNsJEnNGTaSpOZSVaOuYVFI8iDwjwv40cOBnw25nMXOMe8bHPO+YU/G/JKqWtGno2Gzh5JMVtXEqOvYmxzzvsEx7xv21pg9jCZJas6wkSQ1Z9jsuS2jLmAEHPO+wTHvG/bKmD1nI0lqzj0bSVJzhs0eSLI+yZ1JdiQ5f9T1DFOSnyT5QZJbk0x2bYcmuTbJ9u75kK49ST7W/R2+n+Sk0VbfT5JPJ3kgyW3T2uY9xiSbuv7bk2waxVj6mmPMFyf5abetb03yhmnrLujGfGeS06e1L4nPfpJjktyY5I4ktyd5V9c+ttt5F2Me7XauKh8LeADLgB8DxwL7A98Djh91XUMc30+Aw2e0/Tlwfvf6fOCj3es3AH8HBHgVcPOo6+85xt8BTgJuW+gYgUOBu7vnQ7rXh4x6bPMc88XAH8/S9/juc30AsKb7vC9bSp99YCVwUvf6hcBd3bjGdjvvYswj3c7u2SzcOmBHVd1dVU8BVwAbRlxTaxuAy7rXlwFnTmv/bA3cBBycZOUoCpyPqvom8PCM5vmO8XTg2qp6uKoeAa4F1revfmHmGPNcNgBXVNWTVfUPwA4Gn/sl89mvqvuq6jvd638B7gCOYoy38y7GPJe9sp0Nm4U7Crhn2vJOdr1Bl5oCvp5kW5LNXduLquo+GHyggSO69nH6W8x3jOMy9vO6w0afnjqkxJiNOclq4BXAzewj23nGmGGE29mwWbjM0jZOU/teU1UnAWcA5yb5nV30Hfe/Bcw9xnEY+18ALwVOBO4DLunax2bMSX4T+Arw7qr65111naVtXMY80u1s2CzcTuCYactHA/eOqJahq6p7u+cHgKsY7FLfP3V4rHt+oOs+Tn+L+Y5xyY+9qu6vqmeq6lngLxlsaxiTMSfZj8E/up+rqq92zWO9nWcb86i3s2GzcLcAa5OsSbI/sBHYOuKahiLJbyR54dRr4DTgNgbjm5qFswm4unu9FXh7N5PnVcCjU4colqD5jvEa4LQkh3SHJU7r2paMGefX3sxgW8NgzBuTHJBkDbAW+DZL6LOfJMCngDuq6r9OWzW223muMY98O4965sRSfjCYuXIXgxkb7x91PUMc17EMZp58D7h9amzAYcD1wPbu+dCuPcAnur/DD4CJUY+h5zi/wOBwwi8Y/C/unIWMEfj3DE6q7gD+aNTjWsCYL+/G9P3uH5OV0/q/vxvzncAZ09qXxGcfOIXBoZ/vA7d2jzeM83bexZhHup29goAkqTkPo0mSmjNsJEnNGTaSpOYMG0lSc4aNJKk5w0aS1JxhI0lqzrCRFrHuag4f777NLi1Zho20uL2TwX1GThl1IdKeMGykxW09g8uF3DrqQqQ9YdhIi1SS5zO4W+JJwP8ecTnSHjFspMVrLYOw+VFV/WLUxUh7YvmoC5A0pxXAb7FIb7kszYd7NtLi9WIGN8B63rRb+EpLkmEjLUJJljM4V3Mk8D+AZ0ZbkbRnvJ+NJKk592wkSc0ZNpKk5gwbSVJzho0kqTnDRpLUnGEjSWrOsJEkNWfYSJKa+//6Mq33aMve2QAAAABJRU5ErkJggg==\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "cls_new = hyrec_mod.lensed_cl(2500)\n",
+    "cls_std = LambdaCDM.lensed_cl(2500)\n",
+    "\n",
+    "ll_vec_new = cls_new['ell']\n",
+    "ll_vec_std = cls_std['ell']\n",
+    "\n",
+    "plt.figure()\n",
+    "\n",
+    "ax = plt.gca()\n",
+    "\n",
+    "# ax.set_yscale('log')\n",
+    "\n",
+    "plt.plot(ll_vec_new, (cls_new['tt']/cls_std['tt'] - 1.)*100)\n",
+    "# plt.plot(ll_vec_new, (cls_new['ee']/cls_std['ee'] - 1.)*100)\n",
+    "\n",
+    "\n",
+    "# plt.axis([0, 2500, -1e-10, 1e-10])\n",
+    "\n",
+    "plt.xlabel(r'$\\ell$')\n",
+    "plt.ylabel(r'$(C_\\ell - C_\\ell^\\mathrm{ref})/C_\\ell^\\mathrm{ref} [\\%]$')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "np.max(\n",
+    "    new.get_thermodynamics()['x_e'] \n",
+    "    - LambdaCDM.get_thermodynamics()['x_e']\n",
+    ")\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%autoreload\n",
+    "# LambdaCDM_hyrec = Class()\n",
+    "# LambdaCDM_hyrec.set({'recombination': 'hyrec', 'output':'tCl,pCl,lCl,mPk','lensing':'yes','P_k_max_1/Mpc':3.0})\n",
+    "LambdaCDM_newAlpha = Class()\n",
+    "LambdaCDM_newAlpha.set({\n",
+    "    'recombination': 'hyrec', \n",
+    "    'hyrec_Alpha_inf_file': './hyrec/Alpha_inf.dat',\n",
+    "    'output':'tCl,pCl,lCl,mPk',\n",
+    "    'lensing':'yes',\n",
+    "    'P_k_max_1/Mpc':3.0\n",
+    "})\n",
+    "# LambdaCDM_hyrec.compute()\n",
+    "LambdaCDM_newAlpha.compute()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cls_hyrec=LambdaCDM_hyrec.lensed_cl(2500)\n",
+    "cls_newAlpha = LambdaCDM_newAlpha.lensed_cl(2500)\n",
+    "\n",
+    "ll_vec_hyrec = cls_hyrec['ell']\n",
+    "\n",
+    "plt.figure()\n",
+    "\n",
+    "ax = plt.gca()\n",
+    "\n",
+    "# ax.set_yscale('log')\n",
+    "\n",
+    "plt.plot(ll_vec_hyrec, (cls_newAlpha['tt']/cls_hyrec['tt'] - 1.)*100)\n",
+    "plt.plot(ll_vec_hyrec, (cls_newAlpha['ee']/cls_hyrec['ee'] - 1.)*100)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "a = np.loadtxt('../hyrec/Alpha_inf.dat')\n",
+    "a = a/1e5\n",
+    "np.savetxt('../hyrec/Alpha_inf.dat', a)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "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.6.9"
+  },
+  "toc": {
+   "base_numbering": 1,
+   "nav_menu": {},
+   "number_sections": true,
+   "sideBar": true,
+   "skip_h1_title": false,
+   "title_cell": "Table of Contents",
+   "title_sidebar": "Contents",
+   "toc_cell": false,
+   "toc_position": {},
+   "toc_section_display": true,
+   "toc_window_display": false
+  }
+ },
  "nbformat": 4,
  "nbformat_minor": 2
 }