% Reactor for catalytic oxidation of selected VOC 

% Oxidation reactions according to general formula:
% CxHy + (x+y/4) O2 --> x CO2 + (y/2) H2O
% nejdriv vyplnit zadani "CTRL+H" --> "ZADA SE", pak pocitat


clear all; clc; close all; format long;

% Predexponencialny faktor A: [mol/s/Pa/kg_cat]
A_toluen = 14.1;
A_aceton = 28.85; 
A_etanol = 0.45; 
A_metan = 0.21; 

% Aktivacni energie Ea: [J/mol]
Ea_toluen = 75500; 
Ea_aceton = 96700; 
Ea_etanol = 65700; 
Ea_metan = 86400; 

% Stechiometrie oxidacnich reakci

% Metan - CH4 + 2O2 => CO2 + 2 H2O
stech_rmetan_CH4 = 1; stech_rmetan_O2 = 2; stech_rmetan_CO2 = 1; stech_rmetan_H2O = 2;
% Toluen - C7H8 + 9O2 => 7 CO2 + 4 H2O
stech_rtoluen_C7H8 = 1; stech_rtoluen_O2 = 9; stech_rtoluen_CO2 = 7; stech_rtoluen_H2O = 4;
% Etanol - C2H6O + 3O2 => 2 CO2 + 3 H2O
stech_retanol_C2H6O = 1; stech_retanol_O2 = 3; stech_retanol_CO2 = 2; stech_retanol_H2O = 3;
% Aceton - C3H6O + 4O2 => 3 CO2 + 3 H2O
stech_raceton_C3H6O = 1; stech_raceton_O2 = 4; stech_raceton_CO2 = 3; stech_raceton_H2O = 3;

% Molekulove hmotnosti [kg/kmol]
Mw_metan = 16.042;
Mw_toluen = 92.138;
Mw_etanol = 46.068;
Mw_aceton = 58.079;
Mw_N2 = 28.0134;
% Mw_Ar = 39.948;
Mw_CO2 = 44.00995;
Mw_CO = 28.01054;
Mw_O2 = 31.99879;
Mw_H2O = 18.015;

% Pocatecne slozeni [-] - ZADA SE !!!
y0_metan = 0;
y0_toluen = 0.0019;
y0_etanol = 0;
y0_aceton = 0;
% y0_Ar =;
y0_CO2 = 0.0474;
y0_CO = 0;
y0_O2 = 0.136;
y0_H2O = 0.01;
% y0 = 1 =
% y0_metan+y0_toluen+y0_etanol+y0_aceton+y0_N2+y0_Ar+y0_CO2+y0_CO+y0_O2+y0_H2O

% dopocet dusiku do 100 %
y0_ALLbezN2 = y0_metan+y0_toluen+y0_etanol+y0_aceton+y0_CO2+y0_CO+y0_O2+y0_H2O;
y0_N2 =1-y0_ALLbezN2;

% Vypocet pocatecni stredni molekulove hmotnosti smesi [kg/kmol]
Mw0_str = Mw_metan*y0_metan+Mw_toluen*y0_toluen+Mw_etanol*y0_etanol+Mw_aceton*y0_aceton+Mw_CO2*y0_CO2+Mw_CO*y0_CO+Mw_O2*y0_O2+Mw_H2O*y0_H2O+Mw_N2*y0_N2; % [kg/kmol]
%+Mw_Ar*y0_Ar

% Vlastnosti loze a katalyzatoru [-]
w_cat_max = 0.675 ; % [kg_cat] celkova hmotnost katalyzatoru - ZADA SE !!!
dr = 0.150; % [m] prumer reaktoru - ctvercovy prurez
lreactor = 1.5; % [m] vyska celeho reaktoru ZADA SE!!!
area_reactor = dr*dr; % [m2] plocha loze
rho_cat = 600; % [kg/m3] merna hmotnost katalyzatoru - ZDROJ datasheet
lr = w_cat_max/rho_cat/area_reactor; % [m] vyska loze katalyzatoru
porosity_cat = 0.44; % [-] odhad - muze byt pozdeji upresneno
dp_cat = 0.005; % [m] velikost castic katalyzatoru

% Pocatecni (aktualni) podminky
w_cat = 0; % pocatecni hmotnost katalyzatoru
V = 24.3; % [Nm^3/h] normalovy prutok plynu - ZADA SE !!!
Vm = 0.022414; % [Nm^3/mol]  normalni molarni objem
n0_SUM = V/Vm/3600; % [mol/s] pocatecni celkovy latkovy tok
m0_SUM = n0_SUM*Mw0_str/1000; % [kg/s] pocatecni celkovy hmostnostni tok
temp_K = 607.95; %[K] - teplota na zacatku loze - ZADA SE !!! 
temp0_K = temp_K; %[K] - teplota na zacatku loze
temp0_C = temp_K-273.15;
p0 = 99700; %[Pa] - ZADA SE !!!
p = p0; %[Pa] pocatecni tlak

