clear all
clc

%  Toto je základní skript pro výpočet tepelně integrované Carnotovy baterie - model se semiempirickym modelem kompresoru a expanderu
%  
%  
%  Všechny citlivostní analýzi pro teplý model byly vytvořeny modfikací tohoto skriptu
%
%  Pro správné fungování modelu je potřeba nejdříve stáhnout a nainstalovat v matlabu programovací prostředí Python 
%  
%  
%  podrobny navod instalace pythonu i refpropu lze nalezt na teto adrese: https://github.com/usnistgov/REFPROP-wrappers/tree/master/wrappers/MATLAB
%   
%
%  pro propojeni matlabu s refpropem je nutne zadat nasledujici kod : [v,e] = pyversion; system(['"', e, '" -m pip install --user -U ctREFPROP'])
%  
%
%  pro propojeni matlabu s coolpropem je nutne zadat nasledujici kod : [v,e] = pyversion; system([e,' -m pip install --user -U CoolProp'])
%



%         Vypocet semi-empirickeho modelu kompresoru a expanderu je ve forme funkce uveden na konci tohoto skriptu
%         






%%%%%%%%%%%%%% tyto řádky zajištují propojení matlabu a refpropu %%%%%%%%%%

global RP
global MOLSI
global iMass
global iFlag

RP = py.ctREFPROP.ctREFPROP.REFPROPFunctionLibrary('C:\Program Files (x86)\REFPROP');
MOLSI = RP.GETENUMdll(int8(0),'SI WITH C').iEnum;
iMass = int8(0); % 1: molar fractions; 1: mass fractions
iFlag = int8(0); % 0: don't call SATSPLN; 1: call SATSPL
clc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

grafy = 0; %% pokud 1 vytvori vypocet grafy
vysledky = 1; %% pokud 1 vytvori vyskedkove veliciny
tabulky = 1; %% pokud 1 vytvori vysledkovou tabulku jednotlivych stavu

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                       pocatecni vstupy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %% parametry rekauperace
    
    rekuperace_TC = 2; %% pokud je 2 = s rekuperací, 1 = bez rekuperace
    rekuperace_ORC = 2; %% pokud je 2 = s rekuperací, 1 = bez rekuperace 
    PP_rek = 5; %% K

    %% parametry nabíjení
    
    fluid_zdroj = 'WATER';      %  zdroj tepla - charakter
    t_zdroj_in = 70; %% C      teplota zdroje tepla na vstupu do TČ
    t_zdroj_out = 60; %% C      teplota zdroje tepla na výstupu z TČ
    p_zdroj = 0.1; %% Pa        tlak zdroje teopa
    fluid_working ='R1233ZDE';     %% pracovní látka tepelného čerpadla
    PP_cond = 5; %% K           Pinchpoint kondenzátor TČ
    PP_evap = 5; %% K           Pinchpoint výparník TČ
    eta_comp = 0.7;          %  iso s účinnost komprese

    %% parametry vybijení
    
    fluid_working_dis ='R1233ZDE';     %% pracovní látka ORC
    fluid_chlazeni ='AIR';             %% zdroj chlazení ORC
    p_chlazeni = 0.101325; % Pa        %% tlak zdroje chladu
    t_chladivo_in = 20; %% °C          %% teplota zdroje chladu na vstupu
    t_chladivo_out = 30; %% C          %% teplota zdroje chladu na výstupu
    prehrati = 10;                     %% přehřátí pracovní látka TČ    
    eta_expander = 0.70;               %% účinost expanderu
    eta_pump = 0.70;                   %% účinnost nap. čerpadla
    PP_cond_discharg = 5; %% K         %% PP - kondenzátor ORC
    PP_evap_discharg = 5; %% K         %% PP - výparník ORC
     
    %% zasobnik tepla
    
    fluid_zasobnik ='INCOMP::T66';      %% látka uložiště tepla
    t_dis = 80; % C;                    teplota vybitého uložiště
    t_charg = 90; % C                   teplota nabitého uložiště
    p_zasobnik = 0.1; %% MPa;           tlak uložiště tepla
  

%%% od tohoto mista neni nutne upravovat zadny parametr pro spravne fungovani 

    
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% NABIJENI VYPOCET
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% NABIJENI BEZ REKUPERACE

if rekuperace_TC == 1

    t_HS1 = t_zdroj_in;
    t_HS2 = t_zdroj_out;
    h_HS1 = stav(fluid_zdroj,'TP','H',t_HS1,p_zdroj);
    h_HS2 = stav(fluid_zdroj,'TP','H',t_HS2,p_zdroj);

    h_ST1 = py.CoolProp.CoolProp.PropsSI('H','T',t_dis+273.15,'P', p_zasobnik*1000000,fluid_zasobnik)/1000;    
    h_ST3 = py.CoolProp.CoolProp.PropsSI('H','T',t_charg+273.15,'P', p_zasobnik*1000000,fluid_zasobnik)/1000;
    
    t_ST1 = t_dis;
    t_ST3 = t_charg;

    t1 = t_HS2-PP_evap;
    t2 = t1;
    t3 = t2+prehrati;
    p2 = stav(fluid_working,'TQ','P',t2,1);
    p3 = p2;
    p1 = p2;
    h2 = stav(fluid_working,'TQ','H',t2,1);
    h3 = stav(fluid_working,'TP','H',t3,p3);
    s2 = stav(fluid_working,'TQ','S',t2,1);
    s3 = stav(fluid_working,'TP','S',t3,p3);
    t7=t_ST1+PP_cond;
    t_condcharg=t_ST3+PP_cond;

        chyba = 1;

        while abs(chyba) >= 0.01

            if chyba > 0
                t_condcharg=t_condcharg-(abs(chyba)/2);
            else 
                t_condcharg=t_condcharg+(abs(chyba)/2);
            end

            t6=t_condcharg;
            t5=t_condcharg;
            p7 = stav(fluid_working,'TQ','P',t5,1);
            p4 = stav(fluid_working,'TQ','P',t5,1);
            p6 = stav(fluid_working,'TQ','P',t5,1);
            p5 = stav(fluid_working,'TQ','P',t5,1);
            h7 = stav(fluid_working,'TP','H',t7,p7);
            h1=h7; 
            %m_wf_charg=(m_HS*(h_HS1-h_HS2))/(h3-h1);
            [t4, eta_kompresor, m_wf_charg, P_kompresor, KOMPRESOR] = kompresor_fce(t6, t3, t2);
            
            m_HS = (m_wf_charg*(h3-h1))/(h_HS1-h_HS2);
            s1 = stav(fluid_working,'PH','S',p1,h1);
            h5 = stav(fluid_working,'TQ','H',t5,1);
            h6 = stav(fluid_working,'TQ','H',t6,0);
            s5 = stav(fluid_working,'TQ','S',t5,1);
            s6 = stav(fluid_working,'TQ','S',t6,0);
            s7 = stav(fluid_working,'TP','S',t7,p7);
            s4is = s3;
            h4is = stav(fluid_working,'PS','H',p4,s4is);
            %h4 = ((h4is-h3)/eta_comp)+h3;
            %t4 = stav(fluid_working,'PH','T',p4,h4);
            h4 = stav(fluid_working,'PT','H',p4,t4);
            s4 = stav(fluid_working,'HP','S',h4,p4);
            m_ST_charg = (m_wf_charg*(h4-h7))/(h_ST3-h_ST1);
            h_ST2 = ((m_wf_charg*(h5-h7))+(h_ST1*m_ST_charg))/m_ST_charg;
                        
            t_ST2 = py.CoolProp.CoolProp.PropsSI('T','H',h_ST2*1000,'P', p_zasobnik*1000000,fluid_zasobnik)-273.15;
            chyba = t5-(t_ST2+PP_cond);
            t_prehrati = t4-t5;

        end
    
