%-----------------------------------------------------------
%   Geometricke znazorneni principu gradientnich metod
%-----------------------------------------------------------

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 principu gradientnich metod');
disp(tisk);
disp(' ');
clf;
clear;

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);

krok=0.05;
epsilon=4;
mez=-50;
xmin=-2;
xmax=10;
ymin=-4;
ymax=16;

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;
const=11.12;
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 20]);
hold on;
rotate3d on;
view(26,32);


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);


for g=mez:-20:-160
 elipsa3d(-180,xmin:krok:xmax,g,krok,a11,a12,a22,b1,b2,'k-');
 elipsa3d(g,xmin:krok:xmax,g,krok,a11,a12,a22,b1,b2,'k-');
end
elipsa3d(mez,xmin:krok:xmax,mez,krok,a11,a12,a22,b1,b2,'k-');

end












