% fdplanar.m
%	Finite-difference model of diffusion.
%	1-Dimensional planar model.
% 	Impulse of 1 unit deposited at x = 0 at t = 0.
%	REFERENCE calculation extends over full range.
%	TEST calculation extends between x1 and x2. 
%	Boundary conditions at x1 and x2 can be changed
%	and TEST's behavior at boundary can be compared
%	with REFERENCE's behavior at it matched boundary.
clear

%%%
% USER CHOICES
%%%
Uo = 1;        % 1 unit deposited in bin at x = 0 at t = 0;
dx = 0.0010;   % spatial bin size [cm]
xmax = 0.050;  % extent of ±xaxis plot [cm]
ymax = 50;     % extent of yaxis for plot
tmax = 10;     % [s] duration of simulation
pauselength = 1; % [s] duration of pause between timesteps
x1 = -0.010;   % front boundary 
x2 =  0.010;   % rear boundary 
INSULATE_FRONT = 1; % 1 = INSULATED FRONT BOUNDARY0, 0 = MATCHED BOUNDARY
INSULATE_REAR  = 0; % 1 = INSULATED REAR BOUNDARY, 0 = MATCHED BOUNDARY
alpha = 1.3e-3; % diffusivity [cm^2/s]. (this value is for heat in water.)
yymax = 400;
%%%
%%%

%%%%%%%%%%% PROGRAM %%%%%%%%%%%%%%
xmin = -xmax;
x = (xmin:dx:xmax)'; % [cm] position of spatial bins
N = length(x);  %  array length x[1:N] = -xmax:dx:xmax
b = 0.35;     % b must be less than 1/2 to have stable computation
dt = b/alpha*dx^2;  % stable timestep.  b = dt*alpha/dx^4

C = zeros(N,1); %  array of bins for heat C[1:N], calculation with no boundaries.
C(N/2) = Uo/dx;  % put C [J] of energy in zero position bin
CB = zeros(N,1); % array with boundary
CB(N/2) = Uo/dx;  % put C [J] of energy in zero position bin
i1 = round((x1 - xmin)/dx) + 1; % front boundary at x[i1]
i2 = round((x2 - xmin)/dx) + 1; % rear  boundary at x[i2]
if i1 < 2; i1 = 2; end
if i2 > N-1; i2 = N-1; end
x1 = x(i1) - dx/2
x2 = x(i2) + dx/2

figure(1); clf
t = 0;
plot(x,C,'o-','color', 'r')
hold on
plot(x(i1-1:i2+1),CB(i1-1:i2+1),'x-','color', 'b')
plot([x1 x1],[0 yymax],'b--')
plot([x2 x2],[0 yymax],'b--')
set(gca,'fontsize',18)
text(xmax*0.7,0.9*yymax,sprintf('%5.6f s', t),'fontsize',18)
xlabel('position z [cm]')
ylabel('Concentration [units/cm]')
title('Concentration [units/cm] = [units/cm^3][cm^2 area]')
axis([xmin xmax -yymax/4 yymax])
for i=1:N
	plot([x(i) x(i)]+dx/2, [0 yymax],'c--')
	plot([x(i) x(i)]+dx/2, [0 yymax],'c--')
