/* 1-D Shock tube implemented with Roe's Riemann Solver with entropy fix and 
   Hancock's MUSCL Scheme with Limited Slopes.
   
   All parts of the source code are stored in this single file.
   
   This source code was written by Moritz Nadler and is part of a work done for 
   the Department for Theoretical Astrophysics of the University of Tuebingen.
   It is based on a Fortran 90 program written by J. Naber 
   http://www.cwi.nl/ftp/CWIreports/MAS/MAS-E0502.pdf
   
   Manual:
   In the section "global variable and constants" one can alter the adiabatic
   exponent gamma0, the number of space cell, the length of the tube, the number
   of time steps and the end time of the simulation. The step size for time
   and space are automatically determined from these values.   
   In the main function one can select a test case or add another one, set the 
   boundary conditions and the time steps to be written on hard disk (the last 
   or every n-th time step).
   
   Compilation:
   This source code uses only elements of the C++ standard library and therefor
   should work with every C++ compiler but was only tested with the GNU g++ 
   4.0.3 compiler. If the GNU Compiler is used it is highly recommend to switch 
   on the -O3 flag to prevent the program from being very slow.
   Maybe the first statement in the main function 
   std::ios::sync_with_stdio(false);
   is Compiler specific. It can simply be deleted if necessary.

*/

#include <iostream>
using std::cout; using std::cerr; using std::endl; using std::flush;
#include <vector>
using std::vector;
#include <fstream>
using std::ofstream;
#include <cmath>
using std::sqrt; using std::abs;
#include <algorithm>
using std::max; using std::min; using std::min_element;
#include <string>
using std::string;
#include <sstream>
using std::ostringstream;


/* For a more mathematical sound of the types and to make it possible to 
   alter the precision easily */ 
typedef unsigned int Natural; // 32 bit unsigned integer

typedef double Real; /* alter floting point precision here: float 32 bit,
			double 64 bit, long double 96 bit */
typedef vector<vector<Real> > Matrix; // Matrix of real numbers

/* Self-defined "to the power of 2" function because of the lack of integer 
   powers in C/C++ */
inline Real pow2( const Real argument ){ return argument*argument; }


// globle varialbes and constants

// ratio of specific heats gamma0 = c_p / c_v
const Real gamma0 = 1.4;
const Real gamma5 = (gamma0-1.0)/2.0;
const Real gamma6 = 1.0/gamma5;

// space grid
const Natural xN = 200; // Number of space steps
const Real xMax = 1.0; // Length of shocktube
const Real deltaX = xMax/xN; // Length of space step

//time grid
const Natural tN = 3000; // Number of time steps (1500) 
const Real tMax = 0.5; // End time (0.25)
const Real deltaT = tMax/tN; // Lenght of time steps

// relation of delta t and delta x
const Real dtOverdx = deltaT/deltaX;

// the 3 components of the primary state vector w
Matrix rho( xN+2, tN+1 );	//density(x,t)
Matrix u( xN+2, tN+1 );	//velocity(x,t)
Matrix p( xN+2, tN+1 );	//pressure(x,t)

// Cell face state vectors holds w_R and w_L for every cell
vector<Real> rhoCf( 2*(xN+2)+1 );	
vector<Real> uCf( 2*(xN+2)+1 );	
vector<Real> pCf( 2*(xN+2)+1 );

// the 3 components of the conservative Flux state vector	
vector<Real> f1( xN+2 );	
vector<Real> f2( xN+2 );	
vector<Real> f3( xN+2 );

// function to fill state vectors with values of heavyside shape
void fill( const Real& rhoLeft, const Real& rhoRight, const Real& uLeft,
	   const Real& uRight, const Real& pLeft, const Real& pRight ){
	
	for( Natural i = 1; i not_eq xN+1; ++i ){
		if( i < xN/2 ){
			rho[i][0] = rhoLeft;
			u[i][0] = uLeft;
			p[i][0] = pLeft;
		} else {
			rho[i][0] = rhoRight;
			u[i][0] = uRight;
			p[i][0] = pRight;
		}
	}
}