Q_in_char = m_HS*(h_HS1-h_HS2);
Q_out_char = m_ST_charg*(h_ST3-h_ST1);
W_comp_char = m_wf_charg*(h4-h3);
COP=Q_out_char/W_comp_char;
COP_el = Q_out_char/P_kompresor;
m_wf_charg;
m_HS;
m_ST_charg;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% KONEC VYPOCTU BEZ REKUPERACE




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% NABIJENI S REKUPERACÍ

if rekuperace_TC == 2

    t_HS1 = t_zdroj_in;
    t_HS2 = t_zdroj_out;
    h_HS1 = stav(fluid_zdroj,'TP','H',t_HS1,p_zdroj);
    h_HS2 = stav(fluid_zdroj,'TP','H',t_HS2,p_zdroj);

    
    h_ST1 = py.CoolProp.CoolProp.PropsSI('H','T',t_dis+273.15,'P', p_zasobnik*1000000,fluid_zasobnik)/1000;    
    h_ST3 = py.CoolProp.CoolProp.PropsSI('H','T',t_charg+273.15,'P', p_zasobnik*1000000,fluid_zasobnik)/1000;    
    t_ST1 = t_dis;
    t_ST3 = t_charg;

    t1 = t_HS2-PP_evap;
    t2 = t1;
    t3 = t2+prehrati;
    p2 = stav(fluid_working,'TQ','P',t2,1);
    p3 = p2;
    p3rek = p3;
    p1 = p2;
    h2 = stav(fluid_working,'TQ','H',t2,1);
    h3 = stav(fluid_working,'TP','H',t3,p3);
    s2 = stav(fluid_working,'TQ','S',t2,1);
    s3 = stav(fluid_working,'TP','S',t3,p3);
    t7=t_ST1+PP_cond;
    t3rek = t7-PP_rek;
    h3rek = stav(fluid_working,'TP','H',t3rek,p3rek);
    s3rek = stav(fluid_working,'TP','S',t3rek,p3rek);
    t_condcharg=t_ST3+PP_cond;
    Q_rek_nab = h3rek-h3;

    %% iterace na PP

        chyba = 1;

    while abs(chyba) >= 0.01
        if chyba > 0
            t_condcharg=t_condcharg-(abs(chyba)/2);
        else 
            t_condcharg=t_condcharg+(abs(chyba)/2);
        end

        t6=t_condcharg;
        t5=t_condcharg;
        p7 = stav(fluid_working,'TQ','P',t5,1);
        p7rek = p7;
        p4 = stav(fluid_working,'TQ','P',t5,1);
        p6 = stav(fluid_working,'TQ','P',t5,1);
        p5 = stav(fluid_working,'TQ','P',t5,1);
        h7 = stav(fluid_working,'TP','H',t7,p7);
        h7rek = h7-Q_rek_nab;
        h1=h7rek;
        s7rek = stav(fluid_working,'HP','S',h7rek,p7rek);
        t7rek = stav(fluid_working,'HP','T',h7rek,p7rek);

        [t4, eta_kompresor, m_wf_charg, P_kompresor, KOMPRESOR] = kompresor_fce(t6, t3rek, t2);
        
        %m_wf_charg=(m_HS*(h_HS1-h_HS2))/(h3-h1);
        m_HS = (m_wf_charg*(h3-h1))/(h_HS1-h_HS2);
        s1 = stav(fluid_working,'PH','S',p1,h1);
        h5 = stav(fluid_working,'TQ','H',t5,1);
        h6 = stav(fluid_working,'TQ','H',t6,0);
        s5 = stav(fluid_working,'TQ','S',t5,1);
        s6 = stav(fluid_working,'TQ','S',t6,0);
        s7 = stav(fluid_working,'TP','S',t7,p7);
        s4is = s3rek;
        h4is = stav(fluid_working,'PS','H',p4,s4is);
        %h4 = ((h4is-h3rek)/eta_comp)+h3rek;
        %t4 = stav(fluid_working,'PH','T',p4,h4);
        
        h4 = stav(fluid_working,'PT','H',p4,t4);
        s4 = stav(fluid_working,'HP','S',h4,p4);
        m_ST_charg = (m_wf_charg*(h4-h7))/(h_ST3-h_ST1);
        h_ST2 = ((m_wf_charg*(h5-h7))+(h_ST1*m_ST_charg))/m_ST_charg;
        t_ST2 = py.CoolProp.CoolProp.PropsSI('T','H',h_ST2*1000,'P', p_zasobnik*1000000,fluid_zasobnik)-273.15;
        chyba = t5-(t_ST2+PP_cond);
        t_prehrati = t4-t5;

    end

    Q_in_char = m_HS*(h_HS1-h_HS2);
    Q_out_char = m_ST_charg*(h_ST3-h_ST1);
    W_comp_char = m_wf_charg*(h4-h3rek);
    COP=Q_out_char/W_comp_char;
    COP_el = Q_out_char/P_kompresor
    m_wf_charg;
    m_HS;
    m_ST_charg;

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% KONEC VYPOCTU S REKUPERACÍ

















