Question: I have two sine waves, as pictured below. Sometimes they are identical waves, just shifted in the y-axis. Other times the sine waves have different period / amplitude etc.  I need to find the maximum radius circle that can fit between the two lines everywhere along the lines (discretized by the resolution that I generated the lines). What is the best way to do this? Is there a faster way than iteratively creating a circle tangent to the bottom curve and stepping up in radius until intersection with the upper curve, at every point?

Answer: I'm not sure if this a the complete answer, but it may contain a few useful ideas. (There is some math, but I think the code is also important, so I think SO is a reasonable place for this question.)   There may be multiple solutions to the last equation, but in practice scipy.optimize.newton seems to converge to the circle of minimum radius. Also, it's vectorized, so you can solve for several circles at once. import numpy as np from scipy.optimize import newton import matplotlib.pyplot as plt from matplotlib import colors as mcolors colors = list(mcolors.TABLEAU_COLORS.values()) # amplitude, angular frequency, phase, and y-offsets params1 = 1, 1, 0, 0 params2 = 0.5, 0.5, 1, 2 # Sinusoid def y(x, params): A, w, p, y = params return A*np.cos(w*x + p) + y # Derivative of Sinusoid def dy(x, params): A, w, p, y = params return -w*A*np.sin(w*x + p) # Angle of Sinusoid Surface Normal def phi(x, params): return np.arctan2(-1, dy(x, params)) # Radius of circle required given y-coordinates def ry(y1, y2, phi1, phi2): return (y2 - y1)/(np.sin(phi1) + np.sin(phi2)) # Radius of circle required given x-coordinates def rx(x1, x2, phi1, phi2): return (x2 - x1)/(np.cos(phi1) + np.cos(phi2)) # At the solution, the radius of the circle given y-coordinates # must agree with the radius of the circle given x-coordinates def f(x2, x1): y1, y2 = y(x1, params1), y(x2, params2) phi1, phi2 = phi(x1, params1), phi(x2, params2) return rx(x1, x2, phi1, phi2) - ry(y1, y2, phi1, phi2) # Coordinate of center of circle def x0y0(x1, r): y1 = y(x1, params1) phi1 = phi(x1, params1) return (x1 + r*np.cos(phi1), y1 + r*np.sin(phi1)) # Plot Sinusoids t_plot = np.linspace(0, 15, 300) y1 = y(t_plot, params1) y2 = y(t_plot, params2) plt.plot(t_plot, y1, 'k', t_plot, y2, 'k') ax = plt.gca() plt.ylim(-2, 4) plt.axis('equal') # Coordinates on sinusoid 1 x1 = np.asarray([1, 4.5, 7, 11, 13]) y1 = y(x1, params1) # X-coordinates on sinusoid 2 x2 = newton(f, x1, args=(x1,)) y2 = y(x2, params2) # Circle radii and centers phi1 = phi(x1, params1) phi2 = phi(x2, params2) r = rx(x1, x2, phi1, phi2) x0, y0 = x0y0(x1, r) for i in range(len(x1)): color = colors[i] # Mark points on curve plt.plot(x1[i], y1[i], color, marker='*') plt.plot(x2[i], y2[i], color, marker='*') # Plot circles circle = plt.Circle((x0[i], y0[i]), r[i], color=color, fill=False) ax.add_patch(circle)   Two caveats:  scipy.optimize.newton is not guaranteed to converge, and even if it does, it is not guaranteed to converge to the circle with minimum radius (unless it can be proven otherwise...). If it is not reliable for your particular parameters but you like the approach, consider posting a more specific question about how to reliably solve the relevant equations numerically. It may just be a matter of solving for one circle and using that as the guess for the next circle. This code doesn't consider the possibility of the circle intersecting with a second point on the lower curve before it touches the upper curve. I couldn't tell from the problem statement if this was required because the question mentioned "stepping up in radius until intersection with the upper curve". This is easy enough to check, though, due to symmetry of the convex portion of the sinusoid.