I think everyone knows what an ellipse is. It’s an elongated, squished circle that looks kind of like an egg. In mathematical terms, it’s considered a conic section. That means you can slice through a cone and trace an ellipse, like the image below.

Before getting into ellipses, I want to touch on their conic section cousin: the parabola. Parabolas have similar properties and a familiar real-world application.

These big satellite dishes (picture below) were common in the 1980s. They got smaller as technology improved and, in the 90s, shrank to the 2-foot DirecTV-style dishes we still see today.

Technically they are 3-dimensional paraboloids but the same principles apply. Parallel rays are deflected toward a single point, called the focus. Satellite dish manufacturers mount a feedhorn at this point. It collects the amplified signal, filters it, and converts it into an electrical signal. A wire carries that signal to your living room where a satellite receiver decodes it and draws a picture on your TV. Simple!

Notice how all the photons—the little bits of light that carry the signal—hit the focus at the same time. No matter where they strike the parabola, they travel the same distance. A dish collects a signal out of the air like a bowl collects rainwater. Parabolic microphones operate using the same principle.

---

1. The ellipse.

That’s great but we’re here to talk about ellipses.

It turns out that ellipses behave in a similar way. Rather than one focus they have two, and any particle fired from a focus will deflect off the inner wall and hit the other focus. It’s called the reflective property of ellipses. Check out the diagram below to see what I mean.

No matter what direction a particle is fired, it always ends up hitting the other focus. It’s a bank shot that goes in every time. You can imagine the orange tangent line sliding around the perimeter of the ellipse. At every position, the two θ angles will be the same.

And that behavior is what we’ll demonstrate in this post. We’ll fire several particles and animate them to prove that they all impact the focus.

At this point there is a fork in the road. We already know what a particle’s behavior will be after it strikes the ellipse wall: it will change direction and travel toward the focus. So we could simply dictate this behavior and force every particle to move in the direction we know it’s supposed to go. But it will be a lot more interesting to model the deflection pictured above, i.e. calculate the tangent and reflect the particle over the normal line, and let it go wherever it wants to go. I think it’s more fun to define basic rules of motion and see how the system behaves.

It will also require more code, but that’s okay. We’ll step through a main loop to generate each frame of the animation. Particles will be instantiated from a Particle class. This object-oriented approach will make it more difficult to talk through the code. But, as usual, you can find the complete script at the bottom of this post.

---

2. The code.

The formula for an ellipse centered at the origin looks like this:

Where A defines the ellipse’s width along the x-axis and B along the y-axis.

Because we’ll refer to it many times later, let’s create an Ellipse class. It’s not a continuous function—there are two y-values for every value of x—so we need some logic to distinguish between “upper” and “lower.” The distance C from the center (0, 0) to the two foci is defined by C² = A² - B².

class Ellipse:

    def __init__(self, a, b):
        self.A = a
        self.B = b
        self.x = arange(-self.A, self.A + 0.0001, 0.0001)
        self.y_lower, self.y_upper = self.get_y_values()
        self.focus_x_left, self.focus_x_right = self.get_focus_values()
        self.focus_y_left, self.focus_y_right = 0, 0

    def get_y_values(self):
        y_lower = -self.B * sqrt(1 - self.x**2 / self.A**2)
        y_upper = self.B * sqrt(1 - self.x**2 / self.A**2)
        return y_lower, y_upper

    def get_focus_values(self):
        focus_x_left = -sqrt(self.A**2 - self.B**2)
        focus_x_right = sqrt(self.A**2 - self.B**2)
        return focus_x_left, focus_x_right


ellipse = Ellipse(a=5, b=3)

To model a particle’s motion, we need to know its direction (slope) and speed. speed will be the distance a particle travels during each frame of the animation.

We’re going to fire several particles in different directions. Let’s package them into a list so we can easily iterate over them. For now, the list will only contain one particle, but more will be appended later. This approach will allow us to fire particles one after another, rather than all at once.

