clear all
clc


%  Toto je základní skript pro výpočet tepelně integrované Carnotovy baterie - model s konstatními isoentropickými účinnostmi
%  
%  Všechny citlivostní analýzi pro studeny 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 skrze python refpropem je nutne zadat nasledujici kod : [v,e] = pyversion; system(['"', e, '" -m pip install --user -U ctREFPROP'])
%  
%
%  pro propojeni matlabu s skrze python coolpropem je nutne zadat nasledujici kod : [v,e] = pyversion; system([e,' -m pip install --user -U CoolProp'])
%


















%%% 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, jinak zadat 0
vysledky = 1; %% pokud 1 vytvori vypocet text s vysledky, jinak zadat 0
tabulky = 1; %% pokud 1 vytvori vypocet tabulku s jednotlivymi stavt, jinak zadat 0

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





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                       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_working = 'R1233ZDE'; %% pracovní látka TČ
    fluid_chlazeni = 'AIR';   %% zdroj chladu TČ
    p_chlazeni = 0.2; % MPa    tlak zdroje chladu
    t_chladivo_in = 20; %% °C  teplota zdroje chaldu na vstupu
    t_chladivo_out = 30; %% C  teplota zdroje chladu na výstupu
    PP_cond = 5; %% K          PP kondenzátor TČ
    PP_evap = 5; %% K          PP výparník TČ
    eta_comp = 0.7;  %         iso s účinnost komprese

    %% parametry vybijení

    fluid_working_dis = 'R1233ZDE'; %  zdroj tepla - charakter
    fluid_zdroj = 'WATER'; %% C    zdroj tepla
    t_zdroj_in = 70; %% C          teplota zdroje tepla na vstupu
    t_zdroj_out = 60; %% C         teplota zdroje tepla na výstupu
    p_zdroj = 0.1; %% MPa          tlak zdroje tepla
    m_HS = 2;                 %%   hmot tok zdroje tepla v kg/s
    eta_expander = 0.7;       %%   účinnost expandéru
    eta_pump = 0.7;          %%    úč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::MEG-60%';    %% látka uložiště  methylen glykol
    t_dis = -5; % C;                       %% teplota vybitého uložiště
    t_charg = -15; % C;                    %% teplota nabitého uložiště
    p_zasobnik = 0.1;                      %% tlak uložiště chladu
    m_STO_charg = 1;                       %% hmot tok zdroje uložiště do nabijeni
 
%%%  odsdud už není nutné nic přepisovat jen spouštět skript
     
    
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% NABÍJENÍ VÝPOČET
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

if rekuperace_TC == 1

    t_STO1 = t_charg;
    t_STO2 = t_dis;
    h_STO1 = py.CoolProp.CoolProp.PropsSI('H','P',p_zasobnik*(10^6),'T',t_STO1+273.15,fluid_zasobnik)/1000;        
    h_STO2 = py.CoolProp.CoolProp.PropsSI('H','P',p_zasobnik*(10^6),'T',t_STO2+273.15,fluid_zasobnik)/1000;   

    h_COO1 = stav(fluid_chlazeni,'TP','H',t_chladivo_in,p_chlazeni);
    h_COO3 = stav(fluid_chlazeni,'TP','H',t_chladivo_out,p_chlazeni);
    t_COO1 = t_chladivo_in;
    t_COO3 = t_chladivo_out;


    t1 = t_STO1-PP_evap;
    t2 = t1;
    t3 = t_STO2-PP_evap;
    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_COO1+PP_cond;
    t_condcharg=t_COO3+PP_cond;

    QQ = kontrola_kommprese(s3, h3, eta_comp, fluid_working, t_condcharg);

    if QQ < 1 & QQ > 0

        t6 = t_COO3 + PP_cond;
        t5 = t6;
        p7 = stav(fluid_working,'TQ','P',t6,1);
        p4 = p7;
        p6 = p7;
        p5 = p7;
        h7 = stav(fluid_working,'TP','H',t7,p7);
        h1=h7;
        m_wf_charg=(m_STO_charg*(h_STO2-h_STO1))/(h3-h1);
        s1 = stav(fluid_working,'PH','S',p1,h1);
        h6 = stav(fluid_working,'TQ','H',t6,0);
        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);
        q4 = stav(fluid_working,'PH','QMASS',p4,h4)
        s4 = stav(fluid_working,'PH','S',p4,h4);
        m_CL = (m_wf_charg*(h4-h7))/(h_COO3-h_COO1);

    end

    if QQ > 1 | QQ < 0  

        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=t6;
            p7= stav(fluid_working,'TQ','P',t5,1);
            p4=p7;
            p6=p7;
            p5=p7;
            h7 = stav(fluid_working,'TP','H',t7,p7);
            h1=h7;
            m_wf_charg=(m_STO_charg*(h_STO2-h_STO1))/(h3-h1);
            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);
            q4 = stav(fluid_working,'TP','Q',t4,p4);
            s4 = stav(fluid_working,'PH','S',p4,h4);
            m_CL = (m_wf_charg*(h4-h7))/(h_COO3-h_COO1);
            h_COO2 = (m_wf_charg*(h5-h7)+(h_COO1*m_CL))/m_CL;
            t_COO2 = stav(fluid_chlazeni,'PH','T',p_chlazeni,h_COO2);
            chyba = t5-(t_COO2+PP_cond);
            t_prehrati = t4-t5;

        end
    end

  %% konec iterace

    Q_in_char = m_STO_charg*(h_STO2-h_STO1);
    Q_out_char = m_CL*(h_COO3-h_COO1);
    W_comp_char = m_wf_charg*(h4-h3);
    COP=Q_out_char/W_comp_char;
    m_wf_charg;
    m_CL;
    m_STO_charg;

