import femm
import numpy as np
import matplotlib.pyplot as plt

femm.openfemm()

def capacite(L,e,d,R1=100):
    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)
    U = 1
    femm.ei_addconductorprop('v0',U/2,0,1)
    femm.ei_addconductorprop('v1',-U/2,0,1)
    femm.ei_drawrectangle(-d/2-e,-L/2,-d/2,L/2)
    femm.ei_drawrectangle(d/2,-L/2,d/2+e,L/2)
    femm.ei_addblocklabel(0,0)
    femm.ei_selectlabel(0,0)
    femm.ei_setblockprop('air',0,1,0)
    femm.ei_clearselected()
    femm.ei_addblocklabel(-d/2-e/2,0)
    femm.ei_selectlabel(-d/2-e/2,0)
    femm.ei_setblockprop('metal',0,1,0)
    femm.ei_clearselected()
    femm.ei_addblocklabel(d/2+e/2,0)
    femm.ei_selectlabel(d/2+e/2,0)
    femm.ei_setblockprop('metal',0,1,0)
    femm.ei_clearselected()
    femm.ei_selectsegment(-d/2,0)
    femm.ei_selectsegment(-d/2-e,0)
    femm.ei_selectsegment(-d/2-e/2,L/2)
    femm.ei_selectsegment(-d/2-e/2,-L/2)
    femm.ei_setsegmentprop('<None>',0,1,0,0,'v0')
    femm.ei_clearselected()
    femm.ei_selectsegment(d/2,0)
    femm.ei_selectsegment(d/2+e,0)
    femm.ei_selectsegment(d/2+e/2,L/2)
    femm.ei_selectsegment(d/2+e/2,-L/2)
    femm.ei_setsegmentprop('<None>',0,1,0,0,'v1')
    femm.ei_clearselected()
    femm.ei_makeABC(5,R1,0,0,0)
    femm.ei_saveas('condensateurXY.fee')
    femm.ei_analyze()
    femm.ei_loadsolution()
    femm.eo_selectblock(0,0)
    W = femm.eo_blockintegral(0) # stored energie
    femm.eo_clearblock()
    C = 2*W[0]*1e3/U**2 # capacité par mètre de longueur
    eps0= 8.85e-12
    Cinf = eps0*L/d # capacité par mètre de longueur
    print(C/Cinf)
    c1 = C/Cinf
    N = 1000
    y = np.linspace(-L/2,L/2,N)
    Ex1 = np.zeros(N)
    Ex2 = np.zeros(N)
    eps = 1e-3
    for i in range(N):
        E = femm.eo_gete(-d/2+eps,y[i])
        Ex1[i] = E[0]
        E = femm.eo_gete(-d/2-e-eps,y[i])
        Ex2[i] = E[0]
    Q = eps0*(Ex1-Ex2).sum()*(y[1]-y[0])*1e-3
    N = 100
    x = np.linspace(-d/2-e,-d/2,N)
    Ey1 = np.zeros(N)
    Ey2 = np.zeros(N)
    for i in range(N):
        E = femm.eo_gete(x[i],L/2+eps)
        Ey1[i] = E[1]
        E = femm.eo_gete(x[i],-L/2-eps)
        Ey2[i] = E[1]
    Q += eps0*(Ey1-Ey2).sum()*(x[1]-x[0])*1e-3
    C = Q/U
    c2 = C/Cinf
    return c1,c2


L_list = [4,6,8,10,12,14,16,18,22,24,26,30,36,40,45,50] # largeurs des plaques
e = 1 # épaisseur des plaques
d = 2 # distance entre les plaques
C1_list = []
C2_list = []
filename = "capaciteXY-e=1-d=2.txt"
for L in L_list:
    c1,c2 = capacite(L,e,d)
    C1_list.append(c1)
    C2_list.append(c2)
C1 = np.array(C1_list)
C2 = np.array(C2_list)
L = np.array(L_list)
np.savetxt(filename,np.array([L,C1,C2]).T)

plt.figure()
plt.plot(L,C1,"ro",label="C par energie")
plt.plot(L,C2,"bo",label="C par charge")
plt.grid()
plt.xlabel("L")
plt.ylabel("C/Cinf")
plt.show()
plt.xlim(0,20)
plt.ylim(0,2)
plt.legend(loc="upper right")
plt.show()
femm.closefemm()