// Help-functions for averaging
bool absLessThan( const Real& a, const Real& b ){
	return abs(a) < abs(b);
}

Real minmod( const Real& a, const Real& b ){
	if( abs(a) < abs(b) and a*b > 0.0 ){
		return a;
	}
	else if( abs(b) < abs(a) and a*b > 0.0 ){
		return b;
	} else {
		return 0.0;
	}
}

Real maxmod( const Real& a, const Real& b ){
	if( abs(a) > abs(b) and a*b > 0.0 ){
		return a;
	}
	else if( abs(b) > abs(a) and a*b > 0.0 ){
		return b;
	} else {
		return 0.0;
	}
}

// The Limiters
// averaging with algebraic average
void average( vector<Real>& dRho, vector<Real>& dU, vector<Real>& dP,
	      const Natural& t ){
	
	for( Natural i = 1; i not_eq xN+1; ++i ){
		Real dRho1 = rho[i][t] - rho[i-1][t];
		Real dU1 = u[i][t] - u[i-1][t];
		Real dP1 = p[i][t] - p[i-1][t];
		
		Real dRho2 = rho[i+1][t] - rho[i][t];
		Real dU2 = u[i+1][t] - u[i][t];
		Real dP2 = p[i+1][t] - p[i][t];
		
		dRho[i] = (dRho1 + dRho2)/2.0;
		dU[i] = (dU1 + dU2)/2.0;
		dP[i] = (dP1 + dP2)/2.0;
	}
}

// averaging with minmod limiter
void normalMinmod( vector<Real>& dRho, vector<Real>& dU, vector<Real>& dP,
		   const Natural& t ){
	
	for( Natural i = 1; i not_eq xN+1; ++i ){
		Real dRho1 = rho[i][t] - rho[i-1][t];
		Real dU1 = u[i][t] - u[i-1][t];
		Real dP1 = p[i][t] - p[i-1][t];
		
		Real dRho2 = rho[i+1][t] - rho[i][t];
		Real dU2 = u[i+1][t] - u[i][t];
		Real dP2 = p[i+1][t] - p[i][t];
		
		dRho[i] = minmod( dRho1, dRho2 );
		dU[i] = minmod( dU1, dU2 );
		dP[i] = minmod( dP1, dP2 );
	}
}


// averaging with double minmod limiter
void doubleMinmod( vector<Real>& dRho, vector<Real>& dU, vector<Real>& dP,
		   const Natural& t ){
	
	vector<Real> AB(3);
	for( Natural i = 1; i not_eq xN+1; ++i ){
		Real dRho1 = rho[i][t] - rho[i-1][t];
		Real dU1 = u[i][t] - u[i-1][t];
		Real dP1 = p[i][t] - p[i-1][t];
		
		Real dRho2 = rho[i+1][t] - rho[i][t];
		Real dU2 = u[i+1][t] - u[i][t];
		Real dP2 = p[i+1][t] - p[i][t];
				
		
		AB[0] = 0.5*(dRho1 + dRho2);
		AB[1] = 2*dRho1;
		AB[2] = 2*dRho2;
		
		if( dRho1*dRho2 <= 0.0 ){
			dRho[i] = 0.0;
		} else {
			dRho[i] = *min_element( AB.begin(), AB.end(), absLessThan );
		}
		
		
		AB[0] = 0.5*(dU1 + dU2);
		AB[1] = 2.0*dU1;
		AB[2] = 2.0*dU2;
		
		if( dU1*dU2 <= 0.0 ){
			dU[i] = 0.0;
		} else {
			dU[i] = *min_element( AB.begin(), AB.end(), absLessThan );
		}
		
		
		AB[0] = 0.5*(dP1 + dP2);
		AB[1] = 2.0*dP1;
		AB[2] = 2.0*dP2;
		
		if( dP1*dP2 <= 0.0 ){
			dP[i] = 0.0;
		} else {
			dP[i] = *min_element( AB.begin(), AB.end(), absLessThan );
		}	
	}
}

