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

femm.openfemm()

def capacite(R,e,d,R1=100):
    femm.newdocument(1)
    femm.ei_probdef('millimeters', 'axi', 1.0e-8, 1, 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(0,d/2,R,d/2+e)
    femm.ei_drawrectangle(0,-d/2-e,R,-d/2)
    femm.ei_addblocklabel(1,0)
    femm.ei_selectlabel(1,0)
    femm.ei_setblockprop('air',0,1,0)
    femm.ei_clearselected()
    femm.ei_addblocklabel(R/2,-d/2-e/2)
    femm.ei_selectlabel(R/2,-d/2-e/2)
    femm.ei_setblockprop('metal',0,1,0)
    femm.ei_clearselected()
    femm.ei_addblocklabel(R/2,d/2+e/2)
    femm.ei_selectlabel(R/2,d/2+e/2)
    femm.ei_setblockprop('metal',0,1,0)
    femm.ei_clearselected()
    femm.ei_selectsegment(R/2,-d/2)
    femm.ei_selectsegment(R/2,-d/2-e)
    femm.ei_selectsegment(R,-d/2-e/2)
    automesh = 1
    femm.ei_setsegmentprop('<None>',0,automesh,0,0,'v0')
    femm.ei_clearselected()
    femm.ei_selectsegment(R/2,d/2)
    femm.ei_selectsegment(R/2,d/2+e)
    femm.ei_selectsegment(R,d/2+e/2)
    femm.ei_setsegmentprop('<None>',0,automesh,0,0,'v1')
    femm.ei_clearselected()                 
    femm.ei_makeABC(5,R1,0,0,0)
    femm.ei_saveas('condensateurRZ.fee')
    femm.ei_analyze()
    femm.ei_loadsolution()
    femm.eo_selectblock(0,0)
    W = femm.eo_blockintegral(0) # énergie 
    femm.eo_clearblock()       
    C = 2*W[0]/U**2 # capacité 
    eps0= 8.85e-12
    Cinf = eps0*np.pi*R**2/d*1e-3
    c1 = C/Cinf
    print(c1)
    N = 1000
    r = np.linspace(0,R,N)
    Ez1 = np.zeros(N)
    Ez2 = np.zeros(N)
    eps = 1e-2
    for i in range(N):
        E = femm.eo_gete(r[i],d/2-eps)
        Ez1[i] = E[1]
        E = femm.eo_gete(r[i],d/2+e+eps)
        Ez2[i] = E[1]
    Q = eps0*((Ez1-Ez2)*r).sum()*2*np.pi*(r[1]-r[0])*1e-6
    N = 100
    z = np.linspace(-d/2-e,-d/2,N)
    Er = np.zeros(N)
    for i in range(N):
        E = femm.eo_gete(R+eps,z[i])
        Er[i] = E[0]
    Q += eps0*Er.sum()*(z[1]-z[0])*2*np.pi*R*1e-6
    C = Q/U
    c2 = C/Cinf
    return c1,c2

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

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