function [sys,x0]=genpet(T,X,u,flag,Nt,alphas,M0,yselect,Rr)
global coldata dyndata R
%
%Full Petlyuk column
%
%this function connects the different parts of the column
%the model consists of six sections described by gensect.m, one reboiler
%(reboiler.m) and three accumulators (accum.m)
%
%	x		- the states of each section stacked on top of each other
%	u		- the input parameters
%   Rr       - Weight matrix with [0,1] to make the measurement for the
%   combination
if nargin<8
   yselect=0;R=[];
end

nc=length(alphas);
if (length(Nt)==1)
   Nt=Nt*ones(6,1);
end
% Begr
% Splitter tilstander og holdup Nominelle verdier holdup
alphas=ones(7,1)*alphas;
M0W=0.5*M0;		%main column, along wall
M0P=0.5*M0;		%prefractionator
M0R=20*M0;		%reboiler holdup
M0T=20*M0;		%Top accumulator holdup
M0L=M0;			%liquid split accumulator
M0S=M0;			%side product accumulator
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Nominelle strømningsdata
L0=dyndata.L0;			% Reflux rate 
VB=dyndata.V0;			% Boilup rate 
S=dyndata.S0;			% Sideproduct removal rate 
Rl=dyndata.Rl0;			% Liquid split fraction to feed tray side 
Rv=dyndata.Rv0;			% Vapour split fraction to feed tray side 
F=dyndata.F0;			% Feed flow rate
zF=dyndata.zF;	% Feed composition for components 1:nc-1
qF=dyndata.q;		% Feed liquid fraction

% Computing liquid and vapour flows

L3_1=L0*Rl;
L3_4=L0*(1-Rl);
L1=L3_1;
L4=L3_4;
L2=L1+qF*F;
L5=L4-S;
L6=L2+L5;

V6_2=Rv*VB;
V6_5=(1-Rv)*VB;
V2=V6_2;
V5=V6_5;
V1=V2+(1-qF)*F;
V4=V5;
V3=V1+V4;

%Dynamiske data
hydraulics=0;
Tc=1.5;
k1=1/0.063;%0.26
k2=0;
hydpar1=[k1,k2,L1,V1,0];
hydpar2=[k1,k2,L2,V2,0];
hydpar3=[k1,k2,L0,V3,0];
hydpar4=[k1,k2,L4,V4,0];
hydpar5=[k1,k2,L5,V5,0];
hydpar6=[k1,k2,L6,VB,0];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%defining index matrixes
I1=[1:(nc-1)*Nt(1)];
I2=[1:(nc-1)*Nt(2)]+max(I1);
I3=[1:(nc-1)*Nt(3)]+max(I2);
I4=[1:(nc-1)*Nt(4)]+max(I3);
I5=[1:(nc-1)*Nt(5)]+max(I4);
I6=[1:(nc-1)*Nt(6)]+max(I5);
IR=[1:(nc-1)]+max(I6);
IT=[1:(nc-1)]+max(IR);
IL=[1:(nc-1)]+max(IT);
IS=[1:(nc-1)]+max(IL);
if (flag==1)
    x=X(1:104);
    holdup=X(105:end);
    %find derivatives
   %below this if-block
elseif (flag==0)				%initialization
   load Xineq;load mnom;
   nx=length(xineq)+length(holdup);
   nu=11;
   if (yselect==0)
      ny=3;
   elseif (yselect==1)
      ny=51;
   elseif (yselect==2)
      ny=13;
   elseif (yselect==4)
       ny=7;
   end
   sys=[nx;0;ny;nu;0;1];
   x0=[xineq;holdup];
   return
elseif (flag==2)
    sys=[];
    return;
elseif (flag==3)				%Measurment
     holdup=X(105:156);
     x=X(1:104);
    ur0=[zeros(nc,1);u(2)];
	ug0=zeros(2*nc,1);