%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% VYBIJENI VYPOCET
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% VYBÍJENÍ VÝPOČET: BEZ REKUPERACE
if rekuperace_ORC == 1

    t_CL1 = t_chladivo_in;
    h_CL1 = stav(fluid_chlazeni, 'TP', 'H', t_CL1, p_chlazeni);
    t_CL3 = t_chladivo_out;
    h_CL3 = stav(fluid_chlazeni, 'TP', 'H', t_CL3, p_chlazeni);

    t14 = (t_CL3-0.5)+PP_cond_discharg;
    
    chyba3=1;

    while abs(chyba3) > 0.01
        if chyba3 > 0
            t14=t14-(abs(chyba3)/1.1);
        else 
            t14=t14+(abs(chyba3)/1.1);
        end


        t8 = t14;
        p8 = stav(fluid_working_dis, 'TQ', 'P', t8, 0);
        p13 = p8;
        p14 = p8;
        h14 = stav(fluid_working_dis, 'TQ', 'H', t14, 1);
        s14 = stav(fluid_working_dis, 'TQ', 'S', t14, 1);

        chyba2 = 1;
        t10=t_ST1-2;

        while abs(chyba2) > 0.01
            if chyba2 > 0
                t10=t10+((abs(chyba2))/1.1);
            else
                t10=t10-((abs(chyba2))/1.1);
            end

            t11=t10;
            p10 = stav(fluid_working_dis, 'TQ', 'P', t10, 0);
            p9=p10;
            p11=p10;
            p12=p10;
            h11 = stav(fluid_working_dis, 'PQ', 'H', p11, 1);
            s11 = stav(fluid_working_dis, 'PQ', 'S', p11, 1);
            h10 = stav(fluid_working_dis, 'PQ', 'H', p10, 0);
            s10 = stav(fluid_working_dis, 'PQ', 'S', p10, 0);
            t12=t_ST3-PP_evap_discharg;
            h12 = stav(fluid_working_dis, 'TP', 'H', t12, p12);
            h8 = stav(fluid_working_dis, 'PQ', 'H', p8, 0);
            s8 = stav(fluid_working_dis, 'PQ', 'S', p8, 0);
            s9is = s8;   
            h9is = stav(fluid_working_dis, 'PS', 'H', p9, s9is);
            h9 = ((h9is-h8)/eta_comp)+h8;    
            t9 = stav(fluid_working_dis, 'PH', 'T', p9, h9);
            s9 = stav(fluid_working_dis, 'PH', 'S', p9, h9);
            
            [P_expander, eta_expanze, m_wf_discharg, t13, EXPANDER] = expander_fce(t11, t12, t14);
            m_ST_discharg = m_wf_discharg*(h12-h9)/(h_ST3-h_ST1);
            %m_wf_discharg = (m_ST_discharg*(h_ST3-h_ST1))/(h12-h9);

            
            
            
            h_ST4 = ((m_wf_discharg*(h10-h9))+(h_ST1*m_ST_discharg))/m_ST_discharg;
            t_ST4 = py.CoolProp.CoolProp.PropsSI('T','H',h_ST4*1000,'P', p_zasobnik*1000000,fluid_zasobnik)-273.15;
            chyba2=t_ST4-(t10+PP_evap_discharg);            
            
            
        end

        s12= stav(fluid_working_dis, 'TP', 'S', t12, p12);
        s13is=s12;
        h13is = stav(fluid_working_dis, 'PS', 'H', p13, s13is);
        %h13 = h12-((h12-h13is)*eta_expander);
        %t13 = stav(fluid_working_dis, 'PH', 'T', p13, h13);
        h13 = stav(fluid_working_dis, 'PT', 'H', p13, t13);
        s13 = stav(fluid_working_dis, 'PH', 'S', p13, h13);
        m_CL = ((h13-h8)*m_wf_discharg)/(h_CL3-h_CL1);
        h_CL2 = ((m_wf_discharg*(h14-h8))+(h_CL1*m_CL))/m_CL;
        t_CL2 = stav(fluid_chlazeni, 'PH', 'T', p_chlazeni, h_CL2);

        chyba3=(t14-PP_cond_discharg)-t_CL2;

    end

    Q_in_dis=m_wf_discharg*(h12-h9);
    Q_out_dis=m_wf_discharg*(h13-h8);
    W_pump=m_wf_discharg*(h9-h8);
    W_exp=m_wf_discharg*(h12-h13);
    W_cycle=W_exp-W_pump;
    eta_ORC=W_cycle/Q_in_dis;
    eta_ORC_el = (P_expander-W_pump)/Q_in_dis;
    eta_rt=eta_ORC*COP;
    eta_rt_el = eta_ORC_el*COP_el;
    m_ST_discharg;
    m_wf_discharg;
    m_CL;
  
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% KONEC VYBIJENI BEZ REKUPERACE



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% VYBÍJENÍ VÝPOČET S REKUPERACÍ

if rekuperace_ORC == 2

    t_CL1 = t_chladivo_in;
    h_CL1 = stav(fluid_chlazeni, 'TP', 'H', t_CL1, p_chlazeni);
    t_CL3 = t_chladivo_out;
    h_CL3 = stav(fluid_chlazeni, 'TP', 'H', t_CL3, p_chlazeni);

    t14 = (t_CL3-0.5)+PP_cond_discharg;

    
    chyba3=1;

    while abs(chyba3) > 0.01
        if chyba3 > 0
        t14=t14-(abs(chyba3)/1.1);
        else 
        t14=t14+(abs(chyba3)/1.1);
        end

        t8 = t14;
        p8 = stav(fluid_working_dis, 'TQ', 'P', t8, 0);
        p13 = p8;
        p13rek=p8;
        p14=p8;
        h14 = stav(fluid_working_dis, 'TQ', 'H', t14, 1);
        s14 = stav(fluid_working_dis, 'TQ', 'S', t14, 1);

        chyba2 = 1;
        t10=t_ST1-2;

        while abs(chyba2) > 0.01
            if chyba2 > 0
                t10=t10+((abs(chyba2))/1.1);
            else
                t10=t10-((abs(chyba2))/1.1);
            end

            t11=t10;
            p10 = stav(fluid_working_dis, 'TQ', 'P', t10, 0);
            p9=p10;
            p9rek=p10;
            p11=p10;
            p12=p10;
            h11 = stav(fluid_working_dis, 'PQ', 'H', p11, 1);
            s11 = stav(fluid_working_dis, 'PQ', 'S', p11, 1);
            h10 = stav(fluid_working_dis, 'PQ', 'H', p10, 0);
            s10 = stav(fluid_working_dis, 'PQ', 'S', p10, 0);
            t12=t_ST3-PP_evap_discharg;
            h12 = stav(fluid_working_dis, 'TP', 'H', t12, p12);
            h8 = stav(fluid_working_dis, 'PQ', 'H', p8, 0);
            s8 = stav(fluid_working_dis, 'PQ', 'S', p8, 0);
            s9is = s8;
            h9is = stav(fluid_working_dis, 'SP', 'H', s9is, p9);
            h9 = ((h9is-h8)/eta_comp)+h8;
            t9 = stav(fluid_working_dis, 'PH', 'T', p9, h9);
            s9 = stav(fluid_working_dis, 'PH', 'S', p9, h9);
            t13rek = t9+PP_rek;
            h13rek = stav(fluid_working_dis, 'TP', 'H', t13rek, p13rek);
            s12 = stav(fluid_working_dis, 'TP', 'S', t12, p12);
            s13is=s12;
            h13is = stav(fluid_working_dis, 'PS', 'H', p13, s13is);
            [P_expander, eta_expanze, m_wf_discharg, t13, EXPANDER] = expander_fce(t11, t12, t14);            
            m_ST_discharg = m_wf_discharg*(h12-h9)/(h_ST3-h_ST1);
            h13 = stav(fluid_working_dis, 'PT', 'H', p13, t13);
            
            %h13 = h12-((h12-h13is)*eta_expander);
            h9rek = h9+(h13 - h13rek);
            t9rek = stav(fluid_working_dis, 'PH', 'T', p9rek, h9rek);


            
            %m_wf_discharg = (m_ST_discharg*(h_ST3-h_ST1))/(h12-h9rek);
            h_ST4 = ((m_wf_discharg*(h10-h9rek))+(h_ST1*m_ST_discharg))/m_ST_discharg;
            t_ST4 = py.CoolProp.CoolProp.PropsSI('T','H',h_ST4*1000,'P', p_zasobnik*1000000,fluid_zasobnik)-273.15;
            chyba2=t_ST4-(t10+PP_evap_discharg);

        end

        %t13 = stav(fluid_working_dis, 'PH', 'T', p13, h13);
        s13 = stav(fluid_working_dis, 'PH', 'S', p13, h13);
        s13rek = stav(fluid_working_dis, 'PH', 'S', p13rek, h13rek);
        s9rek = stav(fluid_working_dis, 'PH', 'S', p9rek, h9rek);
        m_CL = ((h13rek-h8)*m_wf_discharg)/(h_CL3-h_CL1);
        h_CL2 = ((m_wf_discharg*(h14-h8))+(h_CL1*m_CL))/m_CL;
        t_CL2 = stav(fluid_chlazeni, 'PH', 'T', p_chlazeni, h_CL2);

        chyba3=(t14-PP_cond_discharg)-t_CL2;

    end

    Q_in_dis=m_wf_discharg*(h12-h9rek);
    Q_out_dis=m_wf_discharg*(h13rek-h8);
    W_pump=m_wf_discharg*(h9-h8);
    W_exp=m_wf_discharg*(h12-h13);
    W_cycle=W_exp-W_pump;
    eta_ORC = W_cycle/Q_in_dis;
    eta_ORC_el = (P_expander-W_pump)/Q_in_dis;
    eta_rt=eta_ORC*COP;
    eta_rt_el = eta_ORC_el*COP_el;
    m_ST_discharg;
    m_wf_discharg;
    m_CL;
  
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% KONEC VYBIJENI S REKUPERACI

















