rho_L = interpolate(rho, pos);
rho_R = interpolate(rho, neg);
rho_F = 0.5*(rho_L + rho_R);

U_L = interpolate(U, pos);
U_R = interpolate(U, neg);
Un_L = U_L & mesh.Sf()/mesh.magSf();
Un_R = U_R & mesh.Sf()/mesh.magSf();

p_L = interpolate(p, pos);
p_R = interpolate(p, neg);

H = (rhoE + p)/rho;
H_L = interpolate(H, pos, T.name());
H_R = interpolate(H, neg, T.name());

volScalarField aAster = sqrt(2.0*(thermo.gamma()-1.0)/(thermo.gamma()+1.0)*H);
surfaceScalarField aAster_L = interpolate(aAster, pos, T.name());
surfaceScalarField aAster_R = interpolate(aAster, neg, T.name());
surfaceScalarField a_L = sqr(aAster_L)/max(aAster_L, Un_L);
surfaceScalarField a_R = sqr(aAster_R)/max(aAster_R, -Un_R);
//interface speed of sound
surfaceScalarField a_F = min(a_L, a_R);

Ma_L = Un_L/a_F;
Ma_R = Un_R/a_F;
sqrMaMean = (sqr(Un_L)+sqr(Un_R))/(2.0*sqr(a_F));
scalar MaInf = runTime.controlDict().lookupOrDefault("MaInf",0.01);
MaZero = sqrt(min(1.0, max(sqrMaMean, sqr(MaInf))));
surfaceScalarField fa = MaZero*(2.0 - MaZero);

const scalar beta = 1.0/8.0;
const scalar sigma = 1.0;
const scalar Kp = 0.25;
const scalar Ku = 0.75;
const scalar zero = 0.0;

surfaceScalarField Ma4PlusL = M4(Ma_L, pos, beta);
surfaceScalarField Ma4MinusR = M4(Ma_R, neg, beta);
surfaceScalarField Map = -Kp/fa*max(1.0 - sigma*sqrMaMean, zero)*(p_R - p_L)/(rho_F * sqr(a_F));

Ma_F = Ma4PlusL + Ma4MinusR + Map;

surfaceScalarField alpha = 3.0/16.0*(-4.0+5.0*fa*fa);
surfaceScalarField p5PlusL = p5(Ma_L, pos, alpha);
surfaceScalarField p5MinusR = p5(Ma_R, neg, alpha);
surfaceScalarField pu = -Ku*p5PlusL*p5MinusR*(rho_L + rho_R)*fa*a_F*(Un_R - Un_L);

p_F = p5PlusL*p_L + p5MinusR*p_R + pu;

surfaceScalarField m = a_F*Ma_F*compareLR(rho_L, rho_R, Ma_F, 0.0);
U_F = compareLR(U_L, U_R, m, 0.0);
H_F = compareLR(H_L, H_R, m, 0.0);

Info << "calculating flux..." << endl;

phi = m*mesh.magSf();
phiUp = m*U_F*mesh.magSf() + p_F*mesh.Sf();
phiEp = m*H_F*mesh.magSf();




















