{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e7e545bb",
   "metadata": {},
   "source": [
    "# Build a galaxy model from an interpolated SFH\n",
    "\n",
    "GalaPy provides a non-parametric, interpolated, step-wise SFH model with derived components (like dust/gas mass and metallicity) treated as free parameters.\n",
    "This model is designed to predict the emission from galaxies for which the stellar mass growth history is available (e.g. obtained from hydro-dynamical simulations or with semi-analytical models) or to test the behaviour of exotic and arbitrarily complex SFH shapes.\n",
    "\n",
    "## Load dataset\n",
    "\n",
    "The dataset at hand might look something like"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "monetary-miami",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "tau = numpy.array([\n",
    "    0.0000000e+00, 2.0674900e+08, 4.2752600e+08, 6.6250900e+08,\n",
    "    9.1218400e+08, 1.1766550e+09, 1.4557920e+09, 1.7499790e+09,\n",
    "    2.0585520e+09, 2.3815860e+09, 2.7187860e+09, 3.0688570e+09,\n",
    "    3.4325970e+09, 3.8068280e+09, 4.1934690e+09, 4.5884250e+09,\n",
    "    4.9938270e+09, 5.4046790e+09, 5.8241890e+09, 6.2453270e+09,\n",
    "    6.6728250e+09, 7.0994560e+09, 7.5283510e+09, 7.9562340e+09,\n",
    "    8.3802740e+09, 8.8046430e+09, 9.2189360e+09, 9.6296370e+09,\n",
    "    1.0036022e+10, 1.0428195e+10, 1.0814397e+10, 1.1194359e+10,\n",
    "    1.1555720e+10\n",
    "])\n",
    "sfr = numpy.array([\n",
    "    0.        , 0.        , 0.00258615, 0.04669582, 0.21555495,\n",
    "    0.56162186, 0.60869537, 0.55047424, 0.36031803, 0.27614234,\n",
    "    0.20609573, 0.16298589, 0.2459535 , 0.47935471, 0.57768621,\n",
    "    0.55990713, 0.52003312, 0.48266521, 0.45552792, 0.43251383,\n",
    "    0.40347906, 0.36103059, 0.28409527, 0.22353781, 0.17772765,\n",
    "    0.15988216, 0.1646592 , 0.17408446, 0.19488885, 0.23297054,\n",
    "    0.2524771 , 0.25091753, 0.23612502\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af0aeb2a",
   "metadata": {},
   "source": [
    "where ``tau`` is an array of sampled ages (in units of $[\\text{yr}]$) in the evolution history of the galaxy and ``sfr`` is the star formation rate, $\\psi(\\tau)$ (in units $[M_\\odot/\\text{yr}]$) corresponding to each value in ``ŧau``.\n",
    "\n",
    "This particular data-set has been extracted from a catalogue of galaxies simulated with the semi-analytical model L-GALAXIES2020 in [Parente et al. (2023)](https://doi.org/10.1093/mnras/stad907) in which the autors include new features to the code, namely: (i) a state-of-the-art dust model that adopts the two-size approximation and (ii) a new disc instability criterion that triggers bulge and central black hole growth.\n",
    "\n",
    "Let's fix some further properties of the object simulated by the SAM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0502e5fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# redshift of the source\n",
    "zz = 1.e-4\n",
    "# current age of the object\n",
    "age = tau.max()\n",
    "# average absolute metallicity\n",
    "Zgxy = 6.74e-3\n",
    "# dust mass\n",
    "Mdust = 5.e+6"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0f4218d",
   "metadata": {},
   "source": [
    "## Build galaxy model\n",
    "\n",
    "From the information above we can build a galaxy model with GalaPy, let's first import the class wrapping up all the physical components necessary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "676bffd6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from galapy.Galaxy import GXY"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82011600",
   "metadata": {},
   "source": [
    "We choose a ``parsec22.NTL`` SSP library that along with the emission from stellar atmospheres also includes (i) non-thermal emission due to synchrotron from super-novae and, (ii) nebular thermal emission including lines  (see [Choose the SSP library](https://galapy.readthedocs.io/en/latest/notebooks/choose_ssp_lib.html) for further details).\n",
    "\n",
    "We set the `'interpolated'` model by passing the relevant keyword to the `sfh` argument dictionary and also set some further properties as derived from the SAM and listed above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92e4cdeb",
   "metadata": {},
   "outputs": [],
   "source": [
    "gxy = GXY( \n",
    "    age = age, redshift = zz, \n",
    "    csp = {'ssp_lib':'parsec22.NTL'}, # set the SSP library\n",
    "    sfh = { \n",
    "        'model':'interpolated',       # choose the interpolated model\n",
    "        'tau':tau, 'sfr':sfr,         # pass the dataset simulated with SAM\n",
    "        'Zgxy':Zgxy, 'Mdust':Mdust    # pass the further properties from the SAM\n",
    "    },\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a004f15",
   "metadata": {},
   "source": [
    "The galaxy model above is sufficient to have a simulated emission from the object (and of course further physical properties can be tuned via the free-parameters of the model (e.g. properties of dust), for the full list of tunable free-parameters the user can refer to the following page in this docs: [Free parameters of the galaxy model](https://galapy.readthedocs.io/en/latest/general/free_parameters.html)\n",
    "\n",
    "## Check the interpolation goes as expected\n",
    "\n",
    "Let's check what we are modelling by comparing the input SFH with the interpolated value returned by `gxy.sfh`.\n",
    "First, we can import the ``matplotlib.pyplot`` formatted version present in the ``galapy.analysis.plot`` submodule and some other function for cosmetics:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab543d97",
   "metadata": {},
   "outputs": [],
   "source": [
    "from galapy.analysis.plot import plt, format_axes_ticks"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e84648aa",
   "metadata": {},
   "source": [
    "we also want to model the SFH on a time-grid thinner than the original one, but still containing the original values, we can build one as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f5fcc73",
   "metadata": {},
   "outputs": [],
   "source": [
    "tt = numpy.sort(numpy.append( \n",
    "    tau, # original array\n",
    "    numpy.logspace( 8.0, numpy.log10(age), 256 ) # log-spaced grid of times\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e785830d",
   "metadata": {},
   "source": [
    "The ``tt`` array goes from $10^8\\ \\text{yr}$ up to the age of the object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d28f9fa8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# make a grid of subplots\n",
    "fig, axs = plt.subplots(2,1, sharex=True, \n",
    "                        gridspec_kw={'hspace':0.0}, \n",
    "                        tight_layout=True)\n",
    "\n",
    "#########################################\n",
    "# first sub-plot\n",
    "\n",
    "# fix logarithmic scale for the axes and y-axis label\n",
    "axs[0].set( \n",
    "    #xscale='log', \n",
    "    yscale='log', \n",
    "    ylabel='$\\\\Psi(\\\\tau)\\\\ [M_\\\\odot\\\\;\\\\mathrm{yr}^{-1}]$'\n",
    ")\n",
    "\n",
    "# plot the interpolated model:\n",
    "axs[0].plot(tt, gxy.sfh(tt), \n",
    "            ls='--', color='grey') # cosmetics\n",
    "\n",
    "# plot the SAM input values\n",
    "axs[0].plot(tau, sfr, \n",
    "            ls='none', marker='o', # from here just cosmetics\n",
    "            markerfacecolor='white', \n",
    "            markersize=6, \n",
    "            markeredgewidth=1.75)\n",
    "\n",
    "#########################################\n",
    "# second sub-plot\n",
    "\n",
    "# fix logarithmic scale for the axes and x- and y-axis labels\n",
    "axs[1].set( \n",
    "    #xscale='log', \n",
    "    yscale='log', \n",
    "    xlabel='$\\\\tau\\\\ [\\\\mathrm{yr}]$',\n",
    "    ylabel='$M_\\\\star(\\\\tau)\\\\ [M_\\\\odot]$'\n",
    ")\n",
    "\n",
    "# also plot the evolution of the stellar mass\n",
    "axs[1].plot(tt, gxy.sfh.Mstar(tt, 1000), color='tab:red')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e345a3fe",
   "metadata": {},
   "source": [
    "We can compute the stellar mass of the object at current age with"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ce20603",
   "metadata": {},
   "outputs": [],
   "source": [
    "print( f'{gxy.sfh.Mstar(age):.3e} Msol' )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e3ac447",
   "metadata": {},
   "source": [
    "## Compute emission\n",
    "\n",
    "We can use the ``gxy`` object to simulate a grid of fluxes received from the mock source"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82f63d9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# array of rest-frame wavelengths\n",
    "wave = gxy.wl()\n",
    "# array of fluxes in mJy\n",
    "flux = gxy.get_SED()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74280854",
   "metadata": {},
   "source": [
    "By calling the ``gxy.get_SED()`` method, also the contribution from the different components of the galaxy model is computed and stored internally, it can be accessed from the class attribute ``gxy.components`` which, though, returns luminosities, or by calling the dedicated method to convert it into fluxes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17b57e87",
   "metadata": {},
   "outputs": [],
   "source": [
    "components = gxy.components_to_flux()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c11f5fb0",
   "metadata": {},
   "source": [
    "> Note that the ``components`` object is a dictionary where the keys point out the component short name and\n",
    "> the values are arrays of fluxes in $mJy$ computed on the same grid of wavelengths returned by the ``gxy.wl()`` \n",
    "> method above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f94d27c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "list(components.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "883fa6f8",
   "metadata": {},
   "source": [
    "GalaPy provides some functions for formatted plotting of the SED, they are made available in the ``galapy.analysis.plot`` sub-module:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef52f8f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "from galapy.analysis.plot import sed_layout, sed_components, sed_flux"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8861be0",
   "metadata": {},
   "source": [
    "We can therefore plot the overall SED and components by calling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2882dae6",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,1,figsize=(6,3), constrained_layout=True)\n",
    "\n",
    "# prepares the plot axes\n",
    "ax = sed_layout(gxy.redshift, frame='rest', ax = ax, xlim=(50., 8.e+9), ylim=(1.e-6,1.e+10))\n",
    "\n",
    "# plot the different components\n",
    "_ = sed_components(\n",
    "    wave, components, \n",
    "    redshift=gxy.redshift, frame='rest', \n",
    "    ax=ax\n",
    ")\n",
    "\n",
    "# plot the total flux:\n",
    "_ = sed_flux(\n",
    "    wave, flux, \n",
    "    redshift=gxy.redshift, frame='rest',\n",
    "    ax=ax\n",
    ")\n",
    "\n",
    "# plot a legend\n",
    "ax.legend(ncols=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73e8b68b",
   "metadata": {},
   "source": [
    "## Compute photometry\n",
    "\n",
    "We might want to check what the emission from the SED model above looks like when transmitted through a set of bandpass filters.\n",
    "\n",
    "First let's select some filters from the database, let's say we want the SDSS filters, Herschel's PACS filters and ALMA bands we can extract the correct names and list them by calling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5071d5bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "from galapy.PhotometricSystem import list_filters\n",
    "bands = list_filters('SDSS')+list_filters('Herschel.PACS')+list_filters('ALMA')\n",
    "print(bands)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b20c6d8",
   "metadata": {},
   "source": [
    "We can then build a photometric-system object by instantiating "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb05c957",
   "metadata": {},
   "outputs": [],
   "source": [
    "from galapy.PhotometricSystem import PMS\n",
    "pms = PMS(*bands)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed2c6764",
   "metadata": {},
   "source": [
    "The pivot wavelength associated with each band is contained in the ``lpivot`` attribute:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9162710",
   "metadata": {},
   "outputs": [],
   "source": [
    "lpiv = pms.lpiv\n",
    "lpiv"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "301520f1",
   "metadata": {},
   "source": [
    "which can be plotted using"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f2bdc28",
   "metadata": {},
   "outputs": [],
   "source": [
    "from galapy.analysis.plot import photometric_system as pms_plot\n",
    "fig, ax = plt.subplots(1,1,figsize=(12,3), constrained_layout=True)\n",
    "_ = pms_plot(pms, ax=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a0de0fd",
   "metadata": {},
   "source": [
    "To see the emission in bands we will have to build a ``galapy.Galaxy.PhotoGXY`` model instead of the ``galapy.Galaxy.GXY`` model we have built above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fed3b4a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from galapy.Galaxy import PhotoGXY"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fe10d00",
   "metadata": {},
   "source": [
    "Note that ``PhotoGXY`` inherits from ``GXY``, therefore all the functionalities of the ``gxy`` object built in previous sections of this tutorial.\n",
    "\n",
    "Equivalently to what done in the previous sections, we can build a *photo-galaxy* object by calling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94c718f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "pgxy = PhotoGXY( \n",
    "    pms = pms, # <--- this is the only argument different from the base GXY object \n",
    "    age = age, redshift = zz, \n",
    "    csp = {'ssp_lib':'parsec22.NTL.refined'}, # set the SSP library (in its 'refined' version)\n",
    "    sfh = { \n",
    "        'model':'interpolated',       # choose the interpolated model\n",
    "        'tau':tau, 'sfr':sfr,         # pass the dataset simulated with SAM\n",
    "        'Zgxy':Zgxy, 'Mdust':Mdust    # pass the further properties from the SAM\n",
    "    },\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07f411ec",
   "metadata": {},
   "source": [
    "> **Note** that, since the ALMA bands act on a region of the wavelength space that is otherwise under-sampled,\n",
    "> we have also changed the SSP library, using in this case the refined version of the same model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecba13df",
   "metadata": {},
   "source": [
    "We can compute the transmitted flux in the different bands with"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "067acb7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "pflux = pgxy.photoSED()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58d15f49",
   "metadata": {},
   "source": [
    "and we can plot it together with the full SED flux, being carefull to use observed wavelength this time, and not the rest-frame one (this can be done by calling the function returning the wavelength grid with the argument ``pgxy.wl(obs=True)`` and setting the frame of the plotting functions from ``'rest'`` to either ``'obs'`` or ``'both'`` (as done in the plot below)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7256cb37",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,1,figsize=(6,3), constrained_layout=True)\n",
    "\n",
    "# prepares the plot axes\n",
    "ax = sed_layout(pgxy.redshift, frame='both', ax = ax, xlim=(50., 8.e+9), ylim=(1.e-6,1.e+10))\n",
    "\n",
    "# plot the total flux:\n",
    "_ = sed_flux(\n",
    "    pgxy.wl(obs=True), pgxy.get_SED(), \n",
    "    redshift=pgxy.redshift, frame='both',\n",
    "    ax=ax\n",
    ")\n",
    "\n",
    "# plot the band fluxes here\n",
    "plt.scatter( \n",
    "    lpiv, pflux, \n",
    "    color=plt.cm.plasma(numpy.linspace(0.1,0.9,lpiv.size)), \n",
    "    zorder=2 # otherwise pyplot would put it automatically behind the solid line\n",
    ")\n",
    "\n",
    "# plot a legend\n",
    "ax.legend(ncols=2)"
   ]
  }
 ],
 "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
}