import femm
import numpy as np
import matplotlib.pyplot as plt   
                
femm.openfemm()
femm.newdocument(1)
# 0 : magnetic, 1 : electrostatic, 2 ; heat flow, 3 : current flow
femm.ei_probdef('millimeters', 'planar', 1.0e-8, 1, 30);
# unit, type ('planar' ou 'axi'), precision, depth, minangle
                
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)
                  
L = 8
e = 1
d = 2
R1 = 100        
                  
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)
automesh = 1
femm.ei_setsegmentprop('<None>',0,automesh,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,automesh,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_hidedensityplot()
numcontours = 30
Vmin = -0.5
Vmax = 0.5
femm.eo_showcontourplot(numcontours,Vmin,Vmax)
                    
femm.eo_selectblock(0,0)
W = femm.eo_blockintegral(0) # énergie par mm de longueur
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 = %f"%(C/Cinf))
                     
N = 1000
x = np.linspace(-20,20,N)
Ex = np.zeros(N)
for i in range(N):
    E = femm.eo_gete(x[i],0)
    Ex[i] = E[0]
plt.figure()
plt.plot(x,Ex)
plt.title("y=0")
plt.grid()
plt.xlabel("x")
plt.ylabel("Ex")
plt.savefig("../../../../figures/sciphys/elecmag/condensateur3/Eaxe-L=8-e=1-d=2.png")
                      
print("E inf = %f"%(U/(d*1e-3)))        
                      
N = 1000
y = np.linspace(-L/2,L/2,N)
Ex1 = np.zeros(N)
Ex2 = np.zeros(N)
eps = 1e-3 # décalage par rapport au bord du conducteur
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
plt.figure()
plt.plot(y,Ex1,label="face interne")
plt.plot(y,Ex2,label="face externe")
plt.grid()
plt.xlabel("y")
plt.ylabel("Ex")
plt.legend(loc="upper right")
plt.savefig("../../../../figures/sciphys/elecmag/condensateur3/Eplaque-L=8-e=1-d=2.png")               
                      
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
plt.figure()
plt.plot(x,Ey1,label="bord y = L/2")
plt.plot(x,Ey2,label="bord y = -L/2")
plt.grid()
plt.legend(loc="upper right")
plt.xlabel("x")
plt.ylabel("Ey")
plt.savefig("../../../../figures/sciphys/elecmag/condensateur3/Eplaque2-L=8-e=1-d=2.png")                  
                      
C = Q/U
print("C/Cinf = %f"%(C/Cinf))
                       
plt.show()
femm.closefemm()