function [rho,u,v,p] = ExactSolution(x,y,t)

M = 0;% modifikace oproti puvodnimu zadani (bylo M = 0.5)
alpha_konst = log(2) / 25;
eta_konst = sqrt((x - M * t) ^ 2 + (y - 25) ^ 2);
zeta_konst = sqrt((x - M * t) ^ 2 + (y + 25) ^ 2);

func = @(xi) exp(-(xi .^ 2) / (4 * alpha_konst)) .* sin(t * xi) .* besselj(1, eta_konst * xi) .* xi;
int_1 = quadgk(func, 0, Inf);
func = @(xi) exp(-(xi .^ 2) / (4 * alpha_konst)) .* sin(t * xi) .* besselj(1, zeta_konst * xi) .* xi;
int_2 = quadgk(func, 0, Inf);
u = ((x - M * t) / (2 * alpha_konst * eta_konst)) * int_1 + ((x - M * t) / (2 * alpha_konst * zeta_konst)) * int_2;
v = ((y - 25) / (2 * alpha_konst * eta_konst)) * int_1 + ((y - 25) / (2 * alpha_konst * zeta_konst)) * int_2;
func = @(xi) exp(-(xi .^ 2) / (4 * alpha_konst)) .* cos(t * xi) .* (besselj(0, eta_konst * xi) + besselj(0, zeta_konst * xi)) .* xi;
int_3 = quadgk(func, 0, Inf);
rho = (1 / (2 * alpha_konst)) * int_3;
p = rho;

end