%   [sys,x0]=sfunc(T,x,u,flag,kfactors,M0,hydraulics,Tc)
	yr=reboiler(T,[x(IR);holdup(8*6+1)],ur0,3,alphas(7,:),M0R,hydraulics,Tc);
	%[sys,x0]=sfunc(T,x,u,flag,Ntray,kfactors,M0,hydpar,yselect)
    y6=gensec(T,[x(I6);holdup(8*5+1:8*6)],ug0,3,Nt(6),alphas(6,:),M0,hydpar6,[]);
	y5=gensec(T,[x(I5);holdup(8*4+1:8*5)],ug0,3,Nt(5),alphas(5,:),M0W,hydpar5,[]);
	y4=gensec(T,[x(I4);holdup(8*3+1:8*4)],ug0,3,Nt(4),alphas(4,:),M0W,hydpar4,[]);
	y3=gensec(T,[x(I3);holdup(8*2+1:8*3)],ug0,3,Nt(3),alphas(3,:),M0,hydpar3,[]);
	y2=gensec(T,[x(I2);holdup(8+1:8*2)],ug0,3,Nt(2),alphas(2,:),M0P,hydpar2,[]);
    y1=gensec(T,[x(I1);holdup(1:8)],ug0,3,Nt(1),alphas(1,:),M0P,hydpar1,[]);
    yt=tank(T,[x(IT);holdup(8*6+2)],ur0,3,nc,M0,hydraulics,Tc);
   
     II1=reshape(I1,Nt(1),nc-1);
     II2=reshape(I2,Nt(2),nc-1);
     II3=reshape(I3,Nt(3),nc-1);
     II4=reshape(I4,Nt(4),nc-1);
     II5=reshape(I5,Nt(5),nc-1);
     II6=reshape(I6,Nt(6),nc-1);
     
     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;
     %Calculates the temp.profile, based on Antoines eq.
     if (yselect~=3)
     q=xp(:,1:3)*coldata.alfa';
     pr=coldata.P./q;
     Temp=coldata.B./(log10(pr)-coldata.A)+coldata.C;
     end
    
    if (yselect==0)			%only product composition sys=[xAD;xCB;xBS]
      sys=[y3(nc+2);(1-sum(yr(2:nc)));x(IS(2))];% y=[xA,xC,xB]
    elseif (yselect==1)
    sys=[y3(nc+2);(1-sum(yr(2:nc)));x(IS(2));Temp];% y=[xA,xC,xB,T]
    elseif (yselect==2)% comp, temp, h
    Y=[y3(nc+2);(1-sum(yr(2:nc)));x(IS(2));Temp;yt(3);yr(7)];
    c=coldata.R.*Y;
    cx1=xp(coldata.Rx(1),:);
    cx2=xp(coldata.Rx(2),:);
    sys=nonzeros([c;cx1';cx2']);% y=[xA,xC,xB,T,holdup]
    elseif (yselect==3)
    sys=xp;
    elseif (yselect==4)
    Temp_diff=[Temp(coldata.Rdiff(1,2))-Temp(coldata.Rdiff(1,1));Temp(coldata.Rdiff(2,2))-Temp(coldata.Rdiff(2,1))];
    sys=[y3(nc+2);(1-sum(yr(2:nc)));x(IS(2));Temp_diff;yt(3);yr(7)];    
    end
   return
else
   sys=[];
   return
end

L=u(1);			% Reflux rate 
VB=u(2);			% Boilup rate 
S=u(3);			% Sideproduct removal rate 
Rl=u(4);			% Liquid split fraction to feed tray side 
Rv=u(5);			% Vapour split fraction to feed tray side 
F=	u(6);			% Feed flow rate
zF=u(7:5+nc);	% Feed composition for components 1:nc-1
qF=u(6+nc);		% Feed liquid fraction
D=u(10);
B=u(11);

% holdup=[ones(coldata.Nt,1)*M0P;ones(coldata.Nt,1)*M0P;ones(coldata.Nt,1)*M0;ones(coldata.Nt,1)*M0W;ones(coldata.Nt,1)*M0W;ones(coldata.Nt,1)*M0;M0R;M0T;M0L;M0S];
%         1                    8 9                   16 17                24  
% Væskedynamikk
% Seksjon 3, topp
L3=hydpar3(3)+hydpar3(1)*(holdup(17:24)-M0)-hydpar3(2)*(VB-hydpar3(4));
dm3 = [L-L3(1);L3(1:8-1)-L3(2:8)];
% Seksjon 1, prefraksjonator
L1=hydpar1(3)+hydpar1(1)*(holdup(1:8)-M0P)-hydpar1(2)*(VB-hydpar1(4));
dm1 = [L3(8)*Rl-L1(1);L1(1:8-1)-L1(2:8)];
% Seksjon 2, prefraksjonator
L2=hydpar2(3)+hydpar2(1)*(holdup(9:16)-M0P)-hydpar2(2)*(VB-hydpar2(4));
dm2 = [L1(8)+F*qF-L2(1);L2(1:8-1)-L2(2:8)];
% Sekjson 4, kolonne ved vegg
L4 = hydpar4(3)+hydpar4(1)*(holdup(25:32)-M0W)-hydpar4(2)*(VB-hydpar4(4));
dm4 = [L3(8)*(1-Rl)-L4(1);L4(1:8-1)-L4(2:8)];
% Seksjon 5, kolonne ved vegg
L5 = hydpar5(3)+hydpar5(1)*(holdup(33:40)-M0W)-hydpar5(2)*(VB-hydpar5(4));
dm5 = [L4(8)-S-L5(1);L5(1:8-1)-L5(2:8)];
% Seksjon 6, bunn
L6 = hydpar6(3)+hydpar6(1)*(holdup(41:48)-M0)-hydpar6(2)*(VB-hydpar6(4));
dm6 = [L2(8)+L5(8)-L6(1);L6(1:8-1)-L6(2:8)];

%Computing liquid and vapour flows

% L3_1=L0*Rl;
% L3_4=L0*(1-Rl);
% L1=L3_1;
% L4=L3_4;
% L2=L1+qF*F;
% L5=L4-S;
% L6=L2+L5;

V6_2=Rv*VB;
V6_5=(1-Rv)*VB;
V2=V6_2;
V5=V6_5;
V1=V2+(1-qF)*F;
V4=V5;
V3=V1+V4;

