I'm having trouble simulating multiple balls bouncing in a box.

I've created the following code to simulate a ball bouncing around a box with a rectangular obstacle in the middle.

function ballbox(x, y, vx, vy, rad, col)

%x and y are the initial positions of the ball
%vx is the initial velocity in the x-direction
%vy is the initial velocity in the y-direction
%rad and col define the radius and colour of the ball, respectively.

obswallcond  %defines the positions of the wall and the obstacle

x(1) = x;
y(1) = y;
vx(1) = vx;
vy(1) = vy;

for i = 1:Nt

      x(i+1) = x(i) + vx(i)*dt;
      y(i+1) = y(i) + vy(i)*dt;

      bouncecond %Script that defines when the ball collides with the walls or obstacle

      if corners || obscorners
          vx(i+1) = -vx(i);
          vy(i+1) = -vy(i);

      elseif invertvx
          vx(i+1) = -vx(i);
          vy(i+1) =  vy(i);

      elseif invertvy
          vx(i+1) =  vx(i);
          vy(i+1) = -vy(i);

      else
          vx(i+1) = vx(i);
          vy(i+1) = vy(i);
      end

      scatter(x(i+1), y(i+1), rad, col ,'filled')
      hold on;
      obswallplot
      hold off;
      drawnow
  end
  end

I have the 'corners', 'invertvy' and 'invertvx' commands defined in the 'bouncecond' script. This function works when I only have one ball but I'm unable to simulate more than one ball, simultaneously.

Whenever I try something like,

ballbox(90, -90, 10, 20, 20, 'r');
hold on;
ballbox(0, -80, -10, 20, 20, 'b');
hold off;

the function will bring up the second ball after the simulation for the first one has finished and not simultaneously. I'm only trying to keep the simulation very simple at the minute so I'm not worried about the balls interacting with each other at the minute. Any help would be greatly appreciated.

If it helps, my 'obswallcond' script is

L1 = 100;
L2 = 100;

l1 = 40;
l2 = 40;

Nt = 1000;
dt = 0.1;

'bouncecond' is

hitrwall = x(i+1) ==  L1;
hitlwall = x(i+1) == -L1;
hittwall = y(i+1) ==  L2;
hitbwall = y(i+1) == -L2;

hitrobs = (y(i+1) >= -l2 && y(i+1) <= l2 && x(i+1) ==  l1); 
hitlobs = (y(i+1) >= -l2 && y(i+1) <= l2 && x(i+1) == -l1);
hittobs = (x(i+1) >= -l1 && x(i+1) <= l1 && y(i+1) ==  l2);
hitbobs = (x(i+1) >= -l1 && x(i+1) <= l1 && y(i+1) == -l2);

corners = (hitrwall && hittwall) || (hittwall && hitlwall) || (hitlwall && hitbwall) || (hitbwall && hitrwall);
obscorners = (hitrobs && hittobs) || (hittobs && hitlobs) || (hitlobs && hitbobs) || (hitbobs && hitrobs);
invertvx = hitrwall || hitlwall || hitrobs || hitlobs;
invertvy = hittwall || hitbwall || hittobs || hitbobs;

and 'obswallplot' is

      line([-l1 -l1], [-l2 l2] ,'color', 'k')
      line([-l1 l1],  [l2 l2]  ,'color','k')
      line([l1 l1] ,  [l2 -l2] ,'color','k')
      line([l1 -l1],  [-l2 -l2],'color','k')
      xlim([-L1 L1])
      ylim([-L2 L2])

Answer 1

[Walter Roberson]

1) You need to fix your bouncecond to not test for equality.

For example,

hitrwall = x(i+1) ==  L1;

hitrobs = (y(i+1) >= -l2 && y(i+1) <= l2 && x(i+1) ==  l1);

is wrong because it checks for x(i+1) exactly equal to l1, but x has a non-integer velocity and might be projected to go beyond the limit.

hitrwall = x(i+1) >=  L1;

hitrobs = (y(i+1) >= -l2 && y(i+1) <= l2 && x(i+1) >=  l1);

2) You need to comment what the difference is between L1 and l1 (and similar). For example why isn't hitrobs defined in terms of

(y(i+1) >= -l2 && y(i+1) <= l2 && hitrwall);

3) You need to vectorize ballbox()

For example,

x(1,:) = x;
y(1, :) = y;
vx(1, :) = vx;
vy(1, :) = vy;

for i = 1:Nt

        x(i+1, :) = x(i, :) + vx(i, :)*dt;
        y(i+1, :) = y(i, :) + vy(i, :)*dt;

4) You need to vectorize bouncecond

For example,

hitrobs = (y(i+1, :) >= -l2 & y(i+1, :) <= l2 & x(i+1, :) ==  l1);

except fixed to take into account (1) and (2)

5) As both of those involve testing conditions, you will either need to loop over all of the columns or you will need to use logical indexing. For example,

      if corners || obscorners
          vx(i+1) = -vx(i);
          vy(i+1) = -vy(i);

      elseif invertvx
          vx(i+1) = -vx(i);
          vy(i+1) =  vy(i);

should be more like

      mask = corners | obscorners
      vx(i+1, mask) = -vx(i, mask);
      vy(i+1, mask) = -vy(i, mask);

      mask = invertvx
      vx(i+1, mask) = -vx(i, mask);
      vy(i+1, mask) =  vy(i, mask);

[AllenK]

Ah yeah, sorry. L1 and L2 define the boundaries of the box, i.e the right wall is defined as x = L1 and the base is y = -L2, etc. The inner obstacle is then defined by l1 and l2, where the right of wall of the obstacle is x = l1 as long as -l2 <= y <= l2.

I'm not quite sure why you are setting 'mask = corners | obscorners' and following your advice in step 5 doesn't seem to work.

[Walter Roberson]

You detect hitting the wall and hitting internal obstacles, but you do not adjust the x and y position to be on the proper side of the obstacle. Your total distance traveled in a step should be defined by the velocity and the length of the time step, but part of the time was spent traveling towards the obstacle and part of the time was spent traveling away, so you need to figure out when the bounce occurred and what coordinate the bounce was at, in order to figure out how far you traveled after the bounce until the end of the time step.

Answer 2

[Image Analyst]

ballbox() calls hold off, which, when you're trying to animate your second ball, blows away the plotting from your first ball.

[AllenK]

If I don't have hold off then a trail follows the ball.

[Image Analyst]

You do this:

ballbox(90, -90, 10, 20, 20, 'r'); % Do ball 1
hold on;
ballbox(0, -80, -10, 20, 20, 'b'); % Do ball 2

The first one simulates ball 1. Only after it's done and the it leaves the ballbox function, does the simulation for ball 2 begin. This will not make them animate simultaneously.

[AllenK]

OK, but how would you make them animate simultaneously?

[Walter Roberson]

You need to vectorize like I described.

Or you need to remove the loop from the routine so that it only does one step, and then keep calling it step by step