molecule_diffusion_simulation.m

% Simple 2D Diffusion Simulation
% Created 19 May 2023 by Adam Noel

% Parameters
N = 2000; % Number of molecules
D = 1; % Diffusion coefficient [m^2/s]
dt = 1; % Time step [s]
RX = [9,11,9,11]; % Receiver
numRepeat = 10; % Number of realizations/repeats
rng(1)

% Initialization
tFinal = 200;
numTime = length(0:dt:tFinal);
molecule_count = zeros(numRepeat,numTime);

% Simulation Loop
for n = 1:numRepeat
    molecule_position = zeros(N,2);
    tIndex = 2;
    tCur = 0;
    while tCur < tFinal    
        % Move individual molecules
        for i = 1:N
            molecule_position(i,1) = molecule_position(i,1) + sqrt(2*D*dt)*randn(1); % Update x-coordinate
            molecule_position(i,2) = molecule_position(i,2) + sqrt(2*D*dt)*randn(1); % Update y-coordinate
    
            % Check whether molecule is within RX
            if molecule_position(i,1) > RX(1) && molecule_position(i,1) < RX(2) ...
                    && molecule_position(i,2) > RX(3) && molecule_position(i,2) < RX(4)
                molecule_count(n,tIndex) = molecule_count(n,tIndex) + 1;
            end
        end
    
        % Update simulation time
        tCur = tCur + dt;
        tIndex = tIndex + 1;
    end
end

%% Plot number of molecules observed
t = 0:dt:(tIndex-2)*dt;
figure('Color', 'w');
molecule_watch = mean(molecule_count(:,1:length(t)),1);
plot(t,molecule_watch)
ylabel('Number of Molecules')
xlabel('Time (seconds)')
NExp = N*4./(4*pi*D*t).*exp(-10^2./2./D./t);
hold on
plot(t,NExp, 'LineWidth',2)