import numpy as np
import matplotlib.pyplot as plt

h = 2.5
R = 2.
delta = 0.01

def derivative(Y):
    return np.array([Y[1], ( 1. + Y[1] ** 2 ) / Y[0]]) 

def soap_profile(theta):
    z = - 0.5 * h
    Y = np.array([R, theta])
    Lz = [ z ]
    Lr = [ R ]
    while z < 0.5 * h:
        k1 = delta * derivative(Y)
        k2 = delta * derivative(Y + 0.5 * k1)
        k3 = delta * derivative(Y + 0.5 * k2)
        k4 = delta * derivative(Y + k3)
        Y += ( k1 + k2 + k2 + k3 + k3 + k4 ) / 6.
        z += delta
        Lz.append(z)
        Lr.append(Y[0])
    array_z = np.array(Lz)
    array_rho = np.array(Lr)
    return array_z, array_rho

def deviation_end_radius(theta):
    z, r = soap_profile(theta)
    return ( r[-1] - R ) ** 2
    
import scipy.optimize; res = scipy.optimize.minimize(lambda THETA:deviation_end_radius(THETA[0]), np.array([0.])); theta_opt = res.x[0]

z, r = soap_profile(theta_opt)
plt.plot(r, z, label = 'soap film')
plt.plot(-r, z, label = 'soap film')
plt.plot(np.array([-R, R]), np.array([h / 2., h / 2.]), label = 'supporting ring')
plt.plot(np.array([-R, R]), np.array([- h / 2., - h / 2.]), label = 'supporting ring')
plt.ylabel('z (cm)')
plt.xlabel('x (cm)')
plt.legend()
plt.title('Projection of the soap film in one plane containing the z-axis')
plt.show()