close all
clear all
clc

i=1;
teploty=1:0.1:80;
lambda=[1 1.37 2];

for k=1:length(lambda)
    lam=lambda(k);
    prvnipruchod=1;
    for j=1:length(teploty)
        t=teploty(j);
        cs=1.1;         %měrná tepelná kapacita teoretických suchých spalin = 1,1 [kJ/(kg*K)]
        ca=1;           %měrná tepelná kapacita suchého vzduchuĐ [kJ/(kg*K)]
        a= 1.077;       %součinitel teoretické potřeby vzduchu [-]
        
        x=0.14;         %měrná vlhkost spalin [kg/kg t.s.s.]
        p=98;           %atmosferický tlak [kPa]
        
        %parciální tlak páry na mezi sytosti [kPa] - pro teploty 0 - 80 °C
        pp=10^-3*exp(23.58-4044.2./(235.6+t));
        %[3] Chyský, J., Hemzal,K. a kol.: Technický průvodce "Větrání a klimatizace". Bolit, 1993
        
        %parciální tlak páry v závislosti na měrné vlhkosti spalin
        ppx=p/(0.622/(x*(1+a*(lam-1)))+1); %x=0.14
        
        % měrná vlhkost spalin x pro teplotu pod rosným bodem nebude konstantní
        xn=0.622*(1+a*(lam-1))*pp/(p-pp);
        
        %teplota rosného bodu pro je dána průnikem průběhů těchto entalpií
        EntalpieNadRosnymBodem=(cs + a*(lam-1)*ca)*t+x*(2500+1.84*t);
        EntalpiePodRosnymBodem=(cs + a*(lam-1)*ca)*t+xn*(2500+1.84*t);
        
        if EntalpieNadRosnymBodem <= EntalpiePodRosnymBodem
            %Vztah pro entalpii spalin při teplotě nad rosným bodem
            if (prvnipruchod)
                RosnyBodT(i)=t;
                RosnyBodE(i)=EntalpieNadRosnymBodem;
                prvnipruchod=0;
            end
            Entalpie(i,j)=EntalpieNadRosnymBodem;
            
        else
            %Vztah pro entalpii spalin při teplotě pod rosným bodem
            Entalpie(i,j)=EntalpiePodRosnymBodem;
        end
    end
    plot(teploty,Entalpie(i,:)')
    hold on
    plot(RosnyBodT,RosnyBodE,'o')
    
    i=i+1;
    
end
legend('lambda = 1',['rosny bod = ', num2str(RosnyBodT(1)) ' °C'] , 'lambda = 1.37',['rosny bod = ' num2str(RosnyBodT(2)) ' °C'] , 'lambda = 2',['rosny bod = ' num2str(RosnyBodT(3)) ' °C'])
xlabel('teplota [°C]')
ylabel('entalpie [kJ/kg t.s.s.]')
grid on