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

M = 0.5;
alpha1_konst = log(2) / 9;
alpha2_konst = log(2) / 25;
eta_konst = sqrt((x - M * t) ^ 2 + y ^ 2);

func = @(xi) exp(-(xi .^ 2) / (4 * alpha1_konst)) .* sin(t * xi) .* besselj(1, eta_konst * xi) .* xi;
int_1 = quadgk(func, 0, Inf);
u = ((x - M * t) / (2 * alpha1_konst * eta_konst)) * int_1 + 0.04 * y * exp(-alpha2_konst * ((x - 67 - M * t) ^ 2 + y ^ 2));
v = (y / (2 * alpha1_konst * eta_konst)) * int_1 - 0.04 * (x - 67 - M * t) * exp(-alpha2_konst * ((x - 67 - M * t) ^ 2 + y ^ 2));
func = @(xi) exp(-(xi .^ 2) / (4 * alpha1_konst)) .* cos(t * xi) .* besselj(0, eta_konst * xi) .* xi;
int_2 = quadgk(func, 0, Inf);
p = (1 / (2 * alpha1_konst)) * int_2;
rho = p + 0.1 * exp(-alpha2_konst * ((x - 67 - M * t) ^ 2 + y ^ 2));

end