import numpy as np
from matplotlib.pyplot import *

N = 20
u=np.linspace(2,10,N)
r = np.zeros(N)
theta = np.linspace(0,2*np.pi,1000)
for i in range(N):
    df = (u[i]**2+u[i]*np.cos(theta))/(1+u[i]**2+2*u[i]*np.cos(theta))
    r[i] = df.sum()*(theta[1]-theta[0])/(2*np.pi)

figure() 
plot(u,r)
grid()
xlabel("u")
ylabel("r")
ylim(0,2)
                 

[Y,F] = np.loadtxt("force-R=5.txt",unpack=True)
R=5
figure(figsize=(12,6))
plot(Y/R,F)
grid()
xlabel("Y/R",fontsize=16)
ylabel("Fy",fontsize=16)

lamb = 1e-14
eps0 = 8.85e-12
f = -lamb**2/(2*np.pi*eps0*Y*1e-3)
plot(Y/R,f,"k--")