2.evolution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "147142d6-5c04-4671-bc6e-0473d638307a",
   "metadata": {},
   "source": [
    "# Initialisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d28674b5-d64c-4c19-b1c3-d2d89ed5cf3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib ipympl\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "#import scipy\n",
    "from scipy.special import zeta\n",
    "from time import time\n",
    "from astropy import cosmology as cosmo\n",
    "from astropy import units as u\n",
    "from astropy import constants as c\n",
    "from astropy.visualization import quantity_support\n",
    "quantity_support()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de6ae4b0-00d9-4d12-acb2-8c389a100f31",
   "metadata": {},
   "source": [
    "## Models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6fd15638-bc27-48f5-a6c2-fb5298b04799",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6aaecbfa-9701-44f1-8b6c-bd453ab30f03",
   "metadata": {},
   "outputs": [],
   "source": [
    "from models import flcdm, coasting"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11985058-d2b0-4192-b6ce-2821193bfa27",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82b3fd78-9578-4d44-8fcf-423bf43793de",
   "metadata": {},
   "outputs": [],
   "source": [
    "SN_z, SN_err_z, SN_mu, SN_err_mu = np.loadtxt('data/Pantheon+_Data.dat', usecols=(2, 3, 10, 11), unpack=True, skiprows=1)\n",
    "GRB1 = np.loadtxt('data/GRB_RNN.dat', skiprows=1)\n",
    "GRB2 = np.loadtxt('data/GRB_RNN_2.dat', skiprows=1)\n",
    "GRB_z = np.concatenate([GRB1[:, 0], GRB1[:, 3], GRB1[:, 6], GRB2[:, 0]])\n",
    "GRB_mu = np.concatenate([GRB1[:, 1], GRB1[:, 4], GRB1[:, 7], GRB2[:, 1]])\n",
    "GRB_err_mu = np.concatenate([GRB1[:, 2], GRB1[:, 5], GRB1[:, 8], GRB2[:, 2]])\n",
    "QSO_z, QSO_mu, QSO_err_mu = np.loadtxt('data/table3QSO.dat', usecols=(3, 11, 12), unpack=True, skiprows=1)\n",
    "#QSO_err_mu *= 2\n",
    "\n",
    "def Dl_Mpc(mu):\n",
    "    return 1e-5 * 10**(.2*mu)\n",
    "\n",
    "#SN_Dl = Dl_Mpc(SN_mu)\n",
    "#GRB_Dl = Dl_Mpc(GRB_mu)\n",
    "#QSO_Dl = Dl_Mpc(QSO_mu)\n",
    "\n",
    "\n",
    "all_z = np.concatenate([SN_z, QSO_z, GRB_z])\n",
    "all_mu = np.concatenate([SN_mu, QSO_mu, GRB_mu])\n",
    "all_err_mu = np.concatenate([SN_err_mu, QSO_err_mu, GRB_err_mu])\n",
    "\n",
    "#all_z = SN_z\n",
    "#all_mu = SN_mu\n",
    "#all_err_mu = SN_err_mu\n",
    "#all_z = QSO_z\n",
    "#all_mu = QSO_mu\n",
    "#all_err_mu = QSO_err_mu\n",
    "#all_z = GRB_z\n",
    "#all_mu = GRB_mu\n",
    "#all_err_mu = GRB_err_mu\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9d90afa-5a00-4746-a2da-708796d3fe3a",
   "metadata": {},
   "outputs": [],
   "source": [
    "thetaBAO = np.array([\n",
    "    0.11,19.8,3.26,\n",
    "    0.235,9.06,0.23,\n",
    "    0.365,6.33,0.22,\n",
    "    0.45,4.77,0.17,\n",
    "    0.47,5.02,0.25,\n",
    "    0.49,4.99,0.21,\n",
    "    0.51,4.81,0.17,\n",
    "    0.53,4.29,0.30,\n",
    "    0.55,4.25,0.25,\n",
    "    0.57,4.59,0.36,\n",
    "    0.59,4.39,0.33,\n",
    "    0.61,3.85,0.31,\n",
    "    0.63,3.90,0.43,\n",
    "    0.65,3.55,0.16,\n",
    "    2.225,1.77,0.31])\n",
    "\n",
    "theta_z = thetaBAO[::3]\n",
    "theta_z_err = np.array([0.005,0.035,0.025,0.01,0.005,0.01,0.005,0.005,0.005,0.005,0.005,0.005,0.005,0.01,0.025])\n",
    "theta_angle_deg = thetaBAO[1::3]\n",
    "theta_angle_err = thetaBAO[2::3]\n",
    "#theta_relative_err = theta_angle_err / theta_angle_deg\n",
    "#theta_DM = 180/np.pi/theta_angle_deg"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d85dc01-b250-4988-9a89-310f1bb611c4",
   "metadata": {},
   "source": [
    "# Plot cosmic evolution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed71e385-f51e-4529-997e-202e600006a2",
   "metadata": {},
   "source": [
    "## Distances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "706d02b4-ac11-49ee-9a0a-19761405e5b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig_name = 'distances'\n",
    "plt.close(fig_name)\n",
    "fig = plt.figure(fig_name, figsize=(8, 9))\n",
    "axes = fig.subplots(nrows=3, ncols=1, squeeze=False,\n",
    "                    sharex='col', sharey='row',\n",
    "                    gridspec_kw={'hspace': 0, 'wspace': 0},\n",
    "                   )\n",
    "\n",
    "ax = axes[0, 0]\n",
    "ax.set_ylabel(r'radial comoving distance $D_C$ [Mpc]')\n",
    "ax.set_yscale('log')\n",
    "ax.set_ylim(3e-1, 3e5)\n",
    "\n",
    "ax.plot(coasting.z_t, coasting.Dc_t.to_value(u.Mpc), 'r-', label='coasting')\n",
    "ax.plot(flcdm.z_t, flcdm.Dc_t.to_value(u.Mpc), 'r--', label=r'flat $\\Lambda$CDM')\n",
    "ax.legend()\n",
    "\n",
    "#ax.plot(flcdm.z_t, cosmo.Planck18.comoving_distance(flcdm.z_t).to_value(u.Mpc), 'k+', alpha=.1)\n",
    "#DH = (c.c/coasting.H0).to_value(u.Mpc)\n",
    "#ax.plot(coasting.z_t, np.log(1+coasting.z_t)*DH, 'k+', alpha=.1)\n",
    "\n",
    "\n",
    "ax = axes[1, 0]\n",
    "ax.set_ylabel(r'angular diameter distance $D_A$ [Mpc]')\n",
    "ax.set_yscale('log')\n",
    "ax.set_ylim(3e-1, 3e3)\n",
    "\n",
    "ax.plot(coasting.z_t, coasting.Dm_t.to_value(u.Mpc) / (1+coasting.z_t), 'k-', label='coasting')\n",
    "ax.plot(flcdm.z_t, flcdm.Dm_t.to_value(u.Mpc) / (1+flcdm.z_t), 'k--', label=r'flat $\\Lambda$CDM')\n",
    "ax.legend()\n",
    "\n",
    "#ax.plot(flcdm.z_t, cosmo.Planck18.angular_diameter_distance(flcdm.z_t).to_value(u.Mpc), 'k+', alpha=.1)\n",
    "#z_cmb = np.interp(flcdm.t_cmb, flcdm.t, flcdm.z_t)\n",
    "#print(z_cmb, cosmo.Planck18.angular_diameter_distance(z_cmb).to_value(u.Mpc)*(1+z_cmb))\n",
    "#Ok = .05\n",
    "#DH = (c.c/coasting.H0).to_value(u.Mpc)\n",
    "#ax.plot(coasting.z_t, np.abs(np.sin(np.log((1+coasting.z_t)) * np.sqrt(Ok)) / np.sqrt(Ok))*DH / (1+coasting.z_t), 'k+', alpha=.1)\n",
    "#z_cmb = 1/(coasting.t_cmb * coasting.H0).to_value(u.dimensionless_unscaled)\n",
    "#print(z_cmb, np.abs(np.sin(np.log((1+z_cmb)) * np.sqrt(Ok)) / np.sqrt(Ok))*DH)\n",
    "\n",
    "ax = axes[2, 0]\n",
    "ax.set_ylabel(r'ratio (coasting / flat $\\Lambda$CDM)')\n",
    "#ax.set_yscale('log')\n",
    "ax.set_ylim(.75, 1.25)\n",
    "\n",
    "z = np.logspace(-4, 4, 401)\n",
    "ax.plot(z, coasting.comoving_distance(z) / flcdm.comoving_distance(z), 'r-', label=r'$D_C$')\n",
    "ax.plot(z, coasting.angular_diameter_distance(z) / flcdm.angular_diameter_distance(z), 'k-', label='$D_A$ and $D_L$')\n",
    "ax.legend()\n",
    "\n",
    "\n",
    "ax.set_xlabel('redshift $z$')\n",
    "ax.set_xscale('log')\n",
    "ax.set_xlim(3e-4, 3e3)\n",
    "\n",
    "for ax in axes.flatten():\n",
    "    ax.tick_params(which='both', bottom=True, top=True, left=True, right=True)\n",
    "    ax.tick_params(which='major', direction='inout', length=8, grid_alpha=.3)\n",
    "    ax.tick_params(which='minor', direction='in', length=2, grid_alpha=.1)\n",
    "    ax.grid(True, which='both')\n",
    "\n",
    "fig.set_tight_layout(True)\n",
    "plt.show()\n",
    "plt.savefig('figures/distances.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ff868c5-2589-4833-b14c-e337a41d6288",
   "metadata": {},
   "source": [
    "## Densities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c47e56b5-be58-4d66-a083-8f28d81d22d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig_name = 'evolution'\n",
    "plt.close(fig_name)\n",
    "fig = plt.figure(fig_name, figsize=(8, 8))\n",
    "axes = fig.subplots(nrows=2, ncols=1, squeeze=False,\n",
    "                    sharex='col', sharey='row',\n",
    "                    gridspec_kw={'hspace': 0, 'wspace': 0},\n",
    "                   )\n",
    "alpha = 1\n",
    "\n",
    "ax = axes[0, 0]\n",
    "ax.set_ylabel(r'density $\\rho$ [M$_\\odot$ Mpc$^{-3}$]')\n",
    "\n",
    "def mass_to_E(m):\n",
    "    return((m*c.c**2 * u.Msun/u.Mpc**3).to_value(u.GeV/u.cm**3))\n",
    "def E_to_mass(E):\n",
    "    return((E/c.c**2 * u.GeV/u.cm**3).to_value(u.Msun/u.Mpc**3))\n",
    "secax = ax.secondary_yaxis('right', functions=(mass_to_E, E_to_mass))\n",
    "secax.set_ylabel(r'energy density $\\rho c^2$ [GeV cm$^{-3}$]')\n",
    "\n",
    "\n",
    "ax.set_yscale('log')\n",
    "ax.set_xscale('log')\n",
    "\n",
    "legends = []\n",
    "for model, title, ls, bbox in zip([flcdm, coasting],\n",
    "                                  [r'Flat $\\Lambda$CDM', 'Coasting'],\n",
    "                                  ['--', '-'],\n",
    "                                  [(0.45, 0.5, 0.5, 0.5), (0.72, 0.5, 0.5, 0.5)],\n",
    "                                 ):\n",
    "    t = model.t.to_value(u.yr)\n",
    "    rho_c = model.rho_crit_t.to_value(u.Msun/u.Mpc**3)\n",
    "    rho_de = np.exp(model.ln_Ode_t) * rho_c\n",
    "    rho_m = np.exp(model.ln_Om_t) * rho_c\n",
    "    rho_b = rho_m * model.Ob0/model.Om0\n",
    "    rho_gamma = np.exp(model.ln_Ogamma_t) *rho_c\n",
    "    \n",
    "    l1, = ax.plot(t, rho_de, c='r', ls=ls, alpha=alpha, label='dark energy')\n",
    "    l2, = ax.plot(t, rho_m - rho_b, c='k', ls=ls, alpha=alpha, label='dark matter')\n",
    "    l3, = ax.plot(t, rho_b, c='c', ls=ls, alpha=alpha, label='baryons (atoms)')\n",
    "    l4, = ax.plot(t, rho_gamma, c='y', ls=ls, alpha=alpha, label='photons')\n",
    "    legends.append(ax.legend(handles=[l1, l2, l3, l4], title=title, loc='upper left', bbox_to_anchor=bbox))\n",
    "\n",
    "ax.add_artist(legends[0])\n",
    "ax.set_xlim(3e3, 3e10)\n",
    "ax.set_ylim(3e6, 3e25)\n",
    "ax.set_yticks(np.logspace(7, 25, 19), minor=True, labels='')\n",
    "\n",
    "ax = axes[1, 0]\n",
    "ax.set_ylabel('ionisation fraction $x$')\n",
    "\n",
    "ax.plot(flcdm.t.to_value(u.yr), flcdm.x_ion_t, 'c--', alpha=alpha)\n",
    "ax.plot(coasting.t.to_value(u.yr), coasting.x_ion_t, 'c-', alpha=alpha)\n",
    "\n",
    "ax.plot(flcdm.t_cmb.to_value(u.yr), np.interp(flcdm.t_cmb, flcdm.t, flcdm.x_ion_t), 'cd', fillstyle='none',\n",
    "           label=r'Flat $\\Lambda$CDM: $t_{\\rm CMB} \\sim$'+f'{flcdm.t_cmb.to_value(u.Myr):.3f} Myr')\n",
    "ax.plot(coasting.t_cmb.to_value(u.yr), np.interp(coasting.t_cmb, coasting.t, coasting.x_ion_t), 'cd',\n",
    "        label=r'Coasting: $t_{\\rm CMB} \\sim$'+f'{coasting.t_cmb.to_value(u.Myr):.3f} Myr')\n",
    "\n",
    "l = ax.legend().get_lines()\n",
    "l[0].set_ls('--')\n",
    "l[1].set_ls('-')\n",
    "#h.set_visible(False)\n",
    "\n",
    "ax.set_xscale('log')\n",
    "ax.set_xlabel('cosmic time $t$ [yr]')\n",
    "\n",
    "\n",
    "for ax in axes.flatten():\n",
    "    #ax.tick_params(which='both', bottom=True, top=True, left=True, right=True)\n",
    "    ax.tick_params(which='major', direction='inout', length=8, grid_alpha=.3)\n",
    "    ax.tick_params(which='minor', direction='in', length=2, grid_alpha=.1)\n",
    "    ax.grid(True, which='both')\n",
    "    #ax.axvline(t_cmb_p18, c='r', ls='--')\n",
    "    #ax.axvline(t_cmb_model, c='k', ls='-')\n",
    "fig.set_tight_layout(True)\n",
    "plt.show()\n",
    "plt.savefig('figures/evolution.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61eeb4ea-7521-48bf-a224-b39b83cd61d4",
   "metadata": {},
   "source": [
    "# Compare with observations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8836c9e5-9097-4f4e-8043-2e94c487ddc6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from models import flcdm, coasting\n",
    "\n",
    "def plot_DL_DH_obs(ax, z, mu, err_mu, **kwargs):\n",
    "    #scale = 3e5 * z\n",
    "    scale = z\n",
    "    y = Dl_Mpc(mu) / scale\n",
    "    dy_low = y-Dl_Mpc(mu-err_mu) / scale\n",
    "    dy_high = Dl_Mpc(mu+err_mu) / scale - y\n",
    "    ax.errorbar(z, y, yerr=(dy_low, dy_high), fmt='.', **kwargs)\n",
    "    \n",
    "    \n",
    "fig_name = 'observations'\n",
    "\n",
    "plt.close(fig_name)\n",
    "fig = plt.figure(fig_name, figsize=(8, 9))\n",
    "axes = fig.subplots(nrows=2, ncols=1, squeeze=False,\n",
    "                    sharex=True, sharey='row',\n",
    "                    gridspec_kw={'hspace': 0, 'wspace': 0},\n",
    "                   )\n",
    "\n",
    "ax = axes[0, 0]\n",
    "ax.set_ylabel(r'$D_L(z)/z$ [Mpc]')\n",
    "ax.set_ylim(0, 2.1e4)\n",
    "#ax.set_yscale('log')\n",
    "\n",
    "plot_DL_DH_obs(ax, SN_z, SN_mu, SN_err_mu, c='c', alpha=.05, label='Pantheon+ SN Ia (Brout et al. 2022)')\n",
    "plot_DL_DH_obs(ax, QSO_z, QSO_mu, QSO_err_mu, c='y', alpha=.1, label='QSO (Lusso et al. 2020)')\n",
    "plot_DL_DH_obs(ax, GRB_z, GRB_mu, GRB_err_mu, c='r', alpha=.2, label='GRB (Escanilla-Rivera et al. 2022)')\n",
    "\n",
    "flcdm.fit_mu(SN_z, SN_mu, SN_err_mu)\n",
    "flcdm.fit_mu(QSO_z, QSO_mu, QSO_err_mu)\n",
    "flcdm.fit_mu(GRB_z, GRB_mu, GRB_err_mu)\n",
    "flcdm.fit_mu(all_z, all_mu, all_err_mu)\n",
    "print()\n",
    "coasting.fit_mu(SN_z, SN_mu, SN_err_mu)\n",
    "coasting.fit_mu(QSO_z, QSO_mu, QSO_err_mu)\n",
    "coasting.fit_mu(GRB_z, GRB_mu, GRB_err_mu)\n",
    "coasting.fit_mu(all_z, all_mu, all_err_mu)\n",
    "print()\n",
    "ax.plot(flcdm.z_t, flcdm.Dm_t.to_value(u.Mpc) * (1+flcdm.z_t)/flcdm.z_t, 'k--', label=r'flat $\\Lambda$CDM')\n",
    "ax.plot(coasting.z_t, coasting.Dm_t.to_value(u.Mpc) * (1+coasting.z_t)/coasting.z_t, 'k-', label='coasting')\n",
    "\n",
    "ax.legend()\n",
    "\n",
    "\n",
    "ax = axes[1, 0]\n",
    "ax.set_ylabel(r'$\\theta_{BAO}$ [$^\\circ$]')\n",
    "ax.set_ylim(0, 24)\n",
    "#ax.set_yscale('log')\n",
    "\n",
    "ax.errorbar(theta_z, theta_angle_deg, theta_angle_err, theta_z_err, fmt='b.', label=r'$\\theta_{BAO}$ (Nunes et al. 2020)')\n",
    "\n",
    "ax.plot(flcdm.z_t, (flcdm.r_d / flcdm.Dm_t * (1+flcdm.z_t) * u.rad).to_value(u.deg), 'k--', label=r'flat $\\Lambda$CDM')\n",
    "ax.plot(coasting.z_t, (coasting.r_d / coasting.Dm_t * (1+coasting.z_t) * u.rad).to_value(u.deg), 'k-', label=r'coasting')\n",
    "\n",
    "flcdm.fit_theta_BAO(theta_z, theta_angle_deg, theta_angle_err)\n",
    "coasting.fit_theta_BAO(theta_z, theta_angle_deg, theta_angle_err)\n",
    "\n",
    "ax.legend()\n",
    "\n",
    "\n",
    "ax.set_xlabel('redshift $z$')\n",
    "ax.set_xscale('log')\n",
    "ax.set_xlim(1e-3, 20)\n",
    "#ax.set_xlim(-.2, 8.2)\n",
    "\n",
    "for ax in axes.flatten():\n",
    "    ax.tick_params(which='both', direction='in')\n",
    "    ax.grid(alpha=.5)\n",
    "fig.set_tight_layout(True)\n",
    "plt.show()\n",
    "#plt.savefig('DL.pdf')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27bea391-3649-4179-a435-1bb4fa049bad",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.10"
  },
  "toc-autonumbering": true
 },
 "nbformat": 4,
 "nbformat_minor": 5
}