{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1a342e58",
   "metadata": {},
   "source": [
    "# Astrophysics III: galaxy formation & evolution\n",
    "Moodle: https://go.epfl.ch/PHYS-465\n",
    "\n",
    "## Exercise 11\n",
    "Teaching assistant: Jonathan Petersson <br>\n",
    "Email: jonathan.petersson@epfl.ch"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55040aa8",
   "metadata": {},
   "source": [
    "***"
   ]
  },
  {
   "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 notebook\n",
    "plt.rcParams['font.size'] = 14\n",
    "\n",
    "# Load data:\n",
    "data_dm_winds = np.loadtxt('DMParticles_M3852_winds_z0.dat')\n",
    "data_dm_nomw  = np.loadtxt('DMParticles_M3852_nomw_z0.dat')\n",
    "\n",
    "data_star_winds = np.loadtxt('StarParticles_M3852_winds_z0.dat')\n",
    "data_star_nomw  = np.loadtxt('StarParticles_M3852_nomw_z0.dat')\n",
    "\n",
    "data_gas_winds = np.loadtxt('GasParticles_M3852_winds_z0.dat')\n",
    "data_gas_nomw  = np.loadtxt('GasParticles_M3852_nomw_z0.dat')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e0a369f",
   "metadata": {},
   "source": [
    "## Exercise 11.1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e4ccfc8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Stellar Baryon Conversion Efficiency:\n",
    "Omega_b = 0.04087\n",
    "Omega_m = 0.259 \n",
    "fbar = Omega_b / Omega_m\n",
    "\n",
    "Rvir = 109  # [kpc/h]\n",
    "Rgal = 0.1 * Rvir\n",
    "\n",
    "Mstellar_winds = np.sum(data_star_winds[:,0][np.linalg.norm(data_star_winds[:,1:], axis=1) < Rgal])\n",
    "Mstellar_nomw =  np.sum(data_star_nomw[:,0][np.linalg.norm(data_star_nomw[:,1:], axis=1) < Rgal])\n",
    "\n",
    "Mdm_winds = np.sum(data_dm_winds[:,0][np.linalg.norm(data_dm_winds[:,1:], axis=1) < Rvir])\n",
    "Mdm_nomw = np.sum(data_dm_nomw[:,0][np.linalg.norm(data_dm_nomw[:,1:], axis=1) < Rvir])\n",
    "\n",
    "sbce_winds = Mstellar_winds / (fbar * Mdm_winds)\n",
    "sbce_nomw = Mstellar_nomw / (fbar * Mdm_nomw)\n",
    "\n",
    "print('Momentum-driven winds stellar feedback')\n",
    "print(f'Mhalo: {Mdm_winds}')\n",
    "print(f'sbce: {sbce_winds}\\n')\n",
    "\n",
    "print('Thermal (weak) stellar feedback')\n",
    "print(f'Mhalo: {Mdm_nomw}')\n",
    "print(f'sbce: {sbce_nomw}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e20a91f4",
   "metadata": {},
   "source": [
    "## Exercise 11.2:\n",
    "\n",
    "Kepler's third law states that:\n",
    "\n",
    "\\begin{equation}\n",
    "    \\frac{a^3}{T^2} = \\frac{G(M + m)}{4 \\pi^2}. \n",
    "\\end{equation}\n",
    "\n",
    "For a circular orbit where $a=r$, $T=(2\\pi r)/v_\\mathrm{circ}$ and $M >> m$, Kepler's law gives us that\n",
    "\n",
    "\\begin{equation}\n",
    "    \\frac{r^3}{\\left(\\frac{2\\pi r}{v_\\mathrm{circ}}\\right)^2} = \\frac{GM}{4\\pi^2} \\quad \\Rightarrow \\quad v_\\mathrm{circ} = \\sqrt{\\frac{GM}{r}}. \n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7fea6613",
   "metadata": {},
   "source": [
    "## Exercise 11.3:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44a904e7",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "G = 6.67430e-8   #[cgs]\n",
    "Msol = 1.989e33  #[g]\n",
    "kpc = 3.086e21   #[cm]\n",
    "\n",
    "# Circular velocity:\n",
    "radius = np.linspace(0, 20, 100)\n",
    "\n",
    "vc_winds = {\n",
    "    'dm'    : [],\n",
    "    'star'  : [],\n",
    "    'gas'   : [],\n",
    "    'total' : []\n",
    "}\n",
    "\n",
    "vc_nomw = {\n",
    "    'dm'    : [],\n",
    "    'star'  : [],\n",
    "    'gas'   : [],\n",
    "    'total' : []\n",
    "}\n",
    "\n",
    "for i in range(1, len(radius)):\n",
    "    mass_dm = np.sum(data_dm_winds[:,0][np.linalg.norm(data_dm_winds[:,1:], axis=1) < radius[i]])\n",
    "    mass_star = np.sum(data_star_winds[:,0][np.linalg.norm(data_star_winds[:,1:], axis=1) < radius[i]])\n",
    "    mass_gas = np.sum(data_gas_winds[:,0][np.linalg.norm(data_gas_winds[:,1:], axis=1) < radius[i]])\n",
    "    mass_total = mass_dm + mass_star + mass_gas\n",
    "    \n",
    "    vc_winds['dm'].append(np.sqrt(G * mass_dm * 1e10* Msol / (radius[i] * kpc)) / 1e5)\n",
    "    vc_winds['star'].append(np.sqrt(G * mass_star * 1e10 * Msol / (radius[i] * kpc)) / 1e5)\n",
    "    vc_winds['gas'].append(np.sqrt(G * mass_gas * 1e10 * Msol / (radius[i] * kpc)) / 1e5)\n",
    "    vc_winds['total'].append(np.sqrt(G * mass_total * 1e10 * Msol / (radius[i] * kpc)) / 1e5)\n",
    "    \n",
    "    mass_dm = np.sum(data_dm_nomw[:,0][np.linalg.norm(data_dm_nomw[:,1:], axis=1) < radius[i]])\n",
    "    mass_star = np.sum(data_star_nomw[:,0][np.linalg.norm(data_star_nomw[:,1:], axis=1) < radius[i]])\n",
    "    mass_gas = np.sum(data_gas_nomw[:,0][np.linalg.norm(data_gas_nomw[:,1:], axis=1) < radius[i]])\n",
    "    mass_total = mass_dm + mass_star + mass_gas\n",
    "    \n",
    "    vc_nomw['dm'].append(np.sqrt(G * mass_dm * 1e10* Msol / (radius[i] * kpc)) / 1e5)\n",
    "    vc_nomw['star'].append(np.sqrt(G * mass_star * 1e10 * Msol / (radius[i] * kpc)) / 1e5)\n",
    "    vc_nomw['gas'].append(np.sqrt(G * mass_gas * 1e10 * Msol / (radius[i] * kpc)) / 1e5)\n",
    "    vc_nomw['total'].append(np.sqrt(G * mass_total * 1e10 * Msol / (radius[i] * kpc)) / 1e5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f7f0c8c",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(10, 4), sharey=True, layout='constrained')\n",
    "\n",
    "ax[0].plot(radius[1:], vc_winds['dm'], lw=2, c='blue', alpha=0.5, label='DM')\n",
    "ax[0].plot(radius[1:], vc_winds['star'], lw=2, c='red', alpha=0.5, label='Stars')\n",
    "ax[0].plot(radius[1:], vc_winds['gas'], lw=2, c='orange', alpha=0.5, label='Gas')\n",
    "ax[0].plot(radius[1:], vc_winds['total'], lw=2, c='k', alpha=0.5, label='Total')\n",
    "ax[0].set_xlabel('Radius [kpc]')\n",
    "ax[0].set_ylabel('$v_\\mathrm{circ} \\ [\\mathrm{km \\ s}^{-1}]$')\n",
    "ax[0].set_title('winds')\n",
    "ax[0].legend(frameon=False)\n",
    "\n",
    "ax[1].plot(radius[1:], vc_nomw['dm'], lw=2, c='blue', alpha=0.5,)\n",
    "ax[1].plot(radius[1:], vc_nomw['star'], lw=2, c='red', alpha=0.5)\n",
    "ax[1].plot(radius[1:], vc_nomw['gas'], lw=2, c='orange', alpha=0.5)\n",
    "ax[1].plot(radius[1:], vc_nomw['total'], lw=2, c='k', alpha=0.5)\n",
    "ax[1].set_xlabel('Radius [kpc]')\n",
    "ax[1].set_title('nomw')\n",
    "\n",
    "fig.savefig('ex3.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f8549f5",
   "metadata": {},
   "source": [
    "## Exercise 11.4:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "214e08bb",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 6), layout='constrained')\n",
    "\n",
    "# Avila-Reese et al. (2008):\n",
    "mstellar = np.logspace(9, 11, 100)\n",
    "logvmax = 0.274 * np.log10(mstellar) - 0.639\n",
    "err = 0.058\n",
    "ax.fill_between(np.log10(mstellar), logvmax-err, logvmax+err, color='k', edgecolor='none', alpha=0.2)\n",
    "ax.plot(np.log10(mstellar), logvmax, ls='--', lw=2, c='k', label='Avila-Reese+2008')\n",
    "\n",
    "# Maximum circular velocity & stellar mass:\n",
    "vmax_winds = np.max(vc_winds['total'])\n",
    "mstellar_winds = np.sum(data_star_winds[:,0][np.linalg.norm(data_star_winds[:,1:], axis=1) < Rgal]) * 1e10\n",
    "ax.scatter(np.log10(mstellar_winds), np.log10(vmax_winds), s=300, marker='*', c='blue', label='winds')\n",
    "\n",
    "vmax_nomw = np.max(vc_nomw['total'])\n",
    "mstellar_nomw = np.sum(data_star_nomw[:,0][np.linalg.norm(data_star_nomw[:,1:], axis=1) < Rgal]) * 1e10\n",
    "ax.scatter(np.log10(mstellar_nomw), np.log10(vmax_nomw), s=300, marker='*', c='red', label='nomw')\n",
    "\n",
    "ax.set_xlabel('$\\log_{10}(M_\\star \\ [\\mathrm{M}_\\odot])$')\n",
    "ax.set_ylabel('$\\log_{10}(v_\\mathrm{max} \\ [\\mathrm{km s}^{-1}])$')\n",
    "ax.set_xlim(9, 11)\n",
    "ax.legend(frameon=False)\n",
    "\n",
    "fig.savefig('ex4.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3cd84d31",
   "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
}