import numpy as np

import astropy.units as u
from astropy.modeling.models import Gaussian1D
from astropy import visualization

from matplotlib import pyplot as plt
from roentgen.nuclides import Nuclide

ba133 = Nuclide('Ba', 133)
lines = ba133.get_lines(1 * u.keV, 100 * u.keV, min_intensity=1)
energy_resolution = 0.5 * u.keV

g0 = Gaussian1D(amplitude=lines[0]['intensity'], mean=lines[0]['energy'], stddev=energy_resolution)
for this_line in lines[1:]:
    g0 += Gaussian1D(amplitude=this_line['intensity'], mean=this_line['energy'], stddev=energy_resolution)

energy_ax = np.arange(1, 85, 0.5) * u.keV

with visualization.quantity_support():
    fig, ax = plt.subplots()
    ax.plot(energy_ax, g0(energy_ax))
    ax.set_xlabel(f'Energy')
    ax.vlines(lines['energy'], ymin=[0]*len(lines), ymax=lines['intensity'])
    plt.show()