% File: pipe_1d_tracer.m
% Eli Tziperman and Ian Eisenman
% Solve the 1d advection diffusion equation:
%  d C/dt + u*d C/dx  = kappa d^2 C/dx^2
% plot and animate the solution.
% -----------------------------------------------------------------
close all; clear all; % close figures, clear any stores variables

%run animation (1 or 0)?
animate=1;

% parameters:
kappa=3.0e6;         % cm^2/sec
L=1000.0e5;          % length of pipe, cm
robert_filter_coeff = 0.25;

% period of integration:
year=360.0*86400.0;
T=5.0*year;               % seconds

% advecting velocity
u=2;                    % cm/sec

% finite differencing parameters:
nx=50;               % number of grid points along-pipe (along x direction)
dx=L/(nx-1);         % size of grid cell
dt=0.002*year;       % time step (seconds)
nt=floor(T/dt);      % number of time steps

% Initial conditions:
for i=1:nx
    for n=1:3
        C(i,n)= exp(-(i-nx/2)^2/(nx/8)^2);
    end
end
pipe2d=zeros(nt,nx); %concentration field
time=zeros(nt,1);

% Integrate in time:  n is one of:
% 1 - past    (known)
% 2 - present (known)
% 3 - future  (to be calculated)

for n=1:nt
  
% Loop over all space points and find future value
  for i=2:nx-1
    C_x=(C(i+1,2)-C(i-1,2))/(2*dx);
    C_xx=(C(i+1,1)-2*C(i,1)+C(i-1,1))/(dx*dx);
    C(i,3)=C(i,1)+2*dt*(-u*C_x+kappa*C_xx);
  end
  % cyclic boundary conditions for C at ends of pipe:
  C(1,3)=C(nx-1,3); 
  C(nx,3)=C(2,3);
  % Robert filter:
  C(:,2) = C(:,2)+robert_filter_coeff*(C(:,3) - 2*C(:,2) + C(:,1));
  % prepare for next time step:
  C(:,1)=C(:,2);
  C(:,2)=C(:,3);
  
  % store solution for plotting:
  time(n)=n*dt;
  pipe2d(n,:)=C(:,3)';

end

figure(1) %open figure 1 if it is not already open
clf %clear the figure
% plot concentration at center of pipe as function fo time:
subplot(2,1,1)
plot(time,pipe2d(:,nx/2))
title('C(t)');
xlabel('time (sec)');
ylabel('Concentration at center');

% contour tracer concetration as function of time and space:
subplot(2,1,2)
contour(pipe2d(:,:))
xlabel('distance')
ylabel('time')
title('Concentration')


% Second figure for animation
if animate
    figure(2)
    x=1:size(pipe2d,2);
    h=plot(x,pipe2d(1,:));xlabel('distance'),ylabel('Concentration')
    set(gca,'YLim',[0 1],'XLim',[0 50])
    title('t=1: enter to animate')
    set(h,'EraseMode','xor')
    pause
    % Animation loop
    for i=1:(size(pipe2d,1)/100):size(pipe2d,1)
        drawnow
        y=pipe2d(i,:);
        title(['t=' num2str(i)])
        set(h,'Ydata',y)
        pause(0.1)
    end
end