% Pocatecni latkove toky [mol/s]
n0_metan = n0_SUM*y0_metan;  % [mol/s]
n0_toluen = n0_SUM*y0_toluen; % [mol/s]
n0_etanol = n0_SUM*y0_etanol; % [mol/s]
n0_aceton = n0_SUM*y0_aceton; % [mol/s]
n0_N2 = n0_SUM*y0_N2; % [mol/s]
% n0_Ar = n0_SUM*y0_Ar; % [mol/s]
n0_CO = n0_SUM*y0_CO; % [mol/s]
n0_CO2 = n0_SUM*y0_CO2; % [mol/s]
n0_O2 = n0_SUM*y0_O2; % [mol/s]
n0_H2O = n0_SUM*y0_H2O; % [mol/s]

% pocatecni entalpie a cp smesi plynu
[iCO,iO2,iH2O,iCO2,iN2,iC7H8,iC3H6O,iC2H7O,iCH4] = entalpie(temp_K);
i_SUM = iCO*y0_CO+iO2*y0_O2+iH2O*y0_H2O+iCO2*y0_CO2+iN2*y0_N2+iC7H8*y0_toluen+iC3H6O*y0_aceton+iC2H7O*y0_etanol+iCH4*y0_metan; % [J/mol]
cp_SUM = i_SUM/temp0_C; % [J.mol-1.K-1]

% Veliciny pro integraci
dw_cat = 0.00675; %[kg_cat] integracni krok (krok prirustku hmotnosti katalyzatoru)
vektor_m_cat = 0:dw_cat:w_cat_max;
R = 8.3143; %[J/mol/K] universalni plynova konst.

% Ulozeni pocatecnych promennych pred cyklem
np = 1 ; % indexovanie pro ukladanie hodnot v jednotlivých krocich 

y_metan_1(np) = y0_metan; % [-] - látkove zlomky 
y_toluen_1(np) = y0_toluen;
y_etanol_1(np) = y0_etanol;
y_aceton_1(np) = y0_aceton;
% y_Ar_1(np) = y0_Ar;
y_CO2_1(np) = y0_CO2;
y_CO_1(np) = y0_CO;
y_O2_1(np) = y0_O2;
y_H2O_1(np) = y0_H2O;

temp_K_1(np) = temp0_K; %[K] - teplota na zacatku loze
temp_C_1(np) = temp0_K-273.15; %[°C] - teplota na zacatku loze

p_1(np) = p0; %[Pa]

i_SUM_1(np) = i_SUM; % [J/mol]
cp_SUM_1(np) = cp_SUM; % [J.mol-1.K-1]

Mw_str_1(np) = Mw0_str; % [kg/kmol]

n_metan_1(np) = n0_metan;  % [mol/s] - latkove toky
n_toluen_1(np) = n0_toluen; % [mol/s]
n_etanol_1(np) = n0_etanol; % [mol/s]
n_aceton_1(np) = n0_aceton; % [mol/s]
n_N2_1(np) = n0_N2; % [mol/s]
%n_Ar_1(np) =  n0_Ar; % [mol/s]
n_CO_1(np) = n0_CO; % [mol/s]
n_CO2_1(np) = n0_CO2; % [mol/s]
n_O2_1(np) = n0_O2; % [mol/s]
n_H2O_1(np) = n0_H2O; % [mol/s]

[hustota,viskozita,lambda_spaliny] = hustota_viskozita(temp_K,p,y0_CO2,y0_O2,y0_H2O);
hustota_1(np) = hustota; % [kg/m^3]
viskozita_1(np) = viskozita; % [Pa.s]
lambda_spaliny_1(np) = lambda_spaliny; % [W.m-1.K-1]