end

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






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                ZACATEK NABIJENI S REKUPERACI

if rekuperace_TC == 2

    t_STO1 = t_charg;
    t_STO2 = t_dis;
    h_STO1 = py.CoolProp.CoolProp.PropsSI('H','P',p_zasobnik*(10^6),'T',t_STO1+273.15,fluid_zasobnik)/1000;        
    h_STO2 = py.CoolProp.CoolProp.PropsSI('H','P',p_zasobnik*(10^6),'T',t_STO2+273.15,fluid_zasobnik)/1000;
    

    h_COO1 = stav(fluid_chlazeni,'TP','H',t_chladivo_in,p_chlazeni);
    h_COO3 = stav(fluid_chlazeni,'TP','H',t_chladivo_out,p_chlazeni);
    t_COO1 = t_chladivo_in;
    t_COO3 = t_chladivo_out;
    
    
    t1 = t_STO1-PP_evap;
    t2 = t1;
    t3 = t_STO2-PP_evap;
    t7=t_COO1+PP_cond;
    t3rek = t7 - PP_rek;
    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);
    h3rek = stav(fluid_working,'TP','H',t3rek,p3rek);
    s3rek = stav(fluid_working,'TP','S',t3rek,p3rek);
    s2 = stav(fluid_working,'TQ','S',t2,1);
    s3 = stav(fluid_working,'TP','S',t3,p3);   

    t_condcharg=t_COO3+PP_cond;

    Q_rek_nab = h3rek-h3;

    QQ = kontrola_kommprese(s3, h3, eta_comp, fluid_working, t_condcharg);

