{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1a342e58",
   "metadata": {},
   "source": [
    "# Astrophysics III: galaxy formation & evolution\n",
    "Moodle: https://go.epfl.ch/PHYS-465\n",
    "\n",
    "## Exercise 9\n",
    "Teaching assistant: Jonathan Petersson <br>\n",
    "Email: jonathan.petersson@epfl.ch"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55040aa8",
   "metadata": {},
   "source": [
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4318cec",
   "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",
    "%matplotlib notebook\n",
    "plt.rcParams['font.size'] = 14\n",
    "\n",
    "\n",
    "# Load data:\n",
    "data = np.loadtxt('Gasparticles_M0215_BH_z0.8_corr.dat')\n",
    "mass = data[:,0]    # [1e10 Msol/h]\n",
    "pos  = data[:,1:4]  # [kpc/h]\n",
    "hsml = data[:,4]\n",
    "rho  = data[:,5]    # [(1e10 Msol/h)/((kpc/h)^3)]\n",
    "temp = data[:,6]    # [K]\n",
    "sfr  = data[:,7]    # [Msol / yr]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ed10e8c",
   "metadata": {},
   "source": [
    "## Exercise 9.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0acf0a4b",
   "metadata": {},
   "outputs": [],
   "source": [
    "Rvir = 195         # [kpc/h]\n",
    "Rgal = 0.1 * Rvir  # [kpc/h]\n",
    "Tcut = 1e4         # [K]\n",
    "\n",
    "# Fraction of ALL gas residing within the galaxy:\n",
    "frac_all = np.sum(mass[np.linalg.norm(pos, axis=1) < Rgal]) / np.sum(mass[np.linalg.norm(pos, axis=1) < Rvir])\n",
    "print(f'Fraction of ALL gas residing within the galaxy: {round(frac_all*100, 2)} %')\n",
    "\n",
    "# Fraction of HOT gas residing within the galaxy:\n",
    "frac_hot = np.sum(mass[(np.linalg.norm(pos, axis=1) < Rgal) * (temp > Tcut)]) / np.sum(mass[(np.linalg.norm(pos, axis=1) < Rvir) * (temp > Tcut)])\n",
    "print(f'Fraction of HOT gas residing within the galaxy: {round(frac_hot*100, 2)} %')\n",
    "\n",
    "# Fraction of COLD gas residing within the galaxy:\n",
    "frac_cold = np.sum(mass[(np.linalg.norm(pos, axis=1) < Rgal) * (temp < Tcut)]) / np.sum(mass[(np.linalg.norm(pos, axis=1) < Rvir) * (temp < Tcut)])\n",
    "print(f'Fraction of COLD gas residing within the galaxy: {round(frac_cold*100, 2)} %')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec80269b",
   "metadata": {},
   "source": [
    "## Exercise 9.2a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "244bc847",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mask your data to be within Rgal (0.1*Rvir)\n",
    "x = pos[:,0][np.linalg.norm(pos, axis=1) < Rgal]\n",
    "y = pos[:,1][np.linalg.norm(pos, axis=1) < Rgal]\n",
    "z = pos[:,2][np.linalg.norm(pos, axis=1) < Rgal]\n",
    "t = np.log10(temp[np.linalg.norm(pos, axis=1) < Rgal])\n",
    "\n",
    "# Plot:\n",
    "fig, ax = plt.subplots(1, 3, figsize=(10, 3), sharex=True, sharey=True)\n",
    "\n",
    "ax[0].scatter(x, y, s=2, c=t, alpha=0.5, cmap='plasma', vmin=2.5, vmax=7.5, rasterized=True)\n",
    "ax[0].set_xlabel('$x$ [kpc/h]')\n",
    "ax[0].set_ylabel('$y$ [kpc/h]')\n",
    "ax[0].set_aspect('equal')\n",
    "\n",
    "ax[1].scatter(x, z, s=2, c=t, alpha=0.5, cmap='plasma', vmin=2.5, vmax=7.5, rasterized=True)\n",
    "ax[1].set_xlabel('$x$ [kpc/h]')\n",
    "ax[1].set_ylabel('$z$ [kpc/h]')\n",
    "ax[1].set_aspect('equal')\n",
    "\n",
    "sc = ax[2].scatter(y, z, s=2, c=t, alpha=0.5, cmap='plasma', vmin=2.5, vmax=7.5, rasterized=True)\n",
    "ax[2].set_xlabel('$y$ [kpc/h]')\n",
    "ax[2].set_ylabel('$z$ [kpc/h]')\n",
    "ax[2].set_aspect('equal')\n",
    "\n",
    "divider = make_axes_locatable(ax[2])\n",
    "cax = divider.append_axes('right', size='5%', pad=0.05)\n",
    "cbar = fig.colorbar(sc, cax=cax, label='$\\log_{10}(T \\ [\\mathrm{K}])$')\n",
    "cbar.solids.set(alpha=1)\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig('ex2a.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2361cb9",
   "metadata": {},
   "source": [
    "## Exercise 9.2b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33b9dcaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 3, figsize=(10, 3), sharex=True, sharey=True)\n",
    "\n",
    "alpha = 52\n",
    "cos_alpha = np.cos(alpha*np.pi/180) \n",
    "sin_alpha = np.sin(alpha*np.pi/180)\n",
    "xr = x*cos_alpha - y*sin_alpha \n",
    "yr = x*sin_alpha + y*cos_alpha\n",
    "\n",
    "ax[0].scatter(xr, yr, s=2, c=t, alpha=0.5, cmap='plasma', vmin=2.5, vmax=7.5, rasterized=True)\n",
    "ax[0].set_xlabel('$x$ [kpc/h]')\n",
    "ax[0].set_ylabel('$y$ [kpc/h]')\n",
    "ax[0].set_aspect('equal')\n",
    "\n",
    "ax[1].scatter(xr, z, s=2, c=t, alpha=0.5, cmap='plasma', vmin=2.5, vmax=7.5, rasterized=True)\n",
    "ax[1].set_xlabel('$x$ [kpc/h]')\n",
    "ax[1].set_ylabel('$z$ [kpc/h]')\n",
    "ax[1].set_aspect('equal')\n",
    "\n",
    "sc = ax[2].scatter(yr, z, s=2, c=t, alpha=0.5, cmap='plasma', vmin=2.5, vmax=7.5, rasterized=True)\n",
    "ax[2].set_xlabel('$y$ [kpc/h]')\n",
    "ax[2].set_ylabel('$z$ [kpc/h]')\n",
    "ax[2].set_aspect('equal')\n",
    "\n",
    "divider = make_axes_locatable(ax[2])\n",
    "cax = divider.append_axes('right', size='5%', pad=0.05)\n",
    "cbar = fig.colorbar(sc, cax=cax, label='$\\log_{10}(T \\ [\\mathrm{K}])$')\n",
    "cbar.solids.set(alpha=1)\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig('ex2b.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44f1a163",
   "metadata": {},
   "source": [
    "## Exercise 9.3a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "586d1043",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Convert rho into a number density:\n",
    "n = rho * 1e10 * 1.989e33 / (3.086e21)**3 / (1.67262192 * 1e-24)\n",
    "\n",
    "# Plot:\n",
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "ax.scatter(np.log10(n), np.log10(hsml), s=2, c='blue', alpha=0.2, rasterized=True)\n",
    "ax.set_xlabel('$\\log_{10}(n \\ [(\\mathrm{cm/h})^{-3}])$')\n",
    "ax.set_ylabel('$\\log_{10}(\\mathrm{hsml} \\ [\\mathrm{pc/h}])$')\n",
    "fig.tight_layout()\n",
    "fig.savefig('ex3a.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e43ea50",
   "metadata": {},
   "source": [
    "## Exercise 9.3b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54d15b90",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Plot:\n",
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "\n",
    "h = ax.hist2d(np.log10(n), np.log10(temp), weights=mass, bins=200, norm=mpl.colors.LogNorm(), cmap='plasma', rasterized=True)\n",
    "ax.set_xlabel('$\\log_{10}(n \\ [(\\mathrm{cm/h})^{-3}])$')\n",
    "ax.set_ylabel('$\\log_{10}(T \\ [K])$')\n",
    "fig.colorbar(h[3], ax=ax, label='$M_\\mathrm{Gas} \\ [10^{10} \\ \\mathrm{M}_\\odot/\\mathrm{h}]$')\n",
    "fig.tight_layout()\n",
    "fig.savefig('ex3b.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86f23e12",
   "metadata": {},
   "source": [
    "## Exercise 9.4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09748dbe",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bin your data radially (in x & y):\n",
    "R = np.linalg.norm(pos[:,:2], axis=1)\n",
    "h, bin_edges = np.histogram(R[R < Rvir], bins=200, weights=mass[R < Rvir])\n",
    "bin_mids = (bin_edges[:-1] + bin_edges[1:]) / 2\n",
    "\n",
    "# Convert to the correct unit:\n",
    "h = 1e10 * h / (np.pi * ((bin_edges[1:] * 1e3)**2 - (bin_edges[:-1] * 1e3)**2))\n",
    "\n",
    "# Plot:\n",
    "fig, ax = plt.subplots(1, 2, figsize=(9, 4), sharey=True)\n",
    "ax[0].plot(bin_mids, np.log10(h))\n",
    "ax[0].set_xlabel('Radius [kpc/h]')\n",
    "ax[0].set_ylabel('$\\log_{10}(\\Sigma_\\mathrm{Gas} \\ [\\mathrm{M}_\\odot/\\mathrm{h} \\ (\\mathrm{pc/h})^{-2}]$)')\n",
    "ax[0].axvline(Rvir, c='k', ls=':', lw=2, alpha=0.5)\n",
    "ax[0].axvline(Rgal, c='k', ls='--', lw=2, alpha=0.5)\n",
    "\n",
    "\n",
    "ax[1].plot(np.log10(bin_mids), np.log10(h))\n",
    "ax[1].set_xlabel('$\\log_{10}(\\mathrm{Radius} \\ [\\mathrm{kpc/h}])$')\n",
    "ax[1].axvline(np.log10(Rvir), c='k', ls=':', lw=2, alpha=0.5, label='$R_\\mathrm{vir}$')\n",
    "ax[1].axvline(np.log10(Rgal), c='k', ls='--', lw=2, alpha=0.5, label='$0.1R_\\mathrm{vir}$')\n",
    "ax[1].legend(frameon=False)\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig('ex4.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8557a30",
   "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
}