%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if grafy == 1;
%%                       GENEROVANI GRAFŮ

%%%%%%%%%%%%%%%%%%%%% vytvoreni krivky sytosti pro nabijeni

a=0;
for i=(t1-10):(t6+10)
    a=a+1;
    temp_s1(a) = i;
    s1x0(a) = stav(fluid_working, 'TQ', 'S', i, 0);
    s1x1(a) = stav(fluid_working, 'TQ', 'S', i, 1);
end  

%%%%%%%%%%%%%%%%%%%%% vytvoření křivky sytosti pro vybíjení

b=0;
for j=(t8-10):(t12+10)
    b=b+1;
    temp_s2(b) = j;
    s2x0(b) = stav(fluid_working_dis, 'TQ', 'S', j, 0);
    s2x1(b) = stav(fluid_working_dis, 'TQ', 'S', j, 1);
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                         GRAFY NABIJENI


if rekuperace_TC == 1
    
        S1=[s1 , s2 , s3 , s4 , s5 , s6 , s7 , s1];
        T1=[t1 , t2 , t3 , t4 , t5 , t6 , t7 , t1];
        str1_1={'1' , '2' , '3' , '4' , '5' , '6' , '7' , '1'};

        S_ST1=[s7 , s5 , s4];
        T_ST1=[t_ST1 , t_ST2, t_ST3];
        str1_3={ 'ST_1' , 'ST_2' , 'ST_3'};

    S_HS=[s1 , s3];
    T_HS=[t_HS2, t_HS1];
    str1_2={ 'HS_2' , 'HS_1'};

end

if rekuperace_TC == 2
            
        S1=[s1 , s2 , s3 , s3rek , s4 , s5 , s6 , s7 , s7rek , s1];
        T1=[t1 , t2 , t3 , t3rek , t4 , t5 , t6 , t7 , t7rek , t1];
        str1_1={'1' , '2' , '3' , '3rek' , '4' , '5' , '6' , '7' , '7rek' , '1'};

        S1=[s1 , s2 , s3 , s3rek , s4 , s6 , s7 , s7rek , s1];
        T1=[t1 , t2 , t3 , t3rek , t4 , t6 , t7 , t7rek , t1];
        str1_1={'1' , '2' , '3' , '3rek' , '4' , '6' , '7' , '7rek' , '1'}; 

    S_HS=[s1 , s3];
    T_HS=[t_HS2, t_HS1];
    str1_2={ 'HS_2' , 'HS_1'};
 
        S_ST1=[s7 , s5 , s4];
        T_ST1=[t_ST1 , t_ST2, t_ST3];
        str1_3={ 'ST_1' , 'ST_2' , 'ST_3'};
 
end

    figure(10);
    fig_1 = plot(S1, T1 , 'o-r', S_HS, T_HS, 'o-b', S_ST1, T_ST1, 'o-g', s1x1, temp_s1, 'k', s1x0, temp_s1, 'k');
    fig_1(1).LineWidth = 1.1;
    fig_1(2).LineWidth = 1.1;
    fig_1(3).LineWidth = 1.1;
    fig_1(4).LineWidth = 0.6;
    fig_1(5).LineWidth = 0.6;
    text(S1, T1 , str1_1, 'FontSize', 14);
    text(S_HS, T_HS , str1_2, 'FontSize', 14);
    text(S_ST1, T_ST1 , str1_3, 'FontSize', 14)
    title(['Nabíjení CB', newline,  'pracovní látka: ', fluid_working]);
    legend('TČ', 'HS - zdroj tepla', 'ST - uložiště tepla', 'Location', 'eastoutside');
    xlabel('s [kj/kgK]');
    ylabel('t [°C]');
    grid ON
    set(gcf,'position', [250 250 700 400])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                   GRAFY VYBIJENI

if rekuperace_ORC == 1
    
    S2=[s8 , s9 , s10, s11 , s12 , s13 , s14 , s8];
    T2=[t8 , t9 , t10, t11 , t12 , t13 , t14 , t8];
    str2_1={'8' , '9' , '10' , '11' , '12' , '13' , '14' , '8'};

    S_CL=[s8 , s14, s13];
    T_CL=[t_CL1, t_CL2 , t_CL3];
    str2_2={ 'CL_1' , 'CL_2' , 'CL_3' };

    S_ST2=[s9 , s10 , s12];
    T_ST2=[t_ST1 , t_ST4, t_ST3];
    str2_3={ 'ST_1' , 'ST_4' , 'ST_3' };        
        
end



if rekuperace_ORC == 2

    
    S2=[s8 , s9 , s9rek , s10, s11 , s12 , s13 , s13rek , s14 , s8];
    T2=[t8 , t9 , t9rek , t10, t11 , t12 , t13 , t13rek , t14 , t8];
    str2_1={'8' , '9' , '9rek' , '10' , '11' , '12' , '13' , '13rek' , '14' , '8'};

    S_CL=[s8 , s14, s13rek];
    T_CL=[t_CL1, t_CL2 , t_CL3];
    str2_2={ 'CL_1' , 'CL_2' , 'CL_3' };

    S_ST2=[s9rek , s10 , s12];
    T_ST2=[t_ST1 , t_ST4, t_ST3];
    str2_3={ 'ST_1' , 'ST_4' , 'ST_3' };

end