class Particle:

    def __init__(self, initial_slope, speed):
        self.x = ellipse.focus_x_left
        self.y = ellipse.focus_y_left
        self.movement_history = []
        self.slope = initial_slope
        self.collision_point = self.get_collision_point()
        self.speed = speed
        self.speed_x, self.speed_y = self.get_speed_components()
        self.phase = 1
    
    [...]

particle = Particle(initial_slope=-10, speed=0.1)

particle_list = [particle]

Inside the __init__ you’ll notice a couple other methods.

get_collision_point checks where exactly the particle will collide with the ellipse wall. Calculating this value during initialization makes it much easier to do collision detection. We only have to check the particle’s distance from that exact point. When the distance is near zero, a collision has happened.

To solve for this point, we have to set the particle’s position equation equal to the ellipse equation, and find what values of x and y satisfy the expression. This requires an obnoxious quadratic solution and I’ll happily say I used Wolfram Alpha.

Then we plug the x value back into the y=mx+b particle equation to solve for y. Technically, this expression should include a ± symbol. But if we agree to only fire particles to the right, we can ignore it. I think that’s a fair compromise.

class Particle:

    [...]

    def get_collision_point(self):
        y_intercept = self.y - self.slope * self.x
        x_val = (sqrt(ellipse.A**4 * ellipse.B**2 * self.slope**2 + ellipse.A**2 * ellipse.B**4 - ellipse.A**2 * ellipse.B**2 * y_intercept**2) - ellipse.A**2 * y_intercept * self.slope) / (ellipse.A**2 * self.slope**2 + ellipse.B**2)
        y_val = self.slope * x_val + y_intercept
        return (x_val, y_val)

    [...]

The __init__ also references get_speed_components. This method returns x- and y-components of the speed vector. They tell us exactly how the particle will move during each frame.

class Particle:

    [...]

    def get_speed_components(self):
        angle = atan2(self.collision_point[1] - self.y, self.collision_point[0] - self.x)
        speed_x = self.speed * cos(angle)
        speed_y = self.speed * sin(angle)
        return speed_x, speed_y

    [...]

Particle objects have a movement_history attribute which keeps track of all their past locations. We’ll use it to draw a tail on each particle—sort of like the old NHL glow puck but hopefully less hated by viewers.

The phase attribute will assist in collision detection.

• Phase 1 — The particle was recently fired. It’s moving toward the ellipse wall.
• Phase 2 — The particle has recently deflected off the ellipse wall. It’s moving toward the focus.
• Phase 3 — The particle has impacted the focus and its journey is over.

With our ellipse and particle instances in hand, we can jump back to the main code.

Additional particles will be spawned according to the dictionary below. For example, at frame 20 we’ll launch a particle with a slope of -1.

additional_particles = {20: -1, 40: -0.5, 60: -0.1, 80: 0.2, 100: 0.4, 120: 1.2, 140: 7}

Each iteration of this large for loop will generate one frame of the animation.

Begin by creating a figure and an Axes instance on which to draw. Then plot the ellipse and foci.

plt.style.use("wollen_ellipse.mplstyle")

for frame_counter in range(250):

    print(frame_counter)

    fig, ax = plt.subplots()

    ax.plot(ellipse.x, ellipse.y_upper, color="#444", linewidth=2.0, zorder=1)
    ax.plot(ellipse.x, ellipse.y_lower, color="#444", linewidth=2.0, zorder=1)

    ax.scatter([ellipse.focus_x_left, ellipse.focus_x_right],
               [ellipse.focus_y_left, ellipse.focus_y_right],
               color="None", s=80, edgecolor="#444", linewidth=1.8, zorder=1)

    [...]

Iterate over each particle within particle_list.

First update particle position, which I’ll detail in a moment. Then use scatter to plot it. Plotting the 10 most recent values of movement_history will create a tail behind each particle. It sounds strange but it will make it easier to see what direction a particle is traveling.

for frame_counter in range(250):

    [...]

    for particle in particle_list:

        particle.update_position()

        ax.scatter([particle.x], [particle.y], color="#4C4", s=50, edgecolor="#444", linewidth=1.0, zorder=3)

        ax.plot([pos[0] for pos in particle.movement_history[-10:]],
                [pos[1] for pos in particle.movement_history[-10:]],
                color="#4C4", linewidth=1.5, zorder=2)

    [...]

