clear;
%% Vstupy:
%K - Konstanta, hodnota okrajové podmínky fí(R1)=K, musí platit K<0:
K=-1; 
%Poloměry mezikruží, musí platit R1 < R2:
R1=3; 
R2=8;
%Vzorkování:
nr=5;
nphi=16;

%% Vytvoření sítě
phi=linspace(0,2*pi,nphi);
r=linspace(R1,R2,nr+1);
[Phigrid,Rgrid] = meshgrid(phi,r);

%% Analytické řešení
K0=K/log(R1/R2);
fi=@(r) ( K0*log(r) - K0*(log(R2)) );
v_r=@(r) K0./r;

fi_analyt=fi(r);
vr_analyt=v_r(r);

%% Výpočet potenciálu a rychlosti v uzlech
%Vektorové pole rychlosti:
v=v_r(Rgrid);    
phi=Phigrid;
%Skalární pole potenciálu:
z=fi(Rgrid);

%% Vykreslení
%Převedení do kartézských souřadnic:
[x,y,z]=pol2cart(Phigrid,Rgrid,z);
[phi_xy,v_xy]=pol2cart(phi,v);

%Vykreslení potenciálu:
figure
surf(x,y,z,'EdgeColor','none')
title('\fontsize{16}\phi')
axis('equal',[-R2-1 R2+1 -R2-1 R2+1])
xlabel('x')
ylabel('y')
colorbar
hold on

%Vykreslení vektorového pole:
figure
quiver(x,y,phi_xy,v_xy,0.5);
grid on
axis('equal',[-R2-1 R2+1 -R2-1 R2+1])
title('Analytické řešení vektorového pole')

%% Výpočet metodou sítí
n=nr;
r=linspace(R1,R2,nr+1);
h=r(2)-r(1);    %krok

%Sestavení matice A a vektoru b soustavy Ax=b:
A(1,1)=2*r(1);
A(1,2)=-(r(1)+h/2);
b(1)=K*(r(1)-h/2);

for i=2:n-2
    A(i,i-1)=(-r(i)+h/2);
    A(i,i)=2*r(i);
    A(i,i+1)=-(r(i)+h/2);    
end
A(n-1,n-2)=(-r(n-1)+h/2);
A(n-1,n-1)=2*r(n-1);
b(n-1)=0;
b=[b(1) zeros(1,n-3) b(n-1)];
A;

%Vyřešení soustavy:
fi_sit=A\b';
%Doplnění o okrajové podmínky:
fi_sit=[K fi_sit' 0];

%Výpočet rychlosti
for i=2:n
v_numer(i)=(fi_sit(i+1)-fi_sit(i-1))/(2*h);
end
vr_analyt;
v_numer=[vr_analyt(1) v_numer(2:end) vr_analyt(end)];

%% Vykreslení numerického výsledku
for i=1:nphi
v_sit(:,i)=v_numer;
end
%Převedení do polárních souřadnic:
[phi_xy2,v_xy2]=pol2cart(phi,v_sit);

%Vykreslení vektorového pole rychlosti:
figure
quiver(x,y,phi_xy2,v_xy2,0.5);
grid on
axis('equal',[-R2-1 R2+1 -R2-1 R2+1])
xlabel('x')
ylabel('y')
title('Numerické řešení vektorového pole')