// averaging with superbee limiter
void superbee( vector<Real>& dRho, vector<Real>& dU, vector<Real>& dP,
	       const Natural& t ){

	for( Natural i = 1; i not_eq xN+1; ++i ){
		Real dRho1 = rho[i][t] - rho[i-1][t];
		Real dU1 = u[i][t] - u[i-1][t];
		Real dP1 = p[i][t] - p[i-1][t];
		
		Real dRho2 = rho[i+1][t] - rho[i][t];
		Real dU2 = u[i+1][t] - u[i][t];
		Real dP2 = p[i+1][t] - p[i][t];	
		
		if( dRho1*dRho2 <= 0.0 ){
			dRho[i] = 0.0;
		} else {
			dRho[i] = minmod( maxmod( dRho1, dRho2), minmod( 2.0*dRho1, 2.0*dRho2 ));
		}
		
		if( dU1*dU2 <= 0.0 ){
			dU[i] = 0.0;
		} else {
			dU[i] = minmod( maxmod( dU1, dU2), minmod( 2.0*dU1, 2.0*dU2 ));
		}
		
		if( dP1*dP2 <= 0.0 ){
			dP[i] = 0.0;
		} else {
			dP[i] = minmod( maxmod( dP1, dP2), minmod( 2.0*dP1, 2.0*dP2 ));
		}
	}		
}

// averaging with Koren's Limiter
void koren( vector<Real>& dRho, vector<Real>& dU, vector<Real>& dP,
	    const Natural& t ){
	
	for( Natural i = 1; i not_eq xN+1; ++i ){
		Real dRho1 = rho[i][t] - rho[i-1][t];
		Real dU1 = u[i][t] - u[i-1][t];
		Real dP1 = p[i][t] - p[i-1][t];
		
		Real dRho2 = rho[i+1][t] - rho[i][t];
		Real dU2 = u[i+1][t] - u[i][t];
		Real dP2 = p[i+1][t] - p[i][t];
		
		vector<Real> AB(3);
		
		AB[0] = 2.0*dRho2;
		AB[1] = 2.0*dRho2/3.0 + dRho1/3.0;
		AB[2] = 2*dRho1;
		
		if( dRho1*dRho2 <= 0 ){
			dRho[i] = 0;
		} else {
			dRho[i] = *min_element( AB.begin(), AB.end(), absLessThan );
		}
		
		
		AB[0] = 2.0*dU2;
		AB[1] = 2.0*dU2/3.0 + dU1/3.0;
		AB[2] = 2.0*dU1;
				
		if( dU1*dU2 <= 0.0 ){
			dU[i] = 0.0;
		} else {
			dU[i] = *min_element( AB.begin(), AB.end(), absLessThan );
		}
		

		AB[0] = 2.0*dP2;
		AB[1] = 2.0*dP2/3.0 + dP1/3.0;
		AB[2] = 2.0*dP1;
		
		if( dP1*dP2 <= 0.0 ){
			dP[i] = 0.0;
		} else {
			dP[i] = *min_element( AB.begin(), AB.end(), absLessThan );
		}
	}
}


// averaging with tunable limiter. parameter = 1 => minmod, par. = 2 => superbee
void tunable( vector<Real>& dRho, vector<Real>& dU, vector<Real>& dP,
	      const Natural& t, const Real& parameter ){
	
	const Real null = 0.0;
	for( Natural i = 1; i not_eq xN+1; ++i ){
		Real dRho1 = rho[i][t] - rho[i-1][t];
		Real dU1 = u[i][t] - u[i-1][t];
		Real dP1 = p[i][t] - p[i-1][t];
		
		Real dRho2 = rho[i+1][t] - rho[i][t];
		Real dU2 = u[i+1][t] - u[i][t];
		Real dP2 = p[i+1][t] - p[i][t];
		
		dRho[i] = max( null, min( parameter*dRho2, max( dRho1, min( dRho2, parameter*dRho1 ))))
			 +min( null, max( parameter*dRho2, min( dRho1, max( dRho2, parameter*dRho1 ))));
		
		dU[i] = max( null, min( parameter*dU2, max( dU1, min( dU2, parameter*dU1 ))))
		       +min( null, max( parameter*dU2, min( dU1, max( dU2, parameter*dU1 ))));
		
		dP[i] = max( null, min( parameter*dP2, max( dP1, min( dP2, parameter*dP1 ))))
		       +min( null, max( parameter*dP2, min( dP1, max( dP2, parameter*dP1 ))));
	}
}


