The following line in the MatLab code is incorrect: current_pw = pw_c;



%% RANDOM WALKTS, BROWNIAN MOTION, AND DIFFUSION


% ————————— parameters
width = 5;                                           % width of pipe in mm
a = 0.01;                                         % lattice constant in mm
pw = linspace(0, 0.9, 10);                        % probability of waiting
num_steps = 1000000;                              % # steps for histograms
num_walks = 1000000;                       % # random walks for histograms

% ————————— calculate lattice size based on width
lattice_size = round(width/a);

% ————————— initialize random walker position
x = 0;
y = lattice_size/2;

% ————————— initialize mean squared displacement variables
msd_x = zeros(1,num_steps);
msd_y = zeros(1,num_steps);

% loop over the different values of pw
for p = 1:length(pw)
% set current pw value
current_pw = pw_c;

% ————————— loop over the # of steps
for step = 2:num_steps+1
%           choice: does walker wait or make a step?
            if rand<=current_pw
%           walker waits
            continue;
            else
%           walker makes a step
%           calculate the step direction
            dx = randi([-1,1]);
            dy = randi([-1,1]);
%           update walker position
            x = x+dx;
            y = y+dy;
%           reflect at the y-direction boundary
            if y<1
                y = 2-y;
            elseif y>lattice_size
                y = 2*lattice_size - y - 1;
            end
%           calculate mean squared displacement in each direction
            msd_x(step) = (x*a)^2;
            msd_y(step) = ((y-(lattice_size/2))*a)^2;
      end
end

% —————————— 1) visualize the random walk trajectories
figure;
plot((0:num_steps)*a, (msd_y*a), 'b.')
title(sprintf('Random Walk Trajectories for pw = %.1f', current_pw));
title('Random Walk Trajectories');
xlabel('x (mm)');
ylabel('y (mm)');
xlim([0,num_steps*a]);

% ————————— 2) mean squared displacement in x-direction
figure;
n_steps = 1:num_steps;
msd_x = cumsum(msd_x) ./ (1:num_steps+1);
plot(n_steps(1:num_steps-1)*a, msd_x(2:num_steps)*a^2, 'b')
title(sprintf('Mean Squared Displacement in the x-Direction for pw = %.1f', current_pw));
xlabel('Number of Steps');
ylabel('Mean Squared Displacement');
xlim([0,num_steps*a]);

% ————————— 3) mean squared displacement in y-direction
figure;
msd_y = cumsum(msd_y) ./ (1:num_steps+1);
plot(n_steps(1:num_steps)*a, msd_y(1:num_steps)*a^2, 'b')
title(sprintf('Mean Squared Displacement in the y-Direction for pw = %.1f', current_pw));
xlabel('Number of Steps');
ylabel('Mean Squared Displacement');
xlim([0, num_steps*a]);

% ————————— 4) histogram of locations in x-direction
figure;
for num_steps_hist = [100,1000,10000,100000,1000000]
    x_hist = zeros(1,lattice_size+1);
    for walk = 1:num_walks
        x=0;
        for step_hist = 1:num_steps_hist
            dx = randi([-1,1]);
            x = x+dx;
            if x<0
                x=0;
            elseif x>lattice_size
                x = lattice_size;
            end
        end
        x_hist(x+1) = x_hist(x+1) + 1;
    end
    subplot(2,3,find([100,1000,10000,100000,1000000] == num_steps_hist));
    bar(0:lattice_size, x_hist / num_walks, 'b');
    title(sprintf('Histogram of x-locations after %d steps for pw = %.1f', num_steps_hist, current_pw));
    xlabel('x');
    ylabel('Probability');
end

% ————————— 5) histogram of locations in y-direction
figure;
for num_steps_hist = [100,1000,10000,100000,1000000]
    y_hist = zeros(1,lattice_size+1);
    for walk = 1:num_walks
        x = 0;
        y = lattice_size/2;
        for step_hist = 1:num_steps_hist
            dx = randi([-1,1]);
            dy = randi([-1,1]);
            x = x+dx;
            y = y+dy;
            if y<1
                y = 2-y;
            elseif y>lattice_size
                y = 2*lattice_size - y - 1;
            end
        end
        y_hist(y+1) = y_hist(y+1) + 1;
    end
    subplot(2,3,find([100,1000,10000,100000,1000000] == num_steps_hist));
    bar(0:lattice_size, y_hist / num_walks, 'b');
    title(sprintf('Histogram of y-locations after %d steps for pw = %.1f', num_steps_hist, current_pw));
    xlabel('y');
    ylabel('Probability');
end
end

% ————————— 6) calculate avg time per step for given diffusion constant
D_desired_media_A = 1e-7; % desired diffusion constant for media A in cm^2/sec
D_desired_media_B = 1e-8; % desired diffusion constant for media B in cm^2/sec
t_media_A = D_desired_media_A / (a^2); % avg time per step for media A in sec
t_media_B = D_desired_media_B / (a^2); % avg time per step for media B in sec

% ————————— 7) calculate diffusion constant in the x-direction
figure;
diffusion_x_media_A = gradient(msd_x*a^2) / (2*t_media_A);
diffusion_x_media_B = gradient(msd_x*a^2) / (2*t_media_B);
plot(n_stepsa, diffusion_x_media_A);
hold on;
plot(n_stepsa, diffusion_x_media_B);
hold off;
title('Diffusion Constant in the x-Direction for Media A and B');
xlabel('Time (mm)');
ylabel('Diffusion Constant');
legend({'Media A', 'Media B'});
xlim([0, num_steps*a]);