% INTEGRACE v cyklu
m = length(vektor_m_cat); % pocet kroku integrace
for i = 2:m % (cyklus else zacina u dat ukladanych pod np=3)
    
    if i < 3; % vypocet pocatecnych hodnot 
        
    % rychlostni konstanty
    k_metan = A_metan*exp(-Ea_metan/R/temp_K);
    k_toluen = A_toluen*exp(-Ea_toluen/R/temp_K);
    k_etanol = A_etanol*exp(-Ea_etanol/R/temp_K);
    k_aceton = A_toluen*exp(-Ea_aceton/R/temp_K);
    
    % reakcni rychlosti % [mol.s-1.kg_cat-1]
    rate_metan = k_metan*y0_metan*p0;
    rate_toluen = k_toluen*y0_toluen*p0;
    rate_etanol = k_etanol*y0_etanol*p0;
    rate_aceton = k_aceton*y0_aceton*p0;
    
    % Ulozeni reakcnich rychlosti % [mol.s-1.kg_cat-1]
    rate_metan_1(np) = rate_metan;
    rate_toluen_1(np) = rate_toluen;
    rate_etanol_1(np) = rate_etanol;
    rate_aceton_1(np) = rate_aceton;

    % dfcdw... zmena rychlosti reakce v danem kroku [mol.s-1.kg_cat-1]
    dfcdw_metan = rate_metan*dw_cat;
    dfcdw_toluen = rate_toluen*dw_cat;
    dfcdw_etanol = rate_etanol*dw_cat;
    dfcdw_aceton = rate_aceton*dw_cat;
    
    % zmena latkových toku vlivem reakcii v danem kroku (neni zde reakce CO!) [mol.s-1.kg_cat-1]
    
    % VOC
    dn_metan = dfcdw_metan*stech_rmetan_CH4;
    dn_toluen = dfcdw_toluen*stech_rtoluen_C7H8;
    dn_etanol = dfcdw_etanol*stech_retanol_C2H6O;
    dn_aceton = dfcdw_aceton*stech_raceton_C3H6O;
    
    % O2
    dn_O2_metan = dfcdw_metan*stech_rmetan_O2;
    dn_O2_toluen = dfcdw_toluen*stech_rtoluen_O2;
    dn_O2_etanol = dfcdw_etanol*stech_retanol_O2;
	dn_O2_aceton = dfcdw_aceton*stech_raceton_O2;
    
    % CO2
    dn_CO2_metan = dfcdw_metan*stech_rmetan_CO2;
    dn_CO2_toluen = dfcdw_toluen*stech_rtoluen_CO2;
    dn_CO2_etanol = dfcdw_etanol*stech_retanol_CO2;
	dn_CO2_aceton = dfcdw_aceton*stech_raceton_CO2;
    
    % H2O
    dn_H2O_metan = dfcdw_metan*stech_rmetan_H2O;
    dn_H2O_toluen = dfcdw_toluen*stech_rtoluen_H2O;
    dn_H2O_etanol = dfcdw_etanol*stech_retanol_H2O;
	dn_H2O_aceton = dfcdw_aceton*stech_raceton_H2O;
    
    % Ulozeni dnVOC
    dn_metan_1(np) = dn_metan;
    dn_toluen_1(np) = dn_toluen;
    dn_etanol_1(np) = dn_etanol;
    dn_aceton_1(np) = dn_aceton;
    
    % prirustek teploty vlivem reakcniho tepla a pokles teploty vlivem tepelnych ztrat
    
    [cp_metan,cp_toluen,cp_etanol,cp_aceton,cp_CO2,cp_O2,cp_H2O,cp_N2,dHr_metan,dHr_toluen,dHr_etanol,dHr_aceton] = reakcne_teplo(temp_K); % vypocet tepelnych kapacit a reakcneho tepla 
    
    dCp = (cp_metan*n0_metan+cp_toluen*n0_toluen+cp_etanol*n0_etanol+cp_aceton*n0_aceton+cp_CO2*n0_CO2+cp_O2*n0_O2+cp_H2O*n0_H2O+cp_N2*n0_N2); % [J.K-1.s-1]
    
    dtempdw = (-dHr_metan*rate_metan-dHr_toluen*rate_toluen-dHr_etanol*rate_etanol-dHr_aceton*rate_aceton)/dCp; % [K.kg_cat-1]
    
        % vypocet parametru pro vypocet dQ a dp
        [hustota,viskozita,lambda_spaliny] = hustota_viskozita(temp_K,p,y0_CO2,y0_O2,y0_H2O);
        V_real =n0_SUM*R*temp_K/p; % [m3/s]
        velocity = V_real/area_reactor ; % [m/s]
        Reynolds = dp_cat*velocity*hustota/viskozita; % [-]
        Mw_str = Mw0_str;
        
        % vypocet tepelnych a teplotnych ztrat (dQ,dT)
        [dQ_ztraty,dT_ztraty] = tepelne_ztraty(hustota,viskozita,Reynolds,dw_cat,rho_cat,area_reactor,cp_SUM,dp_cat,dr,temp_K,lr,Mw_str,m0_SUM,lambda_spaliny,w_cat_max,dtempdw,lreactor);
    
        dQ_ztraty_SUM = 0; % [W] celkove kumulujici se tepelne ztraty
        dQ_ztraty_SUM = dQ_ztraty_SUM+dQ_ztraty; % [W] pricitaji se tepelne ztraty na danem kroku 
        dT_ztraty_SUM = 0; % [K] celkove kumulujici se poklesy teploty tepelnymi ztratami
        dT_ztraty_SUM = dT_ztraty_SUM+dT_ztraty; % [K] pricitaji se poklesy teploty tepelnymi ztratami na danem kroku
    
    temp_K = temp_K+dtempdw*dw_cat-dT_ztraty; % nova teplota zvysena vlivem reakcniho tepla a snizena o tepelne ztraty 
    
    % prirustek tlaku vlivem tlakove ztraty
    Re0 = Reynolds; % [-] - pocatecni Reynolds
    
        % Vypocet modifikovaneho Reynoldsova cisla pro vypocet dp
        G_area = m0_SUM/area_reactor; % [kg.s-1.m-2] - hmotnostni tok vztazen na prut. plochu
        Rem = dp_cat*G_area/viskozita; % [-] - modifikovane Reynoldsove cislo
        
        % Vypocet tlakove ztraty
        [N,fm] = tlakova_ztrata_Levy (Rem);
        theta = 1; % [-] - tvarovy faktor (pro koule = 1)
        dp_l = (fm*2*G_area^2*(1-porosity_cat)^(3-N))/(dp_cat*hustota*theta^(3-N)*porosity_cat^3); % [Pa/m]
        dpdw = dp_l*lr*dw_cat/w_cat_max; % [Pa] - tlakova ztrata v integracnom kroku
        p = p+dpdw; % [Pa] - system je v podtlaku takze +dpdw
        
    
    % prirastok katalyzatoru o krok dw
    w_cat = w_cat+dw_cat; % [kg]
    
    % nove latkove toky % [mol/s] - nezahrnuje reakci CO !!!
    n_metan = n0_metan-dn_metan; % [mol/s]
    n_toluen = n0_toluen-dn_toluen; % [mol/s]
    n_etanol = n0_etanol-dn_etanol; % [mol/s]
    n_aceton = n0_aceton-dn_aceton; % [mol/s]
    n_N2 = n0_N2; % [mol/s] 
    % n_Ar = n0_Ar; % [mol/s]
    n_CO = n0_CO; % [mol/s]
    n_CO2 = n0_CO2+dn_CO2_metan+dn_CO2_toluen+dn_CO2_etanol+dn_CO2_aceton; % [mol/s]
    n_O2 = n0_O2-dn_O2_metan-dn_O2_toluen-dn_O2_etanol-dn_O2_aceton; % [mol/s]
    n_H2O = n0_H2O+dn_H2O_metan+dn_H2O_toluen+dn_H2O_etanol+dn_H2O_aceton; % [mol/s]
    n_SUM = n_metan+n_toluen+n_etanol+n_aceton+n_N2+n_CO+n_CO2+n_O2+n_H2O; % [-]
    
    % nove slozeni smesi (molarni zlomky) % [-] - nezahrnuje reakci CO !!!, % nezahrnuje Ar  
    y_metan = n_metan/n_SUM; % [-]
    y_toluen = n_toluen/n_SUM; % [-]
    y_etanol = n_etanol/n_SUM; % [-]
    y_aceton = n_aceton/n_SUM; % [-]
    y_N2 = n_N2/n_SUM; % [-]
    % y_Ar = n_Ar/n_SUM; % [-];
    y_CO = n_CO/n_SUM; % [-];
    y_CO2 = n_CO2/n_SUM; % [-];
    y_O2 = n_O2/n_SUM; % [-];
    y_H2O = n_H2O/n_SUM; % [-];
    
    % Vypocet stredni molekulove hmotnosti smesi [kg/kmol]
    Mw_str = Mw_metan*y_metan+Mw_toluen*y_toluen+Mw_etanol*y_etanol+Mw_aceton*y_aceton+Mw_CO2*y_CO2+Mw_CO*y_CO+Mw_O2*y_O2+Mw_H2O*y_H2O+Mw_N2*y_N2; % [kg/kmol]
    %+Mw_Ar*y_Ar
    
    % nova entalpie a cp smesi
    [iCO,iO2,iH2O,iCO2,iN2,iC7H8,iC3H6O,iC2H7O,iCH4] = entalpie(temp_K);
    i_SUM = iCO*y_CO+iO2*y_O2+iH2O*y_H2O+iCO2*y_CO2+iN2*y_N2+iC7H8*y_toluen+iC3H6O*y_aceton+iC2H7O*y_etanol+iCH4*y_metan; % [kJ/kmol]
    cp_SUM = i_SUM/(temp_K-273.15); % [J.mol-1.K-1]

    
    % ulozeni dat v danem kroku
    np = np+1;
    temp_K_1(np) = temp0_K; %[K] - teplota 
    temp_C_1(np) = temp0_K-273.15; %[C] - teplota 
    dtempdw_1(np-1) = dtempdw; % [K.kg_cat-1]
    dpdw_1(np-1) = dpdw; % [Pa] - tlakova ztrata v integracnom kroku
    hustota_1(np) = hustota; % [kg/m^3]
    viskozita_1(np) = viskozita; % [Pa.s]
    lambda_spaliny_1(np) = lambda_spaliny; % [W.m-1.K-1]
    velocity_1(np) = velocity; % [m/s]
    Re_1(np) = Reynolds; % [-]
    Rem_1(np) = Rem; % [-]
    p_1(np) = p; % [Pa]
    i_SUM_1 (np) = i_SUM; % [J/mol]
    cp_SUM_1(np) = cp_SUM; % [J.mol-1.K-1]
    Mw_str_1(np) = Mw_str; % [kg/kmol]
    
    dQ_ztraty_1(np-1) = dQ_ztraty; % [W] tepelne ztraty v danem kroku
    dQ_ztraty_SUM_1(np-1) = dQ_ztraty_SUM; % [W] celkove kumulujici se tepelne ztraty
    dT_ztraty_1(np-1) = dT_ztraty; % [K] pokles teploty tepelnymi ztratami v danem kroku
    dT_ztraty_SUM_1(np-1) =  dT_ztraty_SUM; % [K] celkove kumulujici se poklesy teploty tepelnymi ztratami
    
    y_metan_1(np) = y_metan; % [-] (molarni zlomek) 
    y_toluen_1(np) = y_toluen; % [-] (molarni zlomek)
    y_etanol_1(np) = y_etanol; % [-] (molarni zlomek)
    y_aceton_1(np) = y_aceton; % [-] (molarni zlomek)
    y_N2_1(np)= y_N2; % [-] (molarni zlomek)
    % y_Ar_1(np) = y_Ar; % [-] (molarni zlomek)
    y_CO_1(np) = y_CO; % [-] (molarni zlomek)
    y_CO2_1(np) = y_CO2; % [-] (molarni zlomek)
    y_O2_1(np) = y_O2; % [-] (molarni zlomek)
    y_H2O_1(np) = y_H2O; % [-] (molarni zlomek)
    
    n_metan_1(np) = n_metan;  % [mol/s] - latkove toky
    n_toluen_1(np) = n_toluen; % [mol/s]
    n_etanol_1(np) = n_etanol; % [mol/s]
    n_aceton_1(np) = n_aceton; % [mol/s]
    n_N2_1(np) = n_N2; % [mol/s]
    %n_Ar_1(np) =  n_Ar; % [mol/s]
    n_CO_1(np) = n_CO; % [mol/s]
    n_CO2_1(np) = n_CO2; % [mol/s]
    n_O2_1(np) = n_O2; % [mol/s]
    n_H2O_1(np) = n_H2O; % [mol/s]
    
    x_metan_1(np) = (y0_metan-y_metan_1(np))/(y0_metan); % [-] stupen konverze
    x_toluen_1(np) = (y0_toluen-y_toluen_1(np))/(y0_toluen); % [-] stupen konverze
    x_etanol_1(np) = (y0_etanol-y_etanol_1(np))/(y0_etanol); % [-] stupen konverze
    x_aceton_1(np) = (y0_aceton-y_aceton_1(np))/(y0_aceton); % [-] stupen konverze
    
