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

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

%format shortE
exact_solution_nodes = ES(x,y); % exact solution in nodes of refined mesh
error_nodes = abs(u - exact_solution_nodes);
error_max = max(max(error_nodes))
error_mean = mean(mean(error_nodes))
%format short
    
% % an attempt to compute error L^2-norm of computed solution
% disp('Global error L^2 norm:');
% tic;
% order_error1 = 15;
% [r_error1,s_error1,w_error1,N_p_error1] = Cubature2D_Hesthaven(order_error1);
% r_error1 = r_error1'; s_error1 = s_error1'; w_error1 = w_error1';
%    
% [a_error1,b_error1] = rstoab(r_error1,s_error1);
% values_mode = zeros(N_p,N_p_error1);
% temp = 1;
% for i = 0:N
%    for j = 0:(N - i)
%        values_mode(temp,:) = ModalBaseValues2D(a_error1,b_error1,i,j);
%        temp = temp + 1;
%    end
% end  
% values_node = inv(V') * values_mode;
%    
% localL2errors = 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_error1 = 0.5 * (-(r_error1 + s_error1) * P(1,v_1_ID) + (r_error1 + 1) * P(1,v_2_ID) + (s_error1 + 1) * P(1,v_3_ID));
%    y_error1 = 0.5 * (-(r_error1 + s_error1) * P(2,v_1_ID) + (r_error1 + 1) * P(2,v_2_ID) + (s_error1 + 1) * P(2,v_3_ID));
%    coef = u(:,i)';
%    LagPol_values_error1 = coef * values_node;
%    localL2errors(i) = sqrt(J(1,i) * sum(w_error1 .* ((ES(x_error1,y_error1) - LagPol_values_error1) .^ 2)));
% end
% globalL2error = sum(localL2errors)
% toc;
% 
% % an attempt to compute error L^infty-norm of computed solution
% disp('Global error L^p norm for p = 1,infty:');
% tic;
% N_error2 = 20; % 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;
% 
% localL1errors = zeros(1,K);
% %localL2errors = zeros(1,K);
% 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));
%     coef = u(:,i)';
%     LagPol_values_error2 = coef * values_node;
%     localL1errors(i) = mean(abs(ES(x_error2,y_error2) - LagPol_values_error2));
%     %localL2errors(i) = mean((ES(x_error2,y_error2) - LagPol_values_error2) .^ 2);
%     localLinftyerrors(i) = max(abs(ES(x_error2,y_error2) - LagPol_values_error2));
% end
% globalL1error = mean(localL1errors)
% %globalL2error = mean(localL2errors)
% globalLinftyerror = max(localLinftyerrors)
% toc;