disp('Computing error norms of the method:');

tic;
[a,b] = rstoab(r,s);
values_mode = zeros(N_p,N_p);
temp = 1;
for i = 0:N
    for j = 0:(N - i)
        values_mode(temp,:) = ModalBaseValues2D(a,b,i,j);
        temp = temp + 1;
    end
end  
values_node = inv(V') * values_mode;


% [rows_NaN,cols_NaN] = find((x == 0) & (y == 0)); % identification of elements, which have vertex [0,0]
% problem_NaN = 0;

u_exactvalues = zeros(1,N_p);
v_exactvalues = zeros(1,N_p);
p_exactvalues = zeros(1,N_p);
u_errors = zeros(N_p,K);
v_errors = zeros(N_p,K);
p_errors = zeros(N_p,K);
for i = 1:K
    element = E(:,i); v_1_ID = element(1); v_2_ID = element(2); v_3_ID = element(3);
    x_error = 0.5 * (-(r + s) * P(1,v_1_ID) + (r + 1) * P(1,v_2_ID) + (s + 1) * P(1,v_3_ID));
    y_error = 0.5 * (-(r + s) * P(2,v_1_ID) + (r + 1) * P(2,v_2_ID) + (s + 1) * P(2,v_3_ID));
    
%     if ~isempty(find(cols_NaN == i))
%         problem_NaN = 1;
%         row_NaN = find((x_error == 0) & (y_error == 0));
%     end
    
    ETA = sqrt(x_error .^ 2 + y_error .^ 2);
    for j = 1:N_p
        eta = ETA(j);
        func = @(xi) exp(-(xi .^ 2) / (4 * kappa_1)) .* sin(xi * T) .* besselj(1, eta * xi) .* xi;
        comp_int = quadgk(func, 0, Inf, 'MaxIntervalCount', 2000);
        u_exactvalues(j) = comp_int / eta;
        v_exactvalues(j) = comp_int / eta;
        func = @(xi) exp(-(xi .^ 2) / (4 * kappa_1)) .* cos(xi * T) .* besselj(0, eta * xi) .* xi;
        p_exactvalues(j) = quadgk(func, 0, Inf, 'MaxIntervalCount', 2000);
    end
    u_exactvalues = (eps_1 / (2 * kappa_1)) * (x_error .* u_exactvalues);
    v_exactvalues = (eps_1 / (2 * kappa_1)) * (y_error .* v_exactvalues);
    p_exactvalues = (eps_1 / (2 * kappa_1)) * p_exactvalues;
    
    coef = u(:,i)';
%     LagPol_values_error2 = coef * values_node;
%     if problem_NaN == 1 % deleting Inf/NaN entries
%         u_exactvalues(row_NaN) = [];
%         LagPol_values_error2(row_NaN) = [];
%     end
%     u_errors(:,i) = (u_exactvalues - LagPol_values_error2)';
    u_values = coef * values_node;
    u_errors(:,i) = (u_exactvalues - u_values)';
    coef = v(:,i)';
%     LagPol_values_error2 = coef * values_node;
%     if problem_NaN == 1 % deleting Inf/NaN entries
%         v_exactvalues(row_NaN) = [];
%         LagPol_values_error2(row_NaN) = [];
%     end
%     v_errors(:,i) = (v_exactvalues - LagPol_values_error2)';
    v_values = coef * values_node;
    v_errors(:,i) = (v_exactvalues - v_values)';
    coef = p(:,i)';
    p_values = coef * values_node;
%     p_errors(:,i) = (p_exactvalues - LagPol_values_error2)';
    p_errors(:,i) = (p_exactvalues - p_values)';
    
    if (i ~= 0) && (rem(i,500) == 0)
        fprintf(2,'%i of total %i elements analyzed.\n',i,K);
    end
    
%     if problem_NaN == 1 % restoring vector lengths
%         u_exactvalues = zeros(1,N_p);
%         v_exactvalues = zeros(1,N_p);
%         LagPol_values_error2 = zeros(1,N_p);
%     end
%     problem_NaN = 0;
end

% deleting NaN entries
u_errors = u_errors(~isnan(u_errors));
v_errors = v_errors(~isnan(v_errors));
format shortE
u_error_nodes = abs(u_errors);
u_error_max = max(u_error_nodes)
u_error_mean = mean(u_error_nodes)
v_error_nodes = abs(v_errors);
v_error_max = max(v_error_nodes)
v_error_mean = mean(v_error_nodes)
p_error_nodes = abs(p_errors);
p_error_max = max(max(p_error_nodes))
p_error_mean = mean(mean(p_error_nodes))
format short
toc;