% Fil dist3.m 
% Case 1: Konstant holdup
% zf=0.5;[t,x]=ode15s('dist3',[0 10],[zf*ones(1,22) 5 5 0])
% Case 2: PI-regulering topp
% Case 3: PI-regulering midt
% Case 4: P-regulering midt
% zf=0.5;[t,x]=ode15s('dist3',[0 10],[zf*ones(1,22) 0.01 9.99 0])
% Tid for aa naa gitt sammensetning i topp OG bunn, sammenligne anal.
%  se litt paa effekt av PI-tuninger (KC, taui)
% se paa ulike verdier av zf og alpha

function xdot = dist3(t,x)

% Konstant relativ flyktighet
alfa = 2;

% Antall trinn
NT=22;

% Foedetrinn

NF=11;

% Oevrige data

F=0;
qF= 0;
zF=0.5;
V=10;
L0=10; 
D=0;
B=0;
MD=x(NT+1);
MB=x(NT+2);
MK=0.0005;
KC=20;
%KC=0;
XDs=0.95;
XFs=0.5;

% Regulator
%err = x(NT)-XDs; % Case 2
err=x(NF)-XFs; % Case 3,4
LT=L0-KC*err;  % Case 4
L0I=0; taui=0.05;
interr = x(NT+3);
%LT = L0I - KC*err - KC/taui*interr; % Case 2,3
%LT=10; % Case 1

% Damp-vaeske likevekt
i=1:NT-1;    
y(i)=alfa*x(i)./(1+(alfa-1)*x(i));

% Dampstroem 
i=1:NT-1;    
V(i)=V*ones(1,NT-1);
i=NF:NT-1;
V(i)=V(i)+(1-qF)*F;

% Vaeskestroem
i=2:NT;
L(i)=LT*ones(1,NT-1);
i=2:NF;
L(i)=L(i)+qF*F;

% Total kondensator
y(NT)=x(NT);

% Holdup i kolonnen
i=2:NT-1;
M(i)=MK*ones(1,NT-2);

% Holdup i kokeren
MBdot=L(2)-V(1)-B;

% Holdup i kondensatoren
MDdot=V(NT-1)-L(NT)-D;

% Dynamiske ligninger (komponentbalanser)

% Kokeren
Mxdot(1) = L(2)*x(2)-V(1)*y(1)-B*x(1);
xdot(1)=(Mxdot(1)-x(1)*MBdot)/MB;

% Kolonnen
i=2:NF-1;
xdot(i)=(L(i+1).*x(i+1)' - L(i).*x(i)' + V(i-1).*y(i-1) - V(i).*y(i))./M(i);
i=NF+1:NT-1;
xdot(i)=(L(i+1).*x(i+1)' - L(i).*x(i)' + V(i-1).*y(i-1) - V(i).*y(i))./M(i);

% Foedetrinn
xdot(NF) =(L(NF+1)*x(NF+1)+V(NF-1)*y(NF-1)-L(NF)*x(NF)-V(NF)*y(NF)+F*zF)/MK;

% Kondensatoren
Mxdot(NT) = V(NT-1)*y(NT-1)-L(NT)*x(NT)-D*x(NT);
xdot(NT)=(Mxdot(NT)-x(NT)*MDdot)/MD;

% Funksjonen f(t,x) maa returnere en kolonnevektor
xdot(NT+1)=MDdot;
xdot(NT+2)=MBdot;

%Integrator in controller
xdot(NT+3)=err;

xdot=xdot(:);











