% Lax-Friedrichs flux: constant C in every face point - fact that v_0 = 0 used for simplification
C_flux = zeros(N_faces * N_fp,K);
faces = [1:N_fp;(N_fp + 1):(2 * N_fp);(2 * N_fp + 1):(3 * N_fp)]';
for i = 1:N_faces
    face = faces(:,i);
    n_x_face = n_x(face(1),:);
    first_eigenvalue_face = n_x_face * u_0;
    eigenvalues_face = [first_eigenvalue_face; first_eigenvalue_face + a_0; first_eigenvalue_face - a_0];
    C_flux(face,:) = repmat(max(abs(eigenvalues_face)),N_fp,1);
end

% Computing matrices for proper upwind numerical flux
TDplusTinv_flux_cell = cell(N_faces * N_fp,K);
TDminusTinv_flux_cell = cell(N_faces * N_fp,K);

FacePoints_def = [1:N_fp; (N_fp + 1):(2 * N_fp); (2 * N_fp + 1):(3 * N_fp)]';
for k = 1:K
    for f = 1:N_faces
        FacePoints = FacePoints_def(:,f);
        n_1 = n_x(FacePoints(1),k);
        n_2 = n_y(FacePoints(1),k);
        
        first_eigenvalue = n_1 * u_0; % generally n_1 * u_0 + n_2 * v_0
        eigenvalues = [first_eigenvalue, first_eigenvalue, first_eigenvalue + a_0, first_eigenvalue - a_0];
        eigenvaluesPlusID = find(eigenvalues > 0);
        eigenvaluesPlus = zeros(1,4);
        eigenvaluesPlus(eigenvaluesPlusID) = eigenvalues(eigenvaluesPlusID);
        eigenvaluesMinusID = find(eigenvalues < 0);
        eigenvaluesMinus = zeros(1,4);
        eigenvaluesMinus(eigenvaluesMinusID) = eigenvalues(eigenvaluesMinusID);
        Dplus_face = diag(eigenvaluesPlus);
        Dminus_face = diag(eigenvaluesMinus);
        
        n1a0 = n_1 * a_0; n2a0 = n_2 * a_0; rho0a02 = rho_0 * a_0 * a_0;
        T_face = [1, 0, rho_0, rho_0; 0, -n_2, n1a0, -n1a0; 0, n_1, n2a0, -n2a0; 0, 0, rho0a02, rho0a02];
        
        n1_2a0 = n_1 / (2 * a_0); n2_2a0 = n_2 / (2 * a_0); frac = 1 / (2 * rho_0 * a_0 * a_0); 
        Tinv_face = [1, 0, 0, -1 / (a_0 * a_0); 0, -n_2, n_1, 0; 0, n1_2a0, n2_2a0, frac; 0, -n1_2a0, -n2_2a0, frac];
        
        temp_matrix_1 = T_face * Dplus_face * Tinv_face;
        temp_matrix_2 = T_face * Dminus_face * Tinv_face;
        for j = 1:N_fp
            TDplusTinv_flux_cell{FacePoints(j),k} = temp_matrix_1;
            TDminusTinv_flux_cell{FacePoints(j),k} = temp_matrix_2;
        end
    end
end