// Hancock's Predictor Corrector MUSCL Scheme
void muscl( const Natural& limiter, const Real& parameter, const bool& wall,
	    const Natural& t ){
	// Coefficient matrix A
	Matrix aw(3,3);

	// Define differences state vectors
	vector<Real> dRho(xN+2);
	vector<Real> dU(xN+2);
	vector<Real> dP(xN+2);
	
	// Predictor step calculates dRho, dU, dP
	if( limiter == 0 ){
		/* No averaging at all. All elements in dRho, dU, dP are 0.
		   The default constructor of <vector> did this already. 
		   This is equivalent to fist order upwind method. */
	}
	else if( limiter == 1 ){
		average( dRho, dU, dP, t );
	}
	else if( limiter == 2 ){
		normalMinmod( dRho, dU, dP, t );
	}
	else if( limiter == 3 ){
		doubleMinmod( dRho, dU, dP, t );
	}
	else if( limiter == 4 ){
		superbee( dRho, dU, dP, t );
	}
	else if( limiter == 5 ){
		koren( dRho, dU, dP, t );
	}
	else if( limiter == 6 ){
		tunable( dRho, dU, dP, t, parameter );
	}
	
	// Walls
	if( wall == false ){
		dRho[0] = 0.0;		dRho[xN+1] = 0.0;
		dU[0] = 0.0; 		dU[xN+1] = 0.0;
		dP[0] = 0.0; 		dP[xN+1] = 0.0;
	}
	else if( wall == true ){
		dRho[0] = -dRho[1];	dRho[xN+1] = -dRho[xN];
		dU[0] = dU[1]; 		dU[xN+1] = dU[xN];
		dP[0] = -dP[1]; 	dP[xN+1] = -dP[xN];
	}
	
	// Calculate aw
	for( Natural i = 0; i not_eq xN+2; ++i ){
		aw[0][0] = u[i][t];
		aw[0][1] = rho[i][t];
		aw[0][2] = 0.0;
		aw[1][0] = 0.0;
		aw[1][1] = u[i][t];
		aw[1][2] = 1.0/rho[i][t];
		aw[2][0] = 0.0;
		aw[2][1] = gamma0*p[i][t];
		aw[2][2] = u[i][t];
		
		// Predictor values
		Real rhoP = rho[i][t] - 0.5*dtOverdx*(aw[0][0]*dRho[i] 
				    + aw[0][1]*dU[i] + aw[0][2]*dP[i]);
		
		Real uP = u[i][t] - 0.5*dtOverdx*(aw[1][0]*dRho[i] 
				+ aw[1][1]*dU[i] + aw[1][2]*dP[i]);
		
		Real pP = p[i][t] - 0.5*dtOverdx*(aw[2][0]*dRho[i]
				+ aw[2][1]*dU[i] + aw[2][2]*dP[i]);
		
		rhoCf[2*i+1] = rhoP - 0.5*dRho[i];
		uCf[2*i+1] = uP - 0.5*dU[i];
		pCf[2*i+1] = pP - 0.5*dP[i];
		
		rhoCf[2*i+2] = rhoP + 0.5*dRho[i];
		uCf[2*i+2] = uP + 0.5*dU[i];
		pCf[2*i+2] = pP + 0.5*dP[i];		
	}	
}


