g=9.81;

% definice topologie a geometrie vetvene site 
% --------------------------------------------
% xi yi zi (3 sloupce)
xyz=load('xyz.txt'); %souradnice uzlu
% nuc-celkovy pocet uzlu, vcetne centralnich uzlu vetveni
nuc = size(xyz,1);
% conektivita elementu vetveni. 8 sloupcu
% i1 i2 i3       icentr       d1 d2 d3 alfa
% koncove uzly   centr.uzel   prumery  uhel vetveni
cv=load('cv.txt');
nev = size(cv,1);
% conektivita elementu trubka. 3 sloupce
% i1 i2           d
% koncove uzly  prumer
cp=load('cp.txt');
nep = size(cp,1);
% okrajove podminky pro tlaky nebo prutoky. 2 sloupce
% i  pQ
% index uzlu (v absolutni hodnote). i>0 tlak, i<0 prutok
bc=load('bc.txt');
nbc = size(bc,1);
%
% kontroly rozmeru
nu=nuc-nev
neq=3*nev+2*nep+nbc
if neq==2*nu
    disp 'topologie je zadana spravne (pocty rovnic souhlasi)'
else
    disp 'mas nekde chybu v zadani topologie - nesouhlasi pocty rovnic'
    return
end
% prenos do vektoru
x=xyz(1:end,1);
y=xyz(1:end,2);
z=xyz(1:end,3);
iuv=cv(:,1:3);
iuc=cv(:,4);
dv=cv(:,5:7);
alpha=cv(:,8);
iup=cp(:,1:2);
dp=cp(:,3);
ib=bc(:,1);
bpq=bc(:,2);
%
% Vysledky budou ulozeny do nasledujicich matic
% pnew(nu), qnew(nu) tlaky a prutoky v uzlech
% rev(nev,1:3)  Reynoldsova cisla ve vetvich rozbocek
% rep(nep)      Reynoldsova cisla v primych trubkach
% ezv(nev,1)    celkovy ztratovy vykon [W] vetveni 
% ezv(nev,2)    vykon odpovidajici jen treci ztrate
% ezp(nep,1)    celkovy ztratovy vykon trubky
% ezp(nep,2)    vykon trubky jen dle Hagen Poisseuile
%
% =======================================================================
% soustava rovnic s matici a(neq,neq), vektor neznamych pq, vektor prave
% strany b.
% pq(1)   = p(1)
% pq(2)   = q(1)
% pq(3)   = p(2)
% ...
% pq(2*i-1) = p(i)
% pq(2*i)   = q(i)
% ...
% rovnice jsou usporadany takto:
% prvnich nbc radku zadane okrajove podminky (tlaky a prutoky)
% nasledujici radky (nbc+1,...,nbc+3.nev) popisuji trojice rovnic vetveni 
%     q1=q2+q3
%     p1-p2=... Bernouliho rovnice
%     p1-p3=...
% nasledujici radky (nbc+3.nev+1,...,nbc+3.nev+2.nep) popisuji dvojici rovnic pro trubku
%     q1=q2
%     p1-p2=... Bernoulliho rovnice
%
% definice rovnic pro okrajove podminky. Nejprve nulovani matice soustavy i
% prave strany
for i=1:neq
    for j=1:neq
        a(i,j)=0;
    end
    b(i)=0;
    pq(i)=0;
end
% prvnich nbc radku zadane okrajove podminky (tlaky a prutoky)
for i=1:nbc
    ibca=2*abs(ib(i));
    b(i)=bpq(i);
    if ib(i)<0   
        a(i,ibca)=1;
        pq(ibca)=b(i);
    else
        a(i,ibca-1)=1;
        pq(ibca-1)=b(i);
    end
end
%
% kapalina
mju=0.001
rho=1000
%
% vyhledani elementu, kde je zadany prutok. 
iq=0;
notfound=1;
while iq<nbc & notfound
    notfound=ib(iq+1)>0;
    iq=iq+1;
end
if notfound
    disp 'nenasel jsem zadny zadany prutok'
    return
