Mori axion 2023

python/snewpy/models/ccsn.py

@@ -742,6 +742,51 @@ def __new__(cls, filename=None, *, progenitor_mass=None):
     # Populate Docstring with abbreviated param values
     __new__.__doc__ = __new__.__doc__.format(**_param_abbrv)
 
+
+class Mori_2023(_RegistryModel):
+    """Model based on 2D simulations with axionlike particles, K. Mori, T.  Takiwaki, K. Kotake and S. Horiuchi, Phys. Rev. D 108:063027, 2023. All models are based on the non-rotating 20 M_sun solar metallicity progenitor model from S.E. Woolsey and A. Heger, Phys. Rep. 442:269, 2007. Data from private communication.
+        """
+    param = {'axion_mass' : [0, 100, 200],
+             'axion_coupling' : [0, 2, 4, 6, 8, 10, 12, 14, 16, 20] }
+
+    _param_validator = lambda p: \
+        (p['axion_mass'] == 0 and p['axion_coupling'] == 0) or \
+        (p['axion_mass'] == 100 and p['axion_coupling'] in (2,4,10,12,14,16,20)) or \
+        (p['axion_mass'] == 200 and p['axion_coupling'] in (2,4,6,8,10,20))
+
+    _param_abbrv = {
+        'axion_mass': '[0, 100, 200] MeV',
+        'axion_coupling' : '[0..20] 1e-10/GeV' }
+
+    def __new__(cls, axion_mass, axion_coupling):
+        """Model Initialization.
+
+        Parameters
+        ----------
+        axion_mass: int
+            Axion mass in units of MeV. Valid values are {axion_mass}.
+        axion_coupling: int
+            Axion-photon coupling, in units of 1e-10/GeV. Valid values are {axion_coupling}.
+        """
+        # Load from Parameters
+        cls.check_valid_params(cls, axion_mass=axion_mass, axion_coupling=axion_coupling)
+
+        if axion_mass == 0:
+            filename = 't-prof_std.dat'
+        else:
+            filename = f't-prof_{axion_mass}_{axion_coupling}.dat'
+
+        metadata = {
+            'Axion mass': axion_mass * u.MeV,
+            'Axion coupling': axion_coupling * 1e-10/u.GeV,
+            'Progenitor mass': 20*u.Msun }
+
+        return loaders.Mori_2023(filename, metadata)
+
+    # Populate Docstring with abbreviated param values
+    __new__.__doc__ = __new__.__doc__.format(**_param_abbrv)
+
+
 class SNOwGLoBES:
     """A model that does not inherit from SupernovaModel (yet) and imports a group of SNOwGLoBES files."""
 

python/snewpy/models/loaders.py

@@ -702,3 +702,83 @@ def get_initial_spectra(self, t, E, flavors=Flavor, interpolation='linear'):
                 raise ValueError('Unrecognized interpolation type "{}"'.format(interpolation))
 
         return initialspectra
