Add presn.Yoshida_2016 model

doc/nb/presn/Yoshida_2016.ipynb

@@ -0,0 +1,287 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "83d2d170-a5b3-44d9-ad05-762448b08c80",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from astropy import units as u\n",
+    "\n",
+    "#import the snewpy modules\n",
+    "from snewpy.models import presn\n",
+    "from snewpy.flavor_transformation import NoTransformation, AdiabaticMSW, MassHierarchy\n",
+    "from snewpy.neutrino import Flavor"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f2ad86ec-ff7e-4ef6-818e-f244ac8112f4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# set the plot parameters\n",
+    "plt.rc('grid', ls=':')\n",
+    "plt.rc('axes', grid=True)\n",
+    "plt.rc('legend', fontsize=12, loc='upper right')\n",
+    "\n",
+    "# define drawing styles for each flavor\n",
+    "styles = {f: dict(color='C0' if f.is_electron else 'C1',\n",
+    "                  ls='-' if f.is_neutrino else ':',\n",
+    "                  label=f.to_tex()) for f in Flavor}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4f1b6439-6702-417a-a9b6-0ca0d5798f47",
+   "metadata": {},
+   "source": [
+    "## Initialize the model and calculate the flux"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ffff2732-a87a-401c-b85d-a0f0d5598657",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# See what progenitors are available …\n",
+    "print(f\"Available progenitors:\\n{presn.Yoshida_2016.param}\")\n",
+    "\n",
+    "# … and select one of them to initialise the model\n",
+    "model = presn.Yoshida_2016(progenitor_mass=15*u.Msun)\n",
+    "model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4d8e7b15-0a7e-49f3-b6bd-9463e8b371a5",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# 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, 25, 201) * u.MeV\n",
+    "t = np.geomspace(-2*u.day, -1*u.s, 101)\n",
+    "distance = 1*u.kpc"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "efb466d7-1b4c-4949-a983-5134fba3f35b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "flux = model.get_flux(t, E,\n",
+    "                      distance=distance,\n",
+    "                      flavor_xform=NoTransformation())\n",
+    "flux"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f06fc5ed-e39f-4f91-a8c6-09a671a7f6d0",
+   "metadata": {},
+   "source": [
+    "## Plotting the integral neutrino fluence and rates"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "34cd5637-201f-493c-9aa8-10a653d4158b",
+   "metadata": {},
+   "source": [
+    "### Integral neutrino fluence vs. Energy for the last hour before collapse"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4f8f5b53-257d-43b7-b933-415738e8b4a1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#integrate the flux over the last hour before the collapse\n",
+    "fluence = flux.integrate('time', limits=[-1, 0]<<u.hour)\n",
+    "#get the relevant arrays for plotting\n",
+    "x = fluence.energy\n",
+    "y = fluence.array.squeeze()\n",
+    "#plot each flavor\n",
+    "for flv in fluence.flavor:\n",
+    "    plt.plot(x,y[flv], **styles[flv])\n",
+    "#add the legend\n",
+    "plt.legend(ncol=2, loc='lower center')\n",
+    "#adjust the scales\n",
+    "plt.autoscale(tight=True,axis='x')\n",
+    "plt.yscale('log')\n",
+    "plt.ylim(1)\n",
+    "#define the labels\n",
+    "xlabel = 'Neutrino energy'\n",
+    "ylabel = 'Neutrino fluence'\n",
+    "plt.xlabel(f'{xlabel} [{x.unit._repr_latex_()}]')\n",
+    "plt.ylabel(f'{ylabel} [{y.unit._repr_latex_()}]')\n",
+    "#add the plot title\n",
+    "plt.title('Neutrinos in the last 1 hour before collapse')\n",
+    "#add the model parameters\n",
+    "plt.annotate(str(model) + f'\\ndistance : {distance}',\n",
+    "             xy=(0.98,0.98),\n",
+    "             xycoords='axes fraction',\n",
+    "             va='top', ha='right'\n",
+    "            )\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "8347961e",
+   "metadata": {},
+   "source": [
+    "### Integral neutrino rate vs. time"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dba0da41-6061-4c7a-b1fc-4b73f29ea960",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#integrate the flux over all energy bins\n",
+    "rate = flux.integrate('energy')\n",
+    "#get the relevant arrays for plotting\n",
+    "x = rate.time\n",
+    "y = rate.array.squeeze()\n",
+    "#we can also convert to a more convenient unit\n",
+    "y = y.to('1/(hour*m**2)')\n",
+    "#plot each flavor\n",
+    "for flv in rate.flavor:\n",
+    "    plt.plot(x,y[flv], **styles[flv])\n",
+    "#add the legend\n",
+    "plt.legend(ncol=2, loc='lower center')\n",
+    "#adjust the scales\n",
+    "plt.autoscale(tight=True,axis='x')\n",
+    "plt.yscale('log')\n",
+    "#plt.ylim(1)\n",
+    "#define the labels\n",
+    "xlabel = 'Time before collapse'\n",
+    "ylabel = 'Neutrino rate'\n",
+    "plt.xlabel(f'{xlabel} [{x.unit._repr_latex_()}]')\n",
+    "plt.ylabel(f'{ylabel} [{y.unit._repr_latex_()}]')\n",
+    "#add the plot title\n",
+    "plt.title('Integral neutrino rate')\n",
+    "#add the model parameters\n",
+    "plt.annotate(str(model) + f'\\ndistance : {distance}',\n",
+    "             xy=(0.02,0.98),\n",
+    "             xycoords='axes fraction',\n",
+    "             va='top', ha='left'\n",
+    "            )\n",
+    "plt.xscale('symlog')\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "2dfec9e5-2819-4993-b352-819224b65279",
+   "metadata": {},
+   "source": [
+    "## Plot energy spectrum with neutrino oscillations"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "94e4605b-fcee-4d3a-8344-caa4100e3ee1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# 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,25,201) * u.MeV\n",
+    "t = -1*u.s\n",
+    "\n",
+    "xforms = [NoTransformation(), AdiabaticMSW(mh=MassHierarchy.NORMAL), AdiabaticMSW(mh=MassHierarchy.INVERTED)]\n",
+    "labels = ['No transformation', 'AdiabaticMSW NMO', 'AdiabaticMSW IMO']\n",
+    "xform_styles = ['-','--',':']\n",
+    "flavor_colors = ['C0','C1','C2','C3']\n",
+    "\n",
+    "#calculate the flux for different transformations\n",
+    "fluxes  = {xform: model.get_flux(t, E, distance, flavor_xform=xform) for xform in xforms}\n",
+    "#make a big figure\n",
+    "fig, axes = plt.subplots(1,4,figsize=(16,4), sharex=True, sharey=True)\n",
+    "\n",
+    "for flv, ax in zip(Flavor,axes):\n",
+    "    plt.sca(ax)\n",
+    "    for xform, line_style, label in zip(xforms,xform_styles, labels):\n",
+    "        flux = fluxes[xform]\n",
+    "        #get the relevant arrays for plotting\n",
+    "        x = flux.energy\n",
+    "        y = flux.array[flv].squeeze()\n",
+    "        #plot each flavor\n",
+    "        plt.plot(x,y, \n",
+    "                 ls=line_style, \n",
+    "                 color=flavor_colors[flv], \n",
+    "                 label=f'{Flavor(flv).to_tex()} {label}'\n",
+    "                )\n",
+    "        \n",
+    "    #add the legend\n",
+    "    plt.legend(ncol=1, loc='lower left', fontsize=10)\n",
+    "    \n",
+    "    #define the labels\n",
+    "    xlabel = 'Neutrino energy'\n",
+    "    ylabel = 'Neutrino flux'\n",
+    "    plt.xlabel(f'{xlabel} [{x.unit._repr_latex_()}]')\n",
+    "    #add the plot title\n",
+    "\n",
+    "    #add the model parameters\n",
+    "    plt.annotate(str(model) + f'\\ndistance : {distance}',\n",
+    "                 xy=(0.98,0.98),\n",
+    "                 xycoords='axes fraction',\n",
+    "                 va='top', ha='right'\n",
+    "                )\n",
+    "    axes[0].set_ylabel(f'{ylabel} [{y.unit._repr_latex_()}]')\n",
+    "    #adjust the scales\n",
+    "    plt.autoscale(tight=True,axis='x')\n",
+    "    plt.autoscale(tight=False,axis='y')\n",
+    "    plt.yscale('log')\n",
+    "    plt.ylim(1)\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d12483ed-e5f7-48a0-99a7-9f7f1bfdd289",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "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.11.8"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

python/snewpy/models/model_files.yml

@@ -5,7 +5,7 @@
 config:
     - &snewpy "https://github.com/SNEWS2/snewpy/raw/v{snewpy_version}/models/{model}/{filename}"
     - &ccsn_repository "https://github.com/SNEWS2/snewpy-models-ccsn/raw/v0.2/models/{model}/{filename}"
-    - &presn_repository "https://github.com/SNEWS2/snewpy-models-presn/raw/v0.1/models/{model}/{filename}"
+    - &presn_repository "https://github.com/SNEWS2/snewpy-models-presn/raw/v0.2/models/{model}/{filename}"
 
 models:
     ccsn:
@@ -65,3 +65,5 @@ models:
         Kato_2017:
             repository: *presn_repository
 
+        Yoshida_2016:
+            repository: *presn_repository

python/snewpy/models/presn.py

@@ -45,3 +45,17 @@ class Kato_2017(loaders.Kato_2017):
     def __init__(self, progenitor_mass:u.Quantity):
         path=f"pre_collapse/m{progenitor_mass.to_value('Msun'):.0f}"
         super().__init__(path)
+
+
+@RegistryModel(
+    progenitor_mass = [12, 15, 20]<<u.Msun
+)
+class Yoshida_2016(loaders.Yoshida_2016):
+    """Presupernova model based on
+    `Yoshida et al. (2016), PRD 93, 123012 <https://doi.org/10.1103/PhysRevD.93.123012>`_
+
+    Dataset available on `Zenodo <https://zenodo.org/records/3778014>`__
+    """
+    def __init__(self, progenitor_mass:u.Quantity):
+        path=f"t_spc_m{progenitor_mass.to_value('Msun'):.0f}_1.2.txt"
+        super().__init__(path)

python/snewpy/models/presn_loaders.py

@@ -139,3 +139,35 @@ def get_initial_spectra(self, t, E, flavors=Flavor):
         flux = self.interpolated(t, E) / (u.MeV * u.s)
         return {f: flux[f] for f in flavors}
 
+class Yoshida_2016(SupernovaModel):
+    """Set up a presupernova model based on 
+    [Yoshida et al. (2016), PRD 93, 123012.]
+    """
+    def __init__(self, filename):
+        with open(self.request_file(filename)) as f:
+            data = []
+            T = []
+            while (line := f.readline()):
+                if not line: break
+                T += [float(line.split()[1])]
+                data += [[np.loadtxt(f, max_rows=100).flatten() for i in range(4)]]
+        times = np.array(T)
+        energies = np.concatenate([
+                np.linspace(0,10,1001)[1:],
+                np.linspace(10,20,501)[1:]
+            ])
+        dNdEdT = np.stack(data, axis=1)
+        #rearrange flavors from NU_E, NU_E_BAR,_NU_X, NU_X_BAR to current 
+        indices = np.argsort([Flavor.NU_E,Flavor.NU_E_BAR, Flavor.NU_X, Flavor.NU_X_BAR]) 
+        dNdEdT = dNdEdT.take(indices, axis=0)
+
+        self.interpolated = _interp_TE(
+            times, energies, dNdEdT, ax_t=1, ax_e=2
+        )
+        super().__init__(times << u.s, self.metadata)
+
+    def get_initial_spectra(self, t, E, flavors=Flavor):
+        t = np.array(-t.to_value("s"), ndmin=1)
+        E = np.array(E.to_value("MeV"), ndmin=1)
+        flux = self.interpolated(t, E) / (u.MeV * u.s)
+        return {f: flux[f] for f in flavors}

python/snewpy/test/test_presn_rates.py

@@ -14,7 +14,7 @@
 T = np.geomspace(-1*u.hour, -1*u.min,1000)
 E = np.linspace(0,20,100)*u.MeV
 
-@pytest.mark.parametrize('model_class',[presn.Odrzywolek_2010, presn.Kato_2017, presn.Patton_2017])
+@pytest.mark.parametrize('model_class',[presn.Odrzywolek_2010, presn.Kato_2017, presn.Patton_2017, presn.Yoshida_2016])
 @pytest.mark.parametrize('transformation',[AdiabaticMSW(mh=mh) for mh in MassHierarchy])
 @pytest.mark.parametrize('detector', ["wc100kt30prct"])
 def test_presn_rate(model_class, transformation, detector):