PISN & type Ia models - On hold

doc/source/nb/SNOwGLoBES_models.ipynb

@@ -19,7 +19,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from snewpy.models.ccsn import SNOwGLoBES\n",
+    "from snewpy.models.base import SNOwGLoBES\n",
     "\n",
     "import matplotlib as mpl\n",
     "import matplotlib.pyplot as plt\n",

python/snewpy/models/base.py

@@ -3,13 +3,20 @@
 from abc import ABC, abstractmethod
 
 import numpy as np
+
 from astropy import units as u
 from astropy.table import Table
 from astropy.units import UnitTypeError, get_physical_type
 from astropy.units.quantity import Quantity
 from scipy.special import loggamma
-from snewpy._model_downloader import LocalFileLoader
 
+import tarfile
+
+from scipy import interpolate 
+
+import logging
+
+from snewpy._model_downloader import LocalFileLoader
 from snewpy.flavor import ThreeFlavor
 from snewpy.flavor_transformation import NoTransformation
 from functools import wraps
@@ -23,7 +30,7 @@ def _wrapper(self, *arg, **kwargs):
         init(self, *arg, **kwargs)
         check(self)
     return _wrapper
-    
+
 class SupernovaModel(ABC, LocalFileLoader):
     """Base class defining an interface to a supernova model."""
 
@@ -132,10 +139,10 @@ def get_initial_spectra(self, t, E):
             A container with the information about the initial neutrino spectra
         """
         spectra_dict = self._get_initial_spectra_dict(t, E, flavors=ThreeFlavor)
-        initial_spectra =  flux.Container['1/(MeV*s)'].from_dict(spectra_dict, 
-                                                                time=t,
-                                                                energy=E,
-                                                                flavor_scheme=ThreeFlavor)
+        initial_spectra =  flux.Spectrum.from_dict(spectra_dict, 
+                                                   time=t,
+                                                   energy=E,
+                                                   flavor_scheme=ThreeFlavor)
         return initial_spectra
 
     def get_transformed_spectra(self, t, E, flavor_xform):
@@ -155,8 +162,8 @@ def get_transformed_spectra(self, t, E, flavor_xform):
         flux.Container 
             A container with the information of the transformed neutrino spectra
         """
-        initialspectra = self.get_initial_spectra(t, E)
-        transformed_spectra = flavor_xform.apply_to(initialspectra)
+        initial_spectra = self.get_initial_spectra(t, E)
+        transformed_spectra = flavor_xform.apply_to(initial_spectra)
         return transformed_spectra
 
     def get_flux (self, t, E, distance, flavor_xform=NoTransformation()):
@@ -249,7 +256,7 @@ def _get_initial_spectra_dict(self, t, E, flavors=ThreeFlavor):
         #Reshape the Energy array to shape [1,len(E)]
         E = np.expand_dims(E, axis=0)
 
-        initialspectra = {}
+        initial_spectra = {}
 
         # Estimate L(t), <E_nu(t)> and alpha(t). Express all energies in erg.
         E = E.to_value('erg')
@@ -283,6 +290,101 @@ def _get_initial_spectra_dict(self, t, E, flavors=ThreeFlavor):
 
             #remove unnecessary dimensions, if E or t was scalar:
             result = np.squeeze(result)
-            initialspectra[flavor] = result
+            initial_spectra[flavor] = result
+
+        return initial_spectra
+
+class SNOwGLoBES(SupernovaModel):
+    """Subclass for models whose input data is in SNOwBLoBES format"""
+    def  __init__(self, filename, metadata={}):
+        """
+        Parameters
+        ----------
+        tarfilename: str
+            Absolute or relative path to tar archive
+        """
+        # Open the requested filename using the model downloader.
+        datafile = self.request_file(filename)
 
