clear;
%% Vstupy:
%Konstanty :
K=-1;  %hodnota okrajové podmínky fí(R1,phi)=K, musí platit K<0

%musí platit L1,L2<K, pro symetrický příklad L=L2:
L=1;  %hodnota okrajové podmínky fí(r,0)=L
L2=L;  %hodnota okrajové podmínky fí(r,60°)=0 nebo fí(r,60°)=L

%Poloměry mezikruží, musí platit R1<R2:
R1=3; 
R2=6;

%Vzorkování:
nr=10;
nphi=12;

%%
%Diskrezizace - síť:
phi=linspace(0,pi/3,nphi);
r=linspace(R1,R2,nr);
[Phigrid,Rgrid] = meshgrid(phi,r);
d_phi=phi(2)-phi(1);
d_r=r(2)-r(1);
nU=(nr-1)*(nphi-2); %Celkový počet uzlů ve kterých hledáme řešení

%% Funkce pro vyplnění matice soustavy
A=@(x)(-2*x*x/((d_r)^2)-2/(d_phi)^2); 
B1=@(x)(x*x/((d_r)^2)+x/(2*d_r));
B2=@(x)(x*x/((d_r)^2)-x/(2*d_r));
C=1/((d_phi)^2);

%% Metoda sítí
%Sestavení matice A a vektoru b soustavy Ax=b
U=sparse(nU,nU);

%Hlavní diagonála:
for i=1:nr-1
   for j=1:nphi-2
   Ar((nphi-2)*(i-1)+j)=A(r(i+1));
   end
end

%Vedlejší diagonály:
for i=1:nr-2
   for j=1:nphi-2
   Bup((nphi-2)*(i-1)+j)=B1(r(i+1));
   end
end

for i=1:nr-3
   for j=1:nphi-2
   Bdown((nphi-2)*(i-1)+j)=B2(r(i+2));
   end
end
for i=1:nphi-2
Bdown(end+1)=B1(r(nr))+B2(r(nr));
end

for i=1:nr-1
   for j=1:nphi-2
   Cu((nphi-2)*(i-1)+j)=C;
   end
end
for i=1:nr-2
   Cu((nphi-2)*i)=0;
end
Cu=Cu(1:end-1);

%Vlastní sestavení matice:
U=diag(Ar)+diag(Bup,nphi-2)+diag(Bdown,-(nphi-2))+diag(Cu,1)+diag(Cu,-1);

%Vektor pravé strany soustavy:
b=zeros(nU,1);
for i=1:nphi-2
b(i)=b(i)-B1(r(2))*K;    
end
for i=1:nr-1
b((i-1)*(nphi-2)+1)=b((i-1)*(nphi-2)+1)-C*L;
b(i*(nphi-2))=b(i*(nphi-2))-C*L2;  
end
b;

%Vyřešení soustavy:
U_reseni=U\b;

%Doplnění známých hodnot z okrajových podmínek:
U_sit=zeros(nr,nphi);
for i=1:nphi
U_sit(1,i)=K;
end
for i=2:nr
U_sit(i,1)=L;
U_sit(i,nphi)=L2;
end
for i=2:nr
   for j=2:(nphi-1)
       U_sit(i,j)=U_reseni((i-2)*(nphi-2)+j-1);
   end
end

%Výpočet vektorového pole rychlosti:
for i=2:(nr-1)
    for j=2:(nphi-1)
        v_r(i,j)=(U_sit(i+1,j)-U_sit(i-1,j))/(2*d_r);
        v_phi(i,j)=(U_sit(i,j+1)-U_sit(i,j-1))/(2*d_phi*r(i));
    end
end

%Převedení do kartézských souřadnic:
for i=2:(nr-1)
   for j=2:(nphi-1)
      v_x(i,j)=v_r(i,j)*cos(phi(j))-v_phi(i,j)*sin(phi(j));
      v_y(i,j)=v_r(i,j)*sin(phi(j))+v_phi(i,j)*cos(phi(j));
   end
end
v_x=v_x(2:end,2:end);
v_y=v_y(2:end,2:end);

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

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

[Phigrid2,Rgrid2] = meshgrid(phi(2:nphi-1),r(2:nr-1));
[x,y]=pol2cart(Phigrid2,Rgrid2);

%Vykreslení rychlosti:
figure
quiver(x,y,v_x,v_y,1.5);
grid on
axis equal
axis([0 R2+1 -1 R2]);
xlabel('x');
ylabel('y');