figure(2);
fig_2 = plot(S2, T2 , 'o-r', S_CL, T_CL, 'o-b', S_ST2, T_ST2, 'o-g', s2x1, temp_s2, 'k', s2x0, temp_s2, 'k');
fig_2(1).LineWidth = 1.1;
fig_2(2).LineWidth = 1.1;
fig_2(3).LineWidth = 1.1;
fig_2(4).LineWidth = 0.6;
fig_2(5).LineWidth = 0.6;
legend('ORC', 'CL - chlazeni', 'ST - uložiště tepla', 'Location', 'eastoutside');
text(S2, T2 , str2_1, 'FontSize', 14);
text(S_CL, T_CL , str2_2, 'FontSize', 14);
text(S_ST2, T_ST2 , str2_3, 'FontSize', 14);
title(['Vybíjení CB', newline,  'pracovní látka: ', fluid_working_dis]);
xlabel('s [kj/kgK]');
ylabel('t [°C]');
grid ON
set(gcf,'position', [250 250 700 400])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
%%                   KONEC VYPISU GRAFŮ

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if vysledky == 1
%%                    VYPIS VYSLEDKU VYPOCTU

disp([newline,'   VSTUPNÍ PARAMETRY:  ', newline]);
disp(['zdroj tepla: ', fluid_zdroj]);
disp(['teplota zdroje tepla na vstupu do tep. čerp.: ', num2str(t_zdroj_in), ' °C']);
disp(['teplota zdroje tepla na vstupu do tep. čerp.: ', num2str(t_zdroj_out), ' °C']);
disp(['hmotnostní tok zdroje tepla na vstupu do tep. čerp.: ', num2str(m_HS), ' kg/s']);
disp(['pracovní látka nabíjení: ', fluid_working]);
disp(['minimální teplotní rozdíl kondenzátor - nabíjení.: ', num2str(PP_cond), ' °C']);
disp(['minimální teplotní rozdíl výparník - nabíjení.: ', num2str(PP_evap), ' °C']);
disp(['hodnota izoentropické komprese - nabíjení.: ', num2str(eta_comp)]);
disp(['pracovní látka vybíjení: ', fluid_working_dis]);
disp(['zdroj chladu: ', fluid_chlazeni]);
disp(['teplota zdroje chladu na vstupu do ORC.: ', num2str(t_chladivo_in), ' °C']);
disp(['teplota zdroje chladu na výstupu z ORC.: ', num2str(t_chladivo_out), ' °C']);
disp(['hodnota izoentropické expanze - vybíjení.: ', num2str(eta_expander)]);
disp(['hodnota izoentropické komprese - čerpadlo - vybíjení.: ', num2str(eta_pump)]);
disp(['minimální teplotní rozdíl kondenzátor - vybíjení.: ', num2str(PP_cond_discharg), ' °C']);
disp(['minimální teplotní rozdíl výparník - vybíjení.: ', num2str(PP_evap_discharg), ' °C']);
disp(['médium zásobníku tepla.: ', fluid_zasobnik]);
disp(['teplota vybitého zásobníku tepla.: ', num2str(t_dis), ' °C']);
disp(['teplota nabitého zásobníku tepla.: ', num2str(t_charg), ' °C']);


disp([newline, '   NABIJENI VYPOCET   ', newline])
disp(['teplo přivedené do cyklu, Q in = ', num2str(Q_in_char),' kW']);
disp(['teplo odvedené z cyklu je, Q out = ', num2str(Q_out_char),' kW']);
disp(['práce přivedená na kompresoru, W comp = ', num2str(W_comp_char),' kW']);
disp(['topný koeficient tepelného čerpadla, COP = ', num2str(COP)]);
disp(['hmotnostní tok tepleného čerpadla, m wf = ', num2str(m_wf_charg), ' kg/s']);
disp(['hmotnostní tok zdroje tepla, m HS = ', num2str(m_HS), ' kg/s']);
disp(['hmotnostní tok uložiště tepla, m ST = ', num2str(m_ST_charg), ' kg/s']);


disp([newline, '   VYBIJENI VYPOCET   ' newline])
disp(['teplo přivedené do cyklu, Q in = ', num2str(Q_in_dis), ' kW' ]);
disp(['teplo odvedené z cyklu, Q out = ', num2str(Q_out_dis), ' kW' ]);
disp(['práce přivedená na čerpadle, W pump = ', num2str(W_pump), ' kW']);
disp(['práce expandéru, W exp = ', num2str(W_exp), ' kW']);
disp(['čistá práce cyklu, W cycle = ', num2str(W_cycle), ' kW']);
disp(['hmotnostní tok uložiště, m ST =', num2str(m_ST_discharg), ' kg/s']);
disp(['hmotnostní tok pracovní látky, m wf =', num2str(m_wf_discharg), ' kg/s']);
disp(['hmotnostní tok chlazení, m ST =', num2str(m_CL), ' kg/s']);
disp(['účinnost ORC systému = ', num2str((eta_ORC)*100), ' %']);
disp(['účinnost zpětné konverze systému(round trip) = ', num2str((eta_rt)*100), ' %']);
disp(['účinnost zpětné konverze systému(round trip elektric) = ', num2str((eta_rt_el)*100), ' %', newline])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
%%                       KONEC VYPISU VYSLEDKU

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if tabulky == 1
%%                        TABULKA JEDNOTLIVYCH STAVU


if rekuperace_TC == 1

        CHARG=[1, p1 , t1 , h1 , s1 ;
               2, p2 , t2 , h2 , s2 ;
               3, p3 , t3 , h3,  s3 ;
               4, p4 , t4 , h4 , s4 ;
               5, p5 , t5 , h5 , s5 ;
               6, p6 , t6 , h6 , s6 ;
               7, p7 , t7 , h7 , s7 ];     
else     
        CHARG=[1, p1 , t1 , h1 , s1 ;
               2, p2 , t2 , h2 , s2 ;
               3, p3 , t3 , h3,  s3 ;
               3, p3rek, t3rek, h3rek, s3rek ;
               4, p4 , t4 , h4 , s4 ;
               5, p5 , t5 , h5 , s5 ;
               6, p6 , t6 , h6 , s6 ;
               7, p7 , t7 , h7 , s7 ;
               7, p7rek, t7rek, h7rek, s7rek];    
end

CHARGING = [" CHARGING   " , " P [MPa]      " , " t  [°C]         " , " h  [kJ/kg]    " , " s  [kJ/kgK]   " ; 
              CHARG];

HEATSOURCE = [" HEAT SOURCE" ," h [kJ/kg]   " , " t [°C]" ;
              " HS1        " , h_HS1, t_HS1 ;
              " HS2        " , h_HS2, t_HS2 ];
   
STORAGE_CHARGING = [" STORAGE CHARGING" , " h [kJ/kg]  " , " t [°C]" ;
                    " ST1             " , h_ST1 , t_ST1 ;
                    " ST2             " , h_ST2 , t_ST2 ;
                    " ST3             " , h_ST3 , t_ST3 ];          




