% initial.m
% M-file to set the initial conditions, and the numerical parameters
delt = 0.01
tfinal = 1
nplot = 1
nstep = tfinal/(nplot*delt)
%revise tfinal in case the there is round-off error
tfinal = nstep*nplot*delt
delx = 0.02
nx = 101
x(1)=0;
for ix=1:nx-1
	x(ix+1) = x(ix)+delx;
end
cinit(1)=inlet(0.);
for ix=2:nx
%	cinit(ix)=0;
	if x(ix)<0.25
		cinit(ix) = 0.;
	elseif x(ix)>0.75
		cinit(ix) = 0.;
	else
		cinit(ix) = 1.;
	end
end
cinit
display('Initial condition is shown above and plotted; press any key to continue.')
plot(x,cinit)
axis([x(1) x(nx) -.5 1.5])
xlabel('x')
ylabel('c')
title('Initial Concentration')
pause
for it=1:nstep*nplot+1
	time=delt*(it-1);
	y(it) = inlet(time);
	xd(it) = time;
end
display('Inlet concentration vs. time is plotted; press any key to continue.')
plot(xd,y)
xlabel('time')
ylabel('c')
title('Inlet Concentration versus time')
pause
%******************************************************************


% inlet.m
function y=inlet(t)
% enter with t = time
% exit with inlet = inlet concentration at that time
y = 0; %******************************************************************


% param.m
% Set the parameters for Langmuir adsorption
phi = 0.485
Kads = 2
%gamma = 0.1
gamma = 2
%******************************************************************


% macflux.m
% MacCormack method with flux correction
% Adsorption (equilibrium) with Langmuir isotherm
% enter with 
	% x(i), i = 1:nx, delx and nx
	% c(i), i = 1:nx, concentration at time t1
	% nstep, the number of time steps to compute for
	% delt, the timestep
	% parameters phi, K, and gamma
% exit with c(i), i = 1:nx for time t1+nstep*delt

for j=2:nx
	csoln(j) = c(j);
end
% add the ficticious point nx+1
csoln(nx+1)=csoln(nx);
para=delt/delx;
eps=delt*(1-phi)/phi;
for ii=1:nstep
	time = time + delt;
	csoln(1)=inlet(time);
	for ix=1:nx+1
   		pp(ix)=1+gamma*(1-phi)/(phi*(1+Kads*csoln(ix))^2);
	end
	for ix=2:nx
   		cnew(ix)=csoln(ix)-para*.5*(csoln(ix+1)-csoln(ix))*(1/pp(ix+1)+1/pp(ix));
	end
	cnew(1)=csoln(1);
	cnew(nx+1) = cnew(nx);
	for ix=1:nx
  	 	pp(ix)=1+gamma*(1-phi)/(phi*(1+Kads*cnew(ix))^2);
	end
	for ix=2:nx
   		ctemp(ix)=.5*(csoln(ix)+cnew(ix))-para*(cnew(ix)-cnew(ix-1))*(1/pp(ix)+1/pp(ix-1))/4;
	end
	ctemp(1) = csoln(1);

	% The following are the Flux-corrected transport steps.  If you don't want them
	% just put ctemp into csoln as in the next three statements.
%	for ix=1:nx
%		csoln(ix)=ctemp(ix)
%	end

	% csoln is c at n, ctemp is c-squiggle,calculate cnew to be c-hat
	for ix=2:nx
	% use eta = 1/8
   		cnew(ix)=ctemp(ix)+(csoln(ix+1)-2*csoln(ix)+csoln(ix-1))/8;
	end
	cnew(1)=csoln(1);
	% anti-diffusion step, calculate the deltas
	for ix=1:nx
		fcdel(ix)=cnew(ix+1)-cnew(ix);
	end
	fcdel;
	% calcualte the fluxes
	for ix = 1:nx-1
   		delhat=(ctemp(ix+1)-ctemp(ix))/8;
  		% use eta = 1/8 again
  		if (delhat>=0)
			sign=1;
		else
  			sign=-1;			
		end
		if ix>1
   			amin= min([sign*fcdel(ix-1),abs(delhat),sign*fcdel(ix+1)]);
		else
			% effectively sets fcdel(0) = 0
			amin = 0;
		end
  		fc(ix) = sign * max ([0 amin]);
	end
	fc(nx)=0;
	for ix=2:nx
		csoln(ix)=cnew(ix)-(fc(ix)-fc(ix-1));
	end
	csoln(nx+1) = csoln(nx);
% end of flux correction step

end
% return answer to c(j), j=1:n
for j=1:nx
	c(j) = csoln(j);
end
%******************************************************************


% runcode.m

% find initial conditions and set numerical parameters
initial
param

% put initial conditions in the working array and the first plotting array
for i=1:nx
	c(i) = cinit(i);
	cplot(1,i)=c(i);
end
time = 0.;

% calculate for each plot
for i=1:nplot
	macflux
	% put solution into next plotting array
	for j=1:nx
		cplot(i+1,j)=c(j);
	end
end

% plot solutions
if nplot==1
		plot(x,cplot(1,:),'r',x,cplot(2,:),'b')
	elseif nplot==2
		plot(x,cplot(1,:),'r',x,cplot(2,:),'b',x,cplot(3,:),'g')
	elseif nplot==3
		plot(x,cplot(1,:),'r',x,cplot(2,:),'b',x,cplot(3,:),'g',x,cplot(4,:),'k')
	elseif nplot==4
		plot(x,cplot(1,:),'r',x,cplot(2,:),'b',x,cplot(3,:),'g',x,cplot(4,:),'k',x,cplot(5,:),'m')
	else
		plot(x,cplot(1,:),'r',x,cplot(nplot-4,:),'b',x,cplot(nplot-3,:),'g',x,cplot(nplot-2,:),'k',x,cplot(nplot-1,:),'m',x,cplot(nplot,:),'c')
end
	
axis([x(1) x(nx) -.5 1.5])
xlabel('x')
ylabel('c')
title (['Solution with phi = ',num2str(phi),', K = ',num2str(Kads),', gamma = ',num2str(gamma)])