from numpy import sqrt,arccos,pi,linspace

def rangex(co):          
    si = sqrt(1.-co*co)  # sin(theta)
    tmp = 1.+co   
    return si +2.*co*  \
        (tmp*si +sqrt(tmp*(2.-co*co*tmp)) ) 

'''
Program carousel-range.py

At which angle of deflection should a physical 
pendulum released at the top be released/broken 
to achieve maximum horizontal distance from the 
pendulum's axis? How far will it fly w/o friction?

x  = max range in units of pendulum rad. R 
co = cos(theta)  <- search on uniform mesh
'''
x = 0.
N = 1000
dc = 0.05/N    #  step of the search for cosine

for i in range(N):
    co = 0.75 - dc*i   # search in (0.7,0.75]
    xprev = x
    x = rangex(co)
    # print(i,rangex(c))
    if (xprev > x): break

co_max = co + dc/2  # half-step correction
theta = arccos(co_max)*180/pi

print('\n N = ',N,'  max@step', i,  
  '\n range/R = ', rangex(co_max), 
  '\n cos(theta)', co_max, '\n theta[deg] =',theta,'\n')  

# optional part of the code: plotting
# redefine co from scalar to vector
# and use MatPlotlib
import matplotlib.pyplot as plt
co  = linspace(0.7,0.75,400)
angle = arccos(co)*180/pi
x = rangex(co)  # function accepts array arg!
plt.plot(angle,x)
plt.title(r'x($\theta$)')
plt.xlabel(r'$\theta$')
plt.ylabel("X")
plt.scatter(theta,rangex(co_max),s=9,
                color=(1,0,0),linewidth=2)
plt.show()