{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1a342e58",
   "metadata": {},
   "source": [
    "# Astrophysics III: galaxy formation & evolution\n",
    "Moodle: https://go.epfl.ch/PHYS-465\n",
    "\n",
    "## Exercise 13\n",
    "Teaching assistant: Jonathan Petersson <br>\n",
    "Email: jonathan.petersson@epfl.ch"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55040aa8",
   "metadata": {},
   "source": [
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bcc04dfd",
   "metadata": {},
   "source": [
    "## Load data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3eaf8a3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "\n",
    "\n",
    "%matplotlib inline\n",
    "plt.rcParams['font.size'] = 14\n",
    "\n",
    "# Load data:\n",
    "data        = np.loadtxt('Mbh_Mstellar_box2-bao_144_gas_gt1e10_lumflag.dat')\n",
    "data = data[data[:,3] != -99]\n",
    "\n",
    "Mhalo     = data[:,0]\n",
    "Mstellar  = data[:,1]\n",
    "Mbh       = data[:,2]\n",
    "BHAccRate = data[:,3]\n",
    "CentOrSat = data[:,5]\n",
    "Nbh       = data[:,6]\n",
    "SFR       = data[:,7]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfc29864",
   "metadata": {},
   "source": [
    "## Exercise 13.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d0255e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot:\n",
    "fig, ax = plt.subplots(2, 2, figsize=(9, 7))\n",
    "\n",
    "# BH accretion rates:\n",
    "ax[0,0].hist(np.log10(BHAccRate[BHAccRate > 0]), bins=50, alpha=0.5, rasterized=True)\n",
    "ax[0,0].hist(np.log10(BHAccRate[(BHAccRate > 0) * (CentOrSat == 0)]), bins=50, alpha=0.5, rasterized=True)\n",
    "ax[0,0].hist(np.log10(BHAccRate[(BHAccRate > 0) * (CentOrSat == 0) * (Mbh > 1e6)]), bins=50, alpha=0.5, rasterized=True)\n",
    "ax[0,0].set_yscale('log')\n",
    "ax[0,0].set_xlabel('$\\log_{10}(\\dot{M}_\\mathrm{BH} \\ [\\mathrm{M}_\\odot \\ \\mathrm{yr}^{-1}])$')\n",
    "ax[0,0].set_ylabel('Number of BHs')\n",
    "\n",
    "# Eddington accretion rate (cgs units):\n",
    "G         = 6.67e-8\n",
    "mp        = 1.67262192e-24\n",
    "epsilon_r = 0.2\n",
    "c         = 3e10\n",
    "Thomson   = 6.6524e-25\n",
    "EddAccRate = (4 * np.pi * G * Mbh * 1.989e33 * mp) / (epsilon_r * c * Thomson)\n",
    "EddAccRate = EddAccRate / 1.989e33 * 365.25 * 24 * 60 * 60\n",
    "fedd = BHAccRate / EddAccRate\n",
    "\n",
    "ax[0,1].hist(np.log10(fedd[BHAccRate > 0]), bins=50, alpha=0.5, rasterized=True)\n",
    "ax[0,1].hist(np.log10(fedd[(BHAccRate > 0) * (CentOrSat == 0)]), bins=50, alpha=0.5, rasterized=True)\n",
    "ax[0,1].hist(np.log10(fedd[(BHAccRate > 0) * (CentOrSat == 0) * (Mbh > 1e6)]), bins=50, alpha=0.5, rasterized=True)\n",
    "ax[0,1].set_yscale('log')\n",
    "ax[0,1].set_xlabel('$\\log_{10}(f_\\mathrm{Edd})$')\n",
    "ax[0,1].set_ylabel('Number of BHs')\n",
    "\n",
    "# AGN luminosity:\n",
    "Lbol = epsilon_r * BHAccRate * 1.989e33 / (365.25 * 24 * 60 * 60) * c**2\n",
    "#Lbol = epsilon_r / (1 - epsilon_r) * BHAccRate * 1.989e33 / (365.25 * 24 * 60 * 60) * c**2\n",
    "#Ledd = (4 * np.pi * G * Mbh * 1.989e33 * c * mp) / Thomson\n",
    "#Lbol[fedd < 0.1] = 10 * Ledd[fedd < 0.1] * fedd[fedd < 0.1]**2\n",
    "\n",
    "ax[1,0].hist(np.log10(Lbol[(Lbol > 1e41) * (BHAccRate > 0)]), bins=50, alpha=0.5, label='All BHs', rasterized=True)\n",
    "ax[1,0].hist(np.log10(Lbol[(Lbol > 1e41) * (BHAccRate > 0) * (CentOrSat == 0)]), bins=50, alpha=0.5, label='Central BHs', rasterized=True)\n",
    "ax[1,0].hist(np.log10(Lbol[(Lbol > 1e41) * (BHAccRate > 0) * (CentOrSat == 0) * (Mbh > 1e6)]), bins=50, alpha=0.5, label='Central BHs > $10^6$ M$_\\odot$', rasterized=True)\n",
    "ax[1,0].set_yscale('log')\n",
    "ax[1,0].set_xlabel('$\\log_{10}(L_\\mathrm{bol} \\ [\\mathrm{erg \\ s}^{-1}])$')\n",
    "ax[1,0].set_ylabel('Number of BHs')\n",
    "fig.tight_layout()\n",
    "ax[1,0].legend(frameon=False, loc='upper left', bbox_to_anchor=(1, 1))\n",
    "\n",
    "fig.delaxes(ax[1,1])\n",
    "fig.savefig('ex1a.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4da5aacc",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Plot:\n",
    "fig, ax = plt.subplots(1, 2, figsize=(11, 5), sharex=True, sharey=True)\n",
    "\n",
    "h = ax[0].hist2d(np.log10(Mbh[BHAccRate > 0]), np.log10(BHAccRate[BHAccRate > 0]), \n",
    "                 bins=80, norm=mpl.colors.LogNorm( vmin=1, vmax=1e3), cmap='plasma', rasterized=True)\n",
    "h = ax[1].hist2d(np.log10(Mbh[(BHAccRate > 0) * (CentOrSat == 0)]), np.log10(BHAccRate[(BHAccRate > 0) * (CentOrSat == 0)]), \n",
    "                 bins=80, norm=mpl.colors.LogNorm(vmin=1, vmax=1e3), cmap='plasma', rasterized=True)\n",
    "\n",
    "ax[0].set_title('All BHs')\n",
    "ax[1].set_title('Central BHs > $10^6$ M$_\\odot$')\n",
    "ax[0].set_xlabel('$\\log_{10}(M_\\mathrm{BH} \\ [\\mathrm{M}_\\odot])$')\n",
    "ax[1].set_xlabel('$\\log_{10}(M_\\mathrm{BH} \\ [\\mathrm{M}_\\odot])$')\n",
    "ax[0].set_ylabel('$\\log_{10}(\\dot{M}_\\mathrm{BH} \\ [\\mathrm{M}_\\odot \\ \\mathrm{yr}^{-1}])$')\n",
    "divider = make_axes_locatable(ax[1])\n",
    "cax = divider.append_axes('right', size='5%', pad=0.05)\n",
    "fig.colorbar(h[3], cax=cax, label='Number of BHs')\n",
    "\n",
    "EddRatios = [f * (4 * np.pi * G * np.logspace(5, 11) * 1.989e33 * mp) / (epsilon_r * c * Thomson) \n",
    "             / 1.989e33 * 365.25 * 24 * 60 * 60 for f in [1e-7, 1e-5, 1e-3, 1e-1]]\n",
    "for i in range(0, 2):\n",
    "    ax[i].plot(np.log10(np.logspace(5, 11)), np.log10(EddRatios[0]), ls='-.', c='k', alpha=0.5, label='$f_\\mathrm{Edd} = 10^{-7}$')\n",
    "    ax[i].plot(np.log10(np.logspace(5, 11)), np.log10(EddRatios[1]), ls=':', c='k', alpha=0.5, label='$f_\\mathrm{Edd} = 10^{-5}$')\n",
    "    ax[i].plot(np.log10(np.logspace(5, 11)), np.log10(EddRatios[2]), ls='--', c='k', alpha=0.5, label='$f_\\mathrm{Edd} = 10^{-3}$')\n",
    "    ax[i].plot(np.log10(np.logspace(5, 11)), np.log10(EddRatios[3]), ls='-', c='k', alpha=0.5, label='$f_\\mathrm{Edd} = 10^{-1}$')\n",
    "\n",
    "ax[1].legend(frameon=False, loc='lower right')\n",
    "fig.tight_layout()\n",
    "fig.savefig('ex1b.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c9b0c4a",
   "metadata": {},
   "source": [
    "## Exercise 13.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e62328c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot:\n",
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "\n",
    "h = ax.hist2d(np.log10(Mstellar[(CentOrSat == 0) * (Mbh > 1e6)]), np.log10(Mbh[(CentOrSat == 0) * (Mbh > 1e6)]), \n",
    "              bins=80, norm=mpl.colors.LogNorm(), cmap='plasma', rasterized=True)\n",
    "ax.set_xlabel('$\\log_{10}(M_\\star \\ [\\mathrm{M}_\\odot])$')\n",
    "ax.set_ylabel('$\\log_{10}(M_\\mathrm{BH} \\ [\\mathrm{M_\\odot}])$')\n",
    "fig.colorbar(h[3], ax=ax, label='Number of BHs')\n",
    "\n",
    "# Stellar mass bins:\n",
    "MstellarBins = 10**h[1]\n",
    "MstellarBins_mid = (MstellarBins[:-1] + MstellarBins[1:]) / 2\n",
    "\n",
    "# Mean relation:\n",
    "mean = [np.mean(Mbh[(CentOrSat == 0) * (Mbh > 1e6) * (Mstellar > MstellarBins[i]) * (Mstellar < MstellarBins[i+1])]) \n",
    "        for i in range(len(MstellarBins[:-1]))]\n",
    "ax.plot(np.log10(MstellarBins_mid), np.log10(mean), ls=':', lw=2, c='k', label='Mean Relation')\n",
    "\n",
    "# McConell & Ma (2013):\n",
    "log10Mbh = 8.46 + 1.05 * np.log10(MstellarBins_mid / 1e11)\n",
    "ax.plot(np.log10(MstellarBins_mid), log10Mbh, ls='-', lw=2, c='k', label='McConell & Ma (2013)')\n",
    "\n",
    "ax.legend(frameon=False, loc='upper left')\n",
    "fig.tight_layout()\n",
    "fig.savefig('ex2.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72aaace0",
   "metadata": {},
   "source": [
    "## Exercise 13.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d06afbe6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot:\n",
    "fig, ax = plt.subplots(1, 3, figsize=(12, 4))\n",
    "\n",
    "# Masks:\n",
    "mask_tot = (CentOrSat == 0) * (Mbh > 1e6)\n",
    "mask_lum = mask_tot * (Lbol > 1e45)\n",
    "mask_mod = mask_tot * (Lbol > 1e42) * (Lbol < 1e45)\n",
    "\n",
    "# Fractions vs Stellar Mass:\n",
    "bins_stellar = np.linspace(10, 13, 6)\n",
    "bins_stellar_mid = (bins_stellar[:-1] + bins_stellar[1:])/2\n",
    "N_tot, bin_edges = np.histogram(np.log10(Mstellar[mask_tot]), bins=bins_stellar)\n",
    "N_lum, bin_edges = np.histogram(np.log10(Mstellar[mask_lum]), bins=bins_stellar)\n",
    "N_mod, bin_edges = np.histogram(np.log10(Mstellar[mask_mod]), bins=bins_stellar)\n",
    "f_lum = N_lum / N_tot\n",
    "f_mod = N_mod / N_tot\n",
    "\n",
    "ax[0].plot(bins_stellar_mid, f_lum, lw=2, marker='o')\n",
    "ax[0].plot(bins_stellar_mid, f_mod, lw=2, marker='o')\n",
    "ax[0].set_xlabel('$\\log_{10}(M_\\star \\ [\\mathrm{M}_\\odot])$')\n",
    "ax[0].set_ylabel('Fraction')\n",
    "ax[0].set_yscale('log')\n",
    "\n",
    "# Fractions vs SFR:\n",
    "bins_sfr = np.linspace(-1.7, 3, 6)\n",
    "bins_sfr_mid = (bins_sfr[:-1] + bins_sfr[1:])/2\n",
    "N_tot, bin_edges = np.histogram(np.log10(SFR[mask_tot * (SFR > 0)]), bins=bins_sfr)\n",
    "N_lum, bin_edges = np.histogram(np.log10(SFR[mask_lum * (SFR > 0)]), bins=bins_sfr)\n",
    "N_mod, bin_edges = np.histogram(np.log10(SFR[mask_mod * (SFR > 0)]), bins=bins_sfr)\n",
    "f_lum = N_lum / N_tot\n",
    "f_mod = N_mod / N_tot\n",
    "\n",
    "ax[1].plot(bins_sfr_mid, f_lum, lw=2, marker='o')\n",
    "ax[1].plot(bins_sfr_mid, f_mod, lw=2, marker='o')\n",
    "ax[1].set_xlabel('$\\log_{10}(\\mathrm{SFR} \\ [\\mathrm{M}_\\odot \\ \\mathrm{yr}^{-1}])$')\n",
    "ax[1].set_yscale('log')\n",
    "\n",
    "# Fractions vs SFR:\n",
    "sSFR = SFR / Mstellar\n",
    "bins_ssfr = np.linspace(-13, -9, 6)\n",
    "bins_ssfr_mid = (bins_ssfr[:-1] + bins_ssfr[1:])/2\n",
    "N_tot, bin_edges = np.histogram(np.log10(sSFR[mask_tot * (SFR > 0)]), bins=bins_ssfr)\n",
    "N_lum, bin_edges = np.histogram(np.log10(sSFR[mask_lum * (SFR > 0)]), bins=bins_ssfr)\n",
    "N_mod, bin_edges = np.histogram(np.log10(sSFR[mask_mod * (SFR > 0)]), bins=bins_ssfr)\n",
    "f_lum = N_lum / N_tot\n",
    "f_mod = N_mod / N_tot\n",
    "\n",
    "ax[2].plot(bins_ssfr_mid, f_lum, lw=2, marker='o', label='$>10^{45}$')\n",
    "ax[2].plot(bins_ssfr_mid, f_mod, lw=2, marker='o', label='$10^{42}-10^{45}$')\n",
    "ax[2].legend(fontsize=12, title_fontsize=12, frameon=False, title='$L_\\mathrm{AGN} \\ [\\mathrm{erg \\ s}^{-1}]$')\n",
    "ax[2].set_xlabel('$\\log_{10}(\\mathrm{sSFR} \\ [\\mathrm{yr}^{-1}])$')\n",
    "ax[2].set_yscale('log')\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig('ex3.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07ddc7fd",
   "metadata": {},
   "source": [
    "## Exercsie 13.4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "927ffd5e",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Plot:\n",
    "fig, ax = plt.subplots(1, 3, figsize=(12, 4), sharex=True, sharey=True)\n",
    "\n",
    "# Masks:\n",
    "mask_tot = (CentOrSat == 0) * (Mbh > 1e6) * (SFR > 0)\n",
    "mask_lum = mask_tot * (Lbol > 1e45)\n",
    "mask_mod = mask_tot * (Lbol > 1e42) * (Lbol < 1e45)\n",
    "\n",
    "# Number of galaxies / fraction of (moderately) luminous AGN:\n",
    "bins = [np.linspace(10, 13, 50), np.linspace(-1.7, 3, 50)]\n",
    "H_tot, xedges, yedges = np.histogram2d(np.log10(Mstellar[mask_tot]), np.log10(SFR[mask_tot]), bins=bins)\n",
    "H_lum, xedges, yedges = np.histogram2d(np.log10(Mstellar[mask_lum]), np.log10(SFR[mask_lum]), bins=bins)\n",
    "H_mod, xedges, yedges = np.histogram2d(np.log10(Mstellar[mask_mod]), np.log10(SFR[mask_mod]), bins=bins)\n",
    "X, Y = np.meshgrid(xedges, yedges)\n",
    "\n",
    "pc = ax[0].pcolormesh(X, Y, H_tot.T, norm=mpl.colors.LogNorm(), cmap='plasma', rasterized=True)\n",
    "fig.colorbar(pc, ax=ax[0], location='top', pad=0, label='Number of galaxies')\n",
    "ax[0].set_xlabel('$\\log_{10}(M_\\star \\ [\\mathrm{M}_\\odot])$')\n",
    "ax[0].set_ylabel('$\\log_{10}(\\mathrm{SFR} \\ [\\mathrm{M}_\\odot \\ \\mathrm{yr}^{-1}])$')\n",
    "\n",
    "pc = ax[1].pcolormesh(X, Y, H_lum.T/H_tot.T, norm=mpl.colors.LogNorm(vmin=1e-2, vmax=1e0), cmap='plasma', rasterized=True)\n",
    "fig.colorbar(pc, ax=ax[1], location='top', pad=0, label='Fraction of luminous AGN')\n",
    "ax[1].set_xlabel('$\\log_{10}(M_\\star \\ [\\mathrm{M}_\\odot])$')\n",
    "\n",
    "pc = ax[2].pcolormesh(X, Y, H_mod.T/H_tot.T, norm=mpl.colors.LogNorm(vmin=1e-2, vmax=1e0), cmap='plasma', rasterized=True)\n",
    "fig.colorbar(pc, ax=ax[2], location='top', pad=0, label='Fraction of moderately luminous AGN')\n",
    "ax[2].set_xlabel('$\\log_{10}(M_\\star \\ [\\mathrm{M}_\\odot])$')\n",
    "\n",
    "# MS fit & Cano-Díaz et al. (2016):\n",
    "fit = np.polyfit(np.log10(Mstellar[mask_tot * (Mstellar < 1e11) * (SFR > 1e0)]), np.log10(SFR[mask_tot  * (Mstellar < 1e11) * (SFR > 1e0)]), 1)\n",
    "ms = fit[0] * bins[0] + fit[1]\n",
    "for i in range(0, 3):\n",
    "    ax[i].plot(bins[0], ms, c='k', lw=2, ls=':', label='Main Sequence Fit')\n",
    "    ax[i].plot(bins[0], 0.81 * bins[0] - 8.34, c='k', lw=2, label='Cano-Díaz et al. (2016)')\n",
    "\n",
    "ax[2].legend(fontsize=12, loc='lower right')\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig('ex4.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aadd767a",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}