void roeSolver( Natural& n ){

	const Real epsilon = 1E-15;
	// loop over space grid
	for( Natural i = 1; i not_eq xN+2; ++i ){

		Real rho1 = rhoCf[2*i];
		Real rho4 = rhoCf[2*i+1];
		Real u1 = uCf[2*i];
		Real u4 = uCf[2*i+1];
		Real p1 = pCf[2*i];
		Real p4 = pCf[2*i+1];
		
		// turn negative valus of p and rho to nearly 0 ( to epsilon)
		if( rho1 <= 0.0 ){ rho1 = epsilon; }
		if( rho4 <= 0.0 ){ rho4 = epsilon; }
		if( p1 <= 0.0 ){ p1 = epsilon; }
		if( p4 <= 0.0 ){ p4 = epsilon; }
		
		// speed of sound
		Real a1 = sqrt(gamma0*p1/rho1);
		Real a4 = sqrt(gamma0*p4/rho4);
		
		// Enthalpy
		Real H1 = p1/(rho1*(gamma0-1.0)) + p1/rho1 + 0.5*pow2(u1);
		Real H4 = p4/(rho4*(gamma0-1.0)) + p4/rho4 + 0.5*pow2(u4);
		
		// Test for cavity
		if( (u1 + gamma6*a1) < (u4 - gamma6*a4) ){
			cerr << "Program quits because of cavity\n"
			     << "n is at: " << n << endl;
			n = tN; /* this will cause the time loop to quit and
				   continue calculations with the next limiter */
		}
		
		// defining roe averaged state quantities
		Real omega = sqrt(rho4/rho1);
		Real rhoh = omega*rho1;
		Real uh = (u1 + omega*u4) / (1 + omega);
		Real Hh = (H1 + omega*H4) / (1 + omega);
		Real ah = sqrt( (gamma0 - 1)*(Hh - 0.5*pow2(uh)) );
		
		// defining variabes for calculation of flux vector
		Real flux1[3];
		Real flux4[3];
		flux1[0] = rho1*u1;
		flux1[1] = rho1*pow2(u1) + p1;
		flux1[2] = rho1*u1*(p1/(rho1*(gamma0-1)) + p1/rho1 + 0.5*pow2(u1));
		flux4[0] = rho4*u4;
		flux4[1] = rho4*pow2(u4) + p4;
		flux4[2] = rho4*u4*(p4/(rho4*(gamma0-1)) + p4/rho4 + 0.5*pow2(u4));
		
		// Eigenvalues
		Real lambda[3];
		lambda[0] = uh - ah;
		lambda[1] = uh;
		lambda[2] = uh + ah;

		// Eigenvectors
		Real v1[3];
		Real v2[3];
		Real v3[3];
		v1[0] = 1.0;
		v1[1] = uh - ah;
		v1[2] = Hh - uh*ah;
		v2[0] = 1.0;
		v2[1] = uh;
		v2[2] = 0.5*pow2(uh);
		v3[0] = 1.0;
		v3[1] = uh + ah;
		v3[2] = Hh + uh*ah;
		
		//absolute eigenvalues
		Real lambda14 = u4 - a4;
		Real lambda11 = u1 - a1;
		Real lambda34 = u4 + a4;
		Real lambda31 = u1 + a1;

		Real dlam3a = 0.0;
		Real dlam3b = 2*(lambda14 - lambda11);
		Real dlam4a = 0.0;
		Real dlam4b = 2*(lambda34 - lambda31);
		Real dlam1a = ah;
		Real dlam1b = max( dlam3a, dlam3b );
		Real dlam2a = ah;
		Real dlam2b = max( dlam4a, dlam4b );
		
		Real dlambda[3];
		dlambda[0] = min( dlam1a, dlam1b );
		dlambda[1] = 0.0;
		dlambda[2] = min( dlam2a, dlam2b );
		
		Real lambdaMin[3];
		Real lambdaPlus[3];
		for( Natural j = 0; j not_eq 3; ++j ){
			if( lambda[j] <= -dlambda[j] ){
				lambdaMin[j] = lambda[j];
				lambdaPlus[j] = 0.0;
			}
			else if( (lambda[j] > -dlambda[j]) and (lambda[j] < dlambda[j]) ){
				lambdaMin[j] = -pow2(lambda[j] - dlambda[j])/(4*dlambda[j]);
				lambdaPlus[j] = pow2(lambda[j] + dlambda[j])/(4*dlambda[j]);
			}
			else if( lambda[j] >= dlambda[j] ){
				lambdaMin[j] = 0.0;
				lambdaPlus[j] = lambda[j];
			}
		}
		// Junmp relations: rho, u, p
		Real drho = rho4 - rho1;
		Real du = u4 - u1;
		Real dp = p4 - p1;

		// jump vector
		Real dV[3];
		dV[0] =  (dp - rhoh*ah*du) / (2.0*pow2(ah));
		dV[1] = -(dp - pow2(ah)*drho) / pow2(ah);
		dV[2] =  (dp + rhoh*ah*du) / (2.0*pow2(ah));
		// calculate flux 
		if( uh >= 0 ){
			f1[i] = flux1[0] + lambdaMin[0]*dV[0]*v1[0];
			f2[i] = flux1[1] + lambdaMin[0]*dV[0]*v1[1];
			f3[i] = flux1[2] + lambdaMin[0]*dV[0]*v1[2];
		} else {
			f1[i] = flux4[0] - lambdaPlus[2]*dV[2]*v3[0];
			f2[i] = flux4[1] - lambdaPlus[2]*dV[2]*v3[1];
			f3[i] = flux4[2] - lambdaPlus[2]*dV[2]*v3[2];
		}
	}
}

