Script for plotting the cross-check of Yoshida_2016 model

plot_Yoshida_fig5.py

import numpy as np
import matplotlib.pyplot as plt
from astropy import units as u

#import the snewpy modules
from snewpy.models import presn
from snewpy.flavor_transformation import NoTransformation, AdiabaticMSW, MassHierarchy
from snewpy.neutrino import Flavor

# set the plot parameters
plt.rc('grid', ls=':')
plt.rc('axes', grid=True)
plt.rc('legend', fontsize=12, loc='upper right')

# define drawing styles for each flavor
styles = {f: dict(color='C0' if f.is_electron else 'C1',
                  ls='-' if f.is_neutrino else ':',
                  label=f.to_tex()) for f in Flavor}

## Initialize the model and calculate the flux

import numpy as np
import matplotlib.pyplot as plt
from astropy import units as u

#import the snewpy modules
from snewpy.models import presn
from snewpy.flavor_transformation import NoTransformation, AdiabaticMSW, MassHierarchy
from snewpy.neutrino import Flavor
# Energy array and time to compute spectra.
# Note that any convenient units can be used and the calculation will remain internally consistent.
model = presn.Yoshida_2016(progenitor_mass=15*u.Msun)
E = np.linspace(0, 25, 201) * u.MeV
t = np.geomspace(-30*u.day, -1e-3*u.day, 1001)
distance = 1*u.kpc

flux = model.get_flux(t, E,
                      distance=distance,
                      flavor_xform=NoTransformation())

neutrino_rate = flux*(4*np.pi*distance**2).to('m**2')
rate = neutrino_rate.integrate('energy') 

#get the relevant arrays for plotting
x = rate.time << u.day
y = rate.array.squeeze()
#we can also convert to a more convenient unit
y = y.to('1/s')
#plot each flavor
for flv in rate.flavor:
    plt.plot(x,y[flv], **styles[flv])
#add the legend
plt.legend(ncol=2, loc='lower center')
#adjust the scales
plt.autoscale(tight=True,axis='x')
plt.yscale('log')
#plt.ylim(1)
#define the labels
xlabel = 'Time before collapse'
ylabel = 'Neutrino rate'
plt.xlabel(f'{xlabel} [{x.unit._repr_latex_()}]')
plt.ylabel(f'{ylabel} [{y.unit._repr_latex_()}]')
#add the plot title
#plt.title('Integral neutrino rate')
#add the model parameters
plt.annotate(str(model) + f'\ndistance : {distance}',
             xy=(0.02,0.98),
             xycoords='axes fraction',
             va='top', ha='left'
            )
plt.xscale('symlog', linthresh=1e-5, subs=-np.power(10.,np.arange(-4,3,1)))
plt.ylim(1e49,1e53)
plt.show()