else    
    
    % rychlostni konstanty
    k_metan = A_metan*exp(-Ea_metan/R/temp_K);
    k_toluen = A_toluen*exp(-Ea_toluen/R/temp_K);
    k_etanol = A_etanol*exp(-Ea_etanol/R/temp_K);
    k_aceton = A_toluen*exp(-Ea_aceton/R/temp_K);
    
    % reakcni rychlosti % [mol.s-1.kg_cat-1]
    rate_metan = k_metan*y_metan*p;
    rate_toluen = k_toluen*y_toluen*p;
    rate_etanol = k_etanol*y_etanol*p;
    rate_aceton = k_aceton*y_aceton*p;
    
    % Ulozeni reakcnich rychlosti % [mol.s-1.kg_cat-1]
    rate_metan_1(np) = rate_metan;
    rate_toluen_1(np) = rate_toluen;
    rate_etanol_1(np) = rate_etanol;
    rate_aceton_1(np) = rate_aceton;
    
    % dfcdw... zmena rychlosti reakce v danem kroku [mol.s-1.kg_cat-1]
    dfcdw_metan = rate_metan*dw_cat;
    dfcdw_toluen = rate_toluen*dw_cat;
    dfcdw_etanol = rate_etanol*dw_cat;
    dfcdw_aceton = rate_aceton*dw_cat;
    
    % zmena latkových toku vlyvem reakcii v danem kroku (neni zde reakce CO!) [mol.s-1.kg_cat-1]
    
    % VOC
    dn_metan = dfcdw_metan*stech_rmetan_CH4;
    dn_toluen = dfcdw_toluen*stech_rtoluen_C7H8;
    dn_etanol = dfcdw_etanol*stech_retanol_C2H6O;
    dn_aceton = dfcdw_aceton*stech_raceton_C3H6O;
    
    % Ulozeni dnVOC
    dn_metan_1(np) = dn_metan;
    dn_toluen_1(np) = dn_toluen;
    dn_etanol_1(np) = dn_etanol;
    dn_aceton_1(np) = dn_aceton;
    
    % O2
    dn_O2_metan = dfcdw_metan*stech_rmetan_O2;
    dn_O2_toluen = dfcdw_toluen*stech_rtoluen_O2;
    dn_O2_etanol = dfcdw_etanol*stech_retanol_O2;
	dn_O2_aceton = dfcdw_aceton*stech_raceton_O2;
    
    % CO2
    dn_CO2_metan = dfcdw_metan*stech_rmetan_CO2;
    dn_CO2_toluen = dfcdw_toluen*stech_rtoluen_CO2;
    dn_CO2_etanol = dfcdw_etanol*stech_retanol_CO2;
	dn_CO2_aceton = dfcdw_aceton*stech_raceton_CO2;
    
    % H2O
    dn_H2O_metan = dfcdw_metan*stech_rmetan_H2O;
    dn_H2O_toluen = dfcdw_toluen*stech_rtoluen_H2O;
    dn_H2O_etanol = dfcdw_etanol*stech_retanol_H2O;
	dn_H2O_aceton = dfcdw_aceton*stech_raceton_H2O;
    
    % nove latkove toky % [mol/s] - nezahrnuje reakci CO !!! 
    % po zreagovani veskereho VOC ukonci cyklus
    n_metan = n_metan-dn_metan; if n_metan < 0; n_metan = 0; end % [mol/s]
    n_toluen = n_toluen-dn_toluen; if n_toluen < 0; n_toluen = 0; end % [mol/s]
    n_etanol = n_etanol-dn_etanol; if n_etanol < 0; n_etanol = 0; end % [mol/s]
    n_aceton = n_aceton-dn_aceton; if n_aceton < 0; n_aceton = 0; end % [mol/s]
    n_N2 = n_N2; % [mol/s] 
    % n_Ar = n0_Ar; % [mol/s]
    n_CO = n_CO; % [mol/s]
    n_CO2 = n_CO2+dn_CO2_metan+dn_CO2_toluen+dn_CO2_etanol+dn_CO2_aceton; % [mol/s]
    n_O2 = n_O2-dn_O2_metan-dn_O2_toluen-dn_O2_etanol-dn_O2_aceton; % [mol/s]
    n_H2O = n_H2O+dn_H2O_metan+dn_H2O_toluen+dn_H2O_etanol+dn_H2O_aceton; % [mol/s]
    n_SUM = n_metan+n_toluen+n_etanol+n_aceton+n_N2+n_CO+n_CO2+n_O2+n_H2O; % [-]
    
    % nove slozeni smesi (molarni zlomky) % [-] - nezahrnuje reakci CO !!!, % nezahrnuje Ar 
    y_metan = n_metan/n_SUM; % [-]
    y_toluen = n_toluen/n_SUM; % [-]
    y_etanol = n_etanol/n_SUM; % [-]
    y_aceton = n_aceton/n_SUM; % [-]
    y_N2 = n_N2/n_SUM; % [-]
    % y_Ar = n_Ar/n_SUM; % [-];
    y_CO = n_CO/n_SUM; % [-];
    y_CO2 = n_CO2/n_SUM; % [-];
    y_O2 = n_O2/n_SUM; % [-];
    y_H2O = n_H2O/n_SUM; % [-];
    
    % Vypocet stredni molekulove hmotnosti smesi [kg/kmol]
    Mw_str = Mw_metan*y_metan+Mw_toluen*y_toluen+Mw_etanol*y_etanol+Mw_aceton*y_aceton+Mw_CO2*y_CO2+Mw_CO*y_CO+Mw_O2*y_O2+Mw_H2O*y_H2O+Mw_N2*y_N2; % [kg/kmol]
    %+Mw_Ar*y_Ar
    
    % prirustek teploty vlivem reakcniho tepla a pokles teploty vlivem tepelnych ztrat
    [cp_metan,cp_toluen,cp_etanol,cp_aceton,cp_CO2,cp_O2,cp_H2O,cp_N2,dHr_metan,dHr_toluen,dHr_etanol,dHr_aceton] = reakcne_teplo(temp_K); % vypocet tepelnych kapacit a reakcneho tepla 
    
    dCp = (cp_metan*n_metan+cp_toluen*n_toluen+cp_etanol*n_etanol+cp_aceton*n_aceton+cp_CO2*n_CO2+cp_O2*n_O2+cp_H2O*n_H2O+cp_N2*n_N2); % [J.K-1.s-1]
    
    dtempdw = (-dHr_metan*rate_metan-dHr_toluen*rate_toluen-dHr_etanol*rate_etanol-dHr_aceton*rate_aceton)/dCp; % [K.kg_cat-1]
    
            % vypocet parametru pro vypocet dQ a dp
            [hustota,viskozita,lambda_spaliny] = hustota_viskozita(temp_K,p,y0_CO2,y0_O2,y0_H2O);
            V_real =n0_SUM*R*temp_K/p; % [m3/s]
            velocity = V_real/area_reactor ; % [m/s]
            Reynolds = dp_cat*velocity*hustota/viskozita; % [-]
    
            % vypocet tepelnych a teplotnych ztrat (dQ,dT)
            [dQ_ztraty,dT_ztraty] = tepelne_ztraty(hustota,viskozita,Reynolds,dw_cat,rho_cat,area_reactor,cp_SUM,dp_cat,dr,temp_K,lr,Mw_str,m0_SUM,lambda_spaliny,w_cat_max,dtempdw,lreactor);
    
            dQ_ztraty_SUM = dQ_ztraty_SUM+dQ_ztraty; % [W] pricitaji se tepelne ztraty na danem kroku 
            dT_ztraty_SUM = dT_ztraty_SUM+dT_ztraty; % [K] pricitaji se poklesy teploty tepelnymi ztratami na danem kroku
    
            temp_K = temp_K+dtempdw*dw_cat-dT_ztraty; % nova teplota zvysena vlivem reakcniho tepla a snizena o tepelne ztraty
    
    % entalpie smesi plynu a cp
    [iCO,iO2,iH2O,iCO2,iN2,iC7H8,iC3H6O,iC2H7O,iCH4] = entalpie(temp_K);
    i_SUM = iCO*y_CO+iO2*y_O2+iH2O*y_H2O+iCO2*y_CO2+iN2*y_N2+iC7H8*y_toluen+iC3H6O*y_aceton+iC2H7O*y_etanol+iCH4*y_metan; % [J/mol]
    cp_SUM = i_SUM/(temp_K-273.15); % [J.mol-1.K-1]

    % prirustek tlaku vlivem tlakove ztraty
        
        % Vypocet modifikovaneho Reynoldsova cisla pro vypocet dp
        G_area = m0_SUM/area_reactor; % [kg.s-1.m-2] - hmotnostni tok vztazen na prut. plochu
        Rem = dp_cat*G_area/viskozita; % [-] - modifikovane Reynoldsove cislo
        
        % Vypocet tlakove ztraty
        [N,fm] = tlakova_ztrata_Levy (Rem);
        theta = 1; % [-] - tvarovy faktor (pro koule = 1)
        dp_l = (fm*2*G_area^2*(1-porosity_cat)^(3-N))/(dp_cat*hustota*theta^(3-N)*porosity_cat^3); % [Pa/m]
        dpdw = dp_l*lr*dw_cat/w_cat_max; % [Pa] - tlakova ztrata v integracnom kroku
        p = p+dpdw; % [Pa] - system je v podtlaku takze +dpdw
        
    % prirastok katalyzatoru o krok dw
    w_cat = w_cat+dw_cat; % [kg]
    
    % ulozeni dat v danem kroku
    np = np+1;
    temp_K_1(np) = temp_K; %[K] - teplota 
    temp_C_1(np) = temp_K-273.15; %[C] - teplota 
    dtempdw_1(np-1) = dtempdw; % [K.kg_cat-1]
    dpdw_1(np-1) = dpdw; % [Pa] - tlakova ztrata v integracnom kroku
    hustota_1(np) = hustota; % [kg/m^3]
    viskozita_1(np) = viskozita; % [Pa.s]
    lambda_spaliny_1(np) = lambda_spaliny; % [W.m-1.K-1]
    velocity_1(np) = velocity; % [m/s]
    Re_1(np) = Reynolds; % [-]
    Rem_1(np) = Rem; % [-]
    p_1(np) = p; % [Pa]
    i_SUM_1 (np) = i_SUM; % [J/mol]
    cp_SUM_1(np) = cp_SUM; % [J.mol-1.K-1]
    Mw_str_1(np) = Mw_str; % [kg/kmol]
    
    dQ_ztraty_1(np-1) = dQ_ztraty; % [W] tepelne ztraty v danem kroku
    dQ_ztraty_SUM_1(np-1) = dQ_ztraty_SUM; % [W] celkove kumulujici se tepelne ztraty
    dT_ztraty_1(np-1) = dT_ztraty; % [K] pokles teploty tepelnymi ztratami v danem kroku
    dT_ztraty_SUM_1(np-1) =  dT_ztraty_SUM; % [K] celkove kumulujici se poklesy teploty tepelnymi ztratami

    y_metan_1(np) = y_metan; % [-] (molarni zlomek) 
    y_toluen_1(np) = y_toluen; % [-] (molarni zlomek)
    y_etanol_1(np) = y_etanol; % [-] (molarni zlomek)
    y_aceton_1(np) = y_aceton; % [-] (molarni zlomek)
    y_N2_1(np)= y_N2; % [-] (molarni zlomek)
    % y_Ar_1(np) = y_Ar; % [-] (molarni zlomek)
    y_CO_1(np) = y_CO; % [-] (molarni zlomek)
    y_CO2_1(np) = y_CO2; % [-] (molarni zlomek)
    y_O2_1(np) = y_O2; % [-] (molarni zlomek)
    y_H2O_1(np) = y_H2O; % [-] (molarni zlomek)
    
    n_metan_1(np) = n_metan;  % [mol/s] - latkove toky
    n_toluen_1(np) = n_toluen; % [mol/s]
    n_etanol_1(np) = n_etanol; % [mol/s]
    n_aceton_1(np) = n_aceton; % [mol/s]
    n_N2_1(np) = n_N2; % [mol/s]
    %n_Ar_1(np) =  n_Ar; % [mol/s]
    n_CO_1(np) = n_CO; % [mol/s]
    n_CO2_1(np) = n_CO2; % [mol/s]
    n_O2_1(np) = n_O2; % [mol/s]
    n_H2O_1(np) = n_H2O; % [mol/s]
    
    x_metan_1(np) = (y0_metan-y_metan_1(np))/(y0_metan); % [-] stupen konverze
    x_toluen_1(np) = (y0_toluen-y_toluen_1(np))/(y0_toluen); % [-] stupen konverze
    x_etanol_1(np) = (y0_etanol-y_etanol_1(np))/(y0_etanol); % [-] stupen konverze
    x_aceton_1(np) = (y0_aceton-y_aceton_1(np))/(y0_aceton); % [-] stupen konverze
    
    end
