%-----------------------------------------------------------
%     Geometricke znazorneni metody nejvetsiho spadu
%-----------------------------------------------------------
clear;
clc;

is_pde=is_toolbox('partial');

if ~is_pde
  disp(' ');
  disp('Lituji, ale neni k dispozici Partial Differential Equation Toolbox !');
  disp(' '); 
else

tisk=sprintf('Geometricke znazorneni metody nejvetsiho spadu');
disp(tisk);
disp(' ');

krok=0.05;
epsilon=4;

mez=-50;
xmin=-2;
xmax=10;
ymin=-4;
ymax=16;
max_pocet=5;

A=[9 2;2 3];
b=[48,26]';
a11=A(1,1);
a12=A(1,2);
a21=A(2,1);
a22=A(2,2);
b1=b(1);
b2=b(2);
res=A\b;
m=res(1);
n=res(2);

const=10.9;

fi=(atan(2*a12/(a11-a22)))/2;
koef_a=1/(0.5*(a11*cos(fi)^2+2*a12*cos(fi)*sin(fi)+a22*sin(fi)^2))^0.5;
koef_b=1/(0.5*(a11*sin(fi)^2-2*a12*cos(fi)*sin(fi)+a22*cos(fi)^2))^0.5;
koef_c=0.5*a11*m^2+a12*m*n+0.5*a22*n^2-b1*m-b2*n;
koef_a=koef_a*const;
koef_b=koef_b*const;

%pdeellip(res(1),res(2),koef_a,koef_b,fi,'e');
gd=[4,res(1),res(2),koef_a,koef_b,fi]';
sf='e';
%g=decsg(gd,sf);
g=decsg(gd);

[P,E,T]=initmesh(g);
[P,E,T]=refinemesh(g,P,E,T);
%P=jigglemesh(P,E,T);

figure(1);
clf;
set(figure(1),'Units','centimeters');
set(figure(1),'Position',[20 6 20 25]);
hold on;
view(20,18);


T=T(1:3,:);
T=T';
P=P';
FP=zeros(max(size(P)),1);

for i=1:max(size(P))
  FP(i)=0.5*P(i,:)*A*(P(i,:))'-b'*(P(i,:))';
end

obj=trisurf(T,P(:,1),P(:,2),FP);
shading interp;
%colormap(pink);
colormap(copper);

rotate3d on;


x0=9;
y0=0;

% Metoda nejvetsiho spadu
[X,Y,pocet]=ns(x0,y0,A,b);
bod3d(X(1),Y(1),hodnotaF(X(1),Y(1),A,b),'ro'); 
bod3d(X(1),Y(1),-200,'ro'); 
disp('  Stiskni klavesu  ');
pause;
pocet=min(pocet,max_pocet);
for i=1:pocet
 bod3d(X(i),Y(i),hodnotaF(X(i),Y(i),A,b),'ro'); 
 bod3d(X(i),Y(i),-200,'ro'); 
 hladina=hodnotaF(X(i),Y(i),A,b);
 elipsa3d(-200,xmin:krok:xmax,hladina,krok,a11,a12,a22,b1,b2,'k-');
 elipsa3d(hladina,xmin:krok:xmax,hladina,krok,a11,a12,a22,b1,b2,'k-');
 j=1;
 clear rX rY rZ;
 for u=0:0.01:1
   rX(j)=u*X(i)+(1-u)*X(i+1);
   rY(j)=u*Y(i)+(1-u)*Y(i+1);
   rZ(j)=hodnotaF(rX(j),rY(j),A,b);
   j=j+1;
 end  
 pr=plot3(rX,rY,rZ,'b-');
 set(pr,'Color',[0 0 1]);
 set(pr,'LineWidth',1.5)
 
 usecka3d(X(i),Y(i),-200,X(i+1),Y(i+1),-200,'b');
 bod3d(X(i+1),Y(i+1),hodnotaF(X(i+1),Y(i+1),A,b),'ro'); 
 bod3d(X(i+1),Y(i+1),-200,'ro'); 
 disp('  Stiskni klavesu  ');
 pause;
end
bod3d(X(pocet),Y(pocet),hodnotaF(X(pocet),Y(pocet),A,b),'ro'); 
bod3d(X(pocet),Y(pocet),-200,'ro'); 
hladina=hodnotaF(X(pocet),Y(pocet),A,b);
elipsa3d(-200,xmin:krok:xmax,hladina,krok,a11,a12,a22,b1,b2,'k-');
elipsa3d(hladina,xmin:krok:xmax,hladina,krok,a11,a12,a22,b1,b2,'k-');

end







