{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "04bf2ca6",
   "metadata": {},
   "source": [
    "# Mori 2023 2D Models with Axionlike Production\n",
    "\n",
    "Neutrino spectra from a set of 2D simulations with axionlike particle production. The models are described in the paper [Multi-messenger signals of heavy axionlike particles in core-collapse supernovae: two-dimensional simulations](https://arxiv.org/abs/2304.11360) by K. Mori, T. Takiwaki, K. Kotake and S. Horiuchi, Phys. Rev. D 108:063027, 2023.\n",
    "\n",
    "The following models are supported:\n",
    "\n",
    "| Model | Axion mass [MeV] | Coupling g10 | tpb,2000 [ms] | Ediag [1e51 erg] | MPNS [Msun] |\n",
    "| ----- | ---------------- | ------------ | ------------- | ---------------- | ---------- |\n",
    "| Standard | − | 0 | 390 | 0.40 | 1.78 |\n",
    "| (100, 2) | 100 | 2 | 385 | 0.37 | 1.77 |\n",
    "| (100, 4) | 100 | 4 | 362 | 0.34 | 1.76 |\n",
    "| (100, 10) | 100 | 10 | 395 | 0.36 | 1.77 |\n",
    "| (100, 12) | 100 | 12 | 357 | 0.43 | 1.77 |\n",
    "| (100, 14) | 100 | 14 | 360 | 0.44 | 1.77 |\n",
    "| (100, 16) | 100 | 16 | 367 | 0.51 | 1.77 |\n",
    "| (100, 20) | 100 | 20 | 330 | 1.10 | 1.74 |\n",
    "| (200, 2) | 200 | 2 | 374 | 0.45 | 1.77 |\n",
    "| (200, 4) | 200 | 4 | 376 | 0.45 | 1.76 |\n",
    "| (200, 6) | 200 | 6 | 333 | 0.54 | 1.75 |\n",
    "| (200, 8) | 200 | 8 | 323 | 0.94 | 1.74 |\n",
    "| (200, 10) | 200 | 10 | 319 | 1.61 | 1.73 |\n",
    "| (200, 20) | 200 | 20 | 248 | 3.87 | 1.62 |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29ceae1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from snewpy.neutrino import Flavor\n",
    "from snewpy.models.ccsn import Mori_2023\n",
    "\n",
    "from astropy import units as u\n",
    "\n",
    "from scipy.interpolate import PchipInterpolator\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b7566ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "mpl.rc('font', size=16)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66d13225",
   "metadata": {},
   "source": [
    "## Initialize the 2D models\n",
    "\n",
    "Use the `param` property of the model class to see the available parameters. Models are initialized using the `axion_mass` and `axion_coupling` parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f62e690f",
   "metadata": {},
   "outputs": [],
   "source": [
    "Mori_2023.param"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55479166",
   "metadata": {},
   "source": [
    "The model with `axion_mass=0` and `axion_coupling=0` is a standard simulation with no ALP production."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "649f5ae4",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_std = Mori_2023(axion_mass=0, axion_coupling=0)\n",
    "model_std.metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f15de7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize a handful of axion models.\n",
    "models = {}\n",
    "for (am, ac) in ((100*u.MeV, 2e-10/u.GeV), (200*u.MeV,2e-10/u.GeV), (100*u.MeV,10e-10/u.GeV), (100*u.MeV,20e-10/u.GeV), (200*u.MeV,10e-10/u.GeV), (200*u.MeV,20e-10/u.GeV)):\n",
    "    models[(am,ac)] = Mori_2023(axion_mass=am, axion_coupling=ac)\n",
    "\n",
    "models"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3cab5a9c",
   "metadata": {},
   "source": [
    "## Plot Model Luminosities\n",
    "\n",
    "Compare axion model luminosity to the standard 2D simulation.\n",
    "\n",
    "Higher mass models with stronger coupling constants should produce a decrease in neutrino luminosity at all flavors relative to the reference simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a398333",
   "metadata": {},
   "outputs": [],
   "source": [
    "for (m,c), model in models.items():\n",
    "    \n",
    "    fig, axes = plt.subplots(2, 4, figsize=(28, 8), sharex=True,\n",
    "                         gridspec_kw={'height_ratios':[3.2,1], 'hspace':0})\n",
    "    \n",
    "    Lmin,  Lmax  = 1e99, -1e99\n",
    "    dLmin, dLmax = 1e99, -1e99\n",
    "    \n",
    "    for j, (flavor) in enumerate([Flavor.NU_E, Flavor.NU_E_BAR, Flavor.NU_MU, Flavor.NU_MU_BAR]):\n",
    "        ax = axes[0][j]\n",
    "        \n",
    "        ax.plot(model_std.time, model_std.luminosity[flavor]/1e51, 'k', label=r'$20M_\\odot$ reference')  # Report luminosity in [foe/s]\n",
    "        Lmin = np.minimum(Lmin, np.min(model_std.luminosity[flavor].to_value('1e51 erg/s')))\n",
    "        Lmax = np.maximum(Lmax, np.max(model_std.luminosity[flavor].to_value('1e51 erg/s')))\n",
    "        \n",
    "        modlabel = rf\"{flavor.to_tex()}: $m_a=${m.to_string(format='latex_inline')}\\n     $g_{{a\\gamma}}=${c.to_string(format='latex_inline')}\"\n",
    "        ax.plot(model.time, model.luminosity[flavor]/1e51,  # Report luminosity in units foe/s\n",
    "                label=modlabel,\n",
    "                color='C0' if flavor.is_electron else 'C1',\n",
    "                ls='-' if flavor.is_neutrino else ':',\n",
    "                lw=2)\n",
    "        if j==0:\n",
    "            ax.set(ylabel=r'luminosity [$10^{51}$ erg s$^{-1}$]')\n",
    "        \n",
    "        ax.legend(fontsize=12)\n",
    "        ax.set(xlim=(model_std.time[0].to_value('s'), model_std.time[-1].to_value('s')))\n",
    "        \n",
    "        ax = axes[1][j]\n",
    "        tmin = np.maximum(model.time[0], model_std.time[0]).to_value('s')\n",
    "        tmax = np.minimum(model.time[-1], model_std.time[-1]).to_value('s')\n",
    "        times = np.arange(tmin, tmax, 0.001)*u.s\n",
    "\n",
    "        Lstd = PchipInterpolator(model_std.time, model_std.luminosity[flavor].to_value('1e51 erg/s'))\n",
    "        Lstd_t = Lstd(times)\n",
    "        select = Lstd_t != 0\n",
    "        \n",
    "        Lmod = PchipInterpolator(model.time, model.luminosity[flavor].to_value('1e51 erg/s'))\n",
    "        Lmod_t = Lmod(times)\n",
    "        dL = (Lmod_t[select] - Lstd_t[select]) / Lstd_t[select]\n",
    "        \n",
    "        dLmin = np.minimum(dLmin, np.min(dL))\n",
    "        dLmax = np.maximum(dLmax, np.max(dL))\n",
    "\n",
    "        ax.plot(times[select], dL)\n",
    "        if j==0:\n",
    "            ax.set(xlabel='time [s]',\n",
    "                   ylabel=r'$\\Delta L_\\nu/L_\\nu$')\n",
    "            \n",
    "    for j in range(4):\n",
    "        axes[0][j].set(ylim=(Lmin, 1.1*Lmax))\n",
    "        axes[1][j].set(ylim=(dLmin, dLmax))\n",
    "        \n",
    "    fig.suptitle(rf\"Axionlike model: $m_a=${m.to_string(format='latex_inline')}, $g_{{a\\gamma}}=${c.to_string(format='latex_inline')}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "snewpy",
   "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}