end   

%% grafy

% vytvoreni zkraceneho vektoru pro promenne dX.. , ktere maji o 1 prvek min
vektor_m_cat_short = vektor_m_cat(2:end); 

% nove vektory (transformace z kg_cat na %_cat)
    steps = w_cat_max/dw_cat; % počet kroku integrace
    step = dw_cat/w_cat_max; % krok integrace
    vektor_percent_cat = 0:(100/steps):100; 
    vektor_percent_cat_short = vektor_percent_cat(2:end); % vytvoreni zkraceneho vektoru

figure1 = figure;
subplot(2,1,1)
plot(vektor_percent_cat,x_metan_1,'Marker','*','LineStyle',':','Color','red'); hold all;
plot(vektor_percent_cat,x_toluen_1,'Marker','*','LineStyle',':','Color','blue'); 
plot(vektor_percent_cat,x_etanol_1,'Marker','*','LineStyle',':','Color','green'); 
plot(vektor_percent_cat,x_aceton_1,'Marker','*','LineStyle',':','Color','cyan'); 
xlabel({'Katalyzátor [%]'});
ylabel({'Konverze VOC [-]'});
grid on; 
legend1 = legend('metan','toluen','etanol','aceton','location','SouthEast');

subplot(2,1,2)
plot(vektor_percent_cat_short,rate_metan_1,'Marker','*','LineStyle',':','Color','red'); hold all;
plot(vektor_percent_cat_short,rate_toluen_1,'Marker','*','LineStyle',':','Color','blue');
plot(vektor_percent_cat_short,rate_etanol_1,'Marker','*','LineStyle',':','Color','green');
plot(vektor_percent_cat_short,rate_aceton_1,'Marker','*','LineStyle',':','Color','cyan');
xlabel({'Katalyzátor [%]'});
ylabel({'Reakční rychlost [mol.s-1.kg kat-1]'});
grid on; 
legend1 = legend('metan','toluen','etanol','aceton','location','SouthEast');

