function prednaska04
%% Vstupni data
a = 1; b = 1; c = 2;% [m] rozmery konstrukce
k = 1;        % [Nm-1] tuhost pruziny
F = 1;      % [N] pusobici sila
L1 = sqrt(a^2+b^2); % puvodni delka pruziny 1
L2 = sqrt(b^2+c^2); % puvodni delka pruziny 2
Dx = 1e-7; % Diferencni krok

%% Omezeni a nastaveni bariery
r=1e-4;
rho = 2;

%% Newtonova metoda
x = [0;0]; % Pocatecni odhad
i = 0; % Citac iteraci 
alpha = 1; % Pocatecni delka kroku
clc; % Vymazej obrazovku

%%

while( r < 1e4 )
    g = gfun(x); % Pocatecni hodnota gradientu
    disp ( ['r = ' num2str(r)] );
    
while( norm(g) > 1e-6 )
    i = i + 1; % Zvysim citac iteraci 
    H = hfun(x);
    d = -H\g;
    
    if( g'*d > 0 )
        disp( 'G: ' );
        d = -g; 
    else
        disp( 'N: ' );
    end

    alpha = fminbnd(@fun,0,1);
    x = x + alpha*d;
    g = gfun(x); % Aktualizovana hodnota gradientu
    disp( ['iterace = ' num2str(i) ' gradient = ' num2str(norm(g)) ' energie=' num2str(Pi(x)) ' alpha = ' num2str(alpha)] );
end

cfun(x)
r = 10*r;
pause;

end

x


%% Cilova funkce a analyza citlivosti
    function val=Pi(x)  % Potencialni energie systemu pro posun v
        u=x(1); v=x(2); %  
        D1=sqrt((a+u)^2+(b+v)^2)-L1;
        D2=sqrt((c-u)^2+(b+v)^2)-L2;
        val = .5*k*(D1^2+D2^2) + F*v + r*max(0,cfun(x))^2;
    end

    function val=fun(alpha)
        val=Pi(x+alpha*d);
    end

    function val=gfun(x) % Hodnota gradientu objektivni funkce
        val=zeros(2,1);
        val(1)=(Pi(x+Dx*[1;0])-Pi(x))/Dx;
        val(2)=(Pi(x+Dx*[0;1])-Pi(x))/Dx;
    end

    function val=hfun(x) % Hodnota druhe derivace funkce
        val(1,:)=(gfun(x+Dx*[1;0])-gfun(x))/Dx;
        val(2,:)=(gfun(x+Dx*[0;1])-gfun(x))/Dx;
    end


    function val=cfun(x)
        u=x(1); v=x(2);
        val = u*u + v*v - rho*rho; 
    end
end
