Merge pull request #347 from SNEWS2/jpkneller-patch-3

Add `Fischer_2020` model

doc/nb/ccsn/Fischer_2020.ipynb

@@ -0,0 +1,189 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Fischer 2020 Model\n",
+    "\n",
+    "CCSN neutrino model from Tobias Fischer\n",
+    "\n",
+    "The citation is: *Neutrino signal from proto-neutron star evolution: Effects of opacities from charged-current-neutrino interactions and inverse neutron decay*, Fischer, Tobias ; Guo, Gang ; Dzhioev, Alan A. ; Martínez-Pinedo, Gabriel ; Wu, Meng-Ru ; Lohs, Andreas ; Qian, Yong-Zhong, Physical Review C, Volume 101, Issue 2, article id.025804 (https://ui.adsabs.harvard.edu/link_gateway/2020PhRvC.101b5804F/doi:10.1103/PhysRevC.101.025804), [arXiv:1804.10890](https://arxiv.org/abs/1804.10890)."
+   ]
+  },
+  {
+   "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",
+    "from snewpy.neutrino import Flavor\n",
+    "from snewpy.models.ccsn import Fischer_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 `Fischer_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": [
+    "Fischer_2020.param"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In this case, there’s just a single progenitor available; so let’s initialize it. If this is the first time you’re using this progenitor, snewpy will automatically download the required data files for you."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "F2020 = Fischer_2020()\n",
+    "\n",
+    "F2020"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Plot the luminosity of different neutrino flavors for this model. "
+   ]
+  },
+  {
+   "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",
+    "    if flavor.is_electron == True:\n",
+    "        color='C0'\n",
+    "    else:\n",
+    "        color='C1'\n",
+    "    ls='-' if flavor.is_neutrino else ':'\n",
+    "    lw = 2\n",
+    "    \n",
+    "    ax.plot(F2020.time, F2020.luminosity[flavor]/1e51,  # Report luminosity in units foe/s\n",
+    "            label=flavor.to_tex(), color=color, ls=ls, lw=lw)\n",
+    "    \n",
+    "ax.set(xlim=(-0.1, 1), xlabel=r'$t-t_{\\rm bounce}$ [s]')\n",
+    "ax.grid()\n",
+    "ax.legend(loc='upper right', ncol=2, fontsize=18)\n",
+    "ax.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 = F2020.get_initial_spectra(t, E)\n",
+    "ospec_nmo = F2020.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",
+    "        if flavor.is_electron == True:\n",
+    "            color='C0'\n",
+    "        else:\n",
+    "            color='C1'\n",
+    "        ax.plot(E, spec[flavor],\n",
+    "                label=flavor.to_tex(),\n",
+    "                color=color,\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.12.3"
+  },
+  "vscode": {
+   "interpreter": {
+    "hash": "e2528887d751495e023d57d695389d9a04f4c4d2e5866aaf6dc03a1ed45c573e"
+   }
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

python/snewpy/models/ccsn.py

@@ -44,6 +44,18 @@ class method to get a list of all valid combinations and filter it:
 from snewpy.models.registry_model import all_models
 from textwrap import dedent
 
+@legacy_filename_initialization
+@RegistryModel(
+    progenitor_mass= [18 * u.Msun],
+    eos = ['HS(DD2)']
+)
+class Fischer_2020(loaders.Fischer_2020):
+    """Model based on simulations from `Fischer et al. (2020) <https://arxiv.org/abs/1804.10890>`
+    """
+    def __init__(self, progenitor_mass:u.Quantity, eos:str):
+        filename='Fischer_2020.tar.gz'
+        return super().__init__(filename, metadata=self.metadata)
+
 @legacy_filename_initialization
 @RegistryModel(
     progenitor_mass = [13, 20, 30, 50] * u.Msun,

python/snewpy/models/ccsn_loaders.py

@@ -893,3 +893,72 @@ def __init__(self, filename, metadata={}):
         self.filename = os.path.basename(filename)
 
         super().__init__(simtab, metadata)
+
+
+class Fischer_2020(PinchedModel):
+    def __init__(self, filename, metadata={}):
+        """
+        Parameters
+        ----------
+        filename : str
+            Absolute or relative path to file
+        """
+        # Open the requested filename using the model downloader.
+        datafile = self.request_file(filename)
+        self.metadata = metadata
+
+        # Open the requested filename using the model downloader.
+        # datafile = _model_downloader.get_model_data(self.__class__.__name__, filename)
+        # self.filename = os.path.basename(filename)
+
+        simtab = Table()
+
+        tf = tarfile.open(datafile)
+
+        # Open luminosity file
+        with tf.extractfile("luminosity.dat") as Lfile:
+             Ldata = np.genfromtxt(Lfile, skip_header=2)
+
+             simtab['TIME'] = Ldata[:, 0]
+
+             simtab['L_NU_E'] = Ldata[:, 1]
+             simtab['L_NU_E_BAR'] = Ldata[:, 2]
+             simtab['L_NU_X'] = Ldata[:, 3]
+             simtab['L_NU_X_BAR'] = Ldata[:, 4]
+
+             Lfile.close()
+
+        # Open mean energy file
+        with tf.extractfile("menergy.dat") as Efile:
+            Edata = np.genfromtxt(Efile, skip_header=2)
+
+            simtab['E_NU_E'] = Edata[:, 1] << u.MeV
+            simtab['E_NU_E_BAR'] = Edata[:, 2] << u.MeV
+            simtab['E_NU_X'] = Edata[:, 3] << u.MeV
+            simtab['E_NU_X_BAR'] = Edata[:, 4] << u.MeV
+
+            Efile.close()
+
+        # Open rms energy file
+        with tf.extractfile("rmsenergy.dat") as RMSEfile:
+            RMSEdata = np.genfromtxt(RMSEfile, skip_header=2)
+
+            simtab['RMS_NU_E'] = RMSEdata[:, 1] << u.MeV
+            simtab['RMS_NU_E_BAR'] = RMSEdata[:, 2] << u.MeV
+            simtab['RMS_NU_X'] = RMSEdata[:, 3] << u.MeV
+            simtab['RMS_NU_X_BAR'] = RMSEdata[:, 4] << u.MeV
+
+            RMSEfile.close()
+
+        simtab['ALPHA_NU_E'] = (2.0 * simtab['E_NU_E'] ** 2 - simtab['RMS_NU_E'] ** 2) / \
+            (simtab['RMS_NU_E'] ** 2 - simtab['E_NU_E'] ** 2)
+        simtab['ALPHA_NU_E_BAR'] = (2.0 * simtab['E_NU_E_BAR'] ** 2 - simtab['RMS_NU_E_BAR'] ** 2) / \
+            (simtab['RMS_NU_E_BAR'] ** 2 - simtab['E_NU_E_BAR'] ** 2)
+        simtab['ALPHA_NU_X'] = (2.0 * simtab['E_NU_X'] ** 2 - simtab['RMS_NU_X'] ** 2) / \
+            (simtab['RMS_NU_X'] ** 2 - simtab['E_NU_X'] ** 2)
+        simtab['ALPHA_NU_X_BAR'] = (2.0 * simtab['E_NU_X_BAR'] ** 2 - simtab['RMS_NU_X_BAR'] ** 2) / \
+            (simtab['RMS_NU_X_BAR'] ** 2 - simtab['E_NU_X_BAR'] ** 2)
+
+        tf.close()
+
+        super().__init__(simtab, metadata)

python/snewpy/models/model_files.yml

@@ -57,6 +57,9 @@ models:
         Bugli_2021:
             repository: *ccsn_repository
 
+        Fischer_2020:
+            repository: *ccsn_repository
+
     presn:
 
         Odrzywolek_2010:

python/snewpy/test/test_01_registry.py

@@ -9,7 +9,7 @@
                           Sukhbold_2015, Bollig_2016, Walk_2018, \
                           Walk_2019, Fornax_2019, Warren_2020, \
                           Kuroda_2020, Fornax_2021, Zha_2021, \
-                          Fornax_2022, Mori_2023, Bugli_2021
+                          Fornax_2022, Mori_2023, Bugli_2021, Fischer_2020
 
 from astropy import units as u
 
@@ -335,6 +335,24 @@ def test_Fornax_2021(self):
             self.assertEqual(len(f), len(Flavor))
             self.assertEqual(f[Flavor.NU_E].unit, 1/(u.erg * u.s))
 
+    def test_Fischer_2020(self):
+        """
+        Instantiate a set of 'Fischer 2020' models
+        """
+        model = Fischer_2020()
+        self.assertEqual(model.metadata['Progenitor mass'], [18 * u.Msun])
+        self.assertEqual(model.metadata['EOS'], 'HS(DD2)')
+
+        # Check that times are in proper units.
+        t = model.get_time()
+        self.assertTrue(t.unit, u.s)
+
+        # Check that we can compute flux dictionaries.
+        f = model.get_initial_spectra(0*u.s, 10*u.MeV)
+        self.assertEqual(type(f), dict)
+        self.assertEqual(len(f), len(Flavor))
+        self.assertEqual(f[Flavor.NU_E].unit, 1/(u.erg * u.s))
+
     def test_Zha_2021(self):
         """
         Instantiate a set of 'Zha 2021' models

python/snewpy/test/test_02_models.py

@@ -7,7 +7,7 @@
 from snewpy.models.ccsn_loaders import Nakazato_2013, Tamborra_2014, OConnor_2013, OConnor_2015, \
                                   Sukhbold_2015, Bollig_2016, Walk_2018, \
                                   Walk_2019, Fornax_2019, Warren_2020, \
-                                  Kuroda_2020, Fornax_2021, Zha_2021, Bugli_2021
+                                  Kuroda_2020, Fornax_2021, Zha_2021, Bugli_2021, Fischer_2020
 from astropy import units as u
 from snewpy import model_path
 import os
@@ -88,6 +88,27 @@ def test_Bugli_2021(self):
                 self.assertEqual(len(f), len(Flavor))
                 self.assertEqual(f[Flavor.NU_E].unit, 1/(u.erg * u.s))
 
+    def test_Fischer_2020(self):
+        """
+        Instantiate a set of 'Fischer 2020' models
+        """
+        metadata = {
+            'Progenitor mass': [18 * u.Msun],
+            'EOS': ['HS(DD2)'],
+        }
+        mfile = 'Fischer_2020/Fischer_2020.tar.gz'
+        model = Fischer_2020(os.path.join(model_path, mfile), metadata=metadata)
+
+        # Check that times are in proper units.
+        t = model.get_time()
+        self.assertTrue(t.unit, u.s)
+
+        # Check that we can compute flux dictionaries.
+        f = model.get_initial_spectra(0*u.s, 10*u.MeV)
+        self.assertEqual(type(f), dict)
+        self.assertEqual(len(f), len(Flavor))
+        self.assertEqual(f[Flavor.NU_E].unit, 1/(u.erg * u.s))
+
     def test_OConnor_2013(self):
         """
         Instantiate a set of "O'Connor 2015" models