if rekuperace_ORC == 1    
    DISCHARG =[8 , p8 , t8 , h8 , s8 ;
             9 , p9 , t9 , h9,  s9 ;
             10 ,p10 , t10 , h10 , s10 ;
             11 , p11 , t11 , h11 , s11 ;
             12 , p12 , t12 , h12 , s12 ;
             13 , p13 , t13 , h13 , s13 ;
             14 , p14 , t14 , h14 , s14 ];   
else
    DISCHARG =[ 8 , p8 , t8 , h8 , s8 ;
                9 , p9 , t9 , h9,  s9 ;
                91 , p9rek , t9rek , h9rek , s9rek ;
                10 , p10 , t10 , h10 , s10 ;
                11 , p11 , t11 , h11 , s11 ;
                12 , p12 , t12 , h12 , s12 ;
                13 , p13 , t13 , h13 , s13 ;
                131 , p13rek , t13rek , h13rek , s13rek ;            
                14 , p14 , t14 , h14 , s14 ];        
end

DISCHARGING=[" DISCHARGING " , " P [MPa]          " , " t  [°C]          " , " h  [kJ/kg]         " , " s  [kJ/kgK]     " ; 
                 DISCHARG];

STORAGE_DISCHARGING = [" STORAGE DISCHARGING " , " h [kJ/kg]   " , "  t [°C] " ;
                        " ST1   " , h_ST1 , t_ST1 ;
                        " ST2   " , h_ST4 , t_ST4 ;
                        " ST3   " , h_ST3 , t_ST3  ];

COOLING = [" COOLING " , " h [kJ/kg] " , " t [°C] " ;
                        " CL1     " , h_CL1 , t_CL1 ;
                        " CL2     " , h_CL2 , t_CL2 ;
                        " CL3     " , h_CL3 , t_CL3];


disp(CHARGING);
disp(STORAGE_CHARGING);
disp(HEATSOURCE);

disp(DISCHARGING);
disp(STORAGE_DISCHARGING);
disp(COOLING);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
%%                       KONEC TABULKA JEDNOTLIVYCH STAVU



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                        EXPANDER FCE

function[P_sh, eta_exp, mdot, t_7, EXPANDER]=expander_fce(t_vyparovani, t_su, t_kondenzace)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% parametry expanderu A

    kapalina = 'R1233ZDE';
    mdot_nom = 0.12; % g/s
    A_leak = 4.06/(10^6); % mm^2 to m^2
    A_su = 27.43/(10^6); % mm^2 to m^2
    AU_su_nom = 21.06; % W/K
    AU_ex_nom = 34.68; % W/K
    AU_amb = 6.4; % W/K
    BVR = 4.05;
    V_sw = 147.987*0.000001; % cm3 to m3
    tau_loss = 0.47; % N*m
    otacky = 2000/60; % rpm = ot/min to ot/s

    
    p_1 = stav(kapalina,'TQ','P',t_vyparovani,1);
    t_1 = t_su;
    p_6 = stav(kapalina,'TQ','P',t_kondenzace,1);

    h_1 = stav(kapalina,'PT','H',p_1,t_1);
    s_1 = stav(kapalina,'PT','S',p_1,t_1);
    d_1 = stav(kapalina,'PT','D',p_1,t_1);
    v_1 = 1/d_1;

    T_amb = 20; % °C

    %% isoentropicky vypocet

    s_7_is = s_1;
    h_7_is = stav(kapalina,'PS','H',p_6 ,s_7_is);

    chyba2 = 2;
    T_w_odhad = 90; % °C

    while abs(chyba2) > 0.01

        if chyba2 < 0
            T_w_odhad = T_w_odhad-abs(chyba2)/8; %0.1;
        else
            T_w_odhad = T_w_odhad+abs(chyba2)/8; %0.1;
        end

    %% vypocet

        chyba = 1;

        t_2_odhad = t_1-1;

        while abs(chyba) > 0.01

            if chyba > 0
                t_2_odhad = t_2_odhad-((abs(8*chyba))/2);   % 0.001   %*i;
            else
                t_2_odhad = t_2_odhad+((abs(8*chyba))/2);   % 0.001   %*i;
            end

            t_2_odhad_kontrola = t_2_odhad;
            s_2_odhad = s_1;
            p_2_odhad = stav(kapalina,'TS','P',t_2_odhad,s_2_odhad);
            h_2_odhad = stav(kapalina,'TS','H',t_2_odhad,s_2_odhad);
            d_2_odhad = stav(kapalina,'TS','D',t_2_odhad,s_2_odhad);

            Cp_2 = stav(kapalina,'TS','Cp',t_2_odhad,s_2_odhad)*1000;

            mdot = d_2_odhad*A_su*((2*((h_1-h_2_odhad)*1000))^(1/2)); % kg/s

            AU_su = AU_su_nom*((mdot/mdot_nom)^0.8); % W/K

            Q_su = (1-exp((-AU_su)/(mdot*Cp_2)))*mdot*Cp_2*(t_2_odhad-T_w_odhad); % W

            h_3 = ((h_2_odhad*1000)-(Q_su/mdot))/1000; % kj/kg
            p_3 = p_2_odhad;
            t_3 = stav(kapalina,'PH','T',p_3,h_3);
            s_3 = stav(kapalina,'PH','S',p_3,h_3);
            d_3 = stav(kapalina,'PH','D',p_3,h_3);
            v_3 = 1/d_3;

            Cp_3 = stav(kapalina,'PH','Cp',p_3,h_3);
            Cv_3 = stav(kapalina,'PH','Cv',p_3,h_3);

            zeta = Cp_3/Cv_3;

            p_crit_leak = p_3*((2/(zeta+1))^(zeta/(zeta-1)));

            if p_crit_leak > p_6
                p_leak = p_crit_leak;
            else
                p_leak = p_1;
            end

            s_leak = s_3;
            h_leak = stav(kapalina,'PS','H',p_leak,s_leak);
            d_leak = stav(kapalina,'PS','D',p_leak,s_leak);

            mdot_leak = d_leak*A_leak*((2*((h_3-h_leak)*1000))^(1/2));

            mdot_nova = d_3*(V_sw/BVR)*otacky+mdot_leak;

            chyba = mdot_nova-mdot;
            chybicka = chyba;
        end

        s_4 = s_3;
        v_4 = v_3*BVR;
        d_4 = 1/v_4;
        h_4 = stav(kapalina,'SD','H',s_4,d_4);
        p_4 = stav(kapalina,'SD','P',s_4,d_4);
        t_4 = stav(kapalina,'SD','T',s_4,d_4);

        w_exp = ((h_3-h_4)*1000); % J/kg
        w_exp1 = (v_4*((p_4-p_6)*1000000));
        w_exp2 = w_exp+w_exp1;

        h_5 = h_4-w_exp1/1000; % kJ/kg
        p_5 = p_6;
        s_5 = stav(kapalina,'PH','S',p_5,h_5);
        t_5 = stav(kapalina,'PH','T',p_5,h_5);

        h_6 = (((mdot-mdot_leak)*h_5)+(mdot_leak*h_3))/mdot;
        s_6 = stav(kapalina,'PH','S',p_6,h_6);
        t_6 = stav(kapalina,'PH','T',p_6,h_6);
        Cp_6 = stav(kapalina,'PH','Cp',p_6,h_6)*1000; % J/kg/K

        AU_ex = AU_ex_nom*((mdot/mdot_nom)^0.8); % W/K
        Q_ex = (1-exp((-AU_ex)/(mdot*Cp_6)))*mdot*Cp_6*(T_w_odhad-t_6); % W

        h_7 = h_6+(Q_ex/1000)/mdot; 
        p_7 = p_6;
        s_7 = stav(kapalina,'PH','S',p_7,h_7);
        t_7 = stav(kapalina,'PH','T',p_7,h_7);

        P_in = (mdot-mdot_leak)*(((h_3-h_4)*1000)+(v_4*((p_4-p_6)*1000000))); % W
        P_loss = ((2*pi*(otacky*60))/60)*tau_loss; % W
        P_sh = (P_in-P_loss)/1000; % kW

        eta_exp = P_sh/(mdot*(h_1-h_7_is));
        eta_is = (h_7-h_1)/(h_7_is-h_1);
        Q_amb = P_loss + Q_su - Q_ex;
        T_w = (Q_amb/AU_amb)+T_amb;
        chyba2 = T_w-T_w_odhad;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                      VYSLEDKY

    EXPANDER = [" EXPANDER   " , " P [MPa]      " , " t  [°C]         " , " h  [kJ/kg]    " , " s  [kJ/kgK]   " ;
                   1, p_1 , t_1, h_1, s_1;
                   2, p_2_odhad , t_2_odhad, h_2_odhad, s_2_odhad;
                   3, p_3 , t_3, h_3, s_3;
                   4, p_4 , t_4, h_4, s_4;
                   5, p_5 , t_5, h_5, s_5;
                   6, p_6 , t_6, h_6, s_6;
                   7, p_7 , t_7, h_7, s_7];

    T = [t_1, t_2_odhad, t_3, t_4, t_5, t_6, t_7];
    S = [s_1, s_2_odhad, s_3, s_4, s_5, s_6, s_7];

    a=0;
    for i=(20):(182)
        a=a+1;
        temp_s1(a) = i;
        s1x0(a) = stav(kapalina, 'TQ', 'S', i, 0);
        s1x1(a) = stav(kapalina, 'TQ', 'S', i, 1);
    end  
    
        bb1=0;
    for i3 = (t_vyparovani+0.1):(t_vyparovani+50)
        bb1 = bb1+1;
        temp_p1(bb1) = i3;
        s1_p1(bb1) = stav(kapalina, 'TP', 'S', i3, p_1);
    end
    
        bb2=0;
    for i4 = (t_kondenzace+0.1):(t_kondenzace+50)
        bb2 = bb2+1;
        temp_p2(bb2) = i4;
        s7_p7(bb2) = stav(kapalina, 'TP', 'S', i4, p_7);        
    end    
    
    t_4_kond = stav(kapalina, 'PQ', 'T', p_4, 1);
        bb3=0;
    for i5 = (t_4_kond+0.1):(t_4_kond+50)
        bb3 = bb3+1;
        temp_p3(bb3) = i5;
        s4_p4(bb3) = stav(kapalina, 'TP', 'S', i5, p_4);
    end              

    figure(67);
    graf_expanze = plot(S, T, 'o-r', s1x1, temp_s1, 'k', s1_p1, temp_p1, '--b', s7_p7, temp_p2, '--b', s4_p4, temp_p3, '--b');
    graf_expanze(1).LineWidth = 1;
    graf_expanze(2).LineWidth = 0.6;
    graf_expanze(3).LineWidth = 0.6;
    graf_expanze(4).LineWidth = 0.6;
    graf_expanze(5).LineWidth = 0.6;
    
    legend('expanze ', 'horní mezní křivka (x = 1)', 'izobara')
    
    if t_4 < t_5        
        ylim([t_4-5 t_1+5])
    else
        ylim([t_5-5 t_1+5])
    end
    
    xlim([s_3*0.95 s_7*1.05 ])
    str1_1 = {'1' , '2' , '3' , '4' , '5' , '6' , '7'};
    text(S, T , str1_1, 'FontSize', 14);
    title('Graf expanze')
    xlabel('s [kJ/kg K]')
    ylabel('T [°C]')

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                        KOMPRESOR FCE