figure2 = figure;
subplot(3,2,1)
plot(vektor_percent_cat,y_metan_1,'Marker','*','LineStyle',':','Color','red'); hold all;
plot(vektor_percent_cat,y_toluen_1,'Marker','*','LineStyle',':','Color','blue');
plot(vektor_percent_cat,y_etanol_1,'Marker','*','LineStyle',':','Color','green');
plot(vektor_percent_cat,y_aceton_1,'Marker','*','LineStyle',':','Color','cyan');
xlabel({'Katalyzátor [%]'});
ylabel({'Molární zlomek VOC [-]'});
grid on; 
legend2 = legend('metan','toluen','etanol','aceton','location','NorthEast');

subplot(3,2,2)
plot(vektor_percent_cat,y_O2_1,'Marker','*','LineStyle',':','Color',[1 0 0]);
xlabel({'Katalyzátor [%]'});
ylabel({'Koncentrace O2 [-]'});
grid on; 
legend3 = legend('O2','location','NorthEast');

subplot(3,2,3)
plot(vektor_percent_cat,y_CO2_1,'Marker','*','LineStyle',':','Color',[1 0 0]);
xlabel({'Katalyzátor [%]'});
ylabel({'Koncentrace CO2 [-]'});
grid on; 
legend4 = legend('CO2','location','SouthEast');