Particle.update_position is where it starts to get interesting. First, check for a collision during the current frame, then move the particle according to its speed. Once a particle enters phase 3, we can stop recording its position in movement_history.

class Particle:

    [...]

    def update_position(self):
        self.check_collision()
        self.x += self.speed_x
        self.y += self.speed_y
        if self.phase < 3:
            self.movement_history.append((self.x, self.y))
        else:
            self.movement_history.append((None, None))

    [...]

The heaviest math lies in the methods check_collision and get_reflected_point. These methods implement the geometry of the ellipse diagram shown earlier in this post. It’s difficult to walk through every line of code so you’ll have to scroll down a little to follow along.

The if and elif conditions check the distance from a particle to the collision point and focus, respectively.

If a particle is very close to its collision point, we enter the first block. Set the phase attribute to 2 to indicate that we’ll now be traveling toward the focus. Since our frame rate isn’t infinite, we have to cheat a little and nudge the particle to its precise collision point. Then record position in movement_history as usual.

The tangent line slope (derivative), which you can think of as a flat surface where the particle will bounce, is defined this way:

The normal line is perpendicular to the tangent line. To calculate its slope, do the negative inverse, e.g. 4 ➞ -¼.

At this point, we want to reflect the particle’s trajectory over the normal line. Take another look at the ellipse diagram to refresh.

We could reflect any arbitrary point found in movement_history, but let’s instead use the left focus, from which every particle is fired. To reflect a point over a line in the general form ax + by + c = 0, use this formula:

We know the particle is located exactly at the collision point and we know it will travel toward this reflected point, so we can calculate its new angle using atan2. Use sin and cos to get the new x and y speed components. Then the collision is over and it’s back to normal. Each frame updates the particle’s position and nothing notable will happen until it impacts the focus.

When a particle does impact the focus, it enters phase 3 and it should disappear from the plot. I chose to relocate it to the point (999, 999). Out of sight, out of mind.

class Particle:

    [...]

    def check_collision(self):
        if self.phase == 1 and sqrt((self.collision_point[1] - self.y)**2 + (self.collision_point[0] - self.x)**2) < self.speed:
            self.phase = 2
            self.x = self.collision_point[0]
            self.y = self.collision_point[1]
            self.movement_history.append((self.x, self.y))
            tangent_line_slope = -(ellipse.B ** 2 / ellipse.A ** 2) * (self.x / self.y)
            normal_line_slope = -1 / tangent_line_slope
            normal_line_intercept = self.y - normal_line_slope * self.x
            trajectory_point = self.get_reflected_point(normal_line_slope, -1, normal_line_intercept)
            trajectory_angle = atan2(trajectory_point[1] - self.y, trajectory_point[0] - self.x)
            self.speed_x = self.speed * cos(trajectory_angle)
            self.speed_y = self.speed * sin(trajectory_angle)
        elif self.phase == 2 and sqrt((ellipse.focus_y_right - self.y)**2 + (ellipse.focus_x_right - self.x)**2) < self.speed:
            self.phase = 3
            self.x = 999
            self.y = 999

    def get_reflected_point(self, a, b, c):
        point_x, point_y = ellipse.focus_x_left, ellipse.focus_y_left
        d = (a * point_x + b * point_y + c) / (a**2 + b**2)
        x_reflected = point_x - 2 * a * d
        y_reflected = point_y - 2 * b * d
        return x_reflected, y_reflected

Back down to the main loop… We need to take care of a few Matplotlib things before finishing. Set ticks and window limits as a function of A and B, the ellipse parameters.

for frame_counter in range(250):

    [...]

    ticks = range(-ellipse.A, ellipse.A + 1)
    window_limits = (-ellipse.A * 1.1, ellipse.A * 1.1)
    ax.set(xticks=ticks,
           yticks=ticks,
           xlim=window_limits,
           ylim=window_limits)

    [...]

