I am trying to parametrize the path of a brachistochrone cycloid. Most examples show the final position at the origin, but I want to generalize the path. I have tried adapting an approach similar to the one found on the scipython blog. While it is technically accurate to call my generated path a brachistochrone cycloid, it does not have the correct concavity.

I am not sure what I did wrong. Can someone help me identify my mistake?
The output path was obtained by running the following code for a minimal working example.
`import numpy as np
from scipy.optimize import newton
from scipy.integrate import quad
import matplotlib.pyplot as plt
class BrachistochroneCycloidPath():
def __init__(self, initial_position, final_position, g):
super().__init__()
self.initial_position = initial_position
self.final_position = final_position
self.xi = initial_position[0]
self.yi = initial_position[1]
self.xf = final_position[0]
self.yf = final_position[1]
self.g = g
self._initial_theta = None
self._final_theta = None
self._radius = None
self._arc_length = None
self._elapsed_duration = None
@property
def initial_theta(self):
return self._initial_theta
@property
def final_theta(self):
return self._final_theta
@property
def radius(self):
return self._radius
@property
def arc_length(self):
return self._arc_length
@property
def elapsed_duration(self):
return self._elapsed_duration
def initialize_theta_and_radius(self, x0):
initial_theta = 0
f = lambda theta : (self.yf - self.yi) / (self.xf - self.xi) - (1 - np.cos(theta)) / (theta - np.sin(theta))
final_theta = newton(
f,
x0=x0)
radius = (self.yi - self.yf) / (1 - np.cos(final_theta))
self._initial_theta = initial_theta
self._final_theta = final_theta
self._radius = radius
def initialize_arc_length(self):
integral = quad(
self.ds,
self.initial_theta,
self.final_theta)
arc_length = abs(
integral[0])
self._arc_length = arc_length
def initialize_elapsed_duration(self):
f = lambda theta : self.ds(theta) / self.v(theta)
integral = quad(
f,
self.initial_theta,
self.final_theta)
elapsed_duration = abs(
integral[0])
self._elapsed_duration = elapsed_duration
def x(self, theta):
x = self.radius * (theta - np.sin(theta)) + self.xf
return x
def y(self, theta):
y = self.radius * (1 - np.cos(theta)) + self.yf
return y
def dx(self, theta):
dx = self.radius * (1 - np.cos(theta))
return dx
def dy(self, theta):
dy = self.radius * np.sin(theta)
return dy
def ds(self, theta):
# dy_dx = self.dy(theta) / self.dx(theta)
# ds = np.sqrt(
# 1 + np.square(dy_dx))
dx = self.dx(theta)
dy = self.dy(theta)
ds = np.sqrt(
np.square(dx) + np.square(dy))
return ds
def v(self, theta):
v = np.sqrt(2 * self.g * self.y(theta))
return v
if __name__ == "__main__":
## initialize inputs
initial_position = (1, 10)
final_position = (10, 1)
# x0 = np.pi / 2 # RuntimeError
x0 = -1 * np.pi / 2
g = 9.8
## initialize path
brachistochrone_cycloid_path = BrachistochroneCycloidPath(
initial_position=initial_position,
final_position=final_position,
g=g)
brachistochrone_cycloid_path.initialize_theta_and_radius(
x0=x0)
brachistochrone_cycloid_path.initialize_arc_length()
brachistochrone_cycloid_path.initialize_elapsed_duration()
## print-check arc-length and elapsed duration
print("\n\tArc-Length [m]:\n{}\n".format(
brachistochrone_cycloid_path.arc_length))
print("\n\tElapsed Duration [s]:\n{}\n".format(
brachistochrone_cycloid_path.elapsed_duration))
## view
theta = np.linspace(
brachistochrone_cycloid_path.initial_theta, ## 0
brachistochrone_cycloid_path.final_theta, ## found by newton method
100)
path_x = brachistochrone_cycloid_path.x(
theta=theta)
path_y = brachistochrone_cycloid_path.y(
theta=theta)
path_label = "Path of Brachistochrone Cycloid (radius = {:,.4})".format(
brachistochrone_cycloid_path.radius)
positions_x = [
initial_position[0],
final_position[0]]
positions_y = [
initial_position[1],
final_position[1]]
positions_label = "Initial/Final Position"
fig, ax = plt.subplots()
ax.plot(
path_x,
path_y,
color="darkorange",
label=path_label)
ax.scatter(
positions_x,
positions_y,
color="black",
label=positions_label)
fig.subplots_adjust(
bottom=0.2)
fig.legend(
mode="expand",
loc="lower center",
ncol=2)
ax.grid(
True)
plt.show()
plt.close()`