% Three-isothermal CSTRs in series  (Luyben)
%  This prgram simulates a three-isothermal CSTRs
%  Assumptions
%  Levels of three reactor are perfectly controlled
%  Disturbance is a step change in feed concentraion

% initilization
  ttime=0;trand_dat=[];tsamp_t=0;
  %Design parameters /residence time of reactor 1,2,and 3/
   Tau1=2;Tau2=2;Tau3=3;
   Reac_rate=.5;
  % Parameter setting for Euler method
   delta_t=.01; endtime=10;samp_t=.1;
  %Initial state condition
   Ca1=.4;Ca2=.2;Ca3=.08; Ca0=0.8;
   % Disturbance occurs
    Ca0=1.8;
   trand_dat=[ttime Ca1 Ca2 Ca3]
   tsamp_t=samp_t;
  %Simulation Starts 
  while ttime<endtime,
  % Evaluation of derivative
    Ca1_dot=(Ca0-Ca1)/Tau1-Reac_rate*Ca1;
	Ca2_dot=(Ca1-Ca2)/Tau2-Reac_rate*Ca2;
	Ca3_dot=(Ca2-Ca3)/Tau3-Reac_rate*Ca3;
	%
	ttime=ttime+delta_t;
 	Ca1=Ca1+Ca1_dot*delta_t;
	Ca2=Ca2+Ca2_dot*delta_t;
	Ca3=Ca3+Ca3_dot*delta_t;
	
 if ttime>=tsamp_t,
   xtemp=[ttime Ca1 Ca2 Ca3]
   trand_dat=[trand_dat;xtemp];
   tsamp_t=tsamp_t+samp_t;
   end;
  end;
 % this is the end of simulation
  subplot(311);plot(trand_dat(:,1),trand_dat(:,2));
  title('concentration in the reactor 1');
  subplot(312);plot(trand_dat(:,1),trand_dat(:,3));
  title('concentration in the reactor 2');
  subplot(313);plot(trand_dat(:,1),trand_dat(:,4));
  title('concentration in the reactor 3');
  end;  
   
   