+
+
+class Mori_2023(PinchedModel):
+    def __init__(self, filename, metadata={}):
+        """
+        Parameters
+        ----------
+        filename : str
+            Absolute or relative path to file prefix.
+        """
+        # PNS mass table, from Mori+ 2023.
+        mpns = { (0, 0):    1.78,
+                 (100, 2):  1.77,
+                 (100, 4):  1.76,
+                 (100, 10): 1.77,
+                 (100, 12): 1.77,
+                 (100, 14): 1.77,
+                 (100, 16): 1.77,
+                 (100, 20): 1.74,
+                 (200, 2):  1.77,
+                 (200, 4):  1.76,
+                 (200, 6):  1.75,
+                 (200, 8):  1.74,
+                 (200, 10): 1.73,
+                 (200, 20): 1.62 }
+
+        # Set up model metadata
+        self.axion_mass = metadata['Axion mass']
+        self.axion_coupling = metadata['Axion coupling']
+        self.progenitor_mass = metadata['Progenitor mass']
+
+        self.metadata = metadata
+
+        filebase = os.path.splitext(os.path.basename(filename))[0]
+        if 'std' in filebase:
+            am, ac = 0, 0
+        else:
+            am, ac = [int(_) for _ in filebase.split('_')[1:]]
+
+        self.pns_mass = mpns[(am,ac)] * u.Msun
+        self.metadata['PNS mass'] = self.pns_mass
+
+        # Open the requested filename using the model downloader.
+        datafile = _model_downloader.get_model_data(self.__class__.__name__, filename)
+
+        # Read ASCII data.
+        simtab = Table.read(datafile, format='ascii')
+
+        # Remove the first table row, which appears to have zero input.
+        simtab = simtab[simtab['1:t_sim[s]'] > 0]
+
+        # Get grid of model times.
+        simtab['TIME'] = simtab['2:t_pb[s]'] << u.s
+        for j, (f, fkey) in enumerate(zip([Flavor.NU_E, Flavor.NU_E_BAR, Flavor.NU_X], 'ebx')):
+            simtab[f'L_{f.name}'] = simtab[f'{6+j}:Le{fkey}[e/s]'] << u.erg / u.s
+            # Compute the pinch parameter from E_rms and E_avg
+            # <E^2> / <E>^2 = (2+a)/(1+a), where
+            # E_rms^2 = <E^2> - <E>^2.
+            Eavg = simtab[f'{9+j}:Em{fkey}[MeV]']
+            Erms = simtab[f'{12+j}:Er{fkey}[MeV]']
+            E2 = Erms**2 + Eavg**2
+            x = E2 / Eavg**2
+            alpha = (2-x) / (x-1)
+
+            simtab[f'E_{f.name}'] = Eavg << u.MeV
+            simtab[f'E2_{f.name}'] = E2 << u.MeV**2
+            simtab[f'ALPHA_{f.name}'] = alpha
+
+#            simtab[f'E_{f.name}'] = simtab[f'{9+j}:Em{fkey}[MeV]'] << u.MeV
+#            Erms = simtab[f'{12+j}:Er{fkey}[MeV]'] * u.MeV
+#
+#            # Compute the pinch parameter from E_rms and E_avg
+#            simtab[f'E2_{f.name}'] = Erms**2 + simtab[f'E_{f.name}']**2
+#            x = simtab[f'E2_{f.name}'] / simtab[f'E_{f.name}']**2
+#            simtab[f'ALPHA_{f.name}'] = (2-x) / (x-1)
+
+        self.filename = os.path.basename(filename)
+
+        super().__init__(simtab, metadata)
+

python/snewpy/models/model_files.yml

@@ -4,6 +4,7 @@
 
 config:
     - &snewpy "https://github.com/SNEWS2/snewpy/raw/v{snewpy_version}/models/{model}/{filename}"
+    - &ccsnmodels "https://github.com/SNEWS2/snewpy-models-ccsn/raw/main/models/{model}/{filename}"
 
 models:
 
@@ -46,3 +47,6 @@ models:
     Zha_2021: 
         repository: *snewpy
 
+    Mori_2023: 
+        repository: *ccsnmodels
+

python/snewpy/test/test_01_registry.py

@@ -8,7 +8,8 @@
 from snewpy.models.ccsn 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
+                          Kuroda_2020, Fornax_2021, Zha_2021, \
+                          Mori_2023
 
 from astropy import units as u
 
@@ -30,7 +31,8 @@ def test_param_combinations(self):
                   Warren_2020,
                   Kuroda_2020,
                   Fornax_2021,
-                  Zha_2021]
+                  Zha_2021,
+                  Mori_2023]
 
         for model in models:
             for pc in model.get_param_combinations():
@@ -325,3 +327,31 @@ def test_Zha_2021(self):
             self.assertEqual(len(f), len(Flavor))
             self.assertEqual(f[Flavor.NU_E].unit, 1./(u.erg * u.s))
 
+    def test_Mori_2023(self):
+        """
+        Instantiate a set of 'Mori 2023' models
+        """
+        axion_mass_coupling = [ (0,0),
+            (100,2), (100,4), (100,10), (100,12), (100,14), (100,16), (100,20),
+            (200,2), (200,4), (200,6), (200,8), (200,10), (200,20)]
+
+        pns_mass = [1.78, 1.77, 1.76, 1.77, 1.77, 1.77, 1.77, 1.74, 1.77, 1.76, 1.75, 1.74, 1.73, 1.62] * u.Msun
+
+        for ((am, ac), mpns) in zip(axion_mass_coupling, pns_mass):
+            model = Mori_2023(axion_mass=am, axion_coupling=ac)
+
+            self.assertEqual(model.metadata['Axion mass'], float(am)*u.MeV)
+            self.assertEqual(model.metadata['Axion coupling'], float(ac)*1e-10/u.GeV)
+            self.assertEqual(model.metadata['Progenitor mass'], 20*u.Msun)
+            self.assertEqual(model.metadata['PNS mass'], mpns)
+
+            # 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))
+