subplot(3,2,4)
plot(vektor_percent_cat,y_H2O_1,'Marker','*','LineStyle',':','Color',[1 0 0]);
xlabel({'Katalyzátor [%]'});
ylabel({'Koncentrace H2O [-]'});
grid on; 
legend5 = legend('H2O','location','SouthEast');

subplot(3,2,5)
plot(vektor_percent_cat,temp_C_1,'Marker','*','LineStyle',':','Color',[1 0 0]);
xlabel({'Katalyzátor [%]'});
ylabel({'Teplota plynu [°C]'});
grid on; 
legend6 = legend('teplota','location','SouthEast');

subplot(3,2,6)
plot(vektor_percent_cat,p_1,'Marker','*','LineStyle',':','Color',[1 0 0]);
xlabel({'Katalyzátor [%]'});
ylabel({'Tlak plynu [Pa]'});
grid on; 
legend7 = legend('tlak','location','NorthWest');

figure3 = figure;
subplot(3,2,1)
plot(vektor_percent_cat,i_SUM_1,'Marker','*','LineStyle',':','Color','red');
xlabel({'Katalyzátor [%]'});
ylabel({'Entalpie plynu [J/mol]'});
grid on; 

subplot(3,2,2)
plot(vektor_percent_cat,cp_SUM_1,'Marker','*','LineStyle',':','Color','red');
xlabel({'Katalyzátor [%]'});
ylabel({'cp plynu [J.mol-1.K-1]'});
grid on;

subplot(3,2,3)
plot(vektor_percent_cat,Mw_str_1,'Marker','*','LineStyle',':','Color','red');
xlabel({'Katalyzátor [%]'});
ylabel({'Mw  plynu [kg.kmol-1]'});
grid on;

subplot(3,2,4)
plot(vektor_percent_cat,viskozita_1,'Marker','*','LineStyle',':','Color','red');
xlabel({'Katalyzátor [%]'});
ylabel({'Viskozita plynu [Pa.s]'});
grid on;

subplot(3,2,5)
plot(vektor_percent_cat_short,dT_ztraty_SUM_1,'Marker','*','LineStyle',':','Color','red');
xlabel({'Katalyzátor [%]'});
ylabel({'Kumulovaná teplotní ztráta [K]'});
grid on;

subplot(3,2,6)
plot(vektor_percent_cat,lambda_spaliny_1,'Marker','*','LineStyle',':','Color','red');
xlabel({'Katalyzátor [%]'});
ylabel({'Tepelná vodivost plynu [W.m-1.K-1]'});
grid on;

    