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

% an attempt to compute L^infty-error of computed solution
disp('Global error L^p norm for p = 1,infty:');
tic;
N_error2 = 10; % should be chosen based on how fine the grid is
N_p_error2 = (N_error2 * N_error2 + N_error2) / 2;
r_1 = linspace(-1,1,N_error2); % nodes on the face of I on axis x
s_1 = linspace(-1,1,N_error2); % nodes on the face of I on axis y
r_error2 = r_1;
s_error2 = -1 * ones(1,N_error2);
for i = 2:N_error2
    r_old = r_error2;
    s_old = s_error2;
    clear r_error2 s_error2;
    r_temp = r_1(1:(N_error2 - (i - 1)));
    s_temp = s_1(i) * ones(1,N_error2 - (i - 1));
    r_error2 = [r_old, r_temp];
    s_error2 = [s_old, s_temp];
end

[a_error2,b_error2] = rstoab(r_error2,s_error2);
values_mode = zeros(N_p,N_p_error2);
temp = 1;
for i = 0:N
    for j = 0:(N - i)
        values_mode(temp,:) = ModalBaseValues2D(a_error2,b_error2,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_error2);
v_exactvalues = zeros(1,N_p_error2);
p_exactvalues = zeros(1,N_p_error2);
u_localL1errors = zeros(1,K);
v_localL1errors = zeros(1,K);
p_localL1errors = zeros(1,K);
u_localLinftyerrors = zeros(1,K);
v_localLinftyerrors = zeros(1,K);
p_localLinftyerrors = zeros(1,K);
for i = 1:K
    element = E(:,i); v_1_ID = element(1); v_2_ID = element(2); v_3_ID = element(3);
    x_error2 = 0.5 * (-(r_error2 + s_error2) * P(1,v_1_ID) + (r_error2 + 1) * P(1,v_2_ID) + (s_error2 + 1) * P(1,v_3_ID));
    y_error2 = 0.5 * (-(r_error2 + s_error2) * P(2,v_1_ID) + (r_error2 + 1) * P(2,v_2_ID) + (s_error2 + 1) * P(2,v_3_ID));
    
    if ~isempty(find(cols_NaN == i))
        problem_NaN = 1;
        row_NaN = find((x_error2 == 0) & (y_error2 == 0));
    end
    
    ETA = sqrt(x_error2 .^ 2 + y_error2 .^ 2);
    for j = 1:N_p_error2
        eta = ETA(j);
        func = @(z) exp(-(z .^ 2) / (4 * kappa_1)) .* sin(a_0 * z * T) .* besselj(1, eta * z) .* z;
        comp_int = quadgk(func, 0, Inf, 'MaxIntervalCount', 2000);
        u_exactvalues(j) = comp_int / eta;
        v_exactvalues(j) = comp_int / eta;
        func = @(z) exp(-(z .^ 2) / (4 * kappa_1)) .* cos(a_0 * z * T) .* besselj(0, eta * z) .* z;
        p_exactvalues(j) = quadgk(func, 0, Inf, 'MaxIntervalCount', 2000);
    end
    u_exactvalues = (eps_1 / (2 * kappa_1)) * (x_error2 .* u_exactvalues);
    v_exactvalues = (eps_1 / (2 * kappa_1)) * (y_error2 .* 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_localL1errors(i) = mean(abs(u_exactvalues - LagPol_values_error2));
    u_localLinftyerrors(i) = max(abs(u_exactvalues - LagPol_values_error2));
    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_localL1errors(i) = mean(abs(v_exactvalues - LagPol_values_error2));
    v_localLinftyerrors(i) = max(abs(v_exactvalues - LagPol_values_error2));
    coef = p(:,i)';
    LagPol_values_error2 = coef * values_node;
    p_localL1errors(i) = mean(abs(p_exactvalues - LagPol_values_error2));
    p_localLinftyerrors(i) = max(abs(p_exactvalues - LagPol_values_error2));
    
    if (i ~= 0) && (rem(i,1000) == 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_error2);
        v_exactvalues = zeros(1,N_p_error2);
        LagPol_values_error2 = zeros(1,N_p_error2);
    end
    problem_NaN = 0;
end
u_globalL1error = mean(u_localL1errors)
v_globalL1error = mean(v_localL1errors)
p_globalL1error = mean(p_localL1errors)
u_globalLinftyerror = max(u_localLinftyerrors)
v_globalLinftyerror = max(v_localLinftyerrors)
p_globalLinftyerror = max(p_localLinftyerrors)
toc;