end
pause(1)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Time step through computation.
%%%%
t      = 0;
dC     = zeros(N,1);
dCB    = zeros(N,1);
while t < tmax
	t = t + dt;
	denom = 2*alpha*t/dx;
	N4at = 4*alpha*t;
	
	%%%% BOUNDARIES
	xx1 = abs(x1);   % front boundary
	xx2 = abs(x2);   % rear  boundary
	xz2 = abs(x(2)); % front initial bin position
	xN2 = x(N-1);    % rear  final bin position
	xb1 = abs(x(i1)); % front bin position at x(i1)
	xb2 = abs(x(i2)); % rear  bin position at x(i2)
	
	% Matched bc for REFERENCE CASE (RED CURVE)
	% Use analytic solution for planar impulse at origin
	% to select C(1) and C(N) to mimic matched boundaries.
	r = x(1);
	C(1) = exp(-r^2/(4*alpha*t))/(4*pi*alpha*t)^(1/2); % front boundary
	r = x(N);
	C(N) = exp(-r^2/(4*alpha*t))/(4*pi*alpha*t)^(1/2); % rear boundary

	% DEFAULT for TEST CASE (BLUE CURVE)
	% Use analytic solution for planar impulse at origin
	% to select C(i1-1) and C(i2+1) to mimic matched boundaries.
	% This is a cheat since it only works for impulse at origin.
	% You need to devise a local rule to mimic a matched boundary.
	r = x(i1-1);
	CB(i1-1) = exp(-r^2/(4*alpha*t))/(4*pi*alpha*t)^(1/2); % front boundary
	r = x(i2+1);
	CB(i2+1) = exp(-r^2/(4*alpha*t))/(4*pi*alpha*t)^(1/2); % rear boundary
	
	%%%%%%%%%%%%%%%%%%%%%
	% STUDENT's solutions
	% to make the rear boundary condition behave as a matched boundary:
	% if 1 = enable, if 0 = disable
	
	if 0
		% Escample of insulating bc
		CB(i1-1) = CB(i1);
		CB(i2+1) = CB(i2);
	end

	if 0
		% Dan's linear approximation
		CB(i2+1) = 2*CB(i2) - CB(i2-1);
	end
		
	if 0
		% A 2nd order approximation
		Di2   = (CB(i2  ) - CB(i2-1))/dx; % neg
		Di2m1 = (CB(i2-1) - CB(i2-2))/dx; % more neg
		Dib = Di2 - (Di2m1 - Di2); % least neg
		CB(i2+1) = CB(i2) + Dib*dx;
	end
		
	if 0
		% Ted's 2nd order approximation
		% dCB(i) = b*(CB(i+1) - 2*CB(i) + CB(i-1)); 
		CB(i2+1) = 2*CB(i2) - CB(i2-1); % dan's
		u = (CB(i2-2+1) - 2*CB(i2-2) + CB(i2-2-1));
		v = (CB(i2-1+1) - 2*CB(i2-1) + CB(i2-1-1));
		CB(i2+1) = CB(i2+1) + (u - v)/dx;
	end
		
	if 0
		% Yinchu approximation
		CB(i2+1) = 3*CB(i2) -3*CB(i2-1) + CB(i2-2);
	end
		
	if 0
		% Laurel's approximation 
		if CB(i2-1) > 0
			CB(i2+1) = CB(i2)*CB(i2)/CB(i2-1); % <----- !!!!!!!
		else
			CB(i2+1) = 0;
		end
	end
		
		
	%%%%%%%%%%%%%%%%%%%%%
	% End student's solutions for rear boundary.	
		
	%%%% Calculate the movement of heat
	% b = dt*alpha/dx^2
	for i=2:N-1
		dC(i) = b*(C(i+1) - 2*C(i) + C(i-1)); 
	end
	for i=i1:i2
		dCB(i) = b*(CB(i+1) - 2*CB(i) + CB(i-1)); 
	end

	% Move the heat
	for i=2:N-1
		C(i) = C(i) + dC(i); 
	end
	for i=i1:i2
		CB(i) = CB(i) + dCB(i); 
	end

	
figure(1); clf
u = max(max(CB),max(C));
if u<yymax/2
	yymax = u;
end

for i=1:N
	plot([x(i) x(i)]+dx/2, [0 yymax],'c--')
	hold on
	plot([x(i) x(i)]+dx/2, [0 yymax],'c--')
end
plot(x,C,'o-','color', 'r')
%plot(x(i1:i2),CB(i1:i2),'o-','color', 'b')
plot(x(i1-1:i2+1),CB(i1-1:i2+1),'x-','color', 'b')
plot([x1 x1],[0 yymax],'b--')
plot([x2 x2],[0 yymax],'b--')
set(gca,'fontsize',18)
text(xmax*0.7,0.9*yymax,sprintf('%5.6f s', t),'fontsize',18)
xlabel('position z [cm]')
ylabel('Concentration [units/cm]')
title('Concentration [units/cm] = [units/cm^3][cm^2 area]')
plot(x(i1),C(i1),'ro','marker','square','markersize',3)
plot(x(i2),C(i2),'ro','marker','diamond','markersize',3)
axis([xmin xmax -yymax/4 yymax])

%%%%%
% ANALYTIC DIFFUSION THEORY (GREEN CURVE)
%%%%
plot(x, exp(-x.^2/(4*alpha*t))/(4*pi*alpha*t)^(1/2), 'g-')

drawnow
%pause(pauselength)
end % while

