%***************************************************************
%***Simulacion de un tanque con control de nivel  ******************
%*** y un mezclador con control de composicion y nivel ******************
%*** la salidas de ambos estan en una proporcion de  2 a 1 ******************
%***M. Rodriguez 2004                               ******************
%FALTA :  poner la equivalencia de a=B y b=C
%****************************************************************

clear 
global area g SP F_setp Kc taui onoff F_in1 F_in2 area_mix CA_in1 CB_in2 ...
area_salida setp_cb SP_m F_setpm Kc_m taui_m onoff SP_h Kc_h taui_h ratio ... 
k K0 EA T R rho SP_r

i=0;

%*********MIXER*************
F_in1=5;%0.01055*10000;     %m3/min
CA_in1=300;    %kmol/m3
CB_in2=240;    %kmol/m3
area_mix=3;        %m2
area_salida=0.016;  %m2 corresponde a unas 8''
g=9.8;         %m/s2
SP_m=0.4; %fraccion molar de B
Kc_m=2;
taui_m=1;
onoff=1;
F_setpm=10;
SP_h=4;
Kc_h=3;
taui_h=2;
ratio=0.5;


%***********TANQUE************
area=4; %m2
F_setp=4;
SP=0.7;%m
Kc=3;
taui=2;

%***********REACTOR************
k=6;  %
K0=[75000000 8500000];
EA=[78240 45605];
R=8.314;       %
T=568;      %K
rho=450;      %kmol/m3
SP_r=2.1;

y_initial=[1 0 1 2 0.8 0 0 2 1 10 6  1 0]; %altura, conc A, conc B 
t_start = 0;						%	min
t_final = 60;						%	min
tspan = [t_start t_final];
[t,y] = ode15s('proceso',tspan,y_initial);

for i=1:length(y(:,1))
F_inx(i)=F_setp+Kc*(SP-y(i,1))+onoff*Kc/taui*y(i,2);
end
figure
subplot(2,1,1),
plot(t,y(:,1)),
title('altura del tanque'),
xlabel('time (min)'),
ylabel('h (m)')
subplot(2,1,2),
plot(t,F_inx),
title('caudal del tanque'),
xlabel('time (min)'),
ylabel('Fin m3/min')

figure
for i=1:length(y(:,3))
F_outxm(i)=F_setpm-Kc_h*(SP_h-y(i,3))-onoff*Kc_h/taui_h*y(i,7); %m3/min
setp_cb=y(i,4)*SP_m/(1-SP_m);
F_in2x(i)=Kc_m*(setp_cb-y(i,5))+onoff*Kc_m/taui_m*y(i,6);
if F_in2x(i)<0
    F_in2x(i)=0;
end
xb(i)=y(i,5)/(y(i,4)+y(i,5));
end


subplot(2,3,1),
plot(t,y(:,3)),
title('altura del tanque'),
xlabel('time (h)'),
ylabel('h (m)')
subplot(2,3,2),
plot(t,F_outxm),
title('caudal del tanque'),
xlabel('time (h)'),
ylabel('Fout (m)')
subplot(2,3,3),
plot(t,y(:,4)),
title('CA'),
xlabel('time (h)'),
ylabel('CA ')
subplot(2,3,4),
plot(t,y(:,5)),
title('CB'),
xlabel('time (h)'),
ylabel('CB')
subplot(2,3,5),
plot(t,F_in2x),
title('Fin2'),
xlabel('time (h)'),
ylabel('Fin2')
subplot(2,3,6),
plot(t,xb),
title('fraccion molar B'),
xlabel('time (h)'),
ylabel('xb')

figure
subplot(2,3,1),
plot(t,y(:,8)),
title('volumen del tanque'),
xlabel('time (min)'),
ylabel('V (m3)')
subplot(2,3,2),
plot(t,y(:,11)),
title('NC'),
xlabel('time (min)'),
ylabel('NC kmol')
subplot(2,3,3),
plot(t,y(:,9)),
title('NA'),
xlabel('time (min)'),
ylabel('NA kmol')
subplot(2,3,4),
plot(t,y(:,10)),
title('NB'),
xlabel('time (min)'),
ylabel('NB kmol')
subplot(2,3,5),
plot(t,y(:,12)),
title('ND'),
xlabel('time (min)'),
ylabel('ND kmol')

pause
close all