% GENKAIBEL.M %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Kaibel column main file

function [dxdt] = Kaibel_dxdt(t,x)

global dyndata coldata

    L = x(1); VB = x(2);
    S1 = x(3); S2 = x(4);
    Rl = x(5); Rv = x(6);
    F = x(7); zF = x(8:10); qF = x(11);
    D = x(12); B = x(13);
    
x=(x(14:end));
    % Index matrices
    indmat;

    L3O=dyndata.u_nom(1);
    L1O=L3O*dyndata.u_nom(5);
    L2O=L1O+dyndata.u_nom(7)*dyndata.u_nom(11);
    L4O=L3O*(1-dyndata.u_nom(5));
    L7O=L4O-dyndata.u_nom(3);
    L5O=L7O-dyndata.u_nom(4);
    L6O=L2O+L5O;
    
    % Compute liquid flows for each section downwards
    L3 = compute_liquid(3,L,x(im_m3),L3O,dyndata.x_nom(im_m3));
    L1 = compute_liquid(1,Rl * L3(end),x(im_m1),L1O,dyndata.x_nom(im_m1));
    L2 = compute_liquid(2,L1(end) + qF * F,x(im_m2),L2O,dyndata.x_nom(im_m2));
    L4 = compute_liquid(4,(1 - Rl) * L3(end),x(im_m4),L4O,dyndata.x_nom(im_m4));
    L7 = compute_liquid(7,L4(end) - S1,x(im_m7),L7O,dyndata.x_nom(im_m7));
    L5 = compute_liquid(5,L7(end) - S2,x(im_m5),L5O,dyndata.x_nom(im_m5));
    L6 = compute_liquid(6,L2(end) + L5(end),x(im_m6),L6O,dyndata.x_nom(im_m6));
    
    % Compute vapour flows for each section
    V6 = VB*ones((coldata.Nt(6)),1);
    V2 = Rv * VB*ones((coldata.Nt(2)),1);
    V1 = (V2(1) + (1 - qF) * F)*ones((coldata.Nt(1)),1);
    V5 = (1 - Rv) * V6(1)*ones((coldata.Nt(5)),1);
    V7 = V5(1)*ones((coldata.Nt(7)),1);
    V4 = V7(1)*ones((coldata.Nt(4)),1);
    V3 = (V1(1) + V4(1))*ones((coldata.Nt(3)),1);
    
    % Change in holdup
    dm3 = holdup_change(L3,V3,L,[V1(1) + V4(1)]);
    dm1 = holdup_change(L1,V1,Rl * L3(end),[V2(1) + (1 - qF) * F]);
    dm2 = holdup_change(L2,V2,L1(end) + qF * F,[Rv * V6(1)]);
    dm4 = holdup_change(L4,V4,(1 - Rl) * L3(end),[V7(1)]);
    dm7 = holdup_change(L7,V7,L4(end) - S1,[V5(1)]);
    dm5 = holdup_change(L5,V5,L7(end) - S2,[(1 - Rv) * V6(1)]);
    dm6 = holdup_change(L6,V6,L2(end) + L5(end),[VB]);
    dm = [dm1; dm2; dm3; dm4; dm5; dm6; dm7];
    
    % Change in compositions and temperatures
    % Find vapour compositions first
    yr = reboiler(t,[x(im_xr); x(im_Tr); x(im_mr)],[[]; VB],3);
    y1 = gensec(t,[x(im_x1); x(im_m1); x(im_T1)],[],3,L1,V1,coldata.Nt(1));
    y2 = gensec(t,[x(im_x2); x(im_m2); x(im_T2)],[],3,L2,V2,coldata.Nt(2));
    y3 = gensec(t,[x(im_x3); x(im_m3); x(im_T3)],[],3,L3,V3,coldata.Nt(3));
    y4 = gensec(t,[x(im_x4); x(im_m4); x(im_T4)],[],3,L4,V4,coldata.Nt(4));
    y5 = gensec(t,[x(im_x5); x(im_m5); x(im_T5)],[],3,L5,V5,coldata.Nt(5));
    y6 = gensec(t,[x(im_x6); x(im_m6); x(im_T6)],[],3,L6,V6,coldata.Nt(6));
    y7 = gensec(t,[x(im_x7); x(im_m7); x(im_T7)],[],3,L7,V7,coldata.Nt(7));
    


    
    [dxr,dmr,~,dTr] = reboiler(t,[x(im_xr); x(im_Tr); x(im_mr)],[L6(end); y6(7:9); VB; B],1);
   
    U1=[L*(Rl); y3(7:9); V2(1) + (1 - qF) * F; y2(1:3)];
    U2=[L1(end) + qF * F; (L1(end) * y1(7:9) + zF * F -(1-qF)*F*y2(1:3))/L2(1);Rv* VB; y6(1:3)];
    
    [dx1,~,dT1] = gensec(t,[x(im_x1); x(im_m1); x(im_T1)],U1,1,L1,V1,coldata.Nt(1));
    [dx2,~,dT2] = gensec(t,[x(im_x2); x(im_m2); x(im_T2)],U2,1,L2,V2,coldata.Nt(2));
    
    [dx3,~,dT3] = gensec(t,[x(im_x3); x(im_m3); x(im_T3)],[L; x(im_xt); V1(1) + V4(1);(V1(1) * y1(1:3) + V4(1) * y4(1:3))/V3(end)],1,L3,V3,coldata.Nt(3));
    [dx4,~,dT4] = gensec(t,[x(im_x4); x(im_m4); x(im_T4)],[(1 - Rl) * L3(end); x(im_xL); V7(1); y7(1:3)],1,L4,V4,coldata.Nt(4));
    [dx5,~,dT5] = gensec(t,[x(im_x5); x(im_m5); x(im_T5)],[L7(end) - S2; x(im_xS2); (1 - Rv) * V6(1); y6(1:3)],1,L5,V5,coldata.Nt(5));
    [dx6,~,dT6] = gensec(t,[x(im_x6); x(im_m6); x(im_T6)],[L6(1); (L2(end) * y2(7:9) + L5(end) * y5(7:9))/L6(1); VB;yr(1:3)'],1,L6,V6,coldata.Nt(6));
    [dx7,~,dT7] = gensec(t,[x(im_x7); x(im_m7); x(im_T7)],[L4(end) - S1; x(im_xS1); V5(1); y5(1:3)],1,L7,V7,coldata.Nt(7));
    [dxt,dmt] = tank(t,[x(im_xt); x(im_mt)],[L; D; V3(1); y3(1:3)],1);
    [dxL,dmL] = accum(t,[x(im_xL); x(im_mL)],[L3(end); y3(7:9); L1(1) + L4(1)],1);
    [dxS1,dmS1] = accum(t,[x(im_xS1); x(im_mS1)],[L4(end); y4(7:9); L7(1) + S1],1);
    [dxS2,dmS2] = accum(t,[x(im_xS2); x(im_mS2)],[L7(end); y7(7:9); L5(1) + S2],1);
    
    dx = [dx1; dx2; dx3; dx4; dx5; dx6; dx7];
    dT = [dT1; dT2; dT3; dT4; dT5; dT6; dT7];
    
    sys = [dx; dxr; dxt; dxL; dxS1; dxS2; dT; dTr;dm; dmr; dmt; dmL; dmS1; dmS2 ];
    dxdt=sys;
    return;