else
    disp 'nasel jsem zadany prutok v uzlu',iq
end
% iq-index uzlu prvni okrajove podminky, kde je zadan prutok
% hledani elementu s uzlem iq
ie=0;
notfound=1;
while ie<nep & notfound
    ie=ie+1;
    notfound= iq ~= iup(ie,1) & iq ~= iup(ie,2);
end
if notfound
    disp 'nenasel jsem zadny element pipe se zadanym prutokem'
    return
else
    disp 'nasel jsem element pipe kde je zadany prutok',ie
end
% parametry trubky, a vypocet Reynoldsova cisla
d=dp(ie);
q=pq(2*iq);
re=4*q*rho/(3.14*d*mju)
%
% Iterace reseni pq
% =============================================================
iter=0;
normp=1e3;
normq=1e3;
while iter<10 & normp>1e-5 & normq>1e-5
    iter=iter+1;
    pqold=pq;
% generovani rovnic nbc+1,nbc+2,... soustavy pro elementy trubka
% ----------------------------------------------
ieq=nbc;
for ie=1:nev
    for j=1:3
        xu(j)=x(iuv(ie,j));
        yu(j)=y(iuv(ie,j));
        zu(j)=z(iuv(ie,j));
    end
    xc=x(iuc(ie));
    yc=y(iuc(ie));
    zc=z(iuc(ie));
    for j=1:3
        l(j)=sqrt((xu(j)-xc)^2+(yu(j)-yc)^2+(zu(j)-zc)^2);
    % indexy prutoku a tlaku
        iq(j)=2*iuv(ie,j);
        ip(j)=2*iuv(ie,j)-1;
    % prurezy a prutoky, turbulentni viskozita a kappa
        d(j)=dv(ie,j);
        s(j)=pi*dv(ie,j)^2/4;
        q(j)=pq(iq(j));
        v(j)=q(j)/s(j);
        [mjut(j),kappa(j),re(j)]=mjuf(mju,d(j),q(j),rho);
        rev(ie,j)=re(j);
    end
    % rovnice kontinuity
    a(ieq+1,iq(1))=1;
    a(ieq+1,iq(2))=-1;
    a(ieq+1,iq(3))=-1;
    
    % Bernoulliho rovnice 1 - 2, viz dokumentace
    
    a1frict=128/(rho*pi)*mjut(1)*l(1)/d(1)^4;
    a2frict=128/(rho*pi)*mjut(2)*l(2)/d(2)^4;
    
    a1=8/pi^2*(q(2)/d(2)^2+q(1)/d(1)^2)*kappa(1)/d(1)^2;
    a2=8/pi^2*(q(2)/d(2)^2+q(1)/d(1)^2)*kappa(2)/d(2)^2;
 
    zet2=zeta(d(2),d(1),q(2),q(1),alpha(ie));
    a2zeta=300*mju/(rho*pi*d(1)^3)+8*abs(q(1))/(pi^2*d(1)^4)*zet2;

    a(ieq+2,iq(2))=a2+a2frict;
    a(ieq+2,iq(1))=-a1+a1frict+a2zeta;
    a(ieq+2,ip(2))=1/rho;
    a(ieq+2,ip(1))=-1/rho;
    b(ieq+2)=g*(zu(2)-zu(1));
    
    % Bernoulliho rovnice 1 - 3, viz dokumentace
    
    a3frict=128/(rho*pi)*mjut(3)*l(3)/d(3)^4;

    a1=8/pi^2*(q(3)/d(3)^2+q(1)/d(1)^2)*kappa(1)/d(1)^2;
    a3=8/pi^2*(q(3)/d(3)^2+q(1)/d(1)^2)*kappa(3)/d(3)^2;
    
    zet3=zeta(d(3),d(1),q(3),q(1),alpha(ie));
    a3zeta=300*mju/(rho*pi*d(1)^3)+8*abs(q(1))/(pi^2*d(1)^4)*zet3;
    
    a(ieq+3,iq(3))=a3+a3frict;
    a(ieq+3,iq(1))=-a1+a1frict+a3zeta;
    a(ieq+3,ip(3))=1/rho;
    a(ieq+3,ip(1))=-1/rho;
    b(ieq+3)=g*(zu(3)-zu(1));
    
    ezv(ie,2)=(q(1)*a1frict+q(2)*a2frict+q(3)*a3frict)*rho; %vykon ve watech
    ezv(ie,1)=ezv(ie,2)+(a2zeta+a3zeta)*q(1)*rho;

    % disp 'fcoef',a1frict,a2frict,a3frict,a2zeta
    
    ieq=ieq+3;