-        return initialspectra
+        tf = tarfile.open(datafile)
+
+        # Find the "NoOsc" files.
+        datafiles = sorted([f.name for f in tf if '.dat' in f.name])
+        nooscfiles = [df for df in datafiles if 'NoOsc' in df]
+        nooscfiles.sort(key=len)
+
+        # Loop through the NoOsc files and pull out the number fluxes.
+        self.time = []
+        self.energy = None   
+        self.initial_spectra = {}
+        self.interpolation = {}         
+
+        self._flavorkeys = {ThreeFlavor.NU_E: 'NuE',
+                            ThreeFlavor.NU_E_BAR: 'aNuE',
+                            ThreeFlavor.NU_MU: 'NuMu',
+                            ThreeFlavor.NU_MU_BAR: 'aNuMu',
+                            ThreeFlavor.NU_TAU: 'NuTau',
+                            ThreeFlavor.NU_TAU_BAR: 'aNuTau'}        
+
+        for nooscfile in nooscfiles:
+            with tf.extractfile(nooscfile) as f:
+                logging.debug('Reading {}'.format(nooscfile))
+                meta = f.readline()
+                metatext = meta.decode('utf-8')
+                t = float(metatext.split('TBinMid=')[-1].split('sec')[0]) 
+                dt = float(metatext.split('tBinWidth=')[-1].split('s')[0]) 
+                dE = float(metatext.split('eBinWidth=')[-1].split('MeV')[0])
+
+                data = Table.read(f, format='ascii.commented_header', header_start=-1)
+                data.meta['t'] = t
+                data.meta['dt'] = dt
+                data.meta['dE'] = dE
+
+                self.time.append(t)
+
+                if self.energy is None:
+                    self.energy = (data['E(GeV)'].data*1000).tolist()
+
+            for flavor in ThreeFlavor:
+                key = self._flavorkeys[flavor]
+                # convert from flux back to initial spectra: number per /s/erg                
+                spectrum = (data[key].data * (4*np.pi*(u.kpc.to(u.cm,10))**2)/dt/(u.MeV.to(u.erg,dE)) ).tolist() 
+                if flavor in self.initial_spectra:
+                    self.initial_spectra[flavor].append(spectrum)
+                else:
+                    self.initial_spectra[flavor] = [spectrum]      
+
+        for flavor in ThreeFlavor:
+            self.interpolation[flavor] = interpolate.RegularGridInterpolator((self.time, self.energy), self.initial_spectra[flavor], method='cubic')
+
+        self.time *= u.s
+        self.energy *= u.MeV
+
+        self.filename = os.path.basename(filename)            
+
+    def _get_initial_spectra_dict(self, t, E, flavors=ThreeFlavor):
+        """Get neutrino spectra/luminosity curves after oscillation.
+
+        Parameters
+        ----------
+        t : astropy.Quantity
+            Time to evaluate initial spectra.
+        E : astropy.Quantity or ndarray of astropy.Quantity
+            Energies to evaluate the initial spectra.
+        flavors: iterable of snewpy.neutrino.Flavor
+            Return spectra for these flavors only (default: all)
+
+        Returns
+        -------
+        spectra : dict
+            Dictionary of model spectra, keyed by neutrino flavor.
+        """   
+        t = u.Quantity(t, ndmin=1).to(u.s).value
+        E = u.Quantity(E, ndmin=1).to(u.MeV).value
+        tE_grid = np.stack(np.meshgrid(t, E, indexing='ij'), axis=-1)
+
+        initial_spectra = {}
+        for flavor in ThreeFlavor:
+            initial_spectra[flavor] = self.interpolation[flavor](tE_grid) / (u.erg * u.s)
+
+        return initial_spectra

python/snewpy/models/ccsn.py