Let’s try Matplotlib’s built-in LaTeX renderer and place the ellipse equation near the top of the window. Activate rendering by putting $ characters at the beginning and end of the string.

for frame_counter in range(250):

    [...]

    ax.text(x=0, y=ellipse.A * 0.98,
            s=r'$\frac{x^2}{' + f'{ellipse.A}' + r'^2} + \frac{y^2}{'  + f'{ellipse.B}' + r'^2} = 1$',
            size=18, ha="center", va="top")

    [...]

I’m saving each frame inside a “frames” folder. When the script finishes, I’ll open them in GIMP and export an animated GIF.

for frame_counter in range(250):

    [...]
    
    plt.savefig(f"frames/{frame_counter}.png")
    
    [...]

Remember, we’re spawning a new particle every 20 frames. Check if the current frame number exists within the additional_particles dictionary. If it does, instantiate a new Particle object and append it to the list.

Finally, iterate frame_counter. One pass through the for loop is complete.

for frame_counter in range(250):

    [...]

    if frame_counter in additional_particles:
        particle_list.append(Particle(initial_slope=additional_particles[frame_counter], speed=0.1))

    frame_counter += 1

When the script finishes running, you’ll have 250 images inside the “frames” folder. As I said, I’ll use GIMP to save them as an animated GIF, but other online and offline tools will work just as well.

---

3. The output.

There it is. It was a bit of a marathon to get here. Every particle deflects off the ellipse wall and impacts the other focus. You can play around with the slopes in additional_particles and see that the values don’t matter. The end result is the same. It’s worth reiterating: this behavior isn’t programmed into the particles. They’re designed to bounce off the wall and follow their natural trajectory.

If you check the movement_history of the particles, you’ll find they all took 100 frames to reach their destination. It’s analogous to the parabola behavior illustrated earlier. All particles travel the same distance regardless of the direction they’re fired. You can use the Pythagorean Theorem to predict exactly how far particles will travel:

Or more intuitively, imagine a particle fired directly to the right.

We don’t encounter ellipses as often in everyday life but I think they’re just as interesting as parabolas. Elliptical pool tables are one (niche) example of the reflective property. It’s also showcased in whispering galleries, where sound is carried across a room seemingly by magic. Now we know how they work.

---

Download the Matplotlib style.

Full code:

from numpy import sqrt, arange
from math import sin, cos, atan2
import matplotlib.pyplot as plt


class Ellipse:

    def __init__(self, a, b):
        self.A = a
        self.B = b
        self.x = arange(-self.A, self.A + 0.0001, 0.0001)
        self.y_lower, self.y_upper = self.get_y_values()
        self.focus_x_left, self.focus_x_right = self.get_focus_values()
        self.focus_y_left, self.focus_y_right = 0, 0

    def get_y_values(self):
        y_lower = -self.B * sqrt(1 - self.x**2 / self.A**2)
        y_upper = self.B * sqrt(1 - self.x**2 / self.A**2)
        return y_lower, y_upper

    def get_focus_values(self):
        focus_x_left = -sqrt(self.A**2 - self.B**2)
        focus_x_right = sqrt(self.A**2 - self.B**2)
        return focus_x_left, focus_x_right


