function [rho_rhs,u_rhs,v_rhs,p_rhs] = ODEs_RHS_sponge(rho,u,v,p,time_RK)

Globals2D;

% local Lax-Friedrichs flux - strong form - fact that v_0 = 0 used for simplification

% computing fluxes on the boundaries within the elements
rho_flux = (u_0 * rho(vmapM) + rho_0 * u(vmapM)) .* n_x(:) + (rho_0 * v(vmapM)) .* n_y(:);
u_flux = (u_0 * u(vmapM) + (1 / rho_0) * p(vmapM)) .* n_x(:);
v_flux = (u_0 * v(vmapM)) .* n_x(:) + ((1 / rho_0) * p(vmapM)) .* n_y(:);
p_flux = (poisson_const * p_0 * u(vmapM) + u_0 * p(vmapM)) .* n_x(:) + (poisson_const * p_0 * v(vmapM)) .* n_y(:);

% computing numerical fluxes
rho_numerical_flux = ((u_0 * (rho(vmapM) + rho(vmapP)) + rho_0 * (u(vmapM) + u(vmapP))) .* n_x(:) + (rho_0 * (v(vmapM) + v(vmapP))) .* n_y(:) + (rho(vmapM) - rho(vmapP)) .* ((1 - alpha_flux) * C_flux(:))) / 2;
u_numerical_flux = ((u_0 * (u(vmapM) + u(vmapP)) + (1 / rho_0) * (p(vmapM) + p(vmapP))) .* n_x(:) + (u(vmapM) - u(vmapP)) .* ((1 - alpha_flux) * C_flux(:))) / 2;
v_numerical_flux = ((u_0 * (v(vmapM) + v(vmapP))) .* n_x(:) + ((1 / rho_0) * (p(vmapM) + p(vmapP))) .* n_y(:) + (v(vmapM) - v(vmapP)) .* ((1 - alpha_flux) * C_flux(:))) / 2;
p_numerical_flux = ((poisson_const * p_0 * (u(vmapM) + u(vmapP)) + u_0 * (p(vmapM) + p(vmapP))) .* n_x(:) + (poisson_const * p_0 * (v(vmapM) + v(vmapP))) .* n_y(:) + (p(vmapM) - p(vmapP)) .* ((1 - alpha_flux) * C_flux(:))) / 2;

% imposing boundary condition (unlimited domain) - postup Toulorge zalozeny na diagonalizaci matice A = n_x * A_1 + n_y * A_2
%for i = 1:length(mapB)
%    node_position_faces = mapB(i);
%    node_position_total = vmapB(i);
%    variables_vector = [rho(node_position_total); u(node_position_total); v(node_position_total); p(node_position_total)];
%    numerical_flux_face = TDplusTinv_flux_cell{node_position_faces} * variables_vector;
%    rho_numerical_flux(node_position_faces) = numerical_flux_face(1);
%    u_numerical_flux(node_position_faces) = numerical_flux_face(2);
%    v_numerical_flux(node_position_faces) = numerical_flux_face(3);
%    p_numerical_flux(node_position_faces) = numerical_flux_face(4);
%end
% imposing boundary condition (unlimited domain) - postup jako ve skriptu Lukacova, tj. pro outflow predepisujeme pouze p, extrapolujeme rho, u, v a pro inflow predepisujeme rho, u, v a extrapolujeme p
% outflow - p(vmapP) = 0
rho_numerical_flux(mapO) = (rho_0 * u(vmapO)) .* n_x(mapO) + (rho_0 * v(vmapO)) .* n_y(mapO);
u_numerical_flux(mapO) = ((1 / rho_0) * (p(vmapO))) .* n_x(mapO) / 2; 
v_numerical_flux(mapO) = ((1 / rho_0) * (p(vmapO))) .* n_y(mapO) / 2;
p_numerical_flux(mapO) = (poisson_const * p_0 * u(vmapO)) .* n_x(mapO) + (poisson_const * p_0 * v(vmapO)) .* n_y(mapO) + (p(vmapO) .* ((1 - alpha_flux) * C_flux(mapO))) / 2;
% inflow - rho(vmapP) = u(vmapP) = v(vmapP) = 0
rho_numerical_flux(mapI) = ((u_0 * rho(vmapI) + rho_0 * u(vmapI)) .* n_x(mapI) + (rho_0 * v(vmapI)) .* n_y(mapI) + rho(vmapI) .* ((1 - alpha_flux) * C_flux(mapI))) / 2;
u_numerical_flux(mapI) = ((u_0 * u(vmapI) + (1 / rho_0) * (p(vmapI) + p(vmapI))) .* n_x(mapI) + u(vmapI) .* ((1 - alpha_flux) * C_flux(mapI))) / 2;
v_numerical_flux(mapI) = ((u_0 * v(vmapI)) .* n_x(mapI) + ((1 / rho_0) * (p(vmapI) + p(vmapI))) .* n_y(mapI) + v(vmapI) .* ((1 - alpha_flux) * C_flux(mapI))) / 2;
p_numerical_flux(mapI) = ((poisson_const * p_0 * u(vmapI) + u_0 * (p(vmapI) + p(vmapI))) .* n_x(mapI) + (poisson_const * p_0 * v(vmapI)) .* n_y(mapI)) / 2;

% resizing both normal and numerical fluxes to appropiate arrays
rho_flux = reshape(rho_flux,N_faces * N_fp,K);
u_flux = reshape(u_flux,N_faces * N_fp,K);
v_flux = reshape(v_flux,N_faces * N_fp,K);
p_flux = reshape(p_flux,N_faces * N_fp,K);

rho_numerical_flux = reshape(rho_numerical_flux,N_faces * N_fp,K);
u_numerical_flux = reshape(u_numerical_flux,N_faces * N_fp,K);
v_numerical_flux = reshape(v_numerical_flux,N_faces * N_fp,K);
p_numerical_flux = reshape(p_numerical_flux,N_faces * N_fp,K);

% building the right-hand side of the ODE system
[rho_x,rho_y] = Grad2D(rho); [u_x,u_y] = Grad2D(u); [v_x,v_y] = Grad2D(v); [p_x,p_y] = Grad2D(p);
rho_rhs = -u_0 * rho_x - rho_0 * (u_x + v_y) + LIFT * (Fscale .* (rho_flux - rho_numerical_flux)) + rho_RHS(x,y,time_RK) - sigma_sponge .* rho;
u_rhs = -u_0 * u_x - (1 / rho_0) * p_x + LIFT * (Fscale .* (u_flux - u_numerical_flux)) - sigma_sponge .* u;
v_rhs = -u_0 * v_x - (1 / rho_0) * p_y + LIFT * (Fscale .* (v_flux - v_numerical_flux)) - sigma_sponge .* v;
p_rhs = -poisson_const * p_0 * (u_x + v_y) - u_0 * p_x + LIFT * (Fscale .* (p_flux - p_numerical_flux)) + p_RHS(x,y,time_RK) - sigma_sponge .* p;

end