% This function are used only for plotting concentration out of reactor

function c=concprof(t,T)
global mc mh Cp Cpc Q1 Q2 Q3 Npoint Tf cf dm1 dm2 dm3 dHrx u conc p 
% comment: For simplicity we use rx((T(i),c(i-1)) rather than rx((T(i),c(i))
 
 %section1.--Reconstruction of unstable ammonia reactor------%
if t<1200
   p=170;
elseif (t>=1200 & t<=7200)
  p=150;
else t>7200
  p=200;
end;
%------------------------------------------------------------%

%section2.--New steady-state---------------------------------%
%if t>1200
%   p=150;
%else
%   p=200;
%end;
%------------------------------------------------------------%

%section3.--P-controller-------------------------------------% 
 %Kc=3;
%Q1=16.11+Kc*(T(1)-369.4); % positiv feedback, e=y-r
%mc=mh-(Q1+Q2+Q3); 
%------------------------------------------------------------% 

m1=mc+Q1;
m2=m1+Q2;
m3=m2+Q3;

% Conc. in bed 1:
c(1)=cf+dm1*rx(T(1),cf)/m1;
for i=2:10
  c(i)=c(i-1)+dm1*rx(T(i),c(i-1))/m1;
end
% Quench 2: 
cin2=m1/m2*c(10)+Q2/m2*cf;
% Conc. in bed 2
c(11)=cin2+dm2*rx(T(11),cin2)/m2;
for i=12:20
  c(i)=c(i-1)+dm2*rx(T(i),c(i-1))/m2;
end
% Quench 3:
cin3=m2/m3*c(20)+Q3/m3*cf;
% Conc. in bed 3
c(21)=cin3+dm3*rx(T(21),cin3)/m3;
for i=22:30
    c(i)=c(i-1)+dm3*rx(T(i),c(i-1))/m3;
 end
 
