clc;

rho = rho_IC_values;
u = u_IC_values;
v = v_IC_values;
p = p_IC_values;

reps = 100;

rho_numerical_flux_cycle = zeros(1,length(mapW));
u_numerical_flux_cycle = zeros(1,length(mapW));
v_numerical_flux_cycle = zeros(1,length(mapW));
p_numerical_flux_cycle = zeros(1,length(mapW));

rho_numerical_flux_matrix = zeros(1,length(mapW));
u_numerical_flux_matrix = zeros(1,length(mapW));
v_numerical_flux_matrix = zeros(1,length(mapW));
p_numerical_flux_matrix = zeros(1,length(mapW));

tic;
for i = 1:reps
    for j = 1:length(mapW)
        node_position_faces = mapW(j);
        node_position_total = vmapW(j);
        velocities_modified = [u(node_position_total); v(node_position_total)] - (u(node_position_total) * n_x(node_position_faces) + v(node_position_total) * n_y(node_position_faces)) * [n_x(node_position_faces); n_y(node_position_faces)];
        variables_vector = [rho(node_position_total); velocities_modified; p(node_position_total)];
        numerical_flux_face = (n_x(node_position_faces) * A_1 + n_y(node_position_faces) * A_2) * variables_vector;
        rho_numerical_flux_cycle(j) = numerical_flux_face(1);
        u_numerical_flux_cycle(j) = numerical_flux_face(2);
        v_numerical_flux_cycle(j) = numerical_flux_face(3);
        p_numerical_flux_cycle(j) = numerical_flux_face(4);
    end
end
toc;

tic;
for i = 1:reps
    scalar_product = u(vmapW) .* n_x(mapW) + v(vmapW) .* n_y(mapW);
    u_modified = u(vmapW) - scalar_product .* n_x(mapW);
    v_modified = v(vmapW) - scalar_product .* n_y(mapW);
   
    rho_numerical_flux_matrix = (u_0 * rho(vmapW) + u_modified) .* n_x(mapW) + v_modified .* n_y(mapW);
    u_numerical_flux_matrix = (u_0 * u_modified + p(vmapW)) .* n_x(mapW);
    v_numerical_flux_matrix = (u_0 * v_modified) .* n_x(mapW) + p(vmapW) .* n_y(mapW);
    p_numerical_flux_matrix = (u_modified + u_0 * p(vmapW)) .* n_x(mapW) + v_modified .* n_y(mapW);
end
toc;