<jupyter_start><jupyter_text>Flavor TransformationsThis notebook produces the figures in the SNEWPY paper showing the effect of the flavor transformation prescriptions upon the nakazato-shen-z0.004-t_rev100ms-s20.0 model<jupyter_code>import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np

from astropy import units as u 

from snewpy.neutrino import Flavor, MassHierarchy, MixingParameters
from snewpy.models import Nakazato_2013
from snewpy.flavor_transformation import AdiabaticMSW, NonAdiabaticMSWH, \
                                         TwoFlavorDecoherence, ThreeFlavorDecoherence, \
                                         NeutrinoDecay, AdiabaticMSWes, NonAdiabaticMSWes

mpl.rc('font', size=18)
infile = '../../models/Nakazato_2013/nakazato-shen-z0.004-t_rev100ms-s20.0.fits'
model = Nakazato_2013(infile)
model<jupyter_output><empty_output><jupyter_text>Transformed and untransformed energy integrated flux at EarthCompute and plot the flux at Earth for a SN at 10 kpc with no flavor transformation, and with the chosen flavor transformation.<jupyter_code>def plot_total_flux(model, xform_nmo, xform_imo):
    """Plot initial and oscillated neutrino luminosities.
    
    Parameters
    ----------
    model : SupernovaModel
        An input model from a CCSN simulation.
    flav_xform : FlavorTransformation
        A FlavorTransformation subclass; used to create an instance.
    """
    
    energies = np.linspace(0,60,121) * u.MeV    
    d = (10*u.kpc).to('cm').value # distance to SN
        
    times = model.get_time()
    burst_epoch = times <= 0.1*u.s
    accretion_epoch = (times > 0.1*u.s) & (times <= 0.5*u.s)
    cooling_epoch = (times > 0.5*u.s) & (times <= 10*u.s)
    
    ilum = {}
    olum_nmo = {}
    olum_imo = {}
    
    for flavor in Flavor:
        ilum[flavor] = np.zeros(len(times))
        olum_nmo[flavor] = np.zeros(len(times))
        olum_imo[flavor] = np.zeros(len(times))

    # Compute the transformed and untransformed flux at each time.
    for i, t in enumerate(times):
        ispec = model.get_initial_spectra(t, energies)
        ospec_nmo = model.get_transformed_spectra(t, energies, xform_nmo)
        ospec_imo = model.get_transformed_spectra(t, energies, xform_imo)

        for flavor in Flavor:
            for j, E in enumerate(energies):
                ispec[flavor][j] /= (4.*np.pi*d**2)
                ospec_nmo[flavor][j] /= (4.*np.pi*d**2)
                ospec_imo[flavor][j] /= (4.*np.pi*d**2)
        
        for flavor in Flavor:
            ilum[flavor][i] = np.trapz(ispec[flavor].to('1/(erg*s)'), energies.to('erg')).value 
            olum_nmo[flavor][i] = np.trapz(ospec_nmo[flavor].to('1/(erg*s)'), energies.to('erg')).value
            olum_imo[flavor][i] = np.trapz(ospec_imo[flavor].to('1/(erg*s)'), energies.to('erg')).value            
            
    # make the figures 
    fig, axes = plt.subplots(3,3, figsize=(20,12), tight_layout=True)
    
    smax = [0.,0.,0.]
    titles = ['Untransformed', 'Transformed (NMO)', 'Transformed (IMO)']
    for i, spec in enumerate([ilum, olum_nmo, olum_imo]):
        for j, phase in enumerate([burst_epoch, accretion_epoch, cooling_epoch]):
            ax = axes[i,j]
            timeunits = 'ms' if j==0 else 's'
                
            for flavor in Flavor:
                if i == 0:
                    smax[j] = np.maximum(smax[j], 1.1*np.max(spec[flavor][phase]))
                    
                ax.plot(times[phase].to(timeunits),
                        spec[flavor][phase], label=flavor.to_tex(), lw=3,
                        color='C0' if flavor.is_electron else 'C1',
                        ls='-' if flavor.is_neutrino else ':')
            
            ax.set(xlim=(times[phase][0].to(timeunits).value, times[phase][-1].to(timeunits).value),
                   ylim=(0, smax[j]))
            
            if j==0:
               ax.set(ylabel=r'flux [$10^{16}$ erg$^{-1}$ cm$^{-2}$ s$^{-1}$]')
               ax.legend(loc='upper right', ncol=1, fontsize=18)
            if j==1:
                ax.set(title=titles[i])
            if i < 2:
                ax.set(xticklabels=[])
            else:
                ax.set(xlabel='time [{}]'.format(timeunits))
            
            ax.grid(ls=':')

      
    return fig<jupyter_output><empty_output><jupyter_text>Untransformed and Transformed Spectra at EarthCompute and plot the spectra at Earth for a SN at 10 kpc with no flavor transformation, and with the chosen flavor transformation.<jupyter_code>def plot_spectra(model, xform_nmo, xform_imo, t):
    """Plot initial and oscillated neutrino luminosities.
    
    Parameters
    ----------
    model : SupernovaModel
        An input model from a CCSN simulation.
    flav_xform : FlavorTransformation
        A FlavorTransformation subclass; used to create an instance.
    t : astropy.Quantity
        Time to compute the spectrum.
    """

    energies = np.linspace(0,60,121) * u.MeV   
    d = (10*u.kpc).to('cm').value # distance to SN

    #get the spectra
    ispec = model.get_initial_spectra(t, energies) 
    ospec_nmo = model.get_transformed_spectra(t, energies, xform_nmo)
    ospec_imo = model.get_transformed_spectra(t, energies, xform_imo)

    for flavor in Flavor:
        for j, E in enumerate(energies):
            ispec[flavor][j] /= (4.*np.pi*d**2)
            ospec_nmo[flavor][j] /= (4.*np.pi*d**2)
            ospec_imo[flavor][j] /= (4.*np.pi*d**2)
    
    fig, axes = plt.subplots(2,2, figsize=(18,15), sharex=True, sharey=True, tight_layout=True)

    for i, spec in enumerate([ispec, ospec_nmo]):
        axes[0][0].plot(energies, spec[Flavor.NU_E]/1e16, 
                    label='Untransformed '+Flavor.NU_E.to_tex() if i==0 else 'Transformed '+Flavor.NU_E.to_tex(),
                    color='C0', ls='-' if i==0 else ':', lw=2,  alpha=0.7)
        axes[0][0].plot(energies, spec[Flavor.NU_X]/1e16, 
                    label='Untransformed '+Flavor.NU_X.to_tex() if i==0 else 'Transformed '+Flavor.NU_X.to_tex(),
                    color='C1', ls='-' if i==0 else ':', lw=2,  alpha=0.7)

        axes[0][0].set(title='Neutrinos in the NMO: $t = ${:.1f}'.format(t))
        axes[0][0].grid()
        axes[0][0].legend(loc='upper right', ncol=2, fontsize=16)
    
        axes[0][1].plot(energies, spec[Flavor.NU_E_BAR]/1e16, 
                    label='Untransformed '+Flavor.NU_E_BAR.to_tex() if i==0 else 'Transformed '+Flavor.NU_E_BAR.to_tex(),
                    color='C0', ls='-' if i==0 else ':', lw=2,  alpha=0.7)
        axes[0][1].plot(energies, spec[Flavor.NU_X_BAR]/1e16, 
                    label='Untransformed '+Flavor.NU_X_BAR.to_tex() if i==0 else 'Transformed '+Flavor.NU_X_BAR.to_tex(),
                    color='C1', ls='-' if i==0 else ':', lw=2,  alpha=0.7)

        axes[0][1].set(title='Antineutrinos in the NMO: $t = ${:.1f}'.format(t)) 
        axes[0][1].grid()
        axes[0][1].legend(loc='upper right', ncol=2, fontsize=16)    
    
    for i, spec in enumerate([ispec, ospec_imo]):
        axes[1][0].plot(energies, spec[Flavor.NU_E]/1e16, 
                    label='Untransformed '+Flavor.NU_E.to_tex() if i==0 else 'Transformed '+Flavor.NU_E.to_tex(),
                    color='C0', ls='-' if i==0 else ':', lw=2,  alpha=0.7)
        axes[1][0].plot(energies, spec[Flavor.NU_X]/1e16, 
                    label='Untransformed '+Flavor.NU_X.to_tex() if i==0 else 'Transformed '+Flavor.NU_X.to_tex(),
                    color='C1', ls='-' if i==0 else ':', lw=2,  alpha=0.7)

        axes[1][0].set(xlabel=r'$E$ [{}]'.format(energies.unit), title='Neutrinos in the IMO: $t = ${:.1f}'.format(t))
        axes[1][0].grid()
        axes[1][0].legend(loc='upper right', ncol=2, fontsize=16)
    
        axes[1][1].plot(energies, spec[Flavor.NU_E_BAR]/1e16, 
                    label='Untransformed '+Flavor.NU_E_BAR.to_tex() if i==0 else 'Transformed '+Flavor.NU_E_BAR.to_tex(),
                    color='C0', ls='-' if i==0 else ':', lw=2,  alpha=0.7)
        axes[1][1].plot(energies, spec[Flavor.NU_X_BAR]/1e16, 
                    label='Untransformed '+Flavor.NU_X_BAR.to_tex() if i==0 else 'Transformed '+Flavor.NU_X_BAR.to_tex(),
                    color='C1', ls='-' if i==0 else ':', lw=2,  alpha=0.7)

        axes[1][1].set(xlabel=r'$E$ [{}]'.format(energies.unit), title='Antineutrinos in the IMO: $t = ${:.1f}'.format(t))
        axes[1][1].grid()
        axes[1][1].legend(loc='upper right', ncol=2, fontsize=18)    

    ax = axes[0][0]
    ax.set(ylabel=r'flux [$10^{16}$ erg$^{-1}$ cm$^{-2}$ s$^{-1}$]')
    ax = axes[1][0]
    ax.set(ylabel=r'flux [$10^{16}$ erg$^{-1}$ cm$^{-2}$ s$^{-1}$]')
    
    return fig


