Plot an osculating circle for a 5th order polynomial together with intersecting normals to the curve

osculating_circle.py

import math
from matplotlib.patches import Circle
import matplotlib.pyplot as plt
import numpy as np

ts = np.linspace(-1.25, 1.25, 500)
dt = ts[1] - ts[0]
x = np.poly1d([1, 0])
y = np.poly1d([0.1, -0.5, 0.2, 1, -1, 0])

u = 0.1255

f, ax = plt.subplots(figsize=(6, 6))
ax.set_xlim([-5, 5])
ax.set_ylim([-5, 5])
curve, = ax.plot(x(ts[100:]), y(ts[100:]), lw=2.5, color='#1f77b4')
ax.text(x(ts[110])+0.02, y(ts[110])+0.02, 'C', fontsize=20, color='#1f77b4')
point, = ax.plot([x(u)], [y(u)], 'o', ms=6, color='#ff7f0e', zorder=200)
ax.text(x(u)-0.12, y(u)-0.15, 'P', fontsize=20, color='#ff7f0e')


def osculating_circle():
    x1 = x.deriv()
    x2 = x1.deriv()
    y1 = y.deriv()
    y2 = y1.deriv()
    numerator = x1(u)**2 + y1(u)**2
    denominator = x1(u)*y2(u) - x2(u)*y1(u)
    r = abs(numerator ** 1.5 / denominator)
    kx = x(u) - y1(u) * numerator / denominator
    ky = y(u) + x1(u) * numerator / denominator
    circle = ax.add_artist(Circle([kx, ky], r, facecolor=None, edgecolor='#ff7f0e', fill=None, lw=2, zorder=100))
    center, = ax.plot([kx], [ky], 'o', ms=6, color='#ff7f0e', zorder=50)
    ax.text(kx + 0.05, ky - 0.1, 'M', fontsize=20, color='#ff7f0e')
    return circle, center


def tangent(t, nt=60, **kwargs):
    x1 = x.deriv()
    y1 = y.deriv()
    tangent_x = np.poly1d([x1(t), x(t) - t*x1(t)])
    tangent_y = np.poly1d([y1(t), y(t) - t*y1(t)])
    t_region = [t - nt*dt, t + nt*dt]
    return ax.plot(tangent_x(t_region), tangent_y(t_region), color='black', lw=1.0, **kwargs)[0]


def normal(t, length, style, **kwargs):
    x1 = x.deriv()
    y1 = y.deriv()
    x2 = x1.deriv()
    y2 = y1.deriv()
    normal_x = np.poly1d([math.copysign(y1(t), y2(t)), x(t) - t*math.copysign(y1(t), y2(t))])
    normal_y = np.poly1d([math.copysign(x1(t), x2(t)), y(t) - t*math.copysign(x1(t), x2(t))])
    t_region = [t, t + length / math.sqrt(x1(t)**2 + y1(t)**2)]
    if y2(t) < 0:
        t_region = t_region[::-1]
    return ax.plot(normal_x(t_region), normal_y(t_region), style, color='black', lw=1.0, **kwargs)[0]


circle = osculating_circle()[0]

u_tangent = tangent(u, zorder=200)
ax.text(u_tangent.get_xdata()[-1]-0.10, u_tangent.get_ydata()[-1]-0.12, 'T', fontsize=20, color='black')
u_normal = normal(u, 2*circle.radius, style='-', zorder=75)
ax.text(u_normal.get_xdata()[-1]-0.2, u_normal.get_ydata()[-1]-0.1, 'N', fontsize=20, color='black')

dns = [-40, -20, 20, 80]  # step-distance from `u`.
normals = [normal(u + dn*dt, circle.radius, style='--') for dn in dns]

xdata, ydata = u_normal.get_xdata(), u_normal.get_ydata()
u_slope = (ydata[1] - ydata[0]) / (xdata[1] - xdata[0])
u_offset = ydata[0] - u_slope*xdata[0]
for j, n in enumerate(normals):
    xdata, ydata = n.get_xdata(), n.get_ydata()
    slope = (ydata[1] - ydata[0]) / (xdata[1] - xdata[0])
    offset = ydata[0] - slope*xdata[0]
    x_intersect = (u_offset - offset) / (slope - u_slope)
    y_intersect = u_slope * x_intersect + u_offset
    ax.plot([x_intersect], [y_intersect], 'o', color='black', ms=3)
    ax.text(x_intersect+0.02, y_intersect-0.02, rf'$\rm S_{j+1}$', fontsize=10, color='black')
    n.set_xdata([xdata[0], x_intersect])
    n.set_ydata([ydata[0], y_intersect])

for j, (n, (dx, dy)) in enumerate(zip(normals, [(-0.03, 0.06), (-0.042, 0.055), (0.035, 0.0), (0.015, 0.025)])):
    xx, yy = n.get_xdata()[0], n.get_ydata()[0]
    ax.text(xx+dx, yy+dy, rf'$\rm P_{j+1}$', fontsize=10, color='black')

ax.set_xlim([-0.75, 1.65])
ax.set_ylim([-0.6, 1.8])
ax.grid()
ax.set_axis_off()
f.savefig('schmiegkreis.png', bbox_inches='tight', pad_inches=0, dpi=600, frameon=False, transparent=False)
f.savefig('schmiegkreis.svg', bbox_inches='tight', pad_inches=0, dpi=600, frameon=False, transparent=False)
plt.show()