Merge master onto branch.

python/asteria/interactions.py

@@ -354,9 +354,9 @@ def _integrated_XSxE(self, params, E, y ):
               1/4. * y**4 * eps_p**2
             - 1/3. * y**3 * (eps_p*eps_m*self.Me/E + 2*eps_p**2 + eps_p**2 * y_ckov)
             + 1/2. * y**2 * (eps_p**2 + eps_m**2 + ( 2*eps_p**2 + eps_p*eps_m*self.Me/E) * y_ckov)
-            +        y    * (eps_m**2 + eps_p**2 ) * y_ckov )
+            -        y    * (eps_m**2 + eps_p**2 ) * y_ckov )
         return XSxE
-
+    
     def cross_section(self, flavor, e_nu, scale_to_H2O=True):
         """Cross section from Marciano and Parsa, J. Phys. G 29:2969, 2003.
 
@@ -422,7 +422,7 @@ def mean_lepton_energy(self, flavor, e_nu, scale_to_H2O=True):
         cut = Enu > (self.e_ckov - self.Me)
 		# See (12) in source - units cm**2 / MeV
         norm = 1.5e-44         
-        ymax = 1./( 1 + self.Me/( 2*Enu[cut] )  )        
+        y_max = 1./( 1 + self.Me/( 2*Enu[cut] )  )        
         y_ckov = (self.e_ckov - self.Me)/Enu[cut]
 
         # Define flavor-dependent parameters
@@ -438,17 +438,11 @@ def mean_lepton_energy(self, flavor, e_nu, scale_to_H2O=True):
             epsilon_m, epsilon_p = epsilons
 
         lep = np.zeros_like( Enu )
-        params = (norm, epsilon_p, epsilon_m, 0) 
+        params = (norm, epsilon_p, epsilon_m, y_ckov) 
 
 
-        # lep[cut] += self._integrated_XSxE(params, Enu[cut], y_max)
-        # lep[cut] -= self._integrated_XSxE(params, Enu[cut], 0)
-        lep[cut] += norm*Enu[cut]**2 * (
-            ymax**4 * epsilon_p**2 /4.
-            - ymax**3 *(epsilon_p*epsilon_m*self.Me/Enu[cut] + 2*epsilon_p**2)/3
-            + ymax**2 * (epsilon_p**2 + epsilon_m**2)/2
-        )
-
+        lep[cut] += self._integrated_XSxE(params, Enu[cut], y_max)
+        lep[cut] -= self._integrated_XSxE(params, Enu[cut], y_ckov)        
 
 		# Note: Units are MeV * cm**2 
         if scale_to_H2O:

python/asteria/source.py

@@ -46,18 +46,27 @@ def __init__(self, name, model, progenitor_mass, progenitor_distance,
         self.v_energy_pdf =  np.vectorize(self.energy_pdf, excluded=['E'], signature='(1,n),(1,n)->(m,n)' )
 
         # Energy CDF, useful for random energy sampling.
-        self.energy_cdf = lambda a, Ea, E: gdtr(1., a + 1., (a + 1.) * (E / Ea))
-
-    def energy_spectrum_profile(self, a, Ea, E):
-        # Returns the energy PDF discretized to every point in the space
-        # ( (a, Ea), E) Where a and Ea are n-element time-profiles of the SN Model 
-        # parameters and E is an m-element list of neutrino energies.
-        # Modifed vesion of the ufunc energy_pdf, but is compatible with vector input
-        # Very memory intensive, for default dt = 0.0001s on [-1,15] and dE = 0.1 [0,100]
-        # The returned n x m matrix of floats is 1.28 GB.
-        # Could be improved.
-        return np.exp((1 + a) * np.log(1 + a) - loggamma(1 + a) - (1 + a) * np.log(Ea)) * \
-               np.exp(np.outer(a, np.log(E)) - np.outer((1 + a)/Ea,  E)).T 
+        self.energy_cdf = lambda a, Ea, E: gdtr(1., a + 1., (a + 1.) * (E / Ea))            
+
+    def parts_by_index(self, x, n): 
+        """Returns a list of size-n numpy arrays containing indices for the 
+        elements of x, and one size-m array ( with m<n ) if there are remaining 
+        elements of x.
+
+
+        Returns
+        -------
+        i_part : list
+            List of index partitions (partitions are numpy array).
+        """
+        nParts = x.size//n    
+        i_part = [ np.arange( i*n, (i+1)*n ) for i in range(nParts) ]
+
+        if len(i_part)*n != x.size:
+            i_part += [ np.arange( len(i_part)*n, x.size ) ]
+
+        return i_part
+
     def get_time(self):
         """Return source time as numpy array.
 .
@@ -214,8 +223,31 @@ def sample_energies(self, t, E, n=1, flavor=Flavor.nu_e_bar):
 
         return energies
 
-    def photonic_energy_per_vol(self, time, E, flavor, photon_spectrum):
-
+    def photonic_energy_per_vol(self, time, E, flavor, photon_spectrum, n=1000):
+        """Compute the energy deposited in a cubic meter of ice by photons
+        from SN neutrino interactions.
+
+        Parameters
+        ----------
+
+        time : float (units s)
+            Time relative to core bounce.            
+        E : `numpy.ndarray`
+            Sorted grid of neutrino energies
+        flavor : :class:`asteria.neutrino.Flavor`
+            Neutrino flavor.
+        photon_spectrum : `numpy.ndarray` (Units vary, m**2)
+            Grid of the product of lepton cross section with lepton mean energy
+            and lepton path length per MeV, sorted according to parameter E
+        n : int
+            Maximum number of time steps to compute at once. A temporary numpy array
+            of size n x time.size is created and can be very memory inefficient.
+
+        Returns
+        -------
+        E_per_V
+            Energy per m**3 of ice deposited  by neutrinos of requested flavor
+        """
         H2O_in_ice = 3.053e28 # 1 / u.m**3
 
         t = time.to(u.s).value
@@ -230,8 +262,11 @@ def photonic_energy_per_vol(self, time, E, flavor, photon_spectrum):
             flux *= 2
 
         print('Beginning {0} simulation....'.format(flavor._name_), end='')
-        E_per_V =  np.zeros( time.size )
-        E_per_V += np.trapz( self.energy_spectrum(t, Enu, flavor) * phot, Enu.value, axis=0)
+        # The following two lines exploit the fact that astropy quantities will
+        # always return a number when numpy size is called on them, even if it is 1.
+        E_per_V =  np.zeros( time.size ) 
+        for i_part in self.parts_by_index(time, n): # Limits memory usage
+            E_per_V[i_part] += np.trapz( self.energy_spectrum(t[i_part], Enu, flavor) * phot, Enu.value, axis=0)
         E_per_V *= flux * H2O_in_ice / ( 4 * np.pi * dist**2)
         print('Completed')
 
@@ -272,8 +307,8 @@ def initialize(config):
             alpha = sn_data_table['ALPHA_{:s}'.format(fl)]
 
             luminosity[flavor] = InterpolatedUnivariateSpline(time, L, ext=1 )
-            mean_energy[flavor] = InterpolatedUnivariateSpline(time, E, ext=1)
-            pinch[flavor] = InterpolatedUnivariateSpline(time, alpha, ext=1) 
+            mean_energy[flavor] = InterpolatedUnivariateSpline(time, E, ext=1 )
+            pinch[flavor] = InterpolatedUnivariateSpline(time, alpha, ext=1 )  
 
     elif config.source.table.format.lower() == 'ascii':
         # ASCII will be supported! Promise, promise.