The Keplerian billiard

This example shows how it possible to use heyoka.py’s event detection system to handle collisions. The goal is to simulate the motion of a test particle around a massive central object, with the following twists:

• the central object is a sphere against which the test particle will bounce in a perfectly elastic fashion (i.e., no loss of kinetic energy),

• the simulation domain is enclosed in a square box against whose walls the test particle will also elastically bounce.

In other words, we are dealing with a Keplerian billiard.

Let us begin by setting up the dynamics of a one-body problem. We will be using the model.nbody() helper for this purpose:

import heyoka as hy

sys = hy.model.nbody(n = 2, masses = [1., 0.])

Note how we created a 2-body problem with the second mass (the test particle) set to zero. This is not the most efficient way of creating a central-force problem in heyoka.py, but it will suffice for the purpose of this example.

Next, we create the symbolic variables representing the \(\left( x, y \right)\) coordinates of the central object and the test particle. Note that these symbolic variables have already been introduced by the model.nbody() helper, we need to re-create them for use in the event handling functions below:

x0, y0, x1, y1 = hy.make_vars("x_0", "y_0", "x_1", "y_1")

Next, we define a couple of constants representing, respectively, the radius of the central object and the size of the square box:

p_radius = 0.15
box_size = 2.

We are now ready to implement the bouncing behaviour. In order to keep things simple, we will limit ourselves to a 2D formulation of the problem, where we will need to detect and react to the following events:

• the test particle hits one of the 4 walls of the box: if it hits the top/bottom walls, then the \(y\) component of the velocity vector needs to be inverted, if it hits the left/right walls, then the \(x\) component of the velocity needs to be inverted;

• the test particle hits the central object: the component of the velocity vector perpendicular to the central object’s surface at the collision point needs to be inverted.

Because the events cause discontinuous changes in the dynamical state of the system, they need to be implemented as terminal events in heyoka.py’s event detection system.

Let us first proceed to the definitions of the callbacks. These are the functions that will be called when the events trigger. We begin with the callbacks that will be called when the test particle hits the box’s walls:

def cb_left_right(ta, d_sgn):
    ta.state[9] = -ta.state[9]
    return True

def cb_top_bottom(ta, d_sgn):
    ta.state[10] = -ta.state[10]
    return True

The callbacks take as first input argument a reference to the integrator object (don’t worry about the second argument for the time being). When invoked, they will invert, respectively, the \(x\) and \(y\) components of the velocity vector of the test particle, which are stored at indices 9 and 10 in the state vector of the system. Both callbacks return True to indicate that the integration should continue after their invocation.

The callback implementing the bouncing behaviour against the central object is slightly more complicated. Let’s take a look:

import numpy as np

def cb_center(ta, d_sgn):
    # Fetch the xy coordinate vector
    # of the test particle from
    # the state vector.
    xy_pos = ta.state[6:8]
    
    # Construct the unit vector
    # of xy_pos.
    xy_pos_uvec = xy_pos / np.linalg.norm(xy_pos)
    
    # Fetch the xy velocity vector of the
    # test particle from the state
    # vector.
    xy_vel = ta.state[9:11]
    
    # Project xy_vel onto the
    # unit vector.
    vproj = np.dot(xy_vel, xy_pos_uvec)
    
    # Compute the velocity change.
    Dv = -vproj*xy_pos_uvec
    
    # Add the velocity change to
    # the state vector.
    xy_vel += 2*Dv
    
    return True

Because the central object is sitting in the origin and it has a spherical shape, the effect of the test particle colliding onto it is always a flip of the component of the velocity vector across the position vector of the contact point.

We can now proceed to the construction of the terminal event objects:

# Events for bouncing off the walls.
ev_bottom = hy.t_event(y1 + box_size/2, callback = cb_top_bottom, direction=hy.event_direction.negative)
ev_top = hy.t_event(y1 - box_size/2, callback = cb_top_bottom, direction=hy.event_direction.positive)
ev_left = hy.t_event(x1 + box_size/2, callback = cb_left_right, direction=hy.event_direction.negative)
ev_right = hy.t_event(x1 - box_size/2, callback = cb_left_right, direction=hy.event_direction.positive)

# Event for bouncing off the central object.
ev_center = hy.t_event((x1 - x0)**2+(y1 - y0)**2 - (p_radius)**2, callback = cb_center, direction=hy.event_direction.negative)

The first argument passed to the constructor of an event is the event equation. The wall-bouncing event equations are:

\[\begin{split} \begin{cases} y_1 \pm \frac{\mathrm{box\_size}}{2} = 0 \\ x_1 \pm \frac{\mathrm{box\_size}}{2} = 0 \end{cases}. \end{split}\]

These equations will be satisfied when the particle hits one of the walls. The center-bouncing event equation is:

\[ \left(x_1 - x_0 \right)^2 + \left(y_1 - y_0 \right)^2 - \mathrm{p\_radius}^2 = 0. \]

This equation will be satisfied when the particle hits the spherical surface of the central object.

Note how we specified an event direction (positive or negative) for all events. These directions indicate that the collisions must be detected only when the particle hits the walls “from the inside” and the central sphere “from the outside”. For this simple example, specifying the direction of collisions is not strictly necessary, however in more numerically-sensitive examples a double-detection of a collision event may take place after the execution of the callback, which would lead to the test particle exiting the box or penetrating the central sphere. Specifying a direction for collision events avoids this issue because it filters out the second spurious bounce.

We are finally ready to construct the integrator object:

ta = hy.taylor_adaptive(sys, [0.,0.,0.,0.,0.,0.,0.2,0.,0.,0.,3.,0.],
                        t_events=[ev_bottom, ev_top, ev_left, ev_right, ev_center])

The vector of initial conditions \(\left[ 0,0,0,0,0,0,0.2,0,0,0,3,0 \right]\) places the central object in the origin with zero velocity, and the test particle on a planar elliptic orbit around it.

Let us now generate a time grid and integrate over it:

grid = np.linspace(0, 15, 10000)
out = ta.propagate_grid(grid)

We can now produce a fancy animation of the Keplerian billiard:

%%capture

import matplotlib.pyplot as plt
from matplotlib import animation, rc
from IPython.display import HTML

fig = plt.figure(figsize = (12, 12))

fig, ax = plt.subplots()

tc = plt.Circle((0.2 , 0.), 0.01, ec='black', fc='black', zorder=3)
ax.add_artist(tc)
line, = ax.plot([], [], 'k--')

def init():
    ax.add_patch(plt.Rectangle((-1, -1), 2, 2, edgecolor='k', facecolor="none", linestyle = '--'))
    cc = plt.Circle((0. , 0.), p_radius, ec='black', fc='orange', zorder=2)
    ax.add_artist(cc)
    ax.axis('equal')
    ax.set_xlim((-1.1,1.1))
    ax.set_ylim((-1.1,1.1))
    return (tc,)

def animate(i):
    tc.set_center((out[-1][i*40,6], out[-1][i*40,7]))
    line.set_data(out[-1][:i*40,6], out[-1][:i*40,7])
    return tc, line


anim = HTML(animation.FuncAnimation(fig, animate, init_func=init,
                                    frames=250, interval=50, 
                                    blit=True).to_jshtml());

anim

Once Loop Reflect