event-driven simulation

1.  Goal: Simulate the motion of N moving particles that behave according to the laws of elastic collision.

     Hard disc model:

         - Moving particles interact via elastic collisions with each other and walls.

         - Each particle is a disc with known position, velocity, mass, and radius.

         - No other forces.

 

2.  N bouncing balls in the unit square without considering balls colliding with each other

public class Ball
{
    private double rx, ry; // position
    private double vx, vy; // velocity
    private final double radius; // radius
    public Ball(...)
    { /* initialize position and velocity */ }
    public void move(double dt)
    {
        if ((rx + vx*dt < radius) || (rx + vx*dt > 1.0 - radius)) {
            rx = 1.0 - radius - vx * dt +  (1.0 - radius - rx);
            vx = -vx; 
        } else {
            rx = rx + vx * dt;    
        }
        if ((ry + vy*dt < radius) || (ry + vy*dt > 1.0 - radius)) {
            ry = 1.0 - radius - vy * dt +  (1.0 - radius - ry);
            vy = -vy; 
       } else {
            ry = ry + vy*dt; 
       }
    }

    public void draw()
    {
        StdDraw.filledCircle(rx, ry, radius); 
    }
}

public class BouncingBalls
{
    public static void main(String[] args)
    {
        int N = Integer.parseInt(args[0]);
        Ball[] balls = new Ball[N];
        for (int i = 0; i < N; i++)
            balls[i] = new Ball();
        while(true)
        {
            StdDraw.clear();
            for (int i = 0; i < N; i++)
            {
                balls[i].move(0.5);
                balls[i].draw();
            }
            StdDraw.show(50);
        }
    }
}

 

3.  Naive Implementation for handling collision:

    -- Discretize time in quanta of size dt.

    -- Update the position of each particle after every dt units of time, and check for overlaps.

    -- If overlap, roll back the clock to the time of the collision, update the velocities of the colliding particles, and continue the simulation.

 

    Main drawbacks:

    -- ~ N^2 / 2 overlap checks per time quantum.

    -- Simulation is too slow if dt is very small.

    -- May miss collisions if dt is too large. (if colliding particles fail to overlap when we are looking)

 

 

4.  Event-Driven Simulation : 

    Change state only when something happens.

    -- Between collisions, particles move in straight-line trajectories.

    -- Focus only on times when collisions occur.

    -- Maintain PQ of collision events, prioritized by time.

    -- Remove the min = get next collision.

   Collision prediction.

   -- Particle i: radius si, position (rxi, ryi), velocity (vxi, vyi).

   -- Particle j: radius sj, position (rxj, ryj), velocity (vxj, vyj).

   -- Will particles i and j collide? If so, when?

    

 

    Collision resolution.

    -- If collision occurs, update colliding particle(s) according to laws of elastic collisions.

 

 

5.  Collision system :

    Initialization:

    -- Fill PQ with all potential particle-wall collisions.

    -- Fill PQ with all potential particle-particle collisions.

    Main loop:

    -- Delete the impending event from PQ (min priority = t).

    -- If the event has been invalidated, ignore it.

    -- Advance all particles to time t, on a straight-line trajectory.

    -- Update the velocities of the colliding particle(s).

    -- Predict future particle-wall and particle-particle collisions involving the colliding particle(s) and insert events onto PQ.

/**
Conventions.
-- Neither particle null ⇒ particle-particle collision.
-- One particle null ⇒ particle-wall collision.
-- Both particles null ⇒ redraw event.
*/
private class Event implements Comparable<Event>
{
    private double time; // time of event
    private Particle a, b; // particles involved in event
    private int countA, countB; // collision counts for a and b, used to check whether the collsion is still valid
    public Event(double t, Particle a, Particle b) { }
    public int compareTo(Event that)
   {return this.time - that.time; }
    public boolean isValid() { }
}

 

public class Particle
{
    private double rx, ry; // position
    private double vx, vy; // velocity
    private final double radius; // radius
    private final double mass; // mass
    private int count; // number of collisions
    public Particle(...) { }
    public void move(double dt) { }
    public void draw() { }
    public double timeToHit(Particle that) {
        if (this == that) return INFINITY;
        double dx = that.rx - this.rx, dy = that.ry - this.ry;
        double dvx = that.vx - this.vx; dvy = that.vy - this.vy;
        double dvdr = dx*dvx + dy*dvy;
        if( dvdr > 0) return INFINITY;
        double dvdv = dvx*dvx + dvy*dvy;
        double drdr = dx*dx + dy*dy;
        double sigma = this.radius + that.radius;
        double d = (dvdr*dvdr) - dvdv * (drdr - sigma*sigma);
        if (d < 0) return INFINITY;
        return -(dvdr + Math.sqrt(d)) / dvdv;
    }
    public double timeToHitVerticalWall() { }
    public double timeToHitHorizontalWall() { }
    public void bounceOff(Particle that) {
        double dx = that.rx - this.rx, dy = that.ry - this.ry;
        double dvx = that.vx - this.vx, dvy = that.vy - this.vy;
        double dvdr = dx*dvx + dy*dvy;
        double dist = this.radius + that.radius;
        double J = 2 * this.mass * that.mass * dvdr / ((this.mass + that.mass) * dist);
        double Jx = J * dx / dist;
        double Jy = J * dy / dist;
        this.vx += Jx / this.mass;
        this.vy += Jy / this.mass;
        that.vx -= Jx / that.mass;
        that.vy -= Jy / that.mass;
        this.count++;
        that.count++;
    }
    public void bounceOffVerticalWall() { }
    public void bounceOffHorizontalWall() { }
}

 

public class CollisionSystem
{
    private MinPQ<Event> pq; // the priority queue
    private double t = 0.0; // simulation clock time
    private Particle[] particles; // the array of particles
    public CollisionSystem(Particle[] particles) { }
    private void predict(Particle a)
    {
        if (a == null) return;
        for (int i = 0; i < N; i++)
        {
            double dt = a.timeToHit(particles[i]);
            pq.insert(new Event(t + dt, a, particles[i]));
        }
        pq.insert(new Event(t + a.timeToHitVerticalWall() , a, null));
        pq.insert(new Event(t + a.timeToHitHorizontalWall(), null, a));
    }

    private void redraw() { }
    public void simulate() {
        pq = new MinPQ<Event>();
        for(int i = 0; i < N; i++) predict(particles[i]);
        pq.insert(new Event(0, null, null));
        while(!pq.isEmpty())
        {
            Event event = pq.delMin();
            if(!event.isValid()) continue;
            Particle a = event.a;
            Particle b = event.b;
            for(int i = 0; i < N; i++)
                particles[i].move(event.time - t);
            t = event.time;
            if (a != null && b != null) a.bounceOff(b);
            else if (a != null && b == null) a.bounceOffVerticalWall();
            else if (a == null && b != null) b.bounceOffHorizontalWall();
            else if (a == null && b == null) redraw();
            predict(a);
            predict(b);
        }
    }
}