%computing vapour compositions for flows connecting the different sections
ur0=[zeros(nc,1);u(2)];
ug0=zeros(2*nc,1);
%  [sys,x0]=sfunc(T,x,u,flag,kfactors,M0,hydraulics,Tc)
yr=reboiler(T,[x(IR);holdup(8*6+1)],ur0,3,alphas(7,:),M0R,hydraulics,Tc);
%[sys,x0]=sfunc(T,x,u,flag,Ntray,kfactors,M0,hydpar,yselect)
y6=gensec(T,[x(I6);holdup(8*5+1:8*6)],ug0,3,Nt(6),alphas(6,:),M0,hydpar6,[]);
y5=gensec(T,[x(I5);holdup(8*4+1:8*5)],ug0,3,Nt(5),alphas(5,:),M0W,hydpar5,[]);
y4=gensec(T,[x(I4);holdup(8*3+1:8*4)],ug0,3,Nt(4),alphas(4,:),M0W,hydpar4,[]);
y3=gensec(T,[x(I3);holdup(8*2+1:8*3)],ug0,3,Nt(3),alphas(3,:),M0,hydpar3,[]);
y2=gensec(T,[x(I2);holdup(8+1:8*2)],ug0,3,Nt(2),alphas(2,:),M0P,hydpar2,[]);
y1=gensec(T,[x(I1);holdup(1:8)],ug0,3,Nt(1),alphas(1,:),M0P,hydpar1,[]);
%generating input vectors for all sections
%reboiler:[Lin;xin(1:(nc-1));VB]
ur=[L6(8);y6(2:nc);VB;B];
%sections 1-6:[Lin;xin;Vin;yin]
u1=[L3(8)*Rl;x(IL);V1;y2(nc+2:2*nc)];

xi2=(L1(8)*y1(2:nc)+F*zF-(1-qF)*F*y2(nc+2:2*nc))/L2(1);%assumed equilibrium on feed stage
u2=[L1(8)+F*qF;xi2;V6_2;y6(nc+2:2*nc)];

xi3=x(IT);
yi3=(V1*y1(nc+2:2*nc)+V4*y4(nc+2:2*nc))/V3;
u3=[L;xi3;V3;yi3];

u4=[L3(8)*(1-Rl);x(IL);V5;y5(nc+2:2*nc)];

u5=[L4(8)-S;x(IS);V6_5;y6(nc+2:2*nc)];

xi6=(L2(8)*y2(2:nc)+L5(8)*y5(2:nc))/L6(1);
u6=[L6(1);xi6;VB;yr(nc+2:2*nc)];

%Liquid split accumulator
ul=[L3(8);y3(2:nc);L1(1)+L4(1)];
%Top accumulator
ut=[V3;y3(nc+2:2*nc);L;D];
%side-accumulator
us=[L4(8);y4(2:nc);L5(1)+S];   
%computing the derivatives
% [sys,x0]=sfunc(T,x,u,flag,kfactors,M0,hydraulics,Tc)
[dxr,dmr]=reboiler(T,[x(IR);holdup(8*6+1)],ur,1,alphas(7,:),M0R,hydraulics,Tc);
% [sys,x0]=sfunc(T,x,u,flag,Ntray,kfactors,M0,hydpar,yselect)
[dx6]=gensec(T,[x(I6);holdup(8*5+1:8*6)],u6,1,Nt(6),alphas(6,:),M0,hydpar6,[]);
[dx5]=gensec(T,[x(I5);holdup(8*4+1:8*5)],u5,1,Nt(5),alphas(5,:),M0W,hydpar5,[]);
[dx4]=gensec(T,[x(I4);holdup(8*3+1:8*4)],u4,1,Nt(4),alphas(4,:),M0W,hydpar4,[]);
[dx3]=gensec(T,[x(I3);holdup(8*2+1:8*3)],u3,1,Nt(3),alphas(3,:),M0,hydpar3,[]);
[dx2]=gensec(T,[x(I2);holdup(8+1:8*2)],u2,1,Nt(2),alphas(2,:),M0P,hydpar2,[]);
[dx1]=gensec(T,[x(I1);holdup(1:8)],u1,1,Nt(1),alphas(1,:),M0P,hydpar1,[]);
% [sys,x0]=sfunc(T,x,u,flag,nc,M0,hydraulics,Tc)
[dxt,dmt]=tank(T,[x(IT);holdup(8*6+2)],ut,1,nc,M0T,hydraulics,0);
[dxl]=accum(T,[x(IL);holdup(8*6+3)],ul,1,nc,M0L,hydraulics,0);
[dxs]=accum(T,[x(IS);holdup(8*6+4)],us,1,nc,M0S,hydraulics,0);
% Holdup i accum
dmL=L3(8)-L1(1)-L4(1);
dmS=L4(8)-S-L5(1);
% De deriverte
sys=[dx1;dx2;dx3;dx4;dx5;dx6;dxr;dxt;dxl;dxs;dm1;dm2;dm3;dm4;dm5;dm6;dmr;dmt;dmL;dmS];

