SKA1-MID-sd.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Author: Alkistis Pourtsidou, ICG Portsmouth\n",
    "\n",
    "Using part of http://camb.readthedocs.io/en/latest/CAMBdemo.html \n",
    "\n",
    "To run this Jupyter 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": 1,
   "metadata": {
    "collapsed": false
   },
   "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": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import scipy\n",
    "from scipy.interpolate import interp1d\n",
    "from __future__ import division"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "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": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy import integrate\n",
    "from scipy import linalg\n",
    "\n",
    "pi=np.pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Fiducial cosmological parameters\n",
    "c=3e5\n",
    "hubble=0.675\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",
    "\n",
    "nss = 0.9655\n",
    "alphass = 0.0\n",
    "betass = 0.0\n",
    "\n",
    "gamma=0.545\n",
    "\n",
    "kpivot = 0.05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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, nrun=alphass, nrunrun=betass, pivot_scalar=kpivot)\n",
    "pars.set_for_lmax(2500, lens_potential_accuracy=0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#calculate results for these parameters\n",
    "results = camb.get_results(pars)\n",
    "#print pars"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Get matter power spectrum at z=0: P(k,z=0)\n",
    "#Not non-linear corrections couples to smaller scales than you want\n",
    "pars.set_matter_power(redshifts=[0.], kmax=1.0)\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-5, maxkh=0.5, npoints = 800)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Derivatives\n",
    "\n",
    "def dlogPd_dns(kk):\n",
    "    return np.log(kk/kpivot)\n",
    "\n",
    "def dlogPd_dalphas(kk):\n",
    "    return (1./2)*pow(np.log(kk/kpivot),2)\n",
    "\n",
    "def dlogPd_dbetas(kk):\n",
    "    return (1./6)*pow(np.log(kk/kpivot),3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.35  0.45  0.55  0.65  0.75  0.85  0.95  1.05  1.15  1.25  1.35  1.45\n",
      "  1.55  1.65  1.75  1.85  1.95  2.05  2.15  2.25  2.35  2.45  2.55  2.65\n",
      "  2.75  2.85  2.95]\n",
      "ztest = 0.5\n",
      "Number of redshift bins = 27\n"
     ]
    }
   ],
   "source": [
    "#Redshift bins for SKA1-MID (Band 1)\n",
    "\n",
    "zlist = np.arange(0.35,3.05,0.1)\n",
    "ztest = 0.5\n",
    "Nzbins = len(zlist)\n",
    "\n",
    "print zlist\n",
    "print \"ztest =\", ztest\n",
    "print \"Number of redshift bins =\", Nzbins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Spatially flat Universe\n",
    "\n",
    "#Define E(z) = H(z)/H0\n",
    "def Ez(zc):\n",
    "    return np.sqrt(1-om0+om0*pow(1+zc,3))\n",
    "def Hz(zc):\n",
    "    return Ez(zc)*H00\n",
    "#Define the comoving distances\n",
    "def drdz(zp):\n",
    "    return (c/H00)/Ez(zp)\n",
    "def rcom(zc):\n",
    "    return scipy.integrate.romberg(drdz,0,zc)\n",
    "def DA(zc):\n",
    "    return rcom(zc)/(1+zc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.758362943343\n"
     ]
    }
   ],
   "source": [
    "#Define the growth function in LCDM\n",
    "def fg(zz):\n",
    "    omz=om0*pow(1+zz,3)/pow(Ez(zz),2)\n",
    "    return pow(omz,gamma)\n",
    "\n",
    "print fg(ztest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.769389642041\n"
     ]
    }
   ],
   "source": [
    "#Get the growth factor \n",
    "def Dg_dz(zz):\n",
    "    return -fg(zz)/(1+zz)\n",
    "\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": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7725 0.00065875\n"
     ]
    }
   ],
   "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 bHI(ztest), OmHI(ztest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Construct P(k,μ,z) \n",
    "def beta(zc):\n",
    "    return fg(zc)/bHI(zc)\n",
    "\n",
    "def Pkz(kk,mu,zc):    \n",
    "    return pow(1+beta(zc)*mu**2,2)*pow(Dgz(zc),2)*Pkz0(kk)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5 0.1661\n"
     ]
    }
   ],
   "source": [
    "#mean brightness temperature [mK] Mario's fit\n",
    "def Tb(zc):\n",
    "    return 0.0559+0.2324*zc-0.024*pow(zc,2)\n",
    "\n",
    "print ztest, Tb(ztest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5\n"
     ]
    }
   ],
   "source": [
    "#Construct P_HI(k,μ,z), Tb assumed known\n",
    "def PHI(kk,mu,zc):\n",
    "    return pow(bHI(zc),2)*Pkz(kk,mu,zc)\n",
    "\n",
    "print ztest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#noise specs\n",
    "\n",
    "Ndishes=130 + 64\n",
    "Ddish=15.0*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",
    "fsky = 0.7\n",
    "Area=fsky*41253 #deg^2\n",
    "omegatot = Area*pow(pi/180,2)\n",
    "ttotal = 5000*60*60 \n",
    "\n",
    "def Tsys(zc):\n",
    "    Tinst = 25.0*1e3 #mK\n",
    "    fc = 1420.4/(1+zc)\n",
    "    Tsky = 60*pow(300/fc,2.55)*1e3 #mK\n",
    "    return Tinst + Tsky \n",
    "\n",
    "def tobs(zc):\n",
    "    return ttotal*(omegapix(zc)/omegatot)*Ndishes*Nbeams\n",
    "\n",
    "Dzbin = 0.1 #redshift bin width\n",
    "dfpix = 50*1e3 #frequency resolution in Hz\n",
    "midfreq = 1420.4e6 #Hz\n",
    "\n",
    "def dzpix(zc):\n",
    "    return pow(1+zc,2)*dfpix/midfreq\n",
    "def sigpix(zc):\n",
    "    return Tsys(zc)/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,mu,zc):\n",
    "    #add small offset to avoid division by zero\n",
    "    return 1e-20+np.exp(-pow(kk,2)*(1-mu**2)*pow(rcom(zc),2)*pow(thetab(zc),2)/(8*np.log(2)))\n",
    "\n",
    "#note division by Tb^2\n",
    "def Pnoise(kk,mu,zc):\n",
    "    return pow(sigpix(zc),2)*Vpix(zc)*pow(Wsq(kk,mu,zc),-1.)/pow(Tb(zc),2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#survey (bin) volume [Mpc^3]\n",
    "\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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.00279771839403 0.131449794558\n"
     ]
    }
   ],
   "source": [
    "def kmin(zc):\n",
    "    return 2*pi*pow(Vsur(zc),-1/3)\n",
    "\n",
    "kNL = 0.1 #1/Mpc    \n",
    "    \n",
    "def kmax(zc):\n",
    "    return kNL*pow(1+zc,2/(2+nss)) #non-linear cutoff \n",
    "\n",
    "print kmin(ztest),kmax(ztest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#effective volume going in the Fisher matrix\n",
    "def Veff(kk,mu,zc):\n",
    "    return Vsur(zc)*(PHI(kk,mu,zc)/(PHI(kk,mu,zc)+Pnoise(kk,mu,zc)))**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def dF(kk,mu):\n",
    "    return (1./(8*pi*pi))*pow(kk,2)*deriv_i(kk,mu,zc)*deriv_j(kk,mu,zc)*Veff(kk,mu,zc)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#2D integration function\n",
    "def integrate2D(dfun, kgrid, mugrid):\n",
    "    \n",
    "    muint = [scipy.integrate.simps(dfun.T[i], mugrid) for i in range(kgrid.size)]\n",
    "    return scipy.integrate.simps(muint, kgrid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mugrid = np.linspace(-1., 1., 513) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Fisher matrix parameters \n",
    "params = [\"0:ns\",\"1:alphas\",\"2:betas\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#%%time\n",
    "#   Fisher matrix   # \n",
    "\n",
    "Npar = 2 #Npar = 3\n",
    "\n",
    "#create array of zeros\n",
    "s = Npar\n",
    "\n",
    "Fishermat = np.zeros([s,s])\n",
    "\n",
    "for zi in range(0,Nzbins):\n",
    "    zc = zlist[zi]\n",
    "    kgrid = np.linspace(kmin(zc), kmax(zc), 513)\n",
    "    if kpivot < kmin(zc) or kpivot > kmax(zc): raise ValueError(\"kpivot out of range\")\n",
    "    K, MU = np.meshgrid(kgrid, mugrid)\n",
    "    for i in range(0,Npar):  \n",
    "        def deriv_i(kk,mu,zc):\n",
    "            if i==0:  return dlogPd_dns(kk)\n",
    "            elif i==1:  return dlogPd_dalphas(kk)\n",
    "            elif i==2:  return dlogPd_dbetas(kk)\n",
    "            else: print \"out of range\"\n",
    "        for  j in range(0,Npar):\n",
    "            if j>=i:\n",
    "                def deriv_j(kk,mu,zc):\n",
    "                    if j==0:  return dlogPd_dns(kk)\n",
    "                    elif j==1:  return dlogPd_dalphas(kk)\n",
    "                    elif j==2:  return dlogPd_dbetas(kk)\n",
    "                    else: print \"index out of range\" \n",
    "                Fishermat[i][j] += integrate2D(dF(K,MU),kgrid,mugrid)\n",
    "            else: Fishermat[i,j] = Fishermat[j,i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 219943.24188762   26231.04205861]\n",
      " [  26231.04205861   57704.94483958]]\n"
     ]
    }
   ],
   "source": [
    "print Fishermat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  4.80724511e-06  -2.18523818e-06]\n",
      " [ -2.18523818e-06   1.83228851e-05]]\n"
     ]
    }
   ],
   "source": [
    "#print covariance matrix\n",
    "print linalg.inv(Fishermat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.00219254306977\n",
      "0.00428052392954\n"
     ]
    }
   ],
   "source": [
    "#print marginalised uncertainties\n",
    "for i in range(0,Npar):\n",
    "    print np.sqrt(linalg.inv(Fishermat)[i,i])"
   ]
  },
  {
   "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
}