The current code does not reflect the particles off the bottle-neck boundaries. It only plots the boxes and checks if particles are inside, but it does not update their trajectories when they hit the boxes.To make particles reflect off the bottle-neck boundaries, you need to:
- At each simulation step, check if a particle will enter a box
- If a particle hits a box, reverse its velocity component perpendicular to the wall to reflect it.
LR = 200e-9; % Length of region
BR = 100e-9; % Breadth of region
numParticles = 1000;
numSteps = 50; % Increase for visible trajectories
dt = 1e-14; % Time step
% Particle properties
x = rand(numParticles,1) * LR;
y = rand(numParticles,1) * BR;
theta = 2*pi*rand(numParticles,1);
v = 1e5; % m/s, arbitrary constant speed
vx = v * cos(theta);
vy = v * sin(theta);
% Obstacles (bottle-neck rectangles)
boxes = [0.8e-7 0.0e-7 0.4e-7 0.4e-7; % [x y width height]
0.8e-7 0.6e-7 0.4e-7 0.4e-7];
% For plotting trajectories
x_traj = zeros(numParticles, numSteps);
y_traj = zeros(numParticles, numSteps);
for t = 1:numSteps
% Save trajectory
x_traj(:,t) = x;
y_traj(:,t) = y;
% Update positions
x_new = x + vx*dt;
y_new = y + vy*dt;
% Reflect off region boundaries
outX = (x_new<0) | (x_new>LR);
outY = (y_new<0) | (y_new>BR);
vx(outX) = -vx(outX);
vy(outY) = -vy(outY);
% Update positions after boundary reflection
x_new(outX) = x(outX) + vx(outX)*dt;
y_new(outY) = y(outY) + vy(outY)*dt;
% Reflect off boxes
for b = 1:size(boxes,1)
bx = boxes(b,1); by = boxes(b,2); bw = boxes(b,3); bh = boxes(b,4);
inBoxX = (x_new > bx) & (x_new < bx+bw);
inBoxY = (y_new > by) & (y_new < by+bh);
inBox = inBoxX & inBoxY;
% Find which wall was crossed
crossedLeft = (x <= bx) & (x_new > bx) & inBoxY;
crossedRight = (x >= bx+bw)& (x_new < bx+bw)& inBoxY;
crossedBottom = (y <= by) & (y_new > by) & inBoxX;
crossedTop = (y >= by+bh)& (y_new < by+bh)& inBoxX;
% Reflect velocities
vx(crossedLeft | crossedRight) = -vx(crossedLeft | crossedRight);
vy(crossedBottom | crossedTop) = -vy(crossedBottom | crossedTop);
% Update positions after box reflection
x_new(crossedLeft | crossedRight) = x(crossedLeft | crossedRight) + vx(crossedLeft | crossedRight)*dt;
y_new(crossedBottom | crossedTop) = y(crossedBottom | crossedTop) + vy(crossedBottom | crossedTop)*dt;
end
% Update positions for next step
x = x_new;
y = y_new;
end
% Plot
figure;
hold on;
for i = 1:numParticles
plot(x_traj(i,:), y_traj(i,:));
end
for b = 1:size(boxes,1)
rectangle('Position', boxes(b,:), 'FaceColor', [1 1 1], 'EdgeColor', 'k', 'LineWidth', 2);
end
The Figure looks like this:
Hope this helps!