function[t_6, eta_comp, mdot_suc, P_comp, KOMPRESOR] = kompresor_fce(t_kondenzace, t_suc, t_vyparovani);

    %% Parametry kompresoru

    kapalina = 'R1233ZDE';
    mdot_nom = 0.988; % kg/s
    a_tl_1 = 0.272;
    a_tl_2 = 180.3;
    A_leak = 2.9/1000000; % m^2
    AU_su_nom = 66.83/1000; % kWK^-1
    AU_ex_nom = 39.325/1000; % kWK^-1
    b_hl = 1.82/1000; % kW.K^-5/4
    BVR = 3.23;
    V_sw = 0.0005727; % m3
    otacky = 2900/60; % ot/s
    T_amb = 10; % °C
    kin_vis = 0.000046; %% m2/s
    hustota_olej = 992; %% kg/m3
    dyn_vis = kin_vis/hustota_olej;

    %%

    p_kondenzace = stav(kapalina,'TQ','P', t_kondenzace, 1);
    p_vyparovani = stav(kapalina,'TQ','P', t_vyparovani, 1);    

    t_1 = t_suc;
    p_1 = p_vyparovani;
    p_6 = p_kondenzace;

    d_1 = stav(kapalina,'TP','D',t_1,p_1);
    v_1= 1/d_1;
    h_1 = stav(kapalina,'TP','H',t_1,p_1);
    s_1 = stav(kapalina,'TP','S',t_1,p_1);

    s_6_is = s_1;
    h_6_is = stav(kapalina,'PS','H',p_6,s_6_is);

    %% druha iterace

    chyba_temp = 2;
    T_w_ocek = 80;

    while abs(chyba_temp) > 0.01

        if chyba_temp > 0
            T_w_ocek = T_w_ocek+(abs(chyba_temp)/8);
        else
            T_w_ocek = T_w_ocek-(abs(chyba_temp)/8);
        end

        t_3_ocek = t_1+0.4;
        chyba =  5;

        while abs(chyba) >0.01

            if chyba > 0
                t_3_ocek = t_3_ocek + (abs(chyba)/2);
            else
                t_3_ocek = t_3_ocek - (abs(chyba)/2);
            end

            p_3 = p_1;
            h_3 = stav(kapalina,'PT','H',p_3,t_3_ocek);
            s_3 = stav(kapalina,'PT','S',p_3,t_3_ocek);
            d_3 = stav(kapalina,'PT','D',p_3,t_3_ocek);
            v_3 = 1/d_3;

            mdot = (V_sw*otacky)/v_3; % kg/s

            v_4 = v_3/BVR;
            d_4 = 1/v_4;
            s_4 = s_3;
            p_4 = stav(kapalina,'SD','P',s_4,d_4);
            h_4 = stav(kapalina,'SD','H',s_4,d_4);
            t_4 = stav(kapalina,'SD','T',s_4,d_4);

            p_5 = p_6;

            w_in = (h_4-h_3)+(v_4*((p_5-p_4)*1000)); % kW

            h_5 = h_3+(w_in);
            t_5 = stav(kapalina,'PH','T',p_5,h_5);
            s_5 = stav(kapalina,'PH','S',p_5,h_5);
            d_5 = stav(kapalina,'PH','D',p_5,h_5);
            v_5 = 1/d_5;
            Cp_5 = stav(kapalina,'PH','Cp',p_5,h_5);

            AU_ex = (AU_ex_nom)*(((mdot)/mdot_nom)^0.8); % kW/K
            Qdot_ex_56 = (1-exp((-AU_ex)/(mdot*Cp_5)))*mdot*Cp_5*(t_5-T_w_ocek); % kW
            h_6 = h_5 - Qdot_ex_56/mdot;
            s_6 = stav(kapalina,'HP','S',h_6,p_6);
            t_6 = stav(kapalina,'HP','T',h_6,p_6);
            d_6 = stav(kapalina,'HP','D',h_6,p_6);
            v_6 = 1/d_6;
            Cp_6 = stav(kapalina,'PH','Cp',p_6,h_6);
            Cv_6 = stav(kapalina,'PH','Cv',p_6,h_6);
            zeta = Cp_6/Cv_6;

            p_crit_leak = p_6*((2/(zeta+1))^(zeta/(zeta-1)));

            if p_crit_leak > p_1
                p_leak = p_crit_leak;
            else
                p_leak = p_1;
            end    

            s_leak = s_6;
            d_leak = stav(kapalina,'PS','D',p_leak,s_leak);
            h_leak = stav(kapalina,'PS','H',p_leak,s_leak);
            mdot_leak = d_leak*A_leak*((2*((h_6-h_leak)*(10^3)))^(1/2)); % kg/s

            mdot_suc = mdot - mdot_leak; % kg/s

            p_2 = p_1;
            h_2 = ((h_1*mdot_suc)+(h_6*mdot_leak))/(mdot_leak+mdot_suc);
            s_2 = stav(kapalina,'PH','S',p_2,h_2);
            t_2 = stav(kapalina,'PH','T',p_2,h_2);
            d_2 = stav(kapalina,'PH','D',p_2,h_2);
            v_2 = 1/d_2;
            Cp_2 = stav(kapalina,'PH','Cp',p_2,h_2);

            AU_su=(AU_su_nom)*(((mdot)/( mdot_nom))^0.8); % kW/K
            Qdot_su_23 = (1-exp((-AU_su)/((mdot)*(Cp_2))))*mdot*(Cp_2)*(T_w_ocek-t_2); % kW    

            h_3_nova = ((Qdot_su_23)/(mdot))+h_2;
            t_3_nova = stav(kapalina,'PH','T',p_3,h_3_nova);

            chyba = t_3_nova-t_3_ocek;

        end

        P_in = (w_in)*mdot; %% kW
        P_loss_1 = P_in*a_tl_1; %% kW
        P_loss_2 = (a_tl_2*dyn_vis*V_sw*((pi*(otacky*60)/30)^2))/1000; %% kW

        Qdot_amb = (P_loss_1+P_loss_2-Qdot_su_23+Qdot_ex_56); % kW

        T_w_nov = ((Qdot_amb)/(b_hl))^(4/5)+T_amb;

        chyba_temp = T_w_nov-T_w_ocek;

    end

    P_sh = P_in+P_loss_1+P_loss_2; % kW
    eta_comp = (mdot*(h_6_is-h_1))/(P_sh); % -
    P_comp = P_sh; % kW
    
        aa1=0;
    for i2=(20):(140)
        aa1=aa1+1;
        temp_s1(aa1) = i2;
        s1x0(aa1) = stav(kapalina, 'TQ', 'S', i2, 0);
        s1x1(aa1) = stav(kapalina, 'TQ', 'S', i2, 1);
    end   
    
        bb1=0;
    for i3 = (t_vyparovani+0.1):(t_vyparovani+50)
        bb1 = bb1+1;
        temp_p1(bb1) = i3;
        s1_p1(bb1) = stav(kapalina, 'TP', 'S', i3, p_1);

    end
    
        bb2=0;
    for i4 = (t_kondenzace+0.1):(t_kondenzace+50)
        bb2 = bb2+1;
        temp_p2(bb2) = i4;
        s6_p6(bb2) = stav(kapalina, 'TP', 'S', i4, p_6);

    end
    
    t_4_kond = stav(kapalina, 'PQ', 'T', p_4, 1);
        bb3=0;
    for i5 = (t_4_kond+0.1):(t_4_kond+50)
        bb3 = bb3+1;
        temp_p3(bb3) = i5;
        s4_p4(bb3) = stav(kapalina, 'TP', 'S', i5, p_4);

    end
    
    
    KOMPRESOR = [" KOMPRESOR   " , " P [MPa]      " , " t  [°C]         " , " h  [kJ/kg]    " , " s  [kJ/kgK]   " ;
                   1, p_1 , t_1, h_1, s_1;
                   2, p_2 , t_2, h_2, s_2;
                   3, p_3 , t_3_ocek, h_3, s_3;
                   4, p_4 , t_4, h_4, s_4;
                   5, p_5 , t_5, h_5, s_5;
                   6, p_6 , t_6, h_6, s_6];

    S1=[s_1, s_2, s_3, s_4, s_5, s_6];
    T1=[t_1, t_2, t_3_ocek, t_4, t_5, t_6];
    
    opengl('software')
    opengl('save', 'software')

    figure(66);
    str1_1 = {'1' , '2' , '3' , '4' , '5' , '6'};
    graf_komprese = plot(S1, T1 , 'o-r', s1x1 , temp_s1, 'k', s1_p1, temp_p1, '--b', s6_p6, temp_p2, '--b', s4_p4, temp_p3, '--b');
    graf_komprese(1).LineWidth = 1.4;
    graf_komprese(2).LineWidth = 0.5;
    graf_komprese(3).LineWidth = 0.5;
    graf_komprese(4).LineWidth = 0.5;
    graf_komprese(5).LineWidth = 0.5;
    title('Graf komprese');
    xlabel(' s [kJ/Kg K]');

    legend('komprese ', 'horni mezni krivka (x = 1)' , 'izobara');
    ylabel(' T [°C]');
    if t_4 > t_5
        ylim([t_1-15 t_4+15])
    else
        ylim([t_1-15 t_5+15])
    end
    xlim([s_1*0.995 s_6*1.005])
    text(S1, T1 , str1_1, 'FontSize', 14);         

    figure(66);

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





function[vysledek]=stav(fluid, input, output, input1, input2)

global RP
global MOLSI
global iMass
global iFlag

r = RP.REFPROPdll(fluid,input,output,MOLSI,iMass,iFlag,input1,input2,{1.0});
o=double(r.Output);
vysledek=o(1);
end