function priklad_magnetostatika()


mu0 = 4*pi*1e-7;
mur = 5000;
J = 0.25;

sirkaOmegar = 4;
vyskaOmegar = 8;

sirkaOmegaJ = 1;
vyskaOmegaJ = 6;

% diskretizace
h = 1;
[Uzly,Segmenty] = diskretizace_magnetostatika_MHP(sirkaOmegar,vyskaOmegar,h);
% vykresli_sit_MHP(Uzly,[Segmenty{1}],1);
% axis equal
% title('Diskretizace stejnosmerne civky metodou hranicnich prvku')
% pause

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

% sestaveni MHP matic a prave strany
[poleV,poleL,polef,poleg] = MHP_magnetostatika(Uzly,Segmenty,sirkaOmegar,sirkaOmegaJ,vyskaOmegaJ);

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

A = [      poleV{1,1}    ,         -poleV{1,1}          ; ...
      -0.5*Ir+poleL{1,1} , -(1/mur)*[0.5*Ir-poleL{1,1}] ];
  
b = mu0*J * [ -polef{1} ; -poleg{1} ];


% reseni soustavy Gaussovou eliminaci
q = A\b;
q0 = q(1:nr);
qr = q(nr+1:2*nr);


% vykresleni hustoty naboju na hranicich
vykresli_reseni_MHP(Uzly,[Segmenty{1}],q(1:nr));
axis equal
title('Nabojova hustota na rozhrani feromagnetikum-vzduch')
pause

% vypocet a vykresleni elektrostatickeho pole mimo hranice
x = -sirkaOmegar/2+h/5:h/5:sirkaOmegar/2-h/5;
y = -vyskaOmegar/2+h/5:h/5:vyskaOmegar/2-h/5;
[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;
RotUr1 = X;
RotUr2 = X;
RotUr1(:) = gradur(2,:);
RotUr2(:) = -gradur(1,:);
figure
contour(X,Y,Ur)
axis equal
hold on
quiver(X,Y,10*RotUr1,10*RotUr2)
hold on
%
clear Body
x = -2*sirkaOmegar:h/5:2*sirkaOmegar;
y = -vyskaOmegar:h/5:vyskaOmegar;
[X,Y] = meshgrid(x,y);
Body(1,:) = X(:);
Body(2,:) = Y(:);
N = length(X(:));
idx = zeros(1,N);
idx = idx | ( abs(Body(1,:))<(sirkaOmegar/2+h/5) & abs(Body(2,:))<(vyskaOmegar/2+h/5) );
% idx = (abs(Body(1,:))<(sirkaOmegar/2+h/5) & abs(Body(2,:))<(vyskaOmegar/2)+h/5);
Body(:,idx) = 0;
[ur,gradur] = vypocti_MHP_reseni_vne_Omega(Uzly,Segmenty{1},sirkaOmegar,sirkaOmegaJ,vyskaOmegaJ,q0,mu0,J,Body);
ur(idx) = 0;
gradur(:,idx) = 0;
Ur = X;
Ur(:) = ur;
RotUr1 = X;
RotUr2 = X;
RotUr1(:) = gradur(2,:);
RotUr2(:) = -gradur(1,:);
contour(X,Y,Ur)
axis equal
hold on
quiver(X,Y,10*RotUr1,10*RotUr2)
hold on
title('Magnetostaticka indukce civky')
