x_l = -20; x_u = 20; y_l = -20; y_u = 20; T = 10; % computational domain

poisson_const = 1.4;
rho_0 = 1;
u_0 = 0;
v_0 = 0;
p_0 = rho_0 / poisson_const;
a_0 = sqrt(poisson_const * p_0 / rho_0);
M_0 = norm([u_0,v_0]) / a_0;

A_1 = [u_0, rho_0, 0, 0; 0, u_0, 0, 1 / rho_0; 0, 0, u_0, 0; 0, poisson_const * p_0, 0, u_0];
A_2 = [v_0, 0, rho_0, 0; 0, v_0, 0, 0; 0, 0, v_0, 1 / rho_0; 0, 0, poisson_const * p_0, v_0];

rho_RHS = @(x,y,t) 0 * x + 0 * y + 0 * t;
u_RHS = @(x,y,t) 0 * x + 0 * y + 0 * t;
v_RHS = @(x,y,t) 0 * x + 0 * y + 0 * t;
p_RHS = @(x,y,t) 0 * x + 0 * y + 0 * t;

eps_1 = 1; b_1 = 2; kappa_1 = log(2) / (b_1 ^ 2);
rho_IC = @(x,y) eps_1 * exp(-kappa_1 * (x .^ 2 + y .^ 2));
u_IC = @(x,y) 0 * x + 0 * y;
v_IC = @(x,y) 0 * x + 0 * y;
p_IC = @(x,y) eps_1 * exp(-kappa_1 * (x .^ 2 + y .^ 2));

mesh_size_parametr = 1.00;
[P,E] = MshToMatlab2D(['Acoustic_pulse_propagation_h=',sprintf('%0.2f',mesh_size_parametr),'.msh']);
%[P,E] = MshToMatlab2D(['Acoustic_pulse_propagation_PML_h=',sprintf('%0.2f',mesh_size_parametr),'.msh']);
%[P,E] = MshToMatlab2D(['Acoustic_pulse_propagation_sponge_h=',sprintf('%0.2f',mesh_size_parametr),'.msh']);
%[Nv, VX, VY, K, EToV] = MeshReaderGambit2D(['Acoustic_pulse_propagation_h=',sprintf('%0.2f',mesh_size_parametr),'.neu']);
%P = [VX; VY];
%E = EToV';

x_axis_limits = [x_l, x_u]; y_axis_limits = [y_l, y_u];