function programbakalarka
clc;
%funguje pouze v kombinaci s knihovnou XSteam
%vstupy:
p_out=1.15;%bar výstupní tlak
qt=7000;%W tepelný výkon
L=0.326;%m delka potrubí
d=0.009;%m prùmìry
D=0.014;%m
T_in=95;%°C vstupní teplota
Eps=10^-5;
k=1.6*10^-6; %m  prùmìrná hodnota absolutní drsnosti
%_____________________________________
S_p=pi*(D^2-d^2)/4;   %prùtoèná plocha
O=pi*(D+d);           %smáèený obvod
De=4*S_p/O;           %hydraulický prùmìr
h_in=XSteam('h_pT',p_out,T_in)*10^3; %J/kg vstupní entalpie
ro_l=XSteam('rho_pT',p_out,T_in); %kg/m^3 hustota kapaliny
%________________________________________
qd=qt/L;
G=0;
A=zeros(2,2);
for i=1:100         %rùzné prùtoky
    G=G+80;
    mt=G*S_p;
    h_out=(qt/mt)+h_in; %J/kg
    A(1,i)=G;
    A(2,i)=h_out;
    T_out=XSteam('T_ph',p_out,h_out*10^(-3));
    T_sat=XSteam('Tsat_p',p_out); %výstpní teplota
    %_______________________________________________________________
    if  T_out < T_sat  %nedojde k varu
        p_z=0;         %tlaková ztráta zrychlením
        p_grav=9.81*ro_l*L;%tlaková ztráta gravitací
        h_sl=(h_in+h_out)*10^(-3)/2;
        my=XSteam('my_ph',p_out,h_sl); %Pas
        Re=G*De/my;
        f=0.184*(Re^(-0.2));
        Lambda_old=f;
        rozdil=1;
        while (rozdil>Eps)%výpoèet koeficientu tøení
            Lambda_new=(1/(-2*log10((k/(De*3.7))+(2.51/(Re*sqrt(Lambda_old))))))^2;
            rozdil=Lambda_new-Lambda_old;
            Lambda_old=Lambda_new;
        end
        f_lo=Lambda_old;
        p_t=(f_lo*(G^2)*L)/(2*De*ro_l);%tlaková ztráta tøením
        %_____________________________________________________________
        pressuredrop=p_z+p_grav+p_t;
        A(3,i)=pressuredrop;
        A(4,i)=p_z;
        A(5,i)=p_grav;
        A(6,i)=p_t;
        %_________________________________________________________________________
    else if T_out==T_sat %dojde k varu
            h_f=XSteam('hL_T',T_sat)*10^3; %J/kg
            h_g=XSteam('hV_T',T_sat)*10^3; %J/kg
            h_fg=h_g-h_f;
            ro_g=XSteam('rhoV_T',T_sat);    %kg/m^3 hustota
            ro_f=XSteam('rhoL_T',T_sat);
            
            z_sat=(h_f-h_in)*mt/qd;    %souøadnice varu
            h=@(z) h_in+qd*z/mt;
            x_e=@(z) ((h_in-h_f)/(h_fg))+((qd*z)/(mt*h_fg));
            beta=@(z) 1/(1+((1-x_e(z))/x_e(z))*(ro_g/ro_f));
            %==============================================================
            %vsuvka pro vykreslení entalpie, suchosti a podílu páry
            %figure;
            %subplot(3,1,1)
            %fplot(h,[0,L]);
            %title('entalpie');
            %xlabel('L [m]');
            %ylabel('h [-]')
            %subplot(3,1,2);
            %fplot(x_e,[0,L]);
            %title('rovnovážná suchost');
            %xlabel('L [m]');
            %ylabel('x_e [-]')
            %subplot(3,1,3);
            %fplot(beta,[z_osb,L]);
            %xlim([0 L]);
            %title('podíl páry');
            %xlabel('L [m]');
            %ylabel('Beta [-]')
            %==========================================================================
            %____________________________________________________________________
            p_z=(((1-x_e(L))^2/((1-beta(L))*ro_f))+((x_e(L))^2/(beta(L)*ro_g))-(1/ro_l))*G^2;%tlaková ztráta zrychlením
            %nebo
            %ro_out=XSteam('rho_ph',p_out,h_out*10^-3);
            %p_z=(1/ro_out-1/ro_l)*G^2;
            %________________________________________________________________
            xa=(x_e(L)-x_e(0))/L;
            xd=x_e(0)+(ro_g/ro_f)*(1-x_e(0));
            xdd=xa*(1-(ro_g/ro_f));
            p_grav=9.81*(ro_l*z_sat+ro_f*(L-z_sat)-(ro_f-ro_g)*((x_e(0)/xdd)*log((xd+xdd*L)/(xd+xdd*z_sat))+(xa/xdd)*(L-z_sat)-xa*xd/xdd^2*log((xd+xdd*L)/(xd+xdd*z_sat))));
            %tlaková ztráta gravitací
            %_____________________________________________________________
            h_sl=(h_in+h_f)*10^(-3)/2;
            my=XSteam('my_ph',p_out,h_sl); %Pas
            Re=G*De/my;
            f=0.184*(Re^(-0.2));
            Lambda_old=f;
            rozdil=1;
            while (rozdil>Eps)
                Lambda_new=(1/(-2*log10((k/(De*3.7))+(2.51/(Re*sqrt(Lambda_old))))))^2; %friction factor 1 phase (we assume turbulent flow)
                rozdil=Lambda_new-Lambda_old;
                Lambda_old=Lambda_new;
            end
            f_lo=Lambda_old;
            p_t=(f_lo*G^2)/(2*De*ro_l)*(L+((ro_f/ro_g)-1)*(x_e(0)*(L-z_sat)+0.5*xa*(L^2-(z_sat)^2)));
            %tlaková ztráta tøením
            %_____________________________________________________________
            pressuredrop=p_z+p_grav+p_t;
            %p_in=p_out+p_drop
            A(3,i)=pressuredrop;
            A(4,i)=p_z;
            A(5,i)=p_grav;
            A(6,i)=p_t;
            %__________________________________________________________________________
        else if T_out>T_sat %všechna kapalina se vypaøí
                h_f=XSteam('hL_T',T_sat)*10^3; %J/kg
                h_g=XSteam('hV_T',T_sat)*10^3; %J/kg
                h_fg=h_g-h_f;
                ro_g=XSteam('rhoV_T',T_sat);    %kg/m^3 hustota
                ro_f=XSteam('rhoL_T',T_sat);
                ro_vout=XSteam('rho_ph',p_out,h_out*10^(-3));
                z_sat=(h_f-h_in)*mt/qd;%souøadnice varu
                z_v=(h_g-h_in)*mt/qd;%souøadnice, kde zùstane jen pára
                h=@(z) h_in+qd*z/mt;
                x_e=@(z) ((h_in-h_f)/(h_fg))+((qd*z)/(mt*h_fg));
                beta=@(z) 1/(1+((1-x_e(z))/x_e(z))*(ro_g/ro_f));
                %____________________________________________________________________
                p_z=((1/ro_vout)-(1/ro_l))*G^2; %tlaková ztráta zrychlením
                %________________________________________________________________
                xa=(x_e(L)-x_e(0))/L;
                xd=x_e(0)+(ro_g/ro_f)*(1-x_e(0));
                xdd=xa*(1-(ro_g/ro_f));
                p_grav=9.81*(ro_l*z_sat+ro_f*(z_v-z_sat)-(ro_f-ro_g)*((x_e(0)/xdd)*log((xd+xdd*z_v)/(xd+xdd*z_sat))+(xa/xdd)*(z_v-z_sat)-xa*xd/(xdd^2)*log((xd+xdd*z_v)/(xd+xdd*z_sat)))+ro_vout*(L-z_v));
                %tlaková ztráta gravitací
                %_____________________________________________________________
                h_sl=(h_in+h_f)*10^(-3)/2;
                my=XSteam('my_ph',p_out,h_sl); %Pas
                h_sv=(h_g+h_out)*10^(-3)/2;
                my2=XSteam('my_ph',p_out,h_sv); %Pas
                Re=G*De/my;
                Re2=G*De/my2;
                f=0.184*(Re^(-0.2));
                Lambda_old=f;
                rozdil=1;
                while (rozdil>Eps)
                    Lambda_new=(1/(-2*log10((k/(De*3.7))+(2.51/(Re*sqrt(Lambda_old))))))^2; %friction factor 1 phase (we assume turbulent flow)
                    rozdil=Lambda_new-Lambda_old;
                    Lambda_old=Lambda_new;
                end
                f_lo=Lambda_old;
                %___________________________
                f2=0.184*(Re2^(-0.2));
                Lambda2_old=f2;
                rozdil2=1;
                while (rozdil2>Eps)
                    Lambda2_new=(1/(-2*log10((k/(De*3.7))+(2.51/(Re2*sqrt(Lambda2_old))))))^2; %friction factor 1 phase (we assume turbulent flow)
                    rozdil2=Lambda2_new-Lambda2_old;
                    Lambda2_old=Lambda2_new;
                end
                f_vo=Lambda2_old;
                p_t=(f_lo*G^2)/(2*De*ro_l)*(z_v+((ro_f/ro_g)-1)*(x_e(0)*(z_v-z_sat)+0.5*xa*(z_v^2-(z_sat)^2)))+(f_vo*G^2)/(2*De*ro_vout)*(L-z_v);
                %tlaková ztráta tøením
                %_____________________________________________________________
                pressuredrop=p_z+p_grav+p_t;
                A(3,i)=pressuredrop;
                A(4,i)=p_z;
                A(5,i)=p_grav;
                A(6,i)=p_t;
            end
        end
    end
    %============================================================================
    %pøidání clony
    ksi=21;
    p_c=ksi*0.5*G^2/ro_l;%tlaková ztráta clony
    A(7,i)=p_c;
    A(8,i)=A(3,i)+A(7,i);%výpoèet nové tlakové ztráty po pøidání clony
