function priklad_elektrostatika()


epsilon0 = 8.854e-12;
epsilonr = 2;
U = 1; % [Volty]


sirkaOmegar = 2;
vyskaOmegar = 2;
sirkaOmegaPlusMinus = 2;
vyskaOmegaPlusMinus = 6;
x1OmegaPlus = 3;


% diskretizace
h = 0.5;
[Uzly,Segmenty] = diskretizace_elektrostatika_MHP(sirkaOmegar,vyskaOmegar,sirkaOmegaPlusMinus,vyskaOmegaPlusMinus,x1OmegaPlus,h);
vykresli_sit_MHP(Uzly,[Segmenty{1},Segmenty{2},Segmenty{3}],1);
title('Diskretizace deskoveho kondenzatoru metodou hranicnich prvku')

nr = size(Segmenty{1},2);
nPlus = size(Segmenty{2},2);
nMinus = size(Segmenty{3},2);

% sestaveni MHP matic
[poleV,poleL] = MHP_elektrostatika(Uzly,Segmenty);

Ir = eye(size(Segmenty{1},2));

A = [         zeros(nPlus,nr)      ,      poleV{2,1}   ,  poleV{2,2} ,  poleV{2,3} ; ...
              zeros(nMinus,nr)     ,      poleV{3,1}   ,  poleV{3,2} ,  poleV{3,3} ; ...
                 poleV{1,1}        ,     -poleV{1,1}   , -poleV{1,2} , -poleV{1,3} ; ...
      epsilonr*(0.5*Ir+poleL{1,1}) , 0.5*Ir-poleL{1,1} , -poleL{1,2} , -poleL{1,3} ];
      
% sestaveni MHP prave strany
b = [ (U/epsilon0)*ones(nPlus,1) ; -(U/epsilon0)*ones(nMinus,1) ; zeros(nr,1) ; zeros(nr,1) ];


% reseni soustavy Gaussovou eliminaci
q = A\b;
qr = q(1:nr);
qe = q(nr+1:2*nr);
qPlus = q(2*nr+1:2*nr+nPlus);
qMinus = q(2*nr+nPlus+1:2*nr+nPlus+nMinus);


% vykresleni hustoty naboju na hranicich
vykresli_reseni_MHP(Uzly,[Segmenty{1},Segmenty{2},Segmenty{3}],q(nr+1:2*nr+nPlus+nMinus));
title('Nabojova hustota na deskach kondenzatoru')
pause

% vypocet a vykresleni elektrostatickeho pole mimo hranice
x = -0.9:0.5:0.9;
y = -0.9:0.5:0.9;
[X,Y] = meshgrid(x,y);
Body(1,:) = X(:);
Body(2,:) = Y(:);
[ur,gradur] = vypocti_MHP_reseni_v_Omega(Uzly,Segmenty{1},qr,Body);
Ur = X;
Ur(:) = ur;
GradUr1 = X;
GradUr2 = X;
GradUr1(:) = gradur(1,:);
GradUr2(:) = gradur(2,:);
contour(X,Y,Ur)
hold on
quiver(X,Y,GradUr1,GradUr2)
hold on

clear Body
x = -10:0.5:10;
y = -7:0.5:7;
[X,Y] = meshgrid(x,y);
Body(1,:) = X(:);
Body(2,:) = Y(:);
N = size(Body,2);
idx = zeros(1,N);
idx = idx | ( abs(Body(1,:))<sirkaOmegar/2+0.1 & abs(Body(2,:))<vyskaOmegar/2+0.1 );
idx = idx | ( abs(Body(1,:))>x1OmegaPlus-0.1 & abs(Body(1,:))<x1OmegaPlus+sirkaOmegaPlusMinus+0.1 & abs(Body(2,:))<vyskaOmegaPlusMinus/2+0.1 );
Body(:,idx) = 0;
[ur,gradur] = vypocti_MHP_reseni_v_Omega(Uzly,[Segmenty{1},Segmenty{2},Segmenty{3}],[qe;qPlus;qMinus],Body);
ur(idx) = 0;
gradur(:,idx) = 0;
Ur = X;
Ur(:) = ur;
GradUr1 = X;
GradUr2 = X;
GradUr1(:) = gradur(1,:);
GradUr2(:) = gradur(2,:);
contour(X,Y,Ur)
hold on
quiver(X,Y,GradUr1,GradUr2)
hold on
axis image
title('Intenzita elektrostatickeho pole deskoveho kondenzatoru')