int main(){

/* Improves C++ style I/O performance.
   Only useful if C style I/O (e.g. printf) is NOT used */
std::ios::sync_with_stdio(false); 

cout << "Memory for state vectors allocated. tMax is at: " << tMax << endl;

// test case
const Natural testCase = 1; // possible test cases: 1; 2; 3; 4; 5;

// boundery. wall = 1 for walls. wall = 0 for no walls
const bool wall = true;

/* every nth time step will be writen to the output file. If nth = tN only the
   final time step will be writen on harddisk. */
const Natural nth = tN;

const Natural limiterMax = 7; // number of diffenrent limiters
Natural limiter = 0; // coutner to loop over every limiter

// Parameter for tunalbe limiter ( 6 )
Real parameter = 1.0;

// iniciate state vectors (t = 0) dependent on selected test case
if( testCase == 1 ){ // this is Sod's Testcase
	fill( 	1.0, 0.125,	// rhoLeft, rhoRight
		0.0, 0.0,	// uLeft, uRight
		1.0, 0.1 );	// pLeft, pRight
}
else if( testCase == 2 ){
	fill( 	1.0, 1.0,	// rhoLeft, rhoRight
		-2.0, 2.0,	// uLeft, uRight
		0.4, 0.4 );	// pLeft, pRight
}
else if( testCase == 3 ){
	fill( 	1.0, 1.0,	// rhoLeft, rhoRight
		0.0, 0.0,	// uLeft, uRight
		1000.0, 0.01 );	// pLeft, pRight
}
else if( testCase == 4 ){
	fill( 	1.0, 1.0,	// rhoLeft, rhoRight
		0.0, 0.0,	// uLeft, uRight
		0.01, 100.0 );	// pLeft, pRight
}
else if( testCase == 5 ){
	fill( 	5.99924, 5.99242,	// rhoLeft, rhoRight
		19.5975, -6.19633,	// uLeft, uRight
		460.894, 46.0950 );	// pLeft, pRight
} else {
	cout << "no proper test case was selected" << endl;
	return 0;
}

vector<Real> x(xN+1); //holds position of centers of the space grid cells
for( Natural i = 1; i not_eq xN+1; ++i ){
	x[i] = ( i-0.5 )*deltaX; //calculate positions
}

/* the 3 componets of the conservative state vector q necessary for the time update
   because the Riemann solver calculates f(q) NOT f(w). q = ( rho, rho*u, rho*E ) */
vector<Real> q1( xN+2 );
vector<Real> q2( xN+2 );	
vector<Real> q3( xN+2 );

// filenames for the results with different limiters
vector< const char* > filenames(limiterMax);
filenames[0] = "nolimiter";
filenames[1] = "average";
filenames[2] = "normalMinmod";
filenames[3] = "doubleMinmod";
filenames[4] = "superbee";
filenames[5] = "koren";
filenames[6] = "tunable";

cout << "State vectors defined for t = 0. Starting time loop" << endl;

Natural j = 0; // counter for looping over tuning parameter;

/* loop over all averaging methods (limiters) and save the resutls in different
   files. */
while( limiter not_eq limiterMax and j not_eq 5 ){

	if( limiter not_eq 6 ){
		cout << filenames[limiter] << ' ' << flush;
	} else {
		cout << filenames[limiter] << "_p" << parameter << ' ' << flush;
	}
	
	//Mach in Time
	for( Natural n = 1; n not_eq tN+1; ++n ){
		// boundary conditions
		if( wall == false ){
			rho[0][n-1] = rho[1][n-1];
			u[0][n-1] = u[1][n-1];
			p[0][n-1] = p[1][n-1];
			rho[xN+1][n-1] = rho[xN][n-1];
			u[xN+1][n-1] = u[xN][n-1];
			p[xN+1][n-1] = p[xN][n-1];
		}
		else if( wall == true ){	
			rho[0][n-1] = rho[1][n-1];
			rho[xN+1][n-1] = rho[xN][n-1];
			u[0][n-1] = -u[1][n-1];
			u[xN+1][n-1] = -u[xN][n-1];
			p[0][n-1] = p[1][n-1];
			p[xN+1][n-1] = p[xN][n-1];		
		}
		
		// Hancock's Muscl scheme. Updates the cell face vectors
		muscl( limiter, parameter, wall, n-1 );


		// Roe's Riemann Solver. Updates the flux vectors
		roeSolver(n);


		/* update conversaton state vector q with flux vector than
		   convert change in q back to change in primary state vector w */
		for( Natural k = 1; k not_eq xN+1; ++k ){
			// time step of the conservative state vector q
			q1[k] = rho[k][n-1]
				- dtOverdx*( f1[k+1] - f1[k] );
			q2[k] = u[k][n-1] * rho[k][n-1] 
				- dtOverdx*( f2[k+1] - f2[k] );
			q3[k] = p[k][n-1]/(gamma0-1) + 0.5*rho[k][n-1]*pow2(u[k][n-1]) 
				- dtOverdx*( f3[k+1] - f3[k] );
		
			// convert change in q to change in primary state vector w
			rho[k][n] = q1[k];
			u[k][n] = q2[k] / q1[k];
			p[k][n] = (gamma0-1) * (q3[k] - 0.5*pow2(q2[k])/q1[k]);
		}
	}

	/* Plot shocktube at end time with all averageing methods */
	ofstream writefile;
	if( limiter not_eq 6 ){
		writefile.open( filenames[limiter] );
	}
	else if( limiter == 6 ){
		ostringstream pToString; 
		pToString << parameter;
		string tunable = filenames[limiter];
		string parameterValue = pToString.str();
		tunable = tunable + parameterValue;
		writefile.open(tunable.c_str());
		parameter += 0.25;
		--limiter; // do the limiter loop once more
		++j;
	}
	
	// write the results in a file
	// x	density	speed	presure
	//writefile.precision(5); // alter the number of decimals used per value in the output file here
	for(Natural i = 1; i not_eq xN+1; ++i ){
		writefile << x[i]; 
		for(Natural t = 1; t not_eq tN+1; ++t){
			if( t % nth == 0 ){
				writefile << "\t" << rho[i][t] << "\t" << u[i][t] << "\t" << p[i][t];
			}
		}
		writefile << "\n";
	}
	writefile.close();

	++limiter;
}


/* Total Mass inside the box: m(t) should be const. for all t if walls = true
   and decrease if walls = false */
/*
vector<Real> mTotal(tN+1); 
// Calculate mass in box:
for( Natural t = 1; t not_eq tN+1; ++t){
	Real m = 0; 
	for( Natural i = 1; i not_eq xN+1; ++i){
		m += rho[i][t];
	}
	m = m*deltaX;
	mTotal[t] = m;
}
// write m(t) in file
ofstream writeMass("m_t");
for( Natural t = 1; t not_eq tN+1; ++t){
	writeMass << t << "\t" << mTotal[t] << "\n";
}
writeMass.close();

*/
cout << "Output generated.\nProgram finished." << endl;

return 0;
}