end %for prutoky
format short g
A
%==========================================================================
%dosazení køivky èerpadla
X=[0.104 0.313 0.649 0.808 1 1.27 1.59 1.88 1.99 2.18 2.44 2.67];%data z internetu
for i=1:12
    X(i)=X(i)*898/(S_p*3600);
end
Y=[2097000 2112000 2106000 2090000 2058000 1994000 1885000 1756000 1699000 1592000 1420000 1253000];
podiln=0.25;
for i=1:12
    Y(i)=Y(i)*podiln^3;
end
p=polyfit(X,Y,3); % koeficienty polynomu 3. stupne pro x a y
v_aprox=polyval(p,X); % hodnoty polynomu v x
%=========================================================================
%vykreslení grafu
figure
title(['tlaková ztráta pro výkon q=' num2str(qt) 'W']);
hold on
xlabel('hmotnostní prùtok [kg/m^2s]');
ylabel('tlaková ztráta [Pa]')
plot(A(1,:),A(3,:),'b');
plot(A(1,:),A(4,:),'r.');
plot(A(1,:),A(5,:),'y.');
plot(A(1,:),A(6,:),'g.');
plot(X,v_aprox,'m');
plot(A(1,:),A(8,:),'c');
ylim([0,85000]);
legend('\Delta p_{int,1}','\Delta p_Z','\Delta p_G','\Delta p_T','\Delta p_{ext}','\Delta p_{int,2}')
end



