{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Warren 2020 Models\n",
    "\n",
    "Data from \"Constraining properties of the next nearby core-collapse supernova with multi-messenger signals: multi-messenger signals\" by Warren, MacKenzie; Couch, Sean; O'Connor, Evan; Morozova, Viktoriya.\n",
    "\n",
    "1D FLASH simulations with STIR, for alpha_lambda = 1.23, 1.25, and 1.27.  Run with SFHo EOS, M1 with 12 energy groups.\n",
    "\n",
    "For more information on these simulations, see Warren, Couch, O'Connor, & Morozova ([arXiv:1912.03328](https://arxiv.org/abs/1912.03328)) and Couch, Warren, & O'Connor (2020).\n",
    "\n",
    "Includes the multi-messenger data from the STIR simulations. The filename indicates the turbulent mixing parameter a and progenitor mass m of the simulation.  Columns are time [s], shock radius [cm], explosion energy [ergs], electron neutrino mean energy [MeV], electron neutrino rms energy [MeV], electron neutrino luminosity [10^51 ergs/s], electron antineutrino mean energy [MeV], electron antineutrino rms energy [MeV], electron antineutrino luminosity [10^51 ergs/s], x neutrino mean energy [MeV], x neutrino rms energy [MeV], x neutrino luminosity [10^51 ergs/s], gravitational wave frequency from eigenmode analysis of the protoneutron star structure [Hz].  Note that the x neutrino luminosity is for one neutrino flavor - to get the total mu/tau neutrino and antineutrino luminosities requires multiplying this number by 4."
   ]
  },
  {
   "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, MassHierarchy\n",
    "from snewpy.models.ccsn import Warren_2020\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 `Warren_2020` 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": [
    "Warren_2020.param"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quite a lot of choice there! Let’s use the `get_param_combinations` function to get a list of all valid combinations and filter it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This will print a long tuple of combinations:\n",
    "#Warren_2020.get_param_combinations()\n",
    "\n",
    "# To get a more manageable list, let’s filter it to show only the\n",
    "# available progenitors with integer masses and mixing of 1.23:\n",
    "filtered_models = list(params for params in Warren_2020.get_param_combinations()\n",
    "                              if params['turbmixing_param'] == 1.23\n",
    "                              and params['progenitor_mass'].value.is_integer()\n",
    "                              and 10 <= params['progenitor_mass'].value <= 25)\n",
    "print(\"Filtered list of progenitors:\\n\", *filtered_models, sep='\\n')\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": [
    "m10 = Warren_2020(**filtered_models[0])\n",
    "m25 = Warren_2020(**filtered_models[-1])\n",
    "# This is equivalent to:\n",
    "#m10 = Warren_2020(progenitor_mass=10*u.solMass, turbmixing_param=1.23)\n",
    "#m25 = Warren_2020(progenitor_mass=25*u.solMass, turbmixing_param=1.23)\n",
    "\n",
    "m10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let’s plot the luminosity of different neutrino flavors for this model. (Note that the `Warren_2020` 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, axes = plt.subplots(1, 2, figsize=(12, 5), sharex=True, sharey=True, tight_layout=True)\n",
    "\n",
    "for i, model in enumerate([m10, m25]):\n",
    "    ax = axes[i]\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",
    "    ax.set(xlim=(-0.05, 0.65),\n",
    "           xlabel=r'$t-t_{\\rm bounce}$ [s]',\n",
    "           title=r'{}: {} $M_\\odot$'.format(model.metadata['EOS'], model.metadata['Progenitor mass'].value))\n",
    "    ax.grid()\n",
    "    ax.legend(loc='upper right', ncol=2, fontsize=18)\n",
    "\n",
    "axes[0].set(ylabel=r'luminosity [foe s$^{-1}$]');"
   ]
  },
  {
   "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()\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 = 50*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",
    "    for flavor in Flavor:\n",
    "        ax.plot(E, spec[flavor],\n",
    "                label=flavor.to_tex(),\n",
    "                color='C0' if flavor.is_electron else 'C1',\n",
    "                ls='-' if flavor.is_neutrino else ':', lw=2,\n",
    "                alpha=0.7)\n",
    "\n",
    "    ax.set(xlabel=r'$E$ [{}]'.format(E.unit),\n",
    "           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 [erg$^{-1}$ s$^{-1}$]')\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.7.3"
  },
  "vscode": {
   "interpreter": {
    "hash": "e2528887d751495e023d57d695389d9a04f4c4d2e5866aaf6dc03a1ed45c573e"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}