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', 'axi', 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)

R = 8
e = 1
d = 2
R1 = 100

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_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 
femm.eo_clearblock()       
                    
C = 2*W[0]/U**2 # capacité 
eps0= 8.85e-12
Cinf = eps0*np.pi*R**2/d*1e-3 
print("C/Cinf = %f"%(C/Cinf))

  

         
                
N = 1000
z = np.linspace(-20,20,N)
Ez = np.zeros(N)
for i in range(N):
    E = femm.eo_gete(0,z[i])
    Ez[i] = E[1]
plt.figure()
plt.plot(z,Ez)
plt.grid()
plt.xlabel("z")
plt.ylabel("Ez")
plt.savefig("../../../../figures/sciphys/elecmag/condensateur3/Eaxe-R=8-e=1-d=2.png")
                
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
plt.figure()
plt.plot(r,Ez1,label="face interne")
plt.plot(r,Ez2,label="face externe")
plt.grid()
plt.xlabel("r")
plt.ylabel("Ez")
plt.legend(loc="upper left")
plt.savefig("../../../../figures/sciphys/elecmag/condensateur3/Efaces-R=8-e=1-d=2.png")
                
N = 300
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 
plt.figure()
plt.plot(z,Er,label="bord r=R")
plt.grid()
plt.xlabel("z")
plt.ylabel("Er")
plt.legend(loc="upper right")
plt.savefig("../../../../figures/sciphys/elecmag/condensateur3/Ebord-R=8-e=1-d=2.png")
                
C = Q/U
print("C/Cinf = %f"%(C/Cinf))
                
plt.figure()
plt.show()
femm.closefemm()