function p_point = ExactSolution_trapz_vector(x_point,y_point,t)

% Mach number of the mean flow
M = 0.5;

% parameters of the source term
s_epsilon = 0.5;
s_alpha = log(2)/2;
% position of the source
x_s = 0; y_s = 0;
% angular frequency
omega = (2 * pi) / 30;


% integration
int_range = 10;
int_step = 0.005;
x_int = -int_range:int_step:int_range;
n = length(x_int); temp_int = zeros(1,n);
for i = 1:n
    f = exp(-s_alpha * ((x_int(i) - x_s) .^ 2 + (x_int - y_s) .^ 2));
    temp = sqrt((x_point - x_int(i)) .^ 2 + (1 - M ^ 2) * (y_point - x_int) .^ 2);
    H0 = besselh(0, omega * temp / (1 - M ^ 2));
    H1 = besselh(1, omega * temp / (1 - M ^ 2));
    exp_eval = exp(-1i * (omega * (x_point - x_int(i)) * M / (1 - M ^ 2) + omega * t));
    dGdt = H0 .* exp_eval;
    dGdx = (M / (1 - M ^ 2)) * (H0 .* exp_eval) + (-1i * (x_point - x_int(i)) ./ ((1 - M ^ 2) * temp)) .* (H1 .* exp_eval);
    func_eval = f .* (dGdt + M * dGdx);
    
    % singularities of Hankel function
    sing_id = temp == 0;
    func_eval(sing_id) = 0 + 0 * 1i;
    
    temp_int(i) = trapz(x_int,func_eval);
end
p_point = trapz(x_int,temp_int);
p_point = s_epsilon * (omega / (4 * sqrt(1 - M ^ 2))) * p_point;

end