%%
    if QQ < 1 & QQ > 0

        t6 = t_COO3 + PP_cond;
        t5 = t6;
        p7 = stav(fluid_working,'TQ','P',t6,1);
        p7rek = p7;
        p4 = p7;
        p6 = p7;
        p5 = p7;
        h7 = stav(fluid_working,'TP','H',t7,p7);
        h7rek = h7-Q_rek_nab;
        s7rek = stav(fluid_working,'HP','S',h7rek,p7rek);
        t7rek = stav(fluid_working,'HP','T',h7rek,p7rek);
        h1=h7rek;
        m_wf_charg=(m_STO_charg*(h_STO2-h_STO1))/(h3-h1);
        s1 = stav(fluid_working,'PH','S',p1,h1);
        h6 = stav(fluid_working,'TQ','H',t6,0);
        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-h3rek)/eta_comp)+h3rek;
        t4 = stav(fluid_working,'PH','T',p4,h4);
        q4 = stav(fluid_working,'PH','QMASS',p4,h4);
        s4 = stav(fluid_working,'PH','S',p4,h4);
        m_CL = (m_wf_charg*(h4-h7))/(h_COO3-h_COO1);

    end

    if QQ > 1 | QQ < 0  
  
        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=t6;
            p7= stav(fluid_working,'TQ','P',t5,1);
            p7rek = p7;
            p4=p7;
            p6=p7;
            p5=p7;
            h7 = stav(fluid_working,'TP','H',t7,p7);
            h7rek = h7 - Q_rek_nab;
            s7rek = stav(fluid_working,'HP','S',h7rek,p7rek);
            t7rek = stav(fluid_working,'HP','T',h7rek,p7rek);
            h1=h7rek;
            m_wf_charg=(m_STO_charg*(h_STO2-h_STO1))/(h3-h1);
            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);
            q4 = stav(fluid_working,'TP','Q',t4,p4);
            s4 = stav(fluid_working,'PH','S',p4,h4);
            m_CL = (m_wf_charg*(h4-h7))/(h_COO3-h_COO1);
            h_COO2 = (m_wf_charg*(h5-h7)+(h_COO1*m_CL))/m_CL;
            t_COO2 = stav(fluid_chlazeni,'PH','T',p_chlazeni,h_COO2);
            chyba = t5-(t_COO2+PP_cond);
            t_prehrati = t4-t5;
  
        end
    end
    Q_in_char = m_STO_charg*(h_STO2-h_STO1);
    Q_out_char = m_CL*(h_COO3-h_COO1);
    W_comp_char = m_wf_charg*(h4-h3rek);
    COP=Q_out_char/W_comp_char;
    m_wf_charg;
    m_CL;
    m_STO_charg;

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                KONEC NABIJENI S REKUPERACI
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%













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












%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                      ZACATEK VYBIJENI BEZ REKUPERACE
if rekuperace_ORC == 1

    t_HHS1 = t_zdroj_in;
    t_HHS2 = t_zdroj_out;
    h_HHS1 = stav(fluid_zdroj,'TP','H',t_HHS1,p_zdroj);    
    h_HHS2 = stav(fluid_zdroj,'TP','H',t_HHS2,p_zdroj);

    t14 = (t_STO2-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_HHS2-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,'QP','H',1,p11);
            s11 = stav(fluid_working_dis,'QP','S',1,p11);
            h10 = stav(fluid_working_dis,'QP','H',0,p10);
            s10 = stav(fluid_working_dis,'QP','S',0,p10);
            t12=t_HHS1-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);
            m_wf_discharg = (m_HS*(h_HHS1-h_HHS2))/(h12-h9);
            h_HHS3 = ((m_wf_discharg*(h10-h9))+(h_HHS2*m_HS))/m_HS;
            t_HHS3 = stav(fluid_zdroj,'PH','T',p_zdroj,h_HHS3);
            chyba2=t_HHS3-(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);
        s13 = stav(fluid_working_dis,'PH','S',p13,h13);
        m_STO_discharg = ((h13-h8)*m_wf_discharg)/(h_STO2-h_STO1);
        h_STO3 = ((m_wf_discharg*(h14-h8))+(h_STO1*m_STO_discharg))/m_STO_discharg;
                        
        t_STO3 = py.CoolProp.CoolProp.PropsSI('T','P',p_zasobnik*(10^6),'H',h_STO3*(10^3),fluid_zasobnik)-273.15;
        chyba3=(t14-PP_cond_discharg)-t_STO3;
        
    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_rt=eta_ORC*COP;
    m_STO_discharg;
    m_wf_discharg;
    m_HS;
    
end

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                ZACATEK VYBIJENI S REKUPERACI

