function [FH] = Fupwind(UL,UR,gamma)

FL = zeros(size(UL));
FR = zeros(size(UR));
FH = zeros(size(UR));

% Calculate left states and flux
rL  = UL(1);
ruL = UL(2);
rEL = UL(3);

uL  =  ruL/rL;
pL  =  (gamma-1)*(rEL - 0.5*rL*uL^2);

FL(1) = ruL;
FL(2) = ruL*uL + pL;
FL(3) = (rEL  + pL)*uL;

% Calculate right states and flux
rR  = UR(1);
ruR = UR(2);
rER = UR(3);

uR  =  ruR/rR;
pR  =  (gamma-1)*(rER - 0.5*rR*uR^2);

FR(1) = ruR;
FR(2) = ruR*uR + pR;
FR(3) = (rER  + pR)*uR;

% Calculate face states based on a simple averaging 
rH = 0.5*(rL+rR);
uH = 0.5*(uL+uR);
pH = 0.5*(pL+pR);

cH = sqrt(gamma*pH/rH);
HH = 0.5*uH^2 + cH^2/(gamma-1);

% Calculate eigenvalues
lambda(1) = uH - cH;
lambda(2) = uH;
lambda(3) = uH + cH;

% Calculate absolute values of eigenvalues and limit them from zero
epsfix = 0.1;
abslambda = max(abs(lambda),epsfix*cH);

% Calculate right eigenvectors
rvect(:,1) = [1; uH - cH; HH - cH*uH];
rvect(:,2) = [1; uH     ; 0.5*uH^2  ];
rvect(:,3) = [1; uH + cH; HH + cH*uH];

% Calculate wave strengths
dp    = pR - pL;
dr    = rR - rL;
du    = uR - uL;

dw(1) = (dp - rH*cH*du)/(2*cH^2);
dw(2) =  dr - dp/cH^2;
dw(3) = (dp + rH*cH*du)/(2*cH^2);

dwlam = dw.*abslambda;

FH = 0.5*(FR + FL) - 0.5*rvect*dwlam';