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

ES = @(x) IC(x - a_advec * T); % function handle exact solution

%format shortE
exact_solution_nodes = ES(x); % exact solution in nodes of refined mesh
error_nodes = abs(u - exact_solution_nodes);
error_max = max(max(error_nodes)) % maximum error in nodes of refined mesh
error_mean = mean(mean(error_nodes))
%format short

%LagPols_cell = LagrangePolynomialBasis(r);
%LagPols_matrix = zeros(N_p);
%for i = 1:N_p
%    LagPols_matrix(i,:) = LagPols_cell{i};
%end
%% computing error L^2-norm exacty using numerical integration (not Gaussian quadrature)
%disp('Global error L^2 norm using numerical integration:');
%tic;
%localL2erros = zeros(1,K);
%for i = 1:K
%    v_1 = E(1,i); v_2 = E(2,i);
%    coef = u(:,i)';
%    LagPol_element = coef * LagPols_matrix;
%    func = @(r) (ES(P(v_1) + ((r + 1) / 2) * (P(v_2) - P(v_1))) - polyval(LagPol_element,r)) .^ 2;
%    localL2errors(i) = sqrt(J(1,i) * quadgk(func,-1,1));
%end
%globalL2error = sum(localL2errors)
%toc;
%
%% an attempt to compute error L^2-norm of computed solution
%disp('Global error L^2 norm using Gaussian quadrature:');
%tic;
%N_p_error1 = 15;
%[GLnodes_error1,GLweights_error1] = GaussLobattoFormula(N_p_error1);
%LagPols_values_error1 = zeros(N_p,N_p_error1);
%for i = 1:N_p
%    LagPols_values_error1(i,:) = polyval(LagPols_cell{i},GLnodes_error1);
%end
%localL2erros = zeros(1,K);
%for i = 1:K
%    v_1 = E(1,i); v_2 = E(2,i);
%    x_error1 = P(v_1) + 0.5 * (GLnodes_error1 + 1) * (P(v_2) - P(v_1));
%    coef = u(:,i)';
%    LagPol_values_error1 = coef * LagPols_values_error1;
%    localL2errors(i) = sqrt(J(1,i) * sum(GLweights_error1 .* ((ES(x_error1) - LagPol_values_error1) .^ 2)));
%end
%globalL2error = sum(localL2errors)
%toc;
%
%% an attempt to compute error L^1-norm of computed solution
%disp('Global error L^1 norm using numerical integration:');
%tic;
%localL1erros = zeros(1,K);
%for i = 1:K
%    v_1 = E(1,i); v_2 = E(2,i);
%    coef = u(:,i)';
%    LagPol_element = coef * LagPols_matrix;
%    func = @(r) abs(ES(P(v_1) + ((r + 1) / 2) * (P(v_2) - P(v_1))) - polyval(LagPol_element,r));
%    localL1errors(i) = J(1,i) * quadgk(func,-1,1);
%end
%globalL1error = sum(localL1errors)
%toc;
%
%% an attempt to compute error L^infty-norm of computed solution
%disp('Global error L^p norms for p = 1,infty:');
%tic;
%%dr_error2 = -1:0.01:1;
%dr_error2 = linspace(-1,1,N_p); % ve stylu Bajka
%LagPols_values_error2 = zeros(N_p,length(dr_error2));
%for i = 1:N_p
%    LagPols_values_error2(i,:) = polyval(LagPols_cell{i},dr_error2);
%end
%localL1errors = zeros(1,K);
%%localL2errors = zeros(1,K);
%localLinftyerrors = zeros(1,K);
%for i = 1:K
%    v_1 = E(1,i); v_2 = E(2,i);
%    x_error2 = P(v_1) + 0.5 * (dr_error2 + 1) * (P(v_2) - P(v_1));
%    coef = u(:,i)';
%    LagPol_values_error2 = coef * LagPols_values_error2;
%    localL1errors(i) = mean(abs(ES(x_error2) - LagPol_values_error2));
%    %localL2errors(i) = sqrt(mean((ES(x_error2) - LagPol_values_error2) .^ 2));
%    localLinftyerrors(i) = max(abs(ES(x_error2) - LagPol_values_error2));
%end
%globalL1error = mean(localL1errors)
%%globalL2error = sqrt(mean(localL2errors .^ 2))
%globalLinftyerror = max(localLinftyerrors)
%toc;