class Particle:

    def __init__(self, initial_slope, speed):
        self.x = ellipse.focus_x_left
        self.y = ellipse.focus_y_left
        self.movement_history = []
        self.slope = initial_slope
        self.collision_point = self.get_collision_point()
        self.speed = speed
        self.speed_x, self.speed_y = self.get_speed_components()
        self.phase = 1

    def get_collision_point(self):
        y_intercept = self.y - self.slope * self.x
        x_val = (sqrt(ellipse.A**4 * ellipse.B**2 * self.slope**2 + ellipse.A**2 * ellipse.B**4 - ellipse.A**2 * ellipse.B**2 * y_intercept**2) - ellipse.A**2 * y_intercept * self.slope) / (ellipse.A**2 * self.slope**2 + ellipse.B**2)
        y_val = self.slope * x_val + y_intercept
        return (x_val, y_val)

    def get_speed_components(self):
        angle = atan2(self.collision_point[1] - self.y, self.collision_point[0] - self.x)
        speed_x = self.speed * cos(angle)
        speed_y = self.speed * sin(angle)
        return speed_x, speed_y

    def update_position(self):
        self.check_collision()
        self.x += self.speed_x
        self.y += self.speed_y
        if self.phase < 3:
            self.movement_history.append((self.x, self.y))
        else:
            self.movement_history.append((None, None))

    def check_collision(self):
        if self.phase == 1 and sqrt((self.collision_point[1] - self.y)**2 + (self.collision_point[0] - self.x)**2) < self.speed:
            self.phase = 2
            self.x = self.collision_point[0]
            self.y = self.collision_point[1]
            self.movement_history.append((self.x, self.y))
            tangent_line_slope = -(ellipse.B ** 2 / ellipse.A ** 2) * (self.x / self.y)
            normal_line_slope = -1 / tangent_line_slope
            normal_line_intercept = self.y - normal_line_slope * self.x
            trajectory_point = self.get_reflected_point(normal_line_slope, -1, normal_line_intercept)
            trajectory_angle = atan2(trajectory_point[1] - self.y, trajectory_point[0] - self.x)
            self.speed_x = self.speed * cos(trajectory_angle)
            self.speed_y = self.speed * sin(trajectory_angle)
        elif self.phase == 2 and sqrt((ellipse.focus_y_right - self.y)**2 + (ellipse.focus_x_right - self.x)**2) < self.speed:
            self.phase = 3
            self.x = 999
            self.y = 999

    def get_reflected_point(self, a, b, c):
        point_x, point_y = ellipse.focus_x_left, ellipse.focus_y_left
        d = (a * point_x + b * point_y + c) / (a**2 + b**2)
        x_reflected = point_x - 2 * a * d
        y_reflected = point_y - 2 * b * d
        return x_reflected, y_reflected


ellipse = Ellipse(a=5, b=3)

particle = Particle(initial_slope=-10, speed=0.1)

particle_list = [particle]

additional_particles = {20: -1, 40: -0.5, 60: -0.1, 80: 0.2, 100: 0.4, 120: 1.2, 140: 7}

plt.style.use("wollen_ellipse.mplstyle")

for frame_counter in range(250):

    print(frame_counter)

    fig, ax = plt.subplots()

    ax.plot(ellipse.x, ellipse.y_upper, color="#444", linewidth=2.0, zorder=1)
    ax.plot(ellipse.x, ellipse.y_lower, color="#444", linewidth=2.0, zorder=1)

    ax.scatter([ellipse.focus_x_left, ellipse.focus_x_right],
               [ellipse.focus_y_left, ellipse.focus_y_right],
               color="None", s=80, edgecolor="#444", linewidth=1.8, zorder=1)

    for particle in particle_list:

        particle.update_position()

        ax.scatter([particle.x], [particle.y], color="#4C4", s=50, edgecolor="#444", linewidth=1.0, zorder=3)

        ax.plot([pos[0] for pos in particle.movement_history[-10:]],
                [pos[1] for pos in particle.movement_history[-10:]],
                color="#4C4", linewidth=1.5, zorder=2)

    ticks = range(-ellipse.A, ellipse.A + 1)
    window_limits = (-ellipse.A * 1.1, ellipse.A * 1.1)
    ax.set(xticks=ticks,
           yticks=ticks,
           xlim=window_limits,
           ylim=window_limits)

    ax.text(x=0, y=ellipse.A * 0.98,
            s=r'$\frac{x^2}{' + f'{ellipse.A}' + r'^2} + \frac{y^2}{'  + f'{ellipse.B}' + r'^2} = 1$',
            size=18, ha="center", va="top")

    plt.savefig(f"frames/{frame_counter}.png")

    if frame_counter in additional_particles:
        particle_list.append(Particle(initial_slope=additional_particles[frame_counter], speed=0.1))

    frame_counter += 1