%% main function
function [sys Temp] = Petlyuk_NLmodel(t,x,U)

global dyndata coldata
load Hel3
load csel3
%L = x(1); VB =  x(2);
%Rl=x(3);Rv=x(4);
%S1=x(5);
%F=x(6);
%zF = x(7:8); qF = x(9);
%D = x(10); B = x(11);
%x=(x(12:end));
L = U(1); VB =  U(2);
Rl=U(3);Rv=U(4);
S1=U(5);
F=U(6);
zF = U(7:8); qF = U(9);
D = U(10); B = U(11);

%x=(x(12:end));
Rv2=1; %No vapor side draw

% Index matrices
indmat_Yes_S1;
     
     
% Nominal Input Parameters
%L=dyndata.u_nom(1);			 % Reflux rate 
%VB=dyndata.u_nom(2);		 % Boilup rate 
%Rl=dyndata.u_nom(3);		 % Liquid split fraction to feed tray side 
%Rv=dyndata.u_nom(4);	     % Vapour split fraction to feed tray side 
%S=dyndata.u_nom(5);	     % Sideproduct removal rate 
%F=dyndata.u_nom(6);		 % Feed flow rate
%zF=dyndata.u_nom(7:8);	     % Feed composition for components 1:2
%qF=dyndata.u_nom(9);		 % Feed liquid fraction

Tboil=coldata.Tboil;
     
     xp1 =  x(II1);
     xp2 =  x(II2);
     xp3 =  x(II3);
     xp4 =  x(II4);
     xp5 =  x(II5);
     xp6 =  x(II6);
     xp1=[xp1 1-sum(xp1,2)]; 
     xp2=[xp2 1-sum(xp2,2)];
     xp3=[xp3 1-sum(xp3,2)];
     xp4=[xp4 1-sum(xp4,2)];
     xp5=[xp5 1-sum(xp5,2)];
     xp6=[xp6 1-sum(xp6,2)];
     xp=[xp1;xp2;xp3;xp4;xp5;xp6];
     Temp=xp*Tboil';
     %Temp(10,1) = dyndata.temp_nom(10);
     %Temp(38,1) = dyndata.temp_nom(38);
     %Temp(60,1) = dyndata.temp_nom(60);
     
     





     
    

% Compute liquid flows for each section downwards

L21 = compute_liquid(3,L,x(im_m21),dyndata.u_nom(1),dyndata.x_nom(im_m21));
L1 = compute_liquid(1,L21(end)*(Rl),x(im_m1),dyndata.u_nom(1)*(dyndata.u_nom(3)),dyndata.x_nom(im_m1));
L22 = compute_liquid(4,L21(end)*(1-Rl),x(im_m22),dyndata.u_nom(1)*(1-dyndata.u_nom(3)),dyndata.x_nom(im_m22));
L2 = compute_liquid(2,L1(end)+F*qF,x(im_m2),dyndata.u_nom(1)*(dyndata.u_nom(3))+F*qF,dyndata.x_nom(im_m2));
L23 = compute_liquid(5,L22(end)-S1,x(im_m23),dyndata.u_nom(1)*(1-dyndata.u_nom(3))-(dyndata.u_nom(5)),dyndata.x_nom(im_m23));
L24 = compute_liquid(6,L23(end)+L2(end),x(im_m24),...
    dyndata.u_nom(1)*(1-dyndata.u_nom(3))-(dyndata.u_nom(5))+dyndata.u_nom(1)*(dyndata.u_nom(3))+F*qF,dyndata.x_nom(im_m24));

% Compute vapour flows for each section
V24 = VB*ones(coldata.Nt(6),1);
V2 = V24(1)*(Rv)*ones(coldata.Nt(2),1);
V23 = V24(1)*(1-Rv)*ones(coldata.Nt(5),1);
V1 = (V2(1)+F*(1-qF))*ones(coldata.Nt(1),1);
V22 = V23(1)*Rv2*ones(coldata.Nt(4),1);
V21 = (V22(1)+V1(1))*ones(coldata.Nt(3),1);

% Change in holdup
dm21 = holdup_change(L21,V21,L,V22(1)+V1(1));
dm1 = holdup_change(L1,V1,L21(end)*(Rl),V2(1)+F*(1-qF));
dm22 = holdup_change(L22,V22,L21(end)*(1-Rl),V23(1)*Rv2);
dm2 = holdup_change(L2,V2,L1(end)+F*qF,V24(1)*(Rv));
dm23 = holdup_change(L23,V23,L22(end)-S1,V24(1)*(1-Rv));
dm24 = holdup_change(L24,V24,L23(end)+L2(end),VB);

% Change in compositions and temperatures
% Find vapour compositions for first two stages in the section and liquid
% composition of last stage
yr = reboiler(t,[x(im_xr); x(im_mr)],[[]; VB],3);
y1 = gensec(t,[x(im_x1); x(im_m1)],[],3,L1,V1,coldata.Nt(1));
y2 = gensec(t,[x(im_x2); x(im_m2)],[],3,L2,V2,coldata.Nt(2));

y21 = gensec(t,[x(im_x21); x(im_m21)],[],3,L21,V21,coldata.Nt(3));
y22 = gensec(t,[x(im_x22); x(im_m22)],[],3,L22,V22,coldata.Nt(4));
y23 = gensec(t,[x(im_x23); x(im_m23)],[],3,L23,V23,coldata.Nt(5));
y24 = gensec(t,[x(im_x24); x(im_m24)],[],3,L24,V24,coldata.Nt(6));

%find the derivatives
[dxr,dmr] = reboiler(t,[x(im_xr);x(im_mr)],[L24(end); y24(5:6); VB; B],1);
[dxt,dmt] = tank(t,[x(im_xt); x(im_mt)],[L; D; V21(1); y21(1:2)],1);

[dx1] = gensec(t,[x(im_x1); x(im_m1)],[L21(end)*(Rl); y21(5:6); V2(1) + (1 - qF) * F; (F*(1-qF)*zF+V2(1)*y2(1:2))/V1(end)],1,L1,V1,coldata.Nt(1));
[dx2] = gensec(t,[x(im_x2); x(im_m2)],[L1(end) + qF * F; (L1(end) * y1(5:6) + zF * F*qF)/L2(1);(Rv) * V24(1); y24(1:2)],1,L2,V2,coldata.Nt(2));

[dx21] = gensec(t,[x(im_x21); x(im_m21)],[L; y21(1:2); V22(1)+V1(1); (V1(1)*y1(1:2)+V22(1)*y22(1:2))/V21(end)],1,L21,V21,coldata.Nt(3));
[dx22] = gensec(t,[x(im_x22); x(im_m22)],[L21(end)*(1-Rl); y21(5:6);Rv2 * V23(1); y23(1:2)],1,L22,V22,coldata.Nt(4));
[dx23] = gensec(t,[x(im_x23); x(im_m23)],[L22(end)-S1; x(im_xS1);(1-Rv) * V24(1); y24(1:2)],1,L23,V23,coldata.Nt(5));
[dx24] = gensec(t,[x(im_x24); x(im_m24)],[L2(end)+L23(end); (L2(end) * y2(5:6) + L23(end) * y23(5:6))/L24(1);VB; yr'],1,L24,V24,coldata.Nt(6));

[dxS1,dmS1] = accum(t,[x(im_xS1); x(im_mS1)],[L22(end); y22(5:6); L23(1) + S1],1);


dm=  [dm1;dm2;dm21;dm22;dm23;dm24];
dx = [dx1; dx2;dx21;dx22;dx23;dx24];
%sys is a column matrix with all the derivatives of compositions and
%holdups

sys = [dx; dxr; dxt;dxS1;dm; dmr; dmt;dmS1];

     

return;