@@ -457,103 +457,3 @@ def __init__(self, *, Bfield:str, direction:str, rotation:int=None, grav:str=Non
             filename = f'{Bfield}_b12_3d_{rotation}deg_snewpy_{direction}.dat'
         return super().__init__(filename=filename, metadata=self.metadata)
 
-class SNOwGLoBES:
-    """A model that does not inherit from SupernovaModel (yet) and imports a group of SNOwGLoBES files."""
-
-    def __init__(self, tarfilename):
-        """
-        Parameters
-        ----------
-        tarfilename: str
-            Absolute or relative path to tar archive with SNOwGLoBES files.
-        """
-        self.tfname = tarfilename
-        tf = tarfile.open(self.tfname)
-
-        # For now just pull out the "NoOsc" files.
-        datafiles = sorted([f.name for f in tf if '.dat' in f.name])
-        noosc = [df for df in datafiles if 'NoOsc' in df]
-        noosc.sort(key=len)
-
-        # Loop through the noosc files and pull out the number fluxes.
-        self.time = []
-        self.energy = None
-        self.flux = {}
-        self.fmin = 1e99
-        self.fmax = -1e99
-
-        for nooscfile in noosc:
-            with tf.extractfile(nooscfile) as f:
-                logging.debug('Reading {}'.format(nooscfile))
-                meta = f.readline()
-                metatext = meta.decode('utf-8')
-                t = float(metatext.split('TBinMid=')[-1].split('sec')[0])
-                dt = float(metatext.split('tBinWidth=')[-1].split('s')[0])
-                dE = float(metatext.split('eBinWidth=')[-1].split('MeV')[0])
-
-                data = Table.read(f, format='ascii.commented_header', header_start=-1)
-                data.meta['t'] = t
-                data.meta['dt'] = dt
-                data.meta['dE'] = dE
-
-                self.time.append(t)
-                if self.energy is None:
-                    self.energy = (data['E(GeV)'].data*1000).tolist()
-
-            for flavor in ['NuE', 'NuMu', 'NuTau', 'aNuE', 'aNuMu', 'aNuTau']:
-                if flavor in self.flux:
-                    self.flux[flavor].append(data[flavor].data.tolist())
-                else:
-                    self.flux[flavor] = [data[flavor].data.tolist()]
-
-        # We now have a table with rows=times and columns=energies. Transpose
-        # so that rows=energy and cols=time.
-        for k, v in self.flux.items():
-            self.flux[k] = np.transpose(self.flux[k])
-            self.fmin = np.minimum(self.fmin, np.min(self.flux[k]))
-            self.fmax = np.maximum(self.fmax, np.max(self.flux[k]))
-
-    def get_fluence(self, t):
-        """Return the fluence at a given time t.
-
-        Parameters
-        ----------
-        t : float
-            Time in seconds.
-
-        Returns
-        -------
-        fluence : dict
-            A dictionary giving fluence at time t, keyed by flavor.
-        """
-        # Get index of closest element in the array
-        idx = np.abs(np.asarray(self.time) - t).argmin()
-
-        fluence = {}
-        for k, fl in self.flux.items():
-            fluence[k] = fl[:,idx]
-
-        return fluence
-
-class Analytic3Species(PinchedModel):
-    """An analytical model calculating spectra given total luminosity,
-    average energy, and rms or pinch, for each species.
-    """
-
-    param = "There are no input files available for this class. Use `doc/scripts/Analytic.py` in the SNEWPY GitHub repo to create a custom input file."
-
-    def get_param_combinations(cls):
-        print(cls.param)
-        return []
-
-    def __init__(self, filename):
-        """
-        Parameters
-        ----------
-        filename : str
-            Absolute or relative path to file with model data.
-        """
-
-        simtab = Table.read(filename,format='ascii')
-        self.filename = filename
-        super().__init__(simtab, metadata={})

python/snewpy/models/ccsn_loaders.py

@@ -1069,3 +1069,26 @@ def __init__(self, filename, metadata={}):
         tf.close()
 
         super().__init__(simtab, metadata)
+
+class Analytic3Species(PinchedModel):
+    """An analytical model calculating spectra given total luminosity,
+    average energy, and rms or pinch, for each species.
+    """
+
+    param = "There are no repository files available for this class. Use `doc/scripts/Analytic.py` in the SNEWPY GitHub repo to create a custom input file."
+
+    def get_param_combinations(cls):
+        print(cls.param)
+        return []
+
+    def __init__(self, filename):
+        """
+        Parameters
+        ----------
+        filename : str
+            Absolute or relative path to file with model data.
+        """
+
+        simtab = Table.read(filename,format='ascii')
+        self.filename = filename
+        super().__init__(simtab, metadata={})

python/snewpy/models/model_files.yml

@@ -6,7 +6,8 @@ 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.4/models/{model}/{filename}"
     - &presn_repository "https://github.com/SNEWS2/snewpy-models-presn/raw/v0.2/models/{model}/{filename}"
-    - &presn_repository_main "https://github.com/SNEWS2/snewpy-models-presn/raw/main/models/{model}/{filename}"
+    - &presn_repository_main "https://github.com/SNEWS2/snewpy-models-presn/raw/main/models/{model}/{filename}"    
+    - &pisn_typeIa_repository "https://github.com/SNEWS2/snewpy-models-pisn_typeIa/raw/v0.2/models/{model}/{filename}"        
 
 models:
     ccsn:
@@ -82,4 +83,14 @@ models:
             repository: *presn_repository
 
         Myers_2026:
-            repository: *presn_repository_main
+            repository: *presn_repository_main            
+
+    pisn:               
+        Wright_2017:
+            repository: *pisn_typeIa_repository                
+
+    typeIa:
+        Wright_2016:
+            repository: *pisn_typeIa_repository  
+
+

python/snewpy/models/pisn.py

@@ -0,0 +1,36 @@
+# -*- coding: utf-8 -*-
+"""
+The submodule ``snewpy.models.pisn`` contains models of pair-instability supernova neutrino fluxes,
+"""
+import logging
+import os
+
+from astropy import units as u
+
+from snewpy.models import pisn_loaders as loaders
+from .base import SupernovaModel
+
+from snewpy.models.registry_model import RegistryModel, Parameter
+from snewpy.models.registry_model import all_models
+
+@RegistryModel(
+    progenitor_mass = [150, 250] * u.Msun,
+    eos = ['SFHo', 'Helm']    
+)
+class Wright_2017(loaders.Wright_2017):
+    """PISN model described in the paper `Neutrino signal from pair-instability supernovae' by Wright et al.,  
+    Phys. Rev. D96 (2017) 103008  <https://journals.aps.org/prd/abstract/10.1103/PhysRevD.96.103008>`_
+    """
+    def _metadata_from_filename(self, filename:str):
+        metadata = {
+            'Progenitor mass': float(filename.split('_')[1].strip('Msun')) * u.Msun,
+            'EOS': filename.split('_')[2].strip('EOS=')
+        }
+        return metadata    
+
+    def __init__(self, progenitor_mass:u.Quantity, eos:str):
+        filename=f"PISN_{progenitor_mass.to_value('Msun'):.0f}Msun_EOS={eos}_NeutrinoFlux.tar.bz2"
+        super().__init__(filename, self.metadata)
+
+
+

python/snewpy/models/pisn_loaders.py

@@ -0,0 +1,14 @@
+# -*- coding: utf-8 -*-
+"""
+The submodule ``snewpy.models.pisn_loaders`` contains classes to load pair-instability supernova
+models from files stored on disk.
+"""
+
+from snewpy.models.base import SNOwGLoBES
+from snewpy import _model_downloader
+
+class Wright_2017(SNOwGLoBES):
+    """PISN model described in the paper `Neutrino signal from pair-instability supernovae' by Wright et al.,  
+    Phys. Rev. D96 (2017) 103008  <https://journals.aps.org/prd/abstract/10.1103/PhysRevD.96.103008>`
+    """
+    pass