Example of plotting an Archimedean spiral, its osculating circles at a selection of tangent points, and the evolute that converges asymptotically to a circle

archimedean_spirals_evolute.py

import numpy as np
import matplotlib.pyplot as plt

# Define the polar equation r = 2*theta
def polar_to_cartesian(theta):
    r = 2 * theta
    x = r * np.cos(theta)
    y = r * np.sin(theta)
    return x, y

# Calculate the first and second derivatives of x and y with respect to theta
def derivatives(theta):
    r = 2 * theta
    dr_dtheta = 2
    x = r * np.cos(theta)
    y = r * np.sin(theta)
    dx_dtheta = dr_dtheta * np.cos(theta) - r * np.sin(theta)
    dy_dtheta = dr_dtheta * np.sin(theta) + r * np.cos(theta)
    d2x_dtheta2 = -2 * np.sin(theta) - r * np.cos(theta) - 2 * np.sin(theta)
    d2y_dtheta2 = 2 * np.cos(theta) - r * np.sin(theta) + 2 * np.cos(theta)
    return dx_dtheta, dy_dtheta, d2x_dtheta2, d2y_dtheta2

# Calculate the curvature and radius of curvature
def curvature_and_radius(theta):
    dx_dtheta, dy_dtheta, d2x_dtheta2, d2y_dtheta2 = derivatives(theta)
    numerator = np.abs(dx_dtheta * d2y_dtheta2 - dy_dtheta * d2x_dtheta2)
    denominator = (dx_dtheta**2 + dy_dtheta**2)**(3/2)
    curvature = numerator / denominator
    radius = 1 / curvature if curvature != 0 else np.inf
    return curvature, radius

# Calculate the normal vector
def normal_vector(dx_dtheta, dy_dtheta):
    magnitude = np.sqrt(dx_dtheta**2 + dy_dtheta**2)
    nx = -dy_dtheta / magnitude
    ny = dx_dtheta / magnitude
    return nx, ny

# Calculate the center of the osculating circle
def osculating_center(x, y, nx, ny, radius):
    center_x = x + radius * nx
    center_y = y + radius * ny
    return center_x, center_y

max_theta = 2*np.pi
# Generate theta values for the spiral
theta_values = np.linspace(0, max_theta, 1000)
x_values, y_values = polar_to_cartesian(theta_values)

# Calculate the osculating circles
num_circles = 60
theta_circle = np.linspace(0, max_theta, num_circles)
circle_centers = []
circle_radii = []
normal_lines = []
for theta in theta_circle:
    x, y = polar_to_cartesian(theta)
    dx_dtheta, dy_dtheta, d2x_dtheta2, d2y_dtheta2 = derivatives(theta)
    curvature, radius = curvature_and_radius(theta)
    nx, ny = normal_vector(dx_dtheta, dy_dtheta)
    center_x, center_y = osculating_center(x, y, nx, ny, radius)
    circle_centers.append((center_x, center_y))
    circle_radii.append(radius)
    normal_lines.append(((x, center_x), (y, center_y)))

# Asymptotic circle (evolute limit)
asymptotic_radius = 2
asymptotic_center = (0, 0)

# Calculate the point where the spiral crosses the asymptotic circle
theta_cross = 1  # From solving 2*theta = 2
x_cross, y_cross = polar_to_cartesian(theta_cross)


# Plot the spiral, osculating circles, the asymptotic circle, and the crossing point
plt.figure(figsize=(8, 8))
plt.plot(x_values, y_values, label='Archimedean Spiral', color='yellow')
for endpoints in normal_lines:
    plt.plot(endpoints[0],endpoints[1],color='mediumorchid', alpha=0.5)
for center, radius in zip(circle_centers, circle_radii):
    circle = plt.Circle(center, radius, color='cyan', alpha=0.5, fill=False)
    plt.gca().add_patch(circle)
    plt.scatter(*center, color='mediumorchid',alpha=0.5,s=20)


# Plot the asymptotic circle
asymptotic_circle = plt.Circle(asymptotic_center, asymptotic_radius, color='red', fill=False, linestyle='--', label='Asymptotic Circle')
plt.gca().add_patch(asymptotic_circle)


# Plot the crossing point
plt.scatter(x_cross, y_cross, color='blue',alpha=0.5, s=60, label='Crossing Point')

plt.xlabel('x')
plt.ylabel('y')
plt.title('Archimedean Spiral r=2*θ, Osculating Circles, Evolute')
plt.axis('equal')
plt.legend()
plt.grid(True)
plt.show()

Sample plot: