New ccsn model: Takata_2025

.gitignore

@@ -5,6 +5,7 @@ __pycache__/
 build/
 dist/
 *.egg-info/
+.DS_Store
 
 # default directory for downloaded model files
 SNEWPY_models/

python/snewpy/models/ccsn.py

@@ -369,6 +369,57 @@ def __init__(self, axion_mass:u.Quantity, axion_coupling:u.Quantity):
         self.metadata['PNS mass'] = pns_mass*u.Msun
         return super().__init__(filename, self.metadata)
 
+
+@RegistryModel(
+               _param_validator = lambda p: \
+               (p['progenitor_mass'].to_value('Msun') in (11.2, 20, 25) and p['axion_mass'] == 0 and p['axion_coupling'] == 0) or \
+               (p['progenitor_mass'].to_value('Msun') in (11.2, 20, 25) and p['axion_mass'].to_value('MeV') in (40, 100, 150, 200, 300, 600, 400, 800) and p['axion_coupling'].to_value('1e-10/GeV') in (2, 4, 6, 8, 10)),
+
+               progenitor_mass = Parameter(values=[11.2, 20, 25]<<u.Msun,
+                                          description='Progenitor star mass in units of M☉'),
+               axion_mass = Parameter(values=[0, 40, 100, 150, 200, 300, 400, 600, 800]<<u.MeV,
+                                      description='Axion mass in units of MeV'),
+               axion_coupling = Parameter(values=[0, 2, 4, 6, 8, 10]<<(1e-10/u.GeV),
+                                          description='Axion-photon coupling, in units of 1e-10/GeV',
+                                          precision=2 #round to 1e-12/u.GeV
+                                         )
+              )
+class Takata_2025(loaders.Takata_2025):
+    """
+    Model based on 1D simulations with axionlike particles, Takata, T. et al., PhysRevD.111.103028, 2025.
+    Data from private communication.
+    """
+
+    def __init__(self, progenitor_mass:u.Quantity, axion_mass:u.Quantity, axion_coupling:u.Quantity):
+        """ 
+        This simulation is conducted for a progenitor_mass of 11.2, 20, and 25 M_sun.
+        Each of these progenitors are simulated for a set of ALP parameters: ALP mass, m_a [MeV] = (0, 40, 100, 150, 300, 400, 600, 800), and
+        ALP-photon coupling strength, g_a𝛄 [10^-10/GeV] = (0, 4, 6, 8, 10)
+        """
+        # Make sure the axion coupling is converted to units 1e-10/GeV:
+        # axion_coupling = np.round(axion_coupling.to('1e-10/GeV'))
+
+        if axion_mass == 0:
+            if progenitor_mass == 11.2*u.Msun:
+                filename = '11_000_00.dat'
+            elif progenitor_mass == 20*u.Msun:
+                filename = '20_000_00.dat'
+            else: 
+                filename = '25_000_00.dat'
+
+        else:
+            if progenitor_mass == 11.2*u.Msun:
+                filename = f'11_{axion_mass.to_value("MeV"):3g}_{axion_coupling.to_value("1e-10/GeV"):02g}.dat'
+            elif progenitor_mass == 20*u.Msun:
+                filename = f'20_{axion_mass.to_value("MeV"):3g}_{axion_coupling.to_value("1e-10/GeV"):02g}.dat'
+            else: 
+                filename = f'25_{axion_mass.to_value("MeV"):3g}_{axion_coupling.to_value("1e-10/GeV"):02g}.dat'
+
+        self.metadata = {}
+
+        return super().__init__(filename, self.metadata)
+
+
 @RegistryModel(
     Bfield = ['hydro','L1','L2'],
     direction = ['average','equator','north','south'],
@@ -377,7 +428,7 @@ def __init__(self, axion_mass:u.Quantity, axion_coupling:u.Quantity):
     _param_validator = lambda p: (p['Bfield'] == 'hydro' and p['grav'] == None and p['rotation'] == None ) or
         (p['Bfield'] == 'L1' and p['grav'] == None and p['rotation'] in [0,90]) or
         (p['Bfield'] == 'L2' and p['rotation'] == None and p['grav'] in ['A','B'])
-)
+)        
 class Bugli_2021(loaders.Bugli_2021):
     """Model based on `Buggli (2021) <https://arxiv.org/abs/2105.00665>`_.
     """

python/snewpy/models/ccsn_loaders.py

@@ -846,6 +846,47 @@ def __init__(self, filename, metadata={}):
         super().__init__(simtab, metadata)
 
 
+class Takata_2025(PinchedModel):
+    def __init__(self, filename, metadata={}):
+        """
+        Parameters
+        ----------
+        filename: str
+            Absolute or relative path to file prefix.
+
+        """
+
+        # Open the requested filename using the model downloader.\
+        datafile = self.request_file(filename)
+
+        self.metadata = metadata
+
+        #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(["NU_E", "NU_E_BAR", "NU_X"], 'ebx')):
+            simtab[f'L_{f}'] = simtab[f'{6+j}:Le{fkey}[e/s]'] << u.erg / u.s
+            # Compute the pinch parameter from E_rms and E_avg
+            # E_rms^2/<E>^2 = (2+a)/(1+a)
+            Eavg = simtab[f'{9+j}:Em{fkey}[MeV]']
+            Erms = simtab[f'{12+j}:Er{fkey}[MeV]']
+            x = Erms**2 / Eavg**2
+            alpha = (2-x)/(x-1)
+
+            simtab[f'E_{f}'] = Eavg << u.MeV
+            simtab[f'Erms_{f}'] = Erms << u.MeV
+            simtab[f'ALPHA_{f}'] = alpha
+
+        self.filename = os.path.basename(filename)
+
+        super().__init__(simtab, metadata)
+
+
 class Bugli_2021(PinchedModel):
     """Model based on `Buggli (2021) <https://arxiv.org/abs/2105.00665>`_.
     """

python/snewpy/models/model_files.yml

@@ -54,6 +54,9 @@ models:
         Mori_2023: 
             repository: *ccsn_repository
 
+        Takata_2025:
+            repository: *ccsn_repository
+
         Bugli_2021:
             repository: *ccsn_repository