{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "388c21c5",
   "metadata": {},
   "source": [
    "# Inter-Stellar Medium (ISM) tutorial "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16f80465",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from galapy.analysis.plot import plt, format_axes_ticks"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76c45c8f",
   "metadata": {},
   "source": [
    "The `galapy.InterStellarMedium` module defines classes and functions for modelling the gas and dust components impact on the total spectrum of a galaxy.\n",
    "\n",
    "It can be imported via:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dde7b96c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import galapy.InterStellarMedium as gpism"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b94822ca",
   "metadata": {},
   "source": [
    "Or called from the global namespace of the library:\n",
    "\n",
    "```python\n",
    ">>> import galapy as gp\n",
    ">>> gp.InterStellarMedium\n",
    "<module 'galapy.InterStellarMedium' from '/path/to/installed/version/of/galapy/InterStellarMedium.py'>\n",
    "```\n",
    "\n",
    "An object of type `ISM` is instantiated via"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce552361",
   "metadata": {},
   "outputs": [],
   "source": [
    "ism = gpism.ISM()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7560ce02",
   "metadata": {},
   "source": [
    "In `galapy`, the ISM is modelled as a medium constituted of two phases, a **Molecular Cloud** (MC) component, embedding star-forming regions, and a **Diffuse Dust** (DD) component which is roughly considered as a *spherical* region centred at the centre of the galaxy.\n",
    "\n",
    "The emission from stars has to first exit the MCs and then escape the dusty diffuse region. The attenuation from the two regions is therefore not independed as only the light escaping from MCs is then attenuated by the DD region."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40d19afa",
   "metadata": {},
   "source": [
    "Let us define some useful quantities:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0097071",
   "metadata": {},
   "source": [
    "* Wavelength grid:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f52f92fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "ll = np.logspace( 0, 10, 300 )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00b23571",
   "metadata": {},
   "source": [
    "* Luminosity unit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a11a268f",
   "metadata": {},
   "outputs": [],
   "source": [
    "Lsun = 3.828e+33"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "356e2997",
   "metadata": {},
   "source": [
    "The default parameters of the model can be printed by calling the dedicated function of the module:\n",
    "\n",
    "* **Molecular-Clouds**: key = `'mc'`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21237ebb",
   "metadata": {},
   "outputs": [],
   "source": [
    "gpism.ism_build_params(phase='mc')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "739d6895",
   "metadata": {},
   "source": [
    "* **Diffuse-Dust**: key = `dd`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "188760cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "gpism.ism_build_params(phase='dd')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbe36f0b",
   "metadata": {},
   "source": [
    "Notice that the two parameter dictionaries share the key-value pair `'f_MC' = 0.5` which regulates the relative contribution of the two components to the overall attenuation in terms of the fraction contributed by the Molecular-Clouds only (i.e. in this case, $50\\%$).\n",
    "\n",
    "\\[Parameterizing in terms of a single component prevents the user to accidentally set values whose sum is greater than 1\\]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c0a9f0d",
   "metadata": {},
   "source": [
    "Let's go through some of the main quantities that can be retrieved with the `ISM` type.\n",
    "\n",
    "### ISM Emission\n",
    "\n",
    "This quantity is approximated with a grey-body emission depending on the model parameters and on the temperature of the medium. The latter quantity has to be set by the user by calling the dedicated function (the argument is assumed in units of Kelvin):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9762bb2",
   "metadata": {},
   "outputs": [],
   "source": [
    "ism.mc.set_temperature( 30 )\n",
    "ism.dd.set_temperature( 38 )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f4a9d0d",
   "metadata": {},
   "source": [
    "the emission of the two different components at varying wavelength can then be obtained with the function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4eb8cb42",
   "metadata": {},
   "outputs": [],
   "source": [
    "Lmc = ism.mc.emission(ll)\n",
    "Ldd = ism.dd.emission(ll)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a69cf26d",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plt.ylim(35,45)\n",
    "plt.xlabel( '$\\\\log\\\\ \\\\lambda\\\\ [\\\\AA]$' )\n",
    "plt.ylabel( '$\\\\log\\\\ \\\\lambda\\\\; L_\\\\lambda\\\\ [$erg$\\\\; s^{-1}]$' )\n",
    "plt.plot( np.ma.log10( ll ), np.ma.log10( ll * Lsun * Ldd ), label='Diffuse Dust' )\n",
    "plt.plot( np.ma.log10( ll ), np.ma.log10( ll * Lsun * Lmc ), label='Molecular Cloud' )\n",
    "plt.plot( np.ma.log10( ll ), np.ma.log10( ll * Lsun * (Ldd+Lmc) ), label='Total', c='k', ls='--', lw=1.5 )\n",
    "_ = plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab3fc9e2",
   "metadata": {},
   "source": [
    "### Attenuation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6e3c89b",
   "metadata": {},
   "outputs": [],
   "source": [
    "Amc = ism.mc.attenuation(ll)\n",
    "Add = ism.dd.attenuation(ll)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "175abe67",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.xlabel( '$\\\\log\\\\ \\\\lambda\\\\ [\\\\AA]$' )\n",
    "plt.plot( np.ma.log10( ll ), Add, label='Diffuse Dust' )\n",
    "plt.plot( np.ma.log10( ll ), Amc, label='Molecular Cloud' )\n",
    "_ = plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bcac5954",
   "metadata": {},
   "source": [
    "### Extinction\n",
    "\n",
    "We compute this quantity in units of its value in the visible-band "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51e891ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "Extmc = ism.mc.extinction(ll) / ism.mc.A_V()\n",
    "Extdd = ism.dd.extinction(ll) / ism.dd.A_V()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7662964",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.xlim(3,7)\n",
    "plt.ylim(-5,1)\n",
    "plt.xlabel( '$\\\\log\\\\ \\\\lambda\\\\ [\\\\AA]$' )\n",
    "plt.ylabel( '$\\\\log\\\\ A_\\\\lambda/A_V$' )\n",
    "plt.plot( np.ma.log10( ll ), np.ma.log10( Extdd ), label='Diffuse Dust' )\n",
    "plt.plot( np.ma.log10( ll ), np.ma.log10( Extmc ), label='Molecular Cloud' )\n",
    "plt.axvline( np.log10(5500), c='k', ls='--', label='$5500\\\\ \\\\AA$' )\n",
    "plt.axvline( 6, c='k', ls=':', label = '$10^6\\\\ \\\\AA$')\n",
    "_ = plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d49629fa",
   "metadata": {},
   "source": [
    "### Total attenuation\n",
    "\n",
    "As this quantity depends on the time spent by the photons within the molecular cloud through the parameter `tau_esc`, it is a function of time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f0a8815",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_esc = ism.mc.params['tau_esc']\n",
    "tunit = np.logspace( np.log10(.8), np.log10(4), 100 ) # time array in units of tau_esc\n",
    "tt = tunit * t_esc                                    # time array in unit of [yr]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27a3efa5",
   "metadata": {},
   "source": [
    "The function `total_attenuation` of the `ISM` class returns a tuple with 2 matrices, of shape `(len(ll), len(tt))` where `ll` is the wavelength grid and `tt` the time grid, containing the contribute to attenuation from MC-only and the global attenuation itself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0daa7c75",
   "metadata": {},
   "outputs": [],
   "source": [
    "attTotMC, attTot = ism.total_attenuation( ll, tt )\n",
    "attTotMC = attTotMC.reshape(ll.size, tt.size)\n",
    "attTot   = attTot.reshape(ll.size, tt.size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cef70ea9",
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.colors import Normalize\n",
    "fig, ax = plt.subplots( 1,1 )\n",
    "\n",
    "# Plot\n",
    "ax.set_xlim( 4, 7 )\n",
    "ax.set_ylim( -2, 0.3 )\n",
    "ax.set_prop_cycle( color=plt.cm.viridis(np.linspace(0.01,0.99,len(tt))) )\n",
    "ax.set_xlabel( '$\\\\log \\\\lambda\\\\ [\\\\AA]$', fontsize = 16 )\n",
    "#ax.set_ylabel( '$\\\\log\\\\ S_\\\\nu\\\\ [m\\\\,\\\\mathrm{Jy}]$', fontsize = 16 )\n",
    "ax.set_title(\"Total attenuation (MC+DD)\", fontsize=16)\n",
    "for _t, _a in zip( tunit, attTot.T ) :\n",
    "    ax.plot(np.ma.log10( ll ), np.ma.log10( _a ) )\n",
    "        \n",
    "# Color-map\n",
    "cmap = plt.cm.viridis\n",
    "norm = Normalize( vmin=np.min(tunit), vmax=np.max(tunit) ) # normalization\n",
    "cax = fig.add_axes([0.42, 0.25, 0.46, 0.04])               # position\n",
    "cbar = fig.colorbar(plt.cm.ScalarMappable(norm=norm, cmap=cmap),\n",
    "                    cax=cax, \n",
    "                    orientation='horizontal', \n",
    "                    ticks=[1, 2, 3, 4])\n",
    "_ = cbar.ax.set_title( \"$\\\\tau / \\\\tau_\\\\mathrm{esc}$\", fontsize=16)\n",
    "\n",
    "# Formatting\n",
    "format_axes_ticks(fig)"
   ]
  }
 ],
 "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}