# read phonon DOS object from file
dos = nlread('si_phonon_dos.nc', PhononDensityOfStates)[0]

# make list of energies
energies = numpy.arange(0.0, 80.0, 0.2)*meV

# calculate the DOS spectrum with two different methods
dos_t = dos.tetrahedronSpectrum(energies)
dos_g = dos.gaussianSpectrum(energies, broadening = 1.0*meV)

#plot the spectra using pylab
import pylab
pylab.figure()
pylab.plot(energies.inUnitsOf(meV), dos_t.inUnitsOf(meV**-1))
pylab.plot(energies.inUnitsOf(meV), dos_g.inUnitsOf(meV**-1))
pylab.xlabel("Energy (meV)")
pylab.ylabel("Phonon DOS (1/meV)")
pylab.show()