if rekuperace_ORC == 2

    t_HHS1 = t_zdroj_in;
    t_HHS2 = t_zdroj_out;
    h_HHS1 = stav(fluid_zdroj,'TP','H',t_HHS1,p_zdroj);
    h_HHS2 = stav(fluid_zdroj,'TP','H',t_HHS2,p_zdroj);

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

    while abs(chyba3) > 0.1
        if chyba3 > 0
            t14=t14-(abs(chyba3)/2);
        else 
            t14=t14+(abs(chyba3)/2);
        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_HHS2;

        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);
            p12 = p10;
            p9 = p10;
            p9rek = p10;
            p11 = p10;
            t12=t_HHS1-PP_evap_discharg;
            h12 = stav(fluid_working_dis,'TP','H',t12,p12);
            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);
            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);
            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);
            h13 = h12-((h12-h13is)*eta_expander);
            h9rek = h9+(h13 - h13rek);
            t9rek = stav(fluid_working_dis,'PH','T',p9rek,h9rek);
            m_wf_discharg = (m_HS*(h_HHS1-h_HHS2))/(h12-h9rek);
            h_HHS3 = ((m_wf_discharg*(h10-h9rek))+(h_HHS2*m_HS))/m_HS;
            t_HHS3 = stav(fluid_zdroj,'PH','T',p_zdroj,h_HHS3);
            chyba2=t_HHS3-(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_STO_discharg = ((h13rek-h8)*m_wf_discharg)/(h_STO2-h_STO1);
        h_STO3 = ((m_wf_discharg*(h14-h8))+(h_STO1*m_STO_discharg))/m_STO_discharg;
        t_STO3 = py.CoolProp.CoolProp.PropsSI('T','P',p_zasobnik*(10^6),'H',h_STO3*(10^3),fluid_zasobnik)-273.15;
        chyba3=(t14-PP_cond_discharg)-t_STO3;
        
    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_rt=eta_ORC*COP;
    m_STO_discharg;
    m_wf_discharg;
    m_HS;
    
end

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





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

%%%%%%%%%%%%%%%%%%%%%  vytvoření křivek sytosti pro nabíjení

    a=0;
    for i=(t1-5):(t5+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-5):(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

    %%%%%%%%%%%%%%%%%%%%%% vytvoreni grafu pro nabijeni

    if rekuperace_TC == 1    

        if QQ > 1 | QQ < 0 
            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_COO1=[s7 , s5 , s4];
            T_COO1=[t_COO1 , t_COO2, t_COO3];
            str1_3={ 'CL_1' , 'CL_2' , 'CL_3'};        
        end

        if QQ < 1 & QQ > 0 
            S1=[s1 , s2 , s3 , s4 , s6 , s7 , s1];
            T1=[t1 , t2 , t3 , t4 , t6 , t7 , t1];
            str1_1={'1' , '2' , '3' , '4' , '6' , '7' , '1'};

            S_COO1=[s7 , s4];
            T_COO1=[t_COO1 , t_COO3];
            str1_3={ 'CL_1' , 'CL_3'};        
        end
    end

    if rekuperace_TC == 2

        if QQ > 1 | QQ < 0 
            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'};

            S_COO1=[s7 , s5 , s4];
            T_COO1=[t_COO1 , t_COO2, t_COO3];
            str1_3={ 'CL_1' , 'CL_2' , 'CL_3'};
        end    

        if QQ < 1 & QQ > 0 
            S1=[s1 , s2 , s3 , s4 , s6 , s7 , s1];
            T1=[t1 , t2 , t3 , t4 , t6 , t7 , t1];
            str1_1={'1' , '2' , '3' , '3rek' , '4' , '6' , '7' , '7rek' , '1'};

            S_COO1=[s7 , s4];
            T_COO1=[t_COO1 , t_COO3];
            str1_3={ 'CL_1' , 'CL_3'};        
        end
    end

    S_STO1=[s1 , s3];
    T_STO1=[t_STO1, t_STO2];
    str1_2={ 'ST_1' , 'ST_2'};


    figure(1)
    fig_1 = plot(S1, T1 , 'o-r', S_STO1, T_STO1, 'o-b', S_COO1, T_COO1, '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;
    legend('TČ', 'ST - uložiště chladu', 'CL - chlazení', 'Location', 'eastoutside');
    text(S1, T1 , str1_1, 'FontSize', 14);
    text(S_STO1, T_STO1 , str1_2, 'FontSize', 14);
    text(S_COO1, T_COO1 , str1_3, 'FontSize', 14)
    title(['Nabíjení CB', newline,  'pracovní látka: ', fluid_working]);
    xlabel('s [kj/kgK]');
    ylabel('t [°C]');
    grid ON
    set(gcf,'position', [250 250 700 400])

    %%%%%%%%%%%%%%% vytvoreni grafu pro 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_STO2=[s8 , s14, s13];
        T_STO2=[t_STO1, t_STO3 , t_STO2];
        str2_2={ 'ST_1' , 'ST_3' , 'ST_2' };

        S_HHS=[s9 , s10 , s12];
        T_HHS=[t_HHS2 , t_HHS3 , t_HHS1];
        str2_3={ 'HS_2' , 'HS_3' , 'HS_1' };

    else

        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_STO2=[s8 , s14, s13rek];
        T_STO2=[t_STO1, t_STO3 , t_STO2];
        str2_2={ 'ST_1' , 'ST_3' , 'ST_2' };

        S_HHS=[s9rek , s10 , s12];
        T_HHS=[t_HHS2 , t_HHS3 , t_HHS1];
        str2_3={ 'HS_2' , 'HS_3' , 'HS_1' };



    end
    
    figure(2)
    fig_2 = plot(S2, T2 , 'o-r', S_HHS, T_HHS, 'o-b', S_STO2, T_STO2, '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;
    text(S2, T2 , str2_1, 'FontSize', 14);
    text(S_STO2, T_STO2 , str2_2, 'FontSize', 14);
    text(S_HHS, T_HHS , str2_3, 'FontSize', 14);
    title(['Vybíjení CB', newline,  'pracovní látka: ', fluid_working_dis]);
    legend('ORC', 'HS - zdroj tepla', 'ST - uložiště chladu', 'Location', 'eastoutside');
    xlabel('s [kj/kgK]');
    ylabel('t [°C]');
    grid ON
    set(gcf,'position', [250 250 700 400])
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
%%                     KONEC VYPISOVANI GRAFU




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if vysledky == 1
%%                      VYPISOVANI VYSLEDKU    

    disp([newline, '     VSTUPNÍ PARAMETRY:     ', newline]);
    disp(['  NABÍJENÍ - PARAMETRY  ']);
    disp(['pracovní látka nabíjení: ', fluid_working]);
    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(['zdroj chladu tep. čerpadlo: ', fluid_chlazeni]);
    disp(['teplota zdroje chladu na vstupu do tep. čerpadla.: ', num2str(t_chladivo_in), ' °C']);
    disp(['teplota zdroje chladu na výstupu z tep. čerpadla.: ', num2str(t_chladivo_out), ' °C']);
    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(['  VYBÍJENÍ - PARAMETRY  ']);
    disp(['pracovní látka vybíjení: ', fluid_working_dis]);
    disp(['zdroj tepla: ', fluid_zdroj]);
    disp(['teplota zdroje tepla na vstupu do ORC.: ', num2str(t_zdroj_in), ' °C']);
    disp(['teplota zdroje tepla na výstupu z ORC.: ', num2str(t_zdroj_out), ' °C']);
    disp(['hmotnostní tok zdroje tepla na vstupu do ORC.: ', num2str(m_HS), ' kg/s']);
    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([newline, '          NABÍJENÍ:        ', 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 chladiva, m CL = ', num2str(m_CL), ' kg/s']);
    disp(['hmotnostní tok uložiště tepla(chladu), m ST = ', num2str(m_STO_charg), ' kg/s']);

    disp([newline, '          VYBÍJENÍ:        ', 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_STO_discharg), ' kg/s']);
    disp(['hmotnostní tok pracovní látky, m wf =', num2str(m_wf_discharg), ' kg/s']);
    disp(['hmotnostní tok zdroje tepla, m HS =', num2str(m_HS), ' kg/s']);
    disp(['účinnost ORC systému = ', num2str((eta_ORC)*100), ' %']);
    disp(['účinnost zpětné konverze systému(round trip) = ', num2str((eta_rt)*100), ' %']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
%%                     KONEC GENEROVANI VYSLEDKU


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if grafy
%%                     GENEREOVANI TABULEK SE STAVY

    if rekuperace_TC == 1

        if QQ > 1 | QQ < 0    
            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 ];

            COOLING = ["COOLING" , " h [J/kg]   " , "   t [°C]" ;
                       "CL1    " , h_COO1 , t_COO1 ;
                       "CL2    " , h_COO2 , t_COO2 ;
                       "CL3    " , h_COO3 , t_COO3]; 
        else    
            CHARG=[1, p1 , t1 , h1 , s1 ;
                   2, p2 , t2 , h2 , s2 ;
                   3, p3 , t3 , h3,  s3 ;
                   4, p4 , t4 , h4 , s4 ;
                   6, p6 , t6 , h6 , s6 ;
                   7, p7 , t7 , h7 , s7 ];

            COOLING = ["COOLING" , " h [J/kg]   " , "   t [°C]" ;
                        "CL1    " , h_COO1 , t_COO1 ;
                        "CL3    " , h_COO3 , t_COO3];            
        end
    else
        if QQ > 1 | QQ < 0  

            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 ];

            COOLING = ["COOLING" , " h [J/kg]   " , "   t [°C]" ;
               "CL1    " , h_COO1 , t_COO1 ;
               "CL2    " , h_COO2 , t_COO2 ;
               "CL3    " , h_COO3 , t_COO3]; 
        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 ;
                   6, p6 , t6 , h6 , s6 ;
                   7, p7 , t7 , h7 , s7 ;
                   7, p7rek , t7rek , h7rek , s7rek ];

            COOLING = ["COOLING" , " h [J/kg]   " , "   t [°C]" ;
                       "CL1    " , h_COO1 , t_COO1 ;
                       "CL3    " , h_COO3 , t_COO3]; 
        end
    end

    STORAGE_CHARGING = ["STORAGE" , " h [J/kg]   " , "    t [°C]" ;
                        "ST1    " , h_STO1 , t_STO1 ;
                        "ST2    " , h_STO2 , t_STO2 ];    
    
    CHARGING = ["CHARGING" , "P [Pa]          " , "t  [°C]             " , "h  [J/kg]         " , "s  [J/kgK]     " ; 
                CHARG];


    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 ];


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

        STORAGE_DISCHARGING = ["STORAGE" , " h [J/kg]   " , "   t [°C]" ;
                                 "ST1    " , h_STO1 , t_STO1 ;
                                 "ST2    " , h_STO3 , t_STO3 ;
                                 "ST3    " , h_STO2 , t_STO2]; 


        HEAT_SOURCE = ["HEAT SOURCE" , " h [J/kg]   " , "   t [°C]" ;
                         "HS1        " , h_HHS1 , t_HHS1 ;
                         "HS3        " , h_HHS3 , t_HHS3 ;
                         "HS2        " , h_HHS2 , t_HHS2 ];
        
    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 ];               

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

        STORAGE_DISCHARGING = ["STORAGE" , " h [J/kg]   " , "   t [°C]" ;
                               "ST1    " , h_STO1 , t_STO1 ;
                               "ST2    " , h_STO3 , t_STO3 ;
                               "ST3    " , h_STO2 , t_STO2  ]; 


        HEAT_SOURCE = ["HEAT SOURCE" , " h [J/kg]   " , "   t [°C] " ;
                         "HS1        " , h_HHS1 , t_HHS1 ;
                         "HS3        " , h_HHS3 , t_HHS3 ;
                         "HS2        " , h_HHS2 , t_HHS2  ];
    end

disp(CHARGING);
disp(STORAGE_CHARGING);
disp(COOLING);
              
disp(DISCHARGING);
disp(STORAGE_DISCHARGING);
disp(HEAT_SOURCE);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
%%                      KONEC GENEROVANI TABULEK 


if rekuperace_TC == 1

    if QQ < 0 | QQ > 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 ;
               4, p4 , t4 , h4 , s4 ;
               6, p6 , t6 , h6 , s6 ;
               7, p7 , t7 , h7 , s7 ];
    end     
else    
    if QQ < 0 | QQ > 1 
        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];    
    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 ;
               6, p6 , t6 , h6 , s6 ;
               7, p7 , t7 , h7 , s7 ;
               7, p7rek, t7rek, h7rek, s7rek ];    
    end
end

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





function[q4]=kontrola_kommprese(s3, h3, eta_comp, fluid_working, t_condcharg)

t5=t_condcharg;
p5=stav(fluid_working,'TQ','P',t5,1);
s4is = s3;
h4is = stav(fluid_working,'PS','H',p5,s4is);
h4 = ((h4is-h3)/eta_comp)+h3;
q4 = stav(fluid_working,'PH','QMASS',p5, h4);

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


