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

function [sys,x] = genkaibel(t,x,u,flag,yselect)

global coldata dyndata


if (flag == 0)
    % Initialization
    nx = length(dyndata.x_nom);
    nu = 13;
    if (yselect == 0)
        ny = 4;
    elseif (yselect == 1)
        ny = 7 + sum(coldata.Nt);
    elseif (yselect == 2)
        ny = 4;
    elseif (yselect == 3)
        ny = nx;
    end
    
    sys = [nx 0 ny nu 0 0]';
    x = dyndata.x_nom;
    return;
elseif (flag == 1)
    % Derivatives
    L = u(1); VB = u(2);
    S1 = u(3); S2 = u(4);
    Rl = u(5); Rv = u(6);
    F = u(7); zF = u(8:10); qF = u(11);
    D = u(12); B = u(13);
    
    % Index matrices
    indmat;
    
    % Compute liquid flows for each section downwards
    L3 = compute_liquid(3,L,[x(im_m3); x(im_T3)]);
    L1 = compute_liquid(1,Rl * L3(end),[x(im_m1); x(im_T1)]);
    L2 = compute_liquid(2,L1(end) + qF * F,[x(im_m2); x(im_T2)]);
    L4 = compute_liquid(4,(1 - Rl) * L3(end),[x(im_m4); x(im_T4)]);
    L7 = compute_liquid(7,L4(end) - S1,[x(im_m7); x(im_T7)]);
    L5 = compute_liquid(5,L7(end) - S2,[x(im_m5); x(im_T5)]);
    L6 = compute_liquid(6,L2(end) + L5(end),[x(im_m6); x(im_T6)]);
    
    % Compute vapour flows for each section upwards
    V6 = compute_vapour(6,[VB; VB],x(im_T6));
    V2 = compute_vapour(2,[Rv * V6(1); Rv * V6(2)],x(im_T2));
    V1 = compute_vapour(1,[V2(1) + (1 - qF) * F; V2(2)],x(im_T1));
    V5 = compute_vapour(5,[(1 - Rv) * V6(1); (1 - Rv) * V6(2)],x(im_T5));
    V7 = compute_vapour(7,[V5(1); V5(2)],x(im_T7));
    V4 = compute_vapour(4,[V7(1); V7(2)],x(im_T4));
    V3 = compute_vapour(3,[V1(1) + V4(1); V1(2) + V4(2)],x(im_T3));
    
    % Change in holdup
    dm3 = holdup_change(L3,V3,L,[V1(1) + V4(1); V1(2) + V4(2)]);
    dm1 = holdup_change(L1,V1,Rl * L3(end),[V2(1) + (1 - qF) * F; V2(2)]);
    dm2 = holdup_change(L2,V2,L1(end) + qF * F,[Rv * V6(1); Rv * V6(2)]);
    dm4 = holdup_change(L4,V4,(1 - Rl) * L3(end),[V7(1); V7(2)]);
    dm7 = holdup_change(L7,V7,L4(end) - S1,[V5(1); V5(2)]);
    dm5 = holdup_change(L5,V5,L7(end) - S2,[(1 - Rv) * V6(1); (1 - Rv) * V6(2)]);
    dm6 = holdup_change(L6,V6,L2(end) + L5(end),[VB; 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)],[zeros(4,1); VB],3);
    y1 = gensec(t,[x(im_x1); x(im_m1); x(im_T1)],zeros(12,1),3,L1,V1,coldata.Nt(1));
    y2 = gensec(t,[x(im_x2); x(im_m2); x(im_T2)],zeros(12,1),3,L2,V2,coldata.Nt(2));
    y3 = gensec(t,[x(im_x3); x(im_m3); x(im_T3)],zeros(12,1),3,L3,V3,coldata.Nt(3));
    y4 = gensec(t,[x(im_x4); x(im_m4); x(im_T4)],zeros(12,1),3,L4,V4,coldata.Nt(4));
    y5 = gensec(t,[x(im_x5); x(im_m5); x(im_T5)],zeros(12,1),3,L5,V5,coldata.Nt(5));
    y6 = gensec(t,[x(im_x6); x(im_m6); x(im_T6)],zeros(12,1),3,L6,V6,coldata.Nt(6));
    y7 = gensec(t,[x(im_x7); x(im_m7); x(im_T7)],zeros(12,1),3,L7,V7,coldata.Nt(7));
    
    [dxr,dmr,Tr,dTr] = reboiler(t,[x(im_xr); x(im_Tr); x(im_mr)],[L6(end); y6(7:9); VB; B],1);
    [dx1,T1,dT1] = gensec(t,[x(im_x1); x(im_m1); x(im_T1)],[Rl * L3(end); x(im_xL); V2(1) + (1 - qF * F); V2(2); y2(1:6)],1,L1,V1,coldata.Nt(1));
    [dx2,T2,dT2] = gensec(t,[x(im_x2); x(im_m2); x(im_T2)],[L1(end) + qF * F; (L1(end) * y1(7:9) + zF * F - (1 - qF) * F * y2(1:3))/L2(1); Rv * V6(1:2); y6(1:6)],1,L2,V2,coldata.Nt(2));
    [dx3,T3,dT3] = gensec(t,[x(im_x3); x(im_m3); x(im_T3)],[L; x(im_xt); V1(1) + V4(1); V1(2) + V4(2); (V1(1) * y1(1:3) + V4(1) * y4(1:3))/V3(end); (V1(2) * y1(4:6) + V4(2) * y4(4:6))/V3(end)],1,L3,V3,coldata.Nt(3));
    [dx4,T4,dT4] = gensec(t,[x(im_x4); x(im_m4); x(im_T4)],[(1 - Rl) * L3(end); x(im_xL); V7(1:2); y7(1:6)],1,L4,V4,coldata.Nt(4));
    [dx5,T5,dT5] = gensec(t,[x(im_x5); x(im_m5); x(im_T5)],[L7(end) - S2; x(im_xS2); (1 - Rv) * V6(1:2); y6(1:6)],1,L5,V5,coldata.Nt(5));
    [dx6,T6,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; VB; yr(1:3)'; yr(1:3)'],1,L6,V6,coldata.Nt(6));
    [dx7,T7,dT7] = gensec(t,[x(im_x7); x(im_m7); x(im_T7)],[L4(end) - S1; x(im_xS1); V5(1:2); y5(1:6)],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];
    return;
elseif (flag == 2)
    sys = [];
    return;
elseif (flag == 3)
    % Measurements
    
    indmat;

    temperatures = [x(im_T1); x(im_T2); x(im_T3); x(im_T4); x(im_T5); x(im_T6); x(im_T7); x(im_Tr)];
    
    y3 = gensec(t,[x(im_x3); x(im_m3); x(im_T3)],zeros(12,1),3,[],[],coldata.Nt(3));

    if (yselect == 0)
        sys = [
            y3(1)
            1 - sum(x(im_xr))
            x(im_xS1(2))
            x(im_xS2(3))
            ];
        return;
    elseif (yselect == 1)
        sys = [
            y3(1)             % x_D
            1 - sum(x(im_xr)) % x_B
            x(im_xS1(2))      % x_S1
            x(im_xS2(3))      % x_S2
            x(im_mt)
            x(im_mr)
            temperatures
            ];
        return;
    elseif (yselect == 2)
        sys = [
            temperatures(17)
            temperatures(30)
            temperatures(59)
            temperatures(49)
            ];
        return;
    elseif (yselect == 3)
        sys = x;
        return;
    else
        sys = [];
        return;
    end
else
    sys = [];
    return;
end