fig = plot_total_flux(model, AdiabaticMSW(), AdiabaticMSW(mh=MassHierarchy.INVERTED) )
fig.savefig('flux_adiabaticmsw.pdf')
fig = plot_spectra(model, AdiabaticMSW(), AdiabaticMSW(mh=MassHierarchy.INVERTED), 100*u.ms)
fig.savefig('spectra_adiabaticmsw.pdf')
mp_nmo = MixingParameters()
angles_nmo = mp_nmo.get_mixing_angles() + (10*u.deg,)
xf_nmo = AdiabaticMSWes(angles_nmo)

mp_imo = MixingParameters(mh=MassHierarchy.INVERTED)
angles_imo = mp_imo.get_mixing_angles() + (10*u.deg,)
xf_imo = AdiabaticMSWes(angles_imo)

fig = plot_total_flux(model, xf_nmo, xf_imo)
fig.savefig('flux_adiabaticmswes.pdf')
fig = plot_spectra(model, xf_nmo, xf_imo, 100*u.ms)
fig.savefig('spectra_adiabaticmswes.pdf')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15757640067911716, 'fraction_numerical': 0.028968590831918505, 'mean_word_length': 2.651685393258427, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '568195', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3683, 'n_tokens_neox': 3399, 'n_words': 650, 'path': 'doc/nb/FlavorTransformation.ipynb', 'repo_name': 'soso128/snewpy', 'size': '9424'}