clear;
%% Vstupy:
%Konstanty:
K=0; %hodnota okrajové podmínky fí(mí,3/2pi)=K
L=1;  %hodnota okrajové podmínky fí(mí,2pi)=L

%Krajní hodnoty mí, musí platit mi1<mi2:
mi1=0.5; 
mi2=0.9;

%Konstanta a - parametr eliptických souřadnic:
a=10;

%Vzorkování:
nmi=10;
nny=14;
%% Diskrezizace - síť:
ny=linspace(1.5*pi,2*pi,nny);
mi=linspace(mi1,mi2,nmi);
[Phigrid,Rgrid] = meshgrid(ny,mi);
d_ny=ny(2)-ny(1);
d_mi=mi(2)-mi(1);
nU=nmi*(nny-2); %Celkový počet uzlů

%% Konstanty pro vyplnění matice soustavy:
A=-2*(d_ny*d_ny+d_mi*d_mi);
B=d_ny*d_ny;
C=d_mi*d_mi;

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

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

%Vedlejší diagonály:
for i=2:nmi-1
   for j=1:nny-2
   Bup((nny-2)*(i-1)+j)=B;
   end
end
for i=1:nny-2
Bup(i)=2*B;
end

for i=1:nmi-2
   for j=1:nny-2
   Bdown((nny-2)*(i-1)+j)=B;
   end
end
for i=1:nny-2
Bdown(end+1)=2*B;
end

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

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

%Vektor pravé strany soustavy:
b=zeros(nU,1);
for i=1:nmi
b((i-1)*(nny-2)+1)=-C*K;
b(i*(nny-2))=-C*L;  
end
b;

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

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

%Převod do kartézských souřadnic:
for i=1:nmi
    for j=1:nny
        x(i,j)=a*cosh(mi(i))*cos(ny(j));
        y(i,j)=a*sinh(mi(i))*sin(ny(j));
    end
end

%Výpočet vektorového pole rychlosti:
for i=2:(nmi-1)
    for j=2:(nny-1)
        pom=1/(a*sqrt((sinh(mi(i)))^2+(sin(ny(j)))^2));
        v_r(i,j)=pom*(U_sit(i+1,j)-U_sit(i-1,j))/(2*d_mi);
        v_phi(i,j)=pom*(U_sit(i,j+1)-U_sit(i,j-1))/(2*d_ny);
    end
end

%Převod vektorového pole do kartézských souřadnic:
for i=2:(nmi-1)
   for j=2:(nny-1)
      jmenovatel=a*(((sinh(mi(i)))^2)*((cos(ny(j)))^2)+((cosh(mi(i)))^2)*((sin(ny(j)))^2));
      v_x(i,j)=((sinh(mi(i)))*(cos(ny(j)))*v_r(i,j)-(cosh(mi(i)))*(sin(ny(j)))*v_phi(i,j))/jmenovatel;
      v_y(i,j)=((sinh(mi(i)))*(cos(ny(j)))*v_phi(i,j)+(cosh(mi(i)))*(sin(ny(j)))*v_r(i,j))/jmenovatel;
   end
end
v_x=v_x(2:end,2:end);
v_y=v_y(2:end,2:end);

%% Vykreslení výsledků
%Vykreslení potenciálu:
figure
surf(x,y,U_sit,'EdgeColor','none')
title('\fontsize{16}\phi')
axis equal;
xlabel('x');
ylabel('y');
colorbar
hold on

x=x(2:end-1,2:end-1);
y=y(2:end-1,2:end-1);

%Vykreslení rychlosti:
figure
quiver(x,y,v_x,v_y,1);
grid on
axis equal
xlabel('x');
ylabel('y');
