function [sys,x0]=genpet(T,x,u,flag,Nt,alphas,M0,yselect,R)
global coldata
%
%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
%   R       - 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

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



%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)
   %find derivatives
   %below this if-block
elseif (flag==0)				%initialization
   load nyx;
   x=X;%xnomopt;
   nx=length(x);
   nu=9;
   if (yselect==0)
      ny=3;
   elseif (yselect==1)
      ny=51+9+3; % Added the u's as measurments and the compositions for the output streams
   elseif (yselect==2)
      ny=5; % All controlled variables
   elseif (yselect==3)
     ny=48; % All states
   end
   sys=[nx;0;ny;nu;0;1];
   x0=x;
   return
elseif (flag==3)				%Measurment
   ur0=[zeros(nc,1);u(2)];
	ug0=zeros(2*nc,1);

	yr=reboiler(T,x(IR),ur0,3,alphas(7,:),M0R);
	y6=gensec(T,x(I6),ug0,3,Nt(6),alphas(6,:),M0,[]);
	y5=gensec(T,x(I5),ug0,3,Nt(5),alphas(5,:),M0W,[]);
	y4=gensec(T,x(I4),ug0,3,Nt(4),alphas(4,:),M0W,[]);
	y3=gensec(T,x(I3),ug0,3,Nt(3),alphas(3,:),M0,[]);
	y2=gensec(T,x(I2),ug0,3,Nt(2),alphas(2,:),M0P,[]);
	y1=gensec(T,x(I1),ug0,3,Nt(1),alphas(1,:),M0P,[]);
   
     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; % Old calculation
%    q=xp*coldata.alfa';
%    pr=coldata.P./q;
% D=u(2)-u(1);%D
% B=u(6)-u(3)-D;%B

if (yselect~=3)
     Temp=coldata.B./(log10(coldata.P./(xp*coldata.alfa'))-coldata.A)+coldata.C; 	
end
   if (yselect==0)			%only product composition sys=[xAD;xCB;xBS;D;B]
     sys=[y3(nc+2);(1-sum(yr(2:nc)));x(IS(2))];
   elseif (yselect==1) % Add input vector as measurement
    sys=[y3(nc+2);(1-sum(yr(2:nc)));x(IS(2));Temp;u];
   elseif (yselect==2)
     Y=[y3(nc+2);(1-sum(yr(2:nc)));x(IS(2));Temp;u];
     C1=R(1,:)*Y;  % Controlled variable 1
     C2=R(2,:)*Y;  % Controlled variable 2
     sys=[y3(nc+2);(1-sum(yr(2:nc)));x(IS(2));C1;C2];
   elseif (yselect==3)
     sys=xp; % Returnes the states
   end
   return
else
   sys=[];
   return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

L0=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);


%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);

yr=reboiler(T,x(IR),ur0,3,alphas(7,:),M0R);
y6=gensec(T,x(I6),ug0,3,Nt(6),alphas(6,:),M0,[]);
y5=gensec(T,x(I5),ug0,3,Nt(5),alphas(5,:),M0W,[]);
y4=gensec(T,x(I4),ug0,3,Nt(4),alphas(4,:),M0W,[]);
y3=gensec(T,x(I3),ug0,3,Nt(3),alphas(3,:),M0,[]);
y2=gensec(T,x(I2),ug0,3,Nt(2),alphas(2,:),M0P,[]);
y1=gensec(T,x(I1),ug0,3,Nt(1),alphas(1,:),M0P,[]);

%generating input vectors for all sections
%reboiler:[Lin;xin(1:(nc-1));VB]
ur=[L6;y6(2:nc);VB];
%sections 1-6:[Lin;xin;Vin;yin]
u1=[L3_1;x(IL);V1;y2(nc+2:2*nc)];

xi2=(L1*y1(2:nc)+F*zF-(1-qF)*F*y2(nc+2:2*nc))/L2;%assumed equilibrium on feed stage
u2=[L1+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=[L0;xi3;V3;yi3];

u4=[L3_4;x(IL);V5;y5(nc+2:2*nc)];

u5=[L4-S;x(IS);V6_5;y6(nc+2:2*nc)];
xi6=(L2*y2(2:nc)+L5*y5(2:nc))/L6;
u6=[L6;xi6;VB;yr(nc+2:2*nc)];

%Liquid split accumulator
ul=[L0;y3(2:nc)];
%Top accumulator
ut=[V3;y3(nc+2:2*nc)];
%side-accumulator
us=[L4;y4(2:nc)];   
%computing the derivatives
dxr=reboiler(T,x(IR),ur,1,alphas(7,:),M0R);
dx6=gensec(T,x(I6),u6,1,Nt(6),alphas(6,:),M0,[]);
dx5=gensec(T,x(I5),u5,1,Nt(5),alphas(5,:),M0W,[]);
dx4=gensec(T,x(I4),u4,1,Nt(4),alphas(4,:),M0W,[]);
dx3=gensec(T,x(I3),u3,1,Nt(3),alphas(3,:),M0,[]);
dx2=gensec(T,x(I2),u2,1,Nt(2),alphas(2,:),M0P,[]);
dx1=gensec(T,x(I1),u1,1,Nt(1),alphas(1,:),M0P,[]);
dxt=accum(T,x(IT),ut,1,nc,M0T);
dxl=accum(T,x(IL),ul,1,nc,M0L);
dxs=accum(T,x(IS),us,1,nc,M0S);

sys=[dx1;dx2;dx3;dx4;dx5;dx6;dxr;dxt;dxl;dxs];

