f: flux at the interface computed using Rusanov scheme
function rusanov(nx,gamma,qL,qR,f,fL,fR)
ps = Array{Float64}(undef,nx+1)
wavespeed(nx,gamma,qL,qR,ps)