#include <stdio.h>
#include <stdlib.h>
#include <math.h>

// physical parameters
const double gamma = 5./3.; // gas adiabatic index
const double pressure = 1.0; // gas pressure

// simulation parameters
const int nx = 100; // number of cells
const double dx = 0.01; // cell size
const double dt = 0.001; // time step
const double tmax = 1.0; // maximum time

// arrays for storing cell values
double density[nx];
double velocity[nx];
double energy[nx];

// helper function for computing the pressure
double pressure_func(double density, double energy)
{
  return (gamma-1)*density*energy;
}

// main simulation loop
int main(int argc, char **argv)
{
  // initialize density, velocity, and energy
  for (int i = 0; i < nx; i++)
  {
    density[i] = 1.0;
    velocity[i] = 0.0;
    energy[i] = pressure/(gamma-1);
  }

  // time stepping loop
  for (double t = 0; t < tmax; t += dt)
  {
    // compute fluxes at cell interfaces
    double flux_density[nx+1], flux_momentum[nx+1], flux_energy[nx+1];
    for (int i = 0; i < nx+1; i++)
    {
      // left and right cell values
      double density_l = density[i];
      double velocity_l = velocity[i];
      double energy_l = energy[i];
      double pressure_l = pressure_func(density_l, energy_l);

      double density_r = density[i+1];
      double velocity_r = velocity[i+1];
      double energy_r = energy[i+1];
      double pressure_r = pressure_func(density_r, energy_r);

      // Roe average values
      double rho_avg = sqrt(density_l*density_r);
      double vel_avg = (sqrt(density_l)*velocity_l + sqrt(density_r)*velocity_r)/(sqrt(density_l)+sqrt(density_r));
      double h_avg = (sqrt(density_l)*(energy_l+pressure_l/density_l) + sqrt(density_r)*(energy_r+pressure_r/density_r))/(sqrt(density_l)+sqrt(density_r));
      double a_avg = sqrt((gamma-1)*(h_avg-0.5*vel_avg*vel_avg));

      // left and right wave speeds
      double lambda_l = vel_avg - a_avg;
      double lambda_r = vel_avg + a_avg;

      // compute fluxes
      flux_density[i] = rho_avg*vel_avg;
      flux_momentum[i] = rho_avg*vel_avg*vel_avg