clc; clear; close all;

r1 = -1:0.002:1; s1 = r1; N1 = length(r1);
n_points = (N1 * N1 + N1) / 2;
dr = r1;
ds = -1 * ones(1,N1);
for i = 2:N1
    r_old = dr;
    s_old = ds;
    clear dr ds;
    r_temp = r1(1:(N1 - (i - 1)));
    s_temp = s1(i) * ones(1,N1 - (i - 1));
    dr = [r_old, r_temp];
    ds = [s_old, s_temp];
end
[da,db] = rstoab(dr,ds);

N = 10;
N_alpha = 50;
alpha_l = 0;
alpha_u = 2;

repeat = 1;
iterations = 0;

while (repeat == 1)

iterations = iterations + 1;

if (iterations > 1)
    alpha_l = input('Specify lower bound for alpha: ');
    alpha_u = input('Specify upper bound for alpha: ');
end

dalpha = linspace(alpha_l,alpha_u,N_alpha + 1);

LebConst_opt = zeros(1,N_alpha + 1);
for n = 1:(N_alpha + 1)
    alpha_nodes = dalpha(n);

    N_p = ((N + 1) * (N + 2)) / 2;
    lambda1 = zeros(1,N_p); lambda2 = lambda1; lambda3 = lambda1;
    temp = 1;
    for i = 0:N
        for j = 0:(N - i)
            lambda1(temp) = i / N;
            lambda3(temp) = j / N;
            temp = temp + 1;
        end
    end
    lambda2 = 1 - lambda1 - lambda3;

    x_opt = -lambda2 + lambda3;
    y_opt = (-lambda2  - lambda3 + 2 * lambda1) / sqrt(3);

    % jak se tyto vztahy dostaly???
    blend1 = 4 * lambda2 .* lambda3;
    blend2 = 4 * lambda1 .* lambda3;
    blend3 = 4 * lambda1 .* lambda2;

    w_edge1 = Warpfactor(N,lambda3 - lambda2);
    w_edge2 = Warpfactor(N,lambda1 - lambda3);
    w_edge3 = Warpfactor(N,lambda2 - lambda1);

    w_final1 = (1 + (alpha_nodes * lambda1) .^ 2) .* blend1 .* w_edge1;
    w_final2 = (1 + (alpha_nodes * lambda2) .^ 2) .* blend2 .* w_edge2;
    w_final3 = (1 + (alpha_nodes * lambda3) .^ 2) .* blend3 .* w_edge3;

    x_opt = x_opt + 1 * w_final1 + cos(2 * pi / 3) * w_final2 + cos(4 * pi / 3) * w_final3;
    y_opt = y_opt + 0 * w_final1 + sin(2 * pi / 3) * w_final2 + sin(4 * pi / 3) * w_final3;
    
    [r_opt,s_opt] = xytors(x_opt,y_opt);
    values_mode = zeros(N_p,n_points);
    temp = 1;
    for i = 0:N
        for j = 0:(N - i)
            values_mode(temp,:) = ModalBaseValues2D(da,db,i,j);
            temp = temp + 1;
        end
    end
    V_opt = VandermondeMatrix2D(N,r_opt,s_opt);
    values_node_opt = inv(V_opt') * values_mode;
    LebConstPoints_opt = sum(abs(values_node_opt),1);
    LebConst_opt(n) = max(LebConstPoints_opt);
end

plot(dalpha,LebConst_opt);

choice = menu('Do you want to narrow down the search:','Yes','No');
if choice == 2
    repeat = 0;
end

end