from matplotlib.pyplot import *
import numpy as np
import poisson.main
                    

def condensateurXY(p,n,largeur,distance,epaisseur,filename,niter=50,ncycles=15):
    levels = n-p+1 
    laplace=poisson.main.Poisson(p,p,levels,1.0,1.0,x0=-0.5,y0=-0.5)
    laplace.laplacien()
    laplace.dirichlet_borders(0.0)
    xc = 2**(p-1) # centre du domaine
    yc = 2**(p-1)
    U = 1
    laplace.dirichlet_polygon(xc-distance//2,yc-largeur//2,[0,-epaisseur,0,epaisseur],[largeur,0,-largeur,0],U/2)
    laplace.dirichlet_polygon(xc+distance//2,yc-largeur//2,[0,epaisseur,0,-epaisseur],[largeur,0,-largeur,0],-U/2)
    # sommets du conducteur de potentiel U/2 sur le maillage fin, noeuds au dessus de la surface
    x1,y1 = (xc-distance//2)*2**(n-p),(yc-largeur//2)*2**(n-p)
    x2,y2 = x1,(yc+largeur//2)*2**(n-p)+1
    x3,y3 = (xc-distance//2-epaisseur)*2**(n-p)-1,y2
    x4,y4 = x3,y1
    result = laplace.multigrid_full_norm(ncycles,levels,niter,2,1)
    V = laplace.get_array()
    Ex=-laplace.get_derivX()
    Ey=-laplace.get_derivY()
    mask = laplace.get_mask()
    extent = laplace.get_extent()
    x=laplace.get_x()
    dx = x[1]-x[0]
    y = laplace.get_y()
    dy = y[1]-y[0]
    # calcul de la capacité à partir du champ sur le conducteur
    Ex_faceInt = Ex[y1:y2,x1+1]
    Ex_faceExt = Ex[y1:y2,x3-1]
    Q = (Ex_faceInt.sum()-Ex_faceExt.sum())*dy
    Ey_bordHaut = Ey[y2+1,x3:x1]
    Ey_bordBas = Ey[y1-1,x3:x1]
    Q += (Ey_bordHaut-Ey_bordBas).sum()*dx
    C = (Q/largeur*distance) # capacité par rapport à C infini
    # calcul de la capacité à partir de l'énergie
    C2 = (Ex*Ex+Ey*Ey).sum()*dx*dy
    C2 = C2/largeur*distance
    # calcul du champ E infini
    Einf = U/(distance/2**p)    
    if filename: np.savez(filename,V,Ex,Ey,mask,extent,x,y,np.array([C,C2,Einf,y1,y2,x1,x3]))
    return C,C2,Einf
                    

def readXY(filename):
    global V,Ex,Ey,mask,extent,x,y,C,C2,Einf,y1,y2,x1,x3
    npzfiles = np.load(filename)
    V = npzfiles['arr_0']
    Ex = npzfiles['arr_1']
    Ey = npzfiles['arr_2']
    mask = npzfiles['arr_3']
    extent = npzfiles['arr_4']
    x = npzfiles['arr_5']
    y = npzfiles['arr_6']
    A = npzfiles['arr_7']
    C,C2,Einf,y1,y2,x1,x3 = A[0],A[1],A[2],A[3],A[4],A[5],A[6]
    y1 = int(y1)
    y2 = int(y2)
    x1 = int(x1)
    x3 = int(x3)
                    

n=11
p=7
distance = 2 # distance entre les plaques
largeur = 6 # largeur des plaques
epaisseur = 1 #épaisseur des plaques
#condensateurXY(p,n,largeur,distance,epaisseur,"condensateurXY-d=2-L=6-e=1.npz")
                    
 
readXY("condensateurXY-d=2-L=6-e=1.npz")
figure(figsize=(8,8))
contour(V,20,extent=extent)
imshow(~mask,extent=extent,alpha=0.5, cmap=cm.gray)
xlabel("x",fontsize=16)
ylabel("y",fontsize=16)
                    

figure(figsize=(8,8))
contour(V,10,extent=extent)
imshow(~mask,extent=extent,alpha=0.5, cmap=cm.gray)
xlabel("x",fontsize=16)
ylabel("y",fontsize=16)
size = largeur/2**p
xlim(-size,size)
ylim(-size,size)

                    

figure(figsize=(12,6))
plot(x,Ex[2**(n-1),:],'.')
xlabel('x',fontsize=16) 
ylabel('Ex',fontsize=16)
grid()
plot([-0.5,0.5],[Einf,Einf],"r--")
                    
 
figure(figsize=(12,6))
Ex_faceInt = Ex[y1:y2,x1+1]
Ex_faceExt = Ex[y1:y2,x3-1]
plot(y[y1:y2],Ex_faceInt,label="face interne")
plot(y[y1:y2],Ex_faceExt,label="face externe")
legend(loc="upper right")
xlabel("y",fontsize=16)
ylabel("Ex",fontsize=16)
ylim(-100,100)
grid()                  
                      

figure(figsize=(12,6))
Ey_bordHaut = Ey[y2+1,x3:x1]
Ey_bordBas = Ey[y1-1,x3:x1]
plot(x[x3:x1],Ey_bordHaut,label="bord y = L/2")
plot(x[x3:x1],Ey_bordBas,label="bord y = -L/2")
xlabel("x",fontsize=16)
ylabel("Ey",fontsize=16)
legend(loc="upper right")
ylim(-100,100)
grid()                   

n=11
p=7
distance = 2 # distance entre les plaques
epaisseur = 1 #épaisseur des plaques
filename = "capacite-d=2-e=1-n=11-p=7.txt"
calcul = False
if calcul:
    C_list = []
    C2_list = []
    L_list = [4,6,8,10,12,14,16,18,22,24,26,30,36,40] # largeurs des plaques
    for L in L_list:
       C,C2,Einf = condensateurXY(p,n,L,distance,epaisseur,None)
       C_list.append(C)
       C2_list.append(C2)
    C = np.array(C_list)
    C2 = np.array(C2_list)
    L = np.array(L_list)
    np.savetxt(filename,np.array([L,C,C2]).T)
L,C,C2 = np.loadtxt(filename,unpack=True)
figure()
plot(L,C,"b.",label="C par Charge")
plot(L,C,"b-")
plot(L,C2,"r.",label="C par energie")
plot(L,C2,"r-")
legend(loc="upper right")
grid()
xlabel("L",fontsize=16)
ylabel(r"$C/C_{\infty}$",fontsize=16)
xlim(0,40)
ylim(0,2)
                

n=12
p=8
filename = "capacite-d=2-e=1-n=12-p=8.txt"
calcul = False
if calcul:
    C_list = []
    C2_list = []
    for L in L_list:
        C,C2,Einf = condensateurXY(p,n,L,distance,epaisseur,None)
        C_list.append(C)
        C2_list.append(C2)
    C = np.array(C_list)
    C2 = np.array(C2_list)
    L = np.array(L_list)
    np.savetxt(filename,np.array([L,C,C2]).T)
L,C,C2 = np.loadtxt(filename,unpack=True)
figure()
plot(L,C,"b.",label="C par charge")
plot(L,C,"b-")
plot(L,C2,"r.",label="C par energie")
plot(L,C2,"r-")
legend(loc="upper right")
grid()
xlabel("L",fontsize=16)
ylabel(r"$C/C_{\infty}$",fontsize=16)
xlim(0,40)
ylim(0,2)
                

def condensateurZR(p,n,rayon,distance,epaisseur,filename,niter=50,ncycles=20):
    levels = n-p+1
    laplace = poisson.main.PoissonAxial(p+1,p,levels,1.0,0.5,x0=-0.5,y0=0.0)
    laplace.laplacien()
    laplace.dirichlet_borders(0.0)
    zc = 2**p
    U = 1
    laplace.dirichlet_polygon(zc-distance//2,0,[0,-epaisseur,0,epaisseur],[rayon,0,-rayon,0],U/2)
    laplace.dirichlet_polygon(zc+distance//2,0,[0,epaisseur,0,-epaisseur],[rayon,0,-rayon,0],-U/2)
    # sommets du conducteur de potentiel U/2 sur le maillage fin
    z1,r1 = (zc-distance//2)*2**(n-p),0
    z2,r2 = z1,(rayon)*2**(n-p)+1
    z3,r3 = (zc-distance//2-epaisseur)*2**(n-p)-1,r2
    z4,r4 = z3,r1
    result = laplace.multigrid_full_norm(ncycles,levels,niter,2,1)
    V = laplace.get_array(symetry=True)
    Ez=-laplace.get_derivZ()
    Er=-laplace.get_derivR()
    mask = laplace.get_mask(symetry=True)
    extent = laplace.get_extent(symetry=True)
    z = laplace.get_z()
    dz = z[1]-z[0]
    r = laplace.get_r()
    dr = r[1]-r[0]
    #calcul de la capacité à partir du champ sur le conducteur
    Ez_faceInt = Ez[r1:r2,z1+1]
    Ez_faceExt = Ez[r1:r2,z3-1]
    Q = ((Ez_faceInt-Ez_faceExt)*r[r1:r2]).sum()*2*np.pi*dr
    Er_bord = Er[r2+1,z3:z1]
    Q += Er_bord.sum()*dz*2*np.pi*r[r2]
    e = distance/2**(p+1)
    R = rayon/2**(p+1)
    S = np.pi*R*R
    C = Q/S*e 
    Einf = U/(distance/2**(p+1))    
    if filename: np.savez(filename,V,Ez,Er,mask,extent,z,r,np.array([C,Einf,r1,r2,z1,z3]))
    return C,Einf
                

def readZR(filename):
    global V,Ez,Er,mask,extent,z,r,C,Einf,r1,r2,z1,z3
    npzfiles = np.load(filename)
    V = npzfiles['arr_0']
    Ez = npzfiles['arr_1']
    Er = npzfiles['arr_2']
    mask = npzfiles['arr_3']
    extent = npzfiles['arr_4']
    z = npzfiles['arr_5']
    r = npzfiles['arr_6']
    A = npzfiles['arr_7']
    C,Einf,r1,r2,z1,z3 = A[0],A[1],A[2],A[3],A[4],A[5]
    r1 = int(r1)
    r2 = int(r2)
    z1 = int(z1)
    z3 = int(z3)
                    

n = 10
p = 6
distance = 2 # distance entre les plaques
rayon = 3 # rayon des plaques
epaisseur = 1 #épaisseur des plaques
#C,Einf = condensateurZR(p,n,rayon,distance,epaisseur,"condensateurZR-d=2-r=3-e=1-n=10-p=6.npz")
readZR("condensateurZR-d=2-r=3-e=1-n=10-p=6.npz")
figure(figsize=(8,8))
contour(V,40,extent=extent)
imshow(~mask,extent=extent,alpha=0.5, cmap=cm.gray)
xlabel("z",fontsize=16)
ylabel("r",fontsize=16)
                    

figure(figsize=(8,8))
contour(V,10,extent=extent)
imshow(~mask,extent=extent,alpha=0.5, cmap=cm.gray)
xlabel("x",fontsize=16)
ylabel("y",fontsize=16)
size = 2*rayon/2**p
xlim(-size,size)
ylim(-size,size)
                    

figure(figsize=(12,6))
plot(z,Ez[0,:],'.')
xlabel('z',fontsize=16) 
ylabel('Ez',fontsize=16)
grid()
plot([-0.5,0.5],[Einf,Einf],"r--")
                    

figure(figsize=(12,6))
Ez_faceInt = Ez[r1:r2,z1+1]
Ez_faceExt = Ez[r1:r2,z3-1]
plot(r[r1:r2],Ez_faceInt,label="face interne")
plot(r[r1:r2],Ez_faceExt,label="face externe")
legend(loc="upper right")
xlabel("r",fontsize=16)
ylabel("Ez",fontsize=16)
ylim(-100,100)
grid()                  
                      

figure(figsize=(12,6))
Er_bord = Er[r2+1,z3:z1]
plot(z[z3:z1],Er_bord,label="bord r = R")
xlabel("z",fontsize=16)
ylabel("Er",fontsize=16)
legend(loc="upper right")
ylim(-100,100)
grid()                   

n=10
p=6
distance = 2 # distance entre les plaques
epaisseur = 1 #épaisseur des plaques
filename = "capaciteZR-d=2-e=1-n=10-p=6.txt"
calcul = False
if calcul:
    C_list = []
    R_list = [2,4,6,8,10,12,14,16,18,20,25,30,35,40,45,50] # rayons
    for R in R_list:
       C,Einf = condensateurZR(p,n,R,distance,epaisseur,None)
       C_list.append(C)
    C = np.array(C_list)
    R = np.array(R_list)
    np.savetxt(filename,np.array([R,C]).T)
R,C = np.loadtxt(filename,unpack=True)
figure()
plot(R,C,"b.")
grid()
xlabel("R",fontsize=16)
ylabel(r"$C/C_{\infty}$",fontsize=16)
xlim(0,50)
ylim(0,2)
                

n=11
p=7
distance = 2 # distance entre les plaques
epaisseur = 1 #épaisseur des plaques
filename = "capaciteZR-d=2-e=1-n=11-p=7.txt"
calcul = False
if calcul:
    C_list = []
    for R in R_list:
       C,Einf = condensateurZR(p,n,R,distance,epaisseur,None)
       C_list.append(C)
    C = np.array(C_list)
    R = np.array(R_list)
    np.savetxt(filename,np.array([R,C]).T)
R,C = np.loadtxt(filename,unpack=True)
figure()
plot(R,C,"b.")
grid()
xlabel("R",fontsize=16)
ylabel(r"$C/C_{\infty}$",fontsize=16)
xlim(0,50)
ylim(0,2)