{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Nakazato 2013 Models\n",
    "\n",
    "This notebook includes an example of how to read the .fits files located in same directory as this notebook. These files contain neutrino luminosity, mean energy, and energy spectrum pinch parameter (alpha) as functions of time for the CCSN neutrino flavors. These were obtained using models from Nakazato et al., 2013 and 2015. Data are public and taken from their [website](http://asphwww.ph.noda.tus.ac.jp/snn/).\n",
    "\n",
    "The citation for use of the database is: *Supernova Neutrino Light Curves and Spectra for Various Progenitor Stars: From Core Collapse to Proto-neutron Star Cooling*, K. Nakazato, K. Sumiyoshi, H. Suzuki, T. Totani, H. Umeda, and S. Yamada, [Astrophys. J. Supp. 205 (2013) 2](http://dx.doi.org/10.1088/0067-0049/205/1/2), [arXiv:1210.6841](http://arxiv.org/abs/1210.6841).\n",
    "\n",
    "If the BH model with LS220 EOS is used, the citation is: *Spectrum of the Supernova Relic Neutrino Background and Metallicity Evolution of Galaxies*, K. Nakazato, E. Mochida, Y. Niino, and H. Suzuki, [Astrophys. J. 804 (2015) 75](http://dx.doi.org/10.1088/0004-637X/804/1/75), [arXiv:1503.01236](http://arxiv.org/abs/1503.01236).\n",
    "\n",
    "For the BH model with Togashi EOS, the citation is: *Numerical Study of Stellar Core Collapse and Neutrino Emission Using the Nuclear Equation of State Obtained by the Variational Method*, K. Nakazato, K. Sumiyoshi, and H. Togashi, [Publ. Astron. Soc. Jpn. 73 (2021) 639-651](https://doi.org/10.1093/pasj/psab026), [arXiv:2103.14386](http://arxiv.org/abs/2103.14386).\n",
    "\n",
    "These examples use code taken from IceCube's fast supernova monte carlo, [ASTERIA](https://github.com/IceCubeOpenSource/ASTERIA). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from astropy import units as u \n",
    "\n",
    "from snewpy.neutrino import Flavor, MixingParameters\n",
    "from snewpy.models.ccsn import Nakazato_2013\n",
    "from snewpy.flavor_transformation import NoTransformation, AdiabaticMSW, ThreeFlavorDecoherence\n",
    "\n",
    "mpl.rc('font', size=16)\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initialize Models\n",
    "\n",
    "To start, let’s see what progenitors are available for the `Nakazato_2013` model. We can use the `param` property to view all physics parameters and their possible values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nakazato_2013.param"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quite a lot of choice there! However, for the Nakazato model, not all combinations of these parameters are valid. We can use the `get_param_combinations` function to get a list of all valid combinations or to filter it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This will print a long tuple of combinations:\n",
    "#Nakazato_2013.get_param_combinations()\n",
    "\n",
    "# To get a more manageable list, let’s filter it to show only the\n",
    "# available progenitors with a metallicity of z=0.004:\n",
    "highz_models = list(params for params in Nakazato_2013.get_param_combinations() if params['metallicity'] == 0.02)\n",
    "print(\"Progenitors with metallicity z=0.02:\", *highz_models, sep='\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We’ll pick one of these progenitors and initialise it. If this is the first time you’re using this progenitor, snewpy will automatically download the required data file for you."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_params = highz_models[3]\n",
    "model = Nakazato_2013(**model_params)\n",
    "# This is equivalent to:\n",
    "#model = Nakazato_2013(progenitor_mass=20*u.solMass, revival_time=100*u.ms, metallicity=0.004, eos='shen')\n",
    "model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let’s plot the luminosity of different neutrino flavors for this model. (Note that the `Nakazato_2013` simulations didn’t distinguish between $\\nu_x$ and $\\bar{\\nu}_x$, so both have the same luminosity.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, figsize=(8,6), tight_layout=False)\n",
    "\n",
    "for flavor in Flavor:\n",
    "    ax.plot(model.time, model.luminosity[flavor]/1e51, # Report luminosity in units foe/s\n",
    "            label=flavor.to_tex(),\n",
    "            color = 'C0' if flavor.is_electron else 'C1',\n",
    "            ls = '-' if flavor.is_neutrino else ':',\n",
    "            lw = 2 )\n",
    "\n",
    "ax.set(xlim=(-0.05, 0.5),\n",
    "       xlabel=r'$t-t_{\\rm bounce}$ [s]',\n",
    "       ylabel=r'luminosity [foe s$^{-1}$]')\n",
    "ax.grid()\n",
    "ax.legend(loc='upper right', ncol=2, fontsize=18);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initial and Oscillated Spectra\n",
    "\n",
    "Plot the neutrino spectra at the source and after the requested flavor transformation has been applied.\n",
    "\n",
    "### Adiabatic MSW Flavor Transformation: Normal mass ordering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Adiabatic MSW effect. NMO is used by default.\n",
    "xform_nmo = AdiabaticMSW(MixingParameters())\n",
    "\n",
    "# Energy array and time to compute spectra.\n",
    "# Note that any convenient units can be used and the calculation will remain internally consistent.\n",
    "E = np.linspace(0,100,201) * u.MeV\n",
    "t = 100*u.ms\n",
    "\n",
    "ispec = model.get_initial_spectra(t, E)\n",
    "ospec_nmo = model.get_transformed_spectra(t, E, xform_nmo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1,2, figsize=(12,5), sharex=True, sharey=True, tight_layout=True)\n",
    "\n",
    "for i, spec in enumerate([ispec, ospec_nmo]):\n",
    "    ax = axes[i]\n",
    "    plt.sca(ax)\n",
    "    spec.plot('energy')\n",
    "    \n",
    "    ax.set(title='Initial Spectra: $t = ${:.1f}'.format(t) if i==0 else 'Oscillated Spectra: $t = ${:.1f}'.format(t))\n",
    "    ax.grid()\n",
    "    ax.legend(loc='upper right', ncol=2, fontsize=16)\n",
    "\n",
    "ax = axes[0]\n",
    "ax.set(ylabel=r'flux, MeV')\n",
    "\n",
    "fig.tight_layout();"
   ]
  }
 ],
 "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.13.5"
  },
  "vscode": {
   "interpreter": {
    "hash": "e2528887d751495e023d57d695389d9a04f4c4d2e5866aaf6dc03a1ed45c573e"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}