{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "df7238a6",
   "metadata": {},
   "source": [
    "# Choose the Simple Stellar Population library\n",
    "\n",
    "We provide a set of Simple Stellar Population libraries, already formatted for usage by ``galapy`` functions and classes.\n",
    "\n",
    "We can print a list of available libraries by calling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "901a3ade",
   "metadata": {},
   "outputs": [],
   "source": [
    "from galapy.CompositeStellarPopulation import print_ssp_libs\n",
    "print_ssp_libs()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2b62551",
   "metadata": {},
   "source": [
    "with the naming convention: **author.method\\[.imf.filling_schema\\]**\n",
    "\n",
    "* **family**: An achronym for the authors of the SSP library\n",
    "* **method**: if present, shows which method was used to compute the SSP\n",
    "* **imf**: if present shows the initial mass function used to compute the SSP\n",
    "* **filling_schema**: all the SSPs' wavelength domain has been extended to\n",
    "   span from $1\\ \\mathring{A}$ to $10^{10}\\ \\mathring{A}$. \n",
    "   This code provides the strategy used (_not present_ = filled with zeros, \n",
    "   ``extrap`` = extrapolated linearly in the logarithm, \n",
    "   ``refined`` = thinner lambda grid obtained by \n",
    "   linearly-interpolating the ``extend`` equivalent)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b97c0b5b",
   "metadata": {},
   "source": [
    "They are divided into two main families providing different features each. The details are reported in the following sections."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1874167a",
   "metadata": {},
   "source": [
    ">   Note that the different families are divided in flavours all of which have different properties. \n",
    ">   The main property that regulates the impact on the sampling performances is the size of the wavelength grid: \n",
    ">   the smaller the fastest."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64da814d",
   "metadata": {},
   "source": [
    "**Look around:**\n",
    "\n",
    "We can load a SSP table by invoking the dedicated function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33dcfbd7",
   "metadata": {},
   "outputs": [],
   "source": [
    "from galapy.CompositeStellarPopulation import load_SSP_table"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c45f760",
   "metadata": {},
   "source": [
    "which returns a tuple with \n",
    "\n",
    "* the wavelength grid\n",
    "* the age grid\n",
    "* the metallicity grid\n",
    "* the table itself\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2da1fbfc",
   "metadata": {},
   "outputs": [],
   "source": [
    "l, t, Z, L = load_SSP_table('bc03.basel.chab')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4494a0db",
   "metadata": {},
   "source": [
    "Our tables are stored as flattened, contiguous objects, mainly for internal usage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f63a071",
   "metadata": {},
   "outputs": [],
   "source": [
    "L.shape, l.size*t.size*Z.size"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3cbe2e5",
   "metadata": {},
   "source": [
    "We can reshape them into a 3D matrix by calling the dedicated function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d09ee455",
   "metadata": {},
   "outputs": [],
   "source": [
    "from galapy.CompositeStellarPopulation import reshape_SSP_table"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "251a683a",
   "metadata": {},
   "source": [
    "this function takes two arguments:\n",
    "\n",
    "* the ``L`` table\n",
    "* a *shape* tuple consisting of 3 elements: the wavelength grid size, the age grid size and the metallicity grid size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b74a979",
   "metadata": {},
   "outputs": [],
   "source": [
    "L = reshape_SSP_table( L, shape=(l.size, t.size, Z.size) )\n",
    "L.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7ebad9b",
   "metadata": {},
   "source": [
    "## Bruzual & Charlot 2003: ``bc03``\n",
    "\n",
    "the classic and popular [Bruzual & Charlot (2003)](https://academic.oup.com/mnras/article/344/4/1000/968846) in its updated version (v2016). \n",
    "This library provides the continuum luminosity from SSPs for a set of different IMFs, at varying wavelength, age and metallicity.\n",
    "\n",
    "> ``bc03`` SSP tables only include stellar emission (atmospheres) thus, if the user wants to also model further\n",
    "> emission components due to stars, these have to be computed separately.\n",
    "> Namely it is required to build the synchrotron and nebular free-free contributions separately\n",
    "> (see ``galapy.Synchrotron`` and ``galapy.NebularFreeFree``).\n",
    "\n",
    "\n",
    "### General to all the libraries in the family:\n",
    "\n",
    "|           | size | values |\n",
    "|-----------|-----|-----------------------------|\n",
    "| **Time grid** | 221 | irregularly distributed values in the interval $0 < \\tau < 2\\times10^{10}\\ \\text{years}$ |\n",
    "| **Metallicity grid** | 7 | ``[0.0001, 0.0004, 0.004, 0.008, 0.02, 0.05, 0.1]`` |\n",
    "\n",
    "### Specific to the different flavours:\n",
    "\n",
    "* ``basel.chab`` : lower resolution wavelength grid ($N_\\lambda < 3000$), with Chabrier IMF\n",
    "* ``stelib.chab``: higher resolution wavelength grid ($N_\\lambda > 7000$), with Chabrier IMF\n",
    "\n",
    "| flavour | Wavelength grid | Description |\n",
    "|---------|-----------------|-------------|\n",
    "| ``basel.chab`` | 2223 | Extended from the original work to span from $1 \\mathring{A}$ to $10^{10}\\mathring{A}$, the filling method is zeros padding |\n",
    "| ``basel.chab.refined`` | 2837 | Same as above but the density of wavelength grid points has been increased to have at least 128 values per-each order of magnitude |\n",
    "| ``stelib.chab`` | 7325 | Same as ``basel.chab`` but with more resolution in the Optical/UV bands |\n",
    "| ``stelib.chab.extrap`` | 7325 | Same as above but instead of zeros padding, the padding values have been extrapolated linearly in log-space from the last two grid points in the original table |\n",
    "\n",
    "### Tips\n",
    "\n",
    "* the SSP libraries of choice for comparison with other works (as these tables are present in most of the other SED fitting tools)\n",
    "* highly tested and verified\n",
    "* the ``stelib`` flavour includes the tables with the highest wavelength resolution from those included in the default database.\n",
    "* the ``stelib.*.extrap`` flavour introduces extrapolation at long and short wavelength instead of padding with zeros the table.\n",
    "* flavours marked ``refined`` should be used when computing the transmitted band flux in photometric systems including bands from the submm/mm bands and longer (lower frequency). \n",
    "\n",
    "### Plot single SSPs\n",
    "\n",
    "It might be useful to visualise the different SSPs before choosing which one best suits someones needs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e6c28c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from galapy.analysis.plot import plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9de8e1c",
   "metadata": {},
   "source": [
    "(the ``plt`` module in ``galapy.analysis.plot`` is just a ``matplotlib.pyplot`` where some tweaking of the default layout has been applied)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4099f87a",
   "metadata": {},
   "source": [
    "Load the two ``stelib`` flavours, to show the difference between the default zeros-padding and the extrapolated library:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "212b2c1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "l, t, Z, Lbc03zerosp = load_SSP_table('bc03.stelib.chab')\n",
    "Lbc03zerosp = reshape_SSP_table( Lbc03zerosp, shape=(l.size, t.size, Z.size) )\n",
    "*_, Lbc03extrap = load_SSP_table('bc03.stelib.chab.extrap')\n",
    "Lbc03extrap = reshape_SSP_table( Lbc03extrap, shape=(l.size, t.size, Z.size) )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18015bfe",
   "metadata": {},
   "source": [
    "(since the wavelength, age and metallicity grids are all the same for the two cases, we store them just from the first table)\n",
    "\n",
    "We select some different ages from the table picking some indices and check the value from the age-grid ``t``"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad7d1232",
   "metadata": {},
   "outputs": [],
   "source": [
    "tidx = [1, 100, 150, 200]\n",
    "iz = 3\n",
    "print( f'Metallicity: {Z[iz]}')\n",
    "print( 'Ages:')\n",
    "for it in tidx :\n",
    "    print( f'{t[it]:.2e} years' )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ce2fa44",
   "metadata": {},
   "source": [
    "And finally, we plot in 4 different panels the four different ages.\n",
    "Note how the tails of the emissions are different (solid blue line VS dashed green line)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb21a2a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(2,2, figsize=(14,8), sharex=True, sharey=True, \n",
    "                        gridspec_kw={'hspace':0.0, 'wspace':0.0}, \n",
    "                        tight_layout=True)\n",
    "\n",
    "for ii, it in enumerate(tidx) :\n",
    "    jj = int(ii/2)\n",
    "    kk = int(ii%2)\n",
    "    if jj != 0 :\n",
    "        axs[jj,kk].set_xlabel('$\\\\lambda\\\\ [\\\\AA]$')\n",
    "    if kk == 0 :\n",
    "        axs[jj,kk].set_ylabel('$L_\\\\lambda\\\\ [L_\\\\odot\\\\;\\\\AA^{-1}]$')\n",
    "    axs[jj,kk].set_xscale('log')\n",
    "    axs[jj,kk].set_yscale('log')\n",
    "    axs[jj,kk].set_xlim(2.*10**0,10**10)\n",
    "    axs[jj,kk].plot(l, Lbc03zerosp[:,it,3], label='zeros padding')\n",
    "    axs[jj,kk].plot(l, Lbc03extrap[:,it,3], ls='--', label='extrapolated')\n",
    "    if ii == 3 :\n",
    "        axs[jj,kk].legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acfe2468",
   "metadata": {},
   "source": [
    "## PARSEC22: ``parsec22``\n",
    "\n",
    "Produced with the PARSEC code ([Bressan et al., 2012](https://ui.adsabs.harvard.edu/abs/2012MNRAS.427..127B/abstract); [Chen et al., 2014](https://academic.oup.com/mnras/article/444/3/2525/1061098); [Chen et al., 2015](https://academic.oup.com/mnras/article/452/1/1068/1749582)) for a Chabrier IMF and varying ages and metallicities, including emission from dusty AGB stars [Bressan et al., 1998](https://ui.adsabs.harvard.edu/abs/1998A%26A...332..135B/abstract). \n",
    "\n",
    "These libraries come in two flavours, the first one with continuum emission only and the second also including nebular emission.\n",
    "In the former, besides continuum stellar emission, non-thermal synchrotron emission from core-collapse supernovae is also included in each SSP spectrum (see, e.g., [Vega et al., 2008](https://www.aanda.org/articles/aa/pdf/2008/24/aa8883-07.pdf)).\n",
    "    In the latter, on top of the stellar continuum and non-thermal synchrotron, nebular emission is also included, with both free-free continuum and nebular emission (see, e.g., [Maya et al., 2004](https://ui.adsabs.harvard.edu/abs/2004ApJ...600..188M/abstract)), calculated with [CLOUDY](https://gitlab.nublado.org/cloudy/cloudy/-/wikis/home).\n",
    "    \n",
    "> We highlight that, using the PARSEC22 SSP libraries come with the advantage of reducing the total amount of \n",
    "> computations the code has to perform for getting to a final equivalent SED.\n",
    "> Namely, using our custom SSP libraries avoids the need to compute separately the radio stellar emissions.\n",
    "> Furthermore, nebular line emission is currently only available with the PARSEC22 SSP libraries.}\n",
    "\n",
    "### General to all the libraries in the family:\n",
    "\n",
    "|           | size | values |\n",
    "|-----------|-----|-----------------------------|\n",
    "| **Time grid** | 146 | irregularly distributed values in the interval $0 < \\tau < 1.4\\times10^{10}\\ \\text{years}$ |\n",
    "| **Metallicity grid** | 7 | ``[1.e-04, 5.e-04, 1.e-03, 4.e-03, 8.e-03, 2.e-02, 1.e+00]`` |\n",
    "\n",
    "### Specific to the different flavours:\n",
    "\n",
    "* ``NT`` : stellar atmospheres + Non-Thermal emission (i.e. Super-Nova Synchrotron)\n",
    "* ``NTL``: stellar atmospheres + Non-Thermal emission (i.e. Super-Nova Synchrotron) + Nebular emission (i.e. Bremstrhalung + Emission Lines)\n",
    "\n",
    "| flavour | Wavelength grid | Description |\n",
    "|---------|-----------------|-------------|\n",
    "| ``NT`` | 1562 | See above, spanning from $1\\mathring{A}$ to $10^{10}\\mathring{A}$, the missing parts are padded with zeros |\n",
    "| ``NT.refined`` | 2189 | Same as above but the density of wavelength grid points has been increased to have at least 128 values per-each order of magnitude in the wavelength grid |\n",
    "| ``NTL`` | 1562 | Same as ``NT`` but with more resolution in the Optical/UV bands |\n",
    "| ``NTL.refined`` | 2189 | Same as ``NT.refined`` |\n",
    "\n",
    "\n",
    "### Tips\n",
    "\n",
    "* the ``parsec22`` family is our default, it is the most updated and has been developed with a resolution largely sufficient for most of the applications.\n",
    "* currently, using ``parsec22.NTL`` libraries is the only method to account for the nebular contribution both in terms of\n",
    "    - energy balance\n",
    "    - line emission\n",
    "* flavours marked ``refined`` should be used when computing the transmitted band flux in photometric systems including bands from the submm/mm bands and longer (lower frequency).\n",
    "\n",
    "### Some plots\n",
    "\n",
    "Let's check the difference between the ``NT`` and ``NTL`` flavours:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cf55117",
   "metadata": {},
   "outputs": [],
   "source": [
    "l, t, Z, Lp22nt = load_SSP_table('parsec22.NT')\n",
    "Lp22nt = reshape_SSP_table( Lp22nt, shape=(l.size, t.size, Z.size) )\n",
    "*_, Lp22ntl = load_SSP_table('parsec22.NTL')\n",
    "Lp22ntl = reshape_SSP_table( Lp22ntl, shape=(l.size, t.size, Z.size) )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c0bc007",
   "metadata": {},
   "source": [
    "(since the wavelength, age and metallicity grids are all the same for the two cases, we store them just from the first table)\n",
    "\n",
    "We select some different ages from the table picking some indices and check the value from the age-grid ``t``"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a45898a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "tidx = [1, 25, 75, 125]\n",
    "iz = 3\n",
    "print( f'Metallicity: {Z[iz]}')\n",
    "print( 'Ages:')\n",
    "for it in tidx :\n",
    "    print( f'{t[it]:.2e} years' )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bcd128c0",
   "metadata": {},
   "source": [
    "And finally, we plot in 4 different panels the four different ages.\n",
    "Note how the difference is evident just at the youngest ages (solid blue line VS dashed green line in the two upper plots)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "638e9ae2",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(2,2, figsize=(14,8), sharex=True, sharey=True, \n",
    "                        gridspec_kw={'hspace':0.0, 'wspace':0.0}, \n",
    "                        tight_layout=True)\n",
    "\n",
    "for ii, it in enumerate(tidx) :\n",
    "    jj = int(ii/2)\n",
    "    kk = int(ii%2)\n",
    "    if jj != 0 :\n",
    "        axs[jj,kk].set_xlabel('$\\\\lambda\\\\ [\\\\AA]$')\n",
    "    if kk == 0 :\n",
    "        axs[jj,kk].set_ylabel('$L_\\\\lambda\\\\ [L_\\\\odot\\\\;\\\\AA^{-1}]$')\n",
    "    axs[jj,kk].set_xscale('log')\n",
    "    axs[jj,kk].set_yscale('log')\n",
    "    axs[jj,kk].set_xlim(2.*10**0,10**10)\n",
    "    axs[jj,kk].plot(l, Lp22nt[:,it,3], label='non-thermal')\n",
    "    axs[jj,kk].plot(l, Lp22ntl[:,it,3], ls='--', label='non-thermal + nebular')\n",
    "    if ii == 3 :\n",
    "        axs[jj,kk].legend()"
   ]
  }
 ],
 "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
}