Following the method suggested here, I've written a Python script that generates random points along the perimeter of an ellipse. I wanted to compare the "naive" method (which fails to produce an even distribution) with the correct method, which works like this:

1. Numerically compute P, the perimeter of the ellipse.
2. Pick a random number s on [0,P].
3. Determine the angle t associated with a randomly chosen number s on [0,P].
4. Compute the x- and y-coordinates associated with t from the parametric form of the ellipse.

Here is my the code I am hoping to receive feedback on. I ran it through flake8 already.

import numpy as np
import matplotlib.pyplot as plt


def ellipse_arc(a, b, theta, n):
    """Cumulative arc length of ellipse with given dimensions"""

    # Divide the interval [0 , theta] into n steps at regular angles
    t = np.linspace(0, theta, n)

    # Using parametric form of ellipse, compute ellipse coord for each t
    x, y = np.array([a * np.cos(t), b * np.sin(t)])

    # Compute vector distance between each successive point
    x_diffs, y_diffs = x[1:] - x[:-1], y[1:] - y[:-1]

    cumulative_distance = [0]
    c = 0

    # Iterate over the vector distances, cumulating the full arc
    for xd, yd in zip(x_diffs, y_diffs):
        c += np.sqrt(xd**2 + yd**2)
        cumulative_distance.append(c)
    cumulative_distance = np.array(cumulative_distance)

    # Return theta-values, distance cumulated at each theta,
    # and total arc length for convenience
    return t, cumulative_distance, c


def theta_from_arc_length_constructor(a, b, theta=2*np.pi, n=100):
    """
    Inverse arc length function: constructs a function that returns the
    angle associated with a given cumulative arc length for given ellipse."""

    # Get arc length data for this ellipse
    t, cumulative_distance, total_distance = ellipse_arc(a, b, theta, n)

    # Construct the function
    def f(s):
        assert np.all(s <= total_distance), "s out of range"
        # Can invert through interpolation since monotonic increasing
        return np.interp(s, cumulative_distance, t)

    # return f and its domain
    return f, total_distance


def rand_ellipse(a=2, b=0.5, size=50, precision=100):
    """
    Returns uniformly distributed random points from perimeter of ellipse.
    """
    theta_from_arc_length, domain = theta_from_arc_length_constructor(a, b, theta=2*np.pi, n=precision)
    s = np.random.rand(size) * domain
    t = theta_from_arc_length(s)
    x, y = np.array([a * np.cos(t), b * np.sin(t)])
    return x, y


def rand_ellipse_bad(a, b, n):
    """
    Incorrect method of generating points evenly spaced along ellipse perimeter.
    Points cluster around major axis.
    """
    t = np.random.rand(n) * 2 * np.pi
    return np.array([a * np.cos(t), b * np.sin(t)])

And some test visualizations:

np.random.seed(4987)

x1, y1 = rand_ellipse_bad(2, .5, 1000)
x2, y2 = rand_ellipse(2, .5, 1000, 1000)

fig, ax = plt.subplots(2, 1, figsize=(13, 7), sharex=True, sharey=True)
fig.suptitle('Generating random points on perimeter of ellipse', size=18)
ax[0].set_aspect('equal')
ax[1].set_aspect('equal')
ax[0].scatter(x1, y1, marker="+", alpha=0.5, color="crimson")
ax[1].scatter(x2, y2, marker="+", alpha=0.5, color="forestgreen")
ax[0].set_title("Bad method: Points clustered along major axis")
ax[1].set_title("Correct method: Evenly distributed points")

# Plot arc length as function of theta
theta_from_arc_length, domain = theta_from_arc_length_constructor(2, .5, theta=2*np.pi, n=100)
s_plot = np.linspace(0, domain, 100)
t_plot = theta_from_arc_length(s_plot)

fig, ax = plt.subplots(figsize=(7,7), sharex=True, sharey=True)
ax.plot(t_plot, s_plot)
ax.set_xlabel(r'$\theta$')
ax.set_ylabel(r'cumulative arc length')

I would appreciate a general review, but I also have a few specific questions:

1. How are my comments? In my code, I know many of my comments are explaining "how" rather than "why," and my inclination is simply to remove them, since I think I have made my variable and function names clear. But I would appreciate an example of how a comment like # Compute vector distance between each successive point could be rewritten in "why" terms.
2. Have I computed the arc length in the most efficient way possible? I started by generating the points in two lists of x- and y-coordinates, then iterated over these lists to get the distance made at each step, cumulating those distances as I went around.
3. Is it there anything unconventional or inefficient about using theta_from_arc_length_constructor to create the inverse arc-length function? My thinking was that I need to evaluate the inverse of the arc-length function for every s, so I should go ahead and "prepare" this function rather than recalculate the total arc length each time. But doesn't return np.interp(s, cumulative_distance, t) mean that numpy has to do the interpolation every time I call f ? Or does it get "automatically" vectorized because I feed it s as an array later on when creating the graphs?