FisherIM-MK-omHIbHI.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Author: Alkistis Pourtsidou, ICG Portsmouth\n",
    "Using part of http://camb.readthedocs.io/en/latest/CAMBdemo.html \n",
    "\n",
    "To run this notebook you need to have CAMB and the CAMB python package \n",
    "installed. In order to install the CAMB python package on your computer follow\n",
    "the instructions in http://camb.readthedocs.io/en/latest/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import sys, platform, os\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import camb\n",
    "from camb import model, initialpower"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pylab as pl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "font = {'size'   : 16, 'family':'STIXGeneral'}\n",
    "axislabelfontsize='x-large'\n",
    "plt.rc('font', **font)\n",
    "plt.rcParams['legend.fontsize']='medium'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import scipy\n",
    "from scipy.interpolate import interp1d\n",
    "from __future__ import division"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy import integrate\n",
    "from scipy import linalg\n",
    "\n",
    "pi=np.pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Fiducial cosmological parameters\n",
    "c=3e5\n",
    "hubble=0.678\n",
    "omegab=0.022*pow(hubble,-2)\n",
    "omegac=0.119*pow(hubble,-2)\n",
    "om0=omegac+omegab\n",
    "H00=100*hubble\n",
    "Ass=2.14e-9\n",
    "nss = 0.968\n",
    "\n",
    "gamma=0.545\n",
    "\n",
    "print om0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Set up the fiducial cosmology\n",
    "pars = camb.CAMBparams()\n",
    "#Set cosmology\n",
    "pars.set_cosmology(H0=H00, ombh2=omegab*pow(hubble,2), omch2=omegac*pow(hubble,2),omk=0,mnu=0)\n",
    "pars.set_dark_energy() #LCDM (default)\n",
    "pars.InitPower.set_params(ns=nss, r=0, As=Ass)\n",
    "pars.set_for_lmax(2500, lens_potential_accuracy=0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#calculate results for these parameters\n",
    "results = camb.get_results(pars)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Get matter power spectrum at z=0: P(k,z=0)\n",
    "\n",
    "#Not non-linear corrections couples to smaller scales than you want\n",
    "pars.set_matter_power(redshifts=[0.], kmax=2.0)\n",
    "\n",
    "#Non-Linear spectra (Halofit)\n",
    "pars.NonLinear = model.NonLinear_both\n",
    "results.calc_power_spectra(pars)\n",
    "kh, z, pk = results.get_matter_power_spectrum(minkh=1e-4, maxkh=2.0, npoints = 200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Construct P(k,z=0) interpolating function, in units of Mpc (no h)\n",
    "Pkz0 = interp1d(kh*hubble, pk[0]/pow(hubble,3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Redshift bins\n",
    "zlist = np.arange(0.1,1.45,0.1)\n",
    "ztest = zlist[4]\n",
    "Nzbins = len(zlist)\n",
    "\n",
    "print zlist\n",
    "print \"ztest =\", ztest\n",
    "print \"Number of redshift bins =\", Nzbins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Define E(z) = H(z)/H0\n",
    "def Ez(zc):\n",
    "    return np.sqrt(1-om0+om0*pow(1+zc,3))\n",
    "\n",
    "#Define the comoving distance\n",
    "def drdz(zp):\n",
    "    return (c/H00)/Ez(zp)\n",
    "def rcom(zc):\n",
    "    return scipy.integrate.romberg(drdz,0,zc)\n",
    "\n",
    "print rcom(ztest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Define the growth function in LCDM\n",
    "def fg(zz):\n",
    "    omz=om0*pow(1+zz,3)/(om0*pow(1+zz,3)+1-om0)\n",
    "    return pow(omz,gamma)\n",
    "\n",
    "print fg(ztest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Get the growth factor \n",
    "def Dg_dz(zz):\n",
    "    return fg(zz)/(1+zz)\n",
    "def Dgz(zc):\n",
    "    ans = scipy.integrate.romberg(Dg_dz, 0.0, zc)\n",
    "    return np.exp(-ans)\n",
    "\n",
    "print Dgz(ztest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#fiducial bHI from Bull et al 2015\n",
    "def bHI(zc):\n",
    "    return 0.67+0.18*zc+0.05*pow(zc,2)\n",
    "\n",
    "#fiducial OmHI Mario's fit\n",
    "def OmHI(zc):\n",
    "    return 0.00048+0.00039*zc-0.000065*pow(zc,2)\n",
    "\n",
    "print OmHI(ztest), bHI(ztest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Construct P(k,z) \n",
    "def Pkz(kk,zc):    \n",
    "    return pow(Dgz(zc),2)*Pkz0(kk)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#mean brightness temperature Mario's fit\n",
    "def Tb(zc): #in mK\n",
    "    return 0.0559+0.2324*zc-0.024*pow(zc,2)\n",
    "\n",
    "print ztest, Tb(ztest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Construct P_HI(k,z) [mK^2]\n",
    "def PHI(kk,zc):\n",
    "    return pow(Tb(zc),2)*pow(bHI(zc),2)*Pkz(kk,zc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#MeerKAT noise specs\n",
    "Ndishes=64\n",
    "Ddish=13.5*100 #cm\n",
    "Nbeams=1\n",
    "\n",
    "def thetab(zc):\n",
    "    return 21*(1+zc)/Ddish\n",
    "\n",
    "def omegapix(zc):\n",
    "    return 1.13*pow(thetab(zc),2)\n",
    "\n",
    "Area=4000.0 #deg^2\n",
    "omegatot = Area*pow(pi/180,2)\n",
    "ttotal = 4000*60*60*(Area/4000) #4000 hrs for 4000 deg^2\n",
    "\n",
    "def fc(zc):\n",
    "    return 1420.4/(1+zc)\n",
    "\n",
    "def Tsky(zc):\n",
    "    return 60*pow(300/fc(zc),2.55)*1e3\n",
    "\n",
    "Tsyslist = [23.5*1e3,23.0*1e3,23.0*1e3,28.0*1e3,29.0*1e3,30.0*1e3,28.5*1e3,29.5*1e3,\\\n",
    "            31.0*1e3,33.0*1e3,34.0*1e3,35.0*1e3,37.0*1e3,38.0*1e3] #mK\n",
    "\n",
    "def tobs(zc):\n",
    "    return ttotal*(omegapix(zc)/omegatot)*Ndishes*Nbeams\n",
    "\n",
    "Dzbin = 0.1\n",
    "dfpix = 50*1e3 #Hz\n",
    "midfreq = 1420.4e6 #Hz\n",
    "\n",
    "def dzpix(zc):\n",
    "    return pow(1+zc,2)*dfpix/midfreq\n",
    "def sigpix(zc,Tsys):\n",
    "    return Tsys/np.sqrt(dfpix*tobs(zc)) \n",
    "def dVpixdz(zz):    \n",
    "    return c*pow(rcom(zz),2)/(H00*Ez(zz))\n",
    "def Vpix(zc):\n",
    "    return omegapix(zc)*scipy.integrate.romberg(dVpixdz,zc-dzpix(zc)/2,zc+dzpix(zc)/2)\n",
    "\n",
    "def Wsq(kk,zc):\n",
    "    return np.exp(-pow(kk,2)*pow(rcom(zc),2)*pow(thetab(zc),2)/(8*np.log(2)))\n",
    "\n",
    "def Pnoise(kk,zc,Tsys):\n",
    "    return pow(sigpix(zc,Tsys),2)*Vpix(zc)*pow(Wsq(kk,zc),-1.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def kmin(zc):\n",
    "    return 2*pi/np.sqrt(pow(rcom(zc),2)*omegatot)\n",
    "def kmax(zc):\n",
    "    return 0.14*pow(1+zc,2/(2+nss)); #non-linear cutoff Smith et al 2003\n",
    "\n",
    "print kmin(ztest), kmax(ztest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#survey (bin) volume [Mpc^3]\n",
    "def dVsurdz(zz):    \n",
    "    return omegatot*c*pow(rcom(zz),2)/(H00*Ez(zz))\n",
    "    \n",
    "def Vsur(zc):\n",
    "    return scipy.integrate.romberg(dVsurdz,zc-Dzbin/2,zc+Dzbin/2)\n",
    "\n",
    "#effective volume going in the Fisher matrix\n",
    "def Veff(kk,zc):\n",
    "    return Vsur(zc)*(PHI(kk,zc)/(PHI(kk,zc)+Pnoise(kk,zc,Tsys)))**2\n",
    "\n",
    "print \"%.4g\" % Vsur(ztest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Fisher matrix derivatives\n",
    "\n",
    "def dlnP_dlnOmHIbHI(kk,zc):\n",
    "    return 2.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Fisher matrix\n",
    "\n",
    "def dFdk(kk):\n",
    "    return (1./(4*pi*pi))*pow(kk,2)*pow(dlnP_dlnOmHIbHI(kk,zc),2)*Veff(kk,zc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for i in range(0,Nzbins):\n",
    "    zc = zlist[i]\n",
    "    Tsys = Tsyslist[i]+Tsky(zc)\n",
    "    K = np.linspace(kmin(zc), kmax(zc), 100)\n",
    "    dF = dFdk(K)\n",
    "    Fisher = scipy.integrate.simps(dF,K)           \n",
    "    print zc, np.sqrt(1/Fisher)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}