end

%Elementy typu PIPE

for ie=1:nep
    for j=1:2
        xu(j)=x(iup(ie,j));
        yu(j)=y(iup(ie,j));
        zu(j)=z(iup(ie,j));
    end
    lpipe=sqrt((xu(2)-xu(1))^2+(yu(2)-yu(1))^2+(zu(2)-zu(1))^2);
    dpipe=dp(ie);
    
    for j=1:2
    % indexy prutoku a tlaku
        iq(j)=2*iup(ie,j);
        ip(j)=2*iup(ie,j)-1;
    % prutoky, turbulentni viskozita a kappa
        q(j)=pq(iq(j));
        [mjut(j),kappa(j),re(j)]=mjuf(mju,dpipe,q(j),rho);
    end
    rep(ie)=(re(1)+re(2))/2;
    % rovnice kontinuity
    a(ieq+1,iq(1))=1;
    a(ieq+1,iq(2))=-1;
    
    % Bernoulliho rovnice 1 - 2, viz dokumentace. Hagenbachuv faktor dle
    % Chen1972
    mjum=(mjut(1)+mjut(2))/2;
    a0=128*mjum*lpipe/(rho*pi*dpipe^4);
    if rep(ie)~=0
        ahagenbach=(1.2+38/rep(ie))*8*q(1)^2/(pi^2*dpipe^4);
    else
        ahagenbach=1.8*8*q(1)^2/(pi^2*dpipe^4);
    end
    a(ieq+2,iq(1))=a0+abs(ahagenbach);
    a(ieq+2,ip(2))=1/rho;
    a(ieq+2,ip(1))=-1/rho;
    b(ieq+2)=g*(zu(2)-zu(1));

    ezp(ie,2)=q(1)*a0*rho;        %vykon ve watech
    ezp(ie,1)=ezp(ie,2)+ahagenbach*q(1)*rho;

    ieq=ieq+2;
end

% reseni soustavy rovnic a.pg=b Gussovou eliminaci

pq=linsolve(a,b');

normp=0;
normq=0;
nump=0;
numq=0;
for i=1:nu
    pnew(i)=pq(2*i-1);
    qnew(i)=pq(2*i);
    pold(i)=pqold(2*i-1);
    qold(i)=pqold(2*i);
    if pnew(i)~=0
        normp=normp+abs((pnew(i)-pold(i))/pnew(i));
        nump=nump+1;
    end
    if qnew(i)~=0
        normq=normq+abs((qnew(i)-qold(i))/qnew(i));
        numq=numq+1;
    end
end
normp=normp/nump;
normq=normq/numq;
% disp 'normp,normq',iter,normp,normq
end
text=sprintf('pocet iteraci =%d',iter);
disp (text);
powerv=0;
powerp=0;
for i=1:nev
    powerv=powerv+ezv(i,1);
end
for i=1:nep
    powerp=powerp+ezp(i,1);
end
pmin=min(pnew);
pmax=max(pnew);
text=sprintf('celkovy dissipovany vykon [W]=%0.3g (%0.3g+%0.3g)',powerv+powerp,powerv,powerp);
disp (text);
text=sprintf('Rozsah tlaku [Pa] od %0.3g do %0.3g',pmin,pmax);
disp (text);
disp 'Podivej se na rev (Reynolds ve vetveni) a rep (Re v trubkach)'
disp 'Podivej se i na ezv (W ztrat ve vetveni) a ezp (v trubkach)'









    




