使粒子在圆圈内随机移动

Question

我试图使绘制点在已知半径和中心的圆内移动.目前我能够在边界内生成点,但现在我需要让它们移动.

我有以下脚本来生成粒子的初始位置.

function [particle_gen] = generate(n,centre,radius)
    %generates n particles in a circle with specified centre and radius

    angle = 2 * pi * rand(n,1);
    r = radius * sqrt(rand(n,1));
    X = r.*cos(angle) + centre(1);
    Y = r.*sin(angle) + centre(2);
    plot(X,Y,'.k')
end

我想为它们设置动画,使粒子以恒定的速度沿直线行进,直到它们撞到圆形边界并反弹.我需要这一切都发生在同一个情节中.

Answer 1

好的,这个问题有几个方面我将分别解决.

更新绘图对象

至于创建绘图和更新绘图值的后勤,您可以使用MATLAB的句柄图形完成所有这些操作.如果创建绘图对象,则可以使用返回的图形来更新绘图.

plt = plot(x, y, 'k.');
set(plt, 'XData', newx, 'YData', newy)

在您的情况下,newx和newy值将在循环中计算,然后您可以分配它们.通过这样做,您将不会经常创建新的图形对象,使用MATLAB可能会非常昂贵.

由于您正在进行大量计算,因此您还希望drawnow每次都通过循环使用MATLAB来渲染图形.

粒子轨迹

您将需要计算要用于更新绘图的值xnew和ynew值.您可以通过多种方式确定这些轨迹,但我在此提出了一种方法.

1. 应用置换.使用速度和方向,您可以将位移应用于每个粒子.
2. 确定点是否在圆圈之外.如果位移将点放在圆外(这里是基本方程),那么我们需要将它放回圆中.如果圆圈外没有点,则返回步骤1.
3. 对于圆外的每个点,找到其路径与圆的交点.使用此解决方案作为参考,这是一项非常简单的练习.
4. 此时找到圆的切线.同样,这是微不足道的,因为我们可以从圆的中心到圆上的交点画一条线.该线斜率的负倒数是交点处曲线的切线斜率.通过将其与交叉点相结合,我们得到了一条线的等式.
5. 在该切线上反射圆外的点,将其"反弹"回圆圈.可以使用此处提供的等式来执行线上点的这种反射.

示例解决方案

我已经实现了上述步骤.虽然这可能不是计算效率最高的方法,但它提供了所需的结果.

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: https://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: https://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

结果

通过运行以下命令,我能够获得以下结果.

bounce(100, 0.01, 0.2, 5, [0 0]);