This example is given by a niubi person in stackoverflow.com.
Here is his code:
function bounce(npts, vmin, vmax, radius, center) % Initial direction/velocity of the points direction = rand(npts, 1) * 2 *pi; velocity = (rand(npts, 1) * (vmax - vmin)) + vmin; % Create random starting locations within the circle theta = rand(npts, 1) * 2*pi; r = radius * sqrt(rand(npts, 1)); XY = [r .* cos(theta(:)) + center(1), ... r .* sin(theta(:)) + center(2)]; % Initial plot objects hfig = figure('Color', 'w'); hax = axes('Parent', hfig); % Plot the dots as black markers hdots = plot(XY(:,1), XY(:,2), ... 'Parent', hax, ... 'Marker', '.', ... 'Color', 'k', ... 'LineStyle', 'none', ... 'MarkerSize', 12); hold(hax, 'on') axis(hax, 'equal') % Plot the circle as a reference t = linspace(0, 2*pi, 100); plot(radius * cos(t) + center(1), ... radius * sin(t) + center(2)) % Keep simulating until we actually close the window while ishghandle(hfig); % Determine new dot locations [XY, direction] = step(XY, direction, velocity, radius, center); % Update the dot plot to reflect new locations set(hdots, 'XData', XY(:,1), 'YData', XY(:,2)) % Force a redraw drawnow end end function [XYnew, direction] = step(XY, direction, velocity, radius, center) % Compute the next position of the points DX = [cos(direction(:)) .* velocity, ... sin(direction(:)) .* velocity]; XYnew = XY + DX; % Now check that they are all inside circle isOutside = sum(bsxfun(@minus, XYnew, center).^2, 2) > radius^2; % The ones that are outside should "bounce" back into the circle if any(isOutside) orig = XY(isOutside,:); new = XYnew(isOutside,:); delta = -DX(isOutside,:); % Find intersection of this path with the circle % Taken from: http://math.stackexchange.com/a/311956 a = sum(delta.^2, 2); b = sum(2 .* delta .* bsxfun(@minus, orig, center), 2); c = sum(bsxfun(@minus, orig, center).^2, 2) - radius^2; t = (2 * c) ./ (-b + sqrt(b.^2 - 4 .* a .* c)); xintersect = orig(:,1) + delta(:,1) .* t; yintersect = orig(:,2) + delta(:,2) .* t; % Get tangent at this intersection (slope/intercept form) m = - 1 ./ ((yintersect - center(2)) ./ (xintersect - center(1))); b = yintersect - m .* xintersect; % "Reflect" outside points across the tangent line to "bounce" them % Equations from: http://stackoverflow.com/a/3307181/670206 d = (new(:,1) + (new(:,2) - b) .* m) ./ (1 + m.^2); XYnew(isOutside,1) = 2 * d - new(:,1); XYnew(isOutside,2) = 2 .* d .* m - new(:,2) + 2 .* b; % Recompute the direction of the particles that "bounced" direction(isOutside) = atan2(XYnew(isOutside,2) - yintersect, ... XYnew(isOutside,1) - xintersect); end end
After that, by running the following command you are able to obtain the following result:
bounce(100, 0.01, 0.2, 5, [0 0]);