[sys x0]=descompunere_CaCO3(t,x,u,flag)

function [sys x0]=descompunere_CaCO3(t,x,u,flag)
if abs(flag)==1
    
    % Intrari
    Qg=u(1);
    Qs=u(2);
    Qc=u(3);
    Ts0=u(4);
    TgL=u(5);
    xCaCO3_0=u(6);
    xCaO_0=0.1;
    xCO2_L=0.01;
    
    % Iesiri (necunoscute)
    for i=1:30
        xCaCO3(i)=x(i);
        xCaO(i)=x(30+i);
        xCO2(i)=x(60+i);
        Ts(i)=x(90+i);
        Tg(i)=x(120+i);
    end;
    
    % Constantele proceseului
    M_CaCO3=100;                    %[kg/kmol]
    M_CaO=56;
    M_CO2=44;
    R=2*4.185;                      % [J/(kmolK)]
    Ea=0.74*10^6;                   % [kJ/kmol]
    A_CaCO3=1.25*3.5*1.26*10^26;    % factor preexponential, legea Arrhenius [1/s]
    dH_CaCO3=420*4.185;             % [kJ/kg]
    Cpg=0.25*4.185;
    Cps=0.35*4.185;                 % [J/kg]
    Kt=1/3600*440*4.185;            % [kJ/(sm2K)]
    At=300;                         % [m2]
    ro=2830;                        % [kg/m3]
    rog=1.2;            
    raport_ro=1.75;
    
    % Parametrii cuptorului
    hr=6;                           % inaltimea zonei reactie [m]
    d=4.5;                          % diametrul cuptorului [m]
    vg=1.5;                         % viteza gazului [m/s]
    vs=0.4062/3600;                 % viteza solidului [m/s]
    
    % Ecuatii algebrice
    dz=hr/30;                       % grosimea elementului dz
    A=pi*d^2/4;
    
    % Ecuatiile diferentiale
    k=A_CaCO3*exp(-Ea/(R*Ts(1)));
    Hgs=At*Kt*(Tg(1)-Ts(1)) /A;
    Hr=k*ro*xCaCO3(1)/M_CaCO3*dH_CaCO3;
    Hc=Qc*5600*4.185/(hr*A);
    
    dxCaCO3(1)=-vs*(xCaCO3(1)-xCaCO3_0)/dz-xCaCO3(1)*k;
    dxCaO(1)=-vs*(xCaO(1)-xCaO_0)/dz + xCaCO3(1)*k*(M_CaO/M_CaCO3)*raport_ro;
    dxCO2(1)=-vg*(xCO2(1)-xCO2(2))/dz +xCaCO3(1)*k*(M_CO2/M_CaCO3)*ro/rog;
    dTs(1)=-vs*(Ts(1)-Ts0)/dz + vs/(Qs*Cps)*(Hgs-Hr);
    dTg(1)=-vg*(Tg(1)-Tg(2))/dz-vg/(Qs*Cps)*(Hgs-Hc);
    
    for i=2:29
        k=A_CaCO3*exp(-Ea/(R*Ts(i)));
        Hgs=At*Kt*(Tg(i)-Ts(i))/A;
        Hr=k*ro*xCaCO3(i)/M_CaCO3*dH_CaCO3;
        Hc=Qc*5600*4.185/(hr*A);
    
        dxCaCO3(i)=-vs*(xCaCO3(i+1)-xCaCO3(i-1))/(2*dz)-xCaCO3(i)*k;
        dxCaO(i)=-vs*(xCaO(i+1)-xCaO(i-1))/(2*dz)+xCaCO3(i)*k*(M_CaO/M_CaCO3)*raport_ro;
        dxCO2(i)=-vg*(xCO2(i-1)-xCO2(i+1))/(2*dz)+xCaCO3(i)*k*(M_CO2/M_CaCO3)*ro/rog;
        dTs(i)=-vs*(Ts(i+1)-Ts(i-1))/(2*dz) + vs/(Qs*Cps)*(Hgs-Hr);
        dTg(i)=-vg*(Tg(i-1)-Tg(i+1))/(2*dz)-vg/(Qs*Cps)*(Hgs-Hc);
    end;
    
    k=A_CaCO3*exp(-Ea/(R*Ts(30)));
    Hgs=At*Kt*(Tg(30)-Ts(30))/A;
    Hr=k*ro*xCaCO3(30)/M_CaCO3*dH_CaCO3;
    Hc=Qc*5600*4.185/(hr*A);
    
    dxCaCO3(30)=-vs*(xCaCO3(30)-xCaCO3(29))/dz-xCaCO3(30)*k;
    dxCaO(30)=-vs*(xCaO(30)-xCaO(29))/dz+xCaCO3(30)*k*(M_CaO/M_CaCO3)*raport_ro;
    dxCO2(30)=-vg*(xCO2(30)-xCO2_L)/dz+ xCaCO3(30)*k*(M_CO2/M_CaCO3)*ro/rog;
    dTs(30)=-vs*(Ts(30)-Ts(29))/dz+vs/(Qs*Cps)*(Hgs-Hr);
    dTg(30)=-vg*(Tg(30)-TgL)/dz-vg/(Qs*Cps)*(Hgs-Hc);
    
sys=[dxCaCO3 dxCaO dxCO2 dTs dTg];   
elseif flag==3
    sys=x;
elseif flag==0
    sys=[150 0 150 6 0 0];
   %for i=1:30
   %    x1(i)=0.9-(0.9-0.07)*i/30;
   %    x2(i)=0.01+(0.89)*i/30;
   %    x3(i)=0.38-(0.37)*i/30;
   %    x4(i)=1033+100*i/30;
   %    x5(i)=1083+100*i/30;
   % end;
   load('init.mat')
   x0=[x1 x2 x3 x4 x5];
end;