flux computation function
function fluxes(nx,gamma,q,f)
for i = 1:nx+1
p = (gamma-1.0)*(q[i,3]-0.5*q[i,2]*q[i,2]/q[i,1])
f[i,1] = q[i,2]
f[i,2] = q[i,2]*q[i,2]/q[i,1] + p
f[i,3] = q[i,2]*q[i,3]/q[i,1] + p*q[i,2]/q[i,1]
end
end
```
The implementation of WENO-5 reconstruction is similar to Listings 8 and 9. We implemented all Riemann solvers for Euler equations in such a way that we compute smoothness indicators for all equations separately. We do not incur much computational expense since our problem is one-dimensional. However, for higher-dimension usually smoothness indicator is computed only for equation (usually momentum or energy equation) and the same smoothness indicators are used for all equations. The computation of flux using Roe averaging is detailed in Listing 15.

✝ ✆

**Listing 15.** Implementation of Julia function to calculate the interface flux using Riemann solver based on Roe averaging. ✞ ☎

```