f: flux at the interface computed using HLLC scheme
function hllc(nx,gamma,uL,uR,f,fL,fR)
gm = gamma-1.0
Ds = Array{Float64}(undef,3)
Ds[1], Ds[2] = 0.0, 1.0
for i = 1:nx+1