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;

rho_exactvalues = zeros(1,N_p_error2);
u_exactvalues = zeros(1,N_p_error2);
v_exactvalues = zeros(1,N_p_error2);
p_exactvalues = zeros(1,N_p_error2);
rho_localL1errors = zeros(1,K);
u_localL1errors = zeros(1,K);
v_localL1errors = zeros(1,K);
p_localL1errors = zeros(1,K);
rho_localLinftyerrors = 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));
    
    ETA = sqrt((x_error2 - x_a - M_0 * T) .^ 2 + (y_error2 - y_a) .^ 2);
    for j = 1:N_p_error2
        eta = ETA(j);
        func = @(z) exp(-(z .^ 2) / (4 * kappa_1)) .* cos(a_0 * z * T) .* besselj(0, eta * z) .* z;
        comp_int = quadgk(func, 0, Inf, 'MaxIntervalCount', 2000);
        rho_exactvalues(j) = comp_int;
        p_exactvalues(j) = comp_int;
        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;
    end
    rho_exactvalues = (eps_1 / (2 * kappa_1)) * rho_exactvalues + eps_2 * exp(-kappa_2 * ((x_error2 - x_e - M_0 * T) .^ 2 + (y_error2 - y_e) .^ 2));
    u_exactvalues = (eps_1 / (2 * kappa_1)) * ((x_error2 - x_a - M_0 * T) .* u_exactvalues) + eps_3 * ((y_error2 - y_a) .* exp(-kappa_3 * ((x_error2 - x_a - M_0 * T) .^ 2 + (y_error2 - y_a) .^ 2)));
    v_exactvalues = (eps_1 / (2 * kappa_1)) * ((y_error2 - y_a) .* v_exactvalues) - eps_3 * ((x_error2 - x_a - M_0 * T) .* exp(-kappa_3 * ((x_error2 - x_a - M_0 * T) .^ 2 + (y_error2 - y_a) .^ 2)));
    p_exactvalues = (eps_1 / (2 * kappa_1)) * p_exactvalues;
    
    coef = rho(:,i)';
    LagPol_values_error2 = coef * values_node;
    rho_localL1errors(i) = mean(abs(rho_exactvalues - LagPol_values_error2));
    rho_localLinftyerrors(i) = max(abs(rho_exactvalues - LagPol_values_error2));
    coef = u(:,i)';
    LagPol_values_error2 = coef * values_node;
    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;
    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
end
rho_globalL1error = mean(rho_localL1errors)
u_globalL1error = mean(u_localL1errors)
v_globalL1error = mean(v_localL1errors)
p_globalL1error = mean(p_localL1errors)
rho_globalLinftyerror = max(rho_localLinftyerrors)
u_globalLinftyerror = max(u_localLinftyerrors)
v_globalLinftyerror = max(v_localLinftyerrors)
p_globalLinftyerror = max(p_localLinftyerrors)
toc;