Simulating a small circle rolling on another big circle

rolling_circle.py

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib.animation import FuncAnimation

tessellatePieces = 300
assert tessellatePieces % 2 == 0

bigCircleAngles = np.linspace(0.0, 2.0 * np.pi, tessellatePieces, endpoint=False)
smallCircleAngles = np.linspace(0.0, 2.0 * np.pi, tessellatePieces // 2, endpoint=False)

def get_points_on_circle(angles, radius):
    x = radius * np.sin(angles)
    y = radius * np.cos(angles)
    return np.stack([x, y], axis=0)

bigCirclePoints = get_points_on_circle(bigCircleAngles, 2.0)
smallCirclePoints = get_points_on_circle(smallCircleAngles, 1.0)

def get_rotation_matrix(theta):
    return np.asarray(
        [
            [np.cos(theta), -np.sin(theta)],
            [np.sin(theta), np.cos(theta)]
        ],
        np.float
    )


def get_circle_point_by_ind(circlePoints, ind):
    numPoints = circlePoints.shape[1]
    a, b = divmod(ind, numPoints)
    return circlePoints[:, b]


def get_shift_rotation(bigCircIstcInd, smallCircItscInd, smallCirclePointsRef=smallCirclePoints):
    vecBigCirc = get_circle_point_by_ind(bigCirclePoints, bigCircIstcInd + 1) - \
                 get_circle_point_by_ind(bigCirclePoints, bigCircIstcInd)

    vecSmallCirc = get_circle_point_by_ind(smallCirclePointsRef, smallCircItscInd - 1) - \
                 get_circle_point_by_ind(smallCirclePointsRef, smallCircItscInd)

    assert abs(np.linalg.norm(vecSmallCirc) - np.linalg.norm(vecBigCirc)) < 1.0e4

    angleCos = vecBigCirc @ vecSmallCirc / (np.linalg.norm(vecBigCirc) * np.linalg.norm(vecSmallCirc))
    angle = -np.arccos(angleCos)

    shift = get_circle_point_by_ind(bigCirclePoints, bigCircIstcInd) - \
            get_circle_point_by_ind(smallCirclePointsRef, smallCircItscInd)

    return shift, get_rotation_matrix(angle)


def get_transformed_small_circle(itscInd, shift, rotation, smallCirclePointsRef=smallCirclePoints):
    scPoints = smallCirclePointsRef - \
               get_circle_point_by_ind(smallCirclePointsRef, itscInd)[:, np.newaxis]
    scPoints = rotation @ scPoints
    scPoints = scPoints + get_circle_point_by_ind(smallCirclePointsRef, itscInd)[:, np.newaxis]
    scPoints = scPoints + shift[:, np.newaxis]
    return scPoints

# solve for an initial position
halfInd = smallCirclePoints.shape[1] // 2

smallCircleItscInd = halfInd
bigCircleIstcInd = 0
shift, rotation = get_shift_rotation(bigCircleIstcInd, smallCircleItscInd)

# the initial small circle
tfSmallCircle = get_transformed_small_circle(smallCircleItscInd, shift, rotation)

# the roll motion
def roll():
    global tfSmallCircle, shift, rotation, smallCircleItscInd, bigCircleIstcInd
    _shift, _rotation = get_shift_rotation(bigCircleIstcInd + 1, smallCircleItscInd - 1, tfSmallCircle)
    _tfSmallCircle = get_transformed_small_circle(smallCircleItscInd - 1, _shift, _rotation, tfSmallCircle)
    oldSmallCircle = tfSmallCircle.copy()
    # update states
    shift = _shift
    rotation = _rotation
    tfSmallCircle = _tfSmallCircle
    bigCircleIstcInd += 1
    smallCircleItscInd -= 1
    return oldSmallCircle

fig, ax = plt.subplots()
bc, = plt.plot(bigCirclePoints[0, :], bigCirclePoints[1, :], 'b')
sc1, = plt.plot([], [], 'b')
sc2, = plt.plot([], [], 'r')

def ani_init():
    ax.set_xlim(-4.5, 4.5)
    ax.set_ylim(-4.5, 4.5)
    ax.set_aspect(1.0)

def ani_update(frame):
    redLen = 10
    circlePoints = roll()
    sc2.set_data(circlePoints[0, :redLen], circlePoints[1, :redLen])
    sc1.set_data(circlePoints[0, redLen:], circlePoints[1, redLen:])

ani = FuncAnimation(fig, ani_update, frames=np.arange(tessellatePieces), init_func=ani_init, interval=20)

Writer = animation.writers['ffmpeg']
writer = Writer(fps=30, metadata=dict(artist='Alan Xiang'), bitrate=800)
ani.save('circle.mp4', writer)

Detailed description can be found here.