import femm
import numpy as np
import matplotlib.pyplot as plt
femm.openfemm()

def energie(h,R,lamb=1e-14,R1=200):
    femm.newdocument(1)
    femm.ei_probdef('millimeters', 'planar', 1.0e-8, 0, 30);
    femm.ei_addmaterial('air',1,1,0)
    femm.ei_addmaterial('metal',1,1,0)
    femm.ei_addconductorprop('v1',0,lamb,0)
    femm.ei_addconductorprop('v2',0,-lamb,0)
    femm.ei_drawarc(R,h,-R,h,180,1)
    femm.ei_drawarc(-R,h,R,h,180,1)
    femm.ei_drawarc(R,-h,-R,-h,180,1)
    femm.ei_drawarc(-R,-h,R,-h,180,1)
    femm.ei_addblocklabel(0,h)
    femm.ei_addblocklabel(0,-h)
    femm.ei_selectlabel(0,h)
    femm.ei_selectlabel(0,-h)
    femm.ei_setblockprop('metal',0,1,0)
    femm.ei_clearselected()
    femm.ei_addblocklabel(0,0)
    femm.ei_selectlabel(0,0)
    femm.ei_setblockprop('air',0,1,0)
    femm.ei_clearselected()
    femm.ei_selectarcsegment(0,h-R)
    femm.ei_selectarcsegment(0,h+R)
    femm.ei_setarcsegmentprop(180,'<None>',0,0,'v1')
    femm.ei_clearselected()
    femm.ei_selectarcsegment(0,-h-R)
    femm.ei_selectarcsegment(0,-h+R)
    femm.ei_setarcsegmentprop(180,'<None>',0,0,'v2')
    femm.ei_clearselected()
    femm.ei_makeABC(5,R1,0,0,0)
    femm.ei_zoomnatural()
    femm.ei_saveas('deuxCylindres.fee')
    femm.ei_analyze()
    femm.ei_loadsolution()
    femm.eo_selectblock(0,0)
    W = femm.eo_blockintegral(0) # énergie par mm de longueur
    femm.eo_clearblock()
    return W[0]

R = 5
N = 100
Y = np.linspace(11,60,N)
deltaY = Y[1]-Y[0]
W = np.zeros(N)
for i in range(N):
    W[i] = energie(Y[i]/2,R)

np.savetxt("energie-R=5.txt",np.array([Y,W]).T)
plt.figure()
plt.plot(Y,W)
plt.grid()
plt.xlabel("Y")
plt.ylabel("W")

F = np.zeros(N-1)
for i in range(N-1):
    F[i] = -(W[i+1]-W[i])/deltaY
Y = Y[0:N-1]
plt.figure()
plt.plot(Y,F)
plt.grid()
plt.xlabel("Y")
plt.ylabel("F")
np.savetxt("force-R=5.txt",np.array([Y,F]).T)

plt.show()
femm.closefemm()