How to make the balls in this program organize and settle at the top of the box instead of leaving it

Question

zina shadidi 2018-3-16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/388682-how-to-make-the-balls-in-this-program-organize-and-settle-at-the-top-of-the-box-instead-of-leaving

评论： zina shadidi 2018-3-21

function y = atoms(x)
clc
n = 25;  % no. of atoms
dt = 0.1;   %  time
ma=1;
mb=1;
x0=[0 0 0 0];
l=200;
p = rand(n,2)*0.5*l;
v = rand(n,2)+.5;
r=0.5;
for tt=1:5000
      plot(p(:,1),p(:,2),'O', 'MarkerFaceColor', 'b');
      axis([0  l  0  l])
      drawnow;
      for i = 1 : n
          if p(i, 1) < 0
              v(i, 1) =- v(i, 1);
          end
          if p(i, 1) > l
              v(i, 1) = -v(i, 1);
          end
          if p(i, 2) < 0
              v(i, 2) = v(i, 2);
          end
      end
      for i = 1 : n
          for j = i + 1 : n
              d = sqrt((p(i, 1) - p(j, 1)) ^ 2 + (p(i, 2) - p(j, 2)) ^ 2);
              if d < r
                  va = sqrt(v(i,1)^2 + v(i,2)^2);
                  vb = sqrt(v(j,1)^2 + v(j,2)^2);
                  gama = atan((p(i,2)-p(j,2))/(p(i,1)-p(j,1)));
                  alpha = atan2(v(i,1),v(i,2)) - gama;
                  beta  = atan2(v(j,1),v(j,2)) - gama;
                  [x,fval ] = fsolve(@gh,x0);
                  v(i,1)=x(1)*cos(x(3)+gama);
                  v(i,2)=x(1)*sin(x(3)+gama);
                  v(j,1)=x(2)*cos(x(4)+gama);
                  v(j,2)=x(2)*sin(x(4)+gama);
      p(i,1)=p(i,1)+v(i,1)*dt;
      p(i,2)=p(i,2)+v(i,2)*dt;
      p(j,1)=p(j,1)+v(j,1)*dt;
      p(j,2)=p(j,2)+v(j,2)*dt;
      p(i,1)=p(i,1)+v(i,1)*dt;
      p(i,2)=p(i,2)+v(i,2)*dt;
      p(j,1)=p(j,1)+v(j,1)*dt;
      p(j,2)=p(j,2)+v(j,2)*dt;
      p(i,1)=p(i,1)+v(i,1)*dt;
      p(i,2)=p(i,2)+v(i,2)*dt;
      p(j,1)=p(j,1)+v(j,1)*dt;
      p(j,2)=p(j,2)+v(j,2)*dt;
      p(i,1)=p(i,1)+v(i,1)*dt;
      p(i,2)=p(i,2)+v(i,2)*dt;
      p(j,1)=p(j,1)+v(j,1)*dt;
      p(j,2)=p(j,2)+v(j,2)*dt;
      p(i,1)=p(i,1)+v(i,1)*dt;
      p(i,2)=p(i,2)+v(i,2)*dt;
      p(j,1)=p(j,1)+v(j,1)*dt;
      p(j,2)=p(j,2)+v(j,2)*dt;
      p(i,1)=p(i,1)+v(i,1)*dt;
      p(i,2)=p(i,2)+v(i,2)*dt;
      p(j,1)=p(j,1)+v(j,1)*dt;
      p(j,2)=p(j,2)+v(j,2)*dt;
              end
          end
      end
      p=p+v*dt;
end
      function y = gh(x)
          e = 1;
          y = [ma*va*cos(alpha) + mb*vb*cos(beta) - ma*x(1)*cos(x(3)) - mb*x(2)*cos(x(4));
               va*sin(alpha) - x(1)*sin(x(3));
               vb*sin(beta)  - x(2)*sin(x(4));
               x(1)*cos(x(3)) - x(2)*cos(x(4)) - e*(va*cos(alpha) - vb*cos(beta))];
      end
end

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

KSSV 2018-3-16

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/388682-how-to-make-the-balls-in-this-program-organize-and-settle-at-the-top-of-the-box-instead-of-leaving#answer_310290

编辑：KSSV 2018-3-16

在 MATLAB Online 中打开

    function y = atoms()    
    n = 25;  % no. of atoms
    dt = 0.1;   %  time
    ma=1;
    mb=1;
    x0=[0 0 0 0];
    l=200;
    p = rand(n,2)*0.5*l;
    v = rand(n,2)+.5;
    r=0.5;
    for tt=1:5000
        idx = abs(p(:,2)-l)<=1 ;  % added this line         
        p(idx,2) = l ;            % added this line 
        plot(p(:,1),p(:,2),'O', 'MarkerFaceColor', 'b');
        axis([0  l  0  l])
        drawnow;
        for i = 1 : n
            if p(i, 1) < 0
                v(i, 1) =- v(i, 1);
            end
            if p(i, 1) > l
                v(i, 1) = -v(i, 1);
            end
            if p(i, 2) < 0
                v(i, 2) = v(i, 2);
            end
        end
        for i = 1 : n
            for j = i + 1 : n
                d = sqrt((p(i, 1) - p(j, 1)) ^ 2 + (p(i, 2) - p(j, 2)) ^ 2);
                if d < r
                    va = sqrt(v(i,1)^2 + v(i,2)^2);
                    vb = sqrt(v(j,1)^2 + v(j,2)^2);
                    gama = atan((p(i,2)-p(j,2))/(p(i,1)-p(j,1)));
                    alpha = atan2(v(i,1),v(i,2)) - gama;
                    beta  = atan2(v(j,1),v(j,2)) - gama;
                    [x,fval ] = fsolve(@gh,x0);
                    v(i,1)=x(1)*cos(x(3)+gama);
                    v(i,2)=x(1)*sin(x(3)+gama);
                    v(j,1)=x(2)*cos(x(4)+gama);
                    v(j,2)=x(2)*sin(x(4)+gama);
                    p(i,1)=p(i,1)+v(i,1)*dt;
                    p(i,2)=p(i,2)+v(i,2)*dt;
                    p(j,1)=p(j,1)+v(j,1)*dt;
                    p(j,2)=p(j,2)+v(j,2)*dt;
                    p(i,1)=p(i,1)+v(i,1)*dt;
                    p(i,2)=p(i,2)+v(i,2)*dt;
                    p(j,1)=p(j,1)+v(j,1)*dt;
                    p(j,2)=p(j,2)+v(j,2)*dt;
                    p(i,1)=p(i,1)+v(i,1)*dt;
                    p(i,2)=p(i,2)+v(i,2)*dt;
                    p(j,1)=p(j,1)+v(j,1)*dt;
                    p(j,2)=p(j,2)+v(j,2)*dt;
                    p(i,1)=p(i,1)+v(i,1)*dt;
                    p(i,2)=p(i,2)+v(i,2)*dt;
                    p(j,1)=p(j,1)+v(j,1)*dt;
                    p(j,2)=p(j,2)+v(j,2)*dt;
                    p(i,1)=p(i,1)+v(i,1)*dt;
                    p(i,2)=p(i,2)+v(i,2)*dt;
                    p(j,1)=p(j,1)+v(j,1)*dt;
                    p(j,2)=p(j,2)+v(j,2)*dt;
                    p(i,1)=p(i,1)+v(i,1)*dt;
                    p(i,2)=p(i,2)+v(i,2)*dt;
                    p(j,1)=p(j,1)+v(j,1)*dt;
                    p(j,2)=p(j,2)+v(j,2)*dt;
                end
            end
        end
        p=p+v*dt;
    end
        function y = gh(x)
            e = 1;
            y = [ma*va*cos(alpha) + mb*vb*cos(beta) - ma*x(1)*cos(x(3)) - mb*x(2)*cos(x(4));
                va*sin(alpha) - x(1)*sin(x(3));
                vb*sin(beta)  - x(2)*sin(x(4));
                x(1)*cos(x(3)) - x(2)*cos(x(4)) - e*(va*cos(alpha